Source code for openaerostruct.functionals.moment_coefficient
importnumpyasnpimportopenmdao.apiasom
[docs]classMomentCoefficient(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. """definitialize(self):self.options.declare("surfaces",types=list)defsetup(self):forsurfaceinself.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="*")defcompute(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 locationforj,surfaceinenumerate(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 areapanel_chords=(chords[1:]+chords[:-1])*0.5MAC=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 beifsurface["symmetry"]:MAC*=2.0# Get the moment arm acting on each panel, relative to the cgpts=(b_pts[:,1:,:]+b_pts[:,:-1,:])*0.5diff=pts-cg# Compute the moment based on the previously computed moment# arm and the section forcesmoment=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 momentifsurface["symmetry"]:moment[:,0]=0.0moment[:,1]*=2.0moment[:,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 CMifj==0:self.MAC_wing=MACself.S_ref_wing=S_refself.M=M# Output the moment vectoroutputs["M"]=M# Compute the normalized CMoutputs["CM"]=M/(0.5*inputs["rho"]*inputs["v"]**2*inputs["S_ref_total"]*self.MAC_wing)defcompute_partials(self,inputs,partials):cg=inputs["cg"]rho=inputs["rho"]S_ref_total=inputs["S_ref_total"]v=inputs["v"]# Cached valuesM=self.MMAC_wing=self.MAC_wingS_ref_wing=self.S_ref_wing# Scaling factor of one over the dynamic pressure times sum of reference areas times the wing MACfact=1.0/(0.5*rho*v**2*S_ref_total*MAC_wing)partials["CM","rho"]=-M*fact/rhopartials["CM","v"]=-2*M*fact/vpartials["CM","S_ref_total"]=-M*fact/S_ref_totalpartials["CM","cg"][:]=0.0# Loop through each surface.forj,surfaceinenumerate(self.options["surfaces"]):name=surface["name"]nx=surface["mesh"].shape[0]ny=surface["mesh"].shape[1]partials["CM",name+"_sec_forces"][:]=0.0partials["CM",name+"_b_pts"][:]=0.0b_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 derivspanel_chords=(chords[1:]+chords[:-1])*0.5MAC=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 laterdpc_dc=np.zeros((ny-1,ny))idx=np.arange(ny-1)dpc_dc[idx,idx]=0.5dpc_dc[idx,idx+1]=0.5dMAC_dc=(2.0/S_ref)*np.einsum("i,ij",panel_chords*widths,dpc_dc)dMAC_dw=(1.0/S_ref)*panel_chords**2dMAC_dS=-MAC/S_ref# If the surface is symmetric, then the previously computed MAC# is half what it should beifsurface["symmetry"]:MAC*=2.0dMAC_dc*=2.0dMAC_dw*=2.0dMAC_dS*=2.0# Compute the bound vortex(quarter chord) points at mid-panelpts=(b_pts[:,1:,:]+b_pts[:,:-1,:])*0.5# Compute the vectors between the cg and the mid-panel bound vortex pointsdiff=pts-cg# Compute the cross product of the panel bound vortex vectors from cg and the panel forcesc=np.cross(diff,sec_forces,axis=2)# Compute the spanwise moment vector distribution by summing over each resulting columnmoment=np.sum(c,axis=0)# Compute the derviative of the moment vectors(c) with respect to the diff vectors(a) multiplied by -1dcda=np.zeros((3,nx-1,ny-1,3))# Compute the derivative wrt to the first element of diffdcda[0,:,:,1]=sec_forces[:,:,2]dcda[0,:,:,2]=-sec_forces[:,:,1]# Compute the derivative wrt to the second element of diffdcda[1,:,:,0]=-sec_forces[:,:,2]dcda[1,:,:,2]=sec_forces[:,:,0]# Compute the derivative wrt to the third element of diffdcda[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 -1dcdb=np.zeros((3,nx-1,ny-1,3))# Compute the derivative wrt to the first element of sec_forcesdcdb[0,:,:,1]=-diff[:,:,2]dcdb[0,:,:,2]=diff[:,:,1]# Compute the derivative wrt to the second element of sec_forcesdcdb[1,:,:,0]=diff[:,:,2]dcdb[1,:,:,2]=-diff[:,:,0]# Compute the derivative wrt to the third element of sec_forcesdcdb[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 rowspartials["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 rowspartials["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 momentifsurface["symmetry"]:moment[:,0]=0.0moment[:,1]*=2.0moment[:,2]=0.0partials["CM",name+"_sec_forces"][0,:]=0.0partials["CM",name+"_sec_forces"][1,:]*=2.0partials["CM",name+"_sec_forces"][2,:]=0.0partials["CM",name+"_b_pts"][0,:]=0.0partials["CM",name+"_b_pts"][1,:]*=2.0partials["CM",name+"_b_pts"][2,:]=0.0partials["M",name+"_sec_forces"][0,:]=0.0partials["M",name+"_sec_forces"][1,:]*=2.0partials["M",name+"_sec_forces"][2,:]=0.0partials["M",name+"_b_pts"][0,:]=0.0partials["M",name+"_b_pts"][1,:]*=2.0partials["M",name+"_b_pts"][2,:]=0.0dcda[0,:]=0.0dcda[1,:]*=2.0dcda[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 spanwiseM_j=np.sum(moment,axis=0)# For first surface, we need to save the deriv resultsifj==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 surfaceterm=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 derivativesbase_name=namebase_dMAC_dc=dMAC_dcbase_dMAC_dw=dMAC_dw# base_dMAC_dS = dMAC_dSelse:# 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 wingterm=fact/MAC_wingpartials["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))