A review of boundary elements methodologies for elastic and viscoelastic rough contact mechanics Текст научной статьи по специальности «Физика»

УДК 539.3

A review of boundary elements methodologies for elastic and viscoelastic rough contact mechanics

C. Putignano and G. Carbone1

Imperial College London, London, SW7 2AZ, UK 1 Politecnico di Bari, Bari, 70126, Italy

In this paper, we review the main contributions developed by the authors to solve problems involving contact between bodies with a roughness covering several orders of magnitude. A periodic self-equilibrated formulation is introduced to avoid any border problem. Furthermore, when focusing on the steady-state rolling or sliding contact mechanics between viscoelastic materials, we develop a mathematical formulation, parametrically dependent on the rolling/sliding speed, which employs the same computational techniques introduced for elastic rough contacts. In particular, in both cases, we adopt an ad hoc discretization technique enabling us to reach the full convergence.

Keywords: rough contact mechanics, boundary element methodologies, viscoelasticity

1. Introduction

Rough contact mechanics is an issue drawing wide research efforts, which involve analytical, numeric and experimental investigations. Indeed, many technologically important aspects like friction, wear, contact stiffness and electrical contact resistance are strongly influenced by the roughness of the contacting surfaces. Therefore, strong interest in these topics has been due to new technological challenges to be faced up: these include the manufacturing of novel bio-inspired adhesives [ 1-9], the optimization of seals [10-14], the extreme downsizing of mechanical and electrical devices down to the nanoscales [1, 15]. An interesting case is represented by the micro- and nanomechanical systems (MEMS/NEMS), where the ratio surface/volume grows with the device miniaturization and may entail an increasing influence of surface phenomena. In general, in all such applications, it is then very important to determine the contact main features in terms of contact area, load, separation, friction.

Historically, one of the first main problems, which scientists have faced up, deals with determining the real contact area between two rough solids. Indeed, such an issue is so crucial since the Amonton-Coulomb friction law, asserting the direct proportionality between the friction force and the load, is an immediate consequence of the direct proportionality between the real contact area and the ap-

plied load. Interestingly, despite the huge amount of research contributions, no definitive answer has been provided so far. The first attempt is due to Bowden and Tabor [16] who assume that the asperities of the rough solid surfaces are plastically deformed as soon as they are in contact, and, as they adopt a plasticity model relying on a constant yield stress, the authors end up with a direct proportionality between the load and the area of intimate contact. However, such an approach is likely to be inaccurate, given the large variety of many experimental evidences [17, 18] suggesting that elastic deformations have the leading role in contact mechanics between rough bodies. Furthermore, soft materials (e.g. rubber, rubber-based composites and polymers), important for many engineering applications such as tyres and seals, do not undergo plastic deformation at all. As a matter of fact, it has been necessary to develop more complicate theories proposing general solutions able to determine not only the contact area, but also other really important quantities, such as stiffness, local separation and local stresses. Some of these approaches are based on multi-asperity models [19-23]: they consider the surface constituted by asperities—with a certain distribution of radii of curvature and height distribution—which behave as independent Hertzian punches. However, this class of theories, which describe the rough surfaces only by means of the joint probability distribution of heights, slopes and curva-

© Putignano C., Carbone G., 2014

tures, does not account for the mutual interaction between asperities in contact and also strongly simplify the shape of the contact spots which are reduced to Hertzian contacts, although they actually present fractal-like boundaries (see [24]) and may even assume the shape of non-simply connected regions. This approximation has a direct consequence: in contrast with experiments and numerical calculations, multi-asperity contact models predict linearity between applied load and true contact area only for extremely small applied loads [25]. This inaccuracy clearly increases when the full contact conditions [17] are approached. Persson proposed a completely different approach [26, 27]: this model shows that the contact pressure probability distribution is governed by a diffusive process as the magnification at which we observe the interface is increased. This enables the model to account for the interaction between the contact spots and makes it exact in full-contact conditions. In the case of partial contact, the theory is approximate but still predicts linearity between contact area and load, although the slope of the relation is different from the multi-asperity contact theories predictions and from numerically calculated values (see e.g. [28, 29]).

Because of the limitations affecting analytical theories, wide research efforts have been dedicated to numerical methodologies that may let achieve quantitatively and qualitatively accurate results. The variety of the proposed techniques includes finite element methods [30], boundary elements methods [6, 28, 31, 32], molecular dynamics simulations [33-35] and hybrid approaches [36, 37]. Actually, these methodologies have hard problems in getting full numerical convergence and, therefore, the reliability of predicted results is widely debated. This is mainly due to the huge amount of scales (covering more than 6 order of magnitudes) involved in the rough contact: the consequent computational complexity usually prevents from achieving the accuracy required to provide definitive answers about relations between contact area, load and penetration. The problem is even more serious when viscoelasticity is considered. In this case, the material time-dependent behavior forces, in general, to include also the time domain and this really increase the simulation computational cost.

In this paper, the boundary elements methodologies developed by the authors in [29, 31, 38, 39] to face up these challenges are reviewed. Moreover, mathematical formulation is generalized to consider possible boundary conditions different from the half-space assumption. Our specific aim is indeed to enlighten the main points that have to be considered when carrying out a simulation taking into account the roughness. Particular attention is paid to discuss the convergence problems and the procedures implemented to reduce the computational cost. To this purpose, we show the main peculiarities of the rough contacts in terms of contact area and friction, thus underlying the importance of full convergence to obtain accurate results.

2. Mathematical formulation

2.1. Numerical model of rough surfaces

