STABILITY AND CONVERGENCE OF AN HYBRID FINITE VOLUME-FINITE ELEMENT METHOD FOR A MULTIPHASIC INCOMPRESSIBLE FLUID MODEL

. In this paper, we construct a fully discrete numerical scheme for approximating a two-dimensional multiphasic incompressible ﬂuid model, also called the Kazhikhov-Smagulov model. We use a ﬁrst-order time discretization and a splitting in time to allow us the construction of an hybrid scheme which combines a Finite Volume and a Finite Element method. Consequently, at each time step, one only needs to solve two decoupled problems, the ﬁrst one for the density and the second one for the velocity and pressure. We will prove the stability of the scheme and the convergence towards the global in time weak solution of the model.


Introduction.
1.1. The model. The Kazhikhov-Smagulov model, which can be deduced from the compressible Navier-Stokes system, describes the motion of a viscous, incompressible mixture of two fluids having different densities. The mixture is subject to a diffusion effect modeled by the Fick's law, which relates the velocity to the derivatives of the density. We assume that the fluid fills the domain Ω ⊂ R 2 , a bounded open set with sufficiently regular boundary Γ. We denote by n the unit outward normal on the boundary Γ and by [0, T ] the time interval, for T > 0. With the notations Q T = (0, T ) × Ω and Σ = (0, T ) × Γ, we consider the following model in Q T :    ∂ t ρ + div ρu = λ∆ρ, ρ ∂ t u + (u · ∇)u − λ(∇ρ · ∇)u + λ div ρ∇u T − µ∆u + ∇P = ρg, div u = 0. (1) The unknowns are ρ : Q T → R the density of the fluid, u : Q T → R 2 the incompressible velocity field and P : Q T → R the pressure of the fluid (a modified pressure). Moreover, g stands for the gravity acceleration (but it can include further external forces) and the parameters λ > 0 and µ > 0 represent mass diffusion and dynamic viscosity coefficients, respectively (which are assumed to be constant). Given a vector a ∈ R d , we set div a = d i=1 ∂ xi a i ; given a matrix valued function A, we denote div(A) the vector having components d j=1 ∂ xj A ij . This model was derived and analyzed for the first time by Kazhikhov and Smagulov [21].
STABILITY AND CONVERGENCE OF AN HYBRID FV-FE METHOD 431 3. For all φ ∈ C 1 ([0, T ]; V) such that φ(T, .) = 0, one has: T 0 − u, ρ∂ t φ + (ρu − λ∇ρ) · ∇ φ + µ ∇u, ∇φ − λ ρ∇u T , ∇φ dt Moreover, it is convenient to write the variational formulation of the problem. Let us assume that (ρ, u, P ) is a sufficiently regular solution of (1)-(2)- (3). Multiplying the equations (1) respectively by arbitrary test functions (ρ,ū,p) ∈ H 1 (Ω) × H 1 0 (Ω) × L 2 0 (Ω), integrating over Ω and using Green's theorem, adding to the momentum equation the density equation where we chooseρ = 1 2 u ·ū and we integrate by parts the convective and diffusive terms, finally we obtain the following formulation for a.e. t ∈ (0, T ): where we have used the following notations: • b ., ., . , a ., ., . and c ., ., . are the trilinear forms defined by: • d ., . is the bilinear form defined by: ∀v ∈ H 1 0 (Ω), ∀p ∈ L 2 0 (Ω). The trilinear forms verify the following properties of continuity, coercivity and antisymmetric: 1.3. Known results. Many authors treat the mathematical analysis of the Kazhikhov-Smagulov model in three-dimensional domains. We can refer for instance to [21,1,2]. In [21,1], under assumption (4) and if the constants λ, µ, m, M are such that λ < 2µ(M − m) −1 , the authors prove the existence of a weak solution of the problem (1)-(2)-(3), when u 0 ∈ H, ρ 0 ∈ H 1 (Ω) and g ∈ L 2 0, T ; L 2 (Ω) . In [2], the existence and uniqueness of a weak solution of (1)-(2)-(3) is proved under assumption (4) but without any restriction on the constant λ. In [23], under de condition λ/µ small enough, the existence and uniqueness of the global solution is proved in the two-dimensional case. Moreover, it is showed the convergence (as λ → 0) towards a weak solution of the Navier-Stokes system for nonhomogeneous fluids in two-and three-dimensional domains. Recently, in [4,5] the authors prove the existence of a regular (resp. strong) solution of a two-dimensional generalized Kazhikhov-Smagulov model. Concerning the numerical study, there exists few numerical schemes in order to approximate the problem (1). Some adequate choices can be found in [17,18,19,20,3], where the authors propose a fully discrete numerical scheme consists of C 0 finite element spatial approximation for all unknowns (density, velocity and pressure) combined with the backward Euler method in time. In [18] the authors obtain unconditional stability and convergence results for the two-dimensional case, by applying a truncating operator in the terms depending on the density. A conditionally stable and convergent numerical scheme is obtained in [19] for the three-dimensional case. In this work, the authors prove an approximate maximum principle bounding by excess and defect the discrete density with respect to the upper and lower bound of the initial density. Also, they study the asymptotic behavior of the numerical scheme as the diffusion parameter λ goes to zero, obtaining convergence towards a weak solution of the density-dependent Navier-Stokes problem. In [20], under hypotheses of regularity for the data and the exact solution, the authors present optimal error estimates of a linearized fully discrete scheme for the three-dimensional case. For the complete three-dimensional Kazhikhov-Smagulov model with O(λ 2 ) terms added in (1) 2 , the existence of regular solutions and some error estimates are given in [17], by assuming smallness conditions on the data. An extension of the results in [19] for the complete model with O(λ 2 ) terms is obtained in [3], where the treatment of the O(λ 2 ) terms requires special attention. Finally, another numerical scheme is developed in [12] by using a backward Euler scheme together with the method of characteristics for the volume fraction of the denser fluid, and a mixed finite element method in space for velocity and pressure.
The finite volume schemes are widely used for the numerical resolution of linear or non-linear conservation laws (see [13,24,27] for instance). The discrete maximum principle, implying the L ∞ -stability of a numerical scheme, is very important in the study of conservation laws. Indeed, the maximum principle is the first fundamental physical property of the problem that a suitable numerical scheme must faithfully reproduce. In [14], a convergence result for the numerical solution of a nonlinear convection-diffusion problem was investigated. The authors use a combined finite volume-finite element scheme, where the nonlinear convective terms are discretized by a monotone finite volume scheme and the diffusion term is approximated using conforming piecewise linear finite elements. In [15], the authors pursue this approach. The discrete maximum principle is necessary and it requires the use of triangulations of a weakly acute type. Under this assumption, the analysis of the error estimates of this combined finite volume-finite element scheme is achieved. Let us mention also the review paper by Droniou [11] which concerns various finite volume methods for solving diffusion equation on general meshes. The development of an hybrid finite volume-finite element scheme was firstly introduced in [7] and used in [6], in order to compute the numerical solution of the variable density incompressible Navier-Stokes system. Using a splitting in time, this hybrid scheme combines a finite volume method to discretize the density equation and a mixed finite element method to compute the velocity field and the pressure. In [6], the L ∞ -stability was obtained under an explicit CFL condition by introducing a second order finite volume scheme with multi-slope gradient reconstruction in order to solve transport equation on unstructured meshes with local refinements. Recently, the hybrid finite volume-finite element scheme was extended in [8] to the numerical simulation of powder-snow avalanches by solving the two-dimensional model (1).

