Cantilever Wing with Polar Corrections¶

The UVLM's own sectional aerodynamics are those of an inviscid flat plate: $c_l = 2\pi\alpha$, $c_d = c_m = 0$ about the quarter chord. However, we can let each surface's sectional forcing be corrected to match tabulated 2D airfoil polar data at the locally-computed strip angle of attack.

This notebook builds a rectangular cantilever wing, attaches asymmetric-airfoil polar , and compares a coupled static aeroelastic angle-of-attack sweep with and without the correction.

Imports¶

In [ ]:

Copied!





from typing import NamedTuple

import jax
from jax import Array, vmap
from jax import numpy as jnp
from matplotlib import pyplot as plt

from flapjax.aero.data_structures import GridDiscretisation
from flapjax.aero.flowfields import ConstantFlowField
from flapjax.aero.utils import make_rectangular_grid
from flapjax.aero.uvlm import UVLM, PolarFunction
from flapjax.coupled import CoupledAeroelastic
from flapjax.models.cantilever_wing.cantilever_wing import K_CS_DEFAULT, M_CS_DEFAULT
from flapjax.structure import BeamStructure
from flapjax.utils.print_utils import set_verbosity

set_verbosity("warning")  # suppress the per-iteration solver tables printed for each sweep point
from typing import NamedTuple

import jax
from jax import Array, vmap
from jax import numpy as jnp
from matplotlib import pyplot as plt

from flapjax.aero.data_structures import GridDiscretisation
from flapjax.aero.flowfields import ConstantFlowField
from flapjax.aero.utils import make_rectangular_grid
from flapjax.aero.uvlm import UVLM, PolarFunction
from flapjax.coupled import CoupledAeroelastic
from flapjax.models.cantilever_wing.cantilever_wing import K_CS_DEFAULT, M_CS_DEFAULT
from flapjax.structure import BeamStructure
from flapjax.utils.print_utils import set_verbosity

set_verbosity("warning")  # suppress the per-iteration solver tables printed for each sweep point

Synthetic airfoil polar¶

A tabulated polar for a stylised symmetric section with mild stall effects. Note that polar_data may be any pytree; here it is a NamedTuple of arrays, and polar_function linearly interpolates it at each strip's angle of attack (which is passed in radians).

In [ ]:

Copied!





class Polar(NamedTuple):
    alpha_deg: Array
    cl: Array
    cd: Array
    cm: Array


# illustrative polar
naca0012_polar = Polar(
    alpha_deg=jnp.array(
        (
            -20.0, -16.0, -13.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0, 0.0,
            2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 13.0, 16.0, 20.0,
        )
    ),
    cl=jnp.array(
        (
            -0.80, -0.75, -1.15, -1.05, -1.00, -0.88, -0.66, -0.44, -0.22, 0.00,
            0.22, 0.44, 0.66, 0.88, 1.00, 1.05, 1.15, 0.75, 0.80,
        )
    ),
    cd=jnp.array(
        (
            0.28, 0.22, 0.16, 0.14, 0.10, 0.06, 0.035, 0.02, 0.012, 0.008,
            0.012, 0.02, 0.035, 0.06, 0.10, 0.14, 0.16, 0.22, 0.28,
        )
    ),
    cm=jnp.array(
        (
            -0.06, -0.07, -0.05, -0.02, -0.01, -0.005, -0.002, -0.001, 0.00, 0.00,
            0.00, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.07, 0.06,
        )
    ),
)


def lookup_polar(alpha: Array, polar: Polar) -> tuple[Array, Array, Array]:
    # linearly interpolate the tabulated cl, cd, cm at the per-strip angle of attack (radians)
    alpha_deg = jnp.rad2deg(alpha)
    cl = jnp.interp(alpha_deg, polar.alpha_deg, polar.cl)
    cd = jnp.interp(alpha_deg, polar.alpha_deg, polar.cd)
    cm = jnp.interp(alpha_deg, polar.alpha_deg, polar.cm)
    return cl, cd, cm
class Polar(NamedTuple):
    alpha_deg: Array
    cl: Array
    cd: Array
    cm: Array


