Swept Pazy wing flutter analysis¶

Evaluate the aeroelastic eigenvalues (and resulting frequency and damping) for key modes for three variants of the Pazy wing to compare the effect of sweep:

1. Swept Pazy wing with the tip mass placed at the leading edge.
2. Swept Pazy wing with the tip mass placed at the trailing edge.
3. Straight Pazy wing rotated by 10 degrees (kinematic-only sweep — retains the straight-wing structural properties).

Imports

In [ ]:

from copy import deepcopy
from functools import partial
from typing import Literal, cast

import jax
from jax import Array
from jax import numpy as jnp
from matplotlib import pyplot as plt
from matplotlib.colors import Normalize

from flapjax.aero.flowfields import ConstantFlowField
from flapjax.coupled import CoupledAeroelastic
from flapjax.models.pazy.straight.pazy_wing import generate_pazy_wing
from flapjax.models.pazy.swept.swept_pazy_wing import generate_swept_pazy_wing
from flapjax.utils.print_utils import set_verbosity

from copy import deepcopy from functools import partial from typing import Literal, cast import jax from jax import Array from jax import numpy as jnp from matplotlib import pyplot as plt from matplotlib.colors import Normalize from flapjax.aero.flowfields import ConstantFlowField from flapjax.coupled import CoupledAeroelastic from flapjax.models.pazy.straight.pazy_wing import generate_pazy_wing from flapjax.models.pazy.swept.swept_pazy_wing import generate_swept_pazy_wing from flapjax.utils.print_utils import set_verbosity

Discretisation and Parameter Setup

In [ ]:

set_verbosity("silent")  # remove flapjax console prints due to the large number of cases

n_vel = 24  # number of velocity points
m = 12  # wing chordwise discretisation
m_star = 50  # number of wake panels per strip
alpha = jnp.deg2rad(5.0)  # angle of attack shared across all sweep cases

cases: tuple[Literal["swept_LE", "swept_TE", "straight_rotated"], ...] = (
    "swept_LE",
    "swept_TE",
    "straight_rotated",
)
u_inf_mags = jnp.linspace(30.0, 80.0, n_vel)  # velocity sweep

# plot colourmap
cmap = plt.get_cmap("viridis")

# noinspection PyArgumentList
norm = cast(Normalize, Normalize(vmin=float(u_inf_mags.min()), vmax=float(u_inf_mags.max())))

set_verbosity("silent") # remove flapjax console prints due to the large number of cases n_vel = 24 # number of velocity points m = 12 # wing chordwise discretisation m_star = 50 # number of wake panels per strip alpha = jnp.deg2rad(5.0) # angle of attack shared across all sweep cases cases: tuple[Literal["swept_LE", "swept_TE", "straight_rotated"], ...] = ( "swept_LE", "swept_TE", "straight_rotated", ) u_inf_mags = jnp.linspace(30.0, 80.0, n_vel) # velocity sweep # plot colourmap cmap = plt.get_cmap("viridis") # noinspection PyArgumentList norm = cast(Normalize, Normalize(vmin=float(u_inf_mags.min()), vmax=float(u_inf_mags.max())))

Factory for the base wing for a given sweep case. Shared across all velocities within a case.

In [ ]:

def make_wing(case: Literal["swept_LE", "swept_TE", "straight_rotated"]) -> CoupledAeroelastic:
    flowfield = ConstantFlowField(
        u_inf=jnp.array([0.0, 0.0, 0.0]),  # overwritten per velocity in eigs_at_velocity
        rho=1.225,
        relative_motion=True,
    )
    match case:
        case "swept_LE":
            return generate_swept_pazy_wing(
                flowfield=flowfield,
                aoa=alpha,
                tip_mass="LE_CORRECTED",
                m=m,
                m_star=m_star,
                node_multiplier=2,
                sweep_angle=10,
            )
        case "swept_TE":
            return generate_swept_pazy_wing(
                flowfield=flowfield,
                aoa=alpha,
                tip_mass="TE_CORRECTED",
                m=m,
                m_star=m_star,
                node_multiplier=2,
                sweep_angle=10,
            )
        case "straight_rotated":
            return generate_pazy_wing(
                flowfield=flowfield,
                aoa=alpha,
                m=m,
                m_star=m_star,
                node_multiplier=2,
                skin=True,
                variable_disc_wake=True,
                sweep=jnp.deg2rad(10.0),
            )
        case _:
            raise ValueError(f"Unknown sweep case: {case}")

