RIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

NUMERICAL APPROXIMATION
AND IDENTIFICATION
IN A 1 - D PARABOLIC
DEGENERATED NON-LINEAR
DIFFUSION AND TRANSPORT
EQUATION

Guy CHAVENT
Gary COHEN

Rapport de Recherche N° 270

Décembre 1977
# NUMERICAL APPROXIMATION AND IDENTIFICATION IN A 1-D PARABOLIC DEGENERATED NON-LINEAR DIFFUSION AND TRANSPORT EQUATION

G. Chavent*, G. Cohen**

## Résumé :

On s'intéresse à la récupération, par injection d'eau, de l'huile contenue dans un barreau poreux. Losque le débit d'injection d'eau est élevé, le terme de transport est prépondérant dans l'équation : nous donnons un schéma de discrétisation qui reste valable dans le cas limite de l'équation du transport. L'identification des non-linéarités apparaissant dans l'équation d'état est alors effectuée par minimisation d'un critère quadratique d'erreur sur l'observation. Des exemples numériques sont donnés.

## Abstract :

The physical motivation of this study is the recovering of oil from a porous medium by injection of water. For high water injection rate, the transport term is preponderant, thus we give an approximation by discontinuous finite elements, which also remains valid for the pure transport equation. Identification of non-linearities appearing in the state equation is then performed by minimizing an output error criterion. Numerical examples are given.

To appear in the Proceedings of the 8th.IFIP Conference on Optimization Techniques, Springer 1978.

* IRIA/LABORIA et Université de Paris IX, Dauphine 75016 Paris.

** IRIA/LABORIA, Domaine de Voluceau B.P. 105 F - 78150 Le Chesnay
NUMERICAL APPROXIMATION AND IDENTIFICATION
IN A 1-D PARABOLIC DEGENERATED NON-LINEAR
DIFFUSION AND TRANSPORT EQUATION

by

G. Chavent (*) and G. Cohen (**)

ABSTRACT

The physical motivation of this study is the recovering of oil from a porous medium
by injection of water(1). For high water injection rate, the transport term is prepon-
derant, thus we give an approximation by discontinuous finite elements, which also
remains valid for the pure transport equation. Identification of non-linearities
appearing in the state equation is then performed by minimizing an output error
criterion. Numerical examples are given.

1. THE PHYSICAL PROBLEM AND ITS MODELLIZATION.

Let us consider a core sample of homogeneous porous material, which we suppose to
be *initially filled with oil*. We inject (cf. fig. 1) some water through its left
hand side end with a *known rate*, and we *observe the rates of oil and water recover-
ed* through the right hand side end (2):

<figure><img src="image_1.png" /><figcaption>Fig. 1</figcaption></figure>

The water does not appear instanteanously at the right hand side end, but only after
the break-through time (B.T.), and then in growing quantities.

The state of the system is the (reduced) water saturation $u(x,t)$ at point $x$ and
time $t$, ranging between 0 and 1. The saturation is not easily accessible to meas-
urement.

<footnote>(1) This study has been made in the frame of a research contract n° 05/77 between IRIA and the IFP (Institut Français du Pétrole) which we thank for allowing us to publish those numerical results.</footnote>

<footnote>(2) This experience corresponds to the laboratory study of core samples obtained during the drilling of wells in the oil industry ; but the two-phase flow process involved here has other applications, as the study of unsaturated soils for instance.</footnote>

<footnote>(*) IRIA, Domaine de Voluceau, POB. 105, F-78150 Le Chesnay and University of PARIS-IX, 75016 Paris, France.</footnote>

<footnote>(**) IRIA, same address.</footnote>
The state equation is (of [1] [2] [3]) :

