# -*- coding: utf-8 -*-
"""
Module for the structural design of steel members.
A structure of classes and functions are provided to assist the basic design
checks for steel structures, based on EN1993-1-1.
The classes included attempt to categorise various calculations into logical
entities like 'part' 'cross-section' etc.
All the design formulas are defined as standalone functions outside classes so
that they can be easily called independently. The same functions are used
within the classes whenever needed.
"""
import logging
import numpy as np
import sys
import pdb
logging.basicConfig(
level=logging.INFO,
filename='logging.log',
filemode='w',
#format='%(name)s - %(levelname)s - %(message)s'
)
[docs]class Geometry:
"""
Structural element geometry.
Class for the geometric properties of a structural element.
Parameters
----------
cs_sketch : CsSketch object
Cross-section sketch.
length : float
Member's length.
"""
[docs] def __init__(self, cs_sketch, length, thickness):
self.cs_sketch = cs_sketch
self.length = length
self.thickness = thickness
[docs]class CsSketch:
"""
Cross-section geometry.
Parameters
----------
nodes : list
List of points.
elem : list
Element connectivity.
"""
[docs] def __init__(self, nodes, elem):
self.nodes = nodes
self.elem = elem
[docs]class CsProps:
"""
Cross-section properties.
Class for the mass properties of cross-sections. The properties can be calculated using the from_cs_sketch() method.
Parameters
----------
area : float
Cross-sectional area.
xc : float
`x` coordinate of the gravity center.
yc : float
`y` coordinate of the gravity center.
moi_xx : float
Moment of inertia around `x` axis.
moi_yy : float
Moment of inertia around `y` axis.
moi_xy : float
Polar moment of inertia.
theta_principal : float
Rotation of the principal axes.
moi_1 : float
Moment of inertia around the major axis.
moi_2 : float
Moment of inertia around the minor axis.
"""
[docs] def __init__(self,
area=None,
a_eff=None,
xc=None,
yc=None,
moi_xx=None,
moi_yy=None,
moi_xy=None,
theta_principal=None,
moi_1=None,
moi_2=None
):
self.area = area
self.a_eff = a_eff
self.xc = xc
self.yc = yc
self.moi_xx = moi_xx
self.moi_yy = moi_yy
self.moi_xy = moi_xy
self.theta_principal = theta_principal
self.moi_1 = moi_1
self.moi_2 = moi_2
[docs] @classmethod
def from_cs_sketch(cls, cs_sketch):
"""
Cross-section calculator.
Alternative constructor, calculates mass properties of a given sc sketch and returns a CsProps object.
Parameters
----------
cs_sketch : CsSketch object
Notes
-----
"""
nele = len(cs_sketch.elem[0])
node = cs_sketch.elem[0] + cs_sketch.elem[1]
nnode = 0
j = 0
while node:
i = [ii for ii, x in enumerate(node) if x == node[0]]
for ii in sorted(i, reverse=True):
del node[ii]
if len(i) == 2:
j += 1
nnode += 1
# classify the section type (currently not used)
# if j == nele:
# section = 'close' # single cell
# elif j == nele - 1:
# section = 'open' # singly-branched
# else:
# section = 'open' # multi-branched
# Calculate the cs-properties
tt = []
xm = []
ym = []
xd = []
yd = []
side_length = []
for i in range(nele):
sn = cs_sketch.elem[0][i]
fn = cs_sketch.elem[1][i]
# thickness of the element
tt = tt + [cs_sketch.elem[2][i]]
# compute the coordinate of the mid point of the element
xm = xm + [mean_list([cs_sketch.nodes[0][sn], cs_sketch.nodes[0][fn]])]
ym = ym + [mean_list([cs_sketch.nodes[1][sn], cs_sketch.nodes[1][fn]])]
# compute the dimension of the element
xd = xd + [(cs_sketch.nodes[0][fn] - cs_sketch.nodes[0][sn])]
yd = yd + [(cs_sketch.nodes[1][fn] - cs_sketch.nodes[1][sn])]
# compute the length of the element
side_length = side_length + [np.sqrt(xd[i] ** 2 + yd[i] ** 2)]
# calculate cross sectional area
area = sum([a * b for a, b in zip(side_length, tt)])
# compute the centroid
xc = sum([a * b * c for a, b, c in zip(side_length, tt, xm)]) / area
yc = sum([a * b * c for a, b, c in zip(side_length, tt, ym)]) / area
if abs(xc / np.sqrt(area)) < 1e-12:
xc = 0
if abs(yc / np.sqrt(area)) < 1e-12:
yc = 0
# Calculate MOI
moi_xx = sum([sum(a) for a in zip([a ** 2 * b * c / 12 for a, b, c in zip(yd, side_length, tt)],
[(a - yc) ** 2 * b * c for a, b, c in
zip(ym, side_length, tt)])])
moi_yy = sum([sum(a) for a in zip([a ** 2 * b * c / 12 for a, b, c in zip(xd, side_length, tt)],
[(a - xc) ** 2 * b * c for a, b, c in
zip(xm, side_length, tt)])])
moi_xy = sum(
[sum(a) for a in zip([a * b * c * d / 12 for a, b, c, d in zip(xd, yd, side_length, tt)],
[(a - xc) * (b - yc) * c * d for a, b, c, d in
zip(xm, ym, side_length, tt)])])
if abs(moi_xy / area ** 2) < 1e-12:
moi_xy = 0
# Calculate angle of principal axes
if moi_xx == moi_yy:
theta_principal = np.pi / 2
else:
theta_principal = np.arctan(
(-2 * moi_xy) / (moi_xx - moi_yy)) / 2
# Change to centroid principal coordinates
# coord12 = [[a - xc for a in cs_sketch.nodes[0]],
# [a - yc for a in cs_sketch.nodes[1]]]
coord12 = np.array([[np.cos(theta_principal), np.sin(theta_principal)],
[-np.sin(theta_principal), np.cos(theta_principal)]]).dot(cs_sketch.nodes)
# re-calculate cross sectional properties for the centroid
for i in range(nele):
sn = cs_sketch.elem[0][i]
fn = cs_sketch.elem[1][i]
# calculate the coordinate of the mid point of the element
xm = xm + [mean_list([coord12[0][sn], coord12[0][fn]])]
ym = ym + [mean_list([coord12[1][sn], coord12[1][fn]])]
# calculate the dimension of the element
xd = xd + [(coord12[0][fn] - coord12[0][sn])]
yd = yd + [(coord12[1][fn] - coord12[1][sn])]
# calculate the principal moment of inertia
moi_1 = sum([sum(a) for a in zip([a ** 2 * b * c / 12 for a, b, c in zip(yd, side_length, tt)],
[(a - yc) ** 2 * b * c for a, b, c in
zip(ym, side_length, tt)])])
moi_2 = sum([sum(a) for a in zip([a ** 2 * b * c / 12 for a, b, c in zip(xd, side_length, tt)],
[(a - xc) ** 2 * b * c for a, b, c in
zip(xm, side_length, tt)])])
return cls(
area=area,
xc=xc,
yc=yc,
moi_xx=moi_xx,
moi_yy=moi_yy,
moi_xy=moi_xy,
theta_principal=theta_principal,
moi_1=moi_1,
moi_2=moi_2
)
[docs]class CsRHS(CsProps):
"""
"""
[docs] def __init__(
self,
width=None,
height=None,
thick=None,
b_flat=None,
h_flat=None,
hcflat=None,
r_in=None,
r_out=None,
area=None,
a_eff=None,
xc=None,
yc=None,
moi_xx=None,
moi_yy=None,
moi_xy=None,
theta_principal=None,
moi_1=None,
moi_2=None,
w_pl=None,
w_el=None,
n_prc=None,
psi_cs=None,
negative_moment=None,
alfa_flange=None,
psi_flange=None,
alfa_web=None,
psi_web=None,
material=None,
n_rd=None,
m_n_rd=None,
m_pl_rd=None,
m_el_rd=None,
m_rd=None,
m_prc=None,
cs_class=None,
nd_width_web=None,
nd_width_flange=None
):
if sys.version_info[0] < 3:
CsProps.__init__(
self,
area=area,
a_eff=a_eff,
xc=xc,
yc=yc,
moi_xx=moi_xx,
moi_yy=moi_yy,
moi_xy=moi_xy,
theta_principal=theta_principal,
moi_1=moi_1,
moi_2=moi_2
)
else:
super().__init__(
area=area,
a_eff=a_eff,
xc=xc,
yc=yc,
moi_xx=moi_xx,
moi_yy=moi_yy,
moi_xy=moi_xy,
theta_principal=theta_principal,
moi_1=moi_1,
moi_2=moi_2
)
self.width = width
self.height = height
self.thick = thick
self.b_flat = b_flat
self.h_flat = h_flat
self.hcflat = hcflat
self.r_in = r_in
self.r_out = r_out
self.material = material
self.w_pl = w_pl
self.w_el = w_el,
self.n_prc = n_prc
self.psi_cs = psi_cs
self.negative_moment=negative_moment
self.alfa_flange = alfa_flange
self.psi_flange = psi_flange
self.alfa_web = alfa_web
self.psi_web = psi_web
self.cs_class = cs_class
self.nd_width_web = nd_width_web
self.nd_width_flange = nd_width_flange
# I add also resistance values that are calculated during the
# classification process but they could be removed if/when I make a
# separate Part for RHS members which will carry the resistance
# information instead.
# Additional comment: Maybe cross-sectional resistances should indeed
# be in the cs_props object (in this case CsRHS). In that case, ignore
# the comment above.
self.n_rd = n_rd
self.m_n_rd = m_n_rd
self.m_pl_rd = m_pl_rd
self.m_el_rd = m_el_rd
self.m_rd = m_rd
self.m_prc = m_prc
[docs] @classmethod
def from_compression(
cls,
width,
height,
thick,
material=None,
force_gross=False
):
"""
Create a CsRHS for a profile under pure compression.
"""
# Useful vars in shorter names
epsilon = material.epsilon
f_yield = material.f_y_nominal
# Calculate the flat parts and the corner radius of the profile.
h_flat, b_flat, r_out, r_in = cls.calc_flat_width_radius(
height,
width,
thick
)
# Calculate the cross-sectional area
area = cls.calc_area(height, width, thick)
logging.info("Area, A: %.3f" % area)
# Distance from the profile centre to the mid-line of the upper flange.
h2um = (height - thick)/2
# Calculation of the effective resistance under pure compression. This
# is needed in case the cs is effective (class 4) under the combined
# axial/bending in order to calculate the applied load N_Ed =
# n_prc*N_eff,Rd
logging.info("PURE COMPRESSION: calculating auxiliary profile")
cs_pc = cls(
width=width,
height=height,
thick=thick,
b_flat=b_flat,
h_flat=h_flat,
hcflat=h_flat,
r_in=r_in,
r_out=r_out,
area=area,
a_eff=area,
xc=height/2,
yc=width/2,
n_prc=-1,
alfa_flange=1.,
psi_flange=1.,
alfa_web=1,
psi_web=1.,
material=material,
)
logging.info(
"Perform classification for pure-compression"
)
(
flange_class_pc,
b_eff_u_pc,
b_ineff_u_pc,
b_eff_d_pc,
b_ineff_d_pc,
web_class_pc,
h_eff_pc,
h_ineff_pc,
alfa_web_pc,
psi_web_pc,
psi_pc,
cs_class_pc,
) = cs_pc.cs_classification(force_gross=force_gross)
logging.info(
"b_eff_u_pc, b_eff_d_pc, h_eff_pc: %.3f, %.3f, %.3f" % (
b_eff_u_pc,
b_eff_d_pc,
h_eff_pc
)
)
# The effective web area is evenly distributed up and down for the case
# of pure compression.
he1_pc = 0.5*h_eff_pc
he2_pc = 0.5*h_eff_pc
logging.info("he1_pc, he2_pc: %.3f, %.3f" % (he1_pc, he2_pc))
# Calculate the effective area for the case of pure compression
a_eff_pc = area - thick*(b_ineff_u_pc + b_ineff_d_pc + 2*h_ineff_pc)
logging.info("A_{eff,pc}, a_eff_pc: %.3f" % a_eff_pc)
# Calculation of the applied axial load, N_{Ed}, based on the requested
# n_prc. N_{Ed} is calculated as n_prc*N_{Rd}, where N_{Rd} regards to
# the full cross-section or the effective, depending on the class of
# the cross-section under pure compression.
n_rd_pc = a_eff_pc*f_yield/gamma_m(0, country="NO", nominal=False)
logging.info("N_Rd for pure compression, n_rd_pc: %.3f" % n_rd_pc)
logging.info("PURE COMPRESSION: End of auxiliary calculations.")
# Populate the joint object with all the remaining calculated
# properties before returning.
cs_pc.a_eff = a_eff_pc
cs_pc.alfa_web = alfa_web_pc
cs_pc.cs_class = cs_class_pc
cs_pc.hcflat = h_flat
cs_pc.m_el_rd = 0
cs_pc.m_n_rd = 0
cs_pc.psi_cs = 1
cs_pc.psi_web = 1
cs_pc.r_in = r_in
cs_pc.m_prc = 0
cs_pc.m_pl_rd = 0
cs_pc.n_rd = n_rd_pc
return(cs_pc)
[docs] @classmethod
def from_bending(
cls,
width,
height,
thick,
material=None,
negative_moment=None,
force_gross=False
):
# Useful vars in shorter names
epsilon = material.epsilon
f_yield = material.f_y_nominal
# Calculate the flat parts and the corner radius of the profile.
h_flat, b_flat, r_out, r_in = cls.calc_flat_width_radius(
height,
width,
thick
)
# Calculate the cross-sectional area
area = cls.calc_area(height, width, thick)
logging.info("Area, A: %.3f" % area)
# Distance from the profile centre to the mid-line of the upper flange.
h2um = (height - thick)/2
logging.info("PURE BENDING: calculating auxiliary profile")
cs_pb = cls(
width=width,
height=height,
thick=thick,
b_flat=b_flat,
h_flat=h_flat,
hcflat=h_flat/2,
r_in=r_in,
r_out=r_out,
area=area,
#a_eff=area,
xc=height/2,
yc=width/2,
n_prc=0,
#TODO: Cross check that the negative_moment makes no difference
#whether it is True or False.
negative_moment=negative_moment,
alfa_flange=1.,
psi_flange=1.,
alfa_web=0.5,
psi_web=-1.,
material=material,
)
logging.info(
"Perform classification for pure-bending"
)
(
flange_class_pb,
b_eff_u_pb,
b_ineff_u_pb,
b_eff_d_pb,
b_ineff_d_pb,
web_class_pb,
h_eff_pb,
h_ineff_pb,
alfa_web_pb,
psi_web_pb,
psi_pb,
cs_class_pb,
) = cs_pb.cs_classification(force_gross=force_gross)
logging.info(
"b_eff_u_pb, b_eff_d_pb, h_eff_pb: %.3f, %.3f, %.3f" % (
b_eff_u_pb,
b_eff_d_pb,
h_eff_pb
)
)
# Proceed with the moment resistance calculation for the cases of
# elastic or plastic response for the pure bending case.
# Plastic:
w_pl_pb = cls.calc_w_pl(
h_flat,
b_flat,
thick,
h2um,
r_in,
r_out,
)
m_pl_rd_pb = f_yield*w_pl_pb/gamma_m(0, country="NO", nominal=False)
# Elastic:
# The effective web area is distributed 40% 60% of the compression
# zone for the case of pure bending.
he1_pb = 0.4*h_eff_pb
he2_pb = 0.6*h_eff_pb
logging.info("he1_pb, he2_pb: %.3f, %.3f" % (he1_pb, he2_pb))
## Calculate the effective area for the case of pure bending.
#a_eff_pb = area - thick*(b_ineff_u_pb + b_ineff_d_pb + 2*h_ineff_pb)
#logging.info("A_{eff,pb}, a_eff_pc: %.3f" % a_eff_pb)
# Calculate cog_y for the effective cross-section (if any)
logging.info("Calculate cog_y for pure bending")
cog_y_pb = cls.calc_cog_y(
area,
thick,
h_flat/2,
he2_pb,
h_ineff_pb,
h_flat,
h2um,
b_ineff_u_pb,
b_ineff_d_pb
)
logging.info("cog_y_pb: %.3f" % (cog_y_pb))
# Calculate the effective moment of inertia
i_total_pb = cls.rhs_moi(
b_flat,
b_eff_u_pb,
b_eff_d_pb,
h_flat,
cog_y_pb,
he1_pb,
he2_pb,
h_flat/2,
thick
)
logging.info("I_pb: %.3f" % (i_total_pb))
# The extreme fiber of the cs
h_xtr_pb = height/2. + abs(cog_y_pb)
logging.info("h_xtr_pb: %.3f" % (h_xtr_pb))
# Calculate the section modulus for the case of pure bending.
w_el_pb = i_total_pb/h_xtr_pb
logging.info("w_el_pb: %.3f" % (w_el_pb))
# Finally, calculate the elastic moment resistance for the case of pure
# bending
m_el_rd_pb = f_yield*w_el_pb/gamma_m(0, country="NO", nominal=False)
# The cross section resistance under pure bending
if cs_class_pb > 2:
m_rd_pb = m_el_rd_pb
else:
m_rd_pb = m_pl_rd_pb
logging.info("PURE BENDING: End of auxiliary calculations.")
a_eff_pb = area - thick*(2*h_ineff_pb + b_ineff_u_pb + b_ineff_d_pb)
# Populate the joint object with all the remaining calculated
# properties before returning.
cs_pb.a_eff = a_eff_pb
cs_pb.alfa_web = alfa_web_pb
cs_pb.cs_class = cs_class_pb
cs_pb.hcflat = h_flat/2
cs_pb.m_el_rd = m_el_rd_pb
cs_pb.m_n_rd = 1
cs_pb.psi_cs = -1
cs_pb.psi_web = -1
cs_pb.r_in = r_in
cs_pb.m_prc = 1
cs_pb.m_pl_rd = m_pl_rd_pb
cs_pb.n_rd = 0
cs_pb.m_rd = m_rd_pb
cs_pb.w_pl = w_pl_pb
cs_pb.w_el = w_el_pb
return(cs_pb)
[docs] @classmethod
def from_n_prc(
cls,
width,
height,
thick,
material=None,
n_prc=None,
negative_moment=None,
force_gross=False
):
"""
Alternative constructor for CsRHS.
This method should be used for creating new RHS profile objects instead
of defining them directly through __init__. The overall geometric
properties are given (width, height, thickness), the material and the
loading condition. The loading is described by the utilisation ratio of
the axial resistance, n_prc. If the profile is not fully utilised by
axial loads (n_prc either 1 or -1), an additional positive moment is
considered to full utilisation.
Axial utilisation, n_prc, is positive for tension (N_Ed = N_prc*N_t,Rd)
and negative for compression (N_Ed = N_prc*N_c,Rd). Positive bending
moment generated sigma_min (compression) at the upper flange.
Parameters
----------
width : float
height : float
thick : float
material " :obj:`Material`, optional
By default, a generic steel S355 is used, with f_u=510 MPa. The
default is currently disabled so if no metarial is given it
crashes.
n_prc : float in [-1, 1], optional
Axial utilisation ratio. Negative for compression. By default, a
profile in pure compression is considered, n_prc=-1.
negative_moment : bool, optional
Defines the sign of the moment when abs(n_prc) < 1. Default is False
which implies sigma_min (compression) in the upper flange.
force_gross : bool, optional
Treat class 4 cross sections as a class 3. The argument is passed
to the "cs_classification" method. Default is false.
"""
logging.info(
"""
From here, starts the calculation of an RHS cross-section. The
requested input is:
width = %.3f
height = %.3f
thick = %.3f
n_prc = %.3f
""" % (width, height, thick, n_prc)
)
if negative_moment is None:
negative_moment = False
# Useful vars in shorter names
epsilon = material.epsilon
f_yield = material.f_y_nominal
# Calculate the flat parts and the corner radius of the profile.
h_flat, b_flat, r_out, r_in = cls.calc_flat_width_radius(
height,
width,
thick
)
# Calculate the cross-sectional area
area = cls.calc_area(height, width, thick)
logging.info("Area, A: %.3f" % area)
# Distance from the profile centre to the mid-line of the upper flange.
h2um = (height - thick)/2
# Calculation of the effective resistance under pure compression. This
# is needed in case the cs is effective (class 4) under the combined
# axial/bending in order to calculate the applied load N_Ed =
# n_prc*N_eff,Rd
cs_pc = cls.from_compression(
width,
height,
thick,
material=material,
force_gross=force_gross
)
cs_pb = cls.from_bending(
width,
height,
thick,
material=material,
negative_moment=negative_moment,
force_gross=force_gross
)
# The applied axial force is calculated by the requested utilisation
# based on the pure compression resistance.
n_ed = n_prc*cs_pc.n_rd
logging.info("N_{Ed}, n_prc * N_rd for pure compression: %.3f" % n_ed)
# The plastic section modulus is identical to the pure-bending case,
# and it is already calculated.
w_pl = cs_pb.w_pl
logging.info("W_pl: %.3f" % w_pl)
# Plastic moment resistance
m_pl_rd = w_pl*f_yield/gamma_m(0, country="NO", nominal=False)
logging.info("M_pl_Rd: %.3f" % m_pl_rd)
# Initial assumption of cog_y (Centre Of Gravity) on symmetry and pure
# bending and gross cross-section
cog_y = 0
a_eff = area
b_eff_u = b_flat
b_eff_d = b_flat
he1 = 0.0*h_flat
he2 = 0.5*h_flat
htflat = 0.5*h_flat
psi_web = -1.
# Recalculate n_prc. The requested n_prc regards to a pure compression
# case, which might be for the gross cross-section or for an
# effective one, as it is calculated for pure compression. Through the
# iterations, the area is different. Might be gross or effective for
# moment + axial, so the n_prc needs to be recalculated for a fixed
# N_Ed on every iteration. The first calculation is assuming plastic
# response, which if true, the area=a_eff will not change so no loops
# will need to be performed (that's why the calculation is out of the
# loop)
n_rd = a_eff*f_yield/gamma_m(0, country="NO", nominal=False)
n_prc_current = n_ed/n_rd
logging.info("n_prc of the current profile: %.3f" % n_prc_current)
approx = 1.
logging.info("\n ITERATIONS START HERE")
while approx>0.00000001:
# Distance of the new CoG to the extreme fiber
# NOTE: This calculation will not know which direction is the
# h_xtr, to the upper or the lower flange. If needed, I can change
# this with an "if negative_moment" condition.
h_xtr = height/2. + abs(cog_y)
logging.info("xtrm fiber: %.3f" % h_xtr)
# Calculate moment of inertia
logging.info(
" b_flat, b_eff_u, b_eff_d, h_flat, cog_y, he1, he2, htflat,"
"thick: %.3f, %.3f, %.3f, %.3f, %.3f, %.3f, %.3f, %.3f, %.3f"%(
b_flat,
b_eff_u,
b_eff_d,
h_flat,
cog_y,
he1,
he2,
htflat,
thick,
)
)
i_total = cls.rhs_moi(
b_flat,
b_eff_u,
b_eff_d,
h_flat,
cog_y,
he1,
he2,
htflat,
thick,
)
logging.info("Moment of inertia, I: %.3f" % i_total)
# Calculate the elastic section modulus.
w_el = i_total / abs(h_xtr)
logging.info("W_el: %.3f" % w_el)
# Calculate the elastic moment resistance without interaction
#NOTE: this could also be taken from the cs_pb, where it's already
# calculated
m_el_rd = w_el*f_yield/gamma_m(0, country="NO", nominal=False)
logging.info("m_el_rd: %.3f" % m_el_rd)
# Calculate the interaction elastic moment resistance
m_n_rd = (f_yield/gamma_m(0, country="NO", nominal=False) - abs(n_ed)/a_eff)*w_el
logging.info("M_N_el_Rd: %.3f" % m_n_rd)
# Assuming that M_ed = M_N_el_Rd, calculate the moment utilisation
# ratio, normalised against the pure bending moment resistance.
# Note that this calculation stands for the elastic case only
# (m_n_rd is calculated for elastic distribution). If this ends up
# not being true, these values are overwritten by the plastic
# calculations further down the script.
m_prc = m_n_rd/cs_pb.m_rd
# Calculate the extreme compressive and tensile stresses for
# the entire profile due to moment and axial separately
# NOTE: add N*eN to the stress due to moment which occures in case
# of non-symmetric A_eff.
sigma_n = n_ed/a_eff
sigma_c_m = -h_xtr*m_n_rd/i_total
sigma_t_m = (height-h_xtr)*m_n_rd/i_total
logging.info(
(
"Stress due to axial, compr. moment, tens. moment:"
" %.3f, %.3f, %.3f"
) % (sigma_n, sigma_c_m, sigma_t_m)
)
# Calculate the extreme compressive and tensile stresses for
# the web only
sigma_c_m_web = -(h_xtr - r_out)*m_n_rd/i_total
sigma_t_m_web = (height-h_xtr-r_out)*m_n_rd/i_total
logging.info(
(
"Web stress due to compr. moment, tens. moment:"
" %.3f, %.3f"
) % (sigma_c_m_web, sigma_t_m_web)
)
sigma_c_web = n_ed/area + sigma_c_m_web
sigma_t_web = n_ed/area + sigma_t_m_web
# Calculate the distribution factor, psi, for the web
psi_web = sigma_t_web/sigma_c_web
logging.info("psi_web: %.3f" % psi_web)
# Calculate the new tension and compression zone heights of
# the web
# NOTE: Here, I have used abs() for the case that we have a psi
# below 0 (the two ends have opposite sign stress). This means that
# there is no distinction for the orientation of the web, whether
# the tension/compression is up or down. Should this change?
if psi_web < 0:
hcflat = abs(h_flat/(1 - psi_web))
else:
# NOTE: Similarly, when psi is positive, the entire web is
# considered to be under compression. Should I consider the
# case of the entire web being in tension as well. (This case
# should be practically impossible: If the web is entirely
# under tension then both flanges should also be in tension ->
# no classification applies. But in the actual calculations I
# take into account the width of the corner radius. Might there
# be a case where the neutral line is on the corner, leaving
# the entire web is in tension and still have a flange in
# compression???)
hcflat = h_flat
htflat = h_flat - hcflat
logging.info("hcflat, htflat: %3.f, %3.f" % (hcflat, htflat))
# Create an object of the RHS member with the current values
joint = cls(
width=width,
height=height,
thick=thick,
b_flat=b_flat,
h_flat=h_flat,
r_out=r_out,
area=area,
material=material,
n_prc=n_prc_current,
)
(
flange_class,
b_eff_u,
b_ineff_u,
b_eff_d,
b_ineff_d,
web_class,
h_eff,
h_ineff,
alfa_web,
psi_web,
psi_cs,
cs_class
) = joint.cs_classification(force_gross=force_gross)
if cs_class == 4:
logging.info(
"Since the cross section is class 4, there is an issue "
"with the calculation of A_eff and I_eff. Essentially, if "
"the psi of the entire profile is positive, there is an "
"issue with the effective area of the min compression "
"flange. Until this is addressed, the results are not "
"accurate."
)
# Proceed with elastic or plastic design accordingly
logging.info(":::Resistance:::")
if cs_class > 2:
logging.info("---Next iteration OR Start iterations")
logging.info(". . . . . . . . . . . . . . . . . . .")
# Elastic design case
# Calculate the effective area of the cross section
a_eff = area - thick*(2*h_ineff + b_ineff_u + b_ineff_d)
logging.info("A_eff: %.3f" % a_eff)
# Calculate the compressive resistance of the profile with the
# current a_eff.
n_rd = a_eff*f_yield/gamma_m(0, country="NO", nominal=False)
logging.info("n_rd_current: %.3f" % n_rd)
# Calculate the axial utilisation ratio for the current state
# of the effective area.
n_prc_current = n_ed/n_rd
logging.info("n_prc_current: %.3f" % n_prc_current)
# New distance from CoG to web's ineffective area
if psi_web > 0:
he1 = h_eff * 2./(5-psi_web)
he2 = h_eff - he1
else:
he1 = 0.4*h_eff
he2 = 0.6*h_eff
logging.info("he1, he2: %.3f, %.3f" % (he1, he2))
# Calculate the new cog_y
cog_y_new = cls.calc_cog_y(
area,
thick,
htflat,
he2,
h_ineff,
h_flat,
h2um,
b_ineff_u,
b_ineff_d
)
# Calculate the convergence relative to the previous step
approx = abs(cog_y_new - cog_y)
logging.info("Convergence residual: %.3f" % approx)
# Keep the new value of cog_y
cog_y = cog_y_new
else:
# Plastic design case
# No effective area
a_eff = area
# The centre of gravity of the profile is at the centre.
cog_y = height/2.
a_w = (area - 2*width*thick)/area
if a_w > 0.5:
a_w = 0.5
logging.info("a_w: %.3f" % a_w)
m_prc = (1 - abs(n_prc))/(1 - 0.5*a_w)
if m_prc > 1.:
m_prc = 1.
logging.info("m_prc: %.3f" % m_prc)
# Applied moment
m_n_rd = m_pl_rd*m_prc
logging.info("M_N_Rd: %.3f\n" % m_n_rd)
# Populate the joint object with all the remaining calculated
# properties before returning.
joint.a_eff = a_eff
joint.alfa_web = alfa_web
joint.cs_class = cs_class
joint.hcflat = hcflat
joint.m_el_rd = m_el_rd
joint.m_n_rd = m_n_rd
joint.psi_cs = psi_cs
joint.psi_web = psi_web
joint.r_in = r_in
joint.w_el = w_el
joint.w_pl = w_pl
joint.m_prc = m_prc
joint.moi_1 = i_total
joint.moi_xx = i_total
joint.m_pl_rd = m_pl_rd
joint.n_rd = n_rd
# Since the profile is below class 4 and there is no a_off
# approximation, the joint is return immediately and the loop
# is terminated.
return(joint)
# If the profile is class 4, the iterations must continue. Here is
# where the loop ends and returns to the beginning but since there
# are no actions performed when the loop has approximated, it is
# just an indentation step removed.
# After the loop has reached the approximation tolerance, the joint is
# populated with the final values and is returned.
joint.a_eff = a_eff
joint.alfa_web = alfa_web
joint.cs_class = cs_class
joint.hcflat = hcflat
joint.m_el_rd = m_el_rd
joint.m_n_rd = m_n_rd
joint.psi_cs = psi_cs
joint.psi_web = psi_web
joint.r_in = r_in
joint.w_el = w_el
joint.w_pl = w_pl
joint.m_prc = m_prc
joint.moi_1 = i_total
joint.moi_xx = i_total
joint.xc = cog_y
joint.yc = width/2
joint.m_pl_rd = m_pl_rd
joint.n_rd = n_rd
#TODO: add m_prc. This calculation is performed for the case of plastic
# but not for the elastic. It has to be calculated further up in the
# calculations and added to the returned joint here.
return(joint)
[docs] @classmethod
def from_m_prc(
cls,
width,
height,
thick,
material=None,
m_prc=None,
tensile_axial=None,
force_gross=False
):
"""
Alternative constructor for CsRHS.
This method should be used for creating new RHS profile objects instead
of defining them directly through __init__. The overall geometric
properties are given (width, height, thickness), the material and the
moment utilisation.
The profiled is considered to be fully utilised by M+N, so, the stress
distribution assumed for the classification is a result of the
requested moment and an additional axial load, N, large enough to fully
utilise the cross section.
The additional axial load assumed is compressive or tensile, according
to the given input.
Parameters
----------
width : float
height : float
thick : float
material " :obj:`Material`, optional
By default, a generic steel S355 is used, with f_u=510 MPa. The
default is currently disabled so if no metarial is given it
crashes.
m_prc : float, optional
Moment utilisation to the profile. Since the profile is assumed
fully utilised by the combination of M+N, a given |m_prc|<1 implies
an additional axial load. Positive m_prc induce compressive stress
in the upper flange. Default value is 1, which implies a profile
under pure positive bending.
tensile_axial : bool, optional
Defines the sign of the axial when abs(m_prc) < 1. Default is False
which implies compression.
force_gross : bool, optional
Treat class 4 cross sections as a class 3. The argument is passed
to the "cs_classification" method. Default is false.
"""
# Default values
if m_prc is None:
m_prc = 1
if tensile_axial is None:
tensile_axial = False
# To achieve the requested values, iterations of the from_n_prc
# classmethod are used.
# For the initial n_prc estimation, elastic distribution for the gross
# cross-section is assumed.
if tensile_axial:
n_prc = 1 - abs(m_prc)
else:
n_prc = abs(m_prc) - 1
# Create a profile for the initial n_prc estimation
cs = cls.from_n_prc(
width,
height,
thick,
material = material,
n_prc = n_prc,
negative_moment = m_prc<0,
force_gross = force_gross
)
return(cs)
[docs] def cs_classification(self, force_gross=False):
"""
Method performing classification of an RHS profile.
Parameters
----------
force_gross : bool, optional
If True, effective width reduction for class 4 cross-sections will
be disregarded. The returned cs_class will still be 4 but the
effective widths will be the full widths (equivalent to a class 3).
Default is False.
"""
logging.info("Begin classification method")
# Calculate the non-dimensional width of the flanges webs (c/et).
self.nd_width_flange = self.b_flat/(self.thick*self.material.epsilon)
self.nd_width_web = self.h_flat/(self.thick*self.material.epsilon)
logging.info(
"Non dimensional widths [flange, web]: %.3f, %.3f"%(
self.nd_width_flange,
self.nd_width_web
)
)
# Calculate the psi and psi_web (assuming elastic distribution) and
# alfa (assuming plastic distribution) for the given n_prc.
# Psi and psi_web. For psi_web, the ratio of the flange and corner
# height to the total height of the profile are calculated.
h_c_ratio = self.r_out/self.height
if self.n_prc <= 0:
psi = -(2*self.n_prc + 1)
psi_web_num = 2*(self.n_prc - h_c_ratio*(self.n_prc + 1)) + 1
psi_web_denom = 2*h_c_ratio*(self.n_prc + 1) - 1
psi_web = psi_web_num/psi_web_denom
else:
psi = 1/(2*self.n_prc -1)
psi_web_num = 2*h_c_ratio*(self.n_prc - 1) + 1
psi_web_denom = 2*(self.n_prc + h_c_ratio*(self.n_prc - 1)) - 1
psi_web = psi_web_num/psi_web_denom
logging.info("Psi_web: %.3f" % psi_web)
# a_w is calculated as a function of n_prc. First, the ratio of the
# web area to the total cs area is needed, w_prc. w_prc is needed to
# evaluate if the plastic neutral axis is in the web or in one of the
# flanges.
w_prc = self.h_flat/self.area
logging.info("debugz. n_prc: %.3f, w_prc: %.3f" %(self.n_prc, w_prc))
if -1 <= self.n_prc < -w_prc:
alfa_web = 1
elif -w_prc <= self.n_prc <= w_prc:
alfa_web = (w_prc - self.n_prc)/(2*w_prc)
else:
alfa_web = 0
logging.info("alfa_web: %.3f" % alfa_web)
# Web classification
# Avoid classification if the web is entirely in tension.
logging.info(":::Web Classification:::")
if psi_web > 1:
logging.info("The web is entirely in tension, no classification")
web_class = NotImplemented
else:
# Perform the classification
web_class = self.classification(
self.nd_width_web,
alfa_web,
psi_web,
)
logging.info("Web class: %.3f" % web_class)
# Calculate the height of the compression zone of the web.
if psi_web >= 0:
hcflat = self.h_flat
else:
hcflat = self.h_flat/(1 - psi_web)
# Reduce to effective height if needed.
h_eff = hcflat
if web_class == 4 and not(force_gross):
h_eff = self.calc_eff_width(
self.h_flat,
self.thick,
psi_web,
self.material.epsilon,
internal=True
)
# Calculate the ineffective height for a web fully in compression
# or only partially. Here, the distinction is made based on the
# psi_web which could be also positive for a web entirely in
# tension but this case was already escaped just before the
# classification.
# NOTE: he following conditional can be avoided and replaced by
# h_ineff = hcflat - h_eff, because this check is already performed
# a few lines above. Cross-check and implement.
if psi_web >= 0:
h_ineff = self.h_flat - h_eff
else:
h_ineff = hcflat - h_eff
logging.info("h_eff, h_ineff: %.3f, %.3f" % (h_eff, h_ineff))
# Flange classification
# Avoid classification if both flanges are in tension.
logging.info(":::Flange Classification:::")
if psi > 1:
logging.info("Both flanges are in tension, no classification.")
b_eff = NotImplemented
else:
# Perform classification
flange_class = self.classification(
self.nd_width_flange,
1,
1,
)
logging.info("Flange class: %.3f" % flange_class)
b_eff = self.b_flat
if flange_class == 4 and not(force_gross):
b_eff = self.calc_eff_width(
self.b_flat,
self.thick,
1,
self.material.epsilon,
internal=True
)
#temp_del
#b_eff = self.b_flat
b_ineff = self.b_flat - b_eff
logging.info("b_eff, b_ineff: %.3f, %.3f" % (b_eff, b_ineff))
# Decide the profile's class based on the maximum class within all
# parts.
cs_class = max(flange_class, web_class)
logging.info("Cross-section class: %.3f" % cs_class)
# Figure out which flanges are in compression.
# in case of class 4, assign the effective widths to the compression
# flanges.
if cs_class == 4:
# In an elastic distribution, if the compressive axial utilisation
# is more than half the cs's capacity, then both flanges are under
# compression.
if self.n_prc < -0.5:
b_eff_u, b_ineff_u = b_eff, b_ineff
b_eff_d, b_ineff_d = b_eff, b_ineff
logging.info(
"Both flanges in compression and reduced to effective."
)
# Equivalently, if the tensile axial utilisation is greater than
# half of the cs's capacity, both flanges are in tension.
elif self.n_prc > 0.5:
b_eff_u, b_ineff_u = self.b_flat, 0
b_eff_d, b_ineff_d = self.b_flat, 0
logging.info(
"Both flanges in tension, no effective width reduction."
)
# for all cases in-between, only one flange under compression.
# Whether it is the upper or the lower flange that is under
# compression is decided by the direction of the moment.
# (reminder: positive moment -> upper flange in compression)
else:
if not self.negative_moment:
b_eff_u, b_ineff_u = b_eff, b_ineff
b_eff_d, b_ineff_d = self.b_flat, 0
logging.info(
"Upper flange in compression, reduced to effective."
)
else:
b_eff_u, b_ineff_u = self.b_flat, 0
b_eff_d, b_ineff_d = b_eff, b_ineff
logging.info(
"Lower flange in compression, reduced to effective."
)
else:
b_eff_u, b_ineff_u = self.b_flat, 0
b_eff_d, b_ineff_d = self.b_flat, 0
# NOTE:I could also return alfa_flange, alfa_web, psi_flange, psi_web
# so that I can check them in the returned CsRHS object. I did not
# really need this so I left it like that.
return(
flange_class,
b_eff_u,
b_ineff_u,
b_eff_d,
b_ineff_d,
web_class,
h_eff,
h_ineff,
alfa_web,
psi_web,
psi,
cs_class
)
# TODO: merge the following 5 methods with the general classification
# functions
[docs] @staticmethod
def classification(nd_width, alfa, psi):
#def classification(width, thick, material, alfa, psi):
if CsRHS.plastic_allowed(nd_width, alfa):
logging.info(
"Is plastic distribution allowed? - %s" %CsRHS.plastic_allowed(
nd_width,
alfa,
)
)
if CsRHS.class1(nd_width, alfa):
pl_class = 1
else:
pl_class = 2
else:
if CsRHS.class4(nd_width, psi):
pl_class = 4
else:
pl_class = 3
return(pl_class)
[docs] @staticmethod
def plastic_allowed(nd_width, alfa):
"""
Static method for checking if a plate exceeds the upper limit of class
2.
The limits are described in table 5.2 of EN1993-1-1. In that table, the
cases of pure compression and pure bending are treated separately. The
general functions that (rightmost column) are continuous on the pure
bending/pure compression limits so it is sufficient to use only them.
Parameters
----------
nd_width : float
Non-dimensional width of the plate, c/(et).
alfa : float in [0, 1]
Ratio of the compression zone height to the overall height of the
profile.
Returns
-------
bool
"""
if alfa > -0.5:
criterion = 456./(13*alfa - 1)
else:
criterion = 41.5/alfa
logging.info("Crit: %.3f, nd: %.3f" %(criterion, nd_width) )
if nd_width <= criterion:
status = True
else:
status = False
logging.info("Is plastic allowed?: %s" %status)
return status
[docs] @staticmethod
def class1(nd_width, alfa):
"""
Static method for checking if a plate exceeds the upper limit of class
1.
The limits are described in table 5.2 of EN1993-1-1. In that table, the
cases of pure compression and pure bending are treated separately. The
general functions that (rightmost column) are continuous on the pure
bending/pure compression limits so it is sufficient to use only them.
Parameters
----------
nd_width : float
Non-dimensional width of the plate, c/(et).
alfa : float in [0, 1]
Ratio of the compression zone height to the overall height of the
profile.
Returns
-------
bool
"""
status = False
if alfa > -0.5:
criterion = 396./(13*alfa - 1)
else:
criterion = 36./alfa
if nd_width <= criterion:
status = True
return status
[docs] @staticmethod
def class4(nd_width, psi):
"""
Static method for checking if a plate exceeds the lower limit of class
4.
The limits are described in table 5.2 of EN1993-1-1. In that table, the
cases of pure compression and pure bending are treated separately. The
general functions that (rightmost column) are continuous on the pure
bending/pure compression limits so it is sufficient to use only them.
Parameters
----------
nd_width : float
Non-dimensional width of the plate, c/(et).
alfa : float in [0, 1]
Ratio of the compression zone height to the overall height of the
profile.
Returns
-------
bool
"""
status = False
if psi > -1:
criterion = 42 / (0.67 + 0.33*psi)
else:
criterion = 62*(1-psi)*(-psi)**0.5
if nd_width >= criterion:
status = True
return status
[docs] @staticmethod
def calc_eff_width(bbbb, thick, psi, epsilon, internal=True):
# Distinguish between internal and outstanding elements
if internal:
# Calculate k_sigma based on psi
if (psi <= 1) and (psi >= 0):
k_sigma = 8.2 / (1.05 + psi)
elif (psi < 0) and (psi >= -1):
k_sigma = 7.81 - 6.29*psi + 9.78*psi**2
elif (psi < -1) and (psi > -3):
k_sigma = 5.98*(1 - psi)**2
else:
logging.error("Wrong psi value. psi: %.3f" % psi)
return
lamda_p = (bbbb/thick) / (28.4 * epsilon * np.sqrt(k_sigma))
logging.info("lamda_p: %.3f, ", lamda_p)
if lamda_p > 0.673:
logging.info("Reduced to effective width.")
rho = (lamda_p-0.055*(3+psi)) / lamda_p**2
if rho > 1:
rho = 1.
else:
logging.info("No effective width reduction")
rho = 1.
logging.info("rho: %.3f, " % rho)
else:
# Calculate k_sigma based on psi
if (psi <= 1) and (psi >= 0):
k_sigma = 0.578 / (psi + 0.34)
elif (psi < 0) and (psi >= -1):
k_sigma = 1.7 - 5*psi + 17.1*psi**2
else:
logging.error("Wrong psi value. psi: %.3f" % psi)
pass
lamda_p = (bbbb/thick) / (28.4 * epsilon * np.sqrt(k_sigma))
logging.info("lamda_p: %.3f" % lamda_p)
if lamda_p > 0.748:
logging.info("Reduced to effectice width")
rho = (lamda_p - 0.188)/lamda_p**2
if rho > 1:
rho = 1.
else:
logging.info("No effective width reduction")
rho = 1.0
# Compression zone width
if psi < 0:
b_c = bbbb / (1 - psi)
else:
b_c = bbbb
# Effective width
b_eff = rho * b_c
return(b_eff)
[docs] @staticmethod
def rhs_moi(
bflat,
b_eff_u,
b_eff_d,
hflat,
cog_y,
he1,
he2,
h_tension,
thick,
):
# NOTE: Complete the help doc string.
"""
Calculate the moment of inertia for a given RHS.
If the RHS is class 4 and effective cross-section reduction is needed,
the effective regions have to be given.
Parameters
----------
Returns
-------
"""
# Outer radius of the profile corners
r_out, r_in = CsRHS.radius_from_thickness(thick)
# NOTE: The origin of the assumed csys for the cross-section is at the
# centre of symmetry.
# Heights of the flat portion of the web, below and above the cog.
h_below = hflat/2. + cog_y
h_above = hflat/2. - cog_y
# Calculate distances from cog to the mid-lines of the two flanges.
h2bm = h_below + r_out - thick/2
h2am = h_above + r_out - thick/2
logging.debug("h2bm, h2am: %.3f, %.3f" % (h2bm, h2am))
# Height of the effective zone on the lower part of the web (the one
# with lower compression or even tension) and dist to cog
he_l = he2 + h_tension
he_l_2cog = h_below - he_l/2
he_u_2cog = h_above - he1/2
logging.debug("he_l, he_l_2cog, he_u_2cog: %.3f, %.3f, %.3f" % (
he_l,
he_l_2cog,
he_u_2cog
))
# Calculate partial moments of inertia
i_u_flange = (b_eff_u * thick**3)/12. + b_eff_u*thick*h2am**2
logging.debug("I - i_u_flange: %.3f" % i_u_flange)
i_d_flange = (b_eff_d * thick**3)/12. + b_eff_d*thick*h2bm**2
logging.debug("I - i_d_flange: %.3f" % i_d_flange)
i_l_web = (thick*he_l**3)/12. + thick*he_l*he_l_2cog**2
logging.debug("I - i_l_web: %.3f" % i_l_web)
i_u_web = (thick*he1**3)/12. + thick*he1*he_u_2cog**2
logging.debug("I - i_u_web: %.3f" % i_u_web)
i_semitorus = (np.pi/8)*(r_out**4-r_in**4)
logging.debug("I - i_semitorus: %.3f" % i_semitorus)
i_t_c = i_semitorus + (np.pi/2)*(r_out**2 - r_in**2)*h_below**2
logging.debug("I - i_t_c: %.3f" % i_t_c)
i_c_c = i_semitorus + (np.pi/2)*(r_out**2 - r_in**2)*h_above**2
logging.debug("I - i_c_c: %.3f" % i_c_c)
# Total moment of inertia
i_total = i_d_flange +\
i_u_flange +\
2*(i_l_web + i_u_web) +\
i_t_c +\
i_c_c
return(i_total)
[docs] @staticmethod
def radius_from_thickness(thickness):
# Corner radius calculated acc. to SSAB's rhs profile specs (taken from
# document: "STRUCTURAL HOLLOW SECTIONS - dimensions and
# cross-sectional properties-2016-Confetti"
if thickness <= 6.:
r_out = 2 * thickness
elif (thickness > 6.) and (thickness <= 10.):
r_out = 2.5 * thickness
else:
r_out = 3 * thickness
r_in = r_out - thickness
return(r_out, r_in)
[docs] @staticmethod
def calc_flat_width_radius(hhhh, bbbb, thick):
r_out, r_in = CsRHS.radius_from_thickness(thick)
h_flat = hhhh - 2*r_out
b_flat = bbbb - 2*r_out
return(h_flat, b_flat, r_out, r_in)
[docs] @staticmethod
def calc_area(hhhh, bbbb, thick):
h_flat, b_flat, r_out, r_in = CsRHS.calc_flat_width_radius(
hhhh,
bbbb,
thick
)
a_web = 2*thick*h_flat
a_flange = 2*thick*b_flat
a_corners = np.pi * (r_out**2 - r_in**2)
area = a_web + a_flange + a_corners
return(area)
[docs] @staticmethod
def calc_w_pl(
hflat,
bflat,
thick,
h2um,
r_in,
r_out
):
# Calculate the partial section modulii
w_f_p = bflat*thick*h2um
w_w_p = (hflat/2.)*thick*(hflat/4.)
area_c = (np.pi/2)*(r_out**2 - r_in**2)
# Formula for the centre of gravity of a semi-torus taken from
# Roark, pg. 809 (there is an typo on the given formula: the
# r_i^2 in the numerator should be r_i^3. Here, in the script,
# it is correct)
cog_c = 4.*(r_out**3 - r_in**3)/(3*np.pi*(r_out**2 - r_in**2))
w_c_p = area_c*(hflat/2. + cog_c)
# Plastic section modulus
w_pl = 2*(w_f_p + 2*w_w_p + w_c_p)
return(w_pl)
[docs] @staticmethod
def calc_cog_y(
area,
thick,
htflat,
he2,
h_ineff,
h_flat,
h2um,
b_ineff_u,
b_ineff_d
):
# Calculate the gravity centre for the effective cross-section
# (if any). The y-axis is positive pointing on the compression
# side of the acting moment.
h_inef_cog = htflat + he2 + h_ineff/2. - h_flat/2
cog_num = thick*(
-2*h_ineff*h_inef_cog - h2um*(b_ineff_u + b_ineff_d)
)
cog_den = area - thick*(2*h_ineff + b_ineff_u + b_ineff_d)
cog_y = cog_num / cog_den
logging.info(
"h_inef_cog, cog_y: %.3f, %.3f" % (h_inef_cog, cog_y)
)
return(cog_y)
# TODO: check the plastic table values. expand the library
[docs]class Material:
"""
Material properties.
Parameters
----------
e_modulus : float
Modulus of elasticity.
poisson : float
Poisson's ratio.
f_yield : float
Yield stress
plasticity : tuple
Plasticity table (tuple of stress-plastic strain pairs).
By default, no plasticity is considered.
"""
[docs] def __init__(
self,
e_modulus,
poisson,
f_y_nominal,
f_y_real=None,
f_u_nominal=None,
plasticity=None
):
self.e_modulus = e_modulus
self.poisson = poisson
self.f_y_nominal = f_y_nominal
self.f_y_real = f_y_real
self.f_u_nominal = f_u_nominal
self.plasticity = plasticity
self.epsilon = np.sqrt(235. / f_y_nominal)
[docs] @staticmethod
def plastic_table(nominal=None):
"""
Plasticity tables.
Returns a tuple with plastic stress-strain curve values for different
steels given a steel name, e.g 'S355'
Parameters
----------
nominal : string [optional]
Steel name. Default value, 'S355'
Attributes
----------
Notes
-----
References
----------
"""
if nominal == None:
nominal = 'S235'
if nominal == 'S355':
table = (
(355.0, 0.0),
(360.0, 0.015),
(390.0, 0.0228),
(420.0, 0.0315),
(440.0, 0.0393),
(480.0, 0.0614),
(520.0, 0.0926),
(550.0, 0.1328),
(570.0, 0.1746),
(585.0, 0.2216),
(586.0, 1.)
)
if nominal == 'S381':
table = (
(381.1, 0.0),
(391.2, 0.0053),
(404.8, 0.0197),
(418.0, 0.0228),
(444.2, 0.0310),
(499.8, 0.0503),
(539.1, 0.0764),
(562.1, 0.1009),
(584.6, 0.1221),
(594.4, 0.1394),
(596, 1.)
)
if nominal == 'S650':
table = (
(760., 0.0),
(770., 0.022),
(850., 0.075),
(900., 0.1),
(901., 1.)
)
if nominal == 'S700':
table = (
(300., 0.00000),
(400., 4.53e-5),
(450., 8.24e-5),
(500., 1.41e-4),
(550., 2.47e-4),
(600., 4.45e-4),
(630., 6.52e-4),
(660., 1.00e-3),
(700., 2.07e-3),
(720., 3.37e-3),
(750., 9.18e-3),
(770., 1.41e-2),
(790., 2.01e-2),
(800., 2.38e-2),
(820., 3.26e-2),
(840., 4.56e-2),
(850., 5.52e-2),
(860., 7.11e-2),
(865., 9.34e-2),
(870., 1.00)
)
return table
[docs] @classmethod
def from_nominal(cls, nominal_strength=None):
"""
Alternative constructor creating a steel material from a given nominal
strength.
Parameters
----------
nominal_strength : str
Steel quality, given in the form of e.g. "S355"
"""
if nominal_strength is None:
f_yield = 235.
else:
f_yield = float(nominal_strength.replace('S', ''))
plasticity = cls.plastic_table(nominal=nominal_strength)
return cls(210000., 0.3, f_yield, plasticity=plasticity)
[docs]class BCs:
"""Boundary conditions
BC objects are meant to be used within 'part' objects and describe the
restrained degrees of freedom at the two ends.
"""
[docs] def __init__(self, bcs):
self.bcs = bcs
[docs] @classmethod
def from_hinged(cls):
return cls([[1, 1, 1, 0, 0, 0], [1, 1, 1, 0, 0, 0]])
[docs]class CsLoads:
"""Loads of a part
"""
[docs] def __init__(
self,
axial,
shear,
moment
):
self.axial = axial
self.shear = shear
self.moment = moment
[docs]class StructProps:
""" Structural properties of a member."""
[docs] def __init__(
self,
t_classification=None,
p_classification=None,
pc_classification=None,
lmbda_y=None,
lmbda_z=None,
n_cr_plate=None,
sigma_cr_plate=None,
n_pl_rd=None,
n_b_rd_plate=None,
sigma_b_rd_plate=None,
sigma_cr_shell=None,
sigma_cr_shell_new=None,
lenca=None,
lenca_new=None,
n_cr_shell=None,
n_cr_shell_new=None,
sigma_b_rk_shell=None,
sigma_b_rk_shell_new=None,
n_b_rk_shell=None,
n_b_rk_shell_new=None,
sigma_b_rd_shell=None,
sigma_b_rd_shell_new=None,
n_b_rd_shell=None,
n_b_rd_shell_new=None,
n_cr_platec=None,
sigma_cr_platec=None,
n_pl_rdc=None,
sigma_b_rd_platec=None,
n_b_rd_platec=None,
):
"""
Main constructor
Parameters
----------
t_classification : float, optional
Classification of a tube, d/(t^2*e)
p_classification : float, optional
Classification of a plate, c/(t*e)
lmbda_y : float, optional
Flexural slenderness on the strong axis.
lmbda_z : float, optional
Flexural slenderness on the weak axis.
n_pl_rd : float, optional
Plastic axial compression resistance.
n_b_rd_shell : float, optional
Shell buckling resistance
"""
self.t_classification = t_classification
self.p_classification = p_classification
self.pc_classification = pc_classification
self.lmbda_y = lmbda_y
self.lmbda_z = lmbda_z
self.n_cr_plate = n_cr_plate
self.sigma_cr_plate = sigma_cr_plate
self.n_pl_rd = n_pl_rd
self.n_b_rd_plate = n_b_rd_plate
self.sigma_b_rd_plate = sigma_b_rd_plate
self.sigma_cr_shell = sigma_cr_shell
self.sigma_cr_shell_new = sigma_cr_shell_new
self.lenca = lenca
self.lenca_new = lenca_new
self.n_cr_shell = n_cr_shell
self.n_cr_shell_new = n_cr_shell_new
self.sigma_b_rk_shell = sigma_b_rk_shell
self.sigma_b_rk_shell_new = sigma_b_rk_shell_new
self.n_b_rk_shell = n_b_rk_shell
self.n_b_rk_shell_new = n_b_rk_shell_new
self.sigma_b_rd_shell = sigma_b_rd_shell
self.sigma_b_rd_shell_new = sigma_b_rd_shell_new
self.n_b_rd_shell = n_b_rd_shell
self.n_b_rd_shell_new = n_b_rd_shell_new
self.n_cr_platec = n_cr_platec
self.sigma_cr_platec = sigma_cr_platec
self.n_pl_rdc = n_pl_rdc
self.sigma_b_rd_platec = sigma_b_rd_platec
self.n_b_rd_platec = n_b_rd_platec
[docs]class Part:
"""
Structural part.
Class describing a structural part, including geometry, boundary conditions
loads and resistance.
"""
[docs] def __init__(self,
geometry=None,
cs_props=None,
material=None,
struct_props=None,
bc_loads=None
):
"""
Parameters
----------
geometry : Geometry object, optional
cs_props : CsProps object, optional
material : Material object, optional
struct_props : StructProps object, optional
bc_loads: BCs object, optional
"""
self.geometry = geometry
self.cs_props = cs_props
self.material = material
self.bc_loads = bc_loads
self.struct_props = struct_props
[docs]class SimplySupportedPlate:
"""
"""
[docs] def __init__(self, width, thickness, length, f_y, psi=None):
if psi is None:
psi = 1
self.width = width
self.thickness = thickness
self.length = length
self.f_y = f_y
self.psi = psi
self.area = width * thickness
self.plate_class = plate_class(thickness, width, f_y)
self.sigma_cr = sigma_cr_plate(thickness, width, psi=psi)
self.a_eff = calc_a_eff(thickness, width, f_y, psi)
# SIMPLY SUPPORTED PLATE
#TODO: Implement EN50341. Currently the resistance is calculated only for pure compression elements. Add interaction.
[docs]def plate_class(
thickness,
width,
f_yield
):
# Docstring
"""
Plate classification.
Returns the class for a given plate, according to EN1993-1-1.
Currently works for simply supported plates under pure compression.
Parameters
----------
thickness : float
[mm] Plate thickness
width : float
[mm] Plate width
f_yield : float
[MPa] Yield stress
Returns
-------
int
[_] Class number
Notes
-----
To be extended to include the rest of the cases of Table 5.3 [1].
Members under combined axial and bending and outstand members.
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-1: General rules and rules for buildings. Brussels: CEN, 2005
"""
# Convert inputs to floats
width, thickness, f_yield = float(width), float(thickness), float(f_yield)
# Calculate classification
classification = width / (thickness * np.sqrt(235. / f_yield))
if classification <= 33.01:
p_class = 1
elif classification <= 38.01:
p_class = 2
elif classification <= 42.01:
p_class = 3
else:
p_class = 4
# Return value
return p_class
[docs]def plate_class_new(
thickness,
width,
f_yield
):
# Docstring
"""
Plate classification acc.to final draft of EC3-1-1.
Returns the class for a given plate, according to EN1993-1-1.
Currently works for simply supported plates under pure compression.
Parameters
----------
thickness : float
[mm] Plate thickness
width : float
[mm] Plate width
f_yield : float
[MPa] Yield stress
Returns
-------
int
[_] Class number
Notes
-----
To be extended to include the rest of the cases of Table 5.3 [1].
Members under combined axial and bending and outstand members.
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-1: General rules and rules for buildings. Brussels: CEN, 2005
"""
# Convert inputs to floats
width, thickness, f_yield = float(width), float(thickness), float(f_yield)
# Calculate classification
classification = width / (thickness * np.sqrt(235. / f_yield))
if classification <= 28.01:
p_class = 1
elif classification <= 34.01:
p_class = 2
elif classification <= 38.01:
p_class = 3
else:
p_class = 4
# Return value
return p_class
[docs]def sigma_cr_plate(
thickness,
width,
psi=None,
young=210000,
poisson=0.3
):
# Docstring
"""
Critical stress of a plate.
Calculates the critical stress for a simply supported plate.
Parameters
----------
thickness : float
[mm] Plate thickness
width : float
[mm] Plate width
psi : float, optional
[_] Ratio of the min over max stress for a linear distribution,
(sigma_min / sigma_max)
Default = 1, which implies a uniform distribution
young : float, optional
[MPa] Young's modulus of elasticity. Default value is 210000
poisson : float, optional
[_] Poisson ratio, default value is 0.3
Returns
-------
float
[MPa] Plate critical stress
Notes
-----
To be extended to include cantilever plate (outstand members)
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-5: Plated structural elements. Brussels: CEN, 2005.
"""
# Convert inputs to floats
thickness, width = float(thickness), float(width)
# Default value for psi
if psi is None:
psi = 1.
else:
psi = float(psi)
# Calculate kapa_sigma
if psi>0:
k_sigma = 8.2 / (1.05 + psi)
elif 0 >= psi > -1:
k_sigma = 7.81 - 6.29*psi + 9.78*psi**2
else:
k_sigma = 5.98*(1 - psi)**2
# Elastic critical stress acc. to EN3-1-5 Annex A
sigma_e = (np.pi**2 * young * thickness) / (12*(1 - poisson**2) * width**2)
sigma_cr = sigma_e * k_sigma
# Return value
return sigma_cr
[docs]def calc_a_eff(
thickness,
width,
f_yield,
psi=None):
# Docstring
"""
Plastic design resistance of a plate.
Calculates the resistance of a plate according to EN1993-1-1 and
EN1993-1-5. The plate is assumed simply supported.
Parameters
----------
thickness : float
[mm] Plate thickness
width : float
[mm] Plate width
f_yield : float
[MPa] Yield stress
psi : float, optional
[_] Ratio of the min over max stress for a linear distribution,
(sigma_min / sigma_max)
Default = 1, which implies a uniform distribution
Returns
-------
float
[N] Plastic design resistance
Notes
-----
To be extended to include cantilever plate (outstand members)
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-1: General rules and rules for buildings.
Brussels: CEN, 2005.
.. [2] Eurocode 3: Design of steel structures - Part 1-5: Plated structural elements. Brussels: CEN, 2005.
"""
# Convert inputs to floats
thickness, width, f_yield = float(thickness), float(width), float(f_yield)
# Default value for psi
if psi is None:
psi = 1.
else:
psi = float(psi)
# Calculate kapa_sigma
k_sigma = 8.2 / (1.05 + psi)
# Aeff calculation.
# Reduction factor for the effective area of the profile acc. to EC3-1-5
classification = width / (thickness * np.sqrt(235. / f_yield))
lambda_p = classification / (28.4 * np.sqrt(k_sigma))
if lambda_p > 0.673 and plate_class(thickness, width, f_yield) == 4:
rho = (lambda_p - 0.055 * (3 + psi)) / lambda_p ** 2
else:
rho = 1.
# Effective area
a_eff = rho * thickness * width
return(a_eff)
[docs]def calc_a_eff_new(
thickness,
width,
f_yield,
psi=None):
# Docstring
"""
Plastic design resistance of a plate.
Calculates the resistance of a plate according to EN1993-1-1 and
EN1993-1-5. The plate is assumed simply supported.
Parameters
----------
thickness : float
[mm] Plate thickness
width : float
[mm] Plate width
f_yield : float
[MPa] Yield stress
psi : float, optional
[_] Ratio of the min over max stress for a linear distribution,
(sigma_min / sigma_max)
Default = 1, which implies a uniform distribution
Returns
-------
float
[N] Plastic design resistance
Notes
-----
To be extended to include cantilever plate (outstand members)
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-1: General rules and rules for buildings.
Brussels: CEN, 2005.
.. [2] Eurocode 3: Design of steel structures - Part 1-5: Plated structural elements. Brussels: CEN, 2005.
"""
# Convert inputs to floats
thickness, width, f_yield = float(thickness), float(width), float(f_yield)
# Default value for psi
if psi is None:
psi = 1.
else:
psi = float(psi)
# Calculate kapa_sigma
k_sigma = 8.2 / (1.05 + psi)
# Aeff calculation.
# Reduction factor for the effective area of the profile acc. to EC3-1-5
classification = width / (thickness * np.sqrt(235. / f_yield))
lambda_p = classification / (28.4 * np.sqrt(k_sigma))
if lambda_p > 0.673 and plate_class_new(thickness, width, f_yield) == 4:
rho = (lambda_p - 0.055 * (3 + psi)) / lambda_p ** 2
else:
rho = 1.
# Effective area
a_eff = rho * thickness * width
return(a_eff)
# CYLINDRICAL SHELLS
[docs]def tube_class(
thickness,
radius,
f_yield
):
"""
CHS classification.
Returns the class for a given plate, according to EN1993-1-1.
Currently works for simply supported plates under pure compression.
Parameters
----------
thickness : float
[mm] Plate thickness
radius : float
[mm] Tube radius
f_yield : float
[MPa] Yield stress
Returns
-------
int
[_] Class number
Notes
-----
To be extended to include the rest of the cases of Table 5.3 [1].
Members under combined axial and bending and outstand members.
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-1: General rules and rules for buildings. Brussels: CEN, 2005
"""
# Convert inputs to floats
radius, thickness, f_yield = float(radius), float(thickness), float(f_yield)
# Calculate classification
classification = 2 * radius / (thickness * np.sqrt(235. / f_yield)**2)
if classification <= 50.01:
t_class = 1
elif classification <= 70.01:
t_class = 2
elif classification <= 90.01:
t_class = 3
else:
t_class = 4
# Return value
return t_class
[docs]def shell_length_category(
radius,
thickness,
length
):
"""Return the length gategory of a cylinder acc. to EC3-1-6 D.1.2.1"""
omega = length / np.sqrt(radius * thickness)
if 1.7 <= omega <= 0.5 * (radius / thickness):
length_category = 1
elif omega < 1.7:
length_category = 0
else:
length_category = 2
return length_category
[docs]def shell_length_category_new(
radius,
thickness,
length
):
"""Return the length gategory of a cylinder acc. to EC3-1-6 D.1.2.1"""
omega = length / np.sqrt(radius * thickness)
if 1.7 <= omega <= 2.86 * (radius / thickness):
length_category = 1
elif omega < 1.7:
length_category = 0
else:
length_category = 2
return length_category
[docs]def sigma_x_rcr(
thickness,
radius,
length,
kapa_bc=None,
e_modulus=None
):
"""
Critical meridional stress for cylindrical shell.
Calculates the critical load for a cylindrical shell under pure
compression and assumes uniform stress distribution. Calculation
according to EN1993-1-6 [1], Annex D.
Parameters
----------
thickness : float
[mm] Shell thickness
radius : float
[mm] Cylinder radius
length : float
[mm] Cylnder length
Returns
-------
list
List of 2 elements:
a) float, Critical load [N]
b) string, length category
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-6: Strength and stability of shell structures.
Brussels: CEN, 2006.
"""
if kapa_bc is None:
kapa_bc = 1
if e_modulus is None:
e_modulus = 210000.
# Convert inputs to floats
thickness, radius, length = float(thickness), float(radius), float(length)
# Length category
lenca = shell_length_category(radius, thickness, length)
# Elastic critical load acc. to EN3-1-6 Annex D
omega = length / np.sqrt(radius * thickness)
if lenca == 1:
c_x = 1.
elif lenca == 0:
c_x = 1.36 - (1.83 / omega) + (2.07 / omega ** 2)
elif lenca == 2:
# c_x_b is read on table D.1 of EN3-1-5 Annex D acc. to BCs
# BC1 - BC1 is used on the Abaqus models (both ends clamped, see EN3-1-5 table 5.1)
c_x_b = 6.
c_x_n = max((1 + 0.2 * (1 - 2 * omega * thickness / radius) / c_x_b), 0.6)
c_x = c_x_n
# Calculate critical stress, eq. D.2 on EN3-1-5 D.1.2.1-5
sigma_cr = 0.605 * 210000 * c_x * thickness / radius
return sigma_cr
[docs]def sigma_x_rcr_new(
thickness,
radius,
length,
kapa_bc=None,
e_modulus=None
):
"""
Critical meridional stress for cylindrical shell.
Calculates the critical load for a cylindrical shell under pure
compression and assumes uniform stress distribution. Calculation
according to EN1993-1-6 [1], Annex D.
Parameters
----------
thickness : float
[mm] Shell thickness
radius : float
[mm] Cylinder radius
length : float
[mm] Cylnder length
Returns
-------
list
List of 2 elements:
a) float, Critical load [N]
b) string, length category
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-6: Strength and stability of shell structures.
Brussels: CEN, 2006.
"""
if kapa_bc is None:
kapa_bc = 1
if e_modulus is None:
e_modulus = 210000.
# Convert inputs to floats
thickness, radius, length = float(thickness), float(radius), float(length)
lenca = shell_length_category_new(radius, thickness, length)
omega = length / np.sqrt(radius * thickness)
if lenca == 0:
c_x = 1.36 - (1.83 / omega) + (2.07 / omega ** 2)
# Calculate critical stress, eq. D.2 on EN3-1-5 D.1.2.1-5
sigma_cr = 0.605 * 210000 * c_x * thickness / radius
elif lenca == 1:
c_x = 1.
# Calculate critical stress, eq. D.2 on EN3-1-5 D.1.2.1-5
sigma_cr = 0.605 * 210000 * c_x * thickness / radius
elif lenca == 2:
# flex buckling
r_o = radius + thickness / 2.
r_i = radius - thickness / 2.
moi = np.pi * (r_o ** 4 - r_i ** 4) / 4.
area = 2 * np.pi * radius * thickness
sigma_cr = n_cr_flex(length, moi, kapa_bc=kapa_bc, e_modulus=e_modulus) / area
else:
print("Wrong length category.")
# Return value
return sigma_cr
[docs]def n_cr_shell(
thickness,
radius,
length
):
"""
Critical compressive load for cylindrical shell.
Calculates the critical load for a cylindrical shell under pure
compression and assumes uniform stress distribution. Calculation
according to EN1993-1-6 [1], Annex D.
Parameters
----------
thickness : float
[mm] Shell thickness
radius : float
[mm] Cylinder radius
length : float
[mm] Cylnder length
Returns
-------
float
[N] Critical load
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-6: Strength and stability of shell structures.
Brussels: CEN, 2006.
"""
# Convert inputs to floats
thickness, radius, length = float(thickness), float(radius), float(length)
# Elastic critical load acc to EN3-1-6 Annex D
nn_cr_shell = 2 * np.pi * radius * thickness * sigma_x_rcr(thickness, radius, length)
# Return value
return nn_cr_shell
[docs]def n_cr_shell_new(
thickness,
radius,
length
):
"""
Critical compressive load for cylindrical shell.
Calculates the critical load for a cylindrical shell under pure
compression and assumes uniform stress distribution. Calculation
according to EN1993-1-6 [1], Annex D.
Parameters
----------
thickness : float
[mm] Shell thickness
radius : float
[mm] Cylinder radius
length : float
[mm] Cylnder length
Returns
-------
float
[N] Critical load
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-6: Strength and stability of shell structures.
Brussels: CEN, 2006.
"""
# Convert inputs to floats
thickness, radius, length = float(thickness), float(radius), float(length)
# Elastic critical load acc to EN3-1-6 Annex D
nn_cr_shell = 2 * np.pi * radius * thickness * sigma_x_rcr_new(thickness, radius, length)
# Return value
return nn_cr_shell
[docs]def sigma_x_rk(
thickness,
radius,
length,
f_y_k,
fab_quality=None,
flex_kapa=None
):
"""
Meridional design buckling stress.
Calculates the meridional buckling stress for a cylindrical shell
according to EN1993-1-6 [1].
Parameters
----------
thickness : float
[mm] Shell thickness
radius : float
[mm] Cylinder radius
length : float
[mm] Cylnder length
f_y_k : float
[MPa] Characteristic yield strength
fab_quality : str, optional
[_] Fabrication quality class. Accepts: 'fcA', 'fcB', 'fcC'
The three classes correspond to .006, .010 and .016 times the
width of a dimple on the shell.
Default = 'fcA', which implies excelent fabrication
gamma_m1 : int, optional
[_] Partial safety factor
Default = 1.1
Returns
-------
float
[MPa] Meridional buckling stress
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-6: Strength and stability of shell structures.
Brussels: CEN, 2006._
"""
# Default values
if fab_quality is None:
fab_quality = 'fcA'
if flex_kapa is None:
flex_kapa = 1.
# Check class acc. to EC3-1
t_class = tube_class(thickness, radius, f_y_k)
if t_class == 4:
# Fabrication quality class acc. to table D2
if fab_quality is 'fcA':
q_factor = 40.
elif fab_quality is 'fcB':
q_factor = 25.
elif fab_quality is 'fcC':
q_factor = 16.
else:
print('Invalid fabrication class input. Choose between \'fcA\', \'fcB\' and \'fcC\' ')
return
# Critical meridional stress, calculated on separate function
category = shell_length_category(radius, thickness, length)
sigma_cr = sigma_x_rcr(thickness, radius, length)
# Shell slenderness
lmda = np.sqrt(f_y_k / sigma_cr)
delta_w_k = (1. / q_factor) * np.sqrt(radius / thickness) * thickness
alpha = 0.62 / (1 + 1.91 * (delta_w_k / thickness) ** 1.44)
beta = 0.6
eta = 1.
if category == 2:
# For long cylinders, a formula is suggested for lambda, EC3-1-6 D1.2.2(4)
# Currently, the general form is used. to be fixed.
lmda_0 = 0.2
# lmda_0 = 0.2 + 0.1 * (sigma_e_M / sigma_e)
else:
lmda_0 = 0.2
lmda_p = np.sqrt(alpha / (1. - beta))
# Buckling reduction factor, chi
if lmda <= lmda_0:
chi_shell = 1.
elif lmda < lmda_p:
chi_shell = 1. - beta * ((lmda - lmda_0) / (lmda_p - lmda_0)) ** eta
else:
chi_shell = alpha / (lmda ** 2)
# Buckling stress
sigma_rk = chi_shell * f_y_k
else:
sigma_rk = f_y_k
# flex buckling
r_o = radius + thickness / 2.
r_i = radius - thickness / 2.
moi = np.pi * (r_o ** 4 - r_i ** 4) / 4.
area = 2 * np.pi * radius * thickness
lmda = lmbda_flex(length, area, moi, f_y_k, kapa_bc=flex_kapa)
chi = chi_flex(lmda, "c")
sigma_rd = sigma_rk * chi
# Return value
return sigma_rd
[docs]def sigma_x_rk_new(
thickness,
radius,
length,
f_y_k,
fab_quality=None,
flex_kapa=None
):
"""
Meridional characteristic buckling stress.
Calculates the characteristic meridional buckling stress for a cylindrical shell according to EN1993-1-6 [1].
Parameters
----------
thickness : float
[mm] Shell thickness
radius : float
[mm] Cylinder radius
length : float
[mm] Cylnder length
f_y_k : float
[MPa] Characteristic yield strength
fab_quality : str, optional
[_] Fabrication quality class. Accepts: 'fcA', 'fcB', 'fcC'
The three classes correspond to .006, .010 and .016 times the
width of a dimple on the shell.
Default = 'fcA', which implies excelent fabrication
gamma_m1 : int, optional
[_] Partial safety factor
Default = 1.2 (new suggestion from prEN draft)
Returns
-------
float
[MPa] Meridional buckling stress
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-6: Strength and stability of shell structures.
Brussels: CEN, 2006._
"""
# Default values
if fab_quality is None:
fab_quality = 'fcA'
if flex_kapa is None:
flex_kapa = 1.
t_class = tube_class(thickness, radius, f_y_k)
if t_class == 4:
# Fabrication quality class acc. to table D2
if fab_quality is 'fcA':
q_factor = 40.
elif fab_quality is 'fcB':
q_factor = 25.
elif fab_quality is 'fcC':
q_factor = 16.
else:
print('Invalid fabrication class input. Choose between \'fcA\', \'fcB\' and \'fcC\' ')
return
# Critical meridional stress, calculated on separate function
category = shell_length_category_new(radius, thickness, length)
sigma_cr = sigma_x_rcr_new(thickness, radius, length)
# Shell slenderness
lmda = np.sqrt(f_y_k / sigma_cr)
delta_w_k = (1. / q_factor) * np.sqrt(radius / thickness) * thickness
alpha_xG = 0.83
alpha_I = 0.06 + 0.93 / (1 + 2.7 * (delta_w_k / thickness)**0.85)
alpha = alpha_xG * alpha_I
#alpha = 0.62 / (1 + 1.91 * (delta_w_k / thickness) ** 1.44)
beta = 1 - 0.95 / (1 + 1.12 * (delta_w_k / thickness))
eta = 5.8 / (1 + 4.6 * (delta_w_k / thickness)**0.87)
chi_xh = 1.15
lambda_p = np.sqrt(alpha / (1 - beta))
if category == 2:
# For long cylinders, a formula is suggested for lambda under compression/bending , EC3-1-6 D1.2.2(4)
# Pure compression is assumed.
lmda_0 = 0.1 + (0 / 1.)
# lmda_0 = 0.2 + 0.1 * (sigma_e_M / sigma_e)
else:
lmda_0 = 0.2
lmda_p = np.sqrt(alpha / (1. - beta))
# Buckling reduction factor, chi
if lmda <= lmda_0:
chi_shell = chi_xh - (lmda / lmda_0)*(chi_xh - 1)
elif lmda < lmda_p:
chi_shell = 1. - beta * ((lmda - lmda_0) / (lmda_p - lmda_0)) ** eta
else:
chi_shell = alpha / (lmda ** 2)
# Buckling stress
sigma_rk = chi_shell * f_y_k
else:
sigma_rk = f_y_k
# flex buckling
r_o = radius + thickness / 2.
r_i = radius - thickness / 2.
moi = np.pi * (r_o ** 4 - r_i ** 4) / 4.
area = 2 * np.pi * radius * thickness
lmda = lmbda_flex(length, area, moi, f_y_k, kapa_bc=flex_kapa)
chi = chi_flex(lmda, "c")
sigma_rd = sigma_rk * chi
# Return value
return sigma_rd
[docs]def sigma_x_rd(
thickness,
radius,
length,
f_y_k,
fab_quality=None,
gamma_m1=None,
flex_kapa=None
):
"""
Meridional design buckling stress.
Calculates the meridional buckling stress for a cylindrical shell
according to EN 1993-1-6 [1]. Flexural buckling is also checked acc to EN 1993-1-1 [2]
Parameters
----------
thickness : float
[mm] Shell thickness
radius : float
[mm] Cylinder radius
length : float
[mm] Cylnder length
f_y_k : float
[MPa] Characteristic yield strength
fab_quality : str, optional
[_] Fabrication quality class. Accepts: 'fcA', 'fcB', 'fcC'
The three classes correspond to .006, .010 and .016 times the
width of a dimple on the shell.
Default = 'fcA', which implies excelent fabrication
gamma_m1 : int, optional
[_] Partial safety factor
Default = 1.1
Returns
-------
float
[MPa] Meridional buckling stress
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-6: Strength and stability of shell structures.
Brussels: CEN, 2006._
.. [2] Eurocode 3: Design of steel structures - Part 1-1: General rules and rules for buildings.
Brussels: CEN, 2005._
"""
# Default values
if fab_quality is None:
fab_quality = 'fcA'
if gamma_m1 is None:
# The default value 1.1 is as the lowest recommended on EC3-1-6 8.5.2(1) NOTE.
gamma_m1 = 1.1
else:
gamma_m1 = float(gamma_m1)
if flex_kapa is None:
flex_kapa = 1.
sigma_xx_rd = sigma_x_rk(
thickness,
radius,
length,
f_y_k,
fab_quality=fab_quality,
flex_kapa=flex_kapa
) / gamma_m1
return sigma_xx_rd
[docs]def sigma_x_rd_new(
thickness,
radius,
length,
f_y_k,
fab_quality=None,
gamma_m1=None,
flex_kapa=None
):
"""
Meridional design buckling stress.
Calculates the meridional buckling stress for a cylindrical shell
according to the new draft of EN1993-1-6 [1].
Parameters
----------
thickness : float
[mm] Shell thickness
radius : float
[mm] Cylinder radius
length : float
[mm] Cylnder length
f_y_k : float
[MPa] Characteristic yield strength
fab_quality : str, optional
[_] Fabrication quality class. Accepts: 'fcA', 'fcB', 'fcC'
The three classes correspond to .006, .010 and .016 times the
width of a dimple on the shell.
Default = 'fcA', which implies excelent fabrication
gamma_m1 : int, optional
[_] Partial safety factor
Default = 1.1
Returns
-------
float
[MPa] Meridional buckling stress
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-6 draft: Strength and stability of shell structures.
Brussels: CEN, 2006._
"""
# Default values
if fab_quality is None:
fab_quality = 'fcA'
if gamma_m1 is None:
# The default value 1.2 is as the lowest recommended on EC3-1-6 8.5.2(10) NOTE.
gamma_m1 = 1.2
else:
gamma_m1 = float(gamma_m1)
if flex_kapa is None:
flex_kapa = 1.
sigma_xx_rd = sigma_x_rk_new(
thickness,
radius,
length,
f_y_k,
fab_quality=fab_quality,
flex_kapa=flex_kapa
) / gamma_m1
return sigma_xx_rd
# OVERALL BUCKLING
[docs]def n_cr_flex(
length,
moi_y,
kapa_bc=None,
e_modulus=None
):
# Docstring
"""
Euler's critical load.
Calculates the critical load for flexural buckling of a given column.
A single direction is considered. If more directions are required
(e.g the two principal axes), the function has to be called multiple
times. For torsional mode critical load use n_cr_tor(), and for
flexural-torsional critical load use n_cr_flex_tor()
Parameters
----------
length : float
[mm] Column length.
moi_y : float
[mm^4] Moment of inertia.
kapa_bc : float, optional
[_] length correction for the effect of the boundary conditions.
Default = 1, which implies simply supported column.
e_modulus : float, optional
[MPa] Modulus of elasticity.
Default = 210000., typical value for steel.
Returns
-------
float
[N] Critical load.
"""
# default values
if kapa_bc is None:
kapa_bc = 1.
else:
kapa_bc = float(kapa_bc)
if e_modulus is None:
e_modulus = 210000.
else:
e_modulus = float(e_modulus)
# Euler's critical load
nn_cr_flex = (np.pi ** 2) * e_modulus * moi_y / (kapa_bc * length) ** 2
# Return the result
return nn_cr_flex
[docs]def n_cr_tor(
length,
area,
moi_y0,
moi_z0,
moi_torsion,
moi_warp,
y_0=None,
z_0=None,
e_modulus=None,
poisson=None,
):
# Docstring
"""
Torsional elastic critical load
Calculates the torsional elastic critical load for a hinged column.
The input values are refering to the principal axes. For flexural
buckling (Euler cases) use n_cr_flex. For the combined
flexural-torsional modes use n_cr_flex_tor.
Parameters
----------
length : float
[mm] Column length.
area : float
[mm^2] Cross-sectional area.
moi_y0 : float
[mm^4] Moment of inertia around `y`-axis.
`y`-axis on the centre of gravity but not necessarily principal.
moi_z0 : float
[mm^4] Moment of inertia around `z`-axis.
`z`-axis on the centre of gravity but not necessarily principal.
moi_torsion : float
[mm^4] Saint Venant constant.
moi_warp : float
[mm^6] Torsion constant.
y_0 : float, optional
[mm] Distance on `y`-axis of the shear center to the origin.
Default = 0, which implies symmetric profile
z_0 : float, optional
[mm] Distance on `z`-axis of the shear center to the origin.
Default = 0, which implies symmetric profile
e_modulus : float, optional
[MPa] Modulus of elasticity.
Default = 210000., general steel.
poisson : float, optional
[_] Young's modulus of elasticity.
Default = 0.3, general steel.
Returns
-------
float
[N] Flexural-torsional critical load.
Notes
-----
The torsional critical load is calculated as:
.. math:: N_{cr, tor} = {GJ + {\pi^2EI_w\over{L^2}}\over{r^2}}
Where:
:math:`E` : Elasticity modulus
:math:`G` : Shear modulus
:math:`J` : Torsional constant (Saint Venant)
:math:`I_w` : Warping constant
:math:`r^2=(moi_y + moi_z)/A + x_0^2 + y_0^2`
:math:`x_0, y_0` : Shear centre coordinates on the principal coordinate system
References
----------
..[1]N. S. Trahair, Flexural-torsional buckling of structures, vol. 6. CRC Press, 1993.
..[2]NS. Trahair, MA. Bradford, DA. Nethercot, and L. Gardner, The behaviour and design of steel structures to EC3, 4th edition. London; New York: Taylor & Francis, 2008.
"""
# default values
if y_0 is None:
y_0 = 0
else:
y_0 = float(y_0)
if z_0 is None:
z_0 = 0
else:
z_0 = float(z_0)
if e_modulus is None:
e_modulus = 210000.
else:
e_modulus = float(e_modulus)
if poisson is None:
poisson = 0.3
else:
poisson = float(poisson)
# Shear modulus
g_modulus = e_modulus / (2 * (1 + poisson))
# Polar radius of gyration.
i_pol = np.sqrt((moi_y0 + moi_z0) / area)
moi_zero = np.sqrt(i_pol ** 2 + y_0 ** 2 + z_0 ** 2)
# Calculation of critical torsional load.
nn_cr_tor = (1 / moi_zero ** 2) * (g_modulus * moi_torsion + (np.pi ** 2 * e_modulus * moi_warp / length ** 2))
# Return the result
return nn_cr_tor
[docs]def n_cr_flex_tor(
length,
area,
moi_y,
moi_z,
moi_yz,
moi_torsion,
moi_warp,
y_sc=None,
z_sc=None,
e_modulus=None,
poisson=None,
):
# Docstring
"""
Flexural-Torsional elastic critical load
Calculates the critical load for flexural-torsional buckling of a
column with hinged ends. The returned value is the minimum of the
the three flexural-torsional and the indepedent torsional mode, as
dictated in EN1993-1-1 6.3.1.4 [1]. (for further details, see Notes).
Parameters
----------
length : float
[mm] Column length.
area : float
[mm^2] Cross-sectional area.
moi_y : float
[mm^4] Moment of inertia around `y`-axis.
`y`-axis on the centre of gravity but not necessarily principal.
moi_z : float
[mm^4] Moment of inertia around `z`-axis.
`z`-axis on the centre of gravity but not necessarily principal.
moi_yz : float
[mm^4] Product of inertia.
moi_torsion : float
[mm^4] Saint Venant constant.
moi_warp : float
[mm^6] Torsion constant.
y_sc : float, optional
[mm] Distance on `y`-axis of the shear center to the origin.
Default = 0, which implies symmetric profile
z_sc : float, optional
[mm] Distance on `z`-axis of the shear center to the origin.
Default = 0, which implies symmetric profile
e_modulus : float, optional
[MPa] Modulus of elasticity.
Default = 210000., general steel.
poisson : float, optional
[_] Young's modulus of elasticity.
Default = 0.3, general steel.
Returns
-------
float
[N] Flexural-torsional critical load.
Notes
-----
The flexural-torsional critical loads are calculated as a combination
of the three independent overall buckling modes:
i) flexural around the major axis,
ii) flexural around the minor axis,
iii) Torsional buckling (around x-axis).
First, the cs-properties are described on the principal axes. Then
the three independent modes are calculated. The combined
flexural-torsional modes are calculated as the roots of a 3rd order
equation, as given in [1], [2]. The minimum of the torsional and the
three combined modes is returned (the two independent flexural modes
are not considered; for critical load of pure flexural mode use
'n_cr_flex').
References
----------
..[1]N. S. Trahair, Flexural-torsional buckling of structures, vol. 6. CRC Press, 1993.
..[2]NS. Trahair, MA. Bradford, DA. Nethercot, and L. Gardner, The behaviour and design of steel structures to EC3, 4th edition. London; New York: Taylor & Francis, 2008.
"""
# default values
if y_sc is None:
y_sc = 0
else:
y_sc = float(y_sc)
if z_sc is None:
z_sc = 0
else:
z_sc = float(z_sc)
if e_modulus is None:
e_modulus = 210000.
else:
e_modulus = float(e_modulus)
if poisson is None:
poisson = 0.3
else:
poisson = float(poisson)
# Angle of principal axes
if abs(moi_y - moi_z) < 1e-20:
theta = np.pi / 4
else:
theta = -np.arctan((2 * moi_yz) / (moi_y - moi_z)) / 2
# Distance of the rotation centre to the gravity centre on the
# principal axes coordinate system
y_0 = y_sc * np.cos(-theta) - z_sc * np.sin(-theta)
z_0 = z_sc * np.cos(-theta) + y_sc * np.sin(-theta)
# Moment of inertia around principal axes.
moi_y0 = (moi_y + moi_z) / 2 + np.sqrt(((moi_y - moi_z) / 2) ** 2 + moi_yz ** 2)
moi_z0 = (moi_y + moi_z) / 2 - np.sqrt(((moi_y - moi_z) / 2) ** 2 + moi_yz ** 2)
# Polar radius of gyration.
i_pol = np.sqrt((moi_y0 + moi_z0) / area)
moi_zero = np.sqrt(i_pol ** 2 + y_0 ** 2 + z_0 ** 2)
# Independent critical loads for flexural and torsional modes.
n_cr_max = (np.pi ** 2 * e_modulus * moi_y0) / (length ** 2)
n_cr_min = (np.pi ** 2 * e_modulus * moi_z0) / (length ** 2)
n_tor = n_cr_tor(
length,
area,
moi_y0,
moi_z0,
moi_torsion,
moi_warp=moi_warp,
y_0=y_0,
z_0=z_0,
e_modulus=e_modulus,
poisson=poisson
)
# Coefficients of the 3rd order equation for the critical loads
# The equation is in the form aaaa * N ^ 3 - bbbb * N ^ 2 + cccc * N - dddd
aaaa = moi_zero ** 2 - y_0 ** 2 - z_0 ** 2
bbbb = ((n_cr_max + n_cr_min + n_tor) * moi_zero ** 2) - (n_cr_min * y_0 ** 2) - (n_cr_max * z_0 ** 2)
cccc = moi_zero ** 2 * (n_cr_min * n_cr_max) + (n_cr_min * n_tor) + (n_tor * n_cr_max)
dddd = moi_zero ** 2 * n_cr_min * n_cr_max * n_tor
det_3 = (
4 * (-bbbb ** 2 + 3 * aaaa * cccc) ** 3 + (2 * bbbb ** 3 - 9 * aaaa * bbbb * cccc + 27 * aaaa ** 2 * dddd) ** 2
)
if det_3 < 0:
det_3 = -1. * det_3
cf = 1j
else:
cf = 1
# Critical load
# The following n_cr formulas are the roots of the 3rd order equation of the global critical load
n_cr_1 = bbbb / (3. * aaaa) - (2 ** (1. / 3) * (-bbbb ** 2 + 3 * aaaa * cccc)) / \
(3. * aaaa * (2 * bbbb ** 3 - 9 * aaaa * bbbb * cccc + 27 * aaaa ** 2 * dddd + \
(cf * np.sqrt(det_3))) ** (1. / 3)) + (
2 * bbbb ** 3 - 9 * aaaa * bbbb * cccc + 27 * aaaa ** 2 * dddd + \
(cf * np.sqrt(det_3))) ** (1. / 3) / (
3. * 2 ** (1. / 3) * aaaa)
n_cr_2 = bbbb / (3. * aaaa) + ((1 + (0 + 1j) * np.sqrt(3)) * (-bbbb ** 2 + 3 * aaaa * cccc)) / \
(3. * 2 ** (2. / 3) * aaaa * (
2 * bbbb ** 3 - 9 * aaaa * bbbb * cccc + 27 * aaaa ** 2 * dddd + \
(cf * np.sqrt(det_3))) ** (1. / 3)) - ((1 - (0 + 1j) * np.sqrt(3)) * \
(
2 * bbbb ** 3 - 9 * aaaa * bbbb * cccc + 27 * aaaa ** 2 * dddd + \
(cf * np.sqrt(det_3))) ** (1. / 3)) / (
6. * 2 ** (1. / 3) * aaaa)
n_cr_3 = bbbb / (3. * aaaa) + ((1 - (0 + 1j) * np.sqrt(3)) * (-bbbb ** 2 + 3 * aaaa * cccc)) / \
(3. * 2 ** (2. / 3) * aaaa * (
2 * bbbb ** 3 - 9 * aaaa * bbbb * cccc + 27 * aaaa ** 2 * dddd + \
(cf * np.sqrt(det_3))) ** (1. / 3)) - ((1 + (0 + 1j) * np.sqrt(3)) * \
(
2 * bbbb ** 3 - 9 * aaaa * bbbb * cccc + 27 * aaaa ** 2 * dddd + \
(cf * np.sqrt(det_3))) ** (1. / 3)) / (
6. * 2 ** (1. / 3) * aaaa)
# Lowest root is the critical load
nn_cr_flex_tor = min(abs(n_cr_1), abs(n_cr_2), abs(n_cr_3), n_tor)
# Return the critical load
return nn_cr_flex_tor
[docs]def lmbda_flex(
length,
area,
moi_y,
f_yield,
kapa_bc=None,
e_modulus=None,
):
# Docstring
"""
Flexural slenderness.
Calculates the slenderness of a columne under pure compression.
Euler's critical load is used.
Parameters
----------
length : float
[mm] Column length
area : float
[mm^2] Cross section area
moi_y : float
[mm^4] Moment of inertia
kapa_bc : float, optional
[_] length correction for the effect of the boundary conditions.
Default = 1, which implies simply supported column
e_modulus : float, optional
[MPa] Modulus of elasticity
Default = 210000., typical value for steel
f_yield : float, optional
[MPa] yield stress.
Default = 380., brcause this value was used extencively while the
function was being written. To be changed to 235.
Returns
-------
float
[_] Member slenderness
"""
# default values
if kapa_bc is None:
kapa_bc = 1.
else:
kapa_bc = float(kapa_bc)
if e_modulus is None:
e_modulus = 210000.
else:
e_modulus = float(e_modulus)
if f_yield is None:
f_yield = 235.
else:
f_yield = float(f_yield)
# Calculate Euler's critical load
n_cr = n_cr_flex(
length,
moi_y,
e_modulus=e_modulus,
kapa_bc=kapa_bc
)
# Flexural slenderness EN3-1-1 6.3.1.3 (1)
lmbda_flexx = np.sqrt(area * f_yield / n_cr)
# Return the result
return lmbda_flexx
[docs]def imp_factor(b_curve):
# Docstring
"""
Imperfection factor.
Returns the imperfection factor for a given buckling curve.
The values are taken from Table 6.1 of EN1993-1-1 [1]
Parameters
----------
b_curve : {'a0', 'a', 'b', 'c', 'd'}
[_] Name of the buckling curve as obtained from Table 6.2 of [1].
Returns
-------
float
[_] Imperfection factor.
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-1: General rules and
rules for buildings. Brussels: CEN, 2005.
"""
switcher = {
'a0': 0.13,
'a': 0.21,
'b': 0.34,
'c': 0.49,
'd': 0.76,
}
return switcher.get(b_curve, "nothing")
[docs]def chi_flex(
lmda,
#length,
#area,
#moi_y,
#f_yield,
b_curve,
):
# Docstring
"""
Flexural buckling reduction factor.
Claculates the reduction factor, chi, according to EN1993-1-1 6.3.1.2
Parameters
----------
length : float
[mm] Column length
area : float
[mm^2] Cross section area
moi_y : float
[mm^4] Moment of inertia
f_yield : float
[MPa] Yield stress.
b_curve : str
[_] Name of the buckling curve as obtained from Table 6.2 of [1].
Valid options are {'a0', 'a', 'b', 'c', 'd'}
kapa_bc : float, optional
[_] length correction for the effect of the boundary conditions.
Default = 1, which implies simply supported column
Returns
-------
float
[_] Reduction factor.
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-1: General rules and
rules for buildings. Brussels: CEN, 2005.
"""
#if kapa_bc is None:
# kapa_bc = 1.
#lmda = lmbda_flex(
# length,
# area,
# moi_y,
# f_yield,
# kapa_bc=kapa_bc,
# e_modulus=None,
#)
if lmda < 0.2:
chi = 1.
else:
alpha = imp_factor(b_curve)
phi = (1 + alpha * (lmda - 0.2) + lmda ** 2) / 2.
chi = 1 / (phi + np.sqrt(phi ** 2 - lmda ** 2))
if chi > 1.:
chi = 1.
return chi
[docs]def n_b_rd(
length,
area,
moi_y,
f_yield,
b_curve,
kapa_bc=None,
gamma_m1=None
):
# Docstring
"""
Flexural buckling resistance.
Verifies the resistance of a column against flexural buckling
according to EN1993-1-1 6.3.1.1.
Parameters
----------
length : float
[mm] Column length
area : float
[mm^2] Cross section area
moi_y : float
[mm^4] Moment of inertia
f_yield : float
[MPa] Yield stress.
b_curve : str
[_] Name of the buckling curve as obtained from Table 6.2 of [1].
Valid options are: {'a0', 'a', 'b', 'c', 'd'}
kapa_bc : float, optional
[_] Length correction for the effect of the boundary conditions.
Default = 1, which implies simply supported column
gamma_m1 : float, optional
[_] Partial safety factor.
Default = 1.
Returns
-------
float
[N] Buckling resistance.
References
----------
.. [1] Eurocode 3: Design of steel structures - Part 1-1: General rules and
rules for buildings. Brussels: CEN, 2005.
"""
if kapa_bc is None:
kapa_bc = 1.
if gamma_m1 is None:
gamma_m1 = 1.
lmda = lmbda_flex(length, area, moi_y, f_yield, kapa_bc=kapa_bc)
chi = chi_flex(lmda, b_curve)
nn_b_rd = area * f_yield * chi / gamma_m1
return nn_b_rd
# CONNECTIONS
[docs]def bolt_grade2stress(bolt_grade):
# Docstring
"""
Convert bolt grade to yield and ultimate stress.
Standard designation for bolt grade as a decimal is converted to yield and
ultimate stress values in MPa. In the standard bolt grade designation, the
integer part of the number represents the ultimate stress in MPa/100 and
the decimal part is the yield stress as a percentage of the ultimate (e.g
4.6 is f_u = 400, f_y = 400 * 0.6 = 240).
Parameters
----------
bolt_grade : float
Returns
-------
tuple : (f_ultimate, f_yield)
"""
# Calculation using divmod
f_ultimate = 100 * divmod(bolt_grade, 1)[0]
f_yield = round(f_ultimate * divmod(bolt_grade, 1)[1])
# Return values
return f_ultimate, f_yield
[docs]def stress_area(bolt_size, shear_threaded=None):
# Docstring
"""
Shear area of a bolt.
Returns the srea to be used for the calculation of shear resistance of a
bolt, either the gross cross-section of the
bolt (circle area) or the reduced area of the threaded part of the bolt.
Parameters
----------
bolt_size : float
Bolt's diameter.
shear_threaded : bool, optional
Designates if the shear plane is on the threaded portion or not.
Default in False, which implies shearing of the non-threaded portion
Returns
-------
float
Notes
-----
Currently, the threaded area is based on an average reduction of the shank
area. To be changed to analytic formula.
"""
# Default
if shear_threaded is None:
shear_threaded = False
# Calculate area
if shear_threaded:
a_shear = 0.784 * (np.pi * bolt_size ** 2 / 4)
else:
a_shear = np.pi * bolt_size ** 2 / 4
# Return
return a_shear
[docs]def f_v_rd(
bolt_size,
bolt_grade,
shear_threaded=None,
gamma_m2=None
):
"""
Bolt's shear resistance.
Calculates the shear resistance of single bolt for one shear plane as given
in table 3.4 of EC3-1-8.
Parameters
----------
bolt_size : float
Diameter of the non-threaded part (nominal bolt size e.g. M16 = 16)
bolt_grade : float
Bolt grade in standard designation format (see documentation of
bolt_grade2stress())
shear_threaded : bool, optional
Designates if the shear plane is on the threaded portion or not.
Default in False, which implies shearing of the non-threaded portion
gamma_m2 : float, optional
Safety factor.
Default value is 1.25
Returns
-------
float
"""
# Defaults
bolt_size = float(bolt_size)
if shear_threaded is None:
shear_threaded = False
if gamma_m2 is None:
gamma_m2 = 1.25
else:
gamma_m2 = float(gamma_m2)
# av coefficient
if shear_threaded and bolt_grade == (4.6 or 8.6):
a_v = 0.5
else:
a_v = 0.6
# Get ultimate stress for bolt
f_ub = bolt_grade2stress(bolt_grade)[0]
# Shear area
a_shear = stress_area(bolt_size, shear_threaded)
# Shear resistance
ff_v_rd = a_v * f_ub * a_shear / gamma_m2
# Return value
return ff_v_rd
[docs]def bolt_min_dist(d_0):
"""
Minimum bolt spacing.
:param d_0:
:return:
"""
e_1 = 1.2 * d_0
e_2 = 1.2 * d_0
e_3 = 1.5 * d_0
p_1 = 2.2 * d_0
p_2 = 2.4 * d_0
return e_1, e_2, e_3, p_1, p_2
[docs]def f_b_rd(bolt_size, bolt_grade, thickness, steel_grade, f_yield, distances, d_0):
"""
Connection bearing capacity.
Calculates the bearing capacity of a single bolt on a plate. The distances to the plate edges/other bolts are
described
:param bolt_size:
:param bolt_grade:
:param thickness:
:param steel_grade:
:param f_yield:
:param distances:
:param d_0:
:return:
"""
pass
[docs]def f_weld_perp(l_w, alpha, f_y, f_u, gamma_m2=None):
beta = 0.8
if f_y>=275:
beta = 0.85
if f_y>=355:
beta = 0.9
if f_y>=420:
beta = 0.9
if gamma_m2 is None:
gamma_m2 = 1.25
# Resistance for a 90 degree fillet weld under tension
f_rd = min([f_u * l_w * alpha * 2**0.5/ (2 * beta * gamma_m2), 0.9 * f_u * l_w * alpha / gamma_m2])
return f_rd
[docs]def f_weld_paral():
pass
[docs]def bolt2washer(m_bolt):
"""
Washer diameter.
Return the diameter of the washer for a given bolt diameter.
The calculation is based on a function derived from linear regression
on ENXXXXXXX[REF].
Parameters
----------
m_bolt : float
Bolt diameter
"""
d_washer = np.ceil(1.5893 * m_bolt + 5.1071)
return d_washer
[docs]def fabclass_2_umax(fab_class):
"""
Maximum displacement for a given fabrication class acc. to EC3-1-6.
Parameters
----------
fab_class : {"fcA", "fcB", "fcC"}
"""
# Assign imperfection amplitude, u_max acc. to the fabrication class
if fab_class is 'fcA':
u_max = 0.006
elif fab_class is 'fcB':
u_max = 0.010
else:
u_max = 0.016
# Return values
return u_max
[docs]def mean_list(numbers):
"""
Mean value.
Calculate the average for a list of numbers.
Parameters
----------
numbers : list
Attributes
----------
Notes
-----
References
----------
"""
return float(sum(numbers)) / max(len(numbers), 1)
[docs]def test_rhs():
width = 100.
height = 100.
thick = 3.
material = Material(210000, 0.3, 500)
cs = CsRHS.from_n_prc(
width,
height,
thick,
material=material,
n_prc=-0.95,
)
return(cs)
# SAFETY FACTORS
[docs]def gamma_m(m_value, nominal=True, country=None):
"""
Get a safety factor.
The function returns the recommended value of a requested safety factor, as
listed in table 4.1 of EN1993-1-8.
Parameters
----------
m_value : str {'0', '1', '2', '3', '3,ser', '4', '5', '6,ser', '7', 'u'}
Description of the requested safety factor. The input can be given as a
number instead of a string, except for the serviceability cases.
nominal : bool, optional
If False, the function will always return 1. It's useful for providing
the ability to cancel all safety factors throughout calculations.
Default value is True.
country : str, optional
The 2 letter country code for the requested safety factor.
"""
# All the recommended values
if country == "NO":
gamma_table = {
'0':1.05, # Value given in NA
'1':1.05, # Value given in NA
'2':1.25,
'3':1.25,
'3,ser':1.1,
'4':1.00,
'5':1.00,
'6,ser':1.00,
'7':1.1,
'u':1.1
}
else:
gamma_table = {
'0':1.00, # Recommended value from EN1993-1-1
'1':1.00, # Recommended value from EN1993-1-1
'2':1.25,
'3':1.25,
'3,ser':1.1,
'4':1.00,
'5':1.00,
'6,ser':1.00,
'7':1.1,
'u':1.1
}
# Cancel the safety factor if measured values are used
if nominal:
gamma = gamma_table[str(m_value)]
else:
gamma = 1.
return(gamma)
# temp function to test the CsRHS.from_n_prc method.
[docs]def test_csnprc(n_prc):
material_col = Material(
210000.,
0.3,
537.82,
614
)
chk = CsRHS.from_n_prc(
100.1,
100.1,
3.3,
material = material_col,
n_prc = n_prc,
negative_moment = False,
force_gross = False
)
print(chk.m_n_rd, chk.cs_class, chk.n_rd, chk.n_prc, chk.n_rd*chk.n_prc)
return(chk)
[docs]def test_csmprc(m_prc):
material_col = Material(
210000.,
0.3,
537.82,
614
)
chk = CsRHS.from_m_prc(
100.1,
100.1,
3.3,
material = material_col,
m_prc = m_prc,
tensile_axial = False,
force_gross = False
)
print(chk.m_n_rd, chk.cs_class, chk.n_rd, chk.n_prc, chk.n_rd*chk.n_prc)
return(chk)