Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://mailybaev.imec.msu.ru/papers/MarchesinMailybaev2005.pdf Äàòà èçìåíåíèÿ: Mon Aug 22 15:20:42 2005 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 19:45:25 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: space station

Dual-family viscous shock waves in n conservation laws with application to multi-phase flow in porous media
Dan Marchesin and Alexei A. Mailybaev August 22, 2005


Abstract We consider shock waves satisfying the viscous profile criterion in general systems of n conservation laws. We study Si,j dual-family shock waves, which are associated with a pair of characteristic families i and j . We explicitly introduce defining equations relating states and speeds of Si,j shocks, which include the Rankine­Hugoniot conditions and additional equations resulting from the viscous profile requirement. Then we develop a constructive method for finding the general local solution of the defining equations for such shocks and derive formulae for the sensitivity analysis of Si,j shocks under change of problem parameters. All possible structures of solutions of the Riemann problems containing Si,j shocks and classical waves are described. As a physical application, all types of Si,j shocks with i > j are detected and studied in a family of models for multi-phase flow in porous media.

Keywords: Dual-family shock, viscous profile, conservation laws, sensitivity analysis, Riemann problem, multi-phase flow, porous medium

1

Intro duction

In this paper, shock waves in general systems of n conservation laws in one space dimension x are considered. When shock waves are required to possess viscous profiles rather than to satisfy Lax's inequalities, new types of shocks arise. In general, these
This work was supported in part by: CNPq under Grant 301532/2003-6, FINEP under CTPETRO Grant 21.01.0248.00, IM-AGIMB/ IMPA, CAPES under Grant 0722/2003 (PAEP no. 0143/03-00), and President of RF grant No. MK-3317.2004.1 Instituto Nacional de Matem´ atica Pura e Aplicada ­ IMPA, Estrada Dona Castorina, 110, 22460-320 Rio de Janeiro RJ, Brazil. E-mail: marchesi@fluid.impa.br Institute of Mechanics, Moscow State Lomonosov University, Michurinsky pr. 1, 119192 Moscow, Russia. E-mail: mailybaev@imec.msu.ru. Current address : Instituto Nacional de Matem´ atica Pura e Aplicada ­ IMPA, Estrada Dona Castorina, 110, 22460-320 Rio de Janeiro RJ, Brazil.


1


shocks may be associated with the i-th characteristic family on the left and the j -th characteristic family on the right. We call such waves Si,j dual-family shocks. It was shown in [5] (see also [6]) that for i > j the viscous profile requirement provides exactly the number of additional equations (i - j equations) that is necessary to ensure that the number of characteristics emanating from the shock in positive time direction equals the number of independent equations at the shock interface. For systems of two equations, transitional shock waves (i = j + 1) were studied in [4, 12, 14], and novel structures of Riemann solutions resulting from such shocks were described. Shock waves with one or several additional equations for the viscous profile were found in problems of wave propagation in ferromagnetics, composite elastic media, elastic beams, and MHD, see [6], and in three phase flow in porous media they were analyzed for the case S2,1 in [11]. A program for studying stability in the sense of Hadamard of Si,j shocks was presented in [7], where a simple example of S3,1 shock was exhibited; another example of S3,1 shock can be found in [8]. In this paper, we explicitly introduce defining equations that relate states and speeds of Si,j shocks. These equations include Rankine­Hugoniot conditions (basic equations) and additional equations resulting from the viscous profile requirement. Then we develop a constructive method for perturbation analysis of general dualfamily shocks under parameter change, in which relationships between states at opposite sides of the shock and shock speed resulting from perturbations of problem parameters are derived. The role of Si,j shocks in solutions of the Riemann problem is described. It turns out that Si,j shocks with i > j may appear in generic Riemann solutions. The presence of Si,j shocks with i > j + 1 leads to the repetition of separated classical waves of the same characteristic family in a Riemann solution. As a physical application, we consider a flow of multi-phase fluid through porous medium with the quadratic Corey model for relative permeabilities of fluid phases. For the identity viscosity matrix, we find analytically Si,j shocks for any i > j . The method for finding these shocks by continuation for state-dependent viscosity matrices of Corey type models is discussed. The variety and form of dual-family shocks may be important for applications such as oil and gas recovery. The paper is organized as follows. Dual-family shocks are introduced in Section 2. Section 3 contains qualitative study of viscous profile requirement and defining equations. The defining equations relating states and speeds of Si,j shocks are explicitly given in Section 4, and variation of these equations under change of problem parameters is studied. Section 5 discusses structures of Riemann solutions containing Si,j shocks. In Section 6, Si,j shocks are analyzed in multi-phase flow through porous medium. The conclusion summarizes the contribution. Some technical proofs are given in the Appendix.

2

Dual-family sho ck waves
G(U ) F (U ) + = t x x D (U ) 2 U x , t 0, x R (2.1)

We consider systems of partial differential equations of the form


in the vanishing viscosity limit 0. The function representing conserved quantities G(U ) Rn , the flux function F (U ) Rn , and the n â n viscosity matrix D(U ) depend smoothly on the state vector U Rn . Taking = 0 in (2.1) yields a system of n first-order conservation laws G(U ) F (U ) + = 0. t x (2.2)

Real eigenvalues (U ) of the characteristic equation det( F / U - G/ U ) = 0 are the characteristic speeds. Assuming that all the eigenvalues are real and distinct in a region of state space U (the strictly hyperbolic region), we list them in increasing order 1 < 2 < · · · < n . A shock wave is a discontinuous (weak) solution of system (2.2) consisting of a left state U- = limx/t s U (x, t) and a right state U+ = limx/t s U (x, t), where s is the shock speed. A shock wave is considered admissible if there is a traveling wave solution (or viscous profile) U (x, t) = U ( ), = (x - st)/ of system (2.1), which represents the shock in the vanishing viscosity limit 0. Substituting this solution into (2.1) and integrating over , we find that U ( ) is a solution (orbit) of the system of ordinary differential equations D(U )U = F (U ) - F (U- ) - s(G(U ) - G(U- )), (2.3)

"connecting" the left equilibrium U (-) = U- to the right equilibrium U (+) = U+ ; the dot denotes the derivative with respect to . By linearizing equation (2.3) about the equilibria U- and U+ we obtain U = B (U± , s)U, where B (U, s) is the n â n matrix B (U, s) = D-1 (U ) F (U ) - F (U- ) - s(G(U ) - G(U- )) . U (2.5) U ( ) = U ( ) - U± , (2.4)

Let µi (U, s), i = 1, . . . , n be the eigenvalues of the matrix B (U, s) ordered with increasing real parts Re µ1 Re µ2 · · · Re µn . Let us define an Si,j shock as a shock possessing a viscous profile and satisfying the inequalities Re µi-1 (U- , s) < 0 < Re µi (U- , s), Si,j : (2.6) Re µj (U+ , s) < 0 < Re µj +1 (U+ , s) (if i - 1 = 0 or j + 1 = n + 1, the corresponding inequality is disregarded). It is easy to see that µi (U- , s) = 0 and µj (U+ , s) = 0 if s = i (U- ) and s = j (U+ ), respectively. Under rather general conditions (see e.g. [6, 10]), inequalities (2.6) reduce to Si,j :
i-1

