No contact penetration of CONTA178 elements - Bearing simulation #2123

herttim · 2023-06-16T06:41:37Z

herttim
Jun 16, 2023

Hello all,

I would like to calculate the contact penetration of a tapered roller bearing by using CONTA178 elements between the raceways. So that, I´ve created a model with 30 contact elements including a pre-defined gap along the surface of the raceway and applied a constant pressure around the surface of the outer ring to simulate the pre-load of the bearing. The gap should simulate the profile of the logarithmic raceways.
However, while solving the model I didn´t get any connection between the outer ring and rollers. The CONTA178 elements are still open and there isn´t any contact penetration. I´m not sure whats the reason for that, I tired to vary the boundary conditions a bit, but that doesn´t helped.

I would be very happy if someone here can help me.

Thank you,
herrtim

Answered by pmaroneh

Jun 16, 2023

Hi @herttim, your question is related to the modelling approach and usage of MAPDL, and not on a specific coding issue with PyMAPDL. For such type of questions, please contact your local Ansys support provider, or use the Ansys Innovation Space forum. Thank you for your understanding.

View full answer

pmaroneh · 2023-06-16T07:38:47Z

pmaroneh
Jun 16, 2023
Maintainer

Hi @herttim, your question is related to the modelling approach and usage of MAPDL, and not on a specific coding issue with PyMAPDL. For such type of questions, please contact your local Ansys support provider, or use the Ansys Innovation Space forum. Thank you for your understanding.

3 replies

herttim Jun 16, 2023
Author

Hi @pmaroneh, your answer is so helpful, you could have saved this one....
It is a Q&A for PyMAPDL, I am coding in PyMAPDL...almost every question in the Q&A could be posted in the Ansys forum or Python forum as well. So I suggest you to copy your helpful comment and paste it under every single discussion.

I just wanted to find someone who faced the same issue and discuss possible solutions and ideas in PyMAPDL.

Thank you,
herttim

pmaroneh Jun 16, 2023
Maintainer

@herttim Just wanted to let you know what guidance we use from a company perspective. But indeed some people from the community might get back to you.

ChrisHarrold Jun 16, 2023

@herttim, Pernelle's answer to your question is the only answer anyone on this repository would or should provide to this question, as it does not relate at all to the open source project PyMAPDL, and is you looking for someone to try and solve your specific problem. This is not the place to post questions about your specific model that isn't working, that's what AIS and tech support is for. Even the developer forum would be a much more reasonable place for this question.

If there is a specific feature of PyMAPDL that you want to post an issue or question about this would be the correct venue for that.

mikerife · 2023-06-19T16:29:18Z

mikerife
Jun 19, 2023
Maintainer

Hi @herttim
Could you please post the PyMAPDL script (or zip it up and post the zip file) and I will take a look at the model.

@ChrisHarrold @pmaroneh @germa89 with all due respect to you (my coworkers), herttim has asked similar questions in the past and the answers have varied between being a pure Python issue, MAPDL command issue, and using PyMAPDL. We really need to look at the code, and probably run it, before assuming the answer is purely a MAPDL answer.

Mike

1 reply

herttim Jun 20, 2023
Author

@mikerife Thank you for your support!

Please find my PyMAPDL code below:
20220711_TRB_Model_V3.zip

import matplotlib.pyplot as plt
import numpy as np
import math as m
import itertools
import pandas

from ansys.mapdl.core import launch_mapdl

mapdl = launch_mapdl()

mapdl.prep7()

mapdl.units("SI")

mapdl.et(1, "SOLID185")
mapdl.mp("EX", 1, 210e9)
mapdl.mp("DENS", 1, 7800)
mapdl.mp("NUXY", 1, 0.3)
mapdl.type(1)
mapdl.keyopt(1, 2, 3)

