RasmahRasmah

Species Transport Topology Optimization

Overview

A mixer or separator transports a dissolved species through a fluid. By minimizing the total concentration J = ∫c dx subject to inlet/outlet concentrations, SpeciesTransportTopOpt shapes the flow so the species is flushed out quickly — the layout that best transports and mixes it.

using Rasmah
using LinearAlgebra

Flush a species through a cell

Set the species concentration to 1 on the inlet and 0 on the outlet, drive the flow, and let the optimizer carve the transport channel.

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)
nn = size(m.nodes, 2)

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

inlet  = [i for i in 1:nn if m.nodes[1, i] < 1e-6]
outlet = [i for i in 1:nn if m.nodes[1, i] > 0.999]
fixed_c = vcat(inlet, outlet)
gc = vcat(ones(length(inlet)), zeros(length(outlet)))

x, hist = solve(SpeciesTransportTopOpt(), m, fluid;
                fixed_u = fixed_u, loads_u = loads_u, fixed_p = [1], gp = [0.0],
                fixed_c = fixed_c, gc = gc, volfrac = 0.5, max_iter = 20)
([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0  …  0.5656033455006573, 0.3450397507373531, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [0.500303449064355, 0.43842765039990955, 0.3344998943192888, 0.20931268117393506, 0.1681716079818875, 0.1681698163390285, 0.16816704985355577, 0.1681652272723082, 0.16816283539872026, 0.16816068971471013, 0.16815976738578844, 0.16815837665544772, 0.1681570212232541, 0.16815638039994152, 0.16815579955764873, 0.16815551297364845, 0.16815568785446752, 0.16815542166155797, 0.1681554765718901, 0.1681554064680417])

Verification — the total concentration drops

println("total concentration: ", 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))
total concentration: 0.5 → 0.168
decreases: true
fluid fraction ≈ 0.5: 0.648
vtk(m; field = x, fieldname = "fluid density", edges = true)

Next steps