1.4.
Main results. The objective of this work is to study the Kazhikhov-Smagulov model (1) in a fully discrete setting by coupling finite volumes to approximate the density and finite elements to approximate the velocity and pressure. Let h > 0, we denote by T h a partition of Ω composed of conforming and isotropic triangles. We (Ω) the finite elements spaces associated to density, velocity and pressure, respectively. In order to simplify the notations, we restrict our study to the case of a uniform time discretization of [0, T ], but all the results given in this work can be extended without any difficulty to the case of a general time discretization. Let N be a positive integer, then we define ∆t = T /N the time step and (t n = n∆t) N n=0 the partition of [0, T ]. Moreover, we consider the following stability condition: where c 0 > 0 is a constant which is independent of h and ∆t, but depends on the velocity field u ∈ V h . Obviously, (10) is a typical CFL condition often used for the numerical solution of conservation laws (see [24]). Let ρ n h , u n h ∈ W h × V h be the approximations of density and velocity at time t n . We denote by ρ h,∆t and u h,∆t the piecewise constant functions in time taking values ρ n h and u n h on (t n−1 , t n ], respectively. Thus, the following main result will be proved: The outline of the paper is organized as follows. In section 2 we describe the hybrid finite volume-finite element scheme, in particular we introduce the finite volume scheme and its properties in the vertex-based framework. In section 3 we study the stability of the numerical scheme. Then, we deduce the weak and weak * convergences results in section 4. Afterwards, we establish the compactness arguments for the discrete density and velocity that provide the strong convergence results in section 5. Finally, section 6 is devoted to the passage to the limit, concluding the proof of Theorem 1.2. 2. Description of the numerical schemes. This section is devoted to the development of an hybrid finite volume-finite element scheme. We recall that finite volume schemes are widely used for the numerical solution of conservation laws (see [13]), whereas finite element methods are naturally applied to approximate the solution of diffusive problems (i.e. elliptic or parabolic problems, see [10,16]). The idea of combined finite volume and finite element methods was used in various works concerning numerical computation of conservation laws. The idea behind this approach is to match the finite volume discretization of convective terms with the finite element discretization of diffusive terms. For the theoretical analysis of these combined schemes we refer to [15] and references therein. In our case, the hybrid finite volume-finite element scheme combines a finite volume approach to approximate the mass conservation equation and a finite element method to solve the momentum equation and the divergence free condition.
At first, we present the triangular mesh of Ω, the discrete spaces associated to density, velocity and pressure, and we set our assumptions on the discretization of the Kazhikhov-Smagulov model (1). Second, we focus on the development of the finite volume scheme to approximate the convection-diffusion equation, and we introduce the properties of the discrete density. Finally, we define the fully discrete scheme by using finite volume and finite element methods to approximate all unknowns (density, velocity, pressure) of Kazhikhov-Smagulov model.