(U- ) < s < i (U- ),
j +1

j (U+ ) < s <

(U+ ).

(2.7)

For i = j , inequalities (2.7) are the Lax conditions. Thus, an Si,i shock is a classical i-shock. Shocks with i < j are called overcompressive. For i = j + 1 such a 3


shock is termed transitional or undercompressive. The inequalities in the first row of (2.7) coincide with the Lax conditions for the left state U- of an i-shock. Analogously, the inequalities in the second row of (2.7) coincide with the Lax conditions for the right state U+ of a j -shock. Therefore, the Si,j shock can be seen as a dual-family shock wave associated with the i-th characteristic family on the left and with the j -th characteristic family on the right. For characteristic shocks, the shock speed coincides with a characteristic speed at left, right, or both sides. As it was noticed above, the eigenvalues µi (U- , s) or/and µj (U+ , s) vanish in these cases. Thus, we can distinguish three types of characteristic shocks as S S S
- i,j + i,j ± i,j

: µi (U- , s) = 0, Re µj (U+ , s) < 0 < Re µ : Re µ
i-1

j +1

(U+ , s); (2.8)

(U- , s) < 0 < Re µi (U- , s), µj (U+ , s) = 0;

: µi (U- , s) = µj (U+ , s) = 0.

Conditions (2.8) can be written in terms of characteristic speeds as S S S
- i,j + i,j ± i,j

: s = i (U- ), j (U+ ) < s < :
i-1

j +1

(U+ ); (2.9)

(U- ) < s < i (U- ), s = j (U+ );

: s = i (U- ) = j (U+ ).

3

Defining equations

Let us consider Si,j as a point in the space (U- , U+ , s) of the left and right states and speed of shocks. For each i, j , the set of all Si,j shocks can be expected to define a smooth surface Si,j in the space (U- , U+ , s); see [3] for an example of such a surface for transitional shocks and quadratic flow functions. Locally this surface can be given by a maximal rank system of equations. Following [6, 12], we distinguish the basic equations defined by the system of conservation laws (2.2) and additional equations determined by the viscous profile requirement. There are n basic equations, which are the Rankine­Hugoniot conditions H(U- , U+ , s) F (U+ ) - F (U- ) - s(G(U+ ) - G(U- )) = 0 Rn . (3.1)

These equations are obtained from the condition that the shock is a weak solution of (2.2); they also follow from requiring that U+ is an equilibrium of (2.3). Additional nadd equations are denoted by H
add

(U- , U+ , s) = 0 Rn

add

.

(3.2)

The number of additional equations nadd can be determined by considering the intersection of the unstable manifold Mu (U- ) and the stable manifold Ms (U+ ) of the equilibria U± . The viscous profile exists if the intersection of these manifolds is not empty, see Figure 3.1(a). Using inequalities (2.6), we find dim Mu (U- ) = n - i + 1, 4 dim Ms (U+ ) = j. (3.3)


(a) M

U (z) M U-

s

u

U+

(b) U (z)

(c)

M M M
u 0u

s

U-

M

0s

U+

Figure 3.1: A viscous profile: (a) for Si,j shock, (b) and (c) for S

± i,j

shock.

Generically, the manifolds Mu (U- ) and Ms (U+ ) intersect forming a locally structurally stable phase state configuration if dim Mu (U- ) + dim Ms (U+ ) > n, which implies i j . In this case nadd = 0, i.e., there are no additional equations. Since the intersection of the manifolds has dimension dim Mu (U- ) + dim Ms (U+ ) - n = i - j + 1, generically there exists a single viscous profile for classical shocks (i = j ) and an infinite number of viscous profiles for overcompressive shocks (i < j ). In case i > j , the manifolds Mu (U- ) and Ms (U+ ) do not intersect in general. More precisely, if Mu (U- ) and Ms (U+ ) intersect and the intersection is a single orbit, then this is a singular situation (the so-called nontransversal intersection) of codimension i - j [2]. The least degenerate case in this situation occurs when the tangent spaces of Mu (U- ) and Ms (U+ ) have one-dimensional intersection at each point of the connecting orbit (so-called quasi-transverse intersection). Therefore, the number of additional equations (3.2) for Si,j shocks with i > j equals nadd = i - j . - + ± Similar surfaces Si,j , Si,j , and Si,j , under some genericity assumptions for the nature of nonhyperbolic equilibria, can be defined for characteristic dual-family shocks. In these cases, there are one or two equations resulting from the conditions s = (U- ) or/and s = (U+ ). These equations are related to system (2.2) only, so that they supplement the basic equations (3.1). In case s = i (U- ), the orbits of equation (2.3) starting at U- as - generally form a smooth manifold M0u (U- ) of dimension n - i + 1 [16]. This manifold has a boundary, which is the unstable manifold Mu (U- ) 5


of dimension n - i, see Figure 3.1(b). We consider the generic situation when the viscous profile U ( ) does not lie in the boundary. Similarly, in case s = j (U+ ), orbits finishing at U+ as + form a smooth manifold M0s (U+ ) of dimension j . This manifold has a boundary, which is the stable manifold Ms (U+ ) of dimension j - 1, see Figure 3.1(c). We assume that the viscous profile U ( ) does not lie in Ms (U+ ). The dimensions of M0u (U- ) and M0s (U+ ) coincide with the dimensions (3.3) for Si,j shocks. Hence, the number of additional equations remains the same, nadd = i - j . The total number of defining equations provides codimensions of the shock surfaces: ij:
- + codim Si,j = n, codim Si,j = codim Si,j = n + 1, ± codim Si,j = n + 2; - codim Si,j = n + i - j, codim Si,j = codim S ± codim Si,j = n + i - j + 2. + i,j

(3.4)

i>j:

= n + i - j + 1,

(3.5)

