Symmetric and Unsymmetric Beams in Complex Bending

Calculate section properties of two different beams given in examples from ‘Aircraft Structures,’ by Peery. These cases have known results, and the output from sectionproperties can be compared for accuracy. These examples represent a more rigourous ‘proof’ against a ‘real’ problem. Only results that have values in the reference material are tested here.

BibTeX Entry for reference:

@Book{Peery,
    title = {Aircraft Structures},
    author = {David J. Peery},
    organization = {Pensylvania State University},
    publisher = {McGraw-Hill Book Company},
    year = {1950},
    edition = {First},
    ISBN = {978-0486485805}
}

from sectionproperties.pre.library import nastran_sections
from sectionproperties.analysis.section import Section

Example 1 in Sec. 6.2 (Symmetric Bending)

This is a symmetric I-section with no lateral supports, undergoing pure unidirectional cantilever bending. Note that units here are inches, to match the text.

We’ll use a very coarse mesh here, to show a conservative comparison for accuracy. Theoretically, with more discretization, we would capture the real results more accurately.

geometry = nastran_sections.nastran_i(6, 3, 3, 1, 1, 1)
geometry = geometry.shift_section(x_offset=0, y_offset=-3)
geometry = geometry.create_mesh(mesh_sizes=[0.25])
section = Section(geometry)
section.plot_mesh()

<AxesSubplot: title={'center': 'Finite Element Mesh'}>

Perform a geometric analysis on the section, and plot properties We don’t need warping analysis for these simple checks, but sectionproperties needs them before evaluating stress.

section.calculate_geometric_properties()
section.calculate_warping_properties()
section.plot_centroids()

<AxesSubplot: title={'center': 'Centroids'}>

Directly from the example, we know that the 2nd moment of inertia resisting the bending is 43.3 in⁴.

section.section_props.ixx_g

43.3333333333333

From statics, we know the max bending moment on the beam will be 80,000 in-lbs. We can apply this moment to the section, and evaluate stress.

moment = 8e5
stress = section.calculate_stress(Mxx=moment)

Next we can extract the max stress from the section, and let’s go ahead and look at the calculated fringe plot. Refer to the stress example for details.

numerical_result = max(stress.get_stress()[0]["sig_zz"])
stress.plot_stress_zz()

<AxesSubplot: title={'center': 'Stress Contour Plot - $\\sigma_{zz}$'}>

From the book, and simple statics, we know the max stress is 55,427.3 psi.

numerical_result

55384.61538461558

This example is admittedly more simple, but it’s still a nice check for the basics on validity of the package.

print(f"Theoretical Result = 55427 [psi]")
print(f"  Numerical Result = {numerical_result:.0f} [psi]")
acc = (55427 - numerical_result) / 55427
print(f"          Accuracy = {acc:%}")

Theoretical Result = 55427 [psi]
  Numerical Result = 55385 [psi]
          Accuracy = 0.076469%

Example 1 in Sec. 7.2. (Unsymmetric Bending)

Moving on to something a bit more advanced… This is an unsymmetric Z-section with no lateral supports. Note that units here are inches, to match the text.

base_geom = nastran_sections.nastran_zed(4, 2, 8, 12)
base_geom = base_geom.shift_section(-5, -6)
base_geom = base_geom.create_mesh([0.25])
section = Section(base_geom)
section.calculate_geometric_properties()
section.calculate_warping_properties()
section.plot_centroids()

<AxesSubplot: title={'center': 'Centroids'}>

Checking each property against the reference text:

props = section.section_props
print("    Property | Theoretical | Numerical")
print(f"    ixx_g    | {693.3:<12.1f}| {props.ixx_g:<.1f}")
print(f"    iyy_g    | {173.3:<12.1f}| {props.iyy_g:<.1f}")
print(f"    ixy_g    | {-240:<12.1f}| {props.ixy_g:<.1f}")
print(f"    i11_c    | {787:<12.1f}| {props.i11_c:<.1f}")
print(f"    i22_c    | {79.5:<12.1f}| {props.i22_c:<.1f}")

Property | Theoretical | Numerical
ixx_g    | 693.3       | 693.3
iyy_g    | 173.3       | 173.3
ixy_g    | -240.0      | -240.0
i11_c    | 787.0       | 787.2
i22_c    | 79.5        | 79.5

Section properties all look good, so we can move on to some stress analysis. Before we do, we will need a quick function to pull stress at a certain point. This is a bit of a hack! Future versions of sectionproperties will have a much more robust system for getting stress at an arbitrary location. This particular function will work for the locations we need, since we know a node will be there.

from typing import Tuple


def get_node(nodes, coord) -> Tuple[int, tuple]:
    """
    This function will loop over the node list provided,
    finding the index of the coordinates you want.
    Returns the index in the nodes list, and the coords.
    """
    for index, var in enumerate(nodes):
        if all(var == coord):
            return index, var
        else:
            continue

    raise ValueError(f"No node found with coordinates: {coord}")

The load applied in the reference is -100,000 in-lbs about the x-axis, and 10,000 inb-lbs about the y-axis.

stress = section.calculate_stress(Mxx=-1e5, Myy=1e4)

Check stress at location A (see docs page for details)

A = (-5, 4)
text_result = 1210
n, _ = get_node(section.mesh_nodes, A)
numerical_result = stress.get_stress()[0]["sig_zz"][n]
print(text_result, numerical_result)

1210 1210.2272727272666

Check stress at location B (see docs page for details)

B = (-5, 6)
text_result = 580
n, _ = get_node(section.mesh_nodes, B)
numerical_result = stress.get_stress()[0]["sig_zz"][n]
print(text_result, numerical_result)

580 579.545454545451

Check stress at location C (see docs page for details)

C = (1, 6)
text_result = -2384
n, _ = get_node(section.mesh_nodes, C)
numerical_result = stress.get_stress()[0]["sig_zz"][n]
print(text_result, numerical_result)

-2384 -2386.363636363626

Looking at total axial stress over the section.

stress.plot_stress_zz()

<AxesSubplot: title={'center': 'Stress Contour Plot - $\\sigma_{zz}$'}>