Composite material model and Tsai-Wu failure criterion; introduce safety_factor (#444)

* initial commit for tsaiwu wingbox package

* Make geometry spline interpolation mesh independent

* formatting changes

* implemented "useComposite" in the surfdict for the option of using composite materials and made changes to all the relevant scripts along with the calculation of effective stiffness values for composite materials

* changed the surface dictionary in the test files to work with the new keys

* Made the changes with if statements of "useComposite", bug fixes with numofplies and num_failure_criteria, added raise commands for surf dict

* cleaning all the codes, removing unnecessary comments and updating the group descriptions for the failure groups and tsaiwuwingbox

* bug fixes due to the hardcoded 4 plies

* Modify surface dictionary checks

* Rename surface dict variables to follow OAS conventions

* Undo artifact from previous merge

* Fix error message

* correcting the spatial_beam_functionals

* `black -l 120 .`

* `black -l 120 .`

* Fix docstring in failure_ks

* fixing the files flagged for pull request tests

* removing the call for shwoing n2, and fixing formatting bug

* removing the asert.near.equal import statement from test file

* failure files and check_surface_dict modified to use "safety_factor". All test files using yield/SF modified to have a safety_factor key

* removed unused import statements, adding option for the cases where safety factor is absent

* creating and editing the documentation files for OAS documentation

* updating heirarchy for the documentation files

* Add code excerpt to composites doc page

* documemtation for composite wingbox

* updated the documentation page

* changed the safety factor documentation and added the Pareto front to the walkthrough page

* Simplifying some code in failure_ks

* small docs fix

* variable renaming

* More robust if check

* Check against lower case strings

* move functions to structure.utils

* remove unused var

* update integration tests

* rename test file

* moved effective E and G computation to runscript level to be more general

* documentation update

* moved example script from docs folder to example folder

* updated example script

* fixed stiffness matrix transformation and added unit test

* fixed strain transformation in tsai-wu comp and refactoring

* working on Gray2024 benchmark model

* updated benchmark example

* ruff

* docs minor update

* docs update

* docs minor fixes

* update tsaiwu test

* format

* typos in docs

* fixed safety_factor and removed debugging leftovers

* Alter composite transformations

---------

Co-authored-by: Alasdair Gray <alachris@umich.edu>
Co-authored-by: kanekosh <shugok@umich.edu>

Author: 3 people

openaerostruct/docs/advanced_features.rst

@@ -22,3 +22,4 @@ These examples show some advanced features in OpenAeroStruct.
    advanced_features/customizing_prob_setup.rst
    advanced_features/mphys_coupling.rst
    advanced_features/stability_derivatives.rst
+   advanced_features/composites_walkthrough.rst

openaerostruct/docs/advanced_features/composites_walkthrough.rst

@@ -0,0 +1,196 @@
+.. _Composites Walkthrough:
+
+Composite Material
+==================
+
+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.
+
+Model setup
+-----------
+
+First, define the following parameters in the `surface` dictionary:
+
+- ``useComposite`` is a boolean variable that is set to True to enable the composites model mode.
+- ``safety_factor`` is the factor of safety used for determining the Tsai-Wu based failure
+- ``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``.
+
+Next, call a utility function to compute the effective E and G for the composite material.
+The following function adds ``E`` and ``G`` to ``surface``.
+
+.. code-block:: python
+
+    from openaerostruct.structures.utils import compute_composite_stiffness
+    compute_composite_stiffness(surf_dict)
+
+The rest of the model setup is the same as aluminum wing aerostructural problems.
+OpenAeroStruct will compute the failure metric based on the Tsai-Wu failure criterion instead of the von Mises failure criterion when we set ``useComposite`` to True.
+But this is done automatically within ``AerostructPoint``, so you don't need to do anything about it.
+
+Theory
+------
+
+Approximation of the Moduli of Elasticity
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+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 following procedure.
+
+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.
+
+We then transform the stiffness matrix of each ply as a function of the ply angle by
+
+.. math::
+
+    \bar{Q} = T_\sigma Q T_\varepsilon^{-1}
+
+where :math:`T_\sigma`, :math:`T_\varepsilon` are the stress and strain transformation matrices, respectively.
+
+The effective stiffness matrix for the laminate is found using the weighted sum of the stiffness matrices of the plies,
+using their respective ply fraction:
+
+.. math::
+
+    Q_{eff} = \sum_{i=1}^{n} f_i \bar{Q}_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 = \frac{1}{S_{eff_{11}}}\\
+    G = \frac{1}{S_{eff_{66}}}
+
+These moduli of elasticity values are hence used to determine the stiffness matrix of the entire FEM spatial beam model.
+
+Tsai-Wu Failure Criterion
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+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:
+
+.. math::
+
+    \begin{pmatrix}
+    \epsilon_1 \\
+    \epsilon_2 \\
+    \gamma_{12}
+    \end{pmatrix}
+    =
+    [T_\varepsilon]
+    \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**, 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 criterion for each ply.
+The Tsai-Wu failure criterion 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 criterion, 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 calculate 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}{\text{safety factor}}
+
+
+The failure is determined by the value of :math:`\hat{g}_{KS}(\rho, g)` exceeding 0.
+
+
+Example runscript
+-----------------
+
+Here is an example runscript of composite wing aerostructural optimization.
+This roughly follows the setup of "Simple Transonic Wing" by `Gray and Martins 2024 <https://www.researchgate.net/publication/377154425_A_Proposed_Benchmark_Model_for_Practical_Aeroelastic_Optimization_of_Aircraft_Wings>`_.
+
+.. embed-code::
+  ../examples/run_aerostruct_composite_benchmark_wing.py

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,

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,

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.

openaerostruct/docs/user_reference/mesh_surface_dict.rst

@@ -243,9 +243,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

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,

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,

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,