Introducing an Open-Source Simulation Model for Track Rollers Considering Friction

Abstract

Locating bearing track rollers are used, for example, in monorail transport systems to enable relative movement between the rail and the shuttle. Due to the two-point contact, both radial and axial forces can be transmitted simultaneously. Since friction is involved, the state of the art does not provide any calculation rules for the dimensioning and design. The development of a calculation model with sophisticated commercial software brings its difficulties since no plausibility check is possible using existing models. For this reason, a model based on analytical descriptions including the Hertzian and the elastic half space theories is presented in this paper. It bridges the gap between very simple approaches and widely developed commercial software. With this model, the contact forces, friction forces, surface tensions, relative velocities and subsurface stresses can be calculated for both free and driven rolling. The main advantages are that the model is easy to apply, and thus comparisons between different track roller designs can be made quickly.

1. Introduction

Locating bearing track roller guidance systems are used, for example, in monorail transport systems. Due to the essential requirements for the systems, no additional lubricant can be used for the guideways. From a contact mechanics point of view, this is therefore a problem that requires tangential surface stresses to be taken into account [1]. The current state of research on the design of locating bearing track roller guidance systems [2,3] does not offer any satisfactory calculation options for this; thus, the author has already described the problem in an initial paper [4] and then presented a proposal for a static model that takes friction into account [5]. In the design phase, it is very important to support decisions with calculation models [6]. However, only very advanced commercial software can currently be used for this purpose.

For contact mechanics calculation problems, the software CONTACT, which uses Kalker’s theory [7], is particularly suitable [8]. Although such contact mechanics problems can also be solved using finite element methods (FEM), FEM approaches will not be considered further within the scope of this article because of their very high computing times.

The application of CONTACT to the machine element of the fixed bearing track roller requires some effort. First of all, the 2D contour geometry of the bodies must be derived from the CAD (computer-aided design) models in order to use it as an input for the calculation. For each geometry change, new contour geometries have to be derived [9]. Rapid changes to the geometry are therefore hardly applicable for the calculation, making it difficult to carry out parameter studies during the design phase.

With commercial software, the calculation process is often not visible and therefore not comprehensible. Therefore, a solution of the problem using contact mechanics modeling in combination with an open-source code is considered useful. In order to provide some plausibility to the calculated results, the gap between simple models and highly developed software must be closed by a scientifically developed, comprehensible model.

In this paper, a simulation model is presented that calculates the surface and material stresses of a continuously rolling track roller. The model is presented using the example of a selected locating bearing track roller guide (Figure 1, right), whose most important geometrical quantities are shown in

Table 1. Parameter and values of the selected example.

Parameter	${\mathbf{A}}_{\mathbf{T}}$ (mm)	${\mathbf{r}}_{\mathbf{C},\mathbf{T}\mathbf{R}}$ (mm)	${\mathbf{r}}_{\mathbf{C}\mathbf{u}\mathbf{r}\mathbf{v}\mathbf{e}}$ (mm)	${\mathbf{r}}_{\mathbf{C},\mathbf{G}}$ (mm)	$\mathsf{\alpha }$ (°)	E (Nmm)	$\mathsf{\nu }$ (-)	${\mathit{F}}_{\mathbf{r}\mathbf{a}\mathbf{d}}$ (mm)	${\mathit{F}}_{\mathbf{a}\mathbf{x}}$ (mm)	${\mathit{F}}_{\mathbf{M}\mathbf{o}\mathbf{t}\mathbf{o}\mathbf{r}}$ (N)
										Free Rolling	Driven Rolling
Value	21.75	7	220	6	30	210,000	0.33	200	30	0	−30

. The model covers both free rolling and driven rolling and is valid for Hertzian contacts with a high proportion of spin.

The simulation model can be used to calculate both straight and curved sections of the guideway. The radius ${\mathrm{r}}_{\mathrm{Curve}}$ indicates the radius of the curvature of the guideway relative to the center of the contour. A distinction between the inner and outer wheel is necessary.

2. Materials and Methods

