7. Stress Analysis

This example demonstrates how to undertake stress analyses and view their results. We start by importing the necessary modules.

import numpy as np
from rich.pretty import pprint
from concreteproperties.material import Concrete, Steel
from concreteproperties.stress_strain_profile import (
    EurocodeNonLinear,
    RectangularStressBlock,
    SteelElasticPlastic,
)
from sectionproperties.pre.library.concrete_sections import concrete_circular_section
from concreteproperties.concrete_section import ConcreteSection

7.1. Assign Materials

The materials used in this example will be 40 MPa concrete and 500 MPa steel, specified in accordance with AS 3600:2018. The service stress-strain profile for the concrete will be adapted from EC2.

concrete = Concrete(
    name="40 MPa Concrete",
    density=2.4e-6,
    stress_strain_profile=EurocodeNonLinear(
        elastic_modulus=32.8e3,
        ultimate_strain=0.0035,
        compressive_strength=40,
        compressive_strain=0.0023,
        tensile_strength=3.8,
        tension_softening_stiffness=10e3,
    ),
    ultimate_stress_strain_profile=RectangularStressBlock(
        compressive_strength=40,
        alpha=0.79,
        gamma=0.87,
        ultimate_strain=0.003,
    ),
    alpha_squash=0.85,
    flexural_tensile_strength=3.8,
    colour="lightgrey",
)

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

7.2. Create Geometry and Concrete Section

The section in this example is a 600 mm diameter concrete column, reinforced with 10N20 bars with 45 mm cover. The geometry is generated using the sectionproperties concrete sections library.

geom = concrete_circular_section(
    d=600,
    n=32,
    dia=20,
    n_bar=10,
    n_circle=4,
    area_conc=np.pi * 600 * 600 / 4,
    area_bar=310,
    cover=45,
    conc_mat=concrete,
    steel_mat=steel,
)

conc_sec = ConcreteSection(geom)
conc_sec.plot_section()

<AxesSubplot:title={'center':'Reinforced Concrete Section'}>

7.3. Elastic Uncracked Stress Analysis

We perform an uncracked stress analysis using the calculate_uncracked_stress() method. This analysis assumes the entire concrete section acts to resist the applied loading. Note that concreteproperties does not check to see if the provided moment is above the cracking moment, this should be confirmed by performing a cracking analysis.

uncr_stress_res_1 = conc_sec.calculate_uncracked_stress(m_x=50e6)
uncr_stress_res_2 = conc_sec.calculate_uncracked_stress(m_x=25e6, m_y=35e6, n=200e3)

We can visualise the stresses by calling the plot_stress() method.

uncr_stress_res_1.plot_stress()

<AxesSubplot:title={'center':'Stress'}>

uncr_stress_res_2.plot_stress()

<AxesSubplot:title={'center':'Stress'}>

7.4. Elastic Cracked Stress Analysis

A cracked stress analysis first requires a cracking analysis to be performed by calling the calculate_cracked_properties() method. Here we perform a cracking analysis for bending about the x axis.

cracked_res = conc_sec.calculate_cracked_properties(theta=0)
print(f"M_cr = {cracked_res.m_cr / 1e6:.2f} kN.m")
print(f"d_n,c = {cracked_res.d_nc:.2f} mm")

M_cr = 86.32 kN.m
d_n,c = 143.81 mm

Since the cracking moment is 86 kN.m, let’s see what the stresses within the cross-section look like under 150 kN.m bending moment. Note that the moment will be applied about the same axis as used in the cracking analysis.

m_ext = 150e6
cracked_stress_res = conc_sec.calculate_cracked_stress(
    cracked_results=cracked_res, m=m_ext
)
cracked_stress_res.plot_stress()

<AxesSubplot:title={'center':'Stress'}>

Many useful results are stored within the StressResult object. Let’s examine a few of these results.

7.4.1. Concrete Results

Concrete geometries that remain after cracking are stored in the concrete_analysis_sections attribute as a list of AnalysisSection objects. We can visualise this geometry by calling the plot_mesh() method. The StressResult object also stores the stresses at the nodes in each AnalysisSection in the concrete_stresses attribute - these are the stresses you see on the plot_stress() plots. Finally, the concrete_forces attribute stores the resultant force for each AnalysisSection and its lever arm about the neutral axis.