The boundary element methodologies developed to solve the rough contact mechanics can obviously be applied to real surfaces that have been experimentally acquired. For example, in [39], we measure the surface roughness by means of a confocal microscope and, then, by an atomic force microscope; these data are, afterwards, employed in the numerical simulation. However, when a large number of simulations has to be carried out in order to obtain the statistical properties typical of the rough contact mechanics, e.g. the contact area, it is convenient to adopt a different approach and use numerically generated surfaces.

In order to define a reliable numerical methodology, let us summarize now the main statistical properties of rough surface. A rough surface is described by the height profile h(x), where x = (x, y) is a two-dimensional vector in the surface mean plane z = 0. The height profile is a random process and, therefore, its statistical properties are completely determined by the following auto-correlation functions:

<h(x)>, <h(xi)h(x2)>, <h(xi)h(x2)h(X3)>, ..., (1)

where ( > stands for ensemble average. Without any loss of generality, we assume that (h(x)> = 0. In all those cases, where the correlation functions (1) involving an odd number of h-functions vanish, the heights of the rough profile are proved to be described by the Gaussian distribution of probability:

9(h) =

1

$

(2n)1/2 hn

-exp

,2 %

2h2

(2)

where hrms = (h2 >1/2 is the rms amplitude of the rough surface. Under these conditions, the statistical properties of the rough surface are completely described by the second order auto-correlation function, which can be rewritten as:

(h(xi)h(x2)>=(h(0)h(x2 - Xi)>. (3)

Interestingly, when the rough surface is an ergodic process, the quantity (h(x1)h(x2)> can be also calculated by means of the following relation:

(h(0) h(x)> = — f d2 xh(x + x')h(x'), (4)

A

which represents the standard way used to measure the correlation function of the rough substrate. An equivalent procedure can be developed by employing its Fourier transform, i.e. the power spectral density of the rough surface:

C (q) =—2 f d2 x(h(x)h(0))e~iqx. (5)

(2n)

By assuming perfect isotropy for the rough surface, the surface auto-correlation function simply becomes:

( h(0) h(x )> = f d2qC( q)eiqx = 2nf qC (q) J0 (qr )dq, (6)

0

which therefore only depends on the radial distance r from the origin of the reference frame. In particular, the rms

roughness amplitude (h2is then determined by the following expression:

< h 2 ) = 2nJ qC (q )dq. (7) Now, it is noteworthy to observe that, in practical applications, rough surfaces are often self-affine fractals, at least in a certain range Xmin < X < Xmax of wavelengths. In this case, the statistical properties of the surface are invariant under the affine transformation:

x ^Zx, h (8)

where Z is a scaling factor and H is the Hurst exponent related to the fractal dimension Df of the surface through Df = 3-

- H. For a self-affine surface, we have therefore:

< h (0) h(x)) = Z"2H < h (0) h(Zx)) (9) and the power spectral density can be expressed as:

C(q) = Coq-2(H+1), qo < q < ft, (10)

where q0 = 2n/Xmax and q1 = 2n/Xmin. Equation (10) shows that for self-affine rough surfaces the power spectral density decreases as q "2( H+1).

Figure 1 shows a schematic view of the power spectral density of a real surface in a logC - logq diagram. The slope a = -2(H + 1) of the curve in range q0 = 2n/Xmax < < q < q1 = 2n/Xmin enables the calculation of the fractal dimension Df of the surface through the relation Df = = 4 + a/ 2.

Moving from the previous assumptions, the rms roughness amplitude of the self-affine rough surface can be calculated as: q

<h2) = 2nC0 J q"2H-1dq =

H

1

1

a 2 H q0

2 H

a1

C 1

H

a0

2 H

(11)

if a1 >> a0. Equation (11) allows to calculate the quantity

Fig. 1. A schematic view of the power spectral density C(q) of the surface roughness. The smallest wave vector qL is related to the lateral size of the sample L (qL = 2n/L), the long distance rolloff wave vector is q0 = 2rc/Xmax, and the short distance cut-off wave vector is q1 = 2n/Xmin. For q < q0 the power spectral density C(q) is almost constant, wherease in the range q0 < q < q1? C(q) is represented (in the logC - logq diagram) by a straight line. The slope of this straight-line determines the fractal dimensions of the surface

C0 as a function of the (h2 ) :

c H

C0

-aoH < h2 ).

(12)

On the other hand, as we show in the next section, it is of outstanding importance to estimate the quantity < (Vh)2). It is possible to relate this parameter to the power spectral density C(q) of the isotropic rough surface through the relation:

< (Vh)2 ) = 2nJ dqq3C (q ). (13)

In the case of a self-affine rough profile with a short distance cutoff wave vector q1 we get:

<(Vh)2) = 2rcC0 ] dqq1"2H =

q0

= nC0 ( q2-2H q2-2H) _ nC0

1 - H

2

( a

2-2 Hx

-a0 )

1 - H

ai

,2(1-H )

(14)

Thus < (Vh)2) ~ q2(1 H) is strongly dependent on the short wavelength content of the substrate. For a very interesting review on the nature of surface roughness the reader is referred to [26].

Now, we can employ these relations on a practical level to numerically generate rough surfaces. In details, a periodic surface with Fourier components up to the value q1 = = Nq0 is generated:

h( x) = E

hk

ahke

Mm x

(15)

where qh = (kq0, hq_), x = (x, y). Since h(x) is real we . Moreover for randomly rough sur-

l- h,-k

ahk

must have

face the following relation must be satisfied ( ahkalm ) = 0 with l ^ -h, m ^ -k. Now, we can calculate the power spectral density of surface Eq. (15):

C(q) = S <Kk|2)ô(q-qhk), (16)

hk

from which it follows:

C (q hk ) = <iahk I2 )S( 0). (17)