The presented project (available at: https://github.com/janwenzel/Dynamic-Contact-Simulation-for-Locating-Bearing-Track-Rollers, accessed on 10 June 2022, Supplementary Materials) is a simulation model based on known theories of contact mechanics. The entire project is written in Python. The most important functions used are listed in the individual sections of this article and can be viewed in the repository. The documentation of the individual functions is openly accessible via the websites of the toolkits.

The calculation of the subsurface stresses is performed with the elastic half space theory. The contact area is calculated according to the Hertzian theory [10,11]. The Hertzian pressure describes the normal stress distribution at the contact of two bodies of constant radii of curvature. Such contact always results in an elliptical contact surface. In the track roller guidance system used here, constant radii of curvature are present in the contacts, which is why Hertzian theory can be applied. The contact kinematics are derived on the basis of roller and guideway geometry following the procedure of Birkhofer [12].

Additionally, as a post processing step, the Wöhler curve, as described in DIN50100 [13], is used to establish a relationship between stress and duration until damage is reached.

A material database script can be used to retrieve and edit material parameters and to add new materials.

2.1. Calculation of the Subsurface Stresses

For the calculation of the stresses below the surface, the model of the elastic half-space, as described by Johnson [14] is chosen. In this paper the further modified formulation by Nikas [15] is used. Figure 2 shows the relationship between surface boundary conditions and any point below the surface.

Here C is any point within the contact zone S at which the traction vector ($q_x,q_y,q_{\mathrm{z}}$) acts. The goal is to obtain the stresses at any point A below the surface. According to [15] the potential functions are formulated as follows

$$Q_i={\displaystyle \iint }{q}_i\left(\xi ,\eta \right)\Omega d\xi d\eta , \left(i=x,y,z\right)$$

(1)

Thereby $\Omega$ is

$$\Omega =z·ln\left(\rho +z\right)-\rho$$

(2)

The vector $\rho$ is formulated as follows

$$\rho =\sqrt{{\left(\xi -x\right)}^2+{\left(\eta -y\right)}^2+z^2}$$

(3)

The calculation of the stresses is only given here as an example for the normal stresses

$${\sigma }_i=\frac{\mathrm{E}}{1+\mathsf{\nu }}\left[\frac{\mathsf{\nu }}{1-2\mathsf{\nu }}\left(\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}\right)+\frac{\partial u_i}{\partial i}\right],\left(i=x,y,z\right)$$

(4)

The displacements are also again listed for only one component. The other formulations can be taken from [10] as well as from the source code

$$u_x=\frac{1+\mathsf{\nu }}{2\mathsf{\pi }\mathrm{E}}\left[2\frac{{\partial }^2{Q}_x}{\partial z^2}-\frac{{\partial }^2{Q}_z}{\partial x\partial z}+2\mathsf{\nu }\left(\frac{{\partial }^2{Q}_x}{\partial x^2}+\frac{{\partial }^2{Q}_y}{\partial x\partial y}+\frac{{\partial }^2{Q}_z}{\partial x\partial z}\right)-z\left(\frac{{\partial }^3{Q}_x}{\partial x^2\partial z}+\frac{{\partial }^3{Q}_y}{\partial x\partial y\partial z}+\frac{{\partial }^3{Q}_z}{\partial x\partial z^2}\right)\right]$$

(5)

The equivalent stress is calculated according to Distortion Energy Hypothesis (DEH) since Broszeit [16] was already able to show that this is more suitable than, for example, the Alternating Shear Stress Hypothesis when additional shear stresses occur at the surface. This approach was confirmed by Harris [17].

$${\sigma }_{\mathrm{DEH}}=\sqrt{\frac{1}{2}\left({\left({\sigma }_x-{\sigma }_y\right)}^2+{\left({\sigma }_z-{\sigma }_x\right)}^2+{\left({\sigma }_y-{\sigma }_z\right)}^2+6\left({\tau }_{xy}^2+{\tau }_{xz}^2+{\tau }_{yz}^2\right)\right)}$$

(6)

The stresses underneath the surface will be calculated with the formulas as described. The necessary deviations of $\Omega$ are calculated using the toolkit SymPy with the functions diff and simplify and afterwards only the results are included in the calculation script. The potential theory in the implemented form is applicable for arbitrary contact areas.

2.2. Description of the Boundary Conditions on the Contact Surface

