How to find the theta angle (neutral axis angle) given external actions (N, Mx, My)

[ocabitza]

Hi, I’m trying to create a Python program to automate ultimate capacity checks given a certain section, reinforcement, and materials.
From the analysis software I get triplets of axial force, bending moment x, and bending moment y (N, Mx, My). For each triplet I’d like to know if there is a way to directly find the corresponding neutral axis angle θ, so that I can then proceed with the calculation of resisting Mx and My.

So far I’ve found a method that works, but it’s too slow. To summarize:

• given the axial force I obtain a biaxial bending diagram Mx, My;
• for each Mx, My pair of the diagram I compute the corresponding angle α in the Mx–My plane, formed by the unit vector in the Mx direction and the vector from the origin to the point Mx, My. This gives me an array of α values;
• these α have a correspondence to the array of θ values (neutral axis angles) used to evaluate the resisting bending moments;
• by interpolation, using as input the α value corresponding to the external Mx and My, I obtain the related θ value;
• with θ known, I can then compute the resistant values.

This method is a bit laborious but it works. The problem, as I mentioned, is the slowness. I would need to know if there is a way to calculate θ directly.

[github-actions]

Thanks for opening your first issue in concreteproperties 🙌 Pull requests are always welcome 😉

[Agent6-6-6]

@ocabitza, that's really the only way to get a more exact capacity. Of course it won't be truly exact as you're still interpolating for the neutral axis angle, which depending on the cross section shape may not be non-linear compared to the corresponding axial load and load angle. It's only as good as the resolution of the data points you're interpolating between.

The way I've done it previously is decrease the neutral axis angle increment in a biaxial analysis, this gives a more refined Mx/My relationship with more points around the interaction curve. Then just interpolate for the moment capacities based on your load angle from the applied Mx/My moments. Of course this is just a capacity ratio in the horizontal Mx/My plane at a fixed axial load. But saves a number of your analysis steps which don't appear to add any further refinement to the answer (certainly not for everyday design).

Usually you're potentially just interested if your N/Mx/My point is inside the 3D M/N surface. Job done. A ratio in the Mx/My plane give an idea of how close you are to the limiting capacity.

Ideally code is modified at some point to include a method to determine the resulting neutral axis angle from a load angle, axial load. Then knowing neutral axis angle back calculate exact capacities. It's been discussed ages ago in an issue from memory and identified as something that needed to be actioned.

[ocabitza]

Ideally code is modified at some point to include a method to determine the resulting neutral axis angle from a load angle, axial load. Then knowing neutral axis angle back calculate exact capacities. It's been discussed ages ago in an issue from memory and identified as something that needed to be actioned.

@Agent6-6-6, thank you for your reply. I think in the meantime I have created a script that does something similar.
I have created an iterative procedure that refines the points around the initial α (angle in Mx, My plane) in an iterative way. Convergence is achieved when the α of the pair of resistant moments is sufficiently close to the input value (based on applied Mx, My). In this way, the obtained θ value should be correct. The issue of the script's slowness still remains, but at least the accuracy no longer depends on the number of values used to build the biaxial moment domain.

[Agent6-6-6]

Perhaps share what you have and someone may have some ideas on how to speed things up.

[ocabitza]

@Agent6-6-6 Sorry for the delay, here is my code.

The first part of the code (section and material definition) is taken from the Ultimate Bending Capacity example

https://concrete-properties.readthedocs.io/en/stable/examples/ultimate_bending.html

import numpy as np
import matplotlib.pyplot as plt
from sectionproperties.pre.library import rectangular_section, triangular_section

import concreteproperties.stress_strain_profile as ssp
from concreteproperties import (
    Concrete,
    ConcreteSection,
    SteelBar,
    add_bar_rectangular_array,
)
from concreteproperties.post import si_kn_m, si_n_mm

