Friction on a Mesoscopic Scale

• Master equation approach

• Some examples:

  • Steplike distribution

  • Rectangular distribution

  • Singular case

  • Two delta-functions

  • Friction loop

• Generalizations

  • Delay in contact formation

  • Aging of junctions

  • Role of temperature

• Distribution of thresholds

  • Rough surfaces

  • Dry contact

  • Lubricated surfaces

• Stick-slip kinetics

  • Motion equation

  • Elastic instability

  • Slip details

  • Interacting contacts

  • Discussion

• Concerted slips

• Complete set of equations

• Conclusion

Master equation approach

Above we have described a simple mesoscopic model – the Burridge-Knopoff spring-block model, known also as the earthquake (EQ) model – which bridges the gap in scales and describes the main experimental observations, such as stick-slip and smooth sliding, in terms of the properties of local contacts (junctions, asperities, solidified lubricant "bridges") that break at a critical force. The drawback of EQ simulations, however, is that heavy calculations with different parameter sets are required to determine the main features of the model, and it is hard to draw conclusions of general validity. Moreover, almost all studies based on the EQ model assume that all contacts have identical properties for simplicity. However, as we show below, this limit is singular and may lead to qualitatively incorrect conclusions.

Earthquake model. The EQ model describes the interface between the bottom of the solid block and the fixed substrate. It assumes that the interaction occurs through  N  asperities that make contacts with the substrate. Each asperity is characterized by its contact area  A_i  and an elastic constant  k_i,  schematized by an elastic spring, which can be estimated from  k_i ~rc²√A_i,  where  r  is the mass density and  c  is the transverse sound velocity of the material which forms the asperity. When the bottom of the solid block is moved by  X,  the stretching  x_i  of an asperity, i.e. its elastic deformation with respect to its relaxed shape, increases. The force at the contact grows as  f_i=k_ix_i  until it reaches the threshold value  f_siµA_i  at  x_si=f_si /k_i µ √A_i;  at this point the contact rapidly slides, and  f_i  and  x_i  drop to a small value before a contact is formed again.

Let  P_c(x_s)  be the normalized probability distribution of values of the thresholds  x_si  at which contacts break. Firstly we consider the model in the quasi-static limit where inertia effects are neglected (this restriction will be removed below). The distribution  P_c(x)  can be characterized by its average value  x_s  and standard deviation  s_s;  a typical example is the Gaussian distribution  P_c(x)=G(x; x_s,s_s) = [1/s √(2p)] exp[-(x-x_s)²/4s²].

To describe the evolution of the model, we introduce the distribution  Q(x; X)  of the stretchings  x_i  when the bottom of the sliding block is at a position  X.  Let all asperities be initially relaxed or weakly stressed, e.g., let the distribution  Q(x; 0)=Q_ini(x)  be the Gaussian,  Q_ini(x)=G(x; x_ini,s_ini)  with  x_ini=0  and  s_ini<<s_s.  Now, let us adiabatically increase the displacement  X  of the bottom of the top sliding block while the base (the bottom substrate) remains fixed. The sum of the elastic forces exerted on the bottom of the block by the stretched asperities makes up the friction force

The evolution of the system, deduced from the numerical simulation of the EQ model, is presented in Fig.eq2a. It shows that, in the long term, the initial distribution approaches a stationary distribution  Q_s(x)  and the total force  F  becomes independent on  X.  The final distribution is independent of the initial one (an elegant mathematical proof of this statement was given in Z.Farkas, S.R.Dahmen, D.E.Wolf, J.Stat.Mech.: Theory and Experiment P06015 (2005); cond-mat/0502644; the authors considered a simplified version of the EQ model, assuming that every contact keeps its own threshold value  f_si  unchanged after breaking/reforming); the statement is valid for any distribution  P_c(x)  except for the singular case of  P_c(x)=δ(x-x_s).

Master equation. Rather than studying the evolution of the distribution  Q(x; X)  by a simulation of the EQ model, it is possible to describe it analytically. Let us consider a small displacement  DX>0  of the bottom of the sliding block (for  DX<0  see below). It induces a variation of the stretching  x_i  of the asperities which has the same value  DX  for all asperities if the deformation of the bottom surface of the block can be neglected. The displacement  X  leads to three kinds of changes in the distribution  Q(x; X):  first, there is a shift due to the global increase of the stretching of the asperities, second, some contacts break because the stretching exceeds the maximum that they can stand, and third, those broken contacts form again, at a lower stretching, after a slip at the scale of the asperities, which locally reduces the tension within the corresponding asperities. These three contributions can be written as a master equation for  Q(x; X):

Once the distribution  Q(x; X)  is known, we can calculate the friction force  F(X)  using Eq.(eq1). The static friction force corresponds to the maximum of  F(X),  i.e.,  F_s=F(X_s),  where  X_s  is a solution of the equation  F'(X)≡dF(X)/dX=0.  In order to simplify the further consideration, let us assume that  P_c(x)=0  for  x≤0  (this agrees with its physical meaning because, if  x<0,  a positive variation  DX  actually reduces the absolute value of the force on a contact, which does not cause its breaking). Also we choose  R(x)=δ(x),  i.e., we assume that a broken contact sticks again only after a complete relaxation.

The steady-state, or smooth-sliding solution, i.e. the solution of Eq.(eq6) which does not depend on  X,  can easily be found. It can be expressed as

where  Q(x)  is the Heaviside step function,  E_P(x)=e^-U(x),  U(x)=∫₀^x dx P(x),  and  C=∫₀^∞ dx E_P(x).  Note also the following useful relationships:  U'(x)=P(x),  E_P'(x)=-E_P(x)P(x),  J_c(x)=E_P(x),  P_c(x)=P(x)E_P(x)  for  x > 0, and  E_P(x)=1  for  x≤0.

In the general case, let the distribution  P_c(x)  be of bell-like shape with the maximum at  x_s  and the width  s_s.  When  X  shifts for the distance  x_s,  due to the breaking and reforming of contacts with a lower stretching, an initially peaked distribution  Q(x; X)  broadens by the value ~s_s (Fig.eq2). Therefore, any initial distribution tends to the stationary one as  |Q(x; X) - Q_s(x)| µ exp(-X/X^*),  where  X^*~x_s²/s_s.