The track roller is always in equilibrium of the applied forces and torques. In the model presented here, only stationary rolling is considered, which is why the inertia forces in the force and torque equilibria can be neglected.

2.2.1. Force and Torque Equilibria

Figure 3 shows all forces and torques acting on the roller, always pointing in the direction of the positive coordinate axes. The coordinate systems are shown in Figure 4. The aim of the calculation model presented here is to determine the stresses, forces and frictional torques acting in the two contacts for given forces on the roller axis.

The radial force $F_{\mathrm{rad}}$, axial force $F_{\mathrm{ax}}$ and motor force $F_{\mathrm{Motor}}$ are specified here. These external loads cause the tangential forces $F_{\mathrm{T}1}$ and $F_{\mathrm{T}2}$ perpendicular to the running direction and the driving forces $F_{\mathrm{A}1}$ and $F_{\mathrm{A}2}$ parallel to the running direction. In addition, the frictional torques $M_{\mathrm{R}1}$ and $M_{\mathrm{R}2}$ in the contact surfaces, the torque about the rolling axis $M_{{\mathrm{x}}_{\mathrm{R}}}$ and the tilting torques $M_{{\mathrm{y}}_{\mathrm{R}}}$ and $M_{{\mathrm{z}}_{\mathrm{R}}}$ are caused. Furthermore, the lever arms $h_{\mathrm{N}}$ and $h_{\mathrm{T}}$ are required to establish the torque equilibria.

First, the equilibria of forces are described

$$\cos \left(\mathsf{\alpha }\right)·\left(F_{\mathrm{T}1}+F_{\mathrm{T}2}\right)+\sin \left(\mathsf{\alpha }\right)·\left(F_{\mathrm{N}1}-F_{\mathrm{N}2}\right)=-F_{\mathrm{ax}}$$

(7)

$$F_{\mathrm{A}1}+F_{\mathrm{A}2}=-F_{\mathrm{Motor}}$$

(8)

$$\sin \left(\mathsf{\alpha }\right)·\left(-F_{\mathrm{T}1}+F_{\mathrm{T}2}\right)+\cos \left(\mathsf{\alpha }\right)·\left(F_{\mathrm{N}1}+F_{\mathrm{N}2}\right)=-F_{\mathrm{rad}}$$

(9)

Subsequently, the description of the torque equilibria takes place

$$\left(M_{\mathrm{R}1}-M_{\mathrm{R}2}\right)·\sin (\mathsf{\alpha })+\left(F_{\mathrm{A}1}+F_{\mathrm{A}2}\right)·{\mathrm{r}}_{\mathrm{A}}=M_{x_{\mathrm{R}}}$$

(10)

$$h_{\mathrm{N}}·\left(F_{\mathrm{N}1}-F_{\mathrm{N}2}\right)+h_{\mathrm{T}}·\left(F_{\mathrm{T}1}+F_{\mathrm{T}2}\right)=M_{y\mathrm{R}}$$

(11)

$$\left(M_{\mathrm{R}1}+M_{\mathrm{R}2}\right)·\cos \left(\mathsf{\alpha }\right)+\left(F_{\mathrm{A}1}-F_{\mathrm{A}2}\right)·{\mathrm{r}}_{\mathrm{A}}·\sin \left(\mathsf{\alpha }\right)=-M_{z_{\mathrm{R}}}$$

(12)

The following applies to the lever arms

$$h_{\mathrm{N}}={\mathrm{A}}_{\mathrm{T}}·\sin \left(\mathsf{\alpha }\right)$$

(13)

$$h_{\mathrm{T}}={\mathrm{A}}_{\mathrm{T}}·\cos \left(\mathsf{\alpha }\right)-{\mathrm{r}}_{\mathrm{C},\mathrm{G}}$$

(14)

2.2.2. Relative Motion of the Contact Area

To describe the relative velocity in the contact area, two assumptions must be made in advance:

• Due to the high proportion of spin slip, no distinction is made between adhesive and sliding areas. Complete sliding in the entire contact area without micro-slip components is assumed.
• In addition, it is assumed that the contact area occurs in the undeformed guide surface. The assumption is therefore made of an ideally hard guide and a relatively soft track roller.

