First moment of area at given point

[LBSmulders]

I was wondering if it is possible to get the first moments of area (for example Sx(y)) at a given distance from the centroidal axis (or at a given point in the section). I would expect it to be in this package, but could not find it anywhere sadly.

Thanks for the help in advance :D

Kind regards,
Luuk

[robbievanleeuwen]

Hi Luuk, thanks for the question!

Could you elaborate on what you are trying to calculate? Do you want the first moments of the entire section taken about an arbitrary point? Or are you looking for the first moments of area of part of the section? Maybe also some context for how you use this value in design/engineering would help :)

Cheers,
Robbie

[LBSmulders]

Thanks for you quick reply!

I want to calculate the shear stress over the height of the cross-section due to a shear force using this formula:

There is a method in the package get_stress_at_points which works. However, I think due to the mesh mostly, there occurs a (for me) small random 'jump' in the calculated shear stresses over height (see pictures below). Since the Dutch design code I am using uses the previous mentioned formula, I would like to fix it using that and therefore I am wondering if there is an option to get the first moment of area at a given height (y) from the centroid.

The points inserted in the get_stress_at_points is a linspace of 100 points over the height with x=-0.2
vy = 1000
For the plot I divided the results by 1000 (KPa to MPa)

Luuk

[robbievanleeuwen]

Hi Luuk, first of all I would suggest refining your mesh - you should achieve much better results with an appropriately refined mesh. Also, see the example on retrieving stresses, where I do essentially the same thing that you are doing however achieve much better curves due to the improved mesh.

I would also make sure you are retrieving the stress in the direction of the shear force and not the resultant shear stress (see the above example) - I'm happy to review your code and help you arrive at a better solution if you'd like.

Lastly, if I'm not mistaken, the S(y) in the formula above is the first moment of area of the section that is only above y (i.e. delete the rest of the section below), or below y if you are below the centroid. This isn't built into sectionproperties but can be easily obtained with some geometry manipulation:

1. Move the geometry so it's centroid sits at the origin - this is so that your first moments of area are taken about the centroidal axis. This can be done with the align_center() method, setting align_to=None aligns the geometric centroid with the origin. This can also be done with CompoundGeometry objects if you have a composite section.
2. Use the split_section() method to split the geometry above and below the point you're interested in.
3. Calculate the geometric properties of the section above the split to obtain the first moment of area S(y)

Here's a short example that calculates S(y) for a rectangle (for a number of different values of y) and plots it (you can see the similarity with the expected shear stress diagram):

import matplotlib.pyplot as plt

from sectionproperties.analysis.section import Section
from sectionproperties.pre.geometry import Geometry
from sectionproperties.pre.library.primitive_sections import rectangular_section


def first_moment_modified(geom: Geometry, y: float) -> float:
    """Calculates the 'modified' first moment of area of the section above y."""
    top, bot = geom.split_section((0, y), vector=(1, 0))
    top[0].create_mesh(mesh_sizes=0)
    sec = Section(top[0])
    sec.calculate_geometric_properties()
    
    return sec.get_q()[0]

# create a 100 x 50 rectangle
geom = rectangular_section(100, 50).align_center()

# calculate q_y for 0, 1, 2, ..., 49
qx = [first_moment_modified(geom, y) for y in range(0, 49)]

# plot result
_, ax = plt.subplots()
ax.plot(qx, range(0, 49))
plt.show()

[robbievanleeuwen]

Here's a shot at replicating your problem:

import matplotlib.pyplot as plt
import numpy as np
from shapely import Polygon

from sectionproperties.analysis.section import Section
from sectionproperties.pre.geometry import Geometry


points = [
    (0, -1),
    (0, -0.4),
    (0.4, -0.2),
    (1.05, -0.2),
    (1.05, 0),
    (-1.5, 0),
    (-1.5, -0.2),
    (-0.75, -0.2),
    (-0.4, -0.4),
    (-0.4, -0.5),
    (-0.55, -0.8),
    (-1.5, -0.8),
    (-1.5, -1),
]
geom = Geometry(Polygon(points))
geom.create_mesh(mesh_sizes=0.003)
sec = Section(geom)
sec.calculate_geometric_properties()
sec.calculate_warping_properties()

x = [-0.2] * 100
y = np.linspace(-1, 0, 100)
pts = list(zip(x, y))

sigs = sec.get_stress_at_points(pts=pts, vy=1000)
tau_yz = [x[2] for x in sigs if x is not None]

fig, ax = plt.subplots()
ax.plot(tau_yz, y, "kx-")
ax.set_xlabel("y-coordinate")
ax.set_ylabel("Shear Stress")
plt.show()

stress = sec.calculate_stress(vy=1000)
stress.plot_stress(stress="vy_zy", cmap="viridis", normalize=False)

get_stress_at_points() uses linear interpolation within the elements (which are actually quadratic elements), so I suspect we're potentially seeing some numerical artefacts here for poor mesh qualities.

Note that the stresses displayed by plot_stress() have consistent interpolation so will always be correct:

[LBSmulders]

Thank you for your comprehensive explanation :D

I chose for the centering and splitting the cross-section, then a check if the point of interest is above or below the centroid to calculate the first moment of area from the top or bottom part, clever way of calculating it like this!

Great package, I use it a lot in my master thesis project, thank you and keep up the good work!

Luuk

[robbievanleeuwen]

Happy to help! Best of luck with masters thesis!

A bit of trivia - sectionproperties was created while I studied my masters @ TU Delft 😄