Examples

The following examples are located in the sectionproperties.examples package.

Simple Example

The following example calculates the geometric, warping and plastic properties of a 50 mm diameter circle. The circle is discretised with 64 points and a mesh size of 2.5 mm².

The geometry and mesh are plotted, and the mesh information printed to the terminal before the analysis is carried out. Detailed time information is printed to the terminal during the cross-section analysis stage. Once the analysis is complete, the cross-section properties are printed to the terminal. The centroidal axis second moments of area and torsion constant are saved to variables and it is shown that, for a circle, the torsion constant is equal to the sum of the second moments of area:

import sectionproperties.pre.sections as sections
from sectionproperties.analysis.cross_section import CrossSection

# create a 50 diameter circle discretised by 64 points
geometry = sections.CircularSection(d=50, n=64)
geometry.plot_geometry()  # plot the geometry

# create a mesh with a mesh size of 2.5
mesh = geometry.create_mesh(mesh_sizes=[2.5])

section = CrossSection(geometry, mesh)  # create a CrossSection object
section.display_mesh_info()  # display the mesh information
section.plot_mesh()  # plot the generated mesh

# perform a geometric, warping and plastic analysis, displaying the time info
section.calculate_geometric_properties(time_info=True)
section.calculate_warping_properties(time_info=True)
section.calculate_plastic_properties(time_info=True)

# print the results to the terminal
section.display_results()

# get the second moments of area and the torsion constant
(ixx_c, iyy_c, ixy_c) = section.get_ic()
j = section.get_j()

# print the sum of the second moments of area and the torsion constant
print("Ixx + Iyy = {0:.3f}".format(ixx_c + iyy_c))
print("J = {0:.3f}".format(j))

The following plots are generated by the above example:

Circular section geometry.

Mesh generated from the above geometry.

The following is printed to the terminal:

Mesh Statistics:
--2562 nodes
--1247 elements
--1 region

--Calculating geometric section properties...
----completed in 0.765906 seconds---

--Assembing 2562x2562 stiffness matrix and load vector...
----completed in 1.107300 seconds---
--Solving for the warping function using the direct solver...
----completed in 0.023138 seconds---
--Computing the torsion constant...
----completed in 0.000254 seconds---
--Assembling shear function load vectors...
----completed in 1.150968 seconds---
--Solving for the shear functions using the direct solver...
----completed in 0.136840 seconds---
--Assembling shear centre and warping moment integrals...
----completed in 0.688029 seconds---
--Calculating shear centres...
----completed in 0.000054 seconds---
--Assembling shear deformation coefficients...
----completed in 1.184013 seconds---
--Assembling monosymmetry integrals...
----completed in 0.784923 seconds---

--Calculating plastic properties...
----completed in 0.669841 seconds---

Section Properties:
A      = 1.960343e+03
Qx     = 1.900702e-13
Qy     = 3.451461e-12
cx     = 1.760642e-15
cy     = 9.695762e-17
Ixx_g  = 3.058119e+05
Iyy_g  = 3.058119e+05
Ixy_g  = -2.785328e-12
Ixx_c  = 3.058119e+05
Iyy_c  = 3.058119e+05
Ixy_c  = -2.785328e-12
Zxx+   = 1.223248e+04
Zxx-   = 1.223248e+04
Zyy+   = 1.223248e+04
Zyy-   = 1.223248e+04
rx     = 1.248996e+01
ry     = 1.248996e+01
phi    = 0.000000e+00
I11_c  = 3.058119e+05
I22_c  = 3.058119e+05
Z11+   = 1.223248e+04
Z11-   = 1.223248e+04
Z22+   = 1.223248e+04
Z22-   = 1.223248e+04
r11    = 1.248996e+01
r22    = 1.248996e+01
J      = 6.116232e+05
Iw     = 4.700106e-02
x_se   = -8.788834e-06
y_se   = -2.644033e-06
x_st   = -8.788834e-06
y_st   = -2.644033e-06
x1_se  = -8.788834e-06
y2_se  = -2.644033e-06
A_sx   = 1.680296e+03
A_sy   = 1.680296e+03
A_s11  = 1.680296e+03
A_s22  = 1.680296e+03
betax+ = -5.288066e-06
betax- = 5.288066e-06
betay+ = -1.757767e-05
betay- = 1.757767e-05
beta11+= -5.288066e-06
beta11-= 5.288066e-06
beta22+= -1.757767e-05
beta22-= 1.757767e-05
x_pc   = 5.313355e-15
y_pc   = 3.649671e-15
Sxx    = 2.078317e+04
Syy    = 2.078317e+04
SF_xx+ = 1.699016e+00
SF_xx- = 1.699016e+00
SF_yy+ = 1.699016e+00
SF_yy- = 1.699016e+00
x11_pc = 5.313355e-15
y22_pc = 3.649671e-15
S11    = 2.078317e+04
S22    = 2.078317e+04
SF_11+ = 1.699016e+00
SF_11- = 1.699016e+00
SF_22+ = 1.699016e+00
SF_22- = 1.699016e+00

