Source code for concreteproperties.concrete_section

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, SteelStrand
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, moment_centroid: Optional[Tuple[float, float]] = None, geometric_centroid_override: bool = False, ) -> None: """Inits the ConcreteSection class. :param geometry: *sectionproperties* CompoundGeometry object describing the reinforced concrete section :param moment_centroid: If specified, all moments for service and ultimate analyses are calculated about this point. If not specified, all moments are calculated about the gross cross-section centroid, i.e. no material properties applied. :param geometric_centroid_override: If set to True, sets ``moment_centroid`` to the geometric centroid i.e. material properties applied (useful for composite section analysis) """ 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, reinforcement (meshed and lumped) and strand 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] = [] self.strand_geometries: 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) elif isinstance(geom.material, SteelStrand): # SteelStrand materials can only be in PrestressedSection if type(self) == ConcreteSection: raise ValueError( "SteelStrand material detected. Use PrestressedSection instead." ) else: cp_geom = CPGeom(geom=geom.geom, material=geom.material) self.strand_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) # calculate gross properties self.gross_properties = res.GrossProperties() self.calculate_gross_area_properties() # set moment centroid if moment_centroid: self.moment_centroid = moment_centroid else: self.moment_centroid = ( self.gross_properties.cx_gross, self.gross_properties.cy_gross, ) # if moment centroid overriden if geometric_centroid_override: self.moment_centroid = self.gross_properties.cx, self.gross_properties.cy
[docs] def calculate_gross_area_properties(self) -> None: """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] ) self.gross_properties.qx_gross += area * centroid[1] self.gross_properties.qy_gross += area * 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 ) self.gross_properties.cx_gross = ( self.gross_properties.qy_gross / self.gross_properties.total_area ) self.gross_properties.cy_gross = ( self.gross_properties.qx_gross / self.gross_properties.total_area ) # 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 + self.strand_geometries: # 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 get_transformed_gross_properties( self, elastic_modulus: float, ) -> res.TransformedGrossProperties: """Transforms gross section properties given a reference elastic modulus. :param elastic_modulus: Reference elastic modulus :return: Transformed concrete properties object """ return res.TransformedGrossProperties( concrete_properties=self.gross_properties, elastic_modulus=elastic_modulus )
[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: msg = "Analysis failed. Please raise an issue at " msg += "https://github.com/robbievanleeuwen/concrete-properties/issues" raise utils.AnalysisError(msg) # calculate cracked section properties self.cracked_section_properties(cracked_results=cracked_results) 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 cracked_section_properties( self, cracked_results: res.CrackedResults, ) -> None: """Given a list of cracked geometries (stored in ``cracked_results``), determines the cracked section properties and stores in ``cracked_results``. :param cracked_results: Cracked results object with stored cracked geometries """ # reset results cracked_results.reset_results() # 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(cracked_results.theta)) ** 2 + cracked_results.e_ixx_c_cr * (np.cos(cracked_results.theta)) ** 2 - 2 * cracked_results.e_ixy_c_cr * np.sin(cracked_results.theta) * np.cos(cracked_results.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, )
[docs] def moment_curvature_analysis( self, theta: float = 0, n: float = 0, kappa0: float = 0, kappa_inc: float = 1e-7, kappa_mult: float = 2, kappa_inc_max: float = 5e-6, delta_m_min: float = 0.15, delta_m_max: float = 0.3, progress_bar: bool = True, ) -> res.MomentCurvatureResults: r"""Performs a moment curvature analysis given a bending angle ``theta`` and applied axial force ``n``. Analysis continues until a material reaches its ultimate strain. :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) :param n: Axial force :param kappa0: Initial curvature :param kappa_inc: Initial curvature increment :param kappa_mult: Multiplier to apply to the curvature increment ``kappa_inc`` when ``delta_m_max`` is satisfied. When ``delta_m_min`` is satisfied, the inverse of this multipler is applied to ``kappa_inc``. :param kappa_inc_max: Maximum curvature increment :param delta_m_min: Relative change in moment at which to reduce the curvature increment :param delta_m_max: Relative change in moment at which to increase the curvature increment :param progress_bar: If set to True, displays the progress bar :return: Moment curvature results object """ # initialise variables moment_curvature = res.MomentCurvatureResults(theta=theta, n_target=n) # function that performs moment curvature analysis def mcurve(kappa_inc=kappa_inc, progress=None): iter = 0 kappa = kappa0 while not moment_curvature._failure: # calculate adaptive step size for curvature if iter > 2: moment_diff = ( abs(moment_curvature.kappa[-1] - moment_curvature.kappa[-2]) / moment_curvature.kappa[-1] ) if moment_diff <= delta_m_min: kappa_inc *= kappa_mult elif moment_diff >= delta_m_max: kappa_inc *= 1 / kappa_mult # enforce maximum curvature increment if kappa_inc > kappa_inc_max: kappa_inc = kappa_inc_max # update curvature kappa = kappa0 if iter == 0 else moment_curvature.kappa[-1] + kappa_inc # find neutral axis that gives convergence of the axial force try: brentq( f=self.service_normal_force_convergence, a=-0.1, b=0.1, args=(kappa, moment_curvature), full_output=True, disp=False, ) except ValueError: if not moment_curvature._failure: msg = "Analysis failed. Please raise an issue at " msg += "https://github.com/robbievanleeuwen/concrete-properties/issues" raise utils.AnalysisError(msg) m_xy = np.sqrt( moment_curvature._m_x_i**2 + moment_curvature._m_y_i**2 ) if progress: 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) moment_curvature.convergence.append( moment_curvature._failure_convergence ) iter += 1 # find kappa corresponding to failure strain: # curvature before and after failure kappa_a = moment_curvature.kappa[-1] kappa_b = kappa # this method (given a kappa) outputs the failure convergence # (normalised to zero) def failure_kappa(kappa_fail): # given kappa find equilibrium brentq( f=self.service_normal_force_convergence, a=-0.1, b=0.1, args=(kappa_fail, moment_curvature), full_output=True, disp=False, ) return moment_curvature._failure_convergence - 1 if progress: progress.update(task, description="[red]Finding failure curvature...") # find curvature corresponding to failure brentq( f=failure_kappa, a=kappa_a, b=kappa_b, full_output=True, disp=False, ) # save final results m_xy = np.sqrt(moment_curvature._m_x_i**2 + moment_curvature._m_y_i**2) moment_curvature.kappa.append(moment_curvature._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) moment_curvature.convergence.append(moment_curvature._failure_convergence) # create progress bar if progress_bar: # 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, ) mcurve(progress=progress) progress.update( task, description="[bold green]:white_check_mark: M-K diagram generated", ) live.refresh() else: mcurve() return moment_curvature
[docs] def service_normal_force_convergence( self, eps0: 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 eps0: Strain at top fibre :param kappa: Curvature :param moment_curvature: Moment curvature results object :return: Net axial force """ # reset failure moment_curvature._failure = False # get global coordinates of extreme compressive fibre ecf, _ = utils.calculate_extreme_fibre( points=self.compound_geometry.points, 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_service( geom=meshed_geom, theta=moment_curvature.theta, ecf=ecf, eps0=eps0, kappa=kappa, ) meshed_split_geoms.extend(split_geoms) # initialise results n = 0 m_x = 0 m_y = 0 failure_convergence = 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( ecf=ecf, eps0=eps0, theta=moment_curvature.theta, kappa=kappa, centroid=self.moment_centroid, ) 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 # update failure convergence # compression failure failure_convergence = max(max_strain / ult_comp_strain, failure_convergence) # tensile failure (ignore concrete) if not isinstance(meshed_geom, CPGeomConcrete): failure_convergence = max( min_strain / ult_tens_strain, failure_convergence ) # calculate lumped geometry actions for lumped_geom in self.reinf_geometries_lumped + self.strand_geometries: # 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]), ecf=ecf, eps0=eps0, theta=moment_curvature.theta, kappa=kappa, ) # add initial prestress strain if isinstance(lumped_geom.material, SteelStrand): eps_pe = -lumped_geom.material.get_prestress_strain(area=area) strain += eps_pe # 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 # update failure convergence # compression failure failure_convergence = max(strain / ult_comp_strain, failure_convergence) # tensile failure failure_convergence = max(strain / ult_tens_strain, failure_convergence) # 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.moment_centroid[1]) m_y += force * (centroid[0] - self.moment_centroid[0]) moment_curvature._kappa = kappa moment_curvature._n_i = n moment_curvature._m_x_i = m_x moment_curvature._m_y_i = m_y moment_curvature._failure_convergence = failure_convergence # return normal force convergence return n - moment_curvature.n_target
[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:: This calculation is code agnostic and no capacity reduction factors are applied. If design capacities are required, use the applicable ``design_code`` module or consult your local design code on how to treat nominal axial loads in ultimate bending calculations. 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 (nominal axial load) :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: msg = "Analysis failed. The solver could not find a neutral axis that " msg += "satisfies equilibrium. This may be due to an axial force that " msg += "exceeds the tensile or compressive capacity of the cross-section." raise utils.AnalysisError(msg) 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_ultimate( geom=meshed_geom, theta=ultimate_results.theta, point_na=point_na, 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.moment_centroid, ) n += n_sec m_x += m_x_sec m_y += m_y_sec # calculate lumped actions for lumped_geom in self.reinf_geometries_lumped + self.strand_geometries: # 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, ) # add initial prestress strain (N.B. ignore eps_ce) if isinstance(lumped_geom.material, SteelStrand): eps_pe = -lumped_geom.material.get_prestress_strain(area=area) strain += eps_pe # 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.moment_centroid[1]) m_y += force * (centroid[0] - self.moment_centroid[0]) # 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, limits: List[Tuple[str, float]] = [ ("D", 1.0), ("d_n", 1e-6), ], control_points: List[Tuple[str, float]] = [ ("kappa0", 0.0), ("fy", 1.0), ("N", 0.0), ], labels: Optional[List[str]] = None, n_points: int = 24, n_spacing: Optional[int] = None, max_comp: Optional[float] = None, max_comp_labels: Optional[List[str]] = None, progress_bar: bool = True, ) -> res.MomentInteractionResults: r"""Generates a moment interaction diagram given a neutral axis angle ``theta``. ``limits`` and ``control_points`` accept a list of tuples that define points on the moment interaction diagram. The first item in the tuple defines the type of control points, while the second item defines the location of the control point. Types of control points are detailed below: .. admonition:: Control points - ``"D"`` - ratio of neutral axis depth to section depth - ``"d_n"`` - neutral axis depth - ``"fy"`` - yield ratio of the most extreme tensile bar - ``"N"`` - net (nominal) axial force - ``"kappa0"`` - zero curvature compression (N.B second item in tuple is not used) :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) :param limits: List of control points that define the start and end of the interaction diagram. List length must equal two. The default limits range from concrete decompression strain to zero curvature tension. :param control_points: List of additional control points to add to the moment interatction diagram. The default control points include the pure compression point (``kappa0``), the balanced point (``fy=1``) and the pure bending point (``N=0``). Control points may lie outside the limits of the moment interaction diagram as long as equilibrium can be found. :param labels: List of labels to apply to the ``limits`` and ``control_points`` for plotting purposes. The first two values in ``labels`` apply labels to the ``limits``, the remaining values apply labels to the ``control_points``. If a single value is provided, this value will be applied to both ``limits`` and all ``control_points``. The length of ``labels`` must equal ``1`` or ``2 + len(control_points)``. :param n_points: Number of points to compute including and between the ``limits`` of the moment interaction diagram. Generates equally spaced neutral axis depths between the ``limits``. :param n_spacing: If provided, overrides ``n_points`` and generates the moment interaction diagram using ``n_spacing`` equally spaced axial loads. Note that using ``n_spacing`` negatively affects performance, as the neutral axis depth must first be located for each point on the moment interaction diagram. :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 :param progress_bar: If set to True, displays the progress bar :return: Moment interaction results object """ # compute extreme tensile fibre _, d_t = utils.calculate_extreme_fibre( points=self.compound_geometry.points, theta=theta ) # validate limits length if len(limits) != 2: raise ValueError("Length of limits must equal 2.") # get neutral axis depths for limits limits_dn = [] for cp in limits: limits_dn.append(self.decode_d_n(theta=theta, cp=cp, d_t=d_t)) # get neutral axis depths for additional control points add_cp_dn = [] for cp in control_points: add_cp_dn.append(self.decode_d_n(theta=theta, cp=cp, d_t=d_t)) # validate labels length if labels and len(labels) != 1 and len(labels) != 2 + len(control_points): raise ValueError( "Length of labels must be 1 or 2 + number of control points" ) # if one label is provided, generate a list if labels and len(labels) == 1: labels = labels * (len(control_points) + 2) # initialise results mi_results = res.MomentInteractionResults() # generate list of neutral axis depths/axial forces to analyse # if we are spacing by axial force if n_spacing: # get axial force of the limits start_res = self.calculate_ultimate_section_actions( d_n=limits_dn[0], ultimate_results=res.UltimateBendingResults(theta=theta), ) end_res = self.calculate_ultimate_section_actions( d_n=limits_dn[1], ultimate_results=res.UltimateBendingResults(theta=theta), ) # generate list of axial forces analysis_list = np.linspace( start=start_res.n, stop=end_res.n, num=n_spacing, dtype=float ).tolist() else: # check for infinity in limits - this will not work with linspace # for sake of distributing neutral axes let kappa0 ~= 2 * D if limits_dn[0] == inf: start = 2 * d_t else: start = limits_dn[0] if limits_dn[1] == inf: stop = 2 * d_t else: stop = limits_dn[1] # generate list of neutral axes analysis_list = np.linspace( start=start, stop=stop, num=n_points, dtype=float ).tolist() # function that performs moment interaction analysis def micurve(progress=None): # loop through all analysis points for idx, analysis_point in enumerate(analysis_list): # calculate ultimate results: # if we have axial forces if n_spacing: # limits should be calculated based on neutral axis values if idx == 0: ult_res = self.calculate_ultimate_section_actions( d_n=limits_dn[0], ultimate_results=res.UltimateBendingResults(theta=theta), ) elif idx == len(analysis_list) - 1: ult_res = self.calculate_ultimate_section_actions( d_n=limits_dn[1], ultimate_results=res.UltimateBendingResults(theta=theta), ) else: ult_res = self.ultimate_bending_capacity( theta=theta, n=analysis_point ) # if we have neutral axes else: ult_res = self.calculate_ultimate_section_actions( d_n=analysis_point, ultimate_results=res.UltimateBendingResults(theta=theta), ) # add labels for limits if labels: if idx == 0: ult_res.label = labels[0] elif idx == len(analysis_list) - 1: ult_res.label = labels[1] # add ultimate result to moment interactions results mi_results.results.append(ult_res) # update progress if progress: progress.update(task, advance=1) # add control points for idx, d_n in enumerate(add_cp_dn): ult_res = self.calculate_ultimate_section_actions( d_n=d_n, ultimate_results=res.UltimateBendingResults(theta=theta), ) # add label if labels: ult_res.label = labels[idx + 2] # add ultimate result to moment interactions results mi_results.results.append(ult_res) # update progress if progress: progress.update(task, advance=1) # sort results mi_results.sort_results() if progress_bar: # 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=len(analysis_list) + len(control_points), ) micurve(progress=progress) # display finished progress bar progress.update( task, description="[bold green]:white_check_mark: M-N diagram generated", ) live.refresh() else: micurve() # cut diagram at max_comp if max_comp: # check input - if value greater than maximum compression if max_comp > mi_results.results[0].n: msg = f"max_comp={max_comp} is greater than the maximum axial capacity " msg += f"{mi_results.results[0].n}." raise ValueError(msg) # get max_comp point ult_res = self.ultimate_bending_capacity(theta=theta, n=max_comp) # determine which results to delete idx_to_keep = 0 for idx, mi_res in enumerate(mi_results.results): # determine which index is the first to keep if idx_to_keep == 0 and mi_res.n < max_comp: idx_to_keep = idx break # remove points in diagram del mi_results.results[:idx_to_keep] # get labels if max_comp_labels: pt1_label = max_comp_labels[0] pt2_label = max_comp_labels[1] else: pt1_label = None pt2_label = None # add first two points to diagram # (m_max_comp, max_comp) # apply label ult_res.label = pt2_label mi_results.results.insert(0, ult_res) # (0, max_comp) mi_results.results.insert( 0, # insertion index res.UltimateBendingResults( theta=theta, d_n=inf, k_u=0, n=max_comp, m_x=0, m_y=0, m_xy=0, label=pt1_label, ), ) return mi_results
[docs] def biaxial_bending_diagram( self, n: float = 0, n_points: int = 48, progress_bar: bool = True, ) -> 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 :param progress_bar: If set to True, displays the progress bar :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) # function that performs biaxial bending analysis def bbcurve(progress=None): # loop through thetas for theta in theta_list: ultimate_results = self.ultimate_bending_capacity(theta=theta, n=n) bb_results.results.append(ultimate_results) if progress: progress.update(task, advance=1) # add first result to end of list top bb_results.results.append(bb_results.results[0]) if progress_bar: # 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, ) bbcurve(progress=progress) progress.update( task, description="[bold green]:white_check_mark: Biaxial bending diagram generated", ) live.refresh() else: bbcurve() 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, n_target=moment_curvature_results.n_target ) # find neutral axis that gives convergence of the axial force try: eps0, r = brentq( f=self.service_normal_force_convergence, a=-0.1, b=0.1, args=(kappa, mk), full_output=True, disp=False, ) except ValueError: msg = "Analysis failed. Confirm that the supplied moment/curvature is " msg += "within the range of the moment-curvature analysis." raise utils.AnalysisError(msg) # 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 global coordinates of extreme compressive fibre ecf, _ = utils.calculate_extreme_fibre( points=self.compound_geometry.points, 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_service( geom=meshed_geom, theta=theta, ecf=ecf, eps0=eps0, 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( kappa=kappa, ecf=ecf, eps0=eps0, theta=theta, centroid=self.moment_centroid, ) # 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]), ecf=ecf, eps0=eps0, 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.moment_centroid[0], centroid[1] - self.moment_centroid[1], ) ) 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_ultimate( geom=meshed_geom, theta=ultimate_results.theta, point_na=point_na, 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.moment_centroid, ) # 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.moment_centroid[0], centroid[1] - self.moment_centroid[1], ) ) 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 decode_d_n( self, theta: float, cp: Tuple[str, float], d_t: float, ) -> float: r"""Decodes a neutral axis depth given a control point ``cp``. :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) :param cp: Control point to decode :param d_t: Depth to extreme tensile fibre :return: Decoded neutral axis depth """ # 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 inf # control point type not valid else: msg = "First value of control point tuple must be D, d_n, fy, N or " msg += "kappa0." raise ValueError(msg)
[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 and strands for lumped_geom in self.reinf_geometries_lumped + self.strand_geometries: 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