Swept Pazy wing time-domain LCO¶

Perform a time-domain limit-cycle oscillation computation for the swept Pazy wing, plotting the wingtip angle of attack over time.

Imports

In [ ]:

from jax import numpy as jnp
from matplotlib import pyplot as plt

from flapjax.aero.flowfields import ConstantFlowField
from flapjax.models.pazy.swept.swept_pazy_wing import generate_swept_pazy_wing
from flapjax.structure.constraints import SpringDamper
from flapjax.utils.print_utils import set_verbosity

# this is essential for a large number of timesteps as Jupyter/jax.debug.print() causes zmq.error.ZMQError: Too many open files. Not an issue for regular Python scripts
set_verbosity("warning")

from jax import numpy as jnp from matplotlib import pyplot as plt from flapjax.aero.flowfields import ConstantFlowField from flapjax.models.pazy.swept.swept_pazy_wing import generate_swept_pazy_wing from flapjax.structure.constraints import SpringDamper from flapjax.utils.print_utils import set_verbosity # this is essential for a large number of timesteps as Jupyter/jax.debug.print() causes zmq.error.ZMQError: Too many open files. Not an issue for regular Python scripts set_verbosity("warning")

Case parameters: freestream velocity, angle of attack, tip force/moment perturbation, and total physical simulation time.

In [ ]:

# fig 23
q_inf = 960.0
model = "LE_CORRECTED"
alpha = jnp.deg2rad(10.0)  # angle of attack
physical_time = 2.5  # time to run the simulation for, in seconds

u_inf_mag = jnp.sqrt(2.0 / 1.225 * q_inf)  # freestream velocity

# torsional moment to apply to the tip for the initial static solve. This leads to an overprediction in deformation compared to the actual aeroelastic equilbrium, giving the system an initial pertubation
m_tip = 2.0

# fig 23 q_inf = 960.0 model = "LE_CORRECTED" alpha = jnp.deg2rad(10.0) # angle of attack physical_time = 2.5 # time to run the simulation for, in seconds u_inf_mag = jnp.sqrt(2.0 / 1.225 * q_inf) # freestream velocity # torsional moment to apply to the tip for the initial static solve. This leads to an overprediction in deformation compared to the actual aeroelastic equilbrium, giving the system an initial pertubation m_tip = 2.0

Build the wing at the chosen freestream and AoA.

In [ ]:

case = generate_swept_pazy_wing(
    flowfield=ConstantFlowField(
        u_inf=jnp.array((u_inf_mag, 0.0, 0.0)),
        rho=1.225,
        relative_motion=True,
    ),
    tip_mass=model,
    aoa=alpha,
    gravity=jnp.array((0.0, -9.81, 0.0)),
    m=8,
    m_star=80,
    node_multiplier=1,
)

case = generate_swept_pazy_wing( flowfield=ConstantFlowField( u_inf=jnp.array((u_inf_mag, 0.0, 0.0)), rho=1.225, relative_motion=True, ), tip_mass=model, aoa=alpha, gravity=jnp.array((0.0, -9.81, 0.0)), m=8, m_star=80, node_multiplier=1, )

Add stiffness-proportional Rayleigh damping targeting 0.5% structural damping at the LCO frequency (~31 Hz).

In [ ]:

lco_freq_hz = 31.0
struct_damping_ratio = 0.005
case.structure.beta_k = float(struct_damping_ratio / (jnp.pi * lco_freq_hz))

lco_freq_hz = 31.0 struct_damping_ratio = 0.005 case.structure.beta_k = float(struct_damping_ratio / (jnp.pi * lco_freq_hz))

Attach a 6-DOF spring/damper at the root node so the wing is elastically restrained rather than fully clamped.

In [ ]:

# root spring/damper: soft in roll/pitch, stiff in translation and yaw.
# damping applied to all 6 degrees of freedom
k_roll = 5e2
k_pitch = 5e3
k_translation = 1e8
k_yaw = 1e8
root_zeta = 2e-2

k_root = jnp.diag(
    jnp.array((k_translation, k_translation, k_translation, k_pitch, k_roll, k_yaw))
)
case.structure.constraints = {
    "root_spring": SpringDamper(
        node_index=0,
        k=k_root,
        c=root_zeta * k_root,
        hg_ref=case.structure.hg0[0],
    ),
}