Ixx + Iyy = 611623.837
J = 611623.214

Creating a Nastran Section

The following example demonstrates how to create a cross-section defined in a Nastran-based finite element analysis program. The following creates a HAT1 cross-section and calculates the geometric, warping and plastic properties. The HAT1 cross-section is meshed with a maximum elemental area of 0.005.

The geometry and mesh are plotted, and the mesh information printed to the terminal before the analysis is carried out. Detailed time information is printed to the terminal during the cross-section analysis stage. Once the analysis is complete, the cross-section properties are printed to the terminal. The centroidal axis second moments of area and torsion constant are saved to variables and it is shown that, for non-circular sections, the torsion constant is not equal to the sum of the second moments of area:

import sectionproperties.pre.nastran_sections as nsections
from sectionproperties.analysis.cross_section import CrossSection

# create a HAT1 section
geometry = nsections.HAT1Section(DIM1=4.0, DIM2=2.0, DIM3=1.5, DIM4=0.1875, DIM5=0.375)
geometry.plot_geometry()  # plot the geometry

# create a mesh with a maximum elemental area of 0.005
mesh = geometry.create_mesh(mesh_sizes=[0.005])

section = CrossSection(geometry, mesh)  # create a CrossSection object
section.display_mesh_info()  # display the mesh information
section.plot_mesh()  # plot the generated mesh`

# perform a geometric, warping and plastic anaylsis, displaying the time info
section.calculate_geometric_properties(time_info=True)
section.calculate_warping_properties(time_info=True)
section.calculate_plastic_properties(time_info=True)

# print the results to the terminal
section.display_results()

# get the second moments of area and the torsion constant
(ixx_c, iyy_c, ixy_c) = section.get_ic()
j = section.get_j()

# print the sum of the second moments of area and the torsion constant
print("Ixx + Iyy = {0:.3f}".format(ixx_c + iyy_c))
print("J = {0:.3f}".format(j))

The following plots are generated by the above example:

Circular section geometry.

Mesh generated from the above geometry.

The following is printed to the terminal:

Mesh Statistics:
--2038 nodes
--926 elements
--2 regions

--Calculating geometric section properties...
----completed in 0.367074 seconds---

--Assembing 2038x2038 stiffness matrix and load vector...
----completed in 0.515934 seconds---
--Solving for the warping function using the direct solver...
----completed in 0.005604 seconds---
--Computing the torsion constant...
----completed in 0.000104 seconds---
--Assembling shear function load vectors...
----completed in 0.525532 seconds---
--Solving for the shear functions using the direct solver...
----completed in 0.064247 seconds---
--Assembling shear centre and warping moment integrals...
----completed in 0.331969 seconds---
--Calculating shear centres...
----completed in 0.000043 seconds---
--Assembling shear deformation coefficients...
----completed in 0.511631 seconds---
--Assembling monosymmetry integrals...
----completed in 0.389498 seconds---

--Calculating plastic properties...
----completed in 0.131321 seconds---

