Heat conduction

A cooling fin in one dimension

The equation \(-u'' + 400\,u = 0\) on \((0, 0.05)\) with \(u(0) = 300\) and \(u'(L) + 2u(L) = 0\) is Example 5.3.1 of the book. It shows the three ingredients of every problem — kernels, boundary conditions, and a solve — and the recovery of the secondary variable at the fixed end.

# SPDX-License-Identifier: LGPL-2.1-or-later
"""Reddy, Example 5.3.1: a cooling fin, -u'' + 400 u = 0 on (0, 0.05).

u(0) = 300, and at x = L the fin loses heat to the surroundings:
u'(L) + 2 u(L) = 0.  The example prints the nodal temperatures, the heat flow
into the fin at the fixed end (the secondary variable), and the same solution
computed with the finite element method for comparison.
"""

import dualmesh as dm


def solve(num_elements, method="dmcdm"):
    mesh = dm.generate_line_mesh(start=0.0, end=0.05, num_elements=num_elements)
    problem = dm.Problem(mesh, method=method)
    problem.add_variable("temperature")
    problem.add_kernel("Diffusion", variable="temperature")
    problem.add_kernel("Reaction", variable="temperature", coefficient=400.0)
    problem.add_boundary_condition(
        "DirichletBC", variable="temperature", boundary="left", value=300.0
    )
    problem.add_boundary_condition(
        "RobinBC", variable="temperature", boundary="right", transfer_coefficient=2.0
    )
    problem.solve()
    return problem


if __name__ == "__main__":
    for method in ("dmcdm", "fem"):
        problem = solve(5, method)
        values = problem.values("temperature")
        print(f"{method:6s} u = " + " ".join(f"{v:8.3f}" for v in values))
        print(f"{method:6s} Q(0) = {problem.total_reaction('temperature', 'left'):.1f}")
    print("book (DMCDM): 300.00  257.62  225.59  202.64  187.83  180.57, Q(0) = 4817")

Running it prints

dmcdm  u =  300.000  257.618  225.593  202.637  187.827  180.568
dmcdm  Q(0) = 4817.0
fem    u =  300.000  257.567  225.507  202.527  187.703  180.437

which are the values of Table 5.3.1 of the book. Note how little separates the two methods: only the method argument of dualmesh.Problem.

A bus bar in two dimensions

Example 5.4.3 adds internal heat generation and a convective boundary. The same problem is also available as an input file (Input files).

# Reddy, Example 5.4.3: a bus bar with internal heat generation, fixed
# temperatures on the sides, an insulated bottom, and convection on top.
mesh:
  type: rectangle
  x_min: 0.0
  x_max: 0.1
  y_min: 0.0
  y_max: 0.05
  num_x_elements: 10
  num_y_elements: 5

problem:
  method: dmcdm

variables:
  temperature: {initial_condition: 0.0}

kernels:
  conduction:
    type: HeatConduction
    variable: temperature
    thermal_conductivity: 20.0
  heating:
    type: HeatSource
    variable: temperature
    heat_source: 1.0e6

boundary_conditions:
  left:
    type: DirichletBC
    variable: temperature
    boundary: left
    value: 40.0
  right:
    type: DirichletBC
    variable: temperature
    boundary: right
    value: 10.0
  top:
    type: ConvectiveHeatFluxBC
    variable: temperature
    boundary: top
    heat_transfer_coefficient: 75.0
    ambient_temperature: 0.0

executioner:
  type: steady

outputs:
  vtu: bus_bar.vtu
  reactions: [[temperature, left], [temperature, right]]
  point_values: [[temperature, 0.05, 0.0], [temperature, 0.05, 0.05]]

Things worth noticing

• An insulated boundary needs no boundary condition at all: a zero natural condition is the default, which in the dual mesh method means that the boundary face of the control domain simply contributes nothing.

• total_reaction returns the heat flow through a boundary, computed from the equations that the Dirichlet conditions replaced — not by differentiating the solution afterwards.

• A temperature-dependent conductivity \(k = k_0 (1 + k_1 T)\) needs only temperature_polynomial=[1.0, k1] on the HeatConduction kernel; the problem then solves by Newton’s method with the exact Jacobian, or by direct iteration with nonlinear_solver="picard".