diff --git a/openaerostruct/aerodynamics/compressible_states.py b/openaerostruct/aerodynamics/compressible_states.py index 5e9cf29e6..308cc052c 100644 --- a/openaerostruct/aerodynamics/compressible_states.py +++ b/openaerostruct/aerodynamics/compressible_states.py @@ -1,6 +1,7 @@ """ Class definition for CompressibleVLMStates. """ + import openmdao.api as om from openaerostruct.aerodynamics.get_vectors import GetVectors diff --git a/openaerostruct/aerodynamics/mesh_point_forces.py b/openaerostruct/aerodynamics/mesh_point_forces.py index 4699d58c7..0188b25e4 100644 --- a/openaerostruct/aerodynamics/mesh_point_forces.py +++ b/openaerostruct/aerodynamics/mesh_point_forces.py @@ -1,6 +1,7 @@ """ Class definition for the MeshPointForces component. """ + import numpy as np import openmdao.api as om diff --git a/openaerostruct/docs/advanced_features.rst b/openaerostruct/docs/advanced_features.rst index 2beee7262..8f265e425 100644 --- a/openaerostruct/docs/advanced_features.rst +++ b/openaerostruct/docs/advanced_features.rst @@ -20,3 +20,4 @@ These examples show some advanced features in OpenAeroStruct. advanced_features/openconcept_coupling.rst advanced_features/customizing_prob_setup.rst advanced_features/mphys_coupling.rst + advanced_features/composites_walkthrough.rst diff --git a/openaerostruct/docs/advanced_features/composites_walkthrough.rst b/openaerostruct/docs/advanced_features/composites_walkthrough.rst new file mode 100644 index 000000000..151f07c07 --- /dev/null +++ b/openaerostruct/docs/advanced_features/composites_walkthrough.rst @@ -0,0 +1,178 @@ +.. _Composites Walkthrough: + +A walkthrough of the composites model +===================================== + +This page will walk you through the composites model in OpenAeroStruct. +The composites module allows you to define composite material properties and laminate layups for the wing structure. + +.. literalinclude:: ../composite_wingbox_mpt_opt_example.py + :start-after: checkpoint 7 + :end-before: checkpoint 8 + +Here, ``useComposite`` is a boolean variable that is set to True to enable the composites model and ``safety_factor`` is the factor of safety used for determining the Tsai-Wu based failure +criteria of the composite material. The composite material properties are defined using the following variables: + +- ``ply_angles`` is a list of the angles of the plies with respect to the x-axis. +- ``ply_fractions`` is a list of the ply fractions of the plies. (Should be the same length as ``ply_angles``, with the sum of the fractions equal to 1). +- ``E1`` is the modulus of elasticity in the fiber direction. +- ``E2`` is the modulus of elasticity in the transverse direction. +- ``G12`` is the shear modulus. +- ``nu12`` is the Poisson's ratio. +- ``sigma_t1`` is the tensile strength in the fiber direction. +- ``sigma_c1`` is the compressive strength in the fiber direction. +- ``sigma_t2`` is the tensile strength in the transverse direction. +- ``sigma_c2`` is the compressive strength in the transverse direction. +- ``sigma_12max`` is the maximum shear strength. + +.. note:: + The composites failure model doesn't use the ``strength_factor_for_upper_skin`` option from the surface dictionary. + If you want to apply a knockdown factor on the compressive strength to account for buckling, you should scale down the values of ``sigma_c1`` and ``sigma_c2``. + +Currently, the moduli of elasticity of the entire FEM spatial beam model are assumed to be isotropic +in 2D plane so as to not change the entire model and is left for the future works. The values of the +moduli of elasticity are found using The unidirectional ply properties are used to find the stiffness matrix of the plies: + +.. math:: + + Q = \begin{bmatrix} + \frac{E_1}{1-\nu_{12}\nu_{21}} & \frac{\nu_{21}E_2}{1-\nu_{12}\nu_{21}} & 0 \\ + \frac{\nu_{12}E_1}{1-\nu_{12}\nu_{21}} & \frac{E_2}{1-\nu_{12}\nu_{21}} & 0 \\ + 0 & 0 & G_{12} + \end{bmatrix} + +where :math:`E_1` and :math:`E_2` are the moduli of elasticity in the fiber direction and transverse direction, respectively, +:math:`\nu_{12}` and :math:`\nu_{21}` are the Poisson's ratios, and :math:`G_{12}` is the shear modulus. + +The transformed reduced stiffness matrix is found using the following equation: + +.. math:: + + \bar{Q} = T^T Q T + +where :math:`T` is the transformation matrix. The transformation matrix is found using the following equation: + +.. math:: + + T = \begin{bmatrix} + \cos \theta & \sin \theta & 0 \\ + -\sin \theta & \cos \theta & 0 \\ + 0 & 0 & 1 + \end{bmatrix} + +where :math:`\theta` is the angle of the ply with respect to the x-axis. + +The effective reduced stiffness matrix for the laminate is found using the weighted sum of the reduced stiffness matrices of the plies, +using their respective ply fraction constituition: + +.. math:: + + Q_{eff} = \sum_{i=1}^{n} f_i Q_{bar_i} + +where :math:`f_i` is the ply fraction of the :math:`i^{th}` ply. + +The effective compliance matrix is found using the following equation: + +.. math:: + + S_{eff} = Q_{eff}^{-1} + +The effective laminate properties are found using the following equations: + +.. math:: + E_{11} = \frac{1}{S_{eff_{11}}}\\ + G = \frac{1}{S_{eff_{66}}} + +These moduli of elasticility values are hence used to determine the stiffness matrix of the entire FEM spatial beam model. Thereafter, at the 4 critical points in the wingbox (mentioned in the aerostruct-wingbox walkthrough), +the strains are calculated for each of the constituent plies by transforming the strains at the critical points to the laminate coordinate system. This is done using the following equation:\\ +using the transformation: + +.. math:: + + \begin{pmatrix} + \epsilon_1 \\ + \epsilon_2 \\ + \gamma_{12} + \end{pmatrix} + = + [T] + \begin{pmatrix} + \epsilon_x \\ + \epsilon_y \\ + \gamma_{xy} + \end{pmatrix} + +The strains are then used to calculate the stresses in the laminate using the following equation: + +.. math:: + + \begin{pmatrix} + \sigma_1 \\ + \sigma_2 \\ + \tau_{12} + \end{pmatrix} + = + [Q] + \begin{pmatrix} + \epsilon_1 \\ + \epsilon_2 \\ + \gamma_{12} + \end{pmatrix} + +These local axial and shear stresses are then utilized to calculate the value of the **Strength Ratios** for each ply using Equation 5, where the coefficients are defined by: + +.. math:: + + F_{11} = \frac{1}{S_L^{(+)} S_L^{(-)}} \quad \text{and} \quad F_1 = \frac{1}{S_L^{(+)}} - \frac{1}{S_L^{(-)}} + +.. math:: + + F_{22} = \frac{1}{S_T^{(+)} S_T^{(-)}} \quad \text{and} \quad F_2 = \frac{1}{S_T^{(+)}} - \frac{1}{S_T^{(-)}} + +.. math:: + + F_{66} = \frac{1}{2 S_{LT}^{2}} + +where :math:`S_L^{(+)} \text{and} S_L^{(-)}` are the longitudinal strengths in tension and compression respectively, +:math:`S_T^{(+)} \text{and} S_T^{(-)}` are the transverse strengths in tension and compression respectively and +:math:`S_{LT}^{(+)}` is the shear strength of a ply. The strength ratios are then used to calculate the Tsai-Wu based failure criteria for each ply. +The Tsai-Wu failure criteria is given by: + +.. math:: + + F_1 \sigma_1 + F_2 \sigma_2 + F_{11} \sigma_1^2 + F_{22} \sigma_2^2 + F_{66} \tau_{12}^2 = 1 + +In order to implement the safety factor in the Tsai-Wu failure criteria, the equation is re-written as: + +.. math:: + a &= F_1 \sigma_1 + F_2 \sigma_2 \\ + b &= F_{11} \sigma_1^2 + F_{22} \sigma_2^2 + F_{12} \sigma_1 \sigma_2 + +We hence caclulate the **Strength Ratios** using the formula: + +.. math:: + + SR = \frac{1}{2} (a + \sqrt{a^2 + 4 b}) + +The strength ratio values hence calculated for each ply (determined by the length of ``ply_angles``) at each critical point (4 total), +(hence 4 x ``numplies`` strength ratio values for each beam element) for all beam elements are aggregated using a **KS Aggregate** function: + +.. math:: + + \hat{g}_{KS}(\rho, g) = \max_j g_j + \frac{1}{\rho} \ln \left( \sum_{j=1}^{n_g} \exp \left( \rho (g_j - \max_j g_j) \right) \right) + + +where :math:`g` is :math:`\left( \frac{SR}{SR_{\text{lim}}} - 1 \right)` value for each ply and :math:`SR_{\text{lim}}` is defined as: + +.. math:: + + SR_{\text{lim}} = \frac{1}{FOS} + + +The failure is determined by the value of :math:`\hat{g}_{KS}(\rho, g)` exceeding 0. + +The effect of using composites can be seen in the following figure. A Pareto-optimal front is generated for the wingbox model using Isotropic (Almunim) and Orthotropic (Carbon Fiber Reinforced Polymer) materials. + +.. image:: /advanced_features/figs/compositeModelPareto.png + :width: 600 + :align: center diff --git a/openaerostruct/docs/advanced_features/figs/compositeModelPareto.png b/openaerostruct/docs/advanced_features/figs/compositeModelPareto.png new file mode 100644 index 000000000..8e753beac Binary files /dev/null and b/openaerostruct/docs/advanced_features/figs/compositeModelPareto.png differ diff --git a/openaerostruct/docs/advanced_features/scripts/run_vsp_777.py b/openaerostruct/docs/advanced_features/scripts/run_vsp_777.py index bc92533ba..ba3f93daf 100644 --- a/openaerostruct/docs/advanced_features/scripts/run_vsp_777.py +++ b/openaerostruct/docs/advanced_features/scripts/run_vsp_777.py @@ -4,6 +4,7 @@ 777 vsp model avaialable here: http://hangar.openvsp.org/vspfiles/375 """ + import os import numpy as np diff --git a/openaerostruct/docs/advanced_features/scripts/two_part_wing_custom_mesh.py b/openaerostruct/docs/advanced_features/scripts/two_part_wing_custom_mesh.py index 81ef249e4..eb50882b6 100644 --- a/openaerostruct/docs/advanced_features/scripts/two_part_wing_custom_mesh.py +++ b/openaerostruct/docs/advanced_features/scripts/two_part_wing_custom_mesh.py @@ -160,7 +160,8 @@ # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin "wing_weight_ratio": 1.25, diff --git a/openaerostruct/docs/advanced_features/scripts/wingbox_mpt_Q400_example.py b/openaerostruct/docs/advanced_features/scripts/wingbox_mpt_Q400_example.py index 77aa0911f..5126c87c8 100644 --- a/openaerostruct/docs/advanced_features/scripts/wingbox_mpt_Q400_example.py +++ b/openaerostruct/docs/advanced_features/scripts/wingbox_mpt_Q400_example.py @@ -93,7 +93,8 @@ # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin "wing_weight_ratio": 1.25, diff --git a/openaerostruct/docs/aerostructural_wingbox_walkthrough.rst b/openaerostruct/docs/aerostructural_wingbox_walkthrough.rst index 878eda108..cb7f8c246 100644 --- a/openaerostruct/docs/aerostructural_wingbox_walkthrough.rst +++ b/openaerostruct/docs/aerostructural_wingbox_walkthrough.rst @@ -100,10 +100,10 @@ The `c_max_t` value is the location of the maximum thickness of the airfoil alon :end-before: checkpoint 6 Next we provide some information related to structures. -We provide the Young's and shear moduli, as well as the allowable yield stress and density of the material being used (the wingbox model currently assumes that the material is isotropic). -Here we use include a safety factor of 1.5 in the allowable yield stress. +We provide the Young's and shear moduli, as well as the maximum yield stress and density of the material being used. The material can either be assumed to be isotropic (like aluminum) or orthotropic (like carbon fiber reinforced polymer). +Here we use include a safety factor of 1.5 using the `safety_factor` key. The `strength_factor_for_upper_skin` value can be used to adjust the yield strength of the upper skin relative to the lower skin. -For example, if we were using different alloys for the top and bottom skin (e.g., Al7075 vs Al2024) we would provide the allowable yield stress of the lower skin for `yield` and the ratio of the yield strengths of the upper skin to the lower skin for `strength_factor_for_upper_skin`. +For example, if we were using different alloys for the top and bottom skin (e.g., Al7075 vs Al2024) we would provide the maximum yield stress of the lower skin for `yield` and the ratio of the yield strengths of the upper skin to the lower skin for `strength_factor_for_upper_skin`. The `wing_weight_ratio` number is used to estimate the weight of other components not modeled in the wingbox structure (e.g., overlaps, fasteners, etc.). With the `exact_failure_constraint` set to `False`, we aggregate the stress constraints for the FEM elements into a single constraint using the Kreisselmeier--Steinhauser function. This helps reduce computational cost during optimization by replacing a large number of constraints (one for each stress combination for each element) with a single constraint. diff --git a/openaerostruct/docs/composite_wingbox_mpt_opt_example.py b/openaerostruct/docs/composite_wingbox_mpt_opt_example.py new file mode 100644 index 000000000..091f9e1bd --- /dev/null +++ b/openaerostruct/docs/composite_wingbox_mpt_opt_example.py @@ -0,0 +1,364 @@ +# Ignore the #docs checkpoint comments. They are just used to split up the code for the documentation webpage. +# docs checkpoint 0 + +import numpy as np +from openaerostruct.geometry.utils import generate_mesh +from openaerostruct.integration.aerostruct_groups import AerostructGeometry, AerostructPoint +from openaerostruct.structures.wingbox_fuel_vol_delta import WingboxFuelVolDelta +import openmdao.api as om + +# docs checkpoint 1 + +# Provide coordinates for a portion of an airfoil for the wingbox cross-section as an nparray with dtype=complex (to work with the complex-step approximation for derivatives). +# These should be for an airfoil with the chord scaled to 1. +# We use the 10% to 60% portion of the NASA SC2-0612 airfoil for this case +# We use the coordinates available from airfoiltools.com. Using such a large number of coordinates is not necessary. +# The first and last x-coordinates of the upper and lower surfaces must be the same + +# fmt: off +upper_x = np.array([0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6], dtype="complex128") +lower_x = np.array([0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6], dtype="complex128") +upper_y = np.array([ 0.0447, 0.046, 0.0472, 0.0484, 0.0495, 0.0505, 0.0514, 0.0523, 0.0531, 0.0538, 0.0545, 0.0551, 0.0557, 0.0563, 0.0568, 0.0573, 0.0577, 0.0581, 0.0585, 0.0588, 0.0591, 0.0593, 0.0595, 0.0597, 0.0599, 0.06, 0.0601, 0.0602, 0.0602, 0.0602, 0.0602, 0.0602, 0.0601, 0.06, 0.0599, 0.0598, 0.0596, 0.0594, 0.0592, 0.0589, 0.0586, 0.0583, 0.058, 0.0576, 0.0572, 0.0568, 0.0563, 0.0558, 0.0553, 0.0547, 0.0541], dtype="complex128") # noqa: E201, E241 +lower_y = np.array([-0.0447, -0.046, -0.0473, -0.0485, -0.0496, -0.0506, -0.0515, -0.0524, -0.0532, -0.054, -0.0547, -0.0554, -0.056, -0.0565, -0.057, -0.0575, -0.0579, -0.0583, -0.0586, -0.0589, -0.0592, -0.0594, -0.0595, -0.0596, -0.0597, -0.0598, -0.0598, -0.0598, -0.0598, -0.0597, -0.0596, -0.0594, -0.0592, -0.0589, -0.0586, -0.0582, -0.0578, -0.0573, -0.0567, -0.0561, -0.0554, -0.0546, -0.0538, -0.0529, -0.0519, -0.0509, -0.0497, -0.0485, -0.0472, -0.0458, -0.0444], dtype="complex128") +# fmt: on + +# docs checkpoint 2 + +# Create a dictionary to store options about the surface +mesh_dict = { + "num_y": 15, + "num_x": 3, + "wing_type": "uCRM_based", + "symmetry": True, + "chord_cos_spacing": 0, + "span_cos_spacing": 0, + "num_twist_cp": 4, +} + +mesh, twist_cp = generate_mesh(mesh_dict) + +# docs checkpoint 3 + +surf_dict = { + # Wing definition + "name": "wing", # give the surface some name + "symmetry": True, # if True, model only one half of the lifting surface + "S_ref_type": "projected", # how we compute the wing area, + # can be 'wetted' or 'projected' + "mesh": mesh, + "fem_model_type": "wingbox", # 'wingbox' or 'tube' + "data_x_upper": upper_x, + "data_x_lower": lower_x, + "data_y_upper": upper_y, + "data_y_lower": lower_y, + # docs checkpoint 4 + "twist_cp": np.array([4.0, 5.0, 8.0, 9.0]), # [deg] + "spar_thickness_cp": np.array([0.004, 0.005, 0.008, 0.01]), # [m] + "skin_thickness_cp": np.array([0.005, 0.01, 0.015, 0.025]), # [m] + "t_over_c_cp": np.array([0.08, 0.08, 0.10, 0.08]), + "original_wingbox_airfoil_t_over_c": 0.12, + # docs checkpoint 5 + # Aerodynamic deltas. + # These CL0 and CD0 values are added to the CL and CD + # obtained from aerodynamic analysis of the surface to get + # the total CL and CD. + # These CL0 and CD0 values do not vary wrt alpha. + # They can be used to account for things that are not included, such as contributions from the fuselage, camber, etc. + "CL0": 0.0, # CL delta + "CD0": 0.0078, # CD delta + "with_viscous": True, # if true, compute viscous drag + "with_wave": True, # if true, compute wave drag + # Airfoil properties for viscous drag calculation + "k_lam": 0.05, # fraction of chord with laminar + # flow, used for viscous drag + "c_max_t": 0.38, # chordwise location of maximum thickness + # docs checkpoint 6 + # Structural values are based on aluminum 7075 + "E": 73.1e9, # [Pa] Young's modulus + "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) + "yield": 420.0e6, # [Pa] yield stress + "mrho": 2.78e3, # [kg/m^3] material density + "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin + "wing_weight_ratio": 1.25, + "exact_failure_constraint": False, # if false, use KS function + "struct_weight_relief": True, + "distributed_fuel_weight": True, + "n_point_masses": 1, # number of point masses in the system; in this case, the engine (omit option if no point masses) + "fuel_density": 803.0, # [kg/m^3] fuel density (only needed if the fuel-in-wing volume constraint is used) + "Wf_reserve": 15000.0, # [kg] reserve fuel mass + # docs checkpoint 7 + "useComposite": True, + "safety_factor": 1.5, + "ply_angles": [0, 45, -45, 90], + "ply_fractions": [0.10, 0.25, 0.25, 0.40], + "E1": 117.7e9, + "E2": 9.7e9, + "nu12": 0.35, + "G12": 4.8e9, + "sigma_t1": 1648.0e6, + "sigma_c1": 1034.0e6, + "sigma_t2": 64.0e6, + "sigma_c2": 228.0e6, + "sigma_12max": 71.0e6, + # docs checkpoint 8 +} + +surfaces = [surf_dict] + +# docs checkpoint 9 + +# Create the problem and assign the model group +prob = om.Problem() + +# Add problem information as an independent variables component +indep_var_comp = om.IndepVarComp() +indep_var_comp.add_output("Mach_number", val=np.array([0.85, 0.64])) +indep_var_comp.add_output("v", val=np.array([0.85 * 295.07, 0.64 * 340.294]), units="m/s") +indep_var_comp.add_output( + "re", + val=np.array([0.348 * 295.07 * 0.85 * 1.0 / (1.43 * 1e-5), 1.225 * 340.294 * 0.64 * 1.0 / (1.81206 * 1e-5)]), + units="1/m", +) +indep_var_comp.add_output("rho", val=np.array([0.348, 1.225]), units="kg/m**3") +indep_var_comp.add_output("speed_of_sound", val=np.array([295.07, 340.294]), units="m/s") + +# docs checkpoint 10 + +indep_var_comp.add_output("CT", val=0.53 / 3600, units="1/s") +indep_var_comp.add_output("R", val=14.307e6, units="m") +indep_var_comp.add_output("W0_without_point_masses", val=128000 + surf_dict["Wf_reserve"], units="kg") + +# docs checkpoint 11 + +indep_var_comp.add_output("load_factor", val=np.array([1.0, 2.5])) +indep_var_comp.add_output("alpha", val=0.0, units="deg") +indep_var_comp.add_output("alpha_maneuver", val=0.0, units="deg") + +indep_var_comp.add_output("empty_cg", val=np.zeros((3)), units="m") + +# docs checkpoint 12 + +indep_var_comp.add_output("fuel_mass", val=10000.0, units="kg") + +prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + +# docs checkpoint 12.5 + +point_masses = np.array([[10.0e3]]) + +point_mass_locations = np.array([[25, -10.0, 0.0]]) + +indep_var_comp.add_output("point_masses", val=point_masses, units="kg") +indep_var_comp.add_output("point_mass_locations", val=point_mass_locations, units="m") + +# Compute the actual W0 to be used within OAS based on the sum of the point mass and other W0 weight +prob.model.add_subsystem( + "W0_comp", om.ExecComp("W0 = W0_without_point_masses + 2 * sum(point_masses)", units="kg"), promotes=["*"] +) + +# docs checkpoint 13 + +# Loop over each surface in the surfaces list +for surface in surfaces: + # Get the surface name and create a group to contain components + # only for this surface + name = surface["name"] + + aerostruct_group = AerostructGeometry(surface=surface) + + # Add groups to the problem with the name of the surface. + prob.model.add_subsystem(name, aerostruct_group) + +# docs checkpoint 14 + +# Loop through and add a certain number of aerostruct points +for i in range(2): + point_name = "AS_point_{}".format(i) + # Connect the parameters within the model for each aero point + + # Create the aerostruct point group and add it to the model + AS_point = AerostructPoint(surfaces=surfaces, internally_connect_fuelburn=False) + + prob.model.add_subsystem(point_name, AS_point) + + # docs checkpoint 15 + + # Connect flow properties to the analysis point + prob.model.connect("v", point_name + ".v", src_indices=[i]) + prob.model.connect("Mach_number", point_name + ".Mach_number", src_indices=[i]) + prob.model.connect("re", point_name + ".re", src_indices=[i]) + prob.model.connect("rho", point_name + ".rho", src_indices=[i]) + prob.model.connect("CT", point_name + ".CT") + prob.model.connect("R", point_name + ".R") + prob.model.connect("W0", point_name + ".W0") + prob.model.connect("speed_of_sound", point_name + ".speed_of_sound", src_indices=[i]) + prob.model.connect("empty_cg", point_name + ".empty_cg") + prob.model.connect("load_factor", point_name + ".load_factor", src_indices=[i]) + prob.model.connect("fuel_mass", point_name + ".total_perf.L_equals_W.fuelburn") + prob.model.connect("fuel_mass", point_name + ".total_perf.CG.fuelburn") + + # docs checkpoint 16 + + for surface in surfaces: + name = surface["name"] + + if surf_dict["distributed_fuel_weight"]: + prob.model.connect("load_factor", point_name + ".coupled.load_factor", src_indices=[i]) + + com_name = point_name + "." + name + "_perf." + prob.model.connect( + name + ".local_stiff_transformed", point_name + ".coupled." + name + ".local_stiff_transformed" + ) + prob.model.connect(name + ".nodes", point_name + ".coupled." + name + ".nodes") + + # Connect aerodyamic mesh to coupled group mesh + prob.model.connect(name + ".mesh", point_name + ".coupled." + name + ".mesh") + if surf_dict["struct_weight_relief"]: + prob.model.connect(name + ".element_mass", point_name + ".coupled." + name + ".element_mass") + + # Connect performance calculation variables + prob.model.connect(name + ".nodes", com_name + "nodes") + prob.model.connect(name + ".cg_location", point_name + "." + "total_perf." + name + "_cg_location") + prob.model.connect(name + ".structural_mass", point_name + "." + "total_perf." + name + "_structural_mass") + + # Connect wingbox properties to von Mises stress calcs + prob.model.connect(name + ".Qz", com_name + "Qz") + prob.model.connect(name + ".J", com_name + "J") + prob.model.connect(name + ".A_enc", com_name + "A_enc") + prob.model.connect(name + ".htop", com_name + "htop") + prob.model.connect(name + ".hbottom", com_name + "hbottom") + prob.model.connect(name + ".hfront", com_name + "hfront") + prob.model.connect(name + ".hrear", com_name + "hrear") + + prob.model.connect(name + ".spar_thickness", com_name + "spar_thickness") + prob.model.connect(name + ".t_over_c", com_name + "t_over_c") + + coupled_name = point_name + ".coupled." + name + prob.model.connect("point_masses", coupled_name + ".point_masses") + prob.model.connect("point_mass_locations", coupled_name + ".point_mass_locations") + + +# docs checkpoint 17 + +prob.model.connect("alpha", "AS_point_0" + ".alpha") +prob.model.connect("alpha_maneuver", "AS_point_1" + ".alpha") + +# docs checkpoint 18 + +# Here we add the fuel volume constraint componenet to the model +prob.model.add_subsystem("fuel_vol_delta", WingboxFuelVolDelta(surface=surface)) +prob.model.connect("wing.struct_setup.fuel_vols", "fuel_vol_delta.fuel_vols") +prob.model.connect("AS_point_0.fuelburn", "fuel_vol_delta.fuelburn") + +if surf_dict["distributed_fuel_weight"]: + prob.model.connect("wing.struct_setup.fuel_vols", "AS_point_0.coupled.wing.struct_states.fuel_vols") + prob.model.connect("fuel_mass", "AS_point_0.coupled.wing.struct_states.fuel_mass") + + prob.model.connect("wing.struct_setup.fuel_vols", "AS_point_1.coupled.wing.struct_states.fuel_vols") + prob.model.connect("fuel_mass", "AS_point_1.coupled.wing.struct_states.fuel_mass") + +# docs checkpoint 19 + +comp = om.ExecComp("fuel_diff = (fuel_mass - fuelburn) / fuelburn", units="kg") +prob.model.add_subsystem("fuel_diff", comp, promotes_inputs=["fuel_mass"], promotes_outputs=["fuel_diff"]) +prob.model.connect("AS_point_0.fuelburn", "fuel_diff.fuelburn") + +# docs checkpoint 20 + +prob.model.add_objective("AS_point_0.fuelburn", scaler=1e-5) + +prob.model.add_design_var("wing.twist_cp", lower=-15.0, upper=15.0, scaler=0.1) +prob.model.add_design_var("wing.spar_thickness_cp", lower=0.003, upper=0.1, scaler=1e2) +prob.model.add_design_var("wing.skin_thickness_cp", lower=0.003, upper=0.1, scaler=1e2) +prob.model.add_design_var("wing.geometry.t_over_c_cp", lower=0.07, upper=0.2, scaler=10.0) +prob.model.add_design_var("alpha_maneuver", lower=-15.0, upper=15) + +# docs checkpoint 21 + +prob.model.add_constraint("AS_point_0.CL", equals=0.5) + +# docs checkpoint 22 + +prob.model.add_constraint("AS_point_1.L_equals_W", equals=0.0) +prob.model.add_constraint("AS_point_1.wing_perf.failure", upper=0.0) + +# docs checkpoint 23 + +prob.model.add_constraint("fuel_vol_delta.fuel_vol_delta", lower=0.0) + +# docs checkpoint 24 + +prob.model.add_design_var("fuel_mass", lower=0.0, upper=2e5, scaler=1e-5) +prob.model.add_constraint("fuel_diff", equals=0.0) + +# docs checkpoint 25 + +prob.driver = om.ScipyOptimizeDriver() +prob.driver.options["optimizer"] = "SLSQP" +prob.driver.options["tol"] = 1e-2 + +# docs checkpoint 26 + +recorder = om.SqliteRecorder("aerostruct.db") +prob.driver.add_recorder(recorder) + +# We could also just use prob.driver.recording_options['includes']=['*'] here, but for large meshes the database file becomes extremely large. So we just select the variables we need. +prob.driver.recording_options["includes"] = [ + "alpha", + "rho", + "v", + "cg", + "AS_point_1.cg", + "AS_point_0.cg", + "AS_point_0.coupled.wing_loads.loads", + "AS_point_1.coupled.wing_loads.loads", + "AS_point_0.coupled.wing.normals", + "AS_point_1.coupled.wing.normals", + "AS_point_0.coupled.wing.widths", + "AS_point_1.coupled.wing.widths", + "AS_point_0.coupled.aero_states.wing_sec_forces", + "AS_point_1.coupled.aero_states.wing_sec_forces", + "AS_point_0.wing_perf.CL1", + "AS_point_1.wing_perf.CL1", + "AS_point_0.coupled.wing.S_ref", + "AS_point_1.coupled.wing.S_ref", + "wing.geometry.twist", + "wing.mesh", + "wing.skin_thickness", + "wing.spar_thickness", + "wing.t_over_c", + "wing.structural_mass", + "AS_point_0.wing_perf.vonmises", + "AS_point_1.wing_perf.vonmises", + "AS_point_0.coupled.wing.def_mesh", + "AS_point_1.coupled.wing.def_mesh", +] + +prob.driver.recording_options["record_objectives"] = True +prob.driver.recording_options["record_constraints"] = True +prob.driver.recording_options["record_desvars"] = True +prob.driver.recording_options["record_inputs"] = True + +# docs checkpoint 27 + +# Set up the problem +prob.setup() + +# change linear solver for aerostructural coupled adjoint +prob.model.AS_point_0.coupled.linear_solver = om.LinearBlockGS(iprint=0, maxiter=30, use_aitken=True) +prob.model.AS_point_1.coupled.linear_solver = om.LinearBlockGS(iprint=0, maxiter=30, use_aitken=True) + +# om.view_model(prob) + +# prob.check_partials(form='central', compact_print=True) + +prob.run_driver() + +print("The fuel burn value is", prob["AS_point_0.fuelburn"][0], "[kg]") +print( + "The wingbox mass (excluding the wing_weight_ratio) is", + prob["wing.structural_mass"][0] / surf_dict["wing_weight_ratio"], + "[kg]", +) + +# docs checkpoint 28 diff --git a/openaerostruct/docs/user_reference/mesh_surface_dict.rst b/openaerostruct/docs/user_reference/mesh_surface_dict.rst index 9f53ccbcd..b9463feb0 100644 --- a/openaerostruct/docs/user_reference/mesh_surface_dict.rst +++ b/openaerostruct/docs/user_reference/mesh_surface_dict.rst @@ -190,9 +190,13 @@ The surface dict will be provided to Groups, including ``Geometry``, ``AeroPoint - Pa - Shear modulus * - yield - - 420.0e6 / 1.5 + - 420.0e6 - Pa - - Allowable yield stress including the safety factor. + - Maximum yield stress of the material. + * - safety_factor + - 1.5 + - + - Factor of safety for the material. * - mrho - 2.78e3 - kg/m^3 diff --git a/openaerostruct/docs/wingbox_mpt_opt_example.py b/openaerostruct/docs/wingbox_mpt_opt_example.py index 6eab896b4..32b5b68e3 100644 --- a/openaerostruct/docs/wingbox_mpt_opt_example.py +++ b/openaerostruct/docs/wingbox_mpt_opt_example.py @@ -76,7 +76,8 @@ # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin "wing_weight_ratio": 1.25, diff --git a/openaerostruct/examples/black_box_example.py b/openaerostruct/examples/black_box_example.py index fcbbf95c5..aa8d88bd8 100644 --- a/openaerostruct/examples/black_box_example.py +++ b/openaerostruct/examples/black_box_example.py @@ -52,7 +52,8 @@ # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/openaerostruct/examples/rectangular_wing/run_rect_wing.py b/openaerostruct/examples/rectangular_wing/run_rect_wing.py index 906f07de6..7cb562314 100644 --- a/openaerostruct/examples/rectangular_wing/run_rect_wing.py +++ b/openaerostruct/examples/rectangular_wing/run_rect_wing.py @@ -3,6 +3,7 @@ Print out lift and drag coefficient when complete. Check output directory for Tecplot solution files. """ + import numpy as np import openmdao.api as om diff --git a/openaerostruct/examples/run_CRM.py b/openaerostruct/examples/run_CRM.py index ff7fcd47d..f42dc9782 100644 --- a/openaerostruct/examples/run_CRM.py +++ b/openaerostruct/examples/run_CRM.py @@ -42,7 +42,8 @@ # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/openaerostruct/examples/run_aerostruct_uCRM_multipoint.py b/openaerostruct/examples/run_aerostruct_uCRM_multipoint.py index 23195f293..02e4fa298 100644 --- a/openaerostruct/examples/run_aerostruct_uCRM_multipoint.py +++ b/openaerostruct/examples/run_aerostruct_uCRM_multipoint.py @@ -19,7 +19,6 @@ older version of OAS (very slight differences due to numerical errors, etc.) """ - import numpy as np from openaerostruct.geometry.utils import generate_mesh @@ -88,7 +87,8 @@ # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin "wing_weight_ratio": 1.25, diff --git a/openaerostruct/examples/run_scaneagle.py b/openaerostruct/examples/run_scaneagle.py index acb23276d..a4f4e7c81 100644 --- a/openaerostruct/examples/run_scaneagle.py +++ b/openaerostruct/examples/run_scaneagle.py @@ -97,6 +97,7 @@ "E": 85.0e9, "G": 25.0e9, "yield": 350.0e6, + "safety_factor": 1, "mrho": 1.6e3, "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 1.0, # multiplicative factor on the computed structural weight diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index 7eb1278be..3cc820f2b 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -12,10 +12,78 @@ from openaerostruct.structures.tube_group import TubeGroup from openaerostruct.structures.wingbox_group import WingboxGroup from openaerostruct.utils.check_surface_dict import check_surface_dict_keys - +import numpy as np import openmdao.api as om +def TransformationMatrix(theta): + """ + function to find the transformation matrix for a given angle + input: theta (angle in degrees) + output: transformation (T) matrix + """ + theta = theta * np.pi / 180 + c = np.cos(theta) + s = np.sin(theta) + T = np.zeros((3, 3)) + T[0, 0] = c**2 + T[0, 1] = s**2 + T[0, 2] = 2 * s * c + T[1, 0] = s**2 + T[1, 1] = c**2 + T[1, 2] = -2 * s * c + T[2, 0] = -s * c + T[2, 1] = s * c + T[2, 2] = c**2 - s**2 + return T + + +def computeCompositeStiffness(surface): + """ + Function to compute the effective E and G stiffness values for a composite material, + based on the ply_fractions, ply angles and individual fiber and matrix properties. + """ + E1 = surface["E1"] + E2 = surface["E2"] + v12 = surface["nu12"] + G12 = surface["G12"] + v21 = (E2 / E1) * v12 + ply_fractions = surface["ply_fractions"] + ply_angles = surface["ply_angles"] + num_plies = len(ply_fractions) + + # finding the Q matrix + Q = np.zeros((3, 3)) + Q[0, 0] = E1 / (1 - v12 * v21) + Q[0, 1] = v12 * E2 / (1 - v12 * v21) + Q[0, 2] = 0 + Q[1, 0] = v21 * E1 / (1 - v12 * v21) + Q[1, 1] = E2 / (1 - v12 * v21) + Q[1, 2] = 0 + Q[2, 0] = 0 + Q[2, 1] = 0 + Q[2, 2] = G12 + + # finding the Q_bar matrix for each ply in the form of a 3D Array + Q_bar = np.zeros((num_plies, 3, 3)) + Q_bar_eff = np.zeros((3, 3)) + for i in range(num_plies): + theta = ply_angles[i] + T = TransformationMatrix(theta) + Q_bar[i] = np.dot(np.dot(np.linalg.inv(T), Q), T) + Q_bar_eff += ply_fractions[i] * Q_bar[i] + + S_bar_eff = np.linalg.inv(Q_bar_eff) + E_eff = 1 / S_bar_eff[0, 0] + G_eff = 1 / S_bar_eff[2, 2] + + # replacing the values in the surface dictionary + surface["E"] = E_eff + surface["G"] = G_eff + + # no need to return anything as the values are updated in the surface dictionary (call by reference) + + class AerostructGeometry(om.Group): def initialize(self): self.options.declare("surface", types=dict) @@ -56,7 +124,7 @@ def setup(self): promotes_outputs=geom_promotes_out, ) - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": tube_promotes_input = [] tube_promotes_output = ["A", "Iy", "Iz", "J", "radius", "thickness"] if "thickness_cp" in surface.keys() and connect_geom_DVs: @@ -70,9 +138,23 @@ def setup(self): promotes_inputs=tube_promotes_input, promotes_outputs=tube_promotes_output, ) - elif surface["fem_model_type"] == "wingbox": + elif ( + surface["fem_model_type"].lower() == "wingbox" + ): # connections and nomenclature remains the same for both isotropic and composite wingbox wingbox_promotes_in = ["mesh", "t_over_c"] - wingbox_promotes_out = ["A", "Iy", "Iz", "J", "Qz", "A_enc", "A_int", "htop", "hbottom", "hfront", "hrear"] + wingbox_promotes_out = [ + "A", + "Iy", + "Iz", + "J", + "Qz", + "A_enc", + "A_int", + "htop", + "hbottom", + "hfront", + "hrear", + ] if "skin_thickness_cp" in surface.keys() and "spar_thickness_cp" in surface.keys(): wingbox_promotes_in.append("skin_thickness_cp") wingbox_promotes_in.append("spar_thickness_cp") @@ -90,7 +172,7 @@ def setup(self): else: raise NameError("Please select a valid `fem_model_type` from either `tube` or `wingbox`.") - if surface["fem_model_type"] == "wingbox": + if surface["fem_model_type"].lower() == "wingbox": # same for both isotropic and composite wingbox promotes = ["A_int"] else: promotes = [] @@ -171,7 +253,7 @@ def setup(self): promotes_outputs=["CDv", "CDw", "L", "D", "CL1", "CDi", "CD", "CL", "Cl"], ) - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": self.add_subsystem( "struct_funcs", SpatialBeamFunctionals(surface=surface), @@ -179,7 +261,12 @@ def setup(self): promotes_outputs=["thickness_intersects", "vonmises", "failure"], ) - elif surface["fem_model_type"] == "wingbox": + elif surface["fem_model_type"].lower() == "wingbox": + if "useComposite" in surface.keys() and surface["useComposite"]: # using the Composite Wing Box + promotedoutput = "tsaiwu_sr" + else: # using the isotropic Wing Box + promotedoutput = "vonmises" + self.add_subsystem( "struct_funcs", SpatialBeamFunctionals(surface=surface), @@ -195,7 +282,7 @@ def setup(self): "nodes", "disp", ], - promotes_outputs=["vonmises", "failure"], + promotes_outputs=[promotedoutput, "failure"], ) else: raise NameError("Please select a valid `fem_model_type` from either `tube` or `wingbox`.") @@ -225,6 +312,11 @@ def setup(self): for surface in surfaces: name = surface["name"] + # if useComposite is enabled, compute the effective E and G values for the composite material + useComposite = "useComposite" in surface.keys() and surface["useComposite"] + if useComposite: + computeCompositeStiffness(surface) + # Connect the output of the loads component with the FEM # displacement parameter. This links the coupling within the coupled # group that necessitates the subgroup solver. diff --git a/openaerostruct/structures/compute_nodes.py b/openaerostruct/structures/compute_nodes.py index af4a44920..0e18f6e90 100644 --- a/openaerostruct/structures/compute_nodes.py +++ b/openaerostruct/structures/compute_nodes.py @@ -31,7 +31,7 @@ def setup(self): nx = mesh.shape[0] ny = mesh.shape[1] - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": self.fem_origin = surface["fem_origin"] else: y_upper = surface["data_y_upper"] diff --git a/openaerostruct/structures/failure_exact.py b/openaerostruct/structures/failure_exact.py index 199047b4d..b71a7ed1e 100644 --- a/openaerostruct/structures/failure_exact.py +++ b/openaerostruct/structures/failure_exact.py @@ -7,15 +7,19 @@ class FailureExact(om.ExplicitComponent): """ Output individual failure constraints on each FEM element. - Parameters + parameters ---------- - vonmises[ny-1, 2] : numpy array + vonmises : ny-1 x 2 numpy array von Mises stress magnitudes for each FEM element. + tsaiwu_sr : ny-1 x 4 * num_plies numpy array + Tsai-Wu strength ratios for each FEM element (ply at each critical element). Returns ------- - failure[ny-1, 2] : numpy array - Array of failure conditions. Positive if element has failed. + failure : ny-1 x num_failure_criteria array + Array of failure conditions. Positive if element has failed. This entity is defined for either failure criteria, + vonmises or tsaiwu_sr. num_failure_criteria is 2 for tube, 4 for the isotropic wingbox and 4*num_plies for the + composite wingbox. """ @@ -24,19 +28,48 @@ def initialize(self): def setup(self): surface = self.options["surface"] + self.useComposite = "useComposite" in self.options["surface"].keys() and self.options["surface"]["useComposite"] + if self.useComposite: + ply_angles = surface["ply_angles"] + num_plies = len(ply_angles) - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": num_failure_criteria = 2 - elif surface["fem_model_type"] == "wingbox": - num_failure_criteria = 4 + + elif surface["fem_model_type"].lower() == "wingbox": + if self.useComposite: # using the Composite wingbox + num_failure_criteria = 4 * num_plies # 4 critical elements * number of plies + self.srlimit = 1 / surface["safety_factor"] + else: # using the Isotropic wingbox + num_failure_criteria = 4 self.ny = surface["mesh"].shape[1] - self.sigma = surface["yield"] - self.add_input("vonmises", val=np.zeros((self.ny - 1, num_failure_criteria)), units="N/m**2") + if "safety_factor" in self.options["surface"].keys(): + self.safety_factor = surface["safety_factor"] + else: + self.safety_factor = 1 + + self.sigma = surface["yield"] / self.safety_factor + + if self.useComposite: # using the Composite wingbox + self.add_input("tsaiwu_sr", val=np.zeros((self.ny - 1, num_failure_criteria)), units=None) + else: # using the Isotropic structures + self.add_input("vonmises", val=np.zeros((self.ny - 1, num_failure_criteria)), units="N/m**2") + self.add_output("failure", val=np.zeros((self.ny - 1, num_failure_criteria))) - self.declare_partials("failure", "vonmises", val=np.eye(((self.ny - 1) * num_failure_criteria)) / self.sigma) + if self.useComposite: # using the Composite wingbox + self.declare_partials( + "failure", "tsaiwu_sr", val=np.eye(((self.ny - 1) * num_failure_criteria)) / self.srlimit + ) + else: # using the Isotropic structures + self.declare_partials( + "failure", "vonmises", val=np.eye(((self.ny - 1) * num_failure_criteria)) / self.sigma + ) def compute(self, inputs, outputs): - outputs["failure"] = inputs["vonmises"] / self.sigma - 1 + if "vonmises" in inputs: + outputs["failure"] = inputs["vonmises"] / self.sigma - 1 + elif "tsaiwu_sr" in inputs: + outputs["failure"] = inputs["tsaiwu_sr"] / self.srlimit - 1 diff --git a/openaerostruct/structures/failure_ks.py b/openaerostruct/structures/failure_ks.py index db7941b95..eda42ee69 100644 --- a/openaerostruct/structures/failure_ks.py +++ b/openaerostruct/structures/failure_ks.py @@ -21,16 +21,18 @@ class FailureKS(om.ExplicitComponent): parameters ---------- - vonmises[ny-1, 2] : numpy array + vonmises : ny-1 x 2 numpy array von Mises stress magnitudes for each FEM element. + tsaiwu_sr : ny-1 x 4 * num_plies numpy array + Tsai-Wu strength ratios for each FEM element (ply at each critical element). Returns ------- failure : float KS aggregation quantity obtained by combining the failure criteria for each FEM node. Used to simplify the optimization problem by - reducing the number of constraints. - + reducing the number of constraints. This entity is defined for either + failure criteria, vonmises or tsaiwu_sr. """ def initialize(self): @@ -39,55 +41,68 @@ def initialize(self): def setup(self): surface = self.options["surface"] - rho = self.options["rho"] - - if surface["fem_model_type"] == "tube": + self.rho = self.options["rho"] + + if "safety_factor" in self.options["surface"].keys(): + self.safety_factor = surface["safety_factor"] + else: + self.safety_factor = 1 + + self.useComposite = "useComposite" in self.options["surface"].keys() and self.options["surface"]["useComposite"] + if self.useComposite: + self.num_plies = len(surface["ply_angles"]) + self.input_name = "tsaiwu_sr" + self.stress_limit = 1 / self.safety_factor + self.stress_units = None + else: + self.input_name = "vonmises" + self.stress_limit = surface["yield"] / self.safety_factor + self.stress_units = "N/m**2" + + if surface["fem_model_type"].lower() == "tube": num_failure_criteria = 2 - elif surface["fem_model_type"] == "wingbox": - num_failure_criteria = 4 + + elif surface["fem_model_type"].lower() == "wingbox": + if self.useComposite: # using the Composite wingbox + num_failure_criteria = 4 * self.num_plies # 4 critical elements * number of plies + else: # using the Isotropic wingbox + num_failure_criteria = 4 self.ny = surface["mesh"].shape[1] - self.add_input("vonmises", val=np.zeros((self.ny - 1, num_failure_criteria)), units="N/m**2") - self.add_output("failure", val=0.0) + self.add_input(self.input_name, val=np.zeros((self.ny - 1, num_failure_criteria)), units=self.stress_units) - self.sigma = surface["yield"] - self.rho = rho + self.add_output("failure", val=0.0) self.declare_partials("*", "*") def compute(self, inputs, outputs): - sigma = self.sigma - rho = self.rho - vonmises = inputs["vonmises"] - fmax = np.max(vonmises / sigma - 1) + stress_array = inputs[self.input_name] + + fmax = np.max(stress_array / self.stress_limit - 1) nlog, nsum, nexp = np.log, np.sum, np.exp - ks = 1 / rho * nlog(nsum(nexp(rho * (vonmises / sigma - 1 - fmax)))) + ks = 1 / self.rho * nlog(nsum(nexp(self.rho * (stress_array / self.stress_limit - 1 - fmax)))) outputs["failure"] = fmax + ks def compute_partials(self, inputs, partials): - vonmises = inputs["vonmises"] - sigma = self.sigma - rho = self.rho + stress_array = inputs[self.input_name] - # Find the location of the max stress constraint - fmax = np.max(vonmises / sigma - 1) - i, j = np.where((vonmises / sigma - 1) == fmax) + fmax = np.max(stress_array / self.stress_limit - 1) + i, j = np.where((stress_array / self.stress_limit - 1) == fmax) i, j = i[0], j[0] - # Set incoming seed as 1 so we simply get the jacobian entries ksb = 1.0 - # Use results from the AD code to compute the jacobian entries - tempb0 = ksb / (rho * np.sum(np.exp(rho * (vonmises / sigma - fmax - 1)))) - tempb = np.exp(rho * (vonmises / sigma - fmax - 1)) * rho * tempb0 + tempb0 = ksb / (self.rho * np.sum(np.exp(self.rho * (stress_array / self.stress_limit - fmax - 1)))) + tempb = np.exp(self.rho * (stress_array / self.stress_limit - fmax - 1)) * self.rho * tempb0 fmaxb = ksb - np.sum(tempb) - # Populate the entries - derivs = tempb / sigma - derivs[i, j] += fmaxb / sigma + derivs = tempb / self.stress_limit + derivs[i, j] += fmaxb / self.stress_limit - # Reshape and save them to the jac dict - partials["failure", "vonmises"] = derivs.reshape(1, -1) + if self.useComposite: # using the Composite wingbox + partials["failure", "tsaiwu_sr"] = derivs.reshape(1, -1) + else: # using the Isotropic structures + partials["failure", "vonmises"] = derivs.reshape(1, -1) diff --git a/openaerostruct/structures/spatial_beam_functionals.py b/openaerostruct/structures/spatial_beam_functionals.py index cff8cdd76..4062e6f0e 100644 --- a/openaerostruct/structures/spatial_beam_functionals.py +++ b/openaerostruct/structures/spatial_beam_functionals.py @@ -5,6 +5,7 @@ # from openaerostruct.structures.spar_within_wing import SparWithinWing from openaerostruct.structures.vonmises_tube import VonMisesTube from openaerostruct.structures.vonmises_wingbox import VonMisesWingbox +from openaerostruct.structures.tsaiwu_wingbox import TsaiWuWingbox from openaerostruct.structures.non_intersecting_thickness import NonIntersectingThickness from openaerostruct.structures.failure_exact import FailureExact from openaerostruct.structures.failure_ks import FailureKS @@ -25,7 +26,7 @@ def setup(self): # Energy(surface=surface), # promotes=['*']) - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": self.add_subsystem( "thicknessconstraint", NonIntersectingThickness(surface=surface), @@ -39,24 +40,44 @@ def setup(self): promotes_inputs=["radius", "nodes", "disp"], promotes_outputs=["vonmises"], ) - elif surface["fem_model_type"] == "wingbox": - self.add_subsystem( - "vonmises", - VonMisesWingbox(surface=surface), - promotes_inputs=[ - "Qz", - "J", - "A_enc", - "spar_thickness", - "htop", - "hbottom", - "hfront", - "hrear", - "nodes", - "disp", - ], - promotes_outputs=["vonmises"], - ) + + elif surface["fem_model_type"].lower() == "wingbox": + if "useComposite" in surface.keys() and surface["useComposite"]: # using the Composite wingbox + self.add_subsystem( + "tsaiwu_sr", + TsaiWuWingbox(surface=surface), + promotes_inputs=[ + "Qz", + "J", + "A_enc", + "spar_thickness", + "htop", + "hbottom", + "hfront", + "hrear", + "nodes", + "disp", + ], + promotes_outputs=["tsaiwu_sr"], + ) + else: # using the Isotropic wingbox + self.add_subsystem( + "vonmises", + VonMisesWingbox(surface=surface), + promotes_inputs=[ + "Qz", + "J", + "A_enc", + "spar_thickness", + "htop", + "hbottom", + "hfront", + "hrear", + "nodes", + "disp", + ], + promotes_outputs=["vonmises"], + ) else: raise NameError("Please select a valid `fem_model_type` from either `tube` or `wingbox`.") @@ -67,10 +88,24 @@ def setup(self): # promotes=['*']) if surface["exact_failure_constraint"]: + if "useComposite" in surface.keys() and surface["useComposite"]: # using the Composite wingbox + promotedinput = "tsaiwu_sr" + else: # using the Isotropic structures + promotedinput = "vonmises" + self.add_subsystem( - "failure", FailureExact(surface=surface), promotes_inputs=["vonmises"], promotes_outputs=["failure"] + "failure", + FailureExact(surface=surface), + promotes_inputs=[promotedinput], + promotes_outputs=["failure"], ) + else: + if "useComposite" in surface.keys() and surface["useComposite"]: # using the Composite wingbox + promotedinput = "tsaiwu_sr" + else: # using the Isotropic structures + promotedinput = "vonmises" + self.add_subsystem( - "failure", FailureKS(surface=surface), promotes_inputs=["vonmises"], promotes_outputs=["failure"] + "failure", FailureKS(surface=surface), promotes_inputs=[promotedinput], promotes_outputs=["failure"] ) diff --git a/openaerostruct/structures/spatial_beam_setup.py b/openaerostruct/structures/spatial_beam_setup.py index 7c4b490d8..0960210f1 100644 --- a/openaerostruct/structures/spatial_beam_setup.py +++ b/openaerostruct/structures/spatial_beam_setup.py @@ -39,7 +39,7 @@ def setup(self): promotes_outputs=["cg_location"], ) - if surface["fem_model_type"] == "wingbox": + if "wingbox" in surface["fem_model_type"].lower(): self.add_subsystem( "fuel_vol", WingboxFuelVol(surface=surface), diff --git a/openaerostruct/structures/struct_groups.py b/openaerostruct/structures/struct_groups.py index 7dce65f3f..92ec8d0a0 100644 --- a/openaerostruct/structures/struct_groups.py +++ b/openaerostruct/structures/struct_groups.py @@ -20,7 +20,7 @@ def setup(self): "geometry", Geometry(surface=surface), promotes_inputs=[], promotes_outputs=["mesh", "t_over_c"] ) - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": tube_promotes_input = [] tube_promotes_output = ["A", "Iy", "Iz", "J", "radius", "thickness"] if "thickness_cp" in surface.keys(): @@ -34,9 +34,23 @@ def setup(self): promotes_inputs=tube_promotes_input, promotes_outputs=tube_promotes_output, ) - elif surface["fem_model_type"] == "wingbox": + elif ( + surface["fem_model_type"].lower() == "wingbox" + ): # connections and nomenclature remains same for both isotropic and composite wingbox wingbox_promotes_in = ["mesh", "t_over_c"] - wingbox_promotes_out = ["A", "Iy", "Iz", "J", "Qz", "A_enc", "A_int", "htop", "hbottom", "hfront", "hrear"] + wingbox_promotes_out = [ + "A", + "Iy", + "Iz", + "J", + "Qz", + "A_enc", + "A_int", + "htop", + "hbottom", + "hfront", + "hrear", + ] if "skin_thickness_cp" in surface.keys() and "spar_thickness_cp" in surface.keys(): wingbox_promotes_in.append("skin_thickness_cp") wingbox_promotes_in.append("spar_thickness_cp") @@ -54,14 +68,16 @@ def setup(self): else: raise NameError("Please select a valid `fem_model_type` from either `tube` or `wingbox`.") - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": self.add_subsystem( "struct_setup", SpatialBeamSetup(surface=surface), promotes_inputs=["mesh", "A", "Iy", "Iz", "J"], promotes_outputs=["nodes", "local_stiff_transformed", "structural_mass", "cg_location", "element_mass"], ) - else: + elif ( + surface["fem_model_type"].lower() == "wingbox" + ): # connections and nomenclature remains same for both isotropic and composite wingbox self.add_subsystem( "struct_setup", SpatialBeamSetup(surface=surface), @@ -92,14 +108,19 @@ def setup(self): promotes_outputs=["disp"], ) - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": self.add_subsystem( "struct_funcs", SpatialBeamFunctionals(surface=surface), promotes_inputs=["thickness", "radius", "nodes", "disp"], promotes_outputs=["thickness_intersects", "vonmises", "failure"], ) - else: + elif surface["fem_model_type"].lower() == "wingbox": + if "useComposite" in surface.keys() and surface["useComposite"]: # using the Composite wingbox + promotedoutput = "tsaiwu_sr" + else: # using the Isotropic wingbox + promotedoutput = "vonmises" + self.add_subsystem( "struct_funcs", SpatialBeamFunctionals(surface=surface), @@ -115,5 +136,5 @@ def setup(self): "hrear", "nodes", ], - promotes_outputs=["vonmises", "failure"], + promotes_outputs=[promotedoutput, "failure"], ) diff --git a/openaerostruct/structures/tsaiwu_wingbox.py b/openaerostruct/structures/tsaiwu_wingbox.py new file mode 100644 index 000000000..6a947f167 --- /dev/null +++ b/openaerostruct/structures/tsaiwu_wingbox.py @@ -0,0 +1,270 @@ +import numpy as np + +import openmdao.api as om + +from openaerostruct.structures.utils import norm, unit + + +class TsaiWuWingbox(om.ExplicitComponent): + """Compute the Tsai-Wu failure criteria for each element, and return an array of the Strength Ratios (SR) for each ply at each element (4). + The Tsai-Wu failure criterion is a quadratic expression that is used to predict the failure of composite materials. + + Parameters + ---------- + nodes[ny, 3] : numpy array + Flattened array with coordinates for each FEM node. + disp[ny, 6] : numpy array + Displacements of each FEM node. + Qz[ny-1] : numpy array + First moment of area above the neutral axis parallel to the local + z-axis (for each wingbox segment). + J[ny-1] : numpy array + Torsion constants for each wingbox segment. + A_enc[ny-1] : numpy array + Cross-sectional enclosed area (measured using the material midlines) of + each wingbox segment. + spar_thickness[ny-1] : numpy array + Material thicknesses of the front and rear spars 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). + + Returns + ------- + tsaiwu_sr[ny-1, 16] : numpy array + Tsai-Wu Strength Ratios for each FEM element (4 critical elements * number of plies). + + """ + + def initialize(self): + self.options.declare("surface", types=dict) + + def setup(self): + self.surface = surface = self.options["surface"] + self.ply_angles = surface["ply_angles"] + self.num_plies = len(self.ply_angles) + + self.ny = surface["mesh"].shape[1] + + self.add_input("nodes", val=np.zeros((self.ny, 3)), units="m") + self.add_input("disp", val=np.zeros((self.ny, 6)), units="m") + self.add_input("Qz", val=np.zeros((self.ny - 1)), units="m**3") + self.add_input("J", val=np.zeros((self.ny - 1)), units="m**4") + self.add_input("A_enc", val=np.zeros((self.ny - 1)), units="m**2") + self.add_input("spar_thickness", val=np.zeros((self.ny - 1)), units="m") + self.add_input("htop", val=np.zeros((self.ny - 1)), units="m") + self.add_input("hbottom", val=np.zeros((self.ny - 1)), units="m") + self.add_input("hfront", val=np.zeros((self.ny - 1)), units="m") + self.add_input("hrear", val=np.zeros((self.ny - 1)), units="m") + self.add_output("tsaiwu_sr", val=np.zeros((self.ny - 1, self.num_plies * 4))) + + self.E = surface["E"] + self.G = surface["G"] + + self.E1 = surface["E1"] + self.E2 = surface["E2"] + self.G12 = surface["G12"] + self.nu12 = surface["nu12"] + self.nu21 = self.nu12 * self.E2 / self.E1 + + self.sigma_c1 = surface["sigma_c1"] + self.sigma_t1 = surface["sigma_t1"] + self.sigma_c2 = surface["sigma_c2"] + self.sigma_t2 = surface["sigma_t2"] + self.sigma_12max = surface["sigma_12max"] + + self.tssf = surface["strength_factor_for_upper_skin"] + + self.declare_partials("*", "*", method="cs") + + def compute(self, inputs, outputs): + disp = inputs["disp"] + nodes = inputs["nodes"] + A_enc = inputs["A_enc"] + Qy = inputs["Qz"] + J = inputs["J"] + htop = inputs["htop"] + hbottom = inputs["hbottom"] + hfront = inputs["hfront"] + hrear = inputs["hrear"] + spar_thickness = inputs["spar_thickness"] + tsaiwu_sr = outputs["tsaiwu_sr"] + + # Only use complex type for these arrays if we're using cs to check derivs + dtype = type(disp[0, 0]) + T = np.zeros((3, 3), dtype=dtype) + x_gl = np.array([1, 0, 0], dtype=dtype) + + E1 = self.E1 + E2 = self.E2 + G12 = self.G12 + nu_12 = self.nu12 + nu_21 = self.nu21 + + num_elems = self.ny - 1 + for ielem in range(num_elems): + P0 = nodes[ielem, :] + P1 = nodes[ielem + 1, :] + L = norm(P1 - P0) + + x_loc = unit(P1 - P0) + y_loc = unit(np.cross(x_loc, x_gl)) + z_loc = unit(np.cross(x_loc, y_loc)) + + T[0, :] = x_loc + T[1, :] = y_loc + T[2, :] = z_loc + + u0x, u0y, u0z = T.dot(disp[ielem, :3]) + r0x, r0y, r0z = T.dot(disp[ielem, 3:]) + u1x, u1y, u1z = T.dot(disp[ielem + 1, :3]) + r1x, r1y, r1z = T.dot(disp[ielem + 1, 3:]) + + # ============================================================================== + # strain equations + # ============================================================================== + + # this is strain + axial_strain = (u1x - u0x) / L + + # this is torsion strain + torsion_shear_strain = J[ielem] / L * (r1x - r0x) / 2 / spar_thickness[ielem] / A_enc[ielem] + + # this is bending strain for the top skin + top_bending_strain = 1.0 / (L**2) * (6 * u0y + 2 * r0z * L - 6 * u1y + 4 * r1z * L) * htop[ielem] + + # this is bending strain for the bottom skin + bottom_bending_strain = -1.0 / (L**2) * (6 * u0y + 2 * r0z * L - 6 * u1y + 4 * r1z * L) * hbottom[ielem] + + # this is bending strain for the front spar + front_bending_strain = -1.0 / (L**2) * (-6 * u0z + 2 * r0y * L + 6 * u1z + 4 * r1y * L) * hfront[ielem] + + # this is bending strain for the rear spar + rear_bending_strain = 1.0 / (L**2) * (-6 * u0z + 2 * r0y * L + 6 * u1z + 4 * r1y * L) * hrear[ielem] + + # shear strain due to bending + vertical_shear_strain = ( + 1.0 + / (L**3) + * (-12 * u0y - 6 * r0z * L + 12 * u1y - 6 * r1z * L) + * Qy[ielem] + / (2 * spar_thickness[ielem]) + ) + + # The strain combinations for the 4 elements under consideration: (split into epsilonx, epsilony and gammatau) + # Defining the epsilon_elem array for the epsionx, epsiony and gammatau for each element + epsilon_elem = np.zeros((4, 3), dtype=dtype) + + # Element 1: + epsilon_elem[0, 0] = top_bending_strain + rear_bending_strain + axial_strain + epsilon_elem[0, 1] = 0 + epsilon_elem[0, 2] = torsion_shear_strain + + # Element 2: + epsilon_elem[1, 0] = bottom_bending_strain + front_bending_strain + axial_strain + epsilon_elem[1, 1] = 0 + epsilon_elem[1, 2] = torsion_shear_strain + + # Element 3: + epsilon_elem[2, 0] = front_bending_strain + axial_strain + epsilon_elem[2, 1] = 0 + epsilon_elem[2, 2] = torsion_shear_strain + vertical_shear_strain + + # Element 4: + epsilon_elem[3, 0] = rear_bending_strain + axial_strain + epsilon_elem[3, 1] = 0 + epsilon_elem[3, 2] = -torsion_shear_strain + vertical_shear_strain + + # defining the array for ply-orientation angles: + ply_angles = self.ply_angles + num_plies = len(ply_angles) + + # defining the epsilon_elem_ply array for the epsilon1, epsilon2 and gamma12 for each ply + epsilon_elem_ply = np.zeros((4, num_plies, 3), dtype=dtype) + sigma_elem_ply = np.zeros((4, num_plies, 3), dtype=dtype) + + # running a loop over the 4 elements and the number of plies to calculate the epsilon_elem_ply array + for elem_num in range(4): + for ply_num in range(num_plies): + epsilon_elem_ply[elem_num, ply_num, 0] = ( + epsilon_elem[elem_num, 0] * np.cos(np.radians(ply_angles[ply_num])) ** 2 + + epsilon_elem[elem_num, 1] * np.sin(np.radians(ply_angles[ply_num])) ** 2 + + 2 + * epsilon_elem[elem_num, 2] + * np.sin(np.radians(ply_angles[ply_num])) + * np.cos(np.radians(ply_angles[ply_num])) + ) + epsilon_elem_ply[elem_num, ply_num, 1] = ( + epsilon_elem[elem_num, 0] * np.sin(np.radians(ply_angles[ply_num])) ** 2 + + epsilon_elem[elem_num, 1] * np.cos(np.radians(ply_angles[ply_num])) ** 2 + - 2 + * epsilon_elem[elem_num, 2] + * np.sin(np.radians(ply_angles[ply_num])) + * np.cos(np.radians(ply_angles[ply_num])) + ) + epsilon_elem_ply[elem_num, ply_num, 2] = ( + -epsilon_elem[elem_num, 0] + * np.sin(np.radians(ply_angles[ply_num])) + * np.cos(np.radians(ply_angles[ply_num])) + + epsilon_elem[elem_num, 1] + * np.sin(np.radians(ply_angles[ply_num])) + * np.cos(np.radians(ply_angles[ply_num])) + + epsilon_elem[elem_num, 2] + * (np.cos(np.radians(ply_angles[ply_num])) ** 2 - np.sin(np.radians(ply_angles[ply_num])) ** 2) + ) + + # defining the Q matrix for the material: + Q11 = E1 / (1 - nu_12 * nu_21) + Q22 = E2 / (1 - nu_12 * nu_21) + Q12 = nu_12 * E2 / (1 - nu_12 * nu_21) + Q21 = Q12 + Q66 = G12 + + # converting the strains to stresses using strain-stress relations + for elem_num in range(4): + for ply_num in range(num_plies): + epsilon1 = epsilon_elem_ply[elem_num, ply_num, 0] + epsilon2 = epsilon_elem_ply[elem_num, ply_num, 1] + gamma12 = epsilon_elem_ply[elem_num, ply_num, 2] + + sigma1 = Q11 * epsilon1 + Q12 * epsilon2 + sigma2 = Q21 * epsilon1 + Q22 * epsilon2 + sigma12 = Q66 * gamma12 + + # assigning the stresses to the sigma_elem_ply array + sigma_elem_ply[elem_num, ply_num, 0] = sigma1 + sigma_elem_ply[elem_num, ply_num, 1] = sigma2 + sigma_elem_ply[elem_num, ply_num, 2] = sigma12 + + sigma_c1 = self.sigma_c1 + sigma_t1 = self.sigma_t1 + sigma_c2 = self.sigma_c2 + sigma_t2 = self.sigma_t2 + sigma_12max = self.sigma_12max + + # defining the constants for the Tsai-Wu Strength Ratios + F1 = 1 / sigma_t1 - 1 / sigma_c1 + F11 = 1 / (sigma_t1 * sigma_c1) + F2 = 1 / sigma_t2 - 1 / sigma_c2 + F22 = 1 / (sigma_t2 * sigma_c2) + F66 = 1 / (sigma_12max**2) + + # Finding the Tsai-Wu Strength Ratios for each ply in each element and storing them in the tsaiwu_sr array + for elem_num in range(4): + for ply_num in range(num_plies): + a = F1 * sigma_elem_ply[elem_num, ply_num, 0] + F2 * sigma_elem_ply[elem_num, ply_num, 1] + b = ( + F11 * sigma_elem_ply[elem_num, ply_num, 0] ** 2 + + F22 * sigma_elem_ply[elem_num, ply_num, 1] ** 2 + + F66 * sigma_elem_ply[elem_num, ply_num, 2] ** 2 + ) + tsaiwu_sr[ielem, elem_num * num_plies + ply_num] = 0.5 * (a + np.sqrt(a**2 + 4 * b)) diff --git a/openaerostruct/transfer/load_transfer.py b/openaerostruct/transfer/load_transfer.py index 9004cfa50..2d8004347 100644 --- a/openaerostruct/transfer/load_transfer.py +++ b/openaerostruct/transfer/load_transfer.py @@ -34,7 +34,7 @@ def setup(self): self.nx = nx = surface["mesh"].shape[0] self.ny = ny = surface["mesh"].shape[1] - if surface["fem_model_type"] == "tube": + if surface["fem_model_type"].lower() == "tube": self.fem_origin = surface["fem_origin"] else: y_upper = surface["data_y_upper"] diff --git a/openaerostruct/utils/check_surface_dict.py b/openaerostruct/utils/check_surface_dict.py index 17c3a10e2..7bc8d7ae3 100755 --- a/openaerostruct/utils/check_surface_dict.py +++ b/openaerostruct/utils/check_surface_dict.py @@ -43,6 +43,7 @@ def check_surface_dict_keys(surface): "E", "G", "yield", + "safety_factor", "mrho", "fem_origin", "wing_weight_ratio", @@ -64,10 +65,29 @@ def check_surface_dict_keys(surface): "data_y_upper", "data_x_lower", "data_y_lower", + # tsaiwu_wingbox structure + "useComposite", # FFD "mx", "my", ] + # keys that are required when useComposite is True + compositeInputs = [ + "safety_factor", + "ply_angles", + "ply_fractions", + "E1", + "E2", + "nu12", + "G12", + "sigma_t1", + "sigma_c1", + "sigma_t2", + "sigma_c2", + "sigma_12max", + ] + + keys_implemented = list(set(keys_implemented + compositeInputs)) for key in surface.keys(): if key not in keys_implemented: @@ -76,3 +96,34 @@ def check_surface_dict_keys(surface): category=RuntimeWarning, stacklevel=2, ) + + # adding checks for using the composite failure model + # check1: if useComposite is True, then the following keys must be present + useComposite = "useComposite" in surface.keys() and surface["useComposite"] + if useComposite: + for key in compositeInputs: + if key not in surface.keys(): + raise ValueError( + f"{key} not found in surface dict, when `useComposite` is True, the following keys must be present: {compositeInputs}", + ) + + # check2: if useComposite is True, then 'fem_model_type' must be 'wingbox' + if useComposite and surface.get("fem_model_type", "") != "wingbox": + raise ValueError( + "`fem_model_type` must be 'wingbox' when `useComposite` is True", + ) + + # check3: if useComposite is True, then length of ply_angles and ply_fractions must be equal + if useComposite: + if len(surface["ply_angles"]) != len(surface["ply_fractions"]): + raise ValueError( + "Length of `ply_angles` and `ply_fractions` arrays must be equal", + ) + + # check5: if useComposite is True, then the ply fractions should add to 1 + if useComposite: + plyFracSum = sum(surface["ply_fractions"]) + if abs(plyFracSum - 1) > 1e-2: + raise ValueError( + f"Sum of `ply_fractions` ({surface['ply_fractions']}) is {plyFracSum} must be 1.", + ) diff --git a/openaerostruct/utils/plot_wing.py b/openaerostruct/utils/plot_wing.py index 0dbb05176..cb46b9e64 100644 --- a/openaerostruct/utils/plot_wing.py +++ b/openaerostruct/utils/plot_wing.py @@ -10,7 +10,6 @@ """ - import sys import numpy as np from openmdao.recorders.sqlite_reader import SqliteCaseReader diff --git a/openaerostruct/utils/plot_wingbox.py b/openaerostruct/utils/plot_wingbox.py index 22aa10298..465660c21 100644 --- a/openaerostruct/utils/plot_wingbox.py +++ b/openaerostruct/utils/plot_wingbox.py @@ -7,7 +7,6 @@ The usage, for example, if running a case from the `examples` directory would be `python ../utils/plot_wingbox.py aerostruct.db`. """ - import sys import numpy as np from openmdao.recorders.sqlite_reader import SqliteCaseReader diff --git a/openaerostruct/utils/testing.py b/openaerostruct/utils/testing.py index cccd2a9d4..262e7f8e6 100644 --- a/openaerostruct/utils/testing.py +++ b/openaerostruct/utils/testing.py @@ -144,7 +144,8 @@ def get_default_surfaces(): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, @@ -195,7 +196,8 @@ def get_ground_effect_surfaces(): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct.py b/tests/integration_tests/test_aerostruct.py index 27f67ac46..4214a4556 100644 --- a/tests/integration_tests/test_aerostruct.py +++ b/tests/integration_tests/test_aerostruct.py @@ -48,7 +48,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct_analysis.py b/tests/integration_tests/test_aerostruct_analysis.py index 52c407651..98757d8e1 100644 --- a/tests/integration_tests/test_aerostruct_analysis.py +++ b/tests/integration_tests/test_aerostruct_analysis.py @@ -47,7 +47,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct_analysis_Sref.py b/tests/integration_tests/test_aerostruct_analysis_Sref.py index eb2327e66..7a18a31b7 100644 --- a/tests/integration_tests/test_aerostruct_analysis_Sref.py +++ b/tests/integration_tests/test_aerostruct_analysis_Sref.py @@ -47,7 +47,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct_analysis_compressible.py b/tests/integration_tests/test_aerostruct_analysis_compressible.py index 9d12a6da1..b2bc51e2c 100644 --- a/tests/integration_tests/test_aerostruct_analysis_compressible.py +++ b/tests/integration_tests/test_aerostruct_analysis_compressible.py @@ -47,7 +47,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct_engine_thrusts.py b/tests/integration_tests/test_aerostruct_engine_thrusts.py index 11917775a..e4f0f6539 100644 --- a/tests/integration_tests/test_aerostruct_engine_thrusts.py +++ b/tests/integration_tests/test_aerostruct_engine_thrusts.py @@ -48,7 +48,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct_ffd.py b/tests/integration_tests/test_aerostruct_ffd.py index 0ee954db1..96b3decdb 100644 --- a/tests/integration_tests/test_aerostruct_ffd.py +++ b/tests/integration_tests/test_aerostruct_ffd.py @@ -60,7 +60,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct_groundeffect.py b/tests/integration_tests/test_aerostruct_groundeffect.py index 80439e706..1bf3b9249 100644 --- a/tests/integration_tests/test_aerostruct_groundeffect.py +++ b/tests/integration_tests/test_aerostruct_groundeffect.py @@ -49,7 +49,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct_point_loads.py b/tests/integration_tests/test_aerostruct_point_loads.py index c19467614..56d5b1f6e 100644 --- a/tests/integration_tests/test_aerostruct_point_loads.py +++ b/tests/integration_tests/test_aerostruct_point_loads.py @@ -48,7 +48,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_aerostruct_tsaiwuwingbox_analysis.py b/tests/integration_tests/test_aerostruct_tsaiwuwingbox_analysis.py new file mode 100644 index 000000000..748924495 --- /dev/null +++ b/tests/integration_tests/test_aerostruct_tsaiwuwingbox_analysis.py @@ -0,0 +1,425 @@ +import unittest +import numpy as np + +from openaerostruct.geometry.utils import generate_mesh + +from openaerostruct.integration.aerostruct_groups import AerostructGeometry, AerostructPoint + +import openmdao.api as om + + +# Provide coordinates for a portion of an airfoil for the wingbox cross-section as an nparray with dtype=complex (to work with the complex-step approximation for derivatives). +# These should be for an airfoil with the chord scaled to 1. +# We use the 10% to 60% portion of the NASA SC2-0612 airfoil for this case +# We use the coordinates available from airfoiltools.com. Using such a large number of coordinates is not necessary. +# The first and last x-coordinates of the upper and lower skins must be the same + +upper_x = np.array( + [ + 0.1, + 0.11, + 0.12, + 0.13, + 0.14, + 0.15, + 0.16, + 0.17, + 0.18, + 0.19, + 0.2, + 0.21, + 0.22, + 0.23, + 0.24, + 0.25, + 0.26, + 0.27, + 0.28, + 0.29, + 0.3, + 0.31, + 0.32, + 0.33, + 0.34, + 0.35, + 0.36, + 0.37, + 0.38, + 0.39, + 0.4, + 0.41, + 0.42, + 0.43, + 0.44, + 0.45, + 0.46, + 0.47, + 0.48, + 0.49, + 0.5, + 0.51, + 0.52, + 0.53, + 0.54, + 0.55, + 0.56, + 0.57, + 0.58, + 0.59, + 0.6, + ], + dtype="complex128", +) +lower_x = np.array( + [ + 0.1, + 0.11, + 0.12, + 0.13, + 0.14, + 0.15, + 0.16, + 0.17, + 0.18, + 0.19, + 0.2, + 0.21, + 0.22, + 0.23, + 0.24, + 0.25, + 0.26, + 0.27, + 0.28, + 0.29, + 0.3, + 0.31, + 0.32, + 0.33, + 0.34, + 0.35, + 0.36, + 0.37, + 0.38, + 0.39, + 0.4, + 0.41, + 0.42, + 0.43, + 0.44, + 0.45, + 0.46, + 0.47, + 0.48, + 0.49, + 0.5, + 0.51, + 0.52, + 0.53, + 0.54, + 0.55, + 0.56, + 0.57, + 0.58, + 0.59, + 0.6, + ], + dtype="complex128", +) +upper_y = np.array( + [ + 0.0447, + 0.046, + 0.0472, + 0.0484, + 0.0495, + 0.0505, + 0.0514, + 0.0523, + 0.0531, + 0.0538, + 0.0545, + 0.0551, + 0.0557, + 0.0563, + 0.0568, + 0.0573, + 0.0577, + 0.0581, + 0.0585, + 0.0588, + 0.0591, + 0.0593, + 0.0595, + 0.0597, + 0.0599, + 0.06, + 0.0601, + 0.0602, + 0.0602, + 0.0602, + 0.0602, + 0.0602, + 0.0601, + 0.06, + 0.0599, + 0.0598, + 0.0596, + 0.0594, + 0.0592, + 0.0589, + 0.0586, + 0.0583, + 0.058, + 0.0576, + 0.0572, + 0.0568, + 0.0563, + 0.0558, + 0.0553, + 0.0547, + 0.0541, + ], + dtype="complex128", +) +lower_y = np.array( + [ + -0.0447, + -0.046, + -0.0473, + -0.0485, + -0.0496, + -0.0506, + -0.0515, + -0.0524, + -0.0532, + -0.054, + -0.0547, + -0.0554, + -0.056, + -0.0565, + -0.057, + -0.0575, + -0.0579, + -0.0583, + -0.0586, + -0.0589, + -0.0592, + -0.0594, + -0.0595, + -0.0596, + -0.0597, + -0.0598, + -0.0598, + -0.0598, + -0.0598, + -0.0597, + -0.0596, + -0.0594, + -0.0592, + -0.0589, + -0.0586, + -0.0582, + -0.0578, + -0.0573, + -0.0567, + -0.0561, + -0.0554, + -0.0546, + -0.0538, + -0.0529, + -0.0519, + -0.0509, + -0.0497, + -0.0485, + -0.0472, + -0.0458, + -0.0444, + ], + dtype="complex128", +) + + +class Test(unittest.TestCase): + def test(self): + # Create a dictionary to store options about the surface + mesh_dict = { + "num_y": 5, + "num_x": 3, + "wing_type": "CRM", + "symmetry": True, + "num_twist_cp": 6, + "chord_cos_spacing": 0, + "span_cos_spacing": 0, + } + + mesh, twist_cp = generate_mesh(mesh_dict) + + surf_dict = { + # Wing definition + "name": "wing", # name of the surface + "symmetry": True, # if true, model one half of wing + # reflected across the plane y = 0 + "S_ref_type": "wetted", # how we compute the wing area, + # can be 'wetted' or 'projected' + "fem_model_type": "wingbox", # NOTE: testing the Tsai Wu wingbox model + # "spar_thickness_cp": np.ones(4)*0.1, + # "skin_thickness_cp": np.ones(4)*0.1, + "spar_thickness_cp": np.array([0.004, 0.005, 0.005, 0.008, 0.008, 0.01]), # [m] + "skin_thickness_cp": np.array([0.005, 0.01, 0.015, 0.020, 0.025, 0.026]), + # "spar_thickness_cp": np.ones(6) * 0.1, # [m] + # "skin_thickness_cp": np.ones(6) * 0.1, + "twist_cp": np.array([4.0, 5.0, 8.0, 8.0, 8.0, 9.0]), + "mesh": mesh, + "data_x_upper": upper_x, + "data_x_lower": lower_x, + "data_y_upper": upper_y, + "data_y_lower": lower_y, + "strength_factor_for_upper_skin": 1.0, + # Aerodynamic performance of the lifting surface at + # an angle of attack of 0 (alpha=0). + # These CL0 and CD0 values are added to the CL and CD + # obtained from aerodynamic analysis of the surface to get + # the total CL and CD. + # These CL0 and CD0 values do not vary wrt alpha. + "CL0": 0.0, # CL of the surface at alpha=0 + "CD0": 0.0078, # CD of the surface at alpha=0 + # Airfoil properties for viscous drag calculation + "k_lam": 0.05, # percentage of chord with laminar + # flow, used for viscous drag + "t_over_c_cp": np.array([0.08, 0.08, 0.08, 0.10, 0.10, 0.08]), + "original_wingbox_airfoil_t_over_c": 0.12, + "c_max_t": 0.38, # chordwise location of maximum thickness + "with_viscous": True, + "with_wave": True, # if true, compute wave drag + # Structural values are based on aluminum 7075 + "E": 73.1e9, # [Pa] Young's modulus + "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor + "mrho": 1550, # [kg/m^3] material density #NOTE: CFRP density + "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin + # 'fem_origin' : 0.35, # normalized chordwise location of the spar + "wing_weight_ratio": 1.25, + "struct_weight_relief": True, # True to add the weight of the structure to the loads on the structure + "distributed_fuel_weight": False, + # Constraints + "exact_failure_constraint": False, # if false, use KS function + "Wf_reserve": 15000.0, # [kg] reserve fuel mass + "useComposite": True, + "safety_factor": 1.5, + "ply_angles": [0, 45, -45, 90], + "ply_fractions": [0.10, 0.25, 0.25, 0.40], + "E1": 117.7e9, + "E2": 9.7e9, + "nu12": 0.35, + "G12": 4.8e9, + "sigma_t1": 1648.0e6, + "sigma_c1": 1034.0e6, + "sigma_t2": 64.0e6, + "sigma_c2": 228.0e6, + "sigma_12max": 71.0e6, + } + + surfaces = [surf_dict] + + # Create the problem and assign the model group + prob = om.Problem() + + # Add problem information as an independent variables component + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=0.85 * 295.07, units="m/s") + indep_var_comp.add_output("alpha", val=0.0, units="deg") + indep_var_comp.add_output("Mach_number", val=0.85) + indep_var_comp.add_output("re", val=0.348 * 295.07 * 0.85 * 1.0 / (1.43 * 1e-5), units="1/m") + indep_var_comp.add_output("rho", val=0.348, units="kg/m**3") + indep_var_comp.add_output("CT", val=0.53 / 3600, units="1/s") + indep_var_comp.add_output("R", val=14.307e6, units="m") + indep_var_comp.add_output("W0", val=148000 + surf_dict["Wf_reserve"], units="kg") + indep_var_comp.add_output("speed_of_sound", val=295.07, units="m/s") + indep_var_comp.add_output("load_factor", val=1.0) + indep_var_comp.add_output("empty_cg", val=np.zeros((3)), units="m") + + prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + # Loop over each surface in the surfaces list + for surface in surfaces: + # Get the surface name and create a group to contain components + # only for this surface + name = surface["name"] + + aerostruct_group = AerostructGeometry(surface=surface) + + # Add tmp_group to the problem with the name of the surface. + prob.model.add_subsystem(name, aerostruct_group) + + # Loop through and add a certain number of aero points + for i in range(1): + point_name = "AS_point_{}".format(i) + # Connect the parameters within the model for each aero point + + # Create the aero point group and add it to the model + AS_point = AerostructPoint(surfaces=surfaces) + + prob.model.add_subsystem(point_name, AS_point) + + # Connect flow properties to the analysis point + prob.model.connect("v", point_name + ".v") + prob.model.connect("alpha", point_name + ".alpha") + prob.model.connect("Mach_number", point_name + ".Mach_number") + prob.model.connect("re", point_name + ".re") + prob.model.connect("rho", point_name + ".rho") + prob.model.connect("CT", point_name + ".CT") + prob.model.connect("R", point_name + ".R") + prob.model.connect("W0", point_name + ".W0") + prob.model.connect("speed_of_sound", point_name + ".speed_of_sound") + prob.model.connect("empty_cg", point_name + ".empty_cg") + prob.model.connect("load_factor", point_name + ".load_factor") + + for _surface in surfaces: + com_name = point_name + "." + name + "_perf." + prob.model.connect( + name + ".local_stiff_transformed", point_name + ".coupled." + name + ".local_stiff_transformed" + ) + prob.model.connect(name + ".nodes", point_name + ".coupled." + name + ".nodes") + + # Connect aerodyamic mesh to coupled group mesh + prob.model.connect(name + ".mesh", point_name + ".coupled." + name + ".mesh") + prob.model.connect(name + ".element_mass", point_name + ".coupled." + name + ".element_mass") + + # Connect performance calculation variables + prob.model.connect(name + ".nodes", com_name + "nodes") + prob.model.connect(name + ".cg_location", point_name + "." + "total_perf." + name + "_cg_location") + prob.model.connect( + name + ".structural_mass", point_name + "." + "total_perf." + name + "_structural_mass" + ) + + # Connect wingbox properties to von Mises stress calcs + prob.model.connect(name + ".Qz", com_name + "Qz") + prob.model.connect(name + ".J", com_name + "J") + prob.model.connect(name + ".A_enc", com_name + "A_enc") + prob.model.connect(name + ".htop", com_name + "htop") + prob.model.connect(name + ".hbottom", com_name + "hbottom") + prob.model.connect(name + ".hfront", com_name + "hfront") + prob.model.connect(name + ".hrear", com_name + "hrear") + + prob.model.connect(name + ".spar_thickness", com_name + "spar_thickness") + prob.model.connect(name + ".t_over_c", com_name + "t_over_c") + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["tol"] = 1e-9 + + # Set up the problem + prob.setup() + # om.n2(prob) + + prob.run_model() + + print(prob["AS_point_0.fuelburn"][0]) + print(prob["wing.structural_mass"][0] / 1.25) + print("TsaiWu SR vals:", prob["AS_point_0.wing_perf.tsaiwu_sr"]) + print(prob["AS_point_0.wing_perf.failure"]) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/integration_tests/test_aerostruct_tsaiwuwingbox_opt.py b/tests/integration_tests/test_aerostruct_tsaiwuwingbox_opt.py new file mode 100644 index 000000000..ea7ab8896 --- /dev/null +++ b/tests/integration_tests/test_aerostruct_tsaiwuwingbox_opt.py @@ -0,0 +1,445 @@ +import unittest +import numpy as np + +from openaerostruct.geometry.utils import generate_mesh + +from openaerostruct.integration.aerostruct_groups import AerostructGeometry, AerostructPoint + +import openmdao.api as om + + +# Provide coordinates for a portion of an airfoil for the wingbox cross-section as an nparray with dtype=complex (to work with the complex-step approximation for derivatives). +# These should be for an airfoil with the chord scaled to 1. +# We use the 10% to 60% portion of the NASA SC2-0612 airfoil for this case +# We use the coordinates available from airfoiltools.com. Using such a large number of coordinates is not necessary. +# The first and last x-coordinates of the upper and lower skins must be the same + +upper_x = np.array( + [ + 0.1, + 0.11, + 0.12, + 0.13, + 0.14, + 0.15, + 0.16, + 0.17, + 0.18, + 0.19, + 0.2, + 0.21, + 0.22, + 0.23, + 0.24, + 0.25, + 0.26, + 0.27, + 0.28, + 0.29, + 0.3, + 0.31, + 0.32, + 0.33, + 0.34, + 0.35, + 0.36, + 0.37, + 0.38, + 0.39, + 0.4, + 0.41, + 0.42, + 0.43, + 0.44, + 0.45, + 0.46, + 0.47, + 0.48, + 0.49, + 0.5, + 0.51, + 0.52, + 0.53, + 0.54, + 0.55, + 0.56, + 0.57, + 0.58, + 0.59, + 0.6, + ], + dtype="complex128", +) +lower_x = np.array( + [ + 0.1, + 0.11, + 0.12, + 0.13, + 0.14, + 0.15, + 0.16, + 0.17, + 0.18, + 0.19, + 0.2, + 0.21, + 0.22, + 0.23, + 0.24, + 0.25, + 0.26, + 0.27, + 0.28, + 0.29, + 0.3, + 0.31, + 0.32, + 0.33, + 0.34, + 0.35, + 0.36, + 0.37, + 0.38, + 0.39, + 0.4, + 0.41, + 0.42, + 0.43, + 0.44, + 0.45, + 0.46, + 0.47, + 0.48, + 0.49, + 0.5, + 0.51, + 0.52, + 0.53, + 0.54, + 0.55, + 0.56, + 0.57, + 0.58, + 0.59, + 0.6, + ], + dtype="complex128", +) +upper_y = np.array( + [ + 0.0447, + 0.046, + 0.0472, + 0.0484, + 0.0495, + 0.0505, + 0.0514, + 0.0523, + 0.0531, + 0.0538, + 0.0545, + 0.0551, + 0.0557, + 0.0563, + 0.0568, + 0.0573, + 0.0577, + 0.0581, + 0.0585, + 0.0588, + 0.0591, + 0.0593, + 0.0595, + 0.0597, + 0.0599, + 0.06, + 0.0601, + 0.0602, + 0.0602, + 0.0602, + 0.0602, + 0.0602, + 0.0601, + 0.06, + 0.0599, + 0.0598, + 0.0596, + 0.0594, + 0.0592, + 0.0589, + 0.0586, + 0.0583, + 0.058, + 0.0576, + 0.0572, + 0.0568, + 0.0563, + 0.0558, + 0.0553, + 0.0547, + 0.0541, + ], + dtype="complex128", +) +lower_y = np.array( + [ + -0.0447, + -0.046, + -0.0473, + -0.0485, + -0.0496, + -0.0506, + -0.0515, + -0.0524, + -0.0532, + -0.054, + -0.0547, + -0.0554, + -0.056, + -0.0565, + -0.057, + -0.0575, + -0.0579, + -0.0583, + -0.0586, + -0.0589, + -0.0592, + -0.0594, + -0.0595, + -0.0596, + -0.0597, + -0.0598, + -0.0598, + -0.0598, + -0.0598, + -0.0597, + -0.0596, + -0.0594, + -0.0592, + -0.0589, + -0.0586, + -0.0582, + -0.0578, + -0.0573, + -0.0567, + -0.0561, + -0.0554, + -0.0546, + -0.0538, + -0.0529, + -0.0519, + -0.0509, + -0.0497, + -0.0485, + -0.0472, + -0.0458, + -0.0444, + ], + dtype="complex128", +) + + +class Test(unittest.TestCase): + def test(self): + # Create a dictionary to store options about the surface + mesh_dict = { + "num_y": 5, + "num_x": 3, + "wing_type": "CRM", + "symmetry": True, + "num_twist_cp": 6, + "chord_cos_spacing": 0, + "span_cos_spacing": 0, + } + + mesh, twist_cp = generate_mesh(mesh_dict) + + surf_dict = { + # Wing definition + "name": "wing", # name of the surface + "symmetry": True, # if true, model one half of wing + # reflected across the plane y = 0 + "S_ref_type": "wetted", # how we compute the wing area, + # can be 'wetted' or 'projected' + "fem_model_type": "wingbox", # NOTE: testing the Tsai Wu wingbox model + "spar_thickness_cp": np.array([0.004, 0.005, 0.005, 0.008, 0.008, 0.01]), # [m] + "skin_thickness_cp": np.array([0.005, 0.01, 0.015, 0.020, 0.025, 0.026]), + "twist_cp": np.array([4.0, 5.0, 8.0, 8.0, 8.0, 9.0]), + "mesh": mesh, + "data_x_upper": upper_x, + "data_x_lower": lower_x, + "data_y_upper": upper_y, + "data_y_lower": lower_y, + "strength_factor_for_upper_skin": 1.0, + # Aerodynamic performance of the lifting surface at + # an angle of attack of 0 (alpha=0). + # These CL0 and CD0 values are added to the CL and CD + # obtained from aerodynamic analysis of the surface to get + # the total CL and CD. + # These CL0 and CD0 values do not vary wrt alpha. + "CL0": 0.0, # CL of the surface at alpha=0 + "CD0": 0.0078, # CD of the surface at alpha=0 + # Airfoil properties for viscous drag calculation + "k_lam": 0.05, # percentage of chord with laminar + # flow, used for viscous drag + "t_over_c_cp": np.array([0.08, 0.08, 0.08, 0.10, 0.10, 0.08]), + "original_wingbox_airfoil_t_over_c": 0.12, + "c_max_t": 0.38, # chordwise location of maximum thickness + "with_viscous": True, + "with_wave": True, # if true, compute wave drag + # Structural values are based on aluminum 7075 + # "E": 73.1e9, # [Pa] Young's modulus + # "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) + "E": 62.53e9, # [Pa] Young's modulus # skin composites + "G": 29.71e9, # [Pa] shear modulus # skin composites + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor + "mrho": 1550, # [kg/m^3] material density #NOTE: CFRP density + "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin + # 'fem_origin' : 0.35, # normalized chordwise location of the spar + "wing_weight_ratio": 1.25, + "struct_weight_relief": True, # True to add the weight of the structure to the loads on the structure + "distributed_fuel_weight": False, + # Constraints + "exact_failure_constraint": False, # if false, use KS function + "Wf_reserve": 15000.0, # [kg] reserve fuel mass + "span": 58, # [m] wingspan + "useComposite": True, + "safety_factor": 1.5, + "ply_angles": [0, 45, -45, 90], + "ply_fractions": [0.10, 0.25, 0.25, 0.40], + "E1": 117.7e9, + "E2": 9.7e9, + "nu12": 0.35, + "G12": 4.8e9, + "sigma_t1": 1648.0e6, + "sigma_c1": 1034.0e6, + "sigma_t2": 64.0e6, + "sigma_c2": 228.0e6, + "sigma_12max": 71.0e6, + } + + surfaces = [surf_dict] + + # Create the problem and assign the model group + prob = om.Problem() + + # Add problem information as an independent variables component + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=0.85 * 295.07, units="m/s") + indep_var_comp.add_output("alpha", val=0.0, units="deg") + indep_var_comp.add_output("Mach_number", val=0.85) + indep_var_comp.add_output("re", val=0.348 * 295.07 * 0.85 * 1.0 / (1.43 * 1e-5), units="1/m") + indep_var_comp.add_output("rho", val=0.348, units="kg/m**3") + indep_var_comp.add_output("CT", val=0.53 / 3600, units="1/s") + indep_var_comp.add_output("R", val=14.307e6, units="m") + indep_var_comp.add_output("W0", val=148000 + surf_dict["Wf_reserve"], units="kg") + indep_var_comp.add_output("speed_of_sound", val=295.07, units="m/s") + indep_var_comp.add_output("load_factor", val=1.0) + indep_var_comp.add_output("empty_cg", val=np.zeros((3)), units="m") + + prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + # Loop over each surface in the surfaces list + for surface in surfaces: + # Get the surface name and create a group to contain components + # only for this surface + name = surface["name"] + + aerostruct_group = AerostructGeometry(surface=surface) + + # Add tmp_group to the problem with the name of the surface. + prob.model.add_subsystem(name, aerostruct_group) + + # Loop through and add a certain number of aero points + for i in range(1): + point_name = "AS_point_{}".format(i) + # Connect the parameters within the model for each aero point + + # Create the aero point group and add it to the model + AS_point = AerostructPoint(surfaces=surfaces) + + prob.model.add_subsystem(point_name, AS_point) + + # Connect flow properties to the analysis point + prob.model.connect("v", point_name + ".v") + prob.model.connect("alpha", point_name + ".alpha") + prob.model.connect("Mach_number", point_name + ".Mach_number") + prob.model.connect("re", point_name + ".re") + prob.model.connect("rho", point_name + ".rho") + prob.model.connect("CT", point_name + ".CT") + prob.model.connect("R", point_name + ".R") + prob.model.connect("W0", point_name + ".W0") + prob.model.connect("speed_of_sound", point_name + ".speed_of_sound") + prob.model.connect("empty_cg", point_name + ".empty_cg") + prob.model.connect("load_factor", point_name + ".load_factor") + + for _surface in surfaces: + com_name = point_name + "." + name + "_perf." + prob.model.connect( + name + ".local_stiff_transformed", point_name + ".coupled." + name + ".local_stiff_transformed" + ) + prob.model.connect(name + ".nodes", point_name + ".coupled." + name + ".nodes") + + # Connect aerodyamic mesh to coupled group mesh + prob.model.connect(name + ".mesh", point_name + ".coupled." + name + ".mesh") + prob.model.connect(name + ".element_mass", point_name + ".coupled." + name + ".element_mass") + + # Connect performance calculation variables + prob.model.connect(name + ".nodes", com_name + "nodes") + prob.model.connect(name + ".cg_location", point_name + "." + "total_perf." + name + "_cg_location") + prob.model.connect( + name + ".structural_mass", point_name + "." + "total_perf." + name + "_structural_mass" + ) + + # Connect wingbox properties to von Mises stress calcs + prob.model.connect(name + ".Qz", com_name + "Qz") + prob.model.connect(name + ".J", com_name + "J") + prob.model.connect(name + ".A_enc", com_name + "A_enc") + prob.model.connect(name + ".htop", com_name + "htop") + prob.model.connect(name + ".hbottom", com_name + "hbottom") + prob.model.connect(name + ".hfront", com_name + "hfront") + prob.model.connect(name + ".hrear", com_name + "hrear") + + prob.model.connect(name + ".spar_thickness", com_name + "spar_thickness") + prob.model.connect(name + ".t_over_c", com_name + "t_over_c") + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["tol"] = 1e-9 + + prob.model.add_objective("AS_point_0.fuelburn", scaler=1e-5) + + prob.model.add_design_var("wing.twist_cp", lower=-15.0, upper=15.0, scaler=0.1) + prob.model.add_design_var("wing.spar_thickness_cp", lower=0.003, upper=0.1, scaler=1e2) + prob.model.add_design_var("wing.skin_thickness_cp", lower=0.003, upper=0.1, scaler=1e2) + prob.model.add_design_var("wing.geometry.t_over_c_cp", lower=0.07, upper=0.2, scaler=10.0) + prob.model.add_design_var("wing.geometry.span", lower=55.0, upper=60.0, scaler=2e-2) + + prob.model.add_constraint("AS_point_0.CL", equals=0.5) + prob.model.add_constraint("AS_point_0.wing_perf.failure", upper=0.0) + + # recorder = om.SqliteRecorder("tsaiwu_opt.db") + # prob.driver.add_recorder(recorder) + + # # We could also just use prob.driver.recording_options['includes']=['*'] here, but for large meshes the database file becomes extremely large. So we just select the variables we need. + # prob.driver.recording_options["includes"] = ['*'] + + # prob.driver.recording_options["record_objectives"] = True + # prob.driver.recording_options["record_constraints"] = True + # prob.driver.recording_options["record_desvars"] = True + # prob.driver.recording_options["record_inputs"] = True + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["tol"] = 1e-6 + prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + + # Set up the problem + prob.setup() + + prob.run_driver() + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/integration_tests/test_aerostruct_wingbox_+weight_analysis.py b/tests/integration_tests/test_aerostruct_wingbox_+weight_analysis.py index c1be382da..b8ec4ff2c 100644 --- a/tests/integration_tests/test_aerostruct_wingbox_+weight_analysis.py +++ b/tests/integration_tests/test_aerostruct_wingbox_+weight_analysis.py @@ -292,7 +292,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_aerostruct_wingbox_analysis.py b/tests/integration_tests/test_aerostruct_wingbox_analysis.py index b4f96fe50..ae461a8cc 100644 --- a/tests/integration_tests/test_aerostruct_wingbox_analysis.py +++ b/tests/integration_tests/test_aerostruct_wingbox_analysis.py @@ -292,7 +292,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_aerostruct_wingbox_fuel_vol_constraint_opt.py b/tests/integration_tests/test_aerostruct_wingbox_fuel_vol_constraint_opt.py index b6e4d7a77..c492f4b70 100644 --- a/tests/integration_tests/test_aerostruct_wingbox_fuel_vol_constraint_opt.py +++ b/tests/integration_tests/test_aerostruct_wingbox_fuel_vol_constraint_opt.py @@ -297,7 +297,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_aerostruct_wingbox_opt.py b/tests/integration_tests/test_aerostruct_wingbox_opt.py index fbb679151..bca40c5c1 100644 --- a/tests/integration_tests/test_aerostruct_wingbox_opt.py +++ b/tests/integration_tests/test_aerostruct_wingbox_opt.py @@ -292,7 +292,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_aerostruct_wingbox_wave_fuel_vol_constraint_opt.py b/tests/integration_tests/test_aerostruct_wingbox_wave_fuel_vol_constraint_opt.py index 3bb331f00..bbfe4dde7 100644 --- a/tests/integration_tests/test_aerostruct_wingbox_wave_fuel_vol_constraint_opt.py +++ b/tests/integration_tests/test_aerostruct_wingbox_wave_fuel_vol_constraint_opt.py @@ -298,7 +298,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_morphing_aerostruct.py b/tests/integration_tests/test_morphing_aerostruct.py index 38c7bc3f8..017348a3e 100644 --- a/tests/integration_tests/test_morphing_aerostruct.py +++ b/tests/integration_tests/test_morphing_aerostruct.py @@ -55,7 +55,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_multipoint_parallel.py b/tests/integration_tests/test_multipoint_parallel.py index 3fea70afe..36622098e 100644 --- a/tests/integration_tests/test_multipoint_parallel.py +++ b/tests/integration_tests/test_multipoint_parallel.py @@ -110,7 +110,8 @@ def test_multipoint_MPI(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin "wing_weight_ratio": 1.25, diff --git a/tests/integration_tests/test_multipoint_wingbox_aerostruct.py b/tests/integration_tests/test_multipoint_wingbox_aerostruct.py index 64217956f..3df443b90 100644 --- a/tests/integration_tests/test_multipoint_wingbox_aerostruct.py +++ b/tests/integration_tests/test_multipoint_wingbox_aerostruct.py @@ -298,7 +298,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_multipoint_wingbox_derivs.py b/tests/integration_tests/test_multipoint_wingbox_derivs.py index 1b79c22b9..a29203268 100644 --- a/tests/integration_tests/test_multipoint_wingbox_derivs.py +++ b/tests/integration_tests/test_multipoint_wingbox_derivs.py @@ -296,7 +296,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_scaneagle.py b/tests/integration_tests/test_scaneagle.py index b67ba5e45..6f23204f7 100644 --- a/tests/integration_tests/test_scaneagle.py +++ b/tests/integration_tests/test_scaneagle.py @@ -92,6 +92,7 @@ def setUp(self): "E": 85.0e9, "G": 25.0e9, "yield": 350.0e6, + "safety_factor": 1, "mrho": 1.6e3, "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 1.0, # multiplicative factor on the computed structural weight diff --git a/tests/integration_tests/test_simple_rect_AS.py b/tests/integration_tests/test_simple_rect_AS.py index 36763bd16..9a50f968d 100644 --- a/tests/integration_tests/test_simple_rect_AS.py +++ b/tests/integration_tests/test_simple_rect_AS.py @@ -49,7 +49,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 2.0, diff --git a/tests/integration_tests/test_struct.py b/tests/integration_tests/test_struct.py index 0e9afc63b..a8143f0ed 100644 --- a/tests/integration_tests/test_struct.py +++ b/tests/integration_tests/test_struct.py @@ -26,7 +26,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness diff --git a/tests/integration_tests/test_struct_analysis.py b/tests/integration_tests/test_struct_analysis.py index 7c8a5b88a..86fe0d479 100644 --- a/tests/integration_tests/test_struct_analysis.py +++ b/tests/integration_tests/test_struct_analysis.py @@ -26,7 +26,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness diff --git a/tests/integration_tests/test_struct_analysis_2g.py b/tests/integration_tests/test_struct_analysis_2g.py index 28ec53128..d43f0c3d8 100644 --- a/tests/integration_tests/test_struct_analysis_2g.py +++ b/tests/integration_tests/test_struct_analysis_2g.py @@ -25,7 +25,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness diff --git a/tests/integration_tests/test_struct_engine_thrusts.py b/tests/integration_tests/test_struct_engine_thrusts.py index 48db58dde..fbf71b8d6 100644 --- a/tests/integration_tests/test_struct_engine_thrusts.py +++ b/tests/integration_tests/test_struct_engine_thrusts.py @@ -25,7 +25,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.5, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness @@ -83,7 +84,8 @@ def test_multiple_masses(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.5, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness diff --git a/tests/integration_tests/test_struct_no_loads.py b/tests/integration_tests/test_struct_no_loads.py index d1d972cbe..c2d621a80 100644 --- a/tests/integration_tests/test_struct_no_loads.py +++ b/tests/integration_tests/test_struct_no_loads.py @@ -25,7 +25,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness diff --git a/tests/integration_tests/test_struct_point_masses.py b/tests/integration_tests/test_struct_point_masses.py index 7e138165a..f4656e854 100644 --- a/tests/integration_tests/test_struct_point_masses.py +++ b/tests/integration_tests/test_struct_point_masses.py @@ -25,7 +25,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.5, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness @@ -83,7 +84,8 @@ def test_multiple_masses(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.5, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness diff --git a/tests/integration_tests/test_v1_aerostruct_analysis.py b/tests/integration_tests/test_v1_aerostruct_analysis.py index e7fb2dab7..f69bba406 100644 --- a/tests/integration_tests/test_v1_aerostruct_analysis.py +++ b/tests/integration_tests/test_v1_aerostruct_analysis.py @@ -54,7 +54,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 1.0, diff --git a/tests/integration_tests/test_v1_aerostruct_opt.py b/tests/integration_tests/test_v1_aerostruct_opt.py index 486ac1d90..164e08a79 100644 --- a/tests/integration_tests/test_v1_aerostruct_opt.py +++ b/tests/integration_tests/test_v1_aerostruct_opt.py @@ -54,7 +54,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "wing_weight_ratio": 1.0, diff --git a/tests/integration_tests/test_v1_struct_analysis.py b/tests/integration_tests/test_v1_struct_analysis.py index 8f9eb49e2..5ee8770f9 100644 --- a/tests/integration_tests/test_v1_struct_analysis.py +++ b/tests/integration_tests/test_v1_struct_analysis.py @@ -33,7 +33,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness diff --git a/tests/integration_tests/test_v1_struct_opt.py b/tests/integration_tests/test_v1_struct_opt.py index 9149b4443..14ff9b5ad 100644 --- a/tests/integration_tests/test_v1_struct_opt.py +++ b/tests/integration_tests/test_v1_struct_opt.py @@ -33,7 +33,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 70.0e9, # [Pa] Young's modulus of the spar "G": 30.0e9, # [Pa] shear modulus of the spar - "yield": 500.0e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case + "yield": 500.0e6, + "safety_factor": 2.5, # [Pa] yield stress divided by 2.5 for limiting case "mrho": 3.0e3, # [kg/m^3] material density "fem_origin": 0.35, # normalized chordwise location of the spar "t_over_c_cp": np.array([0.15]), # maximum airfoil thickness diff --git a/tests/integration_tests/test_wingbox_analysis.py b/tests/integration_tests/test_wingbox_analysis.py index 67b12c60a..4ad837666 100644 --- a/tests/integration_tests/test_wingbox_analysis.py +++ b/tests/integration_tests/test_wingbox_analysis.py @@ -274,7 +274,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_wingbox_derivs.py b/tests/integration_tests/test_wingbox_derivs.py index 9db48e92c..d5e7380a6 100644 --- a/tests/integration_tests/test_wingbox_derivs.py +++ b/tests/integration_tests/test_wingbox_derivs.py @@ -293,7 +293,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar diff --git a/tests/integration_tests/test_wingbox_distributed_fuel.py b/tests/integration_tests/test_wingbox_distributed_fuel.py index 876f7cc8a..ffe847990 100644 --- a/tests/integration_tests/test_wingbox_distributed_fuel.py +++ b/tests/integration_tests/test_wingbox_distributed_fuel.py @@ -297,7 +297,8 @@ def test(self): # Structural values are based on aluminum 7075 "E": 73.1e9, # [Pa] Young's modulus "G": (73.1e9 / 2 / 1.33), # [Pa] shear modulus (calculated using E and the Poisson's ratio here) - "yield": (420.0e6 / 1.5), # [Pa] allowable yield stress + "yield": 420.0e6, # [Pa] yield stress + "safety_factor": 1.5, # safety factor "mrho": 2.78e3, # [kg/m^3] material density "strength_factor_for_upper_skin": 1.0, # the yield stress is multiplied by this factor for the upper skin # 'fem_origin' : 0.35, # normalized chordwise location of the spar