Pilkey - Symmetric Channel Section#
This example re-creates the numerical example “B.2 Symmetric Channel Section” on page 437 of “Analysis and Design of Elastic Beams” by Walter D. Pilkey.
BibTeX reference:
@book{Pilkey,
author = {Pilkey, Walter D},
booktitle = {Analysis and Design of Elastic Beams},
edition = {First},
isbn = {0471381527},
language = {eng},
publisher = {Wiley},
title = {Analysis and Design of Elastic Beams: Computational Methods},
year = {2002},
}
Problem Description#
A symmetric open-channel section is analysed, with t=1
, h=18
and b=8
, see the figure below. The material properties of the section are taken to be E=210000000
and nu=0.33333
. Note that the elastic modulus plays no role in this analysis as the geometry is homogenous, however its value is included for completeness. Also note that the value for E
used by Pilkey is in kPa
(steel), whereas the problem is defined in inches
- this mixing of units is not an issue due to
the elastic modulus not affecting the geometric results.
Note that sectionproperties
uses an x-y
coordinate system rather than the y-z
system used by Pilkey.
[1]:
from IPython.display import Image
display(Image(filename="images/channel-geom.png"))
We can model the above geometry by generating a shapely
Polygon
from a list of points, then passing this Polygon
to the sectionproperties
Geometry
object.
[2]:
from shapely import Polygon
from sectionproperties.pre import Geometry, Material
t = 1 # channel thickness
h = 18 # distance between centreline of channel flanges
b = 8 # distance from edge of flange to centreline of web
# steel material
mat = Material(
name="Steel",
elastic_modulus=2.1e8,
poissons_ratio=0.33333,
yield_strength=1.0,
density=1.0,
color="lightgrey",
)
# define list of points, starting from lower left hand corner
points = [
(-0.5 * t, -0.5 * h - 0.5 * t),
(b, -0.5 * h - 0.5 * t),
(b, -0.5 * h + 0.5 * t),
(0.5 * t, -0.5 * h + 0.5 * t),
(0.5 * t, 0.5 * h - 0.5 * t),
(b, 0.5 * h - 0.5 * t),
(b, 0.5 * h + 0.5 * t),
(-0.5 * t, 0.5 * h + 0.5 * t),
]
# create shapely Polygon object
poly = Polygon(shell=points)
# create sectionproperties Geometry object
geom = Geometry(geom=poly, material=mat)
# plot geometry
geom.plot_geometry()
[2]:
<Axes: title={'center': 'Cross-Section Geometry'}>
Create mesh and Section
object#
The numerical analysis by Pilkey uses 9-noded quadraliteral elements. The mesh used by Pilkey for this problem can be seen below.
[3]:
display(Image(filename="images/channel-mesh.png"))
We can create a mesh in sectionproperties
using 6-noded triangular elements by defining a maximum triangular element area. In this case we choose mesh_sizes=0.1
and create the resulting Section
object.
[4]:
from sectionproperties.analysis import Section
geom.create_mesh(mesh_sizes=0.1)
sec = Section(geometry=geom)
sec.plot_mesh()
[4]:
<Axes: title={'center': 'Finite Element Mesh'}>
Calculate Cross-Section Properties#
Pilkey reports both geometric and warping properties, as such we conduct both analyses.
[5]:
sec.calculate_geometric_properties()
sec.calculate_warping_properties()
Comparison of Results#
The numerical results obtained by Pilkey is listed in the dictionary below.
[6]:
pilkey = {
"area": 34.0,
"qx": 0.0,
"qy": 63.75,
"cx": 1.875,
"cy": 0.0,
"x_sc": -2.86769,
"y_sc": 0.0,
"x_sct": -2.86759,
"y_sct": 0.0,
"ixx_g": 1787.83333, # N.B text erroneously printed 1781
"iyy_g": 342.83333,
"ixy_g": 0.0,
"ixx_c": 1787.83333,
"iyy_c": 223.30208,
"ixy_c": 0.0,
"zxx": 188.19298,
"zyy": 36.45748,
"rx": 7.25144,
"ry": 2.56275,
"phi": 0.0,
"alpha_x": 3.40789,
"alpha_y": 2.15337,
"alpha_xy": 0.0,
"j": 11.28862,
"gamma": 12763.15184,
}
Most of these results can be directly obtained from sectionproperties
, the only properties that require calculation are the shear coefficients, alpha
. The shear coefficient can be obtained from the shear area as follows:
\(\alpha = \frac{A}{A_s}\)
where \(A\) is the cross-section area and \(A_s\) is the shear area.
We create a similar dictionary for the sectionproperties
results.
[7]:
sectionproperties = {
"area": sec.get_area(),
"qx": sec.get_eq(e_ref=mat)[0],
"qy": sec.get_eq(e_ref=mat)[1],
"cx": sec.get_c()[0],
"cy": sec.get_c()[1],
"x_sc": sec.get_sc()[0],
"y_sc": sec.get_sc()[1],
"x_sct": sec.get_sc_t()[0],
"y_sct": sec.get_sc_t()[1],
"ixx_g": sec.get_eig(e_ref=mat)[0],
"iyy_g": sec.get_eig(e_ref=mat)[1],
"ixy_g": sec.get_eig(e_ref=mat)[2],
"ixx_c": sec.get_eic(e_ref=mat)[0],
"iyy_c": sec.get_eic(e_ref=mat)[1],
"ixy_c": sec.get_eic(e_ref=mat)[2],
"zxx": min(sec.get_ez(e_ref=mat)[:2]),
"zyy": min(sec.get_ez(e_ref=mat)[2:]),
"rx": sec.get_rc()[0],
"ry": sec.get_rc()[1],
"phi": sec.get_phi(),
"alpha_x": sec.get_area() / sec.get_eas(e_ref=mat)[0],
"alpha_y": sec.get_area() / sec.get_eas(e_ref=mat)[1],
"alpha_xy": sec.get_area() / sec.section_props.a_sxy,
"j": sec.get_ej(e_ref=mat),
"gamma": sec.get_egamma(e_ref=mat),
}
The comparison of results is summarised in the table below. Relative error is reported in all cases, except where a value is zero, in which the absolute error is reported.
[8]:
from rich.console import Console
from rich.table import Table
from rich.text import Text
# setup table
table = Table(title="Comparison of Results")
table.add_column("Property", justify="left", style="cyan", no_wrap=True)
table.add_column(Text("Pilkey", justify="center"), justify="right", style="green")
table.add_column(Text("sectionproperties", style="i"), justify="right", style="green")
table.add_column(Text("Error", justify="center"), justify="right", style="green")
# create a row for each property
for key in pilkey:
# get results
p_res = pilkey[key]
sp_res = sectionproperties[key]
# calculate relative error
if p_res != 0:
rel_error = (sp_res - p_res) / p_res
else:
rel_error = sp_res
# print row
table.add_row(key, f"{p_res:.4e}", f"{sp_res:.4e}", f"{rel_error:.2e}")
console = Console()
console.print(table)
Comparison of Results ┏━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Property ┃ Pilkey ┃ sectionproperties ┃ Error ┃ ┡━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━┩ │ area │ 3.4000e+01 │ 3.4000e+01 │ 8.36e-16 │ │ qx │ 0.0000e+00 │ -7.6776e-14 │ -7.68e-14 │ │ qy │ 6.3750e+01 │ 6.3750e+01 │ 8.92e-16 │ │ cx │ 1.8750e+00 │ 1.8750e+00 │ 3.55e-16 │ │ cy │ 0.0000e+00 │ -2.2581e-15 │ -2.26e-15 │ │ x_sc │ -2.8677e+00 │ -2.8677e+00 │ -7.52e-06 │ │ y_sc │ 0.0000e+00 │ 5.1688e-05 │ 5.17e-05 │ │ x_sct │ -2.8676e+00 │ -2.8676e+00 │ -6.97e-06 │ │ y_sct │ 0.0000e+00 │ 5.2920e-05 │ 5.29e-05 │ │ ixx_g │ 1.7878e+03 │ 1.7878e+03 │ 1.86e-09 │ │ iyy_g │ 3.4283e+02 │ 3.4283e+02 │ 9.72e-09 │ │ ixy_g │ 0.0000e+00 │ -6.9368e-13 │ -6.94e-13 │ │ ixx_c │ 1.7878e+03 │ 1.7878e+03 │ 1.86e-09 │ │ iyy_c │ 2.2330e+02 │ 2.2330e+02 │ 1.49e-08 │ │ ixy_c │ 0.0000e+00 │ -5.4973e-13 │ -5.50e-13 │ │ zxx │ 1.8819e+02 │ 1.8819e+02 │ 1.31e-08 │ │ zyy │ 3.6457e+01 │ 3.6457e+01 │ 8.21e-08 │ │ rx │ 7.2514e+00 │ 7.2514e+00 │ -4.63e-07 │ │ ry │ 2.5627e+00 │ 2.5628e+00 │ 1.58e-06 │ │ phi │ 0.0000e+00 │ 0.0000e+00 │ 0.00e+00 │ │ alpha_x │ 3.4079e+00 │ 3.4078e+00 │ -4.08e-05 │ │ alpha_y │ 2.1534e+00 │ 2.1533e+00 │ -1.63e-05 │ │ alpha_xy │ 0.0000e+00 │ 4.7965e-13 │ 4.80e-13 │ │ j │ 1.1289e+01 │ 1.1286e+01 │ -1.90e-04 │ │ gamma │ 1.2763e+04 │ 1.2762e+04 │ -7.16e-05 │ └──────────┴─────────────┴───────────────────┴───────────┘
All results are within acceptable limits. The torsion constant had the largest relative error, however this value is relatively small and acceptable given the differences in element type and mesh.
[9]:
err = (sectionproperties["j"] - pilkey["j"]) / pilkey["j"]
print(f"Torsion Constant Relative Error: {err:.6f}")
Torsion Constant Relative Error: -0.000190