Macroscopic friction from the collective slip of contacting asperities

(1)

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

(2)

(3)

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,19Friction 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.⁴The real contact area between two surfaces, how-

(6)

ever, is generally much smaller than it appears to be at a macroscopic 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.⁴

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,25as 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.¹³Experiments have also shown that similar ‘rupture fronts’ traverse the interface on a macroscopic scale.²These fronts initiate at the loaded edge and propagate at velocities up to the shear wave speed. Behind a front the con-

(7)

tact area is reduced.²⁷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.²⁸

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

Signiﬁcant variations in both the normal and shear stresses are reported 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.³ It is speculated that this redistribution of stresses induces ‘local’

Amontons-Coulomb friction.²⁹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.¹⁸

1 2 3 N

K k F

v ¹

P m

Figure 4: Mechanical model of frictional sliding.

Various numerical studies have been able to reproduce these precursor 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,31The studies cited, however, all apply Amontons-Coulomb friction to the blocks individually: no- tably, they distinguish between kinetic and static coefﬁcients of fric-

(8)

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 consider 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 particularly 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,9In 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.

(9)

Chapter 2

Dynamics of interface

rupture

¹ ² ³ ^N

K k F

v ¹

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 kbetween 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 pk := K/k and pv:= v√

mk/F. We will, like the studies cited in the previous chapter^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^Sand friction forces F_n^F. The equation of motion for the n-th block is there- fore

m¨un= F_n^S+ F_n^F (3)

(10)

where ¨un := d²un/dt² 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₁^S = −(k + K)un+ Kvt + kun+1 and F_N^S =

k(u_n−1− un)instead. F_n^S= k(un−1− un) + k(un+1− un). (4)

The ﬁnal term in equation (3),

F_n^F=: F η(F_n^S/F, ˙un), (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 friction 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.

(11)

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





 F₁^S

... F_N^S







= k







−1 − pk 1

1 −2 1

...

1 −2 1

1 −1











 u1

... uN







, (8)

where pk := K/kis 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₁^S/F, ˙u1) η(F₂^S/F, ˙u2)

... η(F_N^S/F, ˙uN)





 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 =: ucxand t =: tcs, where x and s are nondimensional and uc and tc are some characteristic length scale and time scale for the problem, then equation (11) becomes

muc

F t²_c d²

ds²x =kuc

F κ· x + Kvtc

F se + η. (12)

This result indicates that it might be fruitful to set uc := F/k and t_c :=√

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

(12)

can be rewritten as d²

ds²x = κ· x + pvpkse + η, (13) where pv := v√

mk/F is the other dimensionless parameter we identiﬁed before. Although only the product pvpk of the nondimen- sional parameters appears explicitly in this equation, pk 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=: uc/v, so that the time scale of the problem is expressed with respect to an external variable. Then equation (13) changes to

p²_v d²

dζ²x = κ· x + pkζe + η. (14)

We can go from a solution of equation (13) to a solution of equa- tion (14) by substituting pvsfor ζ. 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 pv vanish from equation (13) and equation (14): as pv 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 pv 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

(13)

is constant for the blocks that are sliding, equation (13) reduces to a system of linear ordinary differential equations similar to ω⁻²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 first consider the case when only the first block is moving, so that xn(s) = 0for n > 1. The first block slips when x^′′₁ = p_kp_vs−1 ≥ 0, which implies s ≥ 1/p^kpv. Its equation of motion is

x^′′₁ =−(1 + pk)x1+ pvpks− 1, (16) with x1= 0initially, in which we substitute ˆs = s − 1/p^vpkto get

x^′′₁ =−(1 + pk)x1+ pvpkˆs. (17) The general solution of this equation is

x1(ˆs) = a cos(√

1 + pks) + b sin(ˆ √

1 + pkˆs) + pvpksˆ 1 + pk

. (18)

The initial condition x1(ˆs = 0) = 0gives a = 0 and x^′₁(ˆs = 0) =

√1 + pkb + pvpk/(1 + pk) = 0gives b = −p^vpk/(1 + pk)^3/2, which implies

x1(s) = pvpks− 1

1 + pk − pvpk

