Dynamics of stick-slip in peeling of an adhesive tape

0.47 MB·

by  Rumi De

Dynamics of stick-slip in peeling of an adhesive tape
Rumi De,1 Anil Maybhate,1 ∗ and G. Ananthakrishna1,2, †
1 Materials Research Centre, Indian Institute of Science, Bangalore-560012, India
2 Centre for Condensed Matter Theory, Indian Institute of Science, Bangalore-560012, India
We investigate the dynamics of peeling of an adhesive tape subjected to a constant pull speed.
5
We derive the equations of motion for the angular speed of the roller tape, the peel angle and the
0
pull force used in earlier investigations using a Lagrangian. Due to the constraint between the
0
pullforce, peel angle and thepeel force, it falls into thecategory of differential-algebraic equations
2
requiring an appropriate algorithm for its numerical solution. Using such a scheme, we show that

n stick-slip jumps emerge in a purely dynamical manner. Our detailed numerical study shows that
a these set of equations exhibit rich dynamics hitherto not reported. In particular, our analysis
J
shows that inertia has considerable influence on the nature of the dynamics. Following studies in

8 the Portevin-Le Chatelier effect, we suggest a phenomenological peel force function which includes
1 the influence of the pull speed. This reproduces the decreasing nature of the rupture force with
the pull speed observed in experiments. This rich dynamics is made transparent by using a set

] of approximations valid in different regimes of the parameter space. The approximate solutions
i
c capturemajorfeaturesoftheexactnumericalsolutionsandalsoproducereasonablyaccuratevalues
s for thevarious quantities of interest.
-
l
r
PACSnumbers: 05.45.Pq,62.20.Mk,68.35.Np
t
m
.
t I. INTRODUCTION sliding [9], and even earthquake dynamics is thought to
a
m result from stick-slip of tectonic plates [10]. Stick-slip
is characterized by the system spending most part of
-
d Peeling is a kind of fracture that has been stud- the time in the stuck state and a short time in the slip
n ied experimentally in the context of adhesion and is a state, and is usually seen in systems subjected to a con-
o technologically important subject. Experimental stud- stantresponsewheretheforcedevelopedinthesystemis
c ies on peeling of an adhesive tape mounted on a cylin- measured by dynamically coupling the system to a mea-
[
drical roll are usually in constant pull speed condition suring device. One common feature of such systems is

2 [1, 2, 3, 4, 5, 6]. More recently, constant load experi- thattheforceexhibits“negativeflowratecharacteristic”
v ments have also been reported [3, 7]. Early studies by (NFRC).Modelswhichattempttoexplainthe dynamics
0 Bikermann[5], Kaeble [6] have attempted to explain the of such systems use the macroscopic phenomenological
6
resultsbyconsideringthesystemasafullyelasticobject. NFRC feature as an input, although the unstable region
1
This is clearly inadequate as it ignores the viscoelastic isnotaccessible. Thisistrueformodelsdealingwiththe
3
0 nature of the glue at the contact surface and therefore dynamicsoftheadhesivetapeaswell. Tothebestofour
4 cannotcapturemanyimportantfeaturesofthedynamics. knowledge,thereisnomicroscopictheorywhichpredicts
0 ThefirstdetailedexperimentalstudyofMaugisandBar- the origin of the NFRC macroscopic law except in the
t/ quins [1] show stick-slip oscillations within a window of case of the PLC effect [11, 12] (see below).
a pull velocity with decreasing amplitude of the pull force
m As there is a considerable similarity between the peel-
as a function of the pull velocity. Further, these authors
ingofanadhesivetapeandthePLCeffect,itisusefulto
- reportthatthepullforceshowssinusoidal,sawtoothand
d consider the similarities in some detail. The PLC effect
highly irregular (chaotic as these authors refer to) wave
n refersto a type ofplastic instability observedwhen sam-
o patterns with increasing velocities. More recently, Gan-
ples of dilute alloys are deformed under constant cross
c dur et al. have carriedout a dynamical time series anal-
head speeds [13]. The effect manifests itself in the form
: ysis of the force waveforms, as well as those of acoustic
v of a series of serrations in a range of applied strain rates
i emission signals and report chaotic force waveforms at
X and temperatures. This feature is much like the peeling
the upper end of the pull velocities [3]. One characteris-
of an adhesive tape. Other features common to these
r ticfeatureofthepeelingprocessisthattheexperimental
a two situations are: abrupt onset of the large amplitude
strainenergyreleaserateshowstwostablebranchessep-
oscillationsatlowapplied velocitieswith a graduallyde-
aratedbyanunstablebranch. Stick-slipbehavioriscom-
creasingtrendandNFRC,whichinthePLCeffectrefers
monlyobservedinanumberofsystemssuchasjerkyflow
to the existence of negative strain rate sensitivity of the
or the Portevin-Le Chatelier (PLC) effect [8], frictional
flow stress. In the case of the PLC effect, the physical
originof the negative strainrate sensitivity is attributed
to the ageingofdislocationsandtheir tearingawayfrom
the cloud of solute atoms. Recently, the origin of the
∗Presentaddress: WeillMedicalCollegeofCornellUniversity,New
York,USA. negative SRS has been explicitly demonstrated as aris-
†Electronicmail: [email protected] ingfromcompetingtimescalesofpinningandunpinning
2
in the Ananthakrishna’s model [11, 12]. In the case of
adhesive tape, the origin of NFRC can be attributed to
ω
theviscoelasticbehaviorofthefluid. (Constantloadand
O l O’
constant load rate experiments are possible in the PLC
α
also.) While simple phenomenological models based on R θ
L
NFRCexplainthegenericfeaturesofthePLCeffect[14],
there appears to be some doubts if the equations of mo- P
tionconventionallyusedinthepresentcaseofpeelingare
adequate to describe the velocity jumps [2, 4]. Indeed,
FIG. 1: Schematic plot of experimental setup
theseequationsofmotionaresingularandposeproblems
in the numerical solutions.
Apart from detailed experimental investigation of the capturedtoareasonableaccuracyusingasetofapproxi-
peeling process, Maugis and Barquins [1], have also con- mationsvalidindifferentregimesoftheparameterspace.
tributed substantially to the understanding of the dy- Even though, our emphasis is on demonstrating the cor-
namics of the peeling process. However, the first dy- rectness of these equations of motion and richness of the
namical analysis is due to Hong and Yue [2] who use an inherent dynamics that capture the qualitative features
“N”shapedfunctiontomimicthedependenceofthepeel of the peeling process, we also attempt to make a com-
force on the rupture speed. They showed that the sys- parison of the experimental results mentioned above to
tem of equations exhibits periodic and chaotic stick-slip the extent possible.
oscillations. However, the jumps in the rupture speed
are introduced externally once the rupture velocity ex-
ceeds the limit of stability [4, 15]. Thus, the stick-slip II. EQUATIONS OF MOTION
oscillations are not obtained as a natural consequence of
the equations of motion. Therefore, in our opinion the For the sake of completeness, we start by considering
results presented in Ref. [2] are the artifacts of the nu- the geometry of the experimental setup shown schemat-
merical procedure followed. Ciccotti et al. [4] interpret ically in Fig. 1. An adhesive roll of radius R is mounted
the stick-slip jumps as catastrophes. Again, the belief on an axis passing through O normal to the paper and
that the jumps in the rupture velocity cannot be ob- is pulled at a constant velocity V by a motor positioned
tained from the equations of motion appears to be the at O′ with a force F acting along PO′. Let the distance
motivation for introducing the action of discrete opera- between O and O′ be l, and that between the contact
tors on the state of the system to interpret the stick-slip point P to O′ be L. The point P moves with a local
jumps [4], though they do not demonstrate the correct- velocityv whichcanundergorapidburstsinthe velocity
nessofsuchaframeworkforthesetofequations. Lastly, during rupture. The force required to peel the tape is
there arenoreportsthat explainthe decreasein the am- usually called the force of adhesion denoted by f. The
plitude of the peel force with increasing pull speed as two measured branches referred to earlier, are those of
observed in experiments. As there is a general consen- the function f in a steady state situation of constant
sus that these equations of motion correctly describe the pulling velocity (i.e., there are no accelerations). The
experimental system, a proper resolution of this question line L makes an angle θ with the tangent at the contact
(on the absence of dynamical jumps in these equations) point P. The point P subtends an angle α at O, with
assumes importance. the horizontal line OO′. We denote the elastic constant
The purpose of this paper is to show that the dynam- ofthe adhesivetape byk,the elasticdisplacementofthe
ics of stick-slip during peeling can be explained using a tape by u, the angular velocity by ω and the moment
differential-algebraicscheme meant for such singular sit- of inertia of the roll by I. The angular velocity itself is
uations [16] anddemonstrate the richdynamics inherent identified by ω = α˙ +v/R. The geometry of the setup
to these equations. In what follows we first derive the gives L cosθ = −l sinα and L sinθ = l cosα−R which
equations of motion (used earlier [2]) by introducing an furthergives,L2 =l2+R2−2lRcosα. Thetotalvelocity
appropriate Lagrangian for the system. Then, we use V atO′ is then made up ofthree contributions[1], given
an algorithm meant to solve differential-algebraic equa- by V =v+u˙ −L˙, which gives
tions [16] and present the results of our simulations for
various parameter values. One of our major findings is v =V +L˙ −u˙ =V −Rcosθ α˙ −u˙. (1)
that inertia has a strong influence on the dynamics. In
Following standard methods in mechanics, it is straight-
addition, following the dynamization scheme similar to
forward to derive the equations of motion for α and ω
the one used in the context of the PLC effect [14], we
byconsidering(α,α˙,u,u˙)asthegeneralizedco-ordinates.
suggestthatthe peelforcedepends onthe appliedveloc-
ThecorrespondingLagrangianofthesystemcanbewrit-
ity. Using this form of peel force leads to the decreasing
ten as
nature of the magnitude of the pull force as a function
of applied velocity. For certain values of the inertia, we I k
find canard type solutions. These numerical results are L(α,α˙,u,u˙)= [ω(α,α˙,u,u˙)]2− u2. (2)
2 2
3
We write the dissipation function as
(a)
R=Φ(v,V)= f(v,V)dv, (3)
Z 400
where f(v,V) physically represents the peel force which )
V
we assume is dependent on rupture speed as well as the v,
pullspeedassumedtobederivablefromapotentialfunc- (
f
tion Φ(v,V). The physical origin of this is due to the
competitionbetweentheinternalrelaxationtimescaleof 200
the viscoelastic fluid and the time scale determined by
the applied velocity [17]. When the applied velocity is 10−3 10−2 10−1 v 100 101
low, there is sufficient time for the viscoelastic fluid to
relax. As we increasethe appliedvelocity,the relaxation
(b)
of the fluid gets increasingly difficult and thus behaves
much like an elastic substance. The effect of competing 102
time scales is well represented by Deborah number [18]
)
which is the ratio of time scale for structural relaxation v
(
G
to the characteristic time scale for deformation. Indeed,
in the studies onHele-Shaw cellwith mud as the viscous 1
10
fluid, one observes a transition from viscous fingering to
viscoelastic fracturing [19] with increasing rate of inva- 10−6 10−4 10−2 v 100 102
sion of the displacing fluid.
As stated in the Introduction, the existing models do
FIG. 2: (a) Plots of f(v,V) as a function of v (x axis in log
not explainthe decreasing amplitude ofpull force. Simi-
scale) for V = 1 (solid curve), V = 2 (dashed curve), V = 4
larfeatureobservedinthePLCserrationshasbeenmod-
(dashed and dotted curve), V = 6 (dotted curve); see Eq.
eled using a scheme referred to as dynamization of the
(14). (b)Experimentalstrain energyrelease rate,G(v)curve
negative strain rate sensitivity (SRS) of the flow stress as in Ref. [1]. [Units of f(v,V) is in N, G(v) in J/m2 and v,
f(ǫ˙P) [14, 20], where ǫ˙p is the plastic strain rate. Based V are in m/s ]
on arguments similar to the preceding paragraph, they
modifythisfunctiontodependontheappliedstrainrate,
ǫ˙ , i.e., the negativeSRS ofthe flow stress is takento be we have used cosθ ≃ −sinα ∼ −α. While Eqs. (6)-
a
f(ǫ˙ ,ǫ˙ ) such that the gap between the maximum and (9) are differential equations, Eq. (10) is an algebraic
P a
the minimum of the function f(ǫ˙ ,ǫ˙ ) decreases with in- constraint necessitating the use of differential-algebraic
p a
creasing ǫ˙ . Following this, we consider f to depend on scheme to obtain the numerical solution [16].
a
V also,inawaythatthegapinf decreasesasafunction The fixed point of Eqs. (6), (7), (9), and (10) is given
of the pull speed V (Fig. 2). by α = 0,ω = V/R,v = V,F = f(V,V). (For numerical
Using the Lagrange equations of motion, solution, in the above equations we have actually used
sinα in place of α.) This point is stable for f′(V,V)>0
d ∂L ∂L ∂R and unstable for f′(V,V) < 0. As V is varied such that
− + = 0, (4)
dt(cid:18)∂α˙(cid:19) ∂α ∂α˙ the sign of f′(V,V) changes from negative to positive
d ∂L ∂L ∂R value, the system undergoes a Hopf bifurcation and a
− + = 0. (5) limit cycle appears. The limit cycles reflect the abrupt
dt(cid:18)∂u˙(cid:19) ∂u ∂u˙
jumps between the two positive slope branches of the
function f(v,V).
we obtain the same set of ordinary differential equations
as in Ref. [2] given by
α˙ = ω−v/R, (6) III. ALGORITHM
Iω˙ = FRcosθ =−FRsinα≃−FRα (7)
The singular nature of these equations becomes clear
F˙ = ku˙ =k(V −v)−k cosθ(ωR−v), (8) if one were to consider the differential form of Eq. (10)
≃ k[V −v+Rαα˙], (9) given by
1
with an algebraic constraint v˙ = F˙ (1−cosθ)+F sinθ θ˙ , (11)
f′(v,V)
h i
F(1−cosθ)−f(v,V)≃F(1+α)−f(v,V)=0. (10) ≃ [F˙(1+α)+Fα˙]/f′, (12)
(The last equation results from the elimination of two where the prime denotes the derivative with respect to
second order equations for α.) In Eqs. (7), (9), and (10) v. Equations (11) with Eq. (6), (7), and (8) [or (9)]
4
constitutethefullsetofevolutionequationsforthevector 300
(α,ω,F,v). However, it is clearly singular at points of (a) B C
extremumoff(v,V),requiringanappropriatenumerical
250
algorithm.
WenotethatEqs.(6),(7),(8),and(10)canbewritten F
as 200
MX˙ =φ(X), (13) A
D
150
whereX=(α,ω,F,v),φisavectorfunctionthatgoverns 10−5 10−3 v 10−1 101
the evolution of X and M is a singular “mass matrix”
[16] given by, 290 (b)
b
c
1 0 0 0 260
0 1 0 0
M = . F
0 0 1 0
230
0 0 0 0
  d
a
Equation (13) is a differential-algebraic equation (DAE) 200
and can be solved using the so called singular pertur- 10−3 10−2 10−1 v 100 101
bation technique [16] in which the singular matrix M is
perturbed by adding a small constant ǫ such that the 0.1
( c )
singularityis removed. The resulting equations canthen
0.05
be solved numerically and the limit solution obtained as
ǫ → 0. We have checked the numerical solutions for ǫ
values rangingfrom10−7 to 10−15 in somecasesandthe α 0
results do not depend on the value of ǫ used as long as
−0.05
it is small. The results presentedbelow, however,arefor
ǫ = 10−7. We have solved Eq. (13) using a standard
−0.1
variable-order solver, MATLAB ODE15S program. 4 t 4.1 4.2
We have parametrized the form of f(v,V) as
f(v,V) = 400v0.35+110v0.15+130e(v/11)−2V1.5 290 (d)
−(415−45V0.4−0.35V2.15)v0.5, (14)
270
to give values of the extremum of the peel velocity that F
mimic the general form of the experimental curves [1].
250
The measured strain energy release rate G(V) from sta-
tionary state measurements is shown in Fig. 2(b). The
230
decreasing nature of the gap between the maximum and
4 t 4.1 4.2
minimum of f(v,V) for increasing V is clear from Fig.
2(a). [The valuesoff(v,V)couldnotbecorrectlydeter- FIG. 3: (a) A typical phase space trajectory in the v −F
mined as G(V) is in J/m2 requiring more details. How- plane for V = 0.4,I = 10−5. The corresponding f(v,V) is
ever,thevalueofFmax isclosertoRef. [2]andthejumps shown by a solid curve. (b) A phase space trajectory in the
in v are similar to those in experiments.] The reason for v−F plane for V =1.0 and I =10−5. (c) A plot of α(t) for
usingtheformgivenbyEq. (14)isthattheeffectsofdy- V =1 and I =10−5. (d) Aplot of F(t)(period 4) for V =2
namizationareeasilyincludedthroughitsdependenceon andI =10−5. (Unitsof v, V are in m/s, F in N,I in kgm2
thepullingvelocitywhilemorecomplicatedtermsarere- and t in s.)
quiredtomimic completelythe experimentalcurve(par-
ticularly the flat portion). However, we stress that the
trend of the results remains unaffected when the actual We have found that transients for some regions of
experimental curve is used except for the magnitude of parameters space take considerable time to die out.
velocity jumps and the force values. The results reported here are obtained after these long
transients are omitted. These equations exhibit rich
dynamics, some even unanticipated. Here we report
IV. RESULTS typical results for two important parameters, namely
the pull velocity V (m/s) and the inertia I (kg m2),
We have studied the dynamics of the system of keeping the elastic constant of the tape k = 1000 N/m,
equations for a wide range of values of the parameters. R = 0.1 m and l = 1 m [2]. The influence of k will also
5
bementionedbriefly. (Henceforth,wedropthe unitsfor
300 (a) b c
the sake of brevity.) We find that the observed jumps f g
of the orbit in the v-F plane occur in a fully dynamical e d
way. Moreimportantly,we findallthe three possibilities
260
namely, the orbit can jump when it approaches the
limit of stability, before or beyond that permitted by F
f(v,V). The dynamics can be broadly classified into
i
low, intermediate and high regimes of inertia. 220 l m
k j
(i) Low inertia. Here also, there are three regimes: a n
180
low, intermediate, and high pull velocity.
10−3 10−2 10−1 v 100 101
(a) Consider keeping inertia I at a low value (say
0.3
I = 10−5) and V also at a low value (say, near the top, (b)

say V = 0.4 ). Here we observe regular saw tooth form
for the pull force F. The phase plot in the F-v plane 0.1
is as shown in Fig. 3(a). The corresponding function α a b k a
f(v,V) is also shown by the continuous curve. We see
−0.1
that the trajectory jumps almost instantaneously from
B to C on reaching the maximum of f(v,V) (or from D
to A when it reaches the minimum). The system spends −0.3
considerablymoretimeonAB comparedtothatonCD. 3.9 4.05 t 4.2 4.35
However,this feature ofjumping ofthe trajectoryatthe
limit of stability is only true for small values of I and 300 (c) b,c
f,g
when V is near the limit of stability. At slightly higher
pull velocity, say V = 1, even for small I, say I = 10−5, d,e
the jumps occur even before reaching the top or bottom 260
( the points B and D) as can be seen from Fig. 3(b) for F
V = 1. The small amplitude high frequency oscillations i
seen in the phase plots [ Fig. 3(a), and 3(b)] on the 220 l,m
branchAB are due to the inertial effect, i.e., finite value
of I. These oscillations are better seen on the α(t) plot
j,k
shown in Fig. 3(c). For these values of parameters, the 180 a n,a
system is aperiodic. 3.9 4.05 t 4.2 4.35
16
b) As we increase V, even as the saw tooth form of
(d)
F is retained, various types of periodic orbits [period 4
12
shown in Fig. 3(d) for V = 2] as well as irregular orbits
are seen. In both cases (periodic as well as chaotic) the v
8
trajectory jumps from high velocity branch (CD) to the
low velocity branch before traversing the entire branch
4
or sometimes going beyond the values permitted by f.
The value of F at which the orbit jumps is different
for different cycles. For I = 10−5, at high velocity, say 0
4 4.1 t 4.2 4.3 4.4
V =4, the phase plot is periodic.
FIG. 4: (a) Phase space trajectory in the v −F plane for
(ii) Intermediate and high inertia. a single cycle for I = 10−2 and V = 1. The corresponding
(a) As the results of small V for intermediate and high f(v,V) is shown by a thick solid curve. (b) Corresponding
inertia are similar, we illustrate the results for I =10−2 plots of α(t), (c) the pull force F(t) (period 8) and (d) the
and V = 1. The v − F phase plot, α,F and v are peel velocity v(t). (Units of v, V are in m/s, F in N, I in kg
shown in Fig. 4(a)-4(d). Consider, Fig. 4(a) showing m2 and t in s.)
a typical phase space trajectory for a single cycle. The
corresponding function f(v,V) is also shown by the
thick continuous curve. We see that the maximum (and than f (f ).] When the trajectory jumps from
max min
minimum) value of F is larger (or smaller) than that AB to CD at the highest value of F for the cycle, the
allowed by f(v,V). [This feature holds when the inertia trajectory stays on CD for a significantly shorter time
is in the intermediate regime also, though the values of comparedtothesmallinertiacase(I =10−5)andjumps
maxima (minima) of F are not significantly larger (less) back to AB well before F has reached the minimum of
6
f(v,V), i.e., ∆F is much smallerthan fmax−fmin. The 450
pull force F cascades down through a series of back and (a)
forth jumps between the two branches till the lowest
value of F for the cycle is reached. Note that F at the
350
point n is less than f . For the sake of clarity, two
min
different portions of the trajectory are marked abcdefg F
and ijklmna corresponding to the top and bottom
regions of the plot. The corresponding points are also 250
identified on the F(t) plot. After reaching n, the orbit
jumps to a on AB, the trajectory decides to move up
all the way till F reaches a maximum value (larger than 150
fmax, the point b) without jumping to the CD branch. 4 t 4.2 4.4
This part of F as a function of time, which is nearly
linearonAB, (i.e., the segmentab)displaysanoticeable 0.6
(b)
sinusoidal modulation. The sinusoidal form is better
seen in α [Fig. 4(b)]. Note that the successive drops in 0.6 k
F are of increasing magnitude. The jumps between the 0.3 a
two branches in the v−F plane are seen as bursts of v α
4.21 4.24
[Fig. 4(d)]. For these values of parameters, the system
0.2
is periodic.

(b)AsweincreaseV,thesinusoidalnatureofF andα
−0.2
becomesmoreclearwithitsrangebecominglargerreach-
ing a nearly sinusoidal at V = 4 for large I. [The range
ab in Fig. 4(c) expands. Compare Fig. 5(a).] The mag- 4 t 4.2 4.4
nitude of ∆F on the CD branch for small V and mod-
eratelyorlargeI,graduallydecreaseswithincreasingV. 450
The magnitude of ∆F itself decreases as I is increased. (c)
In the limit of large V and I, the drops in F and α be- b c
comequitesmallwhicharenowlocatednearthemaxima
350
and minima of these curves. This is shown in Figs. 5(a)
and 5(b). The sinusoidal nature is now obvious even in
F
F(t) unlike for smaller V and I where it is clear only in
α(t) for the low v branch. Note that for V = 4, the na-
250
tureoff(v,4)isnearlyflat. Thisinducescertainchanges
a
in the v−F phase plot that are not apparent in F and k
α. Thejumpsbetweenthetwobranchesarenowconcen-
trated in a dense band at low and high values of F. In 150 d
this case, the maximum (minimum) value of F is signifi- 10−1 v 100 101
cantlylarger(less)thanf (f ). Theserapidjumps
max min
between the branches manifest as jitter at the top and 10
bottom of F and α. (d)
8
Unlike for small V [Fig. 4(a)], the nature of the tra-

jectory in Fig. 5(c) is different. After reaching a critical
6
value of F near the maximum value of F (the point b),
v
the orbit spirals upwards and then descends down till
4
another criticalvalue ofF (the pointc) is reached. Hav-

ing reached c, the orbit monotonically comes down till d
2
where it jumps to the AB branch. Beyond this point, it
again spirals upwards till the point a is reached. There-
0
after,F monotonicallyincreasestillbisreached. There-
4 t 4.2 4.4
gionsabandcdaretheregionswhereF showsanearsinu-
soidalform. Theregionsbcandda arethe regionswhere
FIG. 5: (a) A plot of F(t) for V = 4 and I = 10−2. (b)
the orbit jumps between the branches rapidly. These
Corresponding plot of α(t). The inset shows an expanded
manifest themselves as bursts of v which tend to bunch
plot of decreasing trend of α(t). (c) Corresponding plots of
together almost into a band. [Compare Fig. 4(d) with
phasespacetrajectorythatreflectsthechaoticnatureand(d)
Fig. 5(d).] It is interesting to note that the jumps be- the peel velocity v(t). (v, V are in m/s, F in N, I in kg m2
tween the two branches occur exactly at points where and t in s.)
7
df/dv = 0, even when the maximum (minimum) of F
120
are higher (lower) than that allowed by the stationary
curve f(v,V). The variables are aperiodic for the set
of parameters. The phase plots appear to be generated
80
by an effective f(v,V) that is being cycled. [This visual
feeling is mainly due to the fact that jumps between the _F
_∆
branches still occur at the maximum and minimum of
40
f(v,V).]
The influence of k is generally to increase the range
of the pull force F as can be easily anticipated and to
0
decrease the associated time scale.
0 1 2 3 4 5
It may be desirable to comment on the similarity of pull velocity V
thenatureoftheforcewaveformsdisplayedbythemodel
equations with those seen in experiments. As mentioned FIG.6: Theplotshowsthemeanforcedrop∆F asafunction
intheintroduction,apartfromqualitativestatementson of the pull speed V, for two distinct values of I. The dashed
thewaveformsinRef. [1](suchasperiodic,sawtoothetc., linecorrespondstoI =10−2whilethedottedlinecorresponds
toI =10−5. (v,V areinm/s,F inN,I inkgm2 andtins.)
whichareseeninthemodelaswell),itshouldbestressed
thatthereisapaucityofquantitativecharacterizationof
the waveforms. In this respect, the study by Gandur
et al. [3] fills the gap to some extent. These authors 300
have carried out a dynamical analysis of the time series
for various values of the pull velocities (for a fixed value
of the inertia corresponding to their experimental roller
tapegeometry). Inordertocomparethisresult,wehave
F
calculated the largest Lyapunov exponent for a range of
values of I and V. The region of chaos is in the domain
200
ofsmallpullvelocitiesV whenI issmall. Themaximum
Lyapunovexponentturnsouttoberatherhigh,typically
around7.5bits/sincontrasttothesmallvaluesreported
in Ref. [2]. The large magnitude of the positive expo- 1e-05 0.001 0.1 10
nent in our case can be traced to the large changes in
v
the Jacobian, as df(v,V)/dv varies over several order of
magnitude(∼ 106) as a function of the peeling velocity FIG.7: Aphaseplotofcanardtypeofsolutioninv−F plane
and hence as a function of time. In contrast, Hong et for V =0.4 and I =10−3. (v, V are in m/s, F in N, I in kg
al. use anN shapedcurvewheredf(v,V)/dv isconstant m2 and t in s.)
(and small) on both low and high V branches. However,
these large values of Lyapunov exponents are consistent
with rather high values obtained by Gandur et al. [3] comparisonwithexperimentsrequiresanappropriatecut
from time series analysis of the pull force. We also find in the I−V plane consistent with the experimental val-
chaos for intermediate and high inertia in the region of ues of I and V even where they are given. However, as
high velocities where the value of the Lyapunov expo- the values of I are not provided, full mapping of chaotic
nent is small, typically 0.5. The small value here again solutions is not possible. (We also note that Gandur et
canbetracedtothesmallchangesindf(v,V)/dv athigh al. [3] use a different tape from that used in Ref. [1],
velocities. as is clear from the instability range, leading additional
It must be mentioned that comparison with experi- difficulties in comparison.)
mentsisfurthercomplicatedduetothepresenceofatwo Onequantitativeresultthatcanbecomparedwithex-
parameterfamilyofsolutionsstronglydependentonboth periment is the decreasing trend of the force drop mag-
I and V. Thus, the phase diagram is complicated, i.e., nitude. We have calculated the magnitude of the force
the sequence of solutions encountered in the I-V plane dropsduringstick-slipphaseasafunctionthepullveloc-
as we change V or I or both does not in general display ity V for both low ( I = 10−5) and high (I = 10−2) in-
any specific ordering of periodic and chaotic trajectories ertiacases. Figure6 showsthe monotonicallydecreasing
(see Fig. 1 of Ref. [21]) usually found in the well known trend of average∆F(t) as V is increased, for both small
routes to chaos. (For instance 2n periods should be ob- andlargeI,afeatureobservedinexperiments[1]. These
servedbeforetheoddperiods[22].) Indeed,inourmodel, two distinct behaviors are a result of the dynamization
we findthe oddperiods 3,5,7etc, onincreasingV (orI), of f(v,V) as in Eq. (14).
without seeing all the 2n periods. (These odd periods Finally, as another illustration of the richness of the
also imply chaos at parameter values prior to that cor- dynamics seen in our numerical simulations, we show in
responding to these periods.) In view of this, a correct Fig. 7, a plot of an orbit that sticks to unstable part
8
of the manifold before jumping back to the AB branch. extentofvaluesofF(t)rangebetween185and450much
Suchsolutionsareknownascanards[23]. Thoughcanard beyond f(v,4) whose range is around 300. This feature
type of solutions are rare, we have observed them for is less dominant for small I and small V case.
high values of I and low values of V. In our case, such
solutions are due to the competition of time scale due
to inertia and that due to v. This again illustrates the V. APPROXIMATE ANALYSIS OF THE
influenceofinertiaoftherollonthedynamicsofpeeling. DYNAMICS
It is clear that these equations exhibit rich and com-
plex dynamics. A few of these features are easily un- As the dynamics is described by a coupled set of dif-
derstandable, but others are not. For instance, the saw- ferential equations with an algebraic constraint, the re-
tooth form of F for low inertia and low pull velocity can sultsarenottransparent. Wefirstattempttogetinsight
be explained as resulting from the trajectory sticking to intothe complexdynamics throughsomesimple approx-
stable part of f(v,V) and jumping only when it reaches imations valid in each of the regimes of the parameters.
the limit of stability. For these parameter values, as the Solution of these approximate equations will require ap-
time spent by the system is negligible during the jumps propriate initial values for the relevant variables which
between the branches AB and CD (and vice versa), the willbeprovidedfromtheexactnumericalsolutions. Due
system spends most of the time on the branch AB and tothenatureofapproximations,theresultsareexpected
much less on CD due to its steep nature. Then, from tocaptureonlythe trendandorderofmagnitudesofthe
Eq. (9), it is clear that we should find a sawtooth form effects that are being calculated. But as we will show,
whenever the peel velocity v jumps across the branch to even the numbers obtained match quite closely with the
a value of v larger than the pull velocity V. exact numerical results.
However,severalfeatures exhibited by these systemof Our idea is to capture the dynamics through a single
equations are much too complicated to understand. We equation(asfaraspossibleoratmosttwoasinthehigh
first list the issues that need to be explained. I andV case)byincludingalltherelevanttimescalesand
solve the relevant equation on each branch. For this we
(I) Small I. note that the equation for α and v play a crucial role as
(a) We find high frequency tiny oscillations superposed the inertial contribution appears only through Eqs. (6)
on the linearly increasing F [on the AB branch or and (7) and the time spent by the system is controlled
better seen in the α plot Fig. 3(c)]. This needs to be by the equation for v, Eq. (12). Using Eqs. (6) and (7),
understood. we get
F(t)Rα
(b) The numericalsolutionsshow thatthe influence of α¨ =− −v˙/R. (15)
I
inertia can be important even for small I and small V.
For instance, the jumps between AB and CD branches Thegeneralequationforαcanbewrittendownbyusing
occur even before F reaches the extremum values of f. Eq. (12), in Eq. (15), we get
(II)Forintermediateandhighvaluesofinertia,forlow
V case. FRα [F˙(1+α)+Fα˙]
α¨ = − − , (16)
(a) We observe several relatively small amplitude saw I Rf′
tooth form of F on the descending part of the pull force
FRα [F˙ +Fα˙]
F. Theseappearasasequenceofjumpsbetweenthetwo ≃ − − . (17)
branchesinthev−F planewhichweshallrefertoasthe I Rf′
“jumpingmode”. Aproperestimate ofthemagnitudeof
In obtaining Eq. (17), we have used 1+α ≃ 1 which
∆F is desirable.
is valid except for high I and high V. Further, in
(b) In addition, there appears to be a critical value
most cases, we can drop Rαα˙ as the magnitude of
of F for a given cycle below which the return jumps
this term is small and use F˙ ≃ k(V − v). To be
from AB to CD stop and one observes a monotonically
consistent we use F(t) ≃ F + k(V − v)t. We note
increasing trend in F [ab in Fig. 4(c)]. in
however that even for high I and high V where α is
notsmall,dropping1+αandRαα˙ causesonly10%error.
(III) High I and high V.
(a) The jumps between the branches occur at a very
Case I, small I
high frequency [Fig. 5(c)] and now are located near the
Onthe lowvelocitybranchAB, asv/Ris smallinEq.
extremum values of F and α. But these regions are sep-
(6), we can drop v˙ term in Eq. 15. Thus,
arated by a stretch where the orbit monotonically in-
creases on the AB branch and monotonically decreases Iα¨ ≈−FRα. (18)
on the CD branch. We need to elucidate the underly-
ing causes leading to the switching between the jumping Note that for the low inertia case, sinα≈α approxima-
mode and monotonically increasing or decreasing mode. tionisclearlyjustified[seeEq. (7)]. Usingthisequation,
(b) For large V, say V = 4 and large I (Fig. 5), the we first get an idea of the relevant time scales as I is
9
increased. consideringthe approximationused(i.e.,using themean
F). AbetterestimatecanbeobtainedbyusingEq. (19).
Case a For the high I and V case, Fig. 5 for V = 4 shows
Consider the low velocity branch AB where the small that the wave forms are nearly sinusoidal except for a
amplitude high frequency oscillations are seen on the jitter at the top and bottom. For this case, f(v,4) is
nearly linearly increasing part of F [given by F(t) = nearly flat over the entire range of values of v, with a
F +k(V −v)t, see for instance Fig. 3(b)]. A rough value ∼ 300. Here, even on the AB branch, we can not
min
estimate of this time spent on this branch is obtained ignore the v˙ term in Eq. (15). However, one sees that
by (f − f )/kV ∼ t. Using f ∼ 284 and as v = 0.335 and v = 1.25 which suggest that
max min max min max
f ∼ 200, [from Fig. 3(b)], we get t = 0.084 (com- to the leading order, we could ignore the v˙ term. This
min
pared to the correct value of 0.063 which we shall ob- gives the period T =0.115. From Fig. 5(b), considering
tain soon) which is much larger than the period of the onlythe monotonicallydecreasingpart(ab),the valueof
high frequency oscillation. Thus, we could take the lo- T/2=0.53 read off from the figure compares reasonably
cal value F for the purpose of calculating the period well with this value.
of the high frequency oscillation. Consider the orbit For the CD branch, as v is not small, the v˙/R term
at the lowest value of F for which we can use F ∼ appearsto be importantin Eq. (15). Some idea of when
min
f (v,1) ∼ 200. Then using Eq. (18), the frequency this termis importantcanbe hadbylookingatthe time
min
ν = (FR/I)/2π = 225 for I = 10−5 which gives the scales arising from inertia, namely, FR/I and the coeffi-
periopd of oscillation T = 4.44×10−3. This agrees very cient of the damping term, F/Rf′ in Eq. (17). Consider
wellwiththeexactnumericalvalueT =4.1×10−3. This V = 1 for I = 10−5 and 10−2. The period obtained
frequency decreases when the force reaches the maxi- by assuming the mean value of F = 240 in FR/I gives
mum value F ∼ f (v,1) ∼ 284 to ν = 261 giv- 4×10−3 for I = 10−5 compared to 0.128 for I = 10−2.
max max
ing T = 3.69×10−3 which is again surprisingly close to These numbers can be compared with the time scale
the numerical value 3.72×10−3. In the numerical so- Rf′/F which is 0.01 (where we have used f′ ∼ 25 from
lutions, we find that the period gradually decreases [see numerical simulations for V = 1). This shows that for
Fig. 3(c)]. This feature is also easily recovered by using highinertiathedampingcoefficientF/Rf′ inEq. (17)is
F =F +k(V −v)t. This leads to an additional term important. Wewilldiscussthisissueinmoredetaillater.
min
in the equation of motion for α in Eq. (18),
Case b
Iα¨ =−F(t)Rα/I =−[Fmin+k(V −v)t]Rα/I, (19) Now we focus on the origin of jumps between the
branches. We note that the jumps from CD to AB (or
wheretisthe time requiredforF toreachF starting
max vice versa) occur only when the peel velocity v reaches
from F . Here again the v term can be dropped. If
min a value where f′(v,V) = 0. This also means that the
F wasabsent,theequationhastheAiry’sform. (Note
min timescaleoneachbranch,whetheritspendsonlyashort
that for this case also we could assume F ∼ f
min min timeornot,iscontrolledbytheequationforv. However,
and F ∼ f .) Though this equation does not
max max clearlytheinfluenceofinertianeedstobeincluded. Here
have an exact solution, we note that we could take α
we present an approximate equation for v which is valid
to have a sinusoidal form with 2πν = (FR/I) where
in the various limits of the parameters:
F is treated as a slowly increasing parapmeter. (This as-
sumptionworksquitewell.) Theaboveequationcaptures v˙ = [F˙(t)(1+α(t))+F(t)α˙]/f′, (20)
the essential features of the numerical solution. The nu- [k(V −v)+(F +k(V −v)t)α˙)]
merical solution of Eq. (19) ( as also this representation ≃ in , (21)
f′
) gives the decreasing trend of the small amplitude high
frequencyoscillations. (NotethattheAiryequationitself where the time t is time spent on the branch consid-
gives a decreasing amplitude [24].) ered ( low or high V). In Eq. (21), we have again used
WenotethatEq. (18)isvalidontheABbranchwhere F˙ ≃ k(V −v) and F ≃ F +k(V −v)t with the same
in
v is small even for high inertia and small V case. Thus, approximationused in Eq. (17).
we may be able to recover the gross time scales using Wenowattempttoobtaincorrectestimatesofthetime
this equation. Our numerical results show that as we spentby the orbitoneachbranchstartingwiththe least
increase the inertia, α exhibits a sinusoidal form on the complicatedsituationofthelowinertiaandsmallV. For
AB branch [see Fig. 4(b)], although one full cycle is not this case, on the low velocity branch, one can use the
seen. We note that though the value of α is much larger sinusoidal solution for α, namely α = α sin(2πνt+φ),
in
thanthatforsmallI,wecanstillusetheaboveequations whereφisaphasefactorwhichalsoincludesthecontribu-
[Eq. (18) and (19)]. On this branch F increases from a tion arising from the jump as well and 2πν = (FR/I)
value Fmin ∼ fmin to a maximum Fmax ∼ fmax. For with F ≃ Fin+kVt . Both αin and φ needs tpo be sup-
large I = 10−2 ( and V = 1), we get a rough estimate plied. Alternately, one can use Eq. (18) with Eq. (21)
of the period by using the mean value of F ∼240 in Eq. for which we provide α and α˙ at the point from the
in in
18. This gives a period T = 0.128 which already agrees exact numericalsolutions. We stress that this procedure
satisfactorily with the numerically exact value T =0.11 is not equivalent to solving all the equations,as the only
10
equation we use is Eq. (21) with the form of α already 0.6
determined from the equation for α. [We note here that 16
thoughwe have usedthe sinusoidalformof α alongwith
the initial conditions on α ,φ, it is simpler to supply
in
the initial conditions α ,α˙ and use Eq. (18).] We
in in 0.4 v
note here that f′ is a crucial factor that determines the
time at which the orbit jumps from one branch to the
v
other. Equation (21) needs to be integrated from v to
in
8
v thataredeterminedbythepullingvelocityV,i.e.,the
f 0.2 0 t 0.005
form of f(v,V).
For the low v branch f′ term makes a significant con-
tribution for the time spent by the trajectory on AB.
Indeed, one can obtain the order of magnitude of the
0
time spent by the orbit on AB by using a crude approx-
0 0.02 t 0.04 0.06
imation for f′(v,1) = −230(0.5−v). This can be easily
integrated from v =v ∼0.0188, to v =v ∼0.4 which
in f
FIG. 8: Comparison of approximate solution (continuous)
already gives ∆t=0.075. This number is comparable to
with the numerically exact solution (dotted) of v(t) for I =
thenumericallyexactvalue0.063. Acorrectestimatecan 10−5 andV =1fortheABbranch. Theinsetshowsasimilar
beobtainedbyusingf′fromEq. (14)withthesinusoidal
comparison of v(t) on theCD branch. (v, V are in m/s, I in
formαorEq. (18). (We haveusedFin =211.5fromthe kgm2 and t in s.)
numerical simulations for V =1 and 2πν = [RF(t)/I]
with F = Fin+k(V −v)t.) This gives nearply the exact
numericalvalue of ∆t=0.063. In fact, this solution also of F on the AB branch.
captures the oscillatory growth nature of v quite accu- In this case, as already discussed, the coefficient of α˙,
rately. The approximate form of v(t) (continuous line) namely, F/Rf′ term in Eq. (17) determines the time
along with the numerically exact solution ( dotted line) scale on CD, while on AB, the term FRα/I dominates.
are shown in Fig. 8. Using ∆t in F =F +k(V −v)∆t Thus, the general equation valid for this case is
in
gives ∆F = 63 and F = 274.5 which is in good agree-
FRα Fα˙
ment with the exact numerical value of Fmax = 275. It α¨ ≃− − , (22)
I Rf′(v)
is interesting to note that this value is much less than
fmax =283[seeFig. 3(b)]orequivalently∆F islessthan where we use F = F +k(V −v)t. [Note that we have
fmax−fmin, what is also observed in our exact numeri- dropped k(V −v) terinm from Eq. (17) as this term does
calsimulation. The underlyingmechanismofjumpingof
not have any dependence on α or α˙.]
the orbit before F reaches f also becomes clear from
max We start with the cascading effect. Consider the orbit
the analysis(Fig. 8). We notethatthe magnitudeofthe
whenit is atthe highestvalue ofF =295.6onthe CD
in
oscillatorycomponentinv growstillitreachesvmax per- branch for which we can drop FRα/I term. As f′ is a
mitted by f(v,1). Then, the orbit has to jump to CD.
functionofv,andF alsodependsontime,itappearsthat
Thus, the approximate solution gives an insight into the weneedtousecoupledequationsα˙ =−Fα/Rf′withEq.
cause of the orbit jumping even before F reaches f
max (21). However,the numericalsolutionofthese equations
(for small I).
showthatonecanmakefurtherapproximationbytaking
For the CD branch also, the dominant term is f′. In- f′ to be constant taken at v = 15.54 and F = F , as
in
deed, any reasonable function which has the same geo- the time spent on this branch is very small. The error
metrical form of f shown in Fig. 2 will give good results inusingthis approximationiswithin10%. Indeed, using
for ∆t. Using the correct form of f′, we get ∆t = 0.005 α = −0.0304,α˙ = −160 and numerically integrat-
in in
which is close to the exact result. This again gives cor- ing Eq. (21), along with Eq. (22) from v = 15.54 to
in
rect magnitude of ∆F =72.5. In addition the nature of v =8.7 gives ∆t=1.78×10−3. This compares reason-
f
the v(t) obtained by this approximation is close to the ably with the numerical value of 1.89×10−3. Using this
exact numerical solution shown in the inset of Fig. 8. we get ∆F =19.4 which compares well with the numer-
Case II, intermediate and high I and low V ically exact value 19.8. At this point the orbit jumps to
The most difficult feature of our numerical solutions the low velocity AB branch (to the point e). Thus, as
to understand is the dynamical mechanism leading to a ∆t is small, for all practical purposes, we can ignore the
series of drops in the pull force seen on the descending dependence of f on v and F on t and use α to be an ex-
branchofF(t)forintermediateandhighvaluesofinertia ponentially decreasing function for analytical estimates.
andforarangeofV values. Considerthehighinertiaand These analyticalestimates alreadygive reasonablyaccu-
low V case (say I = 10−2 and V = 1) shown in Fig. 4. rate numbers.
As stated earlier,there aretwo differentissues thatneed On the AB branch, the dominant time scale is deter-
to be understood here. First, the series of small force mined by FRα/I, andwe canuse the approximate sinu-
drops ∆F and second the monotonic increasing nature soidal form in Eq. (21), or Eq. (18) along with Eq. (21)