Fully Conservative Characteristic Methods for Flow and Transport: - - PowerPoint PPT Presentation

Fully Conservative Characteristic Methods for Flow and Transport: Part III, Nonlinear Two-Phase Flow

Todd Arbogast Department of Mathematics

and

Center for Subsurface Modeling, Institute for Computational Engineering and Sciences (ICES) The University of Texas at Austin Chieh-Sen (Jason) Huang, National Sun Yat-sen University (Taiwan) Thomas F. Russell, U.S. National Science Foundation

This work was supported by

• U.S. National Science Foundation
• KAUST through the Academic Excellence Alliance

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 1 The University of Texas at Austin, USA

Outline

• 1. Scalar Conservation Laws and Two-Phase Flow
• 2. Background on Scalar Conservation Laws
• 3. Some Standard Numerical Methods
• 4. A Fully Conservative Streamline Method
• 5. Some Numerical Results
• 6. Summary and Conclusions

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 2 The University of Texas at Austin, USA

Scalar Conservation Laws and Two-Phase Flow

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 3 The University of Texas at Austin, USA

The Scalar Conservation Law Recall the general conservation equation φct + ∇ · (cu) = qc

u(x, t) is the Darcy velocity of the fluid;

φ(x, t) is the porosity of the medium (pore volume/bulk volume); c(x, t) is the concentration of tracer (mass/pore volume); qc(x, t) is an external source or sink of fluid. This equation is a first order hyperbolic partial differential equation for c. It is called a scalar balance law. In the absence of sources and sinks (i.e., qc = 0), we have the scalar conservation law for the transport of tracer particles.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 4 The University of Texas at Austin, USA

Immiscible, Incompressible, Two-phase Flow We consider the general immiscible, incompressible, water-oil system. For each pure phase α, where α = w or o, we have c → 1 (by immiscibility) φ → φsα (effective pore space) The general conservation equation gives the water-oil system φst + ∇ · uw = qw and − φst + ∇ · uo = qo, s = sw is the saturation of water (wetting fluid); 1 − s = so is the saturation of oil (nonwetting fluid);

uα is the Darcy phase velocity (α = w or o);

qα represent wells.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 5 The University of Texas at Austin, USA

Darcy’s Law and Capillary Pressure Darcy’s law governs the phase velocity and pressure (ignoring gravity)

uα = −λαk∇pα,

α = w, o, pα is the phase pressure; k(x) is the absolute permeability; λα(s) is the phase mobility, i.e., relative permeability divided by phase viscosity. The capillary pressure relation pc(s) = po − pw completes the description of the system. This is a system of two nearly hyperbolic flow equations.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 6 The University of Texas at Austin, USA

Rearrangement into Flow and Transport Equations By a rearrangement, we can separate into two systems:

• Elliptic flow or pressure equation;
• (Nearly) hyperbolic transport or saturation equation.

Pressure Equation: The sum of the two equations is ∇ · u = q ≡ qw + qo,

u = −k[λw∇pw + λo∇po] = −k[λ∇pw + λo∇pc] = −k[λ∇po − λw∇pc] u = uw + uo is the total velocity;

λ = λw + λo is the total mobility; Saturation Equation: The water conservation equation φst + ∇ · f(s) = qw where the wetting flux function is

f(s, x) = uw = λw

λ u + kλwλo λ ∇pc

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 7 The University of Texas at Austin, USA

Alternate Saturation Equation However, we could equivalently choose the oil conservation equation φ (1 − s)t + ∇ · F(1 − s) = qo with

F(1 − s, x) = uo = λo

λ u − kλwλo λ ∇pc Key point: Our numerics should not be biased toward either fluid!

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 8 The University of Texas at Austin, USA

A Simplified Model System—1 Assume:

• one space dimension, domain (0, 1);
• appropriate boundary conditions so q = qw = qo = 0;
• φ ≡ 1 for simplicity;
• flux functions depend only on s and not on x (or t);
• z = 1 − s = so.

