How to get the stress matrix of elements in global coordinates, and question on compas_fea.functions.principal_stresses()

[ioannaMitropoulou]

Hello!

I am using compas_fea - abaqus to do analysis of shells, and I am having a bit of trouble making sense of the results that I get. I describe here the process that I follow - I try to keep it short-, in case you can spot where I’ve been confused.

I want to get the stress matrix and from that compute the stress for arbitrary orientations on the shell. To do that, for each element I take the structure.results[step][‘element’][‘sxx’], [‘sxy’], [‘syy’] for the ip1_sp1 output and create matrix A, which in my understanding should be the stress matrix in local coordinates for that element's ip1-sp1.

A = Sxx  Sxy  0,0
    Sxy  Syy  0.0
    0.0  0.0  0.0

I also get the local axis vectors ex, ey, ex from the structure.results[step][‘element’][‘axes’] dict and convert A to global coordinates using the compas_functions

frame_local = Frame(Point(0,0,0), ex, ey)
A_world = np.array(local_to_world_coordinates_numpy(frame_local, A.T)).T

To verify if the result is correct, I find the eigenvectors of A_world, and compare them to the results of compas_fea.function.principal_stresses,

vec1, vec5, pr1, pr5, pmax   = compas_fea.functions.principal_stresses(data, ptype='max', scale=1, rotate=0)

In my understanding, vec1should be the same as the eigenvector that corresponds to the largest eigenvalue of A_world, but this is not the case, the results I get are completely different. I must have misunderstood the meaning of the entries in structure.results dict, but I can't figure out where my mistake is. I would be grateful if you could point me in the right direction.

Also, I have a question regarding the compas_fea.functions.principal_stresses(). I was wondering why you are only using the values of ip1, instead of also considering the other integration points, or some average(?) value corresponding to the center of each side of the element.

Many thanks! Ioanna

[franaudo]

Ciao Ioanna and thanks for your feedback!

Probably the best person to answer you is @andrewliew. In the meantime, maybe could you share a simple script replicating the problem? just a simple structure with a few shell elements showing the different results.

I can also have a look at it in the evening :)

best, Francesco

[ioannaMitropoulou]

Hello, thanks for your reply! I attach here a simple example using the file from the compas-fea example: https://compas.dev/compas_fea/latest/examples/mesh_principal_rhino.html

I have added a gh file that visualizes the results for easier comparison. compas_fea_question.zip

Thanks! Ioanna

[andrewliew]

Hi Ioanna,

You are correct, the stresses that are returned from Abaqus, sxx, syy and sxy for shells, are in the local coordinates ex, ey of the element.

Take note (like you were asking) on the number of ip-sp points in your element, if it is reduced integration for a triangular element (which is all I got around to coding up at the time) you will only have one data point (Gauss/integration point). I had to always convert a quad mesh into a triangular one for the function to work. This is because for full integration and/or quads there will be some more data that you might want to take, average or hack and change the function to exploit. Note also that there is data for section point 5 and 1, which can sometimes be quite different.

We need to step back to understand if there is a relationship between the eigenvectors of A and the principal stresses. Principal stresses are calculated via Mohr circles (as is done in the principal stresses function) and is basically a tensor manipulation process, I am not personally aware of a mathematical relationship between these stress Mohr circles and the eigenvectors of A. That is not to say that there isn't one, but check you really are expecting such a relationship.

The principal stresses function was quite experimental when I was working with it, so use it with some caution, you could check some results with Karamba, as that has some principal stress output ability. For the few times I used it it looked right, but sometimes I had to use the rotate argument to rotate the results by 90 degrees for it to be aligned properly (maybe was a local axis thing that I never got to the bottom of), might want to check that too.

[franaudo]

this might help: https://uclageo.com/SoilMechanicsNotes/Section2.3.php

PS: I will have also time to check what you sent Ioanna on Friday :)

[andrewliew]

Nice find^, good to know there is a clean relationship between the two. If you still have problems especially if the results are out by 90 degrees, try the rotate argument. Also check if your loading is in-plane (where sp1 and sp5 data should be almost identical) and if the shell has out-of-plane loading (sp1 and sp5 can differ).

[ioannaMitropoulou]

Hello, thanks for the clarifications!

