# Radial roller bearings with flexible rings: application to rotor dynamics and extension to multibody simulations

Vladimir Ivannikov, Mikhail Leontiev, Sergey Degtyarev

► To cite this version:

Vladimir Ivannikov, Mikhail Leontiev, Sergey Degtyarev. Radial roller bearings with flexible rings: application to rotor dynamics and extension to multibody simulations. 17th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery (ISROMAC2017), Dec 2017, Maui, United States. hal-03004955

HAL Id: hal-03004955

https://hal.science/hal-03004955

Submitted on 13 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Radial roller bearings with flexible rings: application to rotor dynamics and extension to multibody simulations

Mikhail Leontiev<fnref n="1" />, Vladimir Ivannikov<fnref n="2" />*, Sergey Degtyarev<fnref n="2" />

<figure><img src="image_2.png" /></figure>

ISROMAC 2017

International Symposium on Transport Phenomena and Dynamics of Rotating Machinery

Maui, Hawaii

December 16-21, 2017

**Abstract**

To obtain reliable assessments of the roller bearing durability the credible values of the radial forces are required. Most of the existing bearing models omit structural deformations of rings. This assumption does not hold for modern jet engines, where shafts are commonly hollow and the housing and the rings thicknesses may be of comparable magnitudes. Furthermore, if a considerable negative internal clearance occurs at the bearing at its operating regime and the bearing rings are said to be rigid, the computed radial contact forces may be significantly overestimated. To account for this phenomenon, a simplified deformation pattern of a ring is considered, which implies its discretization into finite segments subjected to a constant tensile/compressive force and a constant bending moment. The model key relations are derived from the energy functional that describes equilibrium of a roller located between two rings. The obtained expressions are used to build a specific 2-nodes bearing finite element, which can be conveniently applied for rotor dynamics problems. A numerical example is presented to demonstrate efficacy and performance of the approach. The results are verified with 3D finite element simulations in NX Nastran package.

**Keywords**

radial roller bearing — flexible rings — rotor dynamics — jet engines — structural deformations

## 1. INTRODUCTION

The problem of the internal bearing forces analysis has been well studied in the last 50 years, since their accurate determination is an important departure point for credible subsequent durability assessment. Since a complete solution of this problem involves various physical phenomena (i.e., elastic deformations, contact interactions, inertia effects, elasto-hydrodynamics), numerous mathematical models, which reproduce one or another aspect of the bearing behavior, have been developed. In the majority of the existing approaches the bearing rings are treated as rigid. However, this assumption is not valid, for instance, for modern jet engines, whose shafts are commonly hollow and the housing and the outer ring thicknesses may be of comparable magnitudes. Compliance of such bearing supports should not be omitted in simulations. Particularly, if a considerable negative internal clearance occurs at the bearing at its operating regime, the computed radial forces may be significantly overestimated when structural deformations are neglected.

In general case the deformed state of the bearing rings is difficult to determine and a spatial finite element (FE) problem with large number of degrees of freedom has to be solved to this end. To handle it correctly, besides the accurate modeling of contact between the interacting bodies, one needs to pay special attention to the kinematic constraints to be enforced onto the rings. Imposition of the kinematic boundary conditions directly to the rings surfaces commonly

turns out to be problematic. In order to circumvent this issue, the rotor components neighbouring to the bearing are added to the three dimensional assembly. The kinematic restraints are easier to define in this case. Such approach is safer and less prone to modeling errors and allows one to expect results of higher quality, but it is computationally expensive and nearly unacceptable for rotor dynamics problems of gas turbine engines with several bearings.

One of the first attempts to perform analysis of the bearing accounting for its structural deformations was made by Jones and Harris [12], who incorporated the outer ring flexibility into the model, although for a very particular case — the support of an idler gear in a planetary transmission. Later Filetti and Rumbarger [7], developing a general method for predicting the influence of structural support upon bearing performance, handled the outer ring as a set of beam elements and used 1D linear springs for the rollers in contact. The assembly of the ring and the rollers itself had two connection nodes with the ground support. This model proved to have good correlation with experimental data, however, a priori knowledge about the attachment points of the outer ring significantly limits the approach applicability. Moreover, as with the formulation of Jones and Harris [12], in the model of Filetti and Rumbarger [7] only one of the bearing rings is deformable. Cheng [5] advanced this approach and added to the finite element assembly not only the outer ring itself, but also the neighbouring regions of the housing. Harris [10]

<footnote n="1">¹Moscow Aviation Institute, Moscow, Russia</footnote>

<footnote n="2">²Alfa-Tranzit Co. Ltd, Moscow, Russia</footnote>

<footnote>*Corresponding author: vladimir@civil.ist.utl.pt</footnote>
has also recently studied the influence of the bearing races compliance, although limiting the discussion to some specific cases of external loading.

The recent model of Cavallaro et al. [4], further extended by Leblanc et al. [13], considers deformations of both bearing rings. Distinct from previous formulations, the problem of determination of the rings deformations was formulated in terms of forces. The model also considered the centrifugal expansion of the rings, the cage rotation and the lubricant presence. The ring out-of-roundness produced by a single line contact load was expressed through a Fourier series expansion and the effects of these contributions were superimposed afterwards to draw the resultant deformed shape of the ring. The Newton-Raphson method is then used to simultaneously solve the equilibrium and geometric equations. However, the approach provides reliable results only for bearings with positive clearance in the operating regime.

