Some remarks on simplified double porosity model of immiscible incompressible two-phase flow

1. Introduction

Naturally fractured reservoirs are characterized by a system of fractures existing within a background rock matrix. The fracture system has a low storage capacity and a high conductivity, while the matrix block system has a conductivity that is low in comparison with that in the fractures. The majority of fluid transport will occur along flow paths through the fissure system, and the relative volume and storage capacity of the porous matrix is much larger than that of the fissure system (type 2 reservoirs in [26]). Multiphase flow in subsurface fractured media offers a particular challenge to numerical modeling, since both the flow through the fracture network and transfer to a relatively stagnant matrix need to be modeled. The fractures have a very strong influence on flow and transport but the discrepancy between the width of the fractures and other dimensions involved makes the inclusion of fractures in a numerical model difficult and costly. Possible applications of such systems are in improved oil recovery in hydrocarbon reservoirs, non-aqueous phase contaminant transport, nuclear waste containment etc. (see, e.g., [12]).

Mathematical models of multiphase flow in fractured-porous media can be categorized into two classes: the dual-continuum models and the discrete-fracture network models. Discrete fracture modeling represents each fracture and the matrix as a geometrically well-defined entities which are represented explicitly in the computational domain. This approach leads to the most accurate and physically realistic model at the expense of highly refined or hybrid grid [21,30].

In dual-continuum model one do not consider specific known fractures that might be included individually in the model but a network of small interconnected fractures with certain degree of regularity. There are two large differences in scale in the fractures and the blocks: fracture width is very small compared to the scale of the domain and fracture permeability is much larger than the permeability of surrounding material. In dual-continuum approach a sort of averaging process is applied in order to obtain simplified description of matrix-fracture system and its interactions. In the case of stagnant flow in the matrix averaging leads to a double porosity model that was first obtained in [11,33], where the model was described as a phenomenological model deduced experimentally. In the framework of the model presented in these papers the fissures are assumed to have a negligible volume with respect to the volume of the whole reservoir (see [11] or [29] Chapter 5).

From the mathematical point of view, the usual double-porosity model (or ε 2 -model) assumes that the width of the fracture containing highly permeable porous media is of the same order as the block size. The related homogenization problem was first studied in [9], and was then revisited in the mathematical literature by many other authors (see, e.g. [1,3,4,14,16–18,20,23,25,27,31], and references therein).

Homogenization of linear and nonlinear parabolic equations, in particular, parabolic equations with high contrast coefficients (see the references above) is a long-standing question. In this paper we deal with periodic homogenization of double-porosity-type problems in a special case of asymptotically small volume fraction of highly conductive part of the medium. In other words, when the permeability coefficient of the matrix part vanishes everywhere in the corresponding domain except a set of asymptotically small measure. The first results on the subject were obtained independently and by different approaches in [13] and [28]. In [28] the highly conductive part of the medium (the fissures system) is modelized by the only one small parameter ε which characterizes the scale of the microstructure. The main feature of the homogenized model in this case is that it does not involve neither a cell problem for the calculation of the additional source term nor a cell problem for the calculation of the global permeability tensor. The further results on the subject are connected with the application of the method of two small parameters ε , δ proposed in [10,19]. In this case the parameter ε describes the periodicity of the fractured-porous medium and δ is the relative thickness or opening of the fracture. The homogenization process then splits in two steps. First, one shows that the corresponding problem admits homogenization as ε → 0 and then in the second step it is necessary to pass to the limit as δ → 0 . Thus the effective system does not depend on ε , δ . We refer here to [5,6] and the references therein. The first result on the homogenization of the two-phase flow in double porosity media was obtained in [22] by the method of two small parameters. Further application of the method of two small parameters can be found in [32] where a homogenized model of the two-phase incompressible non-equilibrium flow in fractured porous media in the framework of Kondaurov’s model (where the mobilities and capillary pressure depend both on the real saturation and a non-equilibrium parameter) was derived. In the global double porosity δ-problem obtained in [22] after passing to the limit as ε → 0 , additional matrix-fracture source terms are present that are given implicitly via solutions of a nonlinear local boundary value problem known as the imbibition equation. The nonlinearity of the imbibition equation causes difficulties in numerical simulations of the model since no analytic expression of the matrix-fracture source term is available. In order to overcome this issue, one can linearize the imbibition equation and then express the matrix-fracture source term explicitly from the linearized equation. In [22] the imbibition equation is linearized by using an appropriate constant, as suggested by [7]. In this paper we present a new, variable and more general linearization of the imbibition equation which gives a new simplified double porosity model. We also display numerical simulations comparing the matrix-fracture exchange term calculated by solving nonlinear imbibition equation to the matrix-fracture exchange terms given by two different linearization procedures.

