## Macroscopic friction from the collective slip of

## contacting asperities

A.E. Boerma Master’s thesis

University of Groningen

Supervisor: prof.dr.ir. E. van der Giessen June 3, 2013

## Contents

1 Introduction 5

2 Dynamics of interface rupture 9

2.1 Equations of motion . . . 9

2.2 Analytic solutions . . . 12

2.3 Displacements in small systems . . . 16

3 Quasistatic simulation of rupture 21 3.1 Validation of quasistatic method . . . 23

3.2 Rupture front propagation . . . 24

3.3 Effect of friction forces on rupture speed . . . 26

3.4 Incorporating bulk elasticity . . . 27

4 Rupture under uniform loading 31 4.1 Almost uniform friction forces . . . 32

4.2 Distribution of friction forces . . . 33

4.3 Internal elasticity . . . 33

5 Contact of rough surfaces 37

6 Conclusions and discussion 41

A Supporting material 48

A.1 Validation of equation (21) . . . 48 A.2 Validation of equation (23) and equation (24) . . . 50 A.3 Source codes . . . 50

### Chapter 1

## Introduction

Constructing the equations of motion for a block sliding on an in-
clined plane subject to kinetic friction, or failing to slide subject to
static friction, is a standard problem in introductory mechanics text-
books.^{11,19}Friction is a reaction force that tends to stop two bodies
in contact from moving relative to one another and on the macro-
*scopic scale its magnitude F satisﬁes Amontons’ law*

*F≤ µP ,* (1)

*where P is the normal force and µ is the coefﬁcient of friction, an*
empirical parameter that depends on the materials in contact. In

Figure 1: A block on a ramp.

the classical theory of friction due to Amontons and Coulomb the coefﬁcient of friction does not depend on the relative velocity of the sliding bodies, other than that there is a distinction between kinetic friction (sliding) and static friction (no sliding), and, counter-intuitively, it does not depend on the area in contact either.

Modern theories of friction often do explain friction as a consequence
of contact area.^{4}The real contact area between two surfaces, how-

ever, is generally much smaller than it appears to be at a macro- scopic scale, because no surface is perfectly ﬂat at all length scales.

If we zoom in on what appears to be a smooth surface, at some length scale we will begin to see protrusions on the surface. We call these protrusions ‘asperities’. The contacting surfaces will only touch at the highest of these asperities, as illustrated in Figure 2.

Bowden and Tabor argued that the friction force will be proportional
*to the contact area A between the asperities so that*

*F* *≤ τA,* (2)

*and identiﬁed the proportionality constant τ as the shear stress nec-*
essary to break a junction between the two surfaces.^{4}

Figure 2: Rough surfaces

in contact. In this report, we investigate how macroscopic friction can arise from the collective behavior of asperities. As two surfaces start sliding relative to one another, the joints that form when asperities are put in contact break. At the same time, new connections are made.

This re-creation of junctions is necessary to balance the normal load.

Most interfaces have a lower coefﬁcient of friction for static friction than for dynamic friction, which implies that the average contact area is smaller when the surfaces are moving.

This process of making and breaking connections at the level of individual asperities manifests itself as on larger length scales as

‘self healing’ shear cracks,^{2,13,25}as illustrated in Figure 3. Junctions
break at the tip of such a crack, but new ones form some time after
the tip of the crack has passed.

Figure 3: Self-healing cracks move along the interface before surfaces start sliding.

Simulations suggest that such self-healing cracks occur at sub-micron
scales prior to the onset of frictional sliding.^{13}Experiments have also
shown that similar ‘rupture fronts’ traverse the interface on a macro-
scopic scale.^{2}These fronts initiate at the loaded edge and propagate
at velocities up to the shear wave speed. Behind a front the con-

tact area is reduced.^{27}Shear cracks can propagate with a speed at
most equal to the Rayleigh wave speed, but the greatest reduction
in contact area occurs after passage of “anomalously slow” fronts:

those of speeds an order of magnitude lower than the Rayleigh wave
speed.^{28}

The transition to sliding is accompanied by such a slow detachment
front. Before the surfaces properly start sliding, a number of ‘pre-
cursor events’ will have taken place. As the loading on a sample is
increased, rupture fronts are created that partially traverse the inter-
face. These fronts move at Rayleigh wave speed and arrest abruptly
before reaching the leading edge of the sample. The rupture di-
rection is not necessarily the same as the subsequent direction of
slip.^{3,29}

Signiﬁcant variations in both the normal and shear stresses are re- ported along nominally smooth and statistically uniform interfaces.

The propagation speed of a rupture front is found to be determined
by the local ratio of normal to shear stress. As a rupture front passes,
the stresses along its path are redistributed. This paves the way for
subsequent detachment fronts. When frictional sliding ﬁnally occurs,
the contact area distribution has already become highly nonuniform.^{3}
It is speculated that this redistribution of stresses induces ‘local’

Amontons-Coulomb friction.^{29}When further macroscopic variations
in the normal load are applied, Maegawa et al. found that the length
along which precursor fronts propagate increases as the normal load
on the edge at which they initiate is decreased.^{18}

1 2 3 *N*

*K* *k*
*F*

*v* ^{1}

*P*
*m*

Figure 4: Mechanical model of frictional sliding.

Various numerical studies have been able to reproduce these pre-
cursor events using a model similar to that illustrated in Figure 4. The
numerical simulations show a large number of ‘small’ rupture events
not visible in experimental results.^{5,18,31}The studies cited, however,
all apply Amontons-Coulomb friction to the blocks individually: no-
tably, they distinguish between kinetic and static coefﬁcients of fric-

tion, whereas the transition from static to kinetic friction is precisely what we are interested in.