In the present work a simplified deformation pattern of the bearing rings is introduced to incorporate structural flexibility. To this end a ring is discretized into finite segments subjected to a constant tensile/compressive force and a constant bending moment. The centrifugal load acting on the roller is also considered, while inner friction and cage interactions are omitted. To derive the roller equilibrium equations and the subsequent bearing resultant force vector along with the corresponding tangent stiffness matrix, an expression for the total energy was constructed for each roller separately, and then this functional was minimized. In addition to the radial deflections of the rings, the formulation also permits their misalignments. Applicability of the approach to bearings operating under different conditions (with positive or negative clearance) is discussed and an alternative, simplified structural FE model is proposed for verification. The formulation efficacy is demonstrated in numerical tests.

## 2. RADIAL BEARING MODEL

### 2.1 Bearing equilibrium

Consider a bearing as an element with two nodes A and B attached to centers of the inner and outer rings, respectively, see Figure 1. We introduce a cartesian system x, y, z, whose origin coincides with node A. Its unitary vectors are $e_{x0}$, $e_{y0}$, $e_{z0}$. Both nodes can move under applied load, however, it is possible to recompute the global deflections into relative values keeping one of the rings fixed. In the present case the inner one is chosen the reference. The relative spatial position of the outer ring is then described by small displacements $u_x$, $u_y$ and rotation angles $\gamma_x$, $\gamma_y$, gathered in the vector

$$ \mathbf{y}_{rel} = [u_x \ u_y \ 0 \ \gamma_x \ \gamma_y \ 0]^T, \quad (1) $$

which is related to the externally applied force vector

$$ \mathbf{F} = [F_x \ F_y \ 0 \ M_x \ M_y \ 0]^T. \quad (2) $$

The proposed model is incapable to carry axial load, thus the deflection $u_z$ is not taken into account. Axial rotation $\gamma_z$ is

<figure><img src="image_3.png" /><figcaption>Figure 1. Bearing coordinate system and permitted nodal displacements and rotations.</figcaption></figure>

also absent due to bearing mechanics. The corresponding force vector components are nullified.

Let us introduce a local coordinate system r, n, z for the j-th roller, such that axis r is aligned in the radial direction, axis z is inherited from the bearing node A and axis n = r × z. In this coordinate system, displacements vector (1) is written as

$$ \mathbf{y}_{local\ j} = [u_r \ \gamma_n]^T = \mathbf{S}_j \mathbf{y}_{rel}, \quad (3) $$

where the rotation matrix

$$ \mathbf{S}_j = \begin{bmatrix} \cos \psi & \sin \psi & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\sin \psi & \cos \psi & 0 \end{bmatrix} \quad (4) $$

was introduced. For statics analysis the position of the j-th roller is constant and defined as $\psi = j\theta$ with $\theta = \frac{2\pi}{N}$, where N is the number of rollers. In transient dynamics problems angle $\psi$ varies according to the rotation phase of the outer ring relatively to the inner one, thus an arbitrary roller may periodically enter and leave the contact zone.

The outer ring cross-section at a particular rolling element j is loaded by the force vector

$$ \mathbf{P}_{local\ j} = [P_r \ M_n]^T \quad (5) $$

and the bearing equilibrium can thus be written as

$$ \mathbf{F} + \sum_{j=1}^{N} \mathbf{S}_j^T \mathbf{P}_{local\ j} = \mathbf{0}. \quad (6) $$

The constructed set of equations is nonlinear and the iterative Newton-Raphson scheme is invoked to determine the unknown vector $\mathbf{y}_{rel}$ for a given force vector $\mathbf{F}$:

$$ \sum_{j=1}^{N} \mathbf{S}_j^T \mathbf{K}_{local\ j} \mathbf{S}_j \Delta \mathbf{y}_{rel} = \mathbf{F} + \sum_{j=1}^{N} \mathbf{S}_j^T \mathbf{P}_{local\ j}, \quad (7) $$

where the local roller stiffness matrix $\mathbf{K}_{local}$ was introduced.
Equilibrium equation (6) does not pose any limitations on the nature of external forces. Similarly, vector $P_{localj}$ represents the resultant force and moment acting on the outer ring at the location of the j-th roller and is independent from the effects taken into account at the rolling element level. The model describing the interaction of the roller and the rings can be specifically altered by the researcher to meet one's needs. The rolling element may have different number of local degrees of freedom, include or neglect inertia and gyroscopic effects, consider Hertzian or non-Hertzian contact formulations, etc.

## 2.2 Deformation pattern of a bearing ring

For simplicity, let us describe the deformed state of the bearing inner ring; analogous considerations hold for the outer one as well. A ring can be treated as an initially curved rod under a set of radial forces, see Figure 2a. In the general case, the applied set of forces is not self-equilibrated, thus to perform static analysis some external kinematic constraints have to be imposed. Since their a priori definition is hardly possible, an approximate shape for the deformed ring will be assumed herein.

A classical strength of materials solution states that in a circular rod subjected to K equally distributed identical radial forces the contribution of the bending strains will vanish in comparison to the tensile ones as K grows [3]. In most bearings, especially used in aviation jet industry, the distribution of radial forces varies relatively smoothly and to a certain degree of accuracy the ring can be treated as a set of segments of finite length, each of which is subjected to a constant compressive load N and a bending moment M. The error introduced by this simplification is particularly low for bearings with negative internal clearance in the operating regime. The assumption is similar to replacing the general deformation scheme of the bearing semicircle depicted in Figure 2c with a pattern that does not reproduce deflections at points C and D, which are negligible in comparison to the vertical displacement at point G.