Remark: The pressure equation is ux = 0, so u is a known constant.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 9 The University of Texas at Austin, USA

A Simplified Model System—2 Model system: Writing both transport equations, we have st + f(s)x = 0 (scalar conservation law for water) zt + F(z)x = 0 (scalar conservation law for oil) s + z = 1 f(s) + F(z) = u (for constant u) Since u is known, this system is redundant. Boundary and initial conditions: For simplicity, assume the phases travel to the right, i.e., u, f, f′, F, and F ′ are nonnegative. s(0, t) = sB(t) and z(0, t) = zB(t) = 1 − sB(t) s(x, 0) = s0(x) and z(x, 0) = z0(x) = 1 − s0(x) where 0 ≤ sB(t) ≤ 1 and 0 ≤ s0(x) ≤ 1.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 10 The University of Texas at Austin, USA

Background on Scalar Conservation Laws

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 11 The University of Texas at Austin, USA

The Linear Case General scalar conservation equation in 1-D: ct +

• f(c)
• x = 0

Linearity: Asume f is linear, meaning there is some v independent of c, such that f(c) = cv Then we have ct + (cv)x = 0 In general, v = v(x, t).

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 12 The University of Texas at Austin, USA

The Case of Constant Velocity If v is constant, with an IC, we have ct + vcx = 0, −∞ < x < ∞, t > 0, c(x, 0) = c0(x), which is solved by c(x, t) = c0(x − vt). If v > 0, we have a wave traveling to the right, of fixed shape.

✲ ✻

x c c0 c(x, t)

✲ t ✲ t ✲ t

Particles simply translate to the right with velocity v. We could have a jump in the IC, which propagates as a contact discontinuity.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 13 The University of Texas at Austin, USA

Inflow Boundary Conditions On a domain 0 < x < ℓ, we have only one BC. If v > 0 is constant, we need to specify c on the inflow side x = 0, but not on the outflow side x = ℓ. That is, ct + vcx = 0, 0 < x < ℓ, t > 0, c(0, t) = cI(t), c(x, 0) = c0(x), which is solved by c(x, t) =

  

c0(x − vt), x − vt ≥ 0, cI(t − x/v), x − vt < 0. ℓ

✲ ✻

x c c0 c(x, t)

✲ ① ✲ ① ✲ ①

× ×

✲ ①

“cI” Particle enters domain at x = 0 and some t > 0 and then transports in time. Particles leave the domain.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 14 The University of Texas at Austin, USA

Streamlines Fluid particles travel along paths called streamlines. Given a starting position x, the path would be (ξ(t), t), where Given a starting position x, the streamline is (ξ(t), t): dξ dt = v

• ξ(t), t
• ,

t > 0, ξ(0) = x. Given v(x, t), this is easily solved.

✲

t

✻

ξ(t) x Space Time In the linear case and when vx = 0, c is constant along streamlines, since d dtc

• ξ(t), t
• = ct + cx ξt = ct + cx v = ct + (cv)x = 0.

If we know where the particle starts, we can simply follow it in time along the streamlines.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 15 The University of Texas at Austin, USA

Streamlines for Constant Velocity If v is constant, then

    

dξ dt = v, t > 0, ξ(0) = x. = ⇒ ξ = x + vt, which are straight lines.

✲ ✻ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡

c0

✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡

cI t x x − vt

t t ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✢

Space Time

• Conclusion. If you know c0 and cI, you know the solution for all time.

We saw that in fact c(x, t) =

  

c0(x − vt), x − vt ≥ 0, cI(t − x/v), x − vt < 0.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 16 The University of Texas at Austin, USA

The Nonlinear Burgers’ Equation In practice, often v depends on c. For example, we have the Burgers’ equation ct + c cx = 0 ⇐ ⇒ ct + (f(c))x = 0, where f(c) = 1

2 c2.

Note that the velocity v(c) = c/2 increases as c increases. Remark: In general f = f(c, x, t), but we will assume f = f(c) only.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 17 The University of Texas at Austin, USA

