Stress extraction on beam 188（Create finite element model directly with E command or Use the L command to create a line and then assign section properties and divide the mesh）

[bismarck1943]

I saw a code about calculating beam stress and deflection in pyansy official tutorial.
https://examples.mapdl.docs.pyansys.com/verif-manual/vm-002-beam_stresses_and_deflections.html#sphx-glr-verif-manual-vm-002-beam-stresses-and-deflections-py

I have calculated the stress at the mid-span of a simply supported beam with a height of 5m and a width of 3m under a uniformly distributed load of 10000.But there is a problem with the calculation result.

According to the stress calculation formula, σ=My/Iz, where M=1/8ql²=1/8 * 10000 * 32 * 32=1280000, Iz=bh³/12=35³/12=31.25, so σ =My/Iz=12800002.5/31.25=102400

But my code calculates to 112640，I think it may be because the number of element meshes created with the E command is not enough.

Then I tried to use the L command to create a line and then assign the section properties, and then perform grid division for calculation, but the calculation result is still problematic.

I wonder if the finite element model is created with the E command and the line is created with the L command, and the section properties are assigned and the grid is divided, the result is different.

Here is the Python code(Create finite element model directly with E command):

from ansys.mapdl.core import launch_mapdl

# start MAPDL and enter the pre-processing routine
mapdl = launch_mapdl()
mapdl.clear()
mapdl.prep7()
print(mapdl)

# Type of analysis: Static.
mapdl.antype("STATIC")

# Element type: BEAM188.
mapdl.et(1, "BEAM188")
mapdl.keyopt(1, 3, 3)  # Cubic shape function

mapdl.mp("EX", 1, 2e11)
mapdl.mp("PRXY", 1, 0.3)
print(mapdl.mplist())


# Parameterization of the cross-section dimensions.
sec_num = 1
b=3
h=5
# Define the beam cross-section.
mapdl.sectype(sec_num, "BEAM", "RECT")
mapdl.secdata(b,h)

# Print the section properties.
print(mapdl.slist())

mapdl.n(1, 0, 0, 0)
mapdl.n(2, 16, 0, 0)
mapdl.n(3, 32, 0, 0)

# Define one node for the orientation of the beam cross-section.
orient_node = mapdl.n(4, 8, 1)

# Print the list of the created nodes.
print(mapdl.nlist())

mapdl.e(1, 2, orient_node)
mapdl.e(2, 3, orient_node)
# Print the list of the created elements.
print(mapdl.elist())

# Display elements with their nodes numbers.
mapdl.eplot(show_node_numbering=True, line_width=5, cpos="xy", font_size=40)


# BC for the beams seats
mapdl.d(1, "UX", lab2="UY",lab3="UZ",lab4="ROTX",lab5="ROTY")
mapdl.d(3, "UY",ab2="UZ",lab3="ROTX",lab4="ROTY")

# Parametrization of the distributed load.
w = 10000

# Application of the surface load to the beam element.
mapdl.sfbeam(1, 1, "PRES", w)
mapdl.sfbeam(2, 1, "PRES", w)
mapdl.finish()

mapdl.run("/SOLU")
out = mapdl.solve()
mapdl.finish()

# Enter the post-processing routine and select the first load step.


# Get the maximum stress at the middle of the beam.
s_eqv_max = mapdl.get_value("secr", 2, "s", "eqv", "max")
print('MAX',s_eqv_max)

mapdl.exit()

The output is 112640 , just 1.1 * 102400
when canceling the command”mapdl.keyopt(1, 3, 3) ”，output is 51200，just 1/2 * 102400
Is this just a coincidence？

Then I try to use the L command to create a line and then assign section properties and divide the mesh.
I first performed calculations in ANSYS, and the calculated result was 102399.68, which was consistent with the theoretical value 1024 calculated by the previous formula

Here is the Ansys code(Use the L command to create a line and then assign section properties and divide the mesh):

/CLEAR
/PREP7
ET,1,BEAM188
MP,EX,1,2E11
MP,PRXY,1,0.3
SECTYPE,1,BEAM,RECT
SECDATA,3,5
SECPLOT,1,1
!SLIST,1,,,FULL
L=32
K,1
K,2,L/2
K,3,L
K,10,L/2,L/2
L,1,2
L,2,3
LATT,1,1,1,,,10,1
LESIZE,ALL,,,400
LMESH,ALL