We wee that the viscous profile admissibility criterion for Si,j shocks plays qualitatively different role in cases i j and i > j . Indeed, for classical (i = j ) and overcompressive (i < j ) shocks, generically the viscous profile is a structurally stable connecting orbit that persists under small perturbations of shock states and speed satisfying the Rankine-Hugoniot conditions (3.1). However, for Si,j shocks with i > j , the existence of viscous profile implies additional relations between states and shock speeds. These relations depend on the form of viscous terms governed by the viscosity matrix D(U ). Notice that surfaces corresponding to characteristic shocks form a part of the Si,j surface boundary. The remaining part of the Si,j boundary is related to bifurcations of the viscous profile, see [13].
Example. Let us consider an Si,j shock with states U- , U+ and speed s in a system of three viscous conservation laws. There are n = 3 defining equations if i j (classical and overcompressive shocks), which are the Rankine­Hugoniot conditions (3.1). These equations can be solved for U+ . As a result, there is a dual-family shock of the same type with an arbitrary left state U- and speed s taken in neighborhoods of U- and s ; the uniquely determined right state U+ is close to U+ , see Figure 3.2(a). If i = j + 1 (transitional shocks S2,1 and S3,2 ), there is one additional equation, which generically can be solved for the speed s. Therefore, for each left state U- in a neighborhood of U- there is a dual-family shock of the same type with the uniquely determined right state U+ and speed s, see Figure 3.2(b). Finally, there are two additional equations for an S3,1 shock. Solving these equations for s, the remaining equation defines a surface S- in the space U- passing through U- . Hence, S3,1 shocks exist only for left states U- chosen on this surface. The speed and right state are given uniquely by the choice of a point U- S- , see Figure 3.2(c).

6


(a)

(b)

(c)

U+* U+ s* s UU-* UU-*

U+* U+ s s* s s* S UU-*
-

U+ U+*

U

Figure 3.2: Viscous profiles of Si,j shocks for three conservation laws: (a) classical and overcompressive shocks, i j , (b) S2,1 and S3,2 shocks (transitional shocks), (c) S3,1 shocks.

4

Sensitivity analysis of S

i,j

sho cks

In this section, we study the local form and perturbation of the shock surface Si,j . In case i j , this information can be obtained directly by variation of the explicit Rankine­Hugoniot conditions (3.1). Thus, we will focus on the case i > j when the viscous profile condition results in additional equations (3.2). Let us determine explicitly the system of additional equations. Consider a viscous profile U ( ) of an Si,j shock (i > j ) with states U± and a speed s. Linearizing equation (2.3) near the solution U ( ) yields V = B (U ( ), s)V , V (-) = V (+) = 0, V ( ) Rn , (4.1)

where the matrix B (U, s) is given in (2.5). The corresponding adjoint linear system takes the form W = -B T (U ( ), s)W, W (-) = W (+) = 0, W ( ) Rn . (4.2)

Let us denote by W the linear space of solutions W ( ) of system (4.2). For any + bounded function X ( ) Rn , solutions W ( ) W have the property - W T (X - B (U ( ), s)X )d = 0. The proof of the following proposition is given in the Appendix. Prop osition 1. For any real and function W ( ) W , the vector W ( ) is orthogonal to both Mu (U- ) and Ms (U+ ) at the point U ( ) of the connecting orbit. In case of quasi-transverse intersection of Mu (U- ) and Ms (U+ ), we have dim W = i - j . Consider a Poincar´ hyperplane P orthogonal to the viscous profile at a fixed point e U P = U ( P ), see Figure 4.1. Under perturbations U- + U- , U+ + U+ , and s + s satisfying the Rankine-Hugoniot conditions (3.1), the stable and unstable manifolds 7


W1
U( z)

W2

U+P U-P
M
u

M

s

U+ U+DU+ +

P

UU-+DU-

Figure 4.1: Perturbation of stable and unstable manifolds near a viscous profile. Mu (U- + U- ) and Ms (U+ + U+ ) change. Assuming that the perturbations are P P small, we can uniquely define the points U- Mu (U- + U- ) and U+ Ms (U+ + P P P U+ ) such that U± P and the vector U+ - U- is orthogonal to both Mu (U- ) P P and Ms (U+ ) at U P . The vector U+ - U- is a measure of the separation between the perturbed manifolds Mu (U- + U- ) and Ms (U+ + U+ ). One can see that the manifolds Mu (U- + U- ) and Ms (U+ + U+ ) intersect, i.e., there is a connecting P P orbit between U- + U- and U+ + U+ if and only if U+ = U- , see Figure 4.1. Let us choose a basis W1 ( ), . . . , Wi-j ( ) of W . Then, we can introduce the function Hadd (U- , U+ , s) with values in Ri-j as H
add

(U- , U+ , s) =

T P P W1 ( P )(U+ - U- ), . . . , W

T i-j

P P ( P )(U+ - U- )

P P ^ = W T ( P )(U+ - U- ),

(4.3)

where

^ W ( ) = [W1 ( ), . . . , W

i-j

( )]

(4.4)

is an n â (i - j ) matrix. Indeed, since the vectors W1 ( P ), . . . , Wi-j ( P ) are linearly independent and orthogonal to both Mu (U- ) and Ms (U+ ) at U P (see Proposition 1), P P the condition U+ = U- is equivalent to Hadd (U- , U+ , s) = 0. The local form of the surface Si,j found by linearizing the defining equations near a point (U- , U+ , s) Si,j is described as follows (see the proof in the Appendix). Theorem 1. The tangent plane (U- , U+ , s) of the manifold Si,j at the point (U- , U+ , s) Si,j is given by the equations H = 0, Hadd = 0, where H = G F -s U U H
add

U+ -
U =U
+

G F -s U U
-1 U

U- - (G+ - G- )s, (4.5)
U =U
-

+

=
-

^ WTD

F G -s U U

d
U =U
-

U

-

+

(4.6)

+
-

- ^ W T DU 1 (GU - G- )d s

8


are the linear parts of the functions H and Hadd evaluated terms of deviations (U- , U+ , s). Here we introduced D(U ( )), FU ( ) = F (U ( )), GU ( ) = G(U ( )), F± = where U ( ) is the viscous profile of the Si,j shock at the

at (U- , U+ , s) and written in the short notations DU ( ) = F (U± ), and G± = G(U± ), initial point (U- , U+ , s).

It is remarkable that expression (4.6) does not depend on the choice of the point U P = U ( P ) (the position of the plane P ) in the definition of the function Hadd (U- , U+ , s), just as in the case of Melnikov integrals for systems of planar differential equations [2]. As a model of a physical system, equation (2.1) typically depends on one or more problem parameters. Under variations of these parameters, the functions G(U ), F (U ), and D(U ) undergo perturbations G(U ), F (U ), and D(U ). If these perturbations are small, the manifold Si,j undergoes a small perturbation. The first order approximation of the perturbed manifold can be determined as follows (see the proof in the Appendix). Theorem 2. Let (U- , U+ , s) Si,j and consider perturbations G(U ), F (U ), D(U ) of the system functions. Then the first order approximation of the perturbed manifold Si,j near the point (U- , U+ , s) is given by the equations H = -F+ + F- + s(G+ - G- ), H
add +

(4.7)

=-
- +

- - ^ W T DU 1 DU DU 1 FU - F- - s(GU - G- ) d

(4.8)
-1 U

+
-

^ WTD

FU - F- - s(GU - G- ) d ,

where H and H

add

are given by expressions (4.5) and (4.6).