# illustrative polar
naca0012_polar = Polar(
    alpha_deg=jnp.array(
        (
            -20.0, -16.0, -13.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0, 0.0,
            2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 13.0, 16.0, 20.0,
        )
    ),
    cl=jnp.array(
        (
            -0.80, -0.75, -1.15, -1.05, -1.00, -0.88, -0.66, -0.44, -0.22, 0.00,
            0.22, 0.44, 0.66, 0.88, 1.00, 1.05, 1.15, 0.75, 0.80,
        )
    ),
    cd=jnp.array(
        (
            0.28, 0.22, 0.16, 0.14, 0.10, 0.06, 0.035, 0.02, 0.012, 0.008,
            0.012, 0.02, 0.035, 0.06, 0.10, 0.14, 0.16, 0.22, 0.28,
        )
    ),
    cm=jnp.array(
        (
            -0.06, -0.07, -0.05, -0.02, -0.01, -0.005, -0.002, -0.001, 0.00, 0.00,
            0.00, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.07, 0.06,
        )
    ),
)


def lookup_polar(alpha: Array, polar: Polar) -> tuple[Array, Array, Array]:
    # linearly interpolate the tabulated cl, cd, cm at the per-strip angle of attack (radians)
    alpha_deg = jnp.rad2deg(alpha)
    cl = jnp.interp(alpha_deg, polar.alpha_deg, polar.cl)
    cd = jnp.interp(alpha_deg, polar.alpha_deg, polar.cd)
    cm = jnp.interp(alpha_deg, polar.alpha_deg, polar.cm)
    return cl, cd, cm

Plot the tabulated polar against the flat-plate lift the UVLM would otherwise assume.

In [ ]:

Copied!





alpha_plot_deg = jnp.linspace(-20.0, 20.0, 200)
cl_flat_plate = 2.0 * jnp.pi * jnp.deg2rad(alpha_plot_deg)

fig, axes = plt.subplots(1, 2, figsize=(10.0, 4.0))

axes[0].plot(alpha_plot_deg, cl_flat_plate, "--", label=r"UVLM flat plate ($2\pi\alpha$)")
axes[0].plot(naca0012_polar.alpha_deg, naca0012_polar.cl, "o-", label="tabulated polar")
axes[0].set_xlabel("Angle of attack [deg]")
axes[0].set_ylabel("$c_l$")
axes[0].legend()

axes[1].plot(naca0012_polar.alpha_deg, naca0012_polar.cd, "o-", color="tab:orange")
axes[1].set_xlabel("Angle of attack [deg]")
axes[1].set_ylabel("$c_d$")

fig.tight_layout()
plt.show()
alpha_plot_deg = jnp.linspace(-20.0, 20.0, 200)
cl_flat_plate = 2.0 * jnp.pi * jnp.deg2rad(alpha_plot_deg)

fig, axes = plt.subplots(1, 2, figsize=(10.0, 4.0))

axes[0].plot(alpha_plot_deg, cl_flat_plate, "--", label=r"UVLM flat plate ($2\pi\alpha$)")
axes[0].plot(naca0012_polar.alpha_deg, naca0012_polar.cl, "o-", label="tabulated polar")
axes[0].set_xlabel("Angle of attack [deg]")
axes[0].set_ylabel("$c_l$")
axes[0].legend()

axes[1].plot(naca0012_polar.alpha_deg, naca0012_polar.cd, "o-", color="tab:orange")
axes[1].set_xlabel("Angle of attack [deg]")
axes[1].set_ylabel("$c_d$")

fig.tight_layout()
plt.show()

Wing model¶

A rectangular cantilever wing is constructed, which is symmetric about the root. The aerodynamic surface is given the polar database and evaluation function directly through UVLM's polar_data / polar_function inputs; passing None for both disables the correction and recovers the UVLM's own flat-plate forcing.

In [ ]:

Copied!





# problem discretisation
m = 16  # chordwise panels
n = 16  # spanwise panels

b_ref = 5.0  # semi-span
c_ref = 1.0  # aerodynamic chord

# freestream speed, m/s
u_mag = 7.0

n_nodes = n + 1
beam_coords = jnp.zeros((n_nodes, 3)).at[:, 1].set(jnp.linspace(0.0, b_ref, n_nodes))
grid = make_rectangular_grid(m, n, c_ref, ea=0.25)


def build_wing_base(
        polar_data: Polar | None, polar_function: PolarFunction | None
) -> CoupledAeroelastic:
    """Build the wing without design variables"""
    conn = jnp.zeros((n, 2), dtype=int)
    conn = conn.at[:, 0].set(jnp.arange(n))
    conn = conn.at[:, 1].set(jnp.arange(1, n_nodes))
    beam = BeamStructure(
        num_nodes=n_nodes,
        connectivity=conn,
        y_vector=jnp.array((0.0, 0.0, 1.0)),
    )

    gd = GridDiscretisation(m=m, n=n, m_star=1)
    uvlm = UVLM(
        grid_shapes=[gd],
        dof_mapping=jnp.arange(n_nodes),
        mirror_point=jnp.zeros(3),
        mirror_normal=jnp.array((0.0, 1.0, 0.0)),
        polar_data=[polar_data],
        polar_function=[polar_function],
    )

    return CoupledAeroelastic(beam, uvlm)