/ESHAPE,1
EPLOT

DK,1, UX,,,,,UY,UZ,ROTX, ROTY
DK,3, UY,,,,,UZ, ROTX, ROTY
SFBEAM,ALL,1,PRES,10000
FINISH

/SOLU
SOLVE

/POST1
*get,MZ,SECR,,S,EQV,MAX

The output is 102399.68, which was consistent with the theoretical value 102400 calculated by the previous formula.

So I think there is nothing wrong with my Ansys code，but when I convert the ansys code into python code, the calculation result is wrong again.

Here is the Python code(Use the L command to create a line and then assign section properties and divide the mesh):

from ansys.mapdl.core import launch_mapdl

# start MAPDL and enter the pre-processing routine
mapdl = launch_mapdl()
mapdl.clear()
mapdl.prep7()
print(mapdl)


# Type of analysis: Static.
mapdl.antype("STATIC")
# Element type: BEAM188.
mapdl.et(1, "BEAM188")

mapdl.mp("EX", 1, 2e11)
mapdl.mp("PRXY", 1, 0.3)
print(mapdl.mplist())

# Parameterization of the cross-section dimensions.
sec_num = 1
b=3
h=5
# Define the beam cross-section.
mapdl.sectype(sec_num, "BEAM", "RECT")
mapdl.secdata(b,h)

# Print the section properties.
print(mapdl.slist())

mapdl.k(1, 0, 0, 0)
mapdl.k(2, 16, 0, 0)
mapdl.k(3, 32, 0, 0)
mapdl.k(100,8,1)

mapdl.l(1, 2)
mapdl.l(2, 3)
mapdl.latt(1, 1, 1, "",100, 1)
mapdl.lesize("ALL",'','',400)
mapdl.lmesh("ALL")

print(mapdl.elist())
# Display elements with their nodes numbers.
mapdl.eplot(show_node_numbering=True, line_width=5, cpos="xy", font_size=40)

# BC for the beams seats
mapdl.d(1, "UX", lab2="UY",lab3="UZ",lab4="ROTX",lab5="ROTY")
mapdl.d(3, "UY",lab2="UZ",lab3="ROTX",lab4="ROTY")

# Parametrization of the distributed load.
w = 10000

# Application of the surface load to the beam element.
mapdl.sfbeam('ALL', 1, "PRES", w)
mapdl.finish()

mapdl.run("/SOLU")
out = mapdl.solve()
mapdl.finish()

# Enter the post-processing routine and select the first load step.
mapdl.post1()

#mapdl.set(1)
# Get the maximum stress at the middle of the beam.
s_eqv_max = mapdl.get('mz',"secr",'', "s", "eqv", "max")
#s_eqv_max = mapdl.get_value("SECR",'', "S", "EQV", "MAX")
print('max',s_eqv_max)

mapdl.exit()

---

There is only "get" command in Ansys, but there are "get" and "get_value" in Pyansy

When I use the code as ansys code: (get)
“s_eqv_max = mapdl.get('mz',"secr",'', "s", "eqv", "max")”，the output is 0.00000000

When I use the code as pyansys tutorial code： (get_value)
“mapdl.set(1),
s_eqv_max = mapdl.get_value("SECR",'', "S", "EQV", "MAX")”，the output is 405513.219

I'm weird about the result calculated is wrong. Maybe there is still a problem in the process of converting the ansys code into python code.

But I finally want to know, how to get the correct beam stress in pyansy, use the E command or the L command, or both.

[mikerife]

Hi @bismarck1943
Distributed parallel solves does not store the results in the MAPDL database after a solve; so we need to read in the results from the result file with the SET command in order to post process. Depending on how you ran the MAPDL batch vs the default PyMAPDL MAPDL instance (defaults to distributed parallel) that could explain getting a 0 for some result.

