RasmahRasmah

Wave Control Topology Optimization

Overview

Elastodynamic cloaking asks: where should material go so that a transient wave barely reaches a protected region? WaveControlTopOpt minimizes the time-integrated wave energy (displacement + velocity) in a target region of an undamped elastic body under a transient load, marched with Newmark integration.

using Rasmah

Protect a region from a pulse

Clamp one face, fire a Gaussian pulse in time at a source, and ask the optimizer to keep the wave energy out of the far side.

m = tetrahedralize_box(1.0, 1.0, 1.0, 3, 3, 3)
nn = size(m.nodes, 2)
mat = Material(youngs_modulus = 1e9, poisson_ratio = 0.3, density = 1000.0)

fixed = Int[]
for i in 1:nn
    m.nodes[1, i] < 1e-6 && append!(fixed, 3i - 2, 3i - 1, 3i)
end
src = argmin(abs.(m.nodes[1, :] .- 0.2))

load_hist = [zeros(3nn) for _ in 1:8]
for k in 1:8
    load_hist[k][3src - 2] = exp(-((k - 3)^2) / 2)     # Gaussian pulse in time
end
target = [3i - 2 for i in 1:nn if m.nodes[1, i] > 0.9]

x, hist = solve(WaveControlTopOpt(), m, mat;
                fixed = fixed, load_history = load_hist, target = target,
                Δt = 0.01, volfrac = 0.5, max_iter = 20)
([0.9710436615552276, 0.5688473870859145, 0.49827188350137935, 0.49827188350137935, 0.5688496618317205, 0.9710903714581862, 0.2042007293342635, 0.47505250914789, 0.48925470981650937, 0.4892461561012995  …  0.5113973860916968, 0.5113957729977687, 0.494496810833131, 0.4876899028315846, 0.48892807137972083, 0.48948005785321747, 0.4965300719296656, 0.49652762202338807, 0.4894806219607384, 0.48892076497831055], [2.2933480200016517e-13, 1.9445312924379173e-13, 1.6803815094574964e-13, 1.4382086514350173e-13, 1.2208569153089767e-13, 1.0295410270463964e-13, 8.638616345670027e-14, 7.222487039143683e-14, 6.024777479015799e-14, 5.020084575859794e-14, 4.1827721172640755e-14, 3.488681841243597e-14, 2.9159860840313606e-14, 2.4453834464452426e-14, 2.0599496378356397e-14, 1.744868317978189e-14, 1.4872079283052982e-14, 1.2758030134970095e-14, 1.1012172311445222e-14, 9.556965177056629e-15])

Verification — the protected energy shrinks

The time-integrated energy in the target region must fall:

println("target energy: ", round(hist[1]; sigdigits = 3), " → ", round(hist[end]; sigdigits = 3))
println("energy decreases: ", hist[end] < hist[1])
println("density in [0,1]: ", all(0 .<= x .<= 1))
target energy: 2.29e-13 → 9.56e-15
energy decreases: true
density in [0,1]: true
vtk(m; field = x, fieldname = "density", edges = true)

Next steps