Theorems 1 and 2 determine all nearby Si,j shock waves, even when problem parameters are changed, using the information on a particular shock and its viscous profile. This method is useful for constructing solutions of conservation laws possessing Si,j shocks, continuation procedures, and parametric analysis. - + ± The characteristic shock waves Si,j , Si,j , and Si,j are studied in the same way. In addition to equations (3.1) and (3.2), (4.3), one should impose conditions ensuring that the shock speed is equal to the corresponding characteristic speed at one or both sides of the shock. This adds one or two equations in Theorems 1 and 2 found by linearizing the equations s = (U- ) and s = (U+ ).

5

Dual-family sho cks in Riemann solutions
the Rieat x = 0: found in changing

The basic initial-value problem for a system of conservation laws (2.2) is mann problem, given by piecewise constant initial data with a single jump U (x, 0) = UL for x < 0 and U (x, 0) = UR for x > 0. The solution is ^ scale-invariant form U (x, t) = U ( ), = x/t, consisting of continuously

9


waves (rarefaction waves), jump discontinuities (shock waves), and separating constant states. Classically, there are n families of rarefaction waves, one for each characteristic speed, which we denote by Ri , i = 1, . . . , n. We require all the shock waves - + ± to have a viscous profile, i.e., there can be Si,j , Si,j , Si,j , and Si,j shocks. The structure of a Riemann solution is given by a sequence of waves wk w1 , w 2 , . . . , w m , (5.1)
- + ± appearing with increasing value of . Here each wave wk {Ri , Si,j , Si,j , Si,j , Si,j } is a rarefaction or shock. The wave wk has left and right states U(k)- and U(k)+ and speeds (k)- < (k)+ for a rarefaction wave and s(k) = (k)- = (k)+ for a shock wave. The left state of the first wave w1 and the right state of the last wave wm are the initial conditions of Riemann problem: U(1)- = UL and U(m)+ = UR . The natural requirements in sequence (5.1) are

U

(k)+

=U

(k+1)-

,



(k)+



(k+1)-

.

(5.2)

Relations (5.2) imply that the right state of the wave wk coincides with the left state of the wave wk+1 , and the right speed of wk is lower or equal to the left speed of wk+1 . If (k)+ < (k+1)- then there is a separating constant state between wk and wk+1 . In this case we will use the notation wk -- wk+1 . If (k)+ = (k+1)- then the waves do not possess a separating constant state. This situation will be denoted by wk | wk+1 . A classical Riemann solution consists of n classical wave groups separated by constant states; each classical wave group consists of adjoining rarefactions Ri and classical shocks Si,i (simply Si ) of the same family. The classical structure R1 -- R2 -- S3 of a Riemann solution in a system of three conservation laws is shown in Figure 5.1(a) using characteristic lines in the space-time plane (shock waves are presented by bold lines and rarefaction waves are given by thin line fans). Conditions (5.2) imply that each pair of subsequent waves wk , wk+1 in a general Riemann solution has one of the following types
{Rj or Si,j } -- {Ri or Si ,j },

j
Ri | {S {S
+ i,j

- i,j

± or Si,j },

± or Si,j } | Rj ,

- + ± where Si,j stands for Si,j , Si,j , Si,j , or Si,j . The most important structures of a Riemann solution are the generic ones: they do not change under perturbations of initial conditions UL , UR , flux function, and viscosity matrix. Generic structures are "full" in the sense that no wave can be added to a sequence without violating conditions (5.3). For example, the sequence R1 -- S3 is not generic since it can be extended to R1 -- S2 -- S3 . Only the shocks with i j may appear in generic structures. Generically, overcompressive shocks (i < j ) bifurcate to a set of waves under perturbations with arbitrarily small amplitudes. Let us assume that all the states of a Riemann solution belong to the region of strict hyperbolicity, where all the characteristic speeds are distinct; viscous profiles of shock waves are not restricted to this region, i.e., they may cross elliptic regions in state space. Then, we can describe generic structures of a Riemann solution (for the case of two conservation laws the following theorem was proved in [12]; for the complete proof see [9]).

10


(a)

(b)
R
2

(c)
R S
3

t
R
1

S2 t S R
3 1

S

3,2

2

t
S
3,1

R

R
3

1

S

2

R

3

x
0 0

x
0

x

Figure 5.1: Riemann solutions: (a) classical, (b) with S

- 3,2

shock, (c) with S

3,1

shock.

Theorem 3. Let (5.1) be the generic structure of a Riemann solution. Then w1 + - {R1 , S1 , S1 }, wm {Rn , Sn , Sn }, and each pair wk , wk+1 has one of the types {Rj or Si,j } -- {Ri or Si R i | {S {S
+ i,j - i,j ± or Si,j } | Rj , ,j

},

i = j + 1, i j , i j , i j, (5.4)

± or Si,j },

i j.
- 3,2

As an example, we list two generic nonclassical Riemann solution structures: R1 -- S2 -- R3 | S R1 -- S2 -- S
3,1

-- R3 ,

(5.5) (5.6)

-- S2 -- R3 .

A distinctive feature of structures (5.5) and (5.6) is that the classical waves R3 and S2 appears twice. Riemann solutions with these structures are shown in Figure 5.1(b,c). We see that Riemann solutions with dual-family shock waves violate the customary classical structure of sequences of n classical wave groups separated by constant states with increasing family number from left to right. The shocks with i > j + 1 introduce a "jump back" capability in this sequence allowing classical waves of (j + 1), . . . , (i - 1)-th characteristic families to appear repeatedly. Moreover, from the theoretical point of view, there is no general bound on the number of separated classical waves or of nonclassical shock waves in a Riemann solution for systems of n > 2 conservation laws. The existence of several separated waves corresponding to the same characteristic family is a property of Riemann solutions that was observed only in [8]. We remark that the generic structures of Riemann solutions described above do not include transitional rarefaction waves, see [4]. These waves are related to characteristic speeds given by multiple eigenvalues. Therefore, transitional rarefactions do not appear in strictly hyperbolic systems.

6

Dual-family sho cks in multi-phase flows in p orous media

Let us consider one-dimensional horizontal flow of n + 1 immiscible fluid phases in a porous medium. The fluids can be, for instance, a mixture of gas or CO2 , water, 11


light oil, and heavy (viscous) oil. We assume that the whole pore space is occupied by the fluids; compressibility, thermal and gravitational effects are neglected. The equations expressing conservation of mass of the i-th phase based on Darcy's law of force is (see e.g. [1]) (si ) + (v fi ) = t x x Kl
i j =i

f

j

pij x

,

i = 1 , . . . , n + 1,

(6.1)

where the constants and K denote the porosity and absolute permeability of the porous medium, and v is the total seepage velocity of the fluid. For phase i, si is the saturation, fi is the fractional flow function, and li is the relative mobility, which can be chosen as the quadratic Corey model (see [1]): li fi = , l s2 li = i , µi l = l1 + · · · + ln+1 , (6.2)

