10. AS3600 Design Code#

This example demonstrates how to work with the AS3600:2018 design code. We start by importing the necessary modules.

[1]:
import numpy as np
import matplotlib.pyplot as plt
from concreteproperties.design_codes.as3600 import AS3600
import concreteproperties.stress_strain_profile as ssp
from sectionproperties.pre.library.concrete_sections import concrete_rectangular_section
from concreteproperties.concrete_section import ConcreteSection
from concreteproperties.results import MomentInteractionResults

10.1. Create Design Code and Materials#

In this example 40 MPa concrete will be used with the default 500N steel.

[2]:
design_code = AS3600()
concrete = design_code.create_concrete_material(compressive_strength=40)
steel = design_code.create_steel_material()

We can confirm the concrete material properties.

[3]:
print(concrete.name)
print(f"Density = {concrete.density} kg/mm^3")
concrete.stress_strain_profile.plot_stress_strain(title="Service Profile")
concrete.ultimate_stress_strain_profile.plot_stress_strain(title="Ultimate Profile")
print(
    f"Concrete Flexural Tensile Strength: {concrete.flexural_tensile_strength:.2f} MPa"
)
40 MPa Concrete (AS 3600:2018)
Density = 2.4e-06 kg/mm^3
../_images/notebooks_as3600_5_1.svg
../_images/notebooks_as3600_5_2.svg
Concrete Flexural Tensile Strength: 3.79 MPa

We can confirm the steel material properties.

[4]:
print(steel.name)
print(f"Density = {steel.density} kg/mm^3")
steel.stress_strain_profile.plot_stress_strain()
500 MPa Steel (AS 3600:2018)
Density = 7.85e-06 kg/mm^3
../_images/notebooks_as3600_7_1.svg
[4]:
<Axes: title={'center': 'Stress-Strain Profile'}, xlabel='Strain', ylabel='Stress'>

10.2. Assign Geometry to Design Code#

This example will analyse a 600D x 450W concrete beam with 5N20s top and bottom. After creating the geometry it must be assigned to the design code object.

[5]:
geom = concrete_rectangular_section(
    b=450,
    d=600,
    dia_top=20,
    n_top=5,
    dia_bot=20,
    n_bot=5,
    n_circle=4,
    cover=30,
    area_top=310,
    area_bot=310,
    conc_mat=concrete,
    steel_mat=steel,
)

conc_sec = ConcreteSection(geom)
conc_sec.plot_section()
design_code.assign_concrete_section(concrete_section=conc_sec)
../_images/notebooks_as3600_9_0.svg

10.3. Area Properties#

Obtaining the area properties is identical to that of a ConcreteSection object. The same can be done for a moment-curvature analysis and stress analyses (not carried out in this example).

[6]:
gross_props = design_code.get_gross_properties()
transformed_props = design_code.get_transformed_gross_properties(
    elastic_modulus=concrete.stress_strain_profile.elastic_modulus
)
cracked_props = design_code.calculate_cracked_properties()

gross_props.print_results()
cracked_props.print_results()
     Gross Concrete Section Properties      
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┓
┃ Property                          Value ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━┩
│ Total Area                 2.700000e+05 │
│ Concrete Area              2.669000e+05 │
│ Lumped Reinforcement Area  3.100000e+03 │
│ Axial Rigidity (EA)        9.374320e+09 │
│ Mass (per unit length)     6.648950e-01 │
│ Perimeter                  2.100000e+03 │
│ E.Qx                       2.812296e+12 │
│ E.Qy                       2.109222e+12 │
│ x-Centroid                 2.250000e+02 │
│ y-Centroid                 3.000000e+02 │
│ x-Centroid (Gross)         2.250000e+02 │
│ y-Centroid (Gross)         3.000000e+02 │
│ E.Ixx_g                    1.144420e+15 │
│ E.Iyy_g                    6.329024e+14 │
│ E.Ixy_g                    6.327666e+14 │
│ E.Ixx_c                    3.007311e+14 │
│ E.Iyy_c                    1.583274e+14 │
│ E.Ixy_c                    7.500000e-01 │
│ E.I11                      3.007311e+14 │
│ E.I22                      1.583274e+14 │
│ Principal Axis Angle       0.000000e+00 │
│ E.Zxx+                     1.002437e+12 │
│ E.Zxx-                     1.002437e+12 │
│ E.Zyy+                     7.036774e+11 │
│ E.Zyy-                     7.036774e+11 │
│ E.Z11+                     1.002437e+12 │
│ E.Z11-                     1.002437e+12 │
│ E.Z22+                     7.036774e+11 │
│ E.Z22-                     7.036774e+11 │
│ Ultimate Concrete Strain   3.000000e-03 │
└───────────────────────────┴──────────────┘
  Cracked Concrete Section   
         Properties          