As a simple model of non-rigid coupling between contacts, we con- sider a string of blocks connected by springs. Similar models in- clude the Frenkel-Kontorova model, which consists of a string of Figure 5: Schematic of

the Frenkel-Kontorova model.

harmonic oscillators that is pulled over a periodic potential, the Tom- linson model, which consists of independent oscillators coupled to

Figure 6: Schematic of the Tomlinson model.

a rigid, moving bulk that is pulled over in a periodic potential, and the Burridge-Knopoff (BK) model, which consists of coupled oscil-

Figure 7: Schematic of the Burridge-Knopoff model.

lators that each show Amontons-Coulomb friction. These models
are frequently employed to model kinetic friction^{7,8,20,21}, and partic-
ularly the combined Frenkel-Kontorova and Tomlinson models (FKT
model) are used to model atomic scale friction processes^{24,32}.
In this report we will study variants of the Burridge-Knopoff model.

One characteristic of this model is a velocity-dependent friction law.

A number of these friction laws are illustrated in Figure 8.^{7,9}In this
report we will use constant, velocity-independent friction. Moreover,
we will not distinguish between static and kinetic friction.

*v*

*F*

Constant Carlson Cartwright

Figure 8: Examples
of velocity-dependent
friction laws.^{7,9} In this
report, we will use
*velocity-independent,*
constant friction.

### Chapter 2

## Dynamics of interface

## rupture

^{1}

^{2}

^{3}

^{N}*K* *k*
*F*

*v* ^{1}

*P*
*m*

Figure 9: Mechanical model of frictional sliding.

*We consider a model as illustrated in Figure 9, of a series of N blocks*
*of mass m. The blocks are internally coupled by springs of stiffness*
*k*between nearest neighbors, and the ﬁrst block is loaded through a
*spring with spring constant K that is pulled with a constant speed v.*

Based on dimensional analysis, we expect the motion of the system
*to be governed by two dimensionless parameters p**k* *:= K/k* and
*p**v**:= v√*

*mk/F*. We will, like the studies cited in the previous chap-
ter^{5,18,31}, simulate the initiation of stick-slip motion of such a system,
but we will not differentiate between static and kinetic friction.

### 2.1 Equations of motion

*Two distinct sets of forces act on the blocks: spring forces F*_{n}^{S}and
*friction forces F*_{n}^{F}*. The equation of motion for the n-th block is there-*
fore

*m¨u**n**= F*_{n}^{S}*+ F*_{n}^{F} (3)

*where ¨u**n* *:= d*^{2}*u**n**/dt*^{2} is the block’s acceleration and where the
spring forces are given by Hooke’s law^{⋆}

*⋆*For the ﬁrst and ﬁnal blocks we
*have F*_{1}^{S} = *−(k + K)u**n*+
*Kvt + ku**n+1* *and F*_{N}^{S} =

*k(u*_{n−1}*− u**n*)instead. *F*_{n}^{S}*= k(u**n**−1**− u**n**) + k(u**n+1**− u**n*). (4)

The ﬁnal term in equation (3),

*F*_{n}^{F}*=: F η(F*_{n}^{S}*/F, ˙u**n*), (5)
*represents friction. In this deﬁnition F is the maximum friction force*
*and η is a function that selects the magnitude and direction of the*
friction force, based on the spring forces acting on the block and its
*current velocity ˙u:*^{⋆}

*⋆*In practice, we test whether

*|a| < ϵ for some small ϵ instead*
*of a = 0.*

*η(a, b) :=*

*−a* *if b = 0 and |a| ≤ 1,*

*− sgn a* *if b = 0 and |a| > 1,*

*− sgn b* *if b ̸= 0,*

(6)

where sgn is the signum function

*sgn(a) :=*

*−1* *if a < 1,*
0 *if a = 0,*
+1 *if a > 1.*

(7)

The minus signs in equation (6) serve as a reminder that the fric- tion force is always opposite in direction to the motion of a block.

For a block that is initially at rest, the friction force exactly cancels
*the spring forces until they exceed a maximum value F . When the*
block is sliding, the friction force acts in the opposite direction to the
motion. If a block comes to rest, the friction force again tends to
cancel the spring force so as to keep the block from moving.

We can combine the spring forces for all blocks into a matrix

*F*_{1}^{S}

...
*F*_{N}^{S}

*= k*

*−1 − p**k* 1

1 *−2* 1

...

1 *−2* 1

1 *−1*

*u*1

...
*u**N*

, (8)

*where p**k* *:= K/k*is one of the dimensionless parameters we iden-
tiﬁed earlier. If we let boldface symbols represent matrices, we can
write this equation more concisely as

**F**^{S}* = k κ· u.* (9)

With the further deﬁnitions

**η :=**

*η*1

*η*2

...
*η**N*

:=

*η(F*_{1}^{S}*/F, ˙u*1)
*η(F*_{2}^{S}*/F, ˙u*2)

...
*η(F*_{N}^{S}*/F, ˙u**N*)

**and e :=**

1 0 ... 0

(10)

we combine the equations of motion for all the blocks,

*m¨ u = kκ· u + Kvte + F η.* (11)

**If we deﬁne u =: u***c*

**x**

*and t =: t*

*c*

*s, where x and s are nondimensional*

*and u*

*c*

*and t*

*c*are some characteristic length scale and time scale for the problem, then equation (11) becomes

*mu**c*

*F t*^{2}_{c}*d*^{2}

*ds*^{2}**x =***ku**c*

*F* **κ****· x +***Kvt**c*

*F* * se + η.* (12)

*This result indicates that it might be fruitful to set u**c* *:= F/k* and
*t** _{c}* :=√

*m/k, so that the fully nondimensional form of equation (11)*

can be rewritten as
*d*^{2}

*ds*^{2}**x = κ****· x + p***v**p**k** se + η,* (13)

