The Pressure Equation

By [ 7 ], we can get the method of solving the pressure equation.

In each block B_ij , U_β, β = L, R, B, T is the out normal component of the total flux across the edges of the block. H is the side length of block. The equations to be solved for P^m and U^m are :

The LCELM is the variational type of the MMOC. We shall explain the MMOC briefly.

5.2.1 MMOC procedure

The MMOC procedure for the waterflood problem is based on introducing a character-istic derivative for the transport part of the saturation equation written in non-divergence form.

( Note that the characteristic direction ϕ depends on x , the saturation, and the fluid velocity.) Thus, equation( 15 ) can be written as

Θ∂S

∂ϕ + divV = (1 − Λ_w)q⁺ and equation( 17 ) can be written as

Θ∂ζ_n,k

∂ϕ = (1 − Λ_w(ζ_n,k))q⁺, x ∈ Ω, t_n,k < t ≤ t_n,k+1 (31) We compute the transport microstep by solving equation( 31 ) with the initial values given by equation( 19 ). As seen in Figure 1, the fundamental concept in the MMOC is the discretization of the characteristic derivative by backwards differencing along the tangent to the characteristic through the point (x, t_n,k+1) back to the time level t_n,k for whatever x−point arise in the quadratures used in the finite element scheme. If we ignore for the moment sources and sink and the boundary, the transported values over a micro-step would be defined by

x_n,k(x) = x − Λ⁰_w(ζ_n,k)(E₁U )(x, t_n,k+1)∆t_st

Φ (32)

ζ_n,k+1(x) = ζ_n,k(x_n,k(x)) (33)

The only reasonable criterion [ 8 ] for conservation of mass globally is that the map ( 32 ) have Jacobian identically one.

Figure 1: The point (x, t_n,k+1) back to the time level t_n,k.

5.2.2 LCELM procedure

Recall that the saturation equation can be written in divergence form as

∇_t,x·

Then the transport equation can be written as

∇_t,x·

followed by the diffusive part given by

Φ∂S

∂t + divV = 0 (36)

Let Γ = Ω × [t_n,k, t_n,k+1]. Let G be a reasonable shaped, simply-connected subset of Ω , and define a subset C = C_n,k(G) of Γ as follows. For each x ∈ ∂G , construct the solution

y(x; t) of the final value problem determined by G ,G and the integer curves ( 37 )-( 38 ); see Figure 2 and Figure 3 for an example C in a single space variable setting. ( We let ∆t_st be sufficiently small. Then the map x → x_n,k is one-to-one, so that this construction can be carried out.) Now, denote the outward normal to ∂C by σ(x, t) and note that it is orthogonal to the vector



Thus, mass is conserved locally in the transport step, as defined in equation( 35 ) or equations( 17 )-( 19 ) above, if

The no-flow boundary condition is handled in a natural way in equation( 40 ), since the integral curves ( 32 )-( 34 ) do not exit Ω in this case. In fact, if x ∈ ∂Ω , then

the integral curve remain in ∂Ω and C has a portion of its lateral surface contained in Γ = Ω × [t_n,k, t_n,k+1] . Hence, no special cares arise for subsets G close to the boundary for these boundary conditions.

Figure 2: The space-time domain C.

Figure 3: The area (G, t_n,k+1) back to the time level t_n,k.

5.3 Diffusive fractional step

We shall apply backward differencing in time over [ t_n, t_n+1 ]

V_n+1= −KD(S_n)(∇S_n+1+ (ρ_w− ρ_o)g∇z P_c⁰(S_n) ), ΦS_n+1− S_n

∆t_s + divVn+1 = 0, with the no-flow boundary condition

V_n+1· −→n = 0.

Thus, the mixed finite element equations take the form

( 1

KD(S_n)V_n+1 , −→n₁) − (S_n+1 , div−→n₁)

= −( 1

P_c⁰(S_n)(ρ_w− ρ_o)g∇z , −→n₁) , −→n₁ ∈ η (ΦS_n+1− S_n

∆t_s , n₂) + (divV_n+1 , n₂) = 0 , n₂ ∈ W.

6 Numerical results

The following data and functions are held fixed for the computational results:

Viscosity µ_w = 0.5 cP µ_o = 10 cP

Density ρ_w = 1 g/cm³ ρ_o = 0.7 g/cm³

Porosity Φ = 0.2

Residual saturations S_rw = 0.2 S_ro = 0.15 Absolute permeability K = 6 mdarcy

Relative permeability K_rw(S) = ^(S−S_(1−S^rw)2

In this study, the geometric domain is 256 × 256 meters and the grid size of the model is 1 × 1 meters. In order to solve the transport equation, we use two kinds of partitions as shown in Figure 4 and Figure 5 to compute the injection and other sections, respectively. The red line in Figure 4 and Figure 5 represents the computational domain.

In the following figures, water is injected at the lower right corner at a uniform rate and a mixture of water and oil produced at the top left corner.

Figure 4: The kind of partition (black line) is used to calculate the whole domain except the injection.

Figure 5: This kind of partition (green line) is used to calculate the injection.

Figure 6: The water saturation at the 40th day.

Figure 7: The water saturation at the 80th day.

Figure 8: The water saturation at the 120th day.

Figure 9: The water saturation at the 160th day.

Figure 10: The water saturation at the 200th day.

Figure 11: The water saturation at the 240th day.

22

References

[1] M. Murad and J. Cushman A multiscale theory of swelling porous media, II:Dial porosity models for consolidation of clay incorporating Phys- iochemical effects, Preprint #287, Center for Applied Mathematics, Purdue University, August 1996.

[2] S. N. Antontsev, On the solvability of boundary value problems for degenerate two-phase porous flow equations, Dinamika Splosnoi Sredy Vyp., 10 (1972) 28-53. In Russian.

[3] G. Chavent, A new formulation of diphasic incompressible flows in porous media, in “ Applications of Methods of Functional Analysis to Problems in Mechanics ”,Lecture Notes Mathematics. 503 (1976) 258- 270,Springer-Verlag, Berlin, New York, (P. Germain and B. Nayroles, eds.)

[4] G. Chavent and J.Jaffre, “Mathematical Models and Finite Elemenes for Reservoir Simulation” , North-Holland, Amsterdam, 1986.

[5] J. Douglas, Jr., Superconvergence in the pressure in the simulation of miscible displacement, SIAM J. Numer. Anal., 22 (1985) 962-969.

[6] J. Douglas, Jr., R. E. Ewing, and M. F. Wheeler, The approximation of the pressure by a mixed method in the simulation of miscible displacement, R.A.I.R.O., Anal. Numer., 17 (1983) 17-33.

[7] J. Douglas, F. Pereira, and L.M. Yeh. A parallelizable method for two-phase flows in naturally fractured reservoirs. Computational Geosciences, 1(3):333–368, 1997.

[8] J. Douglas Jr, C.S. Huang, and F. Pereira. The modified method of characteristics with adjusted advection. To appear in Numerische

23

Mathematik; available as Technical Report #298 , Center for Applied Mathematics, Purdue University, June 1997.

[9] J. Douglas, F. Pereira, and L.M. Yeh. A locally conservative Eulerian- Lagrangian numerical method and its application to nonlinear transport In porous media.