To describe the kinematics, and also for the later presentation of the results, the coordinate systems are first defined (Figure 4, left). The R coordinate system (R-COS) is located on the central roller axis in the symmetry plane of the roller. The C coordinate systems (C-COS) are located in the center of each contact surface. The C1 coordinate system (C1-COS) is located in the positive $x_{\mathrm{R}}$ direction.

The relative velocities in contact are determined on the basis of the geometrically-induced drilling and rolling slip components. The decomposition of the angular velocities is shown in Figure 4, right.

The angular velocity is first formulated as a vector in R-COS

$${\overrightarrow{\omega }}_{\mathrm{Roll}}^{\mathrm{R}}=\left(\begin{array}{c}-{\omega }_{\mathrm{r}}\\ 0\\ 0\end{array}\right)$$

(15)

The same applies to the radius

$${\overrightarrow{r}}_{\mathrm{Guideway}}^{\mathrm{R}}=\left(\begin{array}{c}{x}_{\mathrm{R}}\\ y_{\mathrm{R}}\\ \sqrt{{\mathrm{r}}_{\mathrm{C},\mathrm{G}}^2-x_{\mathrm{R}}^2}-{\mathrm{A}}_{\mathrm{T}}\end{array}\right)$$

(16)

Thus, the velocity of the contact point in the R-COS can be described

$${\overrightarrow{v}}_{\mathrm{Roll}}^{\mathrm{R}}={\overrightarrow{\omega }}_{\mathrm{Roll}}^{\mathrm{R}}\times {\overrightarrow{r}}_{\mathrm{Guideway}}^{\mathrm{R}}=\left(\begin{array}{c}0\\ {\omega }_{\mathrm{r}}·\left(\sqrt{{\mathrm{r}}_{\mathrm{C},\mathrm{G}}{}^2-x_{\mathrm{R}}{}^2}-{\mathrm{A}}_{\mathrm{T}}\right)\\ -{\omega }_{\mathrm{r}}·y_{\mathrm{R}}\end{array}\right)\mathrm{a}$$

(17)

The velocity of the contact point in the C-COS is obtained by rotation around the $y_{\mathrm{R}}$ axis

$${\overrightarrow{v}}_{\mathrm{Roll}}^{\mathrm{C}}=\left(\begin{array}{c}\sin \left(\mathsf{\alpha }\right)\cdot {\omega }_{\mathrm{r}}\cdot y_{\mathrm{R}}\\ {\omega }_{\mathrm{r}}\cdot \left(\sqrt{{\mathrm{r}}_{\mathrm{C},\mathrm{G}}{}^2-x_{\mathrm{R}}{}^2}-{\mathrm{A}}_{\mathrm{T}}\right)\\ -\cos \left(\mathsf{\alpha }\right)\cdot {\omega }_{\mathrm{r}}\cdot y_{\mathrm{R}}\end{array}\right)$$

(18)

To be able to determine the relative velocity, the velocity of the guide in the C-COS must first be described

$${\overrightarrow{v}}_{\mathrm{Guideway} }^{\mathrm{C}}=\left(\begin{array}{c}0\\ {\omega }_{\mathrm{r}}\cdot \left({\mathrm{r}}_{\mathrm{C},\mathrm{G}}\cdot \cos (\mathsf{\alpha })-{\mathrm{A}}_{\mathrm{T}}\right)\\ 0\end{array}\right)\mathrm{x}$$

(19)

Thus the relative velocity can be calculated. Furthermore the R-coordinates are substituted by C coordinates

$${\overrightarrow{v}}_{\mathrm{rel}}^{\mathrm{C}}={\overrightarrow{v}}_{\mathrm{Roll}}^{\mathrm{C}}-{\overrightarrow{v}}_{\mathrm{Guideway}}^{\mathrm{C}}=\left(\begin{array}{c}\sin \left(\mathsf{\alpha }\right)\cdot {\omega }_{\mathrm{r}}\cdot y_{\mathrm{C}}\\ {\omega }_{\mathrm{r}}\cdot \left(\sqrt{{\mathrm{r}}_{\mathrm{C},\mathrm{G}}{}^2-{\left(x_{\mathrm{C}}\cdot \cos \left(\mathsf{\alpha }\right)+{\mathrm{r}}_{\mathrm{C},\mathrm{G}}\cdot \sin \left(\mathsf{\alpha }\right)\right)}^2}-{\mathrm{r}}_{\mathrm{C},\mathrm{G}}\cdot \cos \left(\mathsf{\alpha }\right)\right)\\ -\cos \left(\mathsf{\alpha }\right)\cdot {\omega }_{\mathrm{r}}\cdot y_{\mathrm{C}}\end{array}\right)$$

