Traveling wave solutions in the Burridge-Knopoff model

12 Pages·0.23 MB·English

by  C. B. Muratov

Traveling wave solutions in the Burridge-Knopoff model
C. B. Muratov
Courant Institute of Mathematical Sciences, New York University,
251 Mercer St., New York, NY 10012
(February 9, 2008)
The slider-block Burridge-Knopoff model with the Coulomb friction law is studied as an excitable
medium. It is shown that in the continuum limit the system admits solutions in the form of the
self-sustained shock waves traveling with constant speed which depends only on the amount of the
accumulated stress in front of the wave. For a wide class of initial conditions the behavior of the
system is determined by these shock waves and the dynamics of the system can be expressed in
termsoftheirmotion. Thesolutionsintheformoftheperiodicwavetrainsandsourcesofcounter-
propagating waves are analyzed. It is argued that depending on the initial conditions the system
9
will either tend to synchronizeor exhibit chaotic spatiotemporal behavior.
9
9
PACS Number(s): 83.50.Tq, 91.30.Mv, 83.50.By, 03.40.Kf
1

n
a
J I. INTRODUCTION porting steadily propagating solitary waves in the form

4 of shocks [7–10]. These systems are also of special in-
1 Propagating self-sustained waves (autowaves) and terest because they are used for modeling the dynamics

