Steady-state heat diffusion in 1D cylindrical mesh

MW 230817, 240503

This example will evaluate the analytic solution for a simple steady-state heat transfer problem, and use PyFVTool to solve the same problem by the finite-volume method.

import numpy as np
from numpy import exp
import matplotlib.pyplot as plt

# explicity import all required routines from pyfvtool
from pyfvtool import CylindricalGrid1D
from pyfvtool import CellVariable
from pyfvtool import transientTerm, diffusionTerm
from pyfvtool import harmonicMean
from pyfvtool import constantSourceTerm
from pyfvtool import solvePDE

Problem statement

The steady-state heat equation in 1D cylindrical geometry (total cylinder radius \(R\)), with an internal heat source (strength \(S\)) and Dirichlet (constant-value) outer boundary. This can model, for example the temperature (temperature \(T\) relative to outer temperature) profile in a electrically resistive wire through which a current passes.

\[0 = k \nabla^2 T + S\]

which, under the circumstances mentioned above, should have solution

\[T(r) = \frac{S}{4k} (R^2 - r^2)\]

Do not forget that we are in cylindrical coordinates and that :math:`nabla^2` has the corresponding form.

System parameters

It is possibly wise to not set all parameters to 1 or even integer values for testing purposes but rather use something else.

R = 2.1
k_val = 3.7 # heat transfer coefficient
S_val = 4.2 # source strength
T_outer = 0.0  # (outer) boundary temperature

Finite-volume solution with PyFVTool

Define 1D cylindrical grid with radius \(R\).

Nr = 50
Lr = R

mesh = CylindricalGrid1D(Nr, Lr)

Create coefficients of diffusion and source terms in the form of CellVariables.

k = CellVariable(mesh, k_val) # heat transfer coefficient

S = CellVariable(mesh, S_val)

Create solution variable T, initialize to arbitrary 0.0 value.

Subsequently, specify a boundary condition for this variable: outer wall will be kept at 0.0. (Dirichlet boundary condition).

T = CellVariable(mesh, 0.0)

# switch the right (=outer) boundary to Dirichlet: fixed temperature
T.BCs.right.a = 0.0
T.BCs.right.b = 1.0
T.BCs.right.c = T_outer

The diffusion term requires the face values of the diffusion coefficient

k_face = harmonicMean(k)

The equation is defined as a list of the different matrix equation terms.

eqnterms = [-diffusionTerm(k_face),
             constantSourceTerm(S)]

Go solve. Since no transient term is part of the equation, this will yield directly the steady-state solution.

solvePDE(T, eqnterms);

Retrieve solution temperature profile from solution CellVariable.

rnum, Tnum = T.plotprofile()

Comparison with the analytic solution

r_an = np.linspace(0, R)
T_an = (S_val/(4*k_val))*(R**2 - r_an**2)

plt.plot(rnum, Tnum, 'ro', label='FVM')
plt.plot(r_an, T_an, 'k', label='analytic')
plt.ylabel('rel. temp. $T$')
plt.xlabel('radial position $r$')
plt.legend(); # semicolon avoids output of '<matplotlib.legend.Legend...>' descriptor

Quantitative testing of agreement between FV result and analytic solution

The tolerance is manually tuned. This should help detect code changes that break good agreement with experiment and theory.

# analytic solution evaluated at cell positions
T_an_rnum = (S_val/(4*k_val))*(R**2 - rnum**2)

# normalized error
norm_err = (Tnum-T_an_rnum)/T_an_rnum.max()

plt.plot(rnum, norm_err, 'o')
plt.ylabel('error relative to maximum $T$')
plt.xlabel('radial position $r$');

print(np.abs(norm_err) < 0.0001)

[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True False False False False False False False False False
 False False False  True]

assert np.all(np.abs(norm_err) < 0.001)