Streamlines Fluid particles travel along the streamlines. Given a starting position x, the streamline is (ξ(t), t): dξ dt = v(c(ξ(t), t)) = f

• c(ξ(t), t)
• c(ξ(t), t)

, t > 0, ξ(0) = x. This is not solvable without c.

✲

t

✻

ξ(t) x Space Time Problem: The concentration is not constant along streamlines! In general, d dtc

• ξ(t), t
• = ct + cx ξt = ct + f(c)

c cx = ct +

• f(c)
• x = 0.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 18 The University of Texas at Austin, USA

Characteristics Paths of constant concentration are called characteristics. Given a starting position x, the characteristic is (ξ(t), t): dξ dt = f′ c(ξ(t), t)

• ,

t > 0, ξ(0) = x. This is not solvable without c.

✲

t

✻

ξ(t) x Space Time The concentration is constant along characteristics: d dtc

• ξ(t), t
• = ct + cx ξt = ct + f′(c) cx = ct +
• f(c)
• x = 0.

Therefore, characteristics are straight lines! Given c we can find them! Remark: For linear problems, f(c) = cv, so f(c)/c = f′(c) = v and streamlines and characteristics coincide. Conclusion: If we know c at some point x, we can simply follow it in time along the characteristic. Question: Is this correct? Could anything go wrong?

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 19 The University of Texas at Austin, USA

The Riemann Problem We consider the Riemann Problem ct + (f(c))x = 0 c(x, 0) =

  

cL, x < 0, cR, x > 0, for which the IC has just two values or states. Burgers’ Equation. f(c) = 1 2c2.

• Streamlines: f(c)

c = 1 2c

• Characteristics: f′(c) = c

Thus particles move with half the speed of the characteristics, which are ξ(t) =

  

x + cLt, x < 0, x + cRt, x > 0.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 20 The University of Texas at Austin, USA

Rarefactions Consider cL = 0 and cR = 1. Then the characteristics are ξ(t) =

  

x, x < 0, x + t, x > 0.

✲ ✻ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡

Space Time Question: What happens between the two sets of characteristics? Does a vaccuum open up? The solution will form a rarefaction wave, i.e., a spreading of the solution to smooth out the mass distribution. It has the form c(x, t) =

      

cL, x < f′(cL) t, γ(x/t), f′(cL) t < x < f′(cR) t, cR, f′(cR) t < x,

✲ ✻

x c where f′(γ(s)) = s. For Burgers’ equation, f′(s) = s, so γ(x/t) = x/t.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 21 The University of Texas at Austin, USA

Shocks and the Rankin-Hugoniot Speed Consider cL = 1 and cR = 0. Then the characteristics are ξ(t) =

  

x + t, x < 0, x, x > 0.

✲ ✻ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡

x t

✲ ✻ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✡ ✂ ✂ ✂ ✂ ✂ ✂ ✂ ✂ ✂ ✂ ✂ ✂

x t Question: Does mass build up where the two sets of characteristics collide? The solution develops a shock wave, i.e., a discontinuity to maintain mass conservation. The shock travels at the Rankin-Hugoniot speed σ σ = f(cL) − f(cR) cL − cR and c(x, t) =

  

cL, x < σ t, cR, σ t < x.

✲ ✻

x c For Burgers’ equation, σ = 1/2.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 22 The University of Texas at Austin, USA

Summary of Behavior Conservation laws (without diffusion) exhibit:

• Particles transport along the streamlines
• Information or waves transport along the characteristics
• Concentration is constant along the characteristics
• Characteristics are straight lines
• Smooth spreading or rarefaction waves where characteristics diverge
• Discontinuities that
• simply transport with the flow (a contact discontinuity)
• or form as shocks where characteristics collide
• Remark. If diffusion is present, discontinuities cannot form, but the

solution can be nearly discontinuous, i.e., they have steep or sharp fronts.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 23 The University of Texas at Austin, USA

Changes to the Concentration Variable—1 Burgers’ equation is ct +

1

2 c2 x = 0,