where µi is the viscosity of phase i; for simplicity, all irreducible phase saturations were set to zero. The capillarity pressures pij = pi - pj between the phases i and j are measured experimentally as functions of saturations; here pi and pj are the pressures in phases i and j . Since the fluids occupy the whole available space, the saturations satisfy s1 + · · · + sn + sn
+1

= 1.

(6.3)

Similarly, f1 + · · · + fn + fn+1 = 1. As a consequence, any n saturations describe the state of the fluid. Hence, any of the n + 1 equations in system (6.1) is redundant, and the latter can be reduced to an n equation system in n saturations. As the total seepage velocity v is given by boundary conditions, we assume that it ~ is a positive constant and we set t = (L/v )t and x = Lx, where L is the characteristic ~ length of the system. Dividing both sides by v /L, this change of variables removes v and from the left-hand side of system (6.1). For simplicity of notation, we drop the tildes below. Finally, choosing any n saturations as state variables (e.g. Ui = si , i = 1, . . . , n), we arrive at the dimensionless system U F (U ) + = t x x D(U ) U x , U = (U1 , . . . , Un )T . (6.4)

The components of the vector F (U ) = (F1 , . . . , Fn )T represent the fractional flow functions Fi (U ) = where U
n+1

li , l

li (Ui ) = Ui2 /µi , l(U ) = l1 + · · · + ln + ln+1 , = 1 - U1 - · · · - U

(6.5)

n

(6.6)

is the saturation of remaining (n + 1)-th phase, = K/L2 , and D(U ) is the dimensionless viscosity matrix with components dij (U ) = (L/v )li k=i fk pik . The xj 12


1

U3

II I U
0
umb

U2
1

(U4 =1) U1
1

Figure 6.1: State space for four-phase flows. The system restricted to the shaded planes is reduced to a lower dimensional one. dimensionless parameter is small for large characteristic lengths L, so that many aspects of the asymptotic behaviour of solutions may be described by the system of first order conservation laws obtained by taking = 0. For simplicity, we will use dimensionless viscosities µi = µi /(µ1 + · · · + µn + µn+1 ) in (6.5), so that, dropping ~ the tildes, µ1 + · · · + µn + µn+1 = 1. (6.7) The physical domain for the state vector U is given by the inequalities U1 + · · · + Un 1, Ui 0, i = 1, . . . , n. (6.8)

This domain represents a simplex in state space with n + 1 vertices; each vertex corresponds to single phase fluid with Ui = 1 (i = 1, . . . , n + 1); see Figure 6.1. From (6.5), the n â n Jacobian matrix F / U can be expressed in the form F U = 1 dl1 dln diag ,..., l dU1 dUn 1 - 2 (l1 , . . . , ln )T l dln+1 dln dln+1 dl1 - ,..., - dU1 dUn+1 dUn dUn+1 .

(6.9)

One can check that all the eigenvalues of matrix (6.9) are strictly positive inside the physical domain (6.8), i.e., flows with positive speed v have only positive characteristic speeds. There is a so-called umbilic point U umb , at which dl1 dln dln+1 = ··· = = . dU1 dUn dUn+1 (6.10)

This is the resonance state, where F / U is a multiple of the identity matrix and, dl thus, all the characteristic speeds (eigenvalues) merge to the same value 1 dUii . There l 13


are n + 1 umbilic points at the boundary of the physical domain, which are the vertices of the simplex (6.8). At each vertex, all the characteristic speeds merge to the value zero. One can check that for the flux functions (6.5) there is a unique umbilic point U umb = (µ1 , . . . , µn )T , Un+1 = µn+1 in the interior of physical domain (6.8), at which all the characteristic speeds equal 2. The umbilic point U umb exists and is unique for the more general case when the relative mobilities li (Ui ) are arbitrary functions of the corresponding phase saturations Ui with positive first and second derivatives such that l1 (0) = . . . = ln+1 (0) = 0 and dl1 /dU1 (0) = . . . = dln+1 /dUn+1 (0) = 0. Indeed, the vector (l1 , . . . , ln ) has positive components inside the physical domain. Hence, the matrix (6.9) is a multiple of the identity matrix if and only if the vector dl dl dl1 dl - dUn+1 , . . . , dUn - dUn+1 = 0, which yields equations (6.10). Since dli /dUi are dU1 n n+1 n+1 strictly increasing functions of Ui vanishing at Ui = 0, the equations dl1 /dU1 = · · · = dln /dUn = have a unique solution for > 0, and U1 (), . . . , Un () are increasing functions of . The function dln+1 /dUn+1 , where Un+1 () = 1 - U1 () - · · · - Un (), is a decreasing function of . Additionally, we have dln+1 /dUn+1 > 0 at = 0 (Un+1 = 1), and dln+1 /dUn+1 = 0 at = > 0, where the value of is given by the equation Un+1 () = 0. Hence, there exists a unique umb solving the equation dln+1 /dUn+1 = such that 0 < umb < . This yields a unique umbilic point U umb = (U1 (umb ), . . . , Un (umb ))T , which lies inside the physical region. In order to study dual-family shock waves, we artificially take the identity viscosity matrix D(U ) I . Our main motivation here is to show analytically the existence of all types of Si,j shock waves with -n < i - j < n in the system. This existence provides the evidence that any Si,j shock may appear in the system with a realistic viscosity matrix D(U ).

6.1

Reduced dimension systems

By setting one of the Ui = 0 in equation (6.4), we obtain a reduced system describing n instead of n + 1 phase flows. This system "lives" on the face of the simplex; see Figure 6.1, where the lightly shaded face I corresponds to U4 = 0. Of course, all the results of this section hold for the reduced system. This reduction can be done iteratively until we reach a scalar partial differential equation describing two-phase flow. This system is restricted to one of the edges of the simplex. Notice that this reduction is valid for any physical viscosity matrix D(U ). There are other subsystems of (6.4) with dimension n - 1. Let us consider the simplex of codimension 1: U= µ2 µ1 , , U3 , . . . , U µ1 + µ2 µ1 + µ2
T n

; + U3 + · · · + Un 0; , U3 , . . . , Un 0.
umb

(6.11) For n = 3 this is the shaded plane II containing U shown in Figure 6.1. One can check that the first two equations of system (6.4), (6.5) restricted to plane (6.11) are equivalent (in this case we consider D(U ) I ), and the system possesses two equal

14


characteristic speeds =

2 . (µ1 + µ2 )l

(6.12)

As a result, system (6.4) is reduced to the n - 1 dimensional system with state vector (, U3 , . . . , Un )T and viscosities µ1 + µ2 , µ3 , . . . , µn . This system describes the flow of phases U3 , . . . , Un and the mixture of phases U1 , U2 in the proportion U1 /U2 = µ1 /µ2 that acts as a single phase. Analogous reductions can be done taking any two saturations Ui and Uj instead of U1 and U2 . Repeating such a reduction several times, we can obtain systems of lower dimensions. Each of these systems describes a multi-phase flow. Thus, all the results of this section hold for any of the reduced systems. Notice that the reduced systems "live" on lower dimensional simplices in the interior of physical state space (6.8). Using reduction (6.11) n - 1 times, each time with two state variables, we arrive at the system restricted to the line U = µ1 µn ,..., µ µ
T

