RasmahRasmah

Phase-Field Structural Topology Optimization

Overview

The density method’s 0/1 design is often replaced by a smooth phase field: the design ρ is evolved by an Allen–Cahn reaction–diffusion equation

\[\frac{\partial \rho}{\partial t} = \kappa \nabla^2 \rho - \frac{dJ}{d\rho}\]

so the compliance gradient drives a diffusive field that naturally produces smooth, length-scale-controlled interfaces. PhaseFieldStructuralTopOpt marches this phase-field flow, sharpening the design as Δt advances.

using Rasmah

A cantilever via phase-field evolution

Clamp a bar, pull on the free end, and let the phase field evolve toward a stiff layout.

m = tetrahedralize_box(2.0, 0.5, 0.5, 12, 3, 3)
n = size(m.nodes, 2)
fixed = Int[]
for i in 1:n
    m.nodes[1, i] ≈ 0.0 && append!(fixed, 3i - 2, 3i - 1, 3i)
end
loadnodes = [i for i in 1:n if m.nodes[1, i] ≈ 2.0]
forces = zeros(3n)
for i in loadnodes
    forces[3i - 2] = 1.0 / length(loadnodes)
end

ρ, hist = solve(PhaseFieldStructuralTopOpt(), m, Material(youngs_modulus = 1e3, poisson_ratio = 0.3);
                fixed = fixed, loads = forces, volfrac = 0.5, p = 3.0, κ = 0.3, Δt = 0.5, max_iter = 100)
([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 0.001, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [0.06454549133000999, 0.06306077685192105, 0.061554521639490795, 0.05849190508694262, 0.05129167709015097, 0.04106451535593886, 0.03322418286617327, 0.026739963435324806, 0.023757142211277657, 0.02230669905536346  …  0.018036633522967346, 0.0180366291302308, 0.018036620590171932, 0.01803660398951211, 0.018036571729030462, 0.018036509069779917, 0.018036387492931003, 0.01803615207088875, 0.01803569795805415, 0.018034828486407715])

Verification — compliance drops

The compliance falls as the phase field sharpens into a discrete layout:

println("compliance: ", round(hist[1]; digits = 3), " → ", round(hist[end]; digits = 3))
println("decreases: ", hist[end] < hist[1])
println("density in [0,1]: ", all(0 .<= ρ .<= 1))
compliance: 0.065 → 0.018
decreases: true
density in [0,1]: true
vtk(m; field = ρ, fieldname = "phase-field density", edges = true)

Next steps