diff --git a/src/sectionproperties/analysis/plastic_section.py b/src/sectionproperties/analysis/plastic_section.py index 14cfcc32..61aa6e3e 100644 --- a/src/sectionproperties/analysis/plastic_section.py +++ b/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 diff --git a/tests/analysis/test_plastic.py b/tests/analysis/test_plastic.py index de3fc1c5..ed9b38a4 100644 --- a/tests/analysis/test_plastic.py +++ b/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) diff --git a/tests/validation/test_custom_validation.py b/tests/validation/test_custom_validation.py index 155948b8..07642aa9 100644 --- a/tests/validation/test_custom_validation.py +++ b/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)