The rest of the paper is organized as follows. In Section 2 we present a double porosity model of the two-phase flow in the case of thin fissures, where the opening of the fissure is described by a small parameter δ. The main result of the paper is contained in Section 3 in which we propose a new way to linearize the imbibition equation and express the corresponding matrix-fracture source term in the form of a convolution. In the rest of Section 3 we present the new simplified or fully homogenized double porosity model that is obtained from the double porosity δ-model coupled with the newly linearized imbibition equation as the thickness of the fractures δ tends to 0. In Section 4 we present numerical simulations comparing the matrix-fracture exchange term calculated by solving nonlinear imbibition equation to the matrix-fracture exchange terms given by different linearization procedures.

2. Double porosity model

In this section we present a double porosity model of incompressible two-phase flow derived rigorously by the homogenization theory in [8,15,34]. Namely, we consider the reservoir Ω ⊂ R d , d = 1 , 2 , 3 , of the characteristic length L composed of the matrix blocks with the characteristic length l and highly permeable network of fractures. The block size is small compared to the size of the flow domain, i.e., ε = l / L ≪ 1 is a small parameter which goes to zero. The thickness of the fractures is supposed to be of order l δ , where 0 < ε ≪ δ < 1 is a second small parameter. The porosities of the blocks and the fractures are supposed to be constant and are denoted by Φ m and Φ f respectively. The permeabilities of the blocks and the fractures are highly contrasted. In the derivation of the double porosity model it is assumed that if the fracture porosity is k f , then the matrix porosity is ( l δ ) 2 k m , where k f and k m are of the same order.

If we neglect gravitational segregation the double porosity model obtained by homogenization as ε → 0 (see [8,15,34]) can be written in the form (see [22]): (1) Φ ∗ δ ∂ S f δ ∂ t − div ( K ∗ δ λ w , f ( S f δ ) ∇ P w , f δ ) = 𝒬 w δ − Φ ∗ δ ∂ S f δ ∂ t − div ( K ∗ δ λ n , f ( S f δ ) ∇ P n , f δ ) = 𝒬 n δ , P c , f ( S f δ ) = P n , f δ − P w , f δ , where S f δ , P w , f δ and P n , f δ are wetting phase saturation and pressure and non wetting phase pressure in the fractures, respectively; P c , f , λ w , f and λ n , f are the capillary pressure function and the phase mobility functions in the fractures; Φ ∗ δ and K ∗ δ are the effective porosity and permeability of the matrix-fracture system. The terms 𝒬 w δ and 𝒬 n δ are the wetting phase and the non wetting phase source terms modeling the phase mass transfer from the matrix to the fracture system governed by the capillary imbibition.

In this paper as in [22], we adopt the Warren–Root idealization of the fractured media, i.e., each matrix block is similar to a rectangular one and all the blocks are surrounded by the fractures. Then we introduce the reference cell Y = ( 0 , 1 ) d which is decomposed in matrix and fracture parts, where the matrix part, Y m δ , is an open cube with edge length 1 − δ ( 0 < δ ≪ 1 ), and Y f δ = Y ∖ Y m δ represents the fracture part. The flow domain Ω is assumed to be covered by a pavement of cells l Y .

The effective porosity Φ ∗ δ can be expressed as Φ ∗ δ = Φ f vol ( Y f δ ) vol ( Y m δ ) . Moreover, Φ ∗ δ is of order δ. The effective permeability tensor can be expressed using the solutions of certain cell problems (see [22] for more details) and it is again of order δ.