mapdl.et(2, "CONTA178")
mapdl.mp("EX", 2, 210e9)    # Elastic moduli in Pa (kg/(m*s**2))
mapdl.mp("DENS", 2, 7800)   # Density in kg/m3
mapdl.mp("NUXY", 2, 0.3)    # Poisson's Ratio
mapdl.type(2)
mapdl.keyopt(2, 3, 1)
mapdl.keyopt(2, 4, 1)
mapdl.keyopt(2, 5, 0)
mapdl.keyopt(2, 10, 0)
mapdl.keyopt(2, 12, 1)

# Slicing
anzahl_scheiben = 30
anzahl_scheiben_board = 20

# Geometry
anzahl_rollen = 16

# Inner ring

di = 0.065      #Innen-Durchmesser Lagersitz
bi = 0.041      #Breite Lagersitz
α  = 7.6        #Druckwinkel
β  = 9.6        #Hilfswinkel

o = 0.009238
p = 0.008282
q = 0.003962
r = 0.000671
s = 0.02539
t = 0.003391
v = 0.005083
w = 0.00079 #0.00086
b = 0.006658

# Inner ring - Keypoints

mapdl.k(1, 0, di/2, 0)
mapdl.k(2, 0, di/2+o, 0)
mapdl.k(3, 0, di/2+o, p)
mapdl.k(4, 0, di/2+o-q, p+r)
mapdl.k(5, 0, di/2+o-q+t, p+r+s)
mapdl.k(6, 0, di/2+o-q+t+v, p+r+s-w)
mapdl.k(7, 0, di/2+o-q+t+v, p+r+s+b)
mapdl.k(8, 0, di/2, bi)

mapdl.k(13, 0, 0, 0)
mapdl.k(14, 0, 0, -0.005)
mapdl.k(15, 0, 0, 0.005)

# Inner ring - Lines

li_1 = mapdl.l(1, 2)
mapdl.l(2, 3)
mapdl.l(3, 4)
l1_i = mapdl.l(4, 5)
l1_b = mapdl.l(5, 6)
mapdl.l(6, 7)
mapdl.l(7, 8)
li_8 = mapdl.l(8, 1)

# Inner ring - Surface

a1 = mapdl.a(1, 2, 3, 4, 5, 6, 7, 8)

mapdl.lesize(nl1=li_1, ndiv=3) #Boundary Conditions, Definition der Anzahl von Elementen
mapdl.lesize(nl1=li_8, ndiv=10)

# Outer ring

ba = 0.032      #Breite Lagersitz
c = 0.00732
e = 0.010942
f1 = 0.032
g = 0.00439
da = di+2*(o+c+e)      #Außen-Durchmesser 

# Outer ring - Keypoints

mapdl.k(9, 0, di/2+o+c, 0)
mapdl.k(10, 0, da/2, 0)
mapdl.k(11, 0, da/2, f1)
mapdl.k(12, 0, da/2-g, f1)
mapdl.k(16, 0, da/2-0.00439, f1-0.000102) # Kontaktpunkte Kegelrolle, evtl. über zyl. KOS bestimmen, statt über Winkelfunktion auszurechnen!
mapdl.k(17, 0, da/2-0.00955, f1-0.025245) # Kontaktpunkte Kegelrolle, evtl. über zyl. KOS bestimmen, statt über Winkelfunktion auszurechnen!

# Outer ring - Lines

mapdl.l(9, 10)
l_a_f = mapdl.l(10, 11)
mapdl.l(11, 12)
mapdl.l(12, 16)
l1_a = mapdl.l(16, 17)
mapdl.l(17, 9)

mapdl.lesize(nl1=l1_a, ndiv=anzahl_scheiben-1) #Scheibenmodell Außenring-Lauffläche
mapdl.lesize(nl1=l1_i, ndiv=anzahl_scheiben-1) #Scheibenmodell Innenring-Lauffläche
mapdl.lesize(nl1=l1_b, ndiv=anzahl_scheiben_board-1) #Scheibenmodell Innenring-Bordkontakt

mapdl.lesize(nl1=l_a_f, ndiv=10)

# Outer ring - Surface

a2 = mapdl.a(9, 10, 11, 12, 16, 17)

# Extrudieren - Innenring & Außenring