Consider a circular rod with the centerline radius $R_m$, the cross section area A and the area moment of inertia I. Being subjected to radial external load $Q_r$ its segment of length $l = \theta R_m$, shown in Figure 2b, exhibits radial deflection w and the corresponding axil strain can be computed as

$$ \varepsilon = \frac{\Delta l}{l} = \frac{\theta(R_m - w) - \theta R_m}{\theta R_m} = -\frac{w}{R_m}. \quad (8) $$

Recalling the Hook's law for one dimensional case, i.e., $\sigma = \varepsilon E$, the potential energy of this tensile deformation renders:

$$ T = \int_V \frac{1}{2} E \varepsilon^2 dV = \int_0^{\theta R_m} \frac{1}{2} E A \frac{w^2}{R_m^2} A \, ds = \frac{1}{2} E A \theta \frac{w^2}{R_m}. \quad (9) $$

Besides tension phenomenon, radial deflection w alters the segment curvature radius that engenders pure bending stress state. The corresponding deformation energy is [3]

$$ B = \int_{0}^{\theta R_m} \frac{M^2}{2EI} ds = \frac{1}{2} EI \theta \frac{w^2}{R_m(R_m - w)^2}, \quad (10) $$

<figure><img src="image_4.png" /><figcaption>(a) Common loading of the ring.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Ring segment deformation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Ring deformation patterns.</figcaption></figure>

<caption>Figure 2. Radial deformations of the bearing inner ring.</caption>

where the following well-known relation between the bending moment and the change of the curvature radius was applied:

$$ \frac{1}{R_m - w} - \frac{1}{R_m} = \frac{M}{EI}. \quad (11) $$

Comparing expressions (9) and (10) we may observe that, indeed, $B \ll T$. This fact correlates with the initial hypothesis of the prevalence of tensile strains over bending ones.

## 2.3 Roller equilibrium

To obtain the roller force vector $P_{local}$ and the tangent stiffness matrix $K_{local}$ let us describe the roller equilibrium for a given vector $y_{local}$. Equilibrium of a roller in contact is determined by the following parameters:

* $v_r$ - roller radial displacement,

* $\varphi$ - roller rotation about the local tangent axis n,

* $w_i$ - radial deflection of the inner ring central line,

* $w_e$ - radial deflection of the outer ring central line,

which are gathered in the roller state vector

$$ y_{roll} = [v_r \quad \varphi \quad w_i \quad w_e]^T. \quad (12) $$

The roller axial displacement is restrained. Such a simplification is acceptable for small misalignment angles. Moreover, this axial degree of freedom is normally physically inadmissible due to presence of the cage. The roller degree of freedom $\varphi$ and the corresponding outer ring rotation angles $\gamma_x$ and $\gamma_y$ can be removed from the model if it is not necessary to capture the effect of misalignment.

The main bearing parameters required for definition of its geometry are summarized in Figure 3a. Figure 3b displays
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 3. Initial and deformed configurations of the bearing roller and rings.</figcaption></figure>

the deformed configuration of the bearing in the roller cross section. The roller contact compressions between inner and outer rings at its central section are

$$ \delta_i = -\frac{g}{4} - v_r - w_i \quad \text{and} \quad \delta_e = -\frac{g}{4} - u_r + v_r - w_e, \quad (13) $$

where g is the actual diametral clearance. When dealing with slender bearing supports, one should not use the initial diametral clearance in (13). Mounting effects, differences in the operating temperatures of the rings and their centrifugal expansion can significantly affect value of g. These contributions can be easily accounted for, as done, for instance, by Aramaki et al. [2] for centrifugal expansion of the rings, but are beyond the scope of the present paper.

A roller may have a non-cylindrical profile. In this case its geometry is reproduced approximately by means of the slicing technique [6]. The roller profile is discretized into a finite number of so-called cylindrical slices as shown in Figure 3a. For the k-th slice one should define its radius $b_k$, width $h_k$ and distance $s_k$ between its center and the roller central section. For each k-th slice Hertzian linear contact is implied, the distributed line load per unit length is then evaluated as [1]

$$ q_{ak} = C_a \delta_{ak}^n, \quad (14) $$

where index $a = \{e, i\}$ stands for either outer or inner ring, respectively. If the k-th slice is in contact, $\delta_{ak} > 0$, otherwise $q_{ak} = 0$. The contact compressions at the k-th slice are computed from $\delta_a$ and $\varphi$:

$$ \delta_{ik} = \delta_i - R_b + b_k - \varphi s_k, \quad (15a) $$

$$ \delta_{ek} = \delta_e - R_b + b_k + (\varphi - \gamma_n)s_k. \quad (15b) $$

The contact stiffness coefficient $C_a$ depends on the material properties. If the same material with elasticity modulus E and Poisson ratio $\mu$ is used for manufacturing both rollers and rings, then the following relation holds [1]:

$$ C_a = \frac{\pi EL}{(1 - \mu^2) (7.358L)^n}, \quad (16) $$

where for linear contact the exponent number $n = \frac{10}{9}$. Provided the misalignment angle $\gamma_n$ is restrained, the proposed formulation can be directly applied for radial ball bearings analysis. One simply needs to modify contact relation (14), i.e. alter the contact coefficients $C_a$ and the exponent number $n$ (for ball bearings $n = \frac{3}{2}$), also discretization into slices is no longer required.

