Skip to content

Commit bf4e779

Browse files
kanekoshkanekosh
authored
Added a stability derivative example (#462)
* added stability deriv example * minor
1 parent f48575f commit bf4e779

3 files changed

Lines changed: 200 additions & 0 deletions

File tree

‎openaerostruct/docs/advanced_features.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,3 +21,4 @@ These examples show some advanced features in OpenAeroStruct.
2121
advanced_features/openconcept_coupling.rst
2222
advanced_features/customizing_prob_setup.rst
2323
advanced_features/mphys_coupling.rst
24+
advanced_features/stability_derivatives.rst

‎openaerostruct/docs/advanced_features/stability_derivatives.rst‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
.. _Stability_derivatives:
2+
3+
Stability Derivatives
4+
=====================
5+
6+
OpenAeroStruct can compute stability derivatives and use them for design optimization, but the implementation is a bit tricky.
7+
This example shows how to compute the longitudinal stability derivatives and static margins for an aerodynamic lifting surface, but this example can be extended to aerostructural analysis and design optimization.
8+
9+
A challenge associated with stability derivatives in OpenAeroStruct is that for gradient-based optimization, we need derivatives of stability derivatives w.r.t. design variables (e.g., derivatives of CL_alpha w.r.t. wing sweep).
10+
The derivatives of stability derivatives are second-order derivatives of OpenAeroStruct's outputs (e.g., CL), which cannot be computed analytically in OpenAeroStruct.
11+
In fact, this is an OpenMDAO's limitation because the background theory of OpenMDAO is only applicable to first-order derivatives.
12+
13+
To overcome this limitation, we compute the stability derivatives using finite difference, and then analytically differentiate the finite difference component using OpenMDAO.
14+
This can be implemented in multiple ways, but in this example, we use an approach similar to multipoint optimization, where we create two `AeroPoint` instances---one for the specified flight condition, and the other for the perturbed angle of attack for finite difference derivatives.
15+
The script is similar to the :ref:`Multipoint Optimization` example, but in this example, the flight conditions for two AeroPoints are not independent, and we add a few `ExecComps` to compute perturbed angle of attack, stability derivatives, and static margin.
16+
17+
The following script computes CL_alpha, CM_alpha, and static margin for a swept wing.
18+
You can play around with different sweep angles and see how the static margin changes as a function of the sweep angle.
19+
Note that this example does not solve the trim (CM=0).
20+
21+
22+
.. embed-code::
23+
../examples/stability_derivatives.py

‎openaerostruct/examples/stability_derivatives.py‎

Lines changed: 176 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,176 @@
1+
"""
2+
Example of aerodynamic analysis including stability derivatives (CL_alpha and CM_alpha).
3+
4+
We compute CL_alpha and CM_alpha by finite differencing with respect to alpha.
5+
To do so, we instantiate AeroPoint at alpha and alpha + delta_alpha, and add an ExecComp to compute finite difference.
6+
Finite differencing w.r.t. alpha allows us to compute the derivatives of CL_alpha and CM_alpha w.r.t. design variables.
7+
8+
Note that this example does not trim the aircraft (e.g. CM != 0).
9+
"""
10+
11+
import numpy as np
12+
import openmdao.api as om
13+
14+
from openaerostruct.meshing.mesh_generator import generate_mesh
15+
from openaerostruct.geometry.geometry_group import Geometry
16+
from openaerostruct.aerodynamics.aero_groups import AeroPoint
17+
18+
# Create a dictionary to store options about the surface.
19+
# Here, we setup a simple rectangular mesh with chord = 1 m and span = 10 m.
20+
# We then apply sweep after prob.setup()
21+
mesh_dict = {
22+
"num_y": 15,
23+
"num_x": 3,
24+
"wing_type": "rect",
25+
"root_chord": 1.0, # m
26+
"span": 10.0, # m
27+
"symmetry": True,
28+
"num_twist_cp": 2,
29+
"span_cos_spacing": 0.0,
30+
}
31+
32+
mesh = generate_mesh(mesh_dict)
33+
34+
surf_dict = {
35+
# Wing definition
36+
"name": "wing", # name of the surface
37+
"symmetry": True, # if true, model one half of wing
38+
# reflected across the plane y = 0
39+
"S_ref_type": "wetted", # how we compute the wing area,
40+
# can be 'wetted' or 'projected'
41+
"mesh": mesh,
42+
"twist_cp": np.zeros(mesh_dict["num_twist_cp"]), # twist control points
43+
"sweep": 0.0, # wing sweep angle
44+
# Aerodynamic performance of the lifting surface at
45+
# an angle of attack of 0 (alpha=0).
46+
# These CL0 and CD0 values are added to the CL and CD
47+
# obtained from aerodynamic analysis of the surface to get
48+
# the total CL and CD.
49+
# These CL0 and CD0 values do not vary wrt alpha.
50+
"CL0": 0.0, # CL of the surface at alpha=0
51+
"CD0": 0.015, # CD of the surface at alpha=0
52+
# Airfoil properties for viscous drag calculation
53+
"k_lam": 0.05, # percentage of chord with laminar
54+
# flow, used for viscous drag
55+
"t_over_c_cp": np.array([0.15]), # thickness over chord ratio (NACA0015)
56+
"c_max_t": 0.303, # chordwise location of maximum (NACA0015)
57+
# thickness
58+
"with_viscous": True, # if true, compute viscous drag
59+
"with_wave": False, # if true, compute wave drag
60+
}
61+
62+
surfaces = [surf_dict]
63+
64+
# Although this example only considers cruise, hence single point,
65+
# We still need to AeroPoint instances to finite difference CL and CM w.r.t. alpha
66+
n_points = 2
67+
68+
# Create the problem and the model group
69+
prob = om.Problem()
70+
71+
indep_var_comp = om.IndepVarComp()
72+
indep_var_comp.add_output("v", val=248.136, units="m/s")
73+
indep_var_comp.add_output("alpha", val=5.0, units="deg")
74+
indep_var_comp.add_output("Mach_number", val=0.84)
75+
indep_var_comp.add_output("re", val=1.0e6, units="1/m")
76+
indep_var_comp.add_output("rho", val=0.38, units="kg/m**3")
77+
# Set center of gravity to be 0.5 m from the leading edge.
78+
# Note that the CG location does *not* move as a result of wing shape changes
79+
indep_var_comp.add_output("cg", val=np.array([0.5, 0.0, 0.0]), units="m")
80+
81+
prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"])
82+
83+
# Compute alpha perturbation for finite difference
84+
alpha_FD_stepsize = 1e-4 # deg
85+
alpha_perturb_comp = om.ExecComp(
86+
"alpha_plus_delta = alpha + delta_alpha",
87+
units="deg",
88+
delta_alpha={"val": alpha_FD_stepsize, "constant": True},
89+
)
90+
prob.model.add_subsystem("alpha_for_FD", alpha_perturb_comp, promotes=["*"])
91+
92+
# Loop over each surface and create the geometry groups
93+
for surface in surfaces:
94+
# Get the surface name and create a group to contain components only for this surface.
95+
name = surface["name"]
96+
geom_group = Geometry(surface=surface)
97+
98+
# Add geom_group to the problem with the name of the surface.
99+
prob.model.add_subsystem(name + "_geom", geom_group)
100+
101+
# Create aero analysis point for alpha and alpha_plus_delta
102+
point_names = ["aero_point", "aero_point_FD"]
103+
for i in range(n_points):
104+
# Create the aero point group and add it to the model
105+
aero_group = AeroPoint(surfaces=surfaces)
106+
point_name = point_names[i]
107+
prob.model.add_subsystem(point_name, aero_group)
108+
109+
# Connect flow properties to the analysis point
110+
prob.model.connect("v", point_name + ".v")
111+
prob.model.connect("Mach_number", point_name + ".Mach_number")
112+
prob.model.connect("re", point_name + ".re")
113+
prob.model.connect("rho", point_name + ".rho")
114+
prob.model.connect("cg", point_name + ".cg")
115+
116+
# Connect angle of attack. Use perturbed alpha for the second point for finite difference
117+
alpha_name = "alpha" if i == 0 else "alpha_plus_delta"
118+
prob.model.connect(alpha_name, point_name + ".alpha")
119+
120+
# Connect the parameters within the model for each aero point
121+
for surface in surfaces:
122+
name = surface["name"]
123+
124+
# Connect the mesh from the geometry component to the analysis point
125+
prob.model.connect(name + "_geom.mesh", point_name + "." + name + ".def_mesh")
126+
127+
# Perform the connections with the modified names within the 'aero_states' group.
128+
prob.model.connect(name + "_geom.mesh", point_name + ".aero_states." + name + "_def_mesh")
129+
prob.model.connect(name + "_geom.t_over_c", point_name + "." + name + "_perf." + "t_over_c")
130+
131+
# Compute stability derivatives by finite difference
132+
stabibility_derivatives_comp = om.ExecComp(
133+
["CL_alpha = (CL_FD - CL) / delta_alpha", "CM_alpha = (CM_FD - CM) / delta_alpha"],
134+
delta_alpha={"val": alpha_FD_stepsize, "constant": True},
135+
CL_alpha={"val": 0.0, "units": "1/deg"},
136+
CL_FD={"val": 0.0, "units": None},
137+
CL={"val": 0.0, "units": None},
138+
CM_alpha={"val": np.zeros(3), "units": "1/deg"},
139+
CM_FD={"val": np.zeros(3), "units": None},
140+
CM={"val": np.zeros(3), "units": None},
141+
)
142+
prob.model.add_subsystem("stability_derivs", stabibility_derivatives_comp, promotes_outputs=["*"])
143+
# Connect CL and CM from aero points
144+
prob.model.connect("aero_point.CL", "stability_derivs.CL")
145+
prob.model.connect("aero_point.CM", "stability_derivs.CM")
146+
prob.model.connect("aero_point_FD.CL", "stability_derivs.CL_FD")
147+
prob.model.connect("aero_point_FD.CM", "stability_derivs.CM_FD")
148+
149+
# Compute static margin
150+
static_margin_comp = om.ExecComp(
151+
"static_margin = -CM_alpha / CL_alpha",
152+
CM_alpha={"val": 0.0, "units": "1/deg"}, # for pitching moment
153+
CL_alpha={"val": 0.0, "units": "1/deg"},
154+
static_margin={"val": 0.0, "units": None}, # static margin
155+
)
156+
prob.model.add_subsystem("static_margin", static_margin_comp, promotes_outputs=["*"])
157+
# Connect stability derivatives to static margin component
158+
prob.model.connect("CL_alpha", "static_margin.CL_alpha")
159+
prob.model.connect("CM_alpha", "static_margin.CM_alpha", src_indices=1) # connect pitching moment deriv
160+
161+
# Set up the problem
162+
prob.setup()
163+
164+
# Set sweep angle
165+
prob.set_val("wing_geom.sweep", 10.0, units="deg")
166+
167+
# Run model and print outputs
168+
prob.run_model()
169+
170+
print("Sweep angle =", prob.get_val("wing_geom.sweep", units="deg"), "deg")
171+
print("CL =", prob.get_val("aero_point.CL"))
172+
print("CD =", prob.get_val("aero_point.CD"))
173+
print("CM =", prob.get_val("aero_point.CM"))
174+
print("CL_alpha =", prob.get_val("CL_alpha", units="1/deg"), "1/deg")
175+
print("CM_alpha =", prob.get_val("CM_alpha", units="1/deg"), "1/deg")
176+
print("Static Margin =", prob.get_val("static_margin"))

0 commit comments

Comments
 (0)