Some examples

Example 1: Steplike distribution

In particular, if  x_s=0,  then  Q_s(x)=Q(x)pe^-px  and  F=N‹k_i›/p,  while for the case of  x_s>0  and  p®∞  we obtain that  Q_s(x)=x_s^-1  within the interval  0≤x≤x_s  and  0  outside it, so that  F=¹/₂N‹k_i›x_s.

Example 2: Rectangular distribution

Another simple example which admits the exact solution, is the case of rectangular  P_c(x)  distribution shown in Fig.eq25:

and the "rate"  P(x)  is given by

Example 3: The singular case

In a general case the solution of master equation always approaches the steady-state solution corresponded to smooth sliding. However, there is one exception from this general scenario, when the model admits a periodic solution and  F(X)≠Const even in the  X®∞   limit.  Namely, this is the singular case when all contacts are identical, i.e., all contacts are characterized by the same static threshold  x_s,  P_c(x)=δ(x-x_s),  or  P(x)=0  for  x<x_s  and  P(x)=∞  for  x≥x_s.  In this case we can find (guess) the steady-state solution of the master equation (e.g., it may be checked by direct substitution):

The static friction force takes the minimal value  F_s=¹/₂Nkx_s  for the uniform initial distribution  Q_ini(x)=x_s^-1,  when  F(X)  does not depend on  X,  and the maximal value  F_s=Nkx_s  for the delta-function initial distribution  Q_ini(x)=δ(x-x₀)  with some  0≤x₀<x_s,  when the function  F(X)  has sawtooth shape changing from  0  to  F_s.

Example 4: Two delta-functions

The periodic solution described above exists only for the distribution  P_c(x)  with the single threshold. If the contacts are characterized by more than one threshold value, for example, if one part of contacts has the threshold  x_s1  and another part, the threshold  x_s2≠x_s1  [i.e., P_c(x)  is described by a sum of two delta-functions], then the system will always approach the stationary steady state. This is demonstrated in Fig.eq26, where we compare the system evolution in cases of one-peaked and two-peaked  P_c(x)  distributions. Notice, however, that this statement is valid only for the infinite set of contacts (the number of contacts with each threshold must be infinite) and cannot be applied for the system where, e.g., two tips move over a surface.

Figure eq26: Evolution of the model for two examples of the  Pc(x)  distribution: Left coulomb is for one-peaked distribution of threshold,  Pc(x)=G(x; x1,s)  with  x1=1  and  s=10-3  (i.e., close to the delta-function distribution), and right coulomb is for two-peaked distribution of threshold,  Pc(x)=1/2 [G(x; 0.95x1,s) + G(x; 1.05x1,s)]  (i.e., close to the distribution with two delta-functions). Solid curves show the distribution  Q(x; X)  for incrementally increasing values of  X  with the increment  DX»1,  dotted curves show  Pc(x)  (for  X=0  only). The bottom raw shows the corresponding dependences  F(X)  for  ‹ki›=1 and  K=∞.  The initial distribution  Qini(x)  is Gaussian with  xini=0.1  and  sini=0.025.

Example 5: Friction loop

Above we assumed that the slider moves continuously to the right, or  DX>0.  In the case when the top substrate moves to the left,  DX<0,  equations (eq3), (eq5) and (eq6) must be modified in the following manner:

where for the symmetric case (the forward-backward symmetry) we have to put  P_b(x)=P(-x)  and  R_b(x)=R(-x).

The reason for such a behavior is the irreversibility of the master equation. Equation (eq6) describes the "forward" dynamics when  X  increases, while Eq.(eq67) describes the "backward" dynamics when  X  decreases. Indeed, if the force  f_i  on a given contact approaches and overcomes  f_si,  the contact breaks; but if we now reverse the motion direction, this contact cannot jump back to the value  f_i»f_si;  instead,  f_i  will decrease until it reaches the value  f_i»-f_si .

In a typical tribological experiment, the top substrate moves periodically forward/backward, and in the result, the so-called "friction loop" is obtained. The same loop can easily be calculated with the ME approach as demonstrated in Fig.eq27.

Generalizations

Delay in contact formation

Figure eq28: The dependences  F(X)  for different delay times. The distribution  Pc(x)  is Gaussian with  xs=1  and  ss=0.2,  the initial distribution  Qini(x)  is Gaussian with  xini=0.1  and  sini=0.025

Aging of junctions

Now let us take into account the aging of contacts. Experiments as well as MD simulations indicate that the static friction force grows with the time of stationary contact (the waiting time  t_w,  i.e., the duration of static contact prior to sliding). Thus, if the newborn contacts are characterized by a distribution  P_ci(x)  with the average  x_si  and dispersion  s_si,  then typically  x_s  grows while  s_s  decreases with time, and at  t®∞  the distribution  P_c(x, t)  approaches a distribution  P_cf(x)  with  x_sf >x_si  and  s_sf <s_si.  If we assume that the evolution of  P_c(x)  corresponds to a stochastic process, then it should be described by the Smoluchowsky equation

the "diffusion" parameter  D  describes the rate of aging,  B(x)=dU(x)/dx,  and the "potential"  U(x)  is defined by the final distribution,  P_cf (x) µ exp[-U(x)]  so that

Then, the equation  ¶P_c/¶t = DL_xP_c  naturally leads to the growth of the average static threshold from  x_si  to  x_sf  with the time of stationary contact, as widely assumed in earthquakelike models of friction (e.g., see Eq.(Fs) in previous Chapter). At the same time, the contacts continuously break and reborn when the substrate moves, as described by third and forth terms in Eq.(eq6). Combining both the processes together, we come to the system of three equations:

where  D_X=D/v,  and  v=dX(t)/dt  is the driving velocity. Thus, in the case of high driving  (v®∞,  D_X®0)  we obtain the previous behavior with  P_c(x)=P_ci(x),  in the case of low driving  (v®0,  D_X®∞)  we again observe the same type of behavior but with  P_c(x)=P_cf (x),  and in the case of  0 < D_X < ∞  we have an interplay of two processes – the aging which moves  P_c(x)  to  P_cf (x),  and the breaking/reborn process which returns  P_c(x)  to  P_ci(x)  (see Fig.eq8). 