for idx, an_sec in enumerate(cracked_stress_res.concrete_analysis_sections):
    # Label section and plot section mesh
    print(f"Analysis Section {idx + 1}")
    an_sec.plot_mesh()

    # get concrete results
    sigs_conc = cracked_stress_res.concrete_stresses[idx]
    f_conc = cracked_stress_res.concrete_forces[idx][0]
    d_conc = cracked_stress_res.concrete_forces[idx][1]
    m_conc = f_conc * d_conc

    # print results
    print("Concrete Stresses:")
    pprint(sigs_conc)
    print("Concrete Net Force & Lever Arm:")
    print(f"F_c = {f_conc / 1e3:.2f} kN")
    print(f"d_n = {d_conc:.2f} mm")
    print(f"F_c = {m_conc / 1e6:.2f} kN.m")

Analysis Section 1

Concrete Stresses:

array([5.86929734e-05, 1.17773019e+00, 6.52033244e+00, 1.09048928e+01,
│      1.41629150e+01, 1.61691952e+01, 1.68466332e+01, 1.61691952e+01,
│      1.41629150e+01, 1.09048928e+01, 6.52033244e+00, 1.17773019e+00,
│      5.86929734e-05, 8.88584105e+00, 1.03442622e+01, 8.88584105e+00,
│      7.42741987e+00, 7.42741987e+00, 8.88584105e+00, 1.03442622e+01,
│      8.88584105e+00])

Concrete Net Force & Lever Arm:
F_c = 347.99 kN
d_n = 82.77 mm
F_c = 28.80 kN.m

7.4.2. Steel Results

Steel geometries are stored in the steel_geometries attribute of the StressResult object. Similar to the concrete results, the stresses and forces/lever arms are stored, as well as the strains. Here we print the stress for each bar to the terminal.

from rich.table import Table
from rich.console import Console

# store forces & moments for later
forces = []
moments = []

# create a Rich table for pretty printing
table = Table(title="Steel Results")
table.add_column("Bar No.", justify="center", style="cyan", no_wrap=True)
table.add_column("Location (x, y) (mm)", justify="center", style="green")
table.add_column("Stress (MPa)", justify="center", style="green")
table.add_column("Force (kN)", justify="center", style="green")
table.add_column("Lever Arm (mm)", justify="center", style="green")
table.add_column("Moment (kN.m)", justify="center", style="green")

for idx, steel_geom in enumerate(cracked_stress_res.steel_geometries):
    # get the steel results
    centroid = steel_geom.calculate_centroid()
    stress = cracked_stress_res.steel_stresses[idx]
    strain = cracked_stress_res.steel_strains[idx]
    force, d = cracked_stress_res.steel_forces[idx]

    # calculate the moment each bar creates and store the results
    moment = force * d
    forces.append(force)
    moments.append(moment)

    # print compression or tension
    if strain > 0:
        t_or_c = "C"
    else:
        t_or_c = "T"

    table.add_row(
        f"{idx + 1}",
        f"({centroid[0]:.1f}, {centroid[1]:.1f})",
        f"{abs(stress):.1f} ({t_or_c})",
        f"{abs(force) / 1e3:.1f}",
        f"{d:.1f}",
        f"{moment / 1e6:.2f}",
    )

console = Console()
console.print(table)

                                         Steel Results                                         
┏━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┓
┃ Bar No. ┃ Location (x, y) (mm) ┃ Stress (MPa) ┃ Force (kN) ┃ Lever Arm (mm) ┃ Moment (kN.m) ┃
┡━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━┩
│    1    │    (245.0, -0.0)     │  112.3 (T)   │    34.8    │     -157.2     │     5.47      │
│    2    │    (198.2, 144.0)    │   9.4 (T)    │    2.9     │     -13.1      │     0.04      │
│    3    │    (75.7, 233.0)     │   54.2 (C)   │    16.8    │      75.9      │     1.27      │
│    4    │    (-75.7, 233.0)    │   54.2 (C)   │    16.8    │      75.9      │     1.27      │
│    5    │   (-198.2, 144.0)    │   9.4 (T)    │    2.9     │     -13.1      │     0.04      │
│    6    │    (-245.0, -0.0)    │  112.3 (T)   │    34.8    │     -157.2     │     5.47      │
│    7    │   (-198.2, -144.0)   │  215.1 (T)   │    66.7    │     -301.2     │     20.08     │
│    8    │   (-75.7, -233.0)    │  278.7 (T)   │    86.4    │     -390.2     │     33.71     │
│    9    │    (75.7, -233.0)    │  278.7 (T)   │    86.4    │     -390.2     │     33.71     │
│   10    │   (198.2, -144.0)    │  215.1 (T)   │    66.7    │     -301.2     │     20.08     │
└─────────┴──────────────────────┴──────────────┴────────────┴────────────────┴───────────────┘

