Source code for openaerostruct.aerodynamics.wave_drag
importnumpyasnpimportopenmdao.apiasom
[docs]classWaveDrag(om.ExplicitComponent):""" Compute the wave drag if the with_wave option is True. If not, the CDw is 0. This component exists for each lifting surface. Parameters ---------- Mach_number : float Mach number. widths[ny-1] : numpy array 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. chords[ny] : numpy array The chord length of each mesh slice. This is dimension ny rather than ny-1 which would be expected for chord length of each VLM panel. t_over_c[ny-1] : numpy array The streamwise thickness-to-chord ratio of each VLM panel. Returns ------- CDw : float Wave drag coefficient for the lifting surface computed using equations based on the Korn equation """definitialize(self):self.options.declare("surface",types=dict)self.options.declare("with_wave",types=bool)defsetup(self):self.surface=surface=self.options["surface"]self.with_wave=surface["with_wave"]# Thickness over chord for the airfoilself.ka=0.95# airfoil technology level (for NASA SC airfoil)ny=surface["mesh"].shape[1]self.add_input("Mach_number",val=1.6,tags=["mphys_input"])self.add_input("widths",val=np.ones((ny-1))*0.2,units="m",tags=["mphys_coupling"])self.add_input("lengths_spanwise",val=np.arange((ny-1))+1.0,units="m",tags=["mphys_coupling"])# set to np.arange so that d_CDw_d_chords is nonzeroself.add_input("CL",val=0.33)self.add_input("chords",val=np.ones((ny)),units="m",tags=["mphys_coupling"])self.add_input("t_over_c",val=np.arange((ny-1)),tags=["mphys_input"])self.add_output("CDw",val=0.0)self.declare_partials("CDw","*")self.set_check_partial_options(wrt="*",method="cs",step=1e-50)defcompute(self,inputs,outputs):ifself.with_wave:t_over_c=inputs["t_over_c"]widths=inputs["widths"]cos_sweep=widths/inputs["lengths_spanwise"]M=inputs["Mach_number"]chords=inputs["chords"]CL=inputs["CL"]panel_mid_chords=(chords[:-1]+chords[1:])/2.0panel_areas=panel_mid_chords*widthssum_panel_areas=np.sum(panel_areas)avg_cos_sweep=np.sum(cos_sweep*panel_areas)/sum_panel_areas# weighted average of 1/4 chord sweepavg_t_over_c=np.sum(t_over_c*panel_areas)/sum_panel_areas# weighted average of streamwise t/cMDD=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)ifnp.real(M)>np.real(Mcrit):outputs["CDw"]=20*(M-Mcrit)**4else:outputs["CDw"]=0.0ifself.surface["symmetry"]:outputs["CDw"]*=2else:outputs["CDw"]=0.0defcompute_partials(self,inputs,partials):"""Jacobian for wave drag."""# Explicitly zero out the partials to begin, we can't assume the input partials arrays contain zeros alreadypartials["CDw","CL"][:]=0.0partials["CDw","lengths_spanwise"][:]=0.0partials["CDw","widths"][:]=0.0partials["CDw","Mach_number"][:]=0.0partials["CDw","chords"][:]=0.0partials["CDw","t_over_c"][:]=0.0ifself.with_wave:ny=self.surface["mesh"].shape[1]t_over_c=inputs["t_over_c"]widths=inputs["widths"]lengths_spanwise=inputs["lengths_spanwise"]cos_sweep=widths/lengths_spanwiseM=inputs["Mach_number"]chords=inputs["chords"]CL=inputs["CL"]panel_mid_chords=(chords[:-1]+chords[1:])/2.0panel_areas=panel_mid_chords*widthssum_panel_areas=np.sum(panel_areas)avg_cos_sweep=np.sum(cos_sweep*panel_areas)/sum_panel_areasavg_t_over_c=np.sum(t_over_c*panel_areas)/sum_panel_areasMDD=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)ifnp.real(M)>np.real(Mcrit):dCDwdMDD=-80*(M-Mcrit)**3dMDDdCL=-1.0/(10*avg_cos_sweep**3)dMDDdavg=(-10*self.ka*avg_cos_sweep**2+20*avg_t_over_c*avg_cos_sweep+3*CL)/(10*avg_cos_sweep**4)dMDDdtoc=-1.0/(avg_cos_sweep**2)dtocavgdtoc=panel_areas/sum_panel_areasccos=np.sum(widths*panel_mid_chords)ccos2w=np.sum(panel_mid_chords*widths**2/lengths_spanwise)davgdcos=(2*panel_mid_chords*widths/lengths_spanwise/ccos-panel_mid_chords*ccos2w/ccos**2)dtocdcos=(panel_mid_chords*t_over_c/ccos-panel_mid_chords*np.sum(panel_mid_chords*widths*t_over_c)/ccos**2)davgdw=-1*panel_mid_chords*widths**2/lengths_spanwise**2/ccosdavgdc=widths**2/lengths_spanwise/ccos-widths*ccos2w/ccos**2dtocdc=t_over_c*widths/ccos-widths*np.sum(panel_mid_chords*widths*t_over_c)/ccos**2dcdchords=np.zeros((ny-1,ny))i,j=np.indices(dcdchords.shape)dcdchords[i==j]=0.5dcdchords[i==j-1]=0.5partials["CDw","Mach_number"]=-1*dCDwdMDDpartials["CDw","CL"]=dCDwdMDD*dMDDdCLpartials["CDw","lengths_spanwise"]=dCDwdMDD*dMDDdavg*davgdwpartials["CDw","widths"]=dCDwdMDD*dMDDdavg*davgdcos+dCDwdMDD*dMDDdtoc*dtocdcospartials["CDw","chords"]=dCDwdMDD*dMDDdavg*np.matmul(davgdc,dcdchords)+dCDwdMDD*dMDDdtoc*np.matmul(dtocdc,dcdchords)partials["CDw","t_over_c"]=dCDwdMDD*dMDDdtoc*dtocavgdtocifself.surface["symmetry"]:partials["CDw","CL"][0,:]*=2partials["CDw","lengths_spanwise"][0,:]*=2partials["CDw","widths"][0,:]*=2partials["CDw","Mach_number"][0,:]*=2partials["CDw","chords"][0,:]*=2partials["CDw","t_over_c"][0,:]*=2