RasmahRasmah

Heat Conduction and Thermal Topology Optimization

using Rasmah

Overview

Heat conduction is the scalar twin of linear elasticity: the field is a single temperature T, the material supplies a conductivity k, and the assembled operator is the sparse diffusion matrix K = ∫ k ∇Nᵢ·∇Nⱼ. This tutorial walks the steady solve, then reuses the same boundary-condition syntax to run a SIMP thermal-compliance topology optimization with ThermalComplianceTopOpt.

Geometry, face selection, and meshing

The workflow follows the forward FEM exactly: build the geometry, take its exact BRep, select faces before meshing, then mesh.

# 1. the geometry
part = box(1.0, 1.0, 0.2)

# 2. exact boundary representation + face selection
b = brep(part)
left  = faces_in_box(b, [-0.6, -1.0, -1.0], [-0.4, 1.0, 1.0])
right = faces_in_box(b, [ 0.4, -1.0, -1.0], [ 0.6, 1.0, 1.0])

# 3. mesh — the face sets carry over
mesh = tetrahedralize(b, 8)
TetMesh{Matrix{Float64}, Matrix{Int64}, Rasmah.BRepFaceLabels}([-0.5 -0.375 … 0.375 0.5; -0.5 -0.5 … 0.5 0.5; -0.1 -0.1 … 0.1 0.1], [1 1 … 152 152; 2 11 … 233 234; 11 10 … 234 153; 92 92 … 243 243], Rasmah.BRepFaceLabels(FaceSet[FaceSet([49 22 … 28 24; 58 23 … 29 25; 59 32 … 38 34]), FaceSet([184 215 … 186 227; 193 216 … 187 228; 194 225 … 196 237]), FaceSet([1 3 … 7 88; 2 84 … 8 89; 83 85 … 89 170]), FaceSet([73 161 … 75 74; 154 242 … 76 75; 155 243 … 157 156]), FaceSet([64 37 … 19 145; 73 118 … 100 226; 154 127 … 109 235]), FaceSet([135 108 … 126 45; 144 117 … 207 126; 225 198 … 216 135])]))

Visualize the geometry and the mesh:

vtk(b)

Steady heat conduction

Hold the left face at T = 1 and the right at T = 0 — a temperature gradient across the plate:

T = solve(HeatConduction(), mesh, copper;
          fixed = [left => 1.0, right => 0.0])

println("temperature range: ", round.(extrema(T); digits = 4))
temperature range: (0.0 K, 1.0 K)

The temperature falls linearly from 1 to 0 across the plate. Colour the mesh by the nodal field:

vtk(mesh; field = T, fieldname = "temperature")

Thermal compliance minimization (SIMP + MMA)

ThermalComplianceTopOpt minimizes the thermal compliance qᵀT subject to a volume-fraction constraint, using the same region => value pair syntax. Here a unit flux is injected on the left and the right is held cold:

φ, hist = solve(ThermalComplianceTopOpt(), mesh, copper;
                fixed = right => 0.0, loads = left => 1.0,
                volfrac = 0.5, max_iter = 40)

println("initial thermal compliance = ", round(hist[1]; digits = 4))
println("final thermal compliance   = ", round(hist[end]; digits = 4))
println("final volume fraction      = ", round(sum(φ) / length(φ); digits = 4))
initial thermal compliance = 0.004 m^2·kg·K/s^3
final thermal compliance   = 0.004 m^2·kg·K/s^3
final volume fraction      = 0.5

The density concentrates into conductive paths. Plot the final density field:

vtk(mesh; field = φ, fieldname = "thermal density")

Next steps