mapdl.vrotat(a1, pax1=14, pax2=15, arc=360, nseg=anzahl_rollen) #Hier bei nseg noch Anzahl Rollen einfügen
mapdl.vrotat(a2, pax1=14, pax2=15, arc=360, nseg=anzahl_rollen)

# Roller

f2 = 0.044114
a1 = 0.004275
a2 = 0.00788
m1 = 0.025241
h = 0.006887
a3 = 0.005862
a4 = 0.006707
a5 = 0.001135
a6 = 0.006322

# Roller - Keypoints

mapdl.k(2000000, 0, f2, a2)
mapdl.k(2100000, 0, f2+a1, a2+m1)
mapdl.k(3000000, 0, f2+a3, h)
mapdl.k(3100000, 0, f2+a1+a4, a2+m1-a5)

# Roller - Lines

l1_k = mapdl.l(2000000, 3000000)
l2_k = mapdl.l(3000000, 3100000)
mapdl.l(3100000, 2100000)
l3_k = mapdl.l(2000000, 2100000)

# Roller - Surface

f_a3 = mapdl.a(2000000, 3000000, 3100000, 2100000)

mapdl.lesize(nl1=l2_k, ndiv=anzahl_scheiben-1)

# Extrude - Roller

mapdl.vsel("NONE")
k_rot = mapdl.vrotat(f_a3, pax1=2000000, pax2=2100000, arc=360, nseg=2)

mapdl.csys(1)
mapdl.vgen(itime=anzahl_rollen, nv1="ALL", dy=360/anzahl_rollen, imove=0)
mapdl.vsel("ALL")

# Meshing

mapdl.csys(0)
mapdl.vsweep("ALL")

# Creation of contact elements
# Selection of lines

mapdl.lsel("S", "NDIV", vmin=anzahl_scheiben-1)
#print(mapdl.geometry.lnum)
mapdl.lplot(vtk=True, show_line_numbering=True, background="white", color_lines="black")

polar_rollen = 360/anzahl_rollen
radius_ir = di/2+o-q+0.001
radius_ri = di/2+o-q+0.002
radius_ra = f2+a3+0.001
radius_ar = f2+a6+0.0025

x_ir = []
x_ri = []
x_ra = []
x_ar = []

y_ir = []
y_ri = []
y_ra = []
y_ar = []

for i in np.arange(0, 360+polar_rollen, polar_rollen):
    x_ir.append(radius_ir*np.sin(np.radians(i)))
    x_ri.append(radius_ri*np.sin(np.radians(i)))
    x_ra.append(radius_ra*np.sin(np.radians(i)))
    x_ar.append(radius_ar*np.sin(np.radians(i)))
    y_ir.append(radius_ir*np.cos(np.radians(i)))
    y_ri.append(radius_ri*np.cos(np.radians(i)))
    y_ra.append(radius_ra*np.cos(np.radians(i)))
    y_ar.append(radius_ar*np.cos(np.radians(i)))
    
# Inner ring

ir_nnum = []

