Source code for openaerostruct.structures.section_properties_wingbox
importnumpyasnpimportopenmdao.apiasom
[docs]classSectionPropertiesWingbox(om.ExplicitComponent):""" Compute geometric cross-section properties for the wingbox elements. See Chauhan et al. (https://doi.org/10.1007/978-3-319-97773-7_38) for more. Parameters ---------- streamwise_chords[ny-1] : numpy array Average streamwise chord lengths for each streamwise VLM panel. fem_chords[ny-1] : numpy array Effective chord lengths normal to the FEM elements. fem_twists[ny-1] : numpy array Twist angles in planes normal to the FEM elements. spar_thickness[ny-1] : numpy array Material thicknesses of the front and rear spars for each wingbox segment. skin_thickness[ny-1] : numpy array Material thicknesses of the top and bottom skins for each wingbox segment. t_over_c[ny-1] : numpy array Streamwise thickness-to-chord ratios for each wingbox segment. Returns ------- A[ny-1] : numpy array Cross-sectional area of each wingbox segment. A_enc[ny-1] : numpy array Cross-sectional enclosed area (measured using the material midlines) of each wingbox segment. A_int[ny-1] : numpy array Cross-sectional internal area of each wingbox segment (used for fuel volume). Iy[ny-1] : numpy array Second moment of area about the neutral axis parallel to the local y-axis (for each wingbox segment). Qz[ny-1] : numpy array First moment of area above the neutral axis parallel to the local z-axis (for each wingbox segment). Iz[ny-1] : numpy array Second moment of area about the neutral axis parallel to the local z-axis (for each wingbox segment). J[ny-1] : numpy array Torsion constants for each wingbox segment. htop[ny-1] : numpy array Distance to the point on the top skin that is the farthest away from the local-z neutral axis (for each wingbox segment). hbottom[ny-1] : numpy array Distance to the point on the bottom skin that is the farthest away from the local-z neutral axis (for each wingbox segment). hfront[ny-1] : numpy array Distance to the point on the front spar that is the farthest away from the local-y neutral axis (for each wingbox segment). hrear[ny-1] : numpy array Distance to the point on the rear spar that is the farthest away from the local-y neutral axis (for each wingbox segment). """definitialize(self):self.options.declare("surface",types=dict)defsetup(self):self.surface=surface=self.options["surface"]self.mesh=surface["mesh"]self.ny=self.mesh.shape[1]# original thickness-to-chord ratio of the airfoil provided by the userself.orig_wb_af_t_over_c=surface["original_wingbox_airfoil_t_over_c"]# airfoil coordinates provided by the userself.data_x_upper=surface["data_x_upper"]self.data_x_lower=surface["data_x_lower"]self.data_y_upper=surface["data_y_upper"]self.data_y_lower=surface["data_y_lower"]self.add_input("streamwise_chords",val=np.ones((self.ny-1)),units="m")self.add_input("fem_chords",val=np.ones((self.ny-1)),units="m")self.add_input("fem_twists",val=np.ones((self.ny-1)),units="deg")self.add_input("spar_thickness",val=np.ones((self.ny-1)),units="m")self.add_input("skin_thickness",val=np.ones((self.ny-1)),units="m")self.add_input("t_over_c",val=np.ones((self.ny-1)))self.add_output("A",val=np.ones((self.ny-1)),units="m**2")self.add_output("A_enc",val=np.ones((self.ny-1)),units="m**2")self.add_output("A_int",val=np.ones((self.ny-1)),units="m**2")self.add_output("Iy",val=np.ones((self.ny-1)),units="m**4")self.add_output("Qz",val=np.ones((self.ny-1)),units="m**3")self.add_output("Iz",val=np.ones((self.ny-1)),units="m**4")self.add_output("J",val=np.ones((self.ny-1)),units="m**4")self.add_output("htop",val=np.ones((self.ny-1)),units="m")self.add_output("hbottom",val=np.ones((self.ny-1)),units="m")self.add_output("hfront",val=np.ones((self.ny-1)),units="m")self.add_output("hrear",val=np.ones((self.ny-1)),units="m")self.declare_partials("*","*",method="cs")defcompute(self,inputs,outputs):# NOTE: In the code below, the x- and y-axes correspond to the element# local z- and y-axes, respectively.chord=inputs["fem_chords"]spar_thickness=inputs["spar_thickness"]skin_thickness=inputs["skin_thickness"]t_over_c_original=self.orig_wb_af_t_over_ct_over_c=inputs["t_over_c"]streamwise_chord=inputs["streamwise_chords"]theta=inputs["fem_twists"]# Scale data points with chorddata_x_upper=np.outer(self.data_x_upper,chord)data_y_upper=np.outer(self.data_y_upper,chord)data_x_lower=np.outer(self.data_x_lower,chord)data_y_lower=np.outer(self.data_y_lower,chord)# Scale y-coordinates by t/c design variable which is streamwise t/cdata_y_upper*=t_over_c/t_over_c_original*streamwise_chord/chorddata_y_lower*=t_over_c/t_over_c_original*streamwise_chord/chord# Compute enclosed area for torsion constant# This currently does not change with twist# Also compute internal area for internal volume calculation for fuelx_up_diff=data_x_upper[1:]-data_x_upper[:-1]x_low_diff=data_x_lower[1:]-data_x_lower[:-1]y_up_diff=data_y_upper[1:]-data_y_upper[:-1]y_low_diff=data_y_lower[1:]-data_y_lower[:-1]y_up_add=data_y_upper[1:]+data_y_upper[:-1]y_low_add=data_y_lower[1:]+data_y_lower[:-1]A_enc=(x_up_diff)*(y_up_add-skin_thickness)/2# area above 0 lineA_enc+=(x_low_diff)*(-y_low_add-skin_thickness)/2# area below 0 lineA_int=(x_up_diff)*(y_up_add-2*skin_thickness)/2# area above 0 lineA_int+=(x_low_diff)*(-y_low_add-2*skin_thickness)/2# area below 0 lineA_enc=np.sum(A_enc,axis=0)A_int=np.sum(A_int,axis=0)A_enc-=(data_y_upper[0]-data_y_lower[0])*spar_thickness/2# area of sparsA_enc-=(data_y_upper[-1]-data_y_lower[-1])*spar_thickness/2# area of sparsA_int-=(data_y_upper[0]-data_y_lower[0])*spar_thickness# area of sparsA_int-=(data_y_upper[-1]-data_y_lower[-1])*spar_thickness# area of sparsoutputs["A_enc"]=A_encoutputs["A_int"]=A_int# Compute perimeter to thickness ratio for torsion constant# This currently does not change with twistp_by_t_1=(x_up_diff**2+y_up_diff**2)**0.5/skin_thickness# length / thickness of capsp_by_t_2=(x_low_diff**2+y_low_diff**2)**0.5/skin_thickness# length / thickness of capsp_by_t=np.sum(p_by_t_1+p_by_t_2,axis=0)p_by_t+=(data_y_upper[0]-data_y_lower[0]-skin_thickness)/spar_thickness# length / thickness of sparsp_by_t+=(data_y_upper[-1]-data_y_lower[-1]-skin_thickness)/spar_thickness# length / thickness of spars# Torsion constantJ=4*A_enc**2/p_by_toutputs["J"]=J# Rotate the wingboxrot_mat=np.array([[np.cos(theta),np.sin(theta)],[-np.sin(theta),np.cos(theta)]])data_x_upper_2=rot_mat[0,0]*data_x_upper+rot_mat[0,1]*data_y_upperdata_y_upper_2=rot_mat[1,0]*data_x_upper+rot_mat[1,1]*data_y_upperdata_x_lower_2=rot_mat[0,0]*data_x_lower+rot_mat[0,1]*data_y_lowerdata_y_lower_2=rot_mat[1,0]*data_x_lower+rot_mat[1,1]*data_y_lowerdata_x_upper=data_x_upper_2.copy()data_y_upper=data_y_upper_2.copy()data_x_lower=data_x_lower_2.copy()data_y_lower=data_y_lower_2.copy()x_up_diff=data_x_upper[1:]-data_x_upper[:-1]x_low_diff=data_x_lower[1:]-data_x_lower[:-1]y_up_diff=data_y_upper[1:]-data_y_upper[:-1]y_low_diff=data_y_lower[1:]-data_y_lower[:-1]y_up_add=data_y_upper[1:]+data_y_upper[:-1]y_low_add=data_y_lower[1:]+data_y_lower[:-1]# Compute area moment of inertia about x axis# First compute centroid and areafirst_moment_area_upper=(y_up_add/2-(skin_thickness/2))*skin_thickness*x_up_diffupper_area=skin_thickness*x_up_difffirst_moment_area_lower=(y_low_add/2+(skin_thickness/2))*skin_thickness*x_low_difflower_area=skin_thickness*x_low_difffirst_moment_area_front_spar=((data_y_upper[0]-data_y_lower[0]-2*skin_thickness)*spar_thickness*(data_y_upper[0]+data_y_lower[0])/2)first_moment_area_rear_spar=((data_y_upper[-1]-data_y_lower[-1]-2*skin_thickness)*spar_thickness*(data_y_upper[-1]+data_y_lower[-1])/2)area_spars=((data_y_upper[0]-data_y_lower[0]-2*skin_thickness)+(data_y_upper[-1]-data_y_lower[-1]-2*skin_thickness))*spar_thicknessarea=np.sum(upper_area,axis=0)+np.sum(lower_area,axis=0)+area_sparsoutputs["A"]=areacentroid=(np.sum(first_moment_area_upper,axis=0)+np.sum(first_moment_area_lower,axis=0)+first_moment_area_front_spar+first_moment_area_rear_spar)/area# Then compute area moment of inertia for upward bending# This is calculated using derived analytical expression assuming linear interpolation between airfoil data pointsa=y_up_diff/x_up_diffb=(y_up_diff+skin_thickness)/2x2=x_up_diffI_horiz_1=2*((1.0/12.0*a**3*x2**4+1.0/3.0*a**2*x2**3*b+1.0/2.0*a*x2**2*b**2+1.0/3.0*b**3*x2))I_horiz_2=x2*skin_thickness*(y_up_add/2-skin_thickness/2-centroid)**2I_horiz=np.sum(I_horiz_1+I_horiz_2,axis=0)# Compute area moment of inertia about y axisa=-y_low_diff/x_low_diffb=(-y_low_diff+skin_thickness)/2x2=x_low_diffI_horiz+=np.sum(2*((1.0/12.0*a**3*x2**4+1.0/3.0*a**2*x2**3*b+1.0/2.0*a*x2**2*b**2+1.0/3.0*b**3*x2)),axis=0,)I_horiz+=np.sum(x2*skin_thickness*((-y_low_add)/2-skin_thickness/2+centroid)**2,axis=0)# Contribution from the forward sparI_horiz+=(1.0/12.0*spar_thickness*(data_y_upper[0]-data_y_lower[0]-2*skin_thickness)**3+spar_thickness*(data_y_upper[0]-data_y_lower[0]-2*skin_thickness)*((data_y_upper[0]+data_y_lower[0])/2-centroid)**2)# Contribution from the rear sparI_horiz+=(1.0/12.0*spar_thickness*(data_y_upper[-1]-data_y_lower[-1]-2*skin_thickness)**3+spar_thickness*(data_y_upper[-1]-data_y_lower[-1]-2*skin_thickness)*((data_y_upper[-1]+data_y_lower[-1])/2-centroid)**2)outputs["Iz"]=I_horiz# Compute the Q required for transverse shear stress due to upward bendingQ_upper=np.sum(((y_up_add/2-(skin_thickness/2))-centroid)*skin_thickness*x_up_diff,axis=0)Q_upper+=(data_y_upper[0]-skin_thickness-centroid)**2/2*(spar_thickness)Q_upper+=(data_y_upper[-1]-skin_thickness-centroid)**2/2*(spar_thickness)outputs["Qz"]=Q_upper# Compute area moment of inertia for backward bendingI_vert=0first_moment_area_front=((data_y_upper[0]-data_y_lower[0])*spar_thickness*(data_x_upper[0]+spar_thickness/2))first_moment_area_rear=((data_y_upper[-1]-data_y_lower[-1])*spar_thickness*(data_x_upper[-1]-spar_thickness/2))centroid_Ivert=(first_moment_area_front+first_moment_area_rear)/(((data_y_upper[0]-data_y_lower[0])+(data_y_upper[-1]-data_y_lower[-1]))*spar_thickness)I_vert+=(1.0/12.0*(data_y_upper[0]-data_y_lower[0])*spar_thickness**3+(data_y_upper[0]-data_y_lower[0])*spar_thickness*(centroid_Ivert-(data_x_upper[0]+spar_thickness/2))**2)I_vert+=(1.0/12.0*(data_y_upper[-1]-data_y_lower[-1])*spar_thickness**3+(data_y_upper[-1]-data_y_lower[-1])*spar_thickness*(data_x_upper[-1]-spar_thickness/2-centroid_Ivert)**2)# Add contribution of skinsI_vert+=2*(1.0/12.0*skin_thickness*(data_x_upper[-1]-data_x_upper[0]-2*spar_thickness)**3+skin_thickness*(data_x_upper[-1]-data_x_upper[0]-2*spar_thickness)*(centroid_Ivert-(data_x_upper[-1]+data_x_upper[0])/2)**2)outputs["Iy"]=I_vert# Distances for calculating max bending stresses (KS function used)ks_rho=500.0# Hard coded, see Martins and Poon 2005 for morefmax_upper=np.max(data_y_upper,axis=0)htop=fmax_upper+1/ks_rho*np.log(np.sum(np.exp(ks_rho*(data_y_upper-fmax_upper)),axis=0))-centroidfmax_lower=np.max(-data_y_lower,axis=0)hbottom=(fmax_lower+1/ks_rho*np.log(np.sum(np.exp(ks_rho*(-data_y_lower-fmax_lower)),axis=0))+centroid)hfront=centroid_Ivert-data_x_upper[0]hrear=data_x_upper[-1]-centroid_Ivertoutputs["htop"]=htopoutputs["hbottom"]=hbottomoutputs["hfront"]=hfrontoutputs["hrear"]=hrear