flapjax-model-gen

Generate a flapjax aeroelastic model — structure, mass, aerodynamic mesh, control surfaces, half-models — directly from NASTRAN input files, using pyNastran to parse bulk data. Also ships a companion workflow for matching strip-wise aerodynamic corrections to a NASTRAN-ready twist distribution.

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

Building a flapjax model from NASTRAN inputs

flapjax_model_gen.nastran converts a NASTRAN aircraft model into a fully populated flapjax CoupledAeroelastic, ready for .reference_configuration()/.static_solve(). It expects three input files, matching how these models are actually organised:

• Stiffness model (.bdf, GRID + CBEAM + PBEAM, referencing MAT1) — parse_stiffness_bdf builds node coordinates, connectivity, per-element orientation (y_vector) and per-element 6x6 stiffness.
• Mass model (one or more .nsb files, CONM2 entries only — e.g. one structural, one fuel) — parse_mass_files sums them into one lumped 6x6 mass matrix per node.
• Aerodynamic mesh (.bdf, CAERO1 panels + SPLINE1/SET1) — read_aero_bdf reads it, and parse_spline_surfaces groups CAERO1 panels into continuous surfaces using the model's own SPLINE1/SET1 entries, instead of requiring you to hand-list which CAERO eids and structural GRID ids belong together.

Everything derivable from these files is handled automatically: node topology, stiffness, lumped mass, which CAERO1 panels form one continuous surface (including auto-combining a chordwise split, e.g. a main surface + control surface, and a spanwise split, e.g. inboard/outboard discretisation), each surface's aerodynamic grid rediscretized to sit exactly at the structural node positions, the aero-to-structure dof_mapping, and (for a named control surface) the grid_func that deflects it. Everything else — flight condition, wake truncation, boundary constraints, gravity, mirror symmetry — has no NASTRAN-derivable value and is a plain keyword argument to build_coupled_aeroelastic, mirroring the underlying flapjax constructor it's forwarded to; see that function's own docstring for the full parameter list.

import jax.numpy as jnp
from flapjax.aero import ConstantFlowField
from flapjax_model_gen.nastran import (
    parse_stiffness_bdf, parse_mass_files, read_aero_bdf, parse_spline_surfaces, reorder_by_surfaces,
    SurfaceGridSpec, build_coupled_aeroelastic,
)

beam_model = parse_stiffness_bdf("stiffness.bdf")
aero_model = read_aero_bdf("aero.bdf")
surfaces = parse_spline_surfaces(aero_model)  # one SplineSurface per SPLINE1/SET1 grouping

beam_model = reorder_by_surfaces(beam_model, [s.node_ids for s in surfaces], span_axis=1)
m_lumped, _ = parse_mass_files(
    ["structure.nsb", "fuel.nsb"], node_ids=beam_model.node_ids, node_coords=beam_model.coords
)

wing = SurfaceGridSpec(name="wing", surface=surfaces[0], span_axis=1, m=8, m_star=20)
case = build_coupled_aeroelastic(
    beam_model,
    [wing],
    aero_model,
    m_lumped=m_lumped,
    dt=0.01,
    flowfield=ConstantFlowField(u_inf=jnp.array([30.0, 0.0, 0.0]), rho=1.225, relative_motion=True),
)

Multiple aerodynamic surfaces

A model with more than one surface (wing, tail, fin, ...) just needs more SurfaceGridSpecs — no manual grid_shapes/x0_aero/dof_mapping bookkeeping. Each spec names its surface (an arbitrary label, just used to cross-reference control surfaces below) and its own span_axis — e.g. a wing/horizontal tail spans in y, a vertical tail in z; a single shared axis can't sort both correctly:

surfaces = parse_spline_surfaces(aero_model)
# parse_spline_surfaces's own order is ascending SET1-id, not semantically meaningful -- identify each
# surface by its caero_ids/node_ids (from your model), not by position.
wing_surface, htail_surface, vtail_surface = surfaces

