Swept Pazy wing static deflection sweep¶

Run a grid of swept Pazy cases with varying AoA and dynamic pressure in parallel using JAX's pmap. Only the leading-edge tip-mass variant is exercised, as the LE and TE models share the same beam stiffness.

Set XLA flags to set one host device per case. This must run before JAX is imported anywhere in the kernel, otherwise JAX will initialise the default device count and ignore these flags.

In [ ]:

import os

n_vel = 32  # number of velocities
n_alpha = 6  #number of angles of attack
n_case = n_vel * n_alpha
os.environ["XLA_FLAGS"] = (
        os.environ.get("XLA_FLAGS", "")
        + f" --xla_force_host_platform_device_count={n_case} --xla_cpu_multi_thread_eigen=true"
).strip()

import os n_vel = 32 # number of velocities n_alpha = 6 #number of angles of attack n_case = n_vel * n_alpha os.environ["XLA_FLAGS"] = ( os.environ.get("XLA_FLAGS", "") + f" --xla_force_host_platform_device_count={n_case} --xla_cpu_multi_thread_eigen=true" ).strip()

Imports

In [ ]:

import time
from typing import 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.swept.swept_pazy_wing import generate_swept_pazy_wing
from flapjax.utils.print_utils import set_verbosity

import time from typing import 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.swept.swept_pazy_wing import generate_swept_pazy_wing from flapjax.utils.print_utils import set_verbosity

JAX configuration and verbosity

In [ ]:

set_verbosity("silent")  # suppress console prints as there would be a lot for 1k+ cases

# ensure that the XLA flags have worked
assert jax.device_count() == n_case, (
    f"Expected {n_case} devices, got {jax.device_count()}"
)

set_verbosity("silent") # suppress console prints as there would be a lot for 1k+ cases # ensure that the XLA flags have worked assert jax.device_count() == n_case, ( f"Expected {n_case} devices, got {jax.device_count()}" )

Build the AoA / velocity sweep and stack the individual CoupledAeroelastic cases into a single pytree with a leading device axis.

In [ ]:

rho = 1.225  # freestream density
u_inf_vec = jnp.linspace(1e-7, 80.0, n_vel)  # freestream velocity sweep
alphas = jnp.deg2rad(jnp.array((0.0, 1.0, 3.0, 5.0, 7.0, 10.0)))  # angle of attack sweep

# create a list of all the cases. This means creating an instance of the swept Pazy wing for each condition.
cases = []
for u_inf in u_inf_vec:
    u_inf_mag = u_inf
    for alpha in alphas:
        cases.append(
            generate_swept_pazy_wing(
                flowfield=ConstantFlowField(
                    u_inf=jnp.array(
                        (
                            u_inf_mag,
                            0.0,
                            0.0,
                        )
                    ),
                    rho=rho,
                    relative_motion=True,
                ),
                gravity=jnp.array((0.0, -9.81, 0.0)),
                aoa=alpha,
                m=12,
                node_multiplier=2,
                sweep_angle=10,
                tip_mass="TE_CORRECTED",
            )
        )

# stack the list of cases into a JAX pytree so that we can map across them
stacked_case = jax.tree_util.tree_map(lambda *xs: jnp.stack(xs), *cases)

rho = 1.225 # freestream density u_inf_vec = jnp.linspace(1e-7, 80.0, n_vel) # freestream velocity sweep alphas = jnp.deg2rad(jnp.array((0.0, 1.0, 3.0, 5.0, 7.0, 10.0))) # angle of attack sweep # create a list of all the cases. This means creating an instance of the swept Pazy wing for each condition. cases = [] for u_inf in u_inf_vec: u_inf_mag = u_inf for alpha in alphas: cases.append( generate_swept_pazy_wing( flowfield=ConstantFlowField( u_inf=jnp.array( ( u_inf_mag, 0.0, 0.0, ) ), rho=rho, relative_motion=True, ), gravity=jnp.array((0.0, -9.81, 0.0)), aoa=alpha, m=12, node_multiplier=2, sweep_angle=10, tip_mass="TE_CORRECTED", ) ) # stack the list of cases into a JAX pytree so that we can map across them stacked_case = jax.tree_util.tree_map(lambda *xs: jnp.stack(xs), *cases)

Per-case solve function. Returns the mid-chord tip Z displacement normalised by the semi-span.

In [ ]:

def solve(case_: CoupledAeroelastic) -> tuple[Array, Array]:
    static_sol = case_.static_solve(
        prescribed_dofs=tuple(range(6)),  # prescribe the 6 degress of freedom at the root (clamped)
        horseshoe=True,  # use a horseshoe wake as this means only one wake panel per strip
    )

    # compute the mid-chord deflection from the aerodynamic grid coordinate output zeta_b
    tip_z = 0.5 / 0.55 * (
            static_sol.aero.zeta_b[0][0, -1, 2]
            + static_sol.aero.zeta_b[0][-1, -1, 2]
    )

    tip_norm = static_sol.aero.nc[0][0, -1, :]
    tip_norm /= jnp.linalg.norm(tip_norm)

    tip_alpha = jnp.asin(-tip_norm[0])

    return tip_z, tip_alpha


parallel_func = jax.pmap(solve)  # create a function that maps the solve function across the stacked cases

def solve(case_: CoupledAeroelastic) -> tuple[Array, Array]: static_sol = case_.static_solve( prescribed_dofs=tuple(range(6)), # prescribe the 6 degress of freedom at the root (clamped) horseshoe=True, # use a horseshoe wake as this means only one wake panel per strip ) # compute the mid-chord deflection from the aerodynamic grid coordinate output zeta_b tip_z = 0.5 / 0.55 * ( static_sol.aero.zeta_b[0][0, -1, 2] + static_sol.aero.zeta_b[0][-1, -1, 2] ) tip_norm = static_sol.aero.nc[0][0, -1, :] tip_norm /= jnp.linalg.norm(tip_norm) tip_alpha = jnp.asin(-tip_norm[0]) return tip_z, tip_alpha parallel_func = jax.pmap(solve) # create a function that maps the solve function across the stacked cases

Run the analysis for all cases in parallel and time the execution. The results are reshaped into a grid for plotting.

In [ ]:

t_start = time.time()
tip_z_, tip_alpha_ = parallel_func(stacked_case)  # evaluate the function
jax.block_until_ready([tip_z_, tip_alpha_])
t_end = time.time()

print("Total time: ", t_end - t_start)
print("Time per case: ", (t_end - t_start) / n_case)

tip_z_grid = tip_z_.reshape(n_vel, n_alpha).T  # [n_alpha, n_vel] for plotting
tip_alpha_grid = tip_alpha_.reshape(n_vel, n_alpha).T

# subtract initial discrepancy in z (for an undeformed swept wing with aoa, there will be a nonzero tip z)
tip_z_grid -= tip_z_grid[:, [0]]

t_start = time.time() tip_z_, tip_alpha_ = parallel_func(stacked_case) # evaluate the function jax.block_until_ready([tip_z_, tip_alpha_]) t_end = time.time() print("Total time: ", t_end - t_start) print("Time per case: ", (t_end - t_start) / n_case) tip_z_grid = tip_z_.reshape(n_vel, n_alpha).T # [n_alpha, n_vel] for plotting tip_alpha_grid = tip_alpha_.reshape(n_vel, n_alpha).T # subtract initial discrepancy in z (for an undeformed swept wing with aoa, there will be a nonzero tip z) tip_z_grid -= tip_z_grid[:, [0]]

Data from from "Flutter, Post-Flutter, and Limit-Cycle Oscillations of Very Flexible Swept Wings", Fig. 6a) and 8a)

In [ ]:

tip_z_alpha_0 = jnp.array(((249.99999999999977, -0.002529510961214121),
                           (539.8936170212767, -0.004215851602023535),
                           (957.446808510638, -0.010961214165261413),
                           (1510.6382978723402, 0.001686340640809414),
                           (2175.531914893616, 0.001686340640809414),
                           (2962.765957446808, 0.01180438448566612)))
tip_z_alpha_1 = jnp.array(((239.36170212765978, 0.005059021922428353),
                           (547.8723404255318, 0.012647554806070826),
                           (973.4042553191491, 0.024451939291736946),
                           (1489.3617021276593, 0.021922428330522825),
                           (2159.5744680851058, 0.033726812816188945),
                           (2973.4042553191484, 0.0497470489038786)))
tip_z_alpha_3 = jnp.array(((228.72340425531934, 0.018549747048903886),
                           (537.2340425531913, 0.038785834738617186),
                           (957.446808510638, 0.062394603709949426),
                           (1499.9999999999998, 0.1053962900505902),
                           (2186.1702127659573, 0.1247892074198988),
                           (2946.808510638298, 0.16779089376053968)))