# root spring/damper: soft in roll/pitch, stiff in translation and yaw. # damping applied to all 6 degrees of freedom k_roll = 5e2 k_pitch = 5e3 k_translation = 1e8 k_yaw = 1e8 root_zeta = 2e-2 k_root = jnp.diag( jnp.array((k_translation, k_translation, k_translation, k_pitch, k_roll, k_yaw)) ) case.structure.constraints = { "root_spring": SpringDamper( node_index=0, k=k_root, c=root_zeta * k_root, hg_ref=case.structure.hg0[0], ), }

Time-step size is set by the aerodynamic panelling; derive the number of steps required for the requested physical time.

In [ ]:

dt = case.aero.dt
n_tstep = int(physical_time / dt)
print(f"dt={float(dt):.4f} s, n_tstep={n_tstep}")

dt = case.aero.dt n_tstep = int(physical_time / dt) print(f"dt={float(dt):.4f} s, n_tstep={n_tstep}")

Assemble the tip follower force / moment vector used to perturb the static solution.

In [ ]:

f_ext = jnp.zeros((case.structure.n_nodes, 6)).at[-1, 3].set(m_tip)

f_ext = jnp.zeros((case.structure.n_nodes, 6)).at[-1, 3].set(m_tip)

Static equilibrium with the tip perturbation applied. This provides the initial condition for the time-domain solve.

In [ ]:

# root DOFs are attached to the spring/damper, and so we don't eliminate them in prescribed_dofs
static_sol = case.static_solve(
    prescribed_dofs=(), horseshoe=False, f_ext_follower=f_ext
)

# root DOFs are attached to the spring/damper, and so we don't eliminate them in prescribed_dofs static_sol = case.static_solve( prescribed_dofs=(), horseshoe=False, f_ext_follower=f_ext )

Time-domain dynamic solve starting from the perturbed static equilibrium.

In [ ]:

# initialise the dynamic solve with the static solution; root DOFs are free (spring-restrained)
dynamic_sol = case.dynamic_solve(
    init_case=static_sol, prescribed_dofs=(), n_tstep=n_tstep
)

# initialise the dynamic solve with the static solution; root DOFs are free (spring-restrained) dynamic_sol = case.dynamic_solve( init_case=static_sol, prescribed_dofs=(), n_tstep=n_tstep )

Uncomment to write files for visualisation in ParaView.

In [ ]:

# plot_dirs = dynamic_sol.plot("./out/")

# plot_dirs = dynamic_sol.plot("./out/")

Plot beam tip vertical deflection over time.

In [ ]:

tip_norm = dynamic_sol.aero.nc[0][:, 0, -1, :]  # [n_tstep, 3]
tip_norm /= jnp.linalg.norm(tip_norm, axis=-1, keepdims=True)

tip_alpha = jnp.rad2deg(jnp.asin(-tip_norm[:, 0]))  # [n_tstep]
t = jnp.arange(n_tstep) * dt

t_plot_start = 2.2
i_ts_start = int(t_plot_start // dt)
fig, ax = plt.subplots()
ax.plot(t[i_ts_start:], tip_alpha[i_ts_start:])
ax.set_xlabel("Time, s")
ax.set_ylabel("Tip angle of attack, deg")
ax.set_title(f"q_inf={q_inf}")
plt.show()

tip_norm = dynamic_sol.aero.nc[0][:, 0, -1, :] # [n_tstep, 3] tip_norm /= jnp.linalg.norm(tip_norm, axis=-1, keepdims=True) tip_alpha = jnp.rad2deg(jnp.asin(-tip_norm[:, 0])) # [n_tstep] t = jnp.arange(n_tstep) * dt t_plot_start = 2.2 i_ts_start = int(t_plot_start // dt) fig, ax = plt.subplots() ax.plot(t[i_ts_start:], tip_alpha[i_ts_start:]) ax.set_xlabel("Time, s") ax.set_ylabel("Tip angle of attack, deg") ax.set_title(f"q_inf={q_inf}") plt.show()