Recalling (14), we may write the energy of the contact pair roller-ring as follows:

$$ U_a = \frac{1}{n+1} \sum_{k=1}^{K} C_a \delta_{ak}^{n+1} h_k. \quad (17) $$

The roller is also subjected to the centrifugal inertia force, which is related with the kinetic energy as

$$ V = \frac{1}{2} m \omega_c^2 (R_p + v_r)^2, \quad (18) $$

where $m$ is the roller mass and $\omega_c$ is its rotation speed about the bearing axis. Referring to (9) and (10), we formulate the total energy functional of a roller clamped between the two rings:

$$ \Pi = U_e + T_e + B_e + U_i + T_i + B_i - V. \quad (19) $$

Minimization of (19) delivers the roller equilibrium position and determines the unknown vector $\mathbf{y}_{roll}$. Components of the roller force vector are computed as the first order derivatives of the functional with respect to the components of $\mathbf{y}_{roll}$, i.e.,

$$ \mathbf{P}_{roll} = - \frac{\partial \Pi}{\partial \mathbf{y}_{roll}} = \begin{bmatrix} Q_i - Q_e + F_c \\ M_i - M_e \\ Q_e - N_e \\ Q_i - N_i \end{bmatrix} \quad (20) $$

having

$$ Q_a = \sum_{k=1}^{K} C_a \delta_{ak}^n h_k, \quad (21a) $$
$$M_a = \sum_{k=1}^{K} C_a \delta_{ak}^n h_k s_k, \qquad (21b)$$

$$N_a = \frac{E A_a \theta w_a}{R_{am}} + \frac{E I_a \theta w_a}{(R_{am} \pm w_a)^3}, \qquad (21c)$$

$$F_c = m\omega_c^2 (R_p + v_r), \qquad (21d)$$

In (21c) “+” is used in case of the outer ring geometry ($a = e$) and “-” stands for the inner one ($a = i$). The ring geometric parameters $A_a$, $I_a$ and $R_{am}$ can be easily obtained from Figure 3a. Similarly to (20), we may find analytical expressions for the components of the tangent matrix $K_{roll} = -\frac{\partial \mathbf{P}_{roll}}{\partial \mathbf{y}_{roll}}$, however, these are omitted here due to lack of space.

The roller quasi-static equilibrium condition is written as

$$\mathbf{P}_{roll} = \mathbf{0}. \qquad (22)$$

The system of nonlinear equations (22) can be resolved by means of the Newton-Raphson method. The equilibrium position of each roller is determined independently from the others, thus the iterative procedure is applied only to a small system of 4 equations with 4 unknowns. This positively influences performance and stability of the entire algorithm.

Once the roller equilibrium position is governed, the desirable stiffness matrix $K_{roll}$ appearing in (7) can be obtained. To this end, at first, we vary the energy functional (19) not only over the components of vector $\mathbf{y}_{roll}$, but also over the components of vector $\mathbf{y}_{local}$. This results in appearance of the following extended force vector and stiffness matrix:

$$\mathbf{P}_{full} = - \begin{bmatrix} \frac{\partial \Pi}{\partial \mathbf{y}_{roll}} \\ \frac{\partial \Pi}{\partial \mathbf{y}_{local}} \end{bmatrix} = \begin{bmatrix} \mathbf{P}_{roll} \\ \mathbf{P}_{ext} \end{bmatrix}, \qquad (23)$$

$$\mathbf{K}_{full} = - \begin{bmatrix} \frac{\partial \mathbf{P}_{roll}}{\partial \mathbf{y}_{roll}} & \frac{\partial \mathbf{P}_{roll}}{\partial \mathbf{y}_{local}} \\ \frac{\partial \mathbf{P}_{ext}}{\partial \mathbf{y}_{roll}} & \frac{\partial \mathbf{P}_{ext}}{\partial \mathbf{y}_{local}} \end{bmatrix} = \begin{bmatrix} \mathbf{K}_{roll} & \mathbf{K}_{rr} \\ \mathbf{K}_{rr}^{T} & \mathbf{K}_{ext} \end{bmatrix}. \qquad (24)$$

Notation in (24) reflects symmetry of the resultant matrix (this fact is easy to verify) that directly follows from existence of the functional. Submatrices from (24) are omitted herein due to cumbersome but can be analytically derived. Since the bearing element does not have any external degrees of freedom directly related to the rollers, we perform static condensation procedure, taking the components of vector $\mathbf{y}_{local}$ as master degrees of freedom:

$$\mathbf{P}_{local} = \mathbf{P}_{ext} - \mathbf{K}_{rr}^{T} \mathbf{K}_{roll}^{-1} \mathbf{P}_{roll} = -\frac{\partial \Pi}{\partial \mathbf{y}_{local}} = \begin{bmatrix} Q_e \\ M_e \end{bmatrix}, \qquad (25)$$

$$\mathbf{K}_{local} = \mathbf{K}_{ext} - \mathbf{K}_{rr}^{T} \mathbf{K}_{roll}^{-1} \mathbf{K}_{rr}. \qquad (26)$$

When simplifying (25), we also imply that at the roller equilibrium position identity (22) is fulfilled.

## 2.4 Extension to large displacements and rotations