tip_z_alpha_5 = jnp.array(((234.0425531914891, 0.04215851602023607),
                           (531.9148936170211, 0.08178752107925802),
                           (973.4042553191491, 0.12816188870151773),
                           (1515.9574468085104, 0.1736930860033727),
                           (2159.5744680851058, 0.22934232715008435),
                           (2978.723404255318, 0.2782462057335582)))
tip_z_alpha_7 = jnp.array(((228.72340425531934, 0.062394603709949426),
                           (537.2340425531913, 0.10961214165261385),
                           (962.7659574468087, 0.1661045531197302),
                           (1526.5957446808509, 0.23355817875210794),
                           (2159.5744680851058, 0.306070826306914),
                           (3026.5957446808507, 0.38532883642495785)))
tip_z_alpha_10 = jnp.array(((239.36170212765978, 0.08937605396290049),
                            (537.2340425531913, 0.1661045531197302),
                            (989.3617021276593, 0.25042158516020235),
                            (1462.7659574468082, 0.321247892074199),
                            (2132.978723404255, 0.40809443507588533)))
tip_z_plot_data = (tip_z_alpha_0, tip_z_alpha_1, tip_z_alpha_3, tip_z_alpha_5, tip_z_alpha_7, tip_z_alpha_10)

tip_z_alpha_0 = jnp.array(((249.99999999999977, -0.002529510961214121), (539.8936170212767, -0.004215851602023535), (957.446808510638, -0.010961214165261413), (1510.6382978723402, 0.001686340640809414), (2175.531914893616, 0.001686340640809414), (2962.765957446808, 0.01180438448566612))) tip_z_alpha_1 = jnp.array(((239.36170212765978, 0.005059021922428353), (547.8723404255318, 0.012647554806070826), (973.4042553191491, 0.024451939291736946), (1489.3617021276593, 0.021922428330522825), (2159.5744680851058, 0.033726812816188945), (2973.4042553191484, 0.0497470489038786))) tip_z_alpha_3 = jnp.array(((228.72340425531934, 0.018549747048903886), (537.2340425531913, 0.038785834738617186), (957.446808510638, 0.062394603709949426), (1499.9999999999998, 0.1053962900505902), (2186.1702127659573, 0.1247892074198988), (2946.808510638298, 0.16779089376053968))) tip_z_alpha_5 = jnp.array(((234.0425531914891, 0.04215851602023607), (531.9148936170211, 0.08178752107925802), (973.4042553191491, 0.12816188870151773), (1515.9574468085104, 0.1736930860033727), (2159.5744680851058, 0.22934232715008435), (2978.723404255318, 0.2782462057335582))) tip_z_alpha_7 = jnp.array(((228.72340425531934, 0.062394603709949426), (537.2340425531913, 0.10961214165261385), (962.7659574468087, 0.1661045531197302), (1526.5957446808509, 0.23355817875210794), (2159.5744680851058, 0.306070826306914), (3026.5957446808507, 0.38532883642495785))) tip_z_alpha_10 = jnp.array(((239.36170212765978, 0.08937605396290049), (537.2340425531913, 0.1661045531197302), (989.3617021276593, 0.25042158516020235), (1462.7659574468082, 0.321247892074199), (2132.978723404255, 0.40809443507588533))) tip_z_plot_data = (tip_z_alpha_0, tip_z_alpha_1, tip_z_alpha_3, tip_z_alpha_5, tip_z_alpha_7, tip_z_alpha_10)

In [ ]:

tip_aoa_alpha_0 = jnp.array(((0.04867157792008925, -0.17473096473307415),
                             (250.68857802865128, 0.0280048771417043),
                             (543.1673216529865, 0.02931526577801513),
                             (960.975378421518, 0.09839146674920052),
                             (1515.3034672883311, 0.06055243488928319),
                             (2178.232822364964, 0.14416770977761928),
                             (2966.5439061312445, 0.10737698882675417)))
tip_aoa_alpha_1 = jnp.array(((-0.30326137011735455, 1.0887083187214586),
                             (242.04750096407165, 1.0494715389828126),
                             (551.2542915227838, 0.9970934332057464),
                             (974.6595797521242, 0.9721086898734406),
                             (1489.9680389971654, 1.0147400001747187),
                             (2166.8511610667306, 1.0043317704348862),
                             (2974.660827741301, 0.9676284087264406)))
