RasmahRasmah

Stokes Flow Topology Optimization

Overview

Fluid topology optimization shapes a channel to minimize the energy it wastes. In the Darcy (Borrvall–Petersson) approach the design x ∈ [0,1] is the fluid fraction: x=1 is fluid, x=0 is solid, and a penalized inverse-permeability α(x) damps the flow through solid regions. StokesTopOpt minimizes the dissipated power Φ = ∫ (2μ ε:ε + α u·u) of creeping (Stokes) flow subject to a fluid-volume fraction.

using Rasmah
using LinearAlgebra

A flow cell

Clamp all six faces, drive the interior with a body force, and let the optimizer carve the least-dissipative flow path.

m = tetrahedralize_box(1.0, 1.0, 1.0, 3, 3, 3)
U = FESpace(m, ReferenceFE(3, 2); valuetype = VectorValue{3})
fluid = FluidMaterial(dynamic_viscosity = 1.0, density = 1.0)

fixed_u = Int[]
loads_u = zeros(ndofs(U))
for i in 1:U.nnodes
    x = U.nodes[:, i]
    onb = x[1] < 1e-9 || x[1] > 1 - 1e-9 || x[2] < 1e-9 || x[2] > 1 - 1e-9 ||
          x[3] < 1e-9 || x[3] > 1 - 1e-9
    onb && append!(fixed_u, 3i - 2, 3i - 1, 3i)
    onb || (loads_u[3i - 2] = 1.0)
end

x, hist = solve(StokesTopOpt(), m, fluid;
                fixed_u = fixed_u, loads_u = loads_u, fixed_p = [1], gp = [0.0],
                volfrac = 0.5, max_iter = 20)
([0.6329219519540343, 0.6586996091174869, 0.6116879130471055, 0.6532448843826241, 0.6015725748857199, 0.6895218821193521, 0.47076283129221874, 0.34876981885626807, 0.0, 0.0  …  0.5547563727047596, 0.4707628312922176, 0.34876981885626496, 0.0, 0.6532448843826216, 0.6015725748857161, 0.6895218821193576, 0.6329219519540416, 0.6586996091174808, 0.6116879130471023], [5.446683023079157, 5.255587605094582, 4.9941831756239665, 4.734997608513774, 4.498120322675618, 4.284406512676908, 4.154498481310019, 4.0537561995780464, 4.112280131124962, 4.099972289288176, 4.068524203239332, 4.0873880039459, 4.121870955307585, 4.082272469493773, 4.036054330909749, 4.023081931198799, 4.014712464291705, 4.006313959215445, 3.99687017909127, 3.98916979890058])

Verification — the dissipated power drops

println("dissipated power: ", round(hist[1]; digits = 3), " → ", round(hist[end]; digits = 3))
println("decreases: ", hist[end] < hist[1])
println("fluid fraction ≈ 0.5: ", round(sum(x) / length(x); digits = 3))
dissipated power: 5.447 → 3.989
decreases: true
fluid fraction ≈ 0.5: 0.5
vtk(m; field = x, fieldname = "fluid density", edges = true)

Next steps