Inconsistent Bending Strength Calculation for Rotated Sections

[DanielGMorenaFhecor]

We observed discrepancies in the bending strength results when using the calculate_bending_strength method in a GenericSection instance under different rotation angles. Specifically, we compared the results for the same cross-section at angles theta=0 and theta=π/2. The output of calculate_bending_strength(theta=0).m_z does not match calculate_bending_strength(theta=π/2).m_y, as expected.

Additionally, when the section is rotated by π/2, the calculated bending strength m_y for the rotated section do not match as it should with the non-rotated m_z.

Please refer to the example in Point 4 provided by @MestreCarlos in this file for further details.

Code Example:

# Materials
concrete = ConcreteEC2_2004(25)
reinforcement = ReinforcementEC2_2023(fyk=500, Es=200000, density=7850, ftk=550, epsuk=0.07)

# Create section
poly = Polygon(((0, 0), (350, 0), (350, 500), (0, 500)))
geo = SurfaceGeometry(poly, concrete)
geo = add_reinforcement_line(geo, (50, 450), (300, 450), 20, reinforcement, n=10)
geo = add_reinforcement_line(geo, (50, 50), (300, 50), 20, reinforcement, n=10)
sec = GenericSection(geo, integrator='marin')

# Calculate and display bending strengths
print(f"Mz_max [mkN]={round(sec.section_calculator.calculate_bending_strength(theta=0).m_z/1000**2)}")
print(f"Mz_max [mkN]={round(sec.section_calculator.calculate_bending_strength(theta=math.pi/2).m_y/1000**2)}")

# Rotate the section and recalculate
sec.geometry = sec.geometry.rotate(math.pi/2)
print(f"My_max_rotated [mkN]={round(sec.section_calculator.calculate_bending_strength(theta=0).m_y/1000**2)}")

[DanielGMorenaFhecor]

I've updated the file to provide additional clarification.

Here's a clearer example: It appears that when using Marin's algorithm to integrate the bending moment contributions with respect to the Y-axis (Y-axis as defined in the _marin_integration.py file), some terms may be canceling each other out, leading to (almost) zero values. Is this the expected behavior?

# Materials
concrete = ConcreteEC2_2004(25)
reinforcement = ReinforcementEC2_2023(fyk=500, Es=200000, density=7850, ftk=550, epsuk=0.07)

# Create section
poly = Polygon(((0, 0), (350, 0), (350, 500), (0, 500)))
geo = SurfaceGeometry(poly, concrete)
geo = add_reinforcement_line(geo, (50, 450), (300, 450), 20, reinforcement, n=10)
geo = add_reinforcement_line(geo, (50, 50), (300, 50), 20, reinforcement, n=10)
geo = geo.translate(-175,-250)
sec = GenericSection(geo)

# Draw & print
print("====== SECTION ======")
print("====== Theta = 0 ======")
pprint(sec.section_calculator.calculate_bending_strength(theta=0))
print("====== Theta = pi/2 ======")
pprint(sec.section_calculator.calculate_bending_strength(theta=math.pi/2))
draw_section(sec)

====== SECTION ======
====== Theta = 0 ======
UltimateBendingMomentResults(theta=0,
                             n=-0.0029919976950623095,
                             m_y=-562853427.1074414,
                             m_z=4.4796883128622224e-08, # ALMOST ZERO RESULT
                             chi_y=-3.909551620178221e-05,
                             chi_z=0,
                             eps_a=0.00627387905044555)
====== Theta = pi/2 ======
UltimateBendingMomentResults(theta=1.5707963267948966,
                             n=-0.004536689957603812,
                             m_y=-4.6863929542362246e-08, # ALMOST ZERO RESULT
                             m_z=-255115350.12920046,
                             chi_y=-2.689176391228454e-05,
                             chi_z=0,
                             eps_a=0.001206058684649792)