beam_model = reorder_by_surfaces(
    beam_model,
    [wing_surface.node_ids, htail_surface.node_ids, vtail_surface.node_ids],
    span_axis=[1, 1, 2],
)

specs = [
    SurfaceGridSpec(name="wing", surface=wing_surface, span_axis=1, m=8, m_star=20),
    SurfaceGridSpec(name="htail", surface=htail_surface, span_axis=1, m=6, m_star=20),
    SurfaceGridSpec(name="vtail", surface=vtail_surface, span_axis=2, m=6, m_star=20),
]
case = build_coupled_aeroelastic(beam_model, specs, aero_model, m_lumped=m_lumped, dt=..., flowfield=...)

reorder_by_surfaces only groups/orders the structural side; it's a plain 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, so any surface using one of these is skipped entirely, with a UserWarning, rather than guessed at or left to raise.

Control surfaces

Name a control surface with its spanwise extent and chordwise hinge fraction; build_coupled_aeroelastic locates it on the already-built grid and wires it into a grid_func that exposes the name as a deflection-angle keyword:

from flapjax_model_gen.nastran import ControlSurfaceRequest, nudged_x_over_c

# a hinge fraction rarely lands on a nice multiple of 1/m -- nudged_x_over_c keeps the rest of the chordwise
# discretisation uniform and only adjusts the one station nearest the hinge onto it exactly.
wing = SurfaceGridSpec(name="wing", surface=wing_surface, span_axis=1, x_over_c=nudged_x_over_c(20, 0.73))

case = build_coupled_aeroelastic(
    beam_model,
    [wing],
    aero_model,
    control_surfaces={
        "aileron": ControlSurfaceRequest(surface_name="wing", span_lo=8.0, span_hi=10.0, hinge_frac=0.73),
    },
    m_lumped=m_lumped,
    dt=...,
    flowfield=...,
)

deflected = case.aero.grid_func(case.aero.zeta_b0, aileron=jnp.deg2rad(5.0))

Half (symmetric) models

For a half-model mirrored about a plane (e.g. the aircraft's x-z plane, normal to the global y axis), reduce the full-model pieces first: half_model drops the negative-side structure (keeping the mirror-plane nodes) and halves the stiffness of any element lying entirely in the plane; half_model_lumped_mass does the nodal equivalent for lumped mass; restrict_surface_nodes pares each surface's node list to match. Each surface's own CAERO1 panels are then filtered automatically (dropping the negative/mirror-plane side) via SurfaceGridSpec.mirror_axis:

from flapjax_model_gen.nastran import half_model, half_model_lumped_mass, reindex_lumped_mass, restrict_surface_nodes

full_node_ids = beam_model.node_ids
beam_model = half_model(beam_model, mirror_axis=1)
m_lumped, _ = reindex_lumped_mass(m_lumped, full_node_ids, beam_model.node_ids)
m_lumped = half_model_lumped_mass(m_lumped, beam_model.coords, mirror_axis=1)
surfaces = [restrict_surface_nodes(s, beam_model.node_ids) for s in surfaces]

wing = SurfaceGridSpec(name="wing", surface=surfaces[0], span_axis=1, m=8, mirror_axis=1)
case = build_coupled_aeroelastic(
    beam_model,
    [wing],
    aero_model,
    m_lumped=m_lumped,
    dt=...,
    flowfield=...,
    mirror_point=jnp.zeros(3),
    mirror_normal=jnp.array([0.0, 1.0, 0.0]),  # tells flapjax's UVLM about the symmetry plane
)

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). read_spline_nodes_bdf/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:

from flapjax_model_gen.nastran import read_spline_nodes_bdf, resolve_spline_surfaces

surfaces = resolve_spline_surfaces(surfaces, read_spline_nodes_bdf("spline_nodes.bdf"))

Checking panel/surface groupings visually

