FEM/BEM TECHNIQUES FOR MODELLING OF LOCAL FIELDS IN CONTACT MECHANICS FEM/BEM TECHNIQUES FOR MODELLING OF LOCAL FIELDS IN CONTACT MECHANICS

parts of structure by the use of small number of large FE (finite elements) by effective This contribution contains description of modelling technique of contact problems of bodies with curved surfaces when assumption of infinitesimal displacements can be considered to give sufficient accuracy for both displacement and stress analysis. This assumption considers that the local configuration of the bodies in contact will not be influenced during the deformation process and only dimensions and shape of the contact surfaces and contact pressures will change by the load conditions. Such assumptions enable to reduce the contact problem conside-rably, if the local contact displacement and stress fields are modelled as local fields inside large elements (sub-domains) using superposition of a local Hertz type field and a smooth field modelled in a classical way using large FEM or BEM technique. A technique of obtaining the complex solution of the bodies in contact as a combination of the local contact field defined by the Boussinesq’s solution and the smooth field modelled by multi-domain BEM is described and considered to be a basis for the modelling and design of contact problems containing large gradient displacements.


Introduction
If spherical bodies in contact have much different curvatures at least in one direction, like it is in roller bearings, the contact area is very small and under normal loading the contact pressure may occur in order of several thousand N/ mm 2 . The contact area can be small in both contact surface directions (called point contact) or in one direction only, when the difference of the curvatures in the other direction is small (called line contact).
The distribution of the externally applied load over the individual rolling elements and also the distribution of the contact pressure in the element contact area are important parameters, which influence the static and dynamic bearing capacity and life of the structural parts in contact. The elastic deformation of the adjacent components and the local deformation of the rolling elements and raceways at their contact area are decisive parameters for the distribution of the contact pressure. The very localised deformations and pressures in the flexible structure components make the problem high non-linear and the accurate analysis is necessary in order to obtain reliable results for the design of the whole bearing form and arrangements.
The detail distribution of the load over the individual rolling elements and the pressure distribution in the contact area of each element require very fine element mesh [1], high performance computers and effective software for the numerical modelling of the problem. The pressure distribution between the elements and the raceways follows the Hertz distribution approximation in the elliptic form only at the lower load levels.
Many authors [2 -7] gave the analytic expressions for calculation of stress and displacement fields for contact problems of bodies of revolution using Hertz solution obtained from a half plane for 2D and for a half space for 3D problems. The analytic improvement for inclusion of frictional or tractive loads (Smith-Liu equations) was given in [8].
Several authors (e. g. [9 -11] used the Hertz solution, or more general Boussinesq solution (giving the displacements and stresses in a half space loaded by a concentrated force in the surface plane) to solve the contact problem of two bodies of rotation more generally.
The aim of the paper is to define an effective model for accurate enough stress and deformation analysis of the multibody contact of the spherical parts of structure by the use of small number of large FE (finite elements) by effective

2D problems
All real contacts are three-dimensional (3D), and therefore theoretically demand a solution in the 3D theory of elasticity. But few solutions exist and where it is feasible to approximate the geometry to one or two dimensions, it may be that an exact solution thus facilitated is preferable to an approximate numerical solution to the exact problem. In this chapter, we shall be looking at plane problems. A simple, practical problem is the contact of two rollers of the same length. It is clear that the central portion is under plane strain, whilst ends are under plane stress. A second approximation is that the contacting bodies may be considered to be half-planes. We will first examine the displacement and stress field in the semi-infinite plane (Boussinesq problem) subjected to normal force ( Fig. 1).
modelling of both the local displacement and stress fields. In such technique the global field is modelled in a classical way and the local effects containing large field gradients are defined by the Boussinesq solution, which satisfy the governing equations in all points of corresponding sub-domain and the boundary conditions at least in the local sense.
In the multi-body contact problems the deformation characteristics are computed first for some local sub-domain and the load distribution over the individual elements is modelled using the global model using special roller elements containing the non-linear contact stiffness of both the roller element and the raceways including the gap for corresponding rolling element [12].
The infinitesimal displacements are considered in the model. Under these assumptions, both the contact area and the contact pressure are to be determined from the topology of the undeformed bodies in the contact and from the resulting load conditions of each pair of contacted surfaces. However, the load distribution between the contacting bodies is considered to be result of the global analysis of all contacting bodies taking into account the non-linear contact stiffness, the finite but small rotations and bending and displacement of shafts, housing, etc. obtained from the global analysis.
The stiffness characteristic of the contacting pairs is obtained from a local model and it is defined as 1D contact pair if the contact surface is determined by the curvatures of contacted surfaces and it relates the contact force to normal displacements of contacted bodies. On the other side, it is defined as 2D contact if the position of the resulting force changes along one of the minor axes depending on the mutual rotation of contacting bodies and/or, if the curvature of contact surfaces changes along this axis. However, in each case, it is assumed that the topology of each contact surface does not change with the load and the contact element for the global model is considered to be 1D with the stiffness obtained from the 2D characteristic containing the dependence of the resulting force direction. In this way the contact characteristics can be obtained from the local model and the complex multibody contact can be modelled by large elements, which leads to considerable saving of computational time and computer memory.
We do not take into account the frictional forces in this paper, but their inclusion is simple and straightforward as it is known from the referenced literature. The methodology used in this approach is not restricted to the contact problems of the bodies with curved surfaces, but can be used in some modifications to other problems like cracks inside the bodies or on their surface, small inclusions and inhomogenities, etc., i.e. all problems, in which the local effects cause large gradients in the displacement and stress fields.
where R is radius of the disc and λ and μ are Lame constants for the material of the disc. Choosing the Hertz distribution of the contact force acting in the disc, the stress and displacement fields can be computed in a similar way (e.g. by integration of singular integrals) as given above for the half-plane.
The analysis of the stress state helps identify the critical maximum stresses in the contact region. They reveal vital information concerning the surface fatigue mechanism, plastic deformations and bearing capacity of the structure. Maximum shear stress which is in some depth below the surface and large stress gradients in the whole contact surface and sub-surface area require very fine mesh when they are to be evaluated with If we consider the stress state induced within a body by an arbitrary pressure p(ξ), then the stresses are given by x p d 2 where the integrals are evaluated over the loaded area. The sign of pressure is defined in accordance with p(x 1 ) = σ 22 (x 1 ) and hence contact pressures will always be negative.
Similarly, expressions for the displacements can be found. The contact problem of two cylinders pressed together with their axis parallel is a two dimensional equivalent of Hertzian contact. The following considerations apply: 1. The load is sufficiently small for the contact patch, of width 2a, to be small in comparison with both R 1 and R 2 , which are the radii of curvature of the cylinders. This means that the assumption that the contacted bodies can be approximated by half-planes remains valid. 2. No surface shear tractions must arise, that is, both bodies have the same elastic constants, or the coefficient of friction must vanish.
This enables the relative normal approach of particles within the contact to be approximated by a parabola and gives an elliptic contact pressure distribution where p 0 is its peak value in the middle of the contact patch. The integrals are singular with the singularity in the surface point x 2 = 0, x 1 = ξ. So, appropriate technique has to be used for their evaluation as it is known from the Boundary Integral Equation Method (Boundary Element Method) [13]. Muskhelishvili [14] gives a solution for a circular disc loaded by two equilibrated forces (Fig. 2). ; E E (9) Similarly, the stresses for a concentrated tangential force [15], , , The total stress and displacement fields for distributed load over the elliptical area on the half space boundaries are obtainable by integration over the area in two ways: using the cylindrical coordinates r, φ and x 3 , or, using the cartesian coordinates.
If, in the first case, the cartesian coordinates are chosen so that the main axes of the contact ellipse will be identical with the axes x 1 and x 2 and the half axes of the ellipse are a and b, respectively. Then r and φ are defined as The numerical integration is performed using the Gauss quadrature rule for the r coordinate and by the equidistributed quadrature points in the circumferential direction, φ.
In the second case, the Gaussian quadrature rule is used for both x 1 and x 2 coordinate directions, using constant Gauss coordinate along the longer axis and variable integration limits in the second direction.
No equivalent closed solutions are available for 3D rolling elements like that for 2D cylinder given in the previous section. required accuracy using classical numerical models like FEM or BEM.
Accuracy of the displacement analysis is responsible for the correct computation of load distribution on the rolling elements in the multi-body contact problems like rolling bearings.

3D problems
We will suppose infinitesimal elastic deformation of the bodies in contact. It means that contact pressure will be determined from mutual position of the contact surfaces assumed to be rigid in their local position. However, the local position of the contact surfaces is supposed to be determined by the total deformation of the complex system of elastic bodies, i.e. by bending of shafts, elastic deformation of housing, etc. The pressure distribution is then defined by the form of penetration of the surfaces and by the total force acting in the corresponding pair of contacting bodies. It means that the problem is non-linear and has to be solved iteratively.
We will call the 3D contact to be a point contact if the curvature of the contact surfaces does not change in the contact region. In such case the location of the contact is determined from the total deformation of the contact bodies as mentioned above and the contact surface is an ellipse. On the other side, if the contact surfaces have more general form, then we will call the contact to be a line contact. In such case the pressure distribution in the circumferential direction is defined similarly as it was in the 2D case in eq. (5), but the pressure function in longitudinal direction is not a unique function and it will depend on both geometry and location of resulting contact force for each pair of contacting bodies.
Similarly as in 2D, the basis for the numerical computation of the local contact stresses and displacements is the Boussinesq problem giving the response to a force acting in a point of the boundary of the half space.
The stresses for the problem of a concentrated normal force, P 3 , acting at the origin on the half space, x 3 ≥ 0, are given [15]  = + (16) and the last part of equation (15) with summation over all elements. Setting for tractions from (16) into (17), the following resulting system of discretised equations is obtained

M Tu p MU f f U
or shortly KU p = .
Stresses of increased accuracy and smooth over the hole domain can be obtained from the nodal displacements and traction boundary conditions (for the points near or on the domain boundaries with prescribed static conditions) using T-functions as interpolators (see [19] for details). The convergence rate of both displacement and stresses is 4 th order, when quadratic boundary elements are used in the formulation.
As the contact surfaces have mostly rotational shape, special sub-domains are used to model this kind of problems. The whole integration for obtaining the sub-domain matrices is performed over corresponding sub-domain boundaries (boundary elements) and so, the contact surfaces are easy to design. The contact sub-domain shapes are ring segments for the raceways, or full circles for the cylindrical roller elements for 2D problems and 3D ring segments and full 3D rolling elements for complex 3D problems.
The procedure of obtaining the stiffness characteristics of rolling elements in contact is an inverse procedure defined as follows: A. If the contact is 2D, or 3D with constant curvatures: choose both, the contact area and the pressure distribution according to the Hertz theory on the rolling element and raceways. As the problem is linear, the resulting contact force corresponding to the contact area is obtained from the penetration of both pairs of contact areas in the unloaded (initial) state and from the pressure and corresponding deformation of the contact surfaces. In this way the deformation characteristic is computed for a series of sizes of contact area and corresponding resulting contact force. The deformation characteristic gives the dependence of resulting contact force to the total deformation in direction of the resulting force.
B. If the contact is 3D and the curvature changes in the longitudinal direction of the rolling elements (it is often

Sub-domain formulation for the contact region
Very efficient formulation for modelling of displacement and stress fields in the contact regions can be obtained using Trefftz type reciprocity based elements. Here, the principles well known from the basic BEM formulations (see e.g. [16 and 17]) are used and so, we will give only the basic principles. For more details of the formulations used in the models see [18].
Weak formulation of equilibrium equations can be written using the principle of weighted residuals as where σ ij and b i are the stress tensor and body force vector, respectively. We can choose the weight functions u i * to be Trefftz (T-) displacement fields. The T-functions are the functions which satisfy all governing equations inside the domain and they have to satisfy the homogeneous equilibrium equations (with body force absent) of the same body in this case.
Applying both the integration per parts and the Gauss theorem to this equation two times, the following relation is obtained where t i are tractions acting on the domain boundaries Γ of the domain Ω. This equation expresses the reciprocity of works done by two systems of forces, the one without stars which is looked for, and another, reference state (for which all, displacements, stresses and tractions are known inside and on the domain boundaries), denoted by letters with stars. The reciprocity equations are well known from the BEM formulations where the Kelvin fundamental solutions are used for reference state and lead to the singular integral equations. Using T-polynomials, or Kelvin solutions with the singular points located outside the domain, all integrals are regular. However, for more complex problems, the multidomain formulation is advantageous because it leads to sparse matrices, as it is known from FEM. The displacements between the sub-domains are chosen to be compatible and the inter-domain equilibrium is satisfied in a weak, integral sense as follows where Γ i , Γ t a Γ e are the inter-domain boundaries, the boundaries with prescribed tractions and element (subdomain) boundaries respectively and the letters with a bar denote the fields with prescribed value. B. Compute the displacements in points y (radius R i ) inside the domain from the boundary displacements and tractions by equation (25) , , , .
C. Define patch displacement fields from the nodal displacements in a node on the boundary (the point of interest) and in two neighbour nodal points on the boundary and one inside the domain and from the prescribed boundary tractions in the point of interest using (POI) second order T-polynomials with local origin of the co-ordinates in the POI. For such field it is simple to find corresponding stresses in the POI. Similarly we can proceed, when the stresses in points near the boundary are to be obtained, i.e. using the patches of nodal points which are closest to the POI.
When the displacement and stress field contain large gradients like Hertz contact fields, then the resulting field is split into the local part embodying the large gradients and the smooth part of corresponding field as , , where the indices R, L and S denote the resulting, local and smooth part of corresponding field, respectively. The local field is described in an auxiliary region which can only partially coincide with the real region (Fig. 4). The auxiliary region contains the local field solution, whereas the real region is solved numerically with prescribed smooth boundary conditions and thus, it does not require very fine discretisation.
the shape obtained by the optimised design), or if some parts of both the raceway and the rolling elements are cylindrical and the resulting load is asymmetric, then we can proceed in similar way as in previous case, but we have to choose some asymmetric penetration positions (crossing longitudinal axes of rolling element and raceway) and to obtain the characteristics for several such positions and to interpolate between them for the specific conditions. The characteristics give dependence of four parameters: (1) the resulting force, (2) its position, (3) the total deformation and (4) the angle between the axes. The rolling elements have the shape, which is advantageous to be modelled by BEM (small surface to volume ratio). The single domain formulation using T-functions is obtained directly using equation (14) for this purpose with Boussinesq's or Kelvin's functions with singularities located outside the domain for the weighting functions defining the reciprocal states. As we need to compute the stresses on and near the boundaries of such domains a very suitable procedure is defined as follows: A. Solve the equation (14) in the form (20) for unknown displacements from prescribed tractions in the boundary nodes (located in points on the radius R as shown in Fig. 3 for 2D problem) using the Kelvin's functions with source points y located outside the domain (on radius R o ): where , U y x r r r 16 1  The smooth fields are computed by the models described above, i.e. by the equation (19) for the MD formulation, or by the equation (20) for the single domain solution. When the intensity of the local field is known, then the problem can be split into local and smooth part according to (26) and solved separately for each of them, whereas in the case when its intensity is unknown (e.g. cracks, or other geometrical concentrators), then both the local displacements and traction boundary conditions have to be included into the formulation. As in the last case the local solution satisfies the governing equations in strong sense, the accuracy of the solution will not be destroyed by the large gradients and the fineness of discretisation is dictated by the smooth field part of the solution.

Some numerical results and discussion
The local 2D fields are studied on a cylindrical rolling element (Fig. 5). As we can see from Fig. 5, the smooth displacement field contributes to the total deformation of the cylinder by about the same amount as the local displacements. Note that the rolling contact element in the global model consists of the total displacement of the roller plus the local deformation of the raceways. The local deformation of the raceway is to understand the difference of its total displacements and the smooth field solution [18]. Figures  5e and 5f show the isolines of maximum shear stress which are responsible for plastic deformation and damage arising in some depth below the contact surface. Figure 5f contains the extreme region of this field and it is evident that the models give good results even in the close vicinity of integrand singularities.
The integrals for the evaluation of the displacement and stress fields are singular, but Gauss quadrature formulas give good accuracy in all important points with except of the points on the contact surface where stresses are given by the contact pressure distribution. Generally, the stresses in this area are computed by the technique (T-polynomial interpolation from both the nodal displacements and the boundary tractions) described in the previous section.
The raceways are modelled by the complete ring using similar technique as that for cylinders (equidistant integration points in circumferential direction), or by a ring segment both 2D and 3D problems. Examples of such elements deformed by pure bending are shown in Fig. 6 (2 x 5 quadratic boundary elements [BE] in circumferential direction and 2 x 1 BE in radial direction) and Fig. 7 (4 x 5 quadratic BE in circumferential direction and 2 x 1 BE in radial direction). Each of them defines one sub-domain. The pure bending states gave accurate displacements and stresses (with exception of truncation errors) and served for testing the accuracy of the formulation. Also more general forms can be used for the investigation of the local effects, but the geometry of the raceways should be defined as good as possible, e.g. it has to be modelled by the rotational surfaces.
The elliptical Hertzian contact can be solved using the potential methods of Boussinesq and Cerutti, fully described problem of obtaining the stiffness characteristics of the rolling elements is defined as an inverse problem, i.e. each body in the contact is modelled separately, which further simplifies the solution and increases the efficiency of computational models. The technique can be extended to more flexible bodies in contact when the curvature of the contact surfaces is taken into account by an iterative process. Also in such problems considerable improvements can be achieved in the rate of convergence and computer time and memory requirements reductions, because the changes are mostly not large and so small number of iteration steps can be expected to achieve the convergence, unlike it is by the classical numerical methods where many parameters (nodes) are contained in the contact area with the contact conditions changing in many of them during the iteration process. Further improvements of the numerical treatment of the general 3D contact, as well as the use of presented method to other field problems are the aims of our current research.
Similar technique of the splitting the field containing large gradients into the local field with large gradient and another, smooth enough field, can be used also for other types of problems, e.g. problems with cracks, material inhomogenities, inclusions, etc.
by Love [20] and the complete solution for normal as well as shear tractions is introduced in Hills et al. [5]. Classical numerical way of integration by Gauss quadratures gives good results for some distance from the contact area, but not in very close vicinity of the contact surface.

Conclusions
We have shown that the contact of the bodies with cylindrical surfaces can be very effectively solved when the local fields are split out of the solution and the local problem is modelled using the Boussinesq half space solution. The other field, which is smooth enough, is then modelled by large elements. Special T-elements are introduced for the modelling of these fields in order to obtain high accuracy of the models. The stiffness matrix of the T-elements is obtained by integration over the element boundaries, which enables to conserve the cylindrical form of the contact surfaces in the computational models. This is important for the accurate enough definition of rolling contact elements for the multibody contact problems like the rolling bearings as described in our previous model [15], in which, however, classical FE model using fine mesh was used to model the local fields.
The basic assumption by the modelling of contact problem is that the local topology of the contact region is not changed by the loading and thus the Hertz theory can be applied for each pair of contacting surfaces. Using this assumption, the