┏━━━━━━━━━━━━┳━━━━━━━━━━━━━━┓
┃ Property           Value ┃
┡━━━━━━━━━━━━╇━━━━━━━━━━━━━━┩
│ theta       0.000000e+00 │
│ n           0.000000e+00 │
│ m           0.000000e+00 │
│ m_cr        1.159750e+08 │
│ d_nc        1.239649e+02 │
│ E.A_cr      2.398873e+09 │
│ E.Qx_cr     1.141948e+12 │
│ E.Qy_cr     5.397464e+11 │
│ x-Centroid  2.250000e+02 │
│ y-Centroid  4.760351e+02 │
│ E.Ixx_g_cr  6.137603e+14 │
│ E.Iyy_g_cr  1.620731e+14 │
│ E.Ixy_g_cr  2.569383e+14 │
│ E.Ixx_c_cr  7.015297e+13 │
│ E.Iyy_c_cr  4.063014e+13 │
│ E.Ixy_c_cr  1.250000e-01 │
│ E.Iuu_cr    7.015297e+13 │
│ E.I11_cr    7.015297e+13 │
│ E.I22_cr    4.063014e+13 │
│ phi_cr      0.000000e+00 │
└────────────┴──────────────┘

10.4. Ultimate Bending Capacity#

The factored ultimate bending capacity can be found by calling the ultimate_bending_capacity() method. This method returns a factored and unfactored UltimateBendingResults object, as well as the capacity reduction factor phi.

[7]:
f_ult_res, ult_res, phi = design_code.ultimate_bending_capacity()
print(f"Muo = {ult_res.m_xy / 1e6:.2f} kN.m")
print(f"kuo = {ult_res.k_u:.4f}")
print(f"phi = {phi:.3f}")
print(f"phi.Muo = {f_ult_res.m_xy / 1e6:.2f} kN.m")
Muo = 414.19 kN.m
kuo = 0.0000
phi = 0.850
phi.Muo = 352.07 kN.m

We can also pass an axial load to ultimate_bending_capacity(). This will calculate the factored moment capacity by ensuring the supplied axial load equals the factored axial capacity, i.e. given a design axial force, what is the maximum moment my section can handle?

[8]:
n_star = 1000e3
f_ult_res, ult_res, phi = design_code.ultimate_bending_capacity(n_design=n_star)
print(f"N* = {n_star / 1e3:.0f} kN")
print(f"Nu = {ult_res.n / 1e3:.2f} kN")
print(f"Mu = {ult_res.m_xy / 1e6:.2f} kN.m")
print(f"phi = {phi:.3f}")
print(f"phi.Nu = {f_ult_res.n / 1e3:.0f} kN")
print(f"phi.Mu = {f_ult_res.m_xy / 1e6:.2f} kN.m")
N* = 1000 kN
Nu = 1312.32 kN
Mu = 724.37 kN.m
phi = 0.762
phi.Nu = 1000 kN
phi.Mu = 551.98 kN.m

10.5. Moment Interaction Diagram#

The factored moment interaction diagram can be found by calling the moment_interaction_diagram() method. This method returns a factored and unfactored MomentInteractionResults object.

[9]:
f_mi_res, mi_res, phis = design_code.moment_interaction_diagram(progress_bar=False)

We can compare the factored and unfactored moment interaction diagrams.

[10]:
MomentInteractionResults.plot_multiple_diagrams(
    [f_mi_res, mi_res], ["Factored", "Unfactored"], fmt="-"
)
../_images/notebooks_as3600_19_0.svg
[10]:
<Axes: title={'center': 'Moment Interaction Diagram'}, xlabel='Bending Moment', ylabel='Axial Force'>

Using the list of capacity reduction factors phis, we can visualise how the capacity reduction factor varies as a function of the applied axial load.