plot_caero1_panels plots each CAERO1's quadrilateral outline in 3D (matplotlib), labelled by its element id, and plot_beam_model does the equivalent for the structural side -- both useful for sanity-checking span_axis/mirror_axis choices and parse_spline_surfaces groupings before wiring up the full model:

from flapjax_model_gen.nastran import parse_caero1_panels, plot_beam_model, plot_caero1_panels

ax = plot_caero1_panels(parse_caero1_panels(aero_model), surfaces=surfaces)  # colours each panel by its surface
ax.figure.show()  # or ax.figure.savefig("panels.png")

ax = plot_beam_model(beam_model)
ax.figure.show()

Lower-level building blocks

build_coupled_aeroelastic/SurfaceGridSpec are a convenience layer over functions you can call directly for finer control: build_surface_grid (one surface's aero grid, given its SplineSurface and node coordinates) is itself combine_surface_panels (auto-detect/merge a chordwise split into one full-chord panel, order root-to-tip) + build_local_grid (rediscretize). control_surface_indices computes the m_slice/n_slice/ hinge_axis flapjax's own add_control_surface needs, if you're assembling a grid_func by hand. dof_mapping_for_surfaces resolves a surface's node ids to structural indices, if you're assembling a UVLM by hand rather than through build_coupled_aeroelastic.

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), MAT1, and SPLINE1 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); mass is taken entirely from CONM2 mass files, not PBEAM's NSM or any distributed cross-section mass.
• 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) — only the chordwise discretization (m/x_over_c) is a free choice.
• A control surface's hinge must land exactly on a chordwise station (nudged_x_over_c/hinged_x_over_c build one that does for an arbitrary hinge fraction), and its spanwise extent must land exactly on existing structural nodes — there's no clipping at an arbitrary location on either axis.
• There's no NASTRAN SPC/SPC1 parsing, so boundary constraints (a clamped root, a multibody hinge) are built as flapjax constraint objects by hand and passed to build_coupled_aeroelastic(constraints=...).
• 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.

Matching aerodynamic corrections to a twist distribution

A secondary workflow: given 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), find the NASTRAN- ready twist distribution that best reproduces the corrected forces without needing the correction itself at solve time. Most of the underlying machinery (apply_polar_correction/strip_alpha/project_forcing_to_beam, wired into UVLM via polar_data/polar_function) already lives in flapjax itself — see its models/cantilever_wing/polar_correction.ipynb tutorial. This package adds:

• linear_polar.py — a PolarFunction evaluating cl = cl0 + cla * alpha, cm = cm0 + cma * alpha, cd = 0, so cl0/cla/cm0/cma data drops straight into UVLM(polar_data=..., polar_function=...).
• rigid.py — rotate_hg, rigidly rotating a structure's reference SE(3) frames, for representing a change in angle of attack as a rotation of the geometry rather than the freestream.
• twist_grid.py — make_twisted_grid, a per-spanwise-station generalisation of flapjax's make_rectangular_grid(..., twist=...) — 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.
• twist_match.py — match_twist, a JAX-native Levenberg-Marquardt solver (autodiff Jacobian, jax.lax.while_loop) finding the twist distribution minimising the residual between a candidate model's forces and a target force field.

import jax.numpy as jnp
from flapjax_model_gen import LinearPolar, linear_polar_function, make_twisted_grid, match_twist, rigid_corrected_forces, rotate_hg

# Target forces from the correction data on the rigid flight shape (uvlm_corrected built with
# polar_data=[LinearPolar(cl0, cla, cm0, cma)], polar_function=[linear_polar_function]), at two angles of attack.
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)))

# Match an *uncorrected* model's per-station twist (the free variable) against those targets.
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(structure.hg0, jnp.deg2rad(1.0)))
    return jnp.stack([f0, f1])

result = match_twist(forces_fn, jnp.stack([target_0, target_1]), 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.

Development

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

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