tip_aoa_alpha_3 = jnp.array(((-0.8274168246411762, 2.9704264004622534),
                             (227.6070181919383, 2.890804690941719),
                             (539.6180404323534, 2.7712348478513986),
                             (957.4860006814021, 2.6252575537664935),
                             (1497.9576657111227, 2.3319800970685964),
                             (2188.7795789034913, 2.2813117364646214),
                             (2949.2804718397474, 2.0831060953039415)))
tip_aoa_alpha_5 = jnp.array(((-1.463891305134382, 5.255369785433219),
                             (232.63516658783547, 4.839751949671091),
                             (530.7785330885612, 4.505066212065808),
                             (973.7722594469659, 4.157588585391785),
                             (1517.0069725155342, 3.9449686692316996),
                             (2157.8731269242426, 3.2354743419665057),
                             (2976.918440163286, 2.862799813800014)))
tip_aoa_alpha_7 = jnp.array(((0.7937211168508043, 7.150541190506794),
                             (229.36293896459355, 6.5870491171100545),
                             (535.8703289325075, 6.225519132298083),
                             (964.9627038434322, 5.783893202078151),
                             (1525.0939423853315, 4.912746836659432),
                             (2160.3591213656987, 4.310754297138735),
                             (3029.670942693584, 3.481315730029989)))
tip_aoa_alpha_10 = jnp.array(((-2.7892558115731845, 10.01342836354923),
                              (236.91826544480242, 9.463427053160594),
                              (537.9332550428126, 8.819614396303953),
                              (992.2712030241274, 7.746381143382723),
                              (1466.0515744007466, 6.874847901318999),
                              (2134.9001421459666, 5.708489695977355)))

tip_aoa_plot_data = (tip_aoa_alpha_0, tip_aoa_alpha_1, tip_aoa_alpha_3, tip_aoa_alpha_5, tip_aoa_alpha_7,
                     tip_aoa_alpha_10)

tip_aoa_alpha_0 = jnp.array(((0.04867157792008925, -0.17473096473307415), (250.68857802865128, 0.0280048771417043), (543.1673216529865, 0.02931526577801513), (960.975378421518, 0.09839146674920052), (1515.3034672883311, 0.06055243488928319), (2178.232822364964, 0.14416770977761928), (2966.5439061312445, 0.10737698882675417))) tip_aoa_alpha_1 = jnp.array(((-0.30326137011735455, 1.0887083187214586), (242.04750096407165, 1.0494715389828126), (551.2542915227838, 0.9970934332057464), (974.6595797521242, 0.9721086898734406), (1489.9680389971654, 1.0147400001747187), (2166.8511610667306, 1.0043317704348862), (2974.660827741301, 0.9676284087264406))) tip_aoa_alpha_3 = jnp.array(((-0.8274168246411762, 2.9704264004622534), (227.6070181919383, 2.890804690941719), (539.6180404323534, 2.7712348478513986), (957.4860006814021, 2.6252575537664935), (1497.9576657111227, 2.3319800970685964), (2188.7795789034913, 2.2813117364646214), (2949.2804718397474, 2.0831060953039415))) tip_aoa_alpha_5 = jnp.array(((-1.463891305134382, 5.255369785433219), (232.63516658783547, 4.839751949671091), (530.7785330885612, 4.505066212065808), (973.7722594469659, 4.157588585391785), (1517.0069725155342, 3.9449686692316996), (2157.8731269242426, 3.2354743419665057), (2976.918440163286, 2.862799813800014))) tip_aoa_alpha_7 = jnp.array(((0.7937211168508043, 7.150541190506794), (229.36293896459355, 6.5870491171100545), (535.8703289325075, 6.225519132298083), (964.9627038434322, 5.783893202078151), (1525.0939423853315, 4.912746836659432), (2160.3591213656987, 4.310754297138735), (3029.670942693584, 3.481315730029989))) tip_aoa_alpha_10 = jnp.array(((-2.7892558115731845, 10.01342836354923), (236.91826544480242, 9.463427053160594), (537.9332550428126, 8.819614396303953), (992.2712030241274, 7.746381143382723), (1466.0515744007466, 6.874847901318999), (2134.9001421459666, 5.708489695977355))) tip_aoa_plot_data = (tip_aoa_alpha_0, tip_aoa_alpha_1, tip_aoa_alpha_3, tip_aoa_alpha_5, tip_aoa_alpha_7, tip_aoa_alpha_10)

