RasmahRasmah

MHD Stokes Topology Optimization

Overview

A conducting fluid flowing in a magnetic field experiences Lorentz damping. The inductionless MHD Stokes problem couples the velocity u and pressure p with the damping term −σB₀²u, and MHDStokesTopOpt optimizes the flow layout to minimize the dissipated power. The result is a channel that respects both viscous and magnetic drag.

using Rasmah
using LinearAlgebra

A flow cell in a transverse field

Clamp all six faces and drive the interior with a body force; the optimizer carves the flow path that dissipates the least power under the transverse field B₀.

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)   # no-slip on all faces
    onb || (loads_u[3i - 2] = 1.0)                # body force on the interior
end

x, hist = solve(MHDStokesTopOpt(), m, fluid;
                fixed_u = fixed_u, loads_u = loads_u, fixed_p = [1], gp = [0.0],
                B0 = [0.0, 0.0, 2.0], σ = 1.0, volfrac = 0.5, max_iter = 20)
([0.6995955510274725, 0.7835016540727683, 0.653800788091963, 0.7885265753619852, 0.5859220939082761, 0.8149261678378852, 0.4718568190386477, 0.3505485919007301, 0.0, 0.0  …  0.5598414322444503, 0.47185681903864696, 0.3505485919007298, 0.0, 0.788526575361981, 0.5859220939082739, 0.8149261678378956, 0.6995955510274868, 0.7835016540727563, 0.6538007880919551], [5.443150966325855, 5.252311483596564, 4.989499528958463, 4.731960551940015, 4.495667231827182, 4.283159058673794, 4.149148085074769, 4.0799259078731875, 4.030659040951077, 4.102793320590563, 4.06676310778121, 4.0779158723055255, 4.113692976880646, 4.079446725099967, 4.057841551566778, 4.025657451446996, 4.008216668125859, 3.9937507856434475, 3.9784160617542463, 3.9753543412248264])

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

Next steps