Multiply through by 2c and set ξ = c2, to obtain that 2cct + 2 c

1

2 c2 x = 0

= ⇒

• c2

t +

2

3 c3 x = 0

= ⇒ ξt +

2

3 ξ3/2 x = 0.

This is a conservation law for ξ. Question: Are these two equations essentially the same?

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 24 The University of Texas at Austin, USA

Changes to the Concentration Variable—2 Suppose a shock forms. With what speed does it travel? σc =

1 2c2 L − 1 2c2 R

cL − cR and σξ =

2 3ξ3 L − 2 3ξ3 R

ξL − ξR = 2

3(ξ2 L + ξL ξR + ξ2 R)

= 1

2(cL + cR)

= 2

3(c4 L + c2 L c2 R + c4 R).

We are conserving the wrong quantity, so it travels at the wrong speed!

• Theorem. There are infinetly many solutions to a hyperbolic

conservation law. The physically relevant solution is the entropy solution, which is c = lim

ǫ→0 cǫ

where cǫ,t + (fcǫ)x−ǫcxx

= 0.

diffusion That is, the entropy solution is the one that is stable with respect to random diffusive processes.

• Remark. We always have some physical diffusion. However, it may be

very small and unresolved by a numerical method. Thus, the numerical method must respect the complexity of the nearly hyperbolic equation.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 25 The University of Texas at Austin, USA

Some Standard Numerical Methods

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 26 The University of Texas at Austin, USA

Conservative Methods We generally use a cell centered grid

✲ x

xi−2 xi−1 xi+1

① ① ①

xi

①

xi−1/2 xi+1/2 × × xi−5/2 xi−3/2 xi+3/2 × × × A locally conservative method for ct + (f(c))x = 0 is a method such that cn+1

i

− cn

i

∆t + Fi+1/2 − Fi−1/2 h = 0, where the numerical flux Fi±1/2 ≈ f(c) at xi±1/2 needs to be defined. A simple choice. The following implicit method is unstable! cn+1

i

− cn

i

∆t + f

1

2(cn+1 i+1 + cn+1 i

)

• − f

1

2(cn+1 i

+ cn+1

i−1 )

• h

= 0. Question: If v = f′(c) > 0, which direction is the fluid going? Across the point xi+1/2, what is the usual value of c?

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 27 The University of Texas at Austin, USA

Upstream Weighted Finite Differences A simple stable and implicit finite difference approximation to ct + (f(c))x = ct + f′(c)cx = 0 is cn+1

i

− cn

i

∆t + f

• cn+1

up,i+1/2

• − f
• cn+1

up,i−1/2

• h

= 0, wherein the concentration is upstream weighted, meaning that cup,i+1/2 =

  

ci if vi+1/2 = f′ > 0, ci+1 if vi+1/2 = f′ < 0. xi

①

xi−1/2 xi+1/2 × × This is well defined unless f′ is particularly nasty.

• Theorem. The method is O(∆t + h) accurate.

This is not a very accurate method, and we need to use small ∆t and h.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 28 The University of Texas at Austin, USA

Numerical Diffusion The upstream weighted finite difference method suffers from excessive numerical diffusion. That is, it is stable because it smears sharp features in the solution, so it does not achieve high resolution.

• Example. A shock will be approximated like

✲ x ✻

c True shock Approximation

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 29 The University of Texas at Austin, USA

Godunov’s Method Approximate c as a piecewise constant

✻

c

✲ x

xi−2 xi−1 xi xi+1

① ① ① ① ① ① ① ①

xi−5/2 xi−3/2 xi−1/2 xi+1/2 xi+3/2 × × × × × Godunov’s method is explicit. One advances the solution by solving the Riemann problem for left and right states at each grid line xi+1/2. CFL stability constraint. We do not allow the characteristics from two grid lines to interact, so we must limit the time step size. This limitation is the Courant-Friedrichs-Lewy (CFL) constraint: ∆t ≤ h max |f′(c)|.

• Theorem. If the CFL constraint is satisfied, then Godunov’s method is