for i in range(0, anzahl_rollen, 1):
    mapdl.lsel("ALL")
    mapdl.lsel("S", "NDIV", vmin=anzahl_scheiben-1)
    if abs(round(x_ir[i], 6)) != abs(round(y_ir[i], 6)):
        if x_ir[i] >= 0 and y_ir[i] >= 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ir[i], vmax=y_ir[i]+0.001)
            mapdl.lsel("R", "LOC", "X", vmin=x_ir[i], vmax=x_ir[i]+0.001)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ir_nnum.append(mapdl.mesh.nnum)
        if x_ir[i] >= 0 and y_ir[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ir[i], vmax=y_ir[i]-0.001)
            mapdl.lsel("R", "LOC", "X", vmin=x_ir[i], vmax=x_ir[i]+0.001)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ir_nnum.append(mapdl.mesh.nnum)
        if x_ir[i] < 0 and y_ir[i] >= 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ir[i], vmax=y_ir[i]+0.001)
            mapdl.lsel("R", "LOC", "X", vmin=x_ir[i], vmax=x_ir[i]-0.001)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ir_nnum.append(mapdl.mesh.nnum)
        if x_ir[i] < 0 and y_ir[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ir[i], vmax=y_ir[i]-0.001)
            mapdl.lsel("R", "LOC", "X", vmin=x_ir[i], vmax=x_ir[i]-0.001)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ir_nnum.append(mapdl.mesh.nnum)
    if abs(round(x_ir[i], 6)) == abs(round(y_ir[i], 6)):
        if x_ir[i] > 0 and y_ir[i] > 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ir[i], vmax=y_ir[i]+0.0005)
            mapdl.lsel("R", "LOC", "X", vmin=x_ir[i], vmax=x_ir[i]+0.0005)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ir_nnum.append(mapdl.mesh.nnum)
        if x_ir[i] > 0 and y_ir[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ir[i], vmax=y_ir[i]-0.0005)
            mapdl.lsel("R", "LOC", "X", vmin=x_ir[i], vmax=x_ir[i]+0.0005)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ir_nnum.append(mapdl.mesh.nnum)
        if x_ir[i] < 0 and y_ir[i] > 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ir[i], vmax=y_ir[i]+0.0005)
            mapdl.lsel("R", "LOC", "X", vmin=x_ir[i], vmax=x_ir[i]-0.0005)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ir_nnum.append(mapdl.mesh.nnum)
        if x_ir[i] < 0 and y_ir[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ir[i], vmax=y_ir[i]-0.0005)
            mapdl.lsel("R", "LOC", "X", vmin=x_ir[i], vmax=x_ir[i]-0.0005)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ir_nnum.append(mapdl.mesh.nnum)
         
#print(ir_nnum[0])

# Roller - Inner ring contact

ri_nnum = []

for i in range(0, anzahl_rollen, 1):
    mapdl.lsel("ALL")
    mapdl.lsel("S", "NDIV", vmin=anzahl_scheiben-1)
    if abs(round(x_ri[i], 6)) != abs(round(y_ri[i], 6)):
        if x_ri[i] >= 0 and y_ri[i] >= 0: 
            mapdl.lsel("R", "LOC", "Y", vmin=y_ri[i], vmax=y_ri[i]+0.001)
            mapdl.lsel("R", "LOC", "X", vmin=x_ri[i], vmax=x_ri[i]+0.001)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r+0.000001, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ri_nnum.append(mapdl.mesh.nnum)
        if x_ri[i] >= 0 and y_ri[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ri[i], vmax=y_ri[i]-0.001)
            mapdl.lsel("R", "LOC", "X", vmin=x_ri[i], vmax=x_ri[i]+0.001)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r+0.000001, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ri_nnum.append(mapdl.mesh.nnum)
        if x_ri[i] < 0 and y_ri[i] >= 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ri[i], vmax=y_ri[i]+0.001)
            mapdl.lsel("R", "LOC", "X", vmin=x_ri[i], vmax=x_ri[i]-0.001)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r+0.000001, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ri_nnum.append(mapdl.mesh.nnum)
        if x_ri[i] < 0 and y_ri[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ri[i], vmax=y_ri[i]-0.001)
            mapdl.lsel("R", "LOC", "X", vmin=x_ri[i], vmax=x_ri[i]-0.001)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r+0.000001, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ri_nnum.append(mapdl.mesh.nnum)
    if abs(round(x_ri[i], 6)) == abs(round(y_ri[i], 6)):
        if x_ri[i] > 0 and y_ri[i] > 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ri[i], vmax=y_ri[i]+0.0005)
            mapdl.lsel("R", "LOC", "X", vmin=x_ri[i], vmax=x_ri[i]+0.0005)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ri_nnum.append(mapdl.mesh.nnum)
        if x_ri[i] > 0 and y_ri[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ri[i], vmax=y_ri[i]-0.0005)
            mapdl.lsel("R", "LOC", "X", vmin=x_ri[i], vmax=x_ri[i]+0.0005)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ri_nnum.append(mapdl.mesh.nnum)
        if x_ri[i] < 0 and y_ri[i] > 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ri[i], vmax=y_ri[i]+0.0005)
            mapdl.lsel("R", "LOC", "X", vmin=x_ri[i], vmax=x_ri[i]-0.0005)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ri_nnum.append(mapdl.mesh.nnum)
        if x_ri[i] < 0 and y_ri[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ri[i], vmax=y_ri[i]-0.0005)
            mapdl.lsel("R", "LOC", "X", vmin=x_ri[i], vmax=x_ri[i]-0.0005)
            mapdl.lsel("R", "LOC", "Z", vmin=p+r, vmax=p+r+s/2)
            mapdl.nsll("S", 1)
            ri_nnum.append(mapdl.mesh.nnum)
            