Section Properties:
A    = 2.789062e+00
Qx   = 1.626709e+00
Qy   = -1.424642e-16
cx   = -5.107959e-17
cy   = 5.832458e-01
Ixx_g        = 1.935211e+00
Iyy_g        = 3.233734e+00
Ixy_g        = -1.801944e-16
Ixx_c        = 9.864400e-01
Iyy_c        = 3.233734e+00
Ixy_c        = -9.710278e-17
Zxx+         = 6.962676e-01
Zxx-         = 1.691294e+00
Zyy+         = 1.616867e+00
Zyy-         = 1.616867e+00
rx   = 5.947113e-01
ry   = 1.076770e+00
phi  = -9.000000e+01
I11_c        = 3.233734e+00
I22_c        = 9.864400e-01
Z11+         = 1.616867e+00
Z11-         = 1.616867e+00
Z22+         = 1.691294e+00
Z22-         = 6.962676e-01
r11  = 1.076770e+00
r22  = 5.947113e-01
J    = 9.878443e-01
Iw   = 1.160810e-01
x_se         = 4.822719e-05
y_se         = 4.674792e-01
x_st         = 4.822719e-05
y_st         = 4.674792e-01
x1_se        = 1.157666e-01
y2_se        = 4.822719e-05
A_sx         = 1.648312e+00
A_sy         = 6.979733e-01
A_s11        = 6.979733e-01
A_s22        = 1.648312e+00
betax+       = -2.746928e-01
betax-       = 2.746928e-01
betay+       = 9.645438e-05
betay-       = -9.645438e-05
beta11+      = 9.645438e-05
beta11-      = -9.645438e-05
beta22+      = 2.746928e-01
beta22-      = -2.746928e-01
x_pc         = -5.107959e-17
y_pc         = 3.486328e-01
Sxx  = 1.140530e+00
Syy  = 2.603760e+00
SF_xx+       = 1.638062e+00
SF_xx-       = 6.743533e-01
SF_yy+       = 1.610373e+00
SF_yy-       = 1.610373e+00
x11_pc       = -3.671369e-17
y22_pc       = 3.486328e-01
S11  = 2.603760e+00
S22  = 1.140530e+00
SF_11+       = 1.610374e+00
SF_11-       = 1.610374e+00
SF_22+       = 6.743539e-01
SF_22-       = 1.638064e+00

Ixx + Iyy = 4.220
J = 0.988

Creating Custom Geometry

The following example demonstrates how geometry objects can be created from a list of points, facets, holes and control points. An straight angle section with a plate at its base is created from a list of points and facets. The bottom plate is assigned a separate control point meaning two discrete regions are created. Creating separate regions allows the user to control the mesh size in each region and assign material properties to different regions. The geometry is cleaned to remove the overlapping facet at the junction of the angle and the plate. A geometric, warping and plastic analysis is then carried out.

The geometry and mesh are plotted before the analysis is carried out. Once the analysis is complete, a plot of the various calculated centroids is generated:

import sectionproperties.pre.sections as sections
from sectionproperties.analysis.cross_section import CrossSection

# define parameters for the angle section
a = 1
b = 2
t = 0.1

# build the lists of points, facets, holes and control points
points = [[-t/2, -2*a], [t/2, -2*a], [t/2, -t/2], [a, -t/2], [a, t/2],
          [-t/2, t/2], [-b/2, -2*a], [b/2, -2*a], [b/2, -2*a-t],
          [-b/2, -2*a-t]]
facets = [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0], [6, 7], [7, 8],
          [8, 9], [9, 6]]
holes = []
control_points = [[0, 0], [0, -2*a-t/2]]

# create the custom geometry object
geometry = sections.CustomSection(points, facets, holes, control_points)
geometry.clean_geometry()  # clean the geometry
geometry.plot_geometry()  # plot the geometry

# create the mesh - use a smaller refinement for the angle region
mesh = geometry.create_mesh(mesh_sizes=[0.0005, 0.001])

# create a CrossSection object
section = CrossSection(geometry, mesh)
section.plot_mesh()  # plot the generated mesh

# perform a geometric, warping and plastic analysis
section.calculate_geometric_properties()
section.calculate_warping_properties()
section.calculate_plastic_properties()

# plot the centroids
section.plot_centroids()

The following plots are generated by the above example:

Plot of the generated geometry object.

Mesh generated from the above geometry.

Plot of the centroids and the principal axis.

Creating a Merged Section

The following example demonstrates how to merge multiple geometry objects into a single geometry object. A 150x100x6 RHS is modelled with a solid 50x50 triangular section on its top and a 100x100x6 EA section on its right side. The three geometry objects are merged together using the MergedSection class. The order of the geometry objects in the list that is passed into the constructor of the MergedSection class is important, as this same order relates to specifying mesh sizes and material properties.

Once the geometry has been merged, it is vital to clean the geometry to remove any artefacts that may impede the meshing algorithm. A mesh is created with a mesh size of 2.5 mm² for the RHS (first in section_list), 5 mm² for the triangle (second in section_list) and 3 mm² for the angle (last in section_list).

The geometry and mesh are plotted, and the mesh information printed to the terminal before the analysis is carried out. Detailed time information is printed to the terminal during the cross-section analysis stage. Once the analysis is complete, the centroids are plotted:

import sectionproperties.pre.sections as sections
from sectionproperties.analysis.cross_section import CrossSection

# create a 150x100x6 RHS
rhs = sections.Rhs(d=150, b=100, t=6, r_out=15, n_r=8)

# create a triangular section on top of the RHS
points = [[0, 0], [50, 0], [25, 50]]
facets = [[0, 1], [1, 2], [2, 0]]
holes = []
control_points = [[25, 25]]
triangle = sections.CustomSection(points, facets, holes, control_points,
                                  shift=[25, 150])