def make_wing(case: Literal["swept_LE", "swept_TE", "straight_rotated"]) -> CoupledAeroelastic: flowfield = ConstantFlowField( u_inf=jnp.array([0.0, 0.0, 0.0]), # overwritten per velocity in eigs_at_velocity rho=1.225, relative_motion=True, ) match case: case "swept_LE": return generate_swept_pazy_wing( flowfield=flowfield, aoa=alpha, tip_mass="LE_CORRECTED", m=m, m_star=m_star, node_multiplier=2, sweep_angle=10, ) case "swept_TE": return generate_swept_pazy_wing( flowfield=flowfield, aoa=alpha, tip_mass="TE_CORRECTED", m=m, m_star=m_star, node_multiplier=2, sweep_angle=10, ) case "straight_rotated": return generate_pazy_wing( flowfield=flowfield, aoa=alpha, m=m, m_star=m_star, node_multiplier=2, skin=True, variable_disc_wake=True, sweep=jnp.deg2rad(10.0), ) case _: raise ValueError(f"Unknown sweep case: {case}")

Function to evaluate the eigenvalues for a given wing and velocity. This is mapped across the input velocities, allowing multiple design points to be evaluated in parallel.

In [ ]:

def eigs_at_velocity(u_mag: Array, wing: CoupledAeroelastic) -> Array:
    # this function is mapped across the input velocities, allowing multiple design points to be evaluated in
    # parallel
    inner_wing = deepcopy(wing)

    # set the required design variables for this case
    dt = 0.1 / (m * u_mag)
    inner_wing.set_design_variables(
        coords=inner_wing.structure.x0,
        k_cs=inner_wing.structure.k_cs,
        m_cs=inner_wing.structure.m_cs,
        m_lumped=None,
        dt=dt,
        flowfield=ConstantFlowField(
            u_inf=jnp.array([u_mag, 0.0, 0.0]),
            rho=1.225,
            relative_motion=True,
        ),
        x0_aero=inner_wing.aero.zeta_b0,
        remove_checks=True,
    )

    # obtain static equilibrium
    static_sol = inner_wing.static_solve(
        prescribed_dofs=tuple(range(6)),
        horseshoe=False,
    )
    # linearise solution about static equilibrium
    linear_sol = inner_wing.linearise(
        reference=static_sol, skip_checks=True, batch_size=4, n_struct_modes=50
    )

    # obtain eigenvalies of linearised system
    eigvals = linear_sol.modal(
        n_modes=5,
        plot_eigvals=False,
        freq_range=(0.0, 50.0),
        damp_range=(-jnp.inf, 0.3),
    )
    return eigvals

def eigs_at_velocity(u_mag: Array, wing: CoupledAeroelastic) -> Array: # this function is mapped across the input velocities, allowing multiple design points to be evaluated in # parallel inner_wing = deepcopy(wing) # set the required design variables for this case dt = 0.1 / (m * u_mag) inner_wing.set_design_variables( coords=inner_wing.structure.x0, k_cs=inner_wing.structure.k_cs, m_cs=inner_wing.structure.m_cs, m_lumped=None, dt=dt, flowfield=ConstantFlowField( u_inf=jnp.array([u_mag, 0.0, 0.0]), rho=1.225, relative_motion=True, ), x0_aero=inner_wing.aero.zeta_b0, remove_checks=True, ) # obtain static equilibrium static_sol = inner_wing.static_solve( prescribed_dofs=tuple(range(6)), horseshoe=False, ) # linearise solution about static equilibrium linear_sol = inner_wing.linearise( reference=static_sol, skip_checks=True, batch_size=4, n_struct_modes=50 ) # obtain eigenvalies of linearised system eigvals = linear_sol.modal( n_modes=5, plot_eigvals=False, freq_range=(0.0, 50.0), damp_range=(-jnp.inf, 0.3), ) return eigvals

Approximate flutter onset and offset points by interpolating for where the damping is zero.

In [ ]:

def flutter_points(min_zeta: Array) -> tuple[float | None, float | None]:
    # interpolate to approximate flutter onset/offset for the hump mode (first two crossings)
    crossover_point = jnp.sign(min_zeta[1:]) != jnp.sign(min_zeta[:-1])
    crossing_indices = jnp.where(crossover_point)[0]

    def interp_zero(idx: int) -> float:
        # linearly interpolate between positive and negative damping points to find velocity where it would be zero
        z0, z1 = float(min_zeta[idx]), float(min_zeta[idx + 1])
        u0, u1 = float(u_inf_mags[idx]), float(u_inf_mags[idx + 1])
        return u0 + (0.0 - z0) / (z1 - z0) * (u1 - u0)

    u_onset = interp_zero(int(crossing_indices[0])) if crossing_indices.size >= 1 else None
    u_offset = interp_zero(int(crossing_indices[1])) if crossing_indices.size >= 2 else None
    return u_onset, u_offset

def flutter_points(min_zeta: Array) -> tuple[float | None, float | None]: # interpolate to approximate flutter onset/offset for the hump mode (first two crossings) crossover_point = jnp.sign(min_zeta[1:]) != jnp.sign(min_zeta[:-1]) crossing_indices = jnp.where(crossover_point)[0] def interp_zero(idx: int) -> float: # linearly interpolate between positive and negative damping points to find velocity where it would be zero z0, z1 = float(min_zeta[idx]), float(min_zeta[idx + 1]) u0, u1 = float(u_inf_mags[idx]), float(u_inf_mags[idx + 1]) return u0 + (0.0 - z0) / (z1 - z0) * (u1 - u0) u_onset = interp_zero(int(crossing_indices[0])) if crossing_indices.size >= 1 else None u_offset = interp_zero(int(crossing_indices[1])) if crossing_indices.size >= 2 else None return u_onset, u_offset

Plot eigenvalues for a given sweep case, with colour representing velocity.

In [ ]:

def plot_eigenvalues(evals: Array, case: str) -> None:
    fig, ax = plt.subplots()
    for i in range(len(u_inf_mags)):
        ax.scatter(
            evals[i].real / (2.0 * jnp.pi),
            evals[i].imag / (2.0 * jnp.pi),
            color=cmap(norm(float(u_inf_mags[i]))),
            s=10,
        )
    ax.set_xlim(-6.0, 2.0)
    ax.set_ylim(0.0, 50.0)
    ax.axvline(0.0, color="k", linewidth=0.5)
    ax.set_xlabel("Re(eig) / (2 * pi) [1/s]")
    ax.set_ylabel("Im(eig) / (2 * pi) [1/s]")
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
    fig.colorbar(sm, ax=ax, label="u_inf [m/s]")
    ax.set_title(f"Pazy wing eigenvalues vs velocity, case={case}")
    plt.show()

def plot_eigenvalues(evals: Array, case: str) -> None: fig, ax = plt.subplots() for i in range(len(u_inf_mags)): ax.scatter( evals[i].real / (2.0 * jnp.pi), evals[i].imag / (2.0 * jnp.pi), color=cmap(norm(float(u_inf_mags[i]))), s=10, ) ax.set_xlim(-6.0, 2.0) ax.set_ylim(0.0, 50.0) ax.axvline(0.0, color="k", linewidth=0.5) ax.set_xlabel("Re(eig) / (2 * pi) [1/s]") ax.set_ylabel("Im(eig) / (2 * pi) [1/s]") sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) fig.colorbar(sm, ax=ax, label="u_inf [m/s]") ax.set_title(f"Pazy wing eigenvalues vs velocity, case={case}") plt.show()

