diff --git a/examples/aeroelastic_unsteady_first_order_deformation.py b/examples/aeroelastic_unsteady_first_order_deformation.py new file mode 100644 index 00000000..da40f404 --- /dev/null +++ b/examples/aeroelastic_unsteady_first_order_deformation.py @@ -0,0 +1,356 @@ +"""This is script is an example of how to run Ptera Software's +UnsteadyRingVortexLatticeMethodSolver with a custom Airplane with a non-static +Movement.""" + +# First, import the software's main package. Note that if you wished to import this +# software into another package, you would first install it by running "pip install +# pterasoftware" in your terminal. +import pterasoftware as ps + +# Create an Airplane with our custom geometry. I am going to declare every parameter +# for Airplane, even though most of them have usable default values. This is for +# educational purposes, but keep in mind that it makes the code much longer than it +# needs to be. For details about each parameter, read the detailed class docstring. +# The same caveats apply to the other classes, methods, and functions I call in this +# script. + +# Wing cross section initialization +# offsets for the spacing +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 = [] + +# Initialization loop for our wing cross sections. Here we are defining automatically +# wing cross sections with a variable set of chords. All of the wing cross sections for +# deformation simulation are defined to have num_spanwise_panels=1 (except the wing tip which +# is always None). This is because we deform each strip of wing cross section independently by +# modeling them as torsional springs, and that model only really works if those strips are thin. +# Note that if you want to go thinner for the same base definition, you can increase the number +# of spanwise panels and ensure that in Wing you set the single_step_wing parameter to True, +# which will ensure that the wing is split back up into single strips for deformation. +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="cosine" 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, + ), + ) + ) + +# Primary wing definition. Note that the single_step_wing parameter is set to True, +# which means that the wing will be split into strips for deformation, and each +# strip will be modeled as a torsional spring. +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), + single_step_wing=True, + num_chordwise_panels=6, + chordwise_spacing="uniform", +) + +# Actually generating the airplane. A tail is added to the airplane, but it is not +# split into strips for deformation as currently only the first wing is considered +# for deformation in the codebase. Fututre versions of this feature could allow for +# the deformation of multiple wings. For now, it is convenient to not split the tail +# into single strip wing cross sections as it reduces the number of movement variables +# that need to be defined. +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), + single_step_wing=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, +) + +# The main Wing was defined to have symmetric=True, mirror_only=False, and with a +# symmetry plane offset non-coincident with the Wing's axes yz-plane. Therefore, +# that Wing had type 5 symmetry (see the Wing class documentation for more details on +# symmetry types). Therefore, it was actually split into two Wings, the with the +# second Wing being a reflected version of the first. Therefore, we need to define a +# WingMovement for this reflected Wing. To start, we'll first define the reflected +# main wing's root and tip WingCrossSections' WingCrossSectionMovements. + +# definitions for wing cross section movement 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 + +# A list of wing cross section movements for the main wing +main_wcs_movements_list = [] + +# A list of wing cross section movements for the reflected wing +reflected_wcs_movements_list = [] + +# A loop for defining the movement for the main wing and its reflected counterpart's wing +# cross sections. Each wing cross section has its own AeroelasticWingCrossSectionMovement +# which allows the solver to apply deformation angles at each time step based on the +# aerodynamic loads. +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) + + +# Now define the v-tail's root and tip WingCrossSections' WingCrossSectionMovements. +v_tail_root_wcs_movement = ps.movements.aeroelastic_wing_cross_section_movement.AeroelasticWingCrossSectionMovement( + 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.aeroelastic_wing_cross_section_movement.AeroelasticWingCrossSectionMovement( + 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), +) + +# This dephase parameter is used to make the wing start in a flat position +dephase = 169.0 + +# Now define the main wing's AeroelasticWingMovement, the reflected main wing's +# AeroelasticWingMovement, and the v-tail's AeroelasticWingMovement. +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=(15.0, 0.0, 0.0), + periodAngles_Gs_to_Wn_ixyz=(1.0, 0.0, 0.0), + spacingAngles_Gs_to_Wn_ixyz=("sine", "sine", "sine"), + phaseAngles_Gs_to_Wn_ixyz=(dephase, 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=(15.0, 0.0, 0.0), + periodAngles_Gs_to_Wn_ixyz=(1.0, 0.0, 0.0), + spacingAngles_Gs_to_Wn_ixyz=("sine", "sine", "sine"), + phaseAngles_Gs_to_Wn_ixyz=(dephase, 0.0, 0.0), + ) +) + +v_tail_wing_movement = ps.movements.aeroelastic_wing_movement.AeroelasticWingMovement( + 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), +) + +# Delete the extraneous pointers to the WingCrossSectionMovements, as these are now +# contained within the WingMovements. This is optional, but it can make debugging +# easier. +del v_tail_root_wcs_movement +del v_tail_tip_wcs_movement + +# Now define the example airplane's AeroelasticAirplaneMovement. +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), + ) +) + +# Delete the extraneous pointers to the WingMovements. +del main_wing_movement +del reflected_main_wing_movement +del v_tail_wing_movement + +# Define a new OperatingPoint. +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 +) + +# Define the operating point's AeroelasticOperatingPointMovement. +example_operating_point_movement = ( + ps.movements.aeroelastic_operating_point_movement.AeroelasticOperatingPointMovement( + base_operating_point=example_operating_point, + ampVCg__E=0.0, + periodVCg__E=0.0, + spacingVCg__E="sine", + ) +) + +# Delete the extraneous pointer. +del example_operating_point + +# Define the AeroelasticMovement. This contains the AeroelasticAirplaneMovement and +# the AeroelasticOperatingPointMovement. The delta_time and num_steps must be specified +# explicitly. With a flapping period of 1.0s, 3 cycles at dt=0.03s gives 100 steps. +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, +) + +# Define the AeroelasticUnsteadyProblem. +# The deformation parameters are set here. +# The wing_density, spring_constant and damping_constant are the primary parameters +# you should expect to change. The rest are more for considering numerical issues +# with our integrator and debugging. Plotting the flap cycle can give good data as well. +example_problem = ps.problems.AeroelasticUnsteadyProblem( + movement=example_movement, + wing_density=0.012, + spring_constant=20.0, + damping_constant=1.0, + aero_scaling=1.0, + moment_scaling_factor=1.0, + damping_eps=1e-3, + plot_flap_cycle=False, + custom_spacing_second_derivative=None, +) + +# Define a new solver. We'll create an AeroelasticUnsteadyRingVortexLatticeMethodSolver, +# which requires an AeroelasticUnsteadyProblem. +example_solver = ps.aeroelastic_unsteady_ring_vortex_lattice_method.AeroelasticUnsteadyRingVortexLatticeMethodSolver( + aeroelastic_unsteady_problem=example_problem, +) + +# Delete the extraneous pointer. +del example_problem + +# Run the solver. +example_solver.run( + prescribed_wake=True, +) + +# Call the animate function on the solver. This produces a GIF of the wake being +# shed. The GIF is saved in the same directory as this script. Press "q", +# after orienting the view, to begin the animation. +ps.output.animate( + unsteady_solver=example_solver, + scalar_type="lift", + show_wake_vortices=True, + save=True, +) diff --git a/examples/aeroelastic_unsteady_first_order_plotting.py b/examples/aeroelastic_unsteady_first_order_plotting.py new file mode 100644 index 00000000..18ce9b1c --- /dev/null +++ b/examples/aeroelastic_unsteady_first_order_plotting.py @@ -0,0 +1,389 @@ +"""This script runs the coupled aeroelastic UVLM solver sweeping one parameter +at a time (spring constant k, damping constant b, or wing density) and overlays +the Curve 16 Net Deformation from each run.""" + +import matplotlib.pyplot as plt +import numpy as np + +import pterasoftware as ps + +# The curve index to extract from each Net Deformation result. Curve 16 corresponds +# to the wing-tip spanwise station. +CURVE_INDEX = 16 + +# Default values used when a parameter is not being swept +DEFAULT_K = 1.0 +DEFAULT_B = 40.0 +DEFAULT_DENSITY = 0.012 + +# 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: list[float] = [] +B_VALUES: list[float] = [] +DENSITY_VALUES: list[float] = [0.012, 0.12, 0.3] + + +def run_aeroelastic( + spring_constant: float = DEFAULT_K, + damping_constant: float = DEFAULT_B, + wing_density: float = DEFAULT_DENSITY, +) -> tuple[list, object]: + """Run the aeroelastic solver and return the net deformation data. + + :param spring_constant: The torsional spring stiffness value. + :param damping_constant: The damping constant value. + :param wing_density: The wing density per unit height (kg/m^2). + :return: A tuple of (net_data list, solved problem object). + """ + + # Wing cross section initialization + 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=( + "cosine" 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), + single_step_wing=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), + single_step_wing=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, + ) + + # Wing cross section movement 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.aeroelastic_wing_cross_section_movement.AeroelasticWingCrossSectionMovement( + 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.aeroelastic_wing_cross_section_movement.AeroelasticWingCrossSectionMovement( + 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), + ) + + dephase = 169.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=(15.0, 0.0, 0.0), + periodAngles_Gs_to_Wn_ixyz=(1.0, 0.0, 0.0), + spacingAngles_Gs_to_Wn_ixyz=("sine", "sine", "sine"), + phaseAngles_Gs_to_Wn_ixyz=(dephase, 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=(15.0, 0.0, 0.0), + periodAngles_Gs_to_Wn_ixyz=(1.0, 0.0, 0.0), + spacingAngles_Gs_to_Wn_ixyz=("sine", "sine", "sine"), + phaseAngles_Gs_to_Wn_ixyz=(dephase, 0.0, 0.0), + ) + ) + + v_tail_wing_movement = ( + ps.movements.aeroelastic_wing_movement.AeroelasticWingMovement( + 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.aeroelastic_operating_point_movement.AeroelasticOperatingPointMovement( + 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=spring_constant, + damping_constant=damping_constant, + aero_scaling=1.0, + moment_scaling_factor=1.0, + damping_eps=1e-3, + plot_flap_cycle=False, + custom_spacing_second_derivative=None, + ) + + example_solver = ps.aeroelastic_unsteady_ring_vortex_lattice_method.AeroelasticUnsteadyRingVortexLatticeMethodSolver( + aeroelastic_unsteady_problem=example_problem, + ) + + example_solver.run( + prescribed_wake=True, + ) + + problem = example_solver.unsteady_problem + + # ps.output.animate( + # unsteady_solver=example_solver, + # scalar_type="lift", + # show_wake_vortices=True, + # save=True, + # ) + return problem.net_data, problem + + +# Determine which parameter is being swept +active_sweeps = sum(1 for v in (K_VALUES, B_VALUES, DENSITY_VALUES) if v) +if active_sweeps > 1: + raise ValueError( + "Only one of K_VALUES, B_VALUES, or DENSITY_VALUES should be non-empty." + ) +if active_sweeps == 0: + raise ValueError( + "At least one of K_VALUES, B_VALUES, or DENSITY_VALUES must be non-empty." + ) + +if K_VALUES: + sweep_values = K_VALUES + sweep_name = "Spring Constant" + sweep_symbol = "k" + sweep_kwarg = "spring_constant" +elif B_VALUES: + sweep_values = B_VALUES + sweep_name = "Damping Constant" + sweep_symbol = "b" + sweep_kwarg = "damping_constant" +else: + sweep_values = DENSITY_VALUES + sweep_name = "Wing Density" + sweep_symbol = "ρ" + sweep_kwarg = "wing_density" + +# Run for each swept value and collect Curve 16 of the Net Deformation +results = {} +flap_angle = None +for val in sweep_values: + print(f"Running with {sweep_symbol}={val}...") + net_data, problem = run_aeroelastic(**{sweep_kwarg: val}) + # Extract y-component (torsional angle) for Curve 16 across all time steps + curve_16 = np.array(net_data)[:, CURVE_INDEX, 1].tolist() + results[val] = curve_16 + print(f" Completed {sweep_symbol}={val}, {len(curve_16)} steps") + + # Compute the prescribed flap angle once (it is the same for all runs) + if flap_angle is None: + wm = problem.wing_movement + amp = wm.ampAngles_Gs_to_Wn_ixyz[0] + period = wm.periodAngles_Gs_to_Wn_ixyz[0] + phase = np.deg2rad(wm.phaseAngles_Gs_to_Wn_ixyz[0]) + dt = problem.movement.delta_time + num_steps = len(curve_16) + t = np.arange(num_steps) * dt + flap_angle = (amp * np.sin((2 * np.pi / period) * t + phase)).tolist() + +# Overlay plot of Curve 16 Net Deformation for all swept values +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"Net Deformation (Curve {CURVE_INDEX}) — Varying {sweep_name}") +plt.legend() +plt.grid(True) +filename = f"Net_Deformation_Curve_{CURVE_INDEX}_{sweep_name.replace(' ', '_')}.png" +plt.savefig(filename) +plt.show() diff --git a/pterasoftware/__init__.py b/pterasoftware/__init__.py index 6518a6bc..bb0ee8f8 100644 --- a/pterasoftware/__init__.py +++ b/pterasoftware/__init__.py @@ -12,9 +12,15 @@ **Contains the following modules:** +aeroelastic_unsteady_ring_vortex_lattice_method.py: Contains the +AeroelasticUnsteadyRingVortexLatticeMethodSolver class. + convergence.py: Contains functions for analyzing the convergence of SteadyProblems and UnsteadyProblems. +free_flight_unsteady_ring_vortex_lattice_method.py: Contains the +FreeFlightUnsteadyRingVortexLatticeMethodSolver class. + operating_point.py: Contains the OperatingPoint class. output.py: Contains functions for visualizing geometry and results. @@ -51,12 +57,14 @@ # Lazy imports configuration: modules loaded on first access. _LAZY_MODULES = { + "aeroelastic_unsteady_ring_vortex_lattice_method": "pterasoftware.aeroelastic_unsteady_ring_vortex_lattice_method", "convergence": "pterasoftware.convergence", "output": "pterasoftware.output", "steady_horseshoe_vortex_lattice_method": "pterasoftware.steady_horseshoe_vortex_lattice_method", "steady_ring_vortex_lattice_method": "pterasoftware.steady_ring_vortex_lattice_method", "trim": "pterasoftware.trim", "unsteady_ring_vortex_lattice_method": "pterasoftware.unsteady_ring_vortex_lattice_method", + "free_flight_unsteady_ring_vortex_lattice_method": "pterasoftware.free_flight_unsteady_ring_vortex_lattice_method", } # Lazy callable imports: functions that need special handling. diff --git a/pterasoftware/aeroelastic_unsteady_ring_vortex_lattice_method.py b/pterasoftware/aeroelastic_unsteady_ring_vortex_lattice_method.py new file mode 100644 index 00000000..a9ba40d0 --- /dev/null +++ b/pterasoftware/aeroelastic_unsteady_ring_vortex_lattice_method.py @@ -0,0 +1,217 @@ +"""Contains the AeroelasticUnsteadyRingVortexLatticeMethodSolver class. + +**Contains the following classes:** + +AeroelasticUnsteadyRingVortexLatticeMethodSolver: A subclass of +CoupledUnsteadyRingVortexLatticeMethodSolver that extends the coupled solver with Strip +Leading Edge Point (SLEP) functionality for computing aerodynamic moments about the +strip leading edge. This is used for aeroelastic simulations where wing deformations are +coupled with aerodynamic loads. + +**Contains the following functions:** + +None +""" + +from __future__ import annotations + +import numpy as np + +from . import _functions, problems +from ._coupled_unsteady_ring_vortex_lattice_method import ( + CoupledUnsteadyRingVortexLatticeMethodSolver, +) + + +class AeroelasticUnsteadyRingVortexLatticeMethodSolver( + CoupledUnsteadyRingVortexLatticeMethodSolver +): + """A subclass of CoupledUnsteadyRingVortexLatticeMethodSolver that adds SLEP (Strip + Leading Edge Point) functionality for aeroelastic simulations. + + This solver extends the coupled solver with moment calculations about each panel's + strip leading edge point, which is important for analyzing wing loading and + deformation characteristics relative to the wing root. + + **Key additions over parent CoupledUnsteadyRingVortexLatticeMethodSolver:** + + - Initializes and maintains SLEP index mapping and position arrays - Overrides + _reinitialize_step_arrays_hook() to reinitialize SLEP arrays each step - Overrides + _load_calculation_moment_processing_hook() to compute SLEP-based moments - Computes + bound vortex positions relative to strip leading edge points + """ + + def __init__( + self, + aeroelastic_unsteady_problem: problems.AeroelasticUnsteadyProblem, + ) -> None: + """Initialize the solver for an AeroelasticUnsteadyProblem. + + Sets up the solver infrastructure and initializes SLEP (Strip Leading Edge + Point) related attributes. + + :param aeroelastic_unsteady_problem: The AeroelasticUnsteadyProblem to be + solved. + :return: None + """ + if not isinstance( + aeroelastic_unsteady_problem, problems.AeroelasticUnsteadyProblem + ): + raise TypeError( + "aeroelastic_unsteady_problem must be an " "AeroelasticUnsteadyProblem." + ) + + super().__init__(aeroelastic_unsteady_problem) + + first_steady_problem: problems.SteadyProblem = self._get_steady_problem_at(0) + + # Initialize SLEP (Strip Leading Edge Point) information. For each airplane + # and wing, we track the panel index where each new spanwise strip begins. + # This allows efficient computation of moments about the strip leading edge + # (wing root to tip). + panel_count = 0 + slep_point_indices_list: list[int] = [] + for airplane in first_steady_problem.airplanes: + for wing in airplane.wings: + for wing_cross_section in wing.wing_cross_sections: + # Record the first panel index for this wing cross-section + # (start of strip). + slep_point_indices_list.append(panel_count) + if wing_cross_section.num_spanwise_panels is not None: + panel_count += wing_cross_section.num_spanwise_panels + self.slep_point_indices: np.ndarray = np.array( + slep_point_indices_list, dtype=int + ) + + # The current time step's center bound LineVortex points for the right, + # front, left, and back legs (in the first Airplane's geometry axes, + # relative to the local strip leading edge point). + self.stackCblvpr_GP1_Slep: np.ndarray = np.empty(0, dtype=float) + self.stackCblvpf_GP1_Slep: np.ndarray = np.empty(0, dtype=float) + self.stackCblvpl_GP1_Slep: np.ndarray = np.empty(0, dtype=float) + self.stackCblvpb_GP1_Slep: np.ndarray = np.empty(0, dtype=float) + + # The colocation panel points and the front left panel point (in the first + # Airplane's geometry axes, relative to the local strip leading edge point + # and the first Airplane's CG respectively). + self.stackCpp_GP1_Slep: np.ndarray = np.empty(0, dtype=float) + self.stackFlpp_GP1_CgP1: np.ndarray = np.empty(0, dtype=float) + self.moments_GP1_Slep: np.ndarray = np.empty(0, dtype=float) + self.stack_leading_edge_points: np.ndarray = np.empty(0, dtype=float) + + def _reinitialize_step_arrays_hook(self) -> None: + """Reinitialize SLEP arrays at the start of each time step. + + :return: None + """ + self.stackCblvpr_GP1_Slep = np.zeros((self.num_panels, 3), dtype=float) + self.stackCblvpf_GP1_Slep = np.zeros((self.num_panels, 3), dtype=float) + self.stackCblvpl_GP1_Slep = np.zeros((self.num_panels, 3), dtype=float) + self.stackCblvpb_GP1_Slep = np.zeros((self.num_panels, 3), dtype=float) + self.stackCpp_GP1_Slep = np.zeros((self.num_panels, 3), dtype=float) + self.moments_GP1_Slep = np.zeros((self.num_panels, 3), dtype=float) + self.stackFlpp_GP1_CgP1 = np.zeros((self.num_panels, 3), dtype=float) + + def _load_calculation_moment_processing_hook( + self, + rightLegForces_GP1, + frontLegForces_GP1, + leftLegForces_GP1, + backLegForces_GP1, + unsteady_forces_GP1, + ) -> np.ndarray: + """Override parent to compute moments about both center-of-gravity and SLEP. + + This hook extends the parent class's moment calculation by additionally + computing moments about each panel's Strip Leading Edge Point (SLEP). This is + used for analyzing wing loading and deformation characteristics relative to the + wing root. + + The method: 1. Calls parent's implementation to get CG-based moments 2. Updates + bound vortex positions relative to SLEP points 3. Recalculates all moment + contributions in the SLEP frame 4. Stores SLEP moments in self.moments_GP1_Slep + + :return: moments_GP1_CgP1, a (N,3) ndarray of floats representing the moments + (in the first Airplane's geometry axes, relative to the first Airplane's CG) + on every Panel at the current time step. SLEP moments are stored separately + in self.moments_GP1_Slep. + """ + moments_GP1_CgP1 = super()._load_calculation_moment_processing_hook( + rightLegForces_GP1, + frontLegForces_GP1, + leftLegForces_GP1, + backLegForces_GP1, + unsteady_forces_GP1, + ) + + self._update_bound_vortex_positions_relative_to_slep_points() + + rightLegMoments_GP1_Slep = _functions.numba_1d_explicit_cross( + self.stackCblvpr_GP1_Slep, rightLegForces_GP1 + ) + frontLegMoments_GP1_Slep = _functions.numba_1d_explicit_cross( + self.stackCblvpf_GP1_Slep, frontLegForces_GP1 + ) + leftLegMoments_GP1_Slep = _functions.numba_1d_explicit_cross( + self.stackCblvpl_GP1_Slep, leftLegForces_GP1 + ) + backLegMoments_GP1_Slep = _functions.numba_1d_explicit_cross( + self.stackCblvpb_GP1_Slep, backLegForces_GP1 + ) + + # The unsteady moment is calculated at the collocation point because the + # unsteady force acts on the bound RingVortex, whose center is at the + # collocation point, not at the Panel's centroid. + unsteady_moments_GP1_Slep = _functions.numba_1d_explicit_cross( + self.stackCpp_GP1_Slep, unsteady_forces_GP1 + ) + + self.moments_GP1_Slep = ( + rightLegMoments_GP1_Slep + + frontLegMoments_GP1_Slep + + leftLegMoments_GP1_Slep + + backLegMoments_GP1_Slep + + unsteady_moments_GP1_Slep + ) + + return moments_GP1_CgP1 + + def _update_bound_vortex_positions_relative_to_slep_points(self) -> None: + """Transform bound RingVortex leg center positions from CG-relative to SLEP- + relative. + + For each panel, this method: 1. Gets the front-left panel point (leading edge) + from each panel 2. Maps panels to their corresponding strip's leading edge point + using slep_point_indices 3. Subtracts the SLEP position from all vortex leg + center positions 4. Subtracts the SLEP position from collocation points + + This prepares positions for computing moments about the strip leading edge, + which is important for analyzing local wing loading and deformations. + + :return: None + """ + for panel_num, panel in enumerate(self.panels): + self.stackFlpp_GP1_CgP1[panel_num] = panel.Flpp_GP1_CgP1 + slep_points = self.stackFlpp_GP1_CgP1[self.slep_point_indices] + slep_map = ( + np.searchsorted( + self.slep_point_indices, np.arange(self.num_panels), side="right" + ) + - 1 + ) + self.stack_leading_edge_points = np.array([slep_points[i] for i in slep_map]) + self.stackCblvpr_GP1_Slep = ( + self.stackCblvpr_GP1_CgP1 - self.stack_leading_edge_points + ) + self.stackCblvpf_GP1_Slep = ( + self.stackCblvpf_GP1_CgP1 - self.stack_leading_edge_points + ) + self.stackCblvpl_GP1_Slep = ( + self.stackCblvpl_GP1_CgP1 - self.stack_leading_edge_points + ) + self.stackCblvpb_GP1_Slep = ( + self.stackCblvpb_GP1_CgP1 - self.stack_leading_edge_points + ) + + # Find the collocation point positions relative to the SLEP points. + self.stackCpp_GP1_Slep = self.stackCpp_GP1_CgP1 - self.stack_leading_edge_points diff --git a/pterasoftware/free_flight_unsteady_ring_vortex_lattice_method.py b/pterasoftware/free_flight_unsteady_ring_vortex_lattice_method.py new file mode 100644 index 00000000..57dc9042 --- /dev/null +++ b/pterasoftware/free_flight_unsteady_ring_vortex_lattice_method.py @@ -0,0 +1,50 @@ +"""Contains the FreeFlightUnsteadyRingVortexLatticeMethodSolver class. + +**Contains the following classes:** + +FreeFlightUnsteadyRingVortexLatticeMethodSolver: A stub subclass of +CoupledUnsteadyRingVortexLatticeMethodSolver for free flight simulations. Not yet +implemented. + +**Contains the following functions:** + +None +""" + +from __future__ import annotations + +from . import problems +from ._coupled_unsteady_ring_vortex_lattice_method import ( + CoupledUnsteadyRingVortexLatticeMethodSolver, +) + + +class FreeFlightUnsteadyRingVortexLatticeMethodSolver( + CoupledUnsteadyRingVortexLatticeMethodSolver +): + """A stub subclass of CoupledUnsteadyRingVortexLatticeMethodSolver for free flight + simulations. + + This solver will handle FreeFlightUnsteadyProblems where the operating point is + updated dynamically based on computed aerodynamic forces and moments. + + Not yet implemented. + """ + + def __init__( + self, + free_flight_unsteady_problem: problems.FreeFlightUnsteadyProblem, + ) -> None: + """Initialize the solver for a FreeFlightUnsteadyProblem. + + :param free_flight_unsteady_problem: The FreeFlightUnsteadyProblem to be solved. + :return: None + """ + if not isinstance( + free_flight_unsteady_problem, problems.FreeFlightUnsteadyProblem + ): + raise TypeError( + "free_flight_unsteady_problem must be a " "FreeFlightUnsteadyProblem." + ) + + super().__init__(free_flight_unsteady_problem) diff --git a/pterasoftware/geometry/wing.py b/pterasoftware/geometry/wing.py index 7b0bcd4b..57033ecf 100644 --- a/pterasoftware/geometry/wing.py +++ b/pterasoftware/geometry/wing.py @@ -153,6 +153,7 @@ class Wing: "mirror_only", "symmetryNormal_G", "symmetryPoint_G_Cg", + "single_step_wing", # Set once "_symmetry_type", "_num_spanwise_panels", @@ -190,6 +191,7 @@ def __init__( mirror_only: bool | np.bool_ = False, symmetryNormal_G: None | np.ndarray | Sequence[float | int] = None, symmetryPoint_G_Cg: None | np.ndarray | Sequence[float | int] = None, + single_step_wing: bool | np.bool_ = False, num_chordwise_panels: int = 8, chordwise_spacing: str = "cosine", ) -> None: @@ -250,6 +252,9 @@ def __init__( either are True. For more details on how this parameter interacts with symmetryNormal_G, symmetric, and mirror_only, see the class docstring. The units are meters. The default is None. + :param single_step_wing: Set this to True to have the explode_wing method called + on this Wing during initialization, which will return a NEW Wing where all + panels are broken into single strips for deformation. :param num_chordwise_panels: The number of chordwise panels to be used on this Wing, which must be set to a positive integer. The default is 8. :param chordwise_spacing: The type of spacing between the Wing's chordwise @@ -262,6 +267,11 @@ def __init__( wing_cross_sections = _parameter_validation.non_empty_list_return_list( wing_cross_sections, "wing_cross_sections" ) + self.single_step_wing = _parameter_validation.boolLike_return_bool( + single_step_wing, "single_step_wing" + ) + if single_step_wing: + wing_cross_sections = self.explode_wing(wing_cross_sections) num_wing_cross_sections = len(wing_cross_sections) if num_wing_cross_sections < 2: raise ValueError("wing_cross_sections must contain at least two elements.") @@ -445,6 +455,7 @@ def __deepcopy__(self, memo: dict) -> Wing: # process_wing_symmetry for type 5 symmetry). new_wing.symmetric = self.symmetric new_wing.mirror_only = self.mirror_only + new_wing.single_step_wing = self.single_step_wing # Copy immutable numpy arrays and make them read-only. new_wing._Ler_Gs_Cgs = np.copy(self._Ler_Gs_Cgs) @@ -1503,6 +1514,92 @@ def get_plottable_data( return None + def interpolate_between_wing_cross_sections( + self, + wcs1: wing_cross_section_mod.WingCrossSection, + wcs2: wing_cross_section_mod.WingCrossSection, + first_wcs: bool, + ) -> list[wing_cross_section_mod.WingCrossSection]: + """Wing cross section panels are between the line of wcs1 and wcs2. When + exploding a wing to 1 spanwise panel per cross section, we need to interpolate + the intermediate cross sections. + + :param wcs1: The first WingCrossSection. + :param wcs2: The second WingCrossSection. + :param first_wcs: Whether wcs1 is the first WingCrossSection of the wing. If + True, the method will include a WingCrossSection with the same parameters as + wcs1 in the output list. If False, it will not, since it will have already + been included as the last interpolated WingCrossSection between the previous + pair of WingCrossSections. + :return: A list of WingCrossSections representing the interpolated cross + sections + """ + + interpolated = [] + + if first_wcs: + interpolated.append( + wing_cross_section_mod.WingCrossSection( + num_spanwise_panels=1, + chord=wcs1.chord, + Lp_Wcsp_Lpp=wcs1.Lp_Wcsp_Lpp, + angles_Wcsp_to_Wcs_ixyz=wcs1.angles_Wcsp_to_Wcs_ixyz, + control_surface_symmetry_type=wcs1.control_surface_symmetry_type, + control_surface_hinge_point=wcs1.control_surface_hinge_point, + control_surface_deflection=wcs1.control_surface_deflection, + spanwise_spacing="uniform", + airfoil=wcs1.airfoil, + ) + ) + + N = wcs1.num_spanwise_panels + assert N is not None, "wcs1.num_spanwise_panels must not be None" + + for i in range(N): + t = (i + 1) / N # interpolation parameter between 0 and 1 + + chord = (1 - t) * wcs1.chord + t * wcs2.chord + Lp_Wcsp_Lpp = tuple(np.array(wcs2.Lp_Wcsp_Lpp) / N) + # angles_Wcsp_to_Wcs_ixyz = tuple((1 - t) * np.array(wcs1.angles_Wcsp_to_Wcs_ixyz) + t * np.array(wcs2.angles_Wcsp_to_Wcs_ixyz)) + angles_Wcsp_to_Wcs_ixyz = wcs2.angles_Wcsp_to_Wcs_ixyz / N + is_final_section = wcs2.num_spanwise_panels is None and i == N - 1 + + interpolated.append( + wing_cross_section_mod.WingCrossSection( + num_spanwise_panels=None if is_final_section else 1, + chord=chord, + Lp_Wcsp_Lpp=Lp_Wcsp_Lpp, + angles_Wcsp_to_Wcs_ixyz=angles_Wcsp_to_Wcs_ixyz, + control_surface_symmetry_type=wcs1.control_surface_symmetry_type, + control_surface_hinge_point=wcs1.control_surface_hinge_point, + control_surface_deflection=wcs1.control_surface_deflection, + spanwise_spacing=None if is_final_section else "uniform", + airfoil=wcs1.airfoil, + ) + ) + return interpolated + + def explode_wing( + self, wing_cross_sections: list[wing_cross_section_mod.WingCrossSection] + ) -> list[wing_cross_section_mod.WingCrossSection]: + """Takes a list of WingCrossSections and returns a new list where all cross + sections have num_spanwise_panels = 1. + + :param wing_cross_sections: The list of wing cross sections to explode. + :return: A new list of exploded wing cross sections. + """ + + new_cross_sections = [] + + for i in range(len(wing_cross_sections) - 1): + new_cross_sections.extend( + self.interpolate_between_wing_cross_sections( + wing_cross_sections[i], wing_cross_sections[i + 1], i == 0 + ) + ) + + return new_cross_sections + def _assert_T_not_none(T: np.ndarray | None) -> np.ndarray: """Assert that a transformation matrix is not None and return it. diff --git a/pterasoftware/movements/__init__.py b/pterasoftware/movements/__init__.py index 0f1cfa10..d4328080 100644 --- a/pterasoftware/movements/__init__.py +++ b/pterasoftware/movements/__init__.py @@ -10,8 +10,32 @@ **Contains the following modules:** +aeroelastic_airplane_movement.py: Contains the AeroelasticAirplaneMovement class. + +aeroelastic_movement.py: Contains the AeroelasticMovement class. + +aeroelastic_operating_point_movement.py: Contains the AeroelasticOperatingPointMovement +class. + +aeroelastic_wing_cross_section_movement.py: Contains the +AeroelasticWingCrossSectionMovement class. + +aeroelastic_wing_movement.py: Contains the AeroelasticWingMovement class. + airplane_movement.py: Contains the AirplaneMovement class. +free_flight_airplane_movement.py: Contains the FreeFlightAirplaneMovement class. + +free_flight_movement.py: Contains the FreeFlightMovement class. + +free_flight_operating_point_movement.py: Contains the FreeFlightOperatingPointMovement +class. + +free_flight_wing_cross_section_movement.py: Contains the +FreeFlightWingCrossSectionMovement class. + +free_flight_wing_movement.py: Contains the FreeFlightWingMovement class. + movement.py: Contains the Movement class. operating_point_movement.py: Contains the OperatingPointMovement class. @@ -21,7 +45,17 @@ wing_movement.py: Contains the WingMovement class. """ +import pterasoftware.movements.aeroelastic_airplane_movement +import pterasoftware.movements.aeroelastic_movement +import pterasoftware.movements.aeroelastic_operating_point_movement +import pterasoftware.movements.aeroelastic_wing_cross_section_movement +import pterasoftware.movements.aeroelastic_wing_movement import pterasoftware.movements.airplane_movement +import pterasoftware.movements.free_flight_airplane_movement +import pterasoftware.movements.free_flight_movement +import pterasoftware.movements.free_flight_operating_point_movement +import pterasoftware.movements.free_flight_wing_cross_section_movement +import pterasoftware.movements.free_flight_wing_movement import pterasoftware.movements.movement import pterasoftware.movements.operating_point_movement import pterasoftware.movements.wing_cross_section_movement diff --git a/pterasoftware/movements/aeroelastic_airplane_movement.py b/pterasoftware/movements/aeroelastic_airplane_movement.py index 54f9f799..9ad40f5e 100644 --- a/pterasoftware/movements/aeroelastic_airplane_movement.py +++ b/pterasoftware/movements/aeroelastic_airplane_movement.py @@ -145,7 +145,7 @@ def generate_airplane_at_time_step( self, step: int, delta_time: float | int, - wing_deformation_angles_ixyz: list[np.ndarray] | None = None, + wing_deformation_angles_ixyz: list[np.ndarray | None] | None = None, ) -> geometry.airplane.Airplane: """Creates the Airplane at a single time step, optionally applying structural deformation to each Wing. diff --git a/pterasoftware/movements/aeroelastic_movement.py b/pterasoftware/movements/aeroelastic_movement.py index 2a0db51b..72085ecd 100644 --- a/pterasoftware/movements/aeroelastic_movement.py +++ b/pterasoftware/movements/aeroelastic_movement.py @@ -158,7 +158,7 @@ def generate_airplane_at_time_step( self, airplane_movement_index: int, step: int, - wing_deformation_angles_ixyz: list[np.ndarray] | None = None, + wing_deformation_angles_ixyz: list[np.ndarray | None] | None = None, ) -> geometry.airplane.Airplane: """Creates the Airplane at a single time step for a given AeroelasticAirplaneMovement, applying deformation from the solver's structural diff --git a/pterasoftware/output.py b/pterasoftware/output.py index d4bbe005..f65b6b88 100644 --- a/pterasoftware/output.py +++ b/pterasoftware/output.py @@ -484,7 +484,7 @@ def animate( """ if not isinstance( unsteady_solver, - (unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,), + unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver, ): raise TypeError( "unsteady_solver must be an UnsteadyRingVortexLatticeMethodSolver." diff --git a/pterasoftware/problems.py b/pterasoftware/problems.py index 79f0cbfe..998c58d0 100644 --- a/pterasoftware/problems.py +++ b/pterasoftware/problems.py @@ -13,17 +13,24 @@ from __future__ import annotations -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, cast +import matplotlib.pyplot as plt import numpy as np +from scipy.integrate import solve_ivp from . import _core, _parameter_validation, _transformations, geometry, movements from . import operating_point as operating_point_mod +from .movements import aeroelastic_movement as aeroelastic_movement_mod +from .movements import free_flight_movement as free_flight_movement_mod if TYPE_CHECKING: from ._coupled_unsteady_ring_vortex_lattice_method import ( CoupledUnsteadyRingVortexLatticeMethodSolver, ) + from .aeroelastic_unsteady_ring_vortex_lattice_method import ( + AeroelasticUnsteadyRingVortexLatticeMethodSolver, + ) class SteadyProblem: @@ -336,8 +343,845 @@ def initialize_next_problem( Must be overridden by subclasses to compute the geometry for the next time step based on the solver's results. - :param solver: The CoupledUnsteadyRingVortexLatticeMethodSolver instance - providing aerodynamic data from the current time step. + :param solver: The solver instance providing aerodynamic data from the current + time step. :return: None + :raises NotImplementedError: Always. Subclasses must override this method. """ raise NotImplementedError("Subclasses must implement initialize_next_problem.") + + +class AeroelasticUnsteadyProblem(_CoupledUnsteadyProblem): + """A subclass of _CoupledUnsteadyProblem used to couple aeroelastic wing + deformations with unsteady aerodynamics. + + This class couples aerodynamic loads with wing structural dynamics (spring-mass- + damper system) to simulate aeroelastic deformation. Each time step, wing + deformations are calculated based on the combined effects of aerodynamic moments, + inertial forces, and spring-damper restoring forces. + + **Contains the following methods:** + + calculate_wing_panel_accelerations: Computes panel accelerations from finite + difference of positions. + + calculate_mass_matrix: Generates the mass distribution matrix for wing panels. + + calculate_wing_deformation: Computes cumulative wing deformation for the current + step. + + calculate_spring_moments: Calculates spring-damper moments acting on each spanwise + section. + + calculate_torsional_spring_moment: Solves the torsional spring-damper ODE for a + single span section. + + generate_inertial_torque_function: Creates a torque function from prescribed wing + motion. + + spring_numerical_ode: Numerically integrates the spring-damper differential + equation. + + plot_flap_cycle_curves: Visualizes moment and deformation time histories. + + **Notes:** + + The aeroelastic coupling assumes a torsional spring-mass-damper model for each + spanwise section. Wing motion is prescribed through wing flapping, and aerodynamic + moments from the solver are combined with inertial and spring restoring forces via + ODE integration to produce structural deformations. + """ + + def __init__( + self, + movement: aeroelastic_movement_mod.AeroelasticMovement, + wing_density: float, + spring_constant: float, + damping_constant: float, + aero_scaling: float = 1.0, + moment_scaling_factor: float = 1.0, + damping_eps: float = 1e-3, + plot_flap_cycle: bool = False, + custom_spacing_second_derivative=None, + ) -> None: + """The initialization method. + + Sets up the aeroelastic problem with structural parameters for the torsional + spring-mass-damper model applied to each wing spanwise section. Initializes + storage for aerodynamic loads, deformations, moments, and solver state. + + See _CoupledUnsteadyProblem's initialization method for descriptions of + inherited parameters. + + :param movement: An AeroelasticMovement object containing the prescribed motion + and aerodynamic setup for the aeroelastic simulation. + :param wing_density: The mass per unit span area of the wing (kg/m^2). Used to + distribute wing mass across panels for inertial calculations. + :param spring_constant: The torsional spring stiffness for the spring-mass- + damper model (N*m/rad). Controls the restoring torque opposing deformation. + :param damping_constant: The torsional damping coefficient (N*m*s/rad). Controls + the viscous damping in the spring-mass-damper system. + :param aero_scaling: A scaling factor applied to aerodynamic moments (unitless). + The default is 1.0. Use values less than 1 to reduce aerodynamic influence. + :param moment_scaling_factor: A scaling factor applied to the computed wing + deformation angles (unitless). The default is 1.0. Useful for adjusting the + magnitude of structural response. + :param damping_eps: The critical damping tolerance used for diagnostics + (unitless). The default is 1e-3. This parameter is not currently used in the + solver. + :param plot_flap_cycle: If True, plots time histories of moments and + deformations at the end of the simulation. The default is False. + :param custom_spacing_second_derivative: An optional callable function of time + that returns the second time derivative of a custom wing motion spacing + function. Required if custom (non-sinusoidal) wing motion spacing is used. + The default is None. + :return: None + """ + if not isinstance(movement, aeroelastic_movement_mod.AeroelasticMovement): + raise TypeError("movement must be an AeroelasticMovement.") + + # Generate the initial airplane at step 0 with no deformation. + initial_airplane = movement.generate_airplane_at_time_step( + airplane_movement_index=0, step=0 + ) + + super().__init__( + movement=movement, + initial_airplanes=[initial_airplane], + initial_operating_point=movement.operating_points[0], + ) + + # Store a typed reference for aeroelastic-specific operations. + self._aeroelastic_movement: aeroelastic_movement_mod.AeroelasticMovement = ( + movement + ) + + self.plot_flap_cycle = plot_flap_cycle + self.prev_velocities: list[np.ndarray] = [] + self.positions: list[np.ndarray] = [] + self.net_deformation: np.ndarray = np.zeros((0, 3)) + self.angluar_velocities: np.ndarray = np.zeros((0, 3)) + + # Tunable Parameters + self.wing_density = wing_density # per unit height kg/m^2 + self.moment_scaling_factor = moment_scaling_factor + self.spring_constant = spring_constant + self.damping_constant = damping_constant + self.aero_scaling = aero_scaling + self.damping_eps = damping_eps # critical damping tolerance + + # Permanent parameters + self.step_discards = ( + 5 # number of initial steps to discard for numerical stability + ) + self.spacing = ( + self._aeroelastic_movement.airplane_movements[0] + .wing_movements[0] + .spacingAngles_Gs_to_Wn_ixyz[0] + ) + self.wing_movement = self._aeroelastic_movement.airplane_movements[ + 0 + ].wing_movements[0] + + self.per_step_data: list[np.ndarray] = [] + self.net_data: list[np.ndarray] = [] + self.angluar_velocity_data: list[np.ndarray] = [] + self.per_step_inertial: list[np.ndarray] = [] + self.per_step_aero: list[np.ndarray] = [] + self.per_step_spring: list[np.ndarray] = [] + self.base_wing_positions: np.ndarray = np.zeros(0) + self.flap_points: list[np.ndarray] = [] + + # For custom spacing defined in movement. + self.custom_spacing_second_derivative = custom_spacing_second_derivative + + def calculate_wing_panel_accelerations(self) -> np.ndarray: + """Compute panel accelerations using finite difference of stored positions. + + Calculates second-order accelerations using the finite difference formula: a = + (p[n] - 2*p[n-1] + p[n-2]) / dt^2. + + :return: An (N_chordwise, N_spanwise, 3) ndarray of floats representing panel + center accelerations in the global frame. Returns zeros if fewer than 3 + position snapshots are available. + """ + if len(self.positions) <= 2: + if len(self.positions) == 0: + return np.zeros(1) + return np.zeros_like(self.positions[0]) + dt = self.movement.delta_time + # If given a relatively large dt value, the finite difference calculation can produce + # very large accelerations that cause numerical instability in the spring ODE integration. + # A higher order model may be useful if this is the case. + pos_m1: np.ndarray = self.positions[-1] + pos_m2: np.ndarray = self.positions[-2] + pos_m3: np.ndarray = self.positions[-3] + return np.array((pos_m1 - 2 * pos_m2 + pos_m3) / (dt * dt)) + + def calculate_mass_matrix(self, wing: geometry.wing.Wing) -> np.ndarray: + """Generate the mass distribution matrix for all wing panels. + + Distributes the total spanwise mass (wing_density) across panel areas to form a + panel-by-panel mass matrix. Each panel's mass is proportional to its area times + the specified wing_density. + + :param wing: A Wing object whose panels define the mass distribution. + :return: An (N_chordwise, N_spanwise, 3) ndarray of floats representing the mass + at each panel. The three components are identical (mass scalar replicated + for x, y, z axes). + """ + assert wing.panels is not None + areas = np.array([[panel.area for panel in row] for row in wing.panels]) + return np.repeat(areas[:, :, None], 3, axis=2) * self.wing_density + + def initialize_next_problem( + self, solver: CoupledUnsteadyRingVortexLatticeMethodSolver + ) -> None: + # Circular at module level: the aeroelastic solver imports problems.py. + # Needed at runtime for cast(). + from .aeroelastic_unsteady_ring_vortex_lattice_method import ( + AeroelasticUnsteadyRingVortexLatticeMethodSolver, + ) + + aeroelastic_solver = cast( + AeroelasticUnsteadyRingVortexLatticeMethodSolver, solver + ) + + step = len(self._steady_problems) + deformation_matrices = self.calculate_wing_deformation(aeroelastic_solver, step) + + # Build the per-wing deformation list. The main wing (index 0) and its + # symmetric reflection (index 1) both receive the same deformation angles. + # The mirror mesh generation in the Airplane constructor ensures that + # applying the same angles produces physically symmetric deformation. + # All other wings (e.g. the v-tail) get None (no deformation applied). + num_wings = len(self._aeroelastic_movement.airplane_movements[0].wing_movements) + wing_deformation_angles_ixyz: list[np.ndarray | None] = [None] * num_wings + wing_deformation_angles_ixyz[0] = deformation_matrices + wing_deformation_angles_ixyz[1] = deformation_matrices + + # Generate the deformed airplane at this step. + airplane = self._aeroelastic_movement.generate_airplane_at_time_step( + airplane_movement_index=0, + step=step, + wing_deformation_angles_ixyz=wing_deformation_angles_ixyz, + ) + operating_point = self._aeroelastic_movement.operating_points[step] + + self._steady_problems.append( + SteadyProblem( + airplanes=[airplane], + operating_point=operating_point, + ) + ) + + def calculate_wing_deformation( + self, + solver: AeroelasticUnsteadyRingVortexLatticeMethodSolver, + step: int, + ) -> np.ndarray: + """Compute cumulative wing deformation for the current time step. + + Orchestrates the calculation of inertial moments, spring moments, and cumulative + deformation. Updates internal state and optionally generates plots. + + :param solver: The solver instance providing aerodynamic moment data + (moments_GP1_Slep). + :param step: The current time step index (0-indexed). + :return: An (N_spanwise+1, 3) ndarray of floats representing cumulative + deformation angles at each spanwise station. The y-component (index 1) + contains torsional angles in radians; x and z components are zero. + """ + curr_problem: SteadyProblem = self._steady_problems[-1] + airplane = curr_problem.airplanes[0] + wing: geometry.wing.Wing = airplane.wings[0] + + # Compute panel parameters and mass matrix once + num_chordwise_panels = wing.num_chordwise_panels + num_spanwise_panels = wing.num_spanwise_panels + assert num_spanwise_panels is not None, "num_spanwise_panels must not be None" + num_panels = num_chordwise_panels * num_spanwise_panels + mass_matrix = self.calculate_mass_matrix(wing) + + # Initialize deformation state if needed + if self.net_deformation.size == 0: + self.net_deformation = np.zeros((num_spanwise_panels + 1, 3)) + self.angluar_velocities = np.zeros((num_spanwise_panels + 1, 3)) + + # Extract aerodynamic and inertial moments + aeroMoments_GP1_Slep = self._extract_aero_moments( + solver, num_chordwise_panels, num_spanwise_panels, num_panels + ) + inertial_moments = self._calculate_inertial_moments( + solver, + wing, + mass_matrix, + num_chordwise_panels, + num_spanwise_panels, + num_panels, + ) + + # Calculate spring moments and deformation via ODE integration + thetas, omegas, spring_moments = self.calculate_spring_moments( + num_spanwise_panels=num_spanwise_panels, + wing=wing, + mass_matrix=mass_matrix, + aero_moments=aeroMoments_GP1_Slep, + step=step, + ) + + # Build deformation vector and update state + step_deformation = self._build_deformation_vector(thetas, num_spanwise_panels) + self._apply_moment_updates( + step=step, + step_deformation=step_deformation, + omegas=omegas, + inertial_moments=inertial_moments, + aeroMoments_GP1_Slep=aeroMoments_GP1_Slep, + spring_moments=spring_moments, + ) + + # Plot results at end of simulation if enabled + if self.plot_flap_cycle and step == self.num_steps - 1: + self._plot_aeroelastic_results() + + return self.net_deformation + + def _extract_aero_moments( + self, + solver: AeroelasticUnsteadyRingVortexLatticeMethodSolver, + num_chordwise_panels: int, + num_spanwise_panels: int, + num_panels: int, + ) -> np.ndarray: + """Extract and scale aerodynamic moments from the solver output. + + Uses the strip leading edge points as the reference point for moment + calculations, consistent with the assumption of a torsional spring at the + leading edge. + + :param solver: The solver instance with moments_GP1_Slep data. + :param num_chordwise_panels: Number of chordwise panel rows. + :param num_spanwise_panels: Number of spanwise panel rows. + :param num_panels: Total number of panels (num_chordwise * num_spanwise). + :return: An (N_chordwise, N_spanwise, 3) ndarray of scaled aerodynamic moments + in the global panel frame. + """ + aeroMoments_GP1_Slep = ( + np.array(solver.moments_GP1_Slep[:num_panels]).reshape( + num_chordwise_panels, num_spanwise_panels, 3 + ) + * self.aero_scaling + ) + return aeroMoments_GP1_Slep + + def _calculate_inertial_moments( + self, + solver: AeroelasticUnsteadyRingVortexLatticeMethodSolver, + wing: geometry.wing.Wing, + mass_matrix: np.ndarray, + num_chordwise_panels: int, + num_spanwise_panels: int, + num_panels: int, + ) -> np.ndarray: + """Calculate inertial moments from panel accelerations and mass distribution. + + Computes panel accelerations via finite difference, multiplies by mass to get + forces, then calculates moments about the leading edge reference point using + cross products. + + :param solver: The solver instance providing leading edge point positions. + :param wing: The Wing object containing panel definitions. + :param mass_matrix: An (N_chordwise, N_spanwise, 3) ndarray of panel masses. + :param num_chordwise_panels: Number of chordwise panel rows. + :param num_spanwise_panels: Number of spanwise panel rows. + :param num_panels: Total number of panels (num_chordwise * num_spanwise). + :return: An (N_chordwise, N_spanwise, 3) ndarray of inertial moment vectors. + """ + # Store current panel center positions + assert wing.panels is not None + self.positions.append( + np.array([[panel.Cpp_GP1_CgP1 for panel in row] for row in wing.panels]) + ) + + # Calculate panel accelerations and inertial forces + inertial_forces = self.calculate_wing_panel_accelerations() * mass_matrix + + # Calculate moments about leading edge points via cross product + inertial_moments = np.cross( + self.positions[-1] + - solver.stack_leading_edge_points[:num_panels].reshape( + (num_chordwise_panels, num_spanwise_panels, 3) + ), + inertial_forces, + axis=2, + ) + return np.array(inertial_moments) + + def _build_deformation_vector( + self, thetas: np.ndarray, num_spanwise_panels: int + ) -> np.ndarray: + """Construct the step deformation vector from torsional angles. + + Converts the torsional angles output from the spring-damper ODE (one per + spanwise section) into a full (N_spanwise+1, 3) deformation vector with scaling + applied to the y-component (torsional angle). + + :param thetas: An (N_spanwise+1,) ndarray of torsional angles in radians. + :param num_spanwise_panels: Number of spanwise panel rows. + :return: An (N_spanwise+1, 3) ndarray with zero-valued x and z components and + scaled torsional angles in the y component. + """ + step_deformation = np.array( + [ + np.array( + [ + 0, + thetas[i + 1] * self.moment_scaling_factor, + 0, + ] + ) + for i in range(num_spanwise_panels) + ] + ) + step_deformation = np.insert(step_deformation, 0, np.array([0, 0, 0]), axis=0) + return step_deformation + + def _apply_moment_updates( + self, + step: int, + step_deformation: np.ndarray, + omegas: np.ndarray, + inertial_moments: np.ndarray, + aeroMoments_GP1_Slep: np.ndarray, + spring_moments: np.ndarray, + ) -> None: + """Update internal moment and deformation state arrays. + + Stores per-step moment and deformation data, updates the cumulative net + deformation (with discarding of early unstable steps), and tracks wing + deflection points relative to the undeformed baseline. + + :param step: The current time step index. + :param step_deformation: The (N_spanwise+1, 3) deformation vector for this step. + :param omegas: An (N_spanwise+1,) ndarray of angular velocities. + :param inertial_moments: An (N_chordwise, N_spanwise, 3) ndarray of inertial + moments. + :param aeroMoments_GP1_Slep: An (N_chordwise, N_spanwise, 3) ndarray of aero + moments. + :param spring_moments: An (N_spanwise, 3) ndarray of spring-damper moments. + :return: None + """ + # Update angular velocity state + self.angluar_velocities[:, 1] = omegas + + # Generate the reference (undeformed) airplane at this step to get the + # baseline panel positions for tracking wing deflection. + ref_airplane = self._aeroelastic_movement.generate_airplane_at_time_step( + airplane_movement_index=0, step=step + ) + ref_problem = SteadyProblem( + [ref_airplane], self._aeroelastic_movement.operating_points[step] + ) + undeformed_wing = ref_problem.airplanes[0].wings[0] + assert undeformed_wing.panels is not None + undeformed_positions = np.array( + [[panel.Cpp_GP1_CgP1 for panel in row] for row in undeformed_wing.panels] + ) + if self.base_wing_positions.size == 0: + self.base_wing_positions = np.array(undeformed_positions) + + # Track wing deflection relative to undeformed baseline + self.flap_points.append( + np.array(undeformed_positions) - self.base_wing_positions + ) + + # Store per-step moment components for later analysis/plotting + self.per_step_inertial.append(inertial_moments.copy()) + self.per_step_aero.append(aeroMoments_GP1_Slep.copy()) + self.per_step_spring.append(spring_moments.copy()) + + # Update cumulative deformation (with numerical stability discarding) + # Accounts for numerical instability causing large aerodynamic forces in initial steps + if step > self.step_discards: + self.net_deformation = step_deformation + + # Store deformation and angular velocity history + self.per_step_data.append(step_deformation) + self.net_data.append(self.net_deformation.copy()) + self.angluar_velocity_data.append(self.angluar_velocities.copy()) + + def _plot_aeroelastic_results(self) -> None: + """Generate and display time-history plots of aeroelastic results. + + Creates plots of per-step and cumulative deformations, moment components + (inertial, aerodynamic, spring), and wing deflection points. Useful for + visualizing the aeroelastic coupling behavior. + + :return: None + """ + zero_curve = np.zeros((1, np.array(self.per_step_inertial).shape[0])) + + # Deformation time histories + self.plot_flap_cycle_curves( + np.array(self.per_step_data)[:, :, 1].T.tolist(), "Per Step Deformation" + ) + self.plot_flap_cycle_curves( + np.array(self.net_data)[:, :, 1].T.tolist(), "Net Deformation" + ) + + # Moment component time histories + self.plot_flap_cycle_curves( + np.vstack( + ( + zero_curve, + np.array(self.per_step_inertial)[:, :, :, 2].sum(axis=1).T, + ) + ).tolist(), + "Per Step Inertial Moments", + ) + self.plot_flap_cycle_curves( + np.vstack( + (zero_curve, np.array(self.per_step_aero)[:, :, :, 2].sum(axis=1).T) + ).tolist(), + "Per Step Aero Moments", + ) + self.plot_flap_cycle_curves( + np.vstack( + (zero_curve, np.array(self.per_step_spring)[:, :, 2].sum(axis=1).T) + ).tolist(), + "Per Step Spring Moments", + ) + + # Wing deflection tracking + self.plot_flap_cycle_curves( + np.vstack( + ( + zero_curve, + np.array(self.flap_points)[:, :, :, 2].sum(axis=1).T, + ) + ).tolist(), + "Flap Points Z", + ) + + def calculate_spring_moments( + self, + num_spanwise_panels: int, + wing: geometry.wing.Wing, + mass_matrix: np.ndarray, + aero_moments: np.ndarray, + step: int, + ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Calculate spring-damper moments and angular states for each spanwise section. + + Solves the torsional spring-damper ODE independently for each spanwise section, + accounting for aerodynamic moments, inertial forces, and structural properties. + Uses the parallel axis theorem to compute rotational inertia about the flapping + axis. + + :param num_spanwise_panels: Number of spanwise panel rows in the wing. + :param wing: The Wing object containing geometric and structural definitions. + :param mass_matrix: An (N_chordwise, N_spanwise, 3) ndarray of panel masses. + :param aero_moments: An (N_chordwise, N_spanwise, 3) ndarray of aerodynamic + moments from the aerodynamic solver. + :param step: The current time step index. + :return: A tuple of three ndarrays: - thetas: (N_spanwise+1,) ndarray of + torsional angles (radians) at each station. - omegas: (N_spanwise+1,) + ndarray of angular velocities (rad/s) at each station. - spring_moments: + (N_spanwise, 3) ndarray of spring-damper moment vectors. **Notes:** The + rotational inertia is computed as: I = (1/12)*M*(L^2 + W^2) + M*d^2, where M + is panel mass, L is chord, W is span width, and d is distance from the + flapping axis (computed cumulatively using the parallel axis theorem). + """ + spring_moments = np.zeros((num_spanwise_panels, 3)) + thetas = np.zeros(num_spanwise_panels + 1) + omegas = np.zeros(num_spanwise_panels + 1) + d = 0.0 # distance from flapping axis to panel centroid (computed in half-span increments) + for span_panel in range(num_spanwise_panels): + aero_span_moment = np.sum(aero_moments[:, span_panel, 2]) + theta0: float = 0.0 + omega0: float = 0.0 + if span_panel != 0: + theta0 = self.net_deformation[span_panel][1] + omega0 = self.angluar_velocities[span_panel][1] + + dt = self.movement.delta_time + mass = mass_matrix[:, span_panel, :].sum() + # Equation for rotational inertia of rectangular prism about flapping axis + # Considers two factors, the first is the rotational inertial of a rectangular + # prism about its centroid, the second is the parallel axis theorem to + # account for distance from flapping axis to the panel centroid + L = ( + wing.wing_cross_sections[span_panel].chord + + wing.wing_cross_sections[span_panel + 1].chord + ) / 2 + assert wing.panels is not None + W: float = float(np.linalg.norm(wing.panels[0][span_panel].frontLeg_G)) + d += W / 2 + span_I = 1 / 12 * mass * (L**2 + W**2) + mass * (d**2) + theta, omega, moment = self.calculate_torsional_spring_moment( + dt, + # A potential knob to tweak in representation of the torsional inertia + # I=mass * (wing.wing_cross_sections[span_panel].chord ** 2) / 2, + I=1 / 2 * mass * (L**2), + # I=span_I, + theta0=theta0, + omega0=omega0, + aero_span_moment=aero_span_moment, + step=step, + span_I=span_I, + ) + d += W / 2 + thetas[span_panel + 1] = theta + omegas[span_panel + 1] = omega + spring_moments[span_panel] = np.array([0, moment, 0]) + + return thetas, omegas, spring_moments + + def calculate_torsional_spring_moment( + self, + dt: float, + I: float, + theta0: float, + omega0: float, + aero_span_moment: float, + step: int, + span_I: float, + num_steps: int = 2, + ) -> tuple[float, float, float]: + """Solve the torsional spring-damper ODE for a single wing section. + + Integrates the forced torsional damped harmonic oscillator equation: I*dω/dt = + τ_aero + τ_inertial - k*θ - c*ω + + Returns the angular displacement and velocity at the end of the time step, along + with the spring-damper restoring moment. + + :param dt: The time step duration (seconds). + :param I: The rotational inertia about the flapping axis (kg*m^2). + :param theta0: Initial torsional angle at the start of the time step (radians). + :param omega0: Initial angular velocity at the start of the time step (rad/s). + :param aero_span_moment: The z-component aerodynamic moment summed over + chordwise panels for this spanwise section (N*m). + :param step: The current time step index (used for inertial torque evaluation). + :param span_I: The rotational inertia including parallel axis theorem (kg*m^2). + This is the actual inertia used in the ODE solver. + :param num_steps: Number of time sub-steps for numerical integration. The + default is 2. + :return: A tuple of (theta, omega, spring_moment) where: - theta: Final + torsional angle (radians). - omega: Final angular velocity (rad/s). - + spring_moment: The z-component spring-damper moment τ = -k*θ - c*ω (N*m). + """ + k = self.spring_constant + c = self.damping_constant + t = np.linspace(dt * (step - 1), dt * step, num_steps) + + # Forced numerical integration of the spring-damper ODE + theta, omega = self.spring_numerical_ode( + t, + k, + c, + I, + theta0, + omega0, + aero_span_moment, + self.generate_inertial_torque_function(span_I), + ) + + # Internal spring-damper moment (restoring force from structural springs/dampers) + spring_moment = -k * theta - c * omega + + return theta, omega, spring_moment + + def generate_inertial_torque_function(self, span_I: float): + """Generate the prescribed wing motion inertial torque function. + + Extracts the prescribed flapping motion from the wing_movement definition and + creates a callable inertial torque function τ_inertial = I * d²θ_prescribed/dt². + Supports sinusoidal and custom spacing functions. + + :param span_I: The rotational inertia of the wing span section about the + flapping axis (kg*m^2). + :return: A callable function that accepts time and returns the inertial torque + (N*m) due to the prescribed wing motion acceleration. **Notes:** For + sinusoidal spacing: τ = -I * b^2 * sin(b*t + h) * A, where b = 2π/period, h + = phase, A = amplitude. For custom spacing, requires + custom_spacing_second_derivative to be defined. + """ + amp = self.wing_movement.ampAngles_Gs_to_Wn_ixyz[0] + b = 2 * np.pi / self.wing_movement.periodAngles_Gs_to_Wn_ixyz[0] + h = np.deg2rad(self.wing_movement.phaseAngles_Gs_to_Wn_ixyz[0]) + if self.spacing == "sine": + torque_func = lambda time: -1 * (b**2) * np.sin(b * time + h) * amp * span_I + elif self.spacing == "uniform": + raise ValueError( + "Sawtooth function (uniform spacing) is not differentiable, " + "cannot be used for inertial torque function." + ) + elif callable(self.spacing): + if self.custom_spacing_second_derivative is not None: + torque_func = ( + lambda time: self.custom_spacing_second_derivative(time) * span_I + ) + else: + raise ValueError( + "Custom spacing function provided without second derivative function " + "for inertial torque calculation." + ) + + return torque_func + + def spring_numerical_ode( + self, + t: np.ndarray, + k: float, + c: float, + I: float, + theta0: float, + omega0: float, + aero_torque: float, + inertial_torque_func, + ) -> tuple[float, float]: + """Numerically integrate the torsional spring-damper ODE. + + Solves the second-order forced ODE: I * d²θ/dt² = τ_aero + τ_inertial(t) - k*θ - + c*dθ/dt + + using scipy.integrate.solve_ivp with strict tolerances. + + :param t: A (N,) ndarray of time points for integration evaluation. + :param k: Spring constant (N*m/rad). + :param c: Damping constant (N*m*s/rad). + :param I: Rotational inertia (kg*m^2). This parameter is present for potential + alternative models of inertia. + :param theta0: Initial angular displacement (radians). + :param omega0: Initial angular velocity (rad/s). + :param aero_torque: Constant aerodynamic torque acting on the section (N*m). + :param inertial_torque_func: A callable function of time that returns the + inertial torque from prescribed motion acceleration (N*m). + :return: A tuple of (theta, omega) representing the final angle and angular + velocity at the last time point in t. + """ + + def tau(time: float) -> float: + """Total external torque (aerodynamic + inertial from prescribed motion).""" + return float(aero_torque + inertial_torque_func(time)) + + def ode(time: float, y: np.ndarray) -> np.ndarray: + """ODE system: dθ/dt = ω, dω/dt = (τ - c*ω - k*θ)/I.""" + theta, omega = y + return np.array([omega, (tau(time) - c * omega - k * theta) / I]) + + sol = solve_ivp( + ode, + (t[0], t[-1]), + np.array([theta0, omega0]), + t_eval=t, + rtol=1e-9, + atol=1e-12, + ) + + theta = float(sol.y[0][-1]) + omega = float(sol.y[1][-1]) + + return theta, omega + + def plot_flap_cycle_curves( + self, + data: list, + title: str, + flap_cycle=None, + ) -> None: + """Visualize time histories of moments, deformations, or forces. + + Creates a multi-curve line plot showing moment or deformation values across all + time steps, with optional overlay of a reference flap cycle. + + :param data: A list of lists where each inner list represents a curve to plot. + Values in each curve are plotted against step number. + :param title: The title for the plot and the output PNG filename (spaces + replaced with underscores). + :param flap_cycle: Optional reference curve to overlay on the plot. If provided, + should be a list of values to plot with label "Flap Cycle" in black. The + default is None. + :return: None **Notes:** The plot is saved as a PNG file with the title as the + filename. The plot window is displayed to the user. Figure size is 12x6 + inches at 200 DPI. + """ + plt.figure(figsize=(12, 6), dpi=200) + + for i, curve in enumerate(data): + x = range(len(curve)) + plt.plot(x, curve, label=f"Curve {i}") + if flap_cycle is not None: + plt.plot( + range(len(flap_cycle)), flap_cycle, label=f"Flap Cycle", color="black" + ) + plt.xlabel("Step") + plt.ylabel("Value") + plt.title(title) + plt.legend() + plt.grid(True) + plt.savefig(f"{title.replace(' ', '_')}.png") + plt.show() + + +class FreeFlightUnsteadyProblem(_CoupledUnsteadyProblem): + """A subclass of _CoupledUnsteadyProblem for free flight simulations. + + This class manages the geometry for free flight simulations where the operating + point is updated by the solver at each time step based on the computed aerodynamic + forces and moments. The airplane geometry is precomputed by the FreeFlightMovement, + but the operating point evolves dynamically. + + **Contains the following methods:** + + initialize_next_problem: Initializes the next step's problem with an updated + operating point from the FreeFlightMovement. + + **Contains the following class attributes:** + + None + """ + + __slots__ = ("_free_flight_movement",) + + def __init__( + self, + movement: free_flight_movement_mod.FreeFlightMovement, + ) -> None: + """The initialization method. + + :param movement: A FreeFlightMovement object containing the precomputed airplane + geometry and mutable operating point for the free flight simulation. + :param only_final_results: If True, only calculate forces and moments for the + final motion cycle. Can be a bool or numpy bool and will be converted to + bool internally. The default is False. + :return: None + """ + if not isinstance(movement, free_flight_movement_mod.FreeFlightMovement): + raise TypeError("movement must be a FreeFlightMovement.") + + # Extract the initial airplanes (one per airplane movement, at step 0). + initial_airplanes = [airplane_tuple[0] for airplane_tuple in movement.airplanes] + initial_operating_point = movement.operating_point_movement.operating_points[0] + + super().__init__( + movement=movement, + initial_airplanes=initial_airplanes, + initial_operating_point=initial_operating_point, + ) + + self._free_flight_movement = movement + + def initialize_next_problem( + self, solver: CoupledUnsteadyRingVortexLatticeMethodSolver + ) -> None: + """Initialize the next time step's problem. + + :param solver: The solver instance providing aerodynamic data from the current + time step. + :return: None + :raises NotImplementedError: Always. Free flight solver not yet implemented. + """ + raise NotImplementedError("Free flight solver not yet implemented.") diff --git a/pterasoftware/unsteady_ring_vortex_lattice_method.py b/pterasoftware/unsteady_ring_vortex_lattice_method.py index 258c5b87..189c5e72 100644 --- a/pterasoftware/unsteady_ring_vortex_lattice_method.py +++ b/pterasoftware/unsteady_ring_vortex_lattice_method.py @@ -129,7 +129,8 @@ class UnsteadyRingVortexLatticeMethodSolver: def __init__(self, unsteady_problem: _core.CoreUnsteadyProblem) -> None: """The initialization method. - :param unsteady_problem: The UnsteadyProblem to be solved. + :param unsteady_problem: The UnsteadyProblem (or subclass of + CoreUnsteadyProblem) to be solved. :return: None """ # Guard direct instantiation of the base solver against coupled problems while @@ -658,6 +659,18 @@ def _initialize_panel_vortices(self) -> None: """Calculates the locations of the bound RingVortex vertices for all time steps, and then initializes them. + :return: None + """ + for steady_problem_id, steady_problem in enumerate(self.steady_problems): + # Find the freestream velocity (in the first Airplane's geometry axes, + # observed from the Earth frame) at this time step. + self._initialize_panel_vortex(steady_problem, steady_problem_id) + + def _initialize_panel_vortex( + self, steady_problem: problems.SteadyProblem, steady_problem_id: int + ) -> None: + """Initializes the bound RingVortex for each Panel in the given SteadyProblem. + Every Panel has a RingVortex, which is a quadrangle whose front leg is a LineVortex at the Panel's quarter chord. The left and right legs are LineVortices running along the Panel's left and right legs. If the Panel is not @@ -665,6 +678,10 @@ def _initialize_panel_vortices(self) -> None: the rear Panel's quarter chord. Otherwise, they extend backwards and meet the back LineVortex one quarter chord back from the Panel's back leg. + :param steady_problem: The SteadyProblem for which to initialize the bound + RingVortices. + :param steady_problem_id: The index of the given SteadyProblem in the list of + SteadyProblems. :return: None """ for step in range(self.num_steps): @@ -1050,7 +1067,7 @@ def calculate_solution_velocity( bound_singularity_counts: np.ndarray | None = None, wake_singularity_counts: np.ndarray | None = None, ) -> np.ndarray: - """Finds the fluid velocity (in the first Airplane's geometry axes, observed + """Finds the fluid velocity (in the fiirst Airplane's geometry axes, observed from the Earth frame) at one or more points (in the first Airplane's geometry axes, relative to the first Airplane's CG) due to the freestream velocity and the induced velocity from every RingVortex. @@ -1495,6 +1512,33 @@ def _calculate_loads(self) -> None: + unsteady_forces_GP1 ) + moments_GP1_CgP1 = self._load_calculation_moment_processing_hook( + rightLegForces_GP1, + frontLegForces_GP1, + leftLegForces_GP1, + backLegForces_GP1, + unsteady_forces_GP1, + ) + + # TODO: Transform forces_GP1 and moments_GP1_CgP1 to each Airplane's local + # geometry axes before passing to process_solver_loads. + _functions.process_solver_loads(self, forces_GP1, moments_GP1_CgP1) + + def _load_calculation_moment_processing_hook( + self, + rightLegForces_GP1, + frontLegForces_GP1, + leftLegForces_GP1, + backLegForces_GP1, + unsteady_forces_GP1, + ) -> np.ndarray: + """A hook method for processing the moments calculated in _calculate_loads. This + is added to allow for overriding the moment calculation in a child class. + + :return: moments_GP1_CgP1, a (N,3) ndarray of floats representing the moments + (in the first Airplane's geometry axes, relative to the first Airplane's CG) + on every Panel at the current time step. + """ # Find the moments (in the first Airplane's geometry axes, relative to the # first Airplane's CG) on the Panels' RingVortex's right LineVortex, # front LineVortex, left LineVortex, and back LineVortex. @@ -1529,9 +1573,7 @@ def _calculate_loads(self) -> None: + unsteady_moments_GP1_CgP1 ) - # TODO: Transform forces_GP1 and moments_GP1_CgP1 to each Airplane's local - # geometry axes before passing to process_solver_loads. - _functions.process_solver_loads(self, forces_GP1, moments_GP1_CgP1) + return np.array(moments_GP1_CgP1) def _populate_next_airplanes_wake(self) -> None: """Updates the next time step's Airplanes' wakes.