*where p*

*v*

*:= v√*

*mk/F* is the other dimensionless parameter we
*identiﬁed before. Although only the product p**v**p**k* of the nondimen-
*sional parameters appears explicitly in this equation, p**k* appears by
*itself as a permutation of κ*11.

In the above we have chosen our characteristic time and length scale
*with respect to the ‘internal’ parameters m, F and k. Alternatively,*
*we could have set t = t*^{′}*c**ζwith t*^{′}*c**=: u**c**/v, so that the time scale of*
the problem is expressed with respect to an external variable. Then
equation (13) changes to

*p*^{2}_{v}*d*^{2}

*dζ*^{2}**x = κ****· x + p***k** ζe + η.* (14)

We can go from a solution of equation (13) to a solution of equa-
*tion (14) by substituting p**v**sfor ζ. The formulation in equation (14)*
is useful because it allows us to easily take the quasistatic limit as
*p*_{v}*goes to zero (i.e., as v goes to zero). It is illustrative to see how*
*the different interpretation of the time-like variables s and ζ arises if*
*we let p**v* *vanish from equation (13) and equation (14): as p**v* goes
to zero in equation (13) we get dynamic effects with no change in
*external variables, corresponding to the interpretation of s as an in-*
*ternal time scale, whereas as p**v* goes to zero in equation (14) the
internal dynamics vanish.

### 2.2 Analytic solutions

Although we cannot give a closed-form solution to equation (13),
we can give an analytic solution to the equations of motion when a
*small number of blocks is slipping. During any time interval where η**n*

is constant for the blocks that are sliding, equation (13) reduces to a
*system of linear ordinary differential equations similar to ω*^{−2}*x*^{′′}*+x =*

*cs− 1.** ^{⋆}*The solution to this ODE is

*From here on, primes denote*

^{⋆}*differentiation with respect to s.*

*x(s) = a cos(ωs) + b sin(ωs) + cs− 1,* (15)
which indicates that the displacement of the blocks as a function of
time will have both harmonic and linear components. The magnitude
*aand b of the oscillations can be determined from initial conditions,*
i.e. from the blocks’ initial velocity and position.

We will ﬁrst consider the case when only the ﬁrst block is moving, so
*that x**n**(s) = 0for n > 1. The ﬁrst block slips when x*^{′′}_{1} *= p*_{k}*p*_{v}*s−1 ≥*
0*, which implies s ≥ 1/p*^{k}*p**v*. Its equation of motion is

*x*^{′′}_{1} =*−(1 + p**k**)x*1*+ p**v**p**k**s− 1,* (16)
*with x*1= 0*initially, in which we substitute ˆs = s − 1/p*^{v}*p**k*to get

*x*^{′′}_{1} =*−(1 + p**k**)x*1*+ p**v**p**k*ˆ*s.* (17)
The general solution of this equation is

*x*1(ˆ*s) = a cos(*√

*1 + p**k**s) + b sin(*ˆ √

*1 + p**k*ˆ*s) +* *p**v**p**k**s*ˆ
*1 + p**k*

. (18)

*The initial condition x*1(ˆ*s = 0) = 0gives a = 0 and x*^{′}_{1}(ˆ*s = 0) =*

*√1 + p**k**b + p**v**p**k**/(1 + p**k*) = 0*gives b = −p*^{v}*p**k**/(1 + p**k*)* ^{3/2}*, which
implies

*x*1*(s) =* *p**v**p**k**s− 1*

*1 + p**k* *−* *p**v**p**k*

*(1 + p** _{k}*)

*sin(√*

^{3/2}*1 + p**k**(s− 1/p**v**p**k*)), (19)
*under the assumption that the friction force on that block is η*1=*−1.*

While the harmonic terms will lead to an oscillation, due to the linear term the center of oscillation moves in the positive direction, so that the velocity of the block should never be negative.

For the motion of the ﬁrst block, the ratio of magnitude of the linear
term to that of the harmonic term is*√*

*1 + p**k*, so the two are approx-
*imately of the same magnitude if p**k* *≪ 1 i.e. k ≫ K, but the linear*
*term dominates if k ≪ K. As the effective spring constant of the*
*entire slider consisting of N springs in series is k/N, in the contin-*
*uum limit we must have k ∝ N, so in the continuum limit the linear*
and harmonic effects are of the same magnitude for the ﬁrst block.

*The second block slips when x*^{′′}_{2} *= x*1*(s)−1 ≥ 0, so when x*1*(s)≥ 1.*

*Using the substitution ˆs =√*

*1 + p*_{k}*(s− 1/p**v**p** _{k}*)we get
1 =

*p*

_{v}*p*

_{k}*(1 + p**k*)* ^{3/2}*(sin(ˆ

*s) + ˆs).*(20) There is no closed-form solution to this equation, but numerical meth- ods can ﬁnd a solution quickly enough.

*Next we consider the general case when n blocks are sliding. We*
can obtain the corresponding equation of motion by taking the ﬁrst
*n* *out of the N equations in the system (13). Again, we deﬁne s*0

to be the time when the last slip event occurred and change the
*time-like variable ˆs := s − s*^{0}. We attempt a solution of the form of
equation (15)

**x(ˆ***s) = − ˆκ*

^{−1}*n*

*0*

**· (p(ˆs + s****)e**

*+*

**− 1)**∑*n*
*k=1*

*[c*^{+}_{k}*exp(+iω**k**s) + c*ˆ _{k}* ^{−}*exp(

*−iω*

*k*ˆ

**s)]ξ***, (21)*

_{k}*where p = p**v**p*_{k}**and ξ**_{1} **through ξ**_{n}*are eigenvectors of the n × n*
matrix

**ˆ**
**κ*** _{n}*=

*−1 − p**k* 1