Finally, let’s confirm that the net axial force is indeed zero and the internal moment is equal to the external moment.

# sum of forces
int_force = sum(forces) + f_conc
print(f"Sum of Internal Forces: {int_force / 1e3:.0f} kN")

# sum of moments
int_moment = sum(moments) + m_conc
print(f"Sum of Internal Moments: {int_moment / 1e6:.0f} kN.m")
print(f"External Moment: {m_ext / 1e6:.0f} kN.m")

Sum of Internal Forces: 0 kN
Sum of Internal Moments: 150 kN.m
External Moment: 150 kN.m

7.5. Service Stress Analysis

A service stress analysis displays internal stresses based off the stress-strain profiles provided to the material properties. Before viewing the stresses, a moment curvature analysis must be performed to establish the structural response of the reinforced concrete section.

mk_res = conc_sec.moment_curvature_analysis()

Let’s plot the moment curvature diagram so that we can choose some interesting points to generate stress results.

mk_res.plot_results()

<AxesSubplot:title={'center':'Moment-Curvature'}, xlabel='Curvature', ylabel='Moment'>

Based on the above, we compile a list of moments and curvatures at which we want to visualise the cross-section stresses. We use curvatures as input around the softening region as there are two possible values of moment.

# create a list of inputs (either moments [m] or curvatures [k])
stress_inputs = [
    {"m_or_k": "m", "val": 50e6},
    {"m_or_k": "m", "val": cracked_res.m_cr},
    {"m_or_k": "k", "val": 7.5e-07},
    {"m_or_k": "k", "val": 1.15e-06},
    {"m_or_k": "k", "val": 2.15e-06},
    {"m_or_k": "m", "val": 250e6},
    {"m_or_k": "m", "val": 300e6},
    {"m_or_k": "m", "val": 360e6},
]

# loop through each input
for s_in in stress_inputs:
    # determine if the input is a moment or curvature
    if s_in["m_or_k"] == "m":
        m_stress = s_in["val"]
        service_stress_res = conc_sec.calculate_service_stress(
            moment_curvature_results=mk_res, m=m_stress
        )
    elif s_in["m_or_k"] == "k":
        service_stress_res = conc_sec.calculate_service_stress(
            moment_curvature_results=mk_res, m=None, kappa=s_in["val"]
        )
        m_stress = service_stress_res.sum_moments()

    # create plot title and plot stress
    label = f"Moment = {m_stress / 1e6:.0f} kN.m"
    service_stress_res.plot_stress(title=label)

Viewing the above stress plots in combination with the moment curvature diagram is an insightful way of understanding the bending response of a reinforced concrete section!

7.6. Ultimate Stress Analysis

An ultimate stress analysis displayed the internal stresses at the ultimate bending moment. In this example we examine the ultimate stresses at the pure bending point, the balanced point and the decompression point. To locate these points we first plot the moment interaction diagram.

mi_res = conc_sec.moment_interaction_diagram()
mi_res.plot_diagram()

<AxesSubplot:title={'center':'Moment Interaction Diagram'}, xlabel='Bending Moment', ylabel='Axial Force'>

Prior to plotting the stresses, the ultimate bending capacities must be calculated.

ultimate_res_pure = conc_sec.ultimate_bending_capacity()
ultimate_res_bal = conc_sec.ultimate_bending_capacity(n=3500e3)
ultimate_res_decomp = conc_sec.ultimate_bending_capacity(n=9000e3)

We calculate the ultimate stresses by calling the calculate_ultimate_stress() method.

ultimate_stress_pure = conc_sec.calculate_ultimate_stress(
    ultimate_results=ultimate_res_pure
)
ultimate_stress_bal = conc_sec.calculate_ultimate_stress(
    ultimate_results=ultimate_res_bal
)
ultimate_stress_decomp = conc_sec.calculate_ultimate_stress(
    ultimate_results=ultimate_res_decomp
)

ultimate_stress_pure.plot_stress(title="Pure Bending")
ultimate_stress_bal.plot_stress(title="Balanced")
ultimate_stress_decomp.plot_stress(title="Decompression")

<AxesSubplot:title={'center':'Decompression'}>

As a rectangular stress block is used for the ultimate stress-strain profile, the entire cross-section is not in compression due to the zero stress region above the neutral axis in the stress-strain profile.