from typing import List, Union, Optional, Tuple
import copy
from dataclasses import dataclass, asdict
import warnings
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import matplotlib.cm as cm
import matplotlib.patches as mpatches
from matplotlib.colors import ListedColormap, CenteredNorm
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
from rich.console import Console
from rich.text import Text
from rich.live import Live
from rich.table import Table
from rich.panel import Panel
from rich.progress import Progress
import numpy as np
from scipy.sparse import csc_matrix, coo_matrix, linalg
from scipy.optimize import brentq
import sectionproperties.pre.pre as pre
import sectionproperties.pre.geometry as section_geometry
import sectionproperties.analysis.fea as fea
import sectionproperties.analysis.solver as solver
import sectionproperties.post.post as post
from shapely import Polygon, Point
from shapely.strtree import STRtree
[docs]class Section:
"""Class for structural cross-sections.
Stores the finite element geometry, mesh and material information and provides methods to
compute the cross-section properties. The element type used in this program is the six-noded
quadratic triangular element.
The constructor extracts information from the provided mesh object and creates and stores the
corresponding Tri6 finite element objects.
:param geometry: Cross-section geometry object used to generate the mesh
:type geometry: :class:`~sectionproperties.pre.geometry.Geometry`
:param bool time_info: If set to True, a detailed description of the computation and the time
cost is printed to the terminal for every computation performed.
The following example creates a :class:`~sectionproperties.analysis.section.Section`
object of a 100D x 50W rectangle using a mesh size of 5::
import sectionproperties.pre.library.primitive_sections as primitive_sections
from sectionproperties.analysis.section import Section
geometry = primitive_sections.rectangular_section(d=100, b=50)
geometry.create_mesh(mesh_sizes=[5])
section = Section(geometry)
:cvar elements: List of finite element objects describing the cross-section mesh
:vartype elements: list[:class:`~sectionproperties.analysis.fea.Tri6`]
:cvar int num_nodes: Number of nodes in the finite element mesh
:cvar geometry: Cross-section geometry object used to generate the mesh
:vartype geometry: :class:`~sectionproperties.pre.geometry.Geometry`
:cvar mesh: Mesh dict returned by triangle
:vartype mesh: dict(mesh)
:cvar mesh_nodes: Array of node coordinates from the mesh
:vartype mesh_nodes: :class:`numpy.ndarray`
:cvar mesh_elements: Array of connectivities from the mesh
:vartype mesh_elements: :class:`numpy.ndarray`
:cvar mesh_attributes: Array of attributes from the mesh
:vartype mesh_attributes: :class:`numpy.ndarray`
:cvar materials: List of materials
:type materials: list[:class:`~sectionproperties.pre.pre.Material`]
:cvar material_groups: List of objects containing the elements in each defined material
:type material_groups: list[:class:`~sectionproperties.pre.pre.MaterialGroup`]
:cvar section_props: Class to store calculated section properties
:vartype section_props: :class:`~sectionproperties.analysis.section.SectionProperties`
:raises AssertionError: If the number of materials does not equal the number of regions
:raises ValueError: If geometry does not contain a mesh
"""
def __init__(
self,
geometry: Union[section_geometry.Geometry, section_geometry.CompoundGeometry],
time_info: bool = False,
):
"""Inits the Section class."""
if not hasattr(geometry, "mesh") or not geometry.mesh:
raise ValueError(
"Selected Geometry or CompoundGeometry "
"object does not contain a mesh.\n"
"Try running {geometry}.create_mesh() before adding to "
"a Section object for analysis."
)
self.geometry = geometry
self.time_info = time_info
self.mesh = geometry.mesh
self.materials = []
mesh = self.mesh
if isinstance(self.geometry, section_geometry.CompoundGeometry):
self.materials = [geom.material for geom in self.geometry.geoms]
else:
self.materials = [self.geometry.material]
# extract mesh data
nodes = np.array(mesh["vertices"], dtype=np.dtype(float))
elements = np.array(mesh["triangles"], dtype=np.dtype(int))
attributes = np.array(mesh["triangle_attributes"].T[0], dtype=np.dtype(int))
# swap mid-node order to retain node ordering consistency
elements[:, [3, 4, 5]] = elements[:, [5, 3, 4]]
# save total number of nodes in mesh
self.num_nodes = len(nodes)
# initialise material_sections variable
self.material_groups = []
# if materials are specified, check that the right number of material properties are
# specified and then populate material_groups list
if self.materials:
msg = "Number of materials ({0}), ".format(len(self.materials))
msg += "should match the number of regions ({0}).".format(
max(attributes) + 1
)
assert len(self.materials) == max(attributes) + 1, msg
# add a MaterialGroup object to the material_groups list for each uniquely
# encountered material
for (i, material) in enumerate(self.materials):
# add the first material to the list
if i == 0:
self.material_groups.append(MaterialGroup(material, self.num_nodes))
else:
# if the material hasn't been encountered
if material not in self.materials[:i]:
self.material_groups.append(
MaterialGroup(material, self.num_nodes)
)
self.elements = [] # initialise list holding all element objects
# build the mesh one element at a time
for (i, node_ids) in enumerate(elements):
x1 = nodes[node_ids[0]][0]
y1 = nodes[node_ids[0]][1]
x2 = nodes[node_ids[1]][0]
y2 = nodes[node_ids[1]][1]
x3 = nodes[node_ids[2]][0]
y3 = nodes[node_ids[2]][1]
x4 = nodes[node_ids[3]][0]
y4 = nodes[node_ids[3]][1]
x5 = nodes[node_ids[4]][0]
y5 = nodes[node_ids[4]][1]
x6 = nodes[node_ids[5]][0]
y6 = nodes[node_ids[5]][1]
# create a list containing the vertex and mid-node coordinates
coords = np.array([[x1, x2, x3, x4, x5, x6], [y1, y2, y3, y4, y5, y6]])
# if materials are specified, get the material
if self.materials:
# get attribute index of current element
att_el = attributes[i]
# fetch the material
material = self.materials[att_el]
# if there are no materials specified, use a default material
else: # Should not happen but included as failsafe
material = pre.DEFAULT_MATERIAL
# add tri6 elements to the mesh
new_element = fea.Tri6(i, coords, node_ids, material)
self.elements.append(new_element)
# add element to relevant MaterialGroup
for group in self.material_groups:
if material is group.material:
group.add_element(new_element)
break
# save mesh input
self.mesh = mesh
self.mesh_nodes = nodes
self.mesh_elements = elements
self.mesh_attributes = attributes
# create the search tree
p_mesh = [
Polygon(self.geometry.mesh["vertices"][tri][0:3])
for tri in self.geometry.mesh["triangles"]
]
self.mesh_search_tree = STRtree(p_mesh)
# initialise class storing section properties
self.section_props = SectionProperties()
[docs] def calculate_geometric_properties(self):
"""Calculates the geometric properties of the cross-section and stores them in the
:class:`~sectionproperties.analysis.section.SectionProperties` object contained in
the ``section_props`` class variable.
The following geometric section properties are calculated:
* Cross-sectional area
* Cross-sectional perimeter
* Cross-sectional mass
* Area weighted material properties, composite only :math:`E_{eff}`, :math:`G_{eff}`, :math:`{nu}_{eff}`
* Modulus weighted area (axial rigidity)
* First moments of area
* Second moments of area about the global axis
* Second moments of area about the centroidal axis
* Elastic centroid
* Centroidal section moduli
* Radii of gyration
* Principal axis properties
If materials are specified for the cross-section, the moments of area and section moduli
are elastic modulus weighted.
The following example demonstrates the use of this method::
section = Section(geometry)
section.calculate_geometric_properties()
"""
def calculate_geom(progress=None):
# initialise properties
self.section_props.area = 0
self.section_props.perimeter = 0
self.section_props.mass = 0
self.section_props.ea = 0
self.section_props.ga = 0
self.section_props.qx = 0
self.section_props.qy = 0
self.section_props.ixx_g = 0
self.section_props.iyy_g = 0
self.section_props.ixy_g = 0
# calculate perimeter
self.section_props.perimeter = self.geometry.calculate_perimeter()
if progress is not None:
task = progress.add_task(
description="[red]Calculating geometric properties",
total=len(self.elements),
)
# calculate global geometric properties
for el in self.elements:
(
area,
qx,
qy,
ixx_g,
iyy_g,
ixy_g,
e,
g,
rho,
) = el.geometric_properties()
self.section_props.area += area
self.section_props.mass += area * rho
self.section_props.ea += area * e
self.section_props.ga += area * g
self.section_props.qx += qx * e
self.section_props.qy += qy * e
self.section_props.ixx_g += ixx_g * e
self.section_props.iyy_g += iyy_g * e
self.section_props.ixy_g += ixy_g * e
if progress is not None:
progress.update(task, advance=1)
self.section_props.nu_eff = (
self.section_props.ea / (2 * self.section_props.ga) - 1
)
self.section_props.e_eff = self.section_props.ea / self.section_props.area
self.section_props.g_eff = self.section_props.ga / self.section_props.area
self.section_props.calculate_elastic_centroid()
self.section_props.calculate_centroidal_properties(self.mesh)
if progress is not None:
progress.update(
task,
description="[bold green]:white_check_mark: Geometric analysis complete",
)
# conduct geometric analysis
if self.time_info:
progress = solver.create_progress()
progress_table = Table.grid()
panel = Panel.fit(
progress,
title="[bold]Geometric Analysis",
border_style="red",
padding=(1, 1),
)
progress_table.add_row(panel)
with Live(progress_table, refresh_per_second=10) as live:
calculate_geom(progress=progress)
panel.border_style = "green"
live.refresh()
print()
else:
calculate_geom()
[docs] def calculate_warping_properties(self, solver_type="direct"):
"""Calculates all the warping properties of the cross-section and stores them in the
:class:`~sectionproperties.analysis.section.SectionProperties` object contained in
the ``section_props`` class variable.
:param string solver_type: Solver used for solving systems of linear equations, either
using the *'direct'* method or *'cgs'* iterative method
The following warping section properties are calculated:
* Torsion constant
* Shear centre
* Shear area
* Warping constant
* Monosymmetry constant
If materials are specified, the values calculated for the torsion constant, warping
constant and shear area are elastic modulus weighted.
Note that the geometric properties must be calculated prior to the calculation of the
warping properties::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
:raises RuntimeError: If the geometric properties have not been
calculated prior to calling this method
"""
# check that a geometric analysis has been performed
if None in [
self.section_props.area,
self.section_props.ixx_c,
self.section_props.cx,
]:
err = "Calculate geometric properties before performing a warping analysis."
raise RuntimeError(err)
# check if any geometry are disjoint
if isinstance(self.geometry, section_geometry.CompoundGeometry):
polygons = [sec_geom.geom for sec_geom in self.geometry.geoms]
disjoint_regions = section_geometry.check_geometry_disjoint(polygons)
if disjoint_regions:
warnings.warn(
"\nThe section geometry contains disjoint regions which is invalid for warping analysis.\n"
"Please revise your geometry to ensure there is connectivity between all regions.\n"
"Please see https://sectionproperties.readthedocs.io/en/latest/rst/analysis.html#warping-analysis for more "
"information."
)
# create a new Section with the origin shifted to the centroid for calculation of the
# warping properties such that the Lagrangian multiplier approach can be utilised
warping_section = Section(self.geometry)
# shift the coordinates of each element N.B. the mesh class attribute remains unshifted!
for el in warping_section.elements:
el.coords[0, :] -= self.section_props.cx
el.coords[1, :] -= self.section_props.cy
# entire warping analysis is contained in following definition
def warping_analysis(progress=None):
# assemble stiffness matrix and load vector for warping function
if progress:
task = progress.add_task(
description="[red]Assembling {0}x{0} stiffness matrix".format(
self.num_nodes
),
total=len(warping_section.elements),
)
(k_lg, f_torsion) = warping_section.assemble_torsion(
progress=progress, task=task
)
progress.update(
task,
description="[green]:white_check_mark: {0}x{0} stiffness matrix assembled".format(
self.num_nodes
),
)
progress.update(0, advance=1)
else:
(k_lg, f_torsion) = warping_section.assemble_torsion()
# ILU decomposition of stiffness matrices
def ilu_decomp(progress=None, task=None):
# ILU decomposition on regular stiffness matrix
k_precond = linalg.LinearOperator(
(self.num_nodes, self.num_nodes), linalg.spilu(k).solve
)
if progress is not None:
progress.update(task, advance=1)
# ILU decomposition on Lagrangian stiffness matrix
k_lg_precond = linalg.LinearOperator(
(self.num_nodes + 1, self.num_nodes + 1), linalg.spilu(k_lg).solve
)
if progress is not None:
progress.update(task, advance=1)
return (k_precond, k_lg_precond)
# if the cgs method is used, perform ILU decomposition
if solver_type == "cgs":
if progress:
task = progress.add_task(
description="[red]Performing ILU decomposition",
total=2,
)
(k_precond, k_lg_precond) = ilu_decomp(progress=progress, task=task)
progress.update(
task,
description="[green]:white_check_mark: ILU decomposition complete",
)
else:
(k_precond, k_lg_precond) = ilu_decomp()
# solve for warping function
def solve_warping():
if solver_type == "cgs":
omega = solver.solve_cgs_lagrange(k_lg, f_torsion, m=k_lg_precond)
elif solver_type == "direct":
omega = solver.solve_direct_lagrange(k_lg, f_torsion)
return omega
if progress:
task = progress.add_task(
description="[red]Solving warping function ({0})".format(
solver_type
),
total=1,
)
omega = solve_warping()
progress.update(task, advance=1)
progress.update(
task,
description="[green]:white_check_mark: Warping function solved ({0})".format(
solver_type
),
)
progress.update(0, advance=1)
else:
omega = solve_warping()
# save the warping function
self.section_props.omega = omega
# determine the torsion constant
self.section_props.j = (
self.section_props.ixx_c
+ self.section_props.iyy_c
- omega.dot(k_lg[:-1, :-1].dot(np.transpose(omega)))
)
# assemble shear function load vectors
def assemble_shear_load(progress=None, task=None):
f_psi = np.zeros(self.num_nodes)
f_phi = np.zeros(self.num_nodes)
for el in warping_section.elements:
(f_psi_el, f_phi_el) = el.shear_load_vectors(
self.section_props.ixx_c,
self.section_props.iyy_c,
self.section_props.ixy_c,
self.section_props.nu_eff,
)
f_psi[el.node_ids] += f_psi_el
f_phi[el.node_ids] += f_phi_el
if progress is not None:
progress.update(task, advance=1)
return (f_psi, f_phi)
if progress:
task = progress.add_task(
description="[red]Assembling shear function vectors",
total=len(warping_section.elements),
)
(f_psi, f_phi) = assemble_shear_load(progress=progress, task=task)
progress.update(
task,
description="[green]:white_check_mark: Shear function vectors assembled",
)
progress.update(0, advance=1)
else:
(f_psi, f_phi) = assemble_shear_load()
# solve for shear functions psi and phi
def solve_shear_functions(progress=None, task=None):
if solver_type == "cgs":
psi_shear = solver.solve_cgs_lagrange(k_lg, f_psi, m=k_lg_precond)
if progress is not None:
progress.update(task, advance=1)
phi_shear = solver.solve_cgs_lagrange(k_lg, f_phi, m=k_lg_precond)
if progress is not None:
progress.update(task, advance=1)
elif solver_type == "direct":
psi_shear = solver.solve_direct_lagrange(k_lg, f_psi)
if progress is not None:
progress.update(task, advance=1)
phi_shear = solver.solve_direct_lagrange(k_lg, f_phi)
if progress is not None:
progress.update(task, advance=1)
return (psi_shear, phi_shear)
if progress:
task = progress.add_task(
description="[red]Solving shear functions ({0})".format(
solver_type
),
total=2,
)
(psi_shear, phi_shear) = solve_shear_functions(
progress=progress, task=task
)
progress.update(
task,
description="[green]:white_check_mark: Shear functions solved ({0})".format(
solver_type
),
)
progress.update(0, advance=1)
else:
(psi_shear, phi_shear) = solve_shear_functions()
# save the shear functions
self.section_props.psi_shear = psi_shear
self.section_props.phi_shear = phi_shear
# assemble shear centre and warping moment integrals
def assemble_sc_warping_integrals(progress=None, task=None):
sc_xint = 0
sc_yint = 0
q_omega = 0
i_omega = 0
i_xomega = 0
i_yomega = 0
for el in warping_section.elements:
(
sc_xint_el,
sc_yint_el,
q_omega_el,
i_omega_el,
i_xomega_el,
i_yomega_el,
) = el.shear_warping_integrals(
self.section_props.ixx_c,
self.section_props.iyy_c,
self.section_props.ixy_c,
omega[el.node_ids],
)
sc_xint += sc_xint_el
sc_yint += sc_yint_el
q_omega += q_omega_el
i_omega += i_omega_el
i_xomega += i_xomega_el
i_yomega += i_yomega_el
if progress is not None:
progress.update(task, advance=1)
return (sc_xint, sc_yint, q_omega, i_omega, i_xomega, i_yomega)
if progress:
task = progress.add_task(
description="[red]Assembling shear and warping integrals",
total=len(warping_section.elements),
)
(
sc_xint,
sc_yint,
q_omega,
i_omega,
i_xomega,
i_yomega,
) = assemble_sc_warping_integrals(progress=progress, task=task)
progress.update(
task,
description="[green]:white_check_mark: Shear and warping integrals assembled",
)
progress.update(0, advance=1)
else:
(
sc_xint,
sc_yint,
q_omega,
i_omega,
i_xomega,
i_yomega,
) = assemble_sc_warping_integrals()
# calculate shear centres (elasticity approach)
Delta_s = (
2
* (1 + self.section_props.nu_eff)
* (
self.section_props.ixx_c * self.section_props.iyy_c
- self.section_props.ixy_c**2
)
)
x_se = (1 / Delta_s) * (
(self.section_props.nu_eff / 2 * sc_xint) - f_torsion.dot(phi_shear)
)
y_se = (1 / Delta_s) * (
(self.section_props.nu_eff / 2 * sc_yint) + f_torsion.dot(psi_shear)
)
(x11_se, y22_se) = fea.principal_coordinate(
self.section_props.phi, x_se, y_se
)
# calculate shear centres (Trefftz's approach)
x_st = (
self.section_props.ixy_c * i_xomega
- self.section_props.iyy_c * i_yomega
) / (
self.section_props.ixx_c * self.section_props.iyy_c
- self.section_props.ixy_c**2
)
y_st = (
self.section_props.ixx_c * i_xomega
- self.section_props.ixy_c * i_yomega
) / (
self.section_props.ixx_c * self.section_props.iyy_c
- self.section_props.ixy_c**2
)
# save shear centres
self.section_props.Delta_s = Delta_s
self.section_props.x_se = x_se
self.section_props.y_se = y_se
self.section_props.x11_se = x11_se
self.section_props.y22_se = y22_se
self.section_props.x_st = x_st
self.section_props.y_st = y_st
# calculate warping constant
self.section_props.gamma = (
i_omega
- q_omega**2 / self.section_props.ea
- y_se * i_xomega
+ x_se * i_yomega
)
def assemble_shear_deformation(progress=None, task=None):
# assemble shear deformation coefficients
kappa_x = 0
kappa_y = 0
kappa_xy = 0
for el in warping_section.elements:
(kappa_x_el, kappa_y_el, kappa_xy_el) = el.shear_coefficients(
self.section_props.ixx_c,
self.section_props.iyy_c,
self.section_props.ixy_c,
psi_shear[el.node_ids],
phi_shear[el.node_ids],
self.section_props.nu_eff,
)
kappa_x += kappa_x_el
kappa_y += kappa_y_el
kappa_xy += kappa_xy_el
if progress is not None:
progress.update(task, advance=1)
return (kappa_x, kappa_y, kappa_xy)
if progress:
task = progress.add_task(
description="[red]Assembling shear deformation coefficients",
total=len(warping_section.elements),
)
(kappa_x, kappa_y, kappa_xy) = assemble_shear_deformation(
progress=progress, task=task
)
progress.update(
task,
description="[green]:white_check_mark: Shear deformation coefficients assembled",
)
progress.update(0, advance=1)
else:
(kappa_x, kappa_y, kappa_xy) = assemble_shear_deformation()
# calculate shear areas wrt global axis
self.section_props.A_sx = Delta_s**2 / kappa_x
self.section_props.A_sy = Delta_s**2 / kappa_y
self.section_props.A_sxy = Delta_s**2 / kappa_xy
# calculate shear areas wrt principal bending axis:
alpha_xx = kappa_x * self.section_props.area / Delta_s**2
alpha_yy = kappa_y * self.section_props.area / Delta_s**2
alpha_xy = kappa_xy * self.section_props.area / Delta_s**2
# rotate the tensor by the principal axis angle
phi_rad = self.section_props.phi * np.pi / 180
R = np.array(
[
[np.cos(phi_rad), np.sin(phi_rad)],
[-np.sin(phi_rad), np.cos(phi_rad)],
]
)
rotatedAlpha = R.dot(
np.array([[alpha_xx, alpha_xy], [alpha_xy, alpha_yy]])
).dot(np.transpose(R))
# recalculate the shear area based on the rotated alpha value
self.section_props.A_s11 = self.section_props.area / rotatedAlpha[0, 0]
self.section_props.A_s22 = self.section_props.area / rotatedAlpha[1, 1]
# calculate the monosymmetry consants
def calculate_monosymmetry_integrals(progress=None, task=None):
int_x = 0
int_y = 0
int_11 = 0
int_22 = 0
for el in warping_section.elements:
(
int_x_el,
int_y_el,
int_11_el,
int_22_el,
) = el.monosymmetry_integrals(self.section_props.phi)
int_x += int_x_el
int_y += int_y_el
int_11 += int_11_el
int_22 += int_22_el
if progress is not None:
progress.update(task, advance=1)
return (int_x, int_y, int_11, int_22)
if progress:
task = progress.add_task(
description="[red]Assembling monosymmetry integrals",
total=len(self.elements),
)
(int_x, int_y, int_11, int_22) = calculate_monosymmetry_integrals(
progress=progress, task=task
)
progress.update(
task,
description="[green]:white_check_mark: Monosymmetry integrals assembled",
)
progress.update(0, advance=1)
else:
(int_x, int_y, int_11, int_22) = calculate_monosymmetry_integrals()
# calculate the monosymmetry constants
self.section_props.beta_x_plus = (
-int_x / self.section_props.ixx_c + 2 * self.section_props.y_se
)
self.section_props.beta_x_minus = (
int_x / self.section_props.ixx_c - 2 * self.section_props.y_se
)
self.section_props.beta_y_plus = (
-int_y / self.section_props.iyy_c + 2 * self.section_props.x_se
)
self.section_props.beta_y_minus = (
int_y / self.section_props.iyy_c - 2 * self.section_props.x_se
)
self.section_props.beta_11_plus = (
-int_11 / self.section_props.i11_c + 2 * self.section_props.y22_se
)
self.section_props.beta_11_minus = (
int_11 / self.section_props.i11_c - 2 * self.section_props.y22_se
)
self.section_props.beta_22_plus = (
-int_22 / self.section_props.i22_c + 2 * self.section_props.x11_se
)
self.section_props.beta_22_minus = (
int_22 / self.section_props.i22_c - 2 * self.section_props.x11_se
)
# conduct warping analysis
if self.time_info:
# create warping progress
progress = solver.create_progress()
total_task = progress.add_task(
description="[bold red]Conducting warping analysis...", total=7
)
progress_table = Table.grid()
panel = Panel.fit(
progress,
title="[bold]Warping Analysis",
border_style="red",
padding=(1, 1),
)
progress_table.add_row(panel)
with Live(progress_table, refresh_per_second=10) as live:
warping_analysis(progress)
progress.update(
total_task,
description="[bold green]:white_check_mark: Warping analysis completed",
)
panel.border_style = "green"
live.refresh()
print()
else:
warping_analysis()
[docs] def calculate_frame_properties(self, solver_type="direct"):
"""Calculates and returns the properties required for a frame analysis. The properties are
also stored in the :class:`~sectionproperties.analysis.section.SectionProperties`
object contained in the ``section_props`` class variable.
:param string solver_type: Solver used for solving systems of linear equations, either
using the *'direct'* method or *'cgs'* iterative method
:return: Cross-section properties to be used for a frame analysis *(area, ixx, iyy, ixy, j,
phi)*
:rtype: tuple(float, float, float, float, float, float)
The following section properties are calculated:
* Cross-sectional area *(area)*
* Second moments of area about the centroidal axis *(ixx, iyy, ixy)*
* Torsion constant *(j)*
* Principal axis angle *(phi)*
If materials are specified for the cross-section, the area, second moments of area and
torsion constant are elastic modulus weighted.
The following example demonstrates the use of this method::
section = Section(geometry)
(area, ixx, iyy, ixy, j, phi) = section.calculate_frame_properties()
"""
# initialise geometric properties
self.section_props.area = 0
self.section_props.ea = 0
self.section_props.qx = 0
self.section_props.qy = 0
self.section_props.ixx_g = 0
self.section_props.iyy_g = 0
self.section_props.ixy_g = 0
self.section_props.ixx_c = 0
self.section_props.iyy_c = 0
self.section_props.ixy_c = 0
self.section_props.j = 0
self.section_props.phi = 0
# calculate global geometric properties
for el in self.elements:
(area, qx, qy, ixx_g, iyy_g, ixy_g, e, _, _) = el.geometric_properties()
self.section_props.area += area
self.section_props.ea += area * e
self.section_props.qx += qx * e
self.section_props.qy += qy * e
self.section_props.ixx_g += ixx_g * e
self.section_props.iyy_g += iyy_g * e
self.section_props.ixy_g += ixy_g * e
# calculate elastic centroid location
self.section_props.calculate_elastic_centroid()
# calculate second moments of area about the centroidal xy axis
self.section_props.ixx_c = (
self.section_props.ixx_g
- self.section_props.qx**2 / self.section_props.ea
)
self.section_props.iyy_c = (
self.section_props.iyy_g
- self.section_props.qy**2 / self.section_props.ea
)
self.section_props.ixy_c = (
self.section_props.ixy_g
- self.section_props.qx * self.section_props.qy / self.section_props.ea
)
# calculate the principal axis angle
Delta = (
((self.section_props.ixx_c - self.section_props.iyy_c) / 2) ** 2
+ self.section_props.ixy_c**2
) ** 0.5
i11_c = (self.section_props.ixx_c + self.section_props.iyy_c) / 2 + Delta
# calculate initial principal axis angle
if abs(self.section_props.ixx_c - i11_c) < 1e-12 * i11_c:
self.section_props.phi = 0
else:
self.section_props.phi = (
np.arctan2(self.section_props.ixx_c - i11_c, self.section_props.ixy_c)
* 180
/ np.pi
)
# create a new Section with the origin shifted to the centroid for calculation of
# the warping properties
warping_section = Section(self.geometry)
# shift the coordinates of each element N.B. the mesh class attribute remains unshifted
for el in warping_section.elements:
el.coords[0, :] -= self.section_props.cx
el.coords[1, :] -= self.section_props.cy
(k_lg, f) = warping_section.assemble_torsion()
# if the cgs method is used, perform ILU decomposition
if solver_type == "cgs":
k_lg_precond = linalg.LinearOperator(
(self.num_nodes + 1, self.num_nodes + 1), linalg.spilu(k_lg).solve
)
# solve for warping function
if solver_type == "cgs":
omega = solver.solve_cgs_lagrange(k_lg, f, m=k_lg_precond)
elif solver_type == "direct":
omega = solver.solve_direct_lagrange(k_lg, f)
# calculate the torsion constant
self.section_props.j = (
self.section_props.ixx_c
+ self.section_props.iyy_c
- omega.dot(k_lg[:-1, :-1].dot(np.transpose(omega)))
)
return (
self.section_props.ea,
self.section_props.ixx_c,
self.section_props.iyy_c,
self.section_props.ixy_c,
self.section_props.j,
self.section_props.phi,
)
[docs] def calculate_plastic_properties(self, verbose=False):
"""Calculates the plastic properties of the cross-section and stores them in the
:class:`~sectionproperties.analysis.section.SectionProperties` object contained in
the ``section_props`` class variable.
:param bool verbose: If set to True, the number of iterations required for each plastic
axis is printed to the terminal.
The following warping section properties are calculated:
* Plastic centroid for bending about the centroidal and principal axes
* Plastic section moduli for bending about the centroidal and principal axes
* Shape factors for bending about the centroidal and principal axes
If materials are specified for the cross-section, the plastic section moduli are displayed
as plastic moments (i.e :math:`M_p = f_y S`) and the shape factors are not calculated.
Note that the geometric properties must be calculated before the plastic properties are
calculated::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_plastic_properties()
:raises RuntimeError: If the geometric properties have not been calculated prior to calling
this method
"""
# check that a geometric analysis has been performed
if self.section_props.cx is None:
err = "Calculate geometric properties before performing a plastic analysis."
raise RuntimeError(err)
# check if any geometry are overlapped
if isinstance(self.geometry, section_geometry.CompoundGeometry):
polygons = [sec_geom.geom for sec_geom in self.geometry.geoms]
overlapped_regions = section_geometry.check_geometry_overlaps(polygons)
if overlapped_regions:
warnings.warn(
"\nThe section geometry contains overlapping regions and the area of these overlapped regions "
"will be 'double counted' in the plastic analysis which may result in incorrect values.\n"
"If you do not intend for this double counting to occur, use a subtractive modelling approach "
"to remove the overlapping region.\n"
"Please see https://sectionproperties.readthedocs.io/en/latest/rst/advanced_geom.html for more "
"information."
)
def calc_plastic(progress=None):
plastic_section = PlasticSection(self.geometry)
plastic_section.calculate_plastic_properties(self, verbose, progress)
if self.time_info:
progress = solver.create_progress()
progress_table = Table.grid()
panel = Panel.fit(
progress,
title="[bold]Plastic Analysis",
border_style="red",
padding=(1, 1),
)
progress_table.add_row(panel)
with Live(progress_table, refresh_per_second=10) as live:
calc_plastic(progress=progress)
panel.border_style = "green"
live.refresh()
print()
else:
calc_plastic()
[docs] def calculate_stress(self, N=0, Vx=0, Vy=0, Mxx=0, Myy=0, M11=0, M22=0, Mzz=0):
"""Calculates the cross-section stress resulting from design actions and returns
a :class:`~sectionproperties.analysis.section.StressPost` object allowing
post-processing of the stress results.
:param float N: Axial force
:param float Vx: Shear force acting in the x-direction
:param float Vy: Shear force acting in the y-direction
:param float Mxx: Bending moment about the centroidal xx-axis
:param float Myy: Bending moment about the centroidal yy-axis
:param float M11: Bending moment about the centroidal 11-axis
:param float M22: Bending moment about the centroidal 22-axis
:param float Mzz: Torsion moment about the centroidal zz-axis
:return: Object for post-processing cross-section stresses
:rtype: :class:`~sectionproperties.analysis.section.StressPost`
Note that a geometric analysis must be performed prior to performing a stress
analysis. Further, if the shear force or torsion is non-zero a warping analysis
must also be performed::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(N=1e3, Vy=3e3, Mxx=1e6)
:raises RuntimeError: If a geometric and warping analysis (if required) have not
been performed prior to calling this method
"""
# check that a geometric and warping analysis has been performed
if None in [
self.section_props.area,
self.section_props.ixx_c,
self.section_props.cx,
]:
err = "Perform a geometric analysis before carrying out a stress analysis."
raise RuntimeError(err)
if self.section_props.omega is None and (Vx != 0 or Vy != 0 or Mzz != 0):
err = "Perform a warping analysis before carrying out a stress analysis "
err += "with non-zero shear forces or torsion moment."
raise RuntimeError(err)
def calc_stress(progress=None):
if progress is not None:
task = progress.add_task(
description="[red]Calculating cross-section stresses",
total=len(self.elements),
)
# create stress post object
stress_post = StressPost(self)
# get relevant section properties
ea = self.section_props.ea
cx = self.section_props.cx
cy = self.section_props.cy
ixx = self.section_props.ixx_c
iyy = self.section_props.iyy_c
ixy = self.section_props.ixy_c
i11 = self.section_props.i11_c
i22 = self.section_props.i22_c
phi = self.section_props.phi
j = self.section_props.j
Delta_s = self.section_props.Delta_s
nu = self.section_props.nu_eff
# loop through all material groups
for group in stress_post.material_groups:
# allocate nodal weights vector for nodal averaging
nodal_weights = np.zeros(self.num_nodes)
# loop through all elements in the material group
for el in group.elements:
# get element omega and psi
if self.section_props.omega is None:
omega_el = None
psi_shear_el = None
phi_shear_el = None
else:
omega_el = self.section_props.omega[el.node_ids]
psi_shear_el = self.section_props.psi_shear[el.node_ids]
phi_shear_el = self.section_props.phi_shear[el.node_ids]
(
sig_zz_n_el,
sig_zz_mxx_el,
sig_zz_myy_el,
sig_zz_m11_el,
sig_zz_m22_el,
sig_zx_mzz_el,
sig_zy_mzz_el,
sig_zx_vx_el,
sig_zy_vx_el,
sig_zx_vy_el,
sig_zy_vy_el,
weights,
) = el.element_stress(
N,
Mxx,
Myy,
M11,
M22,
Mzz,
Vx,
Vy,
ea,
cx,
cy,
ixx,
iyy,
ixy,
i11,
i22,
phi,
j,
nu,
omega_el,
psi_shear_el,
phi_shear_el,
Delta_s,
)
# add stresses to global vectors
group.stress_result.sig_zz_n[el.node_ids] += sig_zz_n_el * weights
group.stress_result.sig_zz_mxx[el.node_ids] += (
sig_zz_mxx_el * weights
)
group.stress_result.sig_zz_myy[el.node_ids] += (
sig_zz_myy_el * weights
)
group.stress_result.sig_zz_m11[el.node_ids] += (
sig_zz_m11_el * weights
)
group.stress_result.sig_zz_m22[el.node_ids] += (
sig_zz_m22_el * weights
)
group.stress_result.sig_zx_mzz[el.node_ids] += (
sig_zx_mzz_el * weights
)
group.stress_result.sig_zy_mzz[el.node_ids] += (
sig_zy_mzz_el * weights
)
group.stress_result.sig_zx_vx[el.node_ids] += sig_zx_vx_el * weights
group.stress_result.sig_zy_vx[el.node_ids] += sig_zy_vx_el * weights
group.stress_result.sig_zx_vy[el.node_ids] += sig_zx_vy_el * weights
group.stress_result.sig_zy_vy[el.node_ids] += sig_zy_vy_el * weights
# add nodal weights
nodal_weights[el.node_ids] += weights
if progress is not None:
progress.update(task, advance=1)
# nodal averaging
for (i, weight) in enumerate(nodal_weights):
if weight != 0:
group.stress_result.sig_zz_n[i] *= 1 / weight
group.stress_result.sig_zz_mxx[i] *= 1 / weight
group.stress_result.sig_zz_myy[i] *= 1 / weight
group.stress_result.sig_zz_m11[i] *= 1 / weight
group.stress_result.sig_zz_m22[i] *= 1 / weight
group.stress_result.sig_zx_mzz[i] *= 1 / weight
group.stress_result.sig_zy_mzz[i] *= 1 / weight
group.stress_result.sig_zx_vx[i] *= 1 / weight
group.stress_result.sig_zy_vx[i] *= 1 / weight
group.stress_result.sig_zx_vy[i] *= 1 / weight
group.stress_result.sig_zy_vy[i] *= 1 / weight
# calculate combined stresses
group.stress_result.calculate_combined_stresses()
if progress is not None:
progress.update(
task,
description="[bold green]:white_check_mark: Stress analysis complete",
)
return stress_post
if self.time_info:
progress = solver.create_progress()
progress_table = Table.grid()
panel = Panel.fit(
progress,
title="[bold]Stress Analysis",
border_style="red",
padding=(1, 1),
)
progress_table.add_row(panel)
with Live(progress_table, refresh_per_second=10) as live:
stress_post = calc_stress(progress=progress)
panel.border_style = "green"
live.refresh()
print()
else:
stress_post = calc_stress()
# return the stress_post object
return stress_post
[docs] def assemble_torsion(self, progress=None, task=None):
"""Assembles stiffness matrices to be used for the computation of warping
properties and the torsion load vector (f_torsion). A Lagrangian multiplier
(k_lg) stiffness matrix is returned. The stiffness matrix are assembled using
the sparse COO format and returned in the sparse CSC format.
:return: Lagrangian multiplier stiffness matrix and torsion load
vector *(k_lg, f_torsion)*
:rtype: tuple(:class:`scipy.sparse.csc_matrix`, :class:`numpy.ndarray`)
"""
# initialise variables
N = self.num_nodes # size of matrix
row = [] # list holding row indices
col = [] # list holding column indices
data = [] # list holding stiffness matrix entries
f_torsion = np.zeros(N) # force vector array
# loop through all elements in the mesh
for el in self.elements:
# determine number of nodes in the current element
n = len(el.node_ids)
# calculate the element stiffness matrix and torsion load vector
(k_el, f_el) = el.torsion_properties()
# assemble the torsion load vector
f_torsion[el.node_ids] += f_el
# create row index vector
r = np.repeat(el.node_ids, n)
# create column index vector
c = np.tile(el.node_ids, n)
# flatten element stiffness matrix
k = k_el.flatten()
# add to global arrays
row = np.hstack((row, r))
col = np.hstack((col, c))
data = np.hstack((data, k))
if progress is not None:
progress.update(task, advance=1)
# construct Lagrangian multiplier matrix:
# column vector of ones
row = np.hstack((row, range(N)))
col = np.hstack((col, np.repeat(N, N)))
data = np.hstack((data, np.repeat(1, N)))
# row vector of ones
row = np.hstack((row, np.repeat(N, N)))
col = np.hstack((col, range(N)))
data = np.hstack((data, np.repeat(1, N)))
# zero in bottom right corner
row = np.hstack((row, N))
col = np.hstack((col, N))
data = np.hstack((data, 0))
k_lg = coo_matrix((data, (row, col)), shape=(N + 1, N + 1))
return (csc_matrix(k_lg), f_torsion)
[docs] def plot_mesh(
self,
alpha=0.5,
materials=True,
mask=None,
title="Finite Element Mesh",
**kwargs,
):
r"""Plots the finite element mesh.
:param float alpha: Transparency of the mesh outlines: :math:`0 \leq \alpha \leq 1`
:param bool materials: If set to true shades the elements with the specified material
colours
:param mask: Mask array, of length ``num_nodes``, to mask out triangles
:type mask: list[bool]
:param string title: Plot title
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the mesh generated for the second example
listed under the :class:`~sectionproperties.analysis.section.Section` object
definition::
import sectionproperties.pre.library.primitive_sections as primitive_sections
from sectionproperties.pre.pre import Material
from sectionproperties.analysis.section import Section
steel = Material(
name='Steel', elastic_modulus=200e3, poissons_ratio=0.3, density=7.85e-6,
yield_strength=250, color='grey'
)
timber = Material(
name='Timber', elastic_modulus=8e3, poissons_ratio=0.35, density=6.5e-7,
yield_strength=20, color='burlywood'
)
geom_steel = primitive_sections.rectangular_section(d=50, b=50, material=steel)
geom_timber = primitive_sections.rectangular_section(d=50, b=50, material=timber)
geometry = geom_timber.align_to(geom_steel, on="right") + geom_steel
geometry.create_mesh(mesh_sizes=[10, 5])
section = Section(geometry)
section.plot_mesh(materials=True, alpha=0.5)
.. figure:: ../images/composite_mesh.png
:align: center
:scale: 75 %
Finite element mesh generated by the above example.
"""
# create plot and setup the plot
with post.plotting_context(title=title, **kwargs) as (fig, ax):
# if the material colours are to be displayed
if materials and self.materials is not None:
color_array = []
legend_labels = []
c = [] # Indices of elements for mapping colors
# create an array of finite element colours
for idx, element in enumerate(self.elements):
color_array.append(element.material.color)
c.append(idx)
# create a list of unique material legend entries
for (i, material) in enumerate(self.materials):
# if the material has not be entered yet
if i == 0 or material not in self.materials[0:i]:
# add the material colour and name to the legend list
patch = mpatches.Patch(
color=material.color, label=material.name
)
legend_labels.append(patch)
cmap = ListedColormap(color_array) # custom colormap
# plot the mesh colours
ax.tripcolor(
self.mesh_nodes[:, 0],
self.mesh_nodes[:, 1],
self.mesh_elements[:, 0:3],
c,
cmap=cmap,
)
# display the legend
ax.legend(
loc="center left", bbox_to_anchor=(1, 0.5), handles=legend_labels
)
# plot the mesh
ax.triplot(
self.mesh_nodes[:, 0],
self.mesh_nodes[:, 1],
self.mesh_elements[:, 0:3],
lw=0.5,
color="black",
alpha=alpha,
mask=mask,
)
return ax
[docs] def plot_centroids(self, title="Centroids", **kwargs):
"""Plots the elastic centroid, the shear centre, the plastic centroids and the principal
axis, if they have been calculated, on top of the finite element mesh.
:param string title: Plot title
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example analyses a 200 PFC section and displays a plot of
the centroids::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.channel_section(d=200, b=75, t_f=12, t_w=6, r=12, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
section.calculate_plastic_properties()
section.plot_centroids()
.. figure:: ../images/pfc_centroids.png
:align: center
:scale: 75 %
Plot of the centroids generated by the above example.
The following example analyses a 150x90x12 UA section and displays a plot of the
centroids::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
section.calculate_plastic_properties()
section.plot_centroids()
.. figure:: ../images/angle_centroids.png
:align: center
:scale: 75 %
Plot of the centroids generated by the above example.
"""
# create plot and setup the plot
with post.plotting_context(title=title, **kwargs) as (fig, ax):
# plot the finite element mesh
self.plot_mesh(**dict(kwargs, ax=ax))
# if the elastic centroid has been calculated
if self.section_props.cx is not None:
ax.scatter(
self.section_props.cx,
self.section_props.cy,
edgecolors="r",
facecolors="none",
marker="o",
s=100,
label="Elastic centroid",
)
# if the shear centre has been calculated
if self.section_props.x_se is not None:
(x_s, y_s) = self.get_sc()
ax.scatter(x_s, y_s, c="r", marker="+", s=100, label="Shear centre")
# if the global plastic centroid has been calculated
if self.section_props.x_pc is not None:
(x_pc, y_pc) = self.get_pc()
ax.scatter(
x_pc,
y_pc,
c="r",
marker="x",
s=100,
label="Global plastic centroid",
)
# if the principal plastic centroid has been calculated
if self.section_props.x11_pc is not None:
(x11_pc, y22_pc) = self.get_pc_p()
ax.scatter(
x11_pc,
y22_pc,
edgecolors="r",
facecolors="none",
marker="s",
s=100,
label="Principal plastic centroid",
)
# if the principal axis has been calculated
if self.section_props.phi is not None:
post.draw_principal_axis(
ax,
self.section_props.phi * np.pi / 180,
self.section_props.cx,
self.section_props.cy,
)
# display the legend
ax.legend(loc="center left", bbox_to_anchor=(1, 0.5))
return ax
[docs] def display_mesh_info(self):
"""Prints mesh statistics (number of nodes, elements and regions) to the command window.
The following example displays the mesh statistics for a Tee section merged from two
rectangles::
import sectionproperties.pre.library.primitive_sections as primitive_sections
from sectionproperties.analysis.section import Section
rec1 = primitive_sections.rectangular_section(d=100, b=25)
rec2 = primitive_sections.rectangular_section(d=25, b=100)
rec1 = rec1.shift_section(x_offset=-12.5)
rec2 = rec2.shift_section(x_offset=-50, y_offset=100)
geometry = rec1 + rec2
geometry.create_mesh(mesh_sizes=[5, 2.5])
section = Section(geometry)
section.display_mesh_info()
>>>Mesh Statistics:
>>>--4920 nodes
>>>--2365 elements
>>>--2 regions
"""
console = Console()
console.print("Mesh Statistics:", style="bold underline magenta")
console.print(f"- {self.num_nodes} nodes")
console.print(f"- {len(self.elements)} elements")
regions = max(self.mesh_attributes) + 1
text = f"- {regions} region"
if regions > 1:
text += "s"
console.print(text)
console.print()
[docs] def display_results(self, fmt="8.6e"):
"""Prints the results that have been calculated to the terminal.
:param string fmt: Number formatting string
The following example displays the geometric section properties for a 100D x 50W rectangle
with three digits after the decimal point::
import sectionproperties.pre.library.primitive_sections as primitive_sections
from sectionproperties.analysis.section import Section
geometry = primitive_sections.rectangular_section(d=100, b=50)
geometry.create_mesh(mesh_sizes=[5])
section = Section(geometry)
section.calculate_geometric_properties()
section.display_results(fmt='.3f')
"""
post.print_results(self, fmt)
[docs] def get_area(self):
"""
:return: Cross-section area
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
area = section.get_area()
"""
return self.section_props.area
[docs] def get_perimeter(self):
"""
:return: Cross-section perimeter
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
perimeter = section.get_perimeter()
"""
return self.section_props.perimeter
[docs] def get_mass(self):
"""
:return: Cross-section mass
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
perimeter = section.get_mass()
"""
return self.section_props.mass
[docs] def get_ea(self):
"""
:return: Modulus weighted area (axial rigidity)
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
ea = section.get_ea()
"""
return self.section_props.ea
[docs] def get_q(self):
"""
:return: First moments of area about the global axis *(qx, qy)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(qx, qy) = section.get_q()
"""
return (self.section_props.qx, self.section_props.qy)
[docs] def get_ig(self):
"""
:return: Second moments of area about the global axis *(ixx_g, iyy_g, ixy_g)*
:rtype: tuple(float, float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(ixx_g, iyy_g, ixy_g) = section.get_ig()
"""
return (
self.section_props.ixx_g,
self.section_props.iyy_g,
self.section_props.ixy_g,
)
[docs] def get_c(self):
"""
:return: Elastic centroid *(cx, cy)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(cx, cy) = section.get_c()
"""
return (self.section_props.cx, self.section_props.cy)
[docs] def get_ic(self):
"""
:return: Second moments of area centroidal axis *(ixx_c, iyy_c, ixy_c)*
:rtype: tuple(float, float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(ixx_c, iyy_c, ixy_c) = section.get_ic()
"""
return (
self.section_props.ixx_c,
self.section_props.iyy_c,
self.section_props.ixy_c,
)
[docs] def get_z(self):
"""
:return: Elastic section moduli about the centroidal axis with respect to the top and
bottom fibres *(zxx_plus, zxx_minus, zyy_plus, zyy_minus)*
:rtype: tuple(float, float, float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(zxx_plus, zxx_minus, zyy_plus, zyy_minus) = section.get_z()
"""
return (
self.section_props.zxx_plus,
self.section_props.zxx_minus,
self.section_props.zyy_plus,
self.section_props.zyy_minus,
)
[docs] def get_rc(self):
"""
:return: Radii of gyration about the centroidal axis *(rx, ry)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(rx, ry) = section.get_rc()
"""
return (self.section_props.rx_c, self.section_props.ry_c)
[docs] def get_ip(self):
"""
:return: Second moments of area about the principal axis *(i11_c, i22_c)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(i11_c, i22_c) = section.get_ip()
"""
return (self.section_props.i11_c, self.section_props.i22_c)
[docs] def get_phi(self):
"""
:return: Principal bending axis angle
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
phi = section.get_phi()
"""
return self.section_props.phi
[docs] def get_zp(self):
"""
:return: Elastic section moduli about the principal axis with respect to the top and bottom
fibres *(z11_plus, z11_minus, z22_plus, z22_minus)*
:rtype: tuple(float, float, float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(z11_plus, z11_minus, z22_plus, z22_minus) = section.get_zp()
"""
return (
self.section_props.z11_plus,
self.section_props.z11_minus,
self.section_props.z22_plus,
self.section_props.z22_minus,
)
[docs] def get_rp(self):
"""
:return: Radii of gyration about the principal axis *(r11, r22)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
(r11, r22) = section.get_rp()
"""
return (self.section_props.r11_c, self.section_props.r22_c)
[docs] def get_nu_eff(self):
"""
:return: Effective Poisson's ratio
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
nu_eff = section.get_nu_eff()
"""
return self.section_props.nu_eff
[docs] def get_e_eff(self):
"""
:return: Effective elastic modulus based on area
:rtype: float
::
section = Section(geometry)
section.calculate_warping_properties()
e_eff = section.get_e_eff()
"""
return self.section_props.e_eff
[docs] def get_g_eff(self):
"""
:return: Effective shear modulus based on area
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
g_eff = section.get_g_eff()
"""
return self.section_props.g_eff
[docs] def get_j(self):
"""
:return: St. Venant torsion constant
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
j = section.get_j()
"""
return self.section_props.j
[docs] def get_sc(self):
"""
:return: Centroidal axis shear centre (elasticity approach) *(x_se, y_se)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
(x_se, y_se) = section.get_sc()
"""
if self.section_props.x_se is None:
return (None, None)
else:
# add centroid location to move section back to original location
x_se = self.section_props.x_se + self.section_props.cx
y_se = self.section_props.y_se + self.section_props.cy
return (x_se, y_se)
[docs] def get_sc_p(self):
"""
:return: Principal axis shear centre (elasticity approach) *(x11_se, y22_se)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
(x11_se, y22_se) = section.get_sc_p()
"""
if self.section_props.x11_se is None:
return (None, None)
else:
x11_se = self.section_props.x11_se
y22_se = self.section_props.y22_se
return (x11_se, y22_se)
[docs] def get_sc_t(self):
"""
:return: Centroidal axis shear centre (Trefftz's approach) *(x_st, y_st)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
(x_st, y_st) = section.get_sc_t()
"""
if self.section_props.x_st is None:
return (None, None)
else:
# add centroid location to move section back to original location
x_st = self.section_props.x_st + self.section_props.cx
y_st = self.section_props.y_st + self.section_props.cy
return (x_st, y_st)
[docs] def get_gamma(self):
"""
:return: Warping constant
:rtype: float
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
gamma = section.get_gamma()
"""
return self.section_props.gamma
[docs] def get_As(self):
"""
:return: Shear area for loading about the centroidal axis *(A_sx, A_sy)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
(A_sx, A_sy) = section.get_As()
"""
return (self.section_props.A_sx, self.section_props.A_sy)
[docs] def get_As_p(self):
"""
:return: Shear area for loading about the principal bending axis *(A_s11, A_s22)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
(A_s11, A_s22) = section.get_As_p()
"""
return (self.section_props.A_s11, self.section_props.A_s22)
[docs] def get_beta(self):
"""
:return: Monosymmetry constant for bending about both global axes *(beta_x_plus,
beta_x_minus, beta_y_plus, beta_y_minus)*. The *plus* value relates to the top flange
in compression and the *minus* value relates to the bottom flange in compression.
:rtype: tuple(float, float, float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
(beta_x_plus, beta_x_minus, beta_y_plus, beta_y_minus) = section.get_beta()
"""
return (
self.section_props.beta_x_plus,
self.section_props.beta_x_minus,
self.section_props.beta_y_plus,
self.section_props.beta_y_minus,
)
[docs] def get_beta_p(self):
"""
:return: Monosymmetry constant for bending about both principal axes *(beta_11_plus,
beta_11_minus, beta_22_plus, beta_22_minus)*. The *plus* value relates to the top
flange in compression and the *minus* value relates to the bottom flange in
compression.
:rtype: tuple(float, float, float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
(beta_11_plus, beta_11_minus, beta_22_plus, beta_22_minus) = section.get_beta_p()
"""
return (
self.section_props.beta_11_plus,
self.section_props.beta_11_minus,
self.section_props.beta_22_plus,
self.section_props.beta_22_minus,
)
[docs] def get_pc(self):
"""
:return: Centroidal axis plastic centroid *(x_pc, y_pc)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_plastic_properties()
(x_pc, y_pc) = section.get_pc()
"""
if self.section_props.x_pc is None:
return (None, None)
else:
# add centroid location to move section back to original location
x_pc = self.section_props.x_pc + self.section_props.cx
y_pc = self.section_props.y_pc + self.section_props.cy
return (x_pc, y_pc)
[docs] def get_pc_p(self):
"""
:return: Principal bending axis plastic centroid *(x11_pc, y22_pc)*
:rtype: tuple(float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_plastic_properties()
(x11_pc, y22_pc) = section.get_pc_p()
"""
if self.section_props.x11_pc is None:
return (None, None)
else:
# determine the position of the plastic centroid in the global axis
(x_pc, y_pc) = fea.global_coordinate(
self.section_props.phi,
self.section_props.x11_pc,
self.section_props.y22_pc,
)
# add centroid location to move section back to original location
return (x_pc + self.section_props.cx, y_pc + self.section_props.cy)
[docs] def get_s(self):
"""
:return: Plastic section moduli about the centroidal axis *(sxx, syy)*
:rtype: tuple(float, float)
If material properties have been specified, returns the plastic moment :math:`M_p = f_y S`.
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_plastic_properties()
(sxx, syy) = section.get_s()
"""
return (self.section_props.sxx, self.section_props.syy)
[docs] def get_sp(self):
"""
:return: Plastic section moduli about the principal bending axis *(s11, s22)*
:rtype: tuple(float, float)
If material properties have been specified, returns the plastic moment
:math:`M_p = f_y S`.
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_plastic_properties()
(s11, s22) = section.get_sp()
"""
return (self.section_props.s11, self.section_props.s22)
[docs] def get_sf(self):
"""
:return: Centroidal axis shape factors with respect to the top and bottom fibres
*(sf_xx_plus, sf_xx_minus, sf_yy_plus, sf_yy_minus)*
:rtype: tuple(float, float, float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_plastic_properties()
(sf_xx_plus, sf_xx_minus, sf_yy_plus, sf_yy_minus) = section.get_sf()
"""
return (
self.section_props.sf_xx_plus,
self.section_props.sf_xx_minus,
self.section_props.sf_yy_plus,
self.section_props.sf_yy_minus,
)
[docs] def get_sf_p(self):
"""
:return: Principal bending axis shape factors with respect to the top and bottom fibres
*(sf_11_plus, sf_11_minus, sf_22_plus, sf_22_minus)*
:rtype: tuple(float, float, float, float)
::
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_plastic_properties()
(sf_11_plus, sf_11_minus, sf_22_plus, sf_22_minus) = section.get_sf_p()
"""
return (
self.section_props.sf_11_plus,
self.section_props.sf_11_minus,
self.section_props.sf_22_plus,
self.section_props.sf_22_minus,
)
[docs] def get_stress_at_point(
self,
pt: List[float],
N=0,
Mxx=0,
Myy=0,
M11=0,
M22=0,
Mzz=0,
Vx=0,
Vy=0,
agg_func=np.average,
) -> Tuple[float]:
"""Calculates the stress at a point within an element for given design actions
and returns *(sigma_zz, tau_xz, tau_yz)*
:param pt: The point. A list of the x and y coordinate
:type pt: list[float, float]
:param float N: Axial force
:param float Vx: Shear force acting in the x-direction
:param float Vy: Shear force acting in the y-direction
:param float Mxx: Bending moment about the centroidal xx-axis
:param float Myy: Bending moment about the centroidal yy-axis
:param float M11: Bending moment about the centroidal 11-axis
:param float M22: Bending moment about the centroidal 22-axis
:param float Mzz: Torsion moment about the centroidal zz-axis
:param agg_function: A function that aggregates the stresses if the point is shared by several elements.
If the point, pt, is shared by several elements (e.g. if it is a node or on an edge), the stress
(sigma_zz, tau_xz, tau_yz) are retrieved from each element and combined according to this function.
By default, `numpy.average` is used.
:type agg_function: function, optional
:return: Resultant normal and shear stresses list[(sigma_zz, tau_xz, tau_yz)]. If a point it not in the
section then None is returned.
:rtype: Union[Tuple[float, float, float], None]
"""
sigs = self.get_stress_at_points(
[pt], N, Mxx, Myy, M11, M22, Mzz, Vx, Vy, agg_func
)
return next(sigs)
[docs] def get_stress_at_points(
self,
pts: List[List[float]],
N=0,
Mxx=0,
Myy=0,
M11=0,
M22=0,
Mzz=0,
Vx=0,
Vy=0,
agg_func=np.average,
) -> List[Tuple]:
"""Calculates the stress at a set of points within an element for given design actions
and returns *(sigma_zz, tau_xz, tau_yz)*
:param pts: The points. A list of several x and y coordinates
:type pts: list[list[float, float]]
:param float N: Axial force
:param float Vx: Shear force acting in the x-direction
:param float Vy: Shear force acting in the y-direction
:param float Mxx: Bending moment about the centroidal xx-axis
:param float Myy: Bending moment about the centroidal yy-axis
:param float M11: Bending moment about the centroidal 11-axis
:param float M22: Bending moment about the centroidal 22-axis
:param float Mzz: Torsion moment about the centroidal zz-axis
:param agg_function: A function that aggregates the stresses if the point is shared by several elements.
If the point, pt, is shared by several elements (e.g. if it is a node or on an edge), the stress
(sigma_zz, tau_xz, tau_yz) are retrieved from each element and combined according to this function.
By default, `numpy.average` is used.
:type agg_function: function, optional
:return: Resultant normal and shear stresses list[(sigma_zz, tau_xz, tau_yz)]. If a point it not in the
section then None is returned for that element in the list.
:rtype: List[Union[Tuple[float, float, float], None]]
"""
# ensure warping analysis completed for shear and torsion
if self.section_props.omega is None and (Vx != 0 or Vy != 0 or Mzz != 0):
err = "Perform a warping analysis before carrying out a stress analysis "
err += "with non-zero shear forces or torsion moment."
raise RuntimeError(err)
action = {
"N": N,
"Mxx": Mxx,
"Myy": Myy,
"M11": M11,
"M22": M22,
"Mzz": Mzz,
"Vx": Vx,
"Vy": Vy,
}
sect_prop = {
"ea": self.section_props.ea,
"cx": self.section_props.cx,
"cy": self.section_props.cy,
"ixx": self.section_props.ixx_c,
"iyy": self.section_props.iyy_c,
"ixy": self.section_props.ixy_c,
"i11": self.section_props.i11_c,
"i22": self.section_props.i22_c,
"phi": self.section_props.phi,
"j": self.section_props.j,
"Delta_s": self.section_props.Delta_s,
"nu": self.section_props.nu_eff,
}
for pt in pts:
query_geom = Point(pt)
tri_ids = self.mesh_search_tree.query(query_geom, predicate="intersects")
if len(tri_ids) == 0:
sig = None
elif len(tri_ids) == 1:
tri = self.elements[tri_ids[0]]
if self.section_props.omega is None:
omega_el = None
psi_shear_el = None
phi_shear_el = None
else:
omega_el = self.section_props.omega[tri.node_ids]
psi_shear_el = self.section_props.psi_shear[tri.node_ids]
phi_shear_el = self.section_props.phi_shear[tri.node_ids]
sig = tri.local_element_stress(
p=pt,
**action,
**sect_prop,
omega=omega_el,
psi_shear=psi_shear_el,
phi_shear=phi_shear_el,
)
else:
sigs = []
for idx in tri_ids:
tri = self.elements[idx]
if self.section_props.omega is None:
omega_el = None
psi_shear_el = None
phi_shear_el = None
else:
omega_el = self.section_props.omega[tri.node_ids]
psi_shear_el = self.section_props.psi_shear[tri.node_ids]
phi_shear_el = self.section_props.phi_shear[tri.node_ids]
sigs.append(
tri.local_element_stress(
p=pt,
**action,
**sect_prop,
omega=omega_el,
psi_shear=psi_shear_el,
phi_shear=phi_shear_el,
)
)
sig = (
agg_func([sig[0] for sig in sigs]),
agg_func([sig[1] for sig in sigs]),
agg_func([sig[2] for sig in sigs]),
)
yield sig
[docs]class PlasticSection:
"""Class for the plastic analysis of cross-sections.
Stores the finite element geometry and material information and provides methods to compute the
plastic section properties.
:param section: Section object
:type section: :class:`~sectionproperties.analysis.section.Section`
:cvar geometry: Deep copy of the Section geometry object provided to the constructor
:vartype geometry: :class:`~sectionproperties.pre.geometry.Geometry`
:cvar materials: A list of material properties corresponding to various regions in the geometry
and mesh.
:vartype materials: list[:class:`~sectionproperties.pre.pre.Material`]
:cvar mesh: Mesh dict returned by triangle
:vartype mesh: dict(mesh)
:cvar mesh_nodes: Array of node coordinates from the mesh
:vartype mesh_nodes: :class:`numpy.ndarray`
:cvar mesh_elements: Array of connectivities from the mesh
:vartype mesh_elements: :class:`numpy.ndarray`
:cvar elements: List of finite element objects describing the cross-section mesh
:vartype elements: list[:class:`~sectionproperties.analysis.fea.Tri6`]
:cvar f_top: Current force in the top region
:type f_top: float
:cvar c_top: Centroid of the force in the top region *(c_top_x, c_top_y)*
:type c_top: list[float, float]
:cvar c_bot: Centroid of the force in the bottom region *(c_bot_x, c_bot_y)*
:type c_bot: list[float, float]
"""
def __init__(
self, geom: Union[section_geometry.Geometry, section_geometry.CompoundGeometry]
):
"""Inits the PlasticSection class."""
self.geometry = geom.align_center()
self.geometry.compile_geometry()
# initialize variables to be defined later within calculate_plastic_force
self.c_top = [0.0, 0.0]
self.c_bot = [0.0, 0.0]
self.f_top = 0.0
self.f_bot = 0.0
[docs] def calculate_centroid(self, elements):
"""Calculates the elastic centroid from a list of finite elements.
:param elements: A list of Tri6 finite elements.
:type elements: list[:class:`~sectionproperties.analysis.fea.Tri6`]
:return: A tuple containing the x and y location of the elastic centroid.
:rtype: tuple(float, float)
"""
ea = 0
qx = 0
qy = 0
# General case
# loop through all the geometries
if isinstance(self.geometry, section_geometry.CompoundGeometry):
for geom in self.geometry.geoms:
e = geom.material.elastic_modulus
area = geom.calculate_area()
cx, cy = geom.calculate_centroid()
ea += area * e
qx += cx * e
qy += cy * e
# Special case (just one geometry w/ one material)
else:
e = self.geometry.material.elastic_modulus
area = self.geometry.calculate_area()
cx, cy = self.geometry.calculate_centroid()
ea += area * e
qx += cx * area * e
qy += cy * area * e
return (qy / ea, qx, ea)
[docs] def calculate_plastic_properties(self, section, verbose, progress=None):
"""Calculates the location of the plastic centroid with respect to the centroidal and
principal bending axes, the plastic section moduli and shape factors and stores the results
to the supplied :class:`~sectionproperties.analysis.section.Section` object.
:param section: Cross section object that uses the same geometry and materials
specified in the class constructor
:type section: :class:`~sectionproperties.analysis.section.Section`
:param bool verbose: If set to True, the number of iterations required for each plastic
axis is printed to the terminal.
"""
if progress is not None:
task = progress.add_task(
description="[red]Calculating plastic properties",
total=4,
)
# 1) Calculate plastic properties for centroidal axis
# calculate distances to the extreme fibres
fibres = self.calculate_extreme_fibres(
0
) # Angle = 0 ; returns bounding box Provides upper/lower bound for brentq algo
# 1a) Calculate x-axis plastic centroid
(y_pc, r, f, c_top, c_bot) = self.pc_algorithm(
np.array([1, 0]), fibres[2:], 1, verbose
) # fibres[2:] = ymin, ymax
self.check_convergence(r, "x-axis")
section.section_props.y_pc = y_pc
section.section_props.sxx = f * abs(c_top[1] - c_bot[1])
if verbose:
self.print_verbose(y_pc, r, "x-axis") # Location of axis for each iteration
if progress is not None:
progress.update(task, advance=1)
# 1b) Calculate y-axis plastic centroid
(x_pc, r, f, c_top, c_bot) = self.pc_algorithm(
np.array([0, 1]), fibres[0:2], 2, verbose
) # fibres[0:2] = xmin, xmax
self.check_convergence(r, "y-axis")
section.section_props.x_pc = x_pc
section.section_props.syy = f * abs(c_top[0] - c_bot[0])
if verbose:
self.print_verbose(x_pc, r, "y-axis") # Location of axis for each iteration
if progress is not None:
progress.update(task, advance=1)
# 2) Calculate plastic properties for principal axis
# convert principal axis angle to radians
angle = section.section_props.phi * np.pi / 180
# unit vectors in the axis directions
ux = np.array([np.cos(angle), np.sin(angle)]) # Unit vectors
uy = np.array([-np.sin(angle), np.cos(angle)]) # Unit vectors
# calculate distances to the extreme fibres in the principal axis
fibres = self.calculate_extreme_fibres(section.section_props.phi)
# 2a) Calculate 11-axis plastic centroid
(y22_pc, r, f, c_top, c_bot) = self.pc_algorithm(
ux, fibres[2:], 1, verbose
) # ux
# calculate the centroids in the principal coordinate system
c_top_p = fea.principal_coordinate(
section.section_props.phi, c_top[0], c_top[1]
)
c_bot_p = fea.principal_coordinate(
section.section_props.phi, c_bot[0], c_bot[1]
)
self.check_convergence(r, "11-axis")
section.section_props.y22_pc = y22_pc
section.section_props.s11 = f * abs(c_top_p[1] - c_bot_p[1])
if verbose:
self.print_verbose(y22_pc, r, "11-axis")
if progress is not None:
progress.update(task, advance=1)
# 2b) Calculate 22-axis plastic centroid
(x11_pc, r, f, c_top, c_bot) = self.pc_algorithm(
uy, fibres[0:2], 2, verbose
) # uy
# calculate the centroids in the principal coordinate system
c_top_p = fea.principal_coordinate(
section.section_props.phi, c_top[0], c_top[1]
)
c_bot_p = fea.principal_coordinate(
section.section_props.phi, c_bot[0], c_bot[1]
)
self.check_convergence(r, "22-axis")
section.section_props.x11_pc = x11_pc
section.section_props.s22 = f * abs(c_top_p[0] - c_bot_p[0])
if verbose:
self.print_verbose(x11_pc, r, "22-axis")
if progress is not None:
progress.update(task, advance=1)
# if there are no materials specified, calculate shape factors
if list(set(section.materials)) == [pre.DEFAULT_MATERIAL]:
section.section_props.sf_xx_plus = (
section.section_props.sxx / section.section_props.zxx_plus
)
section.section_props.sf_xx_minus = (
section.section_props.sxx / section.section_props.zxx_minus
)
section.section_props.sf_yy_plus = (
section.section_props.syy / section.section_props.zyy_plus
)
section.section_props.sf_yy_minus = (
section.section_props.syy / section.section_props.zyy_minus
)
section.section_props.sf_11_plus = (
section.section_props.s11 / section.section_props.z11_plus
)
section.section_props.sf_11_minus = (
section.section_props.s11 / section.section_props.z11_minus
)
section.section_props.sf_22_plus = (
section.section_props.s22 / section.section_props.z22_plus
)
section.section_props.sf_22_minus = (
section.section_props.s22 / section.section_props.z22_minus
)
if progress is not None:
progress.update(
task,
description="[bold green]:white_check_mark: Plastic analysis complete",
)
[docs] def check_convergence(self, root_result, axis):
"""Checks that the function solver converged and if not, raises a helpful error.
:param root_result: Result object from the root finder
:type root_result: :class:`scipy.optimize.RootResults`
:param string axis: Axis being considered by the function solver
:raises RuntimeError: If the function solver did not converge
"""
if not root_result.converged:
msg = "Plastic centroid calculation about the {0}".format(axis)
msg += " failed. Contact robbie.vanleeuwen@gmail.com with your"
msg += " analysis parameters. Termination flag: {0}".format(
root_result.flag
)
raise RuntimeError(msg)
[docs] def print_verbose(self, d, root_result, axis):
"""Prints information related to the function solver convergence to the terminal.
:param float d: Location of the plastic centroid axis
:param root_result: Result object from the root finder
:type root_result: :class:`scipy.optimize.RootResults`
:param string axis: Axis being considered by the function solver
"""
msg = "---{0} plastic centroid calculation converged at ".format(axis)
msg += "{0:.5e} in {1} iterations.".format(d, root_result.iterations)
print(msg)
[docs] def calculate_extreme_fibres(self, angle):
"""Calculates the locations of the extreme fibres along and perpendicular to the axis
specified by 'angle' using the elements stored in `self.elements`.
:param float angle: Angle (in degrees) along which to calculate the extreme fibre locations
:return: The location of the extreme fibres parallel (u) and perpendicular (v) to the axis
*(u_min, u_max, v_min, v_max)*
:rtype: tuple(float, float, float, float)
"""
# loop through all nodes in the mesh
for (i, pt) in enumerate(self.geometry.points):
# determine the coordinate of the point wrt the axis
(u, v) = fea.principal_coordinate(angle, pt[0], pt[1])
# initialise min, max variables
if i == 0:
u_min = u
u_max = u
v_min = v
v_max = v
# update the mins and maxes where necessary
u_min = min(u_min, u)
u_max = max(u_max, u)
v_min = min(v_min, v)
v_max = max(v_max, v)
return (u_min, u_max, v_min, v_max)
[docs] def evaluate_force_eq(self, d, u, u_p, verbose):
"""Given a distance *d* from the centroid to an axis (defined by unit vector *u*), creates
a mesh including the new and axis and calculates the force equilibrium. The resultant
force, as a ratio of the total force, is returned.
:param float d: Distance from the centroid to current axis
:param u: Unit vector defining the direction of the axis
:type u: :class:`numpy.ndarray`
:param u_p: Unit vector perpendicular to the direction of the axis
:type u_p: :class:`numpy.ndarray`
:param bool verbose: If set to True, the number of iterations required for each plastic
axis is printed to the terminal.
:return: The force equilibrium norm
:rtype: float
"""
p = np.array(
[d * u_p[0], d * u_p[1]]
) # p finding a point on the axis by scaling the perpendicular
# calculate force equilibrium
(f_top, f_bot) = self.calculate_plastic_force(u, p)
# calculate the force norm
f_norm = (f_top - f_bot) / (f_top + f_bot) # Going down to zero
# print verbose results
if verbose:
print("d = {0}; f_norm = {1}".format(d, f_norm))
# return the force norm
return f_norm # This is the result that is the target for the root finding algorithm
[docs] def pc_algorithm(self, u, dlim, axis, verbose):
"""An algorithm used for solving for the location of the plastic centroid. The algorithm
searches for the location of the axis, defined by unit vector *u* and within the section
depth, that satisfies force equilibrium.
:param u: Unit vector defining the direction of the axis
:type u: :class:`numpy.ndarray`
:param dlim: List [dmax, dmin] containing the distances from the centroid to the extreme
fibres perpendicular to the axis
:type dlim: list[float, float]
:param int axis: The current axis direction: 1 (e.g. x or 11) or 2 (e.g. y or 22)
:param bool verbose: If set to True, the number of iterations required for each plastic
axis is printed to the terminal.
:return: The distance to the plastic centroid axis *d*, the result object *r*, the force in
the top of the section *f_top* and the location of the centroids of the top and bottom
areas *c_top* and *c_bottom*
:rtype: tuple(float, :class:`scipy.optimize.RootResults`, float, list[float, float],
list[float, float])
"""
# calculate vector perpendicular to u
if axis == 1:
u_p = np.array([-u[1], u[0]])
else:
u_p = np.array([u[1], -u[0]])
a = dlim[0] # Upper and lower bound for algorithm
b = dlim[1] # Upper and lower bound for algorithm
(
d,
r,
) = brentq( # d = neutral axis depth, measured from the origin; r is the convegence indicator
self.evaluate_force_eq,
a,
b,
args=(u, u_p, verbose),
full_output=True,
disp=False, # Unit vector w/ unit vector perpendicular
xtol=1e-6,
rtol=1e-6,
)
return (d, r, self.f_top, self.c_top, self.c_bot)
[docs] def calculate_plastic_force(
self, u: np.ndarray, p: np.ndarray
) -> Tuple[float, float]:
"""Sums the forces above and below the axis defined by unit vector *u* and point *p*. Also
returns the force centroid of the forces above and below the axis.
:param elements: A list of Tri6 finite elements.
:type elements: list[:class:`~sectionproperties.analysis.fea.Tri6`]
:param u: Unit vector defining the direction of the axis
:type u: :class:`numpy.ndarray`
:param p: Point on the axis
:type p: :class:`numpy.ndarray`
:return: Force in the top and bottom areas *(f_top, f_bot)*
:rtype: tuple(float, float)
"""
# initialise variables
(f_top, f_bot) = (0, 0)
(a_top, a_bot) = (0, 0)
(qx_top, qx_bot) = (0, 0)
(qy_top, qy_bot) = (0, 0)
top_geoms, bot_geoms = self.geometry.split_section(point_i=p, vector=u)
if top_geoms:
for top_geom in top_geoms:
e = top_geom.material.elastic_modulus
f_y = top_geom.material.yield_strength
area_top = top_geom.calculate_area()
a_top += area_top
cx, cy = top_geom.calculate_centroid()
qx_top += cy * area_top
qy_top += cx * area_top
f_top += f_y * area_top
if bot_geoms:
for bot_geom in bot_geoms:
e = bot_geom.material.elastic_modulus
f_y = bot_geom.material.yield_strength
area_bot = bot_geom.calculate_area()
a_bot += area_bot
cx, cy = bot_geom.calculate_centroid()
qx_bot += cy * area_bot
qy_bot += cx * area_bot
f_bot += f_y * area_bot
try:
self.c_top = [qy_top / a_top, qx_top / a_top]
self.f_top = f_top
except ZeroDivisionError:
self.c_top = [0, 0]
self.f_top = 0
try:
self.c_bot = [qy_bot / a_bot, qx_bot / a_bot]
except ZeroDivisionError:
self.c_bot = [0, 0]
return (f_top, f_bot)
[docs]class StressPost:
"""Class for post-processing finite element stress results.
A StressPost object is created when a stress analysis is carried out and is returned as an
object to allow post-processing of the results. The StressPost object creates a deep copy of
the MaterialGroups within the cross-section to allow the calculation of stresses for each
material. Methods for post-processing the calculated stresses are provided.
:param section: Cross section object for stress calculation
:type section: :class:`~sectionproperties.analysis.section.Section`
:cvar section: Cross section object for stress calculation
:vartype section: :class:`~sectionproperties.analysis.section.Section`
:cvar material_groups: A deep copy of the `section` material groups to allow a new stress
analysis
:vartype material_groups: list[:class:`~sectionproperties.pre.pre.MaterialGroup`]
"""
def __init__(self, section):
"""Inits the StressPost class."""
self.section = section
# make a deep copy of the material groups to the StressPost object such that stress results
# can be saved to a new material group
self.material_groups = copy.deepcopy(section.material_groups)
[docs] def plot_stress_contour(self, sigs, title, cmap, normalize, **kwargs):
"""Plots filled stress contours over the finite element mesh.
:param sigs: List of nodal stress values for each material
:type sigs: list[:class:`numpy.ndarray`]
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axe object
:rtype: :class:`matplotlib.axes`
"""
# create plot and setup the plot
with post.plotting_context(title=title, **kwargs) as (fig, ax):
# set up the colormap
cmap = cm.get_cmap(name=cmap)
# create triangulation
triang = tri.Triangulation(
self.section.mesh_nodes[:, 0],
self.section.mesh_nodes[:, 1],
self.section.mesh_elements[:, 0:3],
)
# determine minimum and maximum stress values for the contour list
sig_min = min([min(x) for x in sigs])
sig_max = max([max(x) for x in sigs])
v = np.linspace(sig_min, sig_max, 15, endpoint=True)
if np.isclose(v[0], v[-1], atol=1e-12):
v = 15
ticks = None
else:
ticks = v
norm = None
if normalize:
norm = CenteredNorm()
# plot the filled contour, looping through the materials
for (i, sig) in enumerate(sigs):
# create and set the mask for the current material
mask_array = np.ones(len(self.section.elements), dtype=bool)
mask_array[self.material_groups[i].el_ids] = False
triang.set_mask(mask_array)
# plot the filled contour
trictr = ax.tricontourf(triang, sig, v, cmap=cmap, norm=norm)
# display the colourbar
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
fig.colorbar(trictr, label="Stress", format="%.4e", ticks=ticks, cax=cax)
# plot the finite element mesh
self.section.plot_mesh(materials=False, **dict(kwargs, ax=ax))
return ax
[docs] def plot_stress_vector(self, sigxs, sigys, title, cmap, normalize, **kwargs):
r"""Plots stress vectors over the finite element mesh.
:param sigxs: List of x-components of the nodal stress values for each material
:type sigxs: list[:class:`numpy.ndarray`]
:param sigys: List of y-components of the nodal stress values for each material
:type sigys: list[:class:`numpy.ndarray`]
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
"""
# create plot and setup the plot
with post.plotting_context(title=title, **kwargs) as (fig, ax):
# set up the colormap
cmap = cm.get_cmap(name=cmap)
# initialise quiver plot list max scale
quiv_list = []
max_scale = 0
norm = None
if normalize:
norm = CenteredNorm()
# plot the vectors
for (i, sigx) in enumerate(sigxs):
sigy = sigys[i]
# scale the colour with respect to the magnitude of the vector
c = np.hypot(sigx, sigy)
quiv = ax.quiver(
self.section.mesh_nodes[:, 0],
self.section.mesh_nodes[:, 1],
sigx,
sigy,
c,
cmap=cmap,
norm=norm,
)
# get the scale and store the max value
quiv._init()
max_scale = max(max_scale, quiv.scale)
quiv_list.append(quiv)
# update the colormap values
if i == 0:
c_min = min(c)
c_max = max(c)
else:
c_min = min(c_min, min(c))
c_max = max(c_max, max(c))
# apply the scale
for quiv_plot in quiv_list:
quiv_plot.scale = max_scale
# apply the colourbar
v1 = np.linspace(c_min, c_max, 15, endpoint=True)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
fig.colorbar(quiv, label="Stress", format="%.4e", ticks=v1, cax=cax)
# plot the finite element mesh
self.section.plot_mesh(materials=False, **dict(kwargs, ax=ax))
return ax
[docs] def get_stress(self):
r"""Returns the stresses within each material belonging to the current
:class:`~sectionproperties.analysis.section.StressPost` object.
:return: A list of dictionaries containing the cross-section stresses for each material.
:rtype: list[dict]
A dictionary is returned for each material in the cross-section, containing the following
keys and values:
* *'Material'*: Material name
* *'sig_zz_n'*: Normal stress :math:`\sigma_{zz,N}` resulting from the axial load :math:`N`
* *'sig_zz_mxx'*: Normal stress :math:`\sigma_{zz,Mxx}` resulting from the bending moment
:math:`M_{xx}`
* *'sig_zz_myy'*: Normal stress :math:`\sigma_{zz,Myy}` resulting from the bending moment
:math:`M_{yy}`
* *'sig_zz_m11'*: Normal stress :math:`\sigma_{zz,M11}` resulting from the bending moment
:math:`M_{11}`
* *'sig_zz_m22'*: Normal stress :math:`\sigma_{zz,M22}` resulting from the bending moment
:math:`M_{22}`
* *'sig_zz_m'*: Normal stress :math:`\sigma_{zz,\Sigma M}` resulting from all bending
moments
* *'sig_zx_mzz'*: *x*-component of the shear stress :math:`\sigma_{zx,Mzz}` resulting from
the torsion moment
* *'sig_zy_mzz'*: *y*-component of the shear stress :math:`\sigma_{zy,Mzz}` resulting from
the torsion moment
* *'sig_zxy_mzz'*: Resultant shear stress :math:`\sigma_{zxy,Mzz}` resulting from the
torsion moment
* *'sig_zx_vx'*: *x*-component of the shear stress :math:`\sigma_{zx,Vx}` resulting from
the shear force :math:`V_{x}`
* *'sig_zy_vx'*: *y*-component of the shear stress :math:`\sigma_{zy,Vx}` resulting from
the shear force :math:`V_{x}`
* *'sig_zxy_vx'*: Resultant shear stress :math:`\sigma_{zxy,Vx}` resulting from the shear
force :math:`V_{x}`
* *'sig_zx_vy'*: *x*-component of the shear stress :math:`\sigma_{zx,Vy}` resulting from
the shear force :math:`V_{y}`
* *'sig_zy_vy'*: *y*-component of the shear stress :math:`\sigma_{zy,Vy}` resulting from
the shear force :math:`V_{y}`
* *'sig_zxy_vy'*: Resultant shear stress :math:`\sigma_{zxy,Vy}` resulting from the shear
force :math:`V_{y}`
* *'sig_zx_v'*: *x*-component of the shear stress :math:`\sigma_{zx,\Sigma V}` resulting
from all shear forces
* *'sig_zy_v'*: *y*-component of the shear stress :math:`\sigma_{zy,\Sigma V}` resulting
from all shear forces
* *'sig_zxy_v'*: Resultant shear stress :math:`\sigma_{zxy,\Sigma V}` resulting from all
shear forces
* *'sig_zz'*: Combined normal stress :math:`\sigma_{zz}` resulting from all actions
* *'sig_zx'*: *x*-component of the shear stress :math:`\sigma_{zx}` resulting from all
actions
* *'sig_zy'*: *y*-component of the shear stress :math:`\sigma_{zy}` resulting from all
actions
* *'sig_zxy'*: Resultant shear stress :math:`\sigma_{zxy}` resulting from all actions
* *'sig_1'*: Major principal stress :math:`\sigma_{1}` resulting from all actions
* *'sig_3'*: Minor principal stress :math:`\sigma_{3}` resulting from all actions
* *'sig_vm'*: von Mises stress :math:`\sigma_{vM}` resulting from all actions
The following example returns stresses for each material within a composite section, note
that a result is generated for each node in the mesh for all materials irrespective of
whether the materials exists at that point or not.
::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(
N=50e3, Mxx=-5e6, M22=2.5e6, Mzz=0.5e6, Vx=10e3, Vy=5e3
)
stresses = stress_post.get_stress()
print("Number of nodes: {0}".format(section.num_nodes))
for stress in stresses:
print('Material: {0}'.format(stress['Material']))
print('List Size: {0}'.format(len(stress['sig_zz_n'])))
print('Normal Stresses: {0}'.format(stress['sig_zz_n']))
print('von Mises Stresses: {0}'.format(stress['sig_vm']))
::
$ Number of nodes: 2465
$ Material: Timber
$ List Size: 2465
$ Normal Stresses: [0.76923077 0.76923077 0.76923077 ... 0.76923077 0.76923077 0.76923077]
$ von Mises Stresses: [7.6394625 5.38571866 3.84784964 ... 3.09532948 3.66992556 2.81976647]
$ Material: Steel
$ List Size: 2465
$ Normal Stresses: [19.23076923 0. 0. ... 0. 0. 0.]
$ von Mises Stresses: [134.78886419 0. 0. ... 0. 0. 0.]
"""
stress = []
for group in self.material_groups:
stress.append(
{
"Material": group.material.name,
"sig_zz_n": group.stress_result.sig_zz_n,
"sig_zz_mxx": group.stress_result.sig_zz_mxx,
"sig_zz_myy": group.stress_result.sig_zz_myy,
"sig_zz_m11": group.stress_result.sig_zz_m11,
"sig_zz_m22": group.stress_result.sig_zz_m22,
"sig_zz_m": group.stress_result.sig_zz_m,
"sig_zx_mzz": group.stress_result.sig_zx_mzz,
"sig_zy_mzz": group.stress_result.sig_zy_mzz,
"sig_zxy_mzz": group.stress_result.sig_zxy_mzz,
"sig_zx_vx": group.stress_result.sig_zx_vx,
"sig_zy_vx": group.stress_result.sig_zy_vx,
"sig_zxy_vx": group.stress_result.sig_zxy_vx,
"sig_zx_vy": group.stress_result.sig_zx_vy,
"sig_zy_vy": group.stress_result.sig_zy_vy,
"sig_zxy_vy": group.stress_result.sig_zxy_vy,
"sig_zx_v": group.stress_result.sig_zx_v,
"sig_zy_v": group.stress_result.sig_zy_v,
"sig_zxy_v": group.stress_result.sig_zxy_v,
"sig_zz": group.stress_result.sig_zz,
"sig_zx": group.stress_result.sig_zx,
"sig_zy": group.stress_result.sig_zy,
"sig_zxy": group.stress_result.sig_zxy,
"sig_1": group.stress_result.sig_1,
"sig_3": group.stress_result.sig_3,
"sig_vm": group.stress_result.sig_vm,
}
)
return stress
[docs] def plot_stress_n_zz(
self,
title=r"Stress Contour Plot - $\sigma_{zz,N}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the normal stress :math:`\sigma_{zz,N}` resulting from the
axial load :math:`N`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the normal stress within a 150x90x12 UA section resulting from
an axial force of 10 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(N=10e3)
stress_post.plot_stress_n_zz()
.. figure:: ../images/stress/stress_n_zz.png
:align: center
:scale: 75 %
Contour plot of the axial stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zz_n)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_mxx_zz(
self,
title=r"Stress Contour Plot - $\sigma_{zz,Mxx}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the normal stress :math:`\sigma_{zz,Mxx}` resulting from the
bending moment :math:`M_{xx}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the normal stress within a 150x90x12 UA section resulting from
a bending moment about the x-axis of 5 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mxx=5e6)
stress_post.plot_stress_mxx_zz()
.. figure:: ../images/stress/stress_mxx_zz.png
:align: center
:scale: 75 %
Contour plot of the bending stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zz_mxx)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_myy_zz(
self,
title=r"Stress Contour Plot - $\sigma_{zz,Myy}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the normal stress :math:`\sigma_{zz,Myy}` resulting from the
bending moment :math:`M_{yy}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the normal stress within a 150x90x12 UA section resulting from
a bending moment about the y-axis of 2 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Myy=2e6)
stress_post.plot_stress_myy_zz()
.. figure:: ../images/stress/stress_myy_zz.png
:align: center
:scale: 75 %
Contour plot of the bending stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zz_myy)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_m11_zz(
self,
title=r"Stress Contour Plot - $\sigma_{zz,M11}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the normal stress :math:`\sigma_{zz,M11}` resulting from the
bending moment :math:`M_{11}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the normal stress within a 150x90x12 UA section resulting from
a bending moment about the 11-axis of 5 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(M11=5e6)
stress_post.plot_stress_m11_zz()
.. figure:: ../images/stress/stress_m11_zz.png
:align: center
:scale: 75 %
Contour plot of the bending stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zz_m11)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_m22_zz(
self,
title=r"Stress Contour Plot - $\sigma_{zz,M22}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the normal stress :math:`\sigma_{zz,M22}` resulting from the
bending moment :math:`M_{22}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the normal stress within a 150x90x12 UA section resulting from
a bending moment about the 22-axis of 2 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(M22=5e6)
stress_post.plot_stress_m22_zz()
.. figure:: ../images/stress/stress_m22_zz.png
:align: center
:scale: 75 %
Contour plot of the bending stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zz_m22)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_m_zz(
self,
title=r"Stress Contour Plot - $\sigma_{zz,\Sigma M}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the normal stress :math:`\sigma_{zz,\Sigma M}` resulting from
all bending moments :math:`M_{xx} + M_{yy} + M_{11} + M_{22}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the normal stress within a 150x90x12 UA section resulting from
a bending moment about the x-axis of 5 kN.m, a bending moment about the y-axis of 2 kN.m
and a bending moment of 3 kN.m about the 11-axis::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mxx=5e6, Myy=2e6, M11=3e6)
stress_post.plot_stress_m_zz()
.. figure:: ../images/stress/stress_m_zz.png
:align: center
:scale: 75 %
Contour plot of the bending stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zz_m)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_mzz_zx(
self,
title=r"Stress Contour Plot - $\sigma_{zx,Mzz}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *x*-component of the shear stress :math:`\sigma_{zx,Mzz}`
resulting from the torsion moment :math:`M_{zz}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the x-component of the shear stress within a 150x90x12 UA
section resulting from a torsion moment of 1 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mzz=1e6)
stress_post.plot_stress_mzz_zx()
.. figure:: ../images/stress/stress_mzz_zx.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zx_mzz)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_mzz_zy(
self,
title=r"Stress Contour Plot - $\sigma_{zy,Mzz}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *y*-component of the shear stress :math:`\sigma_{zy,Mzz}`
resulting from the torsion moment :math:`M_{zz}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the y-component of the shear stress within a 150x90x12 UA
section resulting from a torsion moment of 1 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mzz=1e6)
stress_post.plot_stress_mzz_zy()
.. figure:: ../images/stress/stress_mzz_zy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zy_mzz)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_mzz_zxy(
self,
title=r"Stress Contour Plot - $\sigma_{zxy,Mzz}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the resultant shear stress :math:`\sigma_{zxy,Mzz}` resulting
from the torsion moment :math:`M_{zz}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots a contour of the resultant shear stress within a 150x90x12 UA
section resulting from a torsion moment of 1 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mzz=1e6)
stress_post.plot_stress_mzz_zxy()
.. figure:: ../images/stress/stress_mzz_zxy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zxy_mzz)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_vector_mzz_zxy(
self,
title=r"Stress Vector Plot - $\sigma_{zxy,Mzz}$",
cmap="YlOrBr",
normalize=False,
**kwargs,
):
r"""Produces a vector plot of the resultant shear stress :math:`\sigma_{zxy,Mzz}` resulting
from the torsion moment :math:`M_{zz}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example generates a vector plot of the shear stress within a 150x90x12 UA
section resulting from a torsion moment of 1 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mzz=1e6)
stress_post.plot_vector_mzz_zxy()
.. figure:: ../images/stress/vector_mzz_zxy.png
:align: center
:scale: 75 %
Vector plot of the shear stress.
"""
sigxs = []
sigys = []
for group in self.material_groups:
sigxs.append(group.stress_result.sig_zx_mzz)
sigys.append(group.stress_result.sig_zy_mzz)
return self.plot_stress_vector(sigxs, sigys, title, cmap, normalize, **kwargs)
[docs] def plot_stress_vx_zx(
self,
title=r"Stress Contour Plot - $\sigma_{zx,Vx}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *x*-component of the shear stress :math:`\sigma_{zx,Vx}`
resulting from the shear force :math:`V_{x}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the x-component of the shear stress within a 150x90x12 UA
section resulting from a shear force in the x-direction of 15 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vx=15e3)
stress_post.plot_stress_vx_zx()
.. figure:: ../images/stress/stress_vx_zx.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zx_vx)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_vx_zy(
self,
title=r"Stress Contour Plot - $\sigma_{zy,Vx}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *y*-component of the shear stress :math:`\sigma_{zy,Vx}`
resulting from the shear force :math:`V_{x}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the y-component of the shear stress within a 150x90x12 UA
section resulting from a shear force in the x-direction of 15 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vx=15e3)
stress_post.plot_stress_vx_zy()
.. figure:: ../images/stress/stress_vx_zy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zy_vx)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_vx_zxy(
self,
title=r"Stress Contour Plot - $\sigma_{zz,Myy}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the resultant shear stress :math:`\sigma_{zxy,Vx}` resulting
from the shear force :math:`V_{x}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots a contour of the resultant shear stress within a 150x90x12 UA
section resulting from a shear force in the x-direction of 15 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vx=15e3)
stress_post.plot_stress_vx_zxy()
.. figure:: ../images/stress/stress_vx_zxy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
title = r"Stress Contour Plot - $\sigma_{zxy,Vx}$"
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zxy_vx)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_vector_vx_zxy(
self,
title=r"Stress Vector Plot - $\sigma_{zxy,Vx}$",
cmap="YlOrBr",
normalize=False,
**kwargs,
):
r"""Produces a vector plot of the resultant shear stress :math:`\sigma_{zxy,Vx}` resulting
from the shear force :math:`V_{x}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example generates a vector plot of the shear stress within a 150x90x12 UA
section resulting from a shear force in the x-direction of 15 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vx=15e3)
stress_post.plot_vector_vx_zxy()
.. figure:: ../images/stress/vector_vx_zxy.png
:align: center
:scale: 75 %
Vector plot of the shear stress.
"""
sigxs = []
sigys = []
for group in self.material_groups:
sigxs.append(group.stress_result.sig_zx_vx)
sigys.append(group.stress_result.sig_zy_vx)
return self.plot_stress_vector(sigxs, sigys, title, cmap, normalize, **kwargs)
[docs] def plot_stress_vy_zx(
self,
title=r"Stress Contour Plot - $\sigma_{zx,Vy}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *x*-component of the shear stress :math:`\sigma_{zx,Vy}`
resulting from the shear force :math:`V_{y}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the x-component of the shear stress within a 150x90x12 UA
section resulting from a shear force in the y-direction of 30 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vy=30e3)
stress_post.plot_stress_vy_zx()
.. figure:: ../images/stress/stress_vy_zx.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zx_vy)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_vy_zy(
self,
title=r"Stress Contour Plot - $\sigma_{zy,Vy}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *y*-component of the shear stress :math:`\sigma_{zy,Vy}`
resulting from the shear force :math:`V_{y}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the y-component of the shear stress within a 150x90x12 UA
section resulting from a shear force in the y-direction of 30 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vy=30e3)
stress_post.plot_stress_vy_zy()
.. figure:: ../images/stress/stress_vy_zy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zy_vy)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_vy_zxy(
self,
title=r"Stress Contour Plot - $\sigma_{zxy,Vy}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the resultant shear stress :math:`\sigma_{zxy,Vy}` resulting
from the shear force :math:`V_{y}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots a contour of the resultant shear stress within a 150x90x12 UA
section resulting from a shear force in the y-direction of 30 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vy=30e3)
stress_post.plot_stress_vy_zxy()
.. figure:: ../images/stress/stress_vy_zxy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zxy_vy)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_vector_vy_zxy(
self,
title=r"Stress Vector Plot - $\sigma_{zxy,Vy}$",
cmap="YlOrBr",
normalize=False,
**kwargs,
):
r"""Produces a vector plot of the resultant shear stress :math:`\sigma_{zxy,Vy}` resulting
from the shear force :math:`V_{y}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example generates a vector plot of the shear stress within a 150x90x12 UA
section resulting from a shear force in the y-direction of 30 kN::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vy=30e3)
stress_post.plot_vector_vy_zxy()
.. figure:: ../images/stress/vector_vy_zxy.png
:align: center
:scale: 75 %
Vector plot of the shear stress.
"""
sigxs = []
sigys = []
for group in self.material_groups:
sigxs.append(group.stress_result.sig_zx_vy)
sigys.append(group.stress_result.sig_zy_vy)
return self.plot_stress_vector(sigxs, sigys, title, cmap, normalize, **kwargs)
[docs] def plot_stress_v_zx(
self,
title=r"Stress Contour Plot - $\sigma_{zx,\Sigma V}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *x*-component of the shear stress
:math:`\sigma_{zx,\Sigma V}` resulting from the sum of the applied shear forces
:math:`V_{x} + V_{y}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the x-component of the shear stress within a 150x90x12 UA
section resulting from a shear force of 15 kN in the x-direction and 30 kN in the
y-direction::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vx=15e3, Vy=30e3)
stress_post.plot_stress_v_zx()
.. figure:: ../images/stress/stress_v_zx.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zx_v)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_v_zy(
self,
title=r"Stress Contour Plot - $\sigma_{zy,\Sigma V}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *y*-component of the shear stress
:math:`\sigma_{zy,\Sigma V}` resulting from the sum of the applied shear forces
:math:`V_{x} + V_{y}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the y-component of the shear stress within a 150x90x12 UA
section resulting from a shear force of 15 kN in the x-direction and 30 kN in the
y-direction::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vx=15e3, Vy=30e3)
stress_post.plot_stress_v_zy()
.. figure:: ../images/stress/stress_v_zy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zy_v)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_v_zxy(
self,
title=r"Stress Contour Plot - $\sigma_{zxy,\Sigma V}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the resultant shear stress
:math:`\sigma_{zxy,\Sigma V}` resulting from the sum of the applied shear forces
:math:`V_{x} + V_{y}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots a contour of the resultant shear stress within a 150x90x12 UA
section resulting from a shear force of 15 kN in the x-direction and 30 kN in the
y-direction::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vx=15e3, Vy=30e3)
stress_post.plot_stress_v_zxy()
.. figure:: ../images/stress/stress_v_zxy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zxy_v)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_vector_v_zxy(
self,
title=r"Stress Vector Plot - $\sigma_{zxy,\Sigma V}$",
cmap="YlOrBr",
normalize=False,
**kwargs,
):
r"""Produces a vector plot of the resultant shear stress
:math:`\sigma_{zxy,\Sigma V}` resulting from the sum of the applied shear forces
:math:`V_{x} + V_{y}`.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example generates a vector plot of the shear stress within a 150x90x12 UA
section resulting from a shear force of 15 kN in the x-direction and 30 kN in the
y-direction::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Vx=15e3, Vy=30e3)
stress_post.plot_vector_v_zxy()
.. figure:: ../images/stress/vector_v_zxy.png
:align: center
:scale: 75 %
Vector plot of the shear stress.
"""
sigxs = []
sigys = []
for group in self.material_groups:
sigxs.append(group.stress_result.sig_zx_v)
sigys.append(group.stress_result.sig_zy_v)
return self.plot_stress_vector(sigxs, sigys, title, cmap, normalize, **kwargs)
[docs] def plot_stress_zz(
self,
title=r"Stress Contour Plot - $\sigma_{zz}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the combined normal stress :math:`\sigma_{zz}` resulting from
all actions.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the normal stress within a 150x90x12 UA section resulting from
an axial force of 100 kN, a bending moment about the x-axis of 5 kN.m and a bending moment
about the y-axis of 2 kN.m::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(N=100e3, Mxx=5e6, Myy=2e6)
stress_post.plot_stress_zz()
.. figure:: ../images/stress/stress_zz.png
:align: center
:scale: 75 %
Contour plot of the normal stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zz)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_zx(
self,
title=r"Stress Contour Plot - $\sigma_{zx}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *x*-component of the shear stress :math:`\sigma_{zx}`
resulting from all actions.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the x-component of the shear stress within a 150x90x12 UA
section resulting from a torsion moment of 1 kN.m and a shear force of 30 kN in the
y-direction::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mzz=1e6, Vy=30e3)
stress_post.plot_stress_zx()
.. figure:: ../images/stress/stress_zx.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zx)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_zy(
self,
title=r"Stress Contour Plot - $\sigma_{zy}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the *y*-component of the shear stress :math:`\sigma_{zy}`
resulting from all actions.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the y-component of the shear stress within a 150x90x12 UA
section resulting from a torsion moment of 1 kN.m and a shear force of 30 kN in the
y-direction::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mzz=1e6, Vy=30e3)
stress_post.plot_stress_zy()
.. figure:: ../images/stress/stress_zy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zy)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_zxy(
self,
title=r"Stress Contour Plot - $\sigma_{zxy}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the resultant shear stress :math:`\sigma_{zxy}` resulting
from all actions.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots a contour of the resultant shear stress within a 150x90x12 UA
section resulting from a torsion moment of 1 kN.m and a shear force of 30 kN in the
y-direction::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mzz=1e6, Vy=30e3)
stress_post.plot_stress_zxy()
.. figure:: ../images/stress/stress_zxy.png
:align: center
:scale: 75 %
Contour plot of the shear stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_zxy)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_vector_zxy(
self,
title=r"Stress Vector Plot - $\sigma_{zxy}$",
cmap="YlOrBr",
normalize=False,
**kwargs,
):
r"""Produces a vector plot of the resultant shear stress :math:`\sigma_{zxy}` resulting
from all actions.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example generates a vector plot of the shear stress within a 150x90x12 UA
section resulting from a torsion moment of 1 kN.m and a shear force of 30 kN in the
y-direction::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(Mzz=1e6, Vy=30e3)
stress_post.plot_vector_zxy()
.. figure:: ../images/stress/vector_zxy.png
:align: center
:scale: 75 %
Vector plot of the shear stress.
"""
sigxs = []
sigys = []
for group in self.material_groups:
sigxs.append(group.stress_result.sig_zx)
sigys.append(group.stress_result.sig_zy)
return self.plot_stress_vector(sigxs, sigys, title, cmap, normalize, **kwargs)
[docs] def plot_stress_1(
self,
title=r"Stress Contour Plot - $\sigma_{1}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the major principal stress :math:`\sigma_{1}` resulting from all
actions.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots a contour of the major principal stress within a 150x90x12 UA
section resulting from the following actions:
* :math:`N = 50` kN
* :math:`M_{xx} = -5` kN.m
* :math:`M_{22} = 2.5` kN.m
* :math:`M_{zz} = 1.5` kN.m
* :math:`V_{x} = 10` kN
* :math:`V_{y} = 5` kN
::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
mesh = geometry.create_mesh(mesh_sizes=[2.5])
section = CrossSection(geometry, mesh)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(
N=50e3, Mxx=-5e6, M22=2.5e6, Mzz=0.5e6, Vx=10e3, Vy=5e3
)
stress_post.plot_stress_1()
.. figure:: ../images/stress/stress_1.png
:align: center
:scale: 50 %
Contour plot of the major principal stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_1)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_3(
self,
title=r"Stress Contour Plot - $\sigma_{3}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
r"""Produces a contour plot of the Minor principal stress :math:`\sigma_{3}` resulting from all
actions.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots a contour of the Minor principal stress within a 150x90x12 UA
section resulting from the following actions:
* :math:`N = 50` kN
* :math:`M_{xx} = -5` kN.m
* :math:`M_{22} = 2.5` kN.m
* :math:`M_{zz} = 1.5` kN.m
* :math:`V_{x} = 10` kN
* :math:`V_{y} = 5` kN
::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
mesh = geometry.create_mesh(mesh_sizes=[2.5])
section = CrossSection(geometry, mesh)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(
N=50e3, Mxx=-5e6, M22=2.5e6, Mzz=0.5e6, Vx=10e3, Vy=5e3
)
stress_post.plot_stress_3()
.. figure:: ../images/stress/stress_3.png
:align: center
:scale: 50 %
Contour plot of the minor principal stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_3)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_stress_vm(
self,
title=r"Stress Contour Plot - $\sigma_{vM}$",
cmap="coolwarm",
normalize=True,
**kwargs,
):
"""Produces a contour plot of the von Mises stress :math:`\sigma_{vM}` resulting from all
actions.
:param string title: Plot title
:param string cmap: Matplotlib color map.
:param bool normalize: If set to true, the CenteredNorm is used to scale the colormap.
If set to false, the default linear scaling is used.
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots a contour of the von Mises stress within a 150x90x12 UA section
resulting from the following actions:
* :math:`N = 50` kN
* :math:`M_{xx} = -5` kN.m
* :math:`M_{22} = 2.5` kN.m
* :math:`M_{zz} = 1.5` kN.m
* :math:`V_{x} = 10` kN
* :math:`V_{y} = 5` kN
::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
geometry.create_mesh(mesh_sizes=[20])
section = Section(geometry)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(
N=50e3, Mxx=-5e6, M22=2.5e6, Mzz=0.5e6, Vx=10e3, Vy=5e3
)
stress_post.plot_stress_vm()
.. figure:: ../images/stress/stress_vm.png
:align: center
:scale: 75 %
Contour plot of the von Mises stress.
"""
sigs = []
for group in self.material_groups:
sigs.append(group.stress_result.sig_vm)
return self.plot_stress_contour(sigs, title, cmap, normalize, **kwargs)
[docs] def plot_mohrs_circles(self, x, y, title=None, **kwargs):
r"""Plots Mohr's Circles of the 3D stress state at position x, y
:param float x: x-coordinate of the point to draw Mohr's Circle
:param float y: y-coordinate of the point to draw Mohr's Circle
:param string title: Plot title
:param kwargs: Passed to :func:`~sectionproperties.post.post.plotting_context`
:return: Matplotlib axes object
:rtype: :class:`matplotlib.axes`
The following example plots the Mohr's Circles for the 3D stress state within a 150x90x12
UA section resulting from the following actions:
* :math:`N = 50` kN
* :math:`M_{xx} = -5` kN.m
* :math:`M_{22} = 2.5` kN.m
* :math:`M_{zz} = 1.5` kN.m
* :math:`V_{x} = 10` kN
* :math:`V_{y} = 5` kN
at the point (10, 88.9).
::
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
geometry = steel_sections.angle_section(d=150, b=90, t=12, r_r=10, r_t=5, n_r=8)
mesh = geometry.create_mesh(mesh_sizes=[2.5])
section = Section(geometry, mesh)
section.calculate_geometric_properties()
section.calculate_warping_properties()
stress_post = section.calculate_stress(
N=50e3, Mxx=-5e6, M22=2.5e6, Mzz=0.5e6, Vx=10e3, Vy=5e3
)
stress_post.plot_mohrs_circles(10, 88.9)
.. figure:: ../images/stress/mohrs_circles.png
:align: center
:scale: 100 %
Mohr's Circles of the 3D stress state at (10, 88.9).
"""
pt = (x, y)
nodes = self.section.mesh_nodes
ele = self.section.mesh_elements
triang = tri.Triangulation(nodes[:, 0], nodes[:, 1], ele[:, 0:3])
# Find in which material group the point lies
pt_group = None
for group in self.material_groups:
mask_array = np.ones(len(self.section.elements), dtype=bool)
mask_array[group.el_ids] = False
triang.set_mask(mask_array)
trifinder = triang.get_trifinder()
if trifinder(*pt) != -1:
pt_group = group
triang.set_mask(None)
if pt_group is None:
raise Exception(f"Point {(*pt,)} is not within mesh")
# Assesmble the stress results from the relevant material group
sigma_zz_v = pt_group.stress_result.sig_zz
tau_xz_v = pt_group.stress_result.sig_zx
tau_yz_v = pt_group.stress_result.sig_zy
# Get the interpolators
sigma_zz_interp = tri.LinearTriInterpolator(triang, sigma_zz_v)
tau_xz_interp = tri.LinearTriInterpolator(triang, tau_xz_v)
tau_yz_interp = tri.LinearTriInterpolator(triang, tau_yz_v)
# Get the stresses at the point
sigma_zz = sigma_zz_interp(*pt).item()
tau_xz = tau_xz_interp(*pt).item()
tau_yz = tau_yz_interp(*pt).item()
# Assemble the stress tensor
sigma_xx = 0
sigma_yy = 0
tau_xy = 0
sigma = np.array(
[
[sigma_xx, tau_xy, tau_xz],
[tau_xy, sigma_yy, tau_yz],
[tau_xz, tau_yz, sigma_zz],
]
)
# Solve for the principal stresses using the general approach
s, n = np.linalg.eig(sigma)
sigma_3, sigma_2, sigma_1 = np.sort(s)
s_l = s.tolist()
sigma_1i, sigma_2i, sigma_3i = (
s_l.index(sigma_1),
s_l.index(sigma_2),
s_l.index(sigma_3),
)
n1, n2, n3 = n[:, sigma_1i], n[:, sigma_2i], n[:, sigma_3i]
# The tractions on each plane in cartesian coords wrt principal axes
n_inv = np.linalg.inv(n)
tractions = []
for col in range(3):
ss = n_inv[:, col].T @ np.diag(s) @ n_inv[:, col]
ts = np.sqrt(np.linalg.norm(np.diag(s) @ n_inv[:, col]) ** 2 - ss**2)
tractions.append((ss, ts))
def plot_circle(ax, c, R, col, label=None, fill=None):
circ = plt.Circle(c, R, fill=fill, ec=col, label=label)
ax.add_patch(circ)
ax.set_aspect(1)
ax.autoscale_view()
if title is None:
title = f"Mohr's Circles for 3D Stress State at {(*pt,)}"
# create plot and setup the plot
with post.plotting_context(title=title, **kwargs) as (fig, ax):
plot_circle(
ax,
(0.5 * (sigma_2 + sigma_3), 0),
0.5 * (sigma_2 - sigma_3),
"r",
r"C1: ($\sigma_2$, $\sigma_3$)",
)
plot_circle(
ax,
(0.5 * (sigma_1 + sigma_3), 0),
0.5 * (sigma_1 - sigma_3),
"b",
r"C2: ($\sigma_1$, $\sigma_3$)",
)
plot_circle(
ax,
(0.5 * (sigma_1 + sigma_2), 0),
0.5 * (sigma_1 - sigma_2),
"k",
r"C3: ($\sigma_1$, $\sigma_2$)",
)
for i, plane, col in zip(range(3), ["X", "Y", "Z"], ["r", "b", "k"]):
ax.plot(*tractions[i], f"{col}.", label=rf"{plane}-face")
ax.set_axisbelow(True)
ax.grid(which="both")
ax.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.set_ylabel(r"Shear stress $\tau$ (MPa)")
ax.set_xlabel(r"Direct stress $\sigma$ (MPa)")
ax.xaxis.set_tick_params(bottom=True, top=False, direction="inout")
ax.yaxis.set_tick_params(left=True, right=False, direction="inout")
#
# The following is just to get the labels positioned outside the axes
#
# Store default label positions
x_lbl_pos = ax.xaxis.label.get_position()
y_lbl_pos = ax.yaxis.label.get_position()
# Make spines pass through zero of the other axis
ax.spines["bottom"].set_position("zero")
ax.spines["left"].set_position("zero")
# Now set the coords
ax.xaxis.set_label_coords(*x_lbl_pos)
ax.yaxis.set_label_coords(*y_lbl_pos)
return ax
[docs]class MaterialGroup:
"""Class for storing elements of different materials.
A MaterialGroup object contains the finite element objects for a specified `material`. The
`stress_result` variable provides storage for stresses related each material.
:param material: Material object for the current MaterialGroup
:type material: :class:`~sectionproperties.pre.pre.Material`
:param int num_nods: Number of nodes for the entire cross-section
:cvar material: Material object for the current MaterialGroup
:vartype material: :class:`~sectionproperties.pre.pre.Material`
:cvar stress_result: A StressResult object for saving the stresses of the current material
:vartype stress_result: :class:`~sectionproperties.analysis.section.StressResult`
:cvar elements: A list of finite element objects that are of the current material type
:vartype elements: list[:class:`~sectionproperties.analysis.fea.Tri6`]
:cvar el_ids: A list of the element IDs of the elements that are of the current material type
:vartype el_ids: list[int]
"""
def __init__(self, material, num_nodes):
"""Inits the MaterialGroup class."""
self.material = material
self.stress_result = StressResult(num_nodes)
self.elements = []
self.el_ids = []
[docs] def add_element(self, element):
"""Adds an element and its element ID to the MaterialGroup.
:param element: Element to add to the MaterialGroup
:type element: :class:`~sectionproperties.analysis.fea.Tri6`
"""
# add Tri6 element to the list of elements
self.elements.append(element)
self.el_ids.append(element.el_id)
[docs]@dataclass
class StressResult:
"""Class for storing a stress result.
Provides variables to store the results from a cross-section stress analysis. Also provides a
method to calculate combined stresses.
:param int num_nodes: Number of nodes in the finite element mesh
:cvar sig_zz_n: Normal stress (:math:`\sigma_{zz,N}`) resulting from an axial force
:vartype sig_zz_n: :class:`numpy.ndarray`
:cvar sig_zz_mxx: Normal stress (:math:`\sigma_{zz,Mxx}`) resulting from a bending moment about
the xx-axis
:vartype sig_zz_mxx: :class:`numpy.ndarray`
:cvar sig_zz_myy: Normal stress (:math:`\sigma_{zz,Myy}`) resulting from a bending moment about
the yy-axis
:vartype sig_zz_myy: :class:`numpy.ndarray`
:cvar sig_zz_m11: Normal stress (:math:`\sigma_{zz,M11}`) resulting from a bending moment about
the 11-axis
:vartype sig_zz_m11: :class:`numpy.ndarray`
:cvar sig_zz_m22: Normal stress (:math:`\sigma_{zz,M22}`) resulting from a bending moment about
the 22-axis
:vartype sig_zz_m22: :class:`numpy.ndarray`
:cvar sig_zx_mzz: Shear stress (:math:`\sigma_{zx,Mzz}`) resulting from a torsion moment about
the zz-axis
:vartype sig_zx_mzz: :class:`numpy.ndarray`
:cvar sig_zy_mzz: Shear stress (:math:`\sigma_{zy,Mzz}`) resulting from a torsion moment about
the zz-axis
:vartype sig_zy_mzz: :class:`numpy.ndarray`
:cvar sig_zx_vx: Shear stress (:math:`\sigma_{zx,Vx}`) resulting from a shear force in the
x-direction
:vartype sig_zx_vx: :class:`numpy.ndarray`
:cvar sig_zy_vx: Shear stress (:math:`\sigma_{zy,Vx}`) resulting from a shear force in the
x-direction
:vartype sig_zy_vx: :class:`numpy.ndarray`
:cvar sig_zx_vy: Shear stress (:math:`\sigma_{zx,Vy}`) resulting from a shear force in the
y-direction
:vartype sig_zx_vy: :class:`numpy.ndarray`
:cvar sig_zy_vy: Shear stress (:math:`\sigma_{zy,Vy}`) resulting from a shear force in the
y-direction
:vartype sig_zy_vy: :class:`numpy.ndarray`
:cvar sig_zz_m: Normal stress (:math:`\sigma_{zz,\Sigma M}`) resulting from all bending moments
:vartype sig_zz_m: :class:`numpy.ndarray`
:cvar sig_zxy_mzz: Resultant shear stress (:math:`\sigma_{zxy,Mzz}`) resulting from a torsion
moment in the zz-direction
:vartype sig_zxy_mzz: :class:`numpy.ndarray`
:cvar sig_zxy_vx: Resultant shear stress (:math:`\sigma_{zxy,Vx}`) resulting from a a shear
force in the x-direction
:vartype sig_zxy_vx: :class:`numpy.ndarray`
:cvar sig_zxy_vy: Resultant shear stress (:math:`\sigma_{zxy,Vy}`) resulting from a a shear
force in the y-direction
:vartype sig_zxy_vy: :class:`numpy.ndarray`
:cvar sig_zx_v: Shear stress (:math:`\sigma_{zx,\Sigma V}`) resulting from all shear forces
:vartype sig_zx_v: :class:`numpy.ndarray`
:cvar sig_zy_v: Shear stress (:math:`\sigma_{zy,\Sigma V}`) resulting from all shear forces
:vartype sig_zy_v: :class:`numpy.ndarray`
:cvar sig_zxy_v: Resultant shear stress (:math:`\sigma_{zxy,\Sigma V}`) resulting from all
shear forces
:vartype sig_zxy_v: :class:`numpy.ndarray`
:cvar sig_zz: Combined normal force (:math:`\sigma_{zz}`) resulting from all actions
:vartype sig_zz: :class:`numpy.ndarray`
:cvar sig_zx: Combined shear stress (:math:`\sigma_{zx}`) resulting from all actions
:vartype sig_zx: :class:`numpy.ndarray`
:cvar sig_zy: Combined shear stress (:math:`\sigma_{zy}`) resulting from all actions
:vartype sig_zy: :class:`numpy.ndarray`
:cvar sig_zxy: Combined resultant shear stress (:math:`\sigma_{zxy}`) resulting from all
actions
:vartype sig_zxy: :class:`numpy.ndarray`
:cvar sig_1: Major principal stress (:math:`\sigma_{1}`) resulting from all actions
:vartype sig_1: :class:`numpy.ndarray`
:cvar sig_3: Minor principal stress (:math:`\sigma_{3}`) resulting from all actions
:vartype sig_3: :class:`numpy.ndarray`
:cvar sig_vm: von Mises stress (:math:`\sigma_{VM}`) resulting from all actions
:vartype sig_vm: :class:`numpy.ndarray`
"""
def __init__(self, num_nodes):
"""Inits the StressResult class."""
# allocate stresses arising directly from actions
self.sig_zz_n = np.zeros(num_nodes)
self.sig_zz_mxx = np.zeros(num_nodes)
self.sig_zz_myy = np.zeros(num_nodes)
self.sig_zz_m11 = np.zeros(num_nodes)
self.sig_zz_m22 = np.zeros(num_nodes)
self.sig_zx_mzz = np.zeros(num_nodes)
self.sig_zy_mzz = np.zeros(num_nodes)
self.sig_zx_vx = np.zeros(num_nodes)
self.sig_zy_vx = np.zeros(num_nodes)
self.sig_zx_vy = np.zeros(num_nodes)
self.sig_zy_vy = np.zeros(num_nodes)
# allocate combined stresses
self.sig_zz_m = np.zeros(num_nodes)
self.sig_zxy_mzz = np.zeros(num_nodes)
self.sig_zxy_vx = np.zeros(num_nodes)
self.sig_zxy_vy = np.zeros(num_nodes)
self.sig_zx_v = np.zeros(num_nodes)
self.sig_zy_v = np.zeros(num_nodes)
self.sig_zxy_v = np.zeros(num_nodes)
self.sig_zz = np.zeros(num_nodes)
self.sig_zx = np.zeros(num_nodes)
self.sig_zy = np.zeros(num_nodes)
self.sig_zxy = np.zeros(num_nodes)
self.sig_1 = np.zeros(num_nodes)
self.sig_3 = np.zeros(num_nodes)
self.sig_vm = np.zeros(num_nodes)
[docs] def calculate_combined_stresses(self):
"""Calculates the combined cross-section stresses."""
self.sig_zz_m = (
self.sig_zz_mxx + self.sig_zz_myy + self.sig_zz_m11 + self.sig_zz_m22
)
self.sig_zxy_mzz = (self.sig_zx_mzz**2 + self.sig_zy_mzz**2) ** 0.5
self.sig_zxy_vx = (self.sig_zx_vx**2 + self.sig_zy_vx**2) ** 0.5
self.sig_zxy_vy = (self.sig_zx_vy**2 + self.sig_zy_vy**2) ** 0.5
self.sig_zx_v = self.sig_zx_vx + self.sig_zx_vy
self.sig_zy_v = self.sig_zy_vx + self.sig_zy_vy
self.sig_zxy_v = (self.sig_zx_v**2 + self.sig_zy_v**2) ** 0.5
self.sig_zz = self.sig_zz_n + self.sig_zz_m
self.sig_zx = self.sig_zx_mzz + self.sig_zx_v
self.sig_zy = self.sig_zy_mzz + self.sig_zy_v
self.sig_zxy = (self.sig_zx**2 + self.sig_zy**2) ** 0.5
self.sig_1 = self.sig_zz / 2 + np.sqrt(
(self.sig_zz / 2) ** 2 + self.sig_zxy**2
)
self.sig_3 = self.sig_zz / 2 - np.sqrt(
(self.sig_zz / 2) ** 2 + self.sig_zxy**2
)
self.sig_vm = (self.sig_zz**2 + 3 * self.sig_zxy**2) ** 0.5
[docs]@dataclass
class SectionProperties:
"""Class for storing section properties.
Stores calculated section properties. Also provides methods to calculate section properties
entirely derived from other section properties.
:cvar float area: Cross-sectional area
:cvar float perimeter: Cross-sectional perimeter
:cvar float mass: Cross-sectional mass
:cvar float ea: Modulus weighted area (axial rigidity)
:cvar float ga: Modulus weighted product of shear modulus and area
:cvar float nu_eff: Effective Poisson's ratio
:cvar float e_eff: Effective elastic modulus
:cvar float g_eff: Effective shear modulus
:cvar float qx: First moment of area about the x-axis
:cvar float qy: First moment of area about the y-axis
:cvar float ixx_g: Second moment of area about the global x-axis
:cvar float iyy_g: Second moment of area about the global y-axis
:cvar float ixy_g: Second moment of area about the global xy-axis
:cvar float cx: X coordinate of the elastic centroid
:cvar float cy: Y coordinate of the elastic centroid
:cvar float ixx_c: Second moment of area about the centroidal x-axis
:cvar float iyy_c: Second moment of area about the centroidal y-axis
:cvar float ixy_c: Second moment of area about the centroidal xy-axis
:cvar float zxx_plus: Section modulus about the centroidal x-axis for stresses at the positive
extreme value of y
:cvar float zxx_minus: Section modulus about the centroidal x-axis for stresses at the negative
extreme value of y
:cvar float zyy_plus: Section modulus about the centroidal y-axis for stresses at the positive
extreme value of x
:cvar float zyy_minus: Section modulus about the centroidal y-axis for stresses at the negative
extreme value of x
:cvar float rx_c: Radius of gyration about the centroidal x-axis.
:cvar float ry_c: Radius of gyration about the centroidal y-axis.
:cvar float i11_c: Second moment of area about the centroidal 11-axis
:cvar float i22_c: Second moment of area about the centroidal 22-axis
:cvar float phi: Principal axis angle
:cvar float z11_plus: Section modulus about the principal 11-axis for stresses at the positive
extreme value of the 22-axis
:cvar float z11_minus: Section modulus about the principal 11-axis for stresses at the negative
extreme value of the 22-axis
:cvar float z22_plus: Section modulus about the principal 22-axis for stresses at the positive
extreme value of the 11-axis
:cvar float z22_minus: Section modulus about the principal 22-axis for stresses at the negative
extreme value of the 11-axis
:cvar float r11_c: Radius of gyration about the principal 11-axis.
:cvar float r22_c: Radius of gyration about the principal 22-axis.
:cvar float j: Torsion constant
:cvar omega: Warping function
:vartype omega: :class:`numpy.ndarray`
:cvar psi_shear: Psi shear function
:vartype psi_shear: :class:`numpy.ndarray`
:cvar phi_shear: Phi shear function
:vartype phi_shear: :class:`numpy.ndarray`
:cvar float Delta_s: Shear factor
:cvar float x_se: X coordinate of the shear centre (elasticity approach)
:cvar float y_se: Y coordinate of the shear centre (elasticity approach)
:cvar float x11_se: 11 coordinate of the shear centre (elasticity approach)
:cvar float y22_se: 22 coordinate of the shear centre (elasticity approach)
:cvar float x_st: X coordinate of the shear centre (Trefftz's approach)
:cvar float y_st: Y coordinate of the shear centre (Trefftz's approach)
:cvar float gamma: Warping constant
:cvar float A_sx: Shear area about the x-axis
:cvar float A_sy: Shear area about the y-axis
:cvar float A_sxy: Shear area about the xy-axis
:cvar float A_s11: Shear area about the 11 bending axis
:cvar float A_s22: Shear area about the 22 bending axis
:cvar float beta_x_plus: Monosymmetry constant for bending about the x-axis with the top flange
in compression
:cvar float beta_x_minus: Monosymmetry constant for bending about the x-axis with the bottom
flange in compression
:cvar float beta_y_plus: Monosymmetry constant for bending about the y-axis with the top flange
in compression
:cvar float beta_y_minus: Monosymmetry constant for bending about the y-axis with the bottom
flange in compression
:cvar float beta_11_plus: Monosymmetry constant for bending about the 11-axis with the top
flange in compression
:cvar float beta_11_minus: Monosymmetry constant for bending about the 11-axis with the bottom
flange in compression
:cvar float beta_22_plus: Monosymmetry constant for bending about the 22-axis with the top
flange in compression
:cvar float beta_22_minus: Monosymmetry constant for bending about the 22-axis with the bottom
flange in compression
:cvar float x_pc: X coordinate of the global plastic centroid
:cvar float y_pc: Y coordinate of the global plastic centroid
:cvar float x11_pc: 11 coordinate of the principal plastic centroid
:cvar float y22_pc: 22 coordinate of the principal plastic centroid
:cvar float sxx: Plastic section modulus about the centroidal x-axis
:cvar float syy: Plastic section modulus about the centroidal y-axis
:cvar float sf_xx_plus: Shape factor for bending about the x-axis with respect to the top fibre
:cvar float sf_xx_minus: Shape factor for bending about the x-axis with respect to the bottom
fibre
:cvar float sf_yy_plus: Shape factor for bending about the y-axis with respect to the top fibre
:cvar float sf_yy_minus: Shape factor for bending about the y-axis with respect to the bottom
fibre
:cvar float s11: Plastic section modulus about the 11-axis
:cvar float s22: Plastic section modulus about the 22-axis
:cvar float sf_11_plus: Shape factor for bending about the 11-axis with respect to the top
fibre
:cvar float sf_11_minus: Shape factor for bending about the 11-axis with respect to the bottom
fibre
:cvar float sf_22_plus: Shape factor for bending about the 22-axis with respect to the top
fibre
:cvar float sf_22_minus: Shape factor for bending about the 22-axis with respect to the bottom
fibre
"""
area: Optional[float] = None
perimeter: Optional[float] = None
mass: Optional[float] = None
ea: Optional[float] = None
ga: Optional[float] = None
nu_eff: Optional[float] = None
e_eff: Optional[float] = None
g_eff: Optional[float] = None
qx: Optional[float] = None
qy: Optional[float] = None
ixx_g: Optional[float] = None
iyy_g: Optional[float] = None
ixy_g: Optional[float] = None
cx: Optional[float] = None
cy: Optional[float] = None
ixx_c: Optional[float] = None
iyy_c: Optional[float] = None
ixy_c: Optional[float] = None
zxx_plus: Optional[float] = None
zxx_minus: Optional[float] = None
zyy_plus: Optional[float] = None
zyy_minus: Optional[float] = None
rx_c: Optional[float] = None
ry_c: Optional[float] = None
i11_c: Optional[float] = None
i22_c: Optional[float] = None
phi: Optional[float] = None
z11_plus: Optional[float] = None
z11_minus: Optional[float] = None
z22_plus: Optional[float] = None
z22_minus: Optional[float] = None
r11_c: Optional[float] = None
r22_c: Optional[float] = None
j: Optional[float] = None
omega: Optional[np.ndarray] = None
psi_shear: Optional[np.ndarray] = None
phi_shear: Optional[np.ndarray] = None
Delta_s: Optional[float] = None
x_se: Optional[float] = None
y_se: Optional[float] = None
x11_se: Optional[float] = None
y22_se: Optional[float] = None
x_st: Optional[float] = None
y_st: Optional[float] = None
gamma: Optional[float] = None
A_sx: Optional[float] = None
A_sy: Optional[float] = None
A_sxy: Optional[float] = None
A_s11: Optional[float] = None
A_s22: Optional[float] = None
beta_x_plus: Optional[float] = None
beta_x_minus: Optional[float] = None
beta_y_plus: Optional[float] = None
beta_y_minus: Optional[float] = None
beta_11_plus: Optional[float] = None
beta_11_minus: Optional[float] = None
beta_22_plus: Optional[float] = None
beta_22_minus: Optional[float] = None
x_pc: Optional[float] = None
y_pc: Optional[float] = None
x11_pc: Optional[float] = None
y22_pc: Optional[float] = None
sxx: Optional[float] = None
syy: Optional[float] = None
sf_xx_plus: Optional[float] = None
sf_xx_minus: Optional[float] = None
sf_yy_plus: Optional[float] = None
sf_yy_minus: Optional[float] = None
s11: Optional[float] = None
s22: Optional[float] = None
sf_11_plus: Optional[float] = None
sf_11_minus: Optional[float] = None
sf_22_plus: Optional[float] = None
sf_22_minus: Optional[float] = None
[docs] def asdict(self):
"""
Returns the SectionProperties dataclass object as a dictionary.
"""
return asdict(self)
[docs] def calculate_elastic_centroid(self):
"""Calculates the elastic centroid based on the cross-section area and first moments of
area.
"""
self.cx = self.qy / self.ea
self.cy = self.qx / self.ea
[docs] def calculate_centroidal_properties(self, mesh):
"""Calculates the geometric section properties about the centroidal and principal axes
based on the results about the global axis.
"""
# calculate second moments of area about the centroidal xy axis
self.ixx_c = self.ixx_g - self.qx**2 / self.ea
self.iyy_c = self.iyy_g - self.qy**2 / self.ea
self.ixy_c = self.ixy_g - self.qx * self.qy / self.ea
# calculate section moduli about the centroidal xy axis
nodes = np.array(mesh["vertices"])
xmax = nodes[:, 0].max()
xmin = nodes[:, 0].min()
ymax = nodes[:, 1].max()
ymin = nodes[:, 1].min()
self.zxx_plus = self.ixx_c / abs(ymax - self.cy)
self.zxx_minus = self.ixx_c / abs(ymin - self.cy)
self.zyy_plus = self.iyy_c / abs(xmax - self.cx)
self.zyy_minus = self.iyy_c / abs(xmin - self.cx)
# calculate radii of gyration about centroidal xy axis
self.rx_c = (self.ixx_c / self.ea) ** 0.5
self.ry_c = (self.iyy_c / self.ea) ** 0.5
# calculate principal 2nd moments of area about the centroidal xy axis
Delta = (((self.ixx_c - self.iyy_c) / 2) ** 2 + self.ixy_c**2) ** 0.5
self.i11_c = (self.ixx_c + self.iyy_c) / 2 + Delta
self.i22_c = (self.ixx_c + self.iyy_c) / 2 - Delta
# calculate initial principal axis angle
if abs(self.ixx_c - self.i11_c) < 1e-12 * self.i11_c:
self.phi = 0
else:
self.phi = np.arctan2(self.ixx_c - self.i11_c, self.ixy_c) * 180 / np.pi
# calculate section moduli about the principal axis
for (i, pt) in enumerate(nodes):
x = pt[0] - self.cx
y = pt[1] - self.cy
# determine the coordinate of the point wrt the principal axis
(x1, y2) = fea.principal_coordinate(self.phi, x, y)
# initialise min, max variables
if i == 0:
x1max = x1
x1min = x1
y2max = y2
y2min = y2
# update the mins and maxes where necessary
x1max = max(x1max, x1)
x1min = min(x1min, x1)
y2max = max(y2max, y2)
y2min = min(y2min, y2)
# evaluate principal section moduli
self.z11_plus = self.i11_c / abs(y2max)
self.z11_minus = self.i11_c / abs(y2min)
self.z22_plus = self.i22_c / abs(x1max)
self.z22_minus = self.i22_c / abs(x1min)
# calculate radii of gyration about centroidal principal axis
self.r11_c = (self.i11_c / self.ea) ** 0.5
self.r22_c = (self.i22_c / self.ea) ** 0.5