Learning Forum

[tim.hertinger]

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.

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)))
  

### Selection of Nodes on selected Lines between future contact (Not shown in this post) ###

# 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)

# 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
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()

[Ashish Khemka]

Hi,

 

You can search for the following keywords: 'tapered roller bearing ansys' and you can get a video on tapered roller bearing in  Ansys Mechanical.

 

May be a ds.dat file on this analysis might help.

 

Regards,

Ashish Khemka

[mrife]

Hi Tim

The default of Contact 178 gap is to use the node locations plus the real constant defined GAP as the total contact gap.  This is not appropriate for the way this model is set up (otherwise that is a huge contact gap to close!).  Change the element type 2 keyoption to use just the real constant GAP as the contact gap:

keyopt,2,4,1

Also the ET command was issued twice to define element type 2 (guessing copy/pasted the code block).  The command should only be issued once per element type otherwise could easily run into issues.  So take out the second instance.

Next the "Contact Penetration" commands are post processing commands that should be moved to after the mapdl.set("last") command.

Lastly the PyMAPDL Reader, used via the result=mapdl.result is being deprecated.  In the future either post process with native PyMAPDL commands or use PyDPF.  

Mike

[mrife]

Should have added that since this is a contact analysis then it should be run as non-linear i.e. use the nlgeom command to turn on large deformations. 

[tim.hertinger]

I am a bit confused that there are still no stresses on roller and inner ring raceway. I have updated my model as described above and can obtain the contact penetration now. But the stresses are looking the same. Any reasons for that?

 

 

 

 

[mrife]

 

Because the stresses in the bearings are very small compared to the whole model.  Unselect the housing volumes, then reselect the elements/nodes attached to the remaining volumes, and replot the stress….what happens?

In the image above you can kindof see the bearing stresses as one pixel...or a line of single pixels for the top bearing.  Removing the housing volumes/elements will help see these as well.