Plot of wingtip deflection versus root AoA, colour-coded by angle of attack.

In [ ]:

alphas_deg = jnp.rad2deg(alphas)

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

cmap = plt.get_cmap("viridis")

q_inf_vec = 0.5 * 1.225 * u_inf_vec ** 2

fig, ax = plt.subplots()
for i_aoa, aoa in enumerate(alphas_deg):
    ax.plot(q_inf_vec, tip_z_grid[i_aoa, :], color=cmap(norm(float(aoa))))

for data, aoa in zip(tip_z_plot_data, [0.0, 1.0, 3.0, 5.0, 7.0, 10.0]):
    ax.scatter(data[:, 0], data[:, 1], marker="x", color=cmap(norm(aoa)))

sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
fig.colorbar(sm, ax=ax, label="Root angle of Attack [deg]")
ax.plot([], [], color="b", label="Computational")
ax.scatter([], [], color="b", marker="x", label="Experimental")
ax.legend()
ax.set_xlabel("Dynamic pressure [Pa]")
ax.set_ylabel("Tip Z Displacement [z/b]")
ax.set_title("Swept Pazy (10 deg) Wingtip Z Displacement")
plt.show()

alphas_deg = jnp.rad2deg(alphas) # noinspection PyArgumentList norm = cast(Normalize, Normalize(vmin=float(alphas_deg.min()), vmax=float(alphas_deg.max()))) cmap = plt.get_cmap("viridis") q_inf_vec = 0.5 * 1.225 * u_inf_vec ** 2 fig, ax = plt.subplots() for i_aoa, aoa in enumerate(alphas_deg): ax.plot(q_inf_vec, tip_z_grid[i_aoa, :], color=cmap(norm(float(aoa)))) for data, aoa in zip(tip_z_plot_data, [0.0, 1.0, 3.0, 5.0, 7.0, 10.0]): ax.scatter(data[:, 0], data[:, 1], marker="x", color=cmap(norm(aoa))) sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap) fig.colorbar(sm, ax=ax, label="Root angle of Attack [deg]") ax.plot([], [], color="b", label="Computational") ax.scatter([], [], color="b", marker="x", label="Experimental") ax.legend() ax.set_xlabel("Dynamic pressure [Pa]") ax.set_ylabel("Tip Z Displacement [z/b]") ax.set_title("Swept Pazy (10 deg) Wingtip Z Displacement") plt.show()

Plot of wingtip angle of attack versus root AoA, colour-coded by angle of attack.

In [ ]:

fig, ax = plt.subplots()
for i_aoa, aoa in enumerate(alphas_deg):
    ax.plot(q_inf_vec, jnp.rad2deg(tip_alpha_grid[i_aoa]), color=cmap(norm(float(aoa))))

for data, aoa in zip(tip_aoa_plot_data, [0.0, 1.0, 3.0, 5.0, 7.0, 10.0]):
    ax.scatter(data[:, 0], data[:, 1], marker="x", color=cmap(norm(aoa)))

sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
fig.colorbar(sm, ax=ax, label="Root angle of Attack [deg]")
ax.plot([], [], color="b", label="Computational")
ax.scatter([], [], color="b", marker="x", label="Experimental")
ax.legend()
ax.set_xlabel("Dynamic pressure [Pa]")
ax.set_ylabel("Tip angle of attack [deg]")
ax.set_title("Swept Pazy (10 deg) Wingtip Angle of Attack")
plt.show()

fig, ax = plt.subplots() for i_aoa, aoa in enumerate(alphas_deg): ax.plot(q_inf_vec, jnp.rad2deg(tip_alpha_grid[i_aoa]), color=cmap(norm(float(aoa)))) for data, aoa in zip(tip_aoa_plot_data, [0.0, 1.0, 3.0, 5.0, 7.0, 10.0]): ax.scatter(data[:, 0], data[:, 1], marker="x", color=cmap(norm(aoa))) sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap) fig.colorbar(sm, ax=ax, label="Root angle of Attack [deg]") ax.plot([], [], color="b", label="Computational") ax.scatter([], [], color="b", marker="x", label="Experimental") ax.legend() ax.set_xlabel("Dynamic pressure [Pa]") ax.set_ylabel("Tip angle of attack [deg]") ax.set_title("Swept Pazy (10 deg) Wingtip Angle of Attack") plt.show()