Tutorial: 2-D Steady Heat Conduction around a Cylinder (Dirichlet)

This tutorial walks through the Tests/2D_Heat_Conduction_Cylinder_Dirichlet test case. It demonstrates how to use an embedded-boundary level-set body with an isothermal (Dirichlet) thermal BC to compute steady-state heat conduction in a 2-D domain containing a circular cylinder — the canonical Laplace-equation problem with an analytic solution.

Required build flags: USE_EB=TRUE, USE_TEMP=TRUE, DIM=2.

Physical problem

A unit square domain \([0,1]^2\) contains a circular cylinder of radius \(r = 0.15\) centred at \((0.5,\,0.5)\). The four domain walls are held at \(T = 0\), and the cylinder surface is isothermal at \(T = 1\). There are no volumetric heat sources, and the initial elastic solid velocity is zero throughout.

The governing equation is the steady-state heat equation (Laplace equation):

\[\nabla^2 T = 0\]

with boundary conditions:

\[\begin{split}T &= 0 \quad \text{on all four domain walls} \\ T &= 1 \quad \text{on the cylinder surface } (r = 0.15)\end{split}\]

In the time-dependent MPM formulation, the simulation is advanced until the transient decays and the solution converges to the steady state. The thermal diffusivity is \(\alpha = k / (\rho\,c_p) = 1\), so the characteristic diffusion time across the domain is \(t_\alpha \sim 1\). The simulation is run to \(t = 0.1\), by which point the solution has relaxed to within plotting accuracy of the Laplace steady state.

Analytical solution

The steady-state temperature in polar coordinates centred on the cylinder axis is:

\[T(r, \theta) = \frac{\ln(R_\text{out}/r)}{\ln(R_\text{out}/r_0)}\]

where \(r_0 = 0.15\) is the cylinder radius and \(R_\text{out}\) is an effective outer radius determined by the Dirichlet condition \(T = 0\) on the square boundary. For a circular outer boundary of radius \(R_\text{out}\), this reduces to:

\[T(r) = \frac{\ln(R_\text{out}/r)}{\ln(R_\text{out}/r_0)}\]

Because the outer boundary is square rather than circular, the exact solution is a 2-D series; however, at large separation the logarithmic profile provides an excellent approximation. The post-processing script in PostProcess/validate.py computes the numerical angular average at each radial band and compares it against a reference profile obtained by high-resolution numerical integration of the Laplace equation on the same geometry.

Material properties

Property	Value	Notes
Density \(\rho\)	1.0	uniform
Specific heat \(c_p\)	1.0	uniform; gives \(\alpha = 1\)
Thermal conductivity \(k\)	1.0	uniform
Volumetric heat source \(\dot{q}\)	0.0	no internal source
Initial temperature \(T_0\)	0.0	uniform; below the cylinder BC

Domain and discretisation

The computational domain is the unit square \([0,1]^2\) discretised with a \(128 \times 128\) background grid (mpm.ncells = 128 128 0). One material point per cell is placed throughout the domain excluding the cylinder interior, giving approximately 15 000 active particles. The time step is fixed at \(\Delta t = 10^{-5}\) and the simulation runs to \(t_\text{final} = 0.1\).

The level-set grid is refined by a factor of 2 (cylinder.ls_refinement = 2), giving an effective resolution of \(256 \times 256\) for the signed-distance field, which accurately represents the circular geometry.

Boundary conditions

Domain walls — all four faces use Dirichlet temperature BCs at \(T = 0\) and slip momentum BCs (no momentum exchange with the walls):

mpm.bc_xlo_mom = slip
mpm.bc_xhi_mom = slip
mpm.bc_ylo_mom = slip
mpm.bc_yhi_mom = slip

mpm.bc_xlo_temp          = dirichlet
mpm.bc_xlo_temp.T_wall   = 0.0
mpm.bc_xhi_temp          = dirichlet
mpm.bc_xhi_temp.T_wall   = 0.0
mpm.bc_ylo_temp          = dirichlet
mpm.bc_ylo_temp.T_wall   = 0.0
mpm.bc_yhi_temp          = dirichlet
mpm.bc_yhi_temp.T_wall   = 0.0