#print(ri_nnum[0])

# Roller - Outer ring contact

ra_nnum = []

for i in range(0, anzahl_rollen, 1):
    mapdl.lsel("ALL")
    mapdl.lsel("S", "NDIV", vmin=anzahl_scheiben-1)
    if abs(round(x_ra[i], 6)) != abs(round(y_ra[i], 6)):
        if y_ra[i] >= 0 and x_ra[i] >= 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ra[i], vmax=y_ra[i]+0.0016)
            mapdl.lsel("R", "LOC", "X", vmin=x_ra[i], vmax=x_ra[i]+0.0016)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ra_nnum.append(mapdl.mesh.nnum)
        if x_ra[i] >= 0 and y_ra[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ra[i], vmax=y_ra[i]-0.0016)
            mapdl.lsel("R", "LOC", "X", vmin=x_ra[i], vmax=x_ra[i]+0.0016)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ra_nnum.append(mapdl.mesh.nnum)
        if x_ra[i] < 0 and y_ra[i] >= 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ra[i], vmax=y_ra[i]+0.0016)
            mapdl.lsel("R", "LOC", "X", vmin=x_ra[i], vmax=x_ra[i]-0.0016)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ra_nnum.append(mapdl.mesh.nnum)
        if x_ra[i] < 0 and y_ra[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ra[i], vmax=y_ra[i]-0.0016)
            mapdl.lsel("R", "LOC", "X", vmin=x_ra[i], vmax=x_ra[i]-0.0016)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ra_nnum.append(mapdl.mesh.nnum)
    if abs(round(x_ra[i], 6)) == abs(round(y_ra[i], 6)):
        if x_ra[i] > 0 and y_ra[i] > 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ra[i], vmax=y_ra[i]+0.0014)
            mapdl.lsel("R", "LOC", "X", vmin=x_ra[i], vmax=x_ra[i]+0.0014)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ra_nnum.append(mapdl.mesh.nnum)
        if x_ra[i] > 0 and y_ra[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ra[i], vmax=y_ra[i]-0.0014)
            mapdl.lsel("R", "LOC", "X", vmin=x_ra[i], vmax=x_ra[i]+0.0014)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ra_nnum.append(mapdl.mesh.nnum)
        if x_ra[i] < 0 and y_ra[i] > 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ra[i], vmax=y_ra[i]+0.0014)
            mapdl.lsel("R", "LOC", "X", vmin=x_ra[i], vmax=x_ra[i]-0.0014)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ra_nnum.append(mapdl.mesh.nnum)
        if x_ra[i] < 0 and y_ra[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ra[i], vmax=y_ra[i]-0.0014)
            mapdl.lsel("R", "LOC", "X", vmin=x_ra[i], vmax=x_ra[i]-0.0014)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ra_nnum.append(mapdl.mesh.nnum)
        
#print(ra_nnum)

# Outer ring

ar_nnum = []