The matrix-fracture source terms 𝒬 w δ , 𝒬 n δ are given by: (2) 𝒬 w δ ( x , t ) = def − Φ m vol ( Y m δ ) ∫ Y m δ ∂ S m δ ∂ t ( x , y , t ) d y = − 𝒬 n δ ( x , t ) , where the function S m δ ( x , y , t ) is the matrix block saturation defined for each point x ∈ Ω as a solution of the problem (3) Φ m ∂ S m δ ∂ t − δ 2 k m Δ y β m ( S m δ ) = 0 in Y m δ , S m δ ( x , y , t ) = P ( S f δ ( x , t ) ) on Γ δ , S m δ ( x , y , 0 ) = S m 0 ( x ) in Y m δ , where Γ δ stands for the interface between the matrix and fracture parts of the cell Y; P ( S ) = def P c , m − 1 ( P c , f ( S ) ) and (4) β m ( s ) = def ∫ 0 s α m ( ξ ) d ξ , where α m ( s ) = def λ w , m ( s ) λ n , m ( s ) λ m ( s ) | P c , m ′ ( s ) | . Here λ m ( s ) = λ w , m ( s ) + λ n , m ( s ) . Equation (3) 1 is known as the imbibition equation.

We note that the matrix-fracture source terms 𝒬 w δ , 𝒬 n δ given by (2) can also be calculated as: 𝒬 w δ ( x , t ) = − δ 2 k m vol ( Y m δ ) α m ( P ( S f δ ( x , t ) ) ) ∫ ∂ Y m δ ∇ S m δ · n d y = − 𝒬 n δ ( x , t ) .

3. Linearization of the imbibition equation

The main difficulty in application of the double porosity model for the two-phase flow is the fact that the imbibition equation is nonlinear. The nonlinearity does not allow an analytic solution of the equation and an analytic expression of the corresponding source term. This implies that in the numerical simulation of the double porosity model we have to solve the imbibition equation many times. This problem is especially difficult in the case of thin fractures, where the solution of the imbibition equation is of a boundary layer type and a very refined mesh is needed to resolve it. In order to overcome this difficulty, it is possible to linearize the imbibition equation and then use the linearized equation to express the matrix-fracture transfer term by an analytic expression.

In this section we explore different ways of linearization of the imbibition equation. The goal is to derive a linear equation with constant coefficients that can be solved analytically.

Linearization, as suggested in [7], consists in replacing a nonlinear function α m ( S m δ ) in the term div ( α m ( S m δ ) ∇ S m δ ) by some function α ‾ m such that α ‾ m ≈ α m ( S m δ ) for all values S m δ of interest. The simplest variant of such linearization is replacing α m ( S m δ ) by its mean value, namely, we consider a constant α ‾ m > 0 such that (5) α ‾ m = ∫ 0 1 α m ( s ) d s , and we replace the imbibition equation (3) by its linearized version (6) Φ m ∂ S ˆ m δ ∂ t − δ 2 k m α ‾ m Δ y S ˆ m δ = 0 in Y m δ , S ˆ m δ ( x , y , t ) = P ( S f δ ( x , t ) ) on Γ δ , S ˆ m δ ( x , y , 0 ) = S m 0 ( x ) in Y m δ . Now instead of calculating the matrix-fracture source terms defined in (2) by the solution of the original imbibition equation (3) one can calculate these terms using the linearized imbibition equation (6). In this case we denote them by: 𝒬 ˆ w δ = − 𝒬 ˆ n δ . Due to the linearity, these terms can be expressed as a convolution. Namely, 𝒬 ˆ w δ ( x , t ) = − ∂ ∂ t ∫ 0 t K m δ ( t − u ) ( P ( S f δ ( x , u ) ) − P ( S f 0 ( x ) ) ) d u = − ∫ 0 t K m δ ( t − u ) ∂ ∂ t P ( S f δ ( x , u ) ) d u = − 𝒬 ˆ n δ ( x , t ) , where the kernel K m δ ( t ) can easily be calculated (see, e.g., [2]). Due to simplicity of the domain Y m δ the kernel K m δ can be expanded in a function series, truncated at some point and used in numerical calculations. In this way, it enables to avoid the resolution of the local problem numerically. However, the coupling between the local and global problems is still present through the kernel K m δ ( t ) . We will not give details of this approach since we are interested in the case of thin fractures. In Remark 2 we will decouple the local and the global problem by passing to the limit as δ → 0 and obtain a convolution form of the matrix-fracture source terms with the kernel given explicitly.

