An Iterative Perturbation Method for the Pressure Equation in the Simulation of Miscible Displacement in Porous Media.
The miscible displacement problem in porous media is modeled by a nonlinear coupled system of two partial differential equations: the pressure-velocity equation and the concentration equation. An iterative perturbation procedure is proposed and analyzed for the pressure-velocity equation, which is capable of producing as accurate a velocity approximation as the mixed finite element method, and which requires the solution of symmetric positive definite linear systems.  Only the velocity variable is involved in the linear systems, and the pressure variable is obtained by substitution. Trivially applying perturbation methods can only give an error $O(\epsilon)$, while our iterative scheme can improve the error to $O(\epsilon^m)$ at the $m$th iteration level, where $\epsilon$ is a small positive number.  Thus the convergence rate of our iterative procedure is $O(\epsilon)$, and consequently a small number of iterations is required.  Theoretical convergence analysis and numerical experiments are presented to show the efficiency and accuracy of our method.
Introduction
Miscible displacement occurs, for example, in the tertiary oil-recovery process which can enhance
hydrocarbon recovery in the petroleum reservoir. This process involves the injection of a solvent
at injection wells with the intention of displacing resident oil to production wells. The resident
oil may have been left behind after primary production by reservoir pressure and secondary
production by waterflooding. Since the tertiary process requires expensive chemicals and the
performance of the displacement is not guaranteed, its numerical simulation plays an important
role in determining whether enough additional oil is recovered to make the expense worthwhile
and in optimizing the recovery process of hydrocarbon.
Department of Mathematics/Institute of Applied Mathematics, University of British Columbia, Vancouver,
British Columbia, Canada V6T 1Z2. E-mail address: lin@math.ubc.ca
y Department of Mathematics, Wayne State University, Detroit, Michigan 48202, USA. E-mail address:
yang@math.wayne.edu or dyang@na-net.ornl.gov
Mathematically, miscible displacement in porous media is modeled by a nonlinear coupled
system of two partial differential equations with appropriate boundary and initial conditions.
The pressure-velocity equation is elliptic, while the concentration equation is parabolic but
normally convection-dominated. The concentration equation is derived from the conservation
of mass which involves the Darcy velocity of the fluid mixture, but the pressure variable does
not appear in the concentration equation. Thus a good approximation of the concentration
equation requires accurate solution for the velocity variable. Mixed finite element methods
have been applied [9, 10, 11, 15, 16, 25] to the pressure-velocity equation, which can yield the
velocity one order more accurate than the finite difference or element methods. However, the
finite-dimensional spaces for velocity and pressure need to satisfy the Babuska-Brezzi condition,
and the resulting linear system has a nonpositive definite coefficient matrix. Besides, the number
of degrees of freedom in the linear system doubles that of finite difference or element methods.
Thus efforts have been made to improve the performance of mixed finite element methods
[5, 18, 12].
The purpose of this paper is to propose and analyze an iterative perturbation method for the
pressure-velocity equation, which is capable of producing as accurate velocity approximation as
the mixed finite element method, and requires the solution of symmetric positive definite linear
systems. Only the velocity variable is involved in the linear systems and the pressure variable is
obtained by substitution without solving any linear systems at each iteration level. Aside from
this, the finite element spaces for our velocity and pressure variables need not satisfy the so-called
Babuska-Brezzi condition. The iterative perturbation method is a variant of the augmented
Lagrangian method [18, 5, 6] applied directly to the continuous differential problem, which can
also be viewed as the sequential regularization method [2, 20, 3] applied to stationary problems.
Unlike the augmented Lagrangian method using spectral analysis to discuss the convergence and
its rate for the discretized problems, we use the method of asymptotic expansion directly for
the differential problem following the idea in the sequential regularization method [2, 20]. The
asymptotic method is easier to use for more general and more complicated problems than spectral
analysis since the former applies readily to non-symmetric and infinite-dimensional operators.
Trivially applying perturbation methods (e.g. penalty methods) can only give an error O(ffl),
where ffl is a small positive perturbation parameter. We will prove that our iterative scheme can
improve the error to O(ffl m ) at the m-th iteration level. In other words, the convergence rate
of our iterative procedure is O(ffl). Theoretical convergence analysis and numerical experiments
show that the number of iterations is extremely small, usually two or three.
The organization of the rest of the paper is as follows. In Section 2, we describe our iterative
perturbation method for the time-discretized problem at the partial differential equation level.
In Section 3 we show its convergence and the rate of convergence. Then in Section 4 we give a
fully-discretized version for the iterative perturbation method for the pressure-velocity equation
and a Galerkin method for the concentration equation. Finally, in Section 5 we present numerical
examples to demonstrate the effectiveness and accuracy of our method, and in Section 6 we give
some remarks and future directions.
2 The Iterative Perturbation Method
Consider the miscible displacement of one incompressible fluid by another in a porous reservoir
\Omega ae R 2 over a time period denote the pressure and Darcy velocity
of the fluid mixture, and c the concentration of the invading fluid. Then the mathematical
model is a coupled nonlinear system of partial differential equations
2\Omega \Theta J;
(1)
div
2\Omega \Theta J;
(2)
OE
@c
@t
2\Omega \Theta J;
with the boundary conditions
and initial condition
c) is the mobility tensor of the fluid mixture, fl(x; c) and d(x) are the gravity
and vertical coordinate, q is the imposed external rates of flow, OE(x) is the porosity of the rock,
D(x; t) is the coefficient of molecular diffusion and mechanical dispersion of one fluid into the
other, c) is a known linear function of c representing sources, and - is the exterior
normal to the boundary
We assume that the mobility is symmetric and positive definite. For existence of p, we
assume that the mean value of q is zero and for uniqueness we impose that p have mean value
zero.
In recent years much attention has been devoted to the numerical simulation of this problem.
In this paper we are interested in solving the velocity-pressure equation (1)-(2) using an iterative
perturbation method (IPM) for the time-discretized problem. The iterative perturbation method
considered herein is a variant of augmented Lagrangian methods [18, 5] and sequential
regularization methods [2, 20]. We will analyze the IPM using the technique presented in [2, 20],
since the analysis gives the convergence of the iterative procedure and its convergence rate at
the same time, and is applicable to non-symmetric problems.
After a time discretization, we obtain the following system for u and p at current time step:
2\Omega \Theta J;
div
2\Omega \Theta J;
where ~ c is an approximation of c assumed to be known here. Taking ffl to be a small positive
number, we replace the system (1), (2), (4) by the following iterative perturbation method
such that
a(~c)
2\Omega \Theta J;
and
2\Omega \Theta J;
where the initial guess p 0 is required to satisfy the zero mean value property:
R
p s has mean value zero from (11) and (10). We note that, if taking p 0 j 0, each iteration is a
kind of penalty method (probably without its meaning in the optimization context).
Our iterative perturbation method (IPM) has the following salient features:
1. We solve a small system (9)-(10) for velocity u, and get pressure p from (11) directly. We
will show that the accuracy of such a method is O(ffl s ) at the s-th iteration level. Note
that the system (9)-(10) is well-posed since, unlike the usual penalty method, we need not
take ffl very small.
2. The velocity-pressure equation was recently solved by the mixed finite element method
[9, 10, 11, 15, 16, 25], in which the resulting linear system has a non-positive definite
coefficient matrix, and the discrete spaces for u and p need to satisfy the Babuska-Brezzi
condition. While in the IPM procedure, the system (9)-(10) leads to a symmetric positive
definite coefficient matrix. Since u and p are obtained from equations (9)-(10) and (11)
separately, compatibility conditions between the discrete spaces for u and p are not needed.
3. When the standard finite element method [7, 8, 13, 14, 17, 23, 24] is applied for the pressure
equation, the velocity need to be obtained by finite differencing the pressure variable,
which gives less accuracy. The velocity in our method is obtained directly without finite
differencing. Note that the accuracy of the approximate velocity is important, since the
concentration equation involves the velocity only.
4. We will see that in xx4 and 5 the discrete version of our IPM scheme (9)-(11) gives the
same accuracy for the velocity as the mixed method and requires the solution of well-conditioned
linear systems like Galerkin methods. Note that our numerical experiments
will show that a few (usually two or three) iterations are enough for the IPM scheme.
We now describe some notation to be used throughout the rest of the paper. As usual, we
use kuk p and kukH m to denote the standard norms in the Sobolev spaces L
respectively, where 1 the subscript p in the
norm notation when 2. In addition, we define the divergence space
with the following two norms:
We shall denote by (\Delta; \Delta) and h\Delta; \Deltai the inner product and duality
in\Omega and on \Gamma, respectively.
For a normed linear space B with norm k \Delta kB and a sufficiently regular function
we define
ff
2 and kgk L 1
If [ff; we simplify the notation as kgk respectively. We shall also
denote by M and K generic constants, which may be different at different occurrences. We use
M for continuous problems and K for discrete problems.
3 Convergence Analysis
Before stating our convergence theorem, we first give two lemmas.
There exists a unique solution fu; pg, where u 2 H(div) and p 2 H
1(\Omega\Gamma , to the
problem
2\Omega \Theta J;
div
2\Omega \Theta J;
where p and q have mean value zero. Furthermore, there exits a constant M such that the
following estimate holds:
Proof: Substituting (12) into (13) and (14) we see that p satisfies a Poisson equation with
Neumann boundary condition
a @p
@-
Noting the zero mean value of p and using the standard results of the Poisson equation, we
obtain the uniqueness and existence of p and the estimate
An application of the trace inequality[19]
leads to the inequality (15) for p. Then the existence, uniqueness and the estimate (15) for u
directly.
Writing (9)-(11) as
2\Omega \Theta J;
div
2\Omega \Theta J;
we can obtain a similar estimate for u s and p s from Lemma 1:
Lemma 2 For the solution of the following system:
a(~c)
2\Omega \Theta J;
we have the stability estimate:
Proof: Multiplying both sides of the first equation by u and applying Green's formula
leads to (20).
We are now ready to describe our convergence theorem and prove it using the method of
asymptotic expansion.
Theorem 1 Let fu; pg be the solution of system (7), (8), and (4) and fu s the solution of
(9)-(11). Then there exists a constant M , independent of s and ffl, such that
Proof: We first consider the case
in (9). Comparing the coefficients of like powers of ffl, we thus have
5div
5div
where the boundary condition (10). The equation (22) has infinitely
many solutions in general. We should choose u 10 not only to satisfy (22), but also to ensure
that the solution of (23) exists. A choice of u 10 is the exact solution u of (7), (8), and (4), i.e.
where p is the exact solution of (7), (8), and (4), and satisfying zero mean value condition.
Hence, from Lemma 1, we have has the following form
5 div u
The equation (28) suggests that we choose u 11 and a corresponding p 11 (with zero mean value)
to satisfy
div
According to Lemma 1, u 11 and p 11 exist and have the bound:
Generally, assuming that u 1(i\Gamma1) and p 1(i\Gamma1) have been found for i - 2, we choose u 1i and p 1i
(with zero mean value) satisfying
div
Applying Lemma 1 and (32), we obtain that all u 1i and p 1i exist uniquely and
Next we estimate the remainder of the asymptotic expansion to the m\Gammath power of ffl. Denote
the partial sum and the reminder by
-
Then, from (25)-(27), (29)-(31), and (33)-(35), w 1m satisfies, for m - 1,
Then, using (32), (36) and Lemma 2, we obtain
Noting that using (36), we thus have proved the estimate (21) for u 1 . On the other
hand, (11) can be rewritten as
div
Using (26), (30), (34) and (39), we have
. By taking m - 2, this proves (21) for
Now we look at the second iteration
Note that (40) gives us a series expansion for p 1 . Plugging the expansions of u 2 and p 1 into
and comparing like powers of ffl we obtain
div
5 div u
5 div u
Again, required to satisfy the boundary condition (10). As in the case of
choose u 20 j u and thus have
5 div u
This suggests that we construct u 21 and p 21 to satisfy
div
where p 21 has mean value zero. Obviously u is the solution of (44)-(46).
In general, similar to the case of (with zero mean value) to satisfy
div
m. By the same procedure as in the case of we obtain the error equations
similar to (37) and (38) with an addition of a remainder term 5(p
1m ) on the right-hand
side. simulating the proof of Lemma 2 and using (40) (i.e. kp
we have
where the remainder w 2m satisfies the same estimate (39) as w 1m . Noting u 20 j u and u
we obtain (21) for u 2 . Then, using (11), (40), (45), (48) and the estimate of w 2m , we conclude
that
. By taking m - 3, this gives (21) for p 2 since p 21 j 0.
We can repeat this procedure, and by induction, conclude the results of the theorem.
Remark 2 It is not difficult to see that are a solution to a problem in the
form of (12)-(14). Then using the estimates (15) we may obtain a better error estimates for p:
From Theorem 1 we see that the convergence rate of our iterative scheme (9)-(11) is O(M ffl)
M . This implies that the number of iterations needed to achieve a prescribed accuracy
is very small. The fast convergence of our method makes it dramatically different from penalty
methods.
4 The Galerkin Approximation
In this section, we approximate the velocity-pressure iterative perturbation scheme (9)-(10) and
the concentration equation (3) by using the standard Galerkin method.
4.1 The Approximation Scheme
\Gammag and Wg. The variational
form of (9)-(10) can be written into the following: find
such that
(a
The weak form of the concentration equation (3) reads: find
1(\Omega\Gamma such that
@t
For h u ? 0 and an integer k - 0, with respect to the velocity-pressure equation, we introduce
finite element spaces W h ae W and Y associated with a quasi-regular
triangulation
of\Omega into triangles or rectangles of diameter less than h u . Similarly, we
denote by Z h ae H 1
(\Omega\Gamma the finite-dimensional space for the concentration equation with the grid
size h c and approximation index l. Assume that the following approximation properties hold:
z h 2Z h
where K is a constant. The space W h can be taken to be the vector part of the Raviart-Thomas
[22] space of index k, or Brezzi-Douglas-Marini [4] space of index k + 1.
Given a partition of J;
g. Let fp be fp; u; cg and its approximation at
time level t n . We define our approximation scheme at time t n by the following.
Step 1: Given C n , find fU follows. Take the initial guess -
iteratively -
h such that
(a
where q h is a smoothed source and sink function. Let U
U s and P
Step 2: When U n known, find C n+1 2 Z h such that
\Deltat n
Note that in step 1 of the scheme, the initial guess could be more efficiently taken as -
Our numerical experiments will show that the number of
iterations can usually be taken to be or 3 for a range of perturbation parameter
to 10 \Gamma5 .
For convection-dominated displacement problems, sharp fluid interfaces move along characteristic
or near-characteristic directions. Thus the modified method of characteristics [8, 11, 15,
16, 23, 24, 26] may be applied to treat the concentration equation. In this case, the accurate
Darcy velocity computed using the IPM helps to determine the characteristic direction more
precisely.
4.2
For our error analysis below we shall make use of the following elliptic projections of fu; cg :
(a
(div (R u
such that
(D(u)r(R
and q is the right-hand side function of equation (2) [9]. This - function
is chosen to assure coercivity of the projection form. Then, it can be shown that [9, 10]
ck
ck
From (56) and the definition of Y and Y h , for Y we can also find fR
We now start our error analysis by combining (52), (58) and (61) to get
(a
Then we see that
In view of (63), (65) and (68) we have
Presuming that kp
u , we thus have (for
On the other hand, following the way of getting (67) and combining (52), (11), (58) and (59),
we get
(a
It is not difficult to see that p s presuming p Hence we can find R p
satisfying (66). Since -
h we can take w 0 2 W h such that div w
according to [1, Lemma 3.3]
(assuming\Omega is a bounded, Lipschitz-continuous, and connected
domain), we have
Hence we take in (71) and obtain
In view of (66), (69) and (72), we have
Thus the error estimates (70) and (74) are reduced to bounding C \Gamma Rc. We denote that
Combining (53), (60), (62) and taking z = e n+1 as a test function we obtain the error
equation
\Deltat n
re n+1 )
\Deltat n
\Deltat n
re
We now estimate the right-hand side terms of the error equation (75). For T 1 we have by
Since
c \Deltat n kc t k 2
we see that
\Deltat
Assume that D(u) is Lipschitz [9], by (63), (64), (65), and (70) we have
c kck 2
where j is a small positive number. Obviously,
We rewrite T 5 in the following form that
re
is bounded, from (63), (64), (65) and (70) we obtain
c kck 2
Note that g is a linear function of c, we have
c kck 2
Finally, we see that the left-hand side of (75) dominates2\Deltat n
positive number.
Substituting (76)-(82) into (75) and choosing j small we obtain the following error inequality
\Deltat
where
c kck 2
Applying Gronwall's lemma to (84) we see that
\Deltat n
We are now ready to state our main results for scheme (58)-(60).
Theorem 2 Let fp; u; cg be the solution to Problem (3)-(2), and fP; U; Cg the solution to the
scheme (58)-(60) with s iterations at each time step. Then there exists a constant K, for \Deltat
sufficiently small,
c kck L
\Deltat n kC
c kck L
This theorem tells us that for sufficiently small perturbation parameter ffl, the error estimates
for the velocity, pressure and concentration are optimal.
U s
R
U s
U s

Figure

1. One element with velocity on each edge and pressure at the center.
5 Numerical Experiments
In this section, we present some numerical examples to show how well our iterative scheme
performs, and how the parameter ffl affects the number of iterations and accuracy. For simplicity,
we will just consider the pressure-velocity equation, since the concentration equation has been
analyzed previously [9, 10, 11, 14, 15, 16, 17, 23, 24, 25].
Consider the elliptic problem with Neumann boundary condition
div
where
\Omega is a square and \Gamma its boundary. More general
domain\Omega will not present technical
problems. For simplicity, we also assume that the coefficient a is a scalar.
The approximation scheme takes the form: Find U s 2 W h , for
(a
Partition the
domain\Omega into a set of squares of side length h. We take the space W h to be
the vector part of the Raviart-Thomas [22] space of index 0. Thus
where P k is the set of one variable polynomials of order less than or equal to k. Consequently,
the approximate pressure P s lies in the space of piecewise constants. Partitioning the domain
into triangles or rectangles, or applying higher order approximation polynomials can be treated
analogously.
Let U s
ff denote the constant value of the flux in positive x or y-direction on the edge ff;
L, R, B, T (representing left, right, bottom, and top, respectively), of each element. See Figure
1. Consider w to be the basis function (on the standard reference square). Applying
the trapezoidal rule to (88) we
where q ff is the value of q at the middle point of edge ff. Similarly, letting
simplifying we have the following linear system for each element.
R U s
Note the equation (89) has the discrete version on each element:
From equations (91)-(94) we can easily form the element stiffness matrix. Then assembling
all the element matrices and taking into account the boundary condition we obtain the stiffness
matrix. The force vector can be gotten in an analogous way.
To test the convergence, the following stopping criterion will be adopted
where k is a positive integer, and k \Delta k1 denotes the discrete L 1 norm. When this criterion is
satisfied, the iterative process is stopped and the solutions at iteration s are adopted.
We now apply the iterative procedure (91)-(95) to some test problems on Sun SPARCstation
IPC with computation in C++ data type double. We will use the stopping criterion (96) with
pressure errors are measured for iterates fU against the exact
solutions under the L 1 norm, although in the stopping criterion (96) errors are measured for
iterates only. The initial guesses are always chosen to be zero, so all errors are 1.00 before the
iterative procedure starts.
Example 1: Let velocity u and pressure p satisfy
div
@ The true solutions for velocity u and pressure p are given
by
sin(-x)
sin(-y)
The pressure p and external flow rate div u are chosen in such a way that they both have
mean value zero.
Example 2: Let velocity u and pressure p satisfy the problem
div
1. The true solutions for
velocity u and pressure p are given by
Example 3: Let velocity u and pressure p satisfy the nonhomogeneous problem
div
y
10. The function g is chosen such that
the true solutions for velocity u and pressure p are given by
For Example 1, the results with (uniform) grid sizes 1and 1are shown in Tables 1 through
4. Although more than 2 iterations are required for the iterative procedure to stop when
or 10 \Gamma2 , the approximate velocity and pressure at iteration 2 are accurate enough. For Example
2, the results with grid sizes 1and 1are shown in Tables 5 and 6.
The results of Example 3 with grid sizes 1and 1are shown in Tables 7 and 8. However,
the pressure is less accurate than velocity in this example. This might be caused by the fact
errors are calculated relative to the exact solutions and shown in the L 1 norm. Initial guesses
are always chosen to be zero, so all relative errors are 1.00 at iteration 0.
iteration velocity pressure velocity pressure velocity pressure

Table

2: Numerical results for Example 1 with grid size = 1. Both velocity errors and pressure
errors are calculated relative to the exact solutions and shown in the L 1 norm. All runs stop
after 2 iterations except for the case \Gamma6 in which 1 iteration is required.
iteration velocity pressure velocity pressure velocity pressure
that the pressure iterates fP s g converge to the true pressure p up to a constant, since p has
mean value zero while fP s g do not.
From

Tables

through 8 we conclude that our iterative method performs as well as the theory
predicts. In particular, it can achieve good accuracy for velocity, while the linear systems solved
are symmetric and positive definite. Also, the computational work of our method is much smaller
than that of mixed methods, since the number of iterations required is usually very small.

Table

3: Numerical results for Example 1 with grid size = 1. Both velocity errors and pressure
errors are calculated relative to the exact solutions and shown in the L 1 norm.
iteration velocity pressure velocity pressure velocity pressure
errors are calculated relative to the exact solutions and shown in the L 1 norm. All runs stop
after 2 iterations except for the case \Gamma6 in which 1 iteration is required.
iteration velocity pressure velocity pressure velocity pressure

Table

5: Numerical results for Example 2 with grid size = 1. Both velocity errors and pressure
errors are calculated relative to the exact solutions and shown in the L 1 norm.
iteration velocity pressure velocity pressure velocity pressure

Table

Numerical results for Example 2 with grid size = 1. Both velocity errors and pressure
errors are calculated relative to the exact solutions and shown in the L 1 norm.
iteration velocity pressure velocity pressure velocity pressure

Table

7: Numerical results for Example 3 with grid size = 1. Both velocity errors and pressure
errors are calculated relative to the exact solutions and shown in the L 1 norm.
iteration velocity pressure velocity pressure velocity pressure
iteration velocity pressure velocity pressure velocity pressure
6 Concluding Remarks
We have proposed an iterative procedure for the pressure-velocity equation in the numerical
simulation of miscible displacement in porous media. This procedure is first analyzed at the
differential level and then discretized by finite element methods. Theoretical analysis and numerical
experiments show that this procedure converges at the rate of O(ffl), where ffl is a small
positive number. The fast convergence rate and the ease of choosing relaxation parameter make
our iterative method different from penalty methods and Uzawa's algorithm. Indeed, our numerical
experiments show that two or three iterations are usually enough for a variety of problems,
corresponding to
Compared with mixed finite element methods, the discrete version of our scheme can provide
the same accurate approximations for velocity and pressure, which is crucial in reservoir problems
since velocity is intimately involved in the concentration equation. However, in contrast to mixed
finite element methods, our scheme requires only the solution of symmetric and positive definite
linear systems which have a smaller number of degrees of freedom corresponding to the velocity
variable. Since our method can completely decouple the velocity and pressure variables, the
so-called Babuska-Brezzi condition is not needed in constructing the finite dimensional spaces
for velocity and pressure.
In view of the advantages of our iterative method, we can conclude that it can lead to great
savings in computer memory and small execution time of the numerical algorithm. Also, it
has the capability of effectively dealing with heterogeneous and anisotropic media in which the
permeability tensor may be non-diagonal, rapidly-varying and even discontinuous. In [27], one
of the authors conducted numerical experiments for the miscible misplacement using the IPM
for the pressure-velocity equation and a modified characteristic method for the concentration
equation.
Finally, we point out that the iterative procedure presented in this paper can be applied
easily to three-dimensional problems.

Acknowledgement

: The authors would like to thank Professors Uri Ascher and Jim Dou-
glas, Jr. for their inspiration and support in the pursuit of this research.



--R

Decomposition of vector spaces and application to the Stokes
Sequential regularization methods for higher index DAEs with constraint singularities: Linear index-2 case
Sequential regularization methods for nonlinear higher index DAEs

Mixed and Hybrid Finite Element Methods
Discontinuous upwinding and mixed finite elements for two-phase flows in reservoir simulation




On the approximation of miscible displacement in porous media by a method of characteristics combined with a mixed method
Inexact and preconditioned Uzawa algorithms for saddle point problems
The Mathematics of Reservoir Simulation
Efficient time-stepping methods for miscible displacement problems in porous media
Simulation of miscible displacement using mixed methods and a modified method of characteristics
Convergence analysis of an approximation of miscible displacement in porous media by mixed finite elements and a modified method of characteristics
Galerkin Methods for miscible displacement problems in porous media
Applications to the Numerical Solution of Boundary-Value Problems
Finite Element Methods for Navier-Stokes Equations
A sequential regularization method for time-dependent incompressible Navier-Stokes equations
Regularization methods for differential equations and their numerical solution
A mixed finite element method for second order elliptic problems
Finite elements with characteristic finite element method for a miscible displacement problem
Time stepping along characteristics with incomplete iteration for a Galerkin approximation of miscible displacement in porous media
Mixed methods with dynamic finite element spaces for miscible displacement in porous media
A characteristic mixed method with dynamic finite element space for convection-dominated diffusion problems
Numerical simulation of miscible displacement in porous media using an iterative perturbation algorithm combined with a modified method of characteristics.
--TR
miscible displacement;iterative method;perturbation method;galerkin method;flow in porous media