Heat Conduction: Conductivity Tree Benchmark

Description

This tutorial demonstrates heat conduction topology optimization using the classic “conductivity tree” benchmark from Bendsøe & Sigmund’s Topology Optimization textbook (§1.3). The problem minimizes thermal compliance (temperature work) subject to a volume constraint.

The setup: distributed heat flux enters the top edge, the bottom edge is held at T = 0 (heat sink), and the sides are insulated. The optimal design forms a branching tree structure that efficiently conducts heat from the top to the bottom boundary — analogous to structural force flow paths.

Setup

using TopOpt

Define the Problem

The HeatTree problem sets up the standard benchmark with heat flux on top and fixed temperature at the bottom:

nels = (80, 40)  # 80×40 rectangular mesh
k = 1.0          # thermal conductivity
q = 1.0          # heat flux on top boundary
problem = HeatTree(Val{:Linear}, nels, (1.0, 1.0), k; q=q)

The boundary conditions:

Top edge: Distributed heat flux q (Neumann BC)
Bottom edge: Fixed temperature T = 0 (Dirichlet BC)
Left/Right edges: Insulated (zero heat flux)

Parameter Settings

V = 0.3       # volume fraction (30% high-conductivity material)
xmin = 1e-3   # minimum density
rmin = 2.0    # density filter radius

x0 = fill(V, getncells(problem))  # start from uniform 30% density
solver = FEASolver(DirectSolver, problem; xmin=xmin, penalty=PowerPenaltyFun(1.0))
comp = ThermalComplianceFun(solver)
filter = DensityFilterFun(solver; rmin=rmin)

The 30% volume constraint forces the optimizer to strategically place high- conductivity material for maximum heat dissipation.

Objective and Constraint

obj = x -> comp(filter(PseudoDensities(x)))  # thermal compliance
volfrac = VolumeFun(solver)
constr = x -> volfrac(filter(PseudoDensities(x))) - V  # volume constraint

Thermal compliance measures the work done by the heat flux against the temperature field: J = QᵀT. Minimizing it drives heat efficiently toward the cold boundary.

Optimization Setup

model = Model(obj)
addvar!(model, zeros(length(x0)), ones(length(x0)))
add_ineq_constraint!(model, constr)
alg = MMA87()
convcriteria = Nonconvex.KKTCriteria()

Continuation Strategy

A fixed penalty often leaves intermediate (gray) densities. We use a continuation strategy that ramps the penalty exponent through 1.0, 2.0, and 3.0, with exponentially decaying KKT tolerances so each subproblem is solved accurately before the penalty increases:

ps = [1.0, 2.0, 3.0]
tols = [1e-2, 1e-3, 1e-4]

x = x0

Run Optimization

Each continuation step:

Set the penalty exponent p
Set the convergence tolerance tol
Solve the subproblem with MMA87
Warm-start the next step from the current solution

for j in eachindex(ps)
    global x, convcriteria, res
    p = ps[j]
    tol = tols[j]
    setpenalty!(solver, p)
    options = MMAOptions(; tol=Tolerance(; kkt=tol), maxiter=1000, convcriteria)
    res = optimize(model, alg, x; options)
    global x = res.minimizer
    println("Step $j: penalty=$p, tol=$tol, obj=$(obj(x))")
end

[ Info:   iter       obj      Δobj  violation  kkt_residual  
[ Info:      0   1.1e+04       Inf   1.7e-16   1.3e+01
[ Info:      1   1.1e+04   4.5e-06   0.0e+00   4.5e-06
[ Info:      2   1.1e+04   1.3e-09   0.0e+00   4.5e-06
Step 1: penalty=1.0, tol=0.01, obj=10641.83572115921
[ Info:   iter       obj      Δobj  violation  kkt_residual  
[ Info:      0   3.5e+04       Inf   0.0e+00   8.4e+01
[ Info:      1   3.5e+04   2.5e-05   0.0e+00   4.5e-06
[ Info:      2   3.5e+04   7.8e-10   0.0e+00   4.5e-06
Step 2: penalty=2.0, tol=0.001, obj=35199.648008018674
[ Info:   iter       obj      Δobj  violation  kkt_residual  
[ Info:      0   1.1e+05       Inf   0.0e+00   4.0e+02
[ Info:      1   1.1e+05   5.5e-03   0.0e+00   5.5e-03
[ Info:      2   1.1e+05   3.2e-09   0.0e+00   5.5e-03
[ Info:      3   1.1e+05   2.2e-09   0.0e+00   5.5e-03
[ Info:      4   1.1e+05   6.7e-04   0.0e+00   4.9e-03
[ Info:      5   1.1e+05   3.7e-03   0.0e+00   8.6e-03
[ Info:      6   1.1e+05   5.3e-03   0.0e+00   3.3e-03
[ Info:      7   1.1e+05   9.4e-04   0.0e+00   2.3e-03
[ Info:      8   1.1e+05   2.3e-03   0.0e+00   4.9e-06
Step 3: penalty=3.0, tol=0.0001, obj=114396.02474086873

The MMA87 optimizer iteratively updates the density field to minimize thermal compliance while satisfying the volume constraint. Ramping the penalty drives the final design toward a crisp 0/1 distribution.

Results

@show obj(res.minimizer)
@show constr(res.minimizer)

obj(res.minimizer) = 114396.02474086873
constr(res.minimizer) = -2.4707458301520546e-12

-2.4707458301520546e-12

The final design should show the characteristic branching “tree” pattern that efficiently conducts heat from the top to the bottom boundary.

Visualization

using CairoMakie
fig = visualize(problem; topology=res.minimizer)

Precompiling packages...
   5921.8 ms  ✓ QuartoNotebookWorkerMakieExt (serial)
  1 dependency successfully precompiled in 6 seconds
Precompiling packages...
   5071.3 ms  ✓ QuartoNotebookWorkerCairoMakieExt (serial)
  1 dependency successfully precompiled in 5 seconds

Figure 1: Conductivity tree: optimal heat conduction paths from top flux to bottom sink

The visualization reveals the tree-like structure — high-conductivity material forms branching paths similar to how structural members carry loads in compliance minimization. This analogy between heat conduction and structural mechanics is a fundamental insight in topology optimization.