In the result, the friction force  F  now depends on the sliding velocity as shown in Fig.eq9.

 Below we show that, if the sliding block is not rigid,  K<∞,  the system will demonstrate either stick-slip or smooth sliding depending on the distribution  P_c(x). In this case the effect of contact's aging must lead to the transition from stick-slip to smooth sliding with the increase of the sliding velocity. The approach predicts, however, that this transition at  v=v_c  is abrupt (first-order) contrary to what we observed in the EQ model with identical contacts. Moreover, the system should exhibit hysteresis: if the velocity decreases starting from the smooth-sliding regime, the latter will survive till much lower velocities v << v_c.

Note: in the case of contact of rough surfaces, the physical mechanisms of contact aging, according to Baumberger and Caroli (2006), are the following: the first (and more important) mechanism is due to geometrical aging, or the increase of contact area at the asperity, and the second, due to structural aging, or restructuring of the contact. The value of the parameter  D  in Eq.(eq30) may be estimated from experiments which show that in most cases the average static threshold  m_s=F_s/F_load  grows logarithmically with the waiting time  t_w,  dm_s/d(ln t_w) »10^-2.

Role of temperature

and the function  H(x)  in Eq.(eq39) is determined by

On the other hand, for "stiff", or "strong" contacts, when  DE(0;x_s) >> k_BT,  we have to substitute in Eq.(eq40) for the activation energy the expression  DE(x; x_s) = DE(0; x_s) (1-x/x_s)^3/2,  and also to renormalize the prefactor  w,  w ® w (1-x/x_s)^1/2.  In this case the function  H(x)  is to be calculated as

The contribution (eq42) or (eq45) to the rate  P_T(x)  leads to appearance of a low-x tail. Its height grows with temperature as well as with decreasing of the driving velocity  v  as demonstrated in Fig.eq15. The increase of temperature leads to a shift of the distribution  Q(x)  to lower values, so that the friction force decreases when  T  grows. The effect is the larger, the lower is the sliding velocity as shown in Fig.eq16. In the limit  v®0,  all contacts will finally break if  T>0,  so that  Q_s(x)®δ(x)  and  F_k®0;  in this limit we have "smooth sliding" corresponded to creep motion of contacts.

Figure eq15: The rate  PT(x)  for soft contacts at different temperatures  T=0 to 1 for a fixed velocity  w/kv=1, and (inset) at different velocities,  w/kv=0 to 10, for a fixed temperature  T=0.3.  Dotted curve shows the distribution  Pc(x)  (Gaussian with  xs=1  and  ss=0.05).	Figure eq16: The steady-state distribution  Qs(x)  for  T=0.3  and different velocities  (w/kv=0.3, 1, 3 and 10) for soft contacts. The distribution  Pc(x)  is Gaussian with  xs=1  and  ss=0.05.

Figure eq17: (a) Dependence of the kinetic friction force  Fk  in the steady state on temperature for different velocities  v=0.01, 0.1, 1 and 10.  (b) Fk(v)  for different temperatures  T=0.01, 0.1, 0.3 and 1; inset shows the same in log-log scale. Solid curves are for soft contacts, and broken curves, for stiff contacts. The distribution  Pc(x)  is Gaussian with  xs=1  and  s=0.05;  w/k=1.	Figure eq19: The dependences  F(X)  for short  X.  (a)  F(X)  for different temperatures  T=0, 0.03, 0.1, 0.3 and 1 at the fixed velocity  w/kv=1;  (b)  F(X)  for different velocities  v,  w/kv=0, 0.3, 1, 3 and 10 at the fixed temperature  T=0.3.

Nonzero temperature influences on the dynamics of approaching to the steady state (see Fig.eq18 and Fig.eq19). The higher is temperature and/or the lower is velocity, the lower is the static friction force  F_s  determined by the first maximum of  F(X),  and the faster  F(X)  approaches to the steady-state smooth sliding. Also it is important to consider the first cycle of the  F(X)  dependence, which defines the lowest value of  F'(X).  The higher is temperature, the lower is the extremum of  F'(X)  (see Fig.eq19a), so that the larger is the interval of model parameters where the smooth sliding regime operates. At the same time, the higher is driving velocity (Fig.eq19b), the smaller is the lowest value of  F'(X),  so that the narrower is the interval of model parameters where the smooth sliding regime operates (see below). Thus, we come to a surprising conclusion that at  T>0  the decreasing of pulling velocity may lead to the transition from stick-slip to smooth sliding, i.e., the scenario just opposite to the conventional one.

Figure eq18: The dependences  F(X)  for  T=0.3  at different driving velocities  v,  w/kv=0, 0.3, 1, 3 and 10. The distribution  Pc(x)  is Gaussian with  xs=1  and  ss=0.05,  the initial distribution  Qini(x)  is Gaussian with  xini=0.1  and  sini=0.025

The described behavior of  F_k  on  T  and  v  qualitatively agrees with the tip-based experiments. Surprisingly, the dependences obtained within the ME approach, perfectly agree with the experimental ones for a "model" tribological system, where the kinetic friction of the driven lattice of quantized magnetic vortexes in high-temperature cuprate superconductors was studied (e.g., compare inset of Fig.eq17b with Fig.3 of Maeda et.al. Int.J.Mod.Phys. B 19 (2005) 463).

Distribution of thresholds

Let us discuss now physical origins of the distribution  P_c(x).  First of all, it can be coupled with the distribution  P_cf (f_s)  of (static) friction force thresholds of the contacts. If a given contact has an area  A,  then it is characterized by the force threshold  f_sµA  and the elastic constant  k~rc²√A. The displacement threshold for the given contact is  x_s=f_s/k, so that  f_s µ x_s²,  or  df_s/dx_s µ x_s.  Thus, using the relationship  P_c(x_s) dx_s = P_cf (f_s) df_s,  we obtain

Rough surfaces

