RasmahRasmah

Contact-Constrained Topology Optimization

Overview

Many designs touch other bodies — a compliant gripper on a workpiece, a valve seat, a snap-fit. When contact is part of the objective, the optimizer must account for the non-penetration constraint. Rasmah’s ContactTopOpt minimizes compliance under a unilateral contact obstacle a·x ≤ b (normal a, offset b) enforced with a penalty stiffness kpen, so the optimized layout respects the obstacle while carrying its load.

using Rasmah

A cantilever pressing onto an obstacle

Clamp a cube, pull it downward, and place a rigid floor (normal +z, offset 0) beneath it. The design must brace against the obstacle as well as the clamp.

m = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2)
n = size(m.nodes, 2)
fixed = Int[]
for i in 1:n
    m.nodes[1, i] < 1e-6 && append!(fixed, 3i - 2, 3i - 1, 3i)
end
loads = zeros(3n)
for i in 1:n
    m.nodes[1, i] > 0.999 && (loads[3i] = -1e6)     # pull in −z
end

obstacle = ([0.0, 0.0, 1.0], 0.0)                        # rigid floor z = 0

x, hist = solve(ContactTopOpt(), m, Material(youngs_modulus = 1e9, poisson_ratio = 0.3);
                fixed = fixed, loads = loads, obstacle = obstacle, kpen = 1e4,
                volfrac = 0.5, max_iter = 10)
([0.5810054820018875, 0.7433861733157766, 0.001, 0.001, 0.5140384296283357, 0.963346847762727, 0.46046048444464294, 0.3458227683472827, 0.012174276204474295, 0.4847354349205757  …  1.0, 1.0, 0.8625438777534439, 0.001, 0.40847646071210864, 0.7098226302811973, 1.0, 0.6872456273164467, 0.7649504406441124, 0.001], [2.9688902533179275e6, 2.7541155336179e6, 2.517629867134984e6, 2.2922885597548955e6, 2.0396033941263803e6, 1.7417197317914283e6, 1.5832827476267451e6, 1.4821460342943594e6, 1.3441671521933856e6, 1.1658862360028364e6])

Verification — compliance drops

The objective must decrease, and the density field stays in [0,1]:

println("compliance: ", round(hist[1]; digits = 3), " → ", round(hist[end]; digits = 3))
println("density in [0,1]: ", all(0 .<= x .<= 1))
compliance: 2.968890253e6 → 1.165886236e6
density in [0,1]: true
vtk(m; field = x, fieldname = "density (contact)", edges = true)

Next steps