RasmahRasmah

Bioheat: the Pennes Equation

Overview

Living tissue is heated from within: blood perfusion carries heat away and metabolism generates it. The Pennes bioheat equation (Pennes 1948) adds these two effects to ordinary conduction,

\[-k \nabla^2 T + \omega \rho_b c_b (T - T_a) = Q_m\]

where k is the tissue conductivity, ω the perfusion rate, ρ_b c_b the blood heat capacity, T_a the arterial temperature and Q_m the metabolic heat. The perfusion term is a volumetric source/sink proportional to T − T_a — unlike pure conduction, temperature now decays exponentially with a characteristic length.

using Rasmah

The analytic solution

For a 1D slab with T(0) = T₀ and an insulated far end, the steady solution is

\[T(x) = T_a + (T_0 - T_a)\frac{\cosh(m(L-x))}{\cosh(mL)}, \qquad m = \sqrt{\frac{\omega\rho_b c_b}{k}}\]

k, P, Ta, T0, L = 1.0, 4.0, 0.0, 1.0, 1.0      # P = ω ρ_b c_b
mpar = sqrt(P / k)
Texact(x) = Ta + (T0 - Ta) * cosh(mpar * (L - x)) / cosh(mpar * L)
println("T(L)  = ", round(Texact(L); digits = 4))
println("T(L/2) = ", round(Texact(L / 2); digits = 4))
T(L)  = 0.2658
T(L/2) = 0.4102

FEM solution on a slab

Solve the steady Pennes equation on a [0,L]×[0,1]×[0,1] bar with T = T₀ on the x = 0 face and natural (zero-flux) conditions elsewhere, which reduces to the 1D solution above.

m = tetrahedralize_box(1.0, 1.0, 1.0, 32, 2, 2)
n = size(m.nodes, 2)
fixed = [i for i in 1:n if m.nodes[1, i] < 1e-12]
g = fill(T0, length(fixed))

# ω ρ_b c_b = 4, T_a = 0, no metabolic heat
T = solve(Bioheat(), m, Material(youngs_modulus = 1e3, poisson_ratio = 0.3, thermal_conductivity = k);
          fixed = fixed, g = g, ω = 1.0, ρb = 1.0, cb = 4.0, Ta = 0.0, Qm = 0.0)
297-element QuantityVector{Float64}:
 1.0 K
 0.9420315437516167 K
 0.8876916085385923 K
 0.8367724011256851 K
 0.7890791704175067 K
 0.7444294693676308 K
 0.7026524642407874 K
 0.6635882884516597 K
 0.6270874384069982 K
 0.5930102089653653 K
 ⋮
 0.29690138016821277 K
 0.2887331141256054 K
 0.2816894624576991 K
 0.2757414445409309 K
 0.27086434946382243 K
 0.26703763288048443 K
 0.2642448308504241 K
 0.2624734904272918 K
 0.2617151168529074 K

Verification

Compare the FEM temperature at the x = L end against the closed form:

iL = argmin(abs.(m.nodes[1, :] .- L))
println("T_FEM(L)  = ", round(T[iL]; digits = 4))
println("T_exact(L)= ", round(Texact(L); digits = 4))
T_FEM(L)  = 0.2698 K
T_exact(L)= 0.2658
vtk(m; field = T, fieldname = "temperature T", edges = true)

Next steps