well defined, stable, and has error O(∆t + h). Moreover, it has minimal numerical diffusion among the locally conservative explicit methods.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 30 The University of Texas at Austin, USA

A Fully Conservative Streamline Method

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 31 The University of Texas at Austin, USA

The Nonlinear Problem Model system: st + f(s)x = 0 (scalar conservation law for water) zt + F(z)x = 0 (scalar conservation law for oil) s + z = 1, f(s) + F(z) = u. Since u is known, this system is redundant. Difficulties:

• Streamlines are not characteristics.
• The streamlines of the two fluids differ.
• Streamlines follow nonlinear trajectories.
• Shocks and rarefactions may arise.
• A single streamline cannot be determined alone

(it depends on neighboring particles).

✲

tn tn+1

✻

x

ˇ xw ˇ xo

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 32 The University of Texas at Austin, USA

Streamlines The particle trace-back of the point x is ˇ x = ˇ x(x; t). Water streamlines: For water particles,

    

dˇ xw dt = f(s(ˇ xw, t), t) s , tn ≤ t < tn+1 Not f′(s)! ˇ xw(tn+1) = x Oil streamlines: For oil particles,

    

dˇ xo dt = F(z(ˇ xo, t), t) z , tn ≤ t < tn+1 Not F ′(z)! ˇ xo(tn+1) = x

✲

tn tn+1

✻

x

ˇ xw ˇ xo Remarks:

• We use streamlines, not characteristics, since we are interested in the

movement of particles, i.e., mass.

• Particle trajectories are difficult to solve accurately even in 1-D due to

the fact that s is not constant on the streamline.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 33 The University of Texas at Austin, USA

Saturation Changes on a Streamline On the streamline, the saturation is ˇ sw(t) = s(ξ(t), t) This the change is dˇ sw dt = st(ξ, t) + sx(ξ, t) ξ′ = −

• f(s(ξ, t))
• x + sx(ξ, t) ξ′

= −f′(ˇ sw)ˇ sx + ˇ sx f(ˇ sw) ˇ sw Saturation change equations: dˇ sw dt = ˇ sx

f(ˇ

sw) ˇ sw − f′(ˇ sw)

• and

dˇ zo dt = ˇ zx

F(ˇ

zo) ˇ zo − F ′(ˇ zo)

• Remark: We still do not know ˇ

sx = −ˇ zx.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 34 The University of Texas at Austin, USA

Numerical Streamlines We use a Taylor approximation. Given a guess sg ≈ sn+1(x) and δ ≈ sx, we have the numerical water streamlines ˇ xw ≈ x − f(sg) sg ∆t − 1 2sg

f(sg)

sg − f′(sg)

• 2

δ ∆t2 ˇ sn

g ≈ sg −

f(sg)

sg − f′(sg)

• δ ∆t

Similarly, the numerical oil streamlines come from zg = 1 − sg ≈ zn+1(x) and δ ≈ −zx = sx as ˇ xo ≈ x − F(zg) zg ∆t + 1 2zg

F(zg)

zg − F ′(zg)

• 2

δ ∆t2 ˇ zn

g ≈ zg +

F(zg)

zg − F ′(zg)

• δ ∆t

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 35 The University of Texas at Austin, USA

Local Mass Constraints We trace element E back along water and oil streamlines to ˇ Ew and ˇ Eo.

✲ x ✻

tn tn+1 E ˇ Ew ˇ Eo Integrate the conservation equations in space-time over these regions, and use the divergence theorem. Local water and oil mass constraints:

• E sn+1 dx =
• ˇ

Ew

sn dx and

• E zn+1 dx =
• ˇ

Eo

zn dx Local volume constraint: (Sum the two mass constraints) |E| =

• E dx =
• ˇ

Ew

sn dx +

• ˇ

Eo

zn dx

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 36 The University of Texas at Austin, USA

Saturation Approximations We approximate the saturation s in three different ways:

• 1. As a piecewise discontinuous constant sn