[11]:
fig, ax = plt.subplots()
n_list, _ = mi_res.get_results_lists(moment="m_x")  # get list of axial loads
ax.plot(np.array(n_list) / 1e3, phis, "-x")
plt.xlabel("Axial Force [kN]")
plt.ylabel("$\phi$")
plt.grid()
plt.show()
../_images/notebooks_as3600_21_0.svg

We can check to see if a combination of axial force and bending moment lies within the moment interaction diagram using the point_in_diagram() method.

[12]:
# design load cases
n_stars = [4000e3, 5000e3, -500e3, 1000e3]
m_stars = [200e6, 400e6, 100e6, 551e6]
marker_styles = ["x", "+", "o", "*"]
n_cases = len(n_stars)

# plot moment interaction diagram
ax = f_mi_res.plot_diagram(fmt="k-", render=False)

# check to see if combination is within diagram and plot result
for idx in range(n_cases):
    case = f_mi_res.point_in_diagram(n=n_stars[idx], m=m_stars[idx])
    print("Case {num}: {status}".format(num=idx + 1, status="OK" if case else "FAIL"))
    ax.plot(
        m_stars[idx] / 1e6,
        n_stars[idx] / 1e3,
        "k" + marker_styles[idx],
        markersize=10,
        label=f"Case {idx + 1}",
    )

ax.legend()
plt.show()
Case 1: OK
Case 2: FAIL
Case 3: OK
Case 4: OK
../_images/notebooks_as3600_23_1.svg

Let’s compare moment interaction diagrams using a rectangular stress block, a bilinear stress profile and a parabolic stress profile. AS 3600:2018 restricts the maximum stress in the profile to 0.9 * f'c.

[13]:
# bilinear
concrete.ultimate_stress_strain_profile = ssp.BilinearStressStrain(
    compressive_strength=0.9 * 40,
    compressive_strain=0.0015,
    ultimate_strain=0.003,
)
f_mi_res_bil, _, _ = design_code.moment_interaction_diagram(progress_bar=False)

# parabolic
concrete.ultimate_stress_strain_profile = ssp.EurocodeParabolicUltimate(
    compressive_strength=0.9 * 40,
    compressive_strain=0.0015,
    ultimate_strain=0.003,
    n=2,
)
f_mi_res_par, _, _ = design_code.moment_interaction_diagram(progress_bar=False)
[14]:
MomentInteractionResults.plot_multiple_diagrams(
    [f_mi_res, f_mi_res_bil, f_mi_res_par],
    ["Rectangular", "Bilinear", "Parabolic"],
    fmt="-",
)
../_images/notebooks_as3600_26_0.svg
[14]:
<Axes: title={'center': 'Moment Interaction Diagram'}, xlabel='Bending Moment', ylabel='Axial Force'>

10.6. Biaxial Bending Diagram#

We can also compute factored biaxial bending diagrams by calling the biaxial_bending_diagram() method. This method returns a factored BiaxialBendingResults object as well as a list of the capacity reduction factors phis.

[15]:
# reset concrete ultimate profile
concrete.ultimate_stress_strain_profile = ssp.RectangularStressBlock(
    compressive_strength=40,
    alpha=0.79,
    gamma=0.87,
    ultimate_strain=0.003,
)

# create biaxial bending diagram
f_bb_res1, phis1 = design_code.biaxial_bending_diagram(progress_bar=False)
bb_res1 = conc_sec.biaxial_bending_diagram(progress_bar=False)
f_bb_res2, phis2 = design_code.biaxial_bending_diagram(
    n_design=2000e3, progress_bar=False
)
bb_res2 = conc_sec.biaxial_bending_diagram(n=2000e3, progress_bar=False)
[16]:
# plot case 1
ax = f_bb_res1.plot_diagram(fmt="-", render=False)
bb_res1.plot_diagram(fmt="-", ax=ax)
plt.show()
print(f"Average phi = {np.mean(phis1):.3f}")

# plot case 2
ax = f_bb_res2.plot_diagram(fmt="-", render=False)
bb_res2.plot_diagram(fmt="-", ax=ax)
plt.show()
print(f"Average phi = {np.mean(phis2):.3f}")
../_images/notebooks_as3600_29_0.svg
Average phi = 0.850
../_images/notebooks_as3600_29_2.svg
Average phi = 0.618

Note that as the bending angle changes, k_uo and N_ub change, resulting in a constantly changing value for phi.