First, I was using quad elements, which gave me 4 integration points, thus my confusion with them. Now I changed my meshes to triangular and I only get one integration point, so it makes sense. I would suggest if you could also change the example that calculates principal stresses (https://compas.dev/compas_fea/latest/examples/mesh_principal_rhino.html), because that is applied on a quad mesh.

I am now using Karamba to check the results. The results seem to agree (more or less) with the results of compas_fea.functions.principal_stresses (up to a PI/2 or PI rotation), and not with my eigenvectors. However, what confuses me is that the vec-min and vec-max for every sp from compas-fea are parallel. In the following screenshots the arrows are from Karamba, and the lines are from compas-fea (vec max and min on sp1).

on a 2D shell (vec-min and vec-max are identical, so they appear as one line): and on a 3D shell (here vec-min and vec-max are parallel):

In contrast, the results don't match at all with the eigenvectors: So I must have made a mistake there, I'm trying to figure out where.

Regarding the compas-fea vec min and max, in my understanding, they should have an angle of 90 degrees, or is that not the case?

Many thanks! Ioanna

[franaudo]

Hi Ioanna,

This is a sample code for the extraction of the principal stresses and directions from the compas_fea json results file. You can skip the first part (which is just the creation of a simple structure so we don't have to deal with Rhino - which I don't have on this machine).

As you can see I am using quads (which have 4 integration points), so you can get the principal stresses at each integration point or do an average if you want.

"""
Principal stressess extraction from compas_fea results file.
Author(s): Francesco Ranaudo (github.com/franaudo)

This example shows how to extract the principal stresses at an integration point
of an element from the compas_fea results file.
In the first part, a simple planar structure is generated using shell elements.
Once the analysis is complete, the stress components at the integration point of
an element are used to calculate the principal stresses and directions at that
location.
"""

from compas_fea.structure import ElasticIsotropic
from compas_fea.structure import ElementProperties as Properties
from compas_fea.structure import GeneralStep
from compas_fea.structure import GravityLoad
from compas_fea.structure import PinnedDisplacement
from compas_fea.structure import ShellSection
from compas_fea.structure import Structure

import numpy as np
import json
from math import atan2, degrees
from pathlib import Path

## ------------------------------ PART 1 -------------------------------------##
# Create empty Structure object (nxn plate)
folder = 'C:/temp/principal_stresses/'
name = 'principal_stresses'

mdl = Structure(name=name, path=folder)

minX, maxX, minY, maxY = 0., 1., 0., 1.
discretization = 5
x = np.linspace(minX, maxX, int((maxX-minX)*discretization+1))
y = np.linspace(minY, maxY,  int((maxY-minY)*discretization+1))
X, Y = np.meshgrid(x, y)
X = X.reshape((np.prod(X.shape),))
Y = Y.reshape((np.prod(Y.shape),))
Z = np.zeros_like(X)
mdl.add_nodes([node for node in zip(X, Y, Z)])

for j in range(int(maxY*discretization)):
    for i in range(int(maxX*discretization)):
        mdl.add_element(nodes=[int(i+((maxX*discretization+1)*j)), int(i+((maxX*discretization+1)*j)+1), int(i+((maxX*discretization+1)*j)+maxX*discretization+2),
                        int(i+((maxX*discretization+1)*j)+maxX*discretization+1)], type='ShellElement')
mdl.add_set(name='nset_base', type='node', selection=[0, int(
    maxX*discretization), int((maxX*discretization+1)*(maxY*discretization)), int((maxX*discretization+1)*(maxY*discretization+1)-1)])
mdl.add_set(name='elset_shell', type='element', selection=[i for i in mdl.elements.keys()])
mdl.add([ShellSection(name='sec_shell', t=0.05), ])
mdl.add(ElasticIsotropic(name='mat_elastic', E=10*10**9, v=0.3, p=1500))
mdl.add([Properties(name='ep_shell', material='mat_elastic', section='sec_shell', elset='elset_shell'), ])
mdl.add([GravityLoad(name='load_gravity', elements='elset_shell'), ])
mdl.add(PinnedDisplacement(name='disp_pinned', nodes='nset_base'))
mdl.add([
    GeneralStep(name='step_bc', displacements=['disp_pinned']),
    GeneralStep(name='step_loads', loads=['load_gravity']),
])
mdl.steps_order = ['step_bc', 'step_loads']
mdl.summary()
mdl.analyse_and_extract(software='abaqus', fields=['u', 's'])

## ------------------------------ PART 2 -------------------------------------##
# Read the results file
with open(Path(folder).joinpath(name, f"{name}-results.json"), "r") as f:
    results = json.load(f)

# Get the stress values for an element at an integration value
# Note: only 1 section point is considered
element_id = '6'
integration_point = 'ip1'
stress_vector = [results['step_loads']['element'][s][element_id]
                 [f'{integration_point}_sp1'] for s in ['sxx', 'syy', 'sxy']]
stress_matrix = np.array([(stress_vector[0], stress_vector[2]), (stress_vector[2], stress_vector[1])])

# The principal stresses and their directions are computed solving the eigenvalues problem
pricnipal_stresses, principal_directions = np.linalg.eig(stress_matrix)

print('principal stresses: ', pricnipal_stresses)
print('principal axes (basis): ', principal_directions)
print('inclination of the principal axes (w.r.t. World): ', degrees(
    atan2(principal_directions[0, 1], principal_directions[0, 0])))

Below the Max and Min (in-plane) principal stresses directly from Abaqus where you can see that the results match:

I will ASAP create a pull request with this example file and a modification of the compas_fea function using this method (which can be easily extended to the 3d case).

I will close this issue now, but please feel free to reopen it if you have more questions.

Ciao, Francesco

[ioannaMitropoulou]

Great, thanks a lot for the clarifications and the example! Best, Ioanna

[ioannaMitropoulou]

Hello guys, It seems that the changes made by #118 to the file utilities/functions.py in answer to this issue, were afterwards reverted in commit 59b136f88bca5d9a3d8cb2e3493dad3319588175 by @brgcode . Was that intentional?

[tomvanmele]

no this happened by accident. will fix (or if you send a PR, gladly approve it :)

[ioannaMitropoulou]

Hello and thanks for the quick response! I'm looking into it, but unfortunately, the file was replaced completely twice, so I'm not 100% sure what has been changed in the individual commits. I'm not familiar enough with the code to know what to keep.

[franaudo]

no worries, I will have a look tonight :)