for i in range(0, anzahl_rollen, 1):
    mapdl.lsel("ALL")
    mapdl.lsel("S", "NDIV", vmin=anzahl_scheiben-1)
    if abs(round(x_ar[i], 6)) != abs(round(y_ar[i], 6)):
        if y_ar[i] >= 0 and x_ar[i] >= 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ar[i], vmax=y_ar[i]+0.0012)
            mapdl.lsel("R", "LOC", "X", vmin=x_ar[i], vmax=x_ar[i]+0.0012)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ar_nnum.append(mapdl.mesh.nnum)
        if x_ar[i] >= 0 and y_ar[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ar[i], vmax=y_ar[i]-0.0012)
            mapdl.lsel("R", "LOC", "X", vmin=x_ar[i], vmax=x_ar[i]+0.0012)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ar_nnum.append(mapdl.mesh.nnum)
        if x_ar[i] < 0 and y_ar[i] >= 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ar[i], vmax=y_ar[i]+0.0012)
            mapdl.lsel("R", "LOC", "X", vmin=x_ar[i], vmax=x_ar[i]-0.0012)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ar_nnum.append(mapdl.mesh.nnum)
        if x_ar[i] < 0 and y_ar[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ar[i], vmax=y_ar[i]-0.0012)
            mapdl.lsel("R", "LOC", "X", vmin=x_ar[i], vmax=x_ar[i]-0.0012)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ar_nnum.append(mapdl.mesh.nnum)
    if abs(round(x_ra[i], 6)) == abs(round(y_ra[i], 6)):
        if x_ar[i] > 0 and y_ar[i] > 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ar[i], vmax=y_ar[i]+0.0014)
            mapdl.lsel("R", "LOC", "X", vmin=x_ar[i], vmax=x_ar[i]+0.0014)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ar_nnum.append(mapdl.mesh.nnum)
        if x_ar[i] > 0 and y_ar[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ar[i], vmax=y_ar[i]-0.0012)
            mapdl.lsel("R", "LOC", "X", vmin=x_ar[i], vmax=x_ar[i]+0.0012)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ar_nnum.append(mapdl.mesh.nnum)
        if x_ar[i] < 0 and y_ar[i] > 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ar[i], vmax=y_ar[i]+0.0014)
            mapdl.lsel("R", "LOC", "X", vmin=x_ar[i], vmax=x_ar[i]-0.0014)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ar_nnum.append(mapdl.mesh.nnum)
        if x_ar[i] < 0 and y_ar[i] < 0:
            mapdl.lsel("R", "LOC", "Y", vmin=y_ar[i], vmax=y_ar[i]-0.0014)
            mapdl.lsel("R", "LOC", "X", vmin=x_ar[i], vmax=x_ar[i]-0.0014)
            mapdl.lsel("R", "LOC", "Z", vmin=a2+0.000001, vmax=a2+m1/2)
            mapdl.nsll("S", 1)
            ar_nnum.append(mapdl.mesh.nnum)
        
#print(ar_nnum)

# Contact definition
mapdl.vsel("ALL")

# Inner ring - Roller

profil_ir = [0.000012544,	0.000010816,	0.000009216,	0.000007744,	0.0000064,	0.000005184,	0.000004096,	0.000003136,	0.000002304,	0.0000016,	0.000001024,	0.000000576,	2.56E-07,	6.4E-08,	0,	6.4E-08,	2.56E-07,	5.76E-07,	0.000001024,	0.0000016,	2.304E-06,	3.136E-06,	0.000004096,	0.000005184,	0.0000064,	0.000007744,	9.216E-06,	0.000010816,	0.000012544,	0.0000144]

ri_nnum_flip = np.roll(ri_nnum, anzahl_scheiben-1, axis=1)

mapdl.et(2, "CONTA178")
mapdl.type(2)
for a in range(0, anzahl_rollen, 1):
    for i in range(2, anzahl_scheiben+2, 1):
        mapdl.r(i, r1=-210000, r2=profil_ir[i-2])
        mapdl.real(i)
        mapdl.e(ir_nnum[a][i-2], ri_nnum_flip[a][i-2])

# Outer ring - Roller

