RasmahRasmah

Fluid–Structure Interaction Topology Optimization

Overview

A fluid flowing past a deformable structure both loads and is redirected by it. FSITopOpt optimizes the coupled layout — the fluid flow and the elastic solid — through the Darcy interpolation, minimizing the dissipated power of the monolithic fluid + solid system. It takes both a fluid and a solid material.

using Rasmah
using LinearAlgebra

A coupled fluid–solid cell

Clamp the outer faces, drive the interior with a body force, and let the optimizer distribute fluid vs. solid.

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)
solid = Material(youngs_modulus = 1e3, poisson_ratio = 0.3)

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(FSITopOpt(), m, fluid, solid;
                fixed_u = fixed_u, loads_u = loads_u, fixed_p = [1], gp = [0.0],
                volfrac = 0.5, max_iter = 20)
([1.0, 1.0, 0.3329939071885987, 0.33264320351388094, 1.0, 1.0, 0.47741121794757574, 1.0, 0.0, 0.0  …  0.4770452049109532, 0.47741121794757463, 1.0, 0.0, 0.3326432035138782, 1.0, 1.0, 1.0, 1.0, 0.332993907188598], [0.15673044314129952, 0.14861602042422084, 0.14027767604127137, 0.1294703425455196, 0.12069413986972738, 0.11286486444776542, 0.10697999013338694, 0.10479908479175058, 0.10295771670796966, 0.10136398161709487, 0.10011070149708651, 0.09949703769042856, 0.09903634702804809, 0.09898331889263406, 0.09928322755595967, 0.09842527671342242, 0.09833149759463967, 0.09838894404017659, 0.0983304760843842, 0.09826107088804319])

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: 0.157 → 0.098
decreases: true
fluid fraction ≈ 0.5: 0.499
vtk(m; field = x, fieldname = "fluid density", edges = true)

Next steps