# Rotated section
rot_sec = GenericSection(geo.rotate(math.pi/2))
print("====== ROTATED SECTION ======")
print("====== Theta = 0 ======")
pprint(rot_sec.section_calculator.calculate_bending_strength(theta=0))
print("====== Theta = pi/2 ======")
pprint(rot_sec.section_calculator.calculate_bending_strength(theta=math.pi/2))
draw_section(rot_sec)

====== ROTATED SECTION ======
====== Theta = 0 ======
UltimateBendingMomentResults(theta=0,
                             n=-0.004536689957603812,
                             m_y=-255115350.12920046,
                             m_z=2.751732939644625e-08, # ALMOST ZERO RESULT
                             chi_y=-2.689176391228454e-05,
                             chi_z=0,
                             eps_a=0.001206058684649792)
====== Theta = pi/2 ======
UltimateBendingMomentResults(theta=1.5707963267948966,
                             n=-0.0029919976950623095,
                             m_y=-7.926171552343454e-08, # ALMOST ZERO RESULT
                             m_z=-562853427.1074414,
                             chi_y=-3.909551620178221e-05,
                             chi_z=0,
                             eps_a=0.00627387905044555)

We are unsure if we are missing something or if this is an issue with how we are using the package. Any insights would be greatly appreciated. Thanks in advance!

[talledodiego]

Thanks @DanielGMorenaFhecor! I will look at the new example you prepared when I go back from vacation. Before leaving, I tested the Point 4 of @MestreCarlos and is working as expected, as it can be seen from the following code (I copy-paste the cells and the relative outputs in the following).

# imports
from math import pi
from shapely import Polygon
from structuralcodes.materials.concrete import ConcreteEC2_2004
from structuralcodes.materials.reinforcement import ReinforcementEC2_2004
from structuralcodes.geometry import SurfaceGeometry
from structuralcodes.sections import add_reinforcement_line, GenericSection

# Materials
concrete = ConcreteEC2_2004(25)
reinforcemnet = ReinforcementEC2_2004(fyk=500, Es=200000, density=7850, ftk=550, epsuk=0.07)  
# Create section
poly = Polygon(((0, 0), (350, 0), (350, 500), (0, 500)))
geo = SurfaceGeometry(poly, concrete)
geo = add_reinforcement_line(geo, (50, 450), (300, 450), 20, reinforcemnet, n=10)
geo = add_reinforcement_line(geo, (50, 50), (300, 50), 20, reinforcemnet, n=10)
# Without axial force this is actually not needed, but I am doing anyway
geo = geo.translate(-175, -250)
sec = GenericSection(geo)

sec.geometry

res = sec.section_calculator.calculate_bending_strength(theta=0)
print(f"My = {res.m_y * 1e-6:.2f} kNm")
print(f"Mz = {res.m_z * 1e-6:.2f} kNm")

This returns:

My = -562.85 kNm Mz = 0.00 kNm

We rotate the geometry:

sec.geometry = sec.geometry.rotate(pi / 2)
sec.geometry

obtaining the following:

The bending moment strength is now:

res = sec.section_calculator.calculate_bending_strength(theta=pi/2)
print(f"My = {res.m_y * 1e-6:.2f} kNm")
print(f"Mz = {res.m_z * 1e-6:.2f} kNm")

That returns as expected:

My = -0.00 kNm Mz = -562.85 kNm

[DanielGMorenaFhecor]

Thank you, @talledodiego, for taking the time to address these issues, especially while you're on vacation!

For example, when we use theta=0, we're seeing My=-562 kNm, which seems reasonable. However, our primary concern is that Mz is coming out as 0. Shouldn't we be obtaining the bending strength with respect to the Z axis (which, based on our understanding of the coordinate system, is upwards)?

In our example, we can achieve this by rotating the section and re-evaluating with theta=0, which then gives a value of -255 kNm (as shown in the results from my last comment).

Also, the updated example includes some fixes that you applied to Marin's integration algorithm in PR #133 .

It’s possible that we might be misinterpreting the results on our end. Any further clarification you can provide would be greatly appreciated, but please, no rush in answering us while you are on vacation 😉​!