profil_ar = [0.000012544,	0.000010816,	0.000009216,	0.000007744,	0.0000064,	0.000005184,	0.000004096,	0.000003136,	0.000002304,	0.0000016,	0.000001024,	0.000000576,	2.56E-07,	6.4E-08,	0,	6.4E-08,	2.56E-07,	5.76E-07,	0.000001024,	0.0000016,	2.304E-06,	3.136E-06,	0.000004096,	0.000005184,	0.0000064,	0.000007744,	9.216E-06,	0.000010816,	0.000012544,	0.0000144]

ar_nnum_flip = np.flip(ar_nnum, axis=1)
ra_nnum_roll = np.roll(ra_nnum, anzahl_scheiben-1, axis=1)

mapdl.et(2, "CONTA178")
mapdl.type(2)
for a in range(0, anzahl_rollen, 1):
    for i in range(2, anzahl_scheiben+2, 1):
        mapdl.r(i, r1=-210000, r2=profil_ar[i-2])
        mapdl.real(i)
        mapdl.e(ar_nnum_flip[a][i-2], ra_nnum_roll[a][i-2])

mapdl.eplot(vtk=True, show_edges=True, show_axes=True, line_width=1, show_node_numbering=False)

# Boundary conditions
# Innenring
# radial
mapdl.allsel()
mapdl.csys(1)
mapdl.nsel("S", "LOC", "X", vmin=di/2-0.001, vmax=di/2+0.001)
mapdl.d("ALL", "UX")
mapdl.d("ALL", "UY")

# axial
mapdl.allsel()
mapdl.csys(0)
mapdl.asel("S", "LOC", "Z", vmin=p+r+s+b-0.001, vmax=p+r+s+b+0.001)
mapdl.nsla("S", 1)
mapdl.d("ALL", "UZ")

# Outer ring
# axial
mapdl.allsel()
mapdl.csys(0)
mapdl.asel("S", "LOC", "Z", vmin=0, vmax=0.001)
mapdl.csys(1)
mapdl.asel("R", "LOC", "X", vmin=di/2+o+c, vmax=di/2+o+c+e)
mapdl.nsla("S", 1)
mapdl.csys(0)
mapdl.d("ALL", "UZ")

# Roller
mapdl.allsel()
mapdl.csys(0)
mapdl.vsel("S", "LOC", "Z", vmin=p+r+s/2-0.0015, vmax=p+r+s/2+0.0015)
mapdl.csys(1)
mapdl.vsel("S", "LOC", "X", vmin=f2, vmax=f2+0.005)
mapdl.nslv("S", 1)
mapdl.csys(1)
mapdl.d("ALL", "UX")
mapdl.d("ALL", "UZ")

# Forces
mapdl.allsel()
mapdl.csys(1)
mapdl.asel("S", "LOC", "X", vmin=di/2+o+c+e-0.001, vmax=di/2+o+c+e+0.001)
mapdl.sfa("ALL", 1, "PRES", value=2000)

#mapdl.allsel()
#mapdl.csys(0)
#mapdl.lsel("S", "LOC", "Y", vmin=di/2+o+c+e-0.001, vmax=di/2+o+c+e+0.001)
#mapdl.nsll()
#mapdl.sf("ALL", "PRES", 2000)

# Solve
mapdl.run("/SOL")
mapdl.antype(0)
mapdl.outres("ALL", "ALL")
mapdl.run("/STATUS,SOLU")
sol_output = mapdl.solve()
mapdl.finish()
print(sol_output)

# Plot results
result = mapdl.result
result.plot_principal_nodal_stress(
    0,
    "SEQV",
    lighting=False,
    background="w",
    show_edges=True,
    text_color="k",
    add_text=False,)

nnum, stress = result.principal_nodal_stress(0)
von_mises = stress[:, -1]  

max_stress = np.nanmax(von_mises)

# Contact Penetration
# Example to test contact penetration
mapdl.allsel()
mapdl.lsel("S", "LINE", vmin=4)
mapdl.nsll("S", 1)
mapdl.esln("S", 0)
mapdl.esel("R", "ENAM", "", 178)
mapdl.eplot(vtk=True)
contact_penetration_r = mapdl.post_processing.element_values(item='nmisc', comp=3)
print(contact_penetration_r)