The obtained expressions hold for the situation when the inner ring is fixed and the outer exhibits small displacements and rotations. Consider now nodes A and B of the bearing are allowed to have arbitrarily large deflections $\mathbf{u}$ and rotations $\theta$ – the latter is the Euler vector [9]. The element state is then described by means of vector

$$\mathbf{y} = \begin{bmatrix} \mathbf{u}_A \\ \theta_A \\ \mathbf{u}_B \\ \theta_B \end{bmatrix}. \qquad (27)$$

The Euler vector is related to the rotation tensor

$$\mathbf{Q} = \mathbf{I} + h_1(\theta) \mathbf{\Theta} + h_2(\theta) \mathbf{\Theta}^2, \qquad (28)$$

where $\mathbf{I}$ – identity tensor, $\mathbf{\Theta}$ = skew ($\theta$) – skew-symmetric tensor of vector $\theta$ and scalar $\theta = \|\theta\|$. Functions

$$h_1(\theta) = \frac{\sin \theta}{\theta} \quad \text{and} \quad h_2(\theta) = \frac{1}{2} \left( \frac{\sin \theta/2}{\theta/2} \right)^2 \qquad (29)$$

were also introduced in (28). Therefore, $\theta$ = axial ($\mathbf{\Theta}$) is the axial vector of tensor $\mathbf{\Theta}$.

The relative displacement of the outer ring in the coordinate system attached to node A can be evaluated as

$$\xi = [u_x \quad u_y \quad u_z]^T = \mathbf{Q}_A^T \mathbf{u}_{BA} = \mathbf{Q}_A^T (\mathbf{u}_B - \mathbf{u}_A), \qquad (30)$$

whose variation after some vectorial math renders

$$\delta \xi = \delta (\mathbf{Q}_A^T \mathbf{u}_{BA}) = \mathbf{X} \mathbf{\Gamma}_A^T \delta \theta_A + \mathbf{Q}_A^T \delta \mathbf{u}_B - \mathbf{Q}_A^T \delta \mathbf{u}_A, \qquad (31)$$

where $\mathbf{X}$ = skew ($\xi$) and the well-known tangent tensor [9]

$$\mathbf{\Gamma} = \mathbf{I} + h_2(\theta)\mathbf{\Theta} + h_3(\theta)\mathbf{\Theta}^2 \qquad (32)$$

with

$$h_3(\theta) = \frac{1 - h_1(\theta)}{\theta^2} \qquad (33)$$

appears. Operator $\mathbf{\Gamma}$ is typical for description of rotations with the Euler vector. When deriving the first term of (31) we used the following expression [11], valid for any constant vector $\mathbf{t}$:

$$\begin{aligned} \delta \mathbf{Q}^T \mathbf{t} &= -\mathbf{Q}^T \text{skew} (\mathbf{\Gamma} \delta \theta) \mathbf{t} = \mathbf{Q}^T (\mathbf{t} \times (\mathbf{\Gamma} \delta \theta)) \\ &= \mathbf{Q}^T \mathbf{t} \times (\mathbf{Q}^T \mathbf{\Gamma} \delta \theta) = \mathbf{T}' \mathbf{\Gamma}^T \delta \theta, \end{aligned} \qquad (34)$$

