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