Error in computation of bearing capacity for axial loading

[DanielGMorenaFhecor]

When computing the bearing capacity of a cross-section with an axial load N ≠ 0, the results are incorrect at first instance. Comparing these results with established software reveals discrepancies. It works correctly when the center of gravity is at the origin.

See Point 6 in this notebook developed by @MestreCarlos for more information and illustrations.

[talledodiego]

I still need to dig into it and I will do as soon as I am back from vacation. From what you say it works when center of gravity is at the origin. This is an intended behavior: the moments are computed integrating the normal stresses respect to global Y and Z axis. When you have no external axial force it is independent on the coordinates of the center of gravity. When there is external axial force there is also the offset moment due to N * e. It is like having a beam element with a rigid offset at the end.

This is why in our examples generally we translate the section like in this example:

# 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
# because it changes when there is axial force.
geo = geo.translate(-175, -250)
sec = GenericSection(geo)

sec.geometry

In some cases the user may want to account for this rigid offset.

[DanielGMorenaFhecor]

Thanks for your response! I've updated the example to translate the section as you suggested. I'll leave the updated example (Point 6 in this new update) and the results here, and we can discuss them further when you're back.

These results have been updated with your latest fix to Marin's algorithm in PR #133.

# Materials
concrete = ConcreteEC2_2004(25)
reinforcement = 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, 50), (300, 50), 12, reinforcement, n=3)
geo = add_reinforcement_line(geo, (50, 450), (300, 450), 20, reinforcement, n=6)
geo = geo.translate(-175, -250)  # Translate the section
sec = GenericSection(geo)

# Compute bending strength for a list of axial loads
for n_ed in [-3500, -3000, -2500, -2000, -1500, -1000, -500, 0, 500, 900]:
    m_y = sec.section_calculator.calculate_bending_strength(n=n_ed*1e3).m_y / 1e6
    print(f"N={n_ed:.2f} N\t M={m_y:.2f} kNm")

Output:

N=-3500.00 N  M=-207.35 kNm
N=-3000.00 N  M=-289.86 kNm
N=-2500.00 N  M=-343.26 kNm
N=-2000.00 N  M=-369.23 kNm
N=-1500.00 N  M=-340.02 kNm
N=-1000.00 N  M=-262.75 kNm
N=-500.00 N   M=-166.01 kNm
N=0.00 N      M=-65.78 kNm
N=500.00 N    M=35.05 kNm
N=900.00 N    M=118.06 kNm

To summarize, here’s the image with the results from the software we're using (check the table in the middle of the image):