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. Pass include=[...] or exclude=[...] (CAERO element ids, mutually exclusive) to convert only a subset of surfaces from a multi-surface aero model, or to drop an unsupported (non-CAERO1) entry you don't need.
• 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],
)

Multiple aerodynamic surfaces

A model with more than one aerodynamic surface (wing, tail, fin, ...) needs grid_shapes/x0_aero/ dof_mapping as parallel lists, one entry per surface, each pointing at its own subset of the single, shared BeamStructure's nodes (flapjax's UVLM.dof_mapping is (n_surf, )(zeta_n, )). Since a NASTRAN aero model's own SPLINE/SET1 cards already record exactly this "which CAERO panels + which structural GRID ids form one continuous surface" grouping, nastran.parse_spline_surfaces reads it directly instead of requiring it to be hand-listed, and nastran.reorder_by_surfaces/nastran.reorder_nodes reindex a NastranBeamModel's nodes (and elements) to match:

from flapjax_model_gen.nastran import parse_spline_surfaces, reorder_by_surfaces

surfaces = parse_spline_surfaces(aero_model)  # groups CAERO ids sharing a SET1 into one SplineSurface each
# only needed if SET1 lists auxiliary LE/TE spline-node ids rather than real structural GRID ids directly --
# see "Auxiliary spline nodes" below.
# surfaces = resolve_spline_surfaces(surfaces, read_spline_nodes_bdf("spline_nodes.bdf"))
beam_model = reorder_by_surfaces(
    beam_model, [s.node_ids for s in surfaces], span_axis=1
)  # nodes grouped surface-by-surface, root-to-tip within each; unspliced nodes appended at the end

grid_shapes, x0_aero, dof_mapping = [], [], []
index_of = {nid: i for i, nid in enumerate(beam_model.node_ids)}
for surface in surfaces:
    panels = sorted(
        parse_caero1_panels(aero_model, include=surface.caero_ids), key=lambda p: p.le_root[1]
    )
    node_ids = sorted(surface.node_ids, key=lambda nid: beam_model.coords[index_of[nid], 1])
    mapping = jnp.array([index_of[nid] for nid in node_ids])
    grid_shapes.append(GridDiscretisation(m=8, n=len(node_ids) - 1, m_star=20))
    x0_aero.append(build_local_grid(panels, beam_model.coords[mapping], m=8, span_axis=1))
    dof_mapping.append(mapping)

uvlm = UVLM(grid_shapes=grid_shapes, dof_mapping=dof_mapping)

reorder_by_surfaces only groups/orders the structural side; it's a plain nastran.reorder_nodes(beam_model, node_order) under the hood, callable directly if you want to reorder by something other than a spline grouping. Only SPLINE1 is supported (raises NotImplementedError for SPLINE2-5, which pyNastran does parse into an object but this package doesn't handle). SPLINE6/SPLINE7 aren't implemented by pyNastran itself at all -- they land in model.reject_cards, not model.splines -- and a SPLINE7's surface can't be recovered reliably from the raw card either (a blank-CAERO SPLINE7 splines an AELIST of aero boxes that may span multiple CAERO1 panels rather than naming one directly), so any surface using one of these is skipped entirely, with a UserWarning, rather than guessed at or left to raise.

Auxiliary spline nodes

Some models don't splice a SPLINE1 directly to the beam's own structural nodes -- instead the SET1 lists a separate set of auxiliary leading/trailing-edge GRID points (defined in their own file, e.g. spline_nodes.bdf), each rigidly tied to one real structural node via an RBE2 (any PLOTEL entries, used only for plotting, are ignored). nastran.read_spline_nodes_bdf/nastran.resolve_spline_surfaces translate a parse_spline_surfaces result's node ids through those RBE2s into real structural node ids before calling reorder_by_surfaces -- an LE/TE pair rigidly tied to the same station collapses to that one node id, and an id not covered by any RBE2 (i.e. the SET1 already lists a real structural node) passes through unchanged, so it's always safe to call.

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.