Kármán vortex street

State the fluid case

Set the expansion order: 3 for a quick run, 9 to reproduce the reference. The solve is graded, so higher order expansions preserve lower order terms. fluid_model reads the mesh, assembles the FEM, and solves for the steady flow at $\mathrm{Re}_0 = 49.03$.

order = 3  # Change to 9 for the reference calculation.
case = NSE.fluid_model(joinpath(@__DIR__, "cylinder_flow.msh"); Re = 49.03)

Select the Hopf pair

The master modes span the manifold's tangent space: here the Hopf pair, the least-damped conjugate eigenpair of $\mathbf (B_0 + \lambda\, \mathbf B_1) y = \mathbf{0}$. Choosing them is a modelling decision, made here rather than by the backend. case.B is the pair $(\mathbf B_0, \mathbf B_1)$.

spectrum = NSE.solve_hopf_eigenproblem(case.B;
    nev = 40, sigma_re = 3.0, sigma_im = 8.0)

master = [spectrum.hopf_index, spectrum.conjugate_index]
outer = setdiff(eachindex(spectrum.eigenvalues), master)
spectrum.eigenvalues[master] # print the Hopf pair

At $\mathrm{Re}_0 = 49.03$ the pair sits just right of the imaginary axis: marginally unstable, shedding at $\omega_0 = 16.86$ rad/s (St ≈ 0.27). With $r$ the reduced coordinates are $(z_1, \bar z_1, r)$.

Build the MORFE model

build_model turns the case and its spectrum into the physics-independent model and spectral data, to be used by MORFE.jl famous parametrise function.

(; model, spectral, meta) = build_model(case, spectrum;
    master = master, outer = outer, expansion_order = order, scale = 1e-2)
meta.conjugate_permutation # print the conjugate permutation

scale is a mode gauge applied to both eigenvector sides. It changes the scales of $W$ and $R$. meta.conjugate_permutation comes back [2, 1, 3].

Parametrise the invariant manifold

Normal-form style keeps only the resonant monomials, and for a Hopf pair that leaves a single complex equation: the Stuart–Landau equation governing the shedding amplitude.

W, R = parametrise(model, spectral, order;
    resonance = ResonanceConfig(style = :complex_normal_form, tol_relative = 0.1,
        outer_targets = true))
R # print reduced dynamics

$\mathrm{Re}(\lambda_1)$ is the linear growth rate, and $\mathrm{Im}(\lambda_1)$ is the linear frequency of the vortex shedding. At cubic order, $\mathrm{Re}(c_{210}) < 0$ makes the Hopf bifurcation supercritical, which explains why a stable vortex street appears smoothly from $\mathrm{Re}_c$ towards larger Reynolds numbers.

Read the bifurcation diagram off $R$

The first reduced cordinate has dynamics $\dot{z}_1 = R_1(z_1, \overline{z}_1, r)$. We write $z_1$ in polar coordinates $z_1 = \rho e^{i\Omega t}$ such that $\rho$ is the amplitude of the shedding and $\Omega$ its frequency. Complex normal form then yields $\dot\rho = \mathrm{Re}\,R_1(\rho, \rho, r)$ and $\Omega = \mathrm{Im}\,R_1(\rho, \rho, r) / \rho$.

Re₀ = 49.03
to_Re(η) = 1 / (η + 1 / Re₀)
to_η(Re) = 1 / Re - 1 / Re₀

sweep = (parameter = 1, sheet = :primary, amplitudes = range(0, 4; length = 500),
    parameter_range = (to_η(70.0), to_η(48.0)))

branch = normal_form_branch(R; sweep...)
Re_c = to_Re(branch.parameter[1]) # the ρ = 0 end of the branch is the Hopf point

sheet = :primary follows the branch out of the Hopf point rather than returning every root, which from order 9 includes a second sheet at large amplitude.

A physical observable: the lift

The lift is the pressure traction on the cylinder. Pushed through $W$ it becomes a polynomial in $(z_1, \bar z_1, r)$, evaluated around the circular orbit.

l_free, L0 = NSE.lift_functional(case)
L_coeffs, mset_L = NSE.lift_polynomial(W, l_free)
L = DensePolynomial(L_coeffs, mset_L)

θ = range(0, 2π; length = 257)[1:256]

function max_lift(P, ρ, η)
    base = evaluate(P, [0.0im, 0.0im, complex(η)])
    maximum(abs(real(evaluate(P, [ρ * cis(t), ρ * cis(-t), complex(η)]) - base)) for t in θ)
end

max_lift(L, branch.amplitude[end], branch.parameter[end]) # peak lift at the far end

Convergence with expansion order

Truncating $R$ and the projected lift to degree $N$ and reading the branch again gives one curve per order, with no re-solve. Truncating the projected lift rather than $W$ matters: at order 9 that array is ~193 MB, and its projection is a few hundred numbers.

using CairoMakie

D = case.fom.reference_length      # 0.1 m, the reference length in ν = D/Re
strouhal(Ω) = Ω * D / (2π * NSE.U_MEAN)

curves = map(3:2:order) do N
    b = normal_form_branch(restrict_ReducedDynamics_to_degree(R, N); sweep...)
    P = MORFE.Polynomials.restrict_polynomial_to_degree(L, N)
    (; N, b.amplitude, b.frequency, eta = b.parameter, Re = to_Re.(b.parameter),
        St = strouhal.(b.frequency), lift = max_lift.(Ref(P), b.amplitude, b.parameter))
end

fig = Figure(size = (900, 360))
ax_lift = Axis(fig[1, 1]; xlabel = "Re", ylabel = "max |lift|")
ax_st = Axis(fig[1, 2]; xlabel = "Re", ylabel = "Strouhal number")
for c in curves
    lines!(ax_lift, c.Re, c.lift; label = "order $(c.N)")
    lines!(ax_st, c.Re, c.St; label = "order $(c.N)")
end
axislegend(ax_lift; position = :lt)
fig

Save the ROM and the branch

The common saver writes $W$, $R$, their coefficient table and a run summary to the example's results directory. The table is what the example's validate.jl checks against the reference, and the branch written beside it is what the figures above are drawn from.

MORFE.save_rom(joinpath(@__DIR__, "results"), W, R;
    external_system = meta.external_system)

open(joinpath(@__DIR__, "results", "data", "branch.csv"), "w") do io
    println(io, "order,eta,Re,rho,omega,St,max_abs_lift")
    for c in curves, i in eachindex(c.Re)
        println(io, join((c.N, c.eta[i], c.Re[i], c.amplitude[i], c.frequency[i],
            c.St[i], c.lift[i]), ","))
    end
end