flapjax-model-gen

Match strip-wise aerodynamic corrections (cl0, cla, cm0, cma) — defined about an aircraft's flight shape, with no knowledge of local angle of attack or elastic deformation — to a NASTRAN-ready aerodynamic twist distribution, using flapjax as the differentiable solver in the loop.

This is a thin, fast-moving companion package to flapjax, kept separate so it can be iterated on and released independently.

What's here

Most of the workflow this package supports is already implemented in flapjax itself (flapjax.aero.utils.apply_polar_correction / strip_alpha / project_forcing_to_beam, wired into UVLM via its polar_data / polar_function arguments — see flapjax's models/cantilever_wing/polar_correction.ipynb tutorial). This package adds the pieces that were missing:

• linear_polar.py — a PolarFunction (flapjax's (alpha, data) -> (cl, cd, cm) interface) that evaluates an affine correction, cl = cl0 + cla * alpha, cm = cm0 + cma * alpha, cd = 0, so cl0/cla/cm0/cma data can be dropped straight into UVLM(polar_data=..., polar_function=...) exactly like a tabulated polar.
• rigid.py — rigidly rotates a structure's reference SE(3) frames (rotate_hg), for representing a change in angle of attack as a rotation of the geometry, matching the correction data's own convention, rather than rotating the freestream direction.
• twist_grid.py — make_twisted_grid, a per-spanwise-station generalisation of flapjax's make_rectangular_grid(..., twist=...) (which applies one scalar twist to the whole grid). This is the search space the optimizer below works in.
• strip_forces.py — thin helpers around AeroCase.project_forcing_to_beam for pulling global-frame, per-node force/moment vectors out of a solved case, on either the rigid reference shape or a deformed one.
• twist_match.py — match_twist, a JAX-native Levenberg-Marquardt solver (autodiff Jacobian, jax.lax.while_loop) that finds the twist distribution minimising the residual between a candidate model's forces and a target force field. This is the genuinely new piece: no optimizer for this existed in flapjax.

Suggested workflow

This mirrors the four-step process it was built for; wire the pieces together to fit your model rather than treating this as a locked pipeline — in particular, check the assumption flagged in step 3 before trusting the result.

import jax.numpy as jnp
from flapjax.aero import ConstantFlowField
from flapjax_model_gen import (
    LinearPolar, linear_polar_function, rotate_hg, rigid_corrected_forces,
    make_twisted_grid, match_twist,
)

# 1. Target forces directly from the correction data, on the rigid, undeformed flight shape.
#    `uvlm_corrected` is your UVLM/CoupledAeroelastic built with
#    polar_data=[LinearPolar(cl0, cla, cm0, cma)], polar_function=[linear_polar_function].
target_0 = rigid_corrected_forces(uvlm_corrected, hg=structure.hg0)
target_1 = rigid_corrected_forces(uvlm_corrected, hg=rotate_hg(structure.hg0, jnp.deg2rad(1.0)))

# 2. Static deformation under the corrected aero at AoA=0 (full FSI solve).
sol_0 = coupled_corrected.static_solve(prescribed_dofs=range(6), horseshoe=True)

# 3. Rigidly rotate the *converged, deformed* shape and take a single (no-FSI) aero pass at AoA=1.
#    NOTE: apply_polar_correction evaluates cl0/cla/cm0/cma at each strip's *locally computed* AoA, which in a
#    coupled solve includes any elastic twist -- extrapolating a correction that was only derived from rigid
#    sweeps. That's usually the intended reading of "sectional" data, but confirm it before relying on it.
hg_1 = rotate_hg(sol_0.structure.hg, jnp.deg2rad(1.0))
model_1 = rigid_corrected_forces(uvlm_corrected, hg=hg_1)

# 4. Build an *uncorrected* aero model with per-station twist as the free variable, and match it (on the same,
#    fixed deformed geometry from step 2/3) against whichever of the above you're treating as ground truth.
def forces_fn(twist):
    x0_aero = [make_twisted_grid(m, n, chord, ea, twist=twist)]
    uvlm_plain.set_design_variables(..., x0_aero=x0_aero)
    f0 = rigid_corrected_forces(uvlm_plain, hg=structure.hg0)
    f1 = rigid_corrected_forces(uvlm_plain, hg=rotate_hg(sol_0.structure.hg, jnp.deg2rad(1.0)))
    return jnp.stack([f0, f1])

target = jnp.stack([target_0, target_1])
result = match_twist(forces_fn, target, twist0=jnp.zeros(n + 1))

A single per-strip twist angle is a pure shift of local incidence: it can reproduce a cl0-like offset against whatever lift-curve slope the plain panel aerodynamics already has, but it cannot independently fix a mismatched cla, and it can't inject an independent camber-driven cm0 (a flat panel has no camber). Worth checking that the plain UVLM's native cla per strip is already close to the target before trusting the fit — if it isn't, geometry/discretisation needs adjusting, not just twist.

Building a flapjax model from NASTRAN inputs

flapjax_model_gen.nastran converts a NASTRAN aircraft model into flapjax BeamStructure/UVLM inputs, using pyNastran to parse bulk data. It expects three input files, matching how these models are actually organised:

• Aerodynamic mesh (.bdf, CAERO1 panels) — nastran.parse_caero1_panels extracts each panel's leading/ trailing-edge corner points (already resolved to the global frame, including its own CP coordinate system); nastran.build_local_grid then rediscretizes them into a flapjax local grid.
• Mass model (one or more .nsb files, CONM2 entries only — e.g. one structural, one fuel) — nastran.parse_mass_files sums them into one lumped 6x6 mass matrix per node.
• Stiffness model (.bdf, GRID + CBEAM + PBEAM, referencing MAT1) — nastran.parse_stiffness_bdf builds node coordinates, connectivity, per-element orientation (y_vector) and per-element 6x6 stiffness.

import jax.numpy as jnp
from flapjax.aero import GridDiscretisation, UVLM
from flapjax.structure import BeamStructure
from flapjax.coupled import CoupledAeroelastic
from flapjax_model_gen.nastran import (
    parse_stiffness_bdf, parse_mass_files, read_aero_bdf, parse_caero1_panels, build_local_grid,
)

beam_model = parse_stiffness_bdf("stiffness.bdf")
m_lumped, m_lumped_index = parse_mass_files(["structure.nsb", "fuel.nsb"], node_ids=beam_model.node_ids)

structure = BeamStructure(
    num_nodes=len(beam_model.node_ids),
    connectivity=beam_model.connectivity,
    y_vector=beam_model.y_vector,
    k_cs_index=beam_model.k_cs_index,
    m_lumped_index=m_lumped_index,
)

aero_model = read_aero_bdf("aero.bdf")
panels = parse_caero1_panels(aero_model)  # group/order multi-panel surfaces yourself, e.g. sorted by span
local_grid = build_local_grid(panels, beam_model.coords, m=8, span_axis=1)

uvlm = UVLM(
    grid_shapes=[GridDiscretisation(m=8, n=len(beam_model.node_ids) - 1, m_star=20)],
    dof_mapping=jnp.arange(len(beam_model.node_ids)),
)
wing = CoupledAeroelastic(structure, uvlm)
wing.set_design_variables(
    coords=beam_model.coords, k_cs=beam_model.k_cs, m_cs=None, m_lumped=m_lumped,
    dt=..., flowfield=..., x0_aero=[local_grid],
)

What's deliberately out of scope / assumed, each backed by a clear error or warning rather than a silent wrong answer if violated:

• Only CAERO1, CBEAM/PBEAM (not PBEAML/PBCOMP), and MAT1 are supported.
• CBEAM's OFFT must be the default 'GGG', and end offsets (WA/WB) must be zero — no offset beam ends.
• CONM2's CID must be the basic coordinate system (0 or -1).
• A tapered PBEAM (multiple stations) is collapsed to one constant cross-section per element by averaging its station values — flapjax elements are constant-property, so this is an approximation for a genuinely tapered element, not a limitation you can configure around.
• PBEAM's I12 (bend-bend coupling) is not incorporated into k_cs — a warning fires if it's nonzero, since results will be approximate for that element.
• The aerodynamic grid's spanwise stations are placed exactly at the given structural node positions (flapjax's dof_mapping ties each aero column rigidly to one structural node, with no separate spline layer like NASTRAN's beam splines) — build_local_grid doesn't take an independent spanwise panel count.
• The axis convention (CBEAM's orientation vector = flapjax's y_vector directly; PBEAM I1→bending about local z, I2→about local y) was verified against the NASTRAN QRG and cross-checked against flapjax's own beam code, not assumed — see the docstring in nastran/stiffness.py for the derivation.

Development

uv sync --dev
uv run pytest
uv run ruff check src

Tests are deliberately lightweight (geometry/algebra invariants and a synthetic least-squares problem for the optimizer) rather than full aeroelastic solves, so the suite stays fast to iterate against.