Transport & coupled physics

Overview

The physics chapters so far have covered solids, heat, electricity and magnetism, and the main coupled families. This chapter collects the remaining transport and coupled problems — the analyses where a scalar or a field is carried through a body by a physical process, or where two processes act together. They share the same solve(Physics(), …) pattern: name the physics, supply the material (and, for flow, a fluid), give the boundary conditions.

Four of them follow the pattern directly:

• Bioheat — blood perfusion warming tissue (Pennes' equation);
• Poisson–Nernst–Planck — the electric double layer in an electrolyte;
• magnetohydrodynamics — an electrically conducting fluid in a magnetic field;
• time-domain acoustics — a pressure wave marching through time.

Two more are analysis tags rather than a full solve: the J-integral (energy release rate at a crack tip) and the topological derivative (the sensitivity of a functional to a vanishingly small hole).

Bioheat: blood perfusion

Living tissue is not a passive conductor — blood flow carries heat away. Pennes' bioheat equation adds a perfusion term (the blood warms the tissue toward the arterial temperature $T_a$) and a metabolic source $Q_m$:

\[\nabla\cdot(k\nabla T) + \omega\rho_b c_b (T_a - T) + Q_m = 0 ,\]

where $\omega$ is the perfusion rate, $\rho_b c_b$ the blood heat capacity, and $k$ the tissue conductivity (from the material). With $\omega = 0$ it reduces to heat conduction with a source:

Bioheat

m = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4)
left  = nodes_in_box(m, [0.0, 0.0, 0.0], [0.0, 1.0, 1.0])
right = nodes_in_box(m, [1.0, 0.0, 0.0], [1.0, 1.0, 1.0])
mid   = nodes_in_box(m, [0.5, 0.0, 0.0], [0.5, 1.0, 1.0])

tissue = Material(youngs_modulus = 1e9, poisson_ratio = 0.3, thermal_conductivity = 1.0)

T = solve(Bioheat(), m, tissue; Qm = 8.0, fixed = (left => 0.0, right => 0.0))

125-element QuantityVector{Float64}:
 0.0 K
 0.7499999999999982 K
 0.9999999999999973 K
 0.7499999999999982 K
 0.0 K
 0.0 K
 0.7499999999999979 K
 0.9999999999999968 K
 0.7499999999999979 K
 0.0 K
 ⋮
 0.7499999999999976 K
 0.9999999999999963 K
 0.7499999999999974 K
 0.0 K
 0.0 K
 0.7499999999999979 K
 0.9999999999999968 K
 0.7499999999999979 K
 0.0 K

all(isapprox(T[i], 1.0u"K"; atol = 1e-6) for i in mid.ids)

true

With both faces held at $0$ and a uniform metabolic source $Q_m = 8$, the slab develops a parabolic profile whose mid-plane temperature is $Q_m L^{2}/(8k) = 1$ — exactly what the solve returns.

Poisson–Nernst–Planck: the electric double layer

An electrolyte next to a charged surface forms a double layer: ions screen the charge over a distance called the Debye length. In equilibrium the Poisson equation for the potential couples to the Boltzmann distribution of the ions, giving the Poisson–Boltzmann equation. Rasmah solves the equilibrium potential $\phi$:

NernstPlanck

The material supplies the relative permittivity $\varepsilon$; F, RT, and c0 are the (dimensionless) Faraday constant, thermal energy, and background concentration.

Magnetohydrodynamics: a conducting fluid in a field

An electrically conducting fluid moving through a magnetic field feels a Lorentz force, which feeds back into the flow. The inductionless MHD model treats the field as an imposed background $\mathbf{B}_0$ and solves the Stokes flow with the Lorentz body force added:

Magnetohydrodynamics

The fluid supplies the dynamic viscosity $\mu$ and electrical conductivity $\sigma$; the solve returns both the velocity $\mathbf{u}$ and the pressure $p$.

Time-domain acoustics: the wave equation

Sound is a pressure wave obeying

\[\frac{\partial^{2} p}{\partial t^{2}} = c^{2}\,\nabla^{2} p ,\]

with $c$ the sound speed (from the material). Rasmah marches this second-order system forward in time with the Newmark scheme:

TimeAcoustics

wave_solve(M, K, p0, v0, fixed, g, nsteps, dt, f; β = 0.25, γ = 0.5) -> ps

acoustic = Material(youngs_modulus = 1e9, poisson_ratio = 0.3, density = 1.0)
r = solve(TimeAcoustics(), m, acoustic; nsteps = 5, dt = 1e-3)
(length(r.ps), length(r.times))

(6, 6)

The solve returns the pressure field at every time step and the matching time grid. Rigid (sound-hard) walls are the default; pass fixed for a sound-soft ($p = 0$) or prescribed-pressure boundary.

Fracture: the J-integral

A crack grows when it is energetically favourable to do so. The J-integral (Rice 1968) measures the energy release rate: for a mode-I crack of stress intensity $K$ it is $J = K^{2}/E'$, where $E' = E$ (plane stress) or $E/(1-\nu^{2})$ (plane strain). It is path-independent, so it can be evaluated on a contour far from the tip:

JIntegral

j_integral_contour(K, E, ν; plane=:strain, R=1.0, npts=2000)

mode1_stress(K, x, y) -> [σxx, σyy, σxy]

mode1_displacement(K, μ, ν, x, y; plane = :strain) -> (ux, uy)

round(j_integral_contour(1.0, 200e9, 0.3; plane = :strain); sigdigits = 3)

4.55e-12

round((1 - 0.3^2) / 200e9; sigdigits = 3)

4.55e-12

The contour integral around the exact mode-I field recovers $K^{2}/E'$ to machine precision, confirming the identity.

The topological derivative

Where the J-integral asks should this crack grow?, the topological derivative asks a different question: if I punch an infinitesimal hole here, how much does the objective change? For the linear-elastic compliance the answer is the strain-energy density,

\[D_T J(\mathbf{x}) = \boldsymbol{\sigma} : \boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^{\top}\mathbf{C}\boldsymbol{\varepsilon},\]

so a positive value means removing a small volume of material there raises the compliance. This is the shape-sensitivity that steers topological-derivative optimization:

TopologicalDerivative

topological_derivative_elasticity(m, u, C) -> dT

Next steps

These are the last of the built-in physics families. For the fluid mechanics behind the magnetohydrodynamics solve, see the fluid chapter; for the topology-optimization loops that consume the topological derivative and the J-integral, see the Topology optimization chapter.