(20)

Comparisons with findings from the literature and comparative simulations have shown that further assumptions have to be made about the location of the spin center point in order to adjust the relative velocity distribution.

Since the contact surface is assumed to be an ideal plane, the relative velocities in the z-direction have no effect on the stresses induced by friction. For this reason, the relative velocity in z-direction is not considered further.

Assumption 1: Location of the Spin Center

As can be seen from the literature, the spin center point can be located inside or outside the contact, depending on the load case. It is therefore not coupled to the contact center point. According to Kalker’s PhD thesis [18], Lutz and Wernitz already made the assumption that the spin center is located on the contact edge. In the left part of Figure 5 it is schematically shown how the spin center point shifts depending on the load situation in the contact area. The right part of Figure 5 also shows the shifted spin center and the associated change in the theoretical rolling radius.

The displacement of the spin center can be formulated as follows

$${\overrightarrow{v}}_{\mathrm{rel},1}^{\mathrm{C}}=(\begin{array}{c}\sin \left(\mathsf{\alpha }\right)·{\omega }_{\mathrm{r}}·\left(y_{\mathrm{C}}-{\mathit{y}}_{\mathbf{C},\mathbf{0}spin}\right)\\ {\omega }_{\mathrm{r}}·\left(\sqrt{{\mathrm{r}}_{\mathrm{C},\mathrm{G}}{}^2-{\left(x_{\mathrm{C}}·\cos \left(\mathsf{\alpha }\right)+r_{\mathrm{C},\mathrm{G}}-\sin \left(\mathsf{\alpha }\right)\right)}^2}-{\mathrm{r}}_{\mathrm{C},\mathrm{G}}·\cos \left(\mathsf{\alpha }\right)+{\mathit{x}}_{\mathbf{C},\mathbf{0}spin}·\mathbf{sin}\left(\mathsf{\alpha }\right)\right)\end{array})$$

(21)

with

$$\begin{gathered} \mathrm{If}-a<x_{\mathrm{C},0\mathrm{spin}}<a: y_{\mathrm{C},0\mathrm{spin}}=-b·\sqrt{1-{\left(\frac{x_{\mathrm{C},0\mathrm{spin}}}{a}\right)}^2} \\ \mathrm{Else}: y_{\mathrm{C},0\mathrm{spin}}=0 \end{gathered}$$

(22)

and

$$x_{\mathrm{C}1,0\mathrm{spin}}=-x_{\mathrm{C}2,0\mathrm{spin}}=\frac{r_{\mathrm{Delta}0}}{\sin (\mathsf{\alpha })}$$

(23)

The assumption improves the distributions of the relative motion and the shear stresses, but the distributions do not yet match the observations from comparative simulations.

Assumption 2: Non-concentric Distribution of the Relative Motions in the Contact Area

A second assumption ensures even more realistic stress distributions. Assumption 2 was implemented purely on the basis of observations from finite element analysis (Abaqus) and calculations with a commercial contact mechanics software (CONTACT). This assumption is not known in the literature. Assumption 2 offers an empirically determined correction of the results. This correction was determined on the basis of the roller and guideway profiles used in this application. An application of assumption 2 to other contact geometries may therefore lead to errors.

Assumption 2 states that the velocity distribution is not formed purely concentrically around the spin center. Instead, an additional influence of the contact coordinate $x_{\mathrm{C}}$ on the $x$-component of the relative velocity is assumed.

Final description of the relative motion

