Fix plastic calculation for composite sections

src/sectionproperties/analysis/plastic_section.py

@@ -35,13 +35,14 @@ def __init__(
         Args:
             geometry: Section geometry object
         """
+        # move geometry to geometric centroid
         self.geometry = geometry.align_center()
         self.geometry.compile_geometry()
 
         # initialize variables to be defined later within calculate_plastic_force
-        self._c_top = [0.0, 0.0]
-        self._c_bot = [0.0, 0.0]
-        self._f_top = 0.0
+        self._f_list: list[float] = []
+        self._cx_list: list[float] = []
+        self._cy_list: list[float] = []
 
     def calculate_plastic_properties(
         self,
@@ -89,7 +90,15 @@ def calculate_plastic_properties(
 
         self.check_convergence(root_result=r, axis="x-axis")
         section.section_props.y_pc = y_pc
-        section.section_props.sxx = self._f_top * abs(self._c_top[1] - self._c_bot[1])
+
+        # calculate plastic section moduli (plastic moment if composite)
+        section.section_props.sxx = 0.0
+
+        # loop through each geometry force & y-centroid
+        for f, c in zip(self._f_list, self._cy_list):
+            # calculate distance from y-centroid to plastic centroid
+            d_y = abs(c - y_pc)
+            section.section_props.sxx += f * d_y
 
         if verbose:
             self.print_verbose(d=y_pc, root_result=r, axis="x-axis")
@@ -107,7 +116,15 @@ def calculate_plastic_properties(
 
         self.check_convergence(root_result=r, axis="y-axis")
         section.section_props.x_pc = x_pc
-        section.section_props.syy = self._f_top * abs(self._c_top[0] - self._c_bot[0])
+
+        # calculate plastic section moduli (plastic moment if composite)
+        section.section_props.syy = 0.0
+
+        # loop through each geometry force & x-centroid
+        for f, c in zip(self._f_list, self._cx_list):
+            # calculate distance from x-centroid to plastic centroid
+            d_x = abs(c - x_pc)
+            section.section_props.syy += f * d_x
 
         if verbose:
             self.print_verbose(d=x_pc, root_result=r, axis="y-axis")
@@ -137,17 +154,20 @@ def calculate_plastic_properties(
             verbose=verbose,
         )
 
-        # calculate the centroids in the principal coordinate system
-        c_top_p = fea.principal_coordinate(
-            phi=section.section_props.phi, x=self._c_top[0], y=self._c_top[1]
-        )
-        c_bot_p = fea.principal_coordinate(
-            phi=section.section_props.phi, x=self._c_bot[0], y=self._c_bot[1]
-        )
-
         self.check_convergence(root_result=r, axis="11-axis")
         section.section_props.y22_pc = y22_pc
-        section.section_props.s11 = self._f_top * abs(c_top_p[1] - c_bot_p[1])
+
+        # calculate plastic section moduli (plastic moment if composite)
+        section.section_props.s11 = 0.0
+
+        # loop through each geometry force & xy-centroid
+        for f, cx, cy in zip(self._f_list, self._cx_list, self._cy_list):
+            # convert centroid to principal coordinates
+            cen = fea.principal_coordinate(phi=section.section_props.phi, x=cx, y=cy)
+
+            # calculate distance from 22-centroid to plastic centroid
+            d22 = abs(cen[1] - y22_pc)
+            section.section_props.s11 += f * d22
 
         if verbose:
             self.print_verbose(d=y22_pc, root_result=r, axis="11-axis")
@@ -163,17 +183,20 @@ def calculate_plastic_properties(
             verbose=verbose,
         )
 
-        # calculate the centroids in the principal coordinate system
-        c_top_p = fea.principal_coordinate(
-            phi=section.section_props.phi, x=self._c_top[0], y=self._c_top[1]
-        )
-        c_bot_p = fea.principal_coordinate(
-            phi=section.section_props.phi, x=self._c_bot[0], y=self._c_bot[1]
-        )
-
         self.check_convergence(root_result=r, axis="22-axis")
         section.section_props.x11_pc = x11_pc
-        section.section_props.s22 = self._f_top * abs(c_top_p[0] - c_bot_p[0])
+
+        # calculate plastic section moduli (plastic moment if composite)
+        section.section_props.s22 = 0.0
+
+        # loop through each geometry force & xy-centroid
+        for f, cx, cy in zip(self._f_list, self._cx_list, self._cy_list):
+            # convert centroid to principal coordinates
+            cen = fea.principal_coordinate(phi=section.section_props.phi, x=cx, y=cy)
+
+            # calculate distance from 11-centroid to plastic centroid
+            d11 = abs(cen[0] - x11_pc)
+            section.section_props.s22 += f * d11
 
         if verbose:
             self.print_verbose(d=x11_pc, root_result=r, axis="22-axis")
@@ -396,13 +419,13 @@ def calculate_plastic_force(
             p: Point on the axis line
 
         Returns:
-            Force in the above and below the axis line (``f_top``, ``f_bot``)
+            Force in the geometries above and below the axis line (``f_top``, ``f_bot``)
         """
         # initialise variables
         f_top, f_bot = 0.0, 0.0
-        a_top, a_bot = 0.0, 0.0
-        qx_top, qx_bot = 0.0, 0.0
-        qy_top, qy_bot = 0.0, 0.0
+        self._f_list = []
+        self._cx_list = []
+        self._cy_list = []
 
         # split geometry above and below the line
         top_geoms, bot_geoms = self.geometry.split_section(point_i=p, vector=u)
@@ -415,12 +438,14 @@ def calculate_plastic_force(
                 area_top = top_geom.calculate_area()
                 cx, cy = top_geom.calculate_centroid()
 
-                # sum properties
-                a_top += area_top
-                qx_top += cy * area_top
-                qy_top += cx * area_top
+                # sum force
                 f_top += f_y * area_top
 
+                # add force and centroids to list
+                self._f_list.append(f_y * area_top)
+                self._cx_list.append(cx)
+                self._cy_list.append(cy)
+
         if bot_geoms:
             # loop through all bottom geometries
             for bot_geom in bot_geoms:
@@ -429,23 +454,12 @@ def calculate_plastic_force(
                 area_bot = bot_geom.calculate_area()
                 cx, cy = bot_geom.calculate_centroid()
 
-                # sum properties
-                a_bot += area_bot
-                qx_bot += cy * area_bot
-                qy_bot += cx * area_bot
+                # sum force
                 f_bot += f_y * area_bot
 
-        # calculate centroid of force action
-        try:
-            self._c_top = [qy_top / a_top, qx_top / a_top]
-            self._f_top = f_top
-        except ZeroDivisionError:
-            self._c_top = [0.0, 0.0]
-            self._f_top = 0.0
-
-        try:
-            self._c_bot = [qy_bot / a_bot, qx_bot / a_bot]
-        except ZeroDivisionError:
-            self._c_bot = [0.0, 0.0]
+                # add force and centroids to list
+                self._f_list.append(f_y * area_bot)
+                self._cx_list.append(cx)
+                self._cy_list.append(cy)
 
         return f_top, f_bot

tests/analysis/test_plastic.py

@@ -93,7 +93,7 @@ def test_plastic_centroid():
     # create a Section object
     section = Section(geometry=geometry)
 
-    # perform a geometric, warping and plastic analysis
+    # perform a geometric and plastic analysis
     section.calculate_geometric_properties()
     section.calculate_plastic_properties()
 
@@ -125,3 +125,152 @@ def test_plastic_centroid():
     x_pc, y_pc = nested_sec.get_pc()
     assert x_pc == pytest.approx(37.5)
     assert y_pc == pytest.approx(50)
+
+
+def test_plastic_composite_simple():
+    """Tests the plastic properties of a simple composite section."""
+    mat1 = Material(
+        name="mat1",
+        elastic_modulus=10,
+        poissons_ratio=0.3,
+        density=1,
+        yield_strength=50,
+        color="grey",
+    )
+
+    mat2 = Material(
+        name="mat2",
+        elastic_modulus=20,
+        poissons_ratio=0.3,
+        density=1,
+        yield_strength=80,
+        color="black",
+    )
+
+    rect1 = sections.rectangular_section(d=100, b=50, material=mat1)
+    rect2 = sections.rectangular_section(d=100, b=50, material=mat2).shift_section(
+        x_offset=50
+    )
+    geom = rect1 + rect2
+    geom.create_mesh(mesh_sizes=100)
+
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+    sec.calculate_plastic_properties()
+
+    x_pc, y_pc = sec.get_pc()
+    assert y_pc == pytest.approx(50)
+
+    # calculate x_pc
+    f_mat1 = 100 * 50 * 50  # force in mat1
+    x_mat2 = f_mat1 / 80 / 100  # width of mat2 required to equal force in mat1
+    x_left = (50 - x_mat2) / 2  # width of ma2 to left of pc
+    assert x_pc == pytest.approx(x_left + 50)
+
+    mp_xx, mp_yy = sec.get_mp()
+    mp_xx_calc = 50 * 100**2 / 4 * (80 + 50)
+    assert mp_xx == pytest.approx(mp_xx_calc)
+
+    # calculate mp_yy
+    mp_yy_calc = 0
+    mp_yy_calc += f_mat1 * abs(25 - x_pc)  # mat 1 contribution
+    # left mat 2 contribution
+    mp_yy_calc += x_left * 100 * 80 * abs(50 + x_left / 2 - x_pc)
+    # right mat 2 contribution
+    mp_yy_calc += (50 - x_left) * 100 * 80 * abs(100 - (50 - x_left) / 2 - x_pc)
+    assert mp_yy == pytest.approx(mp_yy_calc)
+
+    # check principal calc with rotation 45 deg
+    geom.rotate_section(angle=45)
+    geom.create_mesh(mesh_sizes=100)
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+    sec.calculate_plastic_properties()
+
+    mp_11, mp_22 = sec.get_mp_p()
+    assert mp_11 == pytest.approx(mp_xx_calc)
+    assert mp_22 == pytest.approx(mp_yy_calc)
+
+
+def test_plastic_composite_example():
+    """Tests the composite example on the sectionproperties docs, refer issue #460.
+
+    The steel section is simplified to have no root radii to allow hand calculation
+    comparison.
+    """
+    # create the steel material
+    steel = Material(
+        name="Steel",
+        elastic_modulus=200e3,
+        poissons_ratio=0.3,
+        density=7.85e-6,
+        yield_strength=500,
+        color="grey",
+    )
+
+    # create the timber material
+    timber = Material(
+        name="Timber",
+        elastic_modulus=8e3,
+        poissons_ratio=0.35,
+        yield_strength=20,
+        density=0.78e-6,
+        color="burlywood",
+    )
+
+    # universal steel beam
+    ub = steel_sections.i_section(
+        d=304, b=165, t_f=10.2, t_w=6.1, r=0, n_r=2, material=steel
+    )
+
+    # timber floor panel
+    panel = sections.rectangular_section(d=100, b=600, material=timber)
+    panel = panel.align_center(align_to=ub).align_to(other=ub, on="top")
+
+    # combine geometry
+    geom = ub + panel
+
+    geom.create_mesh(mesh_sizes=[50, 500])
+    sec = Section(geometry=geom)
+    sec.calculate_geometric_properties()
+    sec.calculate_plastic_properties()
+
+    # get calculated values
+    x_pc, y_pc = sec.get_pc()
+    mp_xx, mp_yy = sec.get_mp()
+
+    # hand calcs
+    # 1) x-axis bending
+    # calculate y_pc -> assume centroid in top flange
+    f_timber = 600 * 100 * 20
+    f_web = (304 - 2 * 10.2) * 6.1 * 500
+    f_flange = 165 * 10.2 * 500
+
+    # let y = distance from bottom of top flange to y_pc
+    # A) bot = 165 * y * 500
+    # B) top = 165 * (10.2 - y) * 500
+    # C) (A) - (B) = f_timber - f_web - f_flange  = f_diff
+    # 165 * 500 (2y - 10.2) = f_diff
+    # y = (f_diff / (165 * 500) + 10.2) / 2
+    f_diff = f_timber - f_web - f_flange
+    y = (f_diff / (165 * 500) + 10.2) / 2
+    assert y_pc == pytest.approx(304 - 10.2 + y)
+
+    # calculate mp_xx
+    mp_xx_calc = 0
+    mp_xx_calc += f_timber * abs(304 + 50 - y_pc)  # timber contribution
+    mp_xx_calc += (10.2 - y) * 165 * 500 * abs((10.2 - y) / 2)  # top flange p1
+    mp_xx_calc += y * 165 * 500 * abs(y / 2)  # top flange p2
+    mp_xx_calc += f_web * abs(10.2 + (304 - 2 * 10.2) / 2 - y_pc)  # web contribution
+    mp_xx_calc += f_flange * abs(10.2 / 2 - y_pc)  # bot flange contribution
+    assert mp_xx == pytest.approx(mp_xx_calc)
+
+    # 2) y-axis bending
+    assert x_pc == pytest.approx(165 / 2)
+
+    # calculate mp_yy
+    mp_yy_calc = 0
+    mp_yy_calc += 2 * (10.2 * 165**2 / 4) * 500  # 2 flanges
+    mp_yy_calc += (304 - 10.2 * 2) * 6.1**2 / 4 * 500  # web
+    mp_yy_calc += 100 * 600**2 / 4 * 20  # timber
+    assert mp_yy == pytest.approx(mp_yy_calc)

tests/validation/test_custom_validation.py

@@ -12,6 +12,7 @@
 
 # constants
 tol = 1e-6
+plastic_tol = 1e-5
 warp_tol = 1e-3
 
 
@@ -122,7 +123,7 @@ def test_custom_section_warping(custom_section):
 def test_custom_section_plastic(custom_section):
     """Test custom section plastic properties.
 
-    Warping properties calculated from sectionproperties v3.0.2 with a refined mesh
+    Plastic properties calculated from sectionproperties v3.0.2 with a refined mesh
     [mesh_sizes=0.5].
     """
     sec = custom_section
@@ -131,19 +132,19 @@ def test_custom_section_plastic(custom_section):
     x_pc, y_pc = sec.get_pc()
     x11_pc, y22_pc = sec.get_pc_p()
 
-    check.almost_equal(x_pc, 4.977273e01, rel=tol)
-    check.almost_equal(y_pc, 9.172040e01, rel=tol)
-    check.almost_equal(x11_pc, 5.133714e01, rel=tol)
-    check.almost_equal(y22_pc, 9.158984e01, rel=tol)
-    check.almost_equal(sec.section_props.sxx, 1.531971e05, rel=tol)
-    check.almost_equal(sec.section_props.syy, 1.014943e05, rel=tol)
-    check.almost_equal(sec.section_props.s11, 1.533463e05, rel=tol)
-    check.almost_equal(sec.section_props.s22, 1.015010e05, rel=tol)
-    check.almost_equal(sec.section_props.sf_xx_plus, 8.942884e-01, rel=tol)
-    check.almost_equal(sec.section_props.sf_xx_minus, 1.292703e00, rel=tol)
-    check.almost_equal(sec.section_props.sf_yy_plus, 1.602519e00, rel=tol)
-    check.almost_equal(sec.section_props.sf_yy_minus, 1.567298e00, rel=tol)
-    check.almost_equal(sec.section_props.sf_11_plus, 9.450478e-01, rel=tol)
-    check.almost_equal(sec.section_props.sf_11_minus, 1.341988e00, rel=tol)
-    check.almost_equal(sec.section_props.sf_22_plus, 1.677621e00, rel=tol)
-    check.almost_equal(sec.section_props.sf_22_minus, 1.619711e00, rel=tol)
+    check.almost_equal(x_pc, 4.977273e01, rel=plastic_tol)
+    check.almost_equal(y_pc, 9.172040e01, rel=plastic_tol)
+    check.almost_equal(x11_pc, 5.133714e01, rel=plastic_tol)
+    check.almost_equal(y22_pc, 9.158984e01, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sxx, 1.531971e05, rel=plastic_tol)
+    check.almost_equal(sec.section_props.syy, 1.014943e05, rel=plastic_tol)
+    check.almost_equal(sec.section_props.s11, 1.533463e05, rel=plastic_tol)
+    check.almost_equal(sec.section_props.s22, 1.015010e05, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sf_xx_plus, 8.942884e-01, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sf_xx_minus, 1.292703e00, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sf_yy_plus, 1.602519e00, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sf_yy_minus, 1.567298e00, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sf_11_plus, 9.450478e-01, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sf_11_minus, 1.341988e00, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sf_22_plus, 1.677621e00, rel=plastic_tol)
+    check.almost_equal(sec.section_props.sf_22_minus, 1.619711e00, rel=plastic_tol)