Obtaining stress at specific cross-section location/point

[scote89]

Im interested in getting the cross-section stress at a specific location/point but it does not look like there is a direct function for it right? Any idea how one could achieve this?

Are the mesh points (mesh.points) and the stresses (stresses[0]["sig_zz"]) in the same order? If so, I would imagine having to write an interpolation function between the closest mesh point.

[Spectre5]

This is also something I'd like to see added. I think you could use the few nearest nodes to interpolate it, but a more appropriate way is to find the element the point lies within and then use it's shape functions to calculate the stress from the elements Gauss points.

[scote89]

Here's how I made it work on my end:

from scipy.interpolate import interp2d
# Get stress at specific location
f = interp2d(x=[i[0] for i in mesh.points], y=[i[1] for i in mesh.points], z=stresses[0]["sig_zz"], kind = 'linear', copy=True, bounds_error=True, fill_value=None)
f(-24.0,9.27) # Get stress at point (x,y)

[Czarified]

This is something I ran into in testing. Yes, the stresses are in the same order as the mesh point. If you want stress at a specific node, and you know the coordinates of that node, I wrote a simple function that will return the proper index. This works for getting stress peaks at corners, for example, like my Peery tests.

As Spectre5 pointed out, the most accurate way to get the stress at any arbitrary point would require using the shape function, since the elements sectionproperties uses are not constant-strain triangles. I'd have to dig up my FE book again, but I believe the CTRIA6 have a linear function, so you should may be pretty close the the "actual" value with your linear interpolation.

def get_node(nodes, coord) -> (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

Then, index, _ = get_node(nodes, coord) where coord is a tuple of the (x, y) coordinates you want. Then you can pass index to the get_stress method. See test_Peery.py for example usage (line 134, test_fb_a).

Edit: I was wrong about the CTRIA6 shape function. Updated my comment above.

[robbievanleeuwen]

Hi all,

Just following up on this thread a few months late with a few comments:

For reference, the below is detailed in this document, I have included references for ease of finding the correct section.

• The six-noded triangular elements have quadratic shape functions (3.2.1)
• As per typical FEM practice, the stresses are initially calculated at the six gauss points within the triangular element. For reporting and plotting purposes, the stresses are then extrapolated to the triangular nodes and are averaged with results from adjacent elements (3.4). Section 5 covers the specifics of the stress calculations.
• For an arbitrary point, the 'correct' way would be to first identify which element the point lies within. You would then use quadratic shape function interpolation from the gauss points results (not the nodal results) to calculate the stress at the arbitrary point.

This is definitely something that could be implemented in the future if it's deemed useful? At this stage, it would probably be a post v2 implementation, but let me know your thoughts!

Robbie

[ccaprani]

Just want to note here that the feature added in PR#110 should refer to the API introduced under this feature, when implemented.

[Czarified]

@scote89, Colin's feature will now plot a Mohr's Circle at any point. If you can't work with that, you can place a point in the geometry there, which will place a node in the mesh, and you can pull the stresses with my function above.

If you want the actual code, look at Colin's function source. You basically pull all the data for the stress field, construct an interpolator (he uses a linear interpolator from matplotlib), and solve for whatever combined state you want. The details are frankly beyond me without breaking out my textbooks. Hopefully just having the circles will give you what you want?

[scote89]

Thanks, I tried the following code but I always get a zero value array while the mohr plot gives me some values. I also obtain a 0 value array when I do stress_post.section.material_groups[0].stress_result.sig_zz

def get_stress_at(stress_post,x,y):
    pt = (x,y)
    sect=stress_post.section
    nodes = sect.mesh_nodes
    ele = sect.mesh_elements       
    triang = tri.Triangulation(nodes[:, 0], nodes[:, 1], ele[:, 0:3])

    # Find in which material group the point lies
    pt_group = None
    for group in stress_post.section.material_groups:            
        mask_array = np.ones(len(stress_post.section.elements), dtype=bool)
        mask_array[group.el_ids] = False
        triang.set_mask(mask_array)
        trifinder = triang.get_trifinder()
        if trifinder(*pt) != -1:
            pt_group = group
        triang.set_mask(None)
    
    if pt_group is None:
        raise Exception(f"Point {(*pt,)} is not within mesh")

    # Assesmble the stress results from the relevant material group
    sigma_zz_v = pt_group.stress_result.sig_zz
    tau_xz_v = pt_group.stress_result.sig_zx
    tau_yz_v = pt_group.stress_result.sig_zy

    # Get the interpolators
    sigma_zz_interp = tri.LinearTriInterpolator(triang, sigma_zz_v)
    tau_xz_interp = tri.LinearTriInterpolator(triang, tau_xz_v)
    tau_yz_interp = tri.LinearTriInterpolator(triang, tau_yz_v)

    # Get the stresses at the point
    sigma_zz = sigma_zz_interp(*pt).item()
    tau_xz = tau_xz_interp(*pt).item()
    tau_yz = tau_yz_interp(*pt).item()
    
    print(sigma_zz_v)
    
    return {"sigma_zz":sigma_zz,"tau_xz":tau_xz,"tau_yz":tau_yz}
get_stress_at(stress_post,200,200)['sigma_zz']

[robbievanleeuwen]

Hi @scote89,

There is a small bug in the above code which just fetches the initialised stress results (all zero) instead of the computed stress. Line 10 should read:

for group in stress_post.material_groups:

instead of

for group in stress_post.section.material_groups:

See working gist for the above here.