For this kind of simply supported, uniform pressure loaded Timoshenko beam the displacement is a forth order polynomial along the length of the beam. Whereas the cubic option chosen (shape function) is third order. So it's not that the results are wrong; just inaccurate as we are trying to describe a 4th order polynomial with up to 3rd order term. Just need more elements along the length to get a more accurate answer. Doubling the number of beam elements to 4 total brings the s_eqv_max down to 104,960. Double again to 8 elements and the stress goes down to 103,040 (about 0.6% error). Double again and are down to about 0.1% error.

Mike
p.s. here is the PyMAPDL model I used:

import os
from ansys.mapdl.core import launch_mapdl

path = os.getcwd()
mapdl = launch_mapdl(run_location = path, additional_switches = '-smp')

sec_num = 1
b=3
h=5
w = 10000

mapdl.clear()
mapdl.prep7()

mapdl.et(1, "BEAM188")
mapdl.keyopt(1, 3, 3) # Cubic shape function

mapdl.mp("EX", 1, 2e11)
mapdl.mp("PRXY", 1, 0.3)

mapdl.sectype(sec_num, "BEAM", "RECT")
mapdl.secdata(b,h)

mapdl.k(1, 0, 0, 0)
mapdl.k(2, 32, 0, 0)
mapdl.k(100,8,1)

mapdl.l(1, 2)
mapdl.latt(1, 1, 1, "",100, 1)
mapdl.lesize("ALL",2)
mapdl.lmesh("ALL")

mapdl.finish()
mapdl.slashsolu()

mapdl.antype("STATIC")

node_1 = mapdl.queries.node(0,0,0)
node_2 = mapdl.queries.node(32,0,0)

mapdl.d(node_1, "UX", lab2="UY",lab3="UZ",lab4="ROTX",lab5="ROTY")
mapdl.d(node_2, "UY",lab2="UZ",lab3="ROTX",lab4="ROTY")

mapdl.sfbeam('all', 1, "PRES", w)

out = mapdl.solve()
mapdl.finish()

mapdl.post1()
mapdl.set('last')

s_eqv_max = mapdl.get_value("secr",'',"s","eqv","max")
print('MAX',s_eqv_max)

mapdl.exit()

[bismarck1943]

Thanks！(*￣︶￣)
After reading your code, I found that there is a problem when I define the constraint conditions. When using the L command, it should be the key point of the constraint, and I should use "mapdl.dk" instead of "mapdl.d".
Meshing still uses "mapdl.lesize("ALL",'','',100)"
“mapdl.set(1)”
”s_eqv_max = mapdl.get_value("secr",'',"s","eqv","max")”
the output value is 102401.0234375, which is correct.

But I still have a question. When I use the E command to build a unit like the Pyansys tutorial, will the calculation accuracy have a large deviation? Is there a problem with the tutorial or my understanding is wrong, because from the results, both 51200 and 112400 are just multiples of 102400, and there is no decimal point, so I am curious how it calculated .

https://examples.mapdl.docs.pyansys.com/verif-manual/vm-002-beam_stresses_and_deflections.html#sphx-glr-verif-manual-vm-002-beam-stresses-and-deflections-py

[mikerife]

Hi @bismarck1943 No - it really does not matter if you mesh a line or create the beam elements directly via the E command. Only the elements are sent to the solver. So as long as each method results in the same set of nodes/elements/beam orientation/loads/etc the two methods should return the same results.

Don't over think two numbers that just happen to be multiples of each other. Here it is the mesh density that is driving the accuracy. If we were to run cases of different element sized, for each of the three shape functions, we would get something like this where blue is linear, green is quadratic, and red is cubic:

[bismarck1943]

Thank you so much!!!(๑╹◡╹)ﾉ"""
You are really amazing, I understand what you mean now. That is, when I use the E command, the number of grid divisions depends on the number of elements I have created. When I use the L command, the number of grids depends on the parameters set by LATT. I have verified them separately, and the results of the two modeling methods are consistent.

As you calculated

elements eqv
2 112640
4 104960
8 103040
16 102560

[mikerife]

Hi @bismarck1943 You got it - one small correction is that the LESIZE command is used to size the mesh on a line, not LATT. Mike

[bismarck1943]

Oh yes.It is the LESIZE,my spelling mistake