Complete boundary element formulation for normal and tangential contact problems Текст научной статьи по специальности «Физика»

УДК 539.62

Complete boundary element formulation for normal and tangential contact problems

R. Pohrt and Q. Li

Technische Universität Berlin, Berlin, 10623, Germany

The boundary element method as a numerical tool in contact mechanics is widely used and allows for surface roughness to be investigated with very fine grids. However, for every two grid points, influence coefficients have to be employed for every force-displacement combination. In this paper, we derive the matrixes of influence coefficients for the deformation of an elastic half space, starting from the classical solutions of Boussinesq and Cerruti. We show how to overcome complexity problems by using FFT-based fast convolution. A comprehensive algorithm is given for solving the case of dry Coulomb friction with partial slip. The resulting computer program can be used effectively in iterative schemes also in similar problems, such as mixed lubrication and notably improves the applicability of the boundary element method in contact mechanics.

Keywords: boundary element method, influence matrix, elastic deformation, contact mechanics, dry friction, mixed lubrication

1. Introduction

In many problems in contact mechanics and friction physics, one is primarily concerned with processes and deformations located directly at the surface of the contacting bodies. Non-conforming topographies are brought into contact and deform due to the surface stresses. Roughness features are flattened out or slide relative to the opposite topography. The most important step in solving these problems is to find a physical combination of surface stresses and surface deformations. When the correct state is obtained, the resulting subsurface-stresses can also be evaluated. The analytical solutions of stresses and displacements in an elastic half space under a concentrated force or distributed stress are from Boussinesq (1885) and Cerruti (1882) by use of the theory of potential [1, 2]. For the cases of a single concentrated normal or tangential force, the deformations can be calculated with the known analytical solutions. For some special distributed traction, such as the Hertzian contact and cylindrical indenters, analytical solutions exist as well. In most practical cases, e.g. rough surfaces and for arbitrary load, it is often impossible to derive a direct result by integration or superposition principle, so numerical methods must be applied.

The focus of the solution process on only the surface makes the boundary element method a powerful tool in contact mechanics. In contrast to volume-based formula-

tions such as the finite element analysis, only surface points are represented in the computational grid. Only two dimensions instead of three must be modeled, so a great reduction of grid points is given by design. When analyzing rough surfaces, a fine discretization can be crucial [3, 4].

In the process of analyzing the mathematical model with this method, the first important step is discretizing the equations and obtaining an influence matrix. For contact problems, the normal displacement under normal pressure has been always the point of interest, for example in elastohyd-rodynamic lubrication contacts [5, 6]. In 1929 Love analyzed the effect of uniform normal pressure acting on a rectangle area and gave the normal deflection of a general point on the surface [7]. The deformations in other direction or those resulting from tangential stresses have raised less interest, particularly in boundary element method. Li and Berger have given a solution set for triangular areas [8]. To the best of our knowledge, the discrete equations for rectangular uniform grid have not been published yet. In Sect. 2 of this paper we realize the discretization with flat-roofed elements in rectangular areas to derive the influence matrices for 3D deformation of the elastic surface including tangential loading. We will show that for all cases considered, the influence matrices are fully populated, not sparse. As a result, the computational complexity is high. In Sect. 3, we will show ways to overcome the complexity

© Pohrt R., Li Q., 2014

problem using the FFT-based fast convolution. The method is suitable for parallelization and can indeed be ported to systems such as general purpose graphics processing units. Using these techniques, the inverse problem of finding the stresses for given deflection will be treated in Sect. 4 by applying the conjugate gradient method. Section 5 makes use of both procedures to solve the tangential contact problem with partial slip and Coulomb friction.

2. Equations of elastic deformation, discrete formulations

Let us first consider the case of normal loading. Assume a single point force F acting normal to the free surface of an elastic half space at point (x , y). For the deflections in normal (or z-) direction of a surface spot located at (x, y) the solution by Boussinesq [9] reads

( ) 1 -v 1/7 uz( x y) F •

2nG s

(i)

Here 5 = yj(x - x ')2 + (y - y )2 is the distance between the two points, G is the shear modulus and v is Poisson's ratio of the material. We employ the sign convention often used in contact mechanics that positive normal forces or positive deflections act into the elastic half space. In the case of a distributed normal pressure p(x\ y), this expression must be integrated and for a particular spot (x, y) we obtain 1 -v ff p( x', y)

u„ = -

2nG

JP

- dx'dy'.

(2)

In the same way, we find for the tangential deflections in x and y due to normal loading

1 2v-JJp(x^ y) dx'dy^ (3)

uv = —

4nG

A S

_ 1 2vjj y_i_ p (X, y )dx dy.

uy _ — y 4nG

(4)

As

The solution is radial symmetric, so ux in Eq. (3) and uy in Eq. (4) are similar for normal loading.

For pure tangential loading, the equivalent equations have been found by Cerruti. For a distributed surface traction tx (x , y') acting in x-direction, the deflections read [2]:

1 „ffÎ1 — v (x — x')2 +

:2 JH—+ v

u„ _ -

4 nG

xtx (x', y )dx/dy/,

_ 1 2JJv( x — x/)( y — yr) :

4nG A s / /\ 1 / 1 /

XTx (x , y ' ) dx'dy , 1 — 2v rr x — x

¥

"t x ( x, yf ) dx'dy'

(5)

(6) (7)

4nG'A s-

The case of loading in the y-direction leads to analogous equations. Furthermore Eq. (7) is identical to Eq. (3) except for the sign.

With defining a distributed loading ab in domain A, all Eqs. (2)-(7) can be written in the form

ua(x, y) =

= JJ Sab (x - x , y - y)ab (x ', y)dxdy , (8)

A

where u is the deformation, a direction of deformation and b is direction of load (a, b e {x, y, z}).

With a discrete model as shown in Fig. 1, we are interested in deformation ulJa in discrete points (i, j) only. For each grid point (i, f) we assume a uniform stress distribution Glb in a rectangle around this point with dimensions hx and hy.

We write Eq. (8) in the discrete form

uj J Y

(9)

I J

with

YJ+hy/ 2 xV+hx/2

J _ J J Sab(X — X, Y — Y)dXdY (10)

YJ—hy j2 Xv—hxl2

Here the coefficients KjbJ can be generally expressed as a function of parameters

J _ J(i — i ', J — J, G, v, hx, hy). (11)

So for any case of boundary element method with regular mesh grid, we find the influence coefficient of KJJ , by evaluating Eq. (10). For example, the solution of Love [7] for normal deflections under a uniform normal pressure is obtained by integrating (10) with integrand (2). The solution in this case is

J_l—vx

zz 2nG k ln

m + Vk 2 + m 2 ,, n + Vl2 + n 2 -. +1 ln-

n + V k 2 + n 2

m + V l2 + m 2

. k Wk2 + m2 . l + Vl2 + n2

m ln-, + n ln

l + V l2 + m2

k W k2 + n2

(12)

Fig. 1. Discrete element and uniform stress distribution. Within the discrete element, the normal or tangential stress is assumed to be constant. Its influence on the central spot of surrounding elements is described by the influence coefficients in Eqs. (12)-(18)

(13)

where

k = i- i + 0.5, m = j'- j + 0.5, l = i - i - 0.5, n = j' - j - 0.5. This solution is well known and widely used in tribo-logy applications of the boundary element method. For all contact models that consider only normal stresses and deflections, Eq. (12) is sufficient. For example, it can be used to reproduce the Hertzian solution for parabolic (approximately spherical) indenters [10] or the Pohrt-Popov dependency [11] for randomly rough surfaces.

The corresponding discrete equations for tangential loading or tangential deflections are given in the following section. By applying Eq. (10) to Eq. (3) we can obtain K^^ for the two cases of tangential deflection under normal load, which are symmetric with respect to x and y (radial symmetry):

1 - 2v

kV'1} = _

Kxz + h

4nG m ln

,2 . 2 k + m

,2 . 2 l + m

+ n ln

,2 . 2 l + n

, 2 . 2 k + n

+ hx

m n

k| arctg— _ arctg— |+

k k

,1 ^ n ^ m + 11 arctg — _ arctg—

(14)

Kl 'lJ = .

Kyz

1 _ 2v

4nG

k ln

,2 . 2 k + m

, 2 . 2 k + n

+1 ln

+ h

m| arctg--arctg

m

l2 2 + n

\2 2 + m

l

+

m

( . l . k + n| arctg— _ arctg—

I n n

(15)

From the equivalence of Eqs. (3) and (7), it follows that the influence coefficient of normal deflection under tangential load (x-direction) is given by

} =_KJ'J\ (16)

The coefficient for tangential deflection parallel to loading related to Eq. (5) is in discrete form

1

K

J'1 j _

2nG

hx(1 _v) x

k ln

m + Vk2 + m2 n + 4F+7 6 +1 ln

n + \j k2 + n2

m + V l2 + m2

+ h

ln

k + V k2 + m2

l W l2 + m2

+ n ln

l+¥++n 6

k W k2 + n 2

(17)

and for perpendicular to loading related to Eq. (6) it reads

K1}'1 J =

yx 2nG

h2n2 + h^k2 _*,Jh2m2 + h^k2 +

+ *lhy,m2 + hX-l2 _ ^h2n2 + hll2

h

•y'" ' '"x" \"y" ' '"r" /• (18)

In order to give the complete formulations, equations (19)-(21) give the equivalent formulae for tangential loading in the y-direction

K!; J = 1_2V - h. x 4nG I 2

zy

k ln

,2 . 2 k + m

2 , , 2 n + k

+1 ln

/2 i 2 l + n

2 , /2 m +1

+ hy

, . k / ,

m| arctg--arctg— | +

m m

, . l . k

+ n| arctg--arctg—

n n

(19)

KlJ 'lJ' = 1

yy 2nG

hy (1 _v) x

ln

k + V k2 + m2

l + /m2 +12

+ n ln

l6

k+4~n + k2

+ hx

k ln

m + Vk2 + m2 n+ -\/l2+~«2~ ^ +1 ln

n + -\[k2+n

m + V l2 + m2

K1}'1} = K1}'}

K xy Kyx

(20) (21)

With equations (12)-( 18) the deformation components of a discrete point on the surface under a discrete distributed normal and tangential load can be calculated. The procedure has been employed successfully to the tangential deformation of randomly rough surfaces in [12].

3. Fast methods to obtain surface deflections

When applying Eqs. (12)-(18) with the boundary element method, complexity problems arise quickly. Assume a uniform rectangular grid of N x N points. Each one of the resulting N2 points may have 3 components of surface stress and 3 components of surface deflections. Equation (9) can be written in matrix form

(

= EE

Kxx Kx Kyx Ky

Kzx K

zy

Kx

Ky. Kz

YJiJ'

(22)

with the coefficients as derived in the previous section. As every combination of i, j, i and j' must be evaluated, the overall complexity is O(N4) when calculating the surface deflections from the surface stresses. In the past, different schemes have been presented in order to speed up the evaluation of (22) only for Kjl], but these are also applicable to the other 8 components. Lubrecht and Ioannides [13]

have suggested a multilevel-multi-integration procedure that reduces the complexity to O(N2 log N) by first executing the summation at coarser grids and correcting the resulting errors only in the vicinity of (i, j). A very comprehensible explanation and implementation directives can be found in [14]. Another way of accelerating the evaluation of Eq. (22) is the fast convolution technique [15]. The basic idea is to interpret Eq. (8) as a two-dimensional convolution and doing the same in the Fourier space by means of element-wise multiplication [16]

Uj (x, y) = FFT-1 < FFT (Sab )• FFT (o )= =:

=: FCjb (Ob). (23)

The amount of operations is three times a two-dimensional FFT, O(N2 log N) and an element-wise matrix multiplication O(N2). The overall complexity is thus O(N2 log N).

A comparison of both techniques can be found in [17]. In practice we have found the fast convolution to be by far more favorable. First, the implementation is shorter and less error-prone. Second, highly effective implementations of FFT are available for all computing architectures and the porting to parallel systems such as general purpose graphics processing units is easily done and allows for even more time saving.

4. Solving the inverse problem using conjugate-gradient technique

With the convolution, we now have a fast method at hands for obtaining the surface deflections from the surface stresses. The question remains how to solve the inverse problem, finding the stresses necessary in order to obtain a given distribution of deflections.

Polonsky and Keer [18] introduced a quick algorithm for solving the discrete normal contact problem using the conjugate-gradient method. They used an improved version of Venners multilevel multi-integration implemented into the general conjugate-gradient method to obtain olf from uf, but their algorithm can also be used with the fast convolution method instead. The following algorithm is derived from their idea. It takes as input the deflections u at a subset Is of grid points and finds the stresses within that same domain Is that are necessary for the given deflection. All stresses outside of Is are assumed to vanish. The deflections outside of Is can be calculated as a result of the algorithm, but are not restricted a priori. In the following sections, the procedure will be referred to as FC-l(ua, Is). It can be applied to all nine components of Eq. (22), so we will omit the directional indices and mark the iteration steps instead. In contrast to Polonsky and Keer, we do not alter the computational domain, but use a typical conjugate-gradient algorithm instead.

Applying the usual conjugate-gradient algorithm, we start with an arbitrary stress distribution o0

r0 = Ua - FCab (O0 X d0 = r0. (24)

All steps take place only on domain Is. When using the fast convolution in evaluating FCab (o0) and following, we will thus omit the calculated values outside of Is.

For each step k = 0,1,... of iteration, the following operations are executed:

z = FCjb (dk), (25)

ak = M-. (26)

dkz

Now the stresses are updated

Ok+1 = Ok + akdk >

rk+1 = rk ~akz-A new search direction is found

rk+, rk_

Pk =

rk+1rk+1

rkrk

dk+1 = rk+1 "

(27)

(28)

(29)

(30)

+ Mk.

The loop is stopped, when the norm of the residual is below a predefined value

" |< tol. (31)

rk+1

The result is Ok+1 = Ob = FCab(ua, Is)-

5. Solving problems of partial slip with Coulomb friction using boundary element method

From the previous section we gained two possibilities. First, we can quickly evaluate the surface deflections arising from a given pressure distribution FCab (ob). Second, we can find the stresses necessary to create a given surfaces deflection on domain Is using the CG method

FC~l(ua, Is)-

We will now elaborate on how to solve the tangential contact problem with partial slip under Coulomb friction. As the first step, we will simplify the problem by saying that any deflections other than in the direction of stress are omitted (diagonal version of Eq. (22)). This simplifying assumption is satisfied when v = 1/2.

Assume we have solved the elastic contact problem for normal indentation, with a resulting pressure distribution, see [14, 18, 19]: Pif > ^ j e I^

(32)

Pf= 0, (i, j) g Ic with the contact region Ic. When this contact is forced to move tangentially by distance d, an increasing fraction of the contact area will transition from the stick state to slip. In every point, the Coulomb law of friction is valid, so

o<Mp, (33)

where o is the tangential stress in the direction of movement, ^ is the coefficient of friction.

The main idea to solving the problem is to discriminate the region of slip Is and the region of stick (adherence) Ia Is U Ia = Ic. (34)

If we rearrange all the discrete stresses and deformations into a column vector, we can express Eq. (9) for one particular direction as

u = AO. (35)

We can then sort both the vectors u and o is such a way, that they list first the points of no contact (n), then the points where we assume stick (a) and then slip (s). After this rearrangement Eq. (35) reads

A>o Ao Aso

Aoa

Aos

Asa

Ass

0

(36)

Here ua = (d, d, d,...) is a column vector including only values of d. Also we know Os = MP from Coulomb friction. The only unknown stress oa in this equation can be found by solving

AaaOa = d - AsaMP- (37)

Let us now look at how to implement this.

We start with no slip at all, thus Ia = Ic. On each step of iteration, the following operations are performed: First the current region of stick is found, and the tangential stresses within the slip region are deduced from the local normal pressures

Ia = Ic\ Is, (38)

olf =Mpf, (i, f) e Is

Tf -

= 0, (i, f) g Is.

(39)

All tangential deflection that arises from the slip stresses alone usliponly are obtained via multilevel multi-integration or FC and stored in memory

Aso

MP = FCab (Os). (40)

Usliponly

Asa

Ass

This evaluation will generate deflection values in the whole domain I, not only in the slip regions. Following Eq. (37), we define the additional deflection inside the stick zone uadd that must be due to stick region stresses

(41)

These deflections uadd exist only on Ia and we can find the stresses on Ia using the inverse fast convolution (conjugate-gradient)

Oa = FC- (Uadd, Ia). (42)

If any Oaf >MPif, we need to transfer these points from the stick into the slip region and restart the iteration

Is = Is U{(i, f), Of > MPf}. (43)

Since oslip is zero outside the slip region and oa is undefined (zero in the implementation) anywhere outside the stick region, we can find the combined tangential stresses and calculate the combined deflection arising from both:

Uadd = d Usfiponly, (i' f e If

Fig. 2. Discrete partitioning of the contact area in the elastic contact of a randomly rough fractal indenter. Black areas represent zones of stick, gray areas are in sliding state. Surface generated with the Hurst exponent H = 1, tangential displacement was half the way to complete slip ux/uz =0.5|i(2-v)/(2(1-v)), see Eq. (12). Tangential force was Fj(|MFz) = 0.7572

O = Oslip + Oa

(44)

Uc = FCjb (O). (45)

This last step will give d inside of Ia and some values outside, here we can check for points that are in the slip region but should be sticking

Ia = Ia U{(i, f), uf> d}. (46)

If both Eqs. (43) and (46) do not alter the partitioning of Ic into Is and Ia, then we have found the correct distribution of tangential stresses for both the stick and slip region and the resulting tangential deflections. The resulting tangential force Ft can be found by summing up the stresses

Ft = hh EEof

(47)

With this procedure, the partitioning into stick and slip region and the resulting tangential force can be found for a given tangential displacement d in the dry contact.

We tested the algorithm described here with parabolic indenter shapes and obtained the same results as the classical analytical solution from Mindlin [20, Eqs. (98)-(101)]. Furthermore we applied the algorithm to randomly rough surfaces and obtained congruent solutions. Figure 2 shows a sample for a stick-slip partitioning in the contact of a fractal rough surface.

Starting the algorithm we assumed that the normal and tangential problems are decoupled and the contact area does not change due to tangential loading. This is only exact for v = 1/2. In all other cases, there will generally be normal deflection due to the tangential stresses, equation (7) and these tangential deflections generally demand for the normal contact problem to be solved in each step of the iteration.

In a very similar fashion, the algorithm presented here could be applied to mixed lubrication problems. When two non-conforming bodies are in lubricated contact, there will be a gap filled with lubricant. Inside this zone, there are asperity contacts, which can be treated like the dry friction case, but there are also zones of positive gap width, filled with lubricant. The normal stresses in these points can be found via classical elastohydrodynamic lubrication programs [6]. For the tangential deflection, the stresses of these points are deduced from the Reynolds equation, adding one more partitioning in Eq. (36). We have not implemented this idea yet, but it can be used to determine the local deformations due to macroscopic slippage.

6. Conclusion

For the deflection of an elastic half space under the influence of surface stresses, we discretized the equations at the surface for all nine components. With these discrete formulae, direct formulations of boundary element method are possible for contact problems on rectangular, uniform grids. For the inverse problem of given deflections, we presented a conjugate-gradient formulation using the fast convolution method in order to find surface stresses for various contact problems. The case of partial slip with Coulomb friction is elaborated in more detail and an algorithm is given that can be implemented directly to predict partial slip in arbitrary non-conforming contacts [21]. Furthermore it can be included into other adjacent problems, such as mixed lubrication [22] or transient contacts.

This material is based upon work supported by the Deutsche Forschungsgemeinschaft (DFG, Grant No. PO 810/24-1). Q. Li received support through a scholarship from the China Scholarship Council (CSC).

References

1. Johnson K.L. Contact Mechanics. - Cambridge: Cambrigde University

Press, 2003. - 468 p.

2. Popov V.L. Contact Mechanics and Friction. Physical Principles and Applications. - Berlin: Springer, 2009. - 491 p.

3. Zavarise G., Borri-Brunetto M., Paggi M. On the resolution dependence

of micromechanical contact models // Wear. - 2007. - V. 262. - P. 4254.

4. Pohrt R., Popov V. Contact stiffness of randomly rough surfaces // Sci. Rep. - 2013. - V. 3. - P. 3293.

5. Hamrock B.J., Dowson D. Isothermal elastohydrodynamic lubrication of point contacts. Part I. Theoretical formulation // ASME J. Lubr. Tech. - 1976. - V. 98. - P. 223-229.

6. Spikes H.A. Sixty years of EHL // Lubr. Sci. - 2006. - V 18. - P. 265291.

7. Love H.A.E. The stress produced in a semi-infinite solid by pressure on part of the boundary // Philos. Trans. Roy. Soc. London. - 1929. -V. 228. - P. 377-420.

8. Li J., Berger E.J. A Boussinesq-Cerruti solution set for constant and linear distribution of normal and tangential load over a triangular area // J. Elasticity. - 2001. - V 63. - P. 137-151.

9. Boussinesq V.J. Application des Potentiels l'Etude de l'Equilibre et du. S. l. - Gautier-Villar, 1882.

10. Hertz H. Über die Berührung fester elastischer Körper // J. für die reine und angewandte Mathematik. - 1881. - V. 92. - P. 156-171.

11. PohrtR., Popov V.L., Filippov A.E. Normal contact stiffness of elastic solids with fractal rough surfaces for one- and three-dimensional systems // Phys. Rev. E. - 2012. - V 86. - P. 026710.

12. Grzemba B., Pohrt R., Teidelt E., Popov V.L. Maximum micro-slip in tangential contact of randomly rough self-affine surfaces // Wear. -2014. - V. 309. - No. 1-2. - P. 256-258.

13. Lubrecht A.A., Ioannides E. A fast solution of the dry contact problem and the associated sub-surface stress field, using multilevel techniques // ASME J. Tribology. - 1991. - V 113. - P. 128-133.

14. Venner C.H., Lubrecht A.A. Multilevel Methods in Lubrication. -Amsterdam: Elsevier, 2000. - 400 p.

15. Cho Y.-J., Koo Y.-P., Kim T.-W. A new FFT technique for the analysis of contact pressure and subsurface stress in a semi-infinite solid // KSME Int. J. - 2000. - V. 14/3. - P. 331-337.

16. Stanley H.M., Kato T. An FFT-based method for rough surface contact // J. Tribology. - 1984. - V 119(3). - P. 481-485.

17. Wang W.Z., Wang H., Liu Y.C., Hu Y.Z. A comparative study of the methods for calculation of surface elastic deformation // J. Eng. Tribology: Proc. Inst. of Mechanical Engineers. - 2003. - V. 217. -Part J. - P. 145-153.

18. Polonsky I.A., KeerL.M. A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques // Wear. - 1999. - V. 231. - P. 206-219.

19. Allwooda J., Ciftci H. An incremental solution method for rough contact problems // Wear. - 2004. - V. 258. - P. 1601-1615.

20. Mindlin R. Compliance of elastic bodies in contact // ASME J. Appl. Mech. - 1949. - V. 16. - P. 259-262.

21. Wang F.S., Block J.M., Chen W.W., Martini A., Zhou K., Keer L.M., Wang Q.J. A multilevel model for elastic-plastic contact between a sphere and a flat rough surface // J. Tribology. - 2009. - V. 131. - P. 021409.

22. Zhu D., Wang J. Mixed elastohydrodynamic lubrication in finite roller contacts involving realistic geometry and surface roughness // J. Tribology. - 2012. - V 134. - P. 011504.

Поступила в редакцию 17.02.2014 г.

Ceedeuua 06 aemopax

Pohrt Roman, Dr.-Ing., Technische Universität Berlin, Germany, roman.pohrt@tu-berlin.de Li Qiang, MSc, Technische Universität Berlin, Germany, llqq0108@mailbox.tu-berlin.de