E and zn E = 1 − sn E in each

grid element E.

• These will satisfy the local mass constraint.
• 2. An approximation to the saturation sg = si at the grid point xi.
• These will not satisfy the local mass constraint.
• Used for streamline tracing.
• 3. As a piecewise discontinuous linear ˜

s(x), arising from postprocessing the sE and si.

• These will satisfy the local mass constraint.
• Formally higher order
• Our final and most accurate approximation of s(x).

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 37 The University of Texas at Austin, USA

Saturation Postprocessing On E = (xi−1, xi), define the slopes σi−1 = sE − si−1

1 2(xi − xi−1)

and σi = si − sE

1 2(xi − xi−1)

, and slope limited σE =

      

if σi−1 and σi have opposite signs, σi−1 if not the case above and |σi−1| ≤ |σi|, σi

• therwise.

Then s(x, tn) ≈ ˜ sn(x) = sn

E + σE

• x − 1

2(xi + xi−1)

• for x ∈ E.

xi−1 xi E σE = 0

✈

si−1

✈si

sE σi−1 σi xi−1 xi E

✈

si−1

✈si

sE σi−1= σE σi xi−1 xi E σE = 0

✈

si−1

✈si

sE σi−1 σi= σE

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 38 The University of Texas at Austin, USA

Volume Constraint and Consistency We work from left to right. For each phase tn−1 tn

④

xi−2

④

xi−1

④

xi

④

xi+1 E

④

ˇ xi−2

④

ˇ xi−1

④

ˇ xi ˇ E

✲ ✛

sg, δ ˇ si

?

= ˜ sn−1(ˇ xi) We need to adjust sg and δ so that the volume constraint is met. The system is underdetermined, so we add two consistency conditions. Minimization problem: F(sg, δ) =

• ˜

sn(ˇ xw) − ˇ sn

g

2 +

• ˜

zn(ˇ xo) − ˇ zn

g

2

subject to the volume constraint |E| =

• ˇ

Ew

sn dx +

• ˇ

Eo

zn dx, 0 ≤ sg ≤ 1, and possibly constraints on δ.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 39 The University of Texas at Austin, USA

Time-Step Advance For each grid element, sn+1

E

= 1 |E|

• ˇ

Ew

˜ sn dx Postprocess these to form ˜ sn+1(x). Phase independence: The scheme is unbiased to either phase by

• definition. Because of volume connservation, it is equivalent to advance

zn+1

E

= 1 |E|

• ˇ

Eo

˜ zn dx

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 40 The University of Texas at Austin, USA

Some Numerical Results

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 41 The University of Texas at Austin, USA

1D Buckley-Leverett Problem Governing equations st + ∇ · (f(s)) = 0, x ∈ (0, 1), t ∈ (0, 1), zt + ∇ · (F(z)) = 0, x ∈ (0, 1), t ∈ (0, 1), where s + z = 1, F(z) = 1 − f(s), and the two fractional flows are the same f(σ) = F(σ) = σ2 σ2 + (1 − σ)2.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 42 The University of Texas at Austin, USA

A Simple Buckley-Leverett Test—1

V1 V2

0.2 0.4 0.6 0.7 0.8 0.9 1

V1 V2

0.2 0.4 0.6 0.7 0.8 0.9 1

V1 V2

0.2 0.4 0.6 0.7 0.8 0.9 1

V1 V2

0.2 0.4 0.6 0.7 0.8 0.9 1

V1 V2

0.2 0.4 0.6 0.7 0.8 0.9 1

time = 0.5, at steps 2 dt = 2 h = 0.25, n = 8 Blue -- Cell value Green - Grid value Squares -- uniform bounds on delta Diamonds -- delta_x >0

Notes:

• 1. Shock position and front saturation are assumed.
• 2. δx ∼ sxx =

−f′′′ t2(f′′) > 0 and δ(0) ∼ sx = 1 tf′′(1).

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 43 The University of Texas at Austin, USA

A Simple Buckley-Leverett Test—2 δ versus x

V1 V2