$${\overrightarrow{v}}_{\mathrm{rel},2}^{\mathrm{C}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)=\left(\begin{array}{c}{v}_{\mathrm{rel}, x_{\mathrm{C}}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)\\ v_{\mathrm{rel},y_{\mathrm{C}}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)\end{array}\right)=\left(\begin{array}{c}\sin \left(\mathsf{\alpha }\right)\cdot {\omega }_{\mathrm{r}}\cdot \left(y_{\mathrm{C}}-y_{C,delta}\right)\\ {\omega }_{\mathrm{r}}\left(\sqrt{{\mathrm{r}}_{\mathrm{C},\mathrm{G}}{}^2-{\left(x_{\mathrm{C}}\cdot \cos \left(\mathsf{\alpha }\right)+r_{\mathrm{C},\mathrm{G}}\cdot \sin \left(\mathsf{\alpha }\right)\right)}^2}-{\mathrm{r}}_{\mathrm{C},\mathrm{G}}\cdot \cos \left(\mathsf{\alpha }\right)+x_{\mathrm{C},0\mathrm{spin}}\cdot \sin \left(\mathsf{\alpha }\right)\right)\end{array}\right)$$

(24)

$$\begin{gathered} \mathrm{If} x_{\mathrm{C}}<x_{\mathrm{C},0\mathrm{spin}}: y_{\mathrm{C},\mathrm{delta}}=\frac{\left(x_{\mathrm{C}}+a\right)}{x_{\mathrm{C},0\mathrm{spin}}+a}·y_{\mathrm{C},0\mathrm{spin}} \\ \mathrm{If} x_{\mathrm{C}}>x_{\mathrm{C},0\mathrm{spin}}: y_{\mathrm{C},\mathrm{delta}}=\frac{\left(x_{\mathrm{C}}-a\right)}{x_{\mathrm{C},0\mathrm{spin}}-a}·y_{\mathrm{C},0\mathrm{spin}} \\ \mathrm{Else}: y_{\mathrm{C},\mathrm{delta}}=0 \end{gathered}$$

(25)

2.3. Stresses at the Contact Surface

The stress calculation is described below. The calculation of the pressure as described by Hertz [10]

$$p\left(x_C,y_{\mathrm{C}}\right)=q_{\mathrm{z}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)=p_{\max }·\sqrt{1-{\left(\frac{x_{\mathrm{C}}}{a}\right)}^2-{\left(\frac{y_{\mathrm{C}}}{b}\right)}^2}$$

(26)

Condition for the calculation of the shear stresses with the assumption of total slip in the contact area, wherein $\mu$ is the coefficient of friction

$${\tau }_{\mathrm{abs}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)=\mu ·p\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)$$

(27)

The coefficient of friction could also be replaced by more advanced models since it is often slip-velocity dependent $\mu \left(x_{\mathrm{C}},y_{\mathrm{C}}\right)=f\left({\overrightarrow{v}}_{\mathrm{rel},2}^{\mathrm{C}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)\right)$ [19]. However, this is currently not implemented. The program structure has to be slightly modified for this since the friction coefficient must then be determined iteratively.

Finally, for directional calculation of the surface tangential stresses, the relative velocity directions $v_{\mathrm{rel},\mathrm{dir}, x_{\mathrm{C}}}$ and $v_{\mathrm{rel},\mathrm{dir}, y_{\mathrm{C}}}$ are needed

$$v_{\mathrm{rel},\mathrm{dir}, x_{\mathrm{C}}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)=\frac{v_{\mathrm{rel}, x_{\mathrm{C}}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)}{||{\overrightarrow{v}}_{\mathrm{rel},2}^{\mathrm{C}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)||}$$

(28)

$$v_{\mathrm{rel},\mathrm{dir}, y_{\mathrm{C}}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)=\frac{v_{\mathrm{rel}, y_{\mathrm{C}}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)}{||{\overrightarrow{v}}_{\mathrm{rel},2}^{\mathrm{C}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)||}$$

(29)

This leads to the tangential stresses within the contact zone

