Merge pull request #476 from robbievanleeuwen/feature/yield-moment

Add yield moment calculation

Author: robbievanleeuwen

docs/user_guide.rst

@@ -52,7 +52,7 @@ Cross-Section Analysis
 
   * ☑ Second moments of area
   * ☑ Elastic section moduli
-  * ☐ Yield moment
+  * ☑ Yield moment
   * ☑ Radii of gyration
   * ☑ Plastic centroid
   * ☑ Plastic section moduli
@@ -62,7 +62,7 @@ Cross-Section Analysis
 
   * ☑ Second moments of area
   * ☑ Elastic section moduli
-  * ☐ Yield moment
+  * ☑ Yield moment
   * ☑ Radii of gyration
   * ☑ Plastic centroid
   * ☑ Plastic section moduli

docs/user_guide/results.rst

@@ -99,10 +99,12 @@ Geometric Analysis
     ~sectionproperties.analysis.section.Section.get_c
     ~sectionproperties.analysis.section.Section.get_eic
     ~sectionproperties.analysis.section.Section.get_ez
+    ~sectionproperties.analysis.section.Section.get_my
     ~sectionproperties.analysis.section.Section.get_rc
     ~sectionproperties.analysis.section.Section.get_eip
     ~sectionproperties.analysis.section.Section.get_phi
     ~sectionproperties.analysis.section.Section.get_ezp
+    ~sectionproperties.analysis.section.Section.get_my_p
     ~sectionproperties.analysis.section.Section.get_rp
     ~sectionproperties.analysis.section.Section.get_nu_eff
     ~sectionproperties.analysis.section.Section.get_e_eff

docs/user_guide/theory.rst

@@ -420,6 +420,20 @@ extreme (min. and max.) coordinates of the cross-section in the x and y-directio
   Z_{yy}^+ = \frac{I_{\overline{yy}}}{x_{max} - x_c} \\
   Z_{yy}^- = \frac{I_{\overline{yy}}}{x_c - x_{min}} \\
 
+Yield Moments
+~~~~~~~~~~~~~
+
+The yield moment is defined as the lowest bending moment that causes any point within
+cross-section to reach the yield strength. Note that this implementation is purely
+linear-elastic i.e. uses the linear-elastic modulus and bi-directional yield strength
+only.
+
+``sectionproperties`` applies a unit bending moment, about each axis separately, and
+determines the yield index for each point within the mesh. The yield index is defined as
+the stress divided by the material yield strength. Through this method, a critical yield
+index is determined (i.e. point which will yield first under bending) and the yield
+moment calculated as the inverse of the critical yield index.
+
 .. _label-theory-plastic-section-moduli:
 
 Plastic Section Moduli

src/sectionproperties/analysis/section.py

@@ -202,6 +202,7 @@ def calculate_geometric_properties(self) -> None:
           - Centroidal section moduli
           - Radii of gyration
           - Principal axis properties
