From 2f0d3f9ad41ca5bd9c7fabe1ce6d9399fa48b2cd Mon Sep 17 00:00:00 2001 From: robbievanleeuwen Date: Tue, 2 Aug 2022 17:23:24 +1000 Subject: [PATCH 1/7] Start big restructure --- concreteproperties/analysis_section.py | 178 ++++---- concreteproperties/concrete_section.py | 390 +++++++++--------- concreteproperties/design_codes.py | 2 +- concreteproperties/material.py | 94 +++-- concreteproperties/pre.py | 344 ++++++++++++++- concreteproperties/results.py | 35 +- concreteproperties/stress_strain_profile.py | 36 +- concreteproperties/utils.py | 127 +++--- docs/source/notebooks/area_properties.ipynb | 4 +- .../source/notebooks/cracked_properties.ipynb | 4 +- docs/source/notebooks/moment_curvature.ipynb | 6 +- docs/source/notebooks/ultimate_bending.ipynb | 4 +- 12 files changed, 821 insertions(+), 403 deletions(-) diff --git a/concreteproperties/analysis_section.py b/concreteproperties/analysis_section.py index a0afaf8f..dd0bb20f 100644 --- a/concreteproperties/analysis_section.py +++ b/concreteproperties/analysis_section.py @@ -2,7 +2,7 @@ from dataclasses import dataclass from math import isinf -from typing import TYPE_CHECKING, List, Optional, Tuple +from typing import TYPE_CHECKING, List, Tuple import numpy as np import triangle @@ -10,27 +10,28 @@ import concreteproperties.utils as utils from concreteproperties.post import plotting_context +from concreteproperties.material import Concrete if TYPE_CHECKING: import matplotlib - from sectionproperties.pre.geometry import Geometry - from concreteproperties.material import Concrete + from concreteproperties.material import Material + from concreteproperties.pre import CPGeom class AnalysisSection: - """Class for an analysis section to perform a fast analysis on concrete sections.""" + """Class for an analysis section to perform a fast analysis on meshed sections.""" def __init__( self, - geometry: Geometry, + geometry: CPGeom, ): """Inits the AnalysisSection class. :param geometry: Geometry object """ - self.geometry = geometry + self.material = geometry.material # create simple mesh tri = {} # create tri dictionary @@ -40,6 +41,7 @@ def __init__( if geometry.holes: tri["holes"] = geometry.holes # set holes + # coarse mesh self.mesh = triangle.triangulate(tri, "p") # extract mesh data @@ -51,9 +53,9 @@ def __init__( self.mesh_elements = [] # build elements - self.elements = [] + self.elements: List[Tri3] = [] - for idx, node_ids in enumerate(self.mesh_elements): + for node_ids in self.mesh_elements: x1 = self.mesh_nodes[node_ids[0]][0] y1 = self.mesh_nodes[node_ids[0]][1] x2 = self.mesh_nodes[node_ids[1]][0] @@ -67,10 +69,9 @@ def __init__( # add tri elements to the mesh self.elements.append( Tri3( - el_id=idx, coords=coords, node_ids=node_ids, - conc_material=self.geometry.material, # type: ignore + material=self.material, ) ) @@ -111,22 +112,22 @@ def get_elastic_stress( y = node[1] - cy # axial stress - sig[idx] += n * self.geometry.material.elastic_modulus / e_a + sig[idx] += n * self.material.elastic_modulus / e_a # bending moment - sig[idx] += self.geometry.material.elastic_modulus * ( + sig[idx] += self.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[idx] += self.geometry.material.elastic_modulus * ( + sig[idx] += self.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 ) # initialise section actions - n_conc = 0 - m_x_conc = 0 - m_y_conc = 0 + n_sec = 0 + m_x_sec = 0 + m_y_sec = 0 for el in self.elements: el_n, el_m_x, el_m_y = el.calculate_elastic_actions( @@ -141,28 +142,27 @@ def get_elastic_stress( e_ixy=e_ixy, ) - n_conc += el_n - m_x_conc += el_m_x - m_y_conc += el_m_y + n_sec += el_n + m_x_sec += el_m_x + m_y_sec += el_m_y # calculate point of action - if n_conc == 0: + if n_sec == 0: d_x = 0 d_y = 0 else: - d_x = m_y_conc / n_conc - d_y = m_x_conc / n_conc + d_x = m_y_sec / n_sec + d_y = m_x_sec / n_sec - return sig, n_conc, d_x, d_y + return sig, n_sec, d_x, d_y - def service_stress_analysis( + def service_analysis( self, point_na: Tuple[float, float], - d_n: float, theta: float, kappa: float, centroid: Tuple[float, float], - ) -> Tuple[float, float, float, float]: + ) -> Tuple[float, float, float, float, float]: r"""Performs a service stress analysis on the section. :param point_na: Point on the neutral axis @@ -172,30 +172,37 @@ def service_stress_analysis( :param kappa: Curvature :param centroid: Centroid about which to take moments - :return: Axial force, section moments and max strain + :return: Axial force, section moments and min/max strain """ # initialise section actions - n = 0 - m_x = 0 - m_y = 0 + n_sec = 0 + m_x_sec = 0 + m_y_sec = 0 + min_strain = 0 max_strain = 0 for el in self.elements: - el_n, el_m_x, el_m_y, el_max_strain = el.calculate_service_actions( + ( + el_n, + el_m_x, + el_m_y, + el_min_strain, + el_max_strain, + ) = el.calculate_service_actions( point_na=point_na, - d_n=d_n, theta=theta, kappa=kappa, centroid=centroid, ) + min_strain = min(min_strain, el_min_strain) max_strain = max(max_strain, el_max_strain) - n += el_n - m_x += el_m_x - m_y += el_m_y + n_sec += el_n + m_x_sec += el_m_x + m_y_sec += el_m_y - return n, m_x, m_y, max_strain + return n_sec, m_x_sec, m_y_sec, min_strain, max_strain def get_service_stress( self, @@ -233,30 +240,27 @@ def get_service_stress( ) # get stress at gauss point - sig[idx] = self.geometry.material.stress_strain_profile.get_stress( # type: ignore - strain=strain - ) + sig[idx] = self.material.stress_strain_profile.get_stress(strain=strain) # calculate total force - n, m_x, m_y, _ = self.service_stress_analysis( + n_sec, m_x_sec, m_y_sec, _, _ = self.service_analysis( point_na=point_na, - d_n=d_n, theta=theta, kappa=kappa, centroid=centroid, ) # calculate point of action - if n == 0: + if n_sec == 0: d_x = 0 d_y = 0 else: - d_x = m_y / n - d_y = m_x / n + d_x = m_y_sec / n_sec + d_y = m_x_sec / n_sec - return sig, n, d_x, d_y + return sig, n_sec, d_x, d_y - def ultimate_stress_analysis( + def ultimate_analysis( self, point_na: Tuple[float, float], d_n: float, @@ -277,9 +281,9 @@ def ultimate_stress_analysis( """ # initialise section actions - n = 0 - m_x = 0 - m_y = 0 + n_sec = 0 + m_x_sec = 0 + m_y_sec = 0 for el in self.elements: el_n, el_m_x, el_m_y = el.calculate_ultimate_actions( @@ -290,11 +294,11 @@ def ultimate_stress_analysis( centroid=centroid, ) - n += el_n - m_x += el_m_x - m_y += el_m_y + n_sec += el_n + m_x_sec += el_m_x + m_y_sec += el_m_y - return n, m_x, m_y + return n_sec, m_x_sec, m_y_sec def get_ultimate_stress( self, @@ -336,12 +340,15 @@ def get_ultimate_stress( ) # get stress at gauss point - sig[idx] = self.geometry.material.ultimate_stress_strain_profile.get_stress( # type: ignore - strain=strain - ) + if isinstance(self.material, Concrete): + sig[idx] = self.material.ultimate_stress_strain_profile.get_stress( + strain=strain + ) + else: + sig[idx] = self.material.stress_strain_profile.get_stress(strain=strain) # calculate total force - n, m_x, m_y = self.ultimate_stress_analysis( + n_sec, m_x_sec, m_y_sec = self.ultimate_analysis( point_na=point_na, d_n=d_n, theta=theta, @@ -350,14 +357,14 @@ def get_ultimate_stress( ) # calculate point of action - if n == 0: + if n_sec == 0: d_x = 0 d_y = 0 else: - d_x = m_y / n - d_y = m_x / n + d_x = m_y_sec / n_sec + d_y = m_x_sec / n_sec - return sig, n, d_x, d_y + return sig, n_sec, d_x, d_y def plot_mesh( self, @@ -379,8 +386,8 @@ def plot_mesh( c = [] # Indices of elements for mapping colours # create an array of finite element colours - for idx, element in enumerate(self.elements): - colour_array.append(element.conc_material.colour) + for idx, el in enumerate(self.elements): + colour_array.append(el.material.colour) c.append(idx) cmap = ListedColormap(colour_array) # custom colourmap @@ -421,8 +428,8 @@ def plot_shape( c = [] # Indices of elements for mapping colours # create an array of finite element colours - for idx, element in enumerate(self.elements): - colour_array.append(element.conc_material.colour) + for idx, el in enumerate(self.elements): + colour_array.append(el.material.colour) c.append(idx) cmap = ListedColormap(colour_array) # custom colourmap @@ -441,16 +448,14 @@ def plot_shape( class Tri3: """Class for a three noded linear triangular element. - :param el_id: Unique element id - :param coords: A 2 x 3 array of the coordinates of the tri-3 nodes. + :param coords: A 2 x 3 array of the coordinates of the tri-3 nodes :param node_ids: A list of the global node ids for the current element - :param conc_material: Material object for the current finite element. + :param material: Material object for the current finite element """ - el_id: int coords: np.ndarray node_ids: List[int] - conc_material: Concrete + material: Material def second_moments_of_area( self, @@ -474,19 +479,19 @@ def second_moments_of_area( N, j = utils.shape_function(coords=self.coords, gauss_point=gp) e_ixx += ( - self.conc_material.elastic_modulus + self.material.elastic_modulus * gp[0] * np.dot(N, np.transpose(self.coords[1, :])) ** 2 * j ) e_iyy += ( - self.conc_material.elastic_modulus + self.material.elastic_modulus * gp[0] * np.dot(N, np.transpose(self.coords[0, :])) ** 2 * j ) e_ixy += ( - self.conc_material.elastic_modulus + self.material.elastic_modulus * gp[0] * np.dot(N, np.transpose(self.coords[1, :])) * np.dot(N, np.transpose(self.coords[0, :])) @@ -541,12 +546,12 @@ def calculate_elastic_actions( # axial force force_gp = 0 - force_gp += gp[0] * n * self.conc_material.elastic_modulus / e_a * j + force_gp += gp[0] * n * self.material.elastic_modulus / e_a * j # bending moment force_gp += ( gp[0] - * self.conc_material.elastic_modulus + * self.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 @@ -555,7 +560,7 @@ def calculate_elastic_actions( ) force_gp += ( gp[0] - * self.conc_material.elastic_modulus + * self.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 @@ -573,27 +578,26 @@ def calculate_elastic_actions( def calculate_service_actions( self, point_na: Tuple[float, float], - d_n: float, theta: float, kappa: float, centroid: Tuple[float, float], - ) -> Tuple[float, float, float, float]: + ) -> Tuple[float, float, float, float, float]: r"""Calculates service actions for the current finite element. :param point_na: Point on the neutral axis - :param d_n: Depth of the neutral axis from the extreme compression fibre :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) :param kappa: Curvature :param centroid: Centroid about which to take moments - :return: Axial force, moments and maximum strain + :return: Axial force, moments and min/max strain """ # initialise element results force_e = 0 m_x_e = 0 m_y_e = 0 + min_strain_e = 0 max_strain_e = 0 # get points for 1 point Gaussian integration @@ -615,10 +619,11 @@ def calculate_service_actions( theta=theta, kappa=kappa, ) + min_strain_e = min(min_strain_e, strain) max_strain_e = max(max_strain_e, strain) # get stress at gauss point - stress = self.conc_material.stress_strain_profile.get_stress(strain=strain) + stress = self.material.stress_strain_profile.get_stress(strain=strain) # calculate force (stress * area) force_gp = gp[0] * stress * j @@ -628,7 +633,7 @@ def calculate_service_actions( m_x_e += force_e * (y - centroid[1]) m_y_e += force_e * (x - centroid[0]) - return force_e, m_x_e, m_y_e, max_strain_e + return force_e, m_x_e, m_y_e, min_strain_e, max_strain_e def calculate_ultimate_actions( self, @@ -680,9 +685,12 @@ def calculate_ultimate_actions( ) # get stress at gauss point - stress = self.conc_material.ultimate_stress_strain_profile.get_stress( - strain=strain - ) + if isinstance(self.material, Concrete): + stress = self.material.ultimate_stress_strain_profile.get_stress( + strain=strain + ) + else: + stress = self.material.stress_strain_profile.get_stress(strain=strain) # calculate force (stress * area) force_gp = gp[0] * stress * j diff --git a/concreteproperties/concrete_section.py b/concreteproperties/concrete_section.py index f1946948..2c801070 100644 --- a/concreteproperties/concrete_section.py +++ b/concreteproperties/concrete_section.py @@ -5,7 +5,6 @@ from typing import TYPE_CHECKING, List, Optional, Tuple, Union import matplotlib.patches as mpatches -import matplotlib.pyplot as plt import numpy as np import sectionproperties.pre.geometry as sp_geom from rich.live import Live @@ -17,6 +16,7 @@ from concreteproperties.analysis_section import AnalysisSection from concreteproperties.material import Concrete, Steel from concreteproperties.post import plotting_context +from concreteproperties.pre import CPGeom, CPGeomConcrete if TYPE_CHECKING: import matplotlib @@ -31,38 +31,46 @@ def __init__( ): """Inits the ConcreteSection class. - :param geometry: *sectionproperties* compound geometry object describing the + :param geometry: *sectionproperties* CompoundGeometry object describing the reinforced concrete section """ - self.geometry = geometry - - # sort into concrete and steel geometries - self.concrete_geometries = [] - self.steel_geometries = [] - - for geom in self.geometry.geoms: - if isinstance(geom.material, Concrete): - self.concrete_geometries.append(geom) - if isinstance(geom.material, Steel): - self.steel_geometries.append(geom) - - # validate reinforced concrete input - if len(self.concrete_geometries) == 0 or len(self.steel_geometries) == 0: - raise ValueError( - "geometry must contain both Concrete and Steel geometries." - ) + self.compound_geometry = geometry # check overlapping regions - polygons = [sec_geom.geom for sec_geom in self.geometry.geoms] + 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.ConcreteProperties() + self.gross_properties = res.GrossProperties() # calculate gross area properties self.calculate_gross_area_properties() @@ -72,45 +80,36 @@ def calculate_gross_area_properties( ): """Calculates and stores gross section area properties.""" - # concrete areas - for conc_geom in self.concrete_geometries: + # loop through all geometries + for geom in self.all_geometries: # area and centroid of geometry - area = conc_geom.calculate_area() - centroid = conc_geom.calculate_centroid() + area = geom.calculate_area() + centroid = geom.calculate_centroid() - self.gross_properties.concrete_area += area - self.gross_properties.e_a += area * conc_geom.material.elastic_modulus - self.gross_properties.mass += area * conc_geom.material.density + 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 * conc_geom.material.elastic_modulus * centroid[1] + area * geom.material.elastic_modulus * centroid[1] ) self.gross_properties.e_qy += ( - area * conc_geom.material.elastic_modulus * centroid[0] + area * geom.material.elastic_modulus * centroid[0] ) - # steel area - for steel_geom in self.steel_geometries: - # area and centroid of geometry - area = steel_geom.calculate_area() - centroid = steel_geom.calculate_centroid() + # sum concrete areas + for conc_geom in self.concrete_geometries: + self.gross_properties.concrete_area += conc_geom.calculate_area() - self.gross_properties.steel_area += area - self.gross_properties.e_a += area * steel_geom.material.elastic_modulus - self.gross_properties.mass += area * steel_geom.material.density - self.gross_properties.e_qx += ( - area * steel_geom.material.elastic_modulus * centroid[1] - ) - self.gross_properties.e_qy += ( - area * steel_geom.material.elastic_modulus * centroid[0] - ) + # sum reinforcement meshed areas + for meshed_geom in self.reinf_geometries_meshed: + self.gross_properties.reinf_meshed_area += meshed_geom.calculate_area() - # total area - self.gross_properties.total_area = ( - self.gross_properties.concrete_area + self.gross_properties.steel_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.geometry.calculate_perimeter() + self.gross_properties.perimeter = self.compound_geometry.calculate_perimeter() # centroids self.gross_properties.cx = ( @@ -121,30 +120,30 @@ def calculate_gross_area_properties( ) # global second moments of area - # concrete geometries - for conc_geom in self.concrete_geometries: - conc_sec = AnalysisSection(geometry=conc_geom) + # meshed geometries + for geom in self.meshed_geometries: + sec = AnalysisSection(geometry=geom) - for conc_el in conc_sec.elements: - el_e_ixx_g, el_e_iyy_g, el_e_ixy_g = conc_el.second_moments_of_area() + 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 - # steel geometries - for steel_geom in self.steel_geometries: + # lumped geometries - treat as lumped circles + for geom in self.reinf_geometries_lumped: # area, diameter and centroid of geometry - area = steel_geom.calculate_area() + area = geom.calculate_area() diam = np.sqrt(4 * area / np.pi) - centroid = steel_geom.calculate_centroid() + centroid = geom.calculate_centroid() - self.gross_properties.e_ixx_g += steel_geom.material.elastic_modulus * ( + 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 += steel_geom.material.elastic_modulus * ( + 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 += steel_geom.material.elastic_modulus * ( + self.gross_properties.e_ixy_g += geom.material.elastic_modulus * ( area * centroid[0] * centroid[1] ) @@ -189,7 +188,7 @@ def calculate_gross_area_properties( ) # centroidal section moduli - x_min, x_max, y_min, y_max = self.geometry.calculate_extents() + 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 ) @@ -205,7 +204,7 @@ def calculate_gross_area_properties( # principal section moduli x11_max, x11_min, y22_max, y22_min = utils.calculate_local_extents( - geometry=self.geometry, + geometry=self.compound_geometry, cx=self.gross_properties.cx, cy=self.gross_properties.cy, theta=self.gross_properties.phi, @@ -222,7 +221,7 @@ def calculate_gross_area_properties( for idx, conc_geom in enumerate(self.concrete_geometries): ult_strain = ( - conc_geom.material.ultimate_stress_strain_profile.get_ultimate_strain() + conc_geom.material.ultimate_stress_strain_profile.get_ultimate_compressive_strain() ) if idx == 0: conc_ult_strain = ult_strain @@ -233,10 +232,10 @@ def calculate_gross_area_properties( def get_gross_properties( self, - ) -> res.ConcreteProperties: + ) -> res.GrossProperties: """Returns the gross section properties of the reinforced concrete section. - :return: Concrete properties object + :return: Gross concrete properties object """ return self.gross_properties @@ -260,7 +259,7 @@ def calculate_cracked_properties( self, theta: float = 0, ) -> res.CrackedResults: - r"""Calculates cracked section properties given a neutral axis angle `theta`. + 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`) @@ -273,7 +272,9 @@ def calculate_cracked_properties( # set neutral axis depth limits # depth of neutral axis at extreme tensile fibre - _, d_t = utils.calculate_extreme_fibre(points=self.geometry.points, theta=theta) + _, 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 @@ -312,21 +313,17 @@ def calculate_cracked_properties( # global second moments of area for geom in cracked_results.cracked_geometries: - # if concrete - if isinstance(geom.material, Concrete): - conc_sec = AnalysisSection(geometry=geom) - - for conc_el in conc_sec.elements: - ( - el_e_ixx_g, - el_e_iyy_g, - el_e_ixy_g, - ) = conc_el.second_moments_of_area() + # 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 - - elif isinstance(geom.material, Steel): + # if lumped + else: # area, diameter and centroid of geometry area = geom.calculate_area() diam = np.sqrt(4 * area / np.pi) @@ -391,7 +388,7 @@ def calculate_cracking_moment( self, theta: float, ) -> float: - r"""Calculates the cracking moment given a bending angle `theta`. + 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`) @@ -431,10 +428,10 @@ def calculate_cracking_moment( f_t = conc_geom.material.flexural_tensile_strength m_c_geom = (f_t / conc_geom.material.elastic_modulus) * (e_iuu / d) - # if first geometry, initialise cracking moment + # if we are the first geometry, initialise cracking moment if idx == 0: m_c = m_c_geom - # otherwise take smallest cracking moment + # otherwise take smaller cracking moment else: m_c = min(m_c, m_c_geom) @@ -445,7 +442,7 @@ def cracked_neutral_axis_convergence( d_nc: float, cracked_results: res.CrackedResults, ) -> float: - """Given a trial cracked neutral axis depth `d_nc`, determines the difference + """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 @@ -456,7 +453,7 @@ def cracked_neutral_axis_convergence( # calculate extreme fibre in global coordinates extreme_fibre, d_t = utils.calculate_extreme_fibre( - points=self.geometry.points, theta=cracked_results.theta + points=self.compound_geometry.points, theta=cracked_results.theta ) # validate d_nc input @@ -476,11 +473,10 @@ def cracked_neutral_axis_convergence( ) # split concrete geometries above and below d_nc, discard below - cracked_geoms = [] + cracked_geoms: List[Union[CPGeomConcrete, CPGeom]] = [] for conc_geom in self.concrete_geometries: - top_geoms, _ = utils.split_section( - geometry=conc_geom, + top_geoms, _ = conc_geom.split_section( point=point_na, theta=cracked_results.theta, ) @@ -488,13 +484,13 @@ def cracked_neutral_axis_convergence( # 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 - # add steel geometries to list - cracked_geoms.extend(self.steel_geometries) - - # concrete & steel for geom in cracked_geoms: ea = geom.calculate_area() * geom.material.elastic_modulus centroid = geom.calculate_centroid() @@ -518,10 +514,9 @@ def moment_curvature_analysis( 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`. + r"""Performs a moment curvature analysis given a bending angle ``theta``. - Analysis continues until the steel reaches fracture strain or the concrete - reaches its ultimate strain. + 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`) @@ -538,7 +533,9 @@ def moment_curvature_analysis( # set neutral axis depth limits # depth of neutral axis at extreme tensile fibre - _, d_t = utils.calculate_extreme_fibre(points=self.geometry.points, theta=theta) + _, 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 @@ -551,7 +548,7 @@ def moment_curvature_analysis( total=None, ) - # while there hasn't been a failure in the steel + # while there hasn't been a failure while not moment_curvature._failure: # calculate adaptive step size for curvature if iter > 1: @@ -580,14 +577,6 @@ def moment_curvature_analysis( ) except ValueError: warnings.warn("brentq algorithm failed.") - d_n = 0 - - # calculate moments - self.service_normal_force_convergence( - d_n=d_n, - kappa=kappa, - moment_curvature=moment_curvature, - ) m_xy = np.sqrt( moment_curvature._m_x_i**2 + moment_curvature._m_y_i**2 @@ -636,7 +625,7 @@ def service_normal_force_convergence( # calculate extreme fibre in global coordinates extreme_fibre, d_t = utils.calculate_extreme_fibre( - points=self.geometry.points, theta=moment_curvature.theta + points=self.compound_geometry.points, theta=moment_curvature.theta ) # validate d_n input @@ -650,26 +639,31 @@ def service_normal_force_convergence( extreme_fibre=extreme_fibre, d_n=d_n, theta=moment_curvature.theta ) - # create splits in concrete geometries at points in stress-strain profiles - concrete_split_geoms = utils.split_section_at_strains( - concrete_geometries=self.concrete_geometries, - theta=moment_curvature.theta, - point_na=point_na, - ultimate=False, - kappa=kappa, - ) + # 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 concrete actions - for conc_geom in concrete_split_geoms: - sec = AnalysisSection(geometry=conc_geom) - n_sec, m_x_sec, m_y_sec, max_strain = sec.service_stress_analysis( + # 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, - d_n=d_n, theta=moment_curvature.theta, kappa=kappa, centroid=(self.gross_properties.cx, self.gross_properties.cy), @@ -679,44 +673,58 @@ def service_normal_force_convergence( m_x += m_x_sec m_y += m_y_sec - # check for concrete failure - if ( - max_strain - > conc_geom.material.stress_strain_profile.get_ultimate_strain() # type: ignore + # 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 = conc_geom + moment_curvature.failure_geometry = meshed_geom - # calculate steel actions - for steel_geom in self.steel_geometries: + # calculate lumped geometry actions + for lumped_geom in self.reinf_geometries_lumped: # calculate area and centroid - area = steel_geom.calculate_area() - steel_centroid = steel_geom.calculate_centroid() + area = lumped_geom.calculate_area() + centroid = lumped_geom.calculate_centroid() - # get strain at centroid of steel + # get strain at centroid of lump strain = utils.get_service_strain( - point=(steel_centroid[0], steel_centroid[1]), + point=(centroid[0], centroid[1]), point_na=point_na, theta=moment_curvature.theta, kappa=kappa, ) - # check for steel failure - if ( - abs(strain) - > steel_geom.material.stress_strain_profile.get_ultimate_strain() - ): + # 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 = steel_geom + moment_curvature.failure_geometry = lumped_geom # calculate stress and force - stress = steel_geom.material.stress_strain_profile.get_stress(strain=strain) + stress = lumped_geom.material.stress_strain_profile.get_stress( + strain=strain + ) force = stress * area n += force # calculate moment - m_x += force * (steel_centroid[1] - self.gross_properties.cy) - m_y += force * (steel_centroid[0] - self.gross_properties.cx) + 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 @@ -730,8 +738,8 @@ def ultimate_bending_capacity( 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. + r"""Given a neutral axis angle ``theta`` and an axial force ``n``, calculates + the ultimate bending capacity. :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) @@ -742,7 +750,9 @@ def ultimate_bending_capacity( # set neutral axis depth limits # depth of neutral axis at extreme tensile fibre - _, d_t = utils.calculate_extreme_fibre(points=self.geometry.points, theta=theta) + _, 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 @@ -772,8 +782,9 @@ def ultimate_normal_force_convergence( 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 axial force. + """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 @@ -795,11 +806,11 @@ def calculate_ultimate_section_actions( 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`. + """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 kappa0: If set to true, overwrites d_n and sets zero curvature :param ultimate_results: Ultimate bending results object :return: Ultimate bending results object @@ -810,7 +821,7 @@ def calculate_ultimate_section_actions( # calculate extreme fibre in global coordinates extreme_fibre, _ = utils.calculate_extreme_fibre( - points=self.geometry.points, theta=ultimate_results.theta + points=self.compound_geometry.points, theta=ultimate_results.theta ) # extreme fibre in local coordinates @@ -832,18 +843,27 @@ def calculate_ultimate_section_actions( extreme_fibre=extreme_fibre, d_n=d_n, theta=ultimate_results.theta ) - # create splits in concrete geometries at points in stress-strain profiles + # create splits in meshed geometries at points in stress-strain profiles + meshed_split_geoms: List[Union[CPGeom, CPGeomConcrete]] = [] + if isinf(d_n): - concrete_split_geoms = self.concrete_geometries + meshed_split_geoms = self.meshed_geometries else: - concrete_split_geoms = utils.split_section_at_strains( - concrete_geometries=self.concrete_geometries, - theta=ultimate_results.theta, - point_na=point_na, - ultimate=True, - ultimate_strain=self.gross_properties.conc_ultimate_strain, - d_n=d_n, - ) + 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) + + # # TESTING + # for geom in meshed_split_geoms: + # geom.plot_geometry() # initialise results n = 0 @@ -851,10 +871,10 @@ def calculate_ultimate_section_actions( m_y = 0 k_u = [] - # calculate concrete actions - for conc_geom in concrete_split_geoms: - sec = AnalysisSection(geometry=conc_geom) - n_sec, m_x_sec, m_y_sec = sec.ultimate_stress_analysis( + # 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, @@ -866,13 +886,13 @@ def calculate_ultimate_section_actions( m_x += m_x_sec m_y += m_y_sec - # calculate steel actions - for steel_geom in self.steel_geometries: + # calculate lumped actions + for lumped_geom in self.reinf_geometries_lumped: # calculate area and centroid - area = steel_geom.calculate_area() - centroid = steel_geom.calculate_centroid() + area = lumped_geom.calculate_area() + centroid = lumped_geom.calculate_centroid() - # get strain at centroid of steel + # get strain at centroid of lump if isinf(d_n): strain = self.gross_properties.conc_ultimate_strain else: @@ -885,7 +905,7 @@ def calculate_ultimate_section_actions( ) # calculate stress and force - stress = steel_geom.material.stress_strain_profile.get_stress(strain=strain) + stress = lumped_geom.material.stress_strain_profile.get_stress(strain=strain) force = stress * area n += force @@ -1788,49 +1808,51 @@ def plot_section( plotted_materials = [] legend_labels = [] - # plot concrete geometries - for conc_geom in self.concrete_geometries: - if conc_geom.material not in plotted_materials: + # plot meshed geometries + for meshed_geom in self.meshed_geometries: + if meshed_geom.material not in plotted_materials: patch = mpatches.Patch( - color=conc_geom.material.colour, label=conc_geom.material.name + color=meshed_geom.material.colour, + label=meshed_geom.material.name, ) legend_labels.append(patch) - plotted_materials.append(conc_geom.material) + plotted_materials.append(meshed_geom.material) # TODO - when shapely implements polygon plotting, fix this up - sec = AnalysisSection(geometry=conc_geom) + sec = AnalysisSection(geometry=meshed_geom) if not background: sec.plot_shape(ax=ax) # plot the points and facets - for f in conc_geom.facets: + for f in meshed_geom.facets: if background: fmt = "k-" else: fmt = "ko-" ax.plot( # type: ignore - [conc_geom.points[f[0]][0], conc_geom.points[f[1]][0]], - [conc_geom.points[f[0]][1], conc_geom.points[f[1]][1]], + [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 steel geometries - for steel_geom in self.steel_geometries: - if steel_geom.material not in plotted_materials: + # plot lumped geometries + for lumped_geom in self.reinf_geometries_lumped: + if lumped_geom.material not in plotted_materials: patch = mpatches.Patch( - color=steel_geom.material.colour, label=steel_geom.material.name + color=lumped_geom.material.colour, + label=lumped_geom.material.name, ) legend_labels.append(patch) - plotted_materials.append(steel_geom.material) + plotted_materials.append(lumped_geom.material) # plot the points and facets - coords = list(steel_geom.geom.exterior.coords) + coords = list(lumped_geom.geom.exterior.coords) bar = mpatches.Polygon( - xy=coords, closed=False, color=steel_geom.material.colour + xy=coords, closed=False, color=lumped_geom.material.colour ) ax.add_patch(bar) # type: ignore diff --git a/concreteproperties/design_codes.py b/concreteproperties/design_codes.py index e3fb4597..cf384383 100644 --- a/concreteproperties/design_codes.py +++ b/concreteproperties/design_codes.py @@ -76,7 +76,7 @@ def create_steel_material( def get_gross_properties( self, **kwargs, - ) -> res.ConcreteProperties: + ) -> res.GrossProperties: """Returns the gross section properties of the reinforced concrete section. :param kwargs: Keyword arguments passed to diff --git a/concreteproperties/material.py b/concreteproperties/material.py index 0427e460..12a16cc2 100644 --- a/concreteproperties/material.py +++ b/concreteproperties/material.py @@ -3,15 +3,34 @@ from dataclasses import dataclass, field from typing import TYPE_CHECKING -from concreteproperties.stress_strain_profile import ( - ConcreteServiceProfile, - ConcreteUltimateProfile, - SteelProfile, -) +import concreteproperties.stress_strain_profile as ssp -@dataclass(eq=True) -class Concrete: +@dataclass +class Material: + """Abstract class for a *concreteproperties* material. + + :param name: Material name + :param density: Material density (mass per unit volume) + :param stress_strain_profile: Material stress-strain profile + :param colour: Colour of the material for rendering + :param meshed: If set to True, the material region is meshed; if set to False, the + material region is treated as a lumped circular mass at its centroid + """ + + name: str + density: float + stress_strain_profile: ssp.StressStrainProfile + colour: str + meshed: bool + + def __post_init__(self): + # set elastic modulus + self.elastic_modulus = self.stress_strain_profile.get_elastic_modulus() + + +@dataclass +class Concrete(Material): """Class for a concrete material. :param name: Concrete material name @@ -27,45 +46,64 @@ class Concrete: name: str density: float - stress_strain_profile: ConcreteServiceProfile - ultimate_stress_strain_profile: ConcreteUltimateProfile + stress_strain_profile: ssp.ConcreteServiceProfile + ultimate_stress_strain_profile: ssp.ConcreteUltimateProfile alpha_squash: float flexural_tensile_strength: float colour: str + meshed: bool = field(default=True, init=False) def __post_init__(self): - self.elastic_modulus = self.stress_strain_profile.get_elastic_modulus() + super().__post_init__() - if not isinstance(self.stress_strain_profile, ConcreteServiceProfile): - raise ValueError( - "Concrete stress_strain_profile must be a ConcreteServiceProfile object" - ) + if not isinstance(self.stress_strain_profile, ssp.ConcreteServiceProfile): + msg = "Concrete stress_strain_profile must be a " + msg += "ConcreteServiceProfile object." + raise ValueError(msg) - if not isinstance(self.ultimate_stress_strain_profile, ConcreteUltimateProfile): - raise ValueError( - "Concrete ultimate_stress_strain_profile must be a ConcreteUltimateProfile object" - ) + if not isinstance( + self.ultimate_stress_strain_profile, ssp.ConcreteUltimateProfile + ): + msg = "Concrete ultimate_stress_strain_profile must be a " + msg += "ConcreteUltimateProfile object." + raise ValueError(msg) -@dataclass(eq=True) -class Steel: +@dataclass +class Steel(Material): """Class for a steel material. :param name: Steel material name :param density: Steel density (mass per unit volume) - :param stress_strain_profile: Ultimate steel stress-strain profile + :param stress_strain_profile: Steel stress-strain profile :param colour: Colour of the material for rendering """ name: str density: float - stress_strain_profile: SteelProfile + stress_strain_profile: ssp.StressStrainProfile colour: str + meshed: bool = field(default=True, init=False) + + +@dataclass +class SteelBar(Steel): + """Class for a steel material. + + :param name: Steel bar material name + :param density: Steel bar density (mass per unit volume) + :param stress_strain_profile: Steel bar stress-strain profile + :param colour: Colour of the material for rendering + """ + + name: str + density: float + stress_strain_profile: ssp.StressStrainProfile + colour: str + meshed: bool = field(default=False, init=False) - def __post_init__(self): - self.elastic_modulus = self.stress_strain_profile.get_elastic_modulus() - if not isinstance(self.stress_strain_profile, SteelProfile): - raise ValueError( - "Steel stress_strain_profile must be a SteelProfile object" - ) +@dataclass +class SteelStrand(Steel): + # placeholder + pass diff --git a/concreteproperties/pre.py b/concreteproperties/pre.py index 70062772..e0caea47 100644 --- a/concreteproperties/pre.py +++ b/concreteproperties/pre.py @@ -1,14 +1,352 @@ from __future__ import annotations -from typing import TYPE_CHECKING, Optional, Tuple, Union +from typing import TYPE_CHECKING, List, Tuple, Union +from more_itertools import peekable import numpy as np +from shapely.geometry import Polygon, LineString +from shapely.ops import split import sectionproperties.pre.library.primitive_sections as sp_ps +from sectionproperties.pre.geometry import Geometry + +from concreteproperties.material import Concrete if TYPE_CHECKING: - from sectionproperties.pre.geometry import CompoundGeometry, Geometry + import matplotlib + + from sectionproperties.pre.geometry import CompoundGeometry + + from concreteproperties.material import Material, Steel + + +class CPGeom: + """Watered down implementation of the *sectionproperties* Geometry object, optimised + for *concreteproperties*. + """ + + def __init__( + self, + geom: Polygon, + material: Material, + ): + """Inits the CPGeom class. + + :param geom: Shapely polygon defining the geometry + :param material: Material to apply to the geometry + """ + + # round polygon points and save geometry + self.geom = self.round_geometry(geometry=geom, tol=6) + + # store material + self.material = material + + # create points and facets + self.points, self.facets = self.create_points_and_facets(geometry=self.geom) + + # create holes + self.holes: List[Tuple[float, float]] = [] + + for hole in self.geom.interiors: + hole_polygon = Polygon(hole) + self.holes += tuple(hole_polygon.representative_point().coords) + + def round_geometry( + self, + geometry: Polygon, + tol: int, + ) -> Polygon: + """Rounds the coordinates in ``geometry`` to tolerance ``tol``. + + :param geometry: Geometry to round + :param tol: Number of decimal places to round + + :return: Rounded geometry + """ + + if geometry.exterior: + rounded_exterior = np.round(geometry.exterior.coords, tol) + else: + rounded_exterior = np.array([None]) + + rounded_interiors = [] + for interior in geometry.interiors: + rounded_interiors.append(np.round(interior.coords, tol)) + + if not rounded_exterior.any(): + return Polygon() + else: + return Polygon(rounded_exterior, rounded_interiors) + + def create_points_and_facets( + self, + geometry: Polygon, + ) -> Tuple[List[Tuple[float, float]], List[Tuple[int, int]]]: + """Creates a list of points and facets from a shapely polygon. + + :param geometry: Shapely polygon from which to create points and facets + + :return: Points and facets + """ + + master_count = 0 + points: List[Tuple[float, float]] = [] + facets: List[Tuple[int, int]] = [] + + # perimeter, note in shapely last point == first point + if geometry.exterior: + for coords in list(geometry.exterior.coords[:-1]): + points.append(coords) + master_count += 1 + + facets += self.create_facets(points) + exterior_count = master_count + + # holes + for idx, hole in enumerate(geometry.interiors): + break_count = master_count + int_points = [] + + for coords in hole.coords[:-1]: + int_points.append(coords) + master_count += 1 + + # (idx > 0), (idx < 1) are like a 'step functions' + offset = break_count * (idx > 0) + exterior_count * (idx < 1) + facets += self.create_facets(int_points, offset=offset) + points += int_points + + return points, facets + + def create_facets( + self, points_list: List[Tuple[float, float]], offset: int = 0 + ) -> List[Tuple[int, int]]: + """Generates a list of facets given a list of points and a facet offset. + + :param points_list: List of ordered points to create facets from + :param offset: Facet offset integer + + :return: List of facets + """ + + idx_peeker = peekable([idx + offset for idx, _ in enumerate(points_list)]) + return [(item, idx_peeker.peek(offset)) for item in idx_peeker] + + def calculate_area( + self, + ) -> float: + """Calculates the area of the geometry. + + :return: Geometry area + """ + + return self.geom.area + + def calculate_centroid( + self, + ) -> Tuple[float, float]: + """Calculates the centroid of the geometry. + + :return: Geometry centroid + """ + + return self.geom.centroid.coords[0] + + def calculate_extents( + self, + ) -> Tuple[float, float, float, float]: + """Calculates the minimum and maximum ``x`` and ``y`` values among the points + describing the geometry. + + :return: Extents (``x_min``, ``x_max``, ``y_min``, ``y_max``) + """ + + min_x, min_y, max_x, max_y = self.geom.bounds # type: ignore + + return min_x, max_x, min_y, max_y + + def split_section( + self, + point: Tuple[float, float], + theta: float, + ) -> Tuple[List[CPGeom], List[CPGeom]]: + """Splits the geometry about a line. + + :param point: Point on line + :param theta: Angle line makes with horizontal axis + + :return: Geometries above and below the line + """ + + # generate unit vector + vector = np.cos(theta), np.sin(theta) + + # calculate bounds of geometry + bounds = self.calculate_extents() + + # generate line segment that matches bounds of geometry object + line_seg = self.create_line_segment(point=point, vector=vector, bounds=bounds) + + # check to see if line intersects geometry + if line_seg.intersects(self.geom): + # split geometries + polys = split(geom=self.geom, splitter=line_seg).geoms + else: + polys = [self.geom] + + # sort geometries + top_polys, bot_polys = self.sort_polys(polys=polys, point=point, vector=vector) # type: ignore + + # assign material properties and create cp geometry objects + top_geoms = [ + CPGeomConcrete(geom=poly, material=self.material) + if isinstance(self.material, Concrete) + else CPGeom(geom=poly, material=self.material) + for poly in top_polys + ] + bot_geoms = [ + CPGeomConcrete(geom=poly, material=self.material) + if isinstance(self.material, Concrete) + else CPGeom(geom=poly, material=self.material) + for poly in bot_polys + ] + + # ensure top geoms is in compression + if theta <= np.pi / 2 and theta >= -np.pi / 2: + return top_geoms, bot_geoms + else: + return bot_geoms, top_geoms + + def create_line_segment( + self, + point: Tuple[float, float], + vector: Tuple[float, float], + bounds: Tuple[float, float, float, float], + ) -> LineString: + """Creates a shapely line string defined by a ``point`` and ``vector`` and + bounded by ``bounds``. + + :param point: Point on line + :param vector: Vector defining direction of line + :param bounds: Bounds of the geometry + + :return: Shapely line string + """ + + tol = 1e-6 # distance to displace start of line from bounds + + # not a vertical line + if abs(vector[0]) > 1e-12: + v_ratio = vector[1] / vector[0] + x1 = bounds[0] - tol + x2 = bounds[1] + tol + y1 = v_ratio * (x1 - point[0]) + point[1] + y2 = v_ratio * (x2 - point[0]) + point[1] + + # vertical line + else: + v_ratio = vector[0] / vector[1] + y1 = bounds[2] - tol + y2 = bounds[3] + tol + x1 = v_ratio * (y1 - point[1]) + point[0] + x2 = v_ratio * (y2 - point[1]) + point[0] + + return LineString([(x1, y1), (x2, y2)]) + + def sort_polys( + self, + polys: List[Polygon], + point: Tuple[float, float], + vector: Tuple[float, float], + ) -> Tuple[List[Polygon], List[Polygon]]: + """Sorts polygons that are above and below the line. + + :param polys: Polygons to sort + :param point: Point on line + :param vector: Vector defining direction of line + + :return: Polygons above and below the line + """ + + top_polys: List[Polygon] = [] + bot_polys: List[Polygon] = [] + v_ratio = vector[1] / vector[0] + + for poly in polys: + # get point inside polygon + px, py = poly.representative_point().coords[0] + + # not a vertical line + if abs(vector[0]) > 1e-12: + # get point on line at x-coordinate of representative point + y_line = point[1] + (px - point[0]) * v_ratio + + # if we are below the line + if py < y_line: + bot_polys.append(poly) + # if we are above the line + else: + top_polys.append(poly) + + # vertical line + else: + # if we are to the right of the line + if px < point[0]: + bot_polys.append(poly) + # if we are to the left of the line + else: + top_polys.append(poly) + + return top_polys, bot_polys + + def plot_geometry( + self, + title: str = "Cross-Section Geometry", + **kwargs, + ) -> matplotlib.axes.Axes: # type: ignore + """Plots the geometry. + + :param title: Plot title + :param kwargs: Passed to + :meth:`~sectionproperties.pre.geometry.Geometry.plot_geometry` + + :return: Matplotlib axes object + """ + + return self.to_sp_geom().plot_geometry(title=title, **kwargs) + + def to_sp_geom( + self, + ) -> Geometry: + """Converts self to a *sectionproperties* geometry object. + + :return: *sectionproperties* geometry object + """ + + return Geometry(geom=self.geom, material=self.material) # type: ignore + + +class CPGeomConcrete(CPGeom): + """*concreteproperties* Geometry class for concrete geometries.""" + + def __init__( + self, + geom: Polygon, + material: Concrete, + ): + """Inits the CPGeomConcrete class. + + :param geom: Shapely polygon defining the geometry + :param material: Material to apply to the geometry + """ + + super().__init__( + geom=geom, + material=material, + ) - from concreteproperties.material import Steel + # ensure material is a Concrete object + self.material = material def add_bar( diff --git a/concreteproperties/results.py b/concreteproperties/results.py index eae5eee6..3eda860e 100644 --- a/concreteproperties/results.py +++ b/concreteproperties/results.py @@ -11,7 +11,6 @@ import numpy as np from matplotlib.collections import PatchCollection from matplotlib.colors import CenteredNorm # type: ignore -from mpl_toolkits import mplot3d from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable from rich.console import Console from rich.table import Table @@ -24,14 +23,14 @@ if TYPE_CHECKING: import matplotlib - from sectionproperties.pre.geometry import Geometry from concreteproperties.analysis_section import AnalysisSection from concreteproperties.concrete_section import ConcreteSection + from concreteproperties.pre import CPGeom @dataclass -class ConcreteProperties: +class GrossProperties: """Class for storing gross concrete section properties. All properties with an `e_` preceding the property are multiplied by the elastic @@ -43,7 +42,8 @@ class ConcreteProperties: # section areas total_area: float = 0 concrete_area: float = 0 - steel_area: float = 0 + reinf_meshed_area: float = 0 + reinf_lumped_area: float = 0 e_a: float = 0 # section mass @@ -101,7 +101,14 @@ def print_results( table.add_row("Total Area", "{:>{fmt}}".format(self.total_area, fmt=fmt)) table.add_row("Concrete Area", "{:>{fmt}}".format(self.concrete_area, fmt=fmt)) - table.add_row("Steel Area", "{:>{fmt}}".format(self.steel_area, fmt=fmt)) + table.add_row( + "Meshed Reinforcement Area", + "{:>{fmt}}".format(self.reinf_meshed_area, fmt=fmt), + ) + table.add_row( + "Lumped Reinforcement Area", + "{:>{fmt}}".format(self.reinf_lumped_area, fmt=fmt), + ) table.add_row("Axial Rigidity (EA)", "{:>{fmt}}".format(self.e_a, fmt=fmt)) table.add_row("Mass (per unit length)", "{:>{fmt}}".format(self.mass, fmt=fmt)) table.add_row("Perimeter", "{:>{fmt}}".format(self.perimeter, fmt=fmt)) @@ -143,7 +150,7 @@ class TransformedConcreteProperties: :param elastic_modulus: Reference elastic modulus """ - concrete_properties: ConcreteProperties = field(repr=False) + concrete_properties: GrossProperties = field(repr=False) elastic_modulus: float # area @@ -250,7 +257,7 @@ class CrackedResults: theta: float m_cr: float = 0 d_nc: float = 0 - cracked_geometries: List[Geometry] = field(default_factory=list, repr=False) + cracked_geometries: List[CPGeom] = field(default_factory=list, repr=False) e_a_cr: float = 0 e_qx_cr: float = 0 e_qy_cr: float = 0 @@ -311,8 +318,8 @@ def plot_cracked_geometries( title: str = "Cracked Geometries", **kwargs, ) -> matplotlib.axes.Axes: # type: ignore - """Plots the geometries that remain (are in compression or are steel) after a - cracked analysis. + """Plots the geometries that remain (are in compression or are reinforcement) + after a cracked analysis. :param title: Plot title :param kwargs: Passed to @@ -321,9 +328,9 @@ def plot_cracked_geometries( :return: Matplotlib axes object """ - return CompoundGeometry(self.cracked_geometries).plot_geometry( - title=title, **kwargs - ) + return CompoundGeometry( + [geom.to_sp_geom() for geom in self.cracked_geometries] + ).plot_geometry(title=title, **kwargs) def print_results( self, @@ -410,7 +417,7 @@ class MomentCurvatureResults: m_x: List[float] = field(default_factory=list) m_y: List[float] = field(default_factory=list) m_xy: List[float] = field(default_factory=list) - failure_geometry: Geometry = field(init=False) + failure_geometry: CPGeom = field(init=False) # for analysis _n_i: float = field(default=0, repr=False) @@ -1022,7 +1029,7 @@ class StressResult: concrete_analysis_sections: List[AnalysisSection] concrete_stresses: List[np.ndarray] concrete_forces: List[Tuple[float, float, float]] - steel_geometries: List[Geometry] + steel_geometries: List[CPGeom] steel_stresses: List[float] steel_strains: List[float] steel_forces: List[Tuple[float, float, float]] diff --git a/concreteproperties/stress_strain_profile.py b/concreteproperties/stress_strain_profile.py index b0eaf458..2357e376 100644 --- a/concreteproperties/stress_strain_profile.py +++ b/concreteproperties/stress_strain_profile.py @@ -125,16 +125,26 @@ def get_tensile_strength( return min(self.stresses) - def get_ultimate_strain( + def get_ultimate_compressive_strain( self, ) -> float: - """Returns the largest strain. + """Returns the largest compressive strain. :return: Ultimate strain """ return max(self.strains) + def get_ultimate_tensile_strain( + self, + ) -> float: + """Returns the largest tensile strain. + + :return: Ultimate strain + """ + + return min(self.strains) + def get_unique_strains( self, ) -> List[float]: @@ -173,7 +183,12 @@ def print_properties( "{:>{fmt}}".format(-self.get_tensile_strength(), fmt=fmt), ) table.add_row( - "Ultimate Strain", "{:>{fmt}}".format(self.get_ultimate_strain(), fmt=fmt) + "Ultimate Compressive Strain", + "{:>{fmt}}".format(self.get_ultimate_compressive_strain(), fmt=fmt), + ) + table.add_row( + "Ultimate Tensile Strain", + "{:>{fmt}}".format(self.get_ultimate_tensile_strain(), fmt=fmt), ) console = Console() @@ -238,7 +253,8 @@ def print_properties( "Elastic Modulus", "{:>{fmt}}".format(self.get_elastic_modulus(), fmt=fmt) ) table.add_row( - "Ultimate Strain", "{:>{fmt}}".format(self.get_ultimate_strain(), fmt=fmt) + "Ultimate Compressive Strain", + "{:>{fmt}}".format(self.get_ultimate_compressive_strain(), fmt=fmt), ) console = Console() @@ -277,7 +293,7 @@ def get_tensile_strength( return None - def get_ultimate_strain( + def get_ultimate_compressive_strain( self, ) -> float: """Returns the largest strain. @@ -454,7 +470,7 @@ def get_compressive_strength( return self.compressive_strength - def get_ultimate_strain( + def get_ultimate_compressive_strain( self, ) -> float: """Returns the ultimate strain, or largest compressive strain. @@ -465,7 +481,7 @@ def get_ultimate_strain( try: return self.ultimate_strain # type: ignore except AttributeError: - return super().get_ultimate_strain() + return super().get_ultimate_compressive_strain() def print_properties( self, @@ -485,7 +501,8 @@ def print_properties( "{:>{fmt}}".format(self.get_compressive_strength(), fmt=fmt), ) table.add_row( - "Ultimate Strain", "{:>{fmt}}".format(self.get_ultimate_strain(), fmt=fmt) + "Ultimate Compressive Strain", + "{:>{fmt}}".format(self.get_ultimate_compressive_strain(), fmt=fmt), ) console = Console() console.print(table) @@ -675,7 +692,8 @@ def print_properties( "{:>{fmt}}".format(-self.get_tensile_strength(), fmt=fmt), ) table.add_row( - "Fracture Strain", "{:>{fmt}}".format(self.get_ultimate_strain(), fmt=fmt) + "Fracture Strain", + "{:>{fmt}}".format(self.get_ultimate_tensile_strain(), fmt=fmt), ) console = Console() diff --git a/concreteproperties/utils.py b/concreteproperties/utils.py index 6e6cbee1..5b77220a 100644 --- a/concreteproperties/utils.py +++ b/concreteproperties/utils.py @@ -8,8 +8,10 @@ from rich.text import Text from sectionproperties.pre.geometry import CompoundGeometry +from concreteproperties.pre import CPGeomConcrete + if TYPE_CHECKING: - from sectionproperties.pre.geometry import Geometry + from concreteproperties.pre import CPGeom def get_service_strain( @@ -31,10 +33,10 @@ def get_service_strain( """ # convert point to local coordinates - u, v = global_to_local(theta=theta, x=point[0], y=point[1]) + _, v = global_to_local(theta=theta, x=point[0], y=point[1]) # convert point_na to local coordinates - u_na, v_na = global_to_local(theta=theta, x=point_na[0], y=point_na[1]) + _, v_na = global_to_local(theta=theta, x=point_na[0], y=point_na[1]) # calculate distance between NA and point in `v` direction d = v - v_na @@ -100,24 +102,27 @@ def point_on_neutral_axis( return local_to_global(theta=theta, u=u, v=v) -def split_section_at_strains( - concrete_geometries: List[Geometry], +def split_geom_at_strains( + geom: Union[CPGeom, CPGeomConcrete], theta: float, point_na: Tuple[float, float], ultimate: bool, ultimate_strain: Optional[float] = None, d_n: Optional[float] = None, kappa: Optional[float] = None, -) -> List[Geometry]: - r"""Splits concrete geometries at discontinuities in its stress-strain profile. +) -> Union[List[CPGeom], List[CPGeomConcrete]]: + r"""Splits geometries at discontinuities in its stress-strain profile. - :param concrete_geometries: List of concrete geometries + :param geom: Geometry to split :param theta: Angle (in radians) the neutral axis makes with the horizontal axis (:math:`-\pi \leq \theta \leq \pi`) :param point_na: Point on the neutral axis - :param ultimate: If set to True, uses ultimate stress-strain profile - :param ultimate_strain: Concrete strain at failure - :param d_n: Depth of the neutral axis from the extreme compression fibre + :param ultimate: If set to True, uses ultimate stress-strain profile for concrete + geometries + :param ultimate_strain: Concrete strain at failure (required for ``ultimate=True`` + only) + :param d_n: Depth of the neutral axis from the extreme compression fibre (required + for ``ultimate=True`` only) :param kappa: Curvature :return: List of split geometries @@ -125,82 +130,64 @@ def split_section_at_strains( # handle zero curvature if kappa == 0: - return concrete_geometries + return [geom] # create splits in concrete geometries at points in stress-strain profiles - concrete_split_geoms = [] + split_geoms: Union[List[CPGeom], List[CPGeomConcrete]] = [] - for conc_geom in concrete_geometries: - if ultimate: - strains = ( - conc_geom.material.ultimate_stress_strain_profile.get_unique_strains() # type: ignore - ) + if ultimate and isinstance(geom, CPGeomConcrete): + strains = geom.material.ultimate_stress_strain_profile.get_unique_strains() + else: + strains = geom.material.stress_strain_profile.get_unique_strains() + + # make geom a list of geometries + geom_list = [geom] + + # initialise top_geoms in case of two unique strains + top_geoms = geom_list + + # loop through intermediate points on stress-strain profile + for strain in strains[1:-1]: + # depth to points of *strain* from NA + # ultimate case + if ultimate and ultimate_strain and d_n: + d = strain / ultimate_strain * d_n + # service case + elif kappa: + d = strain / kappa else: - strains = conc_geom.material.stress_strain_profile.get_unique_strains() # type: ignore - - # initialise top_geoms in case of two unique strains - top_geoms = [conc_geom] + raise ValueError("Not enough arguments provided.") - # loop through intermediate points on stress-strain profile - for idx, strain in enumerate(strains[1:-1]): - # depth to point with `strain` from NA - if ultimate: - d = strain / ultimate_strain * d_n - else: - d = strain / kappa + # convert depth to global coordinates + dx, dy = local_to_global(theta=theta, u=0, v=d) - # convert depth to global coordinates - dx, dy = local_to_global(theta=theta, u=0, v=d) + # calculate location of point + pt = point_na[0] + dx, point_na[1] + dy - # calculate location of point with `strain` - pt = point_na[0] + dx, point_na[1] + dy + # make list of geometries that will need to continue to be split after the + # split operation, i.e. those above the split + continuing_geoms = [] - # split concrete geometry (from bottom up) - top_geoms, bot_geoms = split_section( - geometry=conc_geom, + # split concrete geometries (from bottom up) + for g in geom_list: + top_geoms, bot_geoms = g.split_section( point=pt, theta=theta, ) # save bottom geoms - concrete_split_geoms.extend(bot_geoms) - - # continue to split top geoms - conc_geom = CompoundGeometry(geoms=top_geoms) - - # save final top geoms - concrete_split_geoms.extend(top_geoms) + split_geoms.extend(bot_geoms) - return concrete_split_geoms + # save continuing geoms + continuing_geoms.extend(top_geoms) + # update geom_list for next strain + geom_list = continuing_geoms -def split_section( - geometry: Union[Geometry, CompoundGeometry], - point: Tuple[float, float], - theta: float, -) -> Tuple[List[Geometry], List[Geometry]]: - r"""Splits the geometry along a line defined by a `point` and rotation angle - `theta`. - - :param geometry: Geometry to split - :param point: Point at which to split the geometry `(x, y)` - :param theta: Angle (in radians) the neutral axis makes with the horizontal - axis (:math:`-\pi \leq \theta \leq \pi`) + # save final top geoms + split_geoms.extend(top_geoms) - :return: Split geometry above and below the line - """ - - # split the section using the sectionproperties method - top_geoms, bot_geoms = geometry.split_section( - point_i=np.round(point, 8), vector=(np.cos(theta), np.sin(theta)) - ) - - # ensure top geoms is in compression - # sectionproperties definition is based on global coordinate system only - if theta <= np.pi / 2 and theta >= -np.pi / 2: - return top_geoms, bot_geoms - else: - return bot_geoms, top_geoms + return split_geoms def calculate_extreme_fibre( diff --git a/docs/source/notebooks/area_properties.ipynb b/docs/source/notebooks/area_properties.ipynb index d101256d..2fe3bca2 100644 --- a/docs/source/notebooks/area_properties.ipynb +++ b/docs/source/notebooks/area_properties.ipynb @@ -23,7 +23,7 @@ "metadata": {}, "outputs": [], "source": [ - "from concreteproperties.material import Concrete, Steel\n", + "from concreteproperties.material import Concrete, SteelBar\n", "from concreteproperties.stress_strain_profile import (\n", " ConcreteLinear,\n", " RectangularStressBlock,\n", @@ -65,7 +65,7 @@ " colour=\"lightgrey\",\n", ")\n", "\n", - "steel = Steel(\n", + "steel = SteelBar(\n", " name=\"500 MPa Steel\",\n", " density=7.85e-6,\n", " stress_strain_profile=SteelElasticPlastic(\n", diff --git a/docs/source/notebooks/cracked_properties.ipynb b/docs/source/notebooks/cracked_properties.ipynb index 75956ed4..8ff2587f 100644 --- a/docs/source/notebooks/cracked_properties.ipynb +++ b/docs/source/notebooks/cracked_properties.ipynb @@ -24,7 +24,7 @@ "outputs": [], "source": [ "import numpy as np\n", - "from concreteproperties.material import Concrete, Steel\n", + "from concreteproperties.material import Concrete, SteelBar\n", "from concreteproperties.stress_strain_profile import (\n", " ConcreteLinear,\n", " RectangularStressBlock,\n", @@ -66,7 +66,7 @@ " colour=\"lightgrey\",\n", ")\n", "\n", - "steel = Steel(\n", + "steel = SteelBar(\n", " name=\"500 MPa Steel\",\n", " density=7.85e-6,\n", " stress_strain_profile=SteelElasticPlastic(\n", diff --git a/docs/source/notebooks/moment_curvature.ipynb b/docs/source/notebooks/moment_curvature.ipynb index b158f805..27176b49 100644 --- a/docs/source/notebooks/moment_curvature.ipynb +++ b/docs/source/notebooks/moment_curvature.ipynb @@ -24,7 +24,7 @@ "outputs": [], "source": [ "import numpy as np\n", - "from concreteproperties.material import Concrete, Steel\n", + "from concreteproperties.material import Concrete, SteelBar\n", "import concreteproperties.stress_strain_profile as ssp\n", "from sectionproperties.pre.library.primitive_sections import (\n", " rectangular_section,\n", @@ -112,7 +112,7 @@ " conc_nonlinear,\n", "]\n", "\n", - "steel_ep = Steel(\n", + "steel_ep = SteelBar(\n", " name=\"Steel - Elastic-Plastic\",\n", " density=7.85e-6,\n", " stress_strain_profile=ssp.SteelElasticPlastic(\n", @@ -123,7 +123,7 @@ " colour=\"grey\",\n", ")\n", "\n", - "steel_hd = Steel(\n", + "steel_hd = SteelBar(\n", " name=\"Steel - Hardening\",\n", " density=7.85e-6,\n", " stress_strain_profile=ssp.SteelHardening(\n", diff --git a/docs/source/notebooks/ultimate_bending.ipynb b/docs/source/notebooks/ultimate_bending.ipynb index 1d64cf84..98284947 100644 --- a/docs/source/notebooks/ultimate_bending.ipynb +++ b/docs/source/notebooks/ultimate_bending.ipynb @@ -24,7 +24,7 @@ "outputs": [], "source": [ "import numpy as np\n", - "from concreteproperties.material import Concrete, Steel\n", + "from concreteproperties.material import Concrete, SteelBar\n", "import concreteproperties.stress_strain_profile as ssp\n", "import sectionproperties.pre.library.primitive_sections as sp_ps\n", "from concreteproperties.pre import add_bar_rectangular_array\n", @@ -62,7 +62,7 @@ " colour=\"lightgrey\",\n", ")\n", "\n", - "steel = Steel(\n", + "steel = SteelBar(\n", " name=\"500 MPa Steel\",\n", " density=7.85e-6,\n", " stress_strain_profile=ssp.SteelElasticPlastic(\n", From 3ea5b5f59bfbbcd6c3281b3803bcc0944c29f03f Mon Sep 17 00:00:00 2001 From: robbievanleeuwen Date: Tue, 2 Aug 2022 17:52:03 +1000 Subject: [PATCH 2/7] Fix ultimate and moment interaction --- concreteproperties/concrete_section.py | 36 ++++++++++--------- concreteproperties/pre.py | 4 +-- concreteproperties/stress_strain_profile.py | 20 +++++++++++ docs/source/notebooks/biaxial_bending.ipynb | 4 +-- .../source/notebooks/moment_interaction.ipynb | 4 +-- 5 files changed, 46 insertions(+), 22 deletions(-) diff --git a/concreteproperties/concrete_section.py b/concreteproperties/concrete_section.py index 2c801070..5127bbfe 100644 --- a/concreteproperties/concrete_section.py +++ b/concreteproperties/concrete_section.py @@ -741,6 +741,8 @@ def ultimate_bending_capacity( 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 @@ -861,10 +863,6 @@ def calculate_ultimate_section_actions( meshed_split_geoms.extend(split_geoms) - # # TESTING - # for geom in meshed_split_geoms: - # geom.plot_geometry() - # initialise results n = 0 m_x = 0 @@ -905,7 +903,9 @@ def calculate_ultimate_section_actions( ) # calculate stress and force - stress = lumped_geom.material.stress_strain_profile.get_stress(strain=strain) + stress = lumped_geom.material.stress_strain_profile.get_stress( + strain=strain + ) force = stress * area n += force @@ -1011,7 +1011,9 @@ def moment_interaction_diagram( mi_results = res.MomentInteractionResults() # compute extreme tensile fibre - _, d_t = utils.calculate_extreme_fibre(points=self.geometry.points, theta=theta) + _, 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): @@ -1745,8 +1747,9 @@ def extreme_bar( self, theta: float, ) -> Tuple[float, float]: - r"""Given neutral axis angle ``theta``, determines the depth of the furthest bar - from the extreme compressive fibre and also returns its yield strain. + 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`) @@ -1756,19 +1759,20 @@ def extreme_bar( # initialise variables d_ext = 0 - extreme_geom = None # calculate extreme fibre in local coordinates extreme_fibre, _ = utils.calculate_extreme_fibre( - points=self.geometry.points, theta=theta + 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 tensile bar - for steel_geom in self.steel_geometries: - centroid = steel_geom.calculate_centroid() + # 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]) @@ -1778,12 +1782,12 @@ def extreme_bar( if d > d_ext: d_ext = d - extreme_geom = steel_geom + extreme_geom = lumped_geom # calculate yield strain yield_strain = ( - extreme_geom.material.stress_strain_profile.yield_strength # type: ignore - / extreme_geom.material.stress_strain_profile.elastic_modulus # type: ignore + extreme_geom.material.stress_strain_profile.get_yield_strength() + / extreme_geom.material.stress_strain_profile.get_elastic_modulus() ) return d_ext, yield_strain diff --git a/concreteproperties/pre.py b/concreteproperties/pre.py index e0caea47..55af10ea 100644 --- a/concreteproperties/pre.py +++ b/concreteproperties/pre.py @@ -300,8 +300,8 @@ def sort_polys( return top_polys, bot_polys def plot_geometry( - self, - title: str = "Cross-Section Geometry", + self, + title: str = "Cross-Section Geometry", **kwargs, ) -> matplotlib.axes.Axes: # type: ignore """Plots the geometry. diff --git a/concreteproperties/stress_strain_profile.py b/concreteproperties/stress_strain_profile.py index 2357e376..2c7ea7ea 100644 --- a/concreteproperties/stress_strain_profile.py +++ b/concreteproperties/stress_strain_profile.py @@ -125,6 +125,16 @@ def get_tensile_strength( return min(self.stresses) + def get_yield_strength( + self, + ) -> float: + """Returns the yield strength of the stress-strain profile. + + :return: Yield strength + """ + + raise NotImplementedError + def get_ultimate_compressive_strain( self, ) -> float: @@ -668,6 +678,16 @@ def get_elastic_modulus( return self.elastic_modulus + def get_yield_strength( + self, + ) -> float: + """Returns the yield strength of the stress-strain profile. + + :return: Yield strength + """ + + return self.yield_strength + def print_properties( self, fmt: str = "8.6e", diff --git a/docs/source/notebooks/biaxial_bending.ipynb b/docs/source/notebooks/biaxial_bending.ipynb index b797b764..1bad7d7d 100644 --- a/docs/source/notebooks/biaxial_bending.ipynb +++ b/docs/source/notebooks/biaxial_bending.ipynb @@ -24,7 +24,7 @@ "outputs": [], "source": [ "import numpy as np\n", - "from concreteproperties.material import Concrete, Steel\n", + "from concreteproperties.material import Concrete, SteelBar\n", "from concreteproperties.stress_strain_profile import (\n", " ConcreteLinear,\n", " RectangularStressBlock,\n", @@ -67,7 +67,7 @@ " colour=\"lightgrey\",\n", ")\n", "\n", - "steel = Steel(\n", + "steel = SteelBar(\n", " name=\"500 MPa Steel\",\n", " density=7.85e-6,\n", " stress_strain_profile=SteelElasticPlastic(\n", diff --git a/docs/source/notebooks/moment_interaction.ipynb b/docs/source/notebooks/moment_interaction.ipynb index 5fbed945..1bcea160 100644 --- a/docs/source/notebooks/moment_interaction.ipynb +++ b/docs/source/notebooks/moment_interaction.ipynb @@ -24,7 +24,7 @@ "outputs": [], "source": [ "import numpy as np\n", - "from concreteproperties.material import Concrete, Steel\n", + "from concreteproperties.material import Concrete, SteelBar\n", "from concreteproperties.stress_strain_profile import (\n", " ConcreteLinear,\n", " RectangularStressBlock,\n", @@ -66,7 +66,7 @@ " colour=\"lightgrey\",\n", ")\n", "\n", - "steel = Steel(\n", + "steel = SteelBar(\n", " name=\"500 MPa Steel\",\n", " density=7.85e-6,\n", " stress_strain_profile=SteelElasticPlastic(\n", From b6292624ea17412b264bc5492ccf2384e4880966 Mon Sep 17 00:00:00 2001 From: robbievanleeuwen Date: Wed, 3 Aug 2022 00:33:26 +1000 Subject: [PATCH 3/7] Finish overhaul --- concreteproperties/concrete_section.py | 366 ++++++++++-------- concreteproperties/design_codes.py | 35 +- concreteproperties/pre.py | 3 + concreteproperties/results.py | 184 +++++++-- concreteproperties/stress_strain_profile.py | 2 +- .../tests/test_gross_properties.py | 2 +- .../tests/test_moment_interaction.py | 4 +- .../tests/test_reinforced_concrete_basics.py | 35 +- concreteproperties/tests/test_rotation.py | 18 +- .../tests/test_stress_equilibrium.py | 16 +- docs/source/notebooks/stress_analysis.ipynb | 24 +- 11 files changed, 430 insertions(+), 259 deletions(-) diff --git a/concreteproperties/concrete_section.py b/concreteproperties/concrete_section.py index 5127bbfe..f8f2abf2 100644 --- a/concreteproperties/concrete_section.py +++ b/concreteproperties/concrete_section.py @@ -1253,7 +1253,7 @@ def calculate_uncracked_stress( uncracked section. Uses gross area section properties to determine concrete and steel stresses - given an axial force `n`, and bending moments `m_x` and `m_y`. + 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 @@ -1263,12 +1263,16 @@ def calculate_uncracked_stress( """ # initialise stress results - analysis_sections = [] + conc_sections = [] conc_sigs = [] conc_forces = [] - steel_sigs = [] - steel_strains = [] - steel_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 @@ -1288,25 +1292,24 @@ def calculate_uncracked_stress( # point on neutral axis is centroid point_na = (cx, cy) - # split concrete geometries above and below neutral axis - split_conc_geoms = [] + # split meshed geometries above and below neutral axis + split_meshed_geoms = [] - for conc_geom in self.concrete_geometries: - top_geoms, bot_geoms = utils.split_section( - geometry=conc_geom, + for meshed_geom in self.meshed_geometries: + top_geoms, bot_geoms = meshed_geom.split_section( point=point_na, theta=theta, ) - split_conc_geoms.extend(top_geoms) - split_conc_geoms.extend(bot_geoms) + split_meshed_geoms.extend(top_geoms) + split_meshed_geoms.extend(bot_geoms) - # loop through all concrete geometries and calculate stress - for conc_geom in split_conc_geoms: - analysis_section = AnalysisSection(geometry=conc_geom) + # 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_conc, d_x, d_y = analysis_section.get_elastic_stress( + sig, n_sec, d_x, d_y = analysis_section.get_elastic_stress( n=n, m_x=m_x, m_y=m_y, @@ -1317,50 +1320,59 @@ def calculate_uncracked_stress( e_iyy=e_iyy, e_ixy=e_ixy, ) - conc_sigs.append(sig) - conc_forces.append((n_conc, d_x, d_y)) - # save analysis section - analysis_sections.append(analysis_section) + # 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 steel geometries and calculate stress - for steel_geom in self.steel_geometries: - # initialise stress and position of bar + # loop through all lumped geometries and calculate stress + for lumped_geom in self.reinf_geometries_lumped: + # initialise stress and position sig = 0 - centroid = steel_geom.calculate_centroid() + centroid = lumped_geom.calculate_centroid() x = centroid[0] - cx y = centroid[1] - cy # axial stress - sig += n * steel_geom.material.elastic_modulus / e_a + sig += n * lumped_geom.material.elastic_modulus / e_a # bending moment stress - sig += steel_geom.material.elastic_modulus * ( + 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 += steel_geom.material.elastic_modulus * ( + 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 / steel_geom.material.elastic_modulus + strain = sig / lumped_geom.material.elastic_modulus # net force and point of action - n_steel = sig * steel_geom.calculate_area() + n_lumped = sig * lumped_geom.calculate_area() - steel_sigs.append(sig) - steel_strains.append(strain) - steel_forces.append((n_steel, x, y)) + 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=analysis_sections, + concrete_analysis_sections=conc_sections, concrete_stresses=conc_sigs, concrete_forces=conc_forces, - steel_geometries=self.steel_geometries, - steel_stresses=steel_sigs, - steel_strains=steel_strains, - steel_forces=steel_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, ) def calculate_cracked_stress( @@ -1372,9 +1384,9 @@ def calculate_cracked_stress( """Calculates stresses within the reinforced concrete section assuming a cracked section. - Uses cracked area section properties to determine concrete and steel stresses - given an axial force `n` and bending moment `m` about the bending axis stored - in `cracked_results`. + 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 @@ -1384,12 +1396,16 @@ def calculate_cracked_stress( """ # initialise stress results - analysis_sections = [] + conc_sections = [] conc_sigs = [] conc_forces = [] - steel_sigs = [] - steel_strains = [] - steel_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 @@ -1427,26 +1443,13 @@ def calculate_cracked_stress( m_x = sign * np.sqrt(m * m / (1 + 1 / (c * c))) m_y = m_x / c - # depth of neutral axis at extreme tensile fibre - extreme_fibre, d_t = utils.calculate_extreme_fibre( - points=self.geometry.points, theta=theta - ) - - # find point on neutral axis by shifting by d_n - point_na = utils.point_on_neutral_axis( - extreme_fibre=extreme_fibre, d_n=cracked_results.d_nc, theta=theta - ) - - # get principal coordinates of neutral axis - na_local = utils.global_to_local(theta=theta, x=point_na[0], y=point_na[1]) - - # loop through all concrete geometries and calculate stress + # loop through all meshed geometries and calculate stress for geom in cracked_results.cracked_geometries: - if isinstance(geom.material, Concrete): + if geom.material.meshed: analysis_section = AnalysisSection(geometry=geom) # calculate stress, force and point of action - sig, n_conc, d_x, d_y = analysis_section.get_elastic_stress( + sig, n_sec, d_x, d_y = analysis_section.get_elastic_stress( n=n, m_x=m_x, m_y=m_y, @@ -1457,50 +1460,59 @@ def calculate_cracked_stress( e_iyy=e_iyy, e_ixy=e_ixy, ) - conc_sigs.append(sig) - conc_forces.append((n_conc, d_x, d_y)) - # save analysis section - analysis_sections.append(analysis_section) + # 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 steel geometries and calculate stress - for steel_geom in self.steel_geometries: + # 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 = steel_geom.calculate_centroid() + centroid = lumped_geom.calculate_centroid() x = centroid[0] - cx y = centroid[1] - cy # axial stress - sig += n * steel_geom.material.elastic_modulus / e_a + sig += n * lumped_geom.material.elastic_modulus / e_a # bending moment stress - sig += steel_geom.material.elastic_modulus * ( + 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 += steel_geom.material.elastic_modulus * ( + 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 / steel_geom.material.elastic_modulus + strain = sig / lumped_geom.material.elastic_modulus # net force and point of action - n_steel = sig * steel_geom.calculate_area() + n_lumped = sig * lumped_geom.calculate_area() - steel_sigs.append(sig) - steel_strains.append(strain) - steel_forces.append((n_steel, x, y)) + 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=analysis_sections, + concrete_analysis_sections=conc_sections, concrete_stresses=conc_sigs, concrete_forces=conc_forces, - steel_geometries=self.steel_geometries, - steel_stresses=steel_sigs, - steel_strains=steel_strains, - steel_forces=steel_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, ) def calculate_service_stress( @@ -1513,13 +1525,13 @@ def calculate_service_stress( 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, can provided a curvature which overrides the + 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 - plots the stress at the given curvature + calculates the stress at the given curvature :return: Stress results object """ @@ -1537,7 +1549,7 @@ def calculate_service_stress( # set neutral axis depth limits # depth of neutral axis at extreme tensile fibre extreme_fibre, d_t = utils.calculate_extreme_fibre( - points=self.geometry.points, theta=theta + 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 @@ -1559,51 +1571,65 @@ def calculate_service_stress( d_n = 0 # initialise stress results - analysis_sections = [] + conc_sections = [] conc_sigs = [] conc_forces = [] - steel_sigs = [] - steel_strains = [] - steel_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 concrete geometries at points in stress-strain profiles - concrete_split_geoms = utils.split_section_at_strains( - concrete_geometries=self.concrete_geometries, - theta=theta, - point_na=point_na, - ultimate=False, - kappa=kappa, - ) + # create splits in meshed geometries at points in stress-strain profiles + meshed_split_geoms: List[Union[CPGeom, CPGeomConcrete]] = [] - # loop through all concrete geometries and calculate stress - for geom in concrete_split_geoms: - analysis_section = AnalysisSection(geometry=geom) + 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_conc, d_x, d_y = analysis_section.get_service_stress( + 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), ) - conc_sigs.append(sig) - conc_forces.append((n_conc, d_x, d_y)) - # save analysis section - analysis_sections.append(analysis_section) + # 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 steel geometries and calculate stress - for steel_geom in self.steel_geometries: - # get position of bar - centroid = steel_geom.calculate_centroid() + # 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 steel + # get strain at centroid of lump strain = utils.get_service_strain( point=(centroid[0], centroid[1]), point_na=point_na, @@ -1612,28 +1638,32 @@ def calculate_service_stress( ) # calculate stress, force and point of action - sig = steel_geom.material.stress_strain_profile.get_stress(strain=strain) - n_steel = sig * steel_geom.calculate_area() + sig = lumped_geom.material.stress_strain_profile.get_stress(strain=strain) + n_lumped = sig * lumped_geom.calculate_area() - steel_sigs.append(sig) - steel_strains.append(strain) - steel_forces.append( + lumped_reinf_sigs.append(sig) + lumped_reinf_strains.append(strain) + lumped_reinf_forces.append( ( - n_steel, - centroid[0] - -self.gross_properties.cx, - centroid[1] - -self.gross_properties.cy, + 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=analysis_sections, + concrete_analysis_sections=conc_sections, concrete_stresses=conc_sigs, concrete_forces=conc_forces, - steel_geometries=self.steel_geometries, - steel_stresses=steel_sigs, - steel_strains=steel_strains, - steel_forces=steel_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, ) def calculate_ultimate_stress( @@ -1649,7 +1679,7 @@ def calculate_ultimate_stress( # depth of neutral axis at extreme tensile fibre extreme_fibre, _ = utils.calculate_extreme_fibre( - points=self.geometry.points, theta=ultimate_results.theta + points=self.compound_geometry.points, theta=ultimate_results.theta ) # find point on neutral axis by shifting by d_n @@ -1663,48 +1693,62 @@ def calculate_ultimate_stress( ) # initialise stress results for each concrete geometry - analysis_sections = [] + conc_sections = [] conc_sigs = [] conc_forces = [] - steel_sigs = [] - steel_strains = [] - steel_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]] = [] - # create splits in concrete geometries at points in stress-strain profiles if isinf(ultimate_results.d_n): - concrete_split_geoms = self.concrete_geometries + meshed_split_geoms = self.meshed_geometries else: - concrete_split_geoms = utils.split_section_at_strains( - concrete_geometries=self.concrete_geometries, - theta=ultimate_results.theta, - point_na=point_na, - ultimate=True, - ultimate_strain=self.gross_properties.conc_ultimate_strain, - d_n=ultimate_results.d_n, - ) + 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 geom in concrete_split_geoms: - analysis_section = AnalysisSection(geometry=geom) + for meshed_geom in meshed_split_geoms: + analysis_section = AnalysisSection(geometry=meshed_geom) # calculate stress, force and point of action - sig, n_conc, d_x, d_y = analysis_section.get_ultimate_stress( + 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), ) - conc_sigs.append(sig) - conc_forces.append((n_conc, d_x, d_y)) - # save analysis section - analysis_sections.append(analysis_section) + # 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 steel geometries and calculate stress - for steel_geom in self.steel_geometries: - # get position of bar - centroid = steel_geom.calculate_centroid() + # 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 steel if isinf(ultimate_results.d_n): @@ -1719,28 +1763,32 @@ def calculate_ultimate_stress( ) # calculate stress, force and point of action - sig = steel_geom.material.stress_strain_profile.get_stress(strain=strain) - n_steel = sig * steel_geom.calculate_area() + sig = lumped_geom.material.stress_strain_profile.get_stress(strain=strain) + n_lumped = sig * lumped_geom.calculate_area() - steel_sigs.append(sig) - steel_strains.append(strain) - steel_forces.append( + lumped_reinf_sigs.append(sig) + lumped_reinf_strains.append(strain) + lumped_reinf_forces.append( ( - n_steel, + 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=analysis_sections, + concrete_analysis_sections=conc_sections, concrete_stresses=conc_sigs, concrete_forces=conc_forces, - steel_geometries=self.steel_geometries, - steel_stresses=steel_sigs, - steel_strains=steel_strains, - steel_forces=steel_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, ) def extreme_bar( @@ -1854,7 +1902,7 @@ def plot_section( plotted_materials.append(lumped_geom.material) # plot the points and facets - coords = list(lumped_geom.geom.exterior.coords) + coords = list(lumped_geom.geom.exterior.coords) # type: ignore bar = mpatches.Polygon( xy=coords, closed=False, color=lumped_geom.material.colour ) diff --git a/concreteproperties/design_codes.py b/concreteproperties/design_codes.py index cf384383..908296da 100644 --- a/concreteproperties/design_codes.py +++ b/concreteproperties/design_codes.py @@ -2,7 +2,8 @@ from copy import deepcopy from math import inf -from typing import TYPE_CHECKING, List, Optional, Tuple +from multiprocessing.sharedctypes import Value +from typing import TYPE_CHECKING, List, Tuple import numpy as np from rich.live import Live @@ -12,7 +13,7 @@ import concreteproperties.results as res import concreteproperties.stress_strain_profile as ssp import concreteproperties.utils as utils -from concreteproperties.material import Concrete, Steel +from concreteproperties.material import Concrete, SteelBar if TYPE_CHECKING: from concreteproperties.concrete_section import ConcreteSection @@ -60,7 +61,7 @@ def create_steel_material( self, yield_strength: float, colour: str = "grey", - ) -> Steel: + ) -> SteelBar: """Returns a steel material object. List assumptions of material properties here... @@ -232,7 +233,13 @@ def calculate_ultimate_stress( class AS3600(DesignCode): - """Design code class for Australian standard AS 3600:2018.""" + """Design code class for Australian standard AS 3600:2018. + + Note that this design code only supports :class:`~concreteproperties.pre.Concrete` + and :class:`~concreteproperties.pre.SteelBar` material objects. Meshed + :class:`~concreteproperties.pre.Steel` material objects are **not** supported + as this falls under the composite structures design code. + """ def __init__( self, @@ -252,11 +259,15 @@ def assign_concrete_section( 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.steel_geometries: - if steel_geom.material.stress_strain_profile.fracture_strain < 0.05: + 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 @@ -343,7 +354,7 @@ def create_steel_material( yield_strength: float = 500, ductility_class: str = "N", colour: str = "grey", - ) -> Steel: + ) -> SteelBar: r"""Returns a steel material object. | **Material assumptions:** @@ -367,7 +378,7 @@ def create_steel_material( else: raise ValueError("ductility_class must be N or L.") - return Steel( + return SteelBar( name=f"{yield_strength:.0f} MPa Steel (AS 3600:2018)", density=7.85e-6, stress_strain_profile=ssp.SteelElasticPlastic( @@ -406,7 +417,7 @@ def squash_tensile_load( squash_load += force_c # loop through all steel geometries - for steel_geom in self.concrete_section.steel_geometries: + for steel_geom in self.concrete_section.reinf_geometries_lumped: # calculate area and centroid area = steel_geom.calculate_area() @@ -415,7 +426,7 @@ def squash_tensile_load( strain=0.025 ) - force_t = -area * steel_geom.material.stress_strain_profile.yield_strength + force_t = -area * steel_geom.material.stress_strain_profile.get_yield_strength() # add to totals squash_load += force_c @@ -500,10 +511,10 @@ def get_n_ub( # get compressive strain at extreme fibre eps_cu = self.concrete_section.gross_properties.conc_ultimate_strain - # 4) calculate d_n at balanced load + # calculate d_n at balanced load d_nb = d_0 * (eps_cu) / (eps_sy + eps_cu) - # 5) calculate axial force at balanced load + # 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) ) diff --git a/concreteproperties/pre.py b/concreteproperties/pre.py index 55af10ea..20146b45 100644 --- a/concreteproperties/pre.py +++ b/concreteproperties/pre.py @@ -178,6 +178,9 @@ def split_section( :return: Geometries above and below the line """ + # round point + point = np.round(point, 6) + # generate unit vector vector = np.cos(theta), np.sin(theta) diff --git a/concreteproperties/results.py b/concreteproperties/results.py index 3eda860e..b4b2b809 100644 --- a/concreteproperties/results.py +++ b/concreteproperties/results.py @@ -1016,36 +1016,48 @@ class StressResult: :meth:`~concreteproperties.analysis_section.AnalysisSection.plot_shape` :param concrete_stresses: List of concrete stresses at the nodes of each concrete analysis section - :param concrete_forces: List of net forces for each concrete analysis section and its - lever arm (``force``, ``d_x``, ``d_y``) - :param steel_geometries: List of steel geometry objects present in the stress analysis - :param steel_stresses: List of steel stresses for each steel geometry - :param steel_strains: List of steel strains for each steel geometry - :param steel_forces: List of net forces for each steel geometry and its lever arm - (``force``, ``d_x``, ``d_y``) + :param concrete_forces: List of net forces for each concrete analysis section and + its lever arm (``force``, ``d_x``, ``d_y``) + :param meshed_reinforcement_sections: List of meshed reinforcement section objects + present in the stress analysis + :param meshed_reinforcement_stresses: List of meshed reinforcement stresses at the + nodes of each meshed reinforcement analysis section + :param meshed_reinforcement_forces: List of net forces for each meshed reinforcement + analysis section and its lever arm (``force``, ``d_x``, ``d_y``) + :param lumped_reinforcement_geometries: List of lumped reinforcement geometry + objects present in the stress analysis + :param lumped_reinforcement_stresses: List of lumped reinforcement stresses for + each lumped geometry + :param lumped_reinforcement_strains: List of lumped reinforcement strains for each + lumped geometry + :param lumped_reinforcement_forces: List of net forces for each lumped reinforcement + geometry and its lever arm (``force``, ``d_x``, ``d_y``) """ concrete_section: ConcreteSection concrete_analysis_sections: List[AnalysisSection] concrete_stresses: List[np.ndarray] concrete_forces: List[Tuple[float, float, float]] - steel_geometries: List[CPGeom] - steel_stresses: List[float] - steel_strains: List[float] - steel_forces: List[Tuple[float, float, float]] + meshed_reinforcement_sections: List[AnalysisSection] + meshed_reinforcement_stresses: List[np.ndarray] + meshed_reinforcement_forces: List[Tuple[float, float, float]] + lumped_reinforcement_geometries: List[CPGeom] + lumped_reinforcement_stresses: List[float] + lumped_reinforcement_strains: List[float] + lumped_reinforcement_forces: List[Tuple[float, float, float]] def plot_stress( self, title: str = "Stress", conc_cmap: str = "RdGy", - steel_cmap: str = "bwr", + reinf_cmap: str = "bwr", **kwargs, ) -> matplotlib.axes.Axes: # type: ignore """Plots concrete and steel stresses on a concrete section. :param title: Plot title :param conc_cmap: Colour map for the concrete stress - :param steel_cmap: Colour map for the steel stress + :param reinf_cmap: Colour map for the reinforcement stress :param kwargs: Passed to :func:`~concreteproperties.post.plotting_context` :return: Matplotlib axes object @@ -1064,18 +1076,55 @@ def plot_stress( # set up the colormaps cmap_conc = cm.get_cmap(name=conc_cmap) - cmap_steel = cm.get_cmap(name=steel_cmap) + cmap_reinf = cm.get_cmap(name=reinf_cmap) # determine minimum and maximum stress values for the contour list # add tolerance for plotting stress blocks conc_sig_min = min([min(x) for x in self.concrete_stresses]) - 1e-12 conc_sig_max = max([max(x) for x in self.concrete_stresses]) + 1e-12 - steel_sig_min = min(self.steel_stresses) - steel_sig_max = max(self.steel_stresses) + + # if there is meshed reinforcement, calculate min and max + if self.meshed_reinforcement_stresses: + meshed_reinf_sig_min = ( + min([min(x) for x in self.meshed_reinforcement_stresses]) - 1e-12 + ) + meshed_reinf_sig_max = ( + max([max(x) for x in self.meshed_reinforcement_stresses]) + 1e-12 + ) + else: + meshed_reinf_sig_min = None + meshed_reinf_sig_max = None + + # if there is lumped reinforcement, calculate min and max + if self.lumped_reinforcement_stresses: + lumped_reinf_sig_min = min(self.lumped_reinforcement_stresses) + lumped_reinf_sig_max = max(self.lumped_reinforcement_stresses) + else: + lumped_reinf_sig_min = None + lumped_reinf_sig_max = None + + # determine min and max reinforcement stresess + if ( + meshed_reinf_sig_min + and meshed_reinf_sig_max + and lumped_reinf_sig_min + and lumped_reinf_sig_max + ): + reinf_sig_min = min(meshed_reinf_sig_min, lumped_reinf_sig_min) + reinf_sig_max = max(meshed_reinf_sig_max, lumped_reinf_sig_max) + elif meshed_reinf_sig_min and meshed_reinf_sig_max: + reinf_sig_min = meshed_reinf_sig_min + reinf_sig_max = meshed_reinf_sig_max + elif lumped_reinf_sig_min and lumped_reinf_sig_max: + reinf_sig_min = lumped_reinf_sig_min + reinf_sig_max = lumped_reinf_sig_max + else: + reinf_sig_min = 0 + reinf_sig_max = 0 # set up ticks v_conc = np.linspace(conc_sig_min, conc_sig_max, 15, endpoint=True) - v_steel = np.linspace(steel_sig_min, steel_sig_max, 15, endpoint=True) + v_reinf = np.linspace(reinf_sig_min, reinf_sig_max, 15, endpoint=True) if np.isclose(v_conc[0], v_conc[-1], atol=1e-12): v_conc = 15 @@ -1083,12 +1132,12 @@ def plot_stress( else: ticks_conc = v_conc - if np.isclose(v_steel[0], v_steel[-1], atol=1e-12): - ticks_steel = None - steel_tick_same = True + if np.isclose(v_reinf[0], v_reinf[-1], atol=1e-12): + ticks_reinf = None + reinf_tick_same = True else: - ticks_steel = v_steel - steel_tick_same = False + ticks_reinf = v_reinf + reinf_tick_same = False # plot the concrete stresses for idx, sig in enumerate(self.concrete_stresses): @@ -1134,26 +1183,66 @@ def plot_stress( triang, sig, [zero_level], linewidths=1, linestyles="dashed" ) - # plot the steel stresses - steel_patches = [] + # plot the meshed reinforcement stresses + for idx, sig in enumerate(self.meshed_reinforcement_stresses): + # check region has a force + if abs(self.meshed_reinforcement_forces[idx][0]) > 1e-8: + # create triangulation + triang = tri.Triangulation( + self.meshed_reinforcement_sections[idx].mesh_nodes[:, 0], + self.meshed_reinforcement_sections[idx].mesh_nodes[:, 1], + self.meshed_reinforcement_sections[idx].mesh_elements[:, 0:3], # type: ignore + ) + + # plot the filled contour + trictr = fig.axes[0].tricontourf( + triang, sig, v_reinf, cmap=cmap_reinf, norm=CenteredNorm() + ) + + # plot a zero stress contour, supressing warning + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="No contour levels were found within the data range.", + ) + + # set zero stress for neutral axis contour + zero_level = 0 + + if min(sig) > 0: + if min(sig) < 1e-3: + zero_level = min(sig) + 1e-12 + + if max(sig) < 0: + if max(sig) > -1e-3: + zero_level = max(sig) - 1e-12 + + if min(sig) == 0: + zero_level = 1e-12 + + if max(sig) == 0: + zero_level = -1e-12 + + CS = fig.axes[0].tricontour( + triang, sig, [zero_level], linewidths=1, linestyles="dashed" + ) + + # plot the lumped reinforcement stresses + lumped_reinf_patches = [] colours = [] - for idx, sig in enumerate(self.steel_stresses): - steel_patches.append( + for idx, sig in enumerate(self.lumped_reinforcement_stresses): + lumped_reinf_patches.append( mpatches.Polygon( - xy=list( - self.concrete_section.steel_geometries[ - idx - ].geom.exterior.coords - ) + xy=list(self.lumped_reinforcement_geometries[idx].geom.exterior.coords) # type: ignore ) ) colours.append(sig) - patch = PatchCollection(steel_patches, cmap=cmap_steel) + patch = PatchCollection(lumped_reinf_patches, cmap=cmap_reinf) patch.set_array(colours) - if steel_tick_same: - patch.set_clim([0.99 * v_steel[0], 1.01 * v_steel[-1]]) + if reinf_tick_same: + patch.set_clim([0.99 * v_reinf[0], 1.01 * v_reinf[-1]]) fig.axes[0].add_collection(patch) # add the colour bars @@ -1166,9 +1255,9 @@ def plot_stress( ) fig.colorbar( patch, - label="Steel Stress", + label="Reinforcement Stress", format="%.2e", - ticks=ticks_steel, + ticks=ticks_reinf, cax=fig.axes[2], ) @@ -1190,9 +1279,13 @@ def sum_forces( for conc_force in self.concrete_forces: force_sum += conc_force[0] - # sum steel forces - for steel_force in self.steel_forces: - force_sum += steel_force[0] + # sum meshed reinf stresses + for meshed_reinf_force in self.meshed_reinforcement_forces: + force_sum += meshed_reinf_force[0] + + # sum lumped reinf forces + for lumped_reinf_force in self.lumped_reinforcement_forces: + force_sum += lumped_reinf_force[0] return force_sum @@ -1213,10 +1306,15 @@ def sum_moments( moment_sum_x += conc_force[0] * conc_force[2] moment_sum_y += conc_force[0] * conc_force[1] - # sum steel forces - for steel_force in self.steel_forces: - moment_sum_x += steel_force[0] * steel_force[2] - moment_sum_y += steel_force[0] * steel_force[1] + # sum meshed reinf stresses + for meshed_reinf_force in self.meshed_reinforcement_forces: + moment_sum_x += meshed_reinf_force[0] * meshed_reinf_force[2] + moment_sum_y += meshed_reinf_force[0] * meshed_reinf_force[1] + + # sum lumped reinf forces + for lumped_reinf_force in self.lumped_reinforcement_forces: + moment_sum_x += lumped_reinf_force[0] * lumped_reinf_force[2] + moment_sum_y += lumped_reinf_force[0] * lumped_reinf_force[1] moment_sum = np.sqrt(moment_sum_x * moment_sum_x + moment_sum_y * moment_sum_y) diff --git a/concreteproperties/stress_strain_profile.py b/concreteproperties/stress_strain_profile.py index 2c7ea7ea..40d2730d 100644 --- a/concreteproperties/stress_strain_profile.py +++ b/concreteproperties/stress_strain_profile.py @@ -565,7 +565,7 @@ def get_stress( :return: Stress """ - if strain >= self.strains[1] - 1e-12: + if strain >= self.strains[1] - 1e-8: return self.stresses[2] else: return 0 diff --git a/concreteproperties/tests/test_gross_properties.py b/concreteproperties/tests/test_gross_properties.py index c97ab5da..4c63d8c3 100644 --- a/concreteproperties/tests/test_gross_properties.py +++ b/concreteproperties/tests/test_gross_properties.py @@ -8,7 +8,7 @@ def test_rectangle_second_moment_of_area(): b = 100 rect = sp_ps.rectangular_section(d=d, b=b) - sec = AnalysisSection(geometry=rect) + sec = AnalysisSection(geometry=rect) # type: ignore ixx_g = 0 iyy_g = 0 diff --git a/concreteproperties/tests/test_moment_interaction.py b/concreteproperties/tests/test_moment_interaction.py index b7b3d4f8..0a3751d8 100644 --- a/concreteproperties/tests/test_moment_interaction.py +++ b/concreteproperties/tests/test_moment_interaction.py @@ -1,7 +1,7 @@ import numpy as np import pytest from concreteproperties.concrete_section import ConcreteSection -from concreteproperties.material import Concrete, Steel +from concreteproperties.material import Concrete, SteelBar from concreteproperties.stress_strain_profile import ( ConcreteLinear, RectangularStressBlock, @@ -24,7 +24,7 @@ colour="lightgrey", ) -steel = Steel( +steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( diff --git a/concreteproperties/tests/test_reinforced_concrete_basics.py b/concreteproperties/tests/test_reinforced_concrete_basics.py index 6f2947a6..73f41303 100644 --- a/concreteproperties/tests/test_reinforced_concrete_basics.py +++ b/concreteproperties/tests/test_reinforced_concrete_basics.py @@ -2,7 +2,7 @@ import pytest import sectionproperties.pre.library.primitive_sections as sp_ps from concreteproperties.concrete_section import ConcreteSection -from concreteproperties.material import Concrete, Steel +from concreteproperties.material import Concrete, SteelBar from concreteproperties.pre import add_bar, add_bar_rectangular_array from concreteproperties.stress_strain_profile import ( ConcreteLinear, @@ -31,7 +31,7 @@ def test_example_3_1(): colour="lightgrey", ) - steel = Steel( + steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( @@ -72,7 +72,10 @@ def test_example_3_1(): assert pytest.approx(cracked_results.ixx_c_cr, rel=0.01) == 821e6 assert pytest.approx(cracked_results.iuu_cr, rel=0.01) == 821e6 assert pytest.approx(max(cracked_stress.concrete_stresses[0]), rel=0.01) == 15.3 - assert pytest.approx(min(cracked_stress.steel_stresses), rel=0.01) == -213 + assert ( + pytest.approx(min(cracked_stress.lumped_reinforcement_stresses), rel=0.01) + == -213 + ) def test_example_3_2(): @@ -91,7 +94,7 @@ def test_example_3_2(): colour="lightgrey", ) - steel = Steel( + steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( @@ -153,9 +156,17 @@ def test_example_3_2(): conc_stresses.extend(cs) assert pytest.approx(max(conc_stresses), rel=0.01) == 8.1 - assert pytest.approx(max(cracked_stress.steel_stresses), rel=0.02) == 33 - assert pytest.approx(min(cracked_stress.steel_stresses), rel=0.01) == -193 - assert pytest.approx(cracked_stress.steel_stresses[-1], rel=0.01) == -173 + assert ( + pytest.approx(max(cracked_stress.lumped_reinforcement_stresses), rel=0.02) == 33 + ) + assert ( + pytest.approx(min(cracked_stress.lumped_reinforcement_stresses), rel=0.01) + == -193 + ) + assert ( + pytest.approx(cracked_stress.lumped_reinforcement_stresses[-1], rel=0.01) + == -173 + ) def test_example_3_4(): @@ -174,7 +185,7 @@ def test_example_3_4(): colour="lightgrey", ) - steel = Steel( + steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( @@ -238,7 +249,7 @@ def test_example_3_8(): colour="lightgrey", ) - steel = Steel( + steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( @@ -287,7 +298,7 @@ def test_example_3_9(): colour="lightgrey", ) - steel = Steel( + steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( @@ -336,7 +347,7 @@ def test_example_3_11(): colour="lightgrey", ) - steel = Steel( + steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( @@ -388,7 +399,7 @@ def test_example_5_2(): colour="lightgrey", ) - steel = Steel( + steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( diff --git a/concreteproperties/tests/test_rotation.py b/concreteproperties/tests/test_rotation.py index f871d2d3..9983af91 100644 --- a/concreteproperties/tests/test_rotation.py +++ b/concreteproperties/tests/test_rotation.py @@ -1,7 +1,7 @@ import numpy as np import pytest from concreteproperties.concrete_section import ConcreteSection -from concreteproperties.material import Concrete, Steel +from concreteproperties.material import Concrete, SteelBar from concreteproperties.stress_strain_profile import ( ConcreteLinear, RectangularStressBlock, @@ -28,7 +28,7 @@ colour="lightgrey", ) -steel = Steel( +steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( @@ -72,7 +72,7 @@ def test_rotated_gross_properties(theta): pytest.approx(new_gross_results.concrete_area) == ref_gross_results.concrete_area ) - assert pytest.approx(new_gross_results.steel_area) == ref_gross_results.steel_area + assert pytest.approx(new_gross_results.reinf_lumped_area) == ref_gross_results.reinf_lumped_area assert pytest.approx(new_gross_results.e_a) == ref_gross_results.e_a assert pytest.approx(new_gross_results.mass) == ref_gross_results.mass assert pytest.approx(new_gross_results.perimeter) == ref_gross_results.perimeter @@ -147,8 +147,8 @@ def test_rotated_uncracked_stress(theta): assert pytest.approx(cf[0]) == ref_uncr_stress.concrete_forces[i][0] - for idx, sf in enumerate(new_uncr_stress.steel_forces): - assert pytest.approx(sf[0]) == ref_uncr_stress.steel_forces[idx][0] + for idx, sf in enumerate(new_uncr_stress.lumped_reinforcement_forces): + assert pytest.approx(sf[0]) == ref_uncr_stress.lumped_reinforcement_forces[idx][0] # list of normal forces @@ -176,8 +176,8 @@ def test_rotated_cracked_stress(theta): for idx, cf in enumerate(new_cr_stress.concrete_forces): assert pytest.approx(cf[0]) == ref_cr_stress.concrete_forces[idx][0] - for idx, sf in enumerate(new_cr_stress.steel_forces): - assert pytest.approx(sf[0]) == ref_cr_stress.steel_forces[idx][0] + for idx, sf in enumerate(new_cr_stress.lumped_reinforcement_forces): + assert pytest.approx(sf[0]) == ref_cr_stress.lumped_reinforcement_forces[idx][0] # list of normal forces @@ -206,5 +206,5 @@ def test_rotated_ultimate_stress(theta): pytest.approx(cf[0], rel=5e-5) == ref_ult_stress.concrete_forces[idx][0] ) - for idx, sf in enumerate(new_ult_stress.steel_forces): - assert pytest.approx(sf[0], rel=5e-4) == ref_ult_stress.steel_forces[idx][0] + for idx, sf in enumerate(new_ult_stress.lumped_reinforcement_forces): + assert pytest.approx(sf[0], rel=5e-4) == ref_ult_stress.lumped_reinforcement_forces[idx][0] diff --git a/concreteproperties/tests/test_stress_equilibrium.py b/concreteproperties/tests/test_stress_equilibrium.py index cdf54920..e048afa9 100644 --- a/concreteproperties/tests/test_stress_equilibrium.py +++ b/concreteproperties/tests/test_stress_equilibrium.py @@ -2,7 +2,7 @@ import pytest import sectionproperties.pre.library.concrete_sections as sp_cs from concreteproperties.concrete_section import ConcreteSection -from concreteproperties.material import Concrete, Steel +from concreteproperties.material import Concrete, SteelBar from concreteproperties.stress_strain_profile import ( ConcreteLinear, RectangularStressBlock, @@ -28,7 +28,7 @@ colour="lightgrey", ) -steel = Steel( +steel = SteelBar( name="500 MPa Steel", density=7.85e-6, stress_strain_profile=SteelElasticPlastic( @@ -69,7 +69,7 @@ def test_stress_equilibrium_rectangle(theta): # check section equilibirum for uncracked stress uncr_stress = sec.calculate_uncracked_stress(n=nf, m_x=m_star) - assert pytest.approx(uncr_stress.sum_forces(), abs=1e-8) == nf + assert pytest.approx(uncr_stress.sum_forces(), abs=1e-3) == nf assert pytest.approx(uncr_stress.sum_moments()[2], rel=1e-3) == m_star # check section equilibirum for cracked stress @@ -77,7 +77,7 @@ def test_stress_equilibrium_rectangle(theta): cracked_results=cracked, n=nf, m=m_star ) - assert pytest.approx(cr_stress.sum_forces(), abs=1e-8) == nf + assert pytest.approx(cr_stress.sum_forces(), abs=1e-3) == nf assert pytest.approx(cr_stress.sum_moments()[2], rel=5e-3) == m_star # check section equilibirum for ultimate stress @@ -115,13 +115,13 @@ def test_stress_equilibrium_circular(nf): # check section equilibirum for uncracked stress uncr_stress = sec.calculate_uncracked_stress(n=nf, m_x=m_star) - assert pytest.approx(uncr_stress.sum_forces(), abs=1e-8) == nf + assert pytest.approx(uncr_stress.sum_forces(), abs=1e-3) == nf assert pytest.approx(uncr_stress.sum_moments()[2], rel=1e-3) == m_star # check section equilibirum for cracked stress cr_stress = sec.calculate_cracked_stress(cracked_results=cracked, n=nf, m=m_star) - assert pytest.approx(cr_stress.sum_forces(), abs=1e-8) == nf + assert pytest.approx(cr_stress.sum_forces(), abs=1e-3) == nf assert pytest.approx(cr_stress.sum_moments()[2], rel=5e-3) == m_star # check section equilibirum for ultimate stress @@ -163,7 +163,7 @@ def test_stress_equilibrium_tee(theta): # check section equilibirum for uncracked stress uncr_stress = sec.calculate_uncracked_stress(n=nf, m_x=m_star) - assert pytest.approx(uncr_stress.sum_forces(), abs=1e-8) == nf + assert pytest.approx(uncr_stress.sum_forces(), abs=1e-3) == nf assert pytest.approx(uncr_stress.sum_moments()[2], rel=1e-3) == m_star # check section equilibirum for cracked stress @@ -171,7 +171,7 @@ def test_stress_equilibrium_tee(theta): cracked_results=cracked, n=nf, m=m_star ) - assert pytest.approx(cr_stress.sum_forces(), abs=1e-8) == nf + assert pytest.approx(cr_stress.sum_forces(), abs=1e-3) == nf assert pytest.approx(cr_stress.sum_moments()[2], rel=5e-3) == m_star # check section equilibirum for ultimate stress diff --git a/docs/source/notebooks/stress_analysis.ipynb b/docs/source/notebooks/stress_analysis.ipynb index 19a22e7c..5ffef252 100644 --- a/docs/source/notebooks/stress_analysis.ipynb +++ b/docs/source/notebooks/stress_analysis.ipynb @@ -25,7 +25,7 @@ "source": [ "import numpy as np\n", "from rich.pretty import pprint\n", - "from concreteproperties.material import Concrete, Steel\n", + "from concreteproperties.material import Concrete, SteelBar\n", "from concreteproperties.stress_strain_profile import (\n", " EurocodeNonLinear,\n", " RectangularStressBlock,\n", @@ -73,7 +73,7 @@ " colour=\"lightgrey\",\n", ")\n", "\n", - "steel = Steel(\n", + "steel = SteelBar(\n", " name=\"500 MPa Steel\",\n", " density=7.85e-6,\n", " stress_strain_profile=SteelElasticPlastic(\n", @@ -225,7 +225,7 @@ "metadata": {}, "source": [ "### Concrete Results\n", - "Concrete geometries that remain after cracking are stored in the ``concrete_analysis_sections`` attribute as a list of ``AnalysisSection`` objects. We can visualise this geometry by calling the ``plot_mesh()`` method. The ``StressResult`` object also stores the stresses at the nodes in each ``AnalysisSection`` in the ``concrete_stresses`` attribute - these are the stresses you see on the ``plot_stress()`` plots. Finally, the ``concrete_forces`` attribute stores the resultant force for each ``AnalysisSection`` and its lever arm about the neutral axis." + "Concrete geometries that remain after cracking are stored in the ``concrete_analysis_sections`` attribute as a list of ``AnalysisSection`` objects. We can visualise this geometry by calling the ``plot_mesh()`` method. The ``StressResult`` object also stores the stresses at the nodes in each ``AnalysisSection`` in the ``concrete_stresses`` attribute - these are the stresses you see on the ``plot_stress()`` plots. Finally, the ``concrete_forces`` attribute stores the resultant force for each ``AnalysisSection`` and its lever arm to the cracked centroid." ] }, { @@ -266,8 +266,8 @@ "id": "8c3da4fa", "metadata": {}, "source": [ - "### Steel Results\n", - "Steel geometries are stored in the ``steel_geometries`` attribute of the ``StressResult`` object. Similar to the concrete results, the stresses and forces/lever arms are stored, as well as the strains. Here we print the stress for each bar to the terminal." + "### Reinforcement Results\n", + "Lumped reinforcement geometries are stored in the ``lumped_reinforcement_geometries`` attribute of the ``StressResult`` object. Similar to the concrete results, the stresses and forces/lever arms are stored, as well as the strains. Here we print the stress for each lumped bar to the terminal. Note that meshed reinforcement results are also stored, but have not been used for this analysis." ] }, { @@ -285,7 +285,7 @@ "moments_x = []\n", "\n", "# create a Rich table for pretty printing\n", - "table = Table(title=\"Steel Results\")\n", + "table = Table(title=\"Reinforcement Results\")\n", "table.add_column(\"Bar No.\", justify=\"center\", style=\"cyan\", no_wrap=True)\n", "table.add_column(\"Location (x, y) (mm)\", justify=\"center\", style=\"green\")\n", "table.add_column(\"Stress (MPa)\", justify=\"center\", style=\"green\")\n", @@ -293,12 +293,12 @@ "table.add_column(\"Lever Arm (mm)\", justify=\"center\", style=\"green\")\n", "table.add_column(\"Moment (kN.m)\", justify=\"center\", style=\"green\")\n", "\n", - "for idx, steel_geom in enumerate(cracked_stress_res.steel_geometries):\n", - " # get the steel results\n", - " centroid = steel_geom.calculate_centroid()\n", - " stress = cracked_stress_res.steel_stresses[idx]\n", - " strain = cracked_stress_res.steel_strains[idx]\n", - " force, d_x, d_y = cracked_stress_res.steel_forces[idx]\n", + "for idx, reinf_geom in enumerate(cracked_stress_res.lumped_reinforcement_geometries):\n", + " # get the reinforcement results\n", + " centroid = reinf_geom.calculate_centroid()\n", + " stress = cracked_stress_res.lumped_reinforcement_stresses[idx]\n", + " strain = cracked_stress_res.lumped_reinforcement_strains[idx]\n", + " force, d_x, d_y = cracked_stress_res.lumped_reinforcement_forces[idx]\n", "\n", " # calculate the moment each bar creates and store the results\n", " moment_x = force * d_y\n", From 188f35d3496d5d7613ceb22667c0bfbe3bed6d92 Mon Sep 17 00:00:00 2001 From: robbievanleeuwen Date: Wed, 3 Aug 2022 22:42:56 +1000 Subject: [PATCH 4/7] Add composite example; fix reinf stress contours; fix bug in split_geom_at_strains() --- concreteproperties/design_codes.py | 19 +- concreteproperties/pre.py | 10 +- concreteproperties/results.py | 40 +- concreteproperties/tests/test_rotation.py | 20 +- concreteproperties/utils.py | 3 +- docs/source/notebooks/composite_section.ipynb | 527 ++++++++++++++++++ docs/source/notebooks/ultimate_bending.ipynb | 2 +- 7 files changed, 594 insertions(+), 27 deletions(-) create mode 100644 docs/source/notebooks/composite_section.ipynb diff --git a/concreteproperties/design_codes.py b/concreteproperties/design_codes.py index 908296da..21e14fae 100644 --- a/concreteproperties/design_codes.py +++ b/concreteproperties/design_codes.py @@ -234,8 +234,8 @@ def calculate_ultimate_stress( class AS3600(DesignCode): """Design code class for Australian standard AS 3600:2018. - - Note that this design code only supports :class:`~concreteproperties.pre.Concrete` + + Note that this design code only supports :class:`~concreteproperties.pre.Concrete` and :class:`~concreteproperties.pre.SteelBar` material objects. Meshed :class:`~concreteproperties.pre.Steel` material objects are **not** supported as this falls under the composite structures design code. @@ -261,13 +261,20 @@ def assign_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.") + 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: + if ( + abs( + steel_geom.material.stress_strain_profile.get_ultimate_tensile_strain() + ) + < 0.05 + ): self.reinforcement_class = "L" # calculate squash and tensile load @@ -426,7 +433,9 @@ def squash_tensile_load( strain=0.025 ) - force_t = -area * steel_geom.material.stress_strain_profile.get_yield_strength() + force_t = ( + -area * steel_geom.material.stress_strain_profile.get_yield_strength() + ) # add to totals squash_load += force_c diff --git a/concreteproperties/pre.py b/concreteproperties/pre.py index 20146b45..5bba6c1d 100644 --- a/concreteproperties/pre.py +++ b/concreteproperties/pre.py @@ -359,7 +359,7 @@ def add_bar( x: float, y: float, n: int = 4, -) -> Union[Geometry, CompoundGeometry]: +) -> Union[CompoundGeometry]: # type: ignore """Adds a reinforcing bar to a *sectionproperties* geometry. Bars are discretised by four points by default. @@ -392,7 +392,7 @@ def add_bar_rectangular_array( anchor: Tuple[float, float] = (0, 0), exterior_only: bool = False, n: int = 4, -) -> Union[Geometry, CompoundGeometry]: +) -> Union[CompoundGeometry]: # type: ignore """Adds a rectangular array of reinforcing bars to a *sectionproperties* geometry. Bars are discretised by four points by default. @@ -429,7 +429,7 @@ def add_bar_rectangular_array( bar = bar.shift_section(x_offset=x, y_offset=y) geometry = (geometry - bar) + bar - return geometry + return geometry # type: ignore def add_bar_circular_array( @@ -441,7 +441,7 @@ def add_bar_circular_array( theta_0: float = 0, ctr: Tuple[float, float] = (0, 0), n: int = 4, -) -> Union[Geometry, CompoundGeometry]: +) -> Union[CompoundGeometry]: # type: ignore """Adds a circular array of reinforcing bars to a *sectionproperties* geometry. Bars are discretised by four points by default. @@ -469,4 +469,4 @@ def add_bar_circular_array( bar = bar.shift_section(x_offset=x, y_offset=y) geometry = (geometry - bar) + bar - return geometry + return geometry # type: ignore diff --git a/concreteproperties/results.py b/concreteproperties/results.py index b4b2b809..86d7204b 100644 --- a/concreteproperties/results.py +++ b/concreteproperties/results.py @@ -1144,15 +1144,15 @@ def plot_stress( # check region has a force if abs(self.concrete_forces[idx][0]) > 1e-8: # create triangulation - triang = tri.Triangulation( + triang_conc = tri.Triangulation( self.concrete_analysis_sections[idx].mesh_nodes[:, 0], self.concrete_analysis_sections[idx].mesh_nodes[:, 1], self.concrete_analysis_sections[idx].mesh_elements[:, 0:3], # type: ignore ) # plot the filled contour - trictr = fig.axes[0].tricontourf( - triang, sig, v_conc, cmap=cmap_conc, norm=CenteredNorm() + trictr_conc = fig.axes[0].tricontourf( + triang_conc, sig, v_conc, cmap=cmap_conc, norm=CenteredNorm() ) # plot a zero stress contour, supressing warning @@ -1180,23 +1180,29 @@ def plot_stress( zero_level = -1e-12 CS = fig.axes[0].tricontour( - triang, sig, [zero_level], linewidths=1, linestyles="dashed" + triang_conc, + sig, + [zero_level], + linewidths=1, + linestyles="dashed", ) # plot the meshed reinforcement stresses + trictr_reinf = None + for idx, sig in enumerate(self.meshed_reinforcement_stresses): # check region has a force if abs(self.meshed_reinforcement_forces[idx][0]) > 1e-8: # create triangulation - triang = tri.Triangulation( + triang_reinf = tri.Triangulation( self.meshed_reinforcement_sections[idx].mesh_nodes[:, 0], self.meshed_reinforcement_sections[idx].mesh_nodes[:, 1], self.meshed_reinforcement_sections[idx].mesh_elements[:, 0:3], # type: ignore ) # plot the filled contour - trictr = fig.axes[0].tricontourf( - triang, sig, v_reinf, cmap=cmap_reinf, norm=CenteredNorm() + trictr_reinf = fig.axes[0].tricontourf( + triang_reinf, sig, v_reinf, cmap=cmap_reinf, norm=CenteredNorm() ) # plot a zero stress contour, supressing warning @@ -1224,7 +1230,11 @@ def plot_stress( zero_level = -1e-12 CS = fig.axes[0].tricontour( - triang, sig, [zero_level], linewidths=1, linestyles="dashed" + triang_reinf, + sig, + [zero_level], + linewidths=1, + linestyles="dashed", ) # plot the lumped reinforcement stresses @@ -1243,18 +1253,26 @@ def plot_stress( patch.set_array(colours) if reinf_tick_same: patch.set_clim([0.99 * v_reinf[0], 1.01 * v_reinf[-1]]) + else: + patch.set_clim([v_reinf[0], v_reinf[-1]]) fig.axes[0].add_collection(patch) # add the colour bars fig.colorbar( - trictr, # type: ignore + trictr_conc, # type: ignore label="Concrete Stress", format="%.2e", ticks=ticks_conc, cax=fig.axes[1], ) + + if trictr_reinf: + mappable = trictr_reinf + else: + mappable = patch + fig.colorbar( - patch, + mappable, label="Reinforcement Stress", format="%.2e", ticks=ticks_reinf, @@ -1306,7 +1324,7 @@ def sum_moments( moment_sum_x += conc_force[0] * conc_force[2] moment_sum_y += conc_force[0] * conc_force[1] - # sum meshed reinf stresses + # sum meshed reinf stresses for meshed_reinf_force in self.meshed_reinforcement_forces: moment_sum_x += meshed_reinf_force[0] * meshed_reinf_force[2] moment_sum_y += meshed_reinf_force[0] * meshed_reinf_force[1] diff --git a/concreteproperties/tests/test_rotation.py b/concreteproperties/tests/test_rotation.py index 9983af91..81eccf78 100644 --- a/concreteproperties/tests/test_rotation.py +++ b/concreteproperties/tests/test_rotation.py @@ -72,7 +72,10 @@ def test_rotated_gross_properties(theta): pytest.approx(new_gross_results.concrete_area) == ref_gross_results.concrete_area ) - assert pytest.approx(new_gross_results.reinf_lumped_area) == ref_gross_results.reinf_lumped_area + assert ( + pytest.approx(new_gross_results.reinf_lumped_area) + == ref_gross_results.reinf_lumped_area + ) assert pytest.approx(new_gross_results.e_a) == ref_gross_results.e_a assert pytest.approx(new_gross_results.mass) == ref_gross_results.mass assert pytest.approx(new_gross_results.perimeter) == ref_gross_results.perimeter @@ -148,7 +151,10 @@ def test_rotated_uncracked_stress(theta): assert pytest.approx(cf[0]) == ref_uncr_stress.concrete_forces[i][0] for idx, sf in enumerate(new_uncr_stress.lumped_reinforcement_forces): - assert pytest.approx(sf[0]) == ref_uncr_stress.lumped_reinforcement_forces[idx][0] + assert ( + pytest.approx(sf[0]) + == ref_uncr_stress.lumped_reinforcement_forces[idx][0] + ) # list of normal forces @@ -177,7 +183,10 @@ def test_rotated_cracked_stress(theta): assert pytest.approx(cf[0]) == ref_cr_stress.concrete_forces[idx][0] for idx, sf in enumerate(new_cr_stress.lumped_reinforcement_forces): - assert pytest.approx(sf[0]) == ref_cr_stress.lumped_reinforcement_forces[idx][0] + assert ( + pytest.approx(sf[0]) + == ref_cr_stress.lumped_reinforcement_forces[idx][0] + ) # list of normal forces @@ -207,4 +216,7 @@ def test_rotated_ultimate_stress(theta): ) for idx, sf in enumerate(new_ult_stress.lumped_reinforcement_forces): - assert pytest.approx(sf[0], rel=5e-4) == ref_ult_stress.lumped_reinforcement_forces[idx][0] + assert ( + pytest.approx(sf[0], rel=5e-4) + == ref_ult_stress.lumped_reinforcement_forces[idx][0] + ) diff --git a/concreteproperties/utils.py b/concreteproperties/utils.py index 5b77220a..c879c3ae 100644 --- a/concreteproperties/utils.py +++ b/concreteproperties/utils.py @@ -145,6 +145,7 @@ def split_geom_at_strains( # initialise top_geoms in case of two unique strains top_geoms = geom_list + continuing_geoms = [] # loop through intermediate points on stress-strain profile for strain in strains[1:-1]: @@ -185,7 +186,7 @@ def split_geom_at_strains( geom_list = continuing_geoms # save final top geoms - split_geoms.extend(top_geoms) + split_geoms.extend(continuing_geoms) return split_geoms diff --git a/docs/source/notebooks/composite_section.ipynb b/docs/source/notebooks/composite_section.ipynb new file mode 100644 index 00000000..ac2fed36 --- /dev/null +++ b/docs/source/notebooks/composite_section.ipynb @@ -0,0 +1,527 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6efa58a4", + "metadata": {}, + "source": [ + "# Composite Sections" + ] + }, + { + "cell_type": "markdown", + "id": "bb4a755a", + "metadata": {}, + "source": [ + "This example demonstrates how to *concreteproperties* can be used to analyse composite sections. We start by importing the necessary modules.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d12341a2", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from concreteproperties.material import Concrete, Steel, SteelBar\n", + "import concreteproperties.stress_strain_profile as ssp\n", + "import sectionproperties.pre.library.primitive_sections as sp_ps\n", + "import sectionproperties.pre.library.steel_sections as sp_ss\n", + "import sectionproperties.pre.library.concrete_sections as sp_cs\n", + "import concreteproperties.pre as cp_pre\n", + "from concreteproperties.concrete_section import ConcreteSection\n", + "import concreteproperties.results as res" + ] + }, + { + "cell_type": "markdown", + "id": "463c9423", + "metadata": {}, + "source": [ + "## Assign Materials\n", + "The materials used in this example will be 50 MPa concrete with 500 MPa reinforcing steel and 300 MPa structural steel." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a87d1278", + "metadata": {}, + "outputs": [], + "source": [ + "concrete = Concrete(\n", + " name=\"50 MPa Concrete\",\n", + " density=2.4e-6,\n", + " stress_strain_profile=ssp.ConcreteLinearNoTension(\n", + " elastic_modulus=34.8e3,\n", + " ultimate_strain=0.003,\n", + " compressive_strength=0.9 * 50,\n", + " ),\n", + " ultimate_stress_strain_profile=ssp.RectangularStressBlock(\n", + " compressive_strength=50,\n", + " alpha=0.775,\n", + " gamma=0.845,\n", + " ultimate_strain=0.003,\n", + " ),\n", + " alpha_squash=0.85,\n", + " flexural_tensile_strength=4.2,\n", + " colour=\"lightgrey\",\n", + ")\n", + "\n", + "steel_300 = Steel(\n", + " name=\"300 MPa Structural Steel\",\n", + " density=7.85e-6,\n", + " stress_strain_profile=ssp.SteelElasticPlastic(\n", + " yield_strength=300,\n", + " elastic_modulus=200e3,\n", + " fracture_strain=0.05,\n", + " ),\n", + " colour=\"tan\",\n", + ")\n", + "\n", + "steel_bar = SteelBar(\n", + " name=\"500 MPa Steel Bar\",\n", + " density=7.85e-6,\n", + " stress_strain_profile=ssp.SteelElasticPlastic(\n", + " yield_strength=500,\n", + " elastic_modulus=200e3,\n", + " fracture_strain=0.05,\n", + " ),\n", + " colour=\"grey\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "7ee76aca", + "metadata": {}, + "source": [ + "## Square Column with Universal Beam\n", + "First we will analyse a 500 mm x 500 mm square column with 12N20 bars and a central 310UC97 steel column." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7672a60", + "metadata": {}, + "outputs": [], + "source": [ + "# create 500 square concrete\n", + "conc = sp_ps.rectangular_section(d=500, b=500, material=concrete)\n", + "\n", + "# create 310UC97 and centre to column\n", + "uc = sp_ss.i_section(\n", + " d=308,\n", + " b=305,\n", + " t_f=15.4,\n", + " t_w=9.9,\n", + " r=16.5,\n", + " n_r=8,\n", + " material=steel_300,\n", + ").align_center(align_to=conc)\n", + "\n", + "# cut hole in concrete for UC then add UC\n", + "geom = conc - uc + uc\n", + "\n", + "# add 12N20 reinforcing bars\n", + "geom = cp_pre.add_bar_rectangular_array(\n", + " geometry=geom,\n", + " area=310,\n", + " material=steel_bar,\n", + " n_x=4,\n", + " x_s=132,\n", + " n_y=4,\n", + " y_s=132,\n", + " anchor=(52, 52),\n", + " exterior_only=True,\n", + ")\n", + "\n", + "# create concrete section and plot\n", + "conc_sec = ConcreteSection(geom)\n", + "conc_sec.plot_section()" + ] + }, + { + "cell_type": "markdown", + "id": "fd36b69a", + "metadata": {}, + "source": [ + "### Elastic Stress\n", + "Calculate the elastic stress under a bending moment of 100 kN.m." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4ea6715", + "metadata": {}, + "outputs": [], + "source": [ + "el_stress = conc_sec.calculate_uncracked_stress(m_x=100e6)\n", + "el_stress.plot_stress()" + ] + }, + { + "cell_type": "markdown", + "id": "d55583dd", + "metadata": {}, + "source": [ + "### Cracked Stress\n", + "Calculate the cracked stress under a bending moment of 500 kN.m." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91c51681", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "cr_res = conc_sec.calculate_cracked_properties()\n", + "cr_stress = conc_sec.calculate_cracked_stress(cracked_results=cr_res, m=500e6)\n", + "cr_stress.plot_stress()" + ] + }, + { + "cell_type": "markdown", + "id": "a5c1f8bc", + "metadata": {}, + "source": [ + "### Moment Curvature Diagram\n", + "Generate a moment-curvature diagram and plot the stress after yielding of the steel." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0b7e00d", + "metadata": {}, + "outputs": [], + "source": [ + "mk_res = conc_sec.moment_curvature_analysis()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b4dc4669", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "mk_res.plot_results(fmt=\"kx-\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "abd35873", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "serv_stress = conc_sec.calculate_service_stress(\n", + " moment_curvature_results=mk_res, m=None, kappa=2e-5\n", + ")\n", + "serv_stress.plot_stress()" + ] + }, + { + "cell_type": "markdown", + "id": "b03b2bcb", + "metadata": {}, + "source": [ + "### Ultimate Results\n", + "Finally the ultimate bending capacity about the ``x`` and ``y`` axes will be calculated, as well as generating an interaction diagram, a biaxial bending diagram and visualising the ultimate stresses." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5173a385", + "metadata": {}, + "outputs": [], + "source": [ + "ult_res_x = conc_sec.ultimate_bending_capacity()\n", + "ult_res_y = conc_sec.ultimate_bending_capacity(theta=np.pi / 2)\n", + "ult_res_x.print_results()\n", + "ult_res_y.print_results()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f2406b1", + "metadata": {}, + "outputs": [], + "source": [ + "mi_res = conc_sec.moment_interaction_diagram()\n", + "mi_res.plot_diagram()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e735f288", + "metadata": {}, + "outputs": [], + "source": [ + "bb_res = conc_sec.biaxial_bending_diagram()\n", + "bb_res.plot_diagram()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f460ed9", + "metadata": {}, + "outputs": [], + "source": [ + "ult_stress_x = conc_sec.calculate_ultimate_stress(ult_res_x)\n", + "ult_stress_y = conc_sec.calculate_ultimate_stress(ult_res_y)\n", + "ult_stress_x.plot_stress()\n", + "ult_stress_y.plot_stress()" + ] + }, + { + "cell_type": "markdown", + "id": "c49489f9", + "metadata": {}, + "source": [ + "## Concrete Filled Steel Column\n", + "Next we will analyse a 323.9 mm x 12.7 mm steel (350 MPa) circular hollow section filled with concrete and 6N20 bars. The results will be compared to a similarly sized concrete column." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14d9e6e8", + "metadata": {}, + "outputs": [], + "source": [ + "steel_350 = Steel(\n", + " name=\"350 MPa Structural Steel\",\n", + " density=7.85e-6,\n", + " stress_strain_profile=ssp.SteelElasticPlastic(\n", + " yield_strength=350,\n", + " elastic_modulus=200e3,\n", + " fracture_strain=0.05,\n", + " ),\n", + " colour=\"tan\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb581853", + "metadata": {}, + "outputs": [], + "source": [ + "# create outer diameter of steel column\n", + "steel_col = sp_ps.circular_section_by_area(\n", + " area=np.pi * 323.9**2 / 4,\n", + " n=64,\n", + " material=steel_350,\n", + ")\n", + "\n", + "# create inner diameter of steel column, concrete filled\n", + "inner_conc = sp_ps.circular_section_by_area(\n", + " area=np.pi * (323.9 - 2 * 12.7) ** 2 / 4,\n", + " n=64,\n", + " material=concrete,\n", + ")\n", + "\n", + "# create composite geometry\n", + "geom_comp = steel_col - inner_conc + inner_conc\n", + "\n", + "# add reinforcement\n", + "r_bars = 323.9 / 2 - 12.7 - 30 - 10 # 30 mm cover from inside of steel\n", + "\n", + "geom_comp = cp_pre.add_bar_circular_array(\n", + " geometry=geom_comp,\n", + " area=310,\n", + " material=steel_bar,\n", + " n_bar=6,\n", + " r_array=r_bars,\n", + ")\n", + "\n", + "# create concrete section and plot\n", + "conc_sec_comp = ConcreteSection(geom_comp)\n", + "conc_sec_comp.plot_section()\n", + "\n", + "# create 350 diameter column for comparison\n", + "geom_conc = sp_cs.concrete_circular_section(\n", + " d=350,\n", + " n=64,\n", + " dia=20,\n", + " n_bar=6,\n", + " n_circle=4,\n", + " cover=30 + 12, # 30 mm cover + 12 mm tie\n", + " area_conc=np.pi * 350**2 / 4,\n", + " area_bar=310,\n", + " conc_mat=concrete,\n", + " steel_mat=steel_bar,\n", + ")\n", + "\n", + "# create concrete section and plot\n", + "conc_sec_conc = ConcreteSection(geom_conc)\n", + "conc_sec_conc.plot_section()" + ] + }, + { + "cell_type": "markdown", + "id": "d2102f37", + "metadata": {}, + "source": [ + "### Compare Elastic Stresses\n", + "Compare the elastic stresses under a bending moment of 10 kN.m (uncracked)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7d88d9b", + "metadata": {}, + "outputs": [], + "source": [ + "el_stress_comp = conc_sec_comp.calculate_uncracked_stress(m_x=10e6)\n", + "el_stress_conc = conc_sec_conc.calculate_uncracked_stress(m_x=10e6)\n", + "el_stress_comp.plot_stress()\n", + "el_stress_conc.plot_stress()" + ] + }, + { + "cell_type": "markdown", + "id": "f6bd27ac", + "metadata": {}, + "source": [ + "Compare the elastic stresses under a bending moment of 50 kN.m (cracked)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "787f165a", + "metadata": {}, + "outputs": [], + "source": [ + "cr_res_comp = conc_sec_comp.calculate_cracked_properties()\n", + "cr_res_conc = conc_sec_conc.calculate_cracked_properties()\n", + "\n", + "cr_stress_comp = conc_sec_comp.calculate_cracked_stress(\n", + " cracked_results=cr_res_comp, m=50e6\n", + ")\n", + "cr_stress_conc = conc_sec_conc.calculate_cracked_stress(\n", + " cracked_results=cr_res_conc, m=50e6\n", + ")\n", + "cr_stress_comp.plot_stress()\n", + "cr_stress_conc.plot_stress()" + ] + }, + { + "cell_type": "markdown", + "id": "3633d045", + "metadata": {}, + "source": [ + "### Compare Moment Curvature Response\n", + "Calculate the moment curvature response of each column and compare the results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "caecaf08", + "metadata": {}, + "outputs": [], + "source": [ + "mk_res_comp = conc_sec_comp.moment_curvature_analysis()\n", + "mk_res_conc = conc_sec_conc.moment_curvature_analysis()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d7cbc26", + "metadata": {}, + "outputs": [], + "source": [ + "res.MomentCurvatureResults.plot_multiple_results(\n", + " moment_curvature_results=[mk_res_comp, mk_res_conc],\n", + " labels=[\"Composite\", \"Concrete\"],\n", + " fmt=\"-\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "8a02f2f2", + "metadata": {}, + "source": [ + "Although the composite section is significantly stiffer and stronger, the regular concrete section exhibits a more ductile response due to the reduced amount of steel reinforcement." + ] + }, + { + "cell_type": "markdown", + "id": "e867125e", + "metadata": {}, + "source": [ + "### Compare Moment Interaction Results\n", + "Moment interaction diagrams can be produced for each section and compared. Extra points are generated on the tensile side of the composite curve to improve the fidelity of the curve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d2ed0c8", + "metadata": {}, + "outputs": [], + "source": [ + "mi_res_comp = conc_sec_comp.moment_interaction_diagram(\n", + " control_points=[\n", + " (\"kappa0\", 0.0),\n", + " (\"D\", 1.0),\n", + " (\"N\", 0),\n", + " (\"N\", -2e6),\n", + " (\"N\", -3e6),\n", + " (\"d_n\", 1e-06),\n", + " ],\n", + " n_points=[6, 16, 4, 4, 4],\n", + ")\n", + "mi_res_conc = conc_sec_conc.moment_interaction_diagram()\n", + "res.MomentInteractionResults.plot_multiple_diagrams(\n", + " moment_interaction_results=[mi_res_comp, mi_res_conc],\n", + " labels=[\"Composite\", \"Concrete\"],\n", + " fmt=\"x-\",\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/source/notebooks/ultimate_bending.ipynb b/docs/source/notebooks/ultimate_bending.ipynb index 98284947..7a8adb58 100644 --- a/docs/source/notebooks/ultimate_bending.ipynb +++ b/docs/source/notebooks/ultimate_bending.ipynb @@ -37,7 +37,7 @@ "metadata": {}, "source": [ "## Assign Materials\n", - "The materials used in this example will be 50 MPa concrete with 500 MPa, specified in accordance with AS 3600:2018." + "The materials used in this example will be 50 MPa concrete with 500 MPa steel, specified in accordance with AS 3600:2018." ] }, { From 88fee28b72a421c26630f657da77316b4e019ae6 Mon Sep 17 00:00:00 2001 From: robbievanleeuwen Date: Wed, 3 Aug 2022 23:02:57 +1000 Subject: [PATCH 5/7] Fix stress strain profile test --- concreteproperties/tests/test_stress_strain_profile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/concreteproperties/tests/test_stress_strain_profile.py b/concreteproperties/tests/test_stress_strain_profile.py index 521c3d90..864d580b 100644 --- a/concreteproperties/tests/test_stress_strain_profile.py +++ b/concreteproperties/tests/test_stress_strain_profile.py @@ -11,7 +11,7 @@ def test_whitney(): assert pytest.approx(profile.get_stress(-0.001)) == 0 assert pytest.approx(profile.get_stress(-0.003)) == 0 assert pytest.approx(profile.get_stress(-0.1)) == 0 - assert pytest.approx(profile.get_stress(0.00068999)) == 0 + assert pytest.approx(profile.get_stress(0.00068998)) == 0 assert pytest.approx(profile.get_stress(0.00069)) == 0.85 * 40 assert pytest.approx(profile.get_stress(0.001)) == 0.85 * 40 From b75c150d9770b880b999357a0d10e9ba049a339b Mon Sep 17 00:00:00 2001 From: robbievanleeuwen Date: Thu, 4 Aug 2022 23:14:50 +1000 Subject: [PATCH 6/7] Fix ultimate analysis with no lumped reo; update docs --- concreteproperties/concrete_section.py | 8 +- concreteproperties/design_codes.py | 6 +- concreteproperties/material.py | 20 +- concreteproperties/pre.py | 8 +- concreteproperties/results.py | 4 +- docs/source/notebooks/composite_section.ipynb | 177 ++++++++++++++++++ docs/source/rst/analysis.rst | 7 - docs/source/rst/examples.rst | 1 + docs/source/rst/geometry.rst | 35 +++- docs/source/rst/materials.rst | 51 +++-- docs/source/rst/results.rst | 14 +- 11 files changed, 281 insertions(+), 50 deletions(-) diff --git a/concreteproperties/concrete_section.py b/concreteproperties/concrete_section.py index f8f2abf2..ab8a3f2a 100644 --- a/concreteproperties/concrete_section.py +++ b/concreteproperties/concrete_section.py @@ -243,7 +243,7 @@ def get_gross_properties( def get_transformed_gross_properties( self, elastic_modulus: float, - ) -> res.TransformedConcreteProperties: + ) -> res.TransformedGrossProperties: """Transforms gross section properties given a reference elastic modulus. :param elastic_modulus: Reference elastic modulus @@ -251,7 +251,7 @@ def get_transformed_gross_properties( :return: Transformed concrete properties object """ - return res.TransformedConcreteProperties( + return res.TransformedGrossProperties( concrete_properties=self.gross_properties, elastic_modulus=elastic_modulus ) @@ -927,12 +927,14 @@ def calculate_ultimate_section_actions( # save results ultimate_results.d_n = d_n - ultimate_results.k_u = min(k_u) 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 def moment_interaction_diagram( diff --git a/concreteproperties/design_codes.py b/concreteproperties/design_codes.py index 21e14fae..74dd5c0d 100644 --- a/concreteproperties/design_codes.py +++ b/concreteproperties/design_codes.py @@ -62,7 +62,7 @@ def create_steel_material( yield_strength: float, colour: str = "grey", ) -> SteelBar: - """Returns a steel material object. + """Returns a steel bar material object. List assumptions of material properties here... @@ -91,7 +91,7 @@ def get_gross_properties( def get_transformed_gross_properties( self, **kwargs, - ) -> res.TransformedConcreteProperties: + ) -> res.TransformedGrossProperties: """Transforms gross section properties. :param kwargs: Keyword arguments passed to @@ -362,7 +362,7 @@ def create_steel_material( ductility_class: str = "N", colour: str = "grey", ) -> SteelBar: - r"""Returns a steel material object. + r"""Returns a steel bar material object. | **Material assumptions:** | - *Density*: 7850 kg/m\ :sup:`3` diff --git a/concreteproperties/material.py b/concreteproperties/material.py index 12a16cc2..86cfa9d2 100644 --- a/concreteproperties/material.py +++ b/concreteproperties/material.py @@ -8,14 +8,14 @@ @dataclass class Material: - """Abstract class for a *concreteproperties* material. + """Generic class for a *concreteproperties* material. :param name: Material name :param density: Material density (mass per unit volume) :param stress_strain_profile: Material stress-strain profile :param colour: Colour of the material for rendering - :param meshed: If set to True, the material region is meshed; if set to False, the - material region is treated as a lumped circular mass at its centroid + :param meshed: If set to True, the entire material region is meshed; if set to + False, the material region is treated as a lumped circular mass at its centroid """ name: str @@ -71,7 +71,8 @@ def __post_init__(self): @dataclass class Steel(Material): - """Class for a steel material. + """Class for a steel material with the entire region meshed to allow for strain + variation across the section, e.g. structural steel profiles. :param name: Steel material name :param density: Steel density (mass per unit volume) @@ -88,7 +89,8 @@ class Steel(Material): @dataclass class SteelBar(Steel): - """Class for a steel material. + """Class for a steel bar material, treated as a lumped circular mass with a constant + strain. :param name: Steel bar material name :param density: Steel bar density (mass per unit volume) @@ -103,7 +105,7 @@ class SteelBar(Steel): meshed: bool = field(default=False, init=False) -@dataclass -class SteelStrand(Steel): - # placeholder - pass +# @dataclass +# class SteelStrand(Steel): +# # placeholder +# pass diff --git a/concreteproperties/pre.py b/concreteproperties/pre.py index 5bba6c1d..ba500962 100644 --- a/concreteproperties/pre.py +++ b/concreteproperties/pre.py @@ -16,7 +16,7 @@ from sectionproperties.pre.geometry import CompoundGeometry - from concreteproperties.material import Material, Steel + from concreteproperties.material import Material, SteelBar class CPGeom: @@ -355,7 +355,7 @@ def __init__( def add_bar( geometry: Union[Geometry, CompoundGeometry], area: float, - material: Steel, + material: SteelBar, x: float, y: float, n: int = 4, @@ -384,7 +384,7 @@ def add_bar( def add_bar_rectangular_array( geometry: Union[Geometry, CompoundGeometry], area: float, - material: Steel, + material: SteelBar, n_x: int, x_s: float, n_y: int = 1, @@ -435,7 +435,7 @@ def add_bar_rectangular_array( def add_bar_circular_array( geometry: Union[Geometry, CompoundGeometry], area: float, - material: Steel, + material: SteelBar, n_bar: int, r_array: float, theta_0: float = 0, diff --git a/concreteproperties/results.py b/concreteproperties/results.py index 86d7204b..0ecda2e2 100644 --- a/concreteproperties/results.py +++ b/concreteproperties/results.py @@ -143,7 +143,7 @@ def print_results( @dataclass -class TransformedConcreteProperties: +class TransformedGrossProperties: """Class for storing transformed gross concrete section properties. :param concrete_properties: Concrete properties object @@ -597,7 +597,7 @@ def print_results( table.add_row("Bending Angle - theta", "{:>{fmt}}".format(self.theta, fmt=fmt)) table.add_row("Neutral Axis Depth - d_n", "{:>{fmt}}".format(self.d_n, fmt=fmt)) table.add_row( - "Neutral Axis Parameter- k_u", "{:>{fmt}}".format(self.k_u, fmt=fmt) + "Neutral Axis Parameter - k_u", "{:>{fmt}}".format(self.k_u, fmt=fmt) ) table.add_row("Axial Force", "{:>{fmt}}".format(self.n, fmt=fmt)) table.add_row("Bending Capacity - m_x", "{:>{fmt}}".format(self.m_x, fmt=fmt)) diff --git a/docs/source/notebooks/composite_section.ipynb b/docs/source/notebooks/composite_section.ipynb index ac2fed36..03940622 100644 --- a/docs/source/notebooks/composite_section.ipynb +++ b/docs/source/notebooks/composite_section.ipynb @@ -501,6 +501,183 @@ " fmt=\"x-\",\n", ")" ] + }, + { + "cell_type": "markdown", + "id": "31c952e3", + "metadata": {}, + "source": [ + "## Composite Steel Beam\n", + "Finally we will analyse a steel composite beam, consisting of a 530UB92 universal beam compositely attached to a 1500 mm wide x 120 mm deep concrete slab. We start by defining the concrete material (32 MPa) for the slab." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "65577277", + "metadata": {}, + "outputs": [], + "source": [ + "concrete = Concrete(\n", + " name=\"32 MPa Concrete\",\n", + " density=2.4e-6,\n", + " stress_strain_profile=ssp.ConcreteLinearNoTension(\n", + " elastic_modulus=30.1e3,\n", + " ultimate_strain=0.003,\n", + " compressive_strength=0.9 * 32,\n", + " ),\n", + " ultimate_stress_strain_profile=ssp.RectangularStressBlock(\n", + " compressive_strength=32,\n", + " alpha=0.802,\n", + " gamma=0.89,\n", + " ultimate_strain=0.003,\n", + " ),\n", + " alpha_squash=0.85,\n", + " flexural_tensile_strength=3.4,\n", + " colour=\"lightgrey\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "112ad018", + "metadata": {}, + "source": [ + "Next we construct the composite geometry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d08a46da", + "metadata": {}, + "outputs": [], + "source": [ + "# create 530UB92\n", + "ub = sp_ss.i_section(\n", + " d=533,\n", + " b=209,\n", + " t_f=15.6,\n", + " t_w=10.2,\n", + " r=14,\n", + " n_r=8,\n", + " material=steel_300,\n", + ")\n", + "\n", + "# create concrete slab, centre on top of UB\n", + "conc_slab = (\n", + " sp_ps.rectangular_section(\n", + " d=120,\n", + " b=1500,\n", + " material=concrete,\n", + " )\n", + " .align_center(align_to=ub)\n", + " .align_to(other=ub, on=\"top\")\n", + ")\n", + "\n", + "# as there is no overlapping geometry, we can simply add the two sections together\n", + "geom = ub + conc_slab\n", + "\n", + "# the cross-section is rotated 90 degrees to aid viewing the stress plots\n", + "geom = geom.rotate_section(angle=np.pi / 2, use_radians=True)\n", + "\n", + "# create concrete section and plot\n", + "conc_sec = ConcreteSection(geom)\n", + "conc_sec.plot_section()" + ] + }, + { + "cell_type": "markdown", + "id": "77745bac", + "metadata": {}, + "source": [ + "### Moment Curvature Response\n", + "Calculate the moment curvature response of the composite section. The section should show a gradual yielding response as plasticity develops over the depth of the steel beam." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fcbc6333", + "metadata": {}, + "outputs": [], + "source": [ + "mk_res = conc_sec.moment_curvature_analysis(theta=np.pi / 2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f634623", + "metadata": {}, + "outputs": [], + "source": [ + "mk_res.plot_results()" + ] + }, + { + "cell_type": "markdown", + "id": "5890c6c4", + "metadata": {}, + "source": [ + "We can examine the stresses at various point within the moment-curvature response." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b64edebc", + "metadata": {}, + "outputs": [], + "source": [ + "el_stress = conc_sec.calculate_service_stress(moment_curvature_results=mk_res, m=500e6)\n", + "yield_stress = conc_sec.calculate_service_stress(\n", + " moment_curvature_results=mk_res, m=1000e6\n", + ")\n", + "ult_stress = conc_sec.calculate_service_stress(\n", + " moment_curvature_results=mk_res, m=1200e6\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d56e8c7d", + "metadata": {}, + "outputs": [], + "source": [ + "el_stress.plot_stress()\n", + "yield_stress.plot_stress()\n", + "ult_stress.plot_stress()" + ] + }, + { + "cell_type": "markdown", + "id": "b9695c03", + "metadata": {}, + "source": [ + "### Ultimate Bending Capacity\n", + "The ultimate bending capacity can be calculated and compared to the result obtained from the moment curvature analysis." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1c8cd3d", + "metadata": {}, + "outputs": [], + "source": [ + "ult_res = conc_sec.ultimate_bending_capacity(theta=np.pi / 2)\n", + "ult_res.print_results()" + ] + }, + { + "cell_type": "markdown", + "id": "f95ae617", + "metadata": {}, + "source": [ + "The above results match closely with the moment curvature results, as will as confirming that the neutral axis is within the concrete slab. Note that ``k_u`` is not calculated for only meshed sections, it is only computed for lumped reinforcement." + ] } ], "metadata": { diff --git a/docs/source/rst/analysis.rst b/docs/source/rst/analysis.rst index 2baca9c9..b7ca7ac5 100644 --- a/docs/source/rst/analysis.rst +++ b/docs/source/rst/analysis.rst @@ -104,13 +104,6 @@ extreme compressive fibre to the ``ultimate_strain`` parameter (see :ref:`label-conc-ult-profile`) and finding the neutral axis which satisfies the equilibrium of axial forces. -.. note:: - - The bending neutral axis must lie within the cross-section, if a combination of - ``theta`` and ``n`` is provided to the method that results in the neutral axis angle - lying outside the cross-section, the algorithm will not converge and a warning message - displayed. - .. seealso:: For an application of the above, see the example :ref:`/notebooks/ultimate_bending.ipynb`. diff --git a/docs/source/rst/examples.rst b/docs/source/rst/examples.rst index 9f4b403e..d4e141a0 100644 --- a/docs/source/rst/examples.rst +++ b/docs/source/rst/examples.rst @@ -14,4 +14,5 @@ The following examples complement the structured documentation presented thus fa ../notebooks/moment_interaction ../notebooks/biaxial_bending ../notebooks/stress_analysis + ../notebooks/composite_section ../notebooks/design_codes diff --git a/docs/source/rst/geometry.rst b/docs/source/rst/geometry.rst index 01df5d2a..252436d4 100644 --- a/docs/source/rst/geometry.rst +++ b/docs/source/rst/geometry.rst @@ -5,7 +5,7 @@ Geometry concrete cross-sections as :class:`~sectionproperties.pre.geometry.CompoundGeometry` objects. *concreteproperties* expects a :class:`~sectionproperties.pre.geometry.CompoundGeometry` object to contain geometry -information describing both the concrete and steel components of the reinforced conrete +information describing both the concrete and reinforcement components of the cross-section. The concrete portion of the cross-section can be constructed in many ways: @@ -21,6 +21,11 @@ Next, reinforcing bars can be added to the concrete geometry using the methods i :mod:`~concreteproperties.pre` module, and are also included in geometries generated by the Concrete Sections Library. +Alternatively, other sections (such as structural steel sections) can be combined with +concrete geometry to create composite sections. This can be done using geometry creation +methods similar to that of concrete, except that a different material object will be +provided. + Below is a brief overview of the various methods that can be used to create geometries to be used by *concreteproperties*. A more exhaustive overview can be found in the *sectionproperties* |sp_docs_link|. @@ -29,9 +34,9 @@ in the *sectionproperties* |sp_docs_link|. Axis Conventions ---------------- -Reinforced concrete cross-sections in *concreteproperties* are constructed in the +Concrete cross-sections in *concreteproperties* are constructed in the ``x-y`` plane. A key feature of *concreteproperties* is that the analyses are not -restricted to bending about the ``x`` or ``y`` axis. All analysis types (were relevant) +restricted to bending about the ``x`` or ``y`` axis. All analysis types (where relevant) give the user the option of performing the analysis about a rotated ``u-v`` axis, defined by the angle :math:`\theta` (in radians). @@ -306,6 +311,30 @@ geometry. Note how a hole is automatically generated in the closed-off region: geom = beam + slab geom.plot_geometry() +The following code creates a 323.9 mm x 12.7 mm circular hollow section filled with +concrete and demonstrates how composite sections can be created: + +.. plot:: + :include-source: True + :caption: Concrete filled circular hollow section + + import sectionproperties.pre.library.primitive_sections as sp_ps + + conc = None # define your concrete material properties here + steel = None # define your steel material properties here + steel_col = sp_ps.circular_section_by_area( + area=np.pi * 323.9**2 / 4, + n=64, + material=steel, + ) + inner_conc = sp_ps.circular_section_by_area( + area=np.pi * (323.9 - 2 * 12.7) ** 2 / 4, + n=64, + material=conc, + ) + geom = steel_col - inner_conc + inner_conc + geom.plot_geometry() + Adding Reinforcing Bars ----------------------- diff --git a/docs/source/rst/materials.rst b/docs/source/rst/materials.rst index 993e30f9..1322ba40 100644 --- a/docs/source/rst/materials.rst +++ b/docs/source/rst/materials.rst @@ -7,33 +7,61 @@ properties can be used for a single cross-section. For example, higher strength sections can be topped with lower grade in-situ slabs, and high tensile steel can be used in combination with normal grade reinforcing steel. -The structural behaviour of both concrete and steel materials are described by -:ref:`stress-strain-profiles`. +The structural behaviour of materials is described by :ref:`stress-strain-profiles`. .. note:: In *concreteproperties*, a positive sign is given to compressive forces, stresses and strains, while a negative sign is given to tensile forces, stresses and strains. + +Material Classes +---------------- + +*concreteproperties* ships with material objects describing the structural behaviour of +both concrete and steel. The generic :class:`~concreteproperties.material.Material` +class can be used to describe the behaviour of any other material. + +By default, all geometries in *concreteproperties* are meshed to capture strain +variation across the section. However, for smaller geometries (such as reinforcement), +*concreteproperties* can treat the area as having a constant strain with a lumped mass, +which increases the performance of the analysis with almost no loss in fidelity. The +meshing can be switched off by setting the attribute ``meshed=False``. + +The :class:`~concreteproperties.material.SteelBar` class has meshing disabled by default +and should be used when defining steel reinforcement. On the other hand, the +:class:`~concreteproperties.material.Steel` class is meshed by default so should be used +when defining larger sections such as strucutral steel sections used in composite +sections. + + +Material +^^^^^^^^ + +.. autoclass:: concreteproperties.material.Material + :noindex: + + Concrete --------- +^^^^^^^^ .. autoclass:: concreteproperties.material.Concrete :noindex: Steel ------ +^^^^^ .. autoclass:: concreteproperties.material.Steel :noindex: -.. note:: - In the current version of *concreteproperties*, all geometries assigned with the - :class:`~concreteproperties.material.Steel` material are treated as lumped masses. - In *concreteproperties* ``v0.3.0+``, the Materials module will be overhauled allowing - for steel geometries to be treated as structural sections. +SteelBar +^^^^^^^^ + +.. autoclass:: concreteproperties.material.SteelBar + :noindex: + .. _stress-strain-profiles: @@ -45,9 +73,8 @@ service and ultimate analyses. A :class:`~concreteproperties.material.Concrete` requires both a **service** stress-strain profile (calculation of area properties, moment-curvature analysis, elastic and service stress analysis) and an **ultimate** stress-strain profile (ultimate bending capacity, moment interaction diagram, biaxial -bending diagram, ultimate stress analysis). A -:class:`~concreteproperties.material.Steel` object only requires one stress-strain -profile which is used for both service and ultimate analyses. +bending diagram, ultimate stress analysis). All other material objects only requires one +stress-strain profile which is used for both service and ultimate analyses. .. note:: diff --git a/docs/source/rst/results.rst b/docs/source/rst/results.rst index c969c17e..3287009a 100644 --- a/docs/source/rst/results.rst +++ b/docs/source/rst/results.rst @@ -18,12 +18,12 @@ method. .. automethod:: concreteproperties.concrete_section.ConcreteSection.get_gross_properties :noindex: -This method returns a :class:`~concreteproperties.results.ConcreteProperties` object, +This method returns a :class:`~concreteproperties.results.GrossProperties` object, which stores all the calculated section properties as attributes. The gross area properties can be printed to the terminal by calling the -:meth:`~concreteproperties.results.ConcreteProperties.print_results` method. +:meth:`~concreteproperties.results.GrossProperties.print_results` method. -.. autoclass:: concreteproperties.results.ConcreteProperties() +.. autoclass:: concreteproperties.results.GrossProperties() :noindex: :members: @@ -34,13 +34,13 @@ method. .. automethod:: concreteproperties.concrete_section.ConcreteSection.get_transformed_gross_properties :noindex: -This method returns a :class:`~concreteproperties.results.TransformedConcreteProperties` +This method returns a :class:`~concreteproperties.results.TransformedGrossProperties` object, which stores all the calculated transformed section properties as class attributes. The transformed gross area properties can be printed to the terminal by calling the -:meth:`~concreteproperties.results.TransformedConcreteProperties.print_results` method. +:meth:`~concreteproperties.results.TransformedGrossProperties.print_results` method. -.. autoclass:: concreteproperties.results.TransformedConcreteProperties() +.. autoclass:: concreteproperties.results.TransformedGrossProperties() :noindex: :members: @@ -61,7 +61,7 @@ returns a :class:`~concreteproperties.results.CrackedResults` object. :members: Calling -:meth:`~concreteproperties.results.TransformedConcreteProperties.calculate_transformed_properties` +:meth:`~concreteproperties.results.TransformedGrossProperties.calculate_transformed_properties` on a :class:`~concreteproperties.results.CrackedResults` object stores the transformed cracked properties as attributes within the current object. From 41f59366927036ee06891810781744f5ddd72380 Mon Sep 17 00:00:00 2001 From: robbievanleeuwen Date: Thu, 4 Aug 2022 23:27:03 +1000 Subject: [PATCH 7/7] Clean up imports --- concreteproperties/analysis_section.py | 2 +- concreteproperties/concrete_section.py | 12 ++++++------ concreteproperties/material.py | 1 - concreteproperties/post.py | 4 ---- concreteproperties/pre.py | 7 +++---- concreteproperties/stress_strain_profile.py | 2 +- concreteproperties/tests/test_moment_interaction.py | 1 - concreteproperties/utils.py | 3 ++- 8 files changed, 13 insertions(+), 19 deletions(-) diff --git a/concreteproperties/analysis_section.py b/concreteproperties/analysis_section.py index dd0bb20f..66560f26 100644 --- a/concreteproperties/analysis_section.py +++ b/concreteproperties/analysis_section.py @@ -9,8 +9,8 @@ from matplotlib.colors import ListedColormap import concreteproperties.utils as utils -from concreteproperties.post import plotting_context from concreteproperties.material import Concrete +from concreteproperties.post import plotting_context if TYPE_CHECKING: import matplotlib diff --git a/concreteproperties/concrete_section.py b/concreteproperties/concrete_section.py index ab8a3f2a..d2b2ff5f 100644 --- a/concreteproperties/concrete_section.py +++ b/concreteproperties/concrete_section.py @@ -1,7 +1,7 @@ from __future__ import annotations import warnings -from math import inf, nan, isinf +from math import inf, isinf, nan from typing import TYPE_CHECKING, List, Optional, Tuple, Union import matplotlib.patches as mpatches @@ -14,7 +14,7 @@ import concreteproperties.results as res import concreteproperties.utils as utils from concreteproperties.analysis_section import AnalysisSection -from concreteproperties.material import Concrete, Steel +from concreteproperties.material import Concrete from concreteproperties.post import plotting_context from concreteproperties.pre import CPGeom, CPGeomConcrete @@ -1035,7 +1035,7 @@ def decode_d_n(cp): f"Provided d_n {cp[1]:.3f} must be greater than zero." ) return cp[1] - # extreme tensile steel yield ratio + # extreme tensile reinforcement yield ratio elif cp[0] == "fy": # get extreme tensile bar d_ext, eps_sy = self.extreme_bar(theta=theta) @@ -1254,8 +1254,8 @@ def calculate_uncracked_stress( """Calculates stresses within the reinforced concrete section assuming an uncracked section. - Uses gross area section properties to determine concrete and steel stresses - given an axial force ``n``, and bending moments ``m_x`` and ``m_y``. + 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 @@ -1752,7 +1752,7 @@ def calculate_ultimate_stress( # get position of lump centroid = lumped_geom.calculate_centroid() - # get strain at centroid of steel + # get strain at centroid of lump if isinf(ultimate_results.d_n): strain = self.gross_properties.conc_ultimate_strain else: diff --git a/concreteproperties/material.py b/concreteproperties/material.py index 86cfa9d2..03acdd03 100644 --- a/concreteproperties/material.py +++ b/concreteproperties/material.py @@ -1,7 +1,6 @@ from __future__ import annotations from dataclasses import dataclass, field -from typing import TYPE_CHECKING import concreteproperties.stress_strain_profile as ssp diff --git a/concreteproperties/post.py b/concreteproperties/post.py index 5360bd40..6056cdb9 100644 --- a/concreteproperties/post.py +++ b/concreteproperties/post.py @@ -4,14 +4,10 @@ from typing import TYPE_CHECKING, Optional, Tuple, Union import matplotlib.pyplot as plt -import numpy as np if TYPE_CHECKING: import matplotlib - from concreteproperties.analysis_section import AnalysisSection - from concreteproperties.concrete_section import ConcreteSection - @contextlib.contextmanager def plotting_context( diff --git a/concreteproperties/pre.py b/concreteproperties/pre.py index ba500962..2f72a3b7 100644 --- a/concreteproperties/pre.py +++ b/concreteproperties/pre.py @@ -1,19 +1,18 @@ from __future__ import annotations from typing import TYPE_CHECKING, List, Tuple, Union -from more_itertools import peekable import numpy as np -from shapely.geometry import Polygon, LineString -from shapely.ops import split import sectionproperties.pre.library.primitive_sections as sp_ps +from more_itertools import peekable from sectionproperties.pre.geometry import Geometry +from shapely.geometry import LineString, Polygon +from shapely.ops import split from concreteproperties.material import Concrete if TYPE_CHECKING: import matplotlib - from sectionproperties.pre.geometry import CompoundGeometry from concreteproperties.material import Material, SteelBar diff --git a/concreteproperties/stress_strain_profile.py b/concreteproperties/stress_strain_profile.py index 40d2730d..7388633f 100644 --- a/concreteproperties/stress_strain_profile.py +++ b/concreteproperties/stress_strain_profile.py @@ -2,7 +2,7 @@ import warnings from dataclasses import dataclass, field -from typing import TYPE_CHECKING, List, Optional, Union +from typing import TYPE_CHECKING, List, Union import matplotlib.pyplot as plt import numpy as np diff --git a/concreteproperties/tests/test_moment_interaction.py b/concreteproperties/tests/test_moment_interaction.py index 0a3751d8..f4d0ecdb 100644 --- a/concreteproperties/tests/test_moment_interaction.py +++ b/concreteproperties/tests/test_moment_interaction.py @@ -1,5 +1,4 @@ import numpy as np -import pytest from concreteproperties.concrete_section import ConcreteSection from concreteproperties.material import Concrete, SteelBar from concreteproperties.stress_strain_profile import ( diff --git a/concreteproperties/utils.py b/concreteproperties/utils.py index c879c3ae..29f904ff 100644 --- a/concreteproperties/utils.py +++ b/concreteproperties/utils.py @@ -6,11 +6,12 @@ from rich.progress import BarColumn, Progress, ProgressColumn, SpinnerColumn, TextColumn from rich.table import Column from rich.text import Text -from sectionproperties.pre.geometry import CompoundGeometry from concreteproperties.pre import CPGeomConcrete if TYPE_CHECKING: + from sectionproperties.pre.geometry import CompoundGeometry + from concreteproperties.pre import CPGeom