Add a failure criteria for using Composite Materials (Tsai Wu Failure Criteria)

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

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
 

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
 

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

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

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"]
+    numofplies = 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((numofplies, 3, 3))
+    Q_bar_eff = np.zeros((3, 3))
+    for i in range(numofplies):
+        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)
@@ -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"] == "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"] == "wingbox":  # same for both isotropic and composite wingbox
             promotes = ["A_int"]
         else:
             promotes = []
@@ -180,6 +262,11 @@ def setup(self):
             )
 
         elif surface["fem_model_type"] == "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.

openaerostruct/structures/failure_exact.py

@@ -9,13 +9,19 @@ class FailureExact(om.ExplicitComponent):
 
     Parameters
     ----------
+    for Isotropic structures:
     vonmises[ny-1, 2] : numpy array
         von Mises stress magnitudes for each FEM element.
 
+    for Composite wingbox:
+    tsaiwu_sr[ny-1, 4 * numofplies] : 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.
+        Array of failure conditions. Positive if element has failed. This entity is defined for either
+        failure criteria, vonmises or tsaiwu_sr. # TODO: check this
 
     """
 
@@ -24,19 +30,41 @@ def initialize(self):
 
     def setup(self):
         surface = self.options["surface"]
+        ply_angles = surface["ply_angles"]
+        numofplies = len(ply_angles)
 
         if surface["fem_model_type"] == "tube":
             num_failure_criteria = 2
+
         elif surface["fem_model_type"] == "wingbox":
-            num_failure_criteria = 4
+            if "useComposite" in surface.keys() and surface["useComposite"]:  # using the Composite wingbox
+                num_failure_criteria = 4 * numofplies  # 4 critical elements * number of plies
+            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")
+        self.srlimit = 1 / surface["composite_safety_factor"]
+
+        if "useComposite" in surface.keys() and surface["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 "useComposite" in surface.keys() and surface["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

openaerostruct/structures/failure_ks.py

@@ -21,15 +21,21 @@ class FailureKS(om.ExplicitComponent):
 
     parameters
     ----------
+    for Isotropic structures:
     vonmises[ny-1, 2] : numpy array
         von Mises stress magnitudes for each FEM element.
 
+    for Composite wingbox:
+    tsaiwu_sr[ny-1, 4 * numofplies] : 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. # TODO: check this
 
     """
 
@@ -40,17 +46,30 @@ def initialize(self):
     def setup(self):
         surface = self.options["surface"]
         rho = self.options["rho"]
+        self.useComposite = "useComposite" in self.options["surface"].keys() and self.options["surface"]["useComposite"]
+        if self.useComposite:
+            self.numofplies = len(surface["ply_angles"])
 
         if surface["fem_model_type"] == "tube":
             num_failure_criteria = 2
+
         elif surface["fem_model_type"] == "wingbox":
-            num_failure_criteria = 4
+            if self.useComposite:  # using the Composite wingbox
+                num_failure_criteria = 4 * self.numofplies  # 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)
+        if self.useComposite:
+            self.add_input("tsaiwu_sr", val=np.zeros((self.ny - 1, num_failure_criteria)), units=None)
+            self.composite_safety_factor = surface["composite_safety_factor"]
+            self.srlimit = 1 / self.composite_safety_factor
 
+        else:
+            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.sigma = surface["yield"]
         self.rho = rho
 
@@ -59,35 +78,44 @@ def setup(self):
     def compute(self, inputs, outputs):
         sigma = self.sigma
         rho = self.rho
-        vonmises = inputs["vonmises"]
 
-        fmax = np.max(vonmises / sigma - 1)
+        if self.useComposite:  # using the Composite wingbox
+            stress_array = inputs["tsaiwu_sr"]
+            stress_limit = self.srlimit
+        else:  # using the Isotropic structures
+            stress_array = inputs["vonmises"]
+            stress_limit = sigma
+
+        fmax = np.max(stress_array / 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 / rho * nlog(nsum(nexp(rho * (stress_array / stress_limit - 1 - fmax))))
         outputs["failure"] = fmax + ks
 
     def compute_partials(self, inputs, partials):
-        vonmises = inputs["vonmises"]
-        sigma = self.sigma
-        rho = self.rho
-
-        # Find the location of the max stress constraint
-        fmax = np.max(vonmises / sigma - 1)
-        i, j = np.where((vonmises / sigma - 1) == fmax)
+        if self.useComposite:  # using the Composite wingbox
+            stress_array = inputs["tsaiwu_sr"]
+            stress_limit = self.srlimit
+        else:  # using the Isotropic structures
+            stress_array = inputs["vonmises"]
+            stress_limit = self.sigma
+
+        fmax = np.max(stress_array / stress_limit - 1)
+        i, j = np.where((stress_array / 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 / stress_limit - fmax - 1))))
+        tempb = np.exp(self.rho * (stress_array / 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 / stress_limit
+        derivs[i, j] += fmaxb / stress_limit
 
-        # Reshape and save them to the jac dict
-        partials["failure", "vonmises"] = derivs.reshape(1, -1)
+        if (
+            "useComposite" in self.options["surface"].keys() and self.options["surface"]["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)