where $\mathbf{t}' = \mathbf{Q}^T \mathbf{t}$ and $\mathbf{T}' = \text{skew} (\mathbf{t}')$. The useful equivalence [9]

$$\mathbf{Q}^T \mathbf{\Gamma} = \mathbf{\Gamma}^T \qquad (35)$$

was also applied. Similarly to (30), the relative rotation of the inner ring in the coordinate system attached to node A can be determined as

$$\gamma = [\gamma_x \quad \gamma_y \quad \gamma_z]^T = \mathbf{Q}_A^T (\mathbf{e}_{zA} \times \mathbf{e}_{zB}), \qquad (36)$$
<figure><img src="image_9.png" /><figcaption>Figure 4. The rolling element and the neighbouring beams elements of the two rings</figcaption></figure>

where unitary axial vectors of each bearing node at the de-
formed configuration are given by

$$
\mathbf{e}_{zA} = \mathbf{Q}_A \mathbf{e}_{z0} \quad \text{and} \quad \mathbf{e}_{zB} = \mathbf{Q}_B \mathbf{e}_{z0}. \tag{37}
$$

Variation of (36) can be obtained with the aid of (34) as

$$
\begin{equation}
\begin{split}
\delta\boldsymbol{\gamma} = {}& G\boldsymbol{\Gamma}_A^T \delta\theta_A + \boldsymbol{Q}_A^T \boldsymbol{E}_{zB} \boldsymbol{E}_{zA} \boldsymbol{\Gamma}_A \delta\theta_A \\
& - \boldsymbol{Q}_A^T \boldsymbol{E}_{zA} \boldsymbol{E}_{zB} \boldsymbol{\Gamma}_B \delta\theta_B
\end{split}
\tag{38}
\end{equation}
$$

with **G** = skew (γ), **E**<sub>zA</sub> = skew (e<sub>zA</sub>) and **E**<sub>zB</sub> = skew (e<sub>zB</sub>).
Recalling definitions (30) and (36) we can now rewrite (1) as

$$
\mathbf{y}_{rel} = \begin{bmatrix} \xi \\ \gamma \end{bmatrix}, \tag{39}
$$

whose variation after substitution of (31) and (38) can be
formulated in a closed matrix form

$$
\delta \mathbf{y}_{rel} = \mathbf{H} \delta \mathbf{y}, \quad (40)
$$

where the following operator was defined:

$$
\mathbf{H} = \begin{bmatrix} -\mathbf{Q}_A^T & \mathbf{O} & \mathbf{Q}_A^T & \mathbf{X}\mathbf{\Gamma}_A^T \\ \mathbf{O} & \mathbf{G}\mathbf{\Gamma}_A^T + \mathbf{Q}_A^T\mathbf{E}_{zB}\mathbf{E}_{zA}\mathbf{\Gamma}_A & \mathbf{O} & -\mathbf{Q}_A^T\mathbf{E}_{zA}\mathbf{E}_{zB}\mathbf{\Gamma}_B \end{bmatrix} . \quad (41)
$$

Matrix **H**, by definition, allows the developed approach to be
used in multibody dynamics. When the bearing experiences
large motions, operator (41) should be applied to the force
vector and the stiffness matrix of the element:

$$
\mathbf{P}_{\text{tang}} = \mathbf{H}^T \sum_{j=1}^{N} \mathbf{S}_j^T \mathbf{P}_{\text{local } j}, \qquad (42)
$$

$$
\mathbf{K}_{\text{tang}} = \mathbf{H}^T \sum_{j=1}^{N} \mathbf{S}_j^T \mathbf{K}_{\text{local } j} \mathbf{S}_j \mathbf{H}.
$$

The presented multibody dynamics extension possess
some innate limitations. Thus, among the inertia loads, aris-
ing in the bearing, only the centrifugal force acting onto the
roller as a resultant of its relative rotation about the shaft
axis is considered. Coriolis loads and other possible inertia
effects are presently omitted.

2.5 Generalized structural FE model with Dynamics R4

Though the simplified deformation pattern, proposed in Section
2.2, allows to incorporate bearing rings flexibility without

practically any computational overhead, it may provide poorer
results for load cases, when a bearing operates with large
positive clearance and the distribution of radial forces is not
smooth. To provide better assessment of rings deformations
for such cases, a lightweight FE model was built in Dynamics
R4 software [14]. For simplicity, we prohibit misalignment,
but it can be easily added adopting kinematics, described
in Section 2.2. The use of such FE assembly in multibody
dynamics is less straightforward, but still possible.

We discretize the rings with geometrically nonlinear
beam elements. At the location of the *j*-th roller a speci-
ally developed 2-nodes finite element is placed. Basically, the
element uses the same potential (19) to describe the behavi-
our of the roller, though the radial deflections of the rings
central lines are the external degrees of freedom (not the
internal any longer) and related to the spatial displacements
of nodes *G* and *H* as follows:

$$
w_i = \mathbf{u}_G \cdot \mathbf{e}_{r0} \quad \text{and} \quad w_e = \mathbf{u}_H \cdot \mathbf{e}_{r0}, \qquad (43)
$$

The roller radial position $v_r$ is the internal degree of freedom and the inner iterative procedure is required to determine it for given $w_e$ and $w_i$. The contact compressions still can be evaluated with the aid of (13), but one should set $u_r = 0$. The local force vector of such element is written as

$$
\mathbf{P} = \begin{bmatrix} -Q_i \\ Q_e \end{bmatrix}, \qquad (44)
$$

whose components are computed by (21a).

The external load is distributed along the bottom part
of the inner ring. To remove the rigid body motion of the
structure, small stiffness coefficients K<sub>xx</sub> and K<sub>yy</sub>, which
are several orders of magnitude lower than the bending and
tensile stiffnesses of the rings, are added at each node of the
outer ring.

The resultant model is more sophisticated than the one,
presented previously, but should provide better results for
bearings operating with positive clearances.

3. NUMERICAL EXAMPLE

3.1 Static analysis and comparison with FEM

Consider the radial bearing with 24 cylindrical rollers with dimensions R<sub>p</sub> = 75 mm, R<sub>b</sub> = 6.5 mm, L = 16 mm, R<sub>s</sub> = 54 mm, R<sub>h</sub> = 97 mm, b = 30 mm. Radial load F<sub>x</sub> = 5000 N is applied to the inner ring, which rotates at a speed ω<sub>rel</sub> = 12000 rpm. Material elasticity modulus E = 2 · 10<sup>11</sup> MPa, Poisson ratio μ = 0.3 and density ρ = 7850 kg/m<sup>3</sup>.

At first we analyze the bearing with negative clearance
g = -40 µm<fnref n="1" />. If the bearing rings are treated as rigid, the
obtained values of contact radial forces turn out to be 2.5
higher than for the case when structural deformations are
permitted, as exposed in Figure 5a (only the inner ring for-
ces Qᵢ are given to keep the figure readable). In the case of

<footnote n="1">¹To guarantee optimum bearing life, its internal clearance should remain as low as possible or slightly negative at the operating regime. Due to design errors this requirement may be violated that seriously compromises bearing durability.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Distribution of the bearing internal radial forces $Q_i$ for different clearance values.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Model of the bearing with deformable rings in NX Nastran.</figcaption></figure>

positive clearance $g = 10 \ \mu\text{m}$ (see Figure 5b) the influence of rings flexibility on the forces distribution is less crucial, but still apparent. The results obtained by the structural FE model, described in Section 2.5, are also shown in Figure 5 as a reference. For $g = 10 \ \mu\text{m}$ we additionally verified the plots using formulation of Cavallaro et al. [4], which is applicable for bearings with positive clearances. According to these reference curves, the proposed approach based on the considered simplified deformation pattern, though still captures the radial forces reduction, performs worse if the bearing operates with positive clearance (as anticipated in Section 2.5). However, in reality the shaft and the housing are not just hollow cylinders and normally reinforced with some adjacent elements, whose contribution to the resultant stiffness can be hardly estimated analytically. With this additional stiffening it is apparent that the reduction of the radial forces acting on the rolling elements will not be as significant as forecasted by the structural or the Cavallaro model. From this point of view, the results obtained using the simplified deformation pattern hypothesis can be treated as a relatively good first order approximation, particularly useful for design computations. Obviously, it should be further verified by alternative, more accurate, approaches. Parallel use of the

both numerical models, presented in Section 2, leads to a very efficient computational strategy.

The distribution of the radial forces obtained for $g = -40 \ \mu\text{m}$ was also verified by FEM analysis performed with NX Nastran<fnref n="2" />. To this end a 3D model from Figure 6 was used. Surface traction $q_x(\varphi) = q_0 b \cos \varphi$ was distributed along the lower part of the inner ring raceway such that $F_x = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} q_0 b \cos \varphi \, d\varphi$. Kinematic boundary conditions of the outer ring are difficult to determine, so the latter is attached from both ends to long thin cylinders with length 235 mm. The opposite sides of cylinders are clamped. In order to vanish the cylinders contribution to rigidity of the assembly, their elasticity modulus is set $E' = 0,01 E$. The FEM model of the bearing with rigid rings is a lot simpler. In this case we fix the external surface of the outer ring and apply load $F_x$ pointwisely at the inner ring central node, which is connected to the inner raceway surface with rigid elements. For each roller we also set $u_z = 0$ at its gravity center to prevent axial displacements.