,

µ = µ1 + · · · + µn .

(6.13)

Up to multiplication by the constant µi /µ, all n equations in system (6.4) are equivalent to the scalar partial differential equation on the line (6.13): F () 2 + = 2, t x x where 2 F () = , µl 2 - 2µ + µ l() = . µ(1 - µ) (6.14)

(6.15)

One can check that equation (6.14) coincides with the viscous profile equation taken for a scalar conservation law describing two-phase flow, where U = U1 = , U2 = 1 - , µ1 = µ, and µ2 = 1 - µ. According to (6.8), the physical interval for is 0 1. The line (6.13) contains the point U umb at = µ, as well as the vertex U = 0 at = 0. At points with 0 < < µ or µ < 1, there are n - 1 equal characteristic speeds () = 2 . µl (6.16)

There is an (n - 1)-dimensional eigenspace of the matrix F / U corresponding to the multiple eigenvalue (), which consists of the vectors r = (r1 , . . . , rn )T such that r1 + . . . + rn = 0. The remaining characteristic speed is equal to 2(1 - ) ~ () = , µ(1 - µ)l2 (6.17)

with corresponding eigenvector r = (µ1 , . . . , µn )T , which is parallel to the line (6.13). ~

15


6.2

S

n,1

dual-family sho ck waves

In this subsection, we seek for solutions of the viscous profile equation (2.3) that have the form (6.13). Here ( ) is a function of such that ± = (±) and U± = U± . Along the line (6.13), all n equations in system (2.3) are equivalent to the scalar ordinary differential equation = 2 2 - - - s( - - ), µl µl- l- = l(- ), (6.18)

which is the viscous profile equation for two-phase flow described by (6.14). Since the right-hand side of (6.18) vanishes at = + , we find the shock speed s= - + + - 2- + , µ(1 - µ)l- l+ l± = l(± ). (6.19)

Here expression (6.15) for l± was used. One can show that the ordinary differential equation (6.18) has a solution ( ) such that (±) = ± provided that both inequalities (+ - - )(2 + 2- + - 22 + - µ) > 0, - - (6.20) (+ - - )(2 + 2- + - 2- 2 - µ) > 0 + + hold. Conditions (6.20) are equivalent to the inequalities ~ ~ + < s < - , ~ ~ ± = (± ), (6.21)

~ where is the distinct characteristic speed (6.17). Under the conditions (2.7) and (6.21), we see that the shock wave having the viscous profile U is of type Sn,1 exactly if - < s < + , ± = (± ), (6.22)

where () is the characteristic speed (6.16) with multiplicity n - 1. With the use of (6.15) and (6.19), conditions (6.22) reduce to µ+ - µ- + 2- + (µ - + ) > 0, µ+ - µ- - 2- + (µ - - ) > 0. (6.23)

Figure 6.2 shows the region corresponding to the values (- , + ) satisfying inequalities (6.20) and (6.23) for 0 µ 1. The constants - and + taken in this region define left and right states and speeds of Sn,1 shocks by formulae (6.13) and (6.19). One can see that Sn,1 shock waves exist in the system for any µ. If µ 1/2, these shocks can be arbitrarily small, since the boundary of the region contains the point - = + = µ corresponding to the umbilic point. If µ < 1/2, the shock amplitude U+ - U- is bounded away from zero. In all the cases we have - < µ < + , which means that the viscous profile passes through the umbilic point. Since the state variables U1 , . . . , Un can be chosen as any n out of n + 1 saturations s1 , . . . , sn+1 , there are n other lines in state space similar to (6.13) containing the left and right states of Sn,1 shocks. 16


m

1 0.8 0.6 0.4 0.2 0 1 0.9 0.8 0.8 0.7 0.6 0.5 0 0.2 0.6 0.4 1
+

r

r

-

Figure 6.2: Values of (- , + ) corresponding to S

n,1

shocks for different µ.

- + ± Characteristic shocks of types Sn,1 , Sn,1 , and Sn,1 can be found, respectively, when the first, second, or both inequalities (6.20) become equalities. These equalities ~ ~ ~ ~ ~ correspond to s = - , s = + , and s = - = + , respectively, where ± are simple characteristic speeds at different sides of the shock. It turns out that there are no + ± Sn,1 and Sn,1 shocks with viscous profile (6.13). Parameters - and + corresponding - to Sn,1 shocks belong to the front face of the boundary shown darker in Figure 6.2. According to the results of Section 5, a total of n - 1 classical waves may be present at both sides of Sn,1 shock in a Riemann solution. Thus, there are Riemann solutions with Sn,1 shocks for any initial conditions UL and UR taken in certain neighborhoods of the left and right shock states U- and U+ . An example of a structure of Riemann solutions with S3,1 shocks was given in (5.6).

6.3

General Si,j dual-family sho cks
i,j

By reducing system dimension, we can find other types of S consider the case when U1 = · · · = Uk = 0 ,

shocks. First, let us (6.24)

i.e., only the phases Uk+1 , . . . , Un , Un+1 are present. Thus, (6.4) reduces to the system for (n - k + 1)-phase flow. Assume that we found an Si,j shock with speed s in the reduced system. In the full (n + 1)-phase system, this corresponds to a shock for which the first k components of the viscous profile are zero. Due to (6.24), there are k zero characteristic speeds at both sides of the shock. Thus, such a shock will be an Si+k,j +k shock for the full system. The shock speed s is always positive, which follows from the positivity of characteristic speeds. In the previous subsection, we found Sn,1 shocks in a general (n + 1)-phase system. Hence, by using the described recursive relation, we find an Sn,k+1 dual-family shock for any 0 k < n with the viscous profile U ( ) = µn ,..., 0, . . . , 0, µ µ µ
k+1 T

( ), 17

µ=µ

k+1

+ · · · + µn .

(6.25)