0.2 0.4 0.6

• 1
• 0.8
• 0.6
• 0.4
• 0.2

V1 V2

0.2 0.4 0.6

• 1
• 0.8
• 0.6
• 0.4
• 0.2

Squares(red)-- uniform bounds on delta Diamonds(blue) - delta_x>0 Circles(black) -- almost s_x

V1 V2

0.2 0.4 0.6

• 1
• 0.8
• 0.6
• 0.4
• 0.2

Notes:

• 1. δx ∼ sxx =

−f′′′ t2(f′′) > 0 and δ(0) ∼ sx = 1 tf′′(1).

• 2. Results are not too far off even if uniform bounds are given to δ.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 44 The University of Texas at Austin, USA

Shock formation—1 h = 1/20 and ∆t = h = 1.2∆tCFL

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

Triangles and squares show s at (left) t = 0 and t = 0.1, just before the shock forms, (right) t = 0.2 and t = 0.3, after shock formation. Green symbols stand for grid values and blue for average element values.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 45 The University of Texas at Austin, USA

Shock formation—2 h = 1/20 and ∆t = 2h = 2.4∆tCFL

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

The blue squares show the average element s at t = 0.3. For comparison, the red curves are the average element values from ∆t = h = 1/20.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 46 The University of Texas at Austin, USA

Burgers’ Equation st + ∇ · (f(s)) = 0, x ∈ (0, 1), t ∈ (0, 1), zt + ∇ · (F(z)) = 0, x ∈ (0, 1), t ∈ (0, 1), where s + z = 1, F(z) = 1/2 − f(s) and f(s) = s2 2 with the initial and boundary conditions s(0, t) = 1 and s(x, t0) = 0.5 1 + e2k(x−0.5) + 0.5. Notes:

• 1. f(s) + F(z) = 1/2 in this case; otherwise, the oil particle velocity will

go to infinity, which is non-physical.

• 2. s(x, t0) is a smooth approximation of the Heaviside step function

with left state sL = 1 and right state sR = 0.5 at x = 0.5.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 47 The University of Texas at Austin, USA

A Simple Burgers’ Equation Test

0.4 0.6 0.8 0.5 0.75 1 0.4 0.6 0.8 0.5 0.75 1 0.4 0.6 0.8 0.5 0.75 1 0.4 0.6 0.8 0.5 0.75 1

k = 40 for Heavside function Square,

• 100 <= delta <= 100

Blue

• cell values

Green - grid values

0.4 0.6 0.8 0.5 0.75 1 0.4 0.6 0.8 0.5 0.75 1

V1 V2

0.2 0.4 0.6 0.8 1

• 4
• 3
• 2
• 1

V1 V2

0.2 0.4 0.6 0.8 1

• 4
• 3
• 2
• 1

Notes:

• 1. We remove the assumptions on knowing the shock position and

front saturation.

• 2. The shock moves at the right speed, sL + sR

2 = 3 4 (the shock moves 0.075 since ∆t = 0.1).

• 3. Uniform bounds are assigned to δ, and it behaves as sx.

Riemann problems h = 1/20 and ∆t = 2h = 2∆tCFL

0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1

0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0.2 0.4 0.6 0.8 1 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0.2 0.4 0.6 0.8 1 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0.2 0.4 0.6 0.8 1 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0.2 0.4 0.6 0.8 1 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1

Rarefaction (left), shock (center), and both (right). Green squares are grid s values, blue squares are average element values, and red line is true solution.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 49 The University of Texas at Austin, USA

Summary and Conclusions

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 50 The University of Texas at Austin, USA

Summary and Conclusions

• 1. The problem. Streamline tracing not possible.
• 2. The solution. Trace so as to satisfy (two) local mass constraints

and saturation consistency.

• 3. The results.
• Time steps greater than ∆tCFL can be used;
• Accurately approximates shocks and rarefactions;
• Little numerical diffusion on coarse grids;
• Mass locally conserved for both water and oil.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 51 The University of Texas at Austin, USA