From a mesh to a ROM

State the case

Set the expansion order. The cohomological solve is graded, so a degree-N coefficient never depends on a higher degree: orders 3, 5 and 7 are exact truncations of this one order-9 solve rather than three more runs. mechanical_model reads the mesh, the elastic material properties, damping, and boundary information. Pass the finite-element and quadrature orders to mechanical_model to override the defaults (fe_order = 2, quad_order = fe_order + 1).

dirichlet = "Dirichlet" names a facet group in the Gmsh file, for the two clamped beam ends. If your mesh uses a different physical-group name, pass that name instead.

order = 9
case = SVK.mechanical_model(joinpath(@__DIR__, "clamped_clamped_beam.msh");
    material = SVK.SVKMaterial(E = 160e3, ν = 0.22, ρ = 2.32e-3),
    damping = SVK.RayleighDamping(α = 0.0, β = 0.0),
    dirichlet = "Dirichlet")

Build a physics-independent model for MORFE

build_model converts the structural case into MORFE's physics-independent model and computes its spectral data. master = [1] selects the first conjugate vibration-mode pair. meta holds auxiliary backend information.

(; model, spectral, meta) = build_model(case; master = [1], expansion_order = order)
meta.spectrum.eigenvalues[meta.master_indices] # print master eigenvalues

Compute the invariant manifold

parametrise computes the manifold parametrisation W and associated reduced dynamics R up until the chosen order. The ResonanceConfig keyword selects the parametrisation style and the resonance tolerance.

W, R = parametrise(model, spectral, order;
    resonance = ResonanceConfig(style = :complex_normal_form, tol = 0.05))
R # print reduced dynamics

Obtain the backbone

Because the system is conservative, the master eigenvalues form a purely imaginary complex conjugate pair: $\lambda_1 = \overline{\lambda_2} = -\lambda_2$. Master coordinate $z_1$ has dynamics $\dot{z}_1 = R_1(z_1, z_2)$. Complex normal form produces dynamics in polar coordinates, $z_1 = \rho e^{i \theta}$:

The beam is conservative, so $\dot\rho = 0$: every amplitude is a periodic orbit, and $\Omega = \dot\theta$ is an explicit polynomial in $\rho$. The function normal_form_branch traces the backbone $\Omega(\rho)$.

$\rho$ is a reduced coordinate, not a physical quantity. observable_polynomial projects W onto one degree of freedom to recover a physical displacement, and cycle_amplitude reads its amplitude over a cycle. FIG · 01 highlights the y displacement of node 289.

u = observable_polynomial(W, SVK.probe_dof(case, 289, 2)) # y displacement at mid-span

thickness = 10.0                  # beam thickness, mesh length units
ρ = range(0, 85; length = 341)    # modal amplitude; ρ = 85 is ≈ 1.1 × thickness at the probe
ω₀ = abs(imag(meta.spectrum.eigenvalues[meta.master_indices[1]])) # linear eigenfrequency

curves = map(3:2:order) do N
    b = normal_form_branch(restrict_ReducedDynamics_to_degree(R, N); amplitudes = ρ)
    P = MORFE.Polynomials.restrict_polynomial_to_degree(u, N)
    (; N, b.amplitude, b.frequency, dω = 100 .* (b.frequency ./ ω₀ .- 1),
        a = cycle_amplitude.(Ref(P), b.amplitude) ./ thickness)
end

Plot it

The backbone stiffens: at a displacement of about half of its thickness, the frequency has risen roughly 6%. Plotting four truncations together shows where the expansion stops converging, which is the honest way to read the useful range of a ROM.

using CairoMakie
fig = Figure(size = (520, 380))
ax = Axis(fig[1, 1]; xlabel = "Δω / ω₀  [%]", ylabel = "displacement / thickness  [%]")
foreach(c -> lines!(ax, c.dω, 100 .* c.a; label = "order $(c.N)"), curves)
axislegend(ax; position = :lt)
fig

Save the ROM and the backbone

MORFE's common saver, save_rom, writes W, R and a summary to the results directory.

MORFE.save_rom(joinpath(@__DIR__, "results"), W, R; drop_below = 0.0)