# create a 100x100x6 EA on the right of the RHS
angle = sections.AngleSection(d=100, b=100, t=6, r_r=8, r_t=5, n_r=8,
                              shift=[100, 25])

# create a list of the sections to be merged
section_list = [rhs, triangle, angle]

# merge the three sections into one geometry object
geometry = sections.MergedSection(section_list)

# clean the geometry - print cleaning information to the terminal
geometry.clean_geometry(verbose=True)
geometry.plot_geometry()  # plot the geometry

# create a mesh - use a mesh size of 2.5 for the RHS, 5 for the triangle and
# 3 for the angle
mesh = geometry.create_mesh(mesh_sizes=[2.5, 5, 3])

# create a CrossSection object
section = CrossSection(geometry, mesh)
section.display_mesh_info()  # display the mesh information
section.plot_mesh()  # plot the generated mesh

# perform a geometric, warping and plastic analysis, displaying the time info
# and the iteration info for the plastic analysis
section.calculate_geometric_properties(time_info=True)
section.calculate_warping_properties(time_info=True)
section.calculate_plastic_properties(time_info=True, verbose=True)

# plot the centroids
section.plot_centroids()

The following plots are generated by the above example:

Plot of the generated geometry object.

Mesh generated from the above geometry.

Plot of the centroids and the principal axis.

The following is printed to the terminal:

Removed overlapping facets... Rebuilt with points: [30, 67, 93, 32]
Removed overlapping facets... Rebuilt with points: [46, 65, 64, 48]
Mesh Statistics:
--6053 nodes
--2755 elements
--3 regions

--Calculating geometric section properties...
----completed in 1.730845 seconds---

--Assembing 6053x6053 stiffness matrix and load vector...
----completed in 2.793801 seconds---
--Solving for the warping function using the direct solver...
----completed in 0.021323 seconds---
--Computing the torsion constant...
----completed in 0.000316 seconds---
--Assembling shear function load vectors...
----completed in 2.552404 seconds---
--Solving for the shear functions using the direct solver...
----completed in 0.604847 seconds---
--Assembling shear centre and warping moment integrals...
----completed in 1.578075 seconds---
--Calculating shear centres...
----completed in 0.000068 seconds---
--Assembling shear deformation coefficients...
----completed in 2.438405 seconds---

--Calculating plastic properties...
---x-axis plastic centroid calculation converged at 1.66608e+00 in 7 iterations.
---y-axis plastic centroid calculation converged at -5.83761e+00 in 10 iterations.
---11-axis plastic centroid calculation converged at -1.43134e+00 in 7 iterations.
---22-axis plastic centroid calculation converged at -1.21319e+01 in 9 iterations.
----completed in 2.710146 seconds---

Mirroring and Rotating Geometry

The following example demonstrates how geometry objects can be mirrored and rotated. A 200PFC and 150PFC are placed back-to-back by using the mirror_section() method and are rotated counter-clockwise by 30 degrees by using the rotate_section() method. The geometry is cleaned to ensure there are no overlapping facets along the junction between the two PFCs. A geometric, warping and plastic analysis is then carried out.

The geometry and mesh are plotted, and the mesh information printed to the terminal before the analysis is carried out. Detailed time information is printed to the terminal during the cross-section analysis stage and iteration information printed for the plastic analysis. Once the analysis is complete, a plot of the various calculated centroids is generated:

import sectionproperties.pre.sections as sections
from sectionproperties.analysis.cross_section import CrossSection

