RasmahRasmah

Poroelastic Topology Optimization

Overview

Porous structures — filters, bone scaffolds, sound absorbers — carry fluid as well as solid load. Poroelastic topology optimization minimizes the compliance of the coupled Biot system

\[\begin{bmatrix} K_{uu} & -\alpha\,B \\ 0 & k\,K_{pp} \end{bmatrix} \begin{bmatrix} u \\ p \end{bmatrix} = \begin{bmatrix} f \\ q \end{bmatrix}\]

where u is the solid displacement, p the pore pressure, α the Biot coefficient and k the permeability. Optimizing the density x couples both fields, so the design must balance mechanical stiffness against fluid transport.

using Rasmah

A loaded, pressurised block

Clamp the left face of the solid, apply a solid load on the right, and fix the pore pressure on two opposite faces.

m = tetrahedralize_box(1.0, 1.0, 1.0, 3, 3, 3)
n = size(m.nodes, 2)

fixed_u = Int[]
for i in 1:n
    m.nodes[1, i] < 1e-6 && append!(fixed_u, 3i - 2, 3i - 1, 3i)
end
loads_u = zeros(3n)
for i in 1:n
    m.nodes[1, i] > 0.999 && (loads_u[3i - 2] = 1e6)
end
fixed_p = [i for i in 1:n if m.nodes[2, i] < 1e-6 || m.nodes[2, i] > 0.999]

x, hist = solve(PoroelasticTopOpt(), m, Material(youngs_modulus = 1e9, poisson_ratio = 0.3);
                fixed_u = fixed_u, loads_u = loads_u, fixed_p = fixed_p,
                α = 1.0, k = 1.0, volfrac = 0.5, max_iter = 10)
(x = [0.7007554369247262, 0.5876755494686526, 0.44097362031187426, 0.44097362029638587, 0.587675549464797, 0.7007554369170185, 0.6697200260249141, 0.524767209261687, 0.42721504382365305, 0.4272150438175196  …  0.6439077752548416, 0.6439077752624599, 0.5068517819281975, 0.3908453328518009, 0.42184046412292936, 0.556344937500626, 0.6621857955147412, 0.6621857955200424, 0.5563449375019521, 0.4218404641258311], history = [2.1307785360594993e6, 1.480831048222229e6, 2.3256776961592715e6, 2.0108110294654146e6, 1.9355945018142546e6, 1.9053542165114614e6, 1.8748334449193429e6, 1.8378496136283926e6, 1.791854694084322e6, 1.7335421670085972e6])

Verification — compliance drops

The coupled compliance must fall as material redistributes:

println("compliance: ", round(hist[1]; digits = 3), " → ", round(hist[end]; digits = 3))
println("volume fraction ≈ 0.5: ", round(sum(x) / length(x); digits = 3))
compliance: 2.130778536e6 → 1.733542167e6
volume fraction ≈ 0.5: 0.498
vtk(m; field = x, fieldname = "density (poroelastic)", edges = true)

Next steps