Plot damping ratio and damped frequency for a given sweep case, with colour representing velocity.

In [ ]:

def plot_damp_freq(
        zeta: Array,
        freq_d: Array,
        u_onset: float | None,
        u_offset: float | None,
        case: str,
) -> None:
    fig, (ax_z, ax_f) = plt.subplots(1, 2, figsize=(11, 4.5))
    for i in range(len(u_inf_mags)):
        c = cmap(norm(float(u_inf_mags[i])))
        u = float(u_inf_mags[i])
        ax_z.scatter(jnp.full_like(zeta[i], u), -zeta[i], color=c, s=10)
        ax_f.scatter(jnp.full_like(freq_d[i], u), freq_d[i], color=c, s=10)

    ax_z.axhline(0.0, color="k", linewidth=0.5)
    if u_onset is not None:
        ax_z.axvline(
            u_onset,
            color="r",
            linewidth=0.8,
            linestyle="--",
            label=f"onset = {u_onset:.1f} m/s",
        )
    if u_offset is not None:
        ax_z.axvline(
            u_offset,
            color="b",
            linewidth=0.8,
            linestyle="--",
            label=f"offset = {u_offset:.1f} m/s",
        )
    if u_onset is not None or u_offset is not None:
        ax_z.legend(loc="lower right", fontsize=8)
    ax_z.set_xlabel("u_inf [m/s]")
    ax_z.set_ylabel("-Damping ratio")
    ax_z.set_xlim(u_inf_mags.min(), u_inf_mags.max())
    ax_z.set_ylim(-0.1, 0.1)
    ax_z.set_title(f"Damping ratio, case={case}")

    ax_f.set_xlabel("u_inf [m/s]")
    ax_f.set_ylabel("Damped frequency f_d [Hz]")
    ax_f.set_xlim(u_inf_mags.min(), u_inf_mags.max())
    ax_f.set_ylim(0.0, 50.0)
    ax_f.set_title(f"Damped frequency, case={case}")
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
    fig.colorbar(sm, ax=(ax_z, ax_f), label="u_inf [m/s]")
    plt.show()

