Source code for concreteproperties.design_codes.as3600

from __future__ import annotations

from copy import deepcopy
from math import inf
from typing import TYPE_CHECKING, List, Optional, Tuple, Union

import numpy as np
from rich.live import Live
from scipy.interpolate import interp1d
from scipy.optimize import brentq

import concreteproperties.results as res
import concreteproperties.stress_strain_profile as ssp
import concreteproperties.utils as utils
from concreteproperties.design_codes.design_code import DesignCode
from concreteproperties.material import Concrete, SteelBar

if TYPE_CHECKING:
    from concreteproperties.concrete_section import ConcreteSection


[docs]class AS3600(DesignCode): """Design code class for Australian standard AS 3600:2018. .. note:: Note that this design code only supports :class:`~concreteproperties.material.Concrete` and :class:`~concreteproperties.material.SteelBar` material objects. Meshed :class:`~concreteproperties.material.Steel` material objects are **not** supported as this falls under the composite structures design code. """ def __init__(self): """Inits the AS3600 class.""" super().__init__()
[docs] def assign_concrete_section( self, concrete_section: ConcreteSection, ): """Assigns a concrete section to the design code. :param concrete_section: Concrete section object to analyse """ self.concrete_section = concrete_section # check to make sure there are no meshed reinforcement regions if self.concrete_section.reinf_geometries_meshed: raise ValueError( "Meshed reinforcement is not supported in this design code." ) # determine reinforcement class self.reinforcement_class = "N" for steel_geom in self.concrete_section.reinf_geometries_lumped: if ( abs( steel_geom.material.stress_strain_profile.get_ultimate_tensile_strain() ) < 0.05 ): self.reinforcement_class = "L" # calculate squash and tensile load squash, tensile = self.squash_tensile_load() self.squash_load = squash self.tensile_load = tensile
[docs] def create_concrete_material( self, compressive_strength: float, colour: str = "lightgrey", ) -> Concrete: r"""Returns a concrete material object to AS 3600:2018. .. admonition:: Material assumptions - *Density*: 2400 kg/m\ :sup:`3` - *Elastic modulus*: Interpolated from Table 3.1.2 - *Service stress-strain profile*: Linear with no tension, compressive strength at :math:`0.9f'_c` - *Ultimate stress-strain profile*: Rectangular stress block, parameters from Cl. 8.1.3 - *Alpha squash*: From Cl. 10.6.2.2 - *Flexural tensile strength*: From Cl. 3.1.1.3 :param compressive_strength: Characteristic compressive strength of concrete at 28 days in megapascals (MPa) :param colour: Colour of the concrete for rendering :raises ValueError: If ``compressive_strength`` is not between 20 MPa and 100 MPa. :return: Concrete material object """ if compressive_strength < 20 or compressive_strength > 100: raise ValueError("compressive_strength must be between 20 MPa and 100 MPa.") # create concrete name name = f"{compressive_strength:.0f} MPa Concrete (AS 3600:2018)" # calculate elastic modulus fc_list = [20, 25, 32, 40, 50, 65, 80, 100] Ec_list = [24000, 26700, 30100, 32800, 34800, 37400, 39600, 42200] f_Ec = interp1d(fc_list, Ec_list) elastic_modulus = f_Ec(compressive_strength) # calculate stress block parameters alpha = 0.85 - 0.0015 * compressive_strength alpha = max(alpha, 0.67) gamma = 0.97 - 0.0025 * compressive_strength gamma = max(gamma, 0.67) # calculate flexural_tensile_strength flexural_tensile_strength = 0.6 * np.sqrt(compressive_strength) return Concrete( name=name, density=2.4e-6, stress_strain_profile=ssp.ConcreteLinearNoTension( elastic_modulus=elastic_modulus, ultimate_strain=0.003, compressive_strength=0.9 * compressive_strength, ), ultimate_stress_strain_profile=ssp.RectangularStressBlock( compressive_strength=compressive_strength, alpha=alpha, gamma=gamma, ultimate_strain=0.003, ), flexural_tensile_strength=flexural_tensile_strength, colour=colour, )
[docs] def create_steel_material( self, yield_strength: float = 500, ductility_class: str = "N", colour: str = "grey", ) -> SteelBar: r"""Returns a steel bar material object. .. admonition:: Material assumptions - *Density*: 7850 kg/m\ :sup:`3` - *Elastic modulus*: 200000 MPa - *Stress-strain profile*: Elastic-plastic, fracture strain from Table 3.2.1 :param yield_strength: Steel yield strength :param ductility_class: Steel ductility class ("N" or "L") :param colour: Colour of the steel for rendering :raises ValueError: If ``ductility_class`` is not "N" or "L" :return: Steel material object """ if ductility_class == "N": fracture_strain = 0.05 elif ductility_class == "L": fracture_strain = 0.015 else: raise ValueError("ductility_class must be N or L.") return SteelBar( name=f"{yield_strength:.0f} MPa Steel (AS 3600:2018)", density=7.85e-6, stress_strain_profile=ssp.SteelElasticPlastic( yield_strength=yield_strength, elastic_modulus=200e3, fracture_strain=fracture_strain, ), colour=colour, )
[docs] def squash_tensile_load(self) -> Tuple[float, float]: """Calculates the squash and tensile load of the reinforced concrete section. :return: Squash and tensile load """ # initialise the squash load, tensile load and squash moment variables squash_load = 0 tensile_load = 0 # loop through all concrete geometries for conc_geom in self.concrete_section.concrete_geometries: # calculate area area = conc_geom.calculate_area() # calculate alpha_squash comp_strength = ( conc_geom.material.stress_strain_profile.get_compressive_strength() ) if comp_strength: alpha_squash = 1 - 0.003 * comp_strength alpha_squash = min(alpha_squash, 0.85) alpha_squash = max(alpha_squash, 0.72) else: alpha_squash = 1 # calculate compressive force force_c = ( area * alpha_squash * conc_geom.material.ultimate_stress_strain_profile.get_compressive_strength() ) # add to totals squash_load += force_c # loop through all steel geometries for steel_geom in self.concrete_section.reinf_geometries_lumped: # calculate area area = steel_geom.calculate_area() # calculate compressive and tensile force force_c = area * steel_geom.material.stress_strain_profile.get_stress( strain=0.025 ) force_t = ( -area * steel_geom.material.stress_strain_profile.get_yield_strength() ) # add to totals squash_load += force_c tensile_load += force_t return squash_load, tensile_load
[docs] def capacity_reduction_factor( self, n_u: float, n_ub: float, n_uot: float, k_uo: float, phi_0: float, ) -> float: """Returns the AS 3600:2018 capacity reduction factor (Table 2.2.2). ``n_ub`` and ``phi_0`` only required for compression, ``n_uot`` only required for tension. :param n_u: Axial force in member :param n_ub: Axial force at balanced point :param n_uot: Axial force at ultimate tension load :param k_uo: Neutral axis parameter at pure bending :param phi_0: Capacity reduction factor for dominant compression :return: Capacity reduction factor """ # pure bending phi if self.reinforcement_class == "N": phi = 1.24 - 13 * k_uo / 12 phi = min(phi, 0.85) phi = max(phi, 0.65) else: phi = 0.65 # compression if n_u > 0: if n_u >= n_ub: return phi_0 else: return phi_0 + (phi - phi_0) * (1 - n_u / n_ub) # tension else: if self.reinforcement_class == "N": return phi + (0.85 - phi) * (n_u / n_uot) else: return 0.65
[docs] def get_k_uo( self, theta: float, ) -> float: r"""Returns k_uo for the reinforced concrete cross-section given ``theta``. :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) :return: Bending parameter k_uo """ pure_res = self.concrete_section.ultimate_bending_capacity(theta=theta) return pure_res.k_u
[docs] def get_n_ub( self, theta: float, ) -> float: r"""Returns n_ub for the reinforced concrete cross-section given ``theta``. :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) :return: Balanced axial force n_ub """ # get depth to extreme tensile bar and its yield strain d_0, eps_sy = self.concrete_section.extreme_bar(theta=theta) # get compressive strain at extreme fibre eps_cu = self.concrete_section.gross_properties.conc_ultimate_strain # calculate d_n at balanced load d_nb = d_0 * (eps_cu) / (eps_sy + eps_cu) # calculate axial force at balanced load balanced_res = self.concrete_section.calculate_ultimate_section_actions( d_n=d_nb, ultimate_results=res.UltimateBendingResults(theta=theta) ) return balanced_res.n
[docs] def ultimate_bending_capacity( self, theta: float = 0, n_design: float = 0, phi_0: float = 0.6, ) -> Tuple[res.UltimateBendingResults, res.UltimateBendingResults, float]: r"""Calculates the ultimate bending capacity with capacity factors to AS 3600:2018. :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) :param n_design: Design axial force, N* :param phi_0: Compression dominant capacity reduction factor, see Table 2.2.2(d) :return: Factored and unfactored ultimate bending results objects, and capacity reduction factor *(factored_results, unfactored_results, phi)* """ # get parameters to determine phi n_uot = self.tensile_load k_uo = self.get_k_uo(theta=theta) n_ub = self.get_n_ub(theta=theta) # non-linear calculation of phi def non_linear_phi(phi_guess): phi = self.capacity_reduction_factor( n_u=n_design / phi_guess, n_ub=n_ub, n_uot=n_uot, k_uo=k_uo, phi_0=phi_0, ) return phi - phi_guess phi, _ = brentq( f=non_linear_phi, a=phi_0, b=0.85, xtol=1e-3, rtol=1e-6, # type: ignore full_output=True, disp=False, ) # generate basic moment interaction diagram f_mi_res, _, _ = self.moment_interaction_diagram( theta=theta, control_points=[ ("N", 0.0), ], n_points=2, phi_0=phi_0, progress_bar=False, ) # get significant axial loads n_squash = f_mi_res.results[0].n n_decomp = f_mi_res.results[1].n n_tensile = f_mi_res.results[-1].n # DETERMINE where we are on interaction diagram # if we are above the squash load or tensile load if n_design > n_squash: raise utils.AnalysisError( f"N = {n_design} is greater than the squash load, phiNc = {n_squash}." ) elif n_design < n_tensile: raise utils.AnalysisError( f"N = {n_design} is greater than the tensile load, phiNt = {n_tensile}" ) # compression linear interpolation elif n_design > n_decomp: factor = (n_design - n_decomp) / (n_squash - n_decomp) squash = f_mi_res.results[0] decomp = f_mi_res.results[1] ult_res = res.UltimateBendingResults( theta=theta, d_n=inf, k_u=0, n=n_design / phi, m_x=(decomp.m_x + factor * (squash.m_x - decomp.m_x)) / phi, m_y=(decomp.m_y + factor * (squash.m_y - decomp.m_y)) / phi, m_xy=(decomp.m_xy + factor * (squash.m_xy - decomp.m_xy)) / phi, ) # regular calculation elif n_design > 0: ult_res = self.concrete_section.ultimate_bending_capacity( theta=theta, n=n_design / phi ) # tensile linear interpolation else: factor = n_design / n_tensile pure = f_mi_res.results[-2] ult_res = res.UltimateBendingResults( theta=theta, d_n=inf, k_u=0, n=n_design / phi, m_x=(1 - factor) * pure.m_x / phi, m_y=(1 - factor) * pure.m_y / phi, m_xy=(1 - factor) * pure.m_xy / phi, ) # factor ultimate results f_ult_res = deepcopy(ult_res) f_ult_res.n *= phi f_ult_res.m_x *= phi f_ult_res.m_y *= phi f_ult_res.m_xy *= phi return f_ult_res, ult_res, phi
[docs] def moment_interaction_diagram( self, theta: float = 0, limits: List[Tuple[str, float]] = [ ("D", 1.0), ("N", 0.0), ], control_points: List[Tuple[str, float]] = [ ("fy", 1.0), ], labels: Optional[List[str]] = None, n_points: int = 24, n_spacing: Optional[int] = None, phi_0: float = 0.6, progress_bar: bool = True, ) -> Tuple[res.MomentInteractionResults, res.MomentInteractionResults, List[float]]: r"""Generates a moment interaction diagram with capacity factors to AS 3600:2018. See :meth:`concreteproperties.concrete_section.ConcreteSection.moment_interaction_diagram` for allowable control points. .. note:: When providing ``"N"`` to ``limits`` or ``control_points``, ``"N"`` is taken to be the unfactored net (nominal) axial load :math:`N^{*} / \phi`. :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 the pure bending point. :param control_points: List of additional control points to add to the moment interaction diagram. The default control points include the balanced point (``fy=1``). 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 axes 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 phi_0: Compression dominant capacity reduction factor, see Table 2.2.2(d) :param progress_bar: If set to True, displays the progress bar :return: Factored and unfactored moment interaction results objects, and list of capacity reduction factors *(factored_results, unfactored_results, phis)* """ mi_res = self.concrete_section.moment_interaction_diagram( theta=theta, limits=limits, control_points=control_points, labels=labels, n_points=n_points, n_spacing=n_spacing, progress_bar=progress_bar, ) # get theta theta = mi_res.results[0].theta # add squash load mi_res.results.insert( 0, res.UltimateBendingResults( theta=theta, d_n=inf, k_u=0, n=self.squash_load, m_x=0, m_y=0, m_xy=0, ), ) # add tensile load mi_res.results.append( res.UltimateBendingResults( theta=theta, d_n=0, k_u=0, n=self.tensile_load, m_x=0, m_y=0, m_xy=0, ) ) # make a copy of the results to factor f_mi_res = deepcopy(mi_res) # list to store phis phis = [] # get required constants for phi n_uot = self.tensile_load k_uo = self.get_k_uo(theta=theta) n_ub = self.get_n_ub(theta=theta) # factor results for ult_res in f_mi_res.results: phi = self.capacity_reduction_factor( n_u=ult_res.n, n_ub=n_ub, n_uot=n_uot, k_uo=k_uo, phi_0=phi_0 ) ult_res.n *= phi ult_res.m_x *= phi ult_res.m_y *= phi ult_res.m_xy *= phi phis.append(phi) return f_mi_res, mi_res, phis
[docs] def biaxial_bending_diagram( self, n_design: float = 0, n_points: int = 48, phi_0: float = 0.6, progress_bar: bool = True, ) -> Tuple[res.BiaxialBendingResults, List[float]]: """Generates a biaxial bending with capacity factors to AS 3600:2018. :param n_design: Design axial force, N* :param n_points: Number of calculation points :param phi_0: Compression dominant capacity reduction factor, see Table 2.2.2(d) :param progress_bar: If set to True, displays the progress bar :return: Factored biaxial bending results object and list of capacity reduction factors *(factored_results, phis)* """ # initialise results f_bb_res = res.BiaxialBendingResults(n=n_design) phis = [] # 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: # factored capacity f_ult_res, _, phi = self.ultimate_bending_capacity( theta=theta, n_design=n_design, phi_0=phi_0 ) f_bb_res.results.append(f_ult_res) phis.append(phi) if progress: progress.update(task, advance=1) 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() # add first result to end of list top f_bb_res.results.append(f_bb_res.results[0]) phis.append(phis[0]) return f_bb_res, phis