For a contact of two hard rough surfaces (the multi-contact interface, or MCI in notations due to Baumberger and Caroli), the problem reduces to the statistics of asperities. Let the rough surface is characterized by hills of heights  {h_i}  distributed with a probability  P_h(h).  Following the Greenwood and Williamson model of the interface, let us assume that all hills have the spherical shape of the same radius of curvature  r.  When this surface is pressed with another rigid flat surface, which takes a position at the level  h₀,  then the hills of heights  h>h₀  will form contacts, or asperities. If the contacts are elastic (the Hertz contacts), then the contact of height  h  has the compression  (h-h₀),  its area is  pr(h-h₀),  and it bears the normal force  f_l(h)»(4p/3)E^*r^1/2(h-h₀)^3/2,  where  E^*  is the effective Young modulus (if both substrates have the same Young modulus  E,  then  E^*=E/2(1-n²),  where  n  is the Poisson modulus). If we assume that the shear static threshold for the given contact is proportional to the load force,  f_s(h)»mf_l(h),  or

where  m<1  is a constant, then the distribution of static thresholds  P_cf (f_s)  can be coupled with the distribution of asperities heights  P_h(h)  with the help of relation  P_cf (f_s) df_s µ P_h(h) dh,  or  P_cf (f_s) µ (dh/df_s) P_h[h(f_s)],  where  dh/df_s µ f_s^-1/3  according to Eq.(eq48).

For a strong load, when the local stress exceeds the yield threshold  Y,  the contacts begin to deform plastically. When all contacts are plastic, the local pressure on contacts is  p_load=H,  where  H  is the hardness  [ H»3Y  for the spherical geometry of asperities; for metals  H»(10^-3-10^-2)E ~ 10⁹ Pa ]. Then the normal force at the contact is  f_l (h)»pr(h-h₀)H.  Assuming again that  f_s (h)»mf_l(h),  we obtain

so that  P_cf (f_s) µ P_h[h(f_s)].  For example, if the distribution of heights is exponential,  P_h(h)=h₀^-1exp(-h/h₀) Q(h),  where  h₀  is the average height, then for the elastic contacts we obtain (f, x > 0)

while for the plastic contacts

The distribution  P_c(x)  for the elastic and plastic contacts is presented in Fig.eq20. The force  F(X)  (almost) monotonically increases with  X,  approaching the kinetic value  F_k. Thus, in the case of contact of rough surfaces, a relatively large concentration of low-threshold contacts prevents from appearance of the stick-slip motion even for a very soft substrates. The kinetic friction force  F_k  depends on the parameter  B  according to the power law (see Fig.eq21),  F_k µ B^-3/4 µ h₀^3/4  for elastic contacts, and  F_k µ B^-1/2 µ h₀^1/2  for plastic contacts.

Note: Typical values for rough surfaces are of the order  r ~ 10 - 100 mm,  h₀ ~ 1 mm, so that the average size of the contact is  a » rh₀^1/2 ~ 3 - 10 mm. Are the contacts in the plastic or elastic regime, depends on the dimensionless parameter  y=(E/Y)(h₀/r)^1/2:  the former case operates for  y>>1  (as is typical for metals), while the latter, for  y<1  (it corresponds, e.g., to the case of rubber friction); the polymeric glasses belong to an intermediate situation,  y>1,  where only a fraction of contacts is plastic.

Dry contact

Next, let us consider the dry contact of two flat surfaces. In an exotic case when both surfaces have the ideal crystalline structure, we come to the singular case of delta-function distribution  P_c(f).   Such a situation, however, is exceptional. A real surface always consists of domains, which are characterized by different crystalline orientation and even, may be, different structure. MD simulations show a large variation of the friction with relative orientation of the two bare substrates. In order to estimate a shape of the function  P_c(x)  emerging due to domain structure of substrates, let us consider a simple model. Let a domain of the top rigid substrate be crystalline (e.g., of triangular symmetry) with the lattice constant  a  and consists of  N  atoms, while the bottom substrate be also rigid crystalline (e.g., with square symmetry) so that it produces a periodic potential for the motion of the top substrate. For some values of  a  and  N,  the activation energy for domain displacement is high, while for other values, the activation energy is small or even may vanish. Besides, if the domain is rotated on the misfit angle  f,  the activation energy should achieve minima at some angles.

For a fixed misfit angle  f,  one can calculate the total potential energy of the domain  U(X,Y),  where  X  and  Y  are the center of mass coordinates of the domain. The extrema of  U(X,Y)  are determined by the equations  ¶U/¶X=0, ¶U/¶Y=0.  An extremum may correspond either to a minimum  U_m  or to a saddle point  U_s;  then the activation energy is given by  e_a=U_s-U_m.  Assuming that  f_s~e_a/a µ e_a,  we can find the threshold force  f_s  as a function of the misfit angle  f  and the domain size  N.  Then, calculating a histogram of the function  e_a(f),  we obtain the distribution  P_cf (f_s)  if all domains have the same size  N  and all angles are equally presented. Averaging it over different domain sizes  N,  e.g., with a weight function  w(N)=e^-N/‹N›,  where  ‹N›  is the average domain size, we obtain the distribution  P_c(x)  shown in Fig.eq23 (note that it is similar to that of Fig.eq20 characteristic for rough surfaces).

Thus, in this system one also has to expect smooth sliding generally. However, in the estimation we assumed that the domains are rigid, while real substrates are deformable. Also, we supposed that all angles are equally presented, while some angles should have preference due to their lower potential energy. Both these factors should lead to the increase of the threshold values.

Lubricated surfaces: Lifshitz-Slözov coalescence

The dry-friction system considered above, is also exceptional. In a real system, almost always there is a lubricant between the sliding surfaces (called "the third bodies" by tribologists) – either a specially chosen lubricant film, or a grease (oil), or dust, or wear debris produced by sliding, or water or/and a thin layer of hydrocarbons adsorbed from air, etc.

In the conventional melting-freezing mechanism of friction, the lubricant is melted during slip, and solidifies when the motion stops. The solidification process can be described by the Lifshitz-Slözov theory. At the beginning, grains of the solid phase emerge within the liquid lubricant film. Then the grains grow in size according to the law  ‹r› µ t^1/3.  The distribution of grains sizes is described as follows: the number of grains with the radius from  r  to  r+Dr  is equal to  P_LS(r/‹r›)Dr/‹r›,  where the function  P_LS(u)  is zero for  u≥3/2  (i.e., the maximum size of the grains is 1.5‹r›),  while for lower sizes it is given by

