From f79f9cc5994fade2658e161bf0a0cb763fdb0e80 Mon Sep 17 00:00:00 2001 From: Safa Bakhshi Date: Wed, 18 Dec 2024 13:09:28 -0800 Subject: [PATCH] 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 <49300827+kanekosh@users.noreply.github.com> Co-authored-by: Safa Bakhshi --- openaerostruct/__init__.py | 2 +- .../advanced_features/multiple_surfaces.rst | 1 + .../functionals/moment_coefficient.py | 113 +++++++++++++++--- .../test_moment_coefficient.py | 8 +- 4 files changed, 103 insertions(+), 21 deletions(-) diff --git a/openaerostruct/__init__.py b/openaerostruct/__init__.py index 43ce13db0..334087f85 100644 --- a/openaerostruct/__init__.py +++ b/openaerostruct/__init__.py @@ -1 +1 @@ -__version__ = "2.9.0" +__version__ = "2.9.1" diff --git a/openaerostruct/docs/advanced_features/multiple_surfaces.rst b/openaerostruct/docs/advanced_features/multiple_surfaces.rst index fd1842bfe..0c0f07453 100644 --- a/openaerostruct/docs/advanced_features/multiple_surfaces.rst +++ b/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 diff --git a/openaerostruct/functionals/moment_coefficient.py b/openaerostruct/functionals/moment_coefficient.py index cdd7192e5..38fa63845 100644 --- a/openaerostruct/functionals/moment_coefficient.py +++ b/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)) diff --git a/tests/functionals_tests/test_moment_coefficient.py b/tests/functionals_tests/test_moment_coefficient.py index 7224bb35d..849d68d01 100644 --- a/tests/functionals_tests/test_moment_coefficient.py +++ b/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__":