$$
(1) \quad \oint \frac{\partial u}{\partial t} - \frac{\partial}{\partial x} \left( K \frac{\partial u}{\partial x} \alpha(u) \right) + q(t) \frac{\partial}{\partial x} (b(u)) = 0 \quad \text{in } ]0,L[ \times ]0,T[
$$

with the boundary conditions :

$$
(2) \quad \left\{
\begin{array}{l}
\text{Dirichlet at } x = 0: \\
u(0,.) \triangleq u_e \\
u_e \equiv 1
\end{array}
\right.
\qquad (3) \quad
\left\{
\begin{array}{l}
\text{Unilateral condition at } x = L: \\
u(L,.) \triangleq u_s \\
u_s \le 1, \varphi_1 \ge 0, (1-u_s)\varphi_1 = 0
\end{array}
\right.
$$

$$
\left\{
\begin{array}{l}
\text{Unilateral condition at } x = L: \\
u(L,.) \triangleq u_s \\
u_s \le 1, \varphi_1 \ge 0, (1-u_s)\varphi_1 = 0
\end{array}
\right.
$$

and an initial condition :

$$
(4) \quad u(x,0) = u_0(x)
$$

where the water output rate $\varphi_1$ is related to the saturation $u$ by :

$$
(5) \quad \varphi_1(t) = (1+b(u))q - K \frac{\partial}{\partial x} \alpha(u) .
$$

The parameters appearing in (1) to (5) are the porosity $\emptyset$ and the permeability $K$,
that are, in our case, easily determined by preliminary experiences, and the two
functions of saturation $\alpha(u)$ and $b(u)$, related to the relative permeabilities and
capillarity pressure curves that are not precisely known. Their typical forms are gi-
ven in Fig. 2. One may notice that the derivative $\alpha'$ of $\alpha$ vanishes for $u=0$
and $u=1$, so that the diffusion term of equation (1) is degenerated. The mathemati-
cal study of system (1) to (5) under the above conditions, has been made in [4],
where the existence of a solution has been proved.

<figure><img src="image_2.png" /><figcaption>Fig. 2</figcaption></figure>

The identification problem then consists in trying to determine $\alpha$ and/or $b$ from
the measurement $\varphi_{1\text{ob}}$ of the water output rate $\varphi_1(t)$.

This can be achieved in minimizing the observation error functional :

$$
(6) \quad J(\alpha, b) = \int_{0}^{T} (\varphi_{1}(t) - \varphi_{1\text{ob}}(t))^{2} dt = \text{Min}
$$

over a suitable set of functions $\alpha$ and/or $b$. A proof of the existence of a solution
to that optimization problem and a calculus of the Gâteaux-derivative of $J$ may be
found in [3], for a non-degenerate diffusion term and condition (3) being replaced
by a Dirichlet condition.

The choice of a good approximation scheme for the state equations (1) to (5) is a
necessary preliminary to the numerical resolution of the optimization problem (6).
We present in § 2 such an approximation, and in § 3 its use for the optimization
problem.

## 2. NUMERICAL APPROXIMATION OF THE STATE EQUATIONS

Divide $[0,L]$ into $N$ intervals of length $h = \frac{L}{N}$, and set $x_i = ih$ $i = 0,1...N$.
The two ideas which guided us in the choice of the approximation are the following :

i) for high water injection rate $q(t)$, equation (1) behaves quite as a first order hyperbolic equation, and the solution $u$ can present very stiff fronts. So $u$ was approximated by discontinuous finite elements, following the ideas introduced by Lesaint [5] for the linear first order hyperbolic equation :

$$ (7) \qquad \left\{ \begin{aligned} u &\approx u_h \in V_h \\ V_h &= \{ v \in L^2(0,L) \mid |v| \}_{x_i, x_{i+1}}^{x_i, x_{i+1}} \quad \epsilon \in \rho^{k, i=0,1...N-1} \end{aligned} \right. $$

where $\rho^k$ is the set of polynomials in $x$ of degree $\le k$.

ii) we need a good numerical approximation for the water output rate $\phi_1(t)$ (as it appears in the functional (6) to be minimized !), which includes a derivative of $u$. As we have already chosen a discontinuous $u_h$, the only way out is to use a mixed approximation, following the ideas of [6] : besides $u$ we must approximate the quantity $r = -K \frac{\partial}{\partial x}(\alpha(u)) = -K\alpha'(u) \frac{\partial u}{\partial x}$ (and not $\frac{\partial u}{\partial x}$ itself, that is not square integrable as $\alpha'$ vanishes for $u=0$ and $u=1$) :

$$ (8) \qquad \left\{ \begin{aligned} r &= -K \frac{\partial}{\partial x}(\alpha(u)) \approx r_h \in Q_h \\ Q_h &= \{ v \in C([0,L]) \mid |v| \}_{x_i, x_{i+1}}^{x_i, x_{i+1}} \quad \epsilon \in \rho^{k+1}, \quad i=0,1...N-1 \end{aligned} \right. $$

We shall use the straightforward notations :

$$ (8\text{bis}) \qquad b_h = b(u_h), \quad b_{hi}^+ = b(u_h^+(x_i)), \quad b_{hi}^- = b(u_h^-(x_i)), \quad \text{etc...} $$

Following the above guideline, we are led to the following semidiscretization in space of the state equations (1) to (5) :

$$ (9) \quad \left\{ \int_{x_i}^{x_{i+1}} \partial_t \frac{du_h}{dt} v + \int_{x_i}^{x_{i+1}} \frac{\partial r_h}{\partial x} v + q(t) \left\{ (b_{hi}^+ - b_{hi}^-) v(x_i) + \int_{x_i}^{x_{i+1}} \frac{\partial b_h}{\partial x} v \right\} = 0, \quad \forall v \in \rho^k, \forall i = 0,1...N-1, \forall t \in \mathcal{T}[0,T] \right. $$
with

$$
(10) \qquad \bar{u}_{ho} = u_e \qquad \forall t \in [0,T]
$$

and

$$
(11) \quad \int_0^L \frac{r_h s_h}{K} \d y = \int_0^L \alpha_h \frac{\partial s_h}{\partial y} \d y + \alpha(u_e) s_h(0) - \alpha(u_s) s_h(L) \quad \forall s_h \in Q_h, \forall t \in [0,T]
$$

$$
(12) \quad u_e \equiv 1, \quad u_s \le 1, \quad \varphi_{1h} \ge 0, \quad (1-u_s)\varphi_{1h} = 0
$$

$$
(13) \qquad u_h(x,0) = u_{ho}(x)
$$

$$
(14) \qquad \varphi_{1h}(t) = (1+b_{hN}^{-}(t))q(t) + r_{hN}(t) .
$$

$$
\frac{\partial u_k}{\partial t} \quad \text{by} \quad \frac{u_h^{n+1} - u_k^n}{\tau}
$$

The time discretization is simply made by replacing $\frac{\partial u_k}{\partial t}$ by $\frac{u_h^{n+1} - u_k^n}{\tau}$ ($\tau$ being the time step) and using an explicit scheme in (9).

This approximation main features are the following :

. it is valid from the pure parabolic case (q ≡ 0 : imbibition process) to
the pure hyperbolic case (α ≡ 0 : Buckley-Leverett approximation)

. it satisfies exactly the mass balance for each of the two fluids separatel y ; for water for instance, one checks easily that

$$
(15) \left\{
\begin{array}{l}
\displaystyle \frac{d}{dt} \int_{0}^{L} \emptyset u_{h}^{n}(t) dx = \psi_{1h}^{n}(t) - \varphi_{1h}^{n}(t) \quad \text{where} \\
\\
\displaystyle \psi_{1h}^{n}(t) = \left[1 + b_{ho}^{-}(t)\right] q(t) + r_{ho}^{n}(t) = \begin{array}{l} \text{water input rate at } x=0. \\ \text{and } x=L. \end{array}
\end{array}
\right.
$$

. it gives a precise determination of the water break-through time . it is
the first time when u_s becomes equal to 1 .

. for $k=0$, $u_h$ is a piecewise constant function, and the equations (9) to (11) turn out to be identical, for an adequate choice of quadrature formulas, to the usual finite difference scheme. The equations (12)(14) give however a new formulation of the B.C. at $x=L$. We shall refer to that case when using the term "finite differences".

. for $k=1$, $u_h$ is a discontinuous piecewise linear function, and we get new schemes depending on the choice of the quadrature formulas we use. The best choice turned out to be to use the Simpson formula in all integrals appearing in (9) and in (11). We shall refer to that case when using the term "discontinuous finite elements".

Numerical results of the resolution of the state equations :

The data : L = 24cm , Ø = .156 , K = 1.3 10⁻⁸ , σ = 18cm² , q = 4.629 10⁻⁴ cm/s, the functions α and b as in Fig. 2. These data correspond to the rather high water injection rate of 1 cm³/mn, and a stiff water saturation front is expected. The
actual saturations (ranging from .15 to .63) are used in the numerical computations instead of the reduced saturations (ranging from 0 to 1).

The water saturation front : the saturation distribution after 1000 s. are plotted in Fig. 3 for k = 0 and 1 and for a space step size growing from 1/4 to 1 cm.

<figure><img src="image_3.png" /><figcaption>Fig. 3</figcaption></figure>

More instructive is the superposition of the curves for the steps $h = 1/4$ and $h = 1$.

<figure><img src="image_4.png" /><figcaption>Fig. 4: Comparison of the water front obtained for different schemes and step size</figcaption></figure>

One sees that the position of the front is, for discontinuous finite elements, practically independent from the space step size, which is the case, as it is well known, for the finite difference.

The overall process : we have plotted on Fig. 5 the successive positions of the front every 200 s., for finite elements with a step of 1/4 cm.
<figure><img src="image_5.png" /><figcaption>Fig. 5</figcaption></figure>

We have reinforced the two curves at time $t = 1884$ s. and $t = 1800$ s., corresponding respectively to the water break-through time and to the last curve plotted before that time. One sees that the "surface" saturation $v_s$ has grown up from .15 to .63 during the time interval between the two curves.

The water output rate $\phi_1(t)$ and the break through time (B.T.) : we have plotted on Fig. 6, instead of the water output rate $\phi_1(t)$, the output water/oil ratio (WOR) more commonly used by oil engineers :

<figure><img src="image_6.png" /><figcaption>Fig. 6</figcaption></figure>

Here too, the superiority of finite elements over finite difference is apparent : the break-through time is practically independent from the discretization for the finite elements, and the WOR-curves form remains very similar, which is not the case with the finite differences.
The computation time : for the simulation of the flow up to 2400 s, the computation time required on an IRIS 80 CII is :

computer time in
minute and hundred
of minute

<table><thead><tr><th>(h,τ) in (cm,s)</th><th>(1, 4)</th><th>(1/2,2)</th><th>(1/4,1)</th></tr></thead><tbody><tr><td>Finite differences</td><td>.38</td><td>.79</td><td>2.29</td></tr><tr><td>Finite elements</td><td>0.6</td><td>1.71</td><td>5.53</td></tr></tbody></table>

For the same number of unknowns in space, for instance finite differences with h = 1/4 and finite elements with h = 1/2, the maximum admissible time steps τ are about 5s in both cases, so that the corresponding computing time is 0.46 s for finite differences and 0.68 s for finite elements, i.e. an increase of about 50 % which is not prohibitive on account of the gain in accuracy observed.

The limiting cases : we have drawn on Fig. 7 the solution u obtained with the pure parabolic case (left) corresponding to the imbibition process, and the pure hyperbolic case (right). In the latter case, the physical minimum entropy solution is obtained.

<figure><img src="image_7.png" /><figcaption>Fig. 7</figcaption></figure>

CONCLUSION :

The use of discontinuous finite elements, together with the unilateral condition at the right-hand side of the slab, gives a good approximation of the saturation front and of the water output rate, even for large step sizes : it is this scheme that we shall use in the identification procedure.
### 3. NUMERICAL RESOLUTION OF THE IDENTIFICATION PROBLEM.

Because of the unilateral condition in (12), the state equations (9) to (14) are not differentiable. Thus, we approximate them first, by penalizing (cf. [4]), the corresponding variational inequality, this leading to replace the boundary condition (12) by:

$$ (12\text{bis}) \qquad u_e \equiv 1, \qquad \varphi_{1h} - \frac{1}{2\epsilon} [ (u_s - 1)^+ ]^2 = 0 $$

where $\epsilon > 0$ is given "small enough".

We define the -only semi-discretized in space for sake of (relative) simplicity - approximated criterion $J_h$ by

$$ (16) \qquad J_h(\alpha, b) = \int_0^T (\varphi_{1h}(t) - \varphi_{1ob}(t))^2 dt $$

We define then a Lagrangian with the criterion (16) and the state equations (9)(10)(11)(12bis)(13) which leads us, after some meticulous and tedious calculus, to define a discrete adjoint state $p_h(t) \in V_h, \theta_h(t) \in Q_h, t \in [0,T]$; ( think of $\theta_h$ as being an approximation of $-K \frac{\partial p_h}{\partial x}$). With the system of notation (δbis), the adjoint equations are (compare with (9)(10)(11)(12bis) and (13)):

$$ (17) \qquad \left\{ \begin{array}{l} -\int_{x_{i-1}}^x \int_0^T \varphi_h \frac{\partial p_h}{\partial t} v + \int_{x_{i-1}}^x \alpha'_h \frac{\partial \theta_h}{\partial x} v - q(t) \left\{ b'_{hi} (p_{hi}^+ - p_{hi}^-) v(x_i) + \int_{x_{i-1}}^x b'_h \frac{\partial p_h}{\partial x} v \right\} = 0 \\ \forall v \in \rho^k, \quad \forall i=1...N, \quad \forall t \in [0,T] \end{array} \right. $$

$$ (18) \qquad p_{hN}^{+} = p_s \quad \forall t \in [0,T] $$

$$ (19) \qquad \int_0^L \frac{\theta_h s_h}{K} = \int_0^L p_h \frac{\partial s_h}{\partial x} + p_e s_h(0) - p_s s_h(L) \qquad \forall s_h \in Q_h, \forall t \in [0,T] $$

$$ (20) \qquad p_e \equiv 0, \quad \alpha'(u_s)\theta_{hN} - \frac{1}{\epsilon}(u_s - 1)^+[p_s + 2(\varphi_{1h} - \varphi_{1ob})] = 0 \qquad \forall t \in [0,T] $$

$$ (21) \qquad p_h(x,T) = 0, \quad \forall x \in [0,L]. $$

The first order variation $\delta J_k$ of criterion $J_h$ associated to variations $\delta \alpha$ and $\delta b$ of the parameters is then given by:

$$ (22) \qquad \left\{ \begin{array}{l} \delta J_h = -\int_0^T q(t) \left\{ \sum_{i=0}^N \delta b_{hi}^{-} (p_{hi}^{+} - p_{hi}^{-}) + \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} \delta b_h \frac{\partial p_h}{\partial x} \right\} \\ - \int_0^T \left\{ \delta \alpha(u_e) \theta_{ho} - \delta \alpha(u_s) \theta_{hN} + \int_0^L \delta \alpha_h \frac{\partial \theta_h}{\partial x} \right\} \end{array} \right. . $$
The functions $a$ and $b$ are represented as continuous piecewise linear functions of the saturation $u$:

<caption>Fig. 8</caption>

The discretization of the parameter functions $\alpha$ and $b$

<figure><img src="image_8.png" /></figure>

It is then easy to derive from (22) the exact gradient of $J_h$ with respect to the numerical unknowns $\alpha_j$ and $b_j$ of the optimization problem.

The monotonicity ($\alpha$ and $b$ are known to be the increasing functions of $u$) and regularity constraints ($|\alpha''|$ and $|b''|$ bounded, in order to avoid oscillating non linearities) are taken into account simply by penalization of the objective function (16).

A conjugate gradient method is used.

### Numerical results for the identification problem.

**The data** : The same data as in § 2 are used for the state equation, excepted that a shorter slab ($6$ cm) is used, in order to save computation time, as no information is gained before the water break-through time. Discontinuous finite elements ($k=1$) are used, with a space step of $h=1/2$ cm.

The actual saturations (ranging from .15 to .63) are used, and the interval [.15, .63] is divided, as shown in Fig. 8, into 13 intervals for the discretization of $\alpha$ and $b$, so that there are 14 numerical unknown values $\alpha_j$ and/or 14 values $b_j$. The functions $\alpha$ and $b$ of Fig. 2 have been used for the simulation of the observation $\varphi_{10b}(t)$. We shall refer to these as to the "true" parameters. $\alpha^*, b^*$.

**The "no model-error" case** : The observation $\varphi_{10b}(t)$ is computed, using the "true" parameters $\alpha^*$ and $b^*$, with exactly the same numerical model as used for the identification ((9)(10)(11)(12bis)(13)) with $k=1$ and $h = 1/2$) hence $J(\alpha^*, b^*)=0$. The results obtained when we try to identify $b$ from that observation, are shown in Fig. 9, where the water saturation repartition, $u$ after 300s is also shown. Though the minimization of $J_h$ has been very effective (computed and observed water output rates are identical up to the drawing precision), the computed $b$ turns out to be very different from the "true" one, namely for the saturations less than .4, which seems to have absolutely no influence on the observation. In correspondence with that, one observes a discrepancy between the computed and the "true" water front.
<figure><img src="image_9.png" /><figcaption>Fig. 9</figcaption></figure>

The "model-error" case : Assuming the physical laws underlying equations (1) to (5) to be exact, and our approximation scheme to be converging when h→0, a (noiseless) "actual" observation φ₁₀b(t) can be simulated in using the discretized equations with a small h, chosen here equal to 1/4 cm, so that J(α*,b*)≠0. The results for the identification of b are shown on Fig. 10. A small discrepancy between computed and observed water output rate is still visible, but the model error has not affected the computed b, which is very similar to the one shown in Fig. 8.

<figure><img src="image_10.png" /><figcaption>Fig. 10</figcaption></figure>
The identification of $\alpha$ is shown on Fig.11. The initial value of $J$ is very small which means that for the chosen water injection rate $q$, the diffusion parameter $\alpha$ has very little influence on the water output rate $\varphi_{1h}$, i.e. the state equations are nearly of hyperbolic type. The initial discrepancy between $\varphi_{1h}$ and $\varphi_{1ob}$ also results from the model-error, which, in our case, could be as much important as the error $\alpha - \alpha^*$ on the parameter. Nevertheless, the optimization algorithm works, the output error is still reduced, but of course, the computed $\alpha$ does not look at all like $\alpha^*$.

<figure><img src="image_11.png" /><figcaption>Fig. 11</figcaption></figure>

About the simultaneous identifications of $\alpha$ and $b$: The above numerical results, showing the poor influence of $\alpha$ and $b$ on $\varphi_{1h}$, do not encourage the simultaneous identification of those two functions. Engineers, however, did not know condition (12) and used the Dirichlet B.C. :

(12 ter) 	 $u_e \equiv 1$, $u_s \equiv Y(t-t_{BT})$ where Y is the Heavyside function.

and where $t_{BT}$ is the experimental break-through time, which gives a less realistic water output rate $\varphi_{1h}(t)$, having some parasit oscillation just after the B.T. However numerical experimentations (cf. [7]) have shown that this condition artificially increased the influence on $\varphi_{1h}$ of the values of $\alpha$ and $b$ for saturations less than 0.4. The results of simultaneous identification of $\alpha$ and $b$ in that case, all the remaining conditions being unchanged, are shown in Fig. 12.
Conclusion of the numerical results of the identification problem :

The above results show that the proposed optimization procedure works satisfactorily, but that the observation of the water output rate $\phi_{1h}(t)$ is not sufficient, when the realistic B.C. (12) are used, for the determination of $\alpha$ and $b$ over the whole range of saturations. It is suggested to ask the experimentators more observations, for instance the value, along time, of the saturation $u$ at a given point of the slab.

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12</figcaption></figure>

### 3. CONCLUSION

We have given an adequate discretization of a 1-D rather general non-linear degenerated diffusion and transport equation, valid under conditions ranging from the pure parabolic case to the pure hyperbolic case. We have calculated the corresponding adjoint state, and the exact gradient of the approximated output error criterion, and applied it to the resolution of an identification problem. Numerical results have been shown.

### REFERENCES

[1] M. MUSKAT The Flow of Homogeneous Fluids through porous Media, McGrawHill, New-York, (1937).

[2] G. CHAVENT A new formulation of diphasic incompressible flows in porous media, in Lecture Notes in Mathematics, Vol. 503, Springer, (1976).

[3] G. CHAVENT About the identification and modelling of miscible or immiscible displacement in porous media, in Lecture Notes in Control and Information Sciences, Vol. 1, Springer, (1978).

[4] G. CHAVENT Un théorème d'existence dans une inéquation variationnelle parabolique dégénérée modélisant le déplacement d'un fluide par un autre non miscible, LABORIA Report No. 260, IRIA, POBox 105, 78150 Le Chesnay, France.

[5] P. LESAINT and P.A. RAVIART On a finite element method for solving the neutron transport equation. Symposium on Mathematical aspects of finite elements in P.D.E. Math. Research Center, University of Wisconsin, Madison, (1974).

[6] P.A. RAVIART and J.M. THOMAS A mixed finite element method for second order elliptic problems.

[7] G. COHEN Thesis, Paris, (to appear).