+          - Yield moments (composite only)
         """
 
         def calculate_geom(progress: Progress | None = None) -> None:
@@ -259,11 +260,87 @@ def calculate_geom(progress: Progress | None = None) -> None:
             )
             self.section_props.e_eff = self.section_props.ea / self.section_props.area
             self.section_props.g_eff = self.section_props.ga / self.section_props.area
+
+            # calculate derived properties
             self.section_props.calculate_elastic_centroid()
             self.section_props.calculate_centroidal_properties(
                 node_list=self.mesh["vertices"]
             )
 
+            # calculate yield moments
+            self.section_props.my_xx = 0.0
+            self.section_props.my_yy = 0.0
+            self.section_props.my_11 = 0.0
+            self.section_props.my_22 = 0.0
+
+            # calculate yield moments:
+            # 1) loop through each material group and through each element in the group
+            # 2) for each point, calculate the bending stress from a unit bending moment
+            # in each direction (mxx, myy, m11, m22)
+            # 3) from this, calculate the yield index for each point
+            # 4) get the largest yield index and scale the bending moment such that the
+            # yield index is 1
+
+            # initialise the yield indexes
+            yield_index = {
+                "mxx": 0.0,
+                "myy": 0.0,
+                "m11": 0.0,
+                "m22": 0.0,
+            }
+
+            # get useful section properties
+            cx = self.section_props.cx
+            cy = self.section_props.cy
+            phi = self.section_props.phi
+            ixx = self.section_props.ixx_c
+            iyy = self.section_props.iyy_c
+            ixy = self.section_props.ixy_c
+            i11 = self.section_props.i11_c
+            i22 = self.section_props.i22_c
+
+            if ixx is None or iyy is None or ixy is None or i11 is None or i22 is None:
+                msg = "Section properties failed to save."
+                raise RuntimeError(msg)
+
+            # loop through each material group
+            for group in self.material_groups:
+                em = group.material.elastic_modulus
+                fy = group.material.yield_strength
+
+                # loop through each element in the material group
+                for el in group.elements:
+                    # loop through each node in the element
+                    for coord in el.coords.transpose():
+                        # calculate coordinates wrt centroidal & principal axes
+                        x = coord[0] - cx
+                        y = coord[1] - cy
+                        x11, y22 = fea.principal_coordinate(phi=phi, x=x, y=y)
+
+                        # calculate bending stresses due to unit moments
+                        sig_mxx = em * (
+                            -ixy / (ixx * iyy - ixy**2) * x
+                            + iyy / (ixx * iyy - ixy**2) * y
+                        )
+                        sig_myy = em * (
+                            -ixx / (ixx * iyy - ixy**2) * x
+                            + ixy / (ixx * iyy - ixy**2) * y
+                        )
+                        sig_m11 = em / i11 * y22
+                        sig_m22 = -em / i22 * x11
+
+                        # update yield indexes
+                        yield_index["mxx"] = max(yield_index["mxx"], abs(sig_mxx / fy))
+                        yield_index["myy"] = max(yield_index["myy"], abs(sig_myy / fy))
+                        yield_index["m11"] = max(yield_index["m11"], abs(sig_m11 / fy))
+                        yield_index["m22"] = max(yield_index["m22"], abs(sig_m22 / fy))
+
+            # calculate yield moments
+            self.section_props.my_xx = 1.0 / yield_index["mxx"]
+            self.section_props.my_yy = 1.0 / yield_index["myy"]
+            self.section_props.my_11 = 1.0 / yield_index["m11"]
+            self.section_props.my_22 = 1.0 / yield_index["m22"]
+
             if progress and task is not None:
                 msg = "[bold green]:white_check_mark: Geometric analysis complete"
                 progress.update(task_id=task, description=msg)
@@ -1120,7 +1197,7 @@ def calculate_plastic_properties(
 
           - Plastic centroids (centroidal and principal axes)
           - Plastic section moduli (centroidal and principal axes)
-          - Shape factors, non-composite only (centroidal and principal axe)
+          - Shape factors, non-composite only (centroidal and principal axes)
         """
         # check that a geometric analysis has been performed
         if self.section_props.cx is None:
@@ -2240,6 +2317,32 @@ def get_ez(
             self.section_props.zyy_minus / e_ref,
         )
 