concrete = Concrete(
    name="50 MPa Concrete",
    density=2.4e-6,
    stress_strain_profile=ssp.ConcreteLinear(elastic_modulus=34.8e3),
    ultimate_stress_strain_profile=ssp.RectangularStressBlock(
        compressive_strength=50,
        alpha=0.775,
        gamma=0.845,
        ultimate_strain=0.003,
    ),
    flexural_tensile_strength=4.2,
    colour="lightgrey",
)

steel = SteelBar(
    name="500 MPa Steel",
    density=7.85e-6,
    stress_strain_profile=ssp.SteelElasticPlastic(
        yield_strength=500,
        elastic_modulus=200e3,
        fracture_strain=0.05,
    ),
    colour="grey",
)

# construct box by subtracting an inner rectangle from an outer rectangle
outer = rectangular_section(d=1200, b=900, material=concrete)
inner = rectangular_section(d=900, b=600).align_center(align_to=outer)
box = outer - inner

# generate four chamfers
chamfer1 = (
    triangular_section(b=50, h=50, material=concrete)
    .align_to(other=inner, on="left", inner=True)
    .align_to(other=inner, on="bottom", inner=True)
)
chamfer2 = chamfer1.mirror_section(axis="y", mirror_point=(450, 600))
chamfer3 = chamfer1.mirror_section(axis="x", mirror_point=(450, 600))
chamfer4 = chamfer2.mirror_section(axis="x", mirror_point=(450, 600))

# add chamfers to box
geom = box + chamfer1 + chamfer2 + chamfer3 + chamfer4

# add bottom bars
geom = add_bar_rectangular_array(
    geometry=geom,
    area=620,
    material=steel,
    n_x=9,
    x_s=750 / 8,
    anchor=(75, 75),
)

# add top bars
geom = add_bar_rectangular_array(
    geometry=geom,
    area=310,
    material=steel,
    n_x=9,
    x_s=750 / 8,
    anchor=(75, 1125),
)

# add side bars
geom = add_bar_rectangular_array(
    geometry=geom,
    area=200,
    material=steel,
    n_x=2,
    x_s=750,
    n_y=6,
    y_s=150,
    anchor=(75, 225),
)

conc_sec = ConcreteSection(geom)
conc_sec.plot_section()

plt.show()

Here's where my code starts. First I define some functions:

def calculate_MRds(N_Ed: float, Mx_Ed: float, My_Ed: float):
    """
    Calculate resisting bending moments, given internal force components.
    """
    # Convert N from (kN) to (N) and M from (kNm) to (Nmm)
    N_Ed = N_Ed*10**3 #N
    Mx_Ed = Mx_Ed*10**6 #Nmm
    My_Ed = My_Ed*10**6 #Nmm
    
    # Initialise variables
    theta_Ed: float = 0
    Mx_Rd: float = 0
    My_Rd: float = 0
    
    # Initialise necessary variables in case of interpolation  
    # Initialise lists that will contain the alpha and theta values
    Angles = list[float] # Define the type that will be used for lists of angles
    alphas: Angles = [] # Angles in Mx/My plane, measured from Mx axis
    thetas: Angles = [] # Angles of neutral axis, measured from x axis
    
    # Initalise a list that will contain the pairs of resisting moments that make up the domain Mx_Rd, My_Rd
    Moments = list[tuple[float, float]] # Define the type that will be used for lists of bending moments Mx_Rd, My_Rd
    bending_moments: Moments = [] # The list itself
    
    if Mx_Ed != 0 or My_Ed != 0: # Start the calculation only if at least one of Mx_Ed or My_Ed is non-zero 
        if Mx_Ed == 0: # If Mx_Ed = 0, evaluate directly the resisting moment for My_Ed
            Mx_Rd = 0
            if My_Ed > 0:
                theta_Ed = -np.pi/2
                My_Rd = conc_sec.ultimate_bending_capacity(n=N_Ed, theta = theta_Ed).m_xy/1e6
            elif My_Ed < 0:
                theta_Ed = np.pi/2
                My_Rd = -conc_sec.ultimate_bending_capacity(n=N_Ed, theta = theta_Ed).m_xy/1e6
        elif My_Ed == 0: # If My_Ed = 0, evaluate directly the resisting moment for Mx_Ed
            My_Rd = 0
            if Mx_Ed > 0:
                theta_Ed = 0
                Mx_Rd = conc_sec.ultimate_bending_capacity(n=N_Ed, theta = theta_Ed).m_xy/1e6
            elif Mx_Ed < 0:
                theta_Ed = np.pi
                Mx_Rd = -conc_sec.ultimate_bending_capacity(n=N_Ed, theta = theta_Ed).m_xy/1e6
    
        else: # After handling the cases where at least one of the two applied moments is zero, I find theta_Ed by interpolation on alpha_Ed
            if My_Ed > 0: # 0 < alpha < pi (-pi < theta < 0)
                if Mx_Ed > 0: # 0 < alpha < pi/2 (-pi/2 < theta < 0)
                    
                    # Determine the first values of the lists
                    thetas.append(0)
                    bending_moments.append((conc_sec.ultimate_bending_capacity(n=N_Ed, theta = 0).m_xy/1e6, 0))
                    alphas.append(0)
            
                    # Determine the last values of the lists
                    thetas.append(-np.pi/2)
                    bending_moments.append((0, conc_sec.ultimate_bending_capacity(n=N_Ed, theta = -np.pi/2).m_xy/1e6))
                    alphas.append(np.pi/2)
            
                elif Mx_Ed < 0: # pi/2 < alpha < pi (-pi < theta < -pi/2)
                    
                    # Determine the first values of the lists
                    thetas.append(-np.pi/2)
                    bending_moments.append((0, conc_sec.ultimate_bending_capacity(n=N_Ed, theta = -np.pi/2).m_xy/1e6))
                    alphas.append(np.pi/2)
            
                    # Determine the last values of the lists
                    thetas.append(-np.pi)
                    bending_moments.append((conc_sec.ultimate_bending_capacity(n=N_Ed, theta = -np.pi).m_xy/1e6, 0))
                    alphas.append(np.pi)
                    
            elif My_Ed < 0: # -pi < alpha < 0 (0 < theta < pi)
                if Mx_Ed > 0: # -pi/2 < alpha < 0 (0 < theta < pi/2)

                    # Determine the first values of the lists
                    thetas.append(np.pi/2)
                    bending_moments.append((0, conc_sec.ultimate_bending_capacity(n=N_Ed, theta = np.pi/2).m_xy/1e6))
                    alphas.append(-np.pi/2)

                    # Determine the last values of the lists
                    thetas.append(0)
                    bending_moments.append((conc_sec.ultimate_bending_capacity(n=N_Ed, theta = 0).m_xy/1e6, 0))
                    alphas.append(0)
        
                elif Mx_Ed < 0: # -pi < alpha < -pi/2 (pi/2 < theta <pi)

                    # Determine the first values of the lists
                    thetas.append(np.pi)
                    bending_moments.append((conc_sec.ultimate_bending_capacity(n=N_Ed, theta = np.pi).m_xy/1e6, 0))
                    alphas.append(-np.pi)
                    
                    # Determine the last values of the lists
                    thetas.append(np.pi/2)
                    bending_moments.append((0, conc_sec.ultimate_bending_capacity(n=N_Ed, theta = np.pi/2).m_xy/1e6))
                    alphas.append(-np.pi/2)
                    
            theta_Ed, Mx_Rd, My_Rd = _alpha_interp(N_Ed, Mx_Ed, My_Ed, alphas, bending_moments, thetas)
    
    return theta_Ed, Mx_Rd, My_Rd