Here (µk+1 , . . . , µn )T ( )/µ is the viscous profile of the Sn-k,1 shock found for the reduced (n - k + 1)-phase system. These shocks lie on the boundary of the physical domain (6.8), yet they determine certain points (U- , U+ , s) on the shock surfaces Sn,k+1 . By using Theorem 1 in Section 4, we can find all nearby shock waves of the same type. Some of them belong to the interior of the physical domain (6.8) and, thus, they represent Sn,k+1 shocks intrinsic to (n + 1)-phase flow. The second type of dimension reduction corresponds to the plane (6.11). As shown above, in this plane the system reduces to an n - 1 dimensional system describing n phase flows. There are two equal characteristic speeds (6.12) at each point of the plane. Again, let us consider an Si,j shock with speed s in the reduced system. Then, if - < s < + for double characteristic speeds ± , we get an Si+1,j shock in the full (n + 1)-phase system. If + < s < - , this shock becomes an Si,j +1 shock in the full system. Finally, if s < ± or s > ± , we obtain Si,j and Si+1,j +1 shocks, respectively. This reduction can be used repeatedly. For example, Sn,1 shocks previously described are obtained by using n - 1 repeated reductions. The described approach allows finding dual-family shocks Si,j for any 0 < i-j < n. Note that by taking opposite signs in inequalities (6.21), (6.22) (equivalently, in (6.20), (6.23)), we obtain overcompressive shocks S1,n . With these shocks, by using the dimension reduction approach, we can locate overcompressive Si,j shocks for any -n < i - j < 0. The same method can be used for finding particular viscous profiles for these shocks (recall that overcompressive shocks possess an infinite number of viscous profiles).

6.4

Sensitivity analysis and continuation metho d

Let us apply Theorems 1 and 2 for an Sn,1 shock with the viscous profile (6.13). For this purpose, we fix some values of - and + , which determine uniquely U- , U+ , s, and the viscous profile U ( ) = µ1 , . . . , µn ( ). With the use of expression (6.9), µ µ the matrix B (U ( ), s) defined in (2.5) takes the form B (U ( ), s) = 2( ) 2(µ - ( )) -s I + µl(( )) µ(1 - µ) ( ) µl(( ))
2 T

(µ1 , . . . , µn )T (1, . . . , 1). (6.26)

Then, the general solution W ( ) of the adjoint linear system (4.2) is found as W ( ) = w1 wn ,..., µ1 µn
T

( ),

(6.27)

where w1 , . . . , wn are arbitrary constants satisfying the condition w1 + · · · + wn = 0, and the scalar function ( ) is determined by the equation = s- 2( ) µl(( ))
+

,
-

( )d = 1.

(6.28)

The latter equality in (6.28) is the normalization condition. The solution ( ) of (6.28) exist due to relations (6.16) and (6.22). Taking n - 1 linearly independent 18


^ solutions as columns of the n â (n - 1) matrix W ( ), 1/µ1 1/µ1 -1/µ 0 2 ^ ^ ^ 0 -1/µ3 W ( ) = W0 ( ), W0 = . . . . . . 0 0

we obtain · · · . · · · . · · · 1/µ 0 0 . . .
1

. · · · -1/µn

.

(6.29)

Using matrix (6.29) in Theorem 1, we find the differentials H = F - sI U U+ -
U =U
+

F - sI U

U- - (U+ - U- )s,
U =U
-

(6.30)

H Condition H
add

add

=

2- ^T - s W0 U- . µl-
µ1 µ

(6.31) ,...,
µn µ T

= 0 with expression (6.31) yields U- =

- . Then
T

using (6.9) and (6.13) in (6.30), one can check that U+ = µ1 , . . . , µn + and µ µ the expression for s coincides with the linearization of (6.19). This shows that the approximation of Sn,1 shock states and speeds obtained by Theorem 1 agrees with the analytical results. Now let us assume that the viscosity matrix D(U ) I suffers a small variation D(U ), while the functions G(U ) = U and F (U ) remain unchanged. Then the change of states and speeds of Sn,1 shocks can be found using Theorem 2 as H = 0, H
add +

(6.32)
T

=-
-

^T W0 DU

µ1 µn ,..., µ µ

2 2 - - - s( - - ) , µl µl-

(6.33)

where DU = D(U ( )). One can see that if D(U ) = (U )I , where (U ) is a scalar function, the right-hand side of (6.33) vanishes. Indeed, one can check analytically that the states and speeds of Sn,1 shocks remain unchanged for the viscosity matrix D(U ) = (U )I with an arbitrary positive function (U ) (the only change is in equation (6.18) whose left-hand side becomes (U ) ). In general, the perturbation D(U ) of the viscosity matrix changes the parameters of Sn,1 shocks. Since the change of left and right states and speeds can be found using Theorem 2, one can use the continuation method to find Sn,1 shocks in a system with a particular nontrivial D(U ). For this purpose, the matrix D(U ) is changed from I to the required value in a sequence of small steps; at each step the ^ viscous profile and the matrix W ( ) have to be recomputed and used in Theorem 2 for finding approximation for the next step. A similar procedure can be applied to all other types of Si,j shocks. We can expect to find Si,j shocks in the system, at least for viscosity matrices D(U ) close to (U )I with some function (U ) and for small perturbations of flux functions. 19


We note that the umbilic points give rise to elliptic regions in state space under a general perturbation the flux function F (U ). This happens, for instance, in Stone's model for three-phase flow through porous media [17]. Generically, the resonant points in a perturbed system, where all n characteristic merge, will have onedimensional eigenspace (one eigenvector) unlike the umbilic point U umb , which has n eigenvectors. At the same time, all Si,j shocks persist in the system. Hence Si,j shocks are not necessarily related to the existence of umbilic points.

6.5

Physical applications

In many physical applications of system (6.1), the goal is to study the possibility of recovering one phase of the fluid (e.g. oil) on the right side by injecting different fluid phases (e.g. water and steam) on the left. The variety of Si,j shocks that can appear in the system has to be taken into account in analysis of physical behavior, e.g., when predicting oil recovery in petroleum engineering practice. By adopting the method of [7], one should be able to show stability of Si,j dualfamily shocks with i > j as solutions of equation (2.1), at least for shocks of small amplitudes. We expect that this is true also in case of perturbed flux functions and viscosity matrices.

7

Conclusion

We studied a general class of shock waves satisfying the viscous profile admissibility criterion in general systems of n conservation laws. Roughly speaking, shock waves are classified by comparing their speeds with characteristic speeds at opposite sides of the wave. As a result, we are led to consider dual-family shocks Si,j associated with characteristic families i and j at the left and right sides, respectively. We develop a constructive method for analytical and numerical study of such shocks and their perturbations under change of problem parameters. One remarkable feature of Si,j shocks with i > j is that their left and right states and speeds depend on the viscosity matrix. The other is that these shocks dramatically enlarge the possible structures of generic Riemann problem solutions. In Riemann solutions with such shock waves, classical wave groups of the same characteristic family can appear repeatedly from different sides of dual-family shocks, separated by constant states. Thus, classical wave groups in solutions with dual-family shocks does not necessarily follow in increasing family number. As a wide variety of Si,j shocks is found in systems describing multiphase flows through porous media, there is a promising perspective for analyzing and using these shocks in physical applications.

Acknowledgments
The authors thank Johannes Bruining and Marcelo Viana for helpful discussions.

20


8
8.1

App endix
Pro of of Prop osition 1
V = B (U ( ), s)V , (8.1)

Let us define the linear spaces V- and V+ of solutions V ( ) Rn of the system