1 *−2* 1

...

1 *−2*

, (22)

*with the associated eigenvalues λ*1 *through λ**n*. Inserting this into

equation (13) (see the Appendix), we ﬁnd that the eigenfrequen-
*cies ω**k* *of the sliding system are related to the eigenvalues λ**k* of
**ˆ**

**κ**through ω^{2}* _{k}*=

*−λ*

*k*.

**If x**0**= x(ˆ***s = 0)***is the initial position and x**^{′}_{0} **= x*** ^{′}*(ˆ

*s = 0)*the initial

**velocity of the blocks, and Ξ is a matrix whose m-th column is ξ***m*,

*then the coefﬁcients c*

^{±}*can be determined as (see the Appendix)*

_{m}**c**^{+}**+ c**^{−}**= Ξ**^{−1}* · [x(s*0

**) + ˆ**

**κ**

^{−1}*· (ps*0

**e**

*(23)*

**− 1)]****c**

^{+}

**− c**

^{−}**= i Ξ**

^{−1}

**· ω**

^{−1}

**· [x**

^{′}*(s*

_{0}

**) + ˆ**

**κ**

^{−1}*(24) where*

**· e],****c**^{+}=

*c*^{+}_{1}

...
*c*^{+}_{n}

** and c*** ^{−}* =

*c*^{−}_{1}

...
*c*^{−}_{n}

. (25)

The solution (21) that we have now obtained is valid from the moment
*when the n-th block slips to the moment when the n + 1-st block*
*slips, which happens when x**n*= 1, or more explicitly, when

*1 =[0 . . . 0 1] · x(ˆs)*

=**− [0 . . . 0 1] · ˆκ**^{−1}* · (p(ˆs + s*0

**)e**

**− 1)***+ [0 . . . 0 1]·*

∑*n*
*k=1*

*[c*^{+}_{k}*exp(+iω**k**s) + c*ˆ ^{−}* _{k}* exp(

*−iω*

*k*

*ˆ*

**s)]ξ***. (26)*

_{k}Like equation (20), we must use numerical methods to ﬁnd the small- est root for this equation.

Putting all of the above together, we can solve the equations of mo- tion for the entire system by

*1. determining the values of s at which a block slips, i.e.,*

*(a) for the ﬁrst block s = 1/p**k**p**v*,

(b) for the second block from equation (20),
*(c) for the n + 1-st block from equation (26);*

2. determining coefﬁcients from equation (23) and equation (24);

and

3. putting the appropriate numbers into equation (21), of which equation (19) is a special case,

**under the assumption that the friction term η in equation (13) is con-***stant (η**n* =*−1) for blocks that have slipped.*

### 2.3 Displacements in small systems

We have solved equation (13) using velocity-Verlet integration for a
*system of N = 5 blocks until the time when the ﬁfth contact started*
*sliding. The time step was ∆s = 10*^{−4}*/√*

*1 + p**k*,* ^{⋆}* and the initial

*⋆*The reason for this is that
the natural period of the ﬁrst
block when it oscillates freely is
*1/**√*

*1 + p**k*.

**conditions were x = x**^{′}**= 0** *and s = 1/p**v**p**k**− 10∆s, so that the*
ﬁrst block started sliding after ten time steps. The calculation ended
*when x*5*> 10** ^{−4}*.

*The displacements, for p**v* =10^{–1}, 1, 10^{+1}*and p**k* =10^{–1}, 1, 10^{+1},
are plotted in Figure 10 (solid lines).* ^{⋆}*Results for the analytic dynamic

*⋆*In the Appendix (Figures 42
and 41) we provide plots of the
*displacements for p**v* *and p**k*

from 10^{–3}through 10^{+3}.

scheme derived in the previous section are plotted as well (open symbols).

*In the fully numerical results when p**v* *< 1or p**v**p**k* *< 1, the velocity*
of some blocks (plotted in Figure 11) is not strictly positive after they
have started sliding. At various times the velocity of some blocks
becomes zero again, however, these blocks do not ‘stick’ in the
*usual sense: whenever a block’s velocity is x*^{i}_{n}* ^{′}* = 0

*at time step i,*

*it is x*

^{i+1}

_{n}

^{′}*> 0*in the next time step. This ‘stick-slip’ cycle will then repeat for some time until the block starts purely sliding again.

0 50 100 0

2 4 6 8

*s*

*x*

0 5 10

0 2 4 6 8 10

*s*

*x*

0 1 2 3

0 10 20 30

*s*

*x*

0 50 100 150

0 2 4 6 8 10

*s*

*x*

0 5 10

0 2 4 6 8

*s*

*x*

0 1 2 3 4

0 10 20 30

*s*

*x*

0 200 400 600

0 2 4 6 8 10

*s*

*x*

0 20 40 60

0 2 4 6 8 10

*s*

*x*

0 2 4 6 8

0 5 10

*s*

*x*

Figure 10: Displacement
as a function of time for
*p**v*=10^{–1}, 1, 10^{+1}(left to
*right), and p**k* =10^{–1}, 1,
10^{+1}(bottom to top).

0 50 100

0 0.05 0.1 0.15

*s*

*x’*

Numerical Analytical

0 5 10

0 0.5 1 1.5

*s*

*x’*

Numerical Analytical

0 1 2 3

0 5 10 15

*s*

*x’*

Numerical Analytical

0 50 100 150

0 0.05 0.1

*s*

*x’*

Numerical Analytical

0 5 10

0 0.2 0.4 0.6 0.8 1

*s*

*x’*

Numerical Analytical

0 1 2 3 4

0 2 4 6 8 10

*s*

*x’*

Numerical Analytical

0 200 400 600

0 0.01 0.02 0.03 0.04

*s*

*x’*

Numerical Analytical

0 20 40 60

0 0.1 0.2 0.3 0.4

*s*

*x’*

Numerical Analytical

0 2 4 6 8

0 1 2 3

*s*

*x’*

Numerical Analytical

Figure 11: Velocity as a
*function of time for p**v*=
10^{–1}, 1, 10^{+1}(left to right),
*and p**k* = 10^{–1}, 1, 10^{+1}
(bottom to top).

Whether or not this effect occurs does not depend on the time step size, so we do not expect it to be a numerical issue. Most likely, it actually represents very slow sliding. The reason that we do not see true stick-slip is that in our model there is no real transition between the stuck and the sliding state: the maximum magnitude of friction is

*|η**n**| ≤ 1 for both kinetic and static friction. Finally, we never observe*
negative velocities in these results.

Using the analytic scheme we devised in the previous section, we get
results that seem to agree quite well with those from the numerical
integration. Even though the quantitative difference between these
sets of results is very small, there is a qualitative difference in that the
velocities calculated using the analytic scheme do become negative
at some times. This occurs in systems where the numerical results
*show the aforementioned ‘stick-slip’ motion, so when p**v* *< 1* or
*p**v**p**k* *< 1. Negative velocities are at odds with the assumption we*
*made in the previous section, that η**n* =*−1 at all times for a block*
that has slipped, because for negative velocities the friction term
*η** _{n}*= +1is positive.

log_{10}*p** _{v}*
log10

*p*

*k*

Stick−slip

Pure sliding

−5 0 5

−5 0 5

Figure 12: Phase dia- gram showing the tran- sition between ‘stick-slip’

and pure sliding in a sys-
*tem of N = 5 blocks (i.e.,*
four sliding blocks).

In spite of the fact that one of our assumptions in using equation (21) appears to be invalid, the quantitative difference between the two sets of results is negligible. The reason for this is probably that the velocities and the oscillations in them are too small for the dynamic effects to matter much. This notion will be further investigated in the next Chapter.

This threshold for stick-slip versus pure sliding, which appears around
*min(p*_{v}*, p*_{v}*p** _{k}*)

*≈ 1 for N = 5, is dependent on the system size N:*

*when N is increased the solid line in Figure 15 shifts to min(p**v**/p*^{∗}_{v}*,*
*p*_{v}*p** _{k}*)