# create a 200PFC and a 150PFC
pfc1 = sections.PfcSection(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
pfc2 = sections.PfcSection(d=150, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8,
                           shift=[0, 26.5])

# mirror the 200 PFC about the y-axis
pfc1.mirror_section(axis='y', mirror_point=[0, 0])

# merge the pfc sections
geometry = sections.MergedSection([pfc1, pfc2])

# rotate the geometry counter-clockwise by 30 degrees
geometry.rotate_section(angle=30)

# clean the geometry - print cleaning information to the terminal
geometry.clean_geometry(verbose=True)
geometry.plot_geometry()  # plot the geometry

# create a mesh - use a mesh size of 5 for the 200PFC and 4 for the 150PFC
mesh = geometry.create_mesh(mesh_sizes=[5, 4])

# create a CrossSection object
section = CrossSection(geometry, mesh)
section.display_mesh_info()  # display the mesh information
section.plot_mesh()  # plot the generated mesh

# perform a geometric, warping and plastic analysis, displaying the time info
# and the iteration info for the plastic analysis
section.calculate_geometric_properties(time_info=True)
section.calculate_warping_properties(time_info=True)
section.calculate_plastic_properties(time_info=True, verbose=True)

# plot the centroids
section.plot_centroids()

The following plots are generated by the above example:

Plot of the generated geometry object.

Mesh generated from the above geometry.

Plot of the centroids and the principal axis.

The following is printed to the terminal:

Removed overlapping facets... Rebuilt with points: [21, 43, 22, 0]
Mesh Statistics:
--4841 nodes
--2152 elements
--2 regions

--Calculating geometric section properties...
----completed in 1.350236 seconds---

--Assembing 4841x4841 stiffness matrix and load vector...
----completed in 2.002365 seconds---
--Solving for the warping function using the direct solver...
----completed in 0.013307 seconds---
--Computing the torsion constant...
----completed in 0.000222 seconds---
--Assembling shear function load vectors...
----completed in 1.910170 seconds---
--Solving for the shear functions using the direct solver...
----completed in 0.623121 seconds---
--Assembling shear centre and warping moment integrals...
----completed in 1.163591 seconds---
--Calculating shear centres...
----completed in 0.000059 seconds---
--Assembling shear deformation coefficients...
----completed in 1.831169 seconds---

--Calculating plastic properties...
---x-axis plastic centroid calculation converged at 2.77651e+00 in 9 iterations.
---y-axis plastic centroid calculation converged at 3.02247e+00 in 5 iterations.
---11-axis plastic centroid calculation converged at -2.41585e-13 in 3 iterations.
---22-axis plastic centroid calculation converged at 6.10669e-01 in 5 iterations.
----completed in 0.860817 seconds---

Performing a Stress Analysis

The following example demonstrates how a stress analysis can be performed on a cross-section. A 150x100x6 RHS is modelled on its side with a maximum mesh area of 2 mm². The pre-requisite geometric and warping analyses are performed before two separate stress analyses are undertaken. The first combines bending and shear about the x-axis with a torsion moment and the second combines bending and shear about the y-axis with a torsion moment.

After the analysis is performed, various plots of the stresses are generated:

import sectionproperties.pre.sections as sections
from sectionproperties.analysis.cross_section import CrossSection

# create a 150x100x6 RHS on its side
geometry = sections.Rhs(d=100, b=150, t=6, r_out=15, n_r=8)

# create a mesh with a maximum area of 2
mesh = geometry.create_mesh(mesh_sizes=[2])

# create a CrossSection object
section = CrossSection(geometry, mesh)

# perform a geometry and warping analysis
section.calculate_geometric_properties()
section.calculate_warping_properties()

# perform a stress analysis with Mx = 5 kN.m; Vx = 10 kN and Mzz = 3 kN.m
case1 = section.calculate_stress(Mxx=5e6, Vx=10e3, Mzz=3e6)

# perform a stress analysis with My = 15 kN.m; Vy = 30 kN and Mzz = 1.5 kN.m
case2 = section.calculate_stress(Myy=15e6, Vy=30e3, Mzz=1.5e6)

case1.plot_stress_m_zz(pause=False)  # plot the bending stress for case1
case1.plot_vector_mzz_zxy(pause=False)  # plot the torsion vectors for case1
case2.plot_stress_v_zxy(pause=False)  # plot the shear stress for case1
case1.plot_stress_vm(pause=False)  # plot the von mises stress for case1
case2.plot_stress_vm()  # plot the von mises stress for case2

The following plots are generated by the above example:

Contour plot of the bending stress for case 1.

Vector plot of the torsion stress for case 1.

Contour plot of the shear stress for case 2.

Contour plot of the von Mises stress for case 1.

Contour plot of the von Mises stress for case 2.

Creating a Composite Cross-Section

The following example demonstrates how to create a composite cross-section by assigning different material properties to various regions of the mesh. A steel 310UB40.4 is modelled with a 50Dx600W timber panel placed on its top flange.

The geometry and mesh are plotted, and the mesh information printed to the terminal before the analysis is carried out. All types of cross-section analyses are carried out, with an axial force, bending moment and shear force applied during the stress analysis. Once the analysis is complete, the cross-section properties are printed to the terminal and a plot of the centroids and cross-section stresses generated:

import sectionproperties.pre.sections as sections
from sectionproperties.pre.pre import Material
from sectionproperties.analysis.cross_section import CrossSection

# create material properties
steel = Material(name='Steel', elastic_modulus=200e3, poissons_ratio=0.3,
                 yield_strength=500, color='grey')
timber = Material(name='Timber', elastic_modulus=8e3, poissons_ratio=0.35,
                  yield_strength=20, color='burlywood')

# create 310UB40.4
ub = sections.ISection(d=304, b=165, t_f=10.2, t_w=6.1, r=11.4, n_r=8)

# create timber panel on top of the UB
panel = sections.RectangularSection(d=50, b=600, shift=[-217.5, 304])

# merge the two sections into one geometry object
geometry = sections.MergedSection([ub, panel])
geometry.clean_geometry()  # clean the geometry
geometry.plot_geometry()  # plot the geometry

# create a mesh - use a mesh size of 5 for the UB, 20 for the panel
mesh = geometry.create_mesh(mesh_sizes=[5, 20])

# create a CrossSection object - take care to list the materials in the same
# order as entered into the MergedSection
section = CrossSection(geometry, mesh, materials=[steel, timber])
section.display_mesh_info()  # display the mesh information

# plot the mesh with coloured materials and a line transparency of 0.5
section.plot_mesh(materials=True, alpha=0.5)

# perform a geometric, warping and plastic analysis
section.calculate_geometric_properties(time_info=True)
section.calculate_warping_properties(time_info=True)
section.calculate_plastic_properties(time_info=True, verbose=True)

# perform a stress analysis with N = 100 kN, Mxx = 120 kN.m and Vy = 75 kN
stress_post = section.calculate_stress(N=-100e3, Mxx=-120e6, Vy=-75e3,
                                       time_info=True)

# print the results to the terminal
section.display_results()

# plot the centroids
section.plot_centroids()

stress_post.plot_stress_n_zz(pause=False)  # plot the axial stress
stress_post.plot_stress_m_zz(pause=False)  # plot the bending stress
stress_post.plot_stress_v_zxy()  # plot the shear stress

The following plots are generated by the above example:

Plot of the generated geometry object.

Mesh generated from the above geometry.

Plot of the centroids and the principal axis.

Contour plot of the axial stress.

Contour plot of the bending stress.

Contour plot of the shear stress.

The following is printed to the terminal:

Mesh Statistics:
--8972 nodes
--4189 elements
--2 regions

--Calculating geometric section properties...
----completed in 2.619151 seconds---

--Assembing 8972x8972 stiffness matrix and load vector...
----completed in 4.814592 seconds---
--Solving for the warping function using the direct solver...
----completed in 0.032710 seconds---
--Computing the torsion constant...
----completed in 0.000281 seconds---
--Assembling shear function load vectors...
----completed in 3.648590 seconds---
--Solving for the shear functions using the direct solver...
----completed in 0.073731 seconds---
--Assembling shear centre and warping moment integrals...
----completed in 2.288843 seconds---
--Calculating shear centres...
----completed in 0.000064 seconds---
--Assembling shear deformation coefficients...
----completed in 3.597728 seconds---
--Assembling monosymmetry integrals...
----completed in 2.519333 seconds---

--Calculating plastic properties...
d = -185.13088495027134; f_norm = 1.0
d = 168.86911504972866; f_norm = -1.0
d = -8.130884950271337; f_norm = 0.1396051884674814
d = 13.552166872240885; f_norm = 0.0983423820518053
d = 60.60270845168385; f_norm = 0.008805290832496546
d = 64.90008929872263; f_norm = 0.0006273832465500235
d = 65.22746216923525; f_norm = 4.393296044849056e-06
d = 65.22976962543858; f_norm = 2.2112746988930985e-09
d = 65.2298027403234; f_norm = -6.080628821651535e-08
---x-axis plastic centroid calculation converged at 6.52298e+01 in 8 iterations.
d = -300.0; f_norm = -1.0
d = 300.0; f_norm = 1.0
d = 0.0; f_norm = 2.1790636349700628e-16
d = -5e-07; f_norm = -4.7730935851751974e-08
---y-axis plastic centroid calculation converged at 0.00000e+00 in 3 iterations.
d = -185.13088495027134; f_norm = 1.0
d = 168.86911504972866; f_norm = -1.0
d = -8.130884950271337; f_norm = 0.1396051884674814
d = 13.552166872240885; f_norm = 0.0983423820518053
d = 60.60270845168385; f_norm = 0.008805290832496546
d = 64.90008929872263; f_norm = 0.0006273832465500235
d = 65.22746216923525; f_norm = 4.393296044849056e-06
d = 65.22976962543858; f_norm = 2.2112746988930985e-09
d = 65.2298027403234; f_norm = -6.080628821651535e-08
---11-axis plastic centroid calculation converged at 6.52298e+01 in 8 iterations.
d = -300.0; f_norm = -1.0
d = 300.0; f_norm = 1.0
d = 0.0; f_norm = 2.1790636349700628e-16
d = -5e-07; f_norm = -4.7730935851751974e-08
---22-axis plastic centroid calculation converged at 0.00000e+00 in 3 iterations.
----completed in 0.794056 seconds---

--Calculating cross-section stresses...
----completed in 4.240446 seconds---

Section Properties:
A      = 3.521094e+04
E.A    = 1.282187e+09
E.Qx   = 2.373725e+11
E.Qy   = 1.057805e+11
cx     = 8.250000e+01
cy     = 1.851309e+02
E.Ixx_g= 6.740447e+13
E.Iyy_g= 1.745613e+13
E.Ixy_g= 1.958323e+13
E.Ixx_c= 2.345949e+13
E.Iyy_c= 8.729240e+12
E.Ixy_c= -7.421875e-02
E.Zxx+ = 1.389212e+11
E.Zxx- = 1.267184e+11
E.Zyy+ = 2.909747e+10
E.Zyy- = 2.909747e+10
rx     = 1.352644e+02
ry     = 8.251112e+01
phi    = 0.000000e+00
E.I11_c= 2.345949e+13
E.I22_c= 8.729240e+12
E.Z11+ = 1.389212e+11
E.Z11- = 1.267184e+11
E.Z22+ = 2.909747e+10
E.Z22- = 2.909747e+10
r11    = 1.352644e+02
r22    = 8.251112e+01
G.J    = 1.439379e+11
G.Iw   = 2.554353e+16
x_se   = 8.250071e+01
y_se   = 2.863400e+02
x_st   = 8.250070e+01
y_st   = 2.857074e+02
x1_se  = 7.063407e-04
y2_se  = 1.012091e+02
A_sx   = 1.104723e+04
A_sy   = 1.021183e+04
A_s11  = 1.104723e+04
A_s22  = 1.021183e+04
betax+ = 2.039413e+02
betax- = -2.039413e+02
betay+ = 1.412681e-03
betay- = -1.412681e-03
beta11+= 2.039413e+02
beta11-= -2.039413e+02
beta22+= 1.412681e-03
beta22-= -1.412681e-03
x_pc   = 8.250000e+01
y_pc   = 2.503607e+02
M_p,xx = 3.932542e+08
M_p,yy = 1.610673e+08
x11_pc = 8.250000e+01
y22_pc = 2.503607e+02
M_p,11 = 3.932542e+08
M_p,22 = 1.610673e+08

Frame Analysis Example

The following example demonstrates how sectionproperties can be used to calculate the cross-section properties required for a frame analysis. Using this method is preferred over executing a geometric and warping analysis as only variables required for a frame analysis are computed. In this example the torsion constant of a rectangular section is calculated for a number of different mesh sizes and the accuracy of the result compared with the time taken to obtain the solution:

import time
import numpy as np
import matplotlib.pyplot as plt
import sectionproperties.pre.sections as sections
from sectionproperties.analysis.cross_section import CrossSection

# create a rectangular section
geometry = sections.RectangularSection(d=100, b=50)

# create a list of mesh sizes to analyse
mesh_sizes = [1.5, 2, 2.5, 3, 4, 5, 10, 15, 20, 25, 30, 40, 50, 75, 100]
j_calc = []  # list to store torsion constants
t_calc = []  # list to store computation times

# loop through mesh sizes
for mesh_size in mesh_sizes:
    mesh = geometry.create_mesh(mesh_sizes=[mesh_size])  # create mesh
    section = CrossSection(geometry, mesh)  # create a CrossSection object
    start_time = time.time()  # start timing
    # calculate the frame properties
    (_, _, _, _, j, _) = section.calculate_frame_properties()
    t = time.time() - start_time  # stop timing
    t_calc.append(t)  # save the time
    j_calc.append(j)  # save the torsion constant
    # print the result
    str = "Mesh Size: {0}; ".format(mesh_size)
    str += "Solution Time {0:.5f} s; ".format(t)
    str += "Torsion Constant: {0:.12e}".format(j)
    print(str)

correct_val = j_calc[0]  # assume the finest mesh gives the 'correct' value
j_np = np.array(j_calc)  # convert results to a numpy array
error_vals = (j_calc - correct_val) / j_calc * 100  # compute the error

# produce a plot of the accuracy of the torsion constant with computation time
plt.loglog(t_calc[1:], error_vals[1:], 'kx-')
plt.xlabel("Solver Time [s]")
plt.ylabel("Torsion Constant Error [%]")
plt.show()

Plot of the torsion constant as a function of the solution time.

Advanced Examples

The following examples demonstrates how sectionproperties can be used for more academic purposes.

Torsion Constant of a Rectangle

In this example, the aspect ratio of a rectangular section is varied whilst keeping a constant cross-sectional area and the torsion constant calculated. The variation of the torsion constant with the aspect ratio is then plotted:

import numpy as np
import matplotlib.pyplot as plt
import sectionproperties.pre.sections as sections
from sectionproperties.analysis.cross_section import CrossSection

# rectangle dimensions
d_list = []
b_list = np.linspace(0.2, 1, 20)
j_list = []  # list holding torsion constant results

# number of elements for each analysis
n = 500

# loop through all the widths
for b in b_list:
    # calculate d assuming area = 1
    d = 1 / b
    d_list.append(d)

    # compute mesh size
    ms = d * b / n

    # perform a warping analysis on rectangle
    geometry = sections.RectangularSection(d=d, b=b)
    mesh = geometry.create_mesh(mesh_sizes=[ms])
    section = CrossSection(geometry, mesh)
    section.calculate_geometric_properties()
    section.calculate_warping_properties()

    # get the torsion constant
    j = section.get_j()
    print("d/b = {0:.3f}; J = {1:.5e}".format(d/b, j))
    j_list.append(j)

# plot the torsion constant as a function of the aspect ratio
(fig, ax) = plt.subplots()
ax.plot(np.array(d_list) / b_list, j_list, 'kx-')
ax.set_xlabel("Aspect Ratio [d/b]")
ax.set_ylabel("Torsion Constant [J]")
ax.set_title("Rectangular Section Torsion Constant")
plt.show()

Plot of the torsion constant as a function of the aspect ratio.

Mesh Refinement

In this example the convergence of the torsion constant is investigated through an analysis of an I-section. The mesh is refined both by modifying the mesh size and by specifying the number of points making up the root radius. The figure below the example code shows that mesh refinement adjacent to the root radius is a far more efficient method in obtaining fast convergence when compared to reducing the mesh area size for the entire section:

import numpy as np
import matplotlib.pyplot as plt
import sectionproperties.pre.sections as sections
from sectionproperties.analysis.cross_section import CrossSection

# define mesh sizes
mesh_size_list = [50, 20, 10, 5, 3, 2, 1]
nr_list = [4, 8, 12, 16, 20, 24, 32, 64]

# initialise result lists
mesh_results = []
mesh_elements = []
nr_results = []
nr_elements = []

# calculate reference solution
geometry = sections.ISection(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=64)
mesh = geometry.create_mesh(mesh_sizes=[0.5])  # create mesh
section = CrossSection(geometry, mesh)  # create a CrossSection object
section.calculate_geometric_properties()
section.calculate_warping_properties()
j_reference = section.get_j()  # get the torsion constant

# run through mesh_sizes with n_r = 16
for mesh_size in mesh_size_list:
    geometry = sections.ISection(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=16)
    mesh = geometry.create_mesh(mesh_sizes=[mesh_size])  # create mesh
    section = CrossSection(geometry, mesh)  # create a CrossSection object
    section.calculate_geometric_properties()
    section.calculate_warping_properties()

    mesh_elements.append(len(section.elements))
    mesh_results.append(section.get_j())

# run through n_r with mesh_size = 3
for n_r in nr_list:
    geometry = sections.ISection(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=n_r)
    mesh = geometry.create_mesh(mesh_sizes=[3])  # create mesh
    section = CrossSection(geometry, mesh)  # create a CrossSection object
    section.calculate_geometric_properties()
    section.calculate_warping_properties()

    nr_elements.append(len(section.elements))
    nr_results.append(section.get_j())

# convert results to a numpy array
mesh_results = np.array(mesh_results)
nr_results = np.array(nr_results)

# compute the error
mesh_error_vals = (mesh_results - j_reference) / mesh_results * 100
nr_error_vals = (nr_results - j_reference) / nr_results * 100

# plot the results
(fig, ax) = plt.subplots()
ax.loglog(mesh_elements, mesh_error_vals, 'kx-', label='Mesh Size Refinement')
ax.loglog(nr_elements, nr_error_vals, 'rx-', label='Root Radius Refinement')
plt.xlabel("Number of Elements")
plt.ylabel("Torsion Constant Error [%]")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.tight_layout()
plt.show()

Plot of the torsion constant error as a function of number of elements used in the analysis for both general mesh refinement and root radius refinement.