def set_wing_alpha(
        wing: CoupledAeroelastic, alpha_deg: Array | float, *, remove_checks: bool = False
) -> None:
    """Set the design variables on the wing."""
    u_inf = u_mag * jnp.array(
        (jnp.cos(jnp.deg2rad(alpha_deg)), 0.0, jnp.sin(jnp.deg2rad(alpha_deg)))
    )
    flowfield = ConstantFlowField(u_inf=u_inf, rho=1.225, relative_motion=True)

    wing.set_design_variables(
        coords=beam_coords,
        k_cs=K_CS_DEFAULT,
        m_cs=M_CS_DEFAULT,
        m_lumped=None,
        dt=c_ref / (u_mag * m),
        flowfield=flowfield,
        delta_w=None,
        x0_aero=grid,
        remove_checks=remove_checks,
    )


def build_wing(
        polar_data: Polar | None,
        polar_function: PolarFunction | None,
        alpha_deg: Array | float = 0.0,
) -> CoupledAeroelastic:
    r"""
    We map this function across our inputs to build a batch of wings.
    """
    wing = build_wing_base(polar_data, polar_function)
    set_wing_alpha(wing, alpha_deg)
    return wing
# problem discretisation
m = 16  # chordwise panels
n = 16  # spanwise panels

b_ref = 5.0  # semi-span
c_ref = 1.0  # aerodynamic chord

# freestream speed, m/s
u_mag = 7.0

n_nodes = n + 1
beam_coords = jnp.zeros((n_nodes, 3)).at[:, 1].set(jnp.linspace(0.0, b_ref, n_nodes))
grid = make_rectangular_grid(m, n, c_ref, ea=0.25)


def build_wing_base(
        polar_data: Polar | None, polar_function: PolarFunction | None
) -> CoupledAeroelastic:
    """Build the wing without design variables"""
    conn = jnp.zeros((n, 2), dtype=int)
    conn = conn.at[:, 0].set(jnp.arange(n))
    conn = conn.at[:, 1].set(jnp.arange(1, n_nodes))
    beam = BeamStructure(
        num_nodes=n_nodes,
        connectivity=conn,
        y_vector=jnp.array((0.0, 0.0, 1.0)),
    )

    gd = GridDiscretisation(m=m, n=n, m_star=1)
    uvlm = UVLM(
        grid_shapes=[gd],
        dof_mapping=jnp.arange(n_nodes),
        mirror_point=jnp.zeros(3),
        mirror_normal=jnp.array((0.0, 1.0, 0.0)),
        polar_data=[polar_data],
        polar_function=[polar_function],
    )

    return CoupledAeroelastic(beam, uvlm)


def set_wing_alpha(
        wing: CoupledAeroelastic, alpha_deg: Array | float, *, remove_checks: bool = False
) -> None:
    """Set the design variables on the wing."""
    u_inf = u_mag * jnp.array(
        (jnp.cos(jnp.deg2rad(alpha_deg)), 0.0, jnp.sin(jnp.deg2rad(alpha_deg)))
    )
    flowfield = ConstantFlowField(u_inf=u_inf, rho=1.225, relative_motion=True)

    wing.set_design_variables(
        coords=beam_coords,
        k_cs=K_CS_DEFAULT,
        m_cs=M_CS_DEFAULT,
        m_lumped=None,
        dt=c_ref / (u_mag * m),
        flowfield=flowfield,
        delta_w=None,
        x0_aero=grid,
        remove_checks=remove_checks,
    )


def build_wing(
        polar_data: Polar | None,
        polar_function: PolarFunction | None,
        alpha_deg: Array | float = 0.0,
) -> CoupledAeroelastic:
    r"""
    We map this function across our inputs to build a batch of wings.
    """
    wing = build_wing_base(polar_data, polar_function)
    set_wing_alpha(wing, alpha_deg)
    return wing

Angle-of-attack sweep¶

Solve the coupled static aeroelastic problem at a range of angles of attack, with and without the polar correction, and compare the resulting wing-average lift and drag coefficients. Since the planform is rectangular with uniform-width strips, the mean of the per-strip cl / cd over the span is the wing's overall coefficient.