vanishing as - and +, respectively. Similarly, we define the linear spaces W- and W+ of solutions W ( ) Rn of the system W = -B T (U ( ), s)W, (8.2)

vanishing as - and +, respectively. Since the matrix B (U, s) has the eigenvalues µi (U, s), i = 1, . . . , n, and the matrix -B T (U, s) has the eigenvalues -µi (U, s), i = 1, . . . , n satisfying inequalities (2.6), we find dim V- = n - i + 1, dim V+ = j , dim W- = i - 1, and dim W+ = n - j . The set W is the intersection W = W- W+ . Since systems (8.1) and (8.2) are adjoint, the linear space W- is the orthogonal complement of V- at each value of . Indeed, let V ( ) V- and W ( ) W- . Integrating the product W T ( )B (U ( ), s)V ( ) in the interval - < , we obtain T W ( )B (U ( ), s)V ( )d = W T ( )V ( )d
- -

= W ( )V ( ) = W ( )V ( ) +
T

T









-


-
-

W T ( )V ( )d

(8.3)

W T ( )B (U ( ), s)V ( )d .

-

Here, the integration by parts, equations (8.1), (8.2), and the conditions V (-) = W (-) = 0 were used. From (8.3), we obtain that W T ( )V ( ) = 0 (8.4)

for any real . Since the tangent space to the manifold Mu (U- ) at U ( ) is given by the vectors V ( ) V- , the vector W ( ) W is orthogonal to Mu (U- ) at U ( ). Analogously, we can prove that W ( ) is orthogonal to Ms (U+ ) at U ( ). If the manifolds Mu (U- ) and Ms (U+ ) intersect quasi-transversally, we have dim(V- V+ ) = 1. Hence, dim W = dim(W- W+ ) = n - dim(V- V+ ) = n - dim V- + dim V+ - dim(V- V+ ) = i - j. (8.5)

21


8.2

Pro of of Theorems 1 and 2
condiLet us + s . Then

Equations (4.5) and (4.7) are obtained by varying the Rankine­Hugoniot tions (3.1). Let U ( ) be a viscous profile of an Si,j shock with states U± and speed s. take perturbations of the left and right states U± + U± and of the speed s These perturbations are assumed to satisfy the Rankine-Hugoniot conditions. we can write the variational equation for system (2.3) as U = B (U ( ), s)U - D
-1 U

F G -s U U

- U- - DU 1 (GU - G- )s. (8.6) U =U-

The perturbed solution U ( ) + Uu ( ) lying in the unstable manifold Mu (U- + U- ) satisfies the condition Uu (-) = U- . (8.7) Pre-multiplying (8.6) by the transpose of the function W ( ) W and integrating in the interval - < , we find W ( )Uu ( ) = -
-


T









- W T DU 1

F G -s U U

d
U =U
-

U

-

(8.8)

-
-

- W T DU 1 (GU - G- )d s.

Here, we used integration by parts with equations (8.2), (8.7) and the condition W (-) = 0. Analogously, for the perturbed solution U ( ) + Us ( ) lying in the stable manifold Ms (U+ + U+ ), i.e., Us (+) = U+ , we obtain W T ( )Us ( ) =
+


WTD

-1 U

G F -s U U

d
U =U
-

U

-

+

(8.9)

+



- W T DU 1 (GU - G- )d s.

Subtracting (8.8) from (8.9), we get
P P W T ( )(U+ - U- ) = W T ( )(Us ( ) - Uu ( )) +

=
-

WTD

-1 U

F G -s U U

d
U =U-

U

-

(8.10)

+

+
-

- W T DU 1 (GU - G- )d s.

Using expression (8.10) in (4.3), we obtain formula (4.6). Finally, formula (4.8) is derived in the same way from equation (8.6), where the - - - term -DU 1 DU DU 1 FU - F- - s(GU - G- ) + DU 1 FU - F- - s(GU - G- ) resulting from variations of the system functions is added to the right-hand side. 22


References
[1] R. H. Brooks and A. T. Corey, Properties of Porous Media Affecting Fluid Flow, J. Irrig. Drain. Div. 6 (1966), 61­88. [2] J. Guckenheimer and P. Holmes, Nonlinear Oscil lations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 1983. [3] J. M. Hurley and B. J. Plohr, Some effects of viscous terms on Riemann problem solutions, Mat. Contemp. 8 (1995), 203­224. [4] E. L. Isaacson, D. Marchesin, and B. J. Plohr, Transitional waves for conservation laws, SIAM J. Math. Anal. 21 (1990), 837­866. [5] A. G. Kulikovskij, Surfaces of discontinuity separating two perfect media of different properties. Recombination waves in magnetohydrodynamics, Prikl. Mat. Mekh. 32 (1968), 1125­1131 (in Russian). [6] A. G. Kulikovskii, N. V. Pogorelov, and A. Yu. Semenov, Mathematical aspects of numerical solution of hyperbolic systems. Chapman & Hall/CRC, Boca Raton, 2001. [7] T.-P. Liu and K. Zumbrun, On nonlinear stability of general undercompressive viscous shock waves, Comm. Math. Phys. 174 (1995), 319­345. [8] A. A. Mailybaev and D. Marchesin, Dual-family viscous shock waves in systems of conservation laws: a surprising example, in: Proc. Conf. on Analysis, Modeling and Computation of PDE and Multiphase Flow, Stony Brook, NY, 2004, to appear. [9] A. A. Mailybaev and D. Marchesin, Generic solutions of the Riemann problem for n equations, in preparation. [10] A. Ma jda and R. L. Pego, Stable viscosity matrices for systems of conservation laws, J. Differential Equations 56 (1985), 229­262. [11] B. Plohr and D. Marchesin, Wave structure in WAG recovery, SPE 71314, Society of Petroleum Engineers Journal 6 (2001), 209­219. [12] S. Schecter, D. Marchesin, and B. J. Plohr, Structurally stable Riemann solutions, J. Differential Equations 126 (1996), 303­354. [13] S. Schecter, B. J. Plohr, and D. Marchesin, Classification of codimension-one Riemann solutions, J. Dynam. Differential Equations 13 (2001), 523­588. [14] M. Shearer, D. G. Schaeffer, D. Marchesin, and P. L. Paes-Leme, Solution of the Riemann problem for a prototype 2 â 2 system of nonstrictly hyperbolic conservation laws, Arch. Rational Mech. Anal. 97 (1987), 299­320. [15] J. Smoller, Shock Waves and Reaction-Diffusion Equations, Springer, New York, 1983. 23


[16] J. Sotomayor, Generic bifurcations of dynamical systems, in: M. M. Peixoto (Ed.), Dynamical Systems (Proc. Sympos., Univ. Bahia, Salvador, 1971), Academic Press, New York, 1973, pp. 561­582. [17] H. L. Stone, Probability Model for Estimating Three-Phase Relative Permeability, Trans. SPE of AIME, 249 (1970), 214­218.

24