We propose now a new, more general way to linearize the imbibition equation (3). Namely, let us consider the following boundary value problem: (7) Φ m ∂ S ˜ m δ ∂ t − δ 2 k m α ˆ m δ ( x , t ) Δ y S ˜ m δ = 0 in Y m δ , S ˜ m δ ( x , y , t ) = P ( S f δ ( x , t ) ) on Γ δ , S ˜ m δ ( x , y , 0 ) = S m 0 ( x ) in s Y m δ , where the coefficient α ˆ m δ ( x , t ) can be chosen in different ways. One particularly useful choice of the coefficient α ˆ m δ ( x , t ) is the average of the function α m over the range of saturation given by the boundary conditions: α ˆ m δ ( x , t ) = ∫ P ( S min δ ( x , t ) ) P ( S max δ ( x , t ) ) α m ( s ) d s / ( P ( S max δ ( x , t ) ) − P ( S min δ ( x , t ) ) ) (8) where S min δ ( x , t ) = min 0 ⩽ s ⩽ t S f δ ( x , s ) , S max δ ( x , t ) = max 0 ⩽ s ⩽ t S f δ ( x , s ) . Then using the function β m we can express α ˆ m δ as follows: α ˆ m δ ( x , t ) = β m ( P ( S max δ ( x , t ) ) ) − β m ( P ( S min δ ( x , t ) ) ) P ( S max δ ( x , t ) ) − P ( S min δ ( x , t ) ) . Note that in the case of increasing (injection of water in oil saturated media) or decreasing fracture saturation we have a more simple expression for α ˆ m δ : α ˆ m δ ( x , t ) = ∫ P ( S f δ ( x , 0 ) ) P ( S f δ ( x , t ) ) α m ( s ) d s / ( P ( S f δ ( x , t ) ) − P ( S f δ ( x , 0 ) ) ) (9) = β m ( P ( S f δ ( x , t ) ) ) − β m ( P ( S f δ ( x , 0 ) ) ) P ( S f δ ( x , t ) ) − P ( S f δ ( x , 0 ) ) .

With any choice of the function α ˆ m δ ( x , t ) ⩾ 0 problem (7) can be reduced to (6) with α ‾ m = 1 , and thus to an equation with constant coefficients, by the change of the time variable (10) τ x δ ( t ) = def ∫ 0 t α ˆ m δ ( x , s ) d s and passing to a new unknown function S ˆ m δ defined by (11) S ˜ m δ ( x , y , t ) = S ˆ m δ ( x , y , τ x δ ( t ) ) . One can easily check that the function S ˆ m δ is a solution to problem (6) with α ‾ m = 1 and corresponding boundary condition. Notice that the function t ↦ τ x δ ( t ) in (10) is invertible except on the time intervals where α ˆ m δ ( x , t ) = 0 . On these intervals, in the case of α ˆ m δ given by (9), the solution S ˜ m δ of (7) and the boundary function S f δ do not depend on time, making (11) consistent on these intervals and making S ˆ m δ well defined.

From the definition of the matrix-fracture source terms (2), using our linearized imbibition equation (7) and the change of variables (10), (11) we obtain the following (linearized) form of the matrix source term: (12) 𝒬 ˜ w δ ( x , t ) = − α ˆ m δ ( x , t ) Φ m vol ( Y m δ ) ∫ Y m δ ∂ S ˆ m δ ∂ τ ( x , y , τ x δ ( t ) ) d y = α ˆ m δ ( x , t ) 𝒬 ˆ w δ ( x , τ x δ ( t ) ) , where 𝒬 ˆ w δ is the matrix-fracture source term produced by linear imbibition equation (6).