We here use the vmap operator to vectorise the sweep over angles of attack, and jit to compile the solver for speed. This allows up to evalute all angles of attack in a matter of seconds, as we only need one compilation for each of the two cases (with and without polar correction).

In [ ]:

Copied!





def solve_at_alpha(wing: CoupledAeroelastic, alpha_deg: Array) -> tuple[Array, Array]:
    set_wing_alpha(wing, alpha_deg, remove_checks=True)
    sol = wing.static_solve(
        prescribed_dofs=tuple(range(6)), horseshoe=True
    )
    # return average cl and cd over the span
    return jnp.mean(sol.aero.cl[0]), jnp.mean(sol.aero.cd[0])


# angle of attack sweep
alphas_deg = jnp.linspace(-5.0, 15.0, 21)

# build each base wing once, one with and one without polars
wing_shell_baseline = build_wing_base(None, None)
wing_shell_corrected = build_wing_base(naca0012_polar, lookup_polar)

# create solutions of each case for a range of angles of attack
sweep_baseline = jax.jit(vmap(lambda a: solve_at_alpha(wing_shell_baseline, a)))
sweep_corrected = jax.jit(vmap(lambda a: solve_at_alpha(wing_shell_corrected, a)))

cl_baseline, cd_baseline = sweep_baseline(alphas_deg)
cl_corrected, cd_corrected = sweep_corrected(alphas_deg)
def solve_at_alpha(wing: CoupledAeroelastic, alpha_deg: Array) -> tuple[Array, Array]:
    set_wing_alpha(wing, alpha_deg, remove_checks=True)
    sol = wing.static_solve(
        prescribed_dofs=tuple(range(6)), horseshoe=True
    )
    # return average cl and cd over the span
    return jnp.mean(sol.aero.cl[0]), jnp.mean(sol.aero.cd[0])


# angle of attack sweep
alphas_deg = jnp.linspace(-5.0, 15.0, 21)

# build each base wing once, one with and one without polars
wing_shell_baseline = build_wing_base(None, None)
wing_shell_corrected = build_wing_base(naca0012_polar, lookup_polar)

# create solutions of each case for a range of angles of attack
sweep_baseline = jax.jit(vmap(lambda a: solve_at_alpha(wing_shell_baseline, a)))
sweep_corrected = jax.jit(vmap(lambda a: solve_at_alpha(wing_shell_corrected, a)))

cl_baseline, cd_baseline = sweep_baseline(alphas_deg)
cl_corrected, cd_corrected = sweep_corrected(alphas_deg)

The polar-corrected wing follows the tabulated data; the uncorrected UVLM lift grows linearly and has no drag at all. Note that due to finite wing effects, the wing-average lift coefficient is lower than the 2D tabulated value.

In [ ]:

Copied!





fig, axes = plt.subplots(1, 2, figsize=(10.0, 4.0))

axes[0].plot(naca0012_polar.alpha_deg, naca0012_polar.cl, "x", label="tabulated polar")
axes[0].plot(alphas_deg, cl_baseline, "--", label="UVLM (uncorrected)")
axes[0].plot(alphas_deg, cl_corrected, "-", label="UVLM (polar-corrected)")
axes[0].set_xlim(alphas_deg.min(), alphas_deg.max())
axes[0].set_xlabel("Angle of attack [deg]")
axes[0].set_ylabel("$C_L$")

axes[1].plot(naca0012_polar.alpha_deg, naca0012_polar.cd, "x", label="tabulated polar")
axes[1].plot(alphas_deg, cd_baseline, "--", label="UVLM (uncorrected)")
axes[1].plot(alphas_deg, cd_corrected, "-", label="UVLM (polar-corrected)")
axes[1].set_xlim(alphas_deg.min(), alphas_deg.max())
axes[1].set_xlabel("Angle of attack [deg]")
axes[1].set_ylabel("$C_D$")
axes[1].legend()

fig.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 2, figsize=(10.0, 4.0))

axes[0].plot(naca0012_polar.alpha_deg, naca0012_polar.cl, "x", label="tabulated polar")
axes[0].plot(alphas_deg, cl_baseline, "--", label="UVLM (uncorrected)")
axes[0].plot(alphas_deg, cl_corrected, "-", label="UVLM (polar-corrected)")
axes[0].set_xlim(alphas_deg.min(), alphas_deg.max())
axes[0].set_xlabel("Angle of attack [deg]")
axes[0].set_ylabel("$C_L$")