def _alpha_interp(N_Ed: float, Mx_Ed: float, My_Ed: float, alphas, bending_moments, thetas):
    '''
    Returns theta_Ed, Mx_Rd, My_Rd after determining the theta_Ed corresponding to alpha_Ed through successive interpolations.
    Convergence is achieved when the absolute value of the difference between alpha_Ed and alpha_Rd is below a certain threshold.
    '''
    # Initialize variable alpha corresponding to the load (angle between the positive Mx axis and the vector [(0, 0), (Mx_Ed, My_Ed)])
    alpha_Ed: float = 0
    alpha_Ed = _get_alpha(Mx_Ed, My_Ed)
    
    delta_alpha: float = 1 # Initialize a variable delta_alpha (which will later be calculated as abs(alpha_Ed - alpha_Rd)) with a value greater than the convergence threshold
    while delta_alpha >= 8.7e-3:
        for i, alpha in enumerate(alphas):
            if alpha > alpha_Ed:
                alpha_0 = alphas[i-1]
                alpha_1 = alpha
                theta_0 = thetas[i-1]
                theta_1 = thetas[i]
                break
        theta_Ed = (theta_1 - theta_0)/(alpha_1 - alpha_0) * (alpha_Ed - alpha_0) + theta_0 # Estimate of theta (corresponding to alpha_Ed through interpolation)
        thetas.insert(i, theta_Ed)
        Mx_Rd = conc_sec.ultimate_bending_capacity(n=N_Ed, theta = theta_Ed).m_x/1e6 # The value is made negative if Mx_Ed<0
        My_Rd = conc_sec.ultimate_bending_capacity(n=N_Ed, theta = theta_Ed).m_y/1e6 # The value is made negative if My_Ed<0
        bending_moments.insert(i, (Mx_Rd, My_Rd))
        alpha_Rd = _get_alpha(Mx_Rd, My_Rd)
        alphas.insert(i, alpha_Rd)
        delta_alpha = abs(alpha_Ed-alpha_Rd)

    return theta_Ed, Mx_Rd, My_Rd

def _get_alpha(mx: float, my: float):
    '''
    Returns the value of the angle alpha (angle between the positive Mx axis and the vector [(0,0), (Mx, My)])
    '''
    alpha: float = 0
    if mx > 0: #Mx > 0 (I e IV Quadrant)
        alpha = np.arctan(my/mx)
    elif my > 0: #Mx < 0, My > 0 (II Quadrant)
        alpha = np.arctan(my/mx) + np.pi
    else: #Mx < 0, My < 0 (III Quadrant)
        alpha = np.arctan(my/mx) - np.pi
    return alpha

Then I call the main function to get the resisting b.m.

# Check of biaxial bending

# Internal force components
N_Ed = 5000 #kN
Mx_Ed = 3000 #kNm
My_Ed = -2000 #kNm

# Calculation of resisting bending moments

theta_Ed, Mx_Rd, My_Rd = calculate_MRds(N_Ed, Mx_Ed, My_Ed)

print(f'Mx_Rd = {round(Mx_Rd,1)}kNm')
print(f'My_Rd = {round(My_Rd,1)}kNm')
print(f'θ = {round(theta_Ed/np.pi*180,1)}deg')

Then some plots are added...

# Plot the ultimate stresses (corresponding to the attainment of M_Rd)
M_Rd = conc_sec.ultimate_bending_capacity(theta = theta_Ed, n=N_Ed*10**3) # Calculate the resisting moment with respect to the local reference system (neutral axis)
stress = conc_sec.calculate_ultimate_stress(M_Rd)
stress.plot_stress(units = si_n_mm)
plt.show()

# Plot biaxial bending moment diagram Mx, My
n_points = 24

bb_res = conc_sec.biaxial_bending_diagram(n = N_Ed*10**3, n_points = n_points) #returns "BiaxialBendingResults" type
diag = bb_res.plot_diagram(render = False, eng = True, units = si_kn_m) #render = False prevents the plot from showing before we are finished (otherwise the point corresponding to applied actions won't show).
diag.plot(
    Mx_Ed,
    My_Ed,
    "ok"
)

plt.show()

# Check that the point corresponding to the load is inside the diagram
print(bb_res.point_in_diagram(m_x = Mx_Ed, m_y = My_Ed))