Since the term 𝒬 ˆ w δ ( x , τ ) can be expressed in the convolution form with the kernel that is easily calculable, then we can exclude the calculation of the solution of the local problem (7) just as in the case of the constant approximation (5), (6).

Remark 1.

Another possible linearization consists in rewriting equation (3) in terms of the function U δ = β m ( S m δ ) with nonlinearity now appearing in the time derivative term. One can then apply the same type of linearizations as above which lead to different fracture-matrix source terms. In the numerical simulations with this kind of linearization we have observed different results but without net improvement over linearization we have presented above, so we do not consider that kind of linearization further.

Remark 2.

In double porosity model given by (1), (2) and (3) there remains a small parameter δ which measures relative thickness of the fractures. By letting this small parameter to zero one can obtain full decoupling of the local and the global problem. In the case of the two-phase flow that problem was studied in [22]; see the references therein for the case of one-phase flow. It is shown in [22] that in the limit δ → 0 the effective fracture equations (1) coupled with the linearized imbibition equation (6) reduces to the system: (13) Φ f ∂ S ∂ t − div ( k ∗ λ w , f ( S ) ∇ P w ) = 𝒬 w − Φ f ∂ S ∂ t − div ( k ∗ λ n , f ( S ) ∇ P n ) = 𝒬 n , P c , f ( S ) = P n − P w , where S, P n and P w are effective variables in the system of fractures, k ∗ = k f ( d − 1 ) / d and source terms 𝒬 w = − 𝒬 n are given by (14) 𝒬 w ( x , t ) = − C m ∂ ∂ t ∫ 0 t P ( S ( x , u ) ) − P ( S ( x , 0 ) ) t − u d u , (15) C m = 2 d Φ m k m α ‾ m / π . We note that in this final reduction the effective porosity is replaced by the fracture porosity, the effective permeability is reduced fracture permeability, and the matrix-fracture source term 𝒬 w δ ( t ) reduces to the convolution expression (14) with the kernel K ( t ) = C m / t . In the model (13), (14) effective permeability and the matrix-fracture source term are given explicitly and there is no need to solve local problems. The local and the global problems are completely decoupled. We should also mention that in the system (1), (6) the porosity, permeability and the matrix-fracture source term are proportional to δ, and equations (13), (14) are obtained after division by δ. For the source terms we therefore have 𝒬 ˆ w δ = δ 𝒬 w + o ( δ ) (see [22] for more details). In the case of variable linearization of the imbibition equation given by (7) we can use the matrix-fracture source representation (12) to compute, at least formally, the asymptotic behavior of 𝒬 ˜ w δ when δ → 0 . To this end we assume that the boundary function S f δ ( x , t ) converges to some S f ( x , t ) leading to convergence α ˆ m δ ( x , t ) → α ˆ m ( x , t ) and τ x δ ( t ) → τ x ( t ) , where obviously (16) α ˆ m ( x , t ) = β m ( P ( S f ( x , t ) ) ) − β m ( P ( S f ( x , 0 ) ) ) P ( S f ( x , t ) ) − P ( S f ( x , 0 ) ) , and (17) τ x ( t ) = def ∫ 0 t α ˆ m ( x , s ) d s . Now formula (12) gives 𝒬 ˜ w δ ( x , t ) = α ˆ m δ ( x , t ) 𝒬 ˆ w δ ( x , τ x δ ( t ) ) , where by asymptotic expansion of the matrix-fracture source term in the case of constant linearization we have 𝒬 ˆ w δ ( x , τ ) = − δ C m ∂ ∂ τ ∫ 0 τ P ( S f ( x , ( τ x ) − 1 ( u ) ) ) − P ( S f ( x , 0 ) ) τ − u d u + o ( δ ) , and C m = 2 d Φ m k m / π . Therefore we get (18) 𝒬 ˜ w δ ( x , t ) = − δ C m ∂ ∂ t ∫ 0 τ x ( t ) P ( S f ( x , ( τ x ) − 1 ( u ) ) ) − P ( S f ( x , 0 ) ) τ x ( t ) − u d u + o ( δ ) . Further change of variables u = τ x ( s ) gives 𝒬 ˜ w δ ( x , t ) = − δ C m ∂ ∂ t ∫ 0 t P ( S f ( x , s ) ) − P ( S f ( x , 0 ) ) τ x ( t ) − τ x ( s ) α ˆ m ( x , s ) d s + o ( δ ) . We can now conclude that in the limit δ → 0 the effective fracture equations (1) coupled with the linearized imbibition equation (7) reduces to the system (13) with the effective matrix-fracture source term 𝒬 w = − 𝒬 n given by (19) 𝒬 w ( x , t ) = − C m ∂ ∂ t ∫ 0 t β m ( P ( S f ( x , s ) ) ) − β m ( P ( S f ( x , 0 ) ) τ x ( t ) − τ x ( s ) d s , where τ x is given by (16), (17) and C m = 2 d Φ m k m / π . We have obtained 𝒬 w in (19) by leaving out small o ( δ ) terms and dividing by δ in (18).

4. Numerical comparison of matrix-fracture exchange terms

In this section we will compare by means of numerical simulation the matrix-fracture exchange term calculated by solving nonlinear imbibition equation (3) to the matrix-fracture exchange terms given by different linearisation procedures. Since we do not consider the whole two-phase flow simulation in this section we will impose an artificial boundary condition on the block which corresponds to water injection into oil saturated media.

In this example we consider porous block Y = ( 0 , 1 − δ ) d and we will consider relative fracture thickness δ = 0.3 , 0.1 , 0.01 , 0.001 . For fracture permeability we will take fixed value k f = 10 − 13 m² and we take k m = k f . This means that the matrix permeability in the scaled matrix block Y is equal to δ 2 k f . Furthermore, we take the matrix porosity Φ m = 0.35 and fluid viscosities μ w = 10 − 3 Pa·s and μ n = 2 · 10 − 3 Pa·s.

Van Genuchten–Mualem model is used to express the capillary pressure functions and the relative permeabilities. We have P c ( S w ) = P r ( S w − 1 / m − 1 ) 1 / n , k r w ( S w ) = S w [ 1 − ( 1 − S w 1 / m ) m ] 2 , k r g ( S w ) = 1 − S w ( 1 − S w 1 / m ) 2 m , where S w is effective wetting phase saturation and the residual saturations are taken to be zero; P r is reference pressure for Van Genuchten law and n > 1 , m ⩾ 0 are such that m = 1 − 1 / n .

Evolution of the fracture saturation is given by (20) S w f ( t ) = 0.05 + min ( t / 10 , 0.9 ) , and is shown on Fig. 1 together with the boundary condition P ( S w f ( t ) ) . On the same figure we also show the capillary pressure functions in the fracture P c , f and in the matrix P c , m as well as the functions α m , β m and P .

OPEN IN VIEWER

Fig. 1.

Functions used in Simulation 1. Van Genuchten parameters are: P r = 1 bar and n = 2 in the matrix and P r = 0.1 bar and n = 2 in the fractures.

For numerical resolution of the imbibition equation we need a grid that is well adapted for resolving the boundary layers that governs the solution. For that purpose we have adopted meshes of Bakhvalov type (see [24]). The adequacy of chosen parameters of the Bachvalov grid is verified in two ways. First, in 1-D we compared numeric solution to the analytic solution which is easy to calculate in the linear cases. Secondly, since the mass transfer term can be calculated by volume integration and by boundary integration we refined the grid up to point where the two methods give results that differ no more than 1%.

OPEN IN VIEWER

Fig. 2.

Evolution of the matrix-fracture exchange term 𝒬 w δ ( t ) divided by δ for different values of δ, for the nonlinear and two linear models. On x-axis is given time in days.

On Fig. 2 we present time evolution of the matrix-fracture transfer term 𝒬 w δ ( t ) for the nonlinear model (3) (denoted by nlin) and for two linear models: one given by (6) (denoted by clin) and model (7) with (9) (denoted by vlin). Figure 2 shows the matrix-fracture transfer term 𝒬 w δ ( t ) divided by δ for different values of fracture thickness δ. It is shown that the approximation to nonlinear matrix-fracture transfer term given by linear model (7) and (9) is much better than approximation given by model (6) and that the quality of the approximation is rather independent of the fracture size δ. In fact for small δ, expression 𝒬 w δ ( t ) / δ becomes quickly actually independent of δ, which confirms theoretical result in [22].

OPEN IN VIEWER

Fig. 3.

The capillary pressure functions and saturation transfer curve in the case of large difference in the matrix and the fracture. Van Genuchten parameters are: P r = 10 bar and n = 2 in the matrix and P r = 0.1 bar and n = 2 in the fractures.

The matrix-fracture exchange term is strongly influenced by difference between the capillary pressure curves in the matrix and the fracture. If the difference between the two curves is large the saturation transfer function P will have strong derivative near S w = 0 and it will be almost constant in the rest of the domain. This will strongly influence the boundary condition for the imbibition equation. In Fig. 3 we show the case of the van Genuchten capillary pressure functions with the parameters n = 2 and P r = 10 bars in the matrix and with n = 2 and P r = 0.1 bars in the fracture. In the other extreme, where the two capillary pressure functions are mutually equal, the saturation transfer function P is linear.

In Fig. 4 we compare the matrix-fracture exchange term 𝒬 w δ ( t ) in the case of strongly different capillary pressure functions shown in Fig. 3, and the case of equal capillary pressure functions ( P r = 1 bar, n = 2 ) in the matrix and the fracture; in both case δ = 10 − 3 . We see in both cases that better approximation to the nonlinear matrix-fracture exchange term is again given by vlin curve, that is by the model given by (7) and (9). Simpler clin approximation given by model (6) gives in all cases less good approximation.

OPEN IN VIEWER

Fig. 4.

In the left column strong difference in the matrix and fracture P c functions. Parameters are given in Fig. 3; in the right column case of equal P c functions. In the first row δ = 10 − 2 and in the second row δ = 10 − 3 .

Finally we chose an example in which the fracture saturation is not monotone. In that case we need to use the definition of the coefficient α ˆ m δ ( x , t ) given by (8). Keeping all the other parameters as before we change only the boundary conditions which is now given by (21) S w f ( t ) = 0.5 + 0.5 sin ( π t / 5 ) , and the simulation time is, as before, 10 days. The function (21) and corresponding boundary condition P ( S w f ( t ) ) are shown on Fig. 5. The matrix-fracture exchange term 𝒬 w δ ( t ) is shown on Fig. 6. We see that after lost of the monotonicity of the boundary condition the matrix-fracture exchange term given by the variable linearization looses its precision but stays comparable to the constant linearization version of the matrix-fracture exchange term.

OPEN IN VIEWER

Fig. 5.

Non monotone boundary saturation.

OPEN IN VIEWER

Fig. 6.

Matrix-fracture exchange term in the case of non monotone boundary condition and δ = 10 − 3 .

The main difficulty in numerical simulations of the double porosity model (1)-(4) is that the boundary value problem for the nonlinear imbibition equation (3)-(4) needs to be solved many times, as many as there are spatial mesh nodes in the domain Ω. Furthermore, in the case of a fractured porous media with thin fractures the solution of the imbibition equation turns out to be of a boundary layer type and one needs a very refined mesh to resolve it. Up to our knowledge there is no nonlinear simplified model which supports the idea of a linearization. In this work we have proposed a new, more general way to linearize the imbibition equation (3) which appears in the definition of the matrix-fracture transfer source terms in the global double porosity δ-model of incompressible two-phase flow in porous media. After passage to the limit as δ → 0 , we analyze the effective matrix-fracture exchange source term obtained by this new linearization and compare it to the effective matrix-fracture exchange source terms obtained previously by a constant linearization in [22]. Numerical simulations are provided which show that the matrix-fracture exchange term based on the new linearization procedure gives a better approximation of the exact one than the corresponding exchange term obtained earlier by the authors. Thus by introducing the new, variable linearization of the imbibition equation we have also derived the new simplified double porosity model that is more accurate than the simplified double porosity model obtained earlier in [22] by the constant linearization of the imbibition equation. The numerical simulations of the simplified double porosity model are a part of our future research.