Mesh definitions.
The discretization of the Kazhikhov-Smagulov model (1) will be carried out on an unstructured triangular mesh. Let h > 0, we denote by T h a partition of the polygonal domain Ω composed of conforming and isotropic triangles. The triangulation T h is called a basic (or primal) mesh. By h(T ) we denote the length of the longest side of the triangle T ∈ T h , and put We suppose the following assumptions for the triangulations {T h } h>0 of Ω [15,14]: The triangulations T h are of weakly acute type. In other words, the magnitude of all angles of all T ∈ T h is less than or equal to π/2. (A3) The triangulations T h verify the following inverse assumption: where c > 0 is constant independent of h. According to [10, Chap.3, § 3.1, Remark 3.1.3], assumptions (A1) and (A3) imply the existence of a constant c > 0 independent of h, such that where |T | = area of T ∈ T h . Now, let M h = {M i , i ∈ J} be the set of all vertices of the triangulation T h , (J is a suitable index set, of cardinality #J). The set E h of the edges of T h is made of straight segments [M i M j ] joining two vertices of M h . Let us construct the dual mesh C h = {C i , i ∈ J} over the basic mesh T h , which defines a second partition of Ω. The dual finite volume C i associated with each vertex M i ∈ M h is a closed polygon obtained in the following way: we join the barycenter of every triangle T ∈ T h which share the vertex M i with the middle point of every side of T containing M i . If M i ∈ M h ∩ ∂Ω, then we complete the boundary of C i by the segments joining M i with the middle point of boundary sides that contain M i . C i is often called the vertex-based control volume around the node M i . Accordingly, we have Moreover, if we denote |C i | = area of C i ∈ C h , then, we have For i ∈ J, let V(i) = {j ∈ J, C j is a neighbor of C i }. Let i ∈ J and j ∈ V(i), we define by T ij,1 and T ij,2 two neighboring triangles of T h sharing the common edge We denote n ij,1 (resp. n ij,2 ) the unit outward normal to C i along Γ ij,1 (resp. Γ ij,2 ) and |Γ ij,1 | (resp. |Γ ij,2 |) the length of the segment Γ ij,1 (resp. Γ ij,2 ). For Obviously, we have Consequently, there exists a constant c 1 > 0, such that Also, (13), (11) and (16) imply the existence of a constant c 2 > 0, such that 2.2. Discrete spaces. In order to combine the efficiency of the finite volume method for the convection-diffusion equation and the finite element method for the momentum equation, we need to define some discrete spaces associated to the unknowns. Let us define the following discrete spaces for the approximation of the density over the meshes T h and C h : Given a vector (β Mi ) i∈J ∈ R #J , there exists a unique Π h β ∈ W h and a unique L h β ∈ Z h such that As a consequence, there are one-to-one mappings between R #J , W h and Z h . Here, Π h : R #J → W h denotes the Lagrange interpolation operator and L h : C 0 (Ω) → Z h is the so-called lumping operator. Obviously, L h is a continuous linear operator.
For the approximations of velocity and pressure, there are several choices to define the mixed finite element spaces V h ⊂ H 1 0 (Ω) and Y h ⊂ L 2 0 (Ω) verifying the Ladyshenskaya-Babuska-Brezzi (LBB) condition (also called "inf-sup" condition) [16]. Following [7,8], we choose the Taylor-Hood element (P 2 × P 1 ), but others choices are also possibles (the mini-element P 1 -bubble × P 1 for instance). Then Throughout this work, we will suppose the following hypotheses [10,3,18,19]: (H1) Regularity for the data: Let u 0 ∈ V, ρ 0 ∈ H 1 (Ω) with 0 < m ≤ ρ 0 ≤ M in Ω and g ∈ L 2 0, T ; L 2 (Ω) . (H2) The triangulation T h of Ω and the finite elements space W h verify the following inverse inequality: (H3) Inf-sup condition: There exists a constant C > 0 independent of h, such that 2.3. The finite volume scheme. Here, we design the finite volume scheme for solving the convection-diffusion equation (1) 1 on the dual mesh C h . With a view to obtain the numerical scheme corresponding to (7) 1 , we will denote by ., . h (resp. . h ) the approximation of the scalar product (resp. norm) on L 2 (Ω), such that Also, the approximation ., . h can be defined with the aid of a numerical integration using the vertices of T ∈ T h as integration points: In particular, (24) and (25) correspond to the mass lumping technique applied to the mass matrix. Then, there exists constantsĉ 1 ,ĉ 2 > 0 such that ∀ h ∈ (0, h 0 ), For the bilinear form associated to the laplacian, we consider the control volume finite element (CVFE) scheme (see [13]). We denote by (φ i ) i∈J the canonical basis of W h characterized by The following geometrical properties hold: Then, we obtain that −a ii = j =i a ij > 0.
As consequence, ∀ρ, β ∈ W h we have Using assumption (A2) on the triangulation T h , we have which is a necessary condition for guaranteeing the (discrete) maximum principle of the finite volume scheme. Moreover, we will construct an approximation b h of the trilinear form b by using the finite volume approach, according to Feistauer et al. ideas [14,15]. We assume that the velocity field u ∈ V h is a given function which satisfies the divergence free constraint. Then, using (12), (20), Green's theorem and (14), we write for Finally, we characterize the approximation b h as follows. Given u ∈ V h , we define ∀ρ, where we have introduced in (27) an upstream numerical flux G ij l as defined in [7,6]. More specifically, we introduce the cell-averaged velocity u defined for each u ∈ V h : As shown numerically in [7], the choice of the cell-averaged velocity u is a necessary condition in order to ensure the divergence free constraint at the discrete level. Then, for l = 1, 2, i ∈ J, j ∈ V(i), the constant values u ij,l defined on ∂C i are such that With the value of u ij,l at hand, a simple upwind finite volume scheme lead to G ij l ρ 1 , ρ 2 , n ij,l = ρ 1 u ij,l · n ij,l if u ij,l · n ij,l > 0, ρ 2 u ij,l · n ij,l if u ij,l · n ij,l ≤ 0.
Obviously, upwind flux such as (28) restricts the method to first order accuracy. To improve the accuracy, more general fluxes can be designed based on MUSCL strategies. A multi-slope method on general unstructured dual meshes was introduced in [6], where the L ∞ stability for the convection equation was established. The discrete maximum principle can also be proved for the convection-diffusion equation.
As for first order methods, it requires the assumption (A2) on the triangulation T h (see [9]). With the numerical flux (28), some properties are ensured, such as the consistency, conservativity and monotonicity of the numerical flux. Also, G ij l (ρ 1 , ρ 2 , n ij,l ) is locally Lipschitz-continuous with respect to ρ 1 , ρ 2 : Thereafter, we use the following estimate of the approximation b h . Proposition 1. There exists a constant C > 0, such that for ρ, β ∈ W h and Proof. Let ρ, β ∈ W h and u ∈ V h . We write |Γ ji,l |G ji l ρ Mj , ρ Mi , n ji,l .
Using the conservativity of the numerical flux G ij l and the relations Γ ji,l = Γ ij,l , n ij,l = −n ji,l , we obtain |Γ ij,l |G ij l ρ Mi , ρ Mj , n ij,l .
Then, we find that The numerical flux is consistent, then for any constant function ρ, we have Let l = 1, 2, if i ∈ J and j ∈ V(i), then the segment [M i M j ] ∈ E h is the common side of triangles T ij,l ∈ T h , such that Γ ij,l ⊂ T ij,l (see Section 2.1). We have Moreover, the numerical flux G ij l is locally Lipschitz-continuous, then we obtain taking into account that each triangle T ∈ T h appears in the above sum as some T ij,l at most six times. Finally, in virtue of (11), the generalized Hölder inequality and the inverse inequality (23), we arrive at In order to complete the time discretization, we consider an Euler type method, which is implicit with respect to the diffusive term and explicit with respect to the convective term. In conclusion, we define the following finite volume scheme for the approximate solution of (1) 1 .
Initialization: Let ρ 0 h = Π h ρ 0 ∈ W h be approximation of the initial data ρ 0 , with Time step n + 1: Given ρ n h ∈ W h and u n h ∈ V h , find ρ n+1 h ∈ W h such that, for eachρ h ∈ W h : 2.4. Properties of the density. In this paragraph, we are interested to introduce some importants properties of the discrete density computed with the finite volume scheme (30). The existence and uniqueness of the solution ρ n+1 h of (30) follows from the well-known Lax-Milgram theorem. Notice that (30) is equivalent to the following linear problem for the unknown vector (ρ n+1 Mi ) i∈J ∈ R #J : First of all, it is essential to ensure that the previous finite volume scheme preserves the maximum principle. The following proposition clames the L ∞ stability of (30) on unstructured mesh verifying assumption (A2) under an appropriate CFL condition.
Proposition 2. Let u ∈ V h be the velocity field satisfying the divergence free constraint and let the initial density ρ 0 which verifies (4). If ∆t and h satisfy the condition: where c 3 > 0 is a constant independent of h and ∆t, then for each n, 0 ≤ n ≤ N −1, there exists a unique discrete solution ρ n+1 h of finite volume scheme (30) which verifies the pointwise estimates: in Ω.
The proof of Proposition 2 can be found in [14,13]. Evidently, the stability condition (10) and (17) Furthermore, we introduce the discrete Gagliardo-Nirenberg inequality for the density. The proof of this inequality can be found in [18,Lemma 10], using also (26).
Proposition 3. There exists C = C(Ω) > 0 (independent of h) such that, for any ρ h ∈ W h , one has: 2.5. The hybrid finite volume-finite element scheme. We aim to design a fully discrete numerical scheme in order to solve the Kazhikhov-Smagulov model (1) by applying an hybrid scheme which combines finite volume and finite element methods. Using a time splitting procedure, we are able to subdivide the global problem into decoupled sub-problems with respect to ρ and u, P at each time step. Concerning the spatial discretization, we use a finite volume scheme for the convection-diffusion equation for ρ (which was described in Section 2.3), and mixed finite elements for the linearized Navier-Stokes problem. The time discretization is obtained considering a backward Euler type scheme which is implicit with respect to the diffusion terms in both equations and explicit (resp. semi-implicit) with respect to the convective term in the density (resp. velocity) equation. Then, we define the numerical scheme as follows.
where we have used the following notation: In ( Then, using the discrete Laplacian operator ∆ h defined by (33), the finite volume scheme (35) can be rewritten as: In the next section, we shall see that the discrete problem (35)-(36)-(37) is wellposed, that is the existence and uniqueness of a solution holds. Also, we will establish the discrete version of the energy estimate for this hybrid scheme, independent of the discrete parameters.
3. Uniform estimates. In this section energy estimates for the velocity and strong estimates for the density will be obtained, using the discrete Laplacian of the density and involving the scheme (38), in order to prove the global stability of the hybrid scheme (35)-(36)-(37).
First of all, we start to establish some useful inequalities for the hybrid scheme (35)-(36)-(37).
where C 1 , C 2 ,ĉ 1 are positives constants independent of h, ∆t and n.
Proof. At first, we prove the inequality (39). We start takingū h = 2∆t u n+1 h in (36) andP h = P n+1 h in (37): First of all, we compute the first two terms of the above equality (41), using the identity (a − b, 2a) .
Then, using (8) and (9) in the third and fourth terms of (41), we obtain By applying the estimates (32), the Poincaré and Young inequalities, the last term of (42) is estimate as follows: and we get to (39). Next, to prove (40), we consider (38).
h , we find that: Taking into account the estimate (29), using the Young inequality and the discrete Gagliardo-Nirenberg inequality (34), we estimate the right-hand side term I of (43) as follows: Finally, using the Young inequality and the 2D Gagliardo-Nirenberg inequality for the velocity in (44), taking into account (26) in the left-hand side of (43), we obtain (40). Now, by virtue to show the global stability of scheme (35)-(36)-(37), it is easy to obtain, from Theorem 3.1 and from the discrete Gronwall's lemma, the following estimates for the velocity and more regularity for the density.
where C > 0 depends on the data ρ 0 , u 0 , g, λ but is independent of h and ∆t.
At last, as a consequence of Lemma 3.2, we deduce the following result: Under the assumptions of Lemma 3.2, the following estimate holds: where C λ > 0 is independent of h and ∆t.
Proof. This is a direct consequence result of Lemma 3.2 and Proposition 3.
Next, by taking into account the previous estimates given in Lemma 4.2, there exist subsequences (denoted in the same way) with the corresponding weak convergences towards limit functions u, u, ρ, ρ. Thanks to (45), we have the identities of the limits u = u and ρ = ρ. Therefore, we have the following result (for the proof see [ (denoted in the same way) and limit functions u, ρ verifying the following weak convergences, as (h, ∆t) → 0: 5. Strong convergence. As usual in this type of nonlinear system, to obtain the convergence of the hybrid scheme (35)-(36)-(37), we need strong convergence for the approximations in some suitable space in order to make the passage to the limit in the nonlinear terms. To do this, we must get the compactness for the discrete density and velocity.
At first, we establish a time derivative estimate for the discrete density.
Proposition 4. The following estimate holds: where C λ > 0 is independent of h and ∆t, but depends on the data (ρ 0 , u 0 , g, λ).
Proof. We start from (38), and using the Cauchy-Schwarz inequality, the estimate (29), the Sobolev embedding H 1 0 (Ω) ⊂ L 4 (Ω), and (26), we find that Having in mind that ρ n+1 h − ρ n h ∆t ∈ W h ⊂ L 2 (Ω) and using a continuity argument, we deduce from (46) and (26) Finally, by applying the Hölder inequality, the estimates given by Lemma 3.2 and Corollary 1, we get straightforwardly the desired result.
Furthermore, from Lemma 4.2 the discrete density is bounded in L ∞ (Q T ), then one also gets the strong convergence in L p 0, T ; L p (Ω) for any p < ∞. For p = ∞, one can deduce the convergence for at least a subsequence of ρ h,∆t or ρ h,∆t to ρ for a.e. (t, x) ∈ Q T . Now, we are going to estimate the fractional time derivative for the discrete velocity.
Proof. We follow exactly the proof done in [18,Proposition 27], where [ρ n h ] T is replaced by ρ n h for each n = m, . . . m + r and the equation [18, (66)] is replaced by (51) If in the right side of (51) we use the Cauchy-Schwarz inequality, the estimates (29) and (26) Finally, based on the previous result and thanks to a compactness theorems of Aubin-Lions type (see [25,Theorem 5]), we obtain the strong convergences for the velocity.
(52) 6. Passing to the limit. The final step to complete this study is to employ the convergence results obtained throughout this work in order to pass to the limit in the discrete problem. Our goal is to show that the approximate solution ρ n+1 h , u n+1 h , obtained with the aid of the hybrid finite volume-finite element scheme (35)-(36)-(37), converges in some sense to the weak solution ρ, u of the Kazhikhov-Smagulov model (1)-(2)-(3), when the space and time parameters (h, ∆t) tend to zero.
In the case of a convection-diffusion equation, the passage to the limit for the finite volume scheme (35) (as (h, ∆t) → 0) has been done in details by Feistauer et al. (see [14, section E]). But in our case, the velocity in the convective part of (7) 1 depends by (h, ∆t). To handle this difficulty, we use the following strong-weak convergence result.