Due to coalescence of grains, the total number of grains decreases with time as  N(t) µ t^-1.  When the size of a grain exceeds the film thickness  h  (that will occur when  ‹r›(t) ≥ h/3),  such grains will pin the surfaces. Using the relationship  P_cf (f) df µ P_LS(r/‹r›) dr/‹r›,  we obtain that  P_cf (f) µ (dr/df) P_LS(r/‹r›)/‹r›.  Because one grain gives the static threshold proportional to the contacting area,  f_s µ p(r²-h²/4)  for  r>h/2,  we have  dr/df µ f^-1/2.  Combining all things together, we obtain that

where  r(t)=2‹r›(t)/h=at^1/3  (the pinning begins when  r>2/3),  and  B  is determined by the system parameters (elasticity of the contacts, proportionality between the threshold  f_s  and the contact area for a contact/domain/grain, and the thickness of the lubricant film). The distribution (eq53) is shown in Fig.eq24; it should lead to the conventional tribological behavior: the stick-slip motion at low driving velocity and smooth sliding at high velocities.

In a general case, we also have to take into account that the lubricant film may consist of grains (domains) of different orientation or even different structure. Indeed, the proportionality coefficient in the relation  f_s µ p(r²-h²/4)  used above, should depend on the misfit angle between the lubricant domain and the substrate, so that the distribution  P_LS(u)  introduced above, additionally should depend on the misfit angle  f,  P_LS(u; f).  Thus, if there are grains with different orientation, distributed according to a function  R_f(f), then

Stick-slip kinetics

Above we described the rigid motion of the top substrate, when the ME solution always leads to smooth sliding. In order to describe a real dynamics of a tribosystem, we have to take into account the following three issues: (i) the slider inertia, (ii) the elastic instability, and (iii) the delay in contact formation.

Motion equation

Firstly, we have to involve the motion equation for the slider; let it be modeled as a rigid entity (a nonrigid slider will be considered below) of mass  M  (see figure to the right). The slider is driven through a spring of the elastic constant  K  with the velocity  v_d ,  so that the (experimentally measured) spring force is  F_d=K (X_d -X),  X_d=v_dt.  Also the slider experiences the force  -F  from all contacts as described above with the ME or EQ approach. Every contact at the interface moves rigidly with the slider so long as the total force f_i<f_si, where f_si is an upper threshold value. Above that threshold the contact detaches from the slider and slips relative to it for a delay time  t_d  after which it stops and attaches again to the slider; at this moment its stretching is determined from the condition that the total force  f_i  on the contact is equal to the "backward" force  f_b~0.  The force  F  from the interface on the slider is the sum of the forces from the pinned contacts,  f_i^(sub) =k_ix_i  (recall that  x_i  is the stretchings of the contact and  k_i  is its elastic constant), plus additionally the sum of drag forces from the detached/sliding contacts,  f_i^(drag) =-m_ih_i [v-dx_i(t)/dt ],  where  v(t)=dX(t)/dt  is the slider velocity,  m_i  is the mass of the contact and  h_i  is the corresponding damping coefficient (e.g., the "viscosity" of the molten lubricant).

Because the slider is rigid, at any abrupt change of its motion, e.g., at the onset or stopping of motion, it will undergo vibrations with the setup frequency  W_S=(K/M )^1/2,  which cannot decay and therefore will disturb the results. In a real system, these oscillations will be attenuated because of internal friction inside the setup (i.e., due to phonon damping in the slider, when its bottom surface oscillates relative the top one). To incorporate such a damping, we introduce a viscous damping coefficient  h_S  for the slider motion relative its average velocity, so that finally the slider motion equation is

Model parameters. In simulation results presented below we used the EQ algorithm (the ME approach leads to the same results). For the sake of quantitative results, we extracted model parameters from the experiment of Klein (J. Klein, Phys.Rev.Lett. 98 (2007) 056101), where the OMCTS (octamethylcyclotetrasiloxane) lubricant film of the thickness  h=3.5×10^-9 m  (four molecular layers) between two atomically smooth mica surfaces was studied with the surface force balance (SFB) technique. Thus, we took  M=1.47 g and  K=97 N/m, that gives the setup frequency  W_S=257 s^-1  and the period  t_S=2p/W_S=0.0245 sec. Also we used  A=10^-10 m²  for the total contact area and  F_s=18 mN  for the static friction force. The damping coefficient for the drag force may be found from the bulk viscosity of OMCTS; this gives  h=2×10¹¹ s^-1.  Assuming that there are  N  contacts at the sliding interface  (N=4080  in simulation), we can find the parameters for individual contacts:  A_i=A/N  for the contact area,  a_i=√(A_i/p)  for its radius,  m=r_OMCTSA_ih  for its mass,  f_s=F_s/N  for its static threshold, and  a=(A/N)^1/2  for the average distance between neighboring contacts. The contact elastic constant can be found as  k»rc²a_i;  in simulation we took  kN = 2000 N/m (note that it should be  Nk>>K  to have stick-slip). The static thresholds  f_si  for individual contacts take random values from the Gaussian distribution with the mean  f_s  and standard deviation  Df_s;  other contacts' parameters are coupled with  f_si  by  m_i=mf_si/f_s,  k_i=k (f_si /f_s)^1/2  and  h_i=h  (when a contact reborn, its parameters are refreshed). Finally, we used  h_S=0.2W_S ,  v_d=0.1 mm/s,  f_b=0.1f_s,  Df_s=0.01f_s,  and  t_d=5×10^-4 sec  (all the parameters were also varied in wide intervals).

Elastic instability

Secondly, the master equation allows us to compute the force  F(X)  when the bottom of the solid block is displaced by  X,  but actually we don't control  X.  The displacement is caused by a shearing force  F_d=K (X_d -X)  applied on the top of the solid block which displaces the top surface by  X_d.  The total force applied to the bottom of the sliding block, which determines its displacement  X,  is the sum of the applied force and the force from the interface,   F_tot=K (X_d-X) - F(X).  This force can be viewed as deriving from the potential