+    def get_my(self) -> tuple[float, float]:
+        """Returns the yield moment for bending about the centroidal axis.
+
+        This is a composite only property, as such this can only be returned if material
+        properties have been applied to the cross-section.
+
+        Returns:
+            Yield moment for bending about the centroidal ``x`` and ``y`` axes
+            (``my_xx``, ``my_yy``)
+
+        Raises:
+            RuntimeError: If material properties have *not* been applied
+            RuntimeError: If a geometric analysis has not been performed
+        """
+        if not self.is_composite():
+            msg = "Attempting to get a composite only property for a geometric analysis"
+            msg += " (material properties have not been applied). Consider using"
+            msg += " get_z()."
+            raise RuntimeError(msg)
+
+        if self.section_props.my_xx is None or self.section_props.my_yy is None:
+            msg = "Conduct a geometric analysis."
+            raise RuntimeError(msg)
+
+        return (self.section_props.my_xx, self.section_props.my_yy)
+
     def get_rc(self) -> tuple[float, float]:
         """Returns the cross-section centroidal radii of gyration.
 
@@ -2418,6 +2521,32 @@ def get_ezp(
             self.section_props.z22_minus / e_ref,
         )
 
+    def get_my_p(self) -> tuple[float, float]:
+        """Returns the yield moment for bending about the principal axis.
+
+        This is a composite only property, as such this can only be returned if material
+        properties have been applied to the cross-section.
+
+        Returns:
+            Yield moment for bending about the principal ``11`` and ``22`` axes
+            (``my_11``, ``my_22``)
+
+        Raises:
+            RuntimeError: If material properties have *not* been applied
+            RuntimeError: If a geometric analysis has not been performed
+        """
+        if not self.is_composite():
+            msg = "Attempting to get a composite only property for a geometric analysis"
+            msg += " (material properties have not been applied). Consider using"
+            msg += " get_zp()."
+            raise RuntimeError(msg)
+
+        if self.section_props.my_11 is None or self.section_props.my_22 is None:
+            msg = "Conduct a geometric analysis."
+            raise RuntimeError(msg)
+
+        return (self.section_props.my_11, self.section_props.my_22)
+
     def get_rp(self) -> tuple[float, float]:
         """Returns the cross-section principal radii of gyration.
 

src/sectionproperties/post/post.py

@@ -72,6 +72,10 @@ class SectionProperties:
             negative extreme value of the 11-axis
         r11_c: Radius of gyration about the principal 11-axis.
         r22_c: Radius of gyration about the principal 22-axis.
+        my_xx: Yield moment about the x-axis
+        my_yy: Yield moment about the y-axis
+        my_11: Yield moment about the 11-axis
+        my_22: Yield moment about the 22-axis
         j: Torsion constant
         omega: Warping function
         psi_shear: Psi shear function
@@ -165,6 +169,10 @@ class SectionProperties:
     r11_c: float | None = None
     r22_c: float | None = None
     j: float | None = None
+    my_xx: float | None = None
+    my_yy: float | None = None
+    my_11: float | None = None
+    my_22: float | None = None
     omega: npt.NDArray[np.float64] | None = None
     psi_shear: npt.NDArray[np.float64] | None = None
     phi_shear: npt.NDArray[np.float64] | None = None
@@ -278,7 +286,7 @@ def calculate_centroidal_properties(
             else:
                 self.phi = np.arctan2(self.ixx_c - self.i11_c, self.ixy_c) * 180 / np.pi
 
-            # initialise min, max variables TODO: check for `if xxx:` where xxx is float
+            # initialise min, max variables
             if self.phi is not None:
                 x1, y2 = fea.principal_coordinate(
                     phi=self.phi,
@@ -662,6 +670,15 @@ def print_results(
         except RuntimeError:
             pass
 
+    # print cross-section my
+    if is_composite:
+        try:
+            my_xx, my_yy = section.get_my()
+            table.add_row("my_xx", f"{my_xx:>{fmt}}")
+            table.add_row("my_yy-", f"{my_yy:>{fmt}}")
+        except RuntimeError:
+            pass
+
     # print cross-section rc
     try:
         rx, ry = section.get_rc()
@@ -713,6 +730,15 @@ def print_results(
         except RuntimeError:
             pass
 
+    # print cross-section my_p
+    if is_composite:
+        try:
+            my_11, my_22 = section.get_my_p()
+            table.add_row("my_11", f"{my_11:>{fmt}}")
+            table.add_row("my_22-", f"{my_22:>{fmt}}")
+        except RuntimeError:
+            pass
+
     # print cross-section rp
     try:
         r11, r22 = section.get_rp()