$${\overrightarrow{\tau }}^C=\left(\begin{array}{c}{\tau }_{{\mathrm{x}}_{\mathrm{C}}}\\ {\tau }_{{\mathrm{y}}_{\mathrm{C}}}\end{array}\right)=\left(\begin{array}{c}{q}_{\mathrm{x}}\\ q_{\mathrm{y}}\end{array}\right)=\left(\begin{array}{c}-v_{\mathrm{rel},\mathrm{dir}, x_{\mathrm{C}}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)·{\tau }_{\mathrm{abs}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)\\ -v_{\mathrm{rel},\mathrm{dir},y_{\mathrm{C}}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)·{\tau }_{\mathrm{abs}}\left(x_{\mathrm{C}},y_{\mathrm{C}}\right)\end{array}\right)$$

(30)

With the help of the presented equations, the whole contact area can be described. These in turn form the input variables for the subsurface stress calculation.

Calculation Procedure

The actual calculation of the surface parameters is performed by an optimization process. The optimizer (Toolkit: SciPy, scipy.optimize.brentq) iterates until the zeros of the quantity to be optimized are found. The drive torque is used for this purpose since it is uniquely specified. The difference between the calculated drive torque $M_{{\mathrm{x}}_{\mathrm{R}}\mathrm{ist}}$ and the specified drive torque $M_{{\mathrm{x}}_{\mathrm{R}}\mathrm{soll}}$ must be zero

$$M_{{\mathrm{x}}_{\mathrm{R}},\mathrm{Delta}}=M_{{\mathrm{x}}_{\mathrm{R}}\mathrm{ist}}-M_{{\mathrm{x}}_{\mathrm{R}}\mathrm{soll}}$$

(31)

The deviation $r_{\mathrm{Delta}0}$ of the actual rolling radius from the nominal rolling radius ${\mathrm{r}}_{\mathrm{A}}$ is used as input variable for the optimization function. $r_{\mathrm{Delta}0}$ is, thus, the variable to be iterated. The optimization iterations are stopped when the error $M_{{\mathrm{x}}_{\mathrm{R}},\mathrm{Delta}}$ falls below a limit value.

3. Results

With the help of the presented calculation model, different results regarding the contact area and the material stress can be generated. An exemplary output is shown in Figure 6. This can be seen as a guide to the interpretation of the output graphs when using the actual program. The first output is a multi-plot that contains the shear stress distributions of both contact patches and the relative frictional energy distribution. Furthermore, in the tangential stress distribution, the spin center is represented as a black point. In addition, the stress distributions below the surface are plotted in two section planes.

Figure 7 shows the equivalent stress distribution below the surface for a frictionless ($\mu =0.0001, \mu$ cannot be set to 0 due to the formulation of the equations in the current program version) condition. It can be seen that the stress maximum is below the surface, as expected. In Figure 7 below, the same calculation has been carried out for increased friction. As already described by Hamilton and Goodman [20] (only pure sliding in one direction), a shift in the stress field can be detected due to the increased friction. For very small z values ($z\to 0$), inaccuracies occur. Thus, it is recommended to consider the values near the surface only qualitatively.

Since the calculation is performed numerically, a suitable discretization must be set. Higher discretization leads to more accurate results and fewer jumps in the near-surface regions, but the computation time increases. However, the differences below the surface are small, so a lower discretization is sufficient for the regions far from the surface.

Table 2. Discretization, Duration and Accuracy.

Parameter	${\mathit{n}}_{\mathbf{x}\mathbf{C}}$ (-)	${\mathit{n}}_{\mathbf{x}}$ (-)	${\mathit{n}}_{\mathbf{y}}$ (-)	${\mathit{n}}_{\mathbf{z}}$ (-)	${\mathit{t}}_{\mathbf{s}\mathbf{i}\mathbf{m}}$ (min)	${\mathit{\sigma }}_{\mathbf{D}\mathbf{E}\mathbf{H},\mathbf{m}\mathbf{a}\mathbf{x}}$ (MPa)
High	120	120	$x_{\mathrm{n}}·\frac{b}{a}$	31	$\approx 40$	$389.7$
Low	40	20	$x_{\mathrm{n}}·\frac{b}{a}$	12	$\approx 2$	$389.8$

lists the differences between a high and a low discretization. The equivalent stress value differs only marginally, while the calculation duration for the low discretization is significantly lower.

