Fixed incorrect surface moment derivatives and added a moment output …

…to the component (#449)

* Moment deriv comments refactor

* Refactor Moment derivatives. Add comments to explain process. Fix bugs related derivatives wrt reference area, chords, and panel width. Add output for moment with derivatives

* Fixed a derivative. Black and Flake8

* Revert test changes. Fix comment errors. Minor derivative fix.

* Restore complex step tests. Fix incorrect S_ref derivative accumulation for non wing surfaces

* Fix area derivative one more time

* Add new test to properly check moment derivatives

* Removed debugging lines

* Removed comments indicating that there may be problem with the moment derivativesas they have now been resolved. Don't want to confuse people in the future.

* merged tests

* Update multiple_surfaces.rst

* Fix comment grammar and bump version

---------

Co-authored-by: Shugo Kaneko <[email protected]>
Co-authored-by: Safa Bakhshi <[email protected]>

openaerostruct/__init__.py

@@ -1 +1 @@
-__version__ = "2.9.0"
+__version__ = "2.9.1"

openaerostruct/docs/advanced_features/multiple_surfaces.rst

@@ -5,6 +5,7 @@ Multiple Lifting Surfaces
 
 It's easily possible to simulate multiple lifting surfaces simultaneously in OpenAeroStruct.
 The most straightforward example is a wing and a tail for a conventional airplane, as shown below, though OpenAeroStruct can handle any arbitrary collection of lifting surfaces.
+Note that the moment coefficient is always normalized with respect to the MAC of the first surface in the list; therefore, the main lifting surface (i.e., wing) should always be the first entry in the list.
 
 .. literalinclude:: /../../tests/integration_tests/test_multiple_aero_analysis.py
     :start-after: docs checkpoint 0

openaerostruct/functionals/moment_coefficient.py

@@ -5,7 +5,7 @@
 
 class MomentCoefficient(om.ExplicitComponent):
     """
-    Compute the coefficient of moment (CM) for the entire aircraft.
+    Compute the coefficient of moment (CM) and moment (M) for the entire aircraft.
 
     Parameters
     ----------
@@ -34,6 +34,8 @@ class MomentCoefficient(om.ExplicitComponent):
 
     Returns
     -------
+    M[3] : numpy array
+        The moment around the x-, y-, and z-axes at the cg point.
     CM[3] : numpy array
         The coefficient of moment around the x-, y-, and z-axes at the cg point.
     """
@@ -59,6 +61,7 @@ def setup(self):
         self.add_input("S_ref_total", val=1.0, units="m**2", tags=["mphys_input"])
 
         self.add_output("CM", val=np.ones((3)), tags=["mphys_result"])
+        self.add_output("M", val=np.ones((3)), units="N*m", tags=["mphys_result"])
 
         self.declare_partials(of="*", wrt="*")
 
@@ -112,12 +115,15 @@ def compute(self, inputs, outputs):
             # normalize CM
             if j == 0:
                 self.MAC_wing = MAC
+                self.S_ref_wing = S_ref
 
         self.M = M
 
+        # Output the moment vector
+        outputs["M"] = M
+
         # Compute the normalized CM
-        rho = inputs["rho"]
-        outputs["CM"] = M / (0.5 * rho * inputs["v"] ** 2 * inputs["S_ref_total"] * self.MAC_wing)
+        outputs["CM"] = M / (0.5 * inputs["rho"] * inputs["v"] ** 2 * inputs["S_ref_total"] * self.MAC_wing)
 
     def compute_partials(self, inputs, partials):
         cg = inputs["cg"]
@@ -128,12 +134,14 @@ def compute_partials(self, inputs, partials):
         # Cached values
         M = self.M
         MAC_wing = self.MAC_wing
+        S_ref_wing = self.S_ref_wing
 
+        # Scaling factor of one over the dynamic pressure times sum of reference areas times the wing MAC
         fact = 1.0 / (0.5 * rho * v**2 * S_ref_total * MAC_wing)
 
-        partials["CM", "rho"] = -M * fact**2 * 0.5 * v**2 * S_ref_total * MAC_wing
-        partials["CM", "v"] = -M * fact**2 * rho * v * S_ref_total * MAC_wing
-        partials["CM", "S_ref_total"] = -M * fact**2 * 0.5 * rho * v**2 * MAC_wing
+        partials["CM", "rho"] = -M * fact / rho
+        partials["CM", "v"] = -2 * M * fact / v
+        partials["CM", "S_ref_total"] = -M * fact / S_ref_total
 
         partials["CM", "cg"][:] = 0.0
 
@@ -156,7 +164,8 @@ def compute_partials(self, inputs, partials):
             panel_chords = (chords[1:] + chords[:-1]) * 0.5
             MAC = 1.0 / S_ref * np.sum(panel_chords**2 * widths)
 
-            # This transformation is used for multiple derivatives
+            # This produces a bi-diagonal matrix for the derivative of panel_chords with respect to chords
+            # This transformation matrix is further used for multiple derivatives later
             dpc_dc = np.zeros((ny - 1, ny))
             idx = np.arange(ny - 1)
             dpc_dc[idx, idx] = 0.5
@@ -174,34 +183,72 @@ def compute_partials(self, inputs, partials):
                 dMAC_dw *= 2.0
                 dMAC_dS *= 2.0
 
-            # diff derivs
+            # Compute the bound vortex(quarter chord) points at mid-panel
             pts = (b_pts[:, 1:, :] + b_pts[:, :-1, :]) * 0.5
+
+            # Compute the vectors between the cg and the mid-panel bound vortex points
             diff = pts - cg
 
+            # Compute the cross product of the panel bound vortex vectors from cg and the panel forces
             c = np.cross(diff, sec_forces, axis=2)
+
+            # Compute the spanwise moment vector distribution by summing over each resulting column
             moment = np.sum(c, axis=0)
 
+            # Compute the derviative of the moment vectors(c) with respect to the diff vectors(a) multiplied by -1
             dcda = np.zeros((3, nx - 1, ny - 1, 3))
+
+            # Compute the derivative wrt to the first element of diff
             dcda[0, :, :, 1] = sec_forces[:, :, 2]
             dcda[0, :, :, 2] = -sec_forces[:, :, 1]
+
+            # Compute the derivative wrt to the second element of diff
             dcda[1, :, :, 0] = -sec_forces[:, :, 2]
             dcda[1, :, :, 2] = sec_forces[:, :, 0]
+
+            # Compute the derivative wrt to the third element of diff
             dcda[2, :, :, 0] = sec_forces[:, :, 1]
             dcda[2, :, :, 1] = -sec_forces[:, :, 0]
 
+            # Compute the derviative of the moment vectors(c) with respect to the sec_forces vectors(b) multiplied by -1
             dcdb = np.zeros((3, nx - 1, ny - 1, 3))
+
+            # Compute the derivative wrt to the first element of sec_forces
             dcdb[0, :, :, 1] = -diff[:, :, 2]
             dcdb[0, :, :, 2] = diff[:, :, 1]
+
+            # Compute the derivative wrt to the second element of sec_forces
             dcdb[1, :, :, 0] = diff[:, :, 2]
             dcdb[1, :, :, 2] = -diff[:, :, 0]
+
+            # Compute the derivative wrt to the third element of sec_forces
             dcdb[2, :, :, 0] = -diff[:, :, 1]
             dcdb[2, :, :, 1] = diff[:, :, 0]
 
+            # Compute derivative of CM wrt to the sec_forces of the section by reshaping to 3 rows and multiplying by fact.
             partials["CM", name + "_sec_forces"] += dcdb.reshape((3, 3 * (nx - 1) * (ny - 1))) * fact
 
+            # Compute derivative of M wrt to the sec_forces of the section by reshaping to 3 rows
+            partials["M", name + "_sec_forces"] += dcdb.reshape((3, 3 * (nx - 1) * (ny - 1)))
+
+            # Project the derviative of the moment vectors(c) with respect to the diff vectors(a)
+            # onto the derivative of mid-panel chord distribution wrt to the chord distribution giving the derivative of
+            # the moment vectors(c) with respect to the chord distribution(dc_dchord). This works because the diff
+            # vectors are difference between the mid-panel bound vortex(quarter chord) points and cg which is static in this derivative.
+            # The spanwise component of the mid-panel bound vortex(quarter chord) points have the same derivatrive wrt to the chord
+            # distribution as the mid-panel chord distribution does wrt the the chord distribution.
             dc_dchord = np.einsum("ijkl,km->ijml", dcda, dpc_dc)
+
+            # Compute the derivative of CM wrt to the bound vortex points(b_pts) by reshaping dc_dchord to three rows
+            # and multiplying by fact.
             partials["CM", name + "_b_pts"] += dc_dchord.reshape((3, 3 * (nx - 1) * ny)) * fact
 
+            # Compute the derivative of M wrt to the bound vortex points(b_pts) by reshaping dc_dchord to three rows
+            partials["M", name + "_b_pts"] += dc_dchord.reshape((3, 3 * (nx - 1) * ny))
+
+            # Reduce the derivative of the moment vectors(c) with respect to the diff vectors(a) by summing over all
+            # chordwise and spanwise panels(j and k). Reduces to a 3x3 matrix for the whole surface by summing over all
+            # panels.
             dcda = np.einsum("ijkl->il", dcda)
 
             # If the surface is symmetric, set the x- and z-direction moments
@@ -216,28 +263,62 @@ def compute_partials(self, inputs, partials):
                 partials["CM", name + "_b_pts"][0, :] = 0.0
                 partials["CM", name + "_b_pts"][1, :] *= 2.0
                 partials["CM", name + "_b_pts"][2, :] = 0.0
+                partials["M", name + "_sec_forces"][0, :] = 0.0
+                partials["M", name + "_sec_forces"][1, :] *= 2.0
+                partials["M", name + "_sec_forces"][2, :] = 0.0
+                partials["M", name + "_b_pts"][0, :] = 0.0
+                partials["M", name + "_b_pts"][1, :] *= 2.0
+                partials["M", name + "_b_pts"][2, :] = 0.0
                 dcda[0, :] = 0.0
                 dcda[1, :] *= 2.0
                 dcda[2, :] = 0.0
 
+            # Compute the derivative of CM wrt to the cg position which is negative dcda since diff = pts - cg times fact
+            # Accumlate the derivative over each surface as the total moment vector is sum over all surfaces.
             partials["CM", "cg"] -= dcda * fact
 
+            # Compute the derivative of M wrt to the cg position which is negative dcda since diff = pts - cg
+            # Accumlate the derivative over each surface as the total moment vector is sum over all surfaces.
+            partials["M", "cg"] -= dcda
+
+            # Compute the total surface moment vector by summing spanwise
             M_j = np.sum(moment, axis=0)
-            term = fact / MAC
-            partials["CM", name + "_chords"] = -np.outer(M_j * term, dMAC_dc)
-            partials["CM", name + "_widths"] = -np.outer(M_j * term, dMAC_dw)
-            partials["CM", name + "_S_ref"] = -np.outer(M_j, dMAC_dS * term)
 
             # For first surface, we need to save the deriv results
             if j == 0:
+                # Compute a term by dividing fact by MAC. Note that MAC is the mean aerodynamic chord for the surface and
+                # the MAC_wing terms already factored into fact is of the main wing surface
+                term = fact / MAC
+
+                # Compute the derivative of CM wrt to the chord distribution by taking the negative outer product of the
+                # moment vector(M_j) time the term with the derivative of MAC wrt to the chord distribution. We only do
+                # this for the main wing since CM only depends on the MAC of the main wing and the chord distribution of
+                # the main wing is the only chord distribution of all the surfaces that can impact the MAC of the main wing.
+                partials["CM", name + "_chords"] = -np.outer(M_j * term, dMAC_dc)
+
+                # Compute the derivative of CM wrt to the width distribution by taking the negative outer product of the
+                # moment vector(M_j) time the term with the derivative of MAC wrt to the width distribution. We only do
+                # this for the main wing since CM only depends on the MAC of the main wing and the panel width distribution of
+                # the main wing is the only panel width distribution of all the surfaces that can impact the MAC of the main wing.
+                partials["CM", name + "_widths"] = -np.outer(M_j * term, dMAC_dw)
+
+                # Compute the derivative of CM wrt to the surface S_ref by taking the negative outer product of the
+                # moment vector(M_j) time the term with the derivative of MAC wrt to the surface S_ref. The CM depends on
+                # the total references area of all surfaces including the main wing and the MAC of them main wing itself
+                # As result, this derivative has two parts only for the main wing.
+                # partials["CM", name + "_S_ref"] = -np.outer(M_j, dMAC_dS * term)
+                partials["CM", name + "_S_ref"] = np.outer(M_j * fact, (1 / S_ref))
+
+                # Cache the main wing's MAC derivatives
                 base_name = name
                 base_dMAC_dc = dMAC_dc
                 base_dMAC_dw = dMAC_dw
-                base_dMAC_dS = dMAC_dS
-
+                # base_dMAC_dS = dMAC_dS
             else:
                 # Apply this surface's portion of the moment to the MAC_wing term.
-                term = fact / (MAC_wing * MAC)
+                # We need to do this because CM is normalized by the MAC of the main wing
+                term = fact / MAC_wing
                 partials["CM", base_name + "_chords"] -= np.outer(M_j * term, base_dMAC_dc)
                 partials["CM", base_name + "_widths"] -= np.outer(M_j * term, base_dMAC_dw)
-                partials["CM", base_name + "_S_ref"] -= np.outer(M_j, base_dMAC_dS * term)
+                # partials["CM", base_name + "_S_ref"] -= np.outer(M_j, base_dMAC_dS * term)
+                partials["CM", base_name + "_S_ref"] += np.outer(M_j * fact, (1 / S_ref_wing))

tests/functionals_tests/test_moment_coefficient.py

@@ -16,10 +16,8 @@ def test(self):
 
         comp = MomentCoefficient(surfaces=surfaces)
 
-        run_test(self, comp)
+        run_test(self, comp, complex_flag=True, method="cs")
 
-    # This is known to have some issues for sufficiently small values of S_ref_total
-    # There is probably a derivative bug somewhere in the moment_coefficient.py calcs
     def test2(self):
         surfaces = get_default_surfaces()
 
@@ -30,13 +28,15 @@ def test2(self):
         indep_var_comp = om.IndepVarComp()
 
         indep_var_comp.add_output("S_ref_total", val=1e4, units="m**2")
+        indep_var_comp.add_output("cg", val=np.array([-10.0, 10.0, -10.0]), units="m")
 
         group.add_subsystem("moment_calc", comp)
         group.add_subsystem("indep_var_comp", indep_var_comp)
 
         group.connect("indep_var_comp.S_ref_total", "moment_calc.S_ref_total")
+        group.connect("indep_var_comp.cg", "moment_calc.cg")
 
-        run_test(self, group)
+        run_test(self, group, complex_flag=True, method="cs")
 
 
 if __name__ == "__main__":