Source code for openaerostruct.functionals.moment_coefficient

import numpy as np

import openmdao.api as om


[docs] class MomentCoefficient(om.ExplicitComponent): """ Compute the coefficient of moment (CM) and moment (M) for the entire aircraft. Parameters ---------- 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. chords[ny] : numpy array The chordwise length of the entire airfoil following the camber line. S_ref : float The reference area of the lifting surface. sec_forces[nx-1, ny-1, 3] : numpy array Contains the sectional forces acting on each panel. Stored in Fortran order (only relevant with more than one chordwise panel). cg[3] : numpy array The x, y, z coordinates of the center of gravity for the entire aircraft. v : float Freestream air velocity in m/s. rho : float Air density in kg/m^3. S_ref_total : float Total surface area of the aircraft based on the sum of individual surface areas. 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. """ def initialize(self): self.options.declare("surfaces", types=list) def setup(self): for surface in self.options["surfaces"]: name = surface["name"] nx = surface["mesh"].shape[0] ny = surface["mesh"].shape[1] self.add_input(name + "_b_pts", val=np.ones((nx - 1, ny, 3)), units="m", tags=["mphys_coupling"]) self.add_input(name + "_widths", val=np.ones((ny - 1)), units="m", tags=["mphys_coupling"]) self.add_input(name + "_chords", val=np.ones((ny)), units="m", tags=["mphys_coupling"]) self.add_input(name + "_S_ref", val=1.0, units="m**2", tags=["mphys_coupling"]) self.add_input(name + "_sec_forces", val=np.ones((nx - 1, ny - 1, 3)), units="N", tags=["mphys_coupling"]) self.add_input("cg", val=np.ones((3)), units="m", tags=["mphys_input"]) self.add_input("v", val=10.0, units="m/s", tags=["mphys_input"]) self.add_input("rho", val=3.0, units="kg/m**3", tags=["mphys_input"]) 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="*") def compute(self, inputs, outputs): cg = inputs["cg"] M = np.zeros((3)) # Loop through each surface and find its contributions to the moment # of the aircraft based on the section forces and their location for j, surface in enumerate(self.options["surfaces"]): name = surface["name"] b_pts = inputs[name + "_b_pts"] widths = inputs[name + "_widths"] chords = inputs[name + "_chords"] S_ref = inputs[name + "_S_ref"] sec_forces = inputs[name + "_sec_forces"] # Compute the average chord for each panel and then the # mean aerodynamic chord (MAC) based on these chords and the # computed area panel_chords = (chords[1:] + chords[:-1]) * 0.5 MAC = 1.0 / S_ref * np.sum(panel_chords**2 * widths) # If the surface is symmetric, then the previously computed MAC # is half what it should be if surface["symmetry"]: MAC *= 2.0 # Get the moment arm acting on each panel, relative to the cg pts = (b_pts[:, 1:, :] + b_pts[:, :-1, :]) * 0.5 diff = pts - cg # Compute the moment based on the previously computed moment # arm and the section forces moment = np.sum(np.cross(diff, sec_forces, axis=2), axis=0) # If the surface is symmetric, set the x- and z-direction moments # to 0 and double the y-direction moment if surface["symmetry"]: moment[:, 0] = 0.0 moment[:, 1] *= 2.0 moment[:, 2] = 0.0 # Note: a scalar can be factored from a cross product, so I moved the division by MAC # down here for efficiency of calc and derivs. M = M + np.sum(moment, axis=0) # For the first (main) lifting surface, we save the MAC to correctly # 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 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"] rho = inputs["rho"] S_ref_total = inputs["S_ref_total"] v = inputs["v"] # 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 / rho partials["CM", "v"] = -2 * M * fact / v partials["CM", "S_ref_total"] = -M * fact / S_ref_total partials["CM", "cg"][:] = 0.0 # Loop through each surface. for j, surface in enumerate(self.options["surfaces"]): name = surface["name"] nx = surface["mesh"].shape[0] ny = surface["mesh"].shape[1] partials["CM", name + "_sec_forces"][:] = 0.0 partials["CM", name + "_b_pts"][:] = 0.0 b_pts = inputs[name + "_b_pts"] widths = inputs[name + "_widths"] chords = inputs[name + "_chords"] S_ref = inputs[name + "_S_ref"] sec_forces = inputs[name + "_sec_forces"] # MAC derivs panel_chords = (chords[1:] + chords[:-1]) * 0.5 MAC = 1.0 / S_ref * np.sum(panel_chords**2 * widths) # 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 dpc_dc[idx, idx + 1] = 0.5 dMAC_dc = (2.0 / S_ref) * np.einsum("i,ij", panel_chords * widths, dpc_dc) dMAC_dw = (1.0 / S_ref) * panel_chords**2 dMAC_dS = -MAC / S_ref # If the surface is symmetric, then the previously computed MAC # is half what it should be if surface["symmetry"]: MAC *= 2.0 dMAC_dc *= 2.0 dMAC_dw *= 2.0 dMAC_dS *= 2.0 # 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 # to 0 and double the y-direction moment if surface["symmetry"]: moment[:, 0] = 0.0 moment[:, 1] *= 2.0 moment[:, 2] = 0.0 partials["CM", name + "_sec_forces"][0, :] = 0.0 partials["CM", name + "_sec_forces"][1, :] *= 2.0 partials["CM", name + "_sec_forces"][2, :] = 0.0 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) # 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 else: # Apply this surface's portion of the moment to the MAC_wing term. # 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 * fact, (1 / S_ref_wing))