Implement yield moment, update theory, add to display_results(), add more tests

Author: robbievanleeuwen

docs/user_guide/theory.rst

@@ -423,7 +423,16 @@ extreme (min. and max.) coordinates of the cross-section in the x and y-directio
 Yield Moments
 ~~~~~~~~~~~~~
 
-TODO
+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:
 

src/sectionproperties/analysis/section.py

@@ -273,7 +273,73 @@ def calculate_geom(progress: Progress | None = None) -> None:
             self.section_props.my_11 = 0.0
             self.section_props.my_22 = 0.0
 
-            # TODO: calculate yield moments
+            # 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"

src/sectionproperties/post/post.py

@@ -670,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()
@@ -721,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()

tests/analysis/test_yield_moment.py

@@ -1,10 +1,12 @@
 """Tests for the yield moment calculation."""
 
+import numpy as np
 import pytest
 
 from sectionproperties.analysis import Section
+from sectionproperties.post.stress_post import StressPost
 from sectionproperties.pre import Material
-from sectionproperties.pre.library import rectangular_section
+from sectionproperties.pre.library import i_section, rectangular_section
 
 STEEL = Material(
     name="Steel",
@@ -18,7 +20,7 @@
 
 def test_get_without_analysis():
     """Test for raising an error if a geometric analysis has not been performed."""
-    geom = rectangular_section(d=50, b=50)
+    geom = rectangular_section(d=50, b=50, material=STEEL)
     geom.create_mesh(mesh_sizes=0, coarse=True)
     sec = Section(geometry=geom)
 
@@ -52,14 +54,195 @@ def test_rectangle():
 
     my_xx, my_yy = sec.get_my()
     my_11, my_22 = sec.get_my()
+    my_xx_calc = 50 * 100 * 100 / 6 * 500
+    my_yy_calc = 100 * 50 * 50 / 6 * 500
+
+    # compare with theoretical values
+    assert my_xx == pytest.approx(my_xx_calc)
+    assert my_yy == pytest.approx(my_yy_calc)
+    assert my_11 == pytest.approx(my_xx_calc)
+    assert my_22 == pytest.approx(my_yy_calc)
+
+    # compare with stress analysis
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def test_rectangle_rotated():
+    """Test the yield moment of a simple rotated rectangle."""
+    geom = rectangular_section(d=100, b=50, material=STEEL)
+    geom.rotate_section(angle=45.0)
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+
+    my_11, my_22 = sec.get_my()
+    my_xx_calc = 50 * 100 * 100 / 6 * 500
+    my_yy_calc = 100 * 50 * 50 / 6 * 500
+
+    # compare with theoretical values
+    assert my_11 == pytest.approx(my_xx_calc)
+    assert my_22 == pytest.approx(my_yy_calc)
+
+    # compare with stress analysis
+    stress = sec.calculate_stress(m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def test_isection():
+    """Test the yield moment of an isection."""
+    geom = i_section(d=200, b=100, t_f=10, t_w=5, r=12, n_r=8, material=STEEL)
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    ezxx, _, ezyy, _ = sec.get_ez()
+    my_xx_calc = ezxx / 200e3 * 500
+    my_yy_calc = ezyy / 200e3 * 500
 
-    my_xx_calc = 50 * 100 * 100 / 4 * 500
-    my_yy_calc = 100 * 50 * 50 / 4 * 500
+    # compare with theoretical values
+    assert my_xx == pytest.approx(my_xx_calc)
+    assert my_yy == pytest.approx(my_yy_calc)
+    assert my_11 == pytest.approx(my_xx_calc)
+    assert my_22 == pytest.approx(my_yy_calc)
+
+    # compare with stress analysis
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def test_rectangle_composite():
+    """Test the yield moment of a composite rectangular section."""
+    mat1 = Material("a", 2, 0, 2, 1, "b")
+    mat2 = Material("b", 1, 0, 1, 1, "r")
 
+    rect1 = rectangular_section(d=20, b=20, material=mat1)
+    rect2 = rectangular_section(d=20, b=20, material=mat2).align_to(rect1, "top")
+    rect3 = rectangular_section(d=20, b=20, material=mat1).align_to(rect2, "top")
+    geom = rect1 + rect2 + rect3
+    geom.create_mesh(mesh_sizes=0, coarse=True)
+    sec = Section(geom)
+    sec.calculate_geometric_properties()
+
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    eixx_calc = (20 * 20**3 / 12) + 2 * 2 * ((20 * 20**3 / 12) + (20 * 20 * 20**2))
+    eiyy_calc = 5 * (20 * 20**3 / 12)
+    # myield = f * Ieff / y
+    my_xx_calc = 2 * (eixx_calc / 2) / 30
+    my_yy_calc = 1 * eiyy_calc / 10
+
+    # compare with theoretical values
     assert my_xx == pytest.approx(my_xx_calc)
     assert my_yy == pytest.approx(my_yy_calc)
     assert my_11 == pytest.approx(my_xx_calc)
     assert my_22 == pytest.approx(my_yy_calc)
 
+    # compare with stress analysis
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def test_composite_example():
+    """Tests the composite example from the docs, compares to stress analysis."""
+    timber = Material(
+        name="Timber",
+        elastic_modulus=8e3,
+        poissons_ratio=0.35,
+        yield_strength=20,
+        density=0.78e-6,
+        color="burlywood",
+    )
+
+    ub = i_section(d=304, b=165, t_f=10.2, t_w=6.1, r=11.4, n_r=8, material=STEEL)
+    panel = rectangular_section(d=100, b=600, material=timber)
+    panel = panel.align_center(align_to=ub).align_to(other=ub, on="top")
+    geom = ub + panel
+    geom.create_mesh(mesh_sizes=[10, 500])
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+
+    # compare with stress analysis
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def test_yield_internal():
+    """Test where the internal material yields first."""
+    mat1 = Material(
+        name="Material_1",
+        elastic_modulus=0.75,
+        poissons_ratio=0,
+        yield_strength=1,
+        density=1,
+        color="gold",
+    )
+    mat2 = Material(
+        name="Material_2",
+        elastic_modulus=1,
+        poissons_ratio=0,
+        yield_strength=1,
+        density=1,
+        color="blue",
+    )
+    mat3 = Material(
+        name="Material 3",
+        elastic_modulus=3,
+        poissons_ratio=0,
+        yield_strength=1,
+        density=1,
+        color="red",
+    )
+
+    sq1 = rectangular_section(d=100, b=100, material=mat1).align_center()
+    sq2 = rectangular_section(d=75, b=75, material=mat2).align_center()
+    sq3 = rectangular_section(d=50, b=50, material=mat3).align_center()
+    hole = rectangular_section(d=25, b=25).align_center()
+    compound = (sq1 - sq2) + (sq2 - sq3) + (sq3 - hole)
+    compound.create_mesh(10)
+    sec = Section(compound)
+    sec.calculate_geometric_properties()
+
+    # compare with stress analysis
+    my_xx, my_yy = sec.get_my()
+    my_11, my_22 = sec.get_my()
+    stress = sec.calculate_stress(mxx=my_xx, myy=my_yy, m11=my_11, m22=my_22)
+    assert check_yield_index(stress, "sig_zz_mxx") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_myy") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m11") == pytest.approx(1.0)
+    assert check_yield_index(stress, "sig_zz_m22") == pytest.approx(1.0)
+
+
+def check_yield_index(stress: StressPost, key: str) -> float:
+    """Returns the largest yield index.
+
+    Given the output from StressPost object and a dict key representing the type of
+    stress, returns the largest yield index.
+    """
+    yield_index = 0.0
+
+    for idx, mat_dict in enumerate(stress.get_stress()):
+        fy = stress.material_groups[idx].material.yield_strength
+        sigs = mat_dict[key]
+        sig_max = max(np.absolute(sigs))
+
+        yield_index = max(yield_index, sig_max / fy)
 
-# TODO: add more tests
+    return yield_index