RasmahRasmah

Robust and Stress-Constrained Topology Optimization

Overview

Two practical concerns for a manufacturable design: robustness to manufacturing variation (etching/erosion) and stress limits. Rasmah provides both on top of the density method.

  • Robust (three-field) optimization — the filtered density is projected at three thresholds: dilated (η < 0.5, over-etching), blueprint (η = 0.5) and eroded (η > 0.5, under-etching). The objective is the worst-case (smoothed maximum) of the three compliances, so the design stays stiff under either manufacturing error.
  • Stress-constrained optimization — a p-norm of the element von Mises stress is kept below a limit while minimizing compliance.

The part is a machine mounting foot: a 100 mm cube of steel bolted down at its base and carrying a 100 kN press load on its top face (10 MPa spread over the 100 × 100 mm pad). The design question is the same one an engineer asks of any support: how little material can carry that load, and can it do so without overstressing the steel?

Both use the same shape → mesh → boundary conditions → solve(Physics(), …) workflow.

using Rasmah

# 1. the design space — a 100 mm steel mounting foot, 2. meshed
m  = tetrahedralize_box(0.1, 0.1, 0.1, 4, 4, 4)
bs = box_boundary_sets(m, 0.1, 0.1, 0.1)

# 3. boundary conditions — base bolted down, 10 MPa on the top pad (100 kN)
fixed = dofs(m, bs[:zmin])
loads = bs[:zmax] => [0.0, 0.0, -1.0e7]
FaceSet([111 101 … 119 118; 116 102 … 120 119; 117 107 … 125 124]) => [0.0, 0.0, -1.0e7]

Robust (three-field) optimization

RobustProjection(β, Δη) sets the three Heaviside thresholds; the driver minimizes the log-sum-exp of the three compliances subject to the blueprint volume:

x, hist = solve(RobustTopOpt(), m, steel;
                fixed = fixed, loads = loads, volfrac = 0.5, β = 4.0, max_iter = 20)

println("worst-case compliance ", round(hist[1]; digits = 3), " → ", round(hist[end]; digits = 3), " J")
println("final volume fraction = ", round(sum(x) / length(x); digits = 3))

vtk(m; field = x, fieldname = "density", edges = true)
worst-case compliance 2.407 → 1.367 J
final volume fraction = 0.498

The worst-case compliance drops from 2.41 J to 1.37 J: the design became stiffer under the most adverse of the three manufacturing realizations — and it did so while using exactly the 50% volume budget, because the blueprint (design-intent) realization is the one the constraint is applied to. That is the point of the robust formulation: a foot whose etched or over-etched versions are both still stiff.

Stress-constrained optimization

StressTopOpt adds the p-norm von Mises stress constraint σ_PN ≤ σ_lim (the stress is recovered with the full material on the SIMP displacement). The limit here is 500 MPa — a steel-class value, deliberately set high enough that the optimizer still has freedom; in a real design loop it would be the yield strength divided by the code’s factor of safety:

σ_lim = 5.0e8

xs, hists = solve(StressTopOpt(), m, steel;
                  fixed = fixed, loads = loads, volfrac = 0.5, σ_lim = σ_lim, max_iter = 20)

println("compliance ", round(hists[1]; digits = 3), " → ", round(hists[end]; digits = 3), " J")
println("final volume fraction = ", round(sum(xs) / length(xs); digits = 3))

vtk(m; field = xs, fieldname = "density", edges = true)
compliance 3.821 → 1.073 J
final volume fraction = 0.494

The stress constraint changes the design: instead of the lightest possible structure, the optimizer now adds material wherever the stress p-norm approaches the limit. For reference, the same 100 kN load on the full foot (unoptimized, no density scaling) produces a peak von Mises stress of:

u = solve(Elasticity(), m, steel; fixed = fixed, loads = loads)
σvm = von_mises_stress(m, ustrip(u), steel)

println("max von Mises (full material) = ", round(maximum(σvm) / 1e6; digits = 1), " MPa  (limit ", σ_lim / 1e6, " MPa)")
max von Mises (full material) = 11.0 MPa  (limit 500.0 MPa)

11 MPa against a 500 MPa limit — the full foot is nowhere near yielding, which is exactly why the optimizer can afford to remove half its mass; the stress constraint’s job is to make sure the remaining structure does not push any region past the limit as the density field converges. Tighten σ_lim toward the yield strength and the optimizer must start adding material back where the load path concentrates.

vtk(m; field = σvm, fieldname = "von Mises stress")

Next steps