(1 + p_k)^3/2sin(√

1 + pk(s− 1/pvpk)), (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.

(14)

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

1 + pk, so the two are approx- imately of the same magnitude if pk ≪ 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^′′₂ = x1(s)−1 ≥ 0, so when x1(s)≥ 1.

Using the substitution ˆs =√

1 + p_k(s− 1/pvp_k)we get 1 = p_vp_k

(1 + pk)^3/2(sin(ˆs) + ˆs). (20) There is no closed-form solution to this equation, but numerical methods 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 s0

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

x(ˆs) =− ˆκ⁻¹n · (p(ˆs + s0)e− 1) +

∑n k=1

[c⁺_k exp(+iωks) + cˆ _k⁻exp(−iωkˆs)]ξ_k, (21)

where p = pvp_k and ξ₁ through ξ_n are eigenvectors of the n × n matrix

ˆ κ_n=







−1 − pk 1

1 −2 1

...

1 −2





, (22)

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

(15)

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

κthrough ω²_k=−λk.

If x0= x(ˆs = 0)is the initial position and x^′₀ = x^′(ˆs = 0)the initial velocity of the blocks, and Ξ is a matrix whose m-th column is ξm, then the coefﬁcients c^±_mcan be determined as (see the Appendix)

c⁺+ c⁻= Ξ⁻¹· [x(s0) + ˆκ⁻¹· (ps0e− 1)] (23) c⁺− c⁻= i Ξ⁻¹· ω⁻¹· [x^′(s₀) + ˆκ⁻¹· e], (24) where

c⁺=





 c⁺₁

... c⁺_n





 and c⁻ =





 c⁻₁

... 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 xn= 1, or more explicitly, when

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

=− [0 . . . 0 1] · ˆκ⁻¹· (p(ˆs + s0)e− 1) + [0 . . . 0 1]·

∑n k=1

[c⁺_k exp(+iωks) + cˆ ⁻_k exp(−iωks)]ξˆ _k. (26)

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 motion for the entire system by

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

(16)

(a) for the ﬁrst block s = 1/pkpv,

(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⁻⁴/√

1 + pk,^⋆ and the initial

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

1 + pk.

conditions were x = x^′ = 0 and s = 1/pvpk− 10∆s, so that the ﬁrst block started sliding after ten time steps. The calculation ended when x5> 10⁻⁴.

The displacements, for pv =10^–1, 1, 10⁺¹and pk =10^–1, 1, 10⁺¹, 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 pv and pk

from 10^–3through 10⁺³.

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

In the fully numerical results when pv < 1or pvpk < 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ⁱ_n^′ = 0at time step i, it is xⁱ⁺¹_n ^′ > 0in the next time step. This ‘stick-slip’ cycle will then repeat for some time until the block starts purely sliding again.

(17)

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 pv=10^–1, 1, 10⁺¹(left to right), and pk =10^–1, 1, 10⁺¹(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’

0 1 2 3

0 5 10 15

s

x’

0 50 100 150

0 0.05 0.1

s

x’

0 5 10

0 0.2 0.4 0.6 0.8 1

s

x’

0 1 2 3 4

0 2 4 6 8 10

s

x’

0 200 400 600

0 0.01 0.02 0.03 0.04

s

x’

0 20 40 60

0 0.1 0.2 0.3 0.4

s

x’

0 2 4 6 8

0 1 2 3

s

x’

Figure 11: Velocity as a function of time for pv= 10^–1, 1, 10⁺¹(left to right), and pk = 10^–1, 1, 10⁺¹ (bottom to top).

(18)

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 pv < 1 or pvpk < 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₁₀p_v log10pk

Stick−slip

Pure sliding

−5 0 5

Figure 12: Phase dia- gram showing the transition 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_vp_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(pv/p^∗_v, p_vp_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

(19)

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 nondimensional groups scale as

p_k= K k ∝ 1

N and pv =v√ mk

F ∝ N. (27)

This implies that pvp_k is independent of N. Furthermore, both pv∝ N and p^∗_v ∝ N, so the criterion for stick-slip min(pv/p^∗_v, pvpk) < 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

pvpk= 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). ⁻¹⁻⁵ ⁰ ⁵

0 1 2 3 4

log10ω

log₁₀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 pk⋆: the frequencies of the ﬁrst

⋆Recall that the eigenfrequencies are related to the eigenval- ues of κ, which does not de- pend on pv.

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.

(20)

(21)

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(pv, pvpk), whereas the greatest characteristic frequencies are around max(1, pk).

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 pv= 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 + pkζe + η. (29) This is a piece-wise linear system consisting of N equations with 2N + 1degrees of freedom: all ηn, all nnand ζ. The piece-wise behavior is caused by the fact that the absolute value of each element

(22)

of η by deﬁnition satisﬁes |ηⁿ| ≤ 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 − pk 1

1 −2 1

...

1 −2











 x1

x2

... xn





+





 pkζ

0 ... 0





−





 1 1 ... 1





. (30)

This is equivalent to





 1 1 ... 1





=







−1 − pk 1 1

1 −2 1 0

... ...

1 0











 x1

... xn−1

pkζ





+





 0

... 1

−2





x_n. (31)

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





 x1

... xn−1

pkζ0





=







−1 − pk 1 1

1 −2 1 0

... ...

1 0







−1









 1 1 ... 1





−





 0 ... 1

−2











. (32)

The first n − 1 elements of this array are the positions xm(ζ₀)of the first n − 1 blocks at ζ = ζ0, so by setting the final element to xn= 1 we get x(ζ0). For example, the first block slips at ζ = 1/pk, the second block at ζ = 1 + 2/pk and the third at

[ x1

p_kζ₀ ]

=

[−1 − pk 1

1 0

]₋₁[[

1 1 ]

− [

1

−2 ]]

= 3 [

1 1 + p_k

] . (33)

The variable ζ = pvs = tvk/F was introduced as a dimensionless

(23)

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 pv=10^–1, 1, 10⁺¹(left to right), and pk =10^–1, 1, 10⁺¹(bottom to top).

‘time’ in the previous chapter. In the quasistatic limit it is more properly 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 pv, pk =10^–1, 1, 10⁺¹, 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_vsfor ζ.

displacements are plotted for greater ranges of pv and pk, 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 pv≤ 1.

log₁₀p_v log10pk

Quasi-static

Dynamic

10%

1%

−5 0 5

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

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

(24)

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 |∆x1/x₁|, where x1is the displacement of the ﬁrst block at the time when the N-th block slips, as calculated using the dynamic model, and where ∆x1

is the difference between x1 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 pv≪ 1 or pvp_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 pv ∝ N − 1 it seems safe to assume that the condition pv≪ 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

pvpk= vK F^(tot.)

√m^(tot.)

k^(tot.), (34)

which was also the deﬁning parameter in stick-slip versus pure sliding. 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

(25)

Figure 16 shows the number of sliding blocks as a function of time for systems of N = 2, 4, 8, …, 1024 blocks, with pk =1. Both axes are normalized for easy comparison; it bears to mention that 0 on the horizontal axis here corresponds to s = 1/pkp_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⁻² 10⁻¹ 10⁰

s (normalized)

# sliding blocks (norm.) N=1024

N=512 N=256

N=128 N=64

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

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

Using N = 100 and different values of pk, 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 pkincreases, which is illustrated in Figure 18: the solid lines (for pk≥ 1) converge to√

s.

On the other hand, if we let pkdecrease, 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 pk=10^–6, 10^–5, …, 10⁺⁶ (dotted lines pk<1).

This dependence on pk 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 pk ∝ 1/N. This, combined with the ﬁndings in Figure 18, suggests that in the continuum limit

lim

N→∞D(s)= lim^?

pk→0D(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₁₀p_k

log D / log s

N=10 N=100 N=1000

Figure 19: As pk increases, 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 pk = 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 pk ≫ 1, and that there is a transition between log D/ log s =1 and log D/ log s = 1/2 when 10^–2< pk<10⁺³.

(26)

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 pk≪ 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

D

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

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

10^-2 10^-1 10⁰

10^-4 10^-3 10^-2 10^-1 10⁰

D

s (norm.)

Figure 21: As Figure 20, with logarithmic axes.

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

|Fn^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 critical 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

(27)

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, η

0 0.2 0.4 0.6 0.8 1

0 20 40 60 80 100

dD/ds, η

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

(red) as a function of distance along the sliding interface. 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 rupture 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 between neighboring blocks.

(28)

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

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¹⁶.

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(linearly 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-

(29)

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⁰

10^-610^-510^-410^-310^-210^-110⁰

D

s (norm.)

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

(30)

(31)

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 pk = K/k, and we will use pk=1 unless noted otherwise. For the internal elasticity of the slider we consider three cases:

T F1 F2 F3

(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 neighboring 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;

(32)

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²(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²: 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

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 Fn, where 10% of contacts was set to slip immediately (Fn= 0).

In models where there is a stronger coupling between neighboring contacts (1/r, 1/r², 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∑

Fn.

Macroscopic friction from the collective slip of contacting asperities