Here, $n_{\mathrm{xC}}$ is the discretization count on the surface, $n_{\mathrm{x}}$, $n_y$ and $n_{\mathrm{z}}$ are the discretization counts for the subsurface stress calculation; $t_{\mathrm{sim}}$ is the calculation duration for the simulation and ${\sigma }_{\mathrm{DEH},\max }$ is the maximum equivalent stress value.

The results show that the contact surface stresses as well as the subsurface stresses can be calculated within an adequate time of approximately 2 min (Intel Core i5-6200U CPU 2.30 GHz). Additionally, the influence of increased friction can be shown.

4. Discussion

For the discussion of the results, the simulation results from the presented project are compared with an established software for the calculation of contact mechanical and tribological issues. CONTACT is used as the comparison software. The calculation is performed with the help of CONTACTs Module 1 [9]. In contrast to the algorithm of the project presented within this article, where the two contact locations are predefined, CONTACT is able to find the contact patches automatically based on the contour geometry. The two calculations thus differ in this modelling aspect. The same case is calculated with both programs: $F_{\mathrm{rad}}=200 \mathrm{N},F_{\mathrm{ax}}=30 \mathrm{N}, F_{\mathrm{Motor}}=30 \mathrm{N},\mu =0.3$.

In particular, the comparison of the contact forces and maximum stresses shows that the results are very close to each other (Figure 8). The largest deviation is shown by the force $F_{\mathrm{T}2}$. The reason for this difference can be seen in Figure 9. The shear stresses have both positive and negative components. In the solution calculated with CONTACT, the negative shear stress components are larger, which is why the force is smaller in total.

A comparison of the normal and tangential stresses is shown in Figure 9. It can be seen that the basic distributions of the stresses are very similar. As already described, the largest difference is seen in the shear stress perpendicular to the rolling direction (x-direction). This can be explained by the position of the spin center.

With the help of the directional representation of the shear stress distribution (top of Figure 10), it can be clearly seen that the assumptions made are plausible. The position of the spin center can also be found in the CONTACT results; however, it is not located on the edge of the ellipse but is shifted further towards the interior of the contact surface. The main difference is mainly in the relative velocity distribution around the spin center. CONTACT calculates the velocities iteratively and can distinguish between sliding and adhesion areas. Around the spin center, there are only very low to no relative velocities. This difference significantly influences the distribution of the frictional energy.

The comparison of the equivalent stresses is shown in Figure 11. As with the stress distribution at the contact surface, there are only very slight differences here. In particular, the position and the height of the maxima are very close to each other. The weaknesses of the presented program are mainly in the areas close to the contact surface. Inaccuracies occur here.

5. Conclusions

The project presented takes up the existing state of research on the design of track roller guideways and extends it to include existing contact mechanics models so that more in-depth statements can be made on the stress state of the contact pairs, such as the influence of increased friction. The model can be used primarily for basic evaluation and the predetermination of locating bearing track roller guidance systems with two-point contacts. Through the combination of analytical description and open-source code, every step of the calculation, and thus also the results, can be reproduced. The presented calculation model therefore closes the gap in the state of research and can be used as a basis for the interpretation of more in-depth analyses with commercial software. Due to the disclosure of the entire source code, other interested parties can also use only parts of the code for their own projects and questions, such as the calculation of the subsurface stresses.

However, especially when describing the relative velocities in the contact zone, further assumptions are made to correct the simple analytical description. Despite this correction, there are still deviations in comparison to more advanced commercial software, especially in the case of the frictional energy. From the authors’ point of view, the calculation model is thus mainly suitable for surface stress calculation and determination of the equivalent stresses below the surface. To obtain more accurate results for the stress values for small z values, the calculation procedure must be further optimized.

The model can also be used to calculate material pairings with different materials. Here, the comparison with established commercial software showed that, although the shape deviates further from the ideal ellipse with increasing contact areas, there is still very high agreement in the stresses.

Due to the poor data situation in the current state of research, experimental investigations are being sought in order to further evaluate the model. For the development of a universal calculation model for track roller guidance systems, the extension to non-Hertzian contacts is necessary. From the authors’ point of view, this is associated with a high development and programming effort. However, due to the modular structure of the source code, it would be possible to replace the contact area calculation and the relative velocity calculation with more complex models, such as Kalker’s theory.