from __future__ import annotations
import warnings
from math import inf, isinf, nan
from typing import TYPE_CHECKING, List, Optional, Tuple, Union
import matplotlib.patches as mpatches
import numpy as np
import sectionproperties.pre.geometry as sp_geom
from rich.live import Live
from scipy.interpolate import interp1d
from scipy.optimize import brentq
import concreteproperties.results as res
import concreteproperties.utils as utils
from concreteproperties.analysis_section import AnalysisSection
from concreteproperties.material import Concrete
from concreteproperties.post import plotting_context
from concreteproperties.pre import CPGeom, CPGeomConcrete
if TYPE_CHECKING:
import matplotlib
[docs]class ConcreteSection:
"""Class for a reinforced concrete section."""
def __init__(
self,
geometry: sp_geom.CompoundGeometry,
):
"""Inits the ConcreteSection class.
:param geometry: *sectionproperties* CompoundGeometry object describing the
reinforced concrete section
"""
self.compound_geometry = geometry
# check overlapping regions
polygons = [sec_geom.geom for sec_geom in self.compound_geometry.geoms]
overlapped_regions = sp_geom.check_geometry_overlaps(polygons)
if overlapped_regions:
warnings.warn(
"The provided geometry contains overlapping regions, results may be incorrect."
)
# sort into concrete and reinforcement (meshed and lumped) geometries
self.all_geometries: List[Union[CPGeomConcrete, CPGeom]] = []
self.meshed_geometries: List[Union[CPGeomConcrete, CPGeom]] = []
self.concrete_geometries: List[CPGeomConcrete] = []
self.reinf_geometries_meshed: List[CPGeom] = []
self.reinf_geometries_lumped: List[CPGeom] = []
# sort geometry into appropriate list
for geom in self.compound_geometry.geoms:
if isinstance(geom.material, Concrete):
cp_geom = CPGeomConcrete(geom=geom.geom, material=geom.material)
self.concrete_geometries.append(cp_geom)
self.meshed_geometries.append(cp_geom)
else:
cp_geom = CPGeom(geom=geom.geom, material=geom.material) # type: ignore
if cp_geom.material.meshed:
self.reinf_geometries_meshed.append(cp_geom)
self.meshed_geometries.append(cp_geom)
elif not cp_geom.material.meshed:
self.reinf_geometries_lumped.append(cp_geom)
self.all_geometries.append(cp_geom)
# initialise gross properties results class
self.gross_properties = res.GrossProperties()
# calculate gross area properties
self.calculate_gross_area_properties()
[docs] def calculate_gross_area_properties(
self,
):
"""Calculates and stores gross section area properties."""
# loop through all geometries
for geom in self.all_geometries:
# area and centroid of geometry
area = geom.calculate_area()
centroid = geom.calculate_centroid()
self.gross_properties.total_area += area
self.gross_properties.e_a += area * geom.material.elastic_modulus
self.gross_properties.mass += area * geom.material.density
self.gross_properties.e_qx += (
area * geom.material.elastic_modulus * centroid[1]
)
self.gross_properties.e_qy += (
area * geom.material.elastic_modulus * centroid[0]
)
# sum concrete areas
for conc_geom in self.concrete_geometries:
self.gross_properties.concrete_area += conc_geom.calculate_area()
# sum reinforcement meshed areas
for meshed_geom in self.reinf_geometries_meshed:
self.gross_properties.reinf_meshed_area += meshed_geom.calculate_area()
# sum reinforcement lumped areas
for lumped_geom in self.reinf_geometries_lumped:
self.gross_properties.reinf_lumped_area += lumped_geom.calculate_area()
# perimeter
self.gross_properties.perimeter = self.compound_geometry.calculate_perimeter()
# centroids
self.gross_properties.cx = (
self.gross_properties.e_qy / self.gross_properties.e_a
)
self.gross_properties.cy = (
self.gross_properties.e_qx / self.gross_properties.e_a
)
# global second moments of area
# meshed geometries
for geom in self.meshed_geometries:
sec = AnalysisSection(geometry=geom)
for el in sec.elements:
el_e_ixx_g, el_e_iyy_g, el_e_ixy_g = el.second_moments_of_area()
self.gross_properties.e_ixx_g += el_e_ixx_g
self.gross_properties.e_iyy_g += el_e_iyy_g
self.gross_properties.e_ixy_g += el_e_ixy_g
# lumped geometries - treat as lumped circles
for geom in self.reinf_geometries_lumped:
# area, diameter and centroid of geometry
area = geom.calculate_area()
diam = np.sqrt(4 * area / np.pi)
centroid = geom.calculate_centroid()
self.gross_properties.e_ixx_g += geom.material.elastic_modulus * (
np.pi * pow(diam, 4) / 64 + area * centroid[1] * centroid[1]
)
self.gross_properties.e_iyy_g += geom.material.elastic_modulus * (
np.pi * pow(diam, 4) / 64 + area * centroid[0] * centroid[0]
)
self.gross_properties.e_ixy_g += geom.material.elastic_modulus * (
area * centroid[0] * centroid[1]
)
# centroidal second moments of area
self.gross_properties.e_ixx_c = (
self.gross_properties.e_ixx_g
- self.gross_properties.e_qx**2 / self.gross_properties.e_a
)
self.gross_properties.e_iyy_c = (
self.gross_properties.e_iyy_g
- self.gross_properties.e_qy**2 / self.gross_properties.e_a
)
self.gross_properties.e_ixy_c = (
self.gross_properties.e_ixy_g
- self.gross_properties.e_qx
* self.gross_properties.e_qy
/ self.gross_properties.e_a
)
# principal 2nd moments of area about the centroidal xy axis
Delta = (
((self.gross_properties.e_ixx_c - self.gross_properties.e_iyy_c) / 2) ** 2
+ self.gross_properties.e_ixy_c**2
) ** 0.5
self.gross_properties.e_i11 = (
self.gross_properties.e_ixx_c + self.gross_properties.e_iyy_c
) / 2 + Delta
self.gross_properties.e_i22 = (
self.gross_properties.e_ixx_c + self.gross_properties.e_iyy_c
) / 2 - Delta
# principal axis angle
if (
abs(self.gross_properties.e_ixx_c - self.gross_properties.e_i11)
< 1e-12 * self.gross_properties.e_i11
):
self.gross_properties.phi = 0
else:
self.gross_properties.phi = np.arctan2(
self.gross_properties.e_ixx_c - self.gross_properties.e_i11,
self.gross_properties.e_ixy_c,
)
# centroidal section moduli
x_min, x_max, y_min, y_max = self.compound_geometry.calculate_extents()
self.gross_properties.e_zxx_plus = self.gross_properties.e_ixx_c / abs(
y_max - self.gross_properties.cy
)
self.gross_properties.e_zxx_minus = self.gross_properties.e_ixx_c / abs(
y_min - self.gross_properties.cy
)
self.gross_properties.e_zyy_plus = self.gross_properties.e_iyy_c / abs(
x_max - self.gross_properties.cx
)
self.gross_properties.e_zyy_minus = self.gross_properties.e_iyy_c / abs(
x_min - self.gross_properties.cx
)
# principal section moduli
x11_max, x11_min, y22_max, y22_min = utils.calculate_local_extents(
geometry=self.compound_geometry,
cx=self.gross_properties.cx,
cy=self.gross_properties.cy,
theta=self.gross_properties.phi,
)
# evaluate principal section moduli
self.gross_properties.e_z11_plus = self.gross_properties.e_i11 / abs(y22_max)
self.gross_properties.e_z11_minus = self.gross_properties.e_i11 / abs(y22_min)
self.gross_properties.e_z22_plus = self.gross_properties.e_i22 / abs(x11_max)
self.gross_properties.e_z22_minus = self.gross_properties.e_i22 / abs(x11_min)
# store ultimate concrete strain (get smallest from all concrete geometries)
conc_ult_strain = 0
for idx, conc_geom in enumerate(self.concrete_geometries):
ult_strain = (
conc_geom.material.ultimate_stress_strain_profile.get_ultimate_compressive_strain()
)
if idx == 0:
conc_ult_strain = ult_strain
else:
conc_ult_strain = min(conc_ult_strain, ult_strain)
self.gross_properties.conc_ultimate_strain = conc_ult_strain
[docs] def get_gross_properties(
self,
) -> res.GrossProperties:
"""Returns the gross section properties of the reinforced concrete section.
:return: Gross concrete properties object
"""
return self.gross_properties
[docs] def calculate_cracked_properties(
self,
theta: float = 0,
) -> res.CrackedResults:
r"""Calculates cracked section properties given a neutral axis angle ``theta``.
:param theta: Angle (in radians) the neutral axis makes with the horizontal axis
(:math:`-\pi \leq \theta \leq \pi`)
:return: Cracked results object
"""
cracked_results = res.CrackedResults(theta=theta)
cracked_results.m_cr = self.calculate_cracking_moment(theta=theta)
# set neutral axis depth limits
# depth of neutral axis at extreme tensile fibre
_, d_t = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=theta
)
a = 1e-6 * d_t # sufficiently small depth of compressive zone
b = d_t # neutral axis at extreme tensile fibre
# find neutral axis that gives convergence of the the cracked neutral axis
try:
(cracked_results.d_nc, r) = brentq(
f=self.cracked_neutral_axis_convergence,
a=a,
b=b,
args=(cracked_results),
xtol=1e-3,
rtol=1e-6, # type: ignore
full_output=True,
disp=False,
)
except ValueError:
warnings.warn("brentq algorithm failed.")
# calculate cracked section properties
# axial rigidity & first moments of area
for geom in cracked_results.cracked_geometries:
area = geom.calculate_area()
centroid = geom.calculate_centroid()
cracked_results.e_a_cr += area * geom.material.elastic_modulus
cracked_results.e_qx_cr += (
area * geom.material.elastic_modulus * centroid[1]
)
cracked_results.e_qy_cr += (
area * geom.material.elastic_modulus * centroid[0]
)
# centroids
cracked_results.cx = cracked_results.e_qy_cr / cracked_results.e_a_cr
cracked_results.cy = cracked_results.e_qx_cr / cracked_results.e_a_cr
# global second moments of area
for geom in cracked_results.cracked_geometries:
# if meshed
if geom.material.meshed:
sec = AnalysisSection(geometry=geom)
for el in sec.elements:
el_e_ixx_g, el_e_iyy_g, el_e_ixy_g = el.second_moments_of_area()
cracked_results.e_ixx_g_cr += el_e_ixx_g
cracked_results.e_iyy_g_cr += el_e_iyy_g
cracked_results.e_ixy_g_cr += el_e_ixy_g
# if lumped
else:
# area, diameter and centroid of geometry
area = geom.calculate_area()
diam = np.sqrt(4 * area / np.pi)
centroid = geom.calculate_centroid()
cracked_results.e_ixx_g_cr += geom.material.elastic_modulus * (
np.pi * pow(diam, 4) / 64 + area * centroid[1] * centroid[1]
)
cracked_results.e_iyy_g_cr += geom.material.elastic_modulus * (
np.pi * pow(diam, 4) / 64 + area * centroid[0] * centroid[0]
)
cracked_results.e_ixy_g_cr += geom.material.elastic_modulus * (
area * centroid[0] * centroid[1]
)
# centroidal second moments of area
cracked_results.e_ixx_c_cr = (
cracked_results.e_ixx_g_cr
- cracked_results.e_qx_cr**2 / cracked_results.e_a_cr
)
cracked_results.e_iyy_c_cr = (
cracked_results.e_iyy_g_cr
- cracked_results.e_qy_cr**2 / cracked_results.e_a_cr
)
cracked_results.e_ixy_c_cr = (
cracked_results.e_ixy_g_cr
- cracked_results.e_qx_cr * cracked_results.e_qy_cr / cracked_results.e_a_cr
)
cracked_results.e_iuu_cr = (
cracked_results.e_iyy_c_cr * (np.sin(theta)) ** 2
+ cracked_results.e_ixx_c_cr * (np.cos(theta)) ** 2
- 2 * cracked_results.e_ixy_c_cr * np.sin(theta) * np.cos(theta)
)
# principal 2nd moments of area about the centroidal xy axis
Delta = (
((cracked_results.e_ixx_c_cr - cracked_results.e_iyy_c_cr) / 2) ** 2
+ cracked_results.e_ixy_c_cr**2
) ** 0.5
cracked_results.e_i11_cr = (
cracked_results.e_ixx_c_cr + cracked_results.e_iyy_c_cr
) / 2 + Delta
cracked_results.e_i22_cr = (
cracked_results.e_ixx_c_cr + cracked_results.e_iyy_c_cr
) / 2 - Delta
# principal axis angle
if (
abs(cracked_results.e_ixx_c_cr - cracked_results.e_i11_cr)
< 1e-12 * cracked_results.e_i11_cr
):
cracked_results.phi_cr = 0
else:
cracked_results.phi_cr = np.arctan2(
cracked_results.e_ixx_c_cr - cracked_results.e_i11_cr,
cracked_results.e_ixy_c_cr,
)
return cracked_results
[docs] def calculate_cracking_moment(
self,
theta: float,
) -> float:
r"""Calculates the cracking moment given a bending angle ``theta``.
:param theta: Angle (in radians) the neutral axis makes with the
horizontal axis (:math:`-\pi \leq \theta \leq \pi`)
:return: Cracking moment
"""
# get centroidal second moments of area
e_ixx = self.gross_properties.e_ixx_c
e_iyy = self.gross_properties.e_iyy_c
e_ixy = self.gross_properties.e_ixy_c
# determine rotated second moment of area
e_iuu = (
e_iyy * (np.sin(theta)) ** 2
+ e_ixx * (np.cos(theta)) ** 2
- 2 * e_ixy * np.sin(theta) * np.cos(theta)
)
# loop through all concrete geometries to find lowest cracking moment
m_c = 0
for idx, conc_geom in enumerate(self.concrete_geometries):
# get distance from centroid to extreme tensile fibre
d = utils.calculate_max_bending_depth(
points=conc_geom.points,
c_local_v=utils.global_to_local(
theta=theta, x=self.gross_properties.cx, y=self.gross_properties.cy
)[1],
theta=theta,
)
# if no part of the section is in tension, go to next geometry
if d == 0:
continue
# cracking moment for this geometry
f_t = conc_geom.material.flexural_tensile_strength
m_c_geom = (f_t / conc_geom.material.elastic_modulus) * (e_iuu / d)
# if we are the first geometry, initialise cracking moment
if idx == 0:
m_c = m_c_geom
# otherwise take smaller cracking moment
else:
m_c = min(m_c, m_c_geom)
return m_c
[docs] def cracked_neutral_axis_convergence(
self,
d_nc: float,
cracked_results: res.CrackedResults,
) -> float:
"""Given a trial cracked neutral axis depth ``d_nc``, determines the difference
between the first moments of area above and below the trial axis.
:param d_nc: Trial cracked neutral axis
:param cracked_results: Cracked results object
:return: Cracked neutral axis convergence
"""
# calculate extreme fibre in global coordinates
extreme_fibre, d_t = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=cracked_results.theta
)
# validate d_nc input
if d_nc <= 0:
raise ValueError("d_nc must be positive.")
elif d_nc > d_t:
raise ValueError("d_nc must lie within the section, i.e. d_nc <= d_t")
# find point on neutral axis by shifting by d_nc
point_na = utils.point_on_neutral_axis(
extreme_fibre=extreme_fibre, d_n=d_nc, theta=cracked_results.theta
)
# get principal coordinates of neutral axis
na_local = utils.global_to_local(
theta=cracked_results.theta, x=point_na[0], y=point_na[1]
)
# split concrete geometries above and below d_nc, discard below
cracked_geoms: List[Union[CPGeomConcrete, CPGeom]] = []
for conc_geom in self.concrete_geometries:
top_geoms, _ = conc_geom.split_section(
point=point_na,
theta=cracked_results.theta,
)
# save compression geometries
cracked_geoms.extend(top_geoms)
# add reinforcement geometries to list
cracked_geoms.extend(self.reinf_geometries_meshed)
cracked_geoms.extend(self.reinf_geometries_lumped)
# determine moment of area equilibrium about neutral axis
e_qu = 0 # initialise first moment of area
for geom in cracked_geoms:
ea = geom.calculate_area() * geom.material.elastic_modulus
centroid = geom.calculate_centroid()
# convert centroid to local coordinates
_, c_v = utils.global_to_local(
theta=cracked_results.theta, x=centroid[0], y=centroid[1]
)
# calculate first moment of area
e_qu += ea * (c_v - na_local[1])
cracked_results.cracked_geometries = cracked_geoms
return e_qu
[docs] def moment_curvature_analysis(
self,
theta: float = 0,
kappa_inc: float = 1e-7,
delta_m_min: float = 0.15,
delta_m_max: float = 0.3,
) -> res.MomentCurvatureResults:
r"""Performs a moment curvature analysis given a bending angle ``theta``.
Analysis continues until a material reaches its ultimate strain.
:param: Angle (in radians) the neutral axis makes with the horizontal axis
(:math:`-\pi \leq \theta \leq \pi`)
:param kappa_inc: Initial curvature increment
:param delta_m_min: Relative change in moment at which to double step
:param delta_m_max: Relative change in moment at which to halve step
:return: Moment curvature results object
"""
# initialise variables
moment_curvature = res.MomentCurvatureResults(theta=theta)
iter = 0
# set neutral axis depth limits
# depth of neutral axis at extreme tensile fibre
_, d_t = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=theta
)
a = 1e-6 * d_t # sufficiently small depth of compressive zone
b = d_t # neutral axis at extreme tensile fibre
# create progress bar
progress = utils.create_unknown_progress()
with Live(progress, refresh_per_second=10) as live:
task = progress.add_task(
description="[red]Generating M-K diagram",
total=None,
)
# while there hasn't been a failure
while not moment_curvature._failure:
# calculate adaptive step size for curvature
if iter > 1:
moment_diff = (
abs(moment_curvature.kappa[-1] - moment_curvature.kappa[-2])
/ moment_curvature.kappa[-1]
)
if moment_diff <= delta_m_min:
kappa_inc *= 2
elif moment_diff >= delta_m_max:
kappa_inc *= 0.5
kappa = 0 if iter == 0 else moment_curvature.kappa[-1] + kappa_inc
# find neutral axis that gives convergence of the axial force
try:
(d_n, r) = brentq(
f=self.service_normal_force_convergence,
a=a,
b=b,
args=(kappa, moment_curvature),
xtol=1e-3,
rtol=1e-6, # type: ignore
full_output=True,
disp=False,
)
except ValueError:
warnings.warn("brentq algorithm failed.")
m_xy = np.sqrt(
moment_curvature._m_x_i**2 + moment_curvature._m_y_i**2
)
text_update = "[red]Generating M-K diagram: "
text_update += f"M={m_xy:.3e}"
progress.update(task, description=text_update)
# save results
if not moment_curvature._failure:
moment_curvature.kappa.append(kappa)
moment_curvature.n.append(moment_curvature._n_i)
moment_curvature.m_x.append(moment_curvature._m_x_i)
moment_curvature.m_y.append(moment_curvature._m_y_i)
moment_curvature.m_xy.append(m_xy)
iter += 1
progress.update(
task,
description="[bold green]:white_check_mark: M-K diagram generated",
)
live.refresh()
return moment_curvature
[docs] def service_normal_force_convergence(
self,
d_n: float,
kappa: float,
moment_curvature: res.MomentCurvatureResults,
) -> float:
"""Given a neutral axis depth ``d_n`` and curvature ``kappa``, returns the the
net axial force.
:param d_nc: Trial cracked neutral axis
:param kappa: Curvature
:param moment_curvature: Moment curvature results object
:return: Net axial force
"""
# reset failure
moment_curvature._failure = False
# calculate extreme fibre in global coordinates
extreme_fibre, d_t = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=moment_curvature.theta
)
# validate d_n input
if d_n <= 0:
raise ValueError("d_n must be positive.")
elif d_n > d_t:
raise ValueError("d_n must lie within the section, i.e. d_n <= d_t")
# find point on neutral axis by shifting by d_n
point_na = utils.point_on_neutral_axis(
extreme_fibre=extreme_fibre, d_n=d_n, theta=moment_curvature.theta
)
# create splits in meshed geometries at points in stress-strain profiles
meshed_split_geoms: List[Union[CPGeom, CPGeomConcrete]] = []
for meshed_geom in self.meshed_geometries:
split_geoms = utils.split_geom_at_strains(
geom=meshed_geom,
theta=moment_curvature.theta,
point_na=point_na,
ultimate=False,
kappa=kappa,
)
meshed_split_geoms.extend(split_geoms)
# initialise results
n = 0
m_x = 0
m_y = 0
# calculate meshed geometry actions
for meshed_geom in meshed_split_geoms:
sec = AnalysisSection(geometry=meshed_geom)
n_sec, m_x_sec, m_y_sec, min_strain, max_strain = sec.service_analysis(
point_na=point_na,
theta=moment_curvature.theta,
kappa=kappa,
centroid=(self.gross_properties.cx, self.gross_properties.cy),
)
n += n_sec
m_x += m_x_sec
m_y += m_y_sec
# check for failure
ult_comp_strain = (
meshed_geom.material.stress_strain_profile.get_ultimate_compressive_strain()
)
ult_tens_strain = (
meshed_geom.material.stress_strain_profile.get_ultimate_tensile_strain()
)
# don't worry about tension failure in concrete
if max_strain > ult_comp_strain or (
min_strain < ult_tens_strain
and not isinstance(meshed_geom, CPGeomConcrete)
):
moment_curvature._failure = True
moment_curvature.failure_geometry = meshed_geom
# calculate lumped geometry actions
for lumped_geom in self.reinf_geometries_lumped:
# calculate area and centroid
area = lumped_geom.calculate_area()
centroid = lumped_geom.calculate_centroid()
# get strain at centroid of lump
strain = utils.get_service_strain(
point=(centroid[0], centroid[1]),
point_na=point_na,
theta=moment_curvature.theta,
kappa=kappa,
)
# check for failure
ult_comp_strain = (
lumped_geom.material.stress_strain_profile.get_ultimate_compressive_strain()
)
ult_tens_strain = (
lumped_geom.material.stress_strain_profile.get_ultimate_tensile_strain()
)
if strain > ult_comp_strain or strain < ult_tens_strain:
moment_curvature._failure = True
moment_curvature.failure_geometry = lumped_geom
# calculate stress and force
stress = lumped_geom.material.stress_strain_profile.get_stress(
strain=strain
)
force = stress * area
n += force
# calculate moment
m_x += force * (centroid[1] - self.gross_properties.cy)
m_y += force * (centroid[0] - self.gross_properties.cx)
moment_curvature._n_i = n
moment_curvature._m_x_i = m_x
moment_curvature._m_y_i = m_y
# calculate convergence
return n
[docs] def ultimate_bending_capacity(
self,
theta: float = 0,
n: float = 0,
) -> res.UltimateBendingResults:
r"""Given a neutral axis angle ``theta`` and an axial force ``n``, calculates
the ultimate bending capacity.
Note that ``k_u`` is calculated only for lumped (non-meshed) geometries.
:param theta: Angle (in radians) the neutral axis makes with the horizontal axis
(:math:`-\pi \leq \theta \leq \pi`)
:param n: Net axial force
:return: Ultimate bending results object
"""
# set neutral axis depth limits
# depth of neutral axis at extreme tensile fibre
_, d_t = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=theta
)
a = 1e-6 * d_t # sufficiently small depth of compressive zone
b = 6 * d_t # neutral axis at sufficiently large tensile fibre
# initialise ultimate bending results
ultimate_results = res.UltimateBendingResults(theta=theta)
# find neutral axis that gives convergence of the axial force
try:
(d_n, r) = brentq(
f=self.ultimate_normal_force_convergence,
a=a,
b=b,
args=(n, ultimate_results),
xtol=1e-3,
rtol=1e-6, # type: ignore
full_output=True,
disp=False,
)
except ValueError:
warnings.warn("brentq algorithm failed.")
return ultimate_results
[docs] def ultimate_normal_force_convergence(
self,
d_n: float,
n: float,
ultimate_results: res.UltimateBendingResults,
) -> float:
"""Given a neutral axis depth ``d_n`` and neutral axis angle ``theta``,
calculates the difference between the target net axial force ``n`` and the
calculated axial force.
:param d_n: Depth of the neutral axis from the extreme compression fibre
:param n: Net axial force
:param ultimate_results: Ultimate bending results object
:return: Axial force convergence
"""
# calculate convergence
return (
n
- self.calculate_ultimate_section_actions(
d_n=d_n, ultimate_results=ultimate_results
).n
)
[docs] def calculate_ultimate_section_actions(
self,
d_n: float,
ultimate_results: Optional[res.UltimateBendingResults] = None,
) -> res.UltimateBendingResults:
"""Given a neutral axis depth ``d_n`` and neutral axis angle ``theta``,
calculates the resultant bending moments ``m_x``, ``m_y``, ``m_xy`` and the net
axial force ``n``.
:param d_n: Depth of the neutral axis from the extreme compression fibre
:param ultimate_results: Ultimate bending results object
:return: Ultimate bending results object
"""
if ultimate_results is None:
ultimate_results = res.UltimateBendingResults(theta=0)
# calculate extreme fibre in global coordinates
extreme_fibre, _ = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=ultimate_results.theta
)
# extreme fibre in local coordinates
_, ef_v = utils.global_to_local(
theta=ultimate_results.theta,
x=extreme_fibre[0],
y=extreme_fibre[1],
)
# validate d_n input
if d_n <= 0:
raise ValueError("d_n must be positive.")
# find point on neutral axis by shifting by d_n
if isinf(d_n):
point_na = (0, 0)
else:
point_na = utils.point_on_neutral_axis(
extreme_fibre=extreme_fibre, d_n=d_n, theta=ultimate_results.theta
)
# create splits in meshed geometries at points in stress-strain profiles
meshed_split_geoms: List[Union[CPGeom, CPGeomConcrete]] = []
if isinf(d_n):
meshed_split_geoms = self.meshed_geometries
else:
for meshed_geom in self.meshed_geometries:
split_geoms = utils.split_geom_at_strains(
geom=meshed_geom,
theta=ultimate_results.theta,
point_na=point_na,
ultimate=True,
ultimate_strain=self.gross_properties.conc_ultimate_strain,
d_n=d_n,
)
meshed_split_geoms.extend(split_geoms)
# initialise results
n = 0
m_x = 0
m_y = 0
k_u = []
# calculate meshed geometry actions
for meshed_geom in meshed_split_geoms:
sec = AnalysisSection(geometry=meshed_geom)
n_sec, m_x_sec, m_y_sec = sec.ultimate_analysis(
point_na=point_na,
d_n=d_n,
theta=ultimate_results.theta,
ultimate_strain=self.gross_properties.conc_ultimate_strain,
centroid=(self.gross_properties.cx, self.gross_properties.cy),
)
n += n_sec
m_x += m_x_sec
m_y += m_y_sec
# calculate lumped actions
for lumped_geom in self.reinf_geometries_lumped:
# calculate area and centroid
area = lumped_geom.calculate_area()
centroid = lumped_geom.calculate_centroid()
# get strain at centroid of lump
if isinf(d_n):
strain = self.gross_properties.conc_ultimate_strain
else:
strain = utils.get_ultimate_strain(
point=(centroid[0], centroid[1]),
point_na=point_na,
d_n=d_n,
theta=ultimate_results.theta,
ultimate_strain=self.gross_properties.conc_ultimate_strain,
)
# calculate stress and force
stress = lumped_geom.material.stress_strain_profile.get_stress(
strain=strain
)
force = stress * area
n += force
# convert centroid to local coordinates
_, c_v = utils.global_to_local(
theta=ultimate_results.theta, x=centroid[0], y=centroid[1]
)
# calculate moment
m_x += force * (centroid[1] - self.gross_properties.cy)
m_y += force * (centroid[0] - self.gross_properties.cx)
# calculate k_u
d = ef_v - c_v
k_u.append(d_n / d)
# calculate resultant moment
m_xy = np.sqrt(m_x * m_x + m_y * m_y)
# save results
ultimate_results.d_n = d_n
ultimate_results.n = n
ultimate_results.m_x = m_x
ultimate_results.m_y = m_y
ultimate_results.m_xy = m_xy
if k_u:
ultimate_results.k_u = min(k_u)
return ultimate_results
[docs] def moment_interaction_diagram(
self,
theta: float = 0,
control_points: List[Tuple[str, float]] = [
("kappa0", 0.0),
("D", 1.0),
("fy", 1.0),
("N", 0.0),
("d_n", 1e-6),
],
labels: List[Union[str, None]] = [None],
n_points: Union[int, List[int]] = [4, 12, 12, 4],
max_comp: Optional[float] = None,
max_comp_labels: List[Union[str, None]] = [None, None],
) -> res.MomentInteractionResults:
r"""Generates a moment interaction diagram given a neutral axis angle `theta`
and `n_points` calculation points between the decompression case and the pure
bending case.
:param theta: Angle (in radians) the neutral axis makes with the horizontal axis
(:math:`-\pi \leq \theta \leq \pi`)
:param control_points: List of control points over which to generate the
interaction diagram. Each entry in ``control_points`` is a ``Tuple`` with
the first item the type of control point and the second item defining the
location of the control point. Acceptable types of control points are
``"D"`` (ratio of neutral axis depth to section depth), ``"d_n"`` (neutral
axis depth), ``"fy"`` (yield ratio of the most extreme tensile bar), ``"N"``
(axial force) and ``"kappa"`` (zero curvature compression - must be at start
of list, second value in tuple is not used). Control points must be defined
in an order which results in a decreasing neutral axis depth (decreasing
axial force). The default control points define an interaction diagram from
the decompression point to the pure bending point.
:param labels: List of labels to apply to the ``control_points`` for plotting
purposes, length must be the same as the length of ``control_points``. If a
single value is provided, will apply this label to all control points.
:param n_points: Number of neutral axis depths to compute between each control
point. Length must be one less than the length of ``control_points``. If an
integer is provided this will be used between all control points.
:param max_comp: If provided, limits the maximum compressive force in the moment
interaction diagram to ``max_comp``
:param max_comp_labels: Labels to apply to the ``max_comp`` intersection points,
first value is at zero moment, second value is at the intersection with the
interaction diagram
:raises ValueError: If ``control_points``, ``labels`` or ``n_points`` is invalid
:return: Moment interaction results object
"""
# if an integer is provided for n_points, generate a list
if isinstance(n_points, int):
n_points = [n_points] * (len(control_points) - 1)
# if there are no labels provided, generate a list
if len(labels) == 1:
labels = labels * len(control_points)
# validate n_points length
if len(n_points) != len(control_points) - 1:
raise ValueError(
"Length of n_points must be one less than the length of control_points."
)
# validate n_points entries are all longer than 2
for n_pt in n_points:
if n_pt < 3:
raise ValueError("n_points entries must be greater than 2.")
# validate labels length
if len(labels) != len(control_points):
raise ValueError("Length of labels must equal length of control_points.")
# initialise results
mi_results = res.MomentInteractionResults()
# compute extreme tensile fibre
_, d_t = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=theta
)
# function to decode d_n from control_point
def decode_d_n(cp):
# multiple of section depth
if cp[0] == "D":
# check D
if cp[1] <= 0:
raise ValueError(
f"Provided section depth (D) {cp[1]:.3f} must be greater than 0."
)
return cp[1] * d_t
# neutral axis depth
elif cp[0] == "d_n":
# check d_n
if cp[1] <= 0:
raise ValueError(
f"Provided d_n {cp[1]:.3f} must be greater than zero."
)
return cp[1]
# extreme tensile reinforcement yield ratio
elif cp[0] == "fy":
# get extreme tensile bar
d_ext, eps_sy = self.extreme_bar(theta=theta)
# get compressive strain at extreme fibre
eps_cu = self.gross_properties.conc_ultimate_strain
return d_ext * (eps_cu) / (cp[1] * eps_sy + eps_cu)
# provided axial force
elif cp[0] == "N":
ult_res = self.ultimate_bending_capacity(theta=theta, n=cp[1])
return ult_res.d_n
# zero curvature
elif cp[0] == "kappa0":
return 2 * d_t # sufficient depth to capture rectangular block
# control point type not valid
else:
raise ValueError(
"First value of control_point tuple must be D, d_n, fy, N or kappa0."
)
# see if a kappa0 was used
has_kappa0 = False
for cp in control_points:
if cp[0] == "kappa0":
has_kappa0 = True
break
# generate list of neutral axis depths to analyse and list of labels to save
d_n_list = []
label_list = []
start_d_n = 0
end_d_n = 0
for idx, n_pt in enumerate(n_points):
# get netural axis depths from control_points
start_d_n = decode_d_n(control_points[idx])
end_d_n = decode_d_n(control_points[idx + 1])
# generate list of neutral axis depths for this interval
d_n_list.extend(
np.linspace(
start=start_d_n, stop=end_d_n, num=n_pt - 1, endpoint=False
).tolist()
)
# add labels
label_list.append(labels[idx])
label_list.extend([None] * (n_pt - 2))
# add final d_n and label
d_n_list.append(end_d_n)
label_list.append(labels[-1])
# check d_n_list is ordered
if not all(d_n_list[i] >= d_n_list[i + 1] for i in range(len(d_n_list) - 1)):
msg = "control_points must create an ordered list of neutral axes from "
msg += "tensile fibre to compressive fibre."
raise ValueError(msg)
# create progress bar
progress = utils.create_known_progress()
with Live(progress, refresh_per_second=10) as live:
# add progress bar task
task = progress.add_task(
description="[red]Generating M-N diagram",
total=sum(n_points) - len(n_points) + 1,
)
# loop through all neutral axes
for idx, d_n in enumerate(d_n_list):
# calculate ultimate results
if idx == 0 and has_kappa0:
ult_res = self.calculate_ultimate_section_actions(
d_n=inf,
ultimate_results=res.UltimateBendingResults(theta=theta),
)
else:
ult_res = self.calculate_ultimate_section_actions(
d_n=d_n,
ultimate_results=res.UltimateBendingResults(theta=theta),
)
# add label
ult_res.label = label_list[idx]
# add ultimate result to moment interactions results and update progress
mi_results.results.append(ult_res)
progress.update(task, advance=1)
# display finished progress bar
progress.update(
task,
description="[bold green]:white_check_mark: M-N diagram generated",
)
live.refresh()
# cut diagram at max_comp
if max_comp:
# find intersection of max comp with interaction diagram
# and determine which points need to be removed from diagram
x = []
y_mx = []
y_my = []
y_mxy = []
idx_to_keep = 0
for idx, mi_res in enumerate(mi_results.results):
# create coordinates for interpolation
x.append(mi_res.n)
y_mx.append(mi_res.m_x)
y_my.append(mi_res.m_y)
y_mxy.append(mi_res.m_xy)
# determine which index is the first to keep
if idx_to_keep == 0 and mi_res.n < max_comp:
idx_to_keep = idx
# create interpolation function and determine moment which corresponds to
# an axial force of max_comp
f_mx = interp1d(x=x, y=y_mx)
f_my = interp1d(x=x, y=y_my)
f_mxy = interp1d(x=x, y=y_mxy)
m_max_comp_mx = f_mx(max_comp)
m_max_comp_my = f_my(max_comp)
m_max_comp_mxy = f_mxy(max_comp)
# remove points in diagram
del mi_results.results[:idx_to_keep]
# add first two points to diagram
# (m_max_comp, max_comp)
mi_results.results.insert(
0,
res.UltimateBendingResults(
theta=theta,
d_n=nan,
k_u=nan,
n=max_comp,
m_x=m_max_comp_mx,
m_y=m_max_comp_my,
m_xy=m_max_comp_mxy,
label=max_comp_labels[1],
),
)
# (0, max_comp)
mi_results.results.insert(
0,
res.UltimateBendingResults(
theta=theta,
d_n=inf,
k_u=0,
n=max_comp,
m_x=0,
m_y=0,
m_xy=0,
label=max_comp_labels[0],
),
)
return mi_results
[docs] def biaxial_bending_diagram(
self,
n: float = 0,
n_points: int = 48,
) -> res.BiaxialBendingResults:
"""Generates a biaxial bending diagram given a net axial force ``n`` and
``n_points`` calculation points.
:param n: Net axial force
:param n_points: Number of calculation points between the decompression
:return: Biaxial bending results
"""
# initialise results
bb_results = res.BiaxialBendingResults(n=n)
# calculate d_theta
d_theta = 2 * np.pi / n_points
# generate list of thetas
theta_list = np.linspace(start=-np.pi, stop=np.pi - d_theta, num=n_points)
# create progress bar
progress = utils.create_known_progress()
with Live(progress, refresh_per_second=10) as live:
task = progress.add_task(
description="[red]Generating biaxial bending diagram",
total=n_points,
)
# loop through thetas
for theta in theta_list:
ultimate_results = self.ultimate_bending_capacity(theta=theta, n=n)
bb_results.results.append(ultimate_results)
progress.update(task, advance=1)
# add first result to end of list top
bb_results.results.append(bb_results.results[0])
progress.update(
task,
description="[bold green]:white_check_mark: Biaxial bending diagram generated",
)
live.refresh()
return bb_results
[docs] def calculate_uncracked_stress(
self,
n: float = 0,
m_x: float = 0,
m_y: float = 0,
) -> res.StressResult:
"""Calculates stresses within the reinforced concrete section assuming an
uncracked section.
Uses gross area section properties to determine concrete and reinforcement
stresses given an axial force ``n``, and bending moments ``m_x`` and ``m_y``.
:param n: Axial force
:param m_x: Bending moment about the x-axis
:param m_y: Bending moment about the y-axis
:return: Stress results object
"""
# initialise stress results
conc_sections = []
conc_sigs = []
conc_forces = []
meshed_reinf_sections = []
meshed_reinf_sigs = []
meshed_reinf_forces = []
lumped_reinf_geoms = []
lumped_reinf_sigs = []
lumped_reinf_strains = []
lumped_reinf_forces = []
# get uncracked section properties
e_a = self.gross_properties.e_a
cx = self.gross_properties.cx
cy = self.gross_properties.cy
e_ixx = self.gross_properties.e_ixx_c
e_iyy = self.gross_properties.e_iyy_c
e_ixy = self.gross_properties.e_ixy_c
# calculate neutral axis rotation
grad = (e_ixy * m_x - e_ixx * m_y) / (e_iyy * m_x - e_ixy * m_y)
theta = np.arctan2(grad, 1)
if np.isclose(theta, 0):
theta = 0
# point on neutral axis is centroid
point_na = (cx, cy)
# split meshed geometries above and below neutral axis
split_meshed_geoms = []
for meshed_geom in self.meshed_geometries:
top_geoms, bot_geoms = meshed_geom.split_section(
point=point_na,
theta=theta,
)
split_meshed_geoms.extend(top_geoms)
split_meshed_geoms.extend(bot_geoms)
# loop through all meshed geometries and calculate stress
for meshed_geom in split_meshed_geoms:
analysis_section = AnalysisSection(geometry=meshed_geom)
# calculate stress, force and point of action
sig, n_sec, d_x, d_y = analysis_section.get_elastic_stress(
n=n,
m_x=m_x,
m_y=m_y,
e_a=e_a,
cx=cx,
cy=cy,
e_ixx=e_ixx,
e_iyy=e_iyy,
e_ixy=e_ixy,
)
# save results
if isinstance(meshed_geom, CPGeomConcrete):
conc_sigs.append(sig)
conc_forces.append((n_sec, d_x, d_y))
conc_sections.append(analysis_section)
else:
meshed_reinf_sigs.append(sig)
meshed_reinf_forces.append((n_sec, d_x, d_y))
meshed_reinf_sections.append(analysis_section)
# loop through all lumped geometries and calculate stress
for lumped_geom in self.reinf_geometries_lumped:
# initialise stress and position
sig = 0
centroid = lumped_geom.calculate_centroid()
x = centroid[0] - cx
y = centroid[1] - cy
# axial stress
sig += n * lumped_geom.material.elastic_modulus / e_a
# bending moment stress
sig += lumped_geom.material.elastic_modulus * (
-(e_ixy * m_x) / (e_ixx * e_iyy - e_ixy**2) * x
+ (e_iyy * m_x) / (e_ixx * e_iyy - e_ixy**2) * y
)
sig += lumped_geom.material.elastic_modulus * (
+(e_ixx * m_y) / (e_ixx * e_iyy - e_ixy**2) * x
- (e_ixy * m_y) / (e_ixx * e_iyy - e_ixy**2) * y
)
strain = sig / lumped_geom.material.elastic_modulus
# net force and point of action
n_lumped = sig * lumped_geom.calculate_area()
lumped_reinf_sigs.append(sig)
lumped_reinf_strains.append(strain)
lumped_reinf_forces.append((n_lumped, x, y))
lumped_reinf_geoms.append(lumped_geom)
return res.StressResult(
concrete_section=self,
concrete_analysis_sections=conc_sections,
concrete_stresses=conc_sigs,
concrete_forces=conc_forces,
meshed_reinforcement_sections=meshed_reinf_sections,
meshed_reinforcement_stresses=meshed_reinf_sigs,
meshed_reinforcement_forces=meshed_reinf_forces,
lumped_reinforcement_geometries=lumped_reinf_geoms,
lumped_reinforcement_stresses=lumped_reinf_sigs,
lumped_reinforcement_strains=lumped_reinf_strains,
lumped_reinforcement_forces=lumped_reinf_forces,
)
[docs] def calculate_cracked_stress(
self,
cracked_results: res.CrackedResults,
n: float = 0,
m: float = 0,
) -> res.StressResult:
"""Calculates stresses within the reinforced concrete section assuming a cracked
section.
Uses cracked area section properties to determine concrete and reinforcement
stresses given an axial force ``n`` and bending moment ``m`` about the bending
axis stored in ``cracked_results``.
:param cracked_results: Cracked results objects
:param n: Axial force
:param m: Bending moment
:return: Stress results object
"""
# initialise stress results
conc_sections = []
conc_sigs = []
conc_forces = []
meshed_reinf_sections = []
meshed_reinf_sigs = []
meshed_reinf_forces = []
lumped_reinf_geoms = []
lumped_reinf_sigs = []
lumped_reinf_strains = []
lumped_reinf_forces = []
# get cracked section properties
e_a = cracked_results.e_a_cr
cx = cracked_results.cx
cy = cracked_results.cy
e_ixx = cracked_results.e_ixx_c_cr
e_iyy = cracked_results.e_iyy_c_cr
e_ixy = cracked_results.e_ixy_c_cr
# correct small e_ixy sign error
if abs(e_ixy / cracked_results.e_i11_cr) < 1e-12:
e_ixy = 0
# get bending angle
theta = cracked_results.theta
# handle cardinal points (avoid divide by zeros)
tan_theta = np.tan(theta)
with np.errstate(divide="ignore"):
c = (e_ixx - e_ixy * tan_theta) / (e_ixy - e_iyy * tan_theta)
# calculate bending moment about each axis (figure out signs)
sign = 0
if theta <= 0:
if c < 0:
sign = -1
elif c > 0:
sign = 1
else:
if c < 0:
sign = 1
else:
sign = -1
m_x = sign * np.sqrt(m * m / (1 + 1 / (c * c)))
m_y = m_x / c
# loop through all meshed geometries and calculate stress
for geom in cracked_results.cracked_geometries:
if geom.material.meshed:
analysis_section = AnalysisSection(geometry=geom)
# calculate stress, force and point of action
sig, n_sec, d_x, d_y = analysis_section.get_elastic_stress(
n=n,
m_x=m_x,
m_y=m_y,
e_a=e_a,
cx=cx,
cy=cy,
e_ixx=e_ixx,
e_iyy=e_iyy,
e_ixy=e_ixy,
)
# save results
if isinstance(geom, CPGeomConcrete):
conc_sigs.append(sig)
conc_forces.append((n_sec, d_x, d_y))
conc_sections.append(analysis_section)
else:
meshed_reinf_sigs.append(sig)
meshed_reinf_forces.append((n_sec, d_x, d_y))
meshed_reinf_sections.append(analysis_section)
# loop through all lumped geometries and calculate stress
for lumped_geom in self.reinf_geometries_lumped:
# initialise stress and position of bar
sig = 0
centroid = lumped_geom.calculate_centroid()
x = centroid[0] - cx
y = centroid[1] - cy
# axial stress
sig += n * lumped_geom.material.elastic_modulus / e_a
# bending moment stress
sig += lumped_geom.material.elastic_modulus * (
-(e_ixy * m_x) / (e_ixx * e_iyy - e_ixy**2) * x
+ (e_iyy * m_x) / (e_ixx * e_iyy - e_ixy**2) * y
)
sig += lumped_geom.material.elastic_modulus * (
+(e_ixx * m_y) / (e_ixx * e_iyy - e_ixy**2) * x
- (e_ixy * m_y) / (e_ixx * e_iyy - e_ixy**2) * y
)
strain = sig / lumped_geom.material.elastic_modulus
# net force and point of action
n_lumped = sig * lumped_geom.calculate_area()
lumped_reinf_sigs.append(sig)
lumped_reinf_strains.append(strain)
lumped_reinf_forces.append((n_lumped, x, y))
lumped_reinf_geoms.append(lumped_geom)
return res.StressResult(
concrete_section=self,
concrete_analysis_sections=conc_sections,
concrete_stresses=conc_sigs,
concrete_forces=conc_forces,
meshed_reinforcement_sections=meshed_reinf_sections,
meshed_reinforcement_stresses=meshed_reinf_sigs,
meshed_reinforcement_forces=meshed_reinf_forces,
lumped_reinforcement_geometries=lumped_reinf_geoms,
lumped_reinforcement_stresses=lumped_reinf_sigs,
lumped_reinforcement_strains=lumped_reinf_strains,
lumped_reinforcement_forces=lumped_reinf_forces,
)
[docs] def calculate_service_stress(
self,
moment_curvature_results: res.MomentCurvatureResults,
m: float,
kappa: Optional[float] = None,
) -> res.StressResult:
"""Calculates service stresses within the reinforced concrete section.
Uses linear interpolation of the moment-curvature results to determine the
curvature of the section given the user supplied moment, and thus the stresses
within the section. Otherwise, a curvature can be provided which overrides the
supplied moment.
:param moment_curvature_results: Moment-curvature results objects
:param m: Bending moment
:param kappa: Curvature, if provided overrides the supplied bending moment and
calculates the stress at the given curvature
:return: Stress results object
"""
if kappa is None:
# get curvature
kappa = moment_curvature_results.get_curvature(moment=m)
# get theta
theta = moment_curvature_results.theta
# initialise variables
mk = res.MomentCurvatureResults(theta=theta)
# set neutral axis depth limits
# depth of neutral axis at extreme tensile fibre
extreme_fibre, d_t = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=theta
)
a = 1e-6 * d_t # sufficiently small depth of compressive zone
b = d_t # neutral axis at extreme tensile fibre
# find neutral axis that gives convergence of the axial force
try:
(d_n, r) = brentq(
f=self.service_normal_force_convergence,
a=a,
b=b,
args=(kappa, mk),
xtol=1e-3,
rtol=1e-6, # type: ignore
full_output=True,
disp=False,
)
except ValueError:
warnings.warn("brentq algorithm failed.")
d_n = 0
# initialise stress results
conc_sections = []
conc_sigs = []
conc_forces = []
meshed_reinf_sections = []
meshed_reinf_sigs = []
meshed_reinf_forces = []
lumped_reinf_geoms = []
lumped_reinf_sigs = []
lumped_reinf_strains = []
lumped_reinf_forces = []
# find point on neutral axis by shifting by d_n
point_na = utils.point_on_neutral_axis(
extreme_fibre=extreme_fibre, d_n=d_n, theta=theta
)
# create splits in meshed geometries at points in stress-strain profiles
meshed_split_geoms: List[Union[CPGeom, CPGeomConcrete]] = []
for meshed_geom in self.meshed_geometries:
split_geoms = utils.split_geom_at_strains(
geom=meshed_geom,
theta=theta,
point_na=point_na,
ultimate=False,
kappa=kappa,
)
meshed_split_geoms.extend(split_geoms)
# loop through all meshed geometries and calculate stress
for meshed_geom in meshed_split_geoms:
analysis_section = AnalysisSection(geometry=meshed_geom)
# calculate stress, force and point of action
sig, n_sec, d_x, d_y = analysis_section.get_service_stress(
d_n=d_n,
kappa=kappa,
point_na=point_na,
theta=theta,
centroid=(self.gross_properties.cx, self.gross_properties.cy),
)
# save results
if isinstance(meshed_geom, CPGeomConcrete):
conc_sigs.append(sig)
conc_forces.append((n_sec, d_x, d_y))
conc_sections.append(analysis_section)
else:
meshed_reinf_sigs.append(sig)
meshed_reinf_forces.append((n_sec, d_x, d_y))
meshed_reinf_sections.append(analysis_section)
# loop through all lumped geometries and calculate stress
for lumped_geom in self.reinf_geometries_lumped:
# get position of geometry
centroid = lumped_geom.calculate_centroid()
# get strain at centroid of lump
strain = utils.get_service_strain(
point=(centroid[0], centroid[1]),
point_na=point_na,
theta=theta,
kappa=kappa,
)
# calculate stress, force and point of action
sig = lumped_geom.material.stress_strain_profile.get_stress(strain=strain)
n_lumped = sig * lumped_geom.calculate_area()
lumped_reinf_sigs.append(sig)
lumped_reinf_strains.append(strain)
lumped_reinf_forces.append(
(
n_lumped,
centroid[0] - self.gross_properties.cx,
centroid[1] - self.gross_properties.cy,
)
)
lumped_reinf_geoms.append(lumped_geom)
return res.StressResult(
concrete_section=self,
concrete_analysis_sections=conc_sections,
concrete_stresses=conc_sigs,
concrete_forces=conc_forces,
meshed_reinforcement_sections=meshed_reinf_sections,
meshed_reinforcement_stresses=meshed_reinf_sigs,
meshed_reinforcement_forces=meshed_reinf_forces,
lumped_reinforcement_geometries=lumped_reinf_geoms,
lumped_reinforcement_stresses=lumped_reinf_sigs,
lumped_reinforcement_strains=lumped_reinf_strains,
lumped_reinforcement_forces=lumped_reinf_forces,
)
[docs] def calculate_ultimate_stress(
self,
ultimate_results: res.UltimateBendingResults,
) -> res.StressResult:
"""Calculates ultimate stresses within the reinforced concrete section.
:param ultimate_results: Ultimate bending results objects
:return: Stress results object
"""
# depth of neutral axis at extreme tensile fibre
extreme_fibre, _ = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=ultimate_results.theta
)
# find point on neutral axis by shifting by d_n
if isinf(ultimate_results.d_n):
point_na = (0, 0)
else:
point_na = utils.point_on_neutral_axis(
extreme_fibre=extreme_fibre,
d_n=ultimate_results.d_n,
theta=ultimate_results.theta,
)
# initialise stress results for each concrete geometry
conc_sections = []
conc_sigs = []
conc_forces = []
meshed_reinf_sections = []
meshed_reinf_sigs = []
meshed_reinf_forces = []
lumped_reinf_geoms = []
lumped_reinf_sigs = []
lumped_reinf_strains = []
lumped_reinf_forces = []
# create splits in meshed geometries at points in stress-strain profiles
meshed_split_geoms: List[Union[CPGeom, CPGeomConcrete]] = []
if isinf(ultimate_results.d_n):
meshed_split_geoms = self.meshed_geometries
else:
for meshed_geom in self.meshed_geometries:
split_geoms = utils.split_geom_at_strains(
geom=meshed_geom,
theta=ultimate_results.theta,
point_na=point_na,
ultimate=True,
ultimate_strain=self.gross_properties.conc_ultimate_strain,
d_n=ultimate_results.d_n,
)
meshed_split_geoms.extend(split_geoms)
# loop through all concrete geometries and calculate stress
for meshed_geom in meshed_split_geoms:
analysis_section = AnalysisSection(geometry=meshed_geom)
# calculate stress, force and point of action
sig, n_sec, d_x, d_y = analysis_section.get_ultimate_stress(
d_n=ultimate_results.d_n,
point_na=point_na,
theta=ultimate_results.theta,
ultimate_strain=self.gross_properties.conc_ultimate_strain,
centroid=(self.gross_properties.cx, self.gross_properties.cy),
)
# save results
if isinstance(meshed_geom, CPGeomConcrete):
conc_sigs.append(sig)
conc_forces.append((n_sec, d_x, d_y))
conc_sections.append(analysis_section)
else:
meshed_reinf_sigs.append(sig)
meshed_reinf_forces.append((n_sec, d_x, d_y))
meshed_reinf_sections.append(analysis_section)
# loop through all lumped geometries and calculate stress
for lumped_geom in self.reinf_geometries_lumped:
# get position of lump
centroid = lumped_geom.calculate_centroid()
# get strain at centroid of lump
if isinf(ultimate_results.d_n):
strain = self.gross_properties.conc_ultimate_strain
else:
strain = utils.get_ultimate_strain(
point=(centroid[0], centroid[1]),
point_na=point_na,
d_n=ultimate_results.d_n,
theta=ultimate_results.theta,
ultimate_strain=self.gross_properties.conc_ultimate_strain,
)
# calculate stress, force and point of action
sig = lumped_geom.material.stress_strain_profile.get_stress(strain=strain)
n_lumped = sig * lumped_geom.calculate_area()
lumped_reinf_sigs.append(sig)
lumped_reinf_strains.append(strain)
lumped_reinf_forces.append(
(
n_lumped,
centroid[0] - self.gross_properties.cx,
centroid[1] - self.gross_properties.cy,
)
)
lumped_reinf_geoms.append(lumped_geom)
return res.StressResult(
concrete_section=self,
concrete_analysis_sections=conc_sections,
concrete_stresses=conc_sigs,
concrete_forces=conc_forces,
meshed_reinforcement_sections=meshed_reinf_sections,
meshed_reinforcement_stresses=meshed_reinf_sigs,
meshed_reinforcement_forces=meshed_reinf_forces,
lumped_reinforcement_geometries=lumped_reinf_geoms,
lumped_reinforcement_stresses=lumped_reinf_sigs,
lumped_reinforcement_strains=lumped_reinf_strains,
lumped_reinforcement_forces=lumped_reinf_forces,
)
[docs] def extreme_bar(
self,
theta: float,
) -> Tuple[float, float]:
r"""Given neutral axis angle ``theta``, determines the depth of the furthest
lumped reinforcement from the extreme compressive fibre and also returns its
yield strain.
:param theta: Angle (in radians) the neutral axis makes with the horizontal
axis (:math:`-\pi \leq \theta \leq \pi`)
:return: Depth of furthest bar and its yield strain
"""
# initialise variables
d_ext = 0
# calculate extreme fibre in local coordinates
extreme_fibre, _ = utils.calculate_extreme_fibre(
points=self.compound_geometry.points, theta=theta
)
_, ef_v = utils.global_to_local(
theta=theta, x=extreme_fibre[0], y=extreme_fibre[1]
)
# get depth to extreme lumped reinforcement
extreme_geom = self.reinf_geometries_lumped[0]
for lumped_geom in self.reinf_geometries_lumped:
centroid = lumped_geom.calculate_centroid()
# convert centroid to local coordinates
_, c_v = utils.global_to_local(theta=theta, x=centroid[0], y=centroid[1])
# calculate d
d = ef_v - c_v
if d > d_ext:
d_ext = d
extreme_geom = lumped_geom
# calculate yield strain
yield_strain = (
extreme_geom.material.stress_strain_profile.get_yield_strength()
/ extreme_geom.material.stress_strain_profile.get_elastic_modulus()
)
return d_ext, yield_strain
[docs] def plot_section(
self,
title: str = "Reinforced Concrete Section",
background: bool = False,
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots the reinforced concrete section.
:param title: Plot title
:param background: If set to True, uses the plot as a background plot
:param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
:return: Matplotlib axes object
"""
with plotting_context(title=title, aspect=True, **kwargs) as (fig, ax):
# create list of already plotted materials
plotted_materials = []
legend_labels = []
# plot meshed geometries
for meshed_geom in self.meshed_geometries:
if meshed_geom.material not in plotted_materials:
patch = mpatches.Patch(
color=meshed_geom.material.colour,
label=meshed_geom.material.name,
)
legend_labels.append(patch)
plotted_materials.append(meshed_geom.material)
# TODO - when shapely implements polygon plotting, fix this up
sec = AnalysisSection(geometry=meshed_geom)
if not background:
sec.plot_shape(ax=ax)
# plot the points and facets
for f in meshed_geom.facets:
if background:
fmt = "k-"
else:
fmt = "ko-"
ax.plot( # type: ignore
[meshed_geom.points[f[0]][0], meshed_geom.points[f[1]][0]],
[meshed_geom.points[f[0]][1], meshed_geom.points[f[1]][1]],
fmt,
markersize=2,
linewidth=1.5,
)
# plot lumped geometries
for lumped_geom in self.reinf_geometries_lumped:
if lumped_geom.material not in plotted_materials:
patch = mpatches.Patch(
color=lumped_geom.material.colour,
label=lumped_geom.material.name,
)
legend_labels.append(patch)
plotted_materials.append(lumped_geom.material)
# plot the points and facets
coords = list(lumped_geom.geom.exterior.coords) # type: ignore
bar = mpatches.Polygon(
xy=coords, closed=False, color=lumped_geom.material.colour
)
ax.add_patch(bar) # type: ignore
if not background:
ax.legend( # type: ignore
loc="center left", bbox_to_anchor=(1, 0.5), handles=legend_labels
)
return ax