Bounded domain problem for the modified Buckley-Leverett equation
BBOUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION YING WANG AND CHIU-YEN KAO Abstract.
The focus of the present study is the modified Buckley-Leverett(MBL) equation describing two-phase flow in porous media. The MBL equa-tion differs from the classical Buckley-Leverett (BL) equation by including abalanced diffusive-dispersive combination. The dispersive term is a third ordermixed derivatives term, which models the dynamic effects in the pressure dif-ference between the two phases. The classical BL equation gives a monotonewater saturation profile for any Riemann problem; on the contrast, when thedispersive parameter is large enough, the MBL equation delivers non-monotonewater saturation profile for certain Riemann problems as suggested by the ex-perimental observations. In this paper, we first show that the solution of thefinite interval [0 , L ] boundary value problem converges to that of the half-line[0 , + ∞ ) boundary value problem for the MBL equation as L → + ∞ . Thisresult provides a justification for the use of the finite interval boundary valueproblem in numerical studies for the half line problem. Furthermore, we extendthe classical central schemes for the hyperbolic conservation laws to solve theMBL equation which is of pseudo-parabolic type. Numerical results confirmthe existence of non-monotone water saturation profiles consisting of constantstates separated by shocks. Introduction The classical Buckley-Leverett (BL) equation [3] is a simple model for two-phase fluid flow in a porous medium. One application is secondary recovery by water-drive in oil reservoir simulation. In one space dimension the equation has the standard conservation form u t + ( f ( u )) x = 0 in Q = { ( x, t ) : x > , t > } u ( x,
0) = 0 x ∈ (0 , ∞ )(1.1) u (0 , t ) = u B t ∈ [0 , ∞ )with the flux function f ( u ) being defined as f ( u ) = u < , u u + M (1 − u ) ≤ u ≤ , u > . (1.2)In this content, u : ¯ Q → [0 ,
1] denotes the water saturation (e.g. u = 1 means pure water, and u = 0 means pure oil), u B is a constant which indicates water saturation at x = 0, and M > (1.1) is a prototype for conservation laws with convex-concave flux functions. The graph of f ( u ) and f (cid:48) ( u ) with M = 2 is given in Figure 1.1. Date : November 3, 2018.2000
Mathematics Subject Classification.
Key words and phrases. conservation laws, dynamic capillarity, two-phase flows, porous media,shock waves, pseudo-parabolic equations, central schemes.This work is supported in part by NSF Grant DMS-0811003 and an Alfred P. Sloan Fellowship. a r X i v : . [ m a t h . NA ] S e p YING WANG AND CHIU-YEN KAO (a) f ( u ) = u u + M (1 − u ) α u f ( u ) (b) f (cid:48) ( u ) = Mu (1 − u )( u + M (1 − u ) ) f ′ ( u ) Figure 1.1. f ( u ) and f (cid:48) ( u ) with M = 2.Due to the possibility of the existence of shocks in the solution of the hyperbolicconservation laws (1.1), the weak solutions are sought. The function u ∈ L ∞ ( Q ) iscalled a weak solution of the conservation laws (1.1) if (cid:90) Q (cid:26) u ∂φ∂t + f ( u ) ∂φ∂x (cid:27) = 0 for all φ ∈ C ∞ ( Q ) . Notice that the weak solution is not unique. Among the weak solutions, the entropy solution is physically relevant and unique. The weak solution that satisfies Oleinik entropy condition [19] f ( u ) − f ( u l ) u − u l ≥ s ≥ f ( u ) − f ( u r ) u − u r for all u between u l and u r (1.3)is the entropy solution, where u l , u r are the function values to the left and right of the shock respectively, and the shock speed s satisfies Rankine-Hugoniot jump condition [17, 10] (1.4) s = f ( u l ) − f ( u r ) u l − u r . The classical BL equation (1.1) with flux function f ( u ) as given in (1.2) has been well studied (see [14] for an introduction). Let α be the solution of f (cid:48) ( u ) = f ( u ) u , i.e., (1.5) α = (cid:114) MM + 1 . The entropy solution of the classical BL equation can be classified into two cate- gories: (1) If 0 < u B ≤ α , the entropy solution has a single shock at xt = f ( u B ) u B . (2) If α < u B <
1, the entropy solution contains a rarefaction between u B and α for f (cid:48) ( u B ) < xt < f (cid:48) ( α ) and a shock at xt = f ( α ) α . These two types of solutions are shown in Figure 1.2 for M = 2. In either case, the entropy solution of the classical BL equation (1.1) is a non-increasing function of x at any given time t >
0. However, the experiments of two-phase flow in porous medium reveal complex infiltration profiles, which may involve overshoot, i.e., profiles may not be monotone [7]. This suggests the need of modification to the classical BL equation (1.1). OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION 3 (a) u B = 0 . xt u (b) u B = 0 . xt u Figure 1.2.
The entropy solution of the classical BL equation( M = 2 , α = (cid:113) ≈ . < uB = 0 . ≤ α , the solutionconsists of one shock at xt = f ( u B ) u B ; (b) α < uB = 0 . <
1, thesolution consists of a rarefaction between u B and α for f (cid:48) ( u B ) < xt < f (cid:48) ( α ) and a shock at xt = f ( α ) α .To better describe the infiltration profiles, we go back to the origins of (1.1). Let S i be the saturation of water/oil ( i = w, o ) and assume that the medium is completely saturated, i.e. S w + S o = 1. The conservation of mass gives (1.6) φ ∂S i ∂t + ∂q i ∂x = 0where φ is the porosity of the medium (relative volume occupied by the pores) and q i denotes the discharge of water/oil with q w + q o = q , which is assumed to be a constant in space due to the complete saturation assumption. Throughout of this work, we consider it constant in time as well. By Darcy’s law (1.7) q i = − k k ri ( S i ) µ i ∂P i ∂x , i = w, o where k denotes the absolute permeability, k ri is the relative permeability and µ i is the viscosity of water/oil. Instead of considering constant capillary pressure as adopted by the classical BL equation (1.1), Hassanizadeh and Gray [8, 9] have defined the dynamic capillary pressure as (1.8) P c = P o − P w = p c ( S w ) − φτ ∂S w ∂t where p c ( S w ) is the static capillary pressure and τ is a positive constant, and ∂S w ∂t is the dynamic effects. Using Corey [6, 20] expressions with exponent 2, k rw ( S w ) = S w , k ro ( S o ) = S o , rescaling x φq → x and combining (1.6)-(1.8), the single equation for the water saturation u = S w is (1.9) ∂u∂t + ∂∂x (cid:20) u u + M (1 − u ) (cid:21) = − ∂∂x (cid:20) φ q k (1 − u ) u µ w (1 − u ) + µ o u ∂∂x (cid:18) p c ( u ) φ − τ ∂u∂t (cid:19)(cid:21) where M = µ w µ o [22]. Linearizing the right hand side of (1.9) and rescaling the equation as in [21, 20], the modified Buckley-Leverett equation (MBL) is derived YING WANG AND CHIU-YEN KAO as (1.10) ∂u∂t + ∂f ( u ) ∂x = (cid:15) ∂ u∂x + (cid:15) τ ∂ u∂x ∂t where the water fractional flow function f ( u ) is given as in (1.2). Notice that, if P c in (1.8) is taken to be constant, then (1.9) gives the classical BL equation; while if the dispersive parameter τ is taken to be zero, then (1.10) gives the viscous BL equation, which still displays monotone water saturation profile. Thus, in addi- tion to the classical second order viscous term (cid:15)u xx , the MBL equation (1.10) is an extension involving a third order mixed derivative term (cid:15) τ u xxt . Van Dujin et al. [21] showed that the value τ is critical in determining the type of the solution profile. In particular, for certain Riemann problems, the solution profile of (1.10) is not monotone when τ is larger than the threshold value τ ∗ , where τ ∗ was numer- ically determined to be 0.61 [21]. The non-monotonicity of the solution profile is consistent with the experimental observations [7]. The classical BL equation (1.1) is hyperbolic, and the numerical schemes for hyperbolic equations have been well developed (e.g. [14, 15, 4, 5, 18, 12] ). The MBL equation (1.10), however, is pseudo-parabolic, we will illustrate how to extend the central schemes [18, 12, 13] to solve (1.10) numerically. Unlike the finite domain of dependence for the classical BL equation (1.1), the domain of dependence for the MBL equation (1.10) is infinite. This naturally raises the question for the choice of computational domain. To answer this question, we will first study the MBL equation equipped with two types of domains and corresponding boundary conditions. One is the half line boundary value problem u t + ( f ( u )) x = (cid:15)u xx + (cid:15) τ u xxt in Q = { ( x, t ) : x > , t > } u ( x,
0) = u ( x ) x ∈ [0 , ∞ ) u (0 , t ) = g u ( t ) , lim x →∞ u ( x, t ) = 0 t ∈ [0 , ∞ ) u (0) = g u (0) compatibility condition(1.11)and the other one is finite interval boundary value problem v t + ( f ( v )) x = (cid:15)v xx + (cid:15) τ v xxt in (cid:101) Q = { ( x, t ) : x ∈ (0 , L ) , t > } v ( x,
0) = v ( x ) x ∈ [0 , L ] v (0 , t ) = g v ( t ) , v ( L, t ) = h ( t ) t ∈ [0 , ∞ ) v (0) = g v (0) , v ( L ) = h (0) compatibility condition . (1.12)Considering (1.13) u ( x ) = (cid:26) v ( x ) for x ∈ [0 , L ]0 for x ∈ [ L, + ∞ ) , g u ( t ) = g v ( t ) ≡ g ( t ) , h ( t ) ≡ , we will show the relation between the solutions of problems (1.11) and (1.12). To the best knowledge of the authors, there is no such study for MBL equation (1.10). Similar questions were answered for BBM equation [1, 2]. The organization of this paper is as follows. Section 2 will bring forward the exact theory comparing the solutions of (1.11) and (1.12). The difference between the solutions of these two types of problems decays exponentially with respect to the length of the interval L for practically interesting initial profiles. This provides a theoretical justification for the choice of the computational domain. In section OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION 5
3, high order central schemes will be developed for MBL equation in finite interval domain. We provide a detailed derivation on how to extend the central schemes [18, 12] for conservation laws to solve the MBL equation (1.10). The idea of adopting numerical schemes originally designed for hyperbolic equations to pseudo-parabolic equations is not restricted to central type schemes only ([23, 24]). The numerical results in section 4 show that the water saturation profile strongly depends on the dispersive parameter τ value as studied in [21]. For τ > τ ∗ , the MBL equation (1.10) gives non-monotone water saturation profiles for certain Riemann problems as suggested by experimental observations [7]. Section 5 gives the conclusion of the paper and the possible future directions. The half line problem versus the finite interval problem Let u ( x, t ) be the solution to the half line problem (1.11), and let v ( x, t ) be the solution to the finite interval problem (1.12). We consider the natural assumptions (1.13). The goal of this section is to develop an estimate of the difference between u and v on the spatial interval [0 , L ] at a given finite time t . The main result of this section is Theorem 2.1 (The main Theorem) . If u ( x ) satisfies (2.1) u ( x ) = (cid:26) C u x ∈ [0 , L ]0 x > L where L < L and C u , are positive constants, then (cid:107) u ( · , t ) − v ( · , t ) (cid:107) H L,(cid:15),τ ≤ D (cid:15),τ ( t ) e − λL(cid:15) √ τ + D (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ for some < λ < , D (cid:15),τ ( t ) > and D (cid:15),τ ( t ) > , where (cid:107) Y ( · , t ) (cid:107) H L,(cid:15),τ := (cid:115)(cid:90) L Y ( x, t ) + ( (cid:15) √ τ Y x ( x, t )) dx Notice that the initial condition (2.1) we considered is the Riemann problem.
Theorem 2.1 shows that the solution to the half line problem (1.11) can be approx- imated as accurately as one wants by the solution to the finite interval problem (1.12) in the sense that D (cid:15),τ ( t ), D (cid:15),τ ( t ), λL(cid:15) √ τ and λ ( L − L ) (cid:15) √ τ can be controlled. To prove theorem 2.1, we first derive the implicit solution formulae for the half line problem and the finite interval problem in section 2.1 and section 2.2 respec- tively. The implicit solution formulae are in integral form, which are derived by separating the x -derivative from the t -derivative, and formally solving a first order linear ODE in t and a second order non-homogeneous ODE in x . In section 2.3, we use Gronwall’s inequality multiple times to obtain the desired result in theorem Half line problem.
In this section, we derive the implicit solution formula for the half line problem (1.11) (with g u ( t ) = g ( t )). To solve (1.11), we first rewrite (1.11) by separating the x -derivative from the t -derivative, (2.2) (cid:18) I − (cid:15) τ ∂ ∂x (cid:19) (cid:18) u t + 1 (cid:15)τ u (cid:19) = 1 (cid:15)τ u − ( f ( u )) x . YING WANG AND CHIU-YEN KAO
By using integrating factor method, we formally integrate (2.2) over [0 , t ] to obtain (2.3) (cid:18) I − (cid:15) τ ∂ ∂x (cid:19) (cid:16) u − e − t(cid:15)τ u (cid:17) = (cid:90) t (cid:18) (cid:15)τ u − ( f ( u )) x (cid:19) e − t − s(cid:15)τ ds. Furthermore, we let (2.4) A = u − e − t(cid:15)τ u , then (2.3) can be written as (2.5) A (cid:48)(cid:48) − (cid:15) τ A = (cid:90) t (cid:18) − (cid:15) τ u + 1 (cid:15) τ ( f ( u )) x (cid:19) e − t − s(cid:15)τ ds, where (cid:48) = ∂∂x . Notice that (2.5) is a second-order non-homogeneous ODE in x -variable along with the boundary conditions A (0 , t ) = u (0 , t ) − e − t(cid:15)τ u (0) = g ( t ) − e − t(cid:15)τ g (0) ,A ( ∞ , t ) = u ( ∞ , t ) − e − t(cid:15)τ u ( ∞ ) = 0 . (2.6)To solve (2.5), we first solve the corresponding linear homogeneous equation with the non-zero boundary conditions (2.6). We then find a particular solution for the non-homogeneous equation with zero boundary conditions by introducing a Green’s function G ( x, ξ ) and a kernel K ( x, ξ ) for the non-homogeneous terms u and ( f ( u )) x respectively. Combining the solutions for the two non-homogeneous terms and the homogeneous part with boundary conditions, we get the solution for equation (2.5) satisfying the boundary conditions (2.6): A ( x, t ) = − (cid:15) τ (cid:90) t (cid:90) + ∞ G ( x, ξ ) u ( ξ, s ) e − t − s(cid:15)τ dξ ds + 1 (cid:15) τ (cid:90) t (cid:90) + ∞ K ( x, ξ ) f ( u ) e − t − s(cid:15)τ dξ ds + (cid:16) g ( t ) − e − t(cid:15)τ g (0) (cid:17) e − x(cid:15) √ τ (2.7)where the Green’s function G ( x, ξ ) and the kernel K ( x, ξ ) are G ( x, ξ ) = (cid:15) √ τ (cid:16) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:17) , (2.8) K ( x, ξ ) = − ∂G ( x, ξ ) ∂ξ = 12 (cid:16) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:17) . (2.9)To recover the solution for the half line problem (1.11), we refer to the definition of A in (2.4). Thus, the implicit solution formula for the half line problem (1.11) is u ( x, t ) = − (cid:15) τ √ τ (cid:90) t (cid:90) + ∞ (cid:16) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:17) u ( ξ, s ) e − t − s(cid:15)τ dξ ds + 12 (cid:15) τ (cid:90) t (cid:90) + ∞ (cid:16) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:17) f ( u ) e − t − s(cid:15)τ dξ ds + (cid:16) g ( t ) − e − t(cid:15)τ g (0) (cid:17) e − x(cid:15) √ τ + e − t(cid:15)τ u ( x ) . (2.10) OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION 7
Finite interval problem.
The implicit solution for the finite interval problem (1.12) (with g v ( t ) = g ( t )) can be solved in a similar way. The only difference is that the additional boundary condition h ( t ) at x = L in (1.12) gives different boundary conditions for the non-homogeneous ODE in x -variable. Denote (2.11) A L = v − e − t(cid:15)τ v , then it satisfies (2.12) ( A L ) (cid:48)(cid:48) − (cid:15) τ A L = (cid:90) t (cid:18) − (cid:15) τ v + 1 (cid:15) τ ( f ( v )) x (cid:19) e − t − s(cid:15)τ ds where (cid:48) = ∂∂x with the boundary conditions A L (0 , t ) = v (0 , t ) − e − t(cid:15)τ v (0) = g ( t ) − e − t(cid:15)τ g (0) ,A L ( L, t ) = v ( L, t ) − e − t(cid:15)τ v ( L ) = h ( t ) − e − t(cid:15)τ h (0) . These boundary conditions affect both the homogeneous solution and the par- ticular solution of (2.12) as follows A L ( x, t ) = − (cid:15) τ (cid:90) t (cid:90) L G L ( x, ξ ) v ( ξ, s ) e − t − s(cid:15)τ dξ ds + 1 (cid:15) τ (cid:90) t (cid:90) L K L ( x, ξ ) f ( v ) e − t − s(cid:15)τ dξ ds + c ( t ) φ ( x ) + c ( t ) φ ( x )(2.13)where the Green’s function G L ( x, ξ ), the kernel K L ( x, ξ ) and the bases for the homogeneous solutions are (2.14) G L ( x, ξ ) = (cid:15) √ τ e L(cid:15) √ τ − (cid:16) e x + ξ(cid:15) √ τ + e L − ( x + ξ ) (cid:15) √ τ − e | x − ξ | (cid:15) √ τ − e L −| x − ξ | (cid:15) √ τ (cid:17) , K L ( x, ξ ) = − e L(cid:15) √ τ − (cid:16) e x + ξ(cid:15) √ τ − e L − ( x + ξ ) (cid:15) √ τ +sgn( x − ξ ) e | x − ξ | (cid:15) √ τ − sgn( x − ξ ) e L −| x − ξ | (cid:15) √ τ (cid:17) , (2.15) c ( t ) = g ( t ) − e − t(cid:15)τ g (0) , c ( t ) = h ( t ) − e − t(cid:15)τ h (0) , (2.16) φ ( x ) = e L − x(cid:15) √ τ − e − L + x(cid:15) √ τ e L(cid:15) √ τ − e − L(cid:15) √ τ , and φ ( x ) = e x(cid:15) √ τ − e − x(cid:15) √ τ e L(cid:15) √ τ − e − L(cid:15) √ τ . (2.17)Thus, the implicit solution formula for the finite interval problem (1.12) is v ( x, t ) = − (cid:15) τ √ τ ( e L(cid:15) √ τ − (cid:90) t (cid:90) L (cid:16) e x + ξ(cid:15) √ τ + e L − ( x + ξ ) (cid:15) √ τ − e | x − ξ | (cid:15) √ τ − e L −| x − ξ | (cid:15) √ τ (cid:17) v ( ξ, s ) e − t − s(cid:15)τ dξ ds − (cid:15) τ ( e L(cid:15) √ τ − (cid:90) t (cid:90) L (cid:16) e x + ξ(cid:15) √ τ − e L − ( x + ξ ) (cid:15) √ τ + sgn( x − ξ ) e | x − ξ | (cid:15) √ τ − sgn( x − ξ ) e L −| x − ξ | (cid:15) √ τ (cid:17) f ( v ) e − t − s(cid:15)τ dξ ds + c ( t ) φ ( x ) + c ( t ) φ ( x ) + e − t(cid:15)τ v ( x ) . (2.18) YING WANG AND CHIU-YEN KAO
Comparisons.
In this section, we will prove that the solution u ( x, t ) to the half line problem can be approximated as accurately as one wants by the solution v ( x, t ) to the finite interval problem as stated in Theorem 2.1. Due to the difference in the integration domains, we do not use (2.10) and (2.18) directly for the comparison. Instead, we decompose u ( x, t ) ( v ( x, t ) respectively) into two parts: U ( x, t ) and u L ( x, t ) ( V ( x, t ) and v L ( x, t ) respectively), such that U ( x, t ) ( V ( x, t ) respectively) enjoys zero initial condition and boundary conditions at x = 0 and x = L . We estimate the difference between u ( · , t ) and v ( · , t ) by estimating the differences between u L ( · , t ) and v L ( · , t ), U ( · , t ) and V ( · , t ), then applying the triangle inequality. Definitions and lemmas.
To assist the proof of Theorem 2.1 in section 2.3.3,we introduce some new notations in this section. We first decompose u ( x, t ) as sumof two terms U ( x, t ) and u L ( x, t ), such that u ( x, t ) = U ( x, t ) + u L ( x, t ) x ∈ [0 , + ∞ )where (2.19) u L = e − t(cid:15)τ u ( x ) + c ( t ) e − x(cid:15) √ τ + (cid:16) u ( L, t ) − c ( t ) e − L(cid:15) √ τ − e − t(cid:15)τ u ( L ) (cid:17) φ ( x )and c ( t ) and φ ( x ) are given in (2.16) and (2.17) respectively. With this definition, u L takes care of the initial condition u ( x ) and boundary conditions g ( t ) at x = 0 and x = L for u ( x, t ). Then U satisfies an equation slightly different from the equation u satisfies in (1.11): U t − (cid:15)U xx − (cid:15) τ U xxt = (cid:0) u t − (cid:15)u xx − (cid:15) τ u xxt (cid:1) − (cid:0) ( u L ) t − (cid:15) ( u L ) xx − (cid:15) τ ( u L ) xxt (cid:1) = − ( f ( u )) x + 1 (cid:15)τ u L ( x, t )(2.20)In addition, U ( x, t ) has zero initial condition and boundary conditions at x = 0 and x = L , i.e., U ( x,
0) = 0 , U (0 , t ) = 0 , U ( L, t ) = 0 . (2.21)Similarly, for v ( x, t ), let v ( x, t ) = V ( x, t ) + v L ( x, t ) x ∈ [0 , L ]where (2.22) v L = e − t(cid:15)τ v ( x ) + c ( t ) φ ( x ) + c ( t ) φ ( x )and c ( t ), c ( t ) and φ ( x ), φ ( x ) are given in (2.16) and (2.17) respectively. With this definition, v L takes care of the initial condition v ( x ) and boundary conditions g ( t ) and h ( t ) at x = 0 and x = L for v ( x, t ). Then V satisfies an equation slightly different from the equation v satisfies in (1.12): V t − (cid:15)V xx − (cid:15) τ V xxt = − ( f ( v )) x + 1 (cid:15)τ v L ( x, t )(2.23)with V ( x,
0) = 0 , V (0 , t ) = 0 , V ( L, t ) = 0 . (2.24) OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION 9
Since, in the end, we want to study the difference between U ( x, t ) and V ( x, t ), wedefine W ( x, t ) = V ( x, t ) − U ( x, t ) for x ∈ [0 , L ] . Because of (2.20) and (2.23), we have (2.25) W t − (cid:15)W xx − (cid:15) τ W xxt = − ( f ( v ) − f ( u )) x + 1 (cid:15)τ ( v L − u L ) . In lieu of (2.21) and (2.24), W ( x, t ) also has zero initial condition and boundary conditions at x = 0 and x = L , i.e., W ( x,
0) = 0 , W (0 , t ) = 0 , W ( L, t ) = 0 . (2.26)Now, to estimate (cid:107) u − v (cid:107) , we can estimate (cid:107) W (cid:107) = (cid:107) V − U (cid:107) and estimate (cid:107) u L − v L (cid:107) separately. These estimates are done in section 2.3.3. Next, we state the lemmas needed in the proof of Theorem 2.1. The proof of the lemmas can be found in the appendix A and [22]. In all the lemmas, we assume < λ < u ( x ) satisfies u ( x ) = (cid:26) C u x ∈ [0 , L ]0 x > L (2.27)where L < L and C u are positive constants. Notice that the constraint λ ∈ (0 , is crucial in Lemmas 2.3, 2.4. Lemma 2.2. f ( u ) = u u + M (1 − u ) ≤ Du where D = f ( α ) α and α = (cid:113) MM +1 . Lemma 2.3. (i) (cid:82) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − λξ(cid:15) √ τ dξ ≤ (cid:15) √ τ − λ . (ii) (cid:82) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − ξ(cid:15) √ τ dξ ≤ (cid:15) √ τe (1 − λ ) . (iii) (cid:82) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx(cid:15) √ τ | u ( ξ ) | dξ ≤ C u (cid:15) √ τ e λL (cid:15) √ τ . Lemma 2.4. (i) (cid:82) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − λξ(cid:15) √ τ dξ ≤ (cid:15) √ τ − λ . (ii) (cid:82) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − ξ(cid:15) √ τ dξ ≤ (cid:15) √ τ + (cid:15) √ τe (1 − λ ) . (iii) (cid:82) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx(cid:15) √ τ | u ( ξ ) | dξ ≤ C u (cid:15) √ τ e λL (cid:15) √ τ . Lemma 2.5. (i) (cid:12)(cid:12)(cid:12) φ ( x ) − e − x(cid:15) √ τ (cid:12)(cid:12)(cid:12) = e − L(cid:15) √ τ | φ ( x ) | . (ii) | φ ( x ) | ≤ for x ∈ [0 , L ] . (iii) | φ (cid:48) ( x ) | ≤ (cid:15) √ τ if (cid:15) (cid:28) for x ∈ [0 , L ] . Last but not least, the norm that we will use in Theorem 2.1 and its proof is (2.28) (cid:107) Y ( · , t ) (cid:107) H L,(cid:15),τ := (cid:115)(cid:90) L Y ( x, t ) + ( (cid:15) √ τ Y x ( x, t )) dx. A proposition.
In this section, we will give a critical estimate, which is es- sential in the calculation of maximum difference (cid:107) u L ( · , t ) − v L ( · , t ) (cid:107) ∞ in section u L ( x, t ) and v L ( x, t ) given in (2.19) and (2.22) respectively, it is clear that the coefficient u ( L, t ) − c ( t ) e − L(cid:15) √ τ − e − t(cid:15)τ u ( L ) for φ ( x ) appeared in (2.19) needs to be compared with the corresponding coefficient c ( t ) for φ ( x ) appeared in (2.22). We thus define a space-dependent function (2.29) U c ( x, t ) = u ( x, t ) − c ( t ) e − x(cid:15) √ τ − e − t(cid:15)τ u ( x )and establish the following proposition Proposition 2.6. (2.30) | U c ( L, t ) | ≤ a τ ( t ) e bτ t(cid:15)τ e − λL(cid:15) √ τ + c τ t(cid:15)τ e ( bτ − t(cid:15)τ e − λ ( L − L (cid:15) √ τ for some parameter-dependent constants a τ , b τ and c τ . Proof.
Based on the implicit solution formula (2.10) derived in section 2.1, Lemma U c and u given in (2.29), we can get an inequality in terms of U c | U c ( x, t ) | ≤ (cid:15) τ √ τ (cid:20)(cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) | U c ( ξ, s ) | e − t − s(cid:15)τ dξ ds + (cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) | c ( s ) | e − ξ(cid:15) √ τ e − t − s(cid:15)τ dξ ds + (cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) | u ( ξ ) | e − t(cid:15)τ dξ ds (cid:21) + D (cid:15) τ (cid:20)(cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) | U c ( ξ, s ) | e − t − s(cid:15)τ dξ ds + (cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) | c ( s ) | e − ξ(cid:15) √ τ e − t − s(cid:15)τ dξ ds + (cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) | u ( ξ ) | e − t(cid:15)τ dξ ds (cid:21) . (2.31)To show that U c ( x, t ) decays exponentially with respect to x , we pull out an exponential term by writing U c ( x, t ) = e − λx(cid:15) √ τ e − t(cid:15)τ ˜ U ( x, t ), where 0 < λ <
1, such that (2.32) ˜ U ( x, t ) = e λx(cid:15) √ τ e t(cid:15)τ U c ( x, t ) , then (2.31) can be rewritten in terms of ˜ U ( x, t ) as follows (cid:12)(cid:12)(cid:12) ˜ U ( x, t ) (cid:12)(cid:12)(cid:12) ≤ (cid:15) τ √ τ (cid:20)(cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − λξ(cid:15) √ τ (cid:12)(cid:12)(cid:12) ˜ U ( ξ, s ) (cid:12)(cid:12)(cid:12) dξ ds + (cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) | c ( s ) | e λx − ξ(cid:15) √ τ e s(cid:15)τ dξ ds + (cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx(cid:15) √ τ | u ( ξ ) | dξ ds (cid:21) + D (cid:15) τ (cid:20)(cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − λξ(cid:15) √ τ (cid:12)(cid:12)(cid:12) ˜ U ( ξ, s ) (cid:12)(cid:12)(cid:12) dξ ds + (cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) | c ( s ) | e λx − ξ(cid:15) √ τ e s(cid:15)τ dξ ds + (cid:90) t (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx(cid:15) √ τ | u ( ξ ) | dξ ds (cid:21) . (2.33) OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION11
Because of Lemmas 2.3–2.4, we can get the following estimate for (cid:12)(cid:12)(cid:12) ˜ U ( · , t ) (cid:12)(cid:12)(cid:12) ∞ based on (2.33) : (cid:12)(cid:12)(cid:12) ˜ U ( · , t ) (cid:12)(cid:12)(cid:12) ∞ ≤ (cid:15) τ √ τ (cid:20) (cid:15) √ τ − λ (cid:90) t | ˜ U ( · , s ) | ∞ ds + (cid:15) √ τe (1 − λ ) (cid:90) t | c ( s ) | e s(cid:15)τ ds +2 C u (cid:15) √ τ e λL (cid:15) √ τ (cid:90) t ds (cid:21) + D (cid:15) τ (cid:20) (cid:15) √ τ − λ (cid:90) t | ˜ U ( · , s ) | ∞ ds + (cid:15) √ τ (cid:18) e (1 − λ ) (cid:19) (cid:90) t | c ( s ) | e s(cid:15)τ ds +2 C u (cid:15) √ τ e λL (cid:15) √ τ (cid:90) t ds (cid:21) ≤ (cid:90) t b τ (cid:15)τ | ˜ U ( · , s ) | ∞ ds + (cid:90) t ˜ a τ ( s ) (cid:15)τ ds (2.34)where b τ = 1 + D √ τ − λ , ˜ a τ ( t ) = a τ e t(cid:15)τ + c τ e λL (cid:15) √ τ ,a τ = | c ( · ) | ∞ (1 + D √ τ ( e (1 − λ ) + 1))2 e (1 − λ ) , c τ = C u (1 + D √ τ ) . By Gronwall’s inequality, inequality (2.34) gives that (cid:12)(cid:12)(cid:12) ˜ U ( · , t ) (cid:12)(cid:12)(cid:12) ∞ ≤ (cid:90) t ˜ a τ ( t − s ) (cid:15)τ e bτ ( t − s ) (cid:15)τ ds ≤ (cid:18) a τ e t(cid:15)τ + c τ t(cid:15)τ e λL (cid:15) √ τ (cid:19) e bτ t(cid:15)τ Hence | U c ( x, t ) | ≤ (cid:12)(cid:12)(cid:12) ˜ U ( · , t ) (cid:12)(cid:12)(cid:12) ∞ e − λx(cid:15) √ τ e − t(cid:15)τ ≤ (cid:16) a τ e t(cid:15)τ + c τ t(cid:15)τ e λL (cid:15) √ τ (cid:17) e bτ t(cid:15)τ e − λx(cid:15) √ τ e − t(cid:15)τ i.e., U c ( x, t ) decays exponentially with respect to x . In particular, when x = L , we have (2.35) | U c ( L, t ) | ≤ a τ e bτ t(cid:15)τ e − λL(cid:15) √ τ + c τ t(cid:15)τ e ( bτ − t(cid:15)τ e − λ ( L − L (cid:15) √ τ as given in (2.30). (cid:3) Proof of Theorem 2.1.
In this section, we will first find the maximum dif- ference of (cid:107) u L ( · , t ) − v L ( · , t ) (cid:107) ∞ , then we will derive (cid:107) u L ( · , t ) − v L ( · , t ) (cid:107) H L,(cid:15),τ and (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ = (cid:107) U ( · , t ) − V ( · , t ) (cid:107) H L,(cid:15),τ . Combining these two, we will get an estimate for (cid:107) u ( · , t ) − v ( · , t ) (cid:107) H L,(cid:15),τ . Proposition 2.7. If u ( x ) satisfies (2.27), then (cid:107) u L − v L (cid:107) ∞ ≤ E (cid:15),τ ( t ) e − λL(cid:15) √ τ + E (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ where E (cid:15),τ ( t ) = | c ( · ) | ∞ + a τ e bτ t(cid:15)τ and E (cid:15),τ ( t ) = c τ t(cid:15)τ e ( bτ − t(cid:15)τ . Proof.
By the definition of u L and v L given in (2.19) and (2.22) and the assumptionthat u ( x ) = v ( x ) for x ∈ [0 , L ], we can get their difference u L ( x, t ) − v L ( x, t ) = c ( t ) (cid:16) e − x(cid:15) √ τ − φ ( x ) (cid:17) + (cid:16) U c ( L, t ) − h ( t ) + e − t(cid:15)τ h (0) (cid:17) φ ( x )Combining Lemmas 2.5(i), 2.5(ii), inequality (2.35), and h ( t ) ≡
0, we have (2.36) (cid:107) u L ( · , t ) − v L ( · , t ) (cid:107) ∞ ≤ E (cid:15),τ ( t ) e − λL(cid:15) √ τ + E (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ where E (cid:15),τ ( t ) = | c ( · ) | ∞ + a τ e bτ t(cid:15)τ and E (cid:15),τ ( t ) = c τ t(cid:15)τ e ( bτ − t(cid:15)τ . (2.37) (cid:3) Proposition 2.8. If u ( x ) satisfies (2.27), and E (cid:15),τ ( t ) , E (cid:15),τ ( t ) are as in propo-sition 2.7, then (cid:107) u L ( · , t ) − v L ( · , t ) (cid:107) H L,(cid:15),τ ≤ √ L (cid:18) E (cid:15),τ ( t ) e − λL(cid:15) √ τ + E (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ (cid:19) . Proof.
Because of the definition of u L and v L given in (2.19) and (2.22), Lemma (cid:107) ( u L ( · , t ) − v L ( · , t )) x (cid:107) ∞ ≤ | c ( t ) | e − L(cid:15) √ τ | φ (cid:48) ( x ) | + | U c ( L, t ) | | φ (cid:48) ( x ) |≤ (cid:15) √ τ (cid:18) E (cid:15),τ ( t ) e − λL(cid:15) √ τ + E (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ (cid:19) . (2.38)Now, combining (2.36) and (2.38), we obtain that (cid:107) u L ( · , t ) − v L ( · , t ) (cid:107) H L,(cid:15),τ = (cid:115)(cid:90) L | u L − v L | + (cid:12)(cid:12) (cid:15) √ τ ( u L − v L ) x (cid:12)(cid:12) dx ≤√ L (cid:18) E (cid:15),τ ( t ) e − λL(cid:15) √ τ + E (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ (cid:19) . (2.39) (cid:3) Proposition 2.9. If u ( x ) satisfies (2.27), then (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ ≤ γ (cid:15),τ ( t ) e − λL(cid:15) √ τ + γ (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ where the coefficients are given by γ (cid:15),τ ( t ) = e ( M +1)2 t M(cid:15) √ τ (cid:18) ( M + 1) √ τ M + 1 (cid:19) √ L (cid:18) t(cid:15)τ | c ( · ) | ∞ + a τ b τ ( e bτ t(cid:15)τ − (cid:19) γ (cid:15),τ ( t ) = e ( M +1)2 t M(cid:15) √ τ (cid:18) ( M + 1) √ τ M + 1 (cid:19) √ Lc τ · (cid:18) t(cid:15)τ ( b τ − e ( bτ − t(cid:15)τ − b τ − ( e ( bτ − t(cid:15)τ − (cid:19) . (2.40) Proof.
Multiplying the governing equation of W (2.25) by 2 W , integrating over [0 , L ], and using integration by parts, we get ddt (cid:90) L W + ( (cid:15) √ τ W x ) dx = − (cid:15) (cid:90) L W x dx + (cid:90) L W x ( f ( v ) − f ( u )) dx + 2 (cid:15)τ (cid:90) L W ( v L − u L ) dx. OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION13
Therefore, using the norm we defined earlier in (2.28), and f (cid:48) ( u ) ≤ ( M +1) M := C , we have ddt (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ ≤ (cid:90) L | W x || f (cid:48) ( η ) || v − u | dx + 2 √ L(cid:15)τ (cid:107) v L − u L (cid:107) ∞ (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ ≤ C (cid:90) L | W x | ( | W | + (cid:107) v L − u L (cid:107) ∞ ) dx + 2 √ L(cid:15)τ (cid:107) v L − u L (cid:107) ∞ (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ ≤ C(cid:15) √ τ (cid:16) (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ + (cid:107) v L − u L (cid:107) ∞ √ L (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ (cid:17) + 2 √ L(cid:15)τ (cid:107) v L − u L (cid:107) ∞ (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ = 2
C(cid:15) √ τ (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ + (cid:18) C(cid:15) √ τ + 2 (cid:15)τ (cid:19) √ L (cid:107) v L − u L (cid:107) ∞ (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ . Hence, ddt (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ ≤ C(cid:15) √ τ (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ + (cid:18) C(cid:15) √ τ + 1 (cid:15)τ (cid:19) √ L (cid:107) v L − u L (cid:107) ∞ . By Gronwall’s inequality and (2.36) (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ ≤ (cid:90) t (cid:18) C(cid:15) √ τ + 1 (cid:15)τ (cid:19) √ L (cid:107) v L − u L (cid:107) ∞ e C ( t − s ) (cid:15) √ τ ds ≤ e Ct(cid:15) √ τ (cid:18) C(cid:15) √ τ + 1 (cid:15)τ (cid:19) √ L (cid:90) t E (cid:15),τ ( s ) e − λL(cid:15) √ τ + E (cid:15),τ ( s ) e − λ ( L − L (cid:15) √ τ ds ≤ (cid:18) e Ct(cid:15) √ τ (cid:18) C(cid:15) √ τ + 1 (cid:15)τ (cid:19) √ L (cid:90) t E (cid:15),τ ( s ) ds (cid:19) e − λL(cid:15) √ τ + (cid:18) e Ct(cid:15) √ τ (cid:18) C(cid:15) √ τ + 1 (cid:15)τ (cid:19) √ L (cid:90) t E (cid:15),τ ( s ) ds (cid:19) e − λ ( L − L (cid:15) √ τ ≤ e Ct(cid:15) √ τ (cid:18) C(cid:15) √ τ + 1 (cid:15)τ (cid:19) √ L (cid:18) t | c ( · ) | ∞ + a τ (cid:15)τb τ ( e bτ t(cid:15)τ − (cid:19) e − λL(cid:15) √ τ + e Ct(cid:15) √ τ (cid:18) C(cid:15) √ τ + 1 (cid:15)τ (cid:19) √ L c τ (cid:15)τ (cid:18) (cid:15)τb τ − te ( bτ − t(cid:15)τ − ( (cid:15)τb τ − ( e ( bτ − t(cid:15)τ − (cid:19) e − λ ( L − L (cid:15) √ τ . Hence (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ ≤ γ (cid:15),τ ( t ) e − λL(cid:15) √ τ + γ (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ where γ (cid:15),τ ( t ) and γ (cid:15),τ ( t ) are given in (2.40). (cid:3) Now we are in the position to prove the main theorem of this section.
Theorem 2.10. If u ( x ) satisfies u ( x ) = (cid:26) C u x ∈ [0 , L ]0 x > L
04 YING WANG AND CHIU-YEN KAO where L < L and C u , are positive constants, and E (cid:15),τ ( t ) , E (cid:15),τ ( t ) , γ (cid:15),τ ( t ) , γ (cid:15),τ ( t ) are as in (2.37) and (2.40) , then (2.41) (cid:107) u ( · , t ) − v ( · , t ) (cid:107) H L,(cid:15),τ ≤ D (cid:15),τ ( t ) e − λL(cid:15) √ τ + D (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ for some < λ < , and D (cid:15),τ ( t ) = γ (cid:15),τ ( t ) + √ LE (cid:15),τ ( t ) , D (cid:15),τ ( t ) = γ (cid:15),τ ( t ) + √ LE (cid:15),τ ( t ) . Proof of the Main Theorem. (cid:107) u ( · , t ) − v ( · , t ) (cid:107) H L,(cid:15),τ ≤ (cid:107) W ( · , t ) (cid:107) H L,(cid:15),τ + (cid:107) v L ( · , t ) − u L ( · , t ) (cid:107) H L,(cid:15),τ = D (cid:15),τ ( t ) e − λL(cid:15) √ τ + D (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ where D (cid:15),τ ( t ) = γ (cid:15),τ ( t ) + √ LE (cid:15),τ ( t )= e ( M +1)2 t M(cid:15) √ τ (cid:18) ( M + 1) √ τ M + 1 (cid:19) √ L (cid:18) t(cid:15)τ | c ( · ) | ∞ + a τ b τ ( e bτ t(cid:15)τ − (cid:19) + √ L ( | c ( · ) | ∞ + a τ e bτ t(cid:15)τ ) ,D (cid:15),τ ( t ) = γ (cid:15),τ ( t ) + √ LE (cid:15),τ ( t )= e ( M +1)2 t M(cid:15) √ τ (cid:18) ( M + 1) √ τ M + 1 (cid:19) √ Lc τ ·· (cid:18) t(cid:15)τ ( b τ − e ( bτ − t(cid:15)τ − b τ − ( e ( bτ − t(cid:15)τ − (cid:19) + √ Lc τ t(cid:15)τ e ( bτ − t(cid:15)τ . (cid:3) This result gives that (cid:107) u ( · , t ) − v ( · , t ) (cid:107) H L,(cid:15),τ exponentially delays in L . This theorem shows that if λL(cid:15) √ τ and λ ( L − L ) (cid:15) √ τ converge to infinity, then the solution v ( x, t ) of the finite interval problem converges to the solution u ( x, t ) of the half line problem in the sense of (cid:107) · (cid:107) H L,(cid:15),τ . This can be achieved either by letting L → ∞ or (cid:15) →
0. For example, in the extreme case, (cid:15) = 0, the half line problem (1.11) becomes hyperbolic and the domain of dependence is finite, so, certainly, one only need to consider the finite interval problem. This is consistent with the main theorem in the sense that for a fixed final time t , if λL > b τ t and λ ( L − L ) > ( b τ − t , i.e., L > max( b τ tλ , ( b τ − tλ ), then (cid:107) u ( · , t ) − v ( · , t ) (cid:107) H L,(cid:15),τ ≤ D (cid:15),τ ( t ) e − λL(cid:15) √ τ + D (cid:15),τ ( t ) e − λ ( L − L (cid:15) √ τ → (cid:15) →
0. Theorem 2.10 gives a theoretical justification for using the solution of the finite interval problem to approximate the solution of the half line problem with appropriate choice of L and (cid:15) . Hence in the next chapter, the numerical scheme designed to solve the MBL equation (1.10) is given for finite interval problem. OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION15 Numerical schemes
To numerically solve the MBL equation (1.10), We first collect all the terms with time derivative and rewrite MBL equation (1.10) as (3.1) ( u − (cid:15) τ u xx ) t + ( f ( u )) x = (cid:15)u xx . By letting (3.2) w = u − (cid:15) τ u xx ⇐⇒ u = ( I − (cid:15) τ ∂ xx ) − w, MBL equation (3.1) can be written as w t + ( f ( u )) x = (cid:15)u xx . (3.3)Now, the new form of MBL equation (3.3) can be viewed as a PDE in terms of w , and the occurrence of u can be recovered by (3.2). Equation (3.3) can be formally viewed as w t + ( f (( I − (cid:15) τ ∂ xx ) − w )) x = (cid:15) (( I − (cid:15) τ ∂ xx ) − w ) xx , (3.4)which is a balance law in term of w . We adopt numerical schemes originally designed for hyperbolic equations to solve the MBL equation (3.1), which is of pseudo- parabolic type. The local discontinuous Galerkin method has been applied to solve equations involving mixed derivatives u xxt term [23, 24]. To the best knowledge of the authors, the central schemes have not been applied to solve equations of this kind. The main advantage of the central schemes is the simplicity. “the direction of the wind“ is not required to be identified, and hence the field-by-field decomposition can be avoided. In this chapter, we demonstrate how to apply the central schemes to solve the MBL equation (3.1). Second-order schemes.
In this section, we show how to apply the classical second order central schemes [18] originally designed for hyperbolic conservation laws to numerically solve the MBL equation (1.10), which is of pseudo-parabolic type. To solve (3.3), we modify the central scheme given in [18]. As in [18], at each time level, we first reconstruct a piecewise linear approximation of the form L j ( x, t ) = w j ( t ) + ( x − x j ) w (cid:48) j ∆ x , x j − ≤ x ≤ x j + . (3.5)Second-order accuracy is guaranteed if the so-called vector of numerical derivative w (cid:48) j ∆ x , which will be given later, satisfies w (cid:48) j ∆ x = ∂w ( x j , t ) ∂x + O (∆ x ) . (3.6)We denote the staggered piecewise-constant functions ¯ w j + ( t ) as ¯ w j + ( t ) = 1∆ x (cid:90) x j +1 x j w ( x, t ) dx. (3.7) Evolve the piecewise linear interplant (3.5) by integrating (3.3) over [ x j , x j +1 ] × [ t, t + ∆ t ] ¯ w j + ( t + ∆ t ) = ¯ w j + ( t ) − x (cid:34)(cid:90) t +∆ tt f ( u ( x j +1 , s )) ds − (cid:90) t +∆ tt f ( u ( x j , s )) ds (cid:35) + (cid:15) ∆ x (cid:34)(cid:90) t +∆ tt (cid:90) x j +1 x j ∂ u ( x, s ) ∂x dx ds (cid:35) . (3.8)We calculate each term on the right hand side of (3.8) below. For ¯ w j + ( t ), applying the definition of L j ( x, t ) and L j +1 ( x, t ) given in (3.5) to (3.7), we have that ¯ w j + ( t ) = 1∆ x (cid:90) x j + 12 x j L j ( x, t ) dx + 1∆ x (cid:90) x j +1 x j + 12 L j +1 ( x, t ) dx = 12 ( w j ( t ) + w j +1 ( t )) + 18 ( w (cid:48) j − w (cid:48) j +1 ) . (3.9)The middle two integrands can be approximated by the midpoint rule (cid:90) t +∆ tt f ( u ( x j , s )) ds = f ( u ( x j , t + ∆ t t + O (∆ t ) (cid:90) t +∆ tt f ( u ( x j +1 , s )) ds = f ( u ( x j +1 , t + ∆ t t + O (∆ t )(3.10)if the CFL condition λ · max x j ≤ x ≤ x j +1 (cid:12)(cid:12)(cid:12)(cid:12) ∂f ( u ( w ( x, t ))) ∂w (cid:12)(cid:12)(cid:12)(cid:12) < , where λ = ∆ t ∆ x is met. For MBL equation (3.3), we have that at t > u − (cid:15) τ u xx = w, u (0) = w (0) , u ( L ) = w ( L ) . Let v ( x ) = ( L − x ) w (0)+ xw ( L ) L , then u ( x ) = [ Iw ]( x ) = v ( x ) + 1 L (cid:90) L [ w ( y ) − v ( y )] K ( x, y ) dy where K ( x, y ) = ∞ (cid:88) k =1 sin( kπxL ) sin( kπyL )1 + ( kπL ) (cid:15) τ . Hence the eigenvalues for I are λ k = 11 + ( kπL ) (cid:15) τ ≤ , k = 1 , , . . . Therefore, the CFL condition is∆ t ∆ x · max x j ≤ x ≤ x j +1 (cid:12)(cid:12)(cid:12)(cid:12) ∂f ( u ( w ( x, t ))) ∂w (cid:12)(cid:12)(cid:12)(cid:12) = ∆ t ∆ x · max x j ≤ x ≤ x j +1 k =1 , , ... (cid:12)(cid:12)(cid:12)(cid:12) ∂f ( u ( x, t )) ∂u (cid:12)(cid:12)(cid:12)(cid:12) · λ k ≤ ∆ t ∆ x · . < OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION17
In the numerical computations in chapter 4, we chose ∆ t ∆ x = 0 .
1. In (3.10), to estimate u ( · , t + ∆ t )’s, we use Taylor expansion and the conservation law (3.3): w ( x j , t + ∆ t w j ( t ) + ∂w∂t ∆ t O (∆ t )= w j ( t ) + ( (cid:15) ∂ u∂x − ∂f∂x ) ∆ t O (∆ t )= w j ( t ) + ( (cid:15) ∆ x D u j − f (cid:48) j ) λ , (3.11)where D is the discrete central difference operator D u j = u j − − u j + u j +1 ∆ x , and the second-order accuracy is met if f (cid:48) j ∆ x = ∂f ( u ( x j , t )) ∂x + O (∆ x ) . (3.12)The choices for { w (cid:48) j } in (3.6) and { f (cid:48) j } in (3.12) can be found in [18], and we chose w (cid:48) j = M M { ∆ w j + , ∆ w j − } , f (cid:48) j = M M { ∆ f j + , ∆ f j − } (3.13)where M M { x, y } = minmod( x, y ) = (sgn( x )+sgn( y )) · Min( | x | , | y | ) and ∆ w j + = w j +1 − w j . Combining (3.8)-(3.10), we obtain ¯ w j + ( t + ∆ t ) = ¯ w j + ( t ) − λ [ f ( u j +1 ( t + ∆ t − f ( u j ( t + ∆ t (cid:15) ∆ x (cid:34)(cid:90) t +∆ tt (cid:90) x j +1 x j ∂ u ( x, s ) ∂x dx ds (cid:35) . (3.14)Next, we will re-write (3.14) in terms of u . ( u xx ) j + is approximated as( u xx ) j + = 1∆ x (cid:90) x j +1 x j u xx dx = 1∆ x ( u x ( x j +1 , t ) − u x ( x j , t )) , and using the cell averages, it becomes ( u xx ) j + = 1∆ x (cid:18) ¯ u j +3 / − ¯ u j +1 / ∆ x − ¯ u j +1 / − ¯ u j − / ∆ x (cid:19) = ¯ u j +3 / − u j +1 / + ¯ u j − / (∆ x ) = D ¯ u j + . (3.15)Notice that the linear interpolation (similar to (3.5))˜ L j + ( x, t + ∆ t ) = u j + ( t + ∆ t ) + ( x − x j + ) u (cid:48) j + ∆ x for x j ≤ x ≤ x j +1 and the cell average definition (similar to (3.7))¯ u j + ( t + ∆ t ) = 1∆ t (cid:90) x j +1 x j u ( x, t + ∆ t ) dx ensure that ¯ u j + ( t + ∆ t ) = u j + ( t + ∆ t ) , and the convertion between u and w is done using the following relation (3.16) ( I − (cid:15) τ D ) u = w. Hence re-writting (3.14) in terms of u gives the staggered central scheme ( I − (cid:15) τ D ) u j + ( t + ∆ t ) = ( I − (cid:15) τ D )¯ u j + ( t ) − λ [ f ( u j +1 ( t + ∆ t − f ( u j ( t + ∆ t (cid:15) ∆ x (cid:34)(cid:90) t +∆ tt (cid:90) x j +1 x j ∂ u ( x, s ) ∂x dx ds (cid:35) . (3.17)We will focus on the last integral in (3.17). There are many ways to numerically calculate this integral. We will show two ways to do this in the following two subsections, both of them achieve second order accuracy. Trapezoid Scheme.
In this scheme, we use the notion (3.7) and the trapezoid rule to calculate the integral numerically as follows: (cid:90) t +∆ tt (cid:90) x j +1 x j ∂ u ( x, s ) ∂x dx ds = ∆ x (cid:90) t +∆ tt ( u xx ) j + ( s ) ds = ∆ x ∆ t (cid:16) ( u xx ) j + ( t ) + ( u xx ) j + ( t + ∆ t )) (cid:17) (3.18)with O (∆ t ) error. Combining with (3.15) and (3.17), we can get the trapezoid scheme (cid:18) I − ( (cid:15) τ + (cid:15) ∆ t D (cid:19) u j + ( t + ∆ t ) = (cid:18) I − ( (cid:15) τ − (cid:15) ∆ t D (cid:19) ¯ u j + ( t ) − λ (cid:20) f ( u j +1 ( t + ∆ t − f ( u j ( t + ∆ t (cid:21) (3.19)The flow chart of the trapezoid scheme is given in (3.20) (3.20) ¯ w j + ( t ) (3.16) (cid:47) (cid:47) ¯ u j + ( t ) (3.19) (cid:43) (cid:43) (cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88) u j ( t ) (3.16) (cid:47) (cid:47) w j ( t ) (3.9) (cid:52) (cid:52) (cid:104)(cid:104)(cid:104)(cid:104)(cid:104)(cid:104)(cid:104) (3.11) (cid:42) (cid:42) (cid:86)(cid:86)(cid:86)(cid:86)(cid:86) u j + ( t + ∆ t ) w j ( t + ∆ t ) (3.16) (cid:47) (cid:47) u j ( t + ∆ t ) (3.19) (cid:51) (cid:51) (cid:102)(cid:102)(cid:102)(cid:102)(cid:102) Midpoint Scheme.
In this scheme, we use the notion (3.7) and the midpointrule to calculate the integral numerically as follows: (cid:90) t +∆ tt (cid:90) x j +1 x j ∂ u ( x, s ) ∂x dx ds = ∆ x (cid:90) t +∆ tt ( u xx ) j + ( s ) ds = ∆ x ∆ t ( u xx ) j + ( t + ∆ t ( I − (cid:15) τ D ) u j + ( t + ∆ t ) = ¯ w j + ( t ) − λ [ f ( u j +1 ( t + ∆ t − f ( u j ( t + ∆ t (cid:15) ∆ tD ¯ u j + ( t + ∆ t OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION19
The flow chart of the midpoint scheme is given in (3.22) (3.22) ¯ w j + ( t ) (3.21) (cid:43) (cid:43) (cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87)(cid:87) u j ( t ) (3.16) (cid:47) (cid:47) w j ( t ) (3.9) (cid:57) (cid:57) (cid:114)(cid:114)(cid:114)(cid:114)(cid:114)(cid:114)(cid:114)(cid:114)(cid:114)(cid:114) (3.11) (cid:37) (cid:37) (cid:76)(cid:76)(cid:76)(cid:76)(cid:76)(cid:76)(cid:76)(cid:76)(cid:76)(cid:76) ¯ w j + ( t + ∆ t ) (3.16) (cid:47) (cid:47) ¯ u j + ( t + ∆ t ) (3.21) (cid:47) (cid:47) u j + ( t + ∆ t ) w j ( t + ∆ t ) (3.9) (cid:79) (cid:79) (3.16) (cid:47) (cid:47) u j ( t + ∆ t ) (3.21) (cid:55) (cid:55) (cid:111)(cid:111)(cid:111)(cid:111)(cid:111)(cid:111)(cid:111)(cid:111)(cid:111)(cid:111)(cid:111)(cid:111) A third order semi-discrete scheme.
Similarly, we can extend the third order scheme to solve MBL equation (1.10), however, it is more involved. But the third order semi-discrete central scheme proposed in [12] can be extended to solve the MBL equation in a straightforward manner. In order to make the paper self-contained, we include the formulation below. d ¯ w j dt = − H j +1 / ( t ) − H j − / ( t )∆ x + (cid:15)Q j ( t )where ¯ w ( x, t ) denotes the cell average of w ¯ w j ( t ) = 1∆ x (cid:90) x j +1 / x j − / w ( x, t ) dx,H j +1 / ( t ) is the numerical convection flux and Q j ( t ) is a high-order approximation to the diffusion term u xx . H j +1 / ( t ) = f ( u + j +1 / ( t )) + f ( u − j +1 / ( t ))2 − a j +1 / ( t )2 (cid:104) w + j +1 / ( t ) − w − j +1 / ( t ) (cid:105) where u − j +1 / ( t ) , u + j +1 / ( t ) denote the left and right intermediate values of u ( x, t n )at x j +1 / , and their values are converted from the w − j +1 / ( t ) , w + j +1 / ( t ) using (3.2).The way to calculate w − j +1 / ( t ), w + j +1 / ( t ) and a j +1 / ( t ) is w + j +1 / ( t ) = A j +1 − ∆ x B j +1 + (∆ x ) C j +1 ,w − j +1 / ( t ) = A j + ∆ x B j + (∆ x ) C j ,a j +1 / ( t ) = max (cid:26) ∂f∂u ( u − j +1 / ( t )) , ∂f∂u ( u + j +1 / ( t )) (cid:27) , where A j = ¯ w nj − w C
12 ( ¯ w nj +1 − w nj + ¯ w nj − ) ,B j = 1∆ x (cid:20) w R ( ¯ w nj +1 − ¯ w nj ) + w C ¯ w nj +1 − ¯ w nj − w L ( ¯ w nj − ¯ w nj − ) (cid:21) ,C j = 2 w C ¯ w nj − − w nj + ¯ w nj +1 ∆ x ,w i = α i (cid:80) m α m α i = c i ( (cid:15) + IS i ) p , i, m ∈ { C, R, L } c L = c R = 1 / , c C = 1 / , (cid:15) = 10 − , p = 2 ,IS L = ( ¯ w nj − ¯ w nj − ) , IS R = ( ¯ w nj +1 − ¯ w nj ) ,IS C = 133 ( ¯ w nj +1 − w nj + ¯ w nj − ) + 14 ( ¯ w nj +1 − ¯ w nj − ) . The diffusion u xx is approximated using the following fourth-order central differ- encing form Q j ( t ) = − u j − + 16 u j − − u j + 16 u j +1 − u j +2 x . (3.23)The unique feature of this scheme is that the discretization is done in space first, and then the time evolution equation can be solved as a system of ordinary differential equations using any ODE solver of third order or higher. In this paper, we simply use the standard fourth order Runge-Kutta methods. Notice that to achieve the third order accuracy, the linear solver that converts u from w using (3.2) need also to be high order, and (3.23) is used to discretize u xx in our convertion. Computational results
In this section, we show the numerical solutions to the MBL equation u t + ( f ( u )) x = (cid:15)u xx + (cid:15) τ u xxt (4.1)with the initial condition u ( x ) = (cid:26) u B if x = 00 if x > Numerically, it is not practical to solve the half line problem (4.2), and one has to choose an appropriate computational domain. Theorem 2.10 in Chapter 2 provides a theoretical bound for the difference between the solution to the half line problem and that to the finite interval problem. However, the estimate (2.41) in
Theorem 2.10 includes time-dependent parameters D (cid:15),τ ( t ) and D (cid:15),τ ( t ), which cannot be obtained analyticaly. Therefore, we numerically demonstrate how the computational domain size affects the solution. We choose τ = 5, u B = α = (cid:113) and (cid:15) = 0 .
001 as an example here. Figure 4.1 shows the snapshot of the solutions at t = 0 . t = 0 . t = 1 for computational domain [0 , L ] with L = 0 . L = 0 . and L = 1 . In Figure 4.1(a), t = 0 .
1, the leading shock is located at f (¯ u τ =5 )¯ u τ =5 × . . × . . L = 0 . L = 0 . L = 1 .
25 all exceed the leading shock location.
Hence all the three computational domains deliver visually indistinguishable results.
Whereas, in Figure 4.1(b), t = 0 .
5, the leading shock is located at 1 . × . . OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION21 (a) t = 0 . (b) t = 0 . (c) t = 1 Figure 4.1.
Numerical solutions of MBL (4.1) at (a) t = 0 .
1, (b) t = 0 .
5, (c) t = 1 using the trapezoid scheme (3.20). ‘—’, ‘—’, ‘—’ denote the numerical solutions corresponding to computationaldomain [0 , L ] with L = 0 . L = 0 .
75 and L = 1 .
25 respectively.The parameter values are τ = 5, u B = α = (cid:113) , (cid:15) = 0 . x = (cid:15) , ∆ t = 0 . x . L = 0 .
25 is shorter than the computational domain needed to capture this shock, hence the numerical solution halts at x = 0 .
25. On the contrast, L = 0 .
75 and L = .
25 are both large enough to capture this shock front. Similarly, in Figure 4.1(c), t = 1, the leading shock is located at 1.02. L = 0 . < .
02 and L = 0 . < . both result in wrong solution profiles. More specifically, both solutions halt at the boundary of the insufficient computational domain. But L = 1 . > .
02 is large enough to capture the correct solution profile.
In the rest of this chapter, all the computational domains [0 , L ] are thereforechosen based on the principle:
L > leading shock speed × computational time . In addition, numerical solutions for larger L ’s, for example, L = 1 . L = 2 . L = 5, L = 10 are also sought. For all these larger L ’s, the numerical solutions are all consistent with that corresponding to L = 1 .
25 up to t = 1. This confirms that it is not necessary to take L too much larger than leading shock speed × computational time . To validate the order analysis given in chapter 3 for various schemes proposed, we first test the order of our schemes numerically with a smooth initial condition u ( x ) = u B H ( x − , , where H ( x, ξ ) = x < − ξ − (1 + xξ + π sin( πxξ )) if − ξ ≤ x ≤ ξ x > ξ . The final time T = 1 was employed, so that there was no shock created. (cid:15) in the MBL equation (4.1) is taken to be 1, M is taken to be 2, and the computational interval is [ − , L , L , L ∞ order tests of the trapezoid scheme and the third order semi-discrete scheme with different parameter τ value and the initial condition u B are given in Tables 4.1, 4.2. Table 4.1 shows that the trapezoid rule achieved second order accuracy for all the tested cases in L , L , L ∞ sense. Table N (cid:119)(cid:119)(cid:119) u ∆ x − u ∆ x (cid:119)(cid:119)(cid:119) order (cid:119)(cid:119)(cid:119) u ∆ x − u ∆ x (cid:119)(cid:119)(cid:119) order (cid:119)(cid:119)(cid:119) u ∆ x − u ∆ x (cid:119)(cid:119)(cid:119) ∞ order60 7.5416e-03 - 2.5388e-03 - 1.5960e-03 - u B = 0 . τ = 0 . u B = 0 . τ = 1 240 5.5697e-04 1.9488 1.8259e-04 1.9480 1.1283e-04 1.9038480 1.4104e-04 1.9815 4.6109e-05 1.9855 2.8719e-05 1.974060 1.3102e-02 - 4.1784e-03 - 2.2411e-03 - u B = 0 . τ = 5 240 9.6737e-04 1.9039 2.8089e-04 1.9686 1.5667e-04 1.9625480 2.5825e-04 1.9053 7.1250e-05 1.9790 3.9286e-05 1.995660 6.4427e-03 - 2.1578e-03 - 1.1682e-03 - u B = α
120 1.6611e-03 1.9555 5.7775e-04 1.9011 3.6447e-04 1.6804 τ = 0 . u B = α
120 2.0069e-03 1.9185 6.4998e-04 1.8906 3.7650e-04 1.8277 τ = 1 240 5.1832e-04 1.9531 1.6801e-04 1.9519 1.0062e-04 1.9037480 1.3136e-04 1.9803 4.2497e-05 1.9831 2.5599e-05 1.974860 1.1959e-02 - 3.8026e-03 - 1.9938e-03 - u B = α
120 3.2940e-03 1.8602 9.9527e-04 1.9338 5.4231e-04 1.8783 τ = 5 240 8.7736e-04 1.9086 2.5358e-04 1.9727 1.3933e-04 1.9606480 2.3271e-04 1.9146 6.4252e-05 1.9806 3.4967e-05 1.994460 5.7714e-03 - 1.9358e-03 - 1.0481e-03 - u B = 0 .
75 120 1.5035e-03 1.9406 5.1617e-04 1.9070 2.8061e-04 1.9011 τ = 0 . u B = 0 .
75 120 1.8963e-03 1.9213 6.1315e-04 1.8974 3.4013e-03 1.8272 τ = 1 240 4.8284e-04 1.9736 1.5796e-04 1.9567 9.0912e-04 1.9035480 1.2093e-04 1.9974 3.9783e-05 1.9894 2.3121e-05 1.975360 1.1042e-02 - 3.5020e-03 - 1.8299e-03 - u B = 0 .
75 120 3.0287e-03 1.8662 9.1181e-04 1.9414 4.8976e-04 1.9016 τ = 5 240 8.0111e-04 1.9186 2.3118e-04 1.9797 1.2593e-04 1.9595480 2.1076e-04 1.9264 5.8358e-05 1.9860 3.1627e-05 1.9934 Table 4.1.
The accuracy test for the trapezoid scheme for theMBL equation (4.1) with (cid:15) = 1 and M = 2.4.2 shows that the semi-discrete scheme has the order of accuracy greater than 2.5 for all the cases, and exceeds 3 for some cases. This confirms the accuracy study given in sections 3.1.1 and 3.2 respectively. We will now use examples to study the solutions to MBL equation (4.1) using the numerical schemes proposed in chapter 3. We first notice that if we scale t and x as follows ˜ t = t(cid:15) , ˜ x = x(cid:15) , then MBL (4.1) equation can be written in terms of ˜ t and ˜ x as follows u ˜ t + ( f ( u )) ˜ x = u ˜ x ˜ x + τ u ˜ x ˜ x ˜ t . (4.3) OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION23 N (cid:119)(cid:119)(cid:119) u ∆ x − u ∆ x (cid:119)(cid:119)(cid:119) order (cid:119)(cid:119)(cid:119) u ∆ x − u ∆ x (cid:119)(cid:119)(cid:119) order (cid:119)(cid:119)(cid:119) u ∆ x − u ∆ x (cid:119)(cid:119)(cid:119) ∞ order120 2.6992e-03 - 1.1300e-03 - 7.2363e-04 - u B = 0 . τ = 0 . u B = 0 . τ = 1 480 1.2006e-04 2.8606 5.0480e-05 2.8690 4.1985e-05 2.8671960 1.5942e-05 2.9129 6.6663e-06 2.9208 5.1464e-06 3.0282120 3.7573e-03 - 1.2122e-03 - 7.9211e-04 - u B = 0 . τ = 5 480 1.1994e-04 2.6373 3.8434e-05 2.6524 2.5089e-05 2.5857960 1.5565e-05 2.9460 4.9190e-06 2.9660 3.1363e-06 2.9999120 2.1836e-03 - 9.1039e-04 - 5.7219e-04 - u B = α
240 3.2729e-04 2.7381 1.3760e-04 2.7260 8.9550e-05 2.6757 τ = 0 . u B = α
240 7.0517e-04 2.4680 2.9669e-04 2.4656 2.4272e-04 2.5149 τ = 1 480 9.6528e-05 2.8690 4.0354e-05 2.8781 3.3125e-05 2.8733960 1.2890e-05 2.9047 5.3648e-06 2.9111 4.0754e-06 3.0229120 3.0797e-03 - 9.9202e-04 - 6.4456e-04 - u B = α
240 6.1133e-04 2.3328 1.9783e-04 2.3261 1.2277e-04 2.3924 τ = 5 480 9.7351e-05 2.6507 3.1222e-05 2.6637 2.0263e-05 2.5990960 1.2396e-05 2.9733 3.9513e-06 2.9822 2.4962e-06 3.0210120 1.8244e-03 - 7.5548e-04 - 4.6671e-04 - u B = 0 .
75 240 2.7262e-04 2.7425 1.1419e-04 2.7260 7.3299e-05 2.6707 τ = 0 . u B = 0 .
75 240 5.8671e-04 2.4798 2.4585e-04 2.4754 1.9866e-04 2.5304 τ = 1 480 7.9974e-05 2.8750 3.3285e-05 2.8848 2.7033e-05 2.8775960 1.0724e-05 2.8987 4.4466e-06 2.9041 3.3341e-06 3.0193120 2.5902e-03 - 8.3335e-04 - 5.3882e-04 - u B = 0 .
75 240 5.1342e-04 2.3348 1.6611e-04 2.3268 1.0271e-04 2.3913 τ = 5 480 8.1062e-05 2.6630 2.6032e-05 2.6738 1.6813e-05 2.6109960 1.0173e-05 2.9944 3.2662e-06 2.9946 2.0473e-06 3.0377 Table 4.2.
The accuracy test for the third order semi-discretescheme for the MBL equation (4.1) with (cid:15) = 1 and M = 2.The scaled equation (4.3) shows that it is the magnitude of t(cid:15) and x(cid:15) that determine the asymptotic behavior, not t , x , neither (cid:15) alone ([21]). In addition, (4.3) also shows that the dispersive parameter τ denotes the relative importance of the dispersive term u xxt . The bigger τ is, the more dispersive effect (4.1) equation has. This can be seen from the computational results to be shown later in this section. Duijn et al. [21] numerically provided a bifurcation diagram (Figure 4.2) of MBL (4.1) equation as the dispersive parameter τ and the post-shock value u B of the initial condition vary. The solution of (4.1) has been proven to display qualitatively different profiles for parameter values ( τ, u B ) falling in different regimes of the bifurcation diagram. In particular, for every fixed τ value, there are two critical u B τ bifurcation diagram τ * αβ u( τ )u( τ ) Figure 4.2.
The bifurcation diagram of the MBL equation (1.10)with the bifurcation parameters ( τ, u B ). u B values, namely, ¯ u and u . From the bifurcation diagram (Figure 4.2), it is clear that, when τ < τ ∗ , ¯ u = u = α . For a fixed τ value, the solution has three different profiles. (a) If u B ∈ [¯ u, u B for 0 ≤ xt ≤ dfdu ( u B ), a rarefaction wave connection u B to ¯ u for dfdu ( u B ) ≤ xt ≤ dfdu (¯ u ), another plateau value ¯ u for dfdu (¯ u ) < xt < f (¯ u )¯ u , and a shock from ¯ u down to 0 at xt = f (¯ u )¯ u (see Figure 3(a)). (b) If u B ∈ ( u, ¯ u ), the solution contains a plateau value u B for 0 ≤ xt < f (¯ u ) − f ( u B )¯ u − u B , a shock from u B up to ¯ u at xt = f (¯ u ) − f ( u B )¯ u − u B , another plateau value ¯ u for f (¯ u ) − f ( u B )¯ u − u B < xt < f (¯ u )¯ u , and a shock from ¯ u down to 0 at xt = f (¯ u )¯ u (see Figure 3(b)). The solution may exhibit a damped oscillation near u = u B . (c) If u B ∈ (0 , u ], the solution consists a single shock connecting u B and 0 at xt = f ( u B ) u B (see Figure 3(c)). It may exhibit oscillatory behavior near u = u B . Notice that when τ > τ ∗ and u < u B < ¯ u , the solution profiles (3(b)) displays non-monotonicity, which is consistent with the experimental observations ([7]). In the numerical computation we show below, we will therefore test the accuracy and capability of central schemes for different parameter values ( τ and u B ) that fall into various regimes of the bifurcation diagram, and therefore display qualitatively different solution profiles. The numerical experiments were carried out for M = 2, (cid:15) = 0 .
001 and T = 4000 × (cid:15) , i.e. ˜ T = 4000 to get the asymptotic solution profiles, and ∆ x was chosen to be (cid:15) and λ = ∆ t ∆ x was chosen to be 0.1. The scheme used in the computation is the second order Trapezoid scheme as shown in section 3.1.1. The Midpoint scheme delivers similar computational results, hence is omitted here.
The solution profiles at T (blue), ∗ T (green), ∗ T (magenta) and T (black) are OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION25 xt u u B ¯ uu dfdu ( u B ) dfdu ( ¯ u ) f ( ¯ u )¯ u (a) xt u u B u ¯ uf ( ¯ u ) − f ( u B )¯ u − u B f ( ¯ u )¯ u (b) u B xt u u f ( u B ) u B (c) Figure 4.3.
Given a fixed τ , the three qualitatively different so-lution profiles due to different values of u B . In particular, when τ > τ ∗ and u < u B < ¯ u , the solution profiles (Figure 3(b)) dis-plays non-monotonicity, which is consistent with the experimentalobservations ([7]). Figures 3(a), 3(b) and 3(c) are demonstrativefigures.chosen to demonstrate the time evolution of the solutions. The red dashed lines are used to denote the theoretical shock locations and plateau values for comparison purpose. We start with τ >
0. Based on the bifurcation diagram (Figure 4.2), we choose three representative u B values, i.e. u B = 0 . > α , u B = α = (cid:113) MM +1 = (cid:113) (for M = 2) and u B = 0 . < α . For each fixed u B , we choose three representative τ values, i.e. τ = 0 . < τ ∗ ≈ . τ = 1 > τ ∗ with u B = 0 . < u τ =1 < u B = α < ¯ u < u B = 0 .
9, and τ = 5 with u B = 0 . , α, . ∈ [ u τ =5 , ¯ u τ =5 ]. We first use this 9 pairs of ( τ, u B ) values given in Table 4.3 to validate the solution profiles with the demonstrative solution profiles given in Figure 4.3.( τ, u B ) Example 4 Example 5 Example 6Example 1 (0 . , .
9) (1 , .
9) (5 , . . , α ) (1 , α ) (5 , α )Example 3 (0 . , .
75) (1 , .
75) (5 , . Table 4.3. τ, u B ) values with either fixed τ value orfixed u B value used in Examples 1 – 6. Example 1 ( τ, u B ) = (0 . , . , ( τ, u B ) = (1 , . , ( τ, u B ) = (5 , . When u B = 0 . > α is fixed, we increase τ from 0.2 to 1 to 5 (Figure 4(a) , 4(b) , 4(c)), the dispersive effect starts to dominate the solution profile. When τ = 0 . (Figure 4(a)), the solution profile is similar to the classical BL equation solution (see Figure 2(b)), with a rarefaction wave for xt ∈ [ f (cid:48) ( u = 0 . , f (cid:48) ( u = α ) = f (cid:48) ( u = ¯ u τ =0 . )] and a shock from u = α to u = 0 at xt = f (cid:48) ( α ). This corresponds to Figure 3(a) with dfdu (¯ u τ =0 . = α ) = f (¯ u τ =0 . )¯ u τ =0 . = f ( α ) α . When τ = 1 (Figure 4(b)), the rarefaction wave is between xt ∈ [ f (cid:48) ( u = 0 . , f (cid:48) ( u = ¯ u τ =1 )] and the solution remains at the plateau value u = ¯ u τ =1 for xt ∈ [ f (cid:48) ( u = ¯ u τ =1 ) , f (¯ u τ =1 )¯ u τ =1 ] and the shock occurs at xt = f (¯ u τ =1 )¯ u τ =1 . This corresponds to Figure 3(a) with u B = 0 . > ¯ u τ =1 ≈ .
86. When τ = 5 (Figure 4(c)), the solution displays the first shock from u = 0 . to u = ¯ u τ =5 at xt = f (¯ u τ =5 ) − f ( u B )¯ u τ =5 − u B , and then remains at the plateau value u = ¯ u τ =5 for xt ∈ [ f (¯ u τ =5 ) − f ( u B )¯ u τ =5 − u B , f (¯ u τ =5) ¯ u τ =5 ] and the second shocks occurs at xt = f (¯ u τ =5) ¯ u τ =5 . This corresponds to Figure 3(b) with u τ =5 ≈ . < u B = 0 . < ¯ u τ =5 ≈ .
98. Notice that as τ increases, the rarefaction region shrinks and the plateau region enlarges. (a) ( τ, u B ) = (0 . , . u (b) ( τ, u B ) = (1 , . u (c) ( τ, u B ) = (5 , . u (d) ( τ, u B ) = (0 . , α ) u (e) ( τ, u B ) = (1 , α ) u (f) ( τ, u B ) = (5 , α ) u (g) ( τ, u B ) = (0 . , . u (h) ( τ, u B ) = (1 , . u (i) ( τ, u B ) = (5 , . u Figure 4.4.
Numerical solutions to MBL equation with param-eter settings fall in different regimes of the bifurcation diagram(Figure 4.2). The color coding is for different time: T (blue), T (green), T (magenta) and T (black). The results are discussed inexamples 1 – 6. In figures 4(d) – 4(f), α = (cid:113) MM +1 = (cid:113) for M = 2. Example 2 ( τ, u B ) = (0 . , α ) , ( τ, u B ) = (1 , α ) , ( τ, u B ) = (5 , α ). When u B = α is fixed, we increase τ from 0.2 to 1 to 5 (Figure 4(d) , 4(e) , 4(f)), the dispersive effect starts to dominate the solution profile. When τ = 0 .
2, the solution displays one single shock at xt = f ( α ) α . For both τ = 1 and τ = 5, the OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION27 solution has two shocks, one at xt = f (¯ u τ =1( τ =5 respectively) ) − f ( α )¯ u τ =1( τ =5 respectively) − α , and another one at xt = f (¯ u τ =1( τ =5 respectively) )¯ u τ =1( τ =5 respectively) . For both τ = 1 and τ = 5 (Figures 4(e) 4(f)), the solutions correspond to Figure 3(b), which are consistent to the experimental obser- vations. Notice that as τ increases from 1 to 5, i.e., the dispersive effect increases, the inter-shock interval length increases at every fixed time (compare Figure 4(e) with Figure 4(f)). In addition, for fix τ = 1 ( τ = 5 respectively), as time progresses, the inter-shock interval length increases in the linear fashion (see Figure 4(e) (Fig- ure 4(f) respectively) ). Example 3 ( τ, u B ) = (0 . , . , ( τ, u B ) = (1 , . , ( τ, u B ) = (5 , . When u B = 0 . < = α is fixed, we increase τ from 0.2 to 1 to 5 (Figure 4(g) , 4(h) , 4(i)), the dispersive effects starts to dominate the solution profile in the similar fashion as u B = 0 . u B = α . Notice that when τ = 1, since u B = 0 .
75 is very close to u τ =1 , the solution displays oscillation at xt = f ( u B ) u B (Figure 4(h)). If we increase τ further to τ = 5, the dispersive effect is strong enough to create a plateau value at ¯ u ≈ .
98 (see Figure 4(i)).
Example 4 ( τ, u B ) = (0 . , . , ( τ, u B ) = (0 . , α ) , ( τ, u B ) = (0 . , . Now, we fix τ = 0 .
2, decrease u B from 0.9 to α , to 0.75 (Figures4(a) 4(d) 4(g)). If u B > α the solution consists a rarefaction wave connecting u B down to α , then a shock from α to 0, otherwise, the solution consists a single shock from u B down to
0. In all cases, since τ = 0 . < τ ∗ , regardless of the u B value, the solution will not display non-monotone behavior, due to the lack of dispersive effect. Example 5 ( τ, u B ) = (1 , . , ( τ, u B ) = (1 , α ) , ( τ, u B ) = (1 , . Now, we fix τ = 1, decrease u B from 0.9 to α , to 0.75 (Figures4(b) 4(e) 4(h)). If u B = 0 . > ¯ u τ =1 , the solution consists a rarefaction wave connecting u B and ¯ u , and a shock connecting ¯ u down to 0 (Figure 4(b)). Even if u < u B < ¯ u , because τ = 1 > τ ∗ , the solution still has a chance to increase to the plateau value ¯ u as seen in Figure 4(e). But, if u B is too small, for example, u B = 0 . < u , the solution does not increase to ¯ u any more, instead, it consists a single shock connecting u B down to 0 (Figure 4(h)). Example 6 ( τ, u B ) = (5 , . , ( τ, u B ) = (5 , α ) , ( τ, u B ) = (5 , . Now, we fix τ = 5, decrease u B from 0.9 to α , to 0.75 (Figures4(c) 4(f) 4(i)). For all three u B , they are between u τ =5 and ¯ u τ =5 , hence all increase to the plateau value ¯ u τ =5 ≈ .
98 before dropping to 0. Notice that as u B decreases, the inter-shock interval length decreases at every fixed time (compare Figures 4(c), 4(f) and 4(i)). This shows that when the dispersive effect is strong ( τ > τ ∗ ), the bigger u B is, the bigger region the solution stays at the plateau value. Example 7 ( τ, u B ) = (0 , . , ( τ, u B ) = (0 , α ) , ( τ, u B ) = (0 , . We now show the solution profiles for the extreme τ value, i.e. τ = 0 in Figures u B = 0 . u B = α ) and 5(c) ( u B = 0 . of classical BL equation with small diffusion (cid:15)u xx . We compare Figures 5(a), 5(b) and 5(c) with the solution of the classical BL equation given in Figures 2(a) and is that due to the diffusion term in the MBL equation, as shown in Figure 4.5, the solutions do not have sharp edges right at the shock, instead, the solutions smear out a little. If we compare Figures 5(a), 5(b) and 5(c) with Figures 4(a), 4(d) and τ < τ ∗ , solution profile will stay the same for a fixed u B value. (a) ( τ, u B ) = (0 , . u (b) ( τ, u B ) = (0 , α ) u (c) ( τ, u B ) = (0 , . u Figure 4.5.
The numerical solutions of the MBL equation at T= 1 with τ = 0 and different u B values. The results are discussedin example 7. Example 8 ( τ, u B ) = (5 , . , ( τ, u B ) = (5 , . , ( τ, u B ) = (5 , . We also study the solution profiles for u B close to ¯ u . For example, when τ = 5, ¯ u ≈ .
98, we hence choose u B = 0 . u B = 0 . u B = 0 .
97 and solutions are shown in Figure 6(a), 6(b), 6(c). If u B = 0 . > ¯ u τ =5 ≈ .
98, the solution drops to the plateau value ¯ u , then drops to 0 (see Figure 6(a)). If u B = 0 . ≈ ¯ u τ =5 , the solution remains at plateau value ¯ u τ =5 and then drop to 0 (see Figure 6(b)). If u B = 0 . < ¯ u τ =5 , the solution increases to the plateau value ¯ u τ =5 ≈ .
98, then drops to 0. In all cases, the transition from u B to ¯ u τ =5 ≈ .
98 takes very small space. In the majority space, the solution keeps to be the plateau value ¯ u τ =5 ≈ . (a) ( τ, u B ) = (5 , . u (b) ( τ, u B ) = (5 , . u (c) ( τ, u B ) = (5 , . u Figure 4.6.
Numerical solutions to MBL equation with u B closeto ¯ u τ =5 ≈ .
98. The color coding is for different time: T (blue), T (green), T (magenta) and T (black). The results are discussedin example 8. Example 9 ( τ, u B ) = (5 , . τ, u B ) = (5 , . τ, u B ) = (5 , . τ, u B ) = (5 , . τ, u B ) = (5 , . OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION29
In addition, we study the solution profiles for u B close to u . For example, when τ = 5, u ≈ .
68, we hence choose u B = 0 . u B = 0 . u B = 0 . u B = 0 . u B = 0 .
66 and solutions are shown in Figures 7(a), 7(b), 7(c), 7(d), 7(e). As u B decreases crossing u τ =5 ≈ .
68, the solution gradually stops increasing to the plateau value ¯ u τ =5 , and the inter-shock interval length decreases (compare Figures that u B values are too close to u τ =5 . This confirms that even with big dispersive effect (say τ = 5), if u B is too small (e.g. u B < u ), the solution will not exhibit non-monotone behavior. (a) ( τ, u B ) = (5 , . u (b) ( τ, u B ) = (5 , . u (c) ( τ, u B ) = (5 , . u (d) ( τ, u B ) = (5 , . u (e) ( τ, u B ) = (5 , . u Figure 4.7.
Numerical solutions to MBL equation with u B closeto u τ =5 ≈ .
68. The color coding is for different time: T (blue), T (green), T (magenta) and T (black). The results are discussedin example 9. Example 10 ( τ, u B ) = (0 . , . τ, u B ) = (1 , . τ, u B ) = (5 , . We fix u B to be small, and in this example, we take it to be u B = 0 .
6. We vary the τ value, from τ = 0 . < τ ∗ to τ = 1 barely larger than τ ∗ to τ = 5 > τ ∗ . The numerical solutions are given in Figure 8(a), 8(b), 8(c). As τ increases, the post-shock value remains the same, but there will be oscillation generated as τ becomes larger than τ ∗ . Figures 8(d), 8(e) and 8(f) show that as τ increases, the oscillation amplitude increases and oscillates more rounds. Notice that τ is the dispersive parameter, and this means that even for small u B value, different dispersive parameter values still give different dispersive effects, although none can bring the solution to the plateau value ¯ u . Comparing Figures 8(d), 8(e) and 8(f) with Figures 8(g), 8(h) and 8(i), it is clear that the oscillation amplitude remains steady with respect to time. Example 11 (cid:15) = 0 . , (cid:15) = 0 . , (cid:15) = 0 . , (cid:15) = 0 . , (cid:15) = 0 . In this example, we will compare the solution profiles for different (cid:15) values. Fixing (a) ( τ, u B ) = (0 . , . u (b) ( τ, u B ) = (1 , . u (c) ( τ, u B ) = (5 , . u (d) Fig 8(a) zoomed in at T u (e) Fig 8(b) zoomed in at T u (f) Fig 8(c) zoomed in at T u (g) Fig 8(a) zoomed in at T u (h) Fig 8(b) zoomed in at T u (i) Fig 8(c) zoomed in at T u Figure 4.8.
Numerical solutions to MBL equation with small con-stant u B = 0 . τ values. The figures on the secondand third rows are the magnified versions of the first row at t = T and t = T respectively. The color coding is for different time: T (blue), T (green), T (magenta) and T (black). The results arediscussed in examples 10. T = 0 . , ∆ x = 0 . , λ = ∆ t ∆ x = 0 .
1, we show the numerical results in Figure 4.9 for (cid:15) = 0 .
001 (blue), (cid:15) = 0 .
002 (yellow), (cid:15) = 0 .
003 (magenta), (cid:15) = 0 .
004 (green), and (cid:15) = 0 .
005 (black). For the purpose of cross reference, we choose the same nine sets of parameter settings as in examples 1– 6. To assist the observation, the figures in Figure 4.9 are zoomed into the regions where different (cid:15) values introduce different solution profiles. The numerical solutions clearly show that as (cid:15) increases, the numerical solution is smeared out, and the jump location becomes less accurate.
Notice that τ is responsible for the competition between the diffusion and disper- sion, which in turn determines the plateau values. Hence varying (cid:15) value doesn’t affect the plateau location. OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION31 (a) ( τ, u B ) = (0 . , . u (b) ( τ, u B ) = (1 , . u (c) ( τ, u B ) = (5 , . u (d) ( τ, u B ) = (0 . , α ) u (e) ( τ, u B ) = (1 , α ) u (f) ( τ, u B ) = (5 , α ) u (g) ( τ, u B ) = (0 . , . u (h) ( τ, u B ) = (1 , . u (i) ( τ, u B ) = (5 , . u Figure 4.9.
The numerical solutions of MBL equation at T = 0 . (cid:15) = 0 .
001 (blue), (cid:15) = 0 .
002 (yellow), (cid:15) = 0 .
003 (magenta), (cid:15) = 0 .
004 (green), and (cid:15) = 0 .
005 (black). The view windows arezoomed into the regions where different (cid:15) values impose differentsolution profiles. The results are discussed in example 11.5.
Conclusion
We proved that the solution to the infinite domain problem can be approximated by that of the bounded domain problem. This provides a theoretical justification for using finite domain to calculation the numerical solution of the MBL equation (1.10). We also extended the classical central scheme originally designed for the hyperbolic systems to solve the MBL equation, which is of pseudo-parabolic type.
The numerical solutions for qualitatively different parameter values τ and initial conditions u B show that the jump locations are consistent with the theoretical calculation and the plateau heights are consistent with the numerically obtained values given in [21]. In particular, when τ > τ ∗ , for u B ∈ ( u, ¯ u ), the numerical solutions give non-monotone water saturation profiles, which is consistent with the experimental observations. In addition, the order tests show that the proposed second and third order central schemes achieved the desired accuracies. In [22, 20], the two-dimensional space extension of the modified Buckley-Leverett equation has been derived. One of the future directions is to develop high order numerical schemes to solve the two-dimensional MBL equation. Central schemes have been used to solve high dimensional hyperbolic problem and dispersive prob- lem ([11, 16]), which makes it a good candidate for such a task.
Appendix A. Proof of the lemmas
Proof to lemma 2.2.
Let g ( u ) = f ( u ) u = uu + M (1 − u ) , then g (cid:48) ( u ) = M − (1 + M ) u ( u + M (1 − u ) ) > < u < (cid:113) MM +1 = 0 if u = (cid:113) MM +1 < u > (cid:113) MM +1 and hence g ( u ) achieves its maximum at u = (cid:113) MM +1 . Therefore, f ( u ) u = g ( u ) ≤ D , where D = f ( α ) α and α = (cid:113) MM +1 , and in turn, we have that f ( u ) ≤ Du for all ≤ u ≤ (cid:3) Proof to lemma 2.3 (i). (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − λξ(cid:15) √ τ dξ = (cid:15) √ τ − e ( λ − x(cid:15) √ τ λ − ≤ (cid:15) √ τ − λ if λ ∈ (0 , . (cid:3) Proof to lemma 2.3 (ii). (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − ξ(cid:15) √ τ dξ = xe ( λ − x(cid:15) √ τ ≤ (cid:15) √ τe (1 − λ ) if λ ∈ (0 , . (cid:3) Proof to lemma 2.3 (iii).
Based on the assumption on u in (2.27) (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx(cid:15) √ τ | u ( ξ ) | dξ ≤ (cid:90) + ∞ e − | x − ξ | (cid:15) √ τ e λx(cid:15) √ τ | u ( ξ ) | dξ ≤ C u e λx(cid:15) √ τ (cid:90) L e − | x − ξ | (cid:15) √ τ dξ = C u y ( x )(A.1)Calculating y ( x ) with the assumption that λ ∈ (0 , y ( x ) = e λx(cid:15) √ τ (cid:82) L e − | x − ξ | (cid:15) √ τ dξ ≤ (cid:15) √ τ e λx(cid:15) √ τ ≤ (cid:15) √ τ e λL (cid:15) √ τ for x ∈ [0 , L ] e ( λ − x(cid:15) √ τ (cid:82) L e ξ(cid:15) √ τ dξ ≤ (cid:15) √ τ e ( λ − x + L (cid:15) √ τ ≤ (cid:15) √ τ e λL (cid:15) √ τ for x ∈ [ L , + ∞ )Therefore, we get the desired inequality (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ − e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx(cid:15) √ τ | u ( ξ ) | dξ ≤ C u (cid:15) √ τ e λL (cid:15) √ τ . (cid:3) OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION33
Proof to lemma 2.4 (i). (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − λξ(cid:15) √ τ dξ = (cid:15) √ τλ − (cid:16) − λe ( λ − x(cid:15) √ τ − λ − e − x(cid:15) √ τ (cid:17) ≤ (cid:15) √ τ − λ if λ ∈ (0 , . (cid:3) Proof to lemma 2.4 (ii). (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx − ξ(cid:15) √ τ dξ = 2 e ( λ − x(cid:15) √ τ − e ( λ − x(cid:15) √ τ − (cid:15) √ τ + xe ( λ − x(cid:15) √ τ ≤ (cid:15) √ τ + (cid:15) √ τe (1 − λ ) if λ ∈ (0 , . (cid:3) Proof to lemma 2.4 (iii).
Based on the assumption on u in (2.27) (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx(cid:15) √ τ | u ( ξ ) | dξ ≤ C u e λx(cid:15) √ τ (cid:90) L (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) dξ = C u y ( x )(A.2)Calculating y ( x ) with the assumption that λ ∈ (0 , x ∈ [0 , L ] y ( x ) ≤ e ( λ − x(cid:15) √ τ (cid:90) x ( e − ξ(cid:15) √ τ + e ξ(cid:15) √ τ ) dξ + e ( λ +1) x(cid:15) √ τ (cid:90) L x e − ξ(cid:15) √ τ dξ ≤ (cid:15) √ τ e λL (cid:15) √ τ and y ( x ) ≤ e ( λ − x(cid:15) √ τ (cid:90) L ( e − ξ(cid:15) √ τ + e ξ(cid:15) √ τ ) dξ ≤ (cid:15) √ τ e ( λ − x + L (cid:15) √ τ ≤ (cid:15) √ τ e λL (cid:15) √ τ for x ∈ [ L , + ∞ ). Therefore, we get the desired inequality (cid:90) + ∞ (cid:12)(cid:12)(cid:12) e − x + ξ(cid:15) √ τ + sgn( x − ξ ) e − | x − ξ | (cid:15) √ τ (cid:12)(cid:12)(cid:12) e λx(cid:15) √ τ | u ( ξ ) | dξ ≤ C u (cid:15) √ τ e λL (cid:15) √ τ . (cid:3) Proof to lemma 2.5 (i). (cid:12)(cid:12)(cid:12) φ ( x ) − e − x(cid:15) √ τ (cid:12)(cid:12)(cid:12) = e − L(cid:15) √ τ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) e − x(cid:15) √ τ − e x(cid:15) √ τ e L(cid:15) √ τ − e − L(cid:15) √ τ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = e − L(cid:15) √ τ | φ ( x ) | . (cid:3) Proof to lemma 2.5 (ii).
Since φ ( x ) = e x(cid:15) √ τ − e − x(cid:15) √ τ e L(cid:15) √ τ − e − L(cid:15) √ τ , we see that φ (cid:48) ( x ) = (cid:15) √ τ e x(cid:15) √ τ + e − x(cid:15) √ τ e L(cid:15) √ τ − e − L(cid:15) √ τ > φ ( x ) ≤ φ ( L ) = 1 for x ∈ [0 , L ]. (cid:3) Proof to lemma 2.5 (iii). φ (cid:48) ( x ) = (cid:15) √ τ e x(cid:15) √ τ + e − x(cid:15) √ τ e L(cid:15) √ τ − e − L(cid:15) √ τ gives that φ (cid:48)(cid:48) ( x ) = (cid:15) τ φ ( x ) >
0, and hence φ (cid:48) ( x ) ≤ φ (cid:48) ( L ) = (cid:15) √ τ e L(cid:15) √ τ + e − L(cid:15) √ τ e L(cid:15) √ τ − e − L(cid:15) √ τ = (cid:15) √ τ e L(cid:15) √ τ +1 e L(cid:15) √ τ − ≤ (cid:15) √ τ if (cid:15) (cid:28) x ∈ [0 , L ]. (cid:3) Acknowledgments
CYK would like to thank Prof. L.A. Peletier for introducing MBL equation and
Mathematical Biosciences Institute at OSU for the hospitality and support.
References [1] J. L. Bona, H.-Q. Chen, S. M. Sun, and B.-Y. Zhang,
Comparison of quarter-plane and two- point boundary value problems: the BBM-equation , Discrete Contin. Dyn. Syst. (2005), no. 4, 921–940. MR MR2166711 (2006m:35314) [2] J. L. Bona and L.-H. Luo, Initial-boundary value problems for model equations for the prop- agation of long waves , Evolution equations (Baton Rouge, LA, 1992), Lecture Notes in
Pure and Appl. Math., vol. 168, Dekker, New York, 1995, pp. 63, 65–94. MR MR1300420 (95i:35137) [3] S.E. Buckley and M.C. Leverett,
Mechanism of fluid displacement in sands , Petroleum Trans- actions, AIME (1942), 107–116. [4] B. Cockburn, C. Johnson, C.-W. Shu, and E. Tadmor,
Advanced numerical approximation of nonlinear hyperbolic equations , Lecture Notes in Mathematics, vol. 1697, Springer-Verlag,
Berlin, 1998, Papers from the C.I.M.E. Summer School held in Cetraro, June 23–28, 1997,
Edited by Alfio Quarteroni, Fondazione C.I.M.E.. [C.I.M.E. Foundation]. MR MR1729305 (2000h:65004) [5] B. Cockburn, G. E. Karniadakis, and C-W (Eds.) Shu,
Discontinuous galerkin methods: The- ory, computation and applications , Lecture Notes in Computational Science and Engineering, [6] A.T. Corey,
The interrelation between gas and oil relative permeabilities , Producer’s Monthly (1954), no. 1, 38–41. [7] D. A. DiCarlo, Experimental measurements of saturation overshoot on infiltration , Water
Resources Research (2004), 4215.1 – 4215.9. [8] S.M Hassanizadeh and W.G. Gray, Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries , Adv. Water Resour. (1990), 169–186. [9] , Thermodynamic basis of capillary pressure in porous media , Water Resour. Res. (1993), 3389–3405. [10] H. Hugoniot, Propagation des Mouvements dans les Corps et specialement dans les Gaz
Parfaits (in French) , Journal de l’Ecole Polytechnique (1887), 3–97. [11] G-S Jiang and E. Tadmor, Nonoscillatory central schemes for multidimensional hyper- bolic conservation laws , SIAM J. Sci. Comput. (1998), no. 6, 1892–1917 (electronic). MR 1638064 (99f:65128) [12] A. Kurganov and D. Levy,
A third-order semidiscrete central scheme for conservation laws and convection-diffusion equations , SIAM J. Sci. Comput. (2000), no. 4, 1461–1488 (elec- tronic). MR MR1797891 (2001j:65127) [13] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes , Commun. Comput. Phys. (2007), no. 1, 141–163. MR MR2305919 (2007k:35320) [14] R. J. LeVeque, Numerical methods for conservation laws , second ed., Lectures in Mathematics
ETH Z¨urich, Birkh¨auser Verlag, Basel, 1992. MR MR1153252 (92m:65106) [15] ,
Finite volume methods for hyperbolic problems , Cambridge Texts in Applied Math- ematics, Cambridge University Press, Cambridge, 2002. MR MR1925043 (2003h:65001) [16] D. Levy, G. Puppo, and G. Russo,
Compact central WENO schemes for multidimen- sional conservation laws , SIAM J. Sci. Comput. (2000), no. 2, 656–672. MR 1780619 (2001d:65110) [17] W. J. Macquorn Rankine, On the thermodynamic theory of waves of finite longitudinal dis- turbance , Royal Society of London Philosophical Transactions Series I (1870), 277–288.
OUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION35 [18] H. Nessyahu and E. Tadmor,
Nonoscillatory central differencing for hyperbolic conservation laws , J. Comput. Phys. (1990), no. 2, 408–463. MR MR1047564 (91i:65157) [19] O. A. Ole˘ınik, Discontinuous solutions of non-linear differential equations , Uspehi Mat. Nauk (N.S.) (1957), no. 3(75), 3–73. MR MR0094541 (20 [20] C. J. Van Duijn, A. Mikelic, and I.S. Pop, Effective Buckley-Leverett equations by homoge- nization , Progress in industrial mathematics at ECMI (2000), 42–52. [21] C. J. van Duijn, L. A. Peletier, and I. S. Pop,
A new class of entropy solutions of the
Buckley-Leverett equation , SIAM J. Math. Anal. (2007), no. 2, 507–536 (electronic). MR MR2338418 (2008g:35136) [22] Y. Wang,
Central schemes for the modified buckley-leverett equation , Ph.D. thesis, The Ohio
State University, 2010. [23] Y. Xu and C-W. Shu,
A local discontinuous Galerkin method for the Camassa-Holm equation , SIAM J. Numer. Anal. (2008), no. 4, 1998–2021. MR MR2399405 (2009e:65140) [24] , Local discontinuous Galerkin method for the Hunter-Saxton equation and its zero- viscosity and zero-dispersion limits , SIAM J. Sci. Comput. (2008/09), no. 2, 1249–1268. MR MR2466156 (2009k:65187)
Department of Mathematics, The Ohio State University, 231 West 18th Ave, Colum- bus, OH 43210
Current address : School of Mathematics, University of Minnesota, 127 Vincent Hall 206 Church
St SE, Minneapolis, MN 55455
E-mail address : [email protected] Department of Mathematics, The Ohio State University, 231 West 18th Ave, Colum- bus, OH 43210; Mathematical Biosciences Institute, The Ohio State University, 1735
Neil Ave, Columbus, OH 43210
E-mail address : [email protected]@math.ohio-state.edu