RasmahRasmah

Navier–Stokes Topology Optimization

Overview

When inertia matters, the flow is governed by the Navier–Stokes equations rather than Stokes. NavierStokesTopOpt minimizes the dissipated power of steady laminar incompressible flow including the convective term, with the Darcy inverse-permeability interpolation (Gersborg-Hansen et al. 2005) and an optional Reynolds-number continuation that warm-starts the design.

using Rasmah
using LinearAlgebra

A flow cell with inertia

Same setup as the Stokes case, with the convective term active:

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(NavierStokesTopOpt(), m, fluid;
                fixed_u = fixed_u, loads_u = loads_u, fixed_p = [1], gp = [0.0],
                volfrac = 0.5, max_iter = 20)
([0.615428373043191, 0.6396191228702578, 0.5904933725503175, 0.639349901899179, 0.5853544557437249, 0.8620479142130699, 0.4652856210081785, 0.34096289549747855, 0.0, 0.0  …  0.5504610228285669, 0.46524278872772734, 0.34095632899128764, 0.0, 0.6393658589939933, 0.5853700123947728, 0.8620042466390669, 0.6153993475468361, 0.6396153727434832, 0.5904924118779986], [5.446683020883803, 5.255587603379732, 4.994183174977595, 4.734997607787177, 4.498120322562857, 4.284406512931147, 4.154498481514091, 4.053742812092069, 4.112190043415788, 4.10025118735364, 4.066711834878477, 4.084795201398597, 4.11898747090436, 4.081778345122198, 4.042753524939584, 4.023146704244793, 4.014639216724797, 4.00615738741989, 3.9958757405702925, 3.9890450280379603])

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