Cylinder surface — isothermal at \(T = 1\) with a slip momentum BC (no drag on the particles from the cylinder):

eb2.body_names                      = cylinder

cylinder.geom_type                  = sphere
cylinder.sphere_radius              = 0.15
cylinder.sphere_center              = 0.5 0.5 0.0
cylinder.sphere_has_fluid_inside    = false
cylinder.ls_refinement              = 2
cylinder.levelset_mom               = slipwall
cylinder.temp_bc_type               = isothermal
cylinder.lset_T_wall                = 1.0

For a full description of the level-set BC parameters see Embedded Boundaries via Level Sets (USE_EB).

Setting up the test case

The particle file and input file are generated from the JSON configuration in PreProcess/config.json. The key fields are the grid resolution, the particles-per-cell count, and the thermal properties:

{
  "grid": { "nx": 128, "ny": 128 },
  "ppc": [1, 1],
  "bodies": [{
    "temperature": {
      "T": 0.0,
      "spheat": 1.0,
      "thermcond": 1.0,
      "heatsrc": 0.0
    }
  }]
}

To (re)generate the particle and input files, run from the test directory:

cd Tests/2D_Heat_Conduction_Cylinder_Dirichlet/
bash Generate_MPs_Inputfile_Generic.sh

This writes mpm_particles.dat (particle positions and thermal properties) and Inputs_2DHeat_Conduction_Cylinder_Dirichlet.inp one level up.

Building

Navigate to the test directory and build with the required flags:

cd Tests/2D_Heat_Conduction_Cylinder_Dirichlet
# create a temporary copy of the GNUmakefile here
make USE_EB=TRUE USE_TEMP=TRUE DIM=2 -j4

The resulting executable is named ExaGOOP2d.<suffix>.ex (suffix depends on the build environment).

Note

USE_EB=TRUE pulls in the AMReX EB2 library. On some systems this significantly increases compile time; using -j4 or higher parallelises the build across compilation units.

Running

From the test directory:

./ExaGOOP2d.<suffix>.ex Inputs_2DHeat_Conduction_Cylinder_Dirichlet.inp

For a parallel run on four MPI ranks:

mpirun -n 4 ./ExaGOOP2d.<suffix>.ex \
    Inputs_2DHeat_Conduction_Cylinder_Dirichlet.inp

ASCII particle snapshots are written at intervals of mpm.write_output_time = 0.01 to the subdirectory 2D_Heat_Conduction_Cylinder_Dirichlet/<solution folder name set in input file>/matpnt*. A checkpoint is written to 2D_Heat_Conduction_Cylinder_Dirichlet/solution folder name set in input file/chk* at the same frequency.

Post-processing and validation

The validation script reads the final particle snapshot, bins particles by radial distance from the cylinder axis, and compares the angularly-averaged temperature profile against the analytical steady-state:

cd Tests/2D_Heat_Conduction_Cylinder_Dirichlet/
python3 PostProcess/validate.py --time 0.1

A passing run prints the RMS error and PASS. The acceptance criterion is RMS \(< 5 \times 10^{-2}\) (the finite-time run and square-vs-circle outer boundary geometry contribute a background error of order \(10^{-2}\)).

The companion script Plot_Temperature.py generates a colour map of the particle temperature field together with a contour overlay:

python Plot_Temperature.py

Expected result

At \(t = 0.1\) the temperature field has relaxed to a radially symmetric profile that rises from \(T = 0\) at the domain walls to \(T = 1\) at the cylinder surface. The gradient is steepest in the gap between the cylinder and the nearest wall. Representative radially-averaged values are:

Radial distance from axis \(r\)	\(T_\text{avg}(r)\)
0.15 (cylinder surface)	1.000
0.20	≈ 0.82
0.30	≈ 0.57
0.40	≈ 0.32
≥ 0.50 (approaching walls)	→ 0.0

The temperature profile is symmetric about the cylinder axis to within numerical discretisation error.