A necessary condition for smooth sliding is that  X_d  and  X  grow together with  X_d-X=B,  where  B  is a constant that measures the shear strain of the solid block during the sliding. It is determined by the condition  ¶V_eff/¶(X_d-X)=-¶V_eff/¶X=0, which simply means that the total force on the interface vanishes. Smooth sliding also requires this state be stable,

If we start from relaxed asperities, in the early stage of the motion  F(X)  is a growing function of  X  (see Fig.stick3), and then it passes by a maximum when some contacts start to break and reform at lower asperity stress. As a result  F'(X)=dF/dX  becomes negative. Depending on the value of  K  two situations are possible. For large  K  (K>K^*= - max F'(X),  stiff block)  F'(X)  never falls below  -K  and the smooth sliding is a stable steady state. For small  K  (K<K^*,  soft block)  F'(X)  can become smaller than  -K  so that the stability condition (eq10) is no longer valid (this is the well-known elastic instability, e.g., see T.Baumberger and C.Caroli, Advances in Physics 55 (2006) 279). The instability may cause  X_d-X  to change abruptly by a breaking of all contacts and a quick slip of the block before the contacts reform with relaxed asperities. If the process repeats itself, we have the stick-slip motion. The master equation, which gives  F(X),  can be used to compute the period of the stick-slip. However, the existence of a stick slip is not only determined by the stiffness  K  of the block, but also by the distribution  P_c(x)  of the static thresholds and the time of contact reforming, as is described below.

Thus, if the slider is soft enough,  K<K^*,  its motion becomes unstable at the point  X_c,  where  X_c  is the (lowest) solution of the equation  F'(X)=-K  (see Fig.stick3). Let this occur at the moment  t_c.  After the unstable point,  t>t_c,  the motion equation in terms of  DX=X-X_c  and  Dt=t-t_c  is

with the initial condition  DX=0  and  dX/dt=v_d  at  Dt=0.  In the EQ model with the distribution of static thresholds, the function  F(X)  near the point  X_c  can be approximated as  F(X) » F(X_c) + F'(X_c) DX + (1/2)F''(X_c) DX² + (1/6)F'''(X_c) DX³,  so that

Introducing the dimensionless variables  x=K(X-X_c)/F_s  and  t=W_S(t-t_c),  we obtain the following equation for the slider coordinate:

where

Solution of Eq.(stick5) should be found with the initial condition  x(0)=0  and  dx(t)/dt|_t=0 =D₁.  At short times,  t<<1,  the solution has the form  x(t) » D₁t + (1/6)D₁t³.  Thus, as  X(t)  grows at  t>t_c,  the spring force quickly drops during a slip time  t_s»aW_S^-1  in the form

where  a=(6F_sW_S /Kv_d)^1/3  (a»2.21 for the chosen set of parameters). Numerical solution of Eq.(stick5) is shown in Fig.stick4.

On the other hand, if the whole lubricant film melts instantly at  t=t_c,  then for  t>t_c  we have  F(X_c+DX) =Mh d(DX)/dt,  and the motion equation in dimensionless variables takes the form

where  D₀=-F_d(X_c)/F_s.  Its solution for  t <<1  is  x(t) » D₁t + Bt²,  so that the spring force drops as

with  B=(1/2)(D₀-hD₁/W_S).  This gives the slip time  t_s»b W_S^-1  with  b=[2/(D₀-hD₁/W_S)]^1/2.  Numerical solution of Eq.(stick6) is shown in Fig.stick4. One can see that in this case the slip time is much shorter; however, using artificially huge values for  h,  one may make  t_s  as large as desired (see discussion).

Slip details