<footnote n="2">²FEM model of the bearing with positive clearance (for both flexible and rigid rings cases) was not built in NX Nastran within the present work due to its substantive complexity. Discussion of various aspects of such finite element problem itself can become a topic of a separate article.</footnote>
## 3.2 Multibody dynamics simulation

To demonstrate the applicability of the proposed approach for multibody dynamics problems, the model from Figure 7a was analyzed. This assembly was designed in one of our research projects during the development of a 3D model of a planetary gearbox. The inner and outer rings of the bearing, whose parameters are given in Section 3.1, are located at nodes A and B, respectively. The outer ring is connected to the basement C via the spatial rigid link element [9] of length $l = 1$ m. This FE uses Lagrange multipliers to wrap the constraints

$$ 
\begin{aligned}
\mathbf{u}_{BC} &= \mathbf{u}_B - \mathbf{u}_C - (\mathbf{Q}_C \mathbf{e}_{l0} - \mathbf{e}_{l0}) \quad l = \mathbf{o}, \\
\theta_{BC} &= -\frac{1}{2} \text{axial} (\mathbf{Q}_B \mathbf{Q}_C^T) = \mathbf{o}
\end{aligned}
\quad (45) $$

into the functional

$$ \Phi = \lambda_u \cdot \mathbf{u}_{BC} + \lambda_{\theta} \cdot \theta_{BC}, \quad (46) $$

whose minimization, similarly to (23) and (24), leads to the corresponding force vector and tangent matrix. Physically, the unknowns $\lambda_u$ and $\lambda_{\theta}$ are nothing but the reaction forces and moments transmitted through the connecting rod.

All degrees of freedom of node C, except the axial rotation $u_{rz}$, are prohibited. The pointwise mass with $m_B = 2$ kg is located at node A. The basement node C rotates at an axial speed $\omega_C = 50$ rds $\approx 477$ rpm. This value is specially chosen to provide the resultant centrifugal load, acting on the mass at node A, identical to the value of the statically applied radial force from the first test, i.e., $F = m_B l \omega_C^2 = 5000$ N. Independent rotation is also defined at node A. Its value $\omega_A = 12000 + \omega_C = 12477$ rpm is adjusted such that the relative rotation speed of the inner ring would coincide with $\omega_{rel}$, used in the aforementioned statics test. The internal negative clearance $g = -40$ µm is again introduced in the bearing and the rings structural deformations are permitted. Both rotations $\omega_A$ and $\omega_C$ are not applied instantly, but were set linearly increasing as shown in Figure 7b. To advance the solution the implicit Newmark integration scheme [8] is invoked.

Under these circumstances the radial forces arising in the bearing are expected to agree with the ones obtained previously in the quasi-static test. Evolution of the radial forces acting on an arbitrary roller is plotted in Figure 7c. Indeed, the graph for $Q_i$ perfectly correlates with the previously obtained curves: the peaks of the graph coincide with maximum and minimum values of the force distributions from Figure 5a.

## 3.3 Bearing misalignment

To assess the model performance in the case of misalignment, the bearing from Section 3.1 with negative clearance is loaded additionally with constant moment $M_y = 13.7$ Nm. The roller was discretized with 10 identical slices, the same mesh density along the roller length was used in NX Nastran to build contact regions with the underlying solid elements. The applied moment does not change the distribution of the radial

<figure><img src="image_12.png" /><figcaption>Figure 7. Transient analysis of the bearing with flexible rings and negative clearance $g = -40$ µm rotating about the fixed basement.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 8. Distribution of contact forces along the bearing roller in the case of misalignment, $g = -40$ µm.</figcaption></figure>
forces in the bearings (as expected), but rather influences the distribution of the contact pressure along the roller. Such behaviour was also confirmed by the FEM model, where a pair of forces with $||F_z|| = 100$ N was applied at points $R_1$ and $R_2$, as shown in Figure 6, to reproduce the moment load. Figure 8 depicts distribution of the contact pressure along the most loaded roller number 13. The results are in good agreement with NX Nastran. Note the drop in the loads at the ends of the rollers in the 3D FEM simulation: this does not occur in our model as the Hertz line solution has been applied to each slice, distinctly from NX Nastran numerical technique, which solves the non-Hertzian contact problem for finite-length bodies instead. Despite of this the resultant contact forces for both solutions are very close, see Figure 5a.

## 4. CONCLUSION

A model of the radial roller bearing that accounts for structural deformations of its rings was presented. Numerical experiments showed that their contribution is particularly important for those situations when a negative clearance appears at the operating mode, while for regimes with positive clearance this behaviour can be neglected unless the applied loads are extreme (e.g., blade loss situation [15]). The equilibrium position of each rolling element is determined with a set of 4 variables: a radial displacement and an in-plane rotation of the roller and radial deflections of the raceways of each of the two rings. Roller inertia properties are also considered. The model can be implemented in a 2-nodes finite element with 12 degrees of freedom and then integrated practically into any FEM code as it was done for Dynamics R4.

Relative simplicity, high performance and robustness are the most apparent advantages features of the proposed formulation. The internal iterative process that determines the roller equilibrium is computationally efficient due to rolling elements uncoupling and normally converges in 3-4 iterations.

The model currently prohibits axial load to be applied. Addition of flanges and friction into the roller equilibrium along with introduction of axial degree of freedom for the rolling element could alleviate this issue. It is also necessary to capture more complex deformation patterns of rings, especially for configurations when the bending factor becomes essential relatively to the tensile strains. In this case determination of equilibrium for each rolling element can no longer be made independently and the bearing should be handled as a complete assembly of the two rings and the rollers with all internal degrees of freedom influencing each other. A concept of such a model is also proposed in the paper and its further development is in the scope of the authors' future work.

## REFERENCES

[1] S. Andréason. "Load distribution in a taper roller bearing arrangement considering misalignment". In: *Tribology* 6.3 (1973), pp. 84-92.

[2] H. Aramaki, Y. Shoda, Y. Morishita, and T. Sawamoto. "The Performance of Ball Bearings With Silicon Nitride Ceramic Balls in High Speed Spindles for Machine Tools". In: *Journal of Tribology* 110.4 (Oct. 1988), pp. 693-698.

[3] R. Budynas, W. Young, and A. Sadegh. *Roark's Formulas for Stress and Strain, 8th Edition*. McGraw-Hill Education, 2011.

[4] G. Cavallaro, D. Nelias, and F. Bon. "Analysis of High-Speed Intershaft Cylindrical Roller Bearing with Flexible Rings". In: *Tribology Transactions* 48.2 (2005), pp. 154-164.

[5] W. Cheng. "Experimental and Numerical Study of Multibody Contact System with Roller Bearing-Part II: Semi-Finite Element Analysis and Optimal Design of Housing". In: *Tribology Transactions* 39.1 (1996), pp. 166-172.

[6] J. M. de Mul, J. M. Vree, and D. A. Maas. "Equilibrium and Associated Load Distribution in Ball and Roller Bearings Loaded in Five Degrees of Freedom While Neglecting Friction-Part II: Application to Roller Bearings and Experimental Verification". In: *Journal of Tribology* 111.1 (1989), pp. 149-154.

[7] E. G. Filetti and J. H. Rumbarger. "A General Method for Predicting the Influence of Structural Support Upon Rolling Element Bearing Performance". In: *Journal of Lubrication Technology* 92.1 (Jan. 1970), pp. 121-127.

[8] B. Gavrea, D. Negrut, and F. A. Potra. *The Newmark Integration Method for Simulation of Multibody Systems: Analytical Considerations*. 2005.

[9] M. Géradin and A. Cardona. *Flexible multibody dynamics: a finite element approach*. Chichester, England: John Wiley, 2001.

[10] T. A. Harris. *Rolling Bearing Analysis, 4th Edition*. 4th ed. 2000.

[11] V. Ivannikov, C. Tiago, and P. M. Pimenta. "On the boundary conditions of the geometrically nonlinear Kirchhoff-Love shell theory". In: *International Journal of Solids and Structures* 51.18 (2014), pp. 3101-3112.

[12] A. B. Jones and T. A. Harris. "Analysis of a Rolling-Element Idler Gear Bearing Having a Deformable Outer-Race Structure". In: *Journal of Basic Engineering* 85.2 (June 1963), pp. 273-278.

[13] A. Leblanc, D. Nelias, and C. Defaye. "Nonlinear dynamic analysis of cylindrical roller bearing with flexible rings". In: *Journal of Sound and Vibration* 325.1-2 (2009), pp. 145-160.

[14] Alfa-Tranzit Co. Ltd. *Dynamics R4*. http://alfatran.com/. 2005-2017.

[15] Leontiev M. K., Davydov A. V., Degtyarev S. A., and Gladkiy I. L. "To simulation of fan blade out for a high bypass ratio engine". In: *Russian Aeronautics* 57.2 (2014), pp. 154-161.