def plot_damp_freq( zeta: Array, freq_d: Array, u_onset: float | None, u_offset: float | None, case: str, ) -> None: fig, (ax_z, ax_f) = plt.subplots(1, 2, figsize=(11, 4.5)) for i in range(len(u_inf_mags)): c = cmap(norm(float(u_inf_mags[i]))) u = float(u_inf_mags[i]) ax_z.scatter(jnp.full_like(zeta[i], u), -zeta[i], color=c, s=10) ax_f.scatter(jnp.full_like(freq_d[i], u), freq_d[i], color=c, s=10) ax_z.axhline(0.0, color="k", linewidth=0.5) if u_onset is not None: ax_z.axvline( u_onset, color="r", linewidth=0.8, linestyle="--", label=f"onset = {u_onset:.1f} m/s", ) if u_offset is not None: ax_z.axvline( u_offset, color="b", linewidth=0.8, linestyle="--", label=f"offset = {u_offset:.1f} m/s", ) if u_onset is not None or u_offset is not None: ax_z.legend(loc="lower right", fontsize=8) ax_z.set_xlabel("u_inf [m/s]") ax_z.set_ylabel("-Damping ratio") ax_z.set_xlim(u_inf_mags.min(), u_inf_mags.max()) ax_z.set_ylim(-0.1, 0.1) ax_z.set_title(f"Damping ratio, case={case}") ax_f.set_xlabel("u_inf [m/s]") ax_f.set_ylabel("Damped frequency f_d [Hz]") ax_f.set_xlim(u_inf_mags.min(), u_inf_mags.max()) ax_f.set_ylim(0.0, 50.0) ax_f.set_title(f"Damped frequency, case={case}") sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) fig.colorbar(sm, ax=(ax_z, ax_f), label="u_inf [m/s]") plt.show()

Evaluate each sweep case by mapping the eigenvalue function across the velocity sweep. Onset/offset velocities and per-case plots let us compare the effect of sweep between the three variants.

In [ ]:

flutter_summary: dict[str, tuple[float | None, float | None]] = {}

for case_ in cases:  # sweep configuration comparison
    wing_ = make_wing(case=case_)

    # evaluate eigenvalues
    evals_ = jax.lax.map(partial(eigs_at_velocity, wing=wing_), u_inf_mags, batch_size=2)  # [n_vel, n_modes]

    zeta_ = -evals_.real / jnp.abs(evals_)  # damping ratio [n_vel, n_modes]
    freq_d_ = evals_.imag / (2.0 * jnp.pi)  # damped frequency (Hz), [n_vel, n_modes]

    min_zeta_ = zeta_.min(axis=1)  # minimum damping [n_vel]
    u_onset_, u_offset_ = flutter_points(min_zeta=min_zeta_)  # flutter onset and offset velocity
    flutter_summary[case_] = (u_onset_, u_offset_)

    plot_eigenvalues(evals=evals_, case=case_)
    plot_damp_freq(zeta=zeta_, freq_d=freq_d_, u_onset=u_onset_, u_offset=u_offset_, case=case_)

print(f"Flutter boundaries at aoa={float(jnp.rad2deg(alpha)):.1f} deg:")
for case_, (u_onset_, u_offset_) in flutter_summary.items():
    onset_str = f"{u_onset_:.2f} m/s" if u_onset_ is not None else "n/a"
    offset_str = f"{u_offset_:.2f} m/s" if u_offset_ is not None else "n/a"
    print(f"  {case_:>18}: onset = {onset_str:>10}, offset = {offset_str:>10}")

flutter_summary: dict[str, tuple[float | None, float | None]] = {} for case_ in cases: # sweep configuration comparison wing_ = make_wing(case=case_) # evaluate eigenvalues evals_ = jax.lax.map(partial(eigs_at_velocity, wing=wing_), u_inf_mags, batch_size=2) # [n_vel, n_modes] zeta_ = -evals_.real / jnp.abs(evals_) # damping ratio [n_vel, n_modes] freq_d_ = evals_.imag / (2.0 * jnp.pi) # damped frequency (Hz), [n_vel, n_modes] min_zeta_ = zeta_.min(axis=1) # minimum damping [n_vel] u_onset_, u_offset_ = flutter_points(min_zeta=min_zeta_) # flutter onset and offset velocity flutter_summary[case_] = (u_onset_, u_offset_) plot_eigenvalues(evals=evals_, case=case_) plot_damp_freq(zeta=zeta_, freq_d=freq_d_, u_onset=u_onset_, u_offset=u_offset_, case=case_) print(f"Flutter boundaries at aoa={float(jnp.rad2deg(alpha)):.1f} deg:") for case_, (u_onset_, u_offset_) in flutter_summary.items(): onset_str = f"{u_onset_:.2f} m/s" if u_onset_ is not None else "n/a" offset_str = f"{u_offset_:.2f} m/s" if u_offset_ is not None else "n/a" print(f" {case_:>18}: onset = {onset_str:>10}, offset = {offset_str:>10}")