A typical example of stick-slip motion of the EQ model for the chosen set of parameters is presented in Fig.stick2. It is very similar to that observed experimentally, including the large initial stick spike and subsequent stick spikes of smaller amplitude (e.g., compare Fig.stick2 with Figs. 1 and 2 in Klein's paper).

Details of the slip are shown in Fig.stick5. During the slip time  t_s,  the slider velocity grows up to a value  v~W_SF_s /K.  As the slider accelerates, the contacts break, the number of detached contacts grows, and the force  F(t)  on the slider from the contacts drops from  F_d(t_c)  to zero. After the delay time  t_d,  the contacts reborn, and  F(t)  grows as  F(t)=Ct  with the rate  C=Nkvt  back to the value  ~F_d(t).  Therefore, the force  F  oscillates with a period  t_B»F_s /Nkv»K/NkW_S.  Such oscillations exist provided  t_B<t_s,  or  K<Nk.

Thirdly, let us consider the role of delay time. One may think that the elastic instability will result in stick-slip motion. However, the condition  K<K^*  is the necessary but not sufficient one for stick-slip; the second necessary condition is a nonzero delay time,  t_d>0.  This is clearly demonstrated in Fig.stick6: if  t_d=0,  the system still approaches the smooth sliding.

The subsequent (after the slip) motion of the slider depends on a value of the parameter  t_d.  If  t_d>W_S^-1,  the spring force drops to negative values and then oscillates around zero with the frequency  W_S  during the time  t_d  (Fig.stick6c). Otherwise, if  t_d<<W_S^-1,  the spring force drops to a value higher than  F_b  and then begins to grow, oscillating with a high frequency  W_L=(Nk/M)^1/2=1.17×10³ sec^-1  (so that  2p/W_L=5.39×10^-3 sec; see Fig.stick6a and Fig.stick6b). In all cases, the ringing vibrations decay as e^-h_St.

If we neglect ringing oscillations, then the dependence  F_d  versus  v_dt  is "universal". Thus, if the velocity is so small that t_ss >> h_S^-1, i.e., the oscillations are completely damped during the time per stick-slip cycle  t_ss,  the amplitude of stick-slip variation remains the same. But if the driving velocity is so high that  t_ss<h_S^-1,  then the oscillations will disturb the system dynamics and may lead to smooth sliding.

Role of threshold's dispersion. The system behavior (either stick-slip or smooth sliding) is controlled by the dispersion  Df_s.  Indeed,

Thus, if  Df_s  is so small that  K^*>K,  then the motion corresponds to stick-slip; otherwise the smooth sliding regime is achieved. In the stick-slip regime, an increase of  Df_s  leads to the decrease of the period  t_ss  of stick-slips as demonstrated in Fig.stick8. The time  t_s  of the  F(t)  drop during slip also increases with  Df_s,  but this effect is rather small.

Finally, note that the ratio  Df_s/f_s  should decrease with the time of stationary contact due to aging of contacts; namely this aging is responsible for the transition from stick-slip to smooth sliding with the increase of driving velocity.

Role of slider elasticity. Above we have assumed that the slider is a rigid body; in a real setup, however, both the slider mass and its elasticity are distributed through the slider. One may think that in this case, only the most bottom atomic layer of the slider starts to move at the onset of slip. Therefore, a characteristic frequency will be determined by the mass  M_l  of the layer and its elastic constant  K_l;  this could lead to a much higher (atomic scale) frequency. A similar question – which mass, either the total mass  M  or the layer mass  M_l,  defines the time scale of the system? – appeared in the problem of minimal velocity for the atomic-scale smooth sliding. As was proven, when the slider velocity decreases, first the most bottom layer stops; then the stopping wave emerges and takes away the kinetic energy of the slider. In a result, the time scale of that problem is determined by the layer mass  M_l,  and the minimal velocity is of atomic-scale order (m/s).

Now, however, the situation is just opposite. To show this, let us consider the slider as consisting of  N_l  layers, each of the mass  M_l=M/N_l,  elastically coupled by springs of the elastic constant  K_l=KN_l  (see inset in Fig.stick19). Indeed, if we fix the bottom layer and apply a force  F  to the top layer, then the latter will shift on the displacement  DX=∑_l=1^N_l DX_l,  where  DX_l =F/K_l,  so that  DX=F/K  as before. Now let the top layer be driven with the velocity  v_d,  while the bottom layer be in frictional contact with the bottom substrate. The dependence of the elastic force in the slider on time, obtained with simulation for two different values of  N_l,  is shown in Fig.stick19 (note that now we may not use the artificial damping  h_S,  because internal degrees of freedom are included). As seen, the slip kinetics is almost independent on the number of layers  N_l  and is determined by the minimal slider mechanical frequency  W_S.

The frequency  W_S  can be found with the help of elastic theory. Let the slider have a cylinder shape of height  L  and radius  r,  and is characterized by the section  S=pr²,  inertial momentum  I=pr⁴/4,  mass density  r  and Young modulus  E.  If the cylinder foot is fixed and a force  F  is applied to its top, the latter will be shifted on the distance  DX=FL³/3EI  (the problem of bending pivot, see Sec.20, example 3 in the Landau and Lifshitz (LL) textbook). Thus, the effective elastic constant of the slider is  K=3EI/L³.  The minimal frequency of bending vibration of the pivot with one fixed end and one free end, is given by W_S=(3.52/L²)(EI/rS)^1/2  (see Sec. 25, example 6 in LL). Taking  M=rSL,  we finally obtain  W_S»2.03√(K/M).

Interacting contacts

Finally, let us consider the role of interaction between the contacts. A concerted, or synchronized breaking (triggering) of contacts may be studied numerically only within the earthquakelike model, it cannot be included accurately in the ME approach (although it may be incorporated indirectly in a mean-field fashion by renormalization of the distributions  P_c(x)  and  R(x)).

The elastic interaction between contacts separated by a distance  r  is described by the pairwise potential  V(r)=g/r³,  where  g  is a constant. Using  ¶V/¶x=V'(r) ¶r/¶x,  V'(r)=dV/dr=-3g/r⁴  and  ¶r/¶x=x/r,  we obtain that the force acting on the  ith contact from its nearest neighbors (NN), is equal to

The interaction becomes important, when  g/a³~f_sa,  where  a  is the average distance between NN contacts; therefore, it is interesting to check, how the system behavior changes with the dimensionless parameter  x=g/(f_sa⁴).

The interaction between the contacts works roughly in the same way as the dispersion  Df_s:  the stronger is the interaction, the wider is the range of model parameters where stick-slip operates. To demonstrate this effect more clearly, we calculated the system kinetics with increasing interaction  (x=0 to 0.3)  for softer contacts  (kN=200 N/m)  and a larger dispersion  Df_s/f_s=0.3.  The results are presented in Fig.stick11. For these parameters, the system quickly goes to smooth sliding for noninteracting contacts (Fig.stick11a), but demonstrates stick-slip for a strong interaction  x=0.3  (Fig.stick11d).

Details of the slip are shown in Fig.stick12: without interaction, contact's detaching is sequential, while for a strong interaction, all contacts detach simultaneously, and the force  F  abruptly drops to zero. This is clearly demonstrated in Fig.stick12 (compare panels (a) and (d)).

The interaction may lead to avalanches as was observed in the EQ simulation: for  x=0  the contacts break sequentially, one after another, while for  x>0  one contact may stimulate the breaking of nearest contacts, and the average size of avalanche increases with  x  as demonstrated in Fig.stick13. For a large strength of interaction,  x=0.3,  the avalanche occupies the whole system (Fig.stick13d).

Discussion

In tribological community there is a common opinion that the viscosity of a thin confined film is many orders of magnitude higher than the bulk viscosity (P.A.Thompson, M.O.Robbins, G.S.Grest, Israel J.Chem. 35 (1995) 93;  H.-W.Hu, G.A.Carson, S.Granick, Phys.Rev.Lett. 66 (1991) 2758). For example, considering the drop of the spring force during slip, Klein (J.Klein, Phys.Rev.Lett. 98 (2007) 056101) came to the conclusion of "high viscosity" of the molten lubricant film,  h »10⁴h_OMCTS. So large viscosity of a thin film, however, has no reasonable explanation. The cases of a single or two molecular layers confined between two planar surfaces are exceptional in the sense that properties of such films are far from their bulk ones. But already a three-layer film, if it is melted (liquidized) during slips, should exhibit properties not very different from the bulk ones, as follows from almost all MD simulations.

In our model with a distribution of static thresholds, the time  t_s  of  F(t)  drop is almost solely determined by the setup frequency  W_S.  This time remains approximately unchanged even if we increase the model parameter  h  (which corresponds to the film "viscosity") in 10⁵ times as demonstrated in Fig.stick17. 

The conclusion about a huge viscosity of a thin confined film is based on the assumption that the lubricant film is (ideally) homogeneous and, at the onset of slip, the whole film melts simultaneously so that the further slip proceeds with a liquid film. When the film instantly melts, the force  F  on the slider from the interface abruptly drops to the viscous force (which starts from zero being proportional to the slider velocity), while the force from the driving spring remains unchanged at the beginning. Therefore, if the slip is viewed as a uniform, massive lubricant melting event, then the limiting factor would be the lubricant's viscosity, and in that case very large viscosity values for the lubricant's viscosity are required, 10⁴ to 10⁷ times higher than that of the bulk lubricant.

However, this assumption has neither experimental nor theoretical support. From surface science physics it is known that a thin film may hardly be homogeneous on a meso- and even nanoscale, as well as the surfaces themselves are not ideally crystalline on these scales. Generally, a film consists of domains with different orientation or even different structure and, therefore, are characterized by different threshold stress yields. Thus, it is reasonable to assume that the film does not melt and begin to slide all at ones, but different domains start to slide one by one, as in the EQ model. In this case, the decrease of  F  is gradual owing to consequent breaking (melting) of different contacts (domains), and this approach describes the  X(t)  dependences observed experimentally.

The mechanism of slip onset (either sequential breaking of the contacts or the instant melting of the lubricant film) determines details of the slip kinetics which in principle may be resolved with SFB experiments. The order of magnitude of the slip time, however, in both cases is determined by the setup mechanical frequency  W_S.  Indeed, in our model adapted to the SBF setup, the sliding mass  M  is concentrated at the end of the spring, i.e., the slider plus spring form a usual pendulum. A pendulum is characterized by two characteristic times – the inverse of its frequency,  W^-1,  and the inverse of its damping coefficient,  h^-1; whichever of these times is shorter, it determines the system kinetics.

The frequency of ringing vibrations after the slip is either  W_S  (if  F(t)  oscillates around zero in the case of  t_d>W_s^-1)  or  W_L,  i.e., much higher in the case of  t_d<W_s^-1. However, the question is more involved: in the model we assumed that the bottom substrate is rigid and fixed. In a real setup, only the bottom of the base may be fixed, while the base has its own mass  M_B  and elasticity  K_B.  When the sliding stops and the two substrates are pinned together, the whole system of the mass  M_S=M_B+M+Nm  should oscillate with a frequency ~(K_B/M_S)^1/2. Unfortunately, usually we do not know the experimental values for these parameters.

Notice also that, in the case of  dF(v)/dv <0  which emerges because of contact's aging, the regimes of stick-slip and smooth sliding may be separated by a regime of irregular (chaotic) motion due to inertia effects.

Concerted slips

The ME approach is based on the earthquakelike model, which belongs to cellular automaton model and, thus, ignores real dynamics of contacts' motion. Namely, as the force on a given contact  i  reaches the corresponding threshold value, the contact instantly relaxes. In a real situation, however, the relaxation is not instant but takes some time; therefore, the delay time is always nonzero,  t_d>0.  Moreover, if the interaction between the contacts is taken into account, relaxation of neighboring contacts stimulated by the ith contact, also should propagate not instantly but continuously. This leads to appearance of collective (soliton-like) waves in the sliding interface.

...

...

Complete set of equations

Above we considered the bottom plane of the sliding block as a rigid plane, i.e., we ignored a possible variation (or fluctuation) of  F  in the  (x,y)  plane. To develop a full theory of friction on the mesoscale, one has to take into account the elastic deformation of the sliding blocks and the interface by introducing a position dependent distribution  Q[x, X(r)],  where  r={x,y}  denotes the position of an asperity on the interface, and  X(r,t)  describes the local translation of the interface averaged over the mesoscopic scale. The solution of the master equation (or a set of equations if the aging is included) then provides a local force (stress) at the interface  F(r,t).  The master equation has then to be coupled to an equation describing the elastic deformation of the block with the boundary condition at the interface, subjected to the contact forces  F(r,t)  at each point  r.  Besides, one may also take into account a possible heating of the contacts due to sliding. Such a program, of course, may be completed numerically only.

Conclusion

Thus, we described the master equation approach to the earthquakelike model with a continuous distribution of static thresholds, which is much more efficient than simulations and can be solved analytically in cases which are particularly relevant. It provides a deeper understanding of friction analyzed at the mesoscale in terms of the statistical properties of the contacts. This splits the study of friction in independent parts:

(i) the study of the contacts, which needs inputs from the microscopic scale,

(ii) the calculation of the friction force given by the master equation provided the statistical properties of the contacts are known, and

(iii) the incorporation of other aspects such as the interaction between the contacts and their aging.

This approach describes the stick-slip and smooth sliding regimes of tribological systems. The stick-slip motion emerges if and only if the following two conditions are satisfied: first, the slider is soft,  K<K^*, and second, the delay time in contact formation is nonzero,  t_d>0.  If at least one of these conditions is violated, the steady-state motion corresponds to smooth sliding. In stick-slip, the elastic energy stored in the slider during stick when  F(t)  grows, is then partially dissipated during the force dropping, and partially during the following (ringing) oscillations. The transition from stick-slip to smooth sliding with the increase of the driving velocity emerges due to aging of the contacts (although ringing vibrations may also contribute to the transition mechanism).

Note, however, that the ME approach is to be applied to meso- or macroscopic systems, but it cannot be used to explain the AFM/FFM tip-based experiments, except if there is a multi-contact through several tip atoms.

See the articles (pdf files may be found here):

O.M. Braun and M. Peyrard, Phys. Rev. Lett. 100 (2008) 125501 "Modeling friction on a mesoscale: Master equation for the earthquakelike model"; Phys. Rev. E 82 (2010) 036117 "Master equation approach to friction at the mesoscale"; and in preparation

O.M. Braun and E. Tosatti, Europhys. Lett. 88 (2009) 48003 "Kinetics of stick-slip friction in boundary lubrication"; and in preparation

O.M. Braun, I. Barel, and M. Urbakh, Phys. Rev. Lett. 103 (2009) 194301 "Dynamics of transition from static to kinetic friction"

Back to main page or tribology page

Last updated on October 26, 2009 by O.Braun.  Copyright © by O.Braun.  Translated from L^AT_EX by T_TH