*≈ 1, where p*

^{∗}*v*

*≈ N/3. We can contrast this with the expected*scaling of the system parameters. In particular, say we model a

*slider with total mass m*

^{(tot)}

*, linear stiffness k*

^{(tot)}and friction force

*F*

^{(tot)}

*as a series of N + 1 blocks. Each of the blocks has mass*

*m = m*^{(tot)}*/(N + 1)and friction force F = F*^{(tot)}*/(N + 1), and each*
*of the internal springs has spring constant k = Nk*^{(tot)}, so the nondi-
mensional groups scale as

*p** _{k}*=

*K*

*k*

*∝*1

*N* *and p**v* =*v√*
*mk*

*F* *∝ N.* (27)

*This implies that p**v**p*_{k}*is independent of N. Furthermore, both p**v**∝*
*N* *and p*^{∗}_{v}*∝ N, so the criterion for stick-slip min(p**v**/p*^{∗}_{v}*, p**v**p**k**) < 1*
seems to be independent of the system size. In fact, the above
analysis suggests that whether stick-slip occurs in a system of any
size is determined by the macroscopic parameter

*p**v**p**k*= *vK*
*F*^{(tot)}

√
*m*^{(tot)}

*k*^{(tot)}. (28)

*This parameter can be interpreted as a ratio of the rate vK of ap-*
*plied force, versus the resisting force F*^{(tot)} divided by a time scale

√*m*^{(tot)}*/k*^{(tot)}. ^{−1}^{−5} ^{0} ^{5}

0 1 2 3 4

log10ω

log_{10}*p*_{k}

Figure 13: Oscillation
frequencies in a system
*of N = 4 blocks. Solid,*
dashed, dash-dotted
and dotted lines refer to
one through four sliding
blocks.

In Figure 13 we have plotted the eigenfrequencies that occur in the
analytic dynamic model. We can distinguish two sets of frequencies
*that scale differently as a function of p**k**⋆*: the frequencies of the ﬁrst

*⋆*Recall that the eigenfrequen-
cies are related to the eigenval-
**ues of κ, which does not de-***pend on p**v*.

block, which are in the order of*√*

*1 + pk, and the frequencies of*
blocks 2 through 4, which are all in the order of 1.

### Chapter 3

## Quasistatic simulation of rupture

In the one-dimensional nearest-neighbor model introduced in the
*previous chapter, the velocities of blocks are in the order of min(p**v**, p**v**p**k*),
*whereas the greatest characteristic frequencies are around max(1, p**k*).

Both these results tell us something about the dynamics of the model.

In this chapter we will look into which characteristics of our model remain if we ignore those dynamics.

To investigate in what extent the motion of the blocks is governed by dynamic effects, we contrast the results of the previous chapter with

*a quasistatic formulation. We set p**v*= 0in equation (14) to obtain^{⋆}^{⋆}**Recall that x = u k/F ,***where k is an internal spring*
*constant, and ζ = t vk/F is*
a time-like variable.

**0 = κ****· x + p***k** ζe + η.* (29)

*This is a piece-wise linear system consisting of N equations with*

*2N + 1degrees of freedom: all η*

*n*

*, all n*

*n*

*and ζ. The piece-wise be-*havior is caused by the fact that the absolute value of each element

**of η by deﬁnition satisﬁes |η**^{n}*| ≤ 1.*

* We can solve for η at the time when n blocks are sliding by writing*
equation (29) for those blocks as

0 0 ... 0

=

*−1 − p**k* 1

1 *−2* 1

...

1 *−2*

*x*1

*x*2

...
*x**n*

+

*p**k**ζ*

0 ... 0

*−*

1 1 ... 1

. (30)

This is equivalent to

1 1 ... 1

=

*−1 − p**k* 1 1

1 *−2* 1 0

... ...

1 0

*x*1

...
*x**n**−1*

*p**k**ζ*

+

0

... 1

*−2*

*x** _{n}*. (31)

*The n + 1-st block slips when x**n*= 1, so we can determine the time
*ζ*0at which this occurs by calculating

*x*1

...
*x**n**−1*

*p**k**ζ*0

=

*−1 − p**k* 1 1

1 *−2* 1 0

... ...

1 0

*−1*

1 1 ... 1

*−*

0 ... 1

*−2*

. (32)

*The ﬁrst n − 1 elements of this array are the positions x**m**(ζ*_{0})of the
*ﬁrst n − 1 blocks at ζ = ζ*0*, so by setting the ﬁnal element to x**n*= 1
* we get x(ζ*0). For example, the ﬁrst block slips at ζ = 1/p

*k*, the

*second block at ζ = 1 + 2/p*

*k*and the third at

[
*x*1

*p*_{k}*ζ*_{0}
]

=

[*−1 − p**k* 1

1 0

]* _{−1}*[[

1 1 ]

*−*
[

1

*−2*
]]

= 3 [

1
*1 + p*_{k}

] . (33)

*The variable ζ = p**v**s = tvk/F* was introduced as a dimensionless

0 50 100 0

2 4 6 8 10

*s*

*x*

0 5 10

0 2 4 6 8 10

*s*

*x*

0 1 2 3

0 10 20 30

*s*

*x*

0 50 100 150

0 2 4 6 8 10

*s*

*x*

0 5 10 15

0 2 4 6 8 10

*s*

*x*

0 1 2 3 4

0 10 20 30

*s*

*x*

0 200 400 600

0 2 4 6 8 10

*s*

*x*

0 20 40 60

0 2 4 6 8 10

*s*

*x*

0 2 4 6

0 5 10

*s*

*x*

Figure 14: Displacement
as a function of time for
*p**v*=10^{–1}, 1, 10^{+1}(left to
*right), and p**k* =10^{–1}, 1,
10^{+1}(bottom to top).

‘time’ in the previous chapter. In the quasistatic limit it is more prop- erly the position of the loading section, but for ease of discussion we will still refer to it as ‘time’.

### 3.1 Validation of quasistatic method

In Figure 14 we have plotted the displacements so obtained for sys-
*tems of N = 5 blocks (i.e., four sliding blocks) with p**v**, p**k* =10^{–1}, 1,
10^{+1}, along with the displacements for the dynamic model obtained

using the analytic method.* ^{⋆}* In Figures 42 and 41 in the appendix

^{⋆}*We have substituted p*

_{v}*s*

*for ζ.*

*displacements are plotted for greater ranges of p**v* *and p**k*, for this
quasistatic method as well as the numerical and analytic dynamic
methods described in the previous Chapter. Figure 14 suggests
that the quasistatic approximation agrees with the dynamic results
*for p**v**≤ 1.*

log_{10}*p** _{v}*
log10

*p*

*k*

Quasi-static

Dynamic

10%

1%

−5 0 5

−5 0 5

Figure 15: Phase dia- gram showing the tran- sition between ‘stick-slip’

and pure sliding in a sys-
*tem of N = 5 blocks (i.e.,*
23

In Figure 15 (similar to Figure 12 from the previous chapter, which
described the transition between stick-slip and pure sliding), we have
plotted two dashed lines that indicate the approximate error in us-
*ing the quasistatic results. This error is deﬁned as |∆x*1*/x*_{1}*|, where*
*x*1*is the displacement of the ﬁrst block at the time when the N-th*
*block slips, as calculated using the dynamic model, and where ∆x*1

*is the difference between x*1 calculated using the dynamic model
and the corresponding value calculated using the quasistatic model.

*We have used a system of N = 5 blocks here, but the results do*
not appear to depend on the system size (see also Figure 44 in the
*Appendix). In general, the error is small when p**v**≪ 1 or p**v**p*_{k}*≪ 1.*

From the previous chapter we recall that if we keep the character-
*istics m*^{(tot.)}*, F*^{(tot.)} *and k*^{(tot.)} of the entire spring constant, then the
*internal parameters scale as m = m*^{(tot.)}*/N, F = F*^{(tot.)}*/N, and*
*k = (N* *− 1)k*^{(tot.)}*. As p**v* *∝ N − 1 it seems safe to assume that*
*the condition p**v**≪ 1 is irrelevant in the limit of large N. This implies*
that the validity of the quasistatic approximation is also determined
by the macroscopic parameter

*p**v**p**k*= *vK*
*F*^{(tot.)}

√*m*^{(tot.)}

*k*^{(tot.)}, (34)

which was also the deﬁning parameter in stick-slip versus pure slid- ing. In the Appendix, Figure 44, we

### 3.2 Rupture front propagation

The computational simplicity of the quasistatic model allows us to
analyze the evolution of the number of sliding blocks as a function of
time. Because in this one-dimensional system blocks start sliding in
sequence, the number of sliding blocks is a measure of how far the
rupture front has propagated along the interface. We suggest that
this can be compared to the length of global slip precursor events
reported by Rubinstein et al., among others.^{5,29,31}

0 0.5 1

0 0.2 0.4 0.6 0.8 1

*s (normalized)*

# sliding blocks (norm.)

*N=2*
*N=4*
*N=8*
etc.

Figure 16: Rupture length for systems with

2, 4, 8, …, 1024 24 June 3, 2013

Figure 16 shows the number of sliding blocks as a function of time
*for systems of N = 2, 4, 8, …, 1024 blocks, with p**k* =1. Both
axes are normalized for easy comparison; it bears to mention that
*0 on the horizontal axis here corresponds to s = 1/p**k**p** _{v}* (i.e., the

*time at which the ﬁrst block slips) instead of s = 0 as in the rest of*this report. As the number of blocks is increased, the curves tend

*towards √s. This trend is particularly visible if we re-plot Figure 16*

with logarithmic axes (Figure 17). ^{10}^{−3} ^{10}^{−5} ^{10}^{0}

10^{−2}
10^{−1}
10^{0}

*s (normalized)*

# sliding blocks (norm.) *N=1024*

*N=512*
*N=256*

*N=128*
*N=64*

*N=32*
*N=16**N=8**N=4**N=2*

Figure 17: As Figure 16, but with logarithmic axes. Slopes tend towards 1/2.

*Using N = 100 and different values of p**k*, we ﬁnd that there are two
*limiting cases for the rupture front propagation length D as a function*
*of time. As N increases, D(s) tends toward √s, which is the trend*
*we discussed above. The same scaling occurs if p**k*increases, which
*is illustrated in Figure 18: the solid lines (for p**k**≥ 1) converge to√*

*s.*

*On the other hand, if we let p**k**decrease, then D(s) approaches s.*

0 0.5 1

0 20 40 60 80 100

*s (normalized)*

# sliding blocks

Figure 18: Rupture
length for a system with
*N* = 100 blocks, with
*p**k*=10^{–6}, 10^{–5}, …, 10^{+6}
*(dotted lines p**k**<*1).

*This dependence on p**k* occurs in tandem with the dependence on
*N. In the previous Chapter we found that if the macroscopic proper-*
*ties of the slider are kept constant, then p**k* *∝ 1/N. This, combined*
with the ﬁndings in Figure 18, suggests that in the continuum limit

lim

*N**→∞**D(s)*= lim^{?}

*p**k**→0**D(s)= s,*^{?} (35)
whereas the results in Figure 16 suggest

lim

*N**→∞**D(s)*=^{?} *√*

*s.* (36)

−3 −2 −1 0 1 2 3

0.5 0.6 0.7 0.8 0.9 1

log_{10}*p*_{k}

log* D* / log *s*

*N=10*
*N=100*
*N=1000*

Figure 19: *As p**k* in-
creases, the rupture front
scaling transitions from
*D∝ s to D ∝√s.*

*If we look at the scaling (log D/ log s) as a function of the macro-*
*scopic parameter p**k* *= K/k*^{(tot.)}*= N K/k, as illustrated in Figure 19,*
*we can see that the length of the rupture front D scales linearly when*
*p*_{k}*≪ 1, i.e. if the slider is very stiff compared to the loading section,*
*that D ∝√*

*swhen p**k* *≫ 1, and that there is a transition between*
*log D/ log s =1 and log D/ log s = 1/2 when 10*^{–2}*< p**k**<*10^{+3}.

The instantaneous change of the rupture length in response to a
*change in loading conditions is given by D*^{′}*= dD/ds. We recall*
*that our time-like variable s does not actually represent time here;*

rather, it represents the position of the loading section. We will refer
*to dD/ds as the ‘speed’ of the rupture front, for lack of a better term,*
but note that this is not a physical speed per se. In particular, it has
nothing to do with the speed at which cracks propagate along the
interface.

The analysis above implies that the rupture front speed is constant
*when p**k**≪ 1, and that it is otherwise strictly decreasing.*

### 3.3 Effect of friction forces on rupture speed

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

D

s (norm.) η=1,...,N η=N,...,1 η=1,...,1

Figure 20: Rupture front
*size D with different η**n*.

10^{-2}
10^{-1}
10^{0}

10^{-4} 10^{-3} 10^{-2} 10^{-1} 10^{0}

D

s (norm.)

Figure 21: As Figure 20, with logarithmic axes.

Up to this point we have assumed that each block slips at precisely

*|F**n*^{F}*| = F , the same value for each block. In order to eventually*
apply our model to real interfaces, we will now allow for variations
in this critical friction force along the interface. This means that the
* elements of η in equation (29) are no longer restricted to values be-*
tween –1 and +1, and that the equations following equation (29) are
changed accordingly.

We can expect that the rupture front speed will change if the distri-
*bution of η**n* changes. For example, blocks that have a high critical
friction force have a high resistance to sliding, so the rupture front
will arrest at that block.

*We have plotted the rupture length as a function of s (normalized for*
easy comparison) in Figures 20 and 21. In the case where the crit-
ical friction forces are linearly decreasing from the loaded edge, the
*length of the rupture front initially scales as D ∝* *√s, similar to the*
case where the critical friction forces are uniform, but as the normal-
*ized s approaches 1, in the ﬁnal stages of the simulation, the rupture*
front speed increases. In the case where the critical friction forces are

0 0.2 0.4 0.6 0.8 1

0 20 40 60 80 100

dD/ds, η

x (dist. along slider)

0 0.2 0.4 0.6 0.8 1

0 20 40 60 80 100

dD/ds, η

x (dist. along slider)

0 0.2 0.4 0.6 0.8 1

0 20 40 60 80 100

dD/ds, η

x (dist. along slider)

Figure 22: Rupture front
*speed dD/ds (blue) and*
*critical friction force η**n*

(red) as a function of dis- tance along the sliding in- terface. The slider was loaded from the left.

*linearly increasing from the loaded edge, we ﬁnd log D/ log s < 1/2*
instead. This implies that the rupture front length increases quickly
at ﬁrst, but also that its growth decelerates much more rapidly than
in the other two cases.

*In Figure 22 the ‘speed’ dD/ds of the rupture fronts is plotted along*
with the distribution of critical friction forces. This Figure seems to
indicate that there is no one-to-one link between the critical friction
force at one block and the local speed of a rupture front.

.∂

### 3.4 Incorporating bulk elasticity

The model we have considered places a constraint on the way a rupture front can propagate along the interface: because blocks slip sequentially, information can only travel along the interface and the rupture front can become ‘pinned’ at a block. If we apply loading to the left edge of the slider, blocks will slip one by one, from left to right. There is no reason to assume that real interfaces will rup- ture in such a sequential fashion. In reality, stresses that develop along the interface can redistribute through the bulk of both slider and substrate. This means that the rupture front can be felt from a distance, whereas if we have only nearest-neighbor interaction the rupture front cannot be felt over any distance greater than that be- tween neighboring blocks.

Figure 23: Rupture front
size for nearest-neighbor
*model with η = 1, . . . , N*
(solid); nearest-neighbor
*model with η = N, ..., 1*
(dotted); long-range
*model with η = 1, . . . , N*
(dashed) and long-range
*model with η = N, . . . , 1*
(dash-dotted).

### 0 0.2 0.4 0.6 0.8 1

### 0 0.2 0.4 0.6 0.8 1

### D

### s (norm.)

In order to incorporate this long-range interaction we change the
* matrix κ in equation (29) and the equations that follow in such a way*
that there is a spring of stiffness 1 between neighboring blocks, a
spring of stiffness 1/2 between next-nearest neighbors (for example,
blocks 1 and 2), stiffness 1/3 between third-nearest neighbors (for
example, blocks 1 and 4), and so on. If a block is displaced, it exerts
a spring force on all other blocks, and those forces are inversely
proportional to the distance between blocks. This scaling of the
spring forces is similar to the scaling of stresses due to an applied
displacement in plane strain

^{16}.

*In Figure 23 we have plotted the rupture front size D for η**n**= 1, . . . , N*
*(linearly increasing from the loaded edge) and η**n* *= N, . . . , 1*(lin-
early decreasing) for both this model with long-range interaction and
*our previous nearest-neighbor model. We observe again that D(s)*
*displays power-law-like behavior: log D/ log s ≈ 1/3 for the nearest-*

*neighbor model with increasing η**n**; log D/ log s ≈ 1/2 for the nearest-*
*neighbor model with decreasing η**n* and for the long-range model
*with increasing η**n**; and log D/ log s ≈ 2/3 for the long-range model*

*with decreasing η**n*. 10^{-2}

10^{-1}
10^{0}

10^{-6}10^{-5}10^{-4}10^{-3}10^{-2}10^{-1}10^{0}

D

s (norm.)

Figure 24: As Figure 23, but with logarithmic axes.

### Chapter 4

## Rupture under uniform loading

In the previous chapters we applied loading to the slider by pulling it
from one side. Here we will consider a slightly different situation: by
*applying the loading through a rigid ‘bulk’ that is coupled to all blocks*
. This model is illustrated in Figure 25b. The springs between the
*interfacial blocks and the bulk have spring constant p**k* *= K/k, and*
*we will use p**k*=1 unless noted otherwise. For the internal elasticity
of the slider we consider three cases:

*T*
*F*1 *F*2 *F*3

*T*
*F*1 *F*2 *F*3

(a)

(b)

Figure 25: Cartoon of the model considered in the previous chapters (a) and of that considered in this chapter (b).

Horizontal (nearest-
neighbor) springs have
*spring constant k, while*
diagonal springs have
*spring constant K.*

• Nearest-neighbor interaction: springs of stiffness 1 between neighboring blocks along the interface;

*• Spring constants that decrease as 1/r, where r is the distance*
between two blocks: springs of stiffness 1 between neighbor-
ing blocks, stiffness 1/2 between next-nearest neighbors (for
example, blocks 1 and 2), stiffness 1/3 between third-nearest
neighbors (for example, blocks 1 and 4), and so on;

Figure 26: Vertical cross
sections of these plots
indicate which contacts
have slipped (black pix-
els) as a function of
applied load for three
spring models: nearest-
*neighbor (left), 1/r (cen-*

*ter) and 1/r*^{2}(right). Load

0 1/4 1/2 3/4 1

Load

0 1/4 1/2 3/4 1

Load

0 1/4 1/2 3/4 1

*• Spring constants that decrease as 1/r*^{2}: springs of stiffness
1 between neighboring blocks, stiffness 1/4 between next-
nearest neighbors (for example, blocks 1 and 2), stiffness 1/9
between third-nearest neighbors (for example, blocks 1 and
4), and so on;

### 4.1 Almost uniform friction forces

For these cases of elastic coupling along the interface, we have done
*calculations of systems of 500 contacts, where all critical forces η**n*

were set to the same non-zero value, except at 50 contacts where
*we set η**n* = 0, so that these contacts slipped immediately. We
*increased the tangential load T from 0 to 500 in 500 steps. The*
results are shown in Figure 26, where a black pixel at any applied
load indicates that the corresponding contact has slipped. The total
number of slipped contacts as a function of applied load is plotted
in Figure 27.

0 1/4 1/2 3/4 1

0 1/4 1/2 3/4 1

Load

Fraction of slipped contacs

Nearest−neighbor 1/R 1/R2

Figure 27: Number of
slipped contacts as a
function of applied load
*for uniform F**n*, where
10% of contacts was set
*to slip immediately (F**n*=
0).

In models where there is a stronger coupling between neighboring
*contacts (1/r, 1/r*^{2}, and nearest-neighbor, in decreasing order of
interaction distance), we see that we can force a rupture front to
nucleate at these sites. There is no sudden slip of a large number of
*contacts, except when the applied load T is very close to*∑

*F**n*.