mapdl.post1()
mapdl.set("last")
mapdl.allsel()

mapdl.post_processing.plot_nodal_principal_stress("1", smooth_shading=False)

mapdl.exit()

mikerife · 2023-06-20T15:32:37Z

mikerife
Jun 20, 2023
Maintainer

Hi @herttim
Ok there are two issues with the input file. The first is that the PyMAPDL Reader is being deprecated, so we should really stay away from using it in general. It's use starts with:

result = mapdl.result

The PyMAPDL result plotting includes plotting of equivalent stress, so no real need to use the PyMAPDL Reader anyway. The second issue is that the contact penetration post-processing needs to be issued while in /POST1 and after reading in the results. So move these:

# Contact Penetration
# Example to test contact penetration
mapdl.allsel()
mapdl.lsel("S", "LINE", vmin=4)
mapdl.nsll("S", 1)
mapdl.esln("S", 0)
mapdl.esel("R", "ENAM", "", 178)
mapdl.eplot(vtk=True)
contact_penetration_r = mapdl.post_processing.element_values(item='nmisc', comp=3)
print(contact_penetration_r)

To after the

mapdl.set("last")

command.

Lastly the Contact178 elements default to use the defined (by real constant) GAP + the nodes locations as the total contact gap. In this model this results in a gap that is too large and I am pretty sure is not meant in this model. Otherwise the outer housing would almost (or not almost?) need to collapse before that gap was closed. Change the keyopt (4) for the contact elements to define the gap just by the real constant value:

mapdl.keyopt(2,4,1)

Also the second use of the ET command to define element type 2 to be cont178 is redundant (and kind of dangerous to be issues more than once).

Mike

2 replies

mikerife Jun 20, 2023
Maintainer

I should have added that since this is a contact analysis then large deformation should be used (nlgeom command).

herttim Jun 22, 2023
Author

@mikerife Thank you! That helped me a lot.

koubaa · 2023-06-21T16:20:47Z

koubaa
Jun 21, 2023
Maintainer

@germa89 @mikerife This issue suggests possible improvements that could be made to pymapdl:

an option in pymapdl.solve() to do some validation before calling the SOLVE function on the back end (e.g. make sure that nlgeom is on if there are contact elements)
post-processing commands throw a helpful exception if there is no active result set loaded
ET commands can write a warning if elements are already defined with the given id

1 reply

germa89 Jun 22, 2023
Maintainer

I opened an issue because I fear this could be forgotten if left in a discussion comment only.

No contact penetration of CONTA178 elements - Bearing simulation #2123

Uh oh!

Uh oh!

herttim Jun 16, 2023

Replies: 4 comments · 7 replies

Uh oh!

pmaroneh Jun 16, 2023 Maintainer

Uh oh!

herttim Jun 16, 2023 Author

Uh oh!

pmaroneh Jun 16, 2023 Maintainer

Uh oh!

ChrisHarrold Jun 16, 2023

Uh oh!

mikerife Jun 19, 2023 Maintainer

Uh oh!

Uh oh!

herttim Jun 20, 2023 Author

Uh oh!

Uh oh!

mikerife Jun 20, 2023 Maintainer

Uh oh!

mikerife Jun 20, 2023 Maintainer

Uh oh!

herttim Jun 22, 2023 Author

Uh oh!

Uh oh!

koubaa Jun 21, 2023 Maintainer

Uh oh!

germa89 Jun 22, 2023 Maintainer

herttim
Jun 16, 2023

Replies: 4 comments 7 replies

pmaroneh
Jun 16, 2023
Maintainer

herttim Jun 16, 2023
Author

pmaroneh Jun 16, 2023
Maintainer

mikerife
Jun 19, 2023
Maintainer

herttim Jun 20, 2023
Author

mikerife
Jun 20, 2023
Maintainer

mikerife Jun 20, 2023
Maintainer

herttim Jun 22, 2023
Author

koubaa
Jun 21, 2023
Maintainer

germa89 Jun 22, 2023
Maintainer