Fixed viscous drag, moment coeff, and sectional CL bugs (#397)

* fixed viscous drag

* updated test ref values

* version

* black

* renamed variables and fixed moment coeff

* ffd ref value

* CM ref value updated in vsp test

* Updated geometry.py docstrings

Author: kanekosh

openaerostruct/__init__.py

@@ -1 +1 @@
-__version__ = "2.5.1"
+__version__ = "2.5.2"

openaerostruct/aerodynamics/aero_groups.py

@@ -44,7 +44,7 @@ def setup(self):
             self.connect(name + ".widths", name + "_perf.widths")
             self.connect(name + ".chords", name + "_perf.chords")
             self.connect(name + ".lengths", name + "_perf.lengths")
-            self.connect(name + ".cos_sweep", name + "_perf.cos_sweep")
+            self.connect(name + ".lengths_spanwise", name + "_perf.lengths_spanwise")
 
             # Connect S_ref for performance calcs
             self.connect(name + ".S_ref", "total_perf." + name + "_S_ref")

openaerostruct/aerodynamics/functionals.py

@@ -45,14 +45,14 @@ def setup(self):
         self.add_subsystem(
             "viscousdrag",
             ViscousDrag(surface=surface),
-            promotes_inputs=["Mach_number", "re", "widths", "cos_sweep", "lengths", "S_ref", "t_over_c"],
+            promotes_inputs=["Mach_number", "re", "widths", "lengths_spanwise", "lengths", "S_ref", "t_over_c"],
             promotes_outputs=["CDv"],
         )
 
         self.add_subsystem(
             "wavedrag",
             WaveDrag(surface=surface),
-            promotes_inputs=["Mach_number", "cos_sweep", "widths", "CL", "chords", "t_over_c"],
+            promotes_inputs=["Mach_number", "lengths_spanwise", "widths", "CL", "chords", "t_over_c"],
             promotes_outputs=["CDw"],
         )
 

openaerostruct/aerodynamics/geometry.py

@@ -21,7 +21,11 @@ class VLMGeometry(om.ExplicitComponent):
     b_pts[nx-1, ny, 3] : numpy array
         Bound points for the horseshoe vortices, found along the 1/4 chord.
     widths[ny-1] : numpy array
-        The spanwise widths of each individual panel.
+        The widths of each individual panel along y-axis (spanwise direction with zero sweep).
+    lengths_spanwise[ny-1] : numpy array
+        The the length of the quarter-chord line of each panels.
+        This is identical to `widths` if sweep angle is 0.
+        When wing is swept, `lengths_spanwise` is longer than `widths`.
     lengths[ny] : numpy array
         The chordwise length of the entire airfoil following the camber line.
     chords[ny] : numpy array
@@ -49,7 +53,7 @@ def setup(self):
         rng = np.random.default_rng(314)
         self.add_output("b_pts", val=rng.random((nx - 1, ny, 3)), units="m")
         self.add_output("widths", val=np.ones((ny - 1)), units="m")
-        self.add_output("cos_sweep", val=np.zeros((ny - 1)), units="m")
+        self.add_output("lengths_spanwise", val=np.ones((ny - 1)), units="m")
         self.add_output("lengths", val=np.zeros((ny)), units="m")
         self.add_output("chords", val=np.zeros((ny)), units="m")
         self.add_output("normals", val=np.zeros((nx - 1, ny - 1, 3)))
@@ -68,19 +72,19 @@ def setup(self):
         val[size:] = 0.25
         self.declare_partials("b_pts", "def_mesh", rows=rows, cols=cols, val=val)
 
-        # widths
+        # width
         size = ny - 1
         base = np.arange(size)
-        rows = np.tile(base, 12)
-        col = np.tile(3 * base, 6) + np.repeat(np.arange(6), len(base))
+        rows = np.tile(base, 8)
+        col = np.tile(3 * base, 4) + np.repeat([1, 2, 4, 5], len(base))
         cols = np.tile(col, 2) + np.repeat([0, (nx - 1) * ny * 3], len(col))
         self.declare_partials("widths", "def_mesh", rows=rows, cols=cols)
 
-        # cos_sweep
-        rows = np.tile(base, 8)
-        col = np.tile(3 * base, 4) + np.repeat([1, 2, 4, 5], len(base))
+        # length of panel in spanwise direction with sweep
+        rows = np.tile(base, 12)
+        col = np.tile(3 * base, 6) + np.repeat(np.arange(6), len(base))
         cols = np.tile(col, 2) + np.repeat([0, (nx - 1) * ny * 3], len(col))
-        self.declare_partials("cos_sweep", "def_mesh", rows=rows, cols=cols)
+        self.declare_partials("lengths_spanwise", "def_mesh", rows=rows, cols=cols)
 
         # lengths
         size = ny
@@ -116,16 +120,13 @@ def compute(self, inputs, outputs):
         # Compute the bound points at quarter-chord
         b_pts = mesh[:-1, :, :] * 0.75 + mesh[1:, :, :] * 0.25
 
-        # Compute the widths of each panel at the quarter-chord line
+        # Compute the length of the quarter-chord line of each panels
         quarter_chord = 0.25 * mesh[-1] + 0.75 * mesh[0]
-        widths = np.linalg.norm(quarter_chord[1:, :] - quarter_chord[:-1, :], axis=1)
+        lengths_spanwise = np.linalg.norm(quarter_chord[1:, :] - quarter_chord[:-1, :], axis=1)
 
-        # Compute the numerator of the cosine of the sweep angle of each panel
-        # (we need this for the viscous drag dependence on sweep, and we only compute
-        # the numerator because the denominator of the cosine fraction is the width,
-        # which we have already computed. They are combined in the viscous drag
-        # calculation.)
-        cos_sweep = np.linalg.norm(quarter_chord[1:, [1, 2]] - quarter_chord[:-1, [1, 2]], axis=1)
+        # Compute the widths of each panel.
+        # This is a projection of `lengths_spanwise` to the spanwise axis (y-axis)
+        widths = np.linalg.norm(quarter_chord[1:, [1, 2]] - quarter_chord[:-1, [1, 2]], axis=1)
 
         # Compute the length of each chordwise set of mesh points through the camber line.
         dx = mesh[1:, :, 0] - mesh[:-1, :, 0]
@@ -171,7 +172,7 @@ def compute(self, inputs, outputs):
         # Store each array in the outputs dict
         outputs["b_pts"] = b_pts
         outputs["widths"] = widths
-        outputs["cos_sweep"] = cos_sweep
+        outputs["lengths_spanwise"] = lengths_spanwise
         outputs["lengths"] = lengths
         outputs["normals"] = normals
         outputs["S_ref"] = S_ref
@@ -184,18 +185,18 @@ def compute_partials(self, inputs, partials):
         ny = self.ny
         mesh = inputs["def_mesh"]
 
-        # Compute the widths of each panel at the quarter-chord line
+        # Compute the length of the quarter-chord line of each panels
         quarter_chord = 0.25 * mesh[-1] + 0.75 * mesh[0]
-        widths = np.linalg.norm(quarter_chord[1:, :] - quarter_chord[:-1, :], axis=1)
+        lengths_spanwise = np.linalg.norm(quarter_chord[1:, :] - quarter_chord[:-1, :], axis=1)
 
-        # Compute the cosine of the sweep angle of each panel
-        cos_sweep_array = np.linalg.norm(quarter_chord[1:, [1, 2]] - quarter_chord[:-1, [1, 2]], axis=1)
+        # Compute the widths of each panel.
+        widths = np.linalg.norm(quarter_chord[1:, [1, 2]] - quarter_chord[:-1, [1, 2]], axis=1)
 
         delta = np.diff(quarter_chord, axis=0).T
-        d1 = delta / widths
+        d1 = delta / lengths_spanwise
+        partials["lengths_spanwise", "def_mesh"] = np.outer([-0.75, 0.75, -0.25, 0.25], d1.flatten()).flatten()
+        d1 = delta[1:, :] / widths
         partials["widths", "def_mesh"] = np.outer([-0.75, 0.75, -0.25, 0.25], d1.flatten()).flatten()
-        d1 = delta[1:, :] / cos_sweep_array
-        partials["cos_sweep", "def_mesh"] = np.outer([-0.75, 0.75, -0.25, 0.25], d1.flatten()).flatten()
 
         partials["lengths", "def_mesh"][:] = 0.0
         dmesh = np.diff(mesh, axis=0)

openaerostruct/aerodynamics/lift_coeff_2D.py

@@ -45,7 +45,7 @@ def setup(self):
         # Inputs
         self.add_input("alpha", val=3.0, units="deg")
         self.add_input("sec_forces", val=np.ones((self.nx - 1, self.ny - 1, 3)), units="N")
-        self.add_input("widths", val=np.arange((self.ny - 1)) + 1.0, units="m")
+        self.add_input("widths", val=np.ones((self.ny - 1)) * 0.2, units="m")
         self.add_input("chords", val=np.ones((self.ny)), units="m")
         self.add_input("v", val=1.0, units="m/s")
         self.add_input("rho", val=1.0, units="kg/m**3")

openaerostruct/aerodynamics/tests/test_geometry.py

@@ -191,8 +191,8 @@ def test_outputs(self):
         prob.setup()
         prob.run_model()
 
-        assert_near_equal(prob["widths"], np.array([11.95624787, 11.90425878, 11.44086572]), 1e-6)
-        assert_near_equal(prob["cos_sweep"], np.array([9.7938336, 9.79384207, 9.79385053]), 1e-6)
+        assert_near_equal(prob["lengths_spanwise"], np.array([11.95624787, 11.90425878, 11.44086572]), 1e-6)
+        assert_near_equal(prob["widths"], np.array([9.7938336, 9.79384207, 9.79385053]), 1e-6)
         assert_near_equal(prob["S_ref"], np.array([415.02211208]), 1e-6)
         assert_near_equal(prob["chords"], np.array([2.72796, 5.1252628, 7.8891638, 13.6189974]), 1e-6)
         assert_near_equal(prob["lengths"], np.array([2.72796, 5.1252628, 7.8891638, 13.6189974]), 1e-6)

openaerostruct/aerodynamics/tests/test_wave_drag.py

@@ -20,8 +20,8 @@ def test(self):
         indep_var_comp.add_output("t_over_c_cp", val=surface["t_over_c_cp"])
         indep_var_comp.add_output("Mach_number", val=0.95)
         indep_var_comp.add_output("CL", val=0.7)
-        indep_var_comp.add_output("widths", val=np.array([12.14757848, 11.91832712, 11.43730892]), units="m")
-        indep_var_comp.add_output("cos_sweep", val=np.array([10.01555924, 9.80832351, 9.79003729]), units="m")
+        indep_var_comp.add_output("lengths_spanwise", val=np.array([12.14757848, 11.91832712, 11.43730892]), units="m")
+        indep_var_comp.add_output("widths", val=np.array([10.01555924, 9.80832351, 9.79003729]), units="m")
         indep_var_comp.add_output("chords", val=np.array([2.72835132, 5.12528179, 7.88916016, 13.6189974]), units="m")
         group.add_subsystem("indep_var_comp", indep_var_comp, promotes=["*"])
 

openaerostruct/aerodynamics/viscous_drag.py

@@ -55,8 +55,8 @@ def setup(self):
         self.add_input("re", val=5.0e6, units="1/m")
         self.add_input("Mach_number", val=1.6)
         self.add_input("S_ref", val=1.0, units="m**2")
-        self.add_input("cos_sweep", val=np.ones((ny - 1)) * 0.2, units="m")
-        self.add_input("widths", val=np.arange((ny - 1)) + 1.0, units="m")
+        self.add_input("widths", val=np.ones((ny - 1)) * 0.2, units="m")
+        self.add_input("lengths_spanwise", val=np.arange((ny - 1)) + 1.0, units="m")
         self.add_input("lengths", val=np.ones((ny)), units="m")
         self.add_input("t_over_c", val=np.arange((ny - 1)))
         self.add_output("CDv", val=0.0)
@@ -72,7 +72,7 @@ def compute(self, inputs, outputs):
             S_ref = inputs["S_ref"]
             widths = inputs["widths"]
             lengths = inputs["lengths"]
-            cos_sweep = inputs["cos_sweep"] / widths
+            cos_sweep = inputs["widths"] / inputs["lengths_spanwise"]
             t_over_c = inputs["t_over_c"]
 
             # Take panel chord length to be average of its edge lengths
@@ -124,10 +124,9 @@ def compute_partials(self, inputs, partials):
 
             M = inputs["Mach_number"]
             S_ref = inputs["S_ref"]
-
             widths = inputs["widths"]
             lengths = inputs["lengths"]
-            cos_sweep = inputs["cos_sweep"] / widths
+            cos_sweep = widths / inputs["lengths_spanwise"]
 
             # Take panel chord length to be average of its edge lengths
             chords = (lengths[1:] + lengths[:-1]) / 2.0
@@ -186,9 +185,9 @@ def compute_partials(self, inputs, partials):
 
             partials["CDv", "lengths"][0, 1:] += CDv_lengths
             partials["CDv", "lengths"][0, :-1] += CDv_lengths
-            partials["CDv", "widths"][0, :] = d_over_q * FF / S_ref * 0.72
+            partials["CDv", "lengths_spanwise"][0, :] = d_over_q / S_ref * (-0.28 * k_FF * cos_sweep**1.28)
             partials["CDv", "S_ref"] = -D_over_q / S_ref**2
-            partials["CDv", "cos_sweep"][0, :] = 0.28 * k_FF * d_over_q / S_ref / cos_sweep**0.72
+            partials["CDv", "widths"][0, :] = d_over_q / S_ref * (FF + 0.28 * k_FF * cos_sweep**0.28)
             partials["CDv", "t_over_c"] = (
                 d_over_q * widths * 1.34 * M**0.18 * (0.6 / self.c_max_t + 400 * t_over_c**3) * cos_sweep**0.28
             ) / S_ref
@@ -239,9 +238,9 @@ def compute_partials(self, inputs, partials):
 
             if self.surface["symmetry"]:
                 partials["CDv", "lengths"][0, :] *= 2
-                partials["CDv", "widths"][0, :] *= 2
+                partials["CDv", "lengths_spanwise"][0, :] *= 2
                 partials["CDv", "S_ref"] *= 2
-                partials["CDv", "cos_sweep"][0, :] *= 2
+                partials["CDv", "widths"][0, :] *= 2
                 partials["CDv", "Mach_number"][0, :] *= 2
                 partials["CDv", "re"][0, :] *= 2
                 partials["CDv", "t_over_c"][0, :] *= 2

openaerostruct/aerodynamics/wave_drag.py

@@ -12,10 +12,10 @@ class WaveDrag(om.ExplicitComponent):
     ----------
     Mach_number : float
         Mach number.
-    cos_sweep[ny-1] : numpy array
-        The width in the spanwise direction of each VLM panel. This is the numerator of cos(sweep).
     widths[ny-1] : numpy array
-        The actual width of each VLM panel, rotated by the sweep angle. This is the denominator
+        The width in the spanwise direction of each VLM panel. This is the numerator of cos(sweep).
+    lengths_spanwise[ny-1] : numpy array
+        The spanwise length of each VLM panel at 1/4 chord, rotated by the sweep angle. This is the denominator
         of cos(sweep)
     CL : float
         The CL of the lifting surface used for wave drag estimation.
@@ -46,9 +46,9 @@ def setup(self):
         ny = surface["mesh"].shape[1]
 
         self.add_input("Mach_number", val=1.6)
-        self.add_input("cos_sweep", val=np.ones((ny - 1)) * 0.2, units="m")
+        self.add_input("widths", val=np.ones((ny - 1)) * 0.2, units="m")
         self.add_input(
-            "widths", val=np.arange((ny - 1)) + 1.0, units="m"
+            "lengths_spanwise", val=np.arange((ny - 1)) + 1.0, units="m"
         )  # set to np.arange so that d_CDw_d_chords is nonzero
         self.add_input("CL", val=0.33)
         self.add_input("chords", val=np.ones((ny)), units="m")
@@ -62,16 +62,14 @@ def compute(self, inputs, outputs):
         if self.with_wave:
             t_over_c = inputs["t_over_c"]
             widths = inputs["widths"]
-            actual_cos_sweep = inputs["cos_sweep"] / widths
+            cos_sweep = widths / inputs["lengths_spanwise"]
             M = inputs["Mach_number"]
             chords = inputs["chords"]
             CL = inputs["CL"]
 
             mean_chords = (chords[:-1] + chords[1:]) / 2.0
-            panel_areas = mean_chords * inputs["cos_sweep"]
-            avg_cos_sweep = np.sum(actual_cos_sweep * panel_areas) / np.sum(
-                panel_areas
-            )  # weighted average of 1/4 chord sweep
+            panel_areas = mean_chords * widths
+            avg_cos_sweep = np.sum(cos_sweep * panel_areas) / np.sum(panel_areas)  # weighted average of 1/4 chord sweep
             avg_t_over_c = np.sum(t_over_c * panel_areas) / np.sum(panel_areas)  # weighted average of streamwise t/c
             MDD = self.ka / avg_cos_sweep - avg_t_over_c / avg_cos_sweep**2 - CL / (10 * avg_cos_sweep**3)
             Mcrit = MDD - (0.1 / 80.0) ** (1.0 / 3.0)
@@ -92,16 +90,16 @@ def compute_partials(self, inputs, partials):
             ny = self.surface["mesh"].shape[1]
             t_over_c = inputs["t_over_c"]
             widths = inputs["widths"]
-            cos_sweep = inputs["cos_sweep"]
-            actual_cos_sweep = cos_sweep / widths
+            lengths_spanwise = inputs["lengths_spanwise"]
+            cos_sweep = widths / lengths_spanwise
             M = inputs["Mach_number"]
             chords = inputs["chords"]
             CL = inputs["CL"]
 
             chords = (chords[:-1] + chords[1:]) / 2.0
-            panel_areas = chords * inputs["cos_sweep"]
+            panel_areas = chords * widths
             sum_panel_areas = np.sum(panel_areas)
-            avg_cos_sweep = np.sum(actual_cos_sweep * panel_areas) / sum_panel_areas
+            avg_cos_sweep = np.sum(cos_sweep * panel_areas) / sum_panel_areas
             avg_t_over_c = np.sum(t_over_c * panel_areas) / sum_panel_areas
 
             MDD = 0.95 / avg_cos_sweep - avg_t_over_c / avg_cos_sweep**2 - CL / (10 * avg_cos_sweep**3)
@@ -116,14 +114,14 @@ def compute_partials(self, inputs, partials):
                 dMDDdtoc = -1.0 / (avg_cos_sweep**2)
                 dtocavgdtoc = panel_areas / sum_panel_areas
 
-                ccos = np.sum(cos_sweep * chords)
-                ccos2w = np.sum(chords * cos_sweep**2 / widths)
+                ccos = np.sum(widths * chords)
+                ccos2w = np.sum(chords * widths**2 / lengths_spanwise)
 
-                davgdcos = 2 * chords * cos_sweep / widths / ccos - chords * ccos2w / ccos**2
-                dtocdcos = chords * t_over_c / ccos - chords * np.sum(chords * cos_sweep * t_over_c) / ccos**2
-                davgdw = -1 * chords * cos_sweep**2 / widths**2 / ccos
-                davgdc = cos_sweep**2 / widths / ccos - cos_sweep * ccos2w / ccos**2
-                dtocdc = t_over_c * cos_sweep / ccos - cos_sweep * np.sum(chords * cos_sweep * t_over_c) / ccos**2
+                davgdcos = 2 * chords * widths / lengths_spanwise / ccos - chords * ccos2w / ccos**2
+                dtocdcos = chords * t_over_c / ccos - chords * np.sum(chords * widths * t_over_c) / ccos**2
+                davgdw = -1 * chords * widths**2 / lengths_spanwise**2 / ccos
+                davgdc = widths**2 / lengths_spanwise / ccos - widths * ccos2w / ccos**2
+                dtocdc = t_over_c * widths / ccos - widths * np.sum(chords * widths * t_over_c) / ccos**2
 
                 dcdchords = np.zeros((ny - 1, ny))
                 i, j = np.indices(dcdchords.shape)
@@ -132,17 +130,17 @@ def compute_partials(self, inputs, partials):
 
                 partials["CDw", "Mach_number"] = -1 * dCDwdMDD
                 partials["CDw", "CL"] = dCDwdMDD * dMDDdCL
-                partials["CDw", "widths"] = dCDwdMDD * dMDDdavg * davgdw
-                partials["CDw", "cos_sweep"] = dCDwdMDD * dMDDdavg * davgdcos + dCDwdMDD * dMDDdtoc * dtocdcos
+                partials["CDw", "lengths_spanwise"] = dCDwdMDD * dMDDdavg * davgdw
+                partials["CDw", "widths"] = dCDwdMDD * dMDDdavg * davgdcos + dCDwdMDD * dMDDdtoc * dtocdcos
                 partials["CDw", "chords"] = dCDwdMDD * dMDDdavg * np.matmul(
                     davgdc, dcdchords
                 ) + dCDwdMDD * dMDDdtoc * np.matmul(dtocdc, dcdchords)
                 partials["CDw", "t_over_c"] = dCDwdMDD * dMDDdtoc * dtocavgdtoc
 
         if self.surface["symmetry"]:
             partials["CDw", "CL"][0, :] *= 2
+            partials["CDw", "lengths_spanwise"][0, :] *= 2
             partials["CDw", "widths"][0, :] *= 2
-            partials["CDw", "cos_sweep"][0, :] *= 2
             partials["CDw", "Mach_number"][0, :] *= 2
             partials["CDw", "chords"][0, :] *= 2
             partials["CDw", "t_over_c"][0, :] *= 2

openaerostruct/integration/aerostruct_groups.py

@@ -134,7 +134,7 @@ def setup(self):
             "aero_geom",
             VLMGeometry(surface=surface),
             promotes_inputs=["def_mesh"],
-            promotes_outputs=["b_pts", "widths", "cos_sweep", "lengths", "chords", "normals", "S_ref"],
+            promotes_outputs=["b_pts", "widths", "lengths_spanwise", "lengths", "chords", "normals", "S_ref"],
         )
 
         self.linear_solver = om.LinearRunOnce()
@@ -158,7 +158,7 @@ def setup(self):
                 "re",
                 "rho",
                 "widths",
-                "cos_sweep",
+                "lengths_spanwise",
                 "lengths",
                 "S_ref",
                 "sec_forces",
@@ -252,7 +252,7 @@ def setup(self):
             self.connect("coupled." + name + ".widths", name + "_perf.widths")
             # self.connect('coupled.' + name + '.chords', name + '_perf.chords')
             self.connect("coupled." + name + ".lengths", name + "_perf.lengths")
-            self.connect("coupled." + name + ".cos_sweep", name + "_perf.cos_sweep")
+            self.connect("coupled." + name + ".lengths_spanwise", name + "_perf.lengths_spanwise")
 
             # Connect parameters from the 'coupled' group to the total performance group.
             self.connect("coupled." + name + ".S_ref", "total_perf." + name + "_S_ref")