PHYSICA® Physica D 121 (1998) 193-212
ELSEVIER
Plane shear waves under a periodic boundary disturbance in a saturated granular medium Brian T. Hayes *, David G. Schaeffer 1 Department of Mathematics and Centerfor Nonlinear and Complex Systems Duke University, Durham, NC 27708-0320, USA Received 13 October 1997; accepted 27 March 1998 Communicated by R.E Behringer
Abstract We study an initial-boundary value problem for a fully nonlinear system which describes the propagation of plane shear waves in a three-dimensional soil, assuming a hypoplastic flow rule. This model was derived to study liquefaction of soils. For a periodic square-wave velocity disturbance on the boundary, the solution is found to saturate away from the boundary. We calculate the asymptotic state, in terms of the initial and boundary data. This formula exhibits fractal behavior as the shape of the square wave is varied. Our analysis is confirmed by numerical computations. © 1998 Elsevier Science B.V.
Keywords: Initial-boundary value problems; Fully nonlinear partial differential equations; Shock waves; Fractals
1. Introduction In this p a p e r w e study an initial-boundary value p r o b l e m in the quarter-plane {x, t > 0} for the system
3tV = 3xCr,
3tcr = a 3xV + b IOxvl,
(1)
w h e r e the constants a and b satisfy a > b > 0. W e a s s u m e constant initial conditions,
v(x,O)=O
and
cr(x,0)=cr0<0,
(2)
and a periodic square w a v e as b o u n d a r y data, v(0, t)
f(t)
f o r n = 0, 1, 2 . . . . .
{ v_, v+,
n r < t < (n + l - )~)r, (n + 1 -- ;~)~ < t < (n + 1)z,
(3)
In f o r m u l a (3), v+ and v_ are constants with v+ > 0 > v _ ; r is the period; and X c [0, 1]
specifies the fraction o f t i m e spent at v+. * Corresponding author. 1 Research supported by NSF Grant DMS-9504577. 0167-2789/98/$19.00 © 1998 Elsevier Science B .V. All rights reserved Pll SO 1 6 7 - 2 7 8 9 ( 9 8 ) 0 0 1 5 4 - 7
194
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
Region 4
Region 3
I
/ /
Re/gi0.m-25 5 / ~
/
.
! I i
I
i
J
i
!
! ! !
f
t
i
Region 1
X
XI Fig. 1. Schematic representation of the asymptotic solution.
System (1) was derived by Osinov and Gudehus [10] in studying spontaneous liquefaction of soils. Specifically, these equations arise in the description of small-amplitude plane shear waves propagating in a three-dimensional saturated soil, assuming a hypoplastic flow rule. System (1) also applies to purely longitudinal waves in a plastic solid, see [8,9]. The dependent variables v and o- in (1) represent a component of velocity and a component of stress, 2 respectively. The other two components of velocity vanish, and the other components of stress can be determined from knowledge of v (x, t). Although (1) is fully nonlinear, this system reduces to the linear wave equation if the sign of Oxv is definite: Ottl) -- C20xx v = O,
¢2 = a 4- b,
(4)
c 2s = a -- b > 0,
(5)
if ax v _~ 0, and Ott I) - C20xxV = 0,
if OxV < O. In problems where both signs of 0xv are present, waves with four distinct speeds, 4-cf = -4-~/-d 4- b and 4-Cs = ±~/'S - b, may be produced. Osinov and Gudehus [10] solved the quarter-plane problem for (1) in the case of boundary data v(0, t) consisting o f a single negative, trapezoidal pulse. This boundary data generated a kind of traveling-wave solution v ( x , t) propagating into {x > 0}, an expanding trapezoidal pulse w h o s e leading edge m o v e d faster than its trailing edge. In the case of a positive, trapezoidal pulse, the trailing edge of the pulse overtakes its leading edge, resulting in a
complicated set of wave interactions. Gordon et al. [5] showed that for large times, the solution in this case assumed an inverted profile, resembling that of Osinov and Gudehus. 2 Incidentally, as was shown in [10], the sign of cr0 in (2) determines the appropriate sign for b in (1); if ~r0 were positive, b would be negative.
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
195
Stress Profile for the Plane Shear Wave Model
Velocity Profile for the Plane Shear Wave Model
-11
0.5
-12 tl.
0
-13[
-0.5 -14[ N
~ -15
~ -1 it ff
~-16
-1.5
.o
-17 t
-2
-3
(.)
-18[
/
-2.5
50
10 150 x, a=4, b=S, tau=2, [ambda--0.4
-19F 200
250
-20
i
0
x,
i
100 150 a=4, b=S, tau=2, lambda=O.4
i 200
250
(b)
Fig. 2. (a) Numerically computed v(x, t = 73.2) for a = 4, b = 3, r = 2, v+ = 1 = - v _ , ;k = 0.4. (b) Numerically computed (T(x, t = 73.2).
For the time-periodic boundary disturbance we analyze here, the solution of (1) exhibits qualitatively different behavior at large times, in four regions. This is illustrated in Fig. 1. Proceeding from right to left, there are: 1. a quiescent region in which the effect of the boundary disturbance has not yet been felt, 2. a sequence of fast waves of decreasing strength moving to the right, 3. a uniform region in which the stress and velocity are constant, independent of x and t, and 4. a strip of constant width along the boundary with complicated x-dependence and periodic t-dependence. All wave interactions, even prior to the asymptotic regime, are confined to Region 4; reflecting this fact, we shall call Region 4 t h e interaction region. A typical solution, as a function of x, for a fixed (large) time, is graphed in Fig. 2. In this paper we derive a formula for how the uniform state (v~c, cry) in Region 3 depends on the boundary data for v. This formula reflects behavior found by Osinov and Gudehus: i.e., they showed that stress is reduced behind a shear wave produced by a disturbance at the boundary. However, under the repeated disturbances considered here, this stress reduction does not continue indefinitely but rather "saturates"; the uniform state in Region 3 arises in this way. 3 The uniform state ( v ~ , ~r~) varies in a strikingly discontinuous fashion with the boundary data. For example, Fig. 3(a) shows cr~ as a function of )~ for a = 4 and b = 3. This graph has a fractal, scale-invariant character, as m a y be seen in the detail of ~r~ ()~) in Fig. 3(b). Moreover, the graph of ~r~ vs )~ depends strongly on the parameters a and b of (1). Fig. 9(a) of Section 4, for example, shows the graph of a~z vs )~ for a = 1.625 and b -----0.625; this plot bears little resemblance to Fig. 3. The remainder of this paper is divided into three sections plus an appendix. In Section 2, we discuss discontinuous solutions of (1) and show how they can be used to solve (1)-(3). In Section 3, we state and prove our main result, identifying the four different regions and computing the uniform state ( v ~ , crc~) in Region 3. Section 4 3 This saturation is perhaps surprising, in view of the fact that under uniform deformation, hypoplastic constitutive laws exhibit "ratcheting". Indeed, intergranular strain [7] was incorporated into hypoplasticity, precisely to eliminate such unphysical ratcheting. Our constitutive law does not include intergranular strain, yet the solution of (1)-(3) saturates.
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
196
Dependence of the Asymptotic Stress on the Plus Fraction Detail of "sigma vs lambda" Plot
-4
-11.4
-5
l
-7
-11.6
-11.8
-6
-12 .£
~
E
-12.2
-12.4
-12,61
-12 -13
. . 0.2. .
0.1
0.3
i
'
i
0.4 05 0.6 07 lambda, the plus fraction
i
08
i
09
(a)
-%1
0.32
0.34
i 0,36
Plus Fraction (lambda)
i
0.38
0.4
(b) Fig. 3. (a) ~ro~vs ~., for a = 4, b = 3, ~ = 2, and v+ = 1 = - v_. (b) Detail showing the fractal structure of aoo 00.
contains some concluding remarks, while existence and uniqueness of the Riemann solution are discussed in Appendix A.
2. Construction of the solution If ~b is an arbitrary C 2 function of a real variable, then
v(x, t) = (a(x -- cft)
(6)
is a traveling-wave solution of the linear wave equation (4). If, moreover, ¢ i s nondecreasing, then Oxv >_ O, and we may obtain a solution of (1) by defining
(v(x,t))=(qb(x-cft) ) a ( x , t) ] a. - cfqb(x - cft) '
(7)
where a . is an arbitrary constant. If q~ is not differentiable, perhaps not even continuous, but still nondecreasing, then (7) is a weak solution of (1). In particular, suppose ¢ is piecewise constant with one jump; in symbols, ¢(~)
[VR VL
i f ~ >~0, i f ~ < ~0,
(8)
where VL < YR. Then substitution into (7) gives a piecewise constant solution of (1) that jumps across the line
{x -- eft = ~0};
(9)
in other words, the solution has a discontinuity that moves to the right with speed cf. We shall refer to such a discontinuity as a fast positive wave. Likewise, still assuming VL < VR in (8), we may substitute into
V(X, t)
~r(x,t)) = (~,
q~(x -]- eft) q- Cf(b(X -[- cft) ]
(10)
197
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
V
2 3
!
V
2
S_
s.
1
3
1
!
F_ (a)
(b)
At F
U
/
/~SF
\
................. ~'~
UR = (VR,~ R)
......../
/
1
UL
t = t ° .................................................................................
........................................ X=X 0
(c)
-~.
F_/
X
~S+
(d)
Fig. 4. (a) Wave curves in the (v, a)-plane. (b) Riemann solution in the (v, a)-plane. (c) Case with slow negative and fast positive emitted waves. (d) Examples of Riemann solutions containing three waves.
to obtain a fast n e g a t i v e wave. If VL > VR in (8), similar considerations with
cr(x,t)
=
a . ±Cs~b(X ± C s t )
(11)
lead to s l o w waves, either positive or negative according as the m i n u s or plus sign is chosen in (11), respectively. It is c o n v e n i e n t to s u m m a r i z e these solutions in a graphical representation. Consider Fig. 4(a), in which the (v, a ) - p l a n e is divided into four sectors, with UR = (VR, fiR) at their c o m m o n vertex. The dividing lines, called
198
B.T. Hayes, D.G, Schaeffer/Physica D 121 (1998) 193-212
wave curves, have the equations •f'4- = {(V,O'): r7 =O"R 4- Cf(VR -- V), V ~ VR}, S-t- = {(V,O'): cr = o r R-4-Cs(vR- v), v > VR}. If, for example, (v, rr) lies on the curve ~ + , then there is a solution of (1) in which (v, or) and (VR, erR) are separated by a fast positive wave, with (VR, O'R) lying on the right. Similarly ~ _ and S± are associated with fast negative and slow waves, respectively. In all cases, the wave speeds are given by the Rankine-Hugoniot condition (cL [13]): s --
ryR -- o-L UR
-
-
VL
.
(12)
In preparation for our treatment of (1)-(3), we now analyze the pure initial-value problem for system (1) on {t > to} with the initial conditions u(x, to)
[VR / VL
i f x >X0, if x < X0,
and
cr(x, t 0 ) = /crR [ O"L
ifx >x0, if X < X0.
(13)
Such an initial-value problem with piecewise constant initial data is called aRiemannproblem; see [13] for a general discussion. Regarding existence, this problem may be solved by combining one positive and one negative wave, as illustrated in Fig. 4(b). The intermediate state (vi, at) in the figure must be connected to (VL, aL) by a negative wave and to (VR, OR) by a positive wave. Such a connection, in the case of a slow negative and a fast positive wave is indicated in Figs. 4(b) and (c). More generally, the solution consists of slow-slow, slow-fast, fast-fast, or fast-slow waves if (VL, aL) lies in sector number 1, 2, 3 or 4, respectively. Of course, if (VL, rrL) lies on one of the boundaries ~c± or ,5+, then one of the two waves will have zero amplitude. Regarding uniqueness, if one allowed both slow and fast positive waves in the same solution, or both slow and fast negative waves, then (1) with initial data (13) would admit infinitely many solutions! In Fig. 4(d), for example, with UL as before, there are solutions containing three waves. One solution (dashed lines) contains a ~ + wave from UR to UI2, followed by a ,9+ waves from UI2 to UI1, and finally a S - wave from UI1 to UL. A second solution (dash-dot lines) is depicted with a wave sequence of .T'+S_.T'_. In Appendix A, however, we argue that multiple-wave Riemann solutions which contain both slow and fast waves of the same sign are unphysical. Specifically, centered (i.e., x / t self-similar) profiles for a viscous regularization of (1) cannot be constructed. Once these multiple-wave solutions are rejected, there remains the unique solution of (1) and (13), containing precisely one positive and one negative wave. We utilize this solution in solving (1)-(3). Similar considerations apply to the initial/boundary-value problem for system (1) on the quarter plane {x > 0, t > to}, say with initial conditions
V(X, to) = I)R and
tr(x, to) = trR,
(14)
and the boundary data v(0, t) = vB.
(15)
The solution consists of a single positive wave along x = c (t - to), where c equals cf or Cs, according as VB < VR or VB > VR, respectively. The stress trB following the shock is determined from the requirement that (VR, CrR) and (vB, rrB) must be connected by a positive wave: i.e., (vm trB) must lie on the positive wave curve f + U ,5+ through (VR, rrR) in Fig. 4(b).
199
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212 G
x--O
V ----VB
speed = c s
/
v~-vB
(VB,%)
V
(VR,~,R2)< %)~t~ (vB,
~s) speed = - c
(a)
S
co)
3
Fig. 5. (a) Type 2 interaction, shown in (x, t)-plane. (b) The same interaction, shown in (v, cr)-plane. We now solve (1)-(3). The solution will be piecewise constant, with jumps along waves of the various types. 4 To start the construction, consider the solution of (1) with data (14) and (15) where we choose OR = 0'
O'R ~ O'0,
VB -~- V_.
Observe that these initial conditions coincide with (2) and, for 0 < t < (1 - X)r, this boundary condition coincides with (3). Hence, the solution of (1) with these data solves the original problem for 0 < t < (1 - X)r. Thus, in this time interval, the solution of (1)-(3) contains a single positive wave originating at x = t = 0, a fast wave since VB = v _ < 0 = V
R.
Suppose that a piecewise constant solution of (1)-(3) has been constructed for 0 < t < to. More precisely, suppose that for to - 3 < t < to there are jumps in this solution along x = ai -t- ci (t - - to),
i = 1.....
N,
where ci equals -t-cf or 4-Cs. For most values of to the following three conditions hold: 1. The boundary data (3) is continuous at to: i.e., to is not equal to (n - X)r or n r , for n = 1, 2 . . . . 2. All discontinuities of the solution are in the interior of {x > 0}: i.e., ai > 0 for all i. 3. No two shocks coincide: i.e., ai 7& a j for i 7~ j . If all three conditions are satisfied, the solution can be extended to t < to + e, for some positive e, by letting the various waves continue to propagate. At an exceptional time of Type 1, this construction fails near x = 0. However, the local behavior of the solution near x = 0, t = to coincides with the solution of (1), (14) and 15) for appropriate data VR, oR, VB. Thus, at such times one additional positive wave enters the solution: a fast wave if t = n r , or a slow wave if t = (n - )~)r, for n= 1,2,... At an exceptional time of Type 2, the simple construction again fails near x = 0, but again the solution may be continued for a small time by solving (1), (14) and (15) for appropriate data. Fig. 5 shows, for example, that an 4 A similar construction was used in [5]. Incidentally, that paper would be greatly simplified by studying a square pulse rather than a trapezoidal pulse
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
200
incident slow negative wave is reflected as a slow positive wave. Likewise, an incident fast negative wave is reflected as a fast positive wave. The number of waves in the solution is not altered at an exceptional time of Type 2. At an exceptional time of Type 3, the simple construction fails at x -----x0, t = to, where x0 -----ai for at least two different values of i. The solution may be continued by solving (1) and (13) for appropriate initial data. Since at least two waves were incident and at most two are emitted, the number of waves is not increased at such points. Thus, if t < n~:, there are at most 2n waves in the solution; i.e., the waves originating from discontinuities in the boundary data. Remark. Collisions may involve either waves of the opposite family, fast and slow waves of the same family, or both, but the construction for continuing the solution is the same in any event. Proceeding in this way, we may construct a solution of (1) and (3) for all time. Indeed, the numerical solution presented in Fig. 2 was obtained by implementing this procedure.
3. Properties of the solution In this section we verify that the solution of (1) and (3) has the asymptotic structure described in Section 1. This structure depends on the two essential parameters in the problem: viz., the plus fraction, L, and a ratio involving the wave speeds, a =
Cf - - Cs
.
(16)
Cf + C s
The formula for the asymptotic state (voo, o-o~) requires the construction of two sequences, as follows: Let )~0 = L; for j = 1, 2 . . . . . let (a) )~j ---- 1 - (o~-l)~j_l - [o~-l)~j_l]}, (17) (b) Nj = 1 + [ o t - l ~ , j _ l ] , where [- ] denotes the greatest-integer function. Note that 0 < )~j _< 1 and 1 _< Nj <_ Nmax, where Nma x :
1 + [oe-1].
The first sequence, {;~j}, may be characterized as the iterates, starting from )~, of the mapping q~ " [0, 1] --+ [0, 1] given by ¢(~) = 1 -
{ot-l~
-
[~-1~]}.
In words, ¢ operates by expanding the interval by a factor of o~- 1, re-mapping onto [0, 1] by subtracting off the integer part, and finally reflecting about the midpoint. This mapping has much in common with the baker's transformation [3]. The second sequence, {Nj}, which arises from the integer parts that are discarded in constructing {)~j}, appears as coefficients in the formula for the asymptotic state: oo
(a)
vc¢ = v_ -- (v+ -- v_) E
Nj a j,
(18)
j=l (b)
a o o = o-0 -
cfvoo,
In deriving this formula we will need to assume that Xj # 1
for j = 1, 2 . . . .
(19)
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
201
For a given value of or, condition (19) excludes a countable set of values of the plus fraction, )~; as we will see, these are the points at which the graphs of aoc, shown in Fig. 3, are discontinuous. In stating our main result, we define the strength g of a wave as the jump in velocity across it; that is, F=IVR--VL[, where VR and VL are the one-sided limits of the velocity on the right and left, respectively. Note that in the following theorem, the four regions are listed in the opposite order from Section 1. Theorem 3.1.
(a) All wave interactions are confined to the strip {0 < x < x i } , where xi ---- { C f C s / ( C f - Cs)} "C. (b) If (19) holds, then for any fixed x > xb as t ~ e~, the solution tends to the constant, asymptotic value, (voo, o'oo), given by formula (18). (c) The fast positive waves emitted at the boundary at time n r , for n = 0, 1 . . . . . propagate through the interaction region {0 < x < xi}; no other waves escape the interaction region. If Yn denotes the strength, following escape, of the wave along {x = c f ( t - n r ) } , then Fn decays geometrically in n. More precisely, if (19) holds, then (20)
}In ~ (Uq- -- V--)Ol n/Nmax .
(d) For x > c f t , the solution remains in its quiescent, initial state (0, or0). To prove the theorem, we study wave interactions. We now define three types of interactions which occur away from the boundary: - Type L A fast positive wave of strength g overtakes a slow positive wave of strength 8, where 8 < (1 - o e ) - l g .
(21)
- Type II. A slow negative wave and a fast positive wave collide. - Type III. A slow negative wave and a slow positive wave collide.
All interior interactions, i.e. x > 0, belong to one of the above three types. In particular, no fast negative waves arise in the solution of (1)-(3), and every collision of a fast positive wave with a slow positive wave is of Type I, with condition (21) being satisfied.
L e m m a 3.2.
We defer the proof of this lemma until near the end of the argument. To begin, we analyze these three types of interactions. As illustrated in Fig. 6(a), an interaction of Type I results in a slow negative wave and a fast positive wave, of diminished strengths
L e m m a 3.3.
8~ = o t 8
and
F~=y-(1-c~)8,
(22)
respectively. In interactions of Types II and l]I, both waves continue without change in type or strength. Remark. Condition (21) guarantees that F ~ > O.
P r o o f o f L e m m a 3.3. ~,' =
Vl - - V3,
Label the four states adjacent to an interaction of Type I as in Fig. 6(a). Note that 8 =
Vl - - V2,
yr =
V2--
V4,
81 =
V3 - - V 4 -
202
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
----4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.~.
V
3 -C S
4
cf
4
s .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
>.
x
(a)
0~) Fig. 6. (a) Type I collision in the (x, t)-plane. (b) The same collision, shown in the (v, ~r)-plane.
To continue the solution beyond the time of interaction we solve the Riemann problem (1, 13) with data (/)R, 0-R) = (112, 0-2)
and
(VL, 0-L) = (113, 03).
Fig. 6(b) shows the relevant wave curve diagram, and locates the appropriate intermediate state (v4, 0-4) needed to solve the problem. Notice that the intersecting wave curves form a trapezoid. The jump relations are: ~3--0-1 0-4--0-2
=--Cf(V3 --Vl), =--Cf(V4--V2),
0"4--0-3 ="}-Cs(V4--V3), ff2--ffl m---¢s(V2--Vl)-
(23)
Adding the first two equations in (23), then adding the second two equations in (23), and setting their right-hand sides equal to each other gives the result for ~ in (22). From this, the result for yt quickly follows. We omit the analysis of interactions of Types II and III, which is similar but simpler. In these cases the wave-curve diagram is a parallelogram rather than a trapezoid. [] Remark. If constraint (21) is violated, the collision of a fast positive wave with a slow positive wave (no longer a Type I interaction) produces a slow positive wave, as illustrated in Fig. 7; as we show, below, in the proof of Lemma 3.2, this situation never arises.
Recall from Section 2 the three exceptions to ordinary wave propagation: 1. Emission of a new positive wave at a discontinuity of the boundary data; 2. Reflection of a negative wave at the boundary as a positive wave; and 3. Interaction of two waves. According to Lemma 3.3, interactions of Types II and III have no effect on wave types or strengths, while interactions of Type I are in some ways analogous to reflections at the boundary: i.e., an incident slow positive wave is reemitted as a slow negative wave, albeit with its strength attenuated by a factor of or. At both reflections and
203
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
V
-
2
Cs
~ CS :
1 Fig. 7. Wave curves, with an emitted slow positive wave, if condition (21) does not hold.
interactions o f Type I, we regard the reemitted wave as the continuation of the incident wave. With this convention, waves originate only at the boundary. Let us consider the trajectory o f the first slow positive wave, which is emitted at time (1 - )~)r and c o n t i n u e d in this manner. As illustrated in Fig. 8, this wave propagates until it is overtaken b y the fast positive wave emitted at time 3, at which point it gets converted to a slow negative wave; this in turn propagates until it hits the boundary, simply passing through any slow positive and fast positive waves it encounters, and is reflected as a slow positive wave. This reflected wave is overtaken b y the next fast positive wave and thereby converted to a slow negative wave; and so on. Thus, the trajectory o f this wave is an infinite b r o k e n line, say x = X ( t ) , each segment h a v i n g speed -4-cs. A slow wave emitted at a later time, say at (k + 1 - )~)r, undergoes the same sequence of interactions and reflections, except all events are delayed b y k r ; i.e., the wave has trajectory x = X ( t - k r ) . Indeed, since waves originate o n l y at the boundary, {x = X ( t - kv): k = 0, 1, 2 . . . . } is a complete inventory 5 of the slow waves w h i c h occur in constructing the solution of (1)-(3). The sequences {Lj } and {Nj } arise in the quantitative description of the slow wave trajectories. After the slow wave emitted at time (1 - )~0)r collides with the fast wave emitted at time 1:, the slow wave returns to x = 0 at time tl = (1 -t- ot-1)~0)z". Note that N l V < tl < (N1 + 1)v,
5 This is, in fact, an overnumeration if )~j = )~k for some j 5~ k; in this case the sequence {)~j} is periodic, and some of these waves lie on top of one another. Since the strengths of the slow waves are additive, however, we may still consider the waves as separate entities, whose trajectories overlap.
204
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
t
s ~
t 1 = (1+ ~ , / ~ )g NI'I: -C s
3"~
s ~
, ~
~
,
~
~'Cf
2"c
( 1 - ~. ) '17
0
x
Fig. 8. Trajectory of first slow wave; other slow waves not shown•
where N1 = 1 + [ot-1;v0]; by (19) the inequality is strict. Let )Vl = (N1 + 1)r -- tl = (1 - {ot-l~.0 - - [0t-l~.0]}) "C SO that time )q z elapses before the next emission of a fast wave. Repeating the calculation, we see that the reflected wave emerging at time tl returns to x = 0 at time t2 = (N1 + 1 + ~z-1;Vl) r. We define N2 and )~2 analogously. Comparing this process with (17), we conclude that the slow wave emitted at time (1 - )v)r returns to x = 0, after j interactions of Type I, at time
tj=(NI+...+Nj+I--*%j)
r,
j=l,2
.....
Each time it returns to x -----0, however, the slow wave is weaker by a factor of oL. Let us pause to verify Part (a) of Theorem 3.1. The reflection of the slow wave at time emission of a fast wave an interval kj r later. These collide at the point xj-
CfCs
- - z
tj
is followed by the
jr.
c f -- Cs
Since ~.j < 1, all points of interaction lie in the strip {0 _< x < xi}, as claimed. In contrast to the slow waves, the fast (positive) waves are never deflected from their original paths. In order to calculate their strengths, Yn, we need to count the number of interactions of Type I along the fast waves. Every fast wave, except the zeroth wave, suffers at least one such interaction. Specifically, the fast wave emitted at time n r collides with the slow wave emitted at time (n - ~ ) r and, moreover, this is the first interaction for the slow wave. Now, the first reflected slow positive wave appears at time tl, where N l r < tl < ( N 1 ÷ 1)r. Thus, fast waves numbered between 1 and N1 suffer exactly one interaction of Type I, while every later wave suffers at least two. A twice-reflected wave enters the solution at time t2, where N2r < t2 < (N2 + 1)1:. Thus, fast waves numbered
B.T Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
205
between N1 + 1 and N1 + N2 suffer exactly two interactions of Type I, and every later wave suffers at least three. In general, fast waves numbered between N1 + . . . + N j - 1 + 1 and N1 ÷ . . . + N j suffer exactly j interactions of Type I. The initial fast wave has strength Iv_ I. We claim that for 1 < n < N1.
Yn = (v+ - v _ ) a
This follows from substitution into (22) since each of these waves starts with strength v+ - v_, and has one Type I interaction with a slow wave of strength v+ - v_. Thus, each of the fast positive waves 1 . . . . . N1 leaves the interaction region with strength (v+ - v_)ot. Next we claim that Fn=(V+-v-)a
2
forNl+l
Each of these waves has two interactions with slow waves, one of strength v+ - v_ and one of strength (v+ - v _ ) a . Regardless of the order along the fast waves in which the interactions occur, substitution of these data into (22) verifies the claim. In general, ~/n = (v+ - v - ) o t j
forNl+..-+Nj_l+l
< n < NI + . . . + Nj.
Grouping equal terms, we find that Z F n=0
n =
[V_[ +
(V+-
v_)ZNjotJ.
(24)
j=l
The asymptotic state ( v ~ , o-~) is separated from the quiescent state (0, a0) by a sequence o f fast positive waves o f strengths {Fn }- Thus v ~ = - ~n~_0 Fn, and we deduce (18a) on recalling (24). Eq. (18b) follows from (18a) and the j u m p relations. This proves Part (b) of Theorem 3.1. P r o o f o f L e m m a 3.2. T h e process described above - emission of a slow positive wave at x = 0, t = (n - )~)~,
interception by a fast positive wave from x = 0, t = n r , subsequent deflection and undisturbed return to x -- 0, then reflection as a slow positive wave - continues a d infinitum, unless condition (21) is violated. We now show that this condition always holds in solving (1)-(3), using a double-induction argument. The first induction argument is on the fast waves. Along the n = 1 fast wave, there is one Type I interaction, with y = 6 = v+ - v_, so condition (21) holds for this wave. Consider the nth fast positive wave, where N1 + N2 + • .. -4- N j - 1 + 1 < n < N1 + . . . + Nj, and assume that all interactions between fast positive and slow positive waves are of Type I, along previous fast positive waves 1, 2 . . . . . n - 1. The second induction is on Type I interactions along the nth fast positive wave. By repeated application of (22), we observe that after m interactions of Type I, a slow wave has its original strength of v+ - v_ attenuated by a factor o f am. Condition (21) is therefore satisfied at the first Type I interaction along the nth fast positive wave, where Y = v+ - v_ and 6 = (v+ -- v _ ) a m, m ~ [0, j - 1]. Suppose that (21) holds for the first p interactions of Type I, along the fast positive wave; here p _< j - 1. Following the p t h interaction, the fast positive wave's strength has been reduced to (v+ - v_)[ 1 - (1 - o0 (o~il + - . . + a i p ) ] , according to (22), with {in} C {0, 1 . . . . . j - 1}. To verify (21) at the ( p + 1)th interaction, we must show that (Yq- - - V_)Ot m <: (1 - - Ot)--l(uq_ - - V - - ) [ 1 - - (1 - - 01)(O/il -{- • • - -[-
olip)].
Since m ¢( {in}, inequality (25) follows by comparison with the geometric series, and the induction continues.
(25) []
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
206
Part (c) of the theorem is trivial if n = 0. If n > 1, we have from (22) that to the right of the last Type I interaction, the strength o f nth fast positive wave satisfies
}'n=(v+-v-)
1-(1-oe)
d
=(v+-v_)ot
j,
(26)
where j is such that N1 + . . - + Nj-1 + 1 < n < N1 + ... + Nj. Since n <_ jNmax, formula (20) follows on substitution into (26). Part (d) of the theorem follows from the observation that cf is the fastest wave speed for the initial-value problem (1)-(3). This concludes the proof o f Theorem 3.1.
4. Concluding remarks Suppose that for some value of ~., say X = A, we have ).j = 1, for some value of j ; i.e., suppose that condition (19) is violated. Then, following their j t h Type I collision, slow waves return to the boundary precisely when a fast wave is being emitted. Unlike the coincidence of slow waves discussed in footnote 5, the coincidence of a fast wave and a slow wave significantly alters the solution. In particular, Nj jumps by one, as X is varied near A . Moreover, the structure of the solution when 2 = A is different from either one-sided limit. In these latter cases, slow waves continue to interact with fast waves, even after their j t h reflection, while for X = A , the slow wave terminates at the point where the j t h reflection would occur; the fast positive wave is emitted in its place. We close with numerical results which further illustrate properties o f the solution. Recall that Fig. 3 showed the dependence of ac~ on )~, when cf = ~/ff and Cs = 1. In Fig. 9(a), we show aoo vs L, again, but now with cf = 1.5 and Cs = 1. All other parameters in (1)-(3) remain unchanged, including r = 2 and v+ = 1 = - v _ . Fig. 9(a) bears little resemblence to Fig. 3, beyond the fact that it, too, has a fractal structure. The plots of v and a , as functions of x, for a fixed (large) t and cf = 1.5, are quite similar to those for cf = ~/~, in Fig. 2; consequently, we do not display these plots. In Fig. 9(b), however, we present a detail of the interaction region for a(x, t) for cf = 1.5, while in Fig. 9(c), we show the convergence in time of v(x, t) to vo,, f o r x = 6, a value just beyond xi. This latter figure shows that the asymptotic solution is rapidly attained.
Appendix A. Existence and uniqueness of the viscosity solution In this appendix, utilizing a version of (1) augmented by viscosity, we argue that Riemann solutions containing precisely one positive and one negative wave are physical, while those containing additional waves are unphysical. To establish the notation we write (1) as
OtU + F (axU) = 0,
(A.1)
where U=(v,a)
T
and
F=-(Oxa,
a S x v + b [ O x V [ ) T.
(A.2)
The flux F is homogeneous since F(c~V) = a F ( V ) , for a > 0. The Riemann initial data for (A.1) and (A.2) is
U (x, O) = { UL = (VL' aL)' UR ~--- (UR, OR),
where VL, • • •, aR are constants.
x 0,
(A.3)
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
207
Dependence of the AsymptoticStress on the Plus Fraction Detail of the Stress Profile Near the Boundary
-14 -17.2 -14.5
-17.3
to -15
-17.4
U st-W
'q -15.5 -16
2~
g
- t 6.5
~[ ~
-1"/
-178 t
"N -17.5
-18
-17"9 I
6'.t
012
613
ol,
-18"
6'.6 017 61s o10
016
0
lambda, the p]l/sfraction
(a)
x,
a=1.625,13=0.625,lambda=0.15
0,) Convergenceto Asymptotic Velocity in Time 6 -0.2 -0.4 -0.6
~--0.8
ii ~" -1 -1.2 -1.4
~
-1,6 -1.8
1~0
2'0
3'o
4~0
time
50
6'0
70
80
(e) Fig. 9. (a) Complicated dependence of cr~ on )~. (b) Detail of stress profile, featuring the interaction region. (c) Convergence of v(x, t) to vc~ in time.
We shall accept as physical a piecewise-constant solution U (x, t ) of (A. 1) and (A.3) if and only if it can be obtained as the limit of self-similar solutions of a viscosity-augmented equation, as the viscosity tends to zero; 6 that is, iff
U(x, t) = lira U~(x, t), ~--~0+
6 Uniqueness questions for nonlinear hyperbolic equations are notoriously difficult. For systems of hyperbolic conservation laws, G l i m m [4] proved existence in 1965, while uniqueness was not proven until 1995, in a paper by Bressan [1]. The need to reject unphysical shocks is an issue even in the most basic treatment of this subject: j u m p s of one sign are allowed in Bttrgers equation, while those with the "wrong sign" are disallowed, in favor of continuous rarefaction waves; see [13]. All of the entropy criteria, used to select physical shocks, are designed to return information to the hyperbolic model, based on (previously neglected) small-scale physical effects. The vanishing-viscosity criterion - the one best suited to the present problem - is one of several such techniques. In problems where dispersion or capillarity play significant roles, see [6] or [ 12], other regularizations m u s t be employed.
208
B.T Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
where U~ (x, t) satisfies the Dafermos-viscosity 7 regularized version of (A. 1): OtUE + F(OxUe) = e t OxxU~,
e > 0.
(A.4)
The assumption of x/t-similarity reduces (A.4) to a system of ordinary differential equations (ODE). In the course of solving these ODE, we show that weak solutions of (A.1) which contain both fast and slow waves of the same family cannot emerge from limits of (A.4), as e --+ 0. Continuing, we construct explicit solutions of (A.4) which do converge to the (distinguished) Riemann solution used in Section 2. Our argument is based on Shearer and Schaeffer [11], although their flux was nonsmooth only when both components vanish, while ours, given by (A.2), is nonsmooth when just one component, OxV, vanishes. We refer to that paper for some details of the analysis. Setting U~ (x, t) = V, (~), where ~ = x / t , yields a system of two second-order ODE: (A.5)
- ~ V ' -b F ( V ' ) = e V " ,
where I = d/d~. In system (A.5) and below, we suppress the e-dependence of VE. The Riemann initial data (A.3) now become the boundary data: V(~) ---->(UL, O'L),
as ~ --+ --OO
V(~) ---->(VR, OR),
as ~ --+ ~ .
(A.6)
We remark that corresponding to the four constants, OL, • • •, O'R in (A.6), there are four free parameters in system (A.5): arbitrary constants may be added to ~ and to each component of V, and both components of V may be multiplied by X > 0, without altering the system. Substituting W(~) = VI(~) in (A.5) leads to the first-order system (A.7)
-~W + F(W) = eW', while the boundary data (A.6) can be expressed as an integral condition: (x)
f
W ( ~ ) d ~ = UR -- UL.
(A.8)
--00
We now rescale both dependent and independent variables in (A.7) by w07) = e~Z/2eW@),
r1 = ~/E,
in order to obtain the autonomous system d w / d~ = F ( w ) .
(A.9)
Note that the directions of w and W = W are the same. Introducing polar coordinates (p, ¢) in the w-plane, w = p?(¢) = p(cos q~, sin ¢)x, we separate the radial and angular components of (A.9): F(w) = p(R(¢)?(¢) + 0(¢)0(¢)), 7 For existence and uniqueness questions involvinga single shock wave, a traveling-waveansatz, i.e., an ODE involving ~ = x - st, with s being the shock speed, is the appropriate choice. At issue here is not the admissibility of a single shock, but rather that of a combination of shocks with different (but constant) speeds. In this case, a more natura/independent variable is ~ = x/t, with each shock moving with speed ~ = const. The standard viscosity coefficient, eUxx does not admit such x/t-similarity solutions, but a nonstandard choice, introduced by Dafermos [2] and having a factor of t multiplying the viscosity, does allow this scaling.
209
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
d~/~
II
,.i
o
<
."
'
;
'
~- d v / ~
<
IV
0
q~l
'2
g
1~3
'4
2r~
(a)
*4 (b)
Fig. 10. (a) Plot of the right-hand side of (A. 10a), O (~b) v s q~. (b) Radial directions in the W-plane. where 0 = ( - sin q~, cos ~p)T. In this way, we obtain the pair of equations: d~b dr/
J 1 -(1+a+b)cos2q~, |l-(l+a-b)cos2~p,
dp = do pR(d?) --
psin2q~ { 1 + a + b , 2 _ 1 + a -b,
Iq~l < : r / 2 , :r/2
(A.10a)
I~bl < :r/2,
7r/2 < Iq~l < 3:r/2.
(A.10b)
We say that Eq. (A. 10a), which is decoupled from (A.10b), has an equilibrium point at ~, if 69 (q3) = 0. As shown in Fig. 10(a), there are four such points: q~l = tan-1 ~/-a+ b, q~3 = a- + tan-1 ~
q~2 = r e - tan-1 ~/a - b , _ b,
~b4 = - tan-1 ~/a + b.
It was shown in [11] that vectors V in these four directions, called radial directions, solve F ( V ) = s V , with s = R (~bj); moreover, for each radial direction, there are piecewise-constant traveling wave solutions of (A. 1), with speed s. Substituting the points ~b = ~bj, j = 1 . . . . . 4, into formula (A.10b) for R(~p), we find that at ~b = ~bl, s = --cf, and comparing this information with Section 2, we see that 4~1 corresponds to a fast negative wave. Similarly, R ( @ ) = Cs (slow positive wave), R(~b3) :- --Cs (slow negative wave) and R(~b4) = cf (fast positive wave). Thus the picture, in Fig. 10(b), of the radial directions in the (v1(~), a1(~))-plane appears as a reflection through the origin of the (v, a)-wave curves drawn previously in Fig. 4(a). 8 From the sign of O between equilibrium points (cf. the arrows in Fig. 10(a)), we observe that the equilibria ~bl and ~b3, corresponding to slow and fast negative waves, respectively, are unstable (repellers). Meanwhile, q52 and ~b4, corresponding to slow and fast positive waves, respectively, are stable (attractors). For any solution, 8 This reflection arises because we solve the boundary value problem (1.1)-(1.3) from right to left. If, instead, we had drawn the phase plane in Fig. 4(a) relative to UL, the wave curves would exactly match the radial directions in Fig. 10(b).
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
210
therefore, ~b(77) is either a constant, ~bj, or else it traces out a trajectory from a repeller to an adjacent attractor. That is, any trajectory is entirely contained within one of the four sectors of Fig. 10(b). This implies that w(q), and hence W(~) = W(~), must also lie within a single sector - possibly even on the boundary between two sectors. F r o m this analysis of Eq. (A.10a) alone, we can already exclude Riemann solutions of (A.1) which contain both slow and fast waves of the same family. As just observed, W(~) must lie entirely within one of the four sectors shown in Fig. 10(b). By the integral condition (A.8), the vector W(~) must belong to the same sector as UR -
U L.
Notice from Fig. 10(b), that in each sector, one o f the derivatives, dr~ d~ or d a / d ~ , has a fixed sign. Specifically, vt(~) > 0 in Sector I, a t > 0 in Sector II, v t < 0 in Sector III and a t < 0 in Sector IV. Thus, for any similarity solution, one component - either v or a - is monotone; assuming the limit as E - + 0 exists, this monotonicity will be preserved under the limit. By contrast, any Riemann solution containing both slow and fast waves of the same family must be nonmonotone in both v and a . This dual non-monotonicity is clear, since the slow and fast wave curves for each family lie in opposite quadrants of the (v, a ) - p l a n e . Thus, such solutions cannot arise as limits of viscous solutions. To show that the distinguished Riemann solution used in Section 2 does arise as the limit o f viscous solutions, we construct explicit solutions of Eqs. (A. 10). Whereas Shearer and Schaeffer [ 11] were constrained by the generality of their fluxes to use formal asymptotics, the specific choice of flux (A.2) we study here allows for an explicit solution of (A.7). However, the formulas obtained for (p (~), ~b(r/)) are more cumbersome in Sectors II and IV, due to the switch from cf to Cs at kbl = re/2. We therefore construct the solution only in Sector I, leaving other cases to the reader. 9 We find from (A.10a) that q~(r/) = tan -1 (cf tanh(k - c f 0 ) ) ,
(A.11)
where cf = ~ - S + b, as before, and k is an arbitrary constant. Notice that the entire trajectory q~ of (A.11) is monotone decreasing: starting at q~(0) = q~0 in Sector I, 4~ ~ 4~4, as ~ --+ cxD, and q~ - + 4~1, as ~ --+ -cx~. From (A.10b) and (A.11), we have p ( o ) = e f R(O(r/))drl = ~. cosh(k -- cfo)~/1 -1- c 2 tanh2(k - cf~/), where ~ > 0 is a constant of integration. Thus, we arrive at the explicit solution for the derivative of UE, W(~) ----~.e-~2/2E cosh(k - cf~/e)V/1 q- c 2 tanh2(k - cf~/e)? (~b@)).
(A.12)
After some trigonometry, we obtain formulas for the components of the derivative: vt(~) = )~e-~2/2~ cosh(k - cf~/e),
(A.13a)
crt(~) = ~.cfe -~2/2E sinh(k - cf~/e).
(A.13b)
Notice that v t > 0, while tr' changes signs. We now integrate Eqs. (A.13), utilizing the integral constraint (A.8) to find k and )~, in terms o f VR -- VL and ~rR -- erE. We find that (
aR_--~.LL )
k = t a n h - 1 ~ , c f ( v R -- VL),] '
(A.14)
9 The solution for W(~) in Sector III may be obtained from that in Sector I by switching the sign of V and replacing cf with Cs.
211
B.T. Hayes, D.G. S c h a e f f e r l P h y s i c a D 121 (1998) 1 9 3 - 2 1 2 Velocity Profiles for Three Values of epsilon
0
Stress Profiles for Three Values of epsilon
// jl ~
? gJ3
///
? o" ii lz >1-9
!
/
xx
l/
Xx
/i
'/
E
t' ' \x
/ t'' / / "
ii
ii
.!
//
9
¥o
I
L,
g-4
/
,'1.'
] I
i /i i
I'
I i
-2
-6
i
-4
i
-3
-2 xi,
e~ilon
i
v
i
-1
0
1
i
2
= 0.5, dash; 0.1, dash-dot;
0.02.
i
i
3
4
--~s
i
-4
solid
(a)
i
-3
i
i
-2 xi, epsilon
-1
i
t
0
1
= 0.5, d a s h ; O.f, dash-dot;
i
2 0.02,
3 solid
(b) Fig. 11. (a) Viscous velocity profiles with ~ = 0.5, 0.1 and 0.02. (b) Viscous stress profiles with ~ = 0.5, 0.1 and 0.02.
which is well-defined, since ]OR -- OL] ~ cf(vR -- VL) in Sector I. The solution V(~) is then found to be
v(~ ) = VL + ~
VR -- VL coshk
OR - - OL O(~) = OL q- 8 ~ sinhk
ek
e-(Z+Cf)2/2e dz + e-k
e -(z-cf)2/2e dz
ek
e_(Z+Cf)2/2e d z _ e_ k
e_(Z_Cf)2/2e d z
--O~
,
,
(A. 15a)
(A.15b)
--OO
where k is obtained through Eq. (A.14). Even though sinh k appears in the denominator if (A.15b), when k = 0, we must have oR = OL by (A.14), so that this expression is not singular. We claim that as E ~ 0, the integrands in (A.15) converge to multiples of &functions, centered at z = q-cf, and that the integrals converge to piecewise constant functions with jumps at x/t = ~ = -q-cf; that is, as E --+ 0, the solution (A.15) converges to the distinguished solution of the Riemann problem of Section 2. A proof of this claim - based on the sharp peaking of the Gaussian integrands, and analogous to the convergence of the heat kernel - is given in [11]; here we merely illustrate the essential issues graphically. We plot the two components of V(~) from (A.15) in Figs. l l ( a ) and (b), for different values of E. As expected, the smooth profile steepens to resemble the shock solutions of Section 2, as ¢ is decreased. In both figures, we took cf = 2, with UL = (--5, --2) and UR = (0, 0). Note that Cs does not appear in formulas (A.15), which apply to Sector I. Next, in Fig. 12(a), we display Sector I of the (v, o)-plane, 10 corresponding to the above initial data, bounded by dashed lines. For UL = (--5, --2) and UR = (0, 0), the intermediate state is U~ = ( - 2 , 4). The solution of (A.5) and (A.6) with E = 0.5, the solid curve, is practically identical to those in Section 2, except for some slight curvature near UI. (As may be seen by comparison with Fig. 11, motion along the trajectory is very far, indeed, from uniform.) T h e s m o o t h i n g can b e s e e n m o r e clearly in Fig. 12(b), w h i c h is a m a g n i f i c a t i o n o f the p r e v i o u s figure around
U = UI. Note that for ~ = 0.1, the trajectory is essentially a corner in the phase plane, while for larger E, the trajectory d o e s n o t quite reach oi. 10 In the notation of Section 2, this is called sector 3; see footnote 8.
B.T. Hayes, D.G. Schaeffer/Physica D 121 (1998) 193-212
212 Phase Portrait
of Traveling Wave Solution
Detail of Phase Plane near the Intermediate State 4 3.95 I
8.9
~ 4 __
2
d ii
O
3.85 3.8
~ /
~-4
3.71-
3,65 ~ /
-6
-8 -1(
3.75
3.6 3.55
-5
-4 -2 -1 velocity; U _ L = ( - 2 , - 5 ) , U_I=(-2, 4) and U_R=(O,O)
(a)
35~ -2.25
-2'2
-2'15 velocity:
i p i i ' -2 -1.95 -1.~ - t . 3 5 - 1 . s -~.05 epsilon=0,1 (solid), 0.5 (dash), 0,75 (dash-dot)
-2'.1
-1.~5
(b)
Fig. 12. (a) Phase portrait in the (v, cr)-plane, for a viscous profile with e = 0.5. (b) Detail of phase portrait near V(~) = UI, for E = 0.1, 0.5 and 0.75.
Acknowledgements W e are grateful to G. G u d e h u s and V. O s i n o v for stimulating c o n v e r s a t i o n s during w h i c h this p r o b l e m w a s isolated. W e are further i n d e b t e d to V. O s i n o v for c a l l i n g attention to the u n i q u e n e s s i s s u e in s o l v i n g (1) and (13).
References [1] A. Bressan, The unique limit of the Glimm scheme, Arch. Rational Mech. Anal. 130 (1995) 205-230. [2] C.M. Dafermos, Solution of the Riemann problem for a class of hyperbolic systems of conservation laws by the viscosity method, Arch. Rational Mech. Anal. 52 (1973) 1-9. [3] P.G. Drazin, Nonlinear Systems, Cambridge University Press, Cambridge, 1992. [4] J. Glimm, Solutions in the large for nonlinear hyperbolic systems of equations, Comm. Pure Appl. Math. 18 (1965) 95-105. [5] M.S. Gordon, M. Shearer, D.G. Schaeffer, Plane shear waves in a fully saturated granular medium with velocity and stress controlled boundary conditions, Int. J. Non-Linear Mech. 32 (1997) 489-503. [6] D. Jacobs, W.R. McKinney, M. Shearer, Traveling wave solutions of the modified Korteweg-de Vries Burgers equation, J. Diff. Eqs. 16 (1995) 448-467. [7] A. Niemunis, I. Herle, Hypoplastic model for cohesionless soils with elastic strain range, Mech. Cohesive-frictional Mater. 2 (1997) 279-299. [8] V.A. Osinov, Plane waves and dynamic ill-posedness in granular media, in: R.P. Berhinger, J.T. Jenkins (Eds.), Powders and Grains 97, Balkema Press, 1997. [9] V.A. Osinov, Theoretical investigation of large-amplitude waves in granular soils, Soil Dyn. Earthquake Eng. 17 (1998) 13-28. [10] V.A. Osinov, G. Gudehus, Plane shear waves and loss of stability in a saturated granular body, Mech. of Cohesive-Frictional Mater. 1 (1996) 25-44. [11] M. Shearer, D.G. Schaeffer, Fully nonlinear hyperbolic systems of partial differential equations related to hypoplasticity, Comm. PDE 20 (1995) 1133-1153. [12] M. Slemrod, Admissibility criteria for propagating phase boundaries in a van der Waals fluid, Arch. Rational Mech. Anal. 81 (1983) 301-315. [13] J. Smoller, Shock Waves and Reaction-Diffusion Equations, Springer, Berlin, 1983.