axes[1].plot(naca0012_polar.alpha_deg, naca0012_polar.cd, "x", label="tabulated polar")
axes[1].plot(alphas_deg, cd_baseline, "--", label="UVLM (uncorrected)")
axes[1].plot(alphas_deg, cd_corrected, "-", label="UVLM (polar-corrected)")
axes[1].set_xlim(alphas_deg.min(), alphas_deg.max())
axes[1].set_xlabel("Angle of attack [deg]")
axes[1].set_ylabel("$C_D$")
axes[1].legend()

fig.tight_layout()
plt.show()

Spanwise effect at a fixed angle of attack¶

Compare the spanwise sectional lift and the resulting wingtip deflection at a single angle of attack just past the tabulated stall break. This results in the polar-corrected model having reduced deformation.

In [ ]:

Copied!





alpha_showcase = 15.0  # angle of attack where stall occurs

wing_baseline = build_wing(polar_data=None, polar_function=None, alpha_deg=alpha_showcase)
sol_baseline = wing_baseline.static_solve(prescribed_dofs=jnp.arange(6), horseshoe=True)

wing_corrected = build_wing(
    polar_data=naca0012_polar, polar_function=lookup_polar, alpha_deg=alpha_showcase
)
sol_corrected = wing_corrected.static_solve(prescribed_dofs=jnp.arange(6), horseshoe=True)

# mid-chord tip point, z-component
tip_z_baseline = float(sol_baseline.aero.zeta_b[0][0, -1, 2])
tip_z_corrected = float(sol_corrected.aero.zeta_b[0][0, -1, 2])

print(f"Tip deflection at alpha={alpha_showcase:.0f} deg:")
print(f"Uncorrected:     {tip_z_baseline:.3f} m")
print(f"Polar-corrected: {tip_z_corrected:.3f} m")
alpha_showcase = 15.0  # angle of attack where stall occurs

wing_baseline = build_wing(polar_data=None, polar_function=None, alpha_deg=alpha_showcase)
sol_baseline = wing_baseline.static_solve(prescribed_dofs=jnp.arange(6), horseshoe=True)

wing_corrected = build_wing(
    polar_data=naca0012_polar, polar_function=lookup_polar, alpha_deg=alpha_showcase
)
sol_corrected = wing_corrected.static_solve(prescribed_dofs=jnp.arange(6), horseshoe=True)

# mid-chord tip point, z-component
tip_z_baseline = float(sol_baseline.aero.zeta_b[0][0, -1, 2])
tip_z_corrected = float(sol_corrected.aero.zeta_b[0][0, -1, 2])

print(f"Tip deflection at alpha={alpha_showcase:.0f} deg:")
print(f"Uncorrected:     {tip_z_baseline:.3f} m")
print(f"Polar-corrected: {tip_z_corrected:.3f} m")

The per-strip cl differs from the 2D tabulated value shown as a dashed line, since it also reflects the finite wing's spanwise-varying induced downwash and any twist. Note that the case without polar corrections has a notable drop in sectional lift coefficient at the root.

In [ ]:

Copied!





y_nodes = jnp.linspace(0.0, b_ref, n + 1)
y_strip = 0.5 * (y_nodes[:-1] + y_nodes[1:])  # strip midpoints
cl_2d = jnp.interp(alpha_showcase, naca0012_polar.alpha_deg, naca0012_polar.cl)

fig, ax = plt.subplots(figsize=(6.0, 4.0))
ax.plot(y_strip, sol_baseline.aero.cl[0], "--", label="uncorrected")
ax.plot(y_strip, sol_corrected.aero.cl[0], "-", label="polar-corrected")
ax.set_xlabel("Span station [m]")
ax.set_ylabel("Sectional $c_l$")
ax.set_title(rf"Spanwise lift distribution, $\alpha$={alpha_showcase:.0f}$^\circ$")
ax.legend()
fig.tight_layout()
plt.show()
y_nodes = jnp.linspace(0.0, b_ref, n + 1)
y_strip = 0.5 * (y_nodes[:-1] + y_nodes[1:])  # strip midpoints
cl_2d = jnp.interp(alpha_showcase, naca0012_polar.alpha_deg, naca0012_polar.cl)

fig, ax = plt.subplots(figsize=(6.0, 4.0))
ax.plot(y_strip, sol_baseline.aero.cl[0], "--", label="uncorrected")
ax.plot(y_strip, sol_corrected.aero.cl[0], "-", label="polar-corrected")
ax.set_xlabel("Span station [m]")
ax.set_ylabel("Sectional $c_l$")
ax.set_title(rf"Spanwise lift distribution, $\alpha$={alpha_showcase:.0f}$^\circ$")
ax.legend()
fig.tight_layout()
plt.show()