RBE3s causes an error in the plotting functions?

[MrMechanics]

Hi

I am using RBE3s to constrain and apply torque to one of my FE models, and I get some strange plots when I do so. First when I plot the mesh using eplot, I get the same plot as I would in ANSYS. But then after I create a single node with a mass element on it and create an RBE3 between that node and my mesh, the eplot function generates a plot that looks like it is now covered in a triangular mesh.

As far as I can tell the solver runs normally and I can't find that there have been created any new nodes or elements beyond what I myself created, but when I try to plot stresses and displacements of a section of the mesh, the plot is wrong. Specifically where the mesh is cut, the values are set to zero. And I need to see the stresses over the cross-section, so I am hoping maybe there is some small command or something I am using wrong or missing?

I made a small example script which reproduces the problem here:

#------------------------------------
# parameters
#------------------------------------
height = 20
length = 200
width = 5

mesh_size = 2.5

E_mod = 193000.
poisson = 0.27


#------------------------------------
# launch ansys mapdl
#------------------------------------
import time
from ansys.mapdl.core import launch_mapdl
import os
path = os.getcwd()

mapdl = launch_mapdl(run_location=path)


#------------------------------------
# reset mapdl
#------------------------------------
start_time = time.time()
mapdl.clear()
mapdl.prep7()
mapdl.view('1','-1','1','1')


#------------------------------------
# geometry
#------------------------------------
# create block
mapdl.block(0.,length,0.,height,0.,width)


#------------------------------------
# material properties
#------------------------------------
mapdl.mp('EX',102,E_mod)
mapdl.mp('NUXY',102,poisson)


#------------------------------------
# mesh
#------------------------------------
# set element size
mapdl.allsel()
mapdl.esize(mesh_size)

# create mesh
mapdl.allsel()
mapdl.et(102,'SOLID186')
mapdl.type(102)
mapdl.vsweep('ALL')

mapdl.allsel()
mapdl.emodif('ALL','MAT',102)

mapdl.allsel()
print('number of nodes and elements:')
print(len(mapdl.mesh.nnum), len(mapdl.mesh.enum))

# plot the mesh
mapdl.allsel()
mapdl.eplot(vtk=True, style="surface", show_edges=True, background='w',
            show_axes=True, line_width=2, plot_bc=True)

#------------------------------------
# node couplings
#------------------------------------
# make nodes for coupling used to apply
# loads and boundary conditions
mapdl.allsel()
mapdl.n(99999991,length+15,10,2.5)
mapdl.nsel('S','LOC','X',length+15)
mapdl.cm('torque_node','NODE')
mapdl.et(90,'MASS21')
mapdl.real(90)
mapdl.type(90)
mapdl.r(90,1.E-9,1.E-9,1.E-9)
mapdl.e(99999991)

# couple torque node to end of beam
mapdl.allsel()
mapdl.nsel('S','LOC','X',length)
mapdl.rbe3(99999991,'ALL','ALL')

mapdl.allsel()
print('number of nodes and elements (with RBE3s and mass elements):')
print(len(mapdl.mesh.nnum), len(mapdl.mesh.enum))


#------------------------------------
# boundary conditions
#------------------------------------
mapdl.allsel()
mapdl.nsel('S','LOC','X',0)
mapdl.d('ALL','ALL')


#------------------------------------
# loads
#------------------------------------
mapdl.allsel()
mapdl.cmsel('S','torque_node')
mapdl.f('ALL','FY',200)

# plot the mesh
mapdl.allsel()
mapdl.eplot(vtk=True, style="surface", show_edges=True, background='w',
            show_axes=True, line_width=2, plot_bc=True)


end_time = time.time()
print('pre-processing time: ', end_time - start_time, ' seconds')


#------------------------------------
# run solver
#------------------------------------
start_time = time.time()
mapdl.allsel()
mapdl.slashsolu()
mapdl.antype('STATIC')
mapdl.nlgeom('ON')

mapdl.autots('ON')
mapdl.nsubst(10,10,1,'OFF')
mapdl.time(1.)
print('Starting solver...')
mapdl.solve()

end_time = time.time()
print('analysis time: ', end_time - start_time, ' seconds')


#------------------------------------
# post processing
#------------------------------------
# enter post-processor
mapdl.post1()
mapdl.set('LAST')

mapdl.allsel()
mapdl.nsel('S','LOC','X',0,length-100)
mapdl.esln('S')

mapdl.upcoord(1)
# plot VonMises stress (average nodal)
mapdl.post_processing.plot_nodal_eqv_stress(show_edges=True)
# plot displacement
mapdl.post_processing.plot_nodal_displacement(show_edges=True)
mapdl.upcoord(-1)

mapdl.exit()

[mikerife]

Hi @MrMechanics
I think this should be filed as an Issue - I will do that later today. In the meantime unselect the point mass element prior to the element or result plots. Or post process in MAPDL maybe using the mapdl.open_gui() method.

For example when plotting the mesh:

mapdl.eslv()
mapdl.nsle()
mapdl.eplot(vtk=True, style="surface", show_edges=True, background='w',
            show_axes=True, line_width=2, plot_bc=True)

Mike

[MrMechanics]

Hi @mikerife
Unselecting the mass elements did indeed remove the triangular mesh coating the elements, but the plotted values are still wrong for nodal stress and displacement plots. Plotting element stresses is working though, and it looks better without the triangular mesh coating, so I have a workaround with lower resolution there.

Thanks for your help,
Mats

[mikerife]

Hi @MrMechanics Mats
Did you unselect the mass element before plotting results too? Mass21 only stores degree-of-freedom results (as we would expect) and I think this is getting translated as inappropriate zeros getting that to the (VTK) result plots.
Mike

[MrMechanics]

Hi @mikerife

I should be more specific. If you unselect the mass elements, the triangular mesh coating goes away, and the nodal plot values are correct. But if you also unselect other elements because you want to look at a cross-section of your model, the nodal values plotted at that cross-section are set to zero. This is not the case if you plot element values, but plotting element stresses for example plots one average (or minimum or maximum) element stress over the entire element, so the resolution is lower than plotting nodal stresses would be.

Would be great if I could plot element nodal stresses like it is plotted in ANSYS, but I haven't figured out if that is even possible in pymapdl at this time.

Best regards,
Mats

[mikerife]

Hi @MrMechanics Mats
Be sure to select the set of nodes of the selected set of elements. This should give you the expected plot:

mapdl.post1()
mapdl.set('LAST')

mapdl.allsel()
mapdl.nsel('S','LOC','X',0,length-100)
mapdl.esln('S')
mapdl.esel('u','enam','',21)
mapdl.nsle()
mapdl.upcoord(1)
mapdl.post_processing.plot_nodal_eqv_stress(show_edges=True)

Mike

[MrMechanics]

Hi @mikerife
You are correct. That did work.

Thank you again!
best regards,
Mats