"""Demonstrates running Ptera Software's
AeroelasticUnsteadyRingVortexLatticeMethodSolver across a sweep of parameters.
The wing-tip twist histories from each run are overlaid in a plot.
"""
import matplotlib.pyplot as plt
import numpy as np
import pterasoftware as ps
# The element to extract from each entry of the deformation angle time series. The
# example's Wing has 16 spanwise panels, so index 16 corresponds to the wing-tip
# WingCrossSection.
TIP_INDEX = 16
# These are the default values used when a parameter is not being swept.
DEFAULT_K_RAD = 1000.0
DEFAULT_B_RAD = 1000.0
DEFAULT_DENSITY = 6.0
# The prescribed flapping motion's x-component amplitude (degrees), period (seconds),
# and phase (degrees). The two aeroelastic wing movements are constructed with these
# values, and the plotted flap-position overlay is computed from them.
FLAP_AMPLITUDE = 15.0
FLAP_PERIOD = 1.0
FLAP_PHASE = 169.0
# Populate exactly ONE of these lists to sweep that parameter while holding the others
# at their defaults. Leave the other two as empty lists.
K_VALUES_RAD: list[float] = [100.0, 1000.0, 10000.0, 20000.0]
B_VALUES_RAD: list[float] = []
DENSITY_VALUES: list[float] = []
def run_aeroelastic(
spring_constant_rad: float = DEFAULT_K_RAD,
damping_constant_rad: float = DEFAULT_B_RAD,
wing_density: float = DEFAULT_DENSITY,
) -> tuple[list[np.ndarray], ps.problems.AeroelasticUnsteadyProblem]:
"""Run the aeroelastic solver and return the first Wing's deformation angle time
series.
:param spring_constant_rad: The torsional spring stiffness value.
:param damping_constant_rad: The damping constant value.
:param wing_density: The wing density per unit height (kg/m^2).
:return: A tuple of the first Wing's cumulative torsional deformation angle time
series (its entry of listDeformationAnglesYRad_Wcsp_to_Wcs_ixyz) and the solved
AeroelasticUnsteadyProblem.
"""
# Initialize the WingCrossSection parameters.
num_spanwise_panels = 2
Lp_Wcsp_Lpp_Offsets = (0.1, 0.5, 0.0)
cross_section_chords = [1.75, 1.75, 1.75, 1.75, 1.65, 1.55, 1.4, 1.2, 1.0]
wing_cross_sections = []
for i in range(len(cross_section_chords)):
wing_cross_sections.append(
ps.geometry.wing_cross_section.WingCrossSection(
num_spanwise_panels=(
num_spanwise_panels if i < len(cross_section_chords) - 1 else None
),
chord=cross_section_chords[i],
Lp_Wcsp_Lpp=Lp_Wcsp_Lpp_Offsets if i > 0 else (0.0, 0.0, 0.0),
angles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
control_surface_symmetry_type="symmetric",
control_surface_hinge_point=0.75,
control_surface_deflection=0.0,
spanwise_spacing=(
"uniform" if i < len(cross_section_chords) - 1 else None
),
airfoil=ps.geometry.airfoil.Airfoil(
name="naca2412",
outline_A_lp=None,
resample=True,
n_points_per_side=400,
),
)
)
wing_1 = ps.geometry.wing.Wing(
wing_cross_sections=wing_cross_sections,
name="Main Wing",
Ler_Gs_Cgs=(0.0, 0.5, 0.0),
angles_Gs_to_Wn_ixyz=(0.0, 0.0, 0.0),
symmetric=True,
mirror_only=False,
symmetryNormal_G=(0.0, 1.0, 0.0),
symmetryPoint_G_Cg=(0.0, 0.0, 0.0),
explode_into_strips=True,
num_chordwise_panels=6,
chordwise_spacing="uniform",
)
example_airplane = ps.geometry.airplane.Airplane(
wings=[
wing_1,
ps.geometry.wing.Wing(
wing_cross_sections=[
ps.geometry.wing_cross_section.WingCrossSection(
num_spanwise_panels=8,
chord=1.5,
Lp_Wcsp_Lpp=(0.0, 0.0, 0.0),
angles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
control_surface_symmetry_type="symmetric",
control_surface_hinge_point=0.75,
control_surface_deflection=0.0,
spanwise_spacing="uniform",
airfoil=ps.geometry.airfoil.Airfoil(
name="naca0012",
outline_A_lp=None,
resample=True,
n_points_per_side=400,
),
),
ps.geometry.wing_cross_section.WingCrossSection(
num_spanwise_panels=None,
chord=1.0,
Lp_Wcsp_Lpp=(0.5, 2.0, 0.0),
angles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
control_surface_symmetry_type="symmetric",
control_surface_hinge_point=0.75,
control_surface_deflection=0.0,
spanwise_spacing=None,
airfoil=ps.geometry.airfoil.Airfoil(
name="naca0012",
outline_A_lp=None,
resample=True,
n_points_per_side=400,
),
),
],
name="V-Tail",
Ler_Gs_Cgs=(5.0, 0.0, 0.0),
angles_Gs_to_Wn_ixyz=(0.0, -5.0, 0.0),
symmetric=True,
mirror_only=False,
symmetryNormal_G=(0.0, 1.0, 0.0),
symmetryPoint_G_Cg=(0.0, 0.0, 0.0),
explode_into_strips=False,
num_chordwise_panels=6,
chordwise_spacing="uniform",
),
],
name="Example Airplane",
Cg_GP1_CgP1=(0.0, 0.0, 0.0),
weight=0.0,
c_ref=None,
b_ref=None,
)
# Define the WingCrossSectionMovement parameters.
dephase_x = 0.0
period_x = 0.0
amplitude_x = 0.0
dephase_y = 0.0
period_y = 0.0
amplitude_y = 0.0
dephase_z = 0.0
period_z = 0.0
amplitude_z = 0.0
main_wcs_movements_list = []
reflected_wcs_movements_list = []
for i in range(len(example_airplane.wings[0].wing_cross_sections)):
if i == 0:
wcs_movement = ps.movements.aeroelastic_wing_cross_section_movement.AeroelasticWingCrossSectionMovement(
base_wing_cross_section=example_airplane.wings[0].wing_cross_sections[
i
],
)
main_wcs_movements_list.append(wcs_movement)
reflected_wcs_movements_list.append(wcs_movement)
else:
wcs_movement = ps.movements.aeroelastic_wing_cross_section_movement.AeroelasticWingCrossSectionMovement(
base_wing_cross_section=example_airplane.wings[0].wing_cross_sections[
i
],
ampLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
periodLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
spacingLp_Wcsp_Lpp=("sine", "sine", "sine"),
phaseLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
ampAngles_Wcsp_to_Wcs_ixyz=(
amplitude_x,
amplitude_y,
amplitude_z,
),
periodAngles_Wcsp_to_Wcs_ixyz=(period_x, period_y, period_z),
spacingAngles_Wcsp_to_Wcs_ixyz=("sine", "sine", "sine"),
phaseAngles_Wcsp_to_Wcs_ixyz=(dephase_x, dephase_y, dephase_z),
)
main_wcs_movements_list.append(wcs_movement)
reflected_wcs_movements_list.append(wcs_movement)
v_tail_root_wcs_movement = (
ps.movements.wing_cross_section_movement.WingCrossSectionMovement(
base_wing_cross_section=example_airplane.wings[2].wing_cross_sections[0],
ampLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
periodLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
spacingLp_Wcsp_Lpp=("sine", "sine", "sine"),
phaseLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
ampAngles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
periodAngles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
spacingAngles_Wcsp_to_Wcs_ixyz=("sine", "sine", "sine"),
phaseAngles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
)
)
v_tail_tip_wcs_movement = (
ps.movements.wing_cross_section_movement.WingCrossSectionMovement(
base_wing_cross_section=example_airplane.wings[2].wing_cross_sections[1],
ampLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
periodLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
spacingLp_Wcsp_Lpp=("sine", "sine", "sine"),
phaseLp_Wcsp_Lpp=(0.0, 0.0, 0.0),
ampAngles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
periodAngles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
spacingAngles_Wcsp_to_Wcs_ixyz=("sine", "sine", "sine"),
phaseAngles_Wcsp_to_Wcs_ixyz=(0.0, 0.0, 0.0),
)
)
main_wing_movement = ps.movements.aeroelastic_wing_movement.AeroelasticWingMovement(
base_wing=example_airplane.wings[0],
wing_cross_section_movements=main_wcs_movements_list,
ampLer_Gs_Cgs=(0.0, 0.0, 0.0),
periodLer_Gs_Cgs=(0.0, 0.0, 0.0),
spacingLer_Gs_Cgs=("sine", "sine", "sine"),
phaseLer_Gs_Cgs=(0.0, 0.0, 0.0),
ampAngles_Gs_to_Wn_ixyz=(FLAP_AMPLITUDE, 0.0, 0.0),
periodAngles_Gs_to_Wn_ixyz=(FLAP_PERIOD, 0.0, 0.0),
spacingAngles_Gs_to_Wn_ixyz=("sine", "sine", "sine"),
phaseAngles_Gs_to_Wn_ixyz=(FLAP_PHASE, 0.0, 0.0),
)
reflected_main_wing_movement = (
ps.movements.aeroelastic_wing_movement.AeroelasticWingMovement(
base_wing=example_airplane.wings[1],
wing_cross_section_movements=reflected_wcs_movements_list,
ampLer_Gs_Cgs=(0.0, 0.0, 0.0),
periodLer_Gs_Cgs=(0.0, 0.0, 0.0),
spacingLer_Gs_Cgs=("sine", "sine", "sine"),
phaseLer_Gs_Cgs=(0.0, 0.0, 0.0),
ampAngles_Gs_to_Wn_ixyz=(FLAP_AMPLITUDE, 0.0, 0.0),
periodAngles_Gs_to_Wn_ixyz=(FLAP_PERIOD, 0.0, 0.0),
spacingAngles_Gs_to_Wn_ixyz=("sine", "sine", "sine"),
phaseAngles_Gs_to_Wn_ixyz=(FLAP_PHASE, 0.0, 0.0),
)
)
v_tail_wing_movement = ps.movements.wing_movement.WingMovement(
base_wing=example_airplane.wings[2],
wing_cross_section_movements=[
v_tail_root_wcs_movement,
v_tail_tip_wcs_movement,
],
ampLer_Gs_Cgs=(0.0, 0.0, 0.0),
periodLer_Gs_Cgs=(0.0, 0.0, 0.0),
spacingLer_Gs_Cgs=("sine", "sine", "sine"),
phaseLer_Gs_Cgs=(0.0, 0.0, 0.0),
ampAngles_Gs_to_Wn_ixyz=(0.0, 0.0, 0.0),
periodAngles_Gs_to_Wn_ixyz=(0.0, 0.0, 0.0),
spacingAngles_Gs_to_Wn_ixyz=("sine", "sine", "sine"),
phaseAngles_Gs_to_Wn_ixyz=(0.0, 0.0, 0.0),
)
example_airplane_movement = (
ps.movements.aeroelastic_airplane_movement.AeroelasticAirplaneMovement(
base_airplane=example_airplane,
wing_movements=[
main_wing_movement,
reflected_main_wing_movement,
v_tail_wing_movement,
],
ampCg_GP1_CgP1=(0.0, 0.0, 0.0),
periodCg_GP1_CgP1=(0.0, 0.0, 0.0),
spacingCg_GP1_CgP1=("sine", "sine", "sine"),
phaseCg_GP1_CgP1=(0.0, 0.0, 0.0),
)
)
example_operating_point = ps.operating_point.OperatingPoint(
rho=1.225, vCg__E=10.0, alpha=0.0, beta=0.0, externalFX_W=0.0, nu=15.06e-6
)
example_operating_point_movement = (
ps.movements.operating_point_movement.OperatingPointMovement(
base_operating_point=example_operating_point,
ampVCg__E=0.0,
periodVCg__E=0.0,
spacingVCg__E="sine",
)
)
example_movement = ps.movements.aeroelastic_movement.AeroelasticMovement(
airplane_movements=[example_airplane_movement],
operating_point_movement=example_operating_point_movement,
delta_time=0.03,
num_steps=100,
)
example_problem = ps.problems.AeroelasticUnsteadyProblem(
movement=example_movement,
wing_density=wing_density,
spring_constant_rad=spring_constant_rad,
damping_constant_rad=damping_constant_rad,
step_discards=5,
)
example_solver = ps.aeroelastic_unsteady_ring_vortex_lattice_method.AeroelasticUnsteadyRingVortexLatticeMethodSolver(
aeroelastic_unsteady_problem=example_problem,
)
example_solver.run(
prescribed_wake=True,
)
solved_problem = example_solver.unsteady_problem
assert isinstance(solved_problem, ps.problems.AeroelasticUnsteadyProblem)
return solved_problem.listDeformationAnglesYRad_Wcsp_to_Wcs_ixyz[0], solved_problem
# Determine which parameter is being swept.
active_sweeps = sum(1 for v in (K_VALUES_RAD, B_VALUES_RAD, DENSITY_VALUES) if v)
if active_sweeps > 1:
raise ValueError(
"Only one of K_VALUES_RAD, B_VALUES_RAD, or DENSITY_VALUES should be non-empty."
)
if active_sweeps == 0:
raise ValueError(
"At least one of K_VALUES_RAD, B_VALUES_RAD, or DENSITY_VALUES must be "
"non-empty."
)
if K_VALUES_RAD:
sweep_values = K_VALUES_RAD
sweep_name = "Spring Constant"
sweep_symbol = "k"
sweep_kwarg = "spring_constant_rad"
elif B_VALUES_RAD:
sweep_values = B_VALUES_RAD
sweep_name = "Damping Constant"
sweep_symbol = "b"
sweep_kwarg = "damping_constant_rad"
else:
sweep_values = DENSITY_VALUES
sweep_name = "Wing Density"
sweep_symbol = "rho"
sweep_kwarg = "wing_density"
# Run for each swept value and collect the wing-tip twist history.
results = {}
flap_angle = None
for val in sweep_values:
print(f"Running with {sweep_symbol}={val}...")
deformationAnglesYRad_Wcsp_to_Wcs_ixyz, problem = run_aeroelastic(
**{sweep_kwarg: val}
)
# Extract the wing-tip torsional angle from every time series entry (the first entry
# is the pre-solve seed state, and each later entry is one time step's state),
# converting it from radians to degrees.
tip_twist = np.rad2deg(
np.array(deformationAnglesYRad_Wcsp_to_Wcs_ixyz)[:, TIP_INDEX]
).tolist()
results[val] = tip_twist
print(f" Completed {sweep_symbol}={val}, {len(tip_twist)} entries")
# Compute the prescribed flap angle once (it is the same for all runs).
if flap_angle is None:
phase_rad = np.deg2rad(FLAP_PHASE)
dt = problem.movement.delta_time
num_steps = len(tip_twist)
t = np.arange(num_steps) * dt
flap_angle = (
FLAP_AMPLITUDE * np.sin((2 * np.pi / FLAP_PERIOD) * t + phase_rad)
).tolist()
# Create an overlay plot of the wing-tip twist histories for all swept values.
assert flap_angle is not None
plt.figure(figsize=(12, 6), dpi=200)
for val, curve in results.items():
plt.plot(range(len(curve)), curve, label=f"{sweep_symbol}={val}")
plt.plot(
range(len(flap_angle)),
flap_angle,
label="Flap Position",
color="black",
linestyle="--",
)
plt.xlabel("Step")
plt.ylabel("Angle (degrees)")
plt.title(f"Wing-Tip Twist - Varying {sweep_name}")
plt.legend()
plt.grid(True)
filename = f"wing_tip_twist_{sweep_name.lower().replace(' ', '_')}.png"
plt.savefig(filename)
plt.show()