more complex spatiotemporalpatterns arecharacteristic of earthquakes[7–17]. It is clear that systems exhibiting
1
v of excitable media of different nature. A typical exam- stick-slip motion have both necessary ingredients of ex-
3 ple of such a phenomenon is burning of black powder citable systems. The threshold behavior here is due to
0 in a safety fuse. When the fuse is ignited at one end, static friction, which prevents any motion in the system
0 the exothermic reaction releases the heat which is then until some critical amount of stress is accumulated. The
1 spread out by heat diffusion. Thus, the neighboring re- couplingoftheelementsofthe systematdifferentpoints
0 gionsofthe fuse igniteleadingto self-sustainedpropaga- in space is due to the non-locality of elastic stress.
9
tion of the combustion front. The phenomenon of non- Singularperturbationtechniques provedto be veryef-
9
/ attenuated propagation of waves is in fact common for fectiveintreatingproblemsoftravelingwavepropagation
ol a variety of physical, chemical, and biological systems in reaction-diffusionsystems [1–6,18,19]. These methods
s [1–6]. Traveling waves are experimentally observed in usestrongseparationoftimescalesintheproblemtode-
- semiconductorandgasplasma,semiconductorandsuper- compose the dynamics of the system into fast and slow
t
at conductorstructures,combustionsystems,activeoptical motion. Clearly,thissituationisalsorealizedinthemod-
p media,magneticmediaunderillumination,autocatalytic elsofstick-slipmotionwhere(especiallyinthecontextof
: chemical reactions, nerve and heart tissue (see [1–6] and earthquakes) there is a strong separation of time scales
v
i references therein). between fast slipping events and slow accumulation of
X Inorderforself-sustainedwavestobefeasible,thesys- stress. It is therefore advantageous to try to apply these
r tem must possess two basic ingredients. First, the sys- techniques to the problem of stick-slip motion.
a tem must be excitable, that is, there has to be a thresh- In this paper we present a study of the Burridge-
old below which the perturbation of the steady homo- Knopoff slider-block model [11] with the Coulomb fric-
geneous state of the system decays, while perturbations tion law. We will show that for sufficiently slowly spa-
of larger amplitude grow. In the example above, a suf- tially varying displacement variable the dynamics of the
ficient amount of heat is needed to ignite black powder. system is dominated by self-sustained traveling shock
Second, there has to be coupling between the regions waves. We will study the properties of these waves and
of the system at different points in space. In the case reformulatethe dynamics of the systemin terms oftheir
of the safety fuse such a coupling is provided by heat motion.
diffusion leading to spread of the temperature and ig- Our paper is organized as follows. In Sec. II we intro-
nition of powder in front of the combustion zone. Thus, ducethegoverningequationsforthemodelwestudyand
theprototypesystemsexhibitingself-sustainedwavesare discuss the features of the friction law used, in Sec. III
reaction-diffusion systems [1–6]. we construct the solutions in the form of self-sustained
Recently, it was pointed out that an entirely differ- traveling shock waves, in Sec. IV we reformulate the
ent class of systems may be considered as excitable [7]. dynamics of the system in terms of the motion of these
Theseareelasticmediawithfrictionexhibitingstick-slip shock wavesand study generalproperties of the reduced
motion. Both experimental observations and numerical problem, in Sec. V we analyze two different types of
simulations show that such systems are capable of sup- solutions, and in Sec. VI we draw conclusions.
1
II. MODEL f
The Burridge-Knopoff model consists of a one-
∼αv
dimensional array of blocks of mass m resting on a fric-
f
tionalsurface[11]. Theblocksareconnectedtogetherby r
springs with spring constant k and pulled by a loader
c
plate moving with constant speed V via another set of
f
springs having spring constants k (see Fig. 1). Let us s
p
V
a a
v
k k k
p p p FIG.2. Friction force as a function of thevelocity
k k
c c that this model was recently studied numerically in Ref.
m m m [17] (although for the parameters different from those
(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)
considered in this paper).
(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)
Wewouldliketoemphasizethatincontrasttovelocity-
(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)
FIG.1. The Burridge-Knopoff model weakening friction law, for which the kinematic friction
force is a continuous function of the block velocity at
measurethedisplacementXi ofthe i-thblockrelativeto v = 0, in our friction law the force f has a discontinu-
the point ofattachmentof the i-thloaderspring. In this ity at v = 0, reflecting the fact that the coefficient of
case the total force Fi acting on the i-th block is given static friction is typically larger than that of kinematic
by friction. As a result of this difference, the accumulated
stressacceleratestheblockattheonsetofmotion,releas-
Fi =kc(Xi+1+Xi−1 2Xi) kpXi fi, (1) ing potential energy. As we will see below, this provides
− − −
the sustaining force for the traveling waves studied in
wheref istheforceoffriction. Thedynamicsofthesys-
i the following section. Also, as we will show below, this
tem is completely determined by the equation of motion
ensures the existence of the proper continuum limit as
mX¨ = F , provided that the friction law is specified.
i i a 0.
Note that the friction force f is the only nonlinearity →
i Let us now formulate the equations of motion for the
in the equation of motion. Clearly, the dynamics of the
blocks in the case of the friction law introduced above.
system will significantly depend on the particular choice
RecallthatthedisplacementsX aremeasuredrelativeto
i
of the friction law. Recently, a lot of results were pre-
theloaderplate,sotheyaredefinedinthereferenceframe
sentedonthedynamicsoftheBurridge-Knopoffmodelin
movingrelativetothesurfacewithvelocityV. Therefore,
the case of velocity-weakening friction law [9,10,12–16].
when the blocks are at rest, their equation of motion
A characteristic feature of the Burridge-Knopoff model
becomes
with this form of friction is its highly chaotic dynam-
ics that occurs on all length scales down to the smallest X˙ = V. (2)
i
−
length scale a and therefore, the absence of the proper
continuum limit [12]. When the block slides,its equationofmotionchangesto
In contrast to most previous studies, here we adopt
the Coulomb friction law, which is applicable to clean mX¨i =kc(Xi+1+Xi−1−2Xi)−kpXi−fs−α(X˙i+V).
dry surfaces (see, for example, [20]). Namely, we will (3)
characterize the friction force by the value of the static
friction fr and the kinematic friction fs < fr, where fr The transition to sliding (slip) occurs when the force Fi
andfsarepositiveconstants. Thus,theblockwillremain reaches the value of fr:
atrestifF <f . WhenF reachesf theblockstartsto
move (slipsi) andr the frictioin force drrops instantaneously kc(Xi+1+Xi−1−2Xi)−kpXi =fr. (4)
from the value of fr to the value of fs. When the block Inotherwords,theequationofmotionofablockchanges
comestorest(sticks),staticfrictionturnsonagain. Also, from Eq. (2) to Eq. (3) when the condition in Eq. (4) is
we will supplement the friction force with the viscous satisfied. Of course, the block comes back to rest when
friction term fviscous(v) = αv, where v is the velocity of X˙i = V during sliding.
the block and α is a constant (Fig. 2). Note that this Let−usintroducethefollowingdimensionlessquantities
kindoffrictionlawwasconsideredalreadyintheoriginal
paperofBurridgeandKnopoff[11]andisobservedinthe x k k k X +f
′ p ′ p p i s
x = , t =t , u = , (5)
frictionexperimentsandtheiranalogs[21–25]. Alsonote i
a k m − f f
r c r r− s
2
whereaisthedistancebetweentheattachmentpointsof what expresses the fact that the time scale of accumula-
the loader springs (Fig. 1). The dimensionless displace- tion of stress is much longer than that of the motion of
ment variableu is chosenin sucha waythat the system an individual block during sliding.
i
is excitable only at u > 0 [see Eq. (1), in order for slip Finally,letusdiscusstheapplicabilityoftheBurridge-
i
to occur,the force F mustexceed the value f of sliding Knopoff model in the context of real physical systems
i s
friction],while the blockswillalwaysslide foru >1[see exhibiting stick-slip motion. In a realsystem one should
i
Eq. (4)]. replace a one-dimensional array of masses between the
Ifk k ,onecannaturallygotothecontinuumlimit loader plate and the frictional surface by an elastic
c p
bTyhimsaisk≫ianggoaosdubasptpitruotxiiomnaXtioi+n1w+hXeni−t1h−e l2eXngit→h sac2a∂l∂e2xX2of. msioendiuofmthoef cBerutrariindgteh-iKcknnoepsosffh.mAodsetlrawigohutlfdortwhaerrdefoerxetebne-
variation of X is of order 1 in the new variables. Using a two-dimensional array of masses connected by springs
i
the variables in Eq. (5), we rewrite Eqs. (2) and (3), whose one edge is rigidly attached to the loader plate
respectively, as follows: and the other slides on the frictional surface. Naturally,
thethicknessofsuchamediumshouldgreatlyexceedthe
u =v, (6) microscopiclengthscalea. The stick-slipmotionin such
t
a medium is due to surface waves [26]. The dispersion
u =u u 2γ(u v), (7)
tt xx− − t− of these waves is given by ω2 = k2 +(π/2h)2(1+2n)2,
where we introduced dimensionless constants where n is an integer and the transverse speed of sound
was taken to be 1. The main difficulty here is that in-
α V kpm stead of a single displacement variable u one has to deal
γ = , v = , (8)
2 kpm frp−fs wficitahtioanlatrhgaetnisumprbeesrenotfemdobdyetshweiBthurdriiffdegree-nKtnno.pAoffsmimopdlei-l
and dropped the pprimes for simplicity of notation. Note consists of lumping up these modes into a single mode
that in the continuum limit in addition to tracking the withaneffectivestiffnesskp [Eq. (3)]. Itisclearthatthe
motion of individual blocks, one also has to follow the dominant modes will be those with small n, so we must
motionoftheslippoints. Similarly,oneshouldkeeptrack have kp (a/h)2kc kc. Thus, in a physically rele-
∼ ≪
ofthepositionsofthepointsatwhichtheblockscometo vant situation one should consider the continuum limit
rest. The latter are determined by the condition in Eq. ∆x 0 of the Burridge-Knopoffmodel. In short, this is
→
(6) during sliding. merely a simplification of the mathematical handling of
Determining the position of the slip points in the con- the elastic dynamics. A much more important assump-
tinuum limit turns out to be a rather complicated prob- tion is that the friction response to the motion of the
lem, since the positionofthe slipwillstill dependonthe mediumisinstantaneousandhaszerocorrelationlength.
local dynamics of the blocks. In the discretized form, It would be incorrect to think of the “size” of the block
the slip condition [Eq. (4)] in the variables of Eq. (5) a as an atomic distance. Rather, the value of the pa-
becomes rameter a should have a magnitude of the correlation
(memory) length of the cooperative effects responsible
ui+1+ui−1−2ui =u 1, (9) for the surface friction (see, for example, [21] and refer-
(∆x)2 i− ences therein). Then the response of the system can be
considered to be instantaneous if the characteristic ve-
where ∆x = kp/kc. We will get back to this problem locity of the blocks is greater than a/τstick, where τstick
in the following section where we will show that in the is the characteristic sticking time.
p
continuum limit the dynamics of the slip becomes inde-
pendent of ∆x.
Equations(6)–(9)arethe constitutiveequationsthat III. TRAVELING WAVES
will be studied in the rest of the paper. As can be seen
from Eq. (5), we chose time and length scales so that
LetusnowdemonstratethatEqs. (6)–(9)admitsolu-
the speed of sound in the system is equal to 1. The time
tions in the form of traveling shock waves with constant
scaleisdeterminedbytheperiodofoscillationsofaniso-
speed c in the limit of vanishingly small v. But before
lated block due to the loader spring. Note that in the
we do that, let us see what kind of dynamics one would
continuum limit the characteristic speed of a block dur-
observe if there are no slip events and the initial distri-
ingslidingismuchsmallerthanthespeedofsound. This
bution of u is taken to be a sufficiently slowly varying
canbeseenfromEq. (5)ifoneassumesu 1,measures
t ∼ distribution u0(x) (the latter ensures that no slip events
X in the units of Eq. (5), and uses a natural assump-
i will occur at short times). In this case the dynamics is
tion that f f ak . The behavior of the system is
r − s ≪ p trivial: we will have u(x,t) = u0(x)+vt, as long as no
determined only by two parameters: the dimensionless
slipeventsoccur[seeEq. (6)]. Fromthisonecanseethat
dissipationparameterγ whichmeasuresthe effectofvis-
the characteristic time scale of accumulation of stress is
cousfriction,andthe dimensionlessrateofaccumulation v−1 1, when v is small. In other words, for v 1 on
ofstressv. Inthefollowingwewillconsidervtobesmall, ≫ ≪
3
thetimescaleoforder1onewillnotseeanymotionwith 3
2.8
thedistributionofufixedtou (x)+C,whereC isacon-
0
2.6
stant. Inparticular,ifonestartswiththe uniforminitial
2.4
conditions, one will have u=u , where u is some con-
+ + 2.2
stant less than 1. The other possible situation when the
c 2
dynamics of the system becomes trivial is steady creep,
1.8
when all the blocks move together with the loader plate,
1.6
so we simply have u=2γv 0 for v 1 [see Eq. (7)].
1.4
≃ ≪
Consider the system with u = u and v = 0. As was
+ 1.2
noted in the previous section, when u+ < 0 slip events 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
cannot occur, so if one slightly moves a particular block,
u
+
it will quickly move to a new equilibrium position and
FIG. 3. The dependence c(u+). Results of the numerical
no significant changes will occur. A different situation is
solution ofEqs. (A3)–(A7). Dashedlineshowstheresultof
realized for 0 < u < 1 when the system becomes ex-
+ approximation by Eq. (11).
citable. Then, if a single block is slightly moved, it will
accelerate,releasingthe potentialenergyaccumulatedin
u . As a result of this acceleration, the force acting Toactuallycalculatethespeedcofthefront,oneneeds
+
on the adjacent blocks will increase resulting in slipping to solve the system of coupled equations of motion for
of the adjacent block that was at rest. This avalanche- the blocks behind the s-th block. This problem in the
like process will go on, leading to the formation of the limit ∆x 0 is considered in Appendix A. The analy-
→
shock wave that releases the accumulated stress, so that sis of this problem shows that the speed c of the wave
u changes to a new lower value of u−. The situation is uniquely determined by the value of u = u+ in front
hereresemblesagreatdealcombustionofsafetyfusedis- of the wave. Note that similar situation takes place in
cussedintheIntroduction. Thus,for0<u <1asmall the models of combustion (see, for example, [18]) and in
+
perturbationof u may leadto a significantchangeof the reaction-diffusion systems (see, for example, [1–6]). The
state of the system,switching itfrom u=u+ to u=u−. dependence c(u+) found numerically is shown in Fig. 3.
Asinotherexcitablesystems[1–6],itisnaturaltoexpect We would like to emphasize that upon the decrease of
thatsuchaperturbationwilltransformatlongtimestoa ∆x the speed c becomes independent of ∆x, providing a
shockwavetravelingwithconstantvelocityc. Numerical well-defined continuum limit for slip propagation. These
simulations of Eqs. (6) – (9) with v = 0 show that this results are also supported by the direct numerical simu-
indeed happens in our model. lations of Eqs. (6) – (9).
Sinceallthe“nonlinearity”intheproblemiscontained According to Fig. (3), traveling wave solution exists
in the slip condition given by Eq. (9), our system is not for all0<u+ <1. The speed c is alwaysgreaterthan 1,
unlike piecewise-linear model used to study propagation that is, the considered shock waves are supersonic. The
of nerve impulses [27]. For this reason the solution in latterisalsoobservedinothermodelsofstick-slipmotion
the form of the traveling wave can be found exactly for [15,16]. Note, however, that the speed c is the speed of
v = 0 in the continuum limit. For definiteness, we will propagation of the slip andis not relatedwith the actual
considerthewavestravelingfromlefttoright,sothatc> speed of individual blocks in the wave.
0. Because of the reflection symmetry, for each solution As can be seen from Fig. (3), the speed of the travel-
traveling with speed c there is also a solution traveling ing wave diverges as u+ approaches 1. One can get an
with speed c. analytical handle on the dependence c(u+) by expand-
−
Itisclearthatthespeedcofthetravelingwaveshould ing it in the inverse powers of c (see Appendix A). As a
be determined by the motionof the slip point in frontof result, we get the following interpolation formula which
the wave. Let the slip event occur for block s at t = ts. implicitly determines c as a function of u+
Then, for i > s we have ui = u+, so the slip condition 2+c2+4c4 −1
from Eq. (9) becomes
u = 1+ . (11)
+
8c5√c2 1
us−1(ts)=u+−(1−u+)(∆x)2. (10) Thisequationgive(cid:18)sthedepende−nce(cid:19)c(u+)withafewper
cent accuracy.
Since atthe onsetofthe slipu =u anddu /dt=0,at
s + s
timetwewillhaveu =u + (t t )2 ,withu <u , Letus introduceself-similarvariablez =x ct, where
s + s s + −
O − c=c(u ) (Fig. 3). Since the problem possesses transla-
since right after the slip the block is accelerated by an +
(cid:0) (cid:1) tional invariance, we can choose the position of the slip
excess force. This means that at some time t the slip
s+1
point to be at z = 0. Similarly, the stick point z = w,
condition from Eq. (9) will be satisfied for the (s+1)-
−
withw >0,willalsotravelwithspeedcbehindthewave.
th block. According to Eq. (10), this will happen when
Behind the slip the blocks slide according to Eq. (7), so
∆t = t t = (∆x). Therefore, on the time scale
s+1 s
− O for the wave with speed c we will have
of order 1 this will correspond to the motion of the slip
point with the speed c=∆x/∆t= (1). (c2 1)u 2cγu +u=0. (12)
O zz z
− −
4
u u
a) 1 u b) 1 u
+ +
c c
-w
0 z 0 z
u
-
-1 -1
FIG.4. Two typesof traveling wavesolutions: (a) case u+ <uc; (b) case u+ >uc.
The slip condition gives the following initial conditions
for this equation: πγc
u− = u+exp . (16)
u(0)=u , u (0)=0. (13) − − c2(1 γ2) 1!
+ z − −
p
TheanalysisofEq. (12)(seeAppendixB)showsthatthe From this equation one can see that u− < 0, so right
structureofthesolutionchangesqualitativelydepending behind the wave the system is in the state in which no
on whether the value of u+ is smaller or larger than the waves can be further excited. Also, note that when u+
critical value uc > 0, where the value of uc is implicitly approaches the value of uc, the value of u− rapidly ap-
determined by proaches zero.
In the context of earthquakes a shock wave should be
c2(u ) 1 associatedwith an individual earthquake. The total dis-
c
γ = − . (14)
p c(uc) placement A=u+−u− which occurs upon passing of a
wave and which determines the magnitude of the earth-
If we have u+ < uc, the distribution of u will asymp- quake will therefore be
totically approach u− = 0 at z = behind the wave
−∞
[Fig. 4(a)]. πγc
Inotherwords,theblocksbehindthewavenevercome A=u+ 1+exp . (17)
" − c2(1 γ2) 1!#
to rest. Therefore, upon passing of the wave the system − −
shoulddevelopasteadycreepforsmallbutfinitev. Note p
that this will always happen when γ > 1 since in any
case u < 1, so for these values of γ the propagation of
+ IV. FREE BOUNDARY PROBLEM
multiple shock waves becomes impossible (see Sec. V).
Adifferentsituationtakesplaceforu >u . Thenthe
+ c
Sofarwewerestudyingtravelingwavesolutionsinthe
solution has the form shown in Fig. 4(b), and the width
limit v = 0 and u = const in front of the wave. These
w of the travelingwave becomes finite (see Appendix B)
solutions are clearly the fast motions in the systemsince
π(c2 1) they occur on the time scales of order 1. On the other
w = − . (15)
hand, although the accumulation of stress occurs slowly
c2(1 γ2) 1
− − for v 1, it can lead to significant changes in u at long
Thus,inthiscasethepwavepropagatesasa“self-healing” times≪of order v−1 1, so the latter can be associated
≫
pulse (compare with [28]). Note that for γ > 0 we have with the slow motions in the system. Therefore, since
c > c(u ) > 1, so the value of w is bounded from below the fast motions result in the large-scale changes in u,
c
andthewidthofthewaveisalwaysgreateroroforder1. the slip-stick events can be considered as singular per-
This justifies the use of the continuum limit for k k turbations to the slow motions.
c p
≫
and γ 1. Furthermore, when u becomes close to 1, Similar situation is realized in other excitable media
+
∼
the propagation speed c and therefore the width of the where the role of singular perturbation is played by the
wavegrow,sointhiscasethecontinuumlimitisjustified sharp fronts [1–6]. A powerful tool in describing the dy-
even for k k . On the other hand, the width w goes namics of such systems is reduction to free boundary
c p
∼
to zero when u becomes small when γ 1 [see Eq. problem (see, for example, [18]). In this approach one
+
≪
(15)], so for sufficiently small u and γ the continuum obtains the relationship between the speed of the sharp
+
limit will become invalid for the description of traveling frontandslowvariables. Usingthisidea,letusformulate
waves. the dynamics of our system as a free boundary problem.
In the case u+ > uc the value of u = u− behind the Let us assume that the distribution of u varies on the
wave will be (see Appendix B) length scale much greater than 1. Then, introducing an
5
u
adiabatic approximation, we can assume that the speed
c
of the traveling wave is given by the dependence c(u )
+
(Fig. 3) in which now, instead of u one should use an
+
instantaneous value of u right in front of the wave. The
latter will play the role of the slow variable. Thus, the
motion of the shocks can be described by the variables
x
x (t) which give the positions of the i-th shock, and the
i
variables s = 1 which give the direction of their mo-
i
±
tion.
On large length scales the shock will contribute a dis-
continuity to the distribution of u in the limit v 0. It
→
canbeincludedinEq. (6)describingslowmotionsinthe FIG.5. Periodic wave train.
form of a δ-function. Therefore, in the presence of the
shocks Eq. (6) can be rewritten as
The condition of the absence of creep should in fact be
satisfiedfor a wide classof initial conditions. In particu-
u =v c A(u)δ(x x +0s ), (18)
t − i − i i lar, it is easy to see that a wave with u+ >uc will never
Xi be able to reach the points where u<uc if we have
whereAistheamplitudeoftheshockfromEq. (17)eval-
v
uated at u+ = u(xi(t),t), ci is the absolute value of the |ux|< c(u ), (20)
speed of the i-th shock, and the last term in δ-function c
represents the fact that the value of u is evaluated right provided that the regions with u < u are initially at
c
in front of the wave. Equation (18) simply says that at rest. The condition in Eq. (20) simply means that the
the moment t the value of u jumps from u(xi(t),t) at wave with speed c will never catch up with the point at
x = xi(t) to the new value u−(u(xi(t),t)) given by Eq. which u=uc when u grows according to Eq. (6).
(16).
Having now defined the evolution of the slow variable
u,wecanwritedowntheequationofmotionoftheshocks V. SPATIOTEMPORAL PATTERNS
in terms of it
Let us now apply the procedure developed in the pre-
x˙ =s c , (19)
i i i ceding section to spatiotemporalpatterns forming in the
where c is the function of u =u x (t),t (Fig. 3). system under consideration. As in any excitable system
i + i
Equations (18)and (19) arethe basic equations of the [1–6], the basic pattern of this kind is a periodic array
free boundary problem describing(cid:0)the dy(cid:1)namics of the of shocks (wave train) traveling with constant velocity.
system for v 1 and sufficiently slowly varying initial From the results of the preceding section one can see
conditions, as≪suming that no creep occurs in the system that a periodic wave train should consist of sharp fronts
during its evolution. These equations have to be supple- inwhichthe valueofujumps fromu+ tou− followedby
mented with the way of dealing with creation and anni- refractoryregionswhereuslowlyrecoversbacktou+,so
hilation of shocks. From the physical considerations it theresultingpatternissaw-tooth(Fig. 5). Intherefrac-
is clear that a pair of shocks moving toward each other tory region u obeys Eq. (6), which in the frame moving
will annihilate. Similarly, when the value of u reaches 1 togetherwiththewavewithspeedcbecomesuz = v/c,
−
[recallthatourdistributionsofuvarysufficientlyslowly, where z = x ct is the self-similar variable. Therefore,
−
so one can ignore the variation of u in Eq. (9)], a slip thesolutionforuinthe refractoryregionwhichproperly
eventisinitiatedatthepointwhereu=1,creatingapair matches with the back of one shock at z = 0 and the
of counter-propagatingshocks. These features should be front of another at z = L, where L is the period of the
incorporatedinto the free boundary problemandwill be wave train, has the form−u = L−1(z+L)u− L−1zu+,
−
discussed in the following section. Notice that according where
toEq. (9)thedistributionofuforblocksatrestmustbe
cA
a continuously differentiable function in the continuum L= , (21)
v
limit, so the speeds c of the shocks given by Eq. (19)
i
will also be continuously differentiable functions of time. c=c(u )(Fig. 3),andAisgivenbyEq. (17). Equation
+
In writing Eq. (18) we assumed that we always have (21)should in fact be obviousfrom the purely geometric
u x (t),t > u , so that the shock waves are always the considerations. Thus, there exists a family of nonlinear
i c
waves of switching between the blocks with u = u at periodictravelingwavesolutionswhosespeedandperiod
+
(cid:0) (cid:1)
rest in front of the wave to u = u− at rest behind the dependstronglyonamplitude. Similarly,theperiodT of
wave. In other words, we do not consider the possibility thissolutionintime depends onthe amplitude simplyas
of the onset of creep behind the wave. The analysis of T = A/v, so the frequency of the shocks at a particular
creep motion is beyond the scope of the present paper. point as the wave train passes is ω A−1. The latter in
∼
6
u
0.2
0.1
2
v
’a/
x
0
-0.1
0 0.5 1 1.5 2
FIG.6. A source of counter-propagating waves. a/v2
FIG.7. Dependencea′(a) from Eq. (26) at γ =0.3.
fact represents a simple form of the Gutenberg-Richter
scaling law with scaling exponent b = 1 [29]. In other tice that Eq. (24) in general shows how to treat the sin-
words, this exponent would be observed in the system gularity in the free boundary problem (Sec. IV) associ-
under consideration if its dynamics were dominated by atedwiththe creationofa apair ofcounter-propagating
periodic wave trains. It is interesting to note that the shocks.
value of b actually observedfrom the earthquakes gener- The solution for xi(t) in Eq. (24) then allows to cal-
′
ally lies in the range 0.8<b<1.2 [30]. culate the distribution u after the shocks have passed
The above arguments are justified as long as u >
uc, so that no creep develops behind the traveling fr+ont. u′(x,t)≃u−(1)−a′(x−x0)2+vt, (25)
Therefore, according to Eqs. (17), and (21), traveling
′
where a is given by (see Appendix C)
wave trains exist only when L>L , where
min
Lmin =v−1c(uc)uc. (22) a′ = aκ+ 2av(1+κ) , (26)
− v+√v2+8a
Let us now consider another kind of spatiotemporal
patternsthatistypicalofexcitablesystems—thesource u−(1)= exp( πγ/ 1 γ2) [see Eq. (16)], and
− − −
ofcounter-propagatingwaves[1–6]. Theexistenceofsuch
p
asolutionisassociatedwiththefactthatifthereisanin- πγ πγ
κ= 1+ exp . (27)
hmoummogoefnueiosulsodcaistterdibauttxion=oxfu,athternesattssuocmhethpaotinthteinmtaimxie- (cid:18) (1−γ2)3/2(cid:19) − 1−γ2!
0
the value of u(x ) may reach 1, so that the blocks will p
0 From Eq. (27) one can see that the value of κ lies in the
become unstable creating a pair of counter-propagating ′
range 0 < κ < 1. The plot of a as a function of a at a
shocks(Fig. 6). After theseshocksmovedawayfromx ,
0 particular value of γ is presented in Fig. 7.
the distribution of u can once again have a maximum at
From Eq. (26) one can see that there is a maximum
x = x0, so after some time the value of u will increase of u′ at x = 0 (a′ > 0) when the shocks have passed, if
until it reaches 1 and the cycle repeats.
a<a , where
max
These solutions are indeed realized in our system.
Close to x = x0 we can approximate the distribution v2(1+κ)
of u as amax = 2κ2 . (28)
u(x,t) 1 a(x x )2+vt, (23) If this condition is satisfied, after time T =
0
≃ − −
v−1 1+exp πγ/ 1 γ2) the distributionofuwill
[see Eq. (6)], where we alsoassumedthatthe slipoccurs − −
at t = 0. This time-dependent distribution of u gives look(cid:16)exactly(cid:0)like thpe one in (cid:1)E(cid:17)q. (23) which we started
the values of u in front of the shocks at t > 0, after the with,butwithadifferentvalueofa. Therefore,Eq. (26)
slipoccurred,thusdeterminingthevelocityoftheshocks. defines an iterative map for the value of a correspond-
The analysis of the free boundary problem in this case ing to the solutions for u at times nT, where n is the
(Appendix C) then shows that right after the slip the number of iterations. Physically, a source with period T
positions of the shocks will be given by will form at x = 0 creating traveling waves propagating
away from it. Note, however, that even though the pe-
v+√v2+8a riod of the source is a constant that does not depend on
xi x0 √bt, b= , (24) the initial conditions, the solution represents an aperi-
≃ ± 2a
odicprocesssinceaftereachcyclethevalueofachanges.
where the plus and the minus correspond to the shocks The latter is represented in Fig. 7. Observe that we al-
propagatingtotherightandtotheleft,respectively. No- ways have a′ < a. From Fig. 7 one can also see that
7
after many periods the value of a will tend to zero. The want to calculate the dependence of the wave velocity
analysis of Eq. (26) (see also Fig. 7) shows that zero is on the amount of the accumulatedstress u by formally
+
′
the only fixed point of the map a a for any γ. For combining Eq. (12)at z =0 with the continuum version
→
smallenoughvaluesofaEq. (26)canbeapproximatedas of Eq. (9) at z = 0. This would give an incorrect result
a′ a[1 2a(1+κ)/v2]. Itistheneasytoshowthatafter c = 1/√1 u . This means that although the contin-
+
n ≃1 it−erationswe willhavea v2/[2n(1+κ)],where uum appro−ximation may be valid for finding the wave
n
≫ ≃
a is the result of n iterations of the map. Note that profile behind the slip point, one still needs to solve the
n
because of the fact that a becomes smaller and smaller discrete problem at the tip of the wave in order to find
witheachcycle,thedistributionofubecomesflatterand its velocity.
flatter [see Eq. (23)], what means that after a long time In the limit of vanishing loader plate velocity the only
the regions in the neighborhood of the source will tend parameter in the model is the dimensionless coefficient
to synchronize. The characteristic size R of a region in of viscous friction γ. All our analysis was performed for
which the synchronization takes place can be estimated arbitrary values of γ, so it is also applicable to the case
as R a−1/2, so we will have R t1/2. Also note that γ =0,i.e.,whentheviscousfrictionisabsent. Thiscase,
n
forth∼esamereasoniftherearetw∼osourceswithunequal however, has one special feature. The thing is that for
values of a < amax some distance R apart, after time γ = 0 we have uc =0 [see Eq. (14)], so as u approaches
t R2theonewiththesmallervalueofawilloverwhelm zero the speed of the wave approaches 1 and the width
th∼e one with the larger value of a. It is then natural to w goes to zero [see Eq. (15)]. Therefore, for a fixed
expectthatiftheinitialdistributionofuhasallmaxima value of ∆x 1 the continuum limit will no longer be
≪
with a < a , the system will eventually synchronize justifiedforthedescriptionofthewaveswhenubecomes
max
into uniform oscillations in which all the blocks slide at sufficientlyclosetozero. Also,whenγ =0theamplitude
once. A = 2u+ of these waves [see Eq. (17)] will be able to
Different situation is realized when a > a . In this become arbitrarily small.
max
casethedistributionofuformingaftertheshocksmoved According to the results of our analysis, the solution
away from the origin will be a minimum rather than a in the form of the traveling shock wave, as well as a pe-
maximum of u, so the next creation of a pair of shocks riodic wavetrainofa givenperiod, is unique in the limit
will not occur at x = x , but rather at different points of vanishingly small loader plate velocity, and its speed
0
in space. Therefore, it is natural to expect that for suf- is uniquely determined by the value of u in the front of
ficiently random initial conditions one would see shocks the wave. This is in contrast with the results of Langer
created and annihilated at random points in space, cre- and Tang [15] and Myers and Langer [16] for the models
ating a kind of spatiotemporal chaos. with velocity-weakening friction where they find a con-
tinuousfamily ofsuchsolutionsandthereforea selection
problem.
VI. CONCLUSION Notethatthesystemswithvelocity-weakeningfriction
transform to the model considered by us in the limit of
small V , where V is the characteristic velocity of vari-
Let us now summarize the results of our analysis and f f
ation of the friction force. Indeed, if the characteristic
discusstherelationshipofourresultswiththoseobtained
speed of a block during the slip exceeds V , the friction
for other models of stick-slip motion. In the present f
force will drop very rapidly, acting in the same way as
paper we analytically investigated the Burridge-Knopoff
in our model. To obtain the criterion of smallness of
slider-block model with Coulomb friction law, which is
V , we note that the effect of the velocity dependence of
different from the one used in most previous studies f
the frictionis mostpronouncedinthe slipregion,where,
[9,10,12–16]. This difference is expressedinthe fact that
according to Eq. (10), u ∆x= k /k , where we as-
ourfrictionisadiscontinuousfunctionwhenthevelocity t p c
∼
sumed c 1, so u is not near the threshold value u=1.
of the block becomes zero. As a result, the system is ca-
∼ p
In the original variables dX /dt (f f )/√k m, so
pable of propagating ultrasonic traveling solitary waves i r s c
∼ −
the condition V dX /dt gives
in the limit of zero loader plate velocity for any value f i
≪
of accumulated stress above the excitability threshold
f f
u = 0. The existence of such solutions is a character- V r− s. (29)
f
≪ √k m
istic feature of excitable systems [1–6]. c
In the continuum limit the solution in the form of the
This formula shouldin factbe obviousfromthe physical
travelingshockwaveisindependentoftheshort-scalebe-
considerations. Similarly, one should recover our results
havior of the system, except for the precise form of the
for the creep-slip models considered in [7] if the charac-
dependence of the wave’s velocity on the amount of the
teristic velocity of variation of the friction force is small
accumulated stress in front of the wave. The latter in
enough.
fact only weakly depends on the dynamics of the model
In contrast, if the condition opposite to the one in
atshortlengthscalesif the value ofu issufficiently close
Eq. (29) holds, no traveling shock waves can be real-
to the slipping threshold u = 1. Note that one might
ized in the continuum limit. This can be seen from the
8
following argument. In a traveling wave the profile of u f f =pa2(µ µ ). Onecanseethattheconditionof
r s r s
should be described by Eq. (12) in which, in the case of Eq−. (29) with V −=10−5 m/s (see [21])is satisfiedif p=
f
the velocity-weakeningfrictionone shoulddropthe term 2 106 Pa=20 atm, whatgives f f =1.2 10−7 N.
r s
× − ×
2cγu and replace u by u 1 close to the slip point. For this value of p the characteristic sliding velocity will
z
−There we must have u < 0−, so, according to the mod- be X˙ (f f )/ k m = 1.6 cm/s. The displacement
zz r s p
∼ −
ified Eq. (12), we will have c < 1 in the traveling wave. ∆X in a single wave can also be estimated as ∆X =
Onthe otherhand, the propagationofsubsonicwavesin 2(f f )/k = 1p10−8 m. The latter quantity turns
r s p
− ×
the system is impossible since from the dispersion rela- outtobequitesmall. Notethatboth∆X andX˙ increase
tion ω2 =1+k2 for the waves in the absence of friction aspincreases. Also,smallerpressurespwouldbeneeded
wegetthattheirphasevelocityc=ω/k=√1+k−2 >1 for smaller values of h. From [21] one can calculate α=
forallk. Ontheotherhand,numericalsimulationsshow 6.7 10−6 kg/s for this value of p, leading to γ =0.45.
that in the discrete models with velocity-weakening fric- In×thepresentpaperwestudiedperfectlyhomogeneous
tiontravelingshockwavesareindeedobserved. However, systemwithoutanyexternalnoise. Itisclearthatinreal
for ∆x 1 these waves move with speeds very close to systems some degree of noise should be present. It is
≪
the speed of sound and their width is just a few lattice interesting to know what effect the introduction of ran-
spacings (unless one is close to the threshold u = 1, see domnesswillhaveonthebehaviorofthetravelingwaves.
also[15,16]). Oneshouldthereforeexpectthatinthecon- Here one can distinguish two different situations. If one
tinuumlimitonlydiscontinuousshockstravelingwiththe adds some small randomness in the parameters of in-
speed of sound are feasible, in contrast to the situation dividual blocks, such as the friction coefficients, spring
studiedinthepresentpaper. Thisisacrucialdistinction constants and so forth, one would expect the renormal-
of the Burridge-Knopoff model with velocity-weakening ization of the speed of the traveling waves. The other
friction and the one with Coulomb friction law (Fig. 2) kind of noise would result in sudden transitions of the
studied by us. blocks at rest to sliding. This kind of noise can be read-
Aswasalreadydiscussedabove,thepropagationspeed ilyincorporatedintoourfreeboundaryproblembyintro-
of the wave depends on the short-scale dynamics of in- ducing generation of pairs of counter-propagating waves
dividual blocks. Therefore, in order to be able to per- at random points in space. Here the open question is
form numerical simulations of the models under consid- how the frequency of this generation should depend on
erationinthe continuumlimit, onehasto useverysmall the distance to the slipping threshold. It is also inter-
discretization of time. This should make direct simula- esting how the presence of such a noise would affect the
tions of our model rather difficult. On the other hand, wave statistics.
we showed that if the initial conditions vary sufficiently
slowly,thedynamicsofthesystemreducestothemotion
ofindividualshockwaves(seeSec. IV).Therefore,byre- APPENDIX A:
formulating the dynamics of the problem in terms of the
freeboundaryproblemonecandramaticallyincreasethe
Since the variation of u leading to the slip event is of
i
speed of the simulations, thus being able to observe the order ∆x2 and is small for small ∆x, for the blocks near
dynamics for much longer times. This should be impor-
theslippointonecanneglectboththevariationofu and
i
tant in the statistical analysis of our system, especially
the term 2γdu /dt in Eq. (7), so in the limit ∆x
i
in the context of earthquakes (see also [12–14]). − →
0 the right-hand side of this equation simply becomes
Let us see how the results obtained from the studies
u 1 = const. For convenience, let us introduce the
+
ofour modelshould translateto realstick-slipproblems. −
new variables τ and v , such that
i
Consider, for example, a Bristol board of thickness h =
2 mm lying on the surface made of the same material τ2
underexternalpressurepandpulledatthetop[21]. The t=ts+τ∆x, ui =u+ u+(∆x)2 vi+ . (A1)
− 2
coefficientsµ andµ ofthestaticandtheslidingfriction (cid:18) (cid:19)
r s
for these surfaces were found to be µr = 0.37 and µs = Since at τ =0 the s-th block just slipped, we must have
0.31, respectively [21]. The memory length which we [see Eq. (10)]
should use for the distance a between the blocks was
found to be a = 1 µm. A single block in the Burridge- 1 u+
Knopoff model should be associated with a piece of the vs−1(0)= −u , (A2)
+
board of dimensions a a h, so the mass m of the
block should be m = ×2.4 ×10−12 kg, where we used with the initial conditions
the density ρ = 1.2 g/cm3×. From the speed of sound
′
v (0)=0, v (0)=0, (A3)
c = a k /m [Eq. (3)] we can calculate the values of s s
c
k and k = k (πa/2h)2 (see the end of Sec. II). Using
c pp c where the prime denotes the derivative with respect to
c = 4000 m/s, we find k = 3.8 107 N/m, and k =
c × p τ.
2.3 101N/m. Thefrictionforcejumppersingleblockis
× In the limit ∆x 0 Eq. (7) can be written as
→
9
d2v 3 5τ τ2
dτ2i =vi+1+vi−1−2vi, i<s, (A4) vs(1−)1 = 8c4 + 6c3 + 2c2. (A12)
d2v τ2
dτ2i =vi−1−2vi− 2 , i=s. (A5) The knowledge of vs(1−)1 and vs(1) then allows to calculate
the next order correction v(2) by substituting them into
s
Theseequationsessentiallydescribealineararrayofunit
Eq. (A5)
massesconnectedbyspringswithspringconstant1,with
the s-thmassattachedto arigidwallandacteduponby 3τ2 5τ3 τ5 τ6
a time-dependent force −τ2/2. vs(2) = 16c4 + 36c3 − 60c + 360. (A13)
If the wave moves with speed c, after time τ = 1/c,
whichcorrespondsto∆t=∆x/c[seeEq. (A1)],itshould Substituting the obtained solution for v into Eq. (A2)
s
move the distance ∆x to the right. This means that at and using Eq. (A6) with i=s 1, we get
−
timets+1 =ts+∆twemusthaveui−1(ts)=ui(ts+∆t),
u′i−1(ts)=u′i(ts+∆t), or in terms of the new variables 1 + 3 + 5 = 1−u+. (A14)
2c2 8c4 16c6 u
+
1
vi−1(0)=vi(c−1)+ 2c2, (A6) Note that this procedure can be continued to an arbi-
1 traryorderinc−1,althoughthecalculationthenbecomes
vi′−1(0)=vi′(c−1)+ c. (A7) rathercumbersome. Letusonlyquotetheresultthatthe
correction to the left-hand side of Eq. (A14) due to the
NotethattheseconditionsintroduceafeedbackintoEqs. next order terms will be 0.2719/c8.
(A4) and (A5). Also, the solution of Eqs. (A4) should Finally,combiningEq. (A14)withEq. (A9),wearrive
be matched with the continuum profile behind the tip, at the interpolation formula given by Eq. (11).
which satisfies Eq. (12). We should therefore have
k2 APPENDIX B:
vs−k ∼= 2(c2 1), k ≫1. (A8)
−
HerewesolveEq. (12)withtheinitialconditionsgiven
Equations (A3) – (A7) can be solved numerically by
by Eq. (13) for arbitrary γ. Let us introduce the new
an iterative procedure. Then, calculating the value of
variable ξ = z/√c2 1 and a constant γ˜ = cγ/√c2 1,
vs(c−1),usingEq. (A6)asafunctionofc,onecanrelate where c is a function−of u (Fig. 3). Then, Eq. (12)−can
it to u throughEq. (A2). The results of this numerical +
+ be rewritten as
solutionforc>1arepresentedinFig. (3). Theseresults
also show that close to u =0 we have
+ u 2γ˜u +u=0, (B1)
ξξ ξ
−
1
c 1+ u2. (A9) with the same initial conditions written in terms of ξ:
≃ 2 +
u(0)=u , u (0)=0. (B2)
The problem in Eqs. (A3) can be treated analytically + ξ
by expanding its solution in the powers of c−1. Indeed,
Two situations are possible here. If we have γ˜ > 1,
if c is large,for τ <c−1 1 the right-handsides in Eqs.
≪ both roots of the characteristic equation of Eq. (B1) are
(A4) and (A5) can be neglected, so the solution of Eqs.
real, so the solution has the form
(A3) – (A7) for k 0 in the leading order in c−1 will be
≥
k2 kτ u=aeξ(γ˜+√γ˜2−1)+beξ(γ˜−√γ˜2−1). (B3)
v(0) = + . (A10)
s−k 2c2 c
Usingthe initialconditions fromEq. (B2),weobtainfor
Physically,this correspondstothesituationinwhichthe the coefficients a and b
springs between the blocks are absent.
Tocalculatethefirstordercorrectionvs(1),weusevs(0−)1 a= u+ 1 γ˜ , b= u+ 1+ γ˜ .
in Eq. (A5) to obtain 2 − γ˜2 1! 2 γ˜2 1!
− −
p p (B4)
τ2 τ3 τ4
v(1) = + . (A11)
s 4c2 6c − 24
Fromthesolutionweseethatbehindthewavethedistri-
bution of u asymptotically approaches zero. This means
This solution is then used to fix through Eqs. (A6) and
(A7) the initial condition in Eq. (A4) with i = s 1 to that we have u− =0, but the system comes to rest only
− at z = , so the width w of the wave is infinite. The
thenextorderinc−1. Usingthesolutionsvs(0),vs(0−)1,and form of−th∞e solution is shown in Fig. 4(a).
v(0) in Eq. (A4), one can then calculate v(1) :
s−2 s−1
10