For isotropic surfaces we have C(q) = C(q) which simply

gives C (q hk ) = C ( q0V h 2 + k 2 ) and assuming self-affine fractal surface, one obtains:

\- h-1

<1 ahk|2 ) =<1 a11|2 )

'h2+k2

(18)

Hence the quantities < | ahk |2) can be determined once given < | a11|2) and the Hurst exponent of the fractal surface. However, to completely characterize the rough profile, we still need the probability distribution of the quantities ahk. We first observe that the condition <ahkalm) = 0 with l ^ -h, m ^ -k is satisfied if the phases 9hk of the complex quantities ahk are random numbers uniformly distributed between 0 and 2n. Incidentally, the condition a-h _k = ahk also implies that | a-h _k | = | ahk | and that the quantities h-k = -9 hk. So what we need now is only the probability distribution of | ahk Of course there could be several

Section plane

Rigid surface

Fig. 2. An elastic half space in adhesionless contact with a rough periodic rigid surface

choices and the simplest one is to assume that the probability density function of |ahk | is just the Dirac's delta function centered at (| ahk |2 i.e.

P (|ahkl) = 8(|ahk|-(|ahk|2 >1/2). (19)

This choice guarantees also that the random profile h(x) has the Gaussian random distribution.

We observe that the numerical surface so generated is self-affine over the entire power spectrum; in other words, we are considering the wave vector q L coincident with the long distance roll-off wave vector q0. Recent contributions [40, 41] propose a different approach considering qL ^ q0 and, therefore, a surface with a more general power spectral density (see Eq. (1)). Generally speaking, this could be useful to generate surfaces with a higher degree of er-godicity and, consequently, to decrease the number of realizations to be averaged in each numerical simulation.

2.2. Mathematical formulation for elastic materials

We study the adhesionless contact between a periodic numerically generated isotropic randomly rough rigid surface and a linear elastic layer (Fig. 2). Interestingly, we observe that elastic displacement of the layer is the sum of two terms the first is just equal to the average displacement um (z) the second v(x, z) = u(x, z) - um (z), where x is the in-plane position vector, is just the additional displacement caused by the asperities induced deformation of the rough surfaces. Focusing only at the interface z = 0 we observe that the maximum value of the normal displacement v (x) = v z (x, z = 0) is the penetration A of the rough surface into the half space. Furthermore, we notice that, to avoid border effects, a periodic formulation is adopted by considering a periodic domain constituted of square cells D = [(-V 2, V 2) x (-V 2, V 2)].

Now, we can formulate the problem by relating the elastic displacement field at the interface to the interfacial normal stress a(x) = azz (x, z = 0) as [2, 42]:

uz(x) = f d2s0(|x -s\/h)G(x -s)a(s), xe D, (20)

D

where G(x) is the Green function, 0 (| x \/h) is a corrective parameter introduced to account for the slab thickness h:

0 (r/h) =f droS(roh/r)J0(ro), (21)

0

were S (roh/ r) is a correction term which accounts for different constraints or boundary conditions (see [2]) and J0 (ro) being the 0th order Bessel function. We stress that, as should be expected, 0 (r/h) approaches the unit value at relatively low values of r/h and rapidly vanishes as r/h ^ «>.

However, we should observe that, in case of an infinite layer under periodic load conditions, the mean displacement um of the elastic body at the interface is unbounded, and therefore only the term v(x) = uz (x) -um is finite. Since our focus is just on v(x), we need to reformulate the problem in such a way that only the quantities a(x) and v(x) appears. To this end, let us consider that from Eq. (20) one obtains

um =f d2sGma(s).

D

where

Gm = A-2 f d2 xG (x),

D

so that one can write:

v(x) = f d2s0 (| x - s \ /h) L(x - s)a(s), x e D, (22)

D

where L(x) = G(x) - Gm is equal to the elastic displacement at the interface uz (x) caused by a periodically applied self-balanced normal stress distribution a(x) = 8(x) --A-2. Indeed using Eq. (20) one obtains

uz (x) = f d2sG(x - s)[S(s) - A-2] = G(x) - Gm = L(x).

D

To solve the contact problem we control penetration depth A and discretize the domain D in small squares of nonuniform size. The unknown stress acting on each single

square is assumed to be uniformly distributed on it. The discretized version of Eq. (22) becomes:

V = L Cj > (23)

where v{ is the central displacement of each square, c- is the uniform stress on the square and Lu is the elastic response matrix (i.e. the matrix associated to the discretized version of Eq. (22)), caused by the self-balanced applied load shown in Fig. 3.

Lj can be calculated by recalling the Love solution (see [43]) which gives the elastic displacement due to a uniform pressure acting on a rectangular area, and summing up the contribution of each elementary cell D to take into account the periodicity of the problem. Thus Lj = = lj - lm , where:

i -v2 E

= UTS. ¿1« j+d'> x

X ln

% + dj + [(£j + dj )2 + (n j + dj)

nj - dj + [(v. + dj )2 + (n- - d. )2W2

+ (nj + dj )ln

ij j

%j + dj + [(n,j + dj )2 + (tj + dj )2]

-dj yy

2

2^2

+ (tj -d:)ln

+(n

and

lm

X ln

j )ln

1 -v2

> - dj + [(tj. - dj )2 + (n- +dj )¥2

nij - dJ + [(tu - dj )2 + (nij - dj)2]12

n. + dj + [(tj- - dj )2 + (% +dj )2P

'tj - dj + [(tj- - dj )2+(nj - dj^

t« + d, + [(S„ - df )2 + (n

- d ,)2]'/2

(24)

nE

dj X

2

s s +1)x

k +1 + [(h +1)2 + (k +1)2]

2 "i1/2

k -1 + [(h +1)2 + (k -1)2]

2]1/2

+X( k +1) ln +X( h -1) ln +X( k-1)ln

h +1 + [(k +1)2 + (h +1)

,2"12'

vh -1 + [(k +1)2 + (h -1)2]1/2 ' k -1 + [(h -1)2 + (k -1)2]1/2

k +1 + [(h -1)2 + (k +1)2]1/2 h -1 + [(k -1)2 + (h -1)2]1/2

h +1 + [(k -1)2 + (h +1)2]

2]V2

(25)

where dj is the size of the elementary square cell and =1 xj -Xi1 +Xh and nj =1 yj -y1 +Xk.

Equation (23) can be employed to calculate the stresses at the interface if the displacements ut are known. However, our problem belongs to the class of mixed boundary problems since the real contact area has to be determined. Therefore, as we show in Fig. 4, the solution is sought by means of an iterative incremental procedure: (i) fix the displacement A i, (ii) evaluate the contact area as the intersection between the deformed elastic layer, estimated with

Fig. 3. The self-balanced load on the elementary cell utilized to evaluate the elastic response matrix Lj

respect to the elastic solution found previously for the penetration Ai-1 < Ai, and the rigid rough surface, (iii) calculate the displacements in the contact areas as vt = ht -- hmax +A, where hi = h (xi), hmax the maximum height of the rough profile, (iv) solve Eq. (23) to calculate the stress distribution c- in the contact areas, (v) calculate the displacements vt = L.c- out of the contact areas, (vi) update the contact area at every step by eliminating the elements with negative pressure and adding those for which there is penetration.

The inversion of the matrix L:j is done only for those points which belongs to the contact areas, thus enabling a strong reduction of the computational efforts. To invert Eq. (23) we use an iterative method based on the Gauss-Seidel scheme.

2.3. Mathematical formulation for viscoelastic materials

A linear viscoelastic material can be described by the following integral equation [44]:

8(t) = J dTJ(t -t)C(t),

(26)

where 8(t) is the time-dependent strain, c(t) is the corresponding stress (the symbol • stands for the time derivative), and the function J(t) is the creep function which satisfies causality, i.e. J(t < 0) = 0. The most general form of J(t) is [44]:

J(t) = H(t) -i--+J dTC(T)exp(-t/T) , (27)

_ E0 0 _ where H (t) is the Heaviside step function, the real quantity E0 is the rubber elastic modulus of the material at zero frequency, C(t) is a strictly positive function usually defined as the creep (or retardation) spectrum [44, 45], and t is the relaxation time, continuously distributed on the real axis.

To employ Eq. (27) in a numerical procedure, a discrete version is introduced by assuming

C ( T):

s Ck 5( t-

k

'Tk X

Fig. 4. Iterative procedure to calculate the contact area

which gives: J (t ) = H (t )

1

1

J_

E

-X Ck exp(-t/Tk )

k=1

(28)

This form of J(t) corresponds to the general linear vis-coelastic model reported in Fig. 5.

As commonly proposed in the literature [44, 45], it is useful to take the Fourier transform of the equations presented above. In particular, equation (27) can be Fourier transformed to obtain 8(œ) = a(œ)/E(œ) with E(œ) = = [iœJ(œ)]-1. Note that the Fourier transform f (œ) of a function f (t) has been defined as f(œ) = f dtf(t) x x exp(-iœt ). Using Eq. (27), we can easily write the following equation for the viscoelastic complex modulus E (œ) as

Fig. 5. The generalized linear viscoelastic model constituted by a spring in series with n Voigt elements (the latter consists of a Hookean spring in parallel with a Newtonian dashpot)

E ( œ) E

0 0

J dT iœTC (T) 1

l + iœT

f dT-1

C (t) . + i œT

(29)

As shown in Fig. 6, we observe that at "low" frequencies the material is in the "rubbery" region, where E1 = Re E (ro) is relatively small and approximately constant, while the imaginary part E2 = Im E (ro) is very small and usually negligible. In this region, the viscoelastic dissipation can be therefore neglected. At very high frequencies the material is elastically very stiff. In this "glassy" region E1 (ro) is again nearly constant but much larger than in the rubbery region, and E2(ro) becomes small again, i.e. viscoelastic dissipation is negligible. Finally, in the intermediate frequency range (the so called "transition" region), the loss tangent ImE(ro)/Re E(ro) is very large (Fig. 6, b), and, indeed, energy dissipation during the sliding contact is originated by this region of the viscoelastic spectrum. We notice that the linear viscoelastic solid characterization requires a discrete version of Eq. (29). In details, by writing

C (T) = E Ck S( T-Tk),

k

one obtains

1/ E (ro) = 1/ E.+ECk/ (1 + iroT k).

k=1

Once described the reology of linear viscoelastic solids, let us focus on the mathematical formulation deve-

Fig. 6. The real E1 = Re[E(ro)] and the imaginary E2 = Im[E(ro)] parts of the viscoelastic modulus E(ro) of a typical rubber-like material (a); the loss tangent E2(ro)/E1(ro) (b)

loped to solve the viscoelastic contact mechanics. Basically, by recalling the translational invariance and the elas-tic-viscoelastic correspondence principle [45], we may formulate the general linear-viscoelastic contact problem between a rigid indenter and a viscoelastic slab as:

t

t(x, t) =| dxjd2x J(t - T)G(x - x') g(x', t),

(30)

where x is the in-plane position vector, t is the time, u(x, t) is the normal surface displacement of the viscoelastic solid, a(x, t) is the normal interfacial stress, J(t) is given in Eq. (27) and the quantity G(x) is the Green's function. Now, assuming that sliding occurs at constant velocity and neglecting non-uniform temperature effects, one can show that equation (30) can be rewritten in the form:

u(x) = f d2x'G(x - x', v)a(x'). (31)

The authors have shown in [46] that the kernel G(x, v), which depends parametrically on the sliding speed v, has the following form:

G(x, v) = -

1 -v2 n

1

+ J dT C ( t) J dz-

1

0

-0

h

0

0+

| x + vTz | h

exp(-z) S. (32)

| x + v Tz |

For further details regarding the mathematical procedure to obtain G(x, v), the reader is referred to [46].

Here, we observe that, once the Green's function G(x, v) is explicitly given, the viscoelastic problem can be solved by following the same approach already employed by the authors for elastic contact mechanics [29, 31]. This

strategy consists in discretizing the contact domain in N square cells and, then, employing the usual solution scheme for boundary element methodologies. In details, as done before for elastic materials, assuming that in each square cell the normal stress a is constant and equal to ak = = a(Xk), where Xk is the position vector of the center of the square cell Dk, the normal displacement uf = u (Xi) at the center of the i-th square cell can be expressed as:

1 -V2 L

ui =--

n k=1

0

| x - X:

x J d2X'—

D,. | Xi

1

V

1

x0

X',

X + VTZ - Xk

— + J dTC(t) J dz exp(-z):

0+

h

1

f d2X'-

D | X(. + VTz - X' |

(33)

JDk

where the terms

1

J d2X'- :

D | Xi - X'|

1

and J d2X'- ,

D | Xi + VTz - X '|

can be easily calculated by exploiting Love's solution for elastic materials [43], as shown in [31]. In such a way, the problem is reduced to a system of linear equations of the type:

Ui = Lik (v )a k, (34)

where the response matrix Lk (v) parametrically depends on the velocity v. Equation (34) can be easily solved, together with the determination of the real contact area, by employing the iterative procedure previously described.

The tangential friction force Ft needed to slide, at constant velocity v = | v the viscoelastic block against the rough rigid substrate, can be calculated considering that the energy per unit time W provided by the external tangential applied force Ft must balance the energy per unit time dissipated as a consequence of viscoelastic response of the material.

The friction force is, therefore, easily determined as:

W = Ft v = -J d2 xa( x - vt ) = J d2 X ct( X) v -Vu ( X),

du(x - vt) dt

(35)

where ^ is the contact domain. For isotropic surfaces, we can assume, without any loss of generality, v = vi, where i is the unit vector of the X axis, and write the final relation as:

Ft =J d2 X a( X) dX. (36)

The friction coefficient is then calculated as ^ = FjFn, where Fn is the external applied load.

2.4. Adaptive non-uniform discretization technique

An important issue to evaluate the efficiency of a numerical methodology is the capability of converging and, on the other side, the computational costs related to the

methodology. These two points are often opposed, but both related to the adopted discretization methodology. Given the domain D, the classic methodology (see, for example, [47]) relies on a base uniform grid where NxN square elements are always allocated in memory. Such a type of mesh allocates also elements not in contact. Wriggers [48] suggested an alternative method of meshing process using an adaptive uniform grid considering only the elements in contact: in such a way, it is possible to obtain good memory savings. Moving from the approach developed by Carbone and Mangialardi [42] and Carbone et al. [6] for the case of 1D contacts, we improve Wriggers' methodology by introducing a non-uniform adaptive discretization technique: this has a coarse mesh in the inner part of each contact area and a much finer one where the gradients of the interfacial stresses tend to diverge (i.e. at the borders of the contact clusters). Indeed, this procedures based on pure geometric criteria consisting in discretizing the contact areas with the smallest step size at the border and with a following coarsening (with a fixed scale factor 2) as we move toward the inner part of each contact spot. In particular, the smallest step is much smaller than the shortest wavelength of the rigid rough substrate. As we show in the following section, a non-sufficiently small step size can prevent the methodology from obtaining the full numerical methodology.

To elucidate the computational advances related to the proposed meshing scheme, Figure 7 shows the number of elements of discretization of the domain D in terms of the imposed penetration A for a self-affine fractal surface with Hurst coefficient H = 0.8, root mean square roughness <h2}1/2 = 13 |m, q0 = 2n102 and N = 128. The comparison between the different methods of discretization (uniform grid, adaptive uniform grid and adaptive non-uniform

Fig. 7. Number of elements of discretization of the domain D in terms of the imposed penetration A, for a self-affine fractal surface with the Hurst coefficient H = 0.8, root mean square roughness 13 ^m, q0 = 2rc-102 and N = 128. Comparison between different methods of discretization: uniform grid (1), adaptive uniform grid (2) and adaptive non-uniform grid (3)

grid) show how the total number of elements significantly reduces, thus entailing a great benefit in terms of computational costs. For example, with a penetration A = 400 | m the cells allocated in memory are about 7 • 105. Such number would be 2.75 -106 with an adaptive uniform mesh and 4.2 • 106 with a classical uniform grid.

3. Results and discussion

In this section, the boundary element methodologies previously defined are shown fully capable of obtaining reliable results in terms of contact area, contact stiffness, and friction. Furthermore, given the aims of this paper, related to define the main points of a numerical methodology for rough contact mechanics, particular attention is paid to discuss the physical conditions ensuring the real numerical convergence of the methodology.

3.1. Contact area in elastic contacts

Now, let us focus on the contact area between a rough elastic half space and the rigid surface. This issue is really prominent in a large series of phenomena, including friction and wear. Calculations are carried out for self-affine

fractal surfaces with spectral components in the range

2 -1

q0 < q < q1, where q0 = 2n • 10 m , q1 = Nq0 and N number of scales (or wavelengths). For each surface, seven different realizations have been averaged. The half space has been assumed linearly elastic with Young's modulus E = = 0.4 GPa and Poisson's ratio v = 0.45.

Our analysis starts observing that both multi-asperity and Persson's theories suggest that, beside the elastic properties of the material, the rough surface geometry influences the contact area versus load relation only through the quantity m2 = <Vh}2/2:

A K Co (37)

10_ m E *

Ao V2

where c0 is the average pressure at the interface (i.e. c0 = F/A being F the applied load) and k is the slope parameter, which takes the value k = V2tc in the case of multi-asperity contact models and k = yj8/n in the case of Persson's theory. Indeed, our calculations show that neither Persson's theory nor multi-asperity contact theories

Fig. 8. The constant of proportionality k as a function of H, for N= 64 and hrms = 13.26 |m. BGT and Persson predicitons are also plotted for comparison

Fig. 9. The constant of proportionality k as a function of N, for H= 0.85 and hrms = 12.8 ^m. BGT and Persson predicitons are also plotted for comparison

predict the right value of k, as clearly shown in Figs. 8, 9 and 10, where k is plotted as a function of the Hurst exponent H, of the wavelength number N and of the root mean square roughness hrms, respectively. In all cases we find an almost constant value k = 2, which lies in between the two values predicted by multi-asperity models and Persson's theory respectively.

An intermediate value of k was also founded in [30, 35]. In [30] k was significantly affected by H, and also by the surface roughness hrms. However, figures 8, 9 and 10 do not confirm this result, which instead is in (at least partial) agreement with [35], where the focus of the investigation was on the influence of the fractal dimension analyzed with a completely different method based on a molecular dynamics Green's function approach. In general, a wide debate deal with this issue [40]. However, it should be noticed that this difference between our calculations and those presented for example in [30] or in other studies [28], could be due to an incorrect estimation of the probability density function P(a) of the interfacial normal stress a: in all these cases, the probability density function does not vanish as the stress a leans toward zero, although, as analytically proved [27, 49, 50], this must necessarily occur. Our numerical procedure instead predicts a very good proportional decrease of P(a) at small a, as shown in Fig. 11. We observe also that the calculated stress probability density function P (a/ E *) is well fitted both by a double Gaussian distribution ( solid line), as suggested by the Persson's theory [26], and by the Rayleigh distribution, as proposed by Greenwood in [20].

Fig. 11. The probability density function P(a/E*) of the dimen-sionless interfacial normal stress distribution a/ E* for a rough self-affine surface with H = 0.8, hrms = 13.26 ^m, N = 64. The fit of the numerical data is carried out by a double Gaussian distribution (solid lines) [26] and by a Rayleigh distribution (dashed lines) [20]

It should be pointed out that the probability density function P (a/ E *) has strong prominence on the definition of the coefficient k and, more in general, on the convergence of the numerical methodology. In this respect, given a fractal surface with H = 0.8 and N = 64, figure 13 shows the stress probability function P(a/E*) for different values of the discretization parameter M, marking the dimension A) /(M x M) of the smallest cell of the discretized domain. If M is not large enough, the numerical method is unable to compute a correct trend of the probability distribution and, therefore, a correct estimation of the parameter k (see Fig. 12). Therefore, for the method to achieve the numerical convergence a sufficiently fine mesh needs to be defined into the domain of calculation.

This indeed confirm the necessity of having a technique, such the adaptive non uniform mesh, that allows a sufficiently fine discretization of the contact domain. Conversely,

Fig. 10. The constant of proportionality k as a function of hrms, for H = 0.8 and N = 64. BGT and Persson predictions are also plotted for comparison

Fig. 12. Coefficient of proportionality k for different values of the discretization parameter M

P(a/E*)

M

o 128

a 256 v 512 ♦ 1024 □ 2048

Fig. 13. Interfacial stress probability distribution P(g) for a fractal surface with H = 0.8, N = 64 and different mesh refinements

numerical results could be really far from the actual physical solution.

3.2. Viscoelasticfriction

When considering the sliding contact between viscoelastic bodies, viscoelastic friction must be accounted for. In this section, we focus on the contact of a rigid rough fractal surface sliding over a viscoelastic half space (h ^ +<») marked by only one relaxation time. In details, the viscoelastic material is characterized by the following values of E^ = 107 Pa, Ej E0 = 3, and t = 0.01 s. As for the rough surface, we employ self-affine fractal surfaces numerically generated by means of the spectral method described in [31]. These surfaces have spectral components in the range q0 < q < q1, where q0 = 2n/L, the side of the square computational cell is L = 0.01 m, q1 = Nq0 and N number of scales (or wavelengths). In particular, results shown in this section are obtained with N = 64. Indeed, in Fig. 14, we analyze the viscoelastic friction as a function of the dimen-sionless speed £= vt0/L for a fixed normal load Fn = = 0.3 N. We observe a bell-shaped curve that vanishes for very low and very high speeds, i.e. when the solid behaves as an elastic material. In details, for very low speeds, the material shows an elastic behavior with the soft elastic modulus E0; on the other side, for very high speed, the half-space is still elastic, but characterized by the hard vis-coelastic modulus E^.

The boundary element methodology is, therefore, shown capable of carrying out numerical simulations of the viscoelastic friction. Interestingly, in [39], experimental outcomes seem to validate the reliability of these measures. However, also in this case, it should be carefully checked the numerical convergence of the methodology. Indeed, in [39], we show that all the roughness scales contribute to the viscoelastic friction. As a matter of fact, if the smallest

M^Mrnax

Fig. 14. Viscoelastic friction coefficient as a function of the di-mensionless sliding speed £ for a constant normal load P = 0.3 N

scales are not discretized with a step size small enough, the friction evaluation could be strongly inaccurate. In this respect, the proposed boundary element methodology relying on the adaptive non uniform mesh seems to be quite effective in obtaining the full numerical convergence.

4. Conclusions

In this paper, we review the boundary element methodologies developed by the authors to solve the contact mechanics between elastic and viscoelastic bodies. We introduce an integral equation relating the surface displacement and the interfacial stresses. We adopt a self-equilibrated periodic formulation: this enables us to avoid any border problem. Furthermore, in the case of steady-state rolling or sliding contact between viscoelastic bodies, the mathematical formulation, depending parametrically only on the viscoelastic layer speed, allows to formally solve the problem with the same techniques already developed for the elastic case. In particular, either for elastic and viscoelastic materials, we employ an adaptive non-uniform mesh, being coarser in the inner part of each contact cluster and finer at the borders. This entails not only outstanding computational savings, but also enable us to reach the full numerical convergence. Conversely, if the convergence is not reached, usually due to an non-sufficiently small step size, significant errors can be carried out in evaluating important quantities like the contact area and the viscoelastic friction.

References

1. Geim A.K., Dubonos S.V., Gricorieva I.V., Novoselov K.S., Zhukov A.A., Shapoval S. Yu. Microfabricated adhesive mimicking gecko foot-hair // Nature Mater. - 2003. - V. 2. - P. 461-463.

2. Carbone G., Mangialardi L. Adhesion and friction of an elastic halfspace in contact with a slightly wavy rigid surface // J. Mech. Phys. Sol. - 2004. - V. 52. - No. 6. - P. 1267-1287.

3. Carbone G., Pierro E., Gorb S. Origin of the superior adhesive performance of mushroom shaped microstructured surfaces // Soft Matter. -

2011. - V. 7. - No. 12. - P. 5545-5552.

4. Carbone G., Pierro E. Sticky bio-inspired micropillars: Finding the best shape // Small. - 2012. - V. 8. - No. 9. - P. 1449-1454.

5. Carbone G., Pierro E. Effect of interfacial air entrapment on the adhesion of bio-inspired mushroom-shaped micro-pillars // Soft Matter. -

2012. - V. 30. - No. 8. - P. 7904-7908.

6. Carbone G., Scaraggi M., Tartaglino U. Adhesive contact of rough surfaces: comparison between numerical calculations and analytical

theories // Euro. Phys. J. E. Soft Matter. - 2009. - V. 30. - No. 1. -P. 65-74.

7. Carbone G., Mangialardi L., Persson B.N.J. Adhesion between a thin elastic plate, and a hard randomly rough substrate // Phys. Rev. B. -2004. - V. 70. - No. 12. - P. 125-407.

8. Afferrante L., Carbone G., Demelio G., Pugno N. Adhesion of elastic thin films: Double peeling of tapes versus axisymmetric peeling of membranes // Tribol. Lett. - 2013. - V. 52. - P. 439-447.

9. Putignano C., Afferrante L., Carbone G. Peeling stability of pre-stressed

tapes // Submitted to Beilstein J. Nanotechnol. - 2014.

10. Bottiglione F., Carbone G., Mangialardi L., Mantriota G. Leakage mechanism in flat seals // J. Appl. Phys. - 2009. - V. 106. - P. 104902.

11. Bottiglione F., Carbone G., Mantriota G. Fluid leakage in seals: An approach based on percolation theory // Tribol. Int. - 2009. - V. 42. -No. 5. - P. 731-737.

12. Lorenz B., Persson B.N.J. Time-dependent fluid squeeze-out between solids with rough surfaces // Euro. Phys. J. E. - 2010. - V. 32. -No. 3. - P. 281-290.

13. LorenzB., Persson B.N.J. On the dependence of the leak rate of seals on the skewness of the surface height probability distribution // EPL. -2010. - V. 90. - P. 38002.

14. Lorenz B., Persson B.N.J. Leak rate of seals: Effective-medium theory and comparison with experiment // Euro. Phys. J. E. - 2010. - V.31.-No. 2. - P. 159-167.

15. Bhushan B. Springer Handbook of Nanotechnology. - Berlin: Springer, 2004. - 1222 p.

16. Bowden F.P., Tabor D. The friction and lubrication of solids // Proc. R. Soc. A. Math. Phys. Eng. Sci. - 1939. - V. 169. - P. 391-413.

17. Archard J.F., Hirst W. The wear of metals under unlubricated conditions // Proc. R. Soc. A. Math. Phys. Eng. Sci. - 1956. - V. 236. -P. 397-410.

18. Fuller K.N.G., Tabor D. The effect of surface roughness on the adhesion of elastic solids // Proc. R. Soc. Lond. A. - 1975. - V. 345. -1642. - P. 327-342.

19. Greenwood J.A., Williamson J.B.P. Contact of nominally flat surfaces // Proc. R. Soc. Lond. A. - 1966. - V. 295. - P. 300-319.

20. Greenwood J.A. A simplified elliptic model of rough surface contact // Wear. - 2006. - V. 261. - P. 191-200.

21. Greenwood J.A., Putignano C., Ciavarella M. A Greenwood and Williamson theory for line contact // Wear. - 2011. - V. 270. - P. 332334.

22. Bush A. W., Gibson R.D., Thomas T.R. The elastic contact of a rough surface // Wear. - 1975. - V. 35. - P. 87-111.

23. Ciavarella M., Delfine V., Demelio G. A "re-vitalized" Greenwood and Williamson model of elastic contat between fractal surfaces // J. Mech. Phys. Sol. - 2006. - V. 54. - P. 2569-2591.

24. Mandelbrot B.B. The Fractal Geometry of Nature. - New York: W.H. Freeman and Company, 1982. - 170 p.

25. Carbone G., Bottiglione F. Asperity contact theories: Do they predict linearity between contact area and load? // J. Mech. Phys. Sol. -2008. - V. 56. - P. 2555-2572.

26. Persson B.N.J. Theory of rubber friction and contact mechanics // J. Chem. Phys. - 2001. - V. 115. - P. 3840-3861.

27. Persson B.N.J. Contact mechanics for randomly rough surfaces // Surf. Sci. Rep. - 2006. - V. 61. - P. 201-227.

28. Paggi M., Ciaveralla M. The coefficient of proportionality k between real contact area and load, with new asperity models // Wear. -2010. - V. 268. - P. 1020-1029.

29. Putignano C., Afferrante L., Carbone G., Demelio G. The influence of the statistical properties of self-affine surfaces in elastic contact: a numerical investigation // J. Mech. Phys. Sol. - 2012. - V. 60. - No. 5. -P. 973-982.

30. Hyun S., Pei L., Molinari J.-F., Robbins M.O. Finite-element analysis of contact between elastic self-affine surfaces // Phys. Rev. E. -2004. - V. 70. - P. 026117.

31. Putignano C., Afferrante L., Carbone G., Demelio G. A new efficient numerical method for contact mechanics of rough surfaces // Int. J. Sol. Struct. - 2012. - V. 49. - No. 2. - P. 338-343.

32. Pohrt R., Popov V.L. Normal contact stiffness of elastic solids with fractal rough surfaces // Phys. Rev. Lett. - 2012. - V. 108. - P. 104301.

33. Yang C., Persson B.N.J. Molecular dynamics study of contact mechanics: Contact area and interfacial separation from small to full contact // Phys. Rev. Lett. - 2008. - V. 100. - P. 024303.

34. Yang C., Tartaglino U., Persson B.N.J. A multiscale molecular dynamics approach to contact mechanics // Euro. Phys. J. E. Soft Matter. - 2006. - V. 19. - No. 1. - P. 47-58.

35. Campana C., Muser M.H. Contact mechanics of real vs. randomly rough surfaces: A Green's function molecular dynamics study// Euro. Phys. Lett. - 2007. - V. 77. - No. 3. - P. 38005.

36. Hyun S., Robbins M.O. Elastic contact between rough surfaces: Effect of roughness at large and small wavelengths // Tribol. Int. - 2007. -V. 40. - P. 413-1422.

37. Luan B.Q., Hyun S., Molinari J.F., Bernstein N., Robbins M.O. Multiscale modeling of two-dimensional contacts // Phys. Rev. E. -2006. - V. 74. - P. 046710.

38. Scaraggi M., Putignano C., Carbone G. Elastic contact of rough surfaces: A simple criterion to make 2D isotropic roughness equivalent to 1D one // Wear. - 2013. - V. 297. - No. 1-2. - 5. - P. 811-817.

39. Carbone G., Putignano C. Rough viscoelastic sliding contact: theory and experiments // Phys. Rev. E. - 2014. - V. 39. - P. 032408.

40. Yastrebov V.A., Anciaux G., Molinari J.-F. Contact between representative rough surfaces // Phys. Rev. E. Stat. Nonlin. Soft Matter Phys. - 2012. - V. 86. - No. 3. - P. 035601(R).

41. Pastewka L., Prodanov N., Lorenz B., Muser M.H., Robbins M.O., Persson B.N.J. Finite-size scaling in the interfacial stiffness of rough elastic contacts // Phys. Rev. E. - 2013. - V. 87. - P. 062809.

42. Carbone G., Mangialardi L. Analysis of adhesive contact of confined layers by using a Green's function approach // J. Mech. Phys. Sol. - 2008. - V. 56. - No. 2. - P. 684-706.

43. Johnson K.L.J. Contact Mechanics. - Cambridge: Cambridge University Press, 1985. - 452 p.

44. Christensen R.M. Theory of Viscoelasticity. - New York: Academic Press, 1971. - 245 p.

45. Williams M.L., Landel R.F., Ferry J.D. The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids // J. Amer. Chem. Soc. - 1955. - V. 77. - P. 37013707.

46. Carbone G., Putignano C. A novel methodology to predict sliding/ rolling friction in viscoelastic materials: theory and experiments // J. Mech. Phys. Sol. - 2013. - V. 61(8). - P. 1822-1834.

47. Borri Brunetto M., Chiaia B., Ciavarella M. Incipient sliding ofrough surfaces in contact: a multi-scale numerical analysis // Comp. Meth. Appl. Mech. Eng. - 2001. - V. 190. - P. 6053-6073.

48. Wriggers P. Computational Contact Mechanics. - Chichester: Wiley & Sons Ltd, 2002. - 520 p.

49. Persson B.N.J., Bucher F., Chiaia B. Elastic contact between randomly rough surfaces: Comparison of theory with numerical results // Phys. Rev. B. - 2002. - V. 65. - P. 184106.

50. Manners W., Greenwood J.A. Some observations on Persson's diffusion theory of elastic contact // Wear. - 2006. - V. 261. - P. 600-610.

nocTynH^a b pe^aKUHro 10.03.2014 r.

CeedeHua 06 aemopax

Putignano Carmine, Intra European Marie Curie Research Fellow, Imperial College London, UK, c.putignano12@imperial.ac.uk Carbone Giuseppe, Associate Prof., Politecnico di Bari, Italy, giuseppe.carbone@poliba.it