Homogenization of a Biot-Stokes system modeling deformable vuggy porous media
Zhaoqin Huang, Xu Zhou, Tao Huang, Jun Yao, Xiaoguang Wang, Hervé Jourde
11 Homogenization of a Biot-Stokes system modeling deformable vuggy porous media
Zhaoqin Huang , Xu Zhou , Tao Huang , Jun Yao , Xiaoguang Wang , Hervé Jourde
1. Research Center of Multiphase Flow in Porous Media, China University of Petroleum (East China), Qingdao Shandong, 266580, China 2. Laboratoire HydroSciences Montpellier (HSM), UMR 5569 CNRS-IRD-UM, Université de Montpellier, 163 Rue Auguste Broussonnet, Montpellier 34095, France *Corresponding author: [email protected], [email protected]
Abstract
Vugs are small to medium-sized cavities commonly present in rocks. The presence of vugs may have non-trivial impacts on the hydro-mechanical behaviors of the rock. How to effectively quantify and analyze such effects is still an challenging problem. To address the problem, we derive a macroscopic coupled hydro-mechanical model for single-phase viscous fluid flow through a deformable vuggy porous medium. We first model the hydro-mechanical coupling process on a fine scale using Biot’s equations within porous matrix, Stokes equations within the vugs, and an extended Beavers-Joseph-Saffman boundary condition on the porous-fluid interface. Based on the homogenization theory, we then obtain macroscopic Biot’s equations that govern the hydro-mechanical coupling behaviors of vuggy porous media on a large scale. Subsequently, the macroscopic poroelastic coefficients, such as the effective permeability, effective Young’s modulus and effective Biot coefficient, are derived from three cell problems. Finally, several numerical examples are designed to validate the proposed model and to demonstrate the computational procedure for evaluation of the hydro-mechanical behaviors of vuggy porous media.
Keywords
Biot-Stokes system; Homogenization; Vuggy porous media; Poroelasticity; Extended Beavers-Joseph-Saffman boundary condition.
1. Introduction
Vug are defined as visible pores that are significantly larger than adjacent grains or crystals (Lucia, 2007; Huang et al., 2010, 2011), as shown in Figure 1. The presence of vugs can significantly increase the porosity and permeability of the rock (Arbogast and Lehr, 2006; Popov et al., 2009; Karasuyama et al., 2011; Xie et al., 2017). The main challenge is the coupling of porous flow and free-fluid flow (Arbogast and Brunson, 2007; Zhang et al., 2016), which is the key issue for many environmental, biomedical and industrial applications (Mosthaf et al., 2011; Baber et al., 2012; Chen et al., 2014). There are usually two schemes to model the coupling porous-free flow. One is the single domain approach using Stokes-Brinkman equations, the other one is the two domain approach using Darcy-Stokes equations. The former provides a model that can be continuously varied from a Darcy dominated flow to a Stokes dominated flow, so it avoids to introduce extra boundary conditions between porous media and free fluid region (Popov et al., 2009). However, the recent results show that the Stokes-Brinkman equations are only applicable for a special parameter set (Karasuyama et al., 2011), where the effective vugs permeability should be less than four orders of magnitude different from the matrix permeability. Therefore, the Darcy-Stokes system is more used. But one should introduce suitable boundary conditions at the porous-fluid interface (Goyeau et al., 2003) to couple these two equations. The Beavers-Joseph-Saffman (BJS) boundary condition is usually adopted (Beavers and Joseph, 1967; Saffman, 1971; Mikelic and Jäger, 2000; Carraro et al., 2015). It should be note that the mentioned coupling equations are just focus on the fluid flow and ignore the deformation of the porous solid skeleton. porousmatrix vugs0 1 2 3 cm CT scan slice
Figure 1: A typical core of vuggy carbonate reservoirs.
Recently, the interaction between the free fluid and a deformable porous medium has been attracted extensive attention because of its wide range of applications, e.g. ground surface water flow, seabed-wave, blood-vessel interactions and geomechanics of fractured vuggy carbonate reservoirs (Wang et al., 2004; Badia et al., 2009; Zhang et al., 2017). Let us focus on the last application. As illustrated in Figure 1, a vuggy carbonate reservoir is a typical double-porosity medium. Due to the different sizes between matrix pores and vugs, the vuggy porous media are more stress-sensitive than the homogeneous porous media (Jaeger et al., 2007; Zhang et al., 2017). The Biot’s equations are usually used to model the hydro-mechanical coupling effects of porous media (Biot, 1941, 1956; Jaeger et al., 2007; Both et al., 2017).
Murad et al. (2000, 2001) were the first to use Biot-Stokes equations to study the hydro-mechanical behavior of a double-porosity medium, in which a porous cell contains micro pores and the surrounding system of macro pores, void spaces or bulk flow paths (e.g., fractured rock or aggregated soil). And the extended BJS boundary conditions were proposed in the framework of small strains. Following their works, Showalter (2005) proposed a mixed formulation for the resolvent Biot-Stokes equations with extended BJS boundary conditions. Whereafter, some numerical studies of the Biot-Stokes coupling problem have been presented (Badia et al., 2009; Ambartsumyan et al., 2017). Although we can use numerical techniques to get accurate and detailed poroelastic and fluid behaviors on the fine scale (the scale of the periodic region), it will lead to a computationally expensive problem if one focuses on the fine-scale structure on the macroscale (Kouznetsova et al., 2001; Yalchin and Hou, 2009; E, 2011). Therefore, an upscaling description of the fine-scale structural behaviors is needed when we focuses on the macroscale problem (e.g., core analysis and reservoir simulation). Usually, one can get the macroscale constitutive relations through the experiments. This approach, however, is less effective when dealing with complex coupling multi-physical process. An alternative is provided by volume averaging or homogenization techniques (Whitaker, 1999; Auriault et al., 2009), which allow the inclusion of the fine-scale descriptions within the larger scale problem. Murad et al. (2000, 2001) applied the homogenization technique to upscale the Biot-Stokes equations aiming to get the constitutive laws for the double-porosity media. However, they assumed the local Biot number Bi →∞, and then they neglected the interfacial resistance. As a result, the continuity condition of the tangential velocity on the interface is enforced instead of the BJS condition. Lewandowska and Auriault (2013) developed a macroscopic model of hydro-mechanical coupling for the case of a porous medium containing isolated cracks or vugs. The Biot-Stokes equations are applied in their work, in which the continuity conditions of the stress, fluid pressure and velocity are imposed at the porous-fluid interface. More recently, Wan and Eghbalian (2016a, 2016b) studied the hydro-mechanical description of fractured porous media based on micro poromechanics by using the mean-field theory and Mori-Tanaka homogenization scheme. They also used Biot-Stokes equations on the fine scale, but the continuity of the velocity and stress are directly applied instead of the BJS condition. Therefore, a general homogenized result of a Biot-Stokes systems coupling with BJS condition is still an opening problem. In this work, we start from the Biot-Stokes equations coupling with the BJS condition on fine scale to develop a general macroscopic hydro-mechanical model for vuggy porous media via the asymptotic expansions homogenization theory (Sanchez-Palencia, 1980; Auriault et al., 2009; Argilaga et al., 2016). To our best knowledge, the derived macroscopic Biot’s equations are novel and general for deformable vuggy porous media. It is shown that the general structure of Biot’s equations is the same as in the case of homogeneous medium but the effective poroelastic coefficients are modified and have different physical meanings. In Section 2, the Biot-Stokes equations of the problem at the fine scale is presented. In Section 3, the homogenization upscaling process is described in detail. Section 4 gives several numerical examples and discussions.
2. Biot-Stokes equations on the fine scale
The development in this section is similar to that in the papers (Chen et al., 2004; Arbogast and Lehr, 2006). Let Ω be a bounded deformable fractured or vuggy porous medium in R d , d =2, 3, as depicted in Figure 2. And domain Ω is a union of non-overlapping regions Ω m and Ω f , which are poroelastic matrix and free fluid region respectively. Let Γ fm = Ω m ∩ Ω f . The Biot’s equations in Ω m and the Stokes equations in Ω f , and the BJS interface conditions at the porous-fluid interface Γ fm are stated as follows. matrix vug Ω m Ω f Γ fm Figure 2: Schematic of a periodic porous medium with micro vugs. (1) Biot’s equations for poroelastic matrix m m mm p s s s mm p s p f mm s p m p pp σ gσ α σ σ a e ukv v u gv α e u (1) where, σ m and σ s are the total stress tensor of matrix and stress tensor of solid matrix skeleton, respectively. The density of saturated matrix ρ m = ρ f +(1- ) ρ s , here is the porosity of matrix, ρ f and ρ s are the density of fluid and solid matrix skeleton. p p is the matrix pore pressure; α is the Biot coefficient, which is a second rank tensor; a is the fourth rank elastic tensor. u s is the displacement of solid matrix skeleton, and s u denotes the differentiation of u s with respect to time. v m is the Darcy’s velocity, and v p is the average fluid velocity in matrix; k is the permeability tensor of matrix; µ is the fluid viscosity; g is the gravitation acceleration vector; γ is the storage coefficient. The strain tensor e ( u s ) is defined by Ts s s e u u u (2) (2) Stokes equations for free fluid f f ff f f ff f p σ gσ I e vv (3) where, σ f is fluid stress tensor; p f is the fluid pressure. v f is free fluid velocity. The strain tensor e ( v f ) is defined by Tf f f e v v v (4) (3) Interface conditions f d s fmf m fmf f p fmf f s fm onon2 on2 on p p v n v u nσ n σ nn e v nn e v τ v u ττ k τ (5) where, n is the outward unit normal vector directed from Ω f to Ω m ; is the unit tangential vector along the interface. The first equation is the mass balance, which requires that the normal fluid flux must be continuous across the interface. We note that this continuity of flux constrains the normal velocity of the solid matrix skeleton. The second equation represent the conservation of momentum, which requires that the stress of the porous matrix is balanced by the stress of the fluid. The third one is the balance of the normal components of the stress in the fluid phase across the interface. The fourth equation is the BJS condition modeling slip with friction, and β is the BJS slip coefficient determined by experiments. L ε=l/L l
Figure 3: Schematic of a periodic vuggy porous medium
3. Upscaling based on homogenization theory
Now, we consider the homogenization procedure for the above Biot-Stokes system. The important characteristics of this procedure is existence of two different length scales: the fine scale l , which characterizes the typical REV (Representative Elementary Volume) size, and the macro scale L , which characterizes the global variation of external forces and boundary data. Let ε = l / L , and with ε <<1, as illustrated in Figure 3. Here, we also emphasize that the microscale l p << l , which characterizes the matrix pore size. Let the vuggy porous medium Ω has a periodic mesostructured with period Y (cf. Figure 3), where m f Y Y Y , with Y m and Y f being the matrix and free fluid parts, respectively. Define m m f f : , : x x Y x x Y (6) In this work, we only consider a formal asymptotic expansion of the solution in (1)~(5), and the boundary of Ω does not play a role in this expansion. So the above definition can be rewritten as follows m m f f : , : x x Y x x Y (7) Next, we consider the following scaled problem. (1) Biot’s equations in matrix m m mm p s s s mm p s p f mm s p m p pp σ gσ α σ σ a e ukv v u gv α e u (8) (2) Stokes equations in vugs f f f2f f f ff f p σ gσ I e vv (9) (3) Interface conditions f d s fmf m fm2 f f p fmf f s fm onon2 on2 on p p v n v u nσ n σ nn e v nτ e v n v u ττ k τ (10) where, we have scaled both the fluid viscosity and permeability k ε by ε . This is the usual scaling for deriving Darcy’s equation from Stokes equations (Arbogast and Lehr, 2006). Since as ε →0, flow paths (in our case matrix pores and vugs) will constrict, and a corresponding decrease in fluid viscosity is required to maintain the flow rates. Moreover, the permeability O l k , so it have a similar scaling of viscosity. Note that these scaling treatments are implied in the stated scaling of BJS condition. Following the custom of homogenization theory, we assume that any point is described by two coordinates: x describing the global location of the point and Y y describing the location of the point within the ε -cell εY . Consequently, x and y are related by the scaling constant ε : y ~ ε -1 x . By the chain rule, the following relation holds x y (11) where x and y denotes the gradient operators with respect to x and y , respectively. Assuming that the solution to (8)~(10) behave as if it was a function of these two coordinates and that it can be expanded in a power series in terms of ε , the stress, pressure, displacement and velocity are then can be expanded in the following asymptotic form , , , , , , , , t t t t x y x y x y x y (12) where, i denotes the different variables, e.g. stress, pressure, displacement and velocity, which are Y -periodic in y , x , Y y . Then, we shall substitute (12) into (8)~(10), and analyze the resulting equations. (1) Biot’s equations in matrix We now we shall substitute (12) into (8), apply (11), and collect terms with like powers of ε -1 and ε , which yield
0m m0 0 0 0m p s s m0 0 1s s s m0 0 0m p s m0p m00m s m0 1m m m m0 0 11 m p p f m0 1 0 1m m s s y yx yyy yx yx yx y x y
Yp YYYp YYYp p Y σσ α σ a e uσ a e u a e uv v uv α e uσ σ gk v gv v α e u e u
0p m in p Y (13) By the second and fifth equations of (13), we have , , , t p p t u u x x (14) That is, u and p are independent of y . This implies the intuition that the local average of s u and p p does not oscillate. (2) Stokes equations in vugs yx y yyx y p YYp YYY σ σ Iσ σ gσ I e vvv v (15) By the first equation of (15), we have , , , p p t t x σ σ x (16) So σ and p are also independent of y . (3) Interface conditions onononon2 on2 onon yy YYp p YYp p YYY v n v u nσ n σ nσ n σ nn e v nτ e v n v u ττ k τv n v u n (17) where Y fm = Y f ∩ Y m , and the third equation implies that , p p p t x (18) Apply the first three equations of (13), the second equation of (13), the second equation of (17), and equation (18) to see that
1s m1 0 0s s fm : 0 in: : on y yy x
Yp Y a e ua e u n a e u I α n (19) This system forms an elliptic problem in y for u that can be solved in terms of x e u and p . Therefore, by applying the linearity property, u can be represented in terms of x e u and p : , , x t p t u ξ y e u x ζ y I α x (20) where, khi ξ y is a third rank tensor; ζ y is a vector. Substitute (20) into (19), we can get the following two local base cell problems. First, for each k and h ( k , h = 1, 2, 3), we define the vector kh ξ y to be the solution of mfm : 0 in: : on khy ykhy YY a e ξa e ξ n a e n (21) where the tensor , , 1, 2,3 ij ik jh e i j e , with ik being the Kronecker symbol. In the same way, the vector ζ y satisfies mfm : 0 in: on y yy YY a e ζa e ζ n I n (22) Then, applying the second and third equations of (13) and (20), we obtain : : : y y x p σ α a e ζ I α a I e ξ e u (23) Now, we define the total stress as
00 m mT 0f f inin YY σσ σ (24) And let us define the average operator <∙> l l Y y lY y (25) Averaging the total stress σ over Y and using the first equation of (15), (16) and (23), we obtain f m0 0 0 0 0T f m eff eff s : x p σ σ σ α a e u (26) where meff m c meff :: yy α α I a e ζ I αa a I e ξ (27) Integration of the seventh equation of (13) over Y m , and applying the divergence theorem, the periodicity condition, and the fourth equation of (17), we obtain m fm m d d d 0 xY Y Y y A y σ n σ g (28) Again, applying the divergence theorem, the periodicity condition, and integrating the second equation of (15), we see that f fm f d d d 0 xY Y Y y A y σ n σ g (29) By adding the equations (28) and (29), we obtain d d 0 xY Y y y σ g (30) where, the mass density of the bulk material is defined as f fs s inin YY (31) Then, by dividing by volume or area | Y |, (30) implies that x σ g (32) In the derivation of (32), a volume averaging theorem (Whitaker, 1999) has been used together with the periodicity condition, which implies that fm pb x x xY Y
A AY Y σ σ n σ n σ σ (33) The equation (32) is the macroscopic Navier equation for a deformable vuggy porous medium.
In this section, it follows from the sixth and eighth equations of (13), the second, third, fourth equations of (15), and the first, sixth, seventh equations of (17) that y xy y y xy yy p p YYp p YYYp p YY k v gv e v gvv n v u nn e v nτ e v n v u ττ k τ (34) Applying the equations (14) and (18) leads to in0 in2 in0 inon2 on2 on y y y yy yy p p YYp p YYYp p YY k v gv e v gvv n v nn e v nτ e v n v ττ k τ (35) where * 0 0f f s v v u (36) The first and third equations of (35) imply that , ,, ,1, , ,1, , , i i ii i ii i ii i i p g p tp g p tt g p tt g p t x y y xx y y xv x y ω y xv x y ω y x (37) for each i -th dimension ( i =1,…, d ). With i e being the standard Cartesian basis vector in the i -th direction, let il y and il ω y ( l =f, m) be the periodic solution of the following Darcy-Stokes problem in the base cell in0 in2 in0 inon2 on2 on i i iiy i i iy yiyi i i i iy i iy YYYYYYY k ω eω ω eωω n ω nn e ω nτ e ω n ω ττ k τ (38) We now apply the averaging operator in the following sense m f m0 0 * 0T m f m f f i i i iY Yi y y p gY v v v ω ω (39) Then together with the fourth equation of (13) and (36), we see that m f m f0 0 * 0 0 0 0T m f p s f sm f p0 0 0 0 0p f c m s T T stotal porositytotal pore fluid velocity v v v v u v uv v u v u (40) where, p0T v is the fluid velocity through the total pores in fractured or vuggy porous media, including the matrix pores and fractures (or vugs); T is the total porosity of the fractured or vuggy porous media, including double porosity, i.e. matrix porosity and vuggs’ porosity. Now let the matrix k eff be defined by m f eff m f ij i ij jY Y k y yY ω ω (41) Then we can obtain the following macroscopic Darcy’s equation p p x x k v g kv v u g (42) Finally, let us collect the ninth equation of (13) and the fifth equation of (15), which leads to : in0 in x y x yx y p YY v v α e u e uv v (43) Following the definition of total Darcy velocity in (40) and (36), we require the following forms : in: in x y x yx y x p YY v v α e u e uv v I e u (44) where, the relation : x x u I e u has been used. Now integrating the above equations over Y , and applying the average operators, we find fm m f0 1 0 1 0 0T s s s s x x y xY y pY v n u α e u e u I e u (45) Then it follows from a further application of the divergence theorem and periodicity condition for the second term on the left hand side, and use of (23) that effeff m0 0T m c sm 0m m 0m c sm 0m : :: : :: x y xyy xy pp α v α I I α e ξ e uI α e ζ I αα I a e ζ I α e uI α e ζ I α (46) here we used the following relations (Chen et al., 2004) m m m : : y y y I e ξ ξ a e ζ (47)
With this, the macroscopic equations of a deformable vuggy porous medium are given by (26), (32), (42) and (46). And the macroscopic equations can be rewritten as follows: x x x p p p xx σ gσ α a e ukv v u gv u α e u (48) where the corresponding coefficients are defined as following: m feff eff m f meff eff meff eff m c meff m ij i ij j khijkh ijlm ylmijkh khkh kh ij ij yijkhij ij yjk ki ki ka a a eI I eI e I k ω ωa ξα ξζ (49) The formula of (48) are classical Biot’s equations. That is, the macroscopic hydro-mechanical behavior of vuggy porous media is also described by using the Biot’s theory in the framework of small3
With this, the macroscopic equations of a deformable vuggy porous medium are given by (26), (32), (42) and (46). And the macroscopic equations can be rewritten as follows: x x x p p p xx σ gσ α a e ukv v u gv u α e u (48) where the corresponding coefficients are defined as following: m feff eff m f meff eff meff eff m c meff m ij i ij j khijkh ijlm ylmijkh khkh kh ij ij yijkhij ij yjk ki ki ka a a eI I eI e I k ω ωa ξα ξζ (49) The formula of (48) are classical Biot’s equations. That is, the macroscopic hydro-mechanical behavior of vuggy porous media is also described by using the Biot’s theory in the framework of small3 strains. However, the corresponding coefficients are more complex and difficult to evaluate. Here, we summarize the properties of the above coefficients. At first, it can be shown that eff k is positive definite and symmetric (Arbogast and Lehr, 2006). We also note that eff α is symmetric. Second, it is easy to prove that the effective elasticity coefficient eff a is symmetric and positive definite (Chen et al., 2004; Auriault et al., 2009). Finally, the last equation of (49) shows that eff .
4. Numerical experiments and discussions
In this section, several numerical experiments are conducted to verify the validation of the proposed macroscopic equations and to evaluate the macroscopic hydro-mechanical properties. The identification and the use of the macroscopic constitutive equations of the vuggy porous media require the solution of the cell problems (21)-(22) for the mechanical part and cell problem (38) for the fluid flow part. However, even for simple configurations of the elementary cell, those problems cannot be solved analytically and a numerical computation is needed. In this work, all computations are performed by using the finite element method (FEM). In this work, the cell problems of the mechanical part (Eqs. 21 and 22) are solved using classical Galerkin FEM. For the cell problem of the fluid part (Eq 38), we use a mixed FEM in the porous flow and free fluid flow regions. For more details on the said numerical schemes, the reader is referred to (Bower, 2009; Arbogast and Brunson, 2007; Huang et al., 2011).
As shown in Figure 4-a, a periodic medium (10 -3 m × 10 -3 m) with a circular vug is considered, and the porous matrix is assumed to be isotropic and homogeneous (i.e. no vug). We study the impacts of Poisson’s ratio of matrix and volumetric fraction of the vug on effective Biot coefficient α eff with two cases. In Case 1, the matrix is assumed to be non-porous, which is the same as the work of (Lydzba and Shao, 2000), and the Young modulus and Biot coefficient α of matrix are 1 GPa and , respectively. In Case 2, the matrix is porous, and the parameters are given as: Young modulus E , 1 GPa; the Poisson’s ratio ν , 0.3; Porosity of matrix ϕ , 0.126; Biot coefficient α , 0.334 I ; Storage coefficient γ , 0.218 GPa -1 . Due to the symmetries of the porous medium, it yields an isotropic macroscopic model. Herein, we applied the plane strain numerical simulation. At first, we keep the volumetric fraction of vug as a constant (i.e. ϕ c =0.126), and change Poisson’s ratio of the matrix to different values. Then keeping the Poisson’s ratio as a constant ( ν =0.2), we change the volumetric fraction of vug to different values. Figure 4-a and 4-b show the solutions of the cell base problems, i.e. Eqs. (21)-(22) . The α eff obtained for different situations are plotted in Figure 5, and excellent agreements between the current work for non-porous matrix and the references (Lydzba and Shao, 2000) are achieved. Comparing to the non-porous matrix case, the Biot coefficient α eff of the porous matrix cases is increased, which means the hydro-mechanical coupling behavior becomes stronger. At the same time, this enhancement feature becomes less marked with the increase of the Poisson’s ratio and volumetric fraction of the vug. Table 1 shows the comparison of the calculated effective Biot’s coefficients between this work and the references. l r matrix vug (a) (b) (c) Figure 4: (a) Schematic of a porous medium with a circular vug, (b) the solution ξ of base problem (21) corresponding to the x -direction displacement, (c) the solution ζ of base problem (22) corresponding to the total displacement results. Figure 5: Biot coefficient α versus Poisson's ratio (top) and porosity of the vug (bottom). To assess the validity of this work for calculating effective permeability, a same porous medium as above with different vug radius is considered, and the matrix permeability is 10 mD (milli Darcy). The isotropic effective permeability for different vug radius is presented in Table 1, and results of the current B i o t c o e ff i c i e n t α Poisson's ratioLydzba and Shao (2000)Current work for non-porous matrixCurrent work for porous matrix B i o t c o e ff i c i e n t α Porosity of the vug
Lydzba and Shao (2000)Current work for non-porous matrixCurrent work for porous matrix model are compared with those of (Huang et al., 2011) and (Wan and Eghbalian, 2016b). As shown in Table 2, the effective permeability increases obviously with the increasing of vug’s radius, and our numerical results get a good match with that of the reference results. It should be noted that the results of (Wan and Eghbalian, 2016b) are little different from the current work and (Huang et al., 2011). This is because the former used the continuity condition of tangential velocity but the last two used the BJS slip conditions. The impacts of the BJS slip coefficient on the effective permeability have been investigated in our previous work (Huang et al., 2011). Table 1: The comparison of the effective Biot’s coefficients between this work and the references
Biot’s Coefficients Lewandowska and Auriault (2013) Wan and Eghbalian, (2016b) This work α eff γ eff (GPa -1 ) 0.287 - 0.28697 Table 2: The effective permeability of base cell with different vug radius r / L k = k (mD) Wan and Eghbalian (2016b) Huang et al. (2011) This work 0.02 1.00063E+01 1.000628E+01 1.000628E+01 0.1 1.01571E+01 1.015803E+01 1.015810E+01 0.2 1.06283E+01 1.064756E+01 1.064840E+01 0.4 1.25133E+01 1.287336E+01 1.287413E+01 0.6 ― 1.790391E+01 1.790512E+01 0.8 ― 3.104431E+01 3.104633E+01 0.98 ― 1.374829E+02 1.374740E+02 In this subsection, a porous medium with a penny-shaped vug is used to verify the current work. The semi-axes of the vug are 0.2 mm in the x -direction and 0.002 mm in the y -direction, and the parameters are the same as Case 1 in section 4.1. The effective Biot’s coefficients are reported by (Lewandowska and Auriault, 2013; Wan and Eghbalian, 2016b) are compared to this work, which is shown in Table 3. It can be seen that the results of current work are in excellent agreement with those of (Lewandowska and Auriault, 2013), while there is a slight error between those of (Wan and Eghbalian, 2016b). It should be note that the macroscopic model proposed by Wan and Eghbalian is aimed to fractured porous media. When the shape of vugs becomes closer, the closed-form of the poroelastic parameters will agree better with our work. Two calculations were performed to solve problems (21) by applying a macroscopic unit strain along the x -direction and y -direction, respectively. Displacements ξ and ξ are presented in Figure 6. It should be noted that there are small differences between our current work and Lewandowska and Auriault (2013). This is because they used the simplified boundary conditions for the base problem (21) as described in section 1. Although the Biot coefficient of the porous matrix is isotropic, the effective Biot coefficient of the porous medium is anisotropic due to the penny-shaped vug. l matrix r x r y vug Figure 6: (a) Schematic of a porous medium with a penny-shaped vug, (b) the solution ξ of base problem (21) corresponding to the total displacement, (c) the solution ξ of base problem (21) corresponding to the total displacement results. Table 3: The comparison of the effective Biot’s coefficients between this work and the references
Biot’s Coefficients Lewandowska and Auriault (2013) Wan and Eghbalian, (2016b) This work 𝛼 (cid:2915)(cid:2916)(cid:2916)(cid:2869)(cid:2869) ( ξ ) 0.404 0.394 0.4043 𝛼 (cid:2915)(cid:2916)(cid:2916)(cid:2870)(cid:2870) ( ξ ) γ eff (GPa -1 ) 0.298 - 0.2981 The effective elasticity tensor of the vuggy porous medium is compared to that of the porous matrix, as listed below. All the components of effective elastic tensor are decreased due to the presence of a penny-shaped vug, especially in the y -direction (i.e. along the semi-minor axis). a a , e a a In this subsection, the evaluations of the geometric dependences of effective poroelastic parameters are investigated (as shown in Figure 7). The sizes and parameters of the periodic media are the same as the Case 2 in section 4.1, but with Poisson’s ratio 0.25. In addition, the volume fraction of vugs is 0.126 in all examples. As shown in Table 4, the numerical results indicate that the effective parameters are correlated with the vug’s shape. l r matrix vug l matrix r vug l matrix r vug l matrix r vug l matrix r x r y vug Figure 7: Geometry of porous media with different shape vugs (left), the solution ξ of base problem (21) corresponding to the x -direction displacement (middle), and the solution il ω y of base problem (38) corresponding to the x -direction velocity profile (right). Table 4: The comparison of the effective coefficients between different shape vugs Effective properties matrix circular square cross-shape rhombus ellipse γ eff (1/GPa) 0.218 0.28792 0.31138 0.33744 0.33811 0.38030 a (GPa) 1.2 0.82625 0.79061 0.70950 0.56220 0.92757 a (GPa) 1.2 0.82625 0.79061 0.70950 0.86512 0.35898 a (GPa) 0.4 0.26379 0.22629 0.23379 0.25592 0.18956 a (GPa) 0.4 0.23437 0.18863 0.17935 0.15603 0.12851 k (10 -14 m ) For effective Darcy permeability k eff , its directional value is dependent on the length of a vug along with the corresponding direction. And it increases with the increasing length. This implies that the vug connectivity may be the most important factor for effective permeability. For effective elastic Young’s modulus, the second equation of (49) indicates that its value is dependent on the averaging stress under the traction of kh ξ imposed at porous-fluid interface. Therefore, its value is strongly correlated to the interface distribution and position on the medium. The numerical results of a ijkl in Table 4 confirm that. For the effective Biot’s coefficient, the third equation of (49) indicates that its value depends on the porosity of vugs and the traction of kh ξ imposed at porous-fluid interface. In this example, the porosity of vugs is same. As a result, the effective Biot’s coefficient also relies on the shape and distribution of the porous-fluid interface. The parameter γ eff has the similar variation trend to the effective permeability. Actually, t he shape of vug determines the interface between porous matrix and free fluid flow region, which is the key for the analysis. Vug are defined as visible pores that are significantly larger than adjacent grains or crystals (Lucia, 2007). They may be formed through a variety of processes. Most commonly, vugs may form when mineral crystals or fossils inside a rock matrix are later removed through erosion or dissolution processes, leaving behind irregular voids. In this subsection, we will give a simplified 3D model for this typical vug, as shown in Figure 8. There is a tortuous pipe, and the medium part was enlarged due to the erosion process. The volume fraction of the vug is 0.0791. Figure 9 shows the meshes for different cell problems, which are corresponding to solid part and fluid flow part. The numerical results of different components have been shown in Figure 10. The calculated effective elastic parameters and the effective Darcy permeability are listed in Table 5. As illustrated in Figure 8, the tortuous vug actually connect with each other along the x -direction due to the periodicity. As a result, the component of the effective permeability in x -direction is larger than the components in y and z -directions will decrease more. Figure 8: 3D Geometry of a porous medium with a tortuous vug Figure 9: 3D Geometry of a porous medium with a tortuous vug (a) the solution ξ of cell problem (21) (b) the solution ζ of cell problem (22) (c) the solution x il e of cell problem (38) (d) the solution x il e ω of cell problem (38) Figure 10: The numerical results of several components of the solutions of cell problems Table 5: The calculated effective parameters based on the solutions of the cell problems
Effective properties Porous matrix The 3D vuggy porous medium γ eff (1/GPa) 0.218 0.27237 k (m ) -14 -05 k (m ) -14 -14 k (m ) -14 ×10 -14 The calculated effective elasticity tensor is compared to that of the porous matrix, as listed below. All the components of effective elastic tensor are decreased due to the presence of the tortuous vug, especially along the y and z -directions. a a e a a
5. Conclusions
We have systematically derived the macroscopic equations for the mechanical behavior of a deformable vuggy porous medium saturated with a Newtonian fluid via the homogenization theory. The developed formulation presented in this paper is novel since it takes into account the general Beavers-Joseph-Saffman (BJS) interface boundary conditions for Biot-Stokes system. In the case of Biot-Stokes coupling hydro-mechanical system, these macroscopic equations coincide with classical Biot’s equations. And the homogenization approach determines the form of the macroscopic constitutive relationships between variables, and shows how to compute the poroelastic coefficients in these relationships. Several numerical tests have been conducted to demonstrated the calculation procedures. The numerical results indicates that the existence of the vugs have significant impacts on the effective elastic and hydro parameters of the vuggy porous media. In addition, we note that the calculations of the macroscopic properties only depend on the decoupling three cell problems, i.e., two Navier equations for elastic problem and one Darcy-Stokes flow problem with BJS condition. It also show that our homogenized results provide a natural way of modeling realistic deformable vuggy porous media. In this work, only small strain and elastic problem is considered. In the future, the finite deformation analysis will be conducted, and the corresponding elastoplastic problem is also the next topic.
Acknowledgement
This work was supported by National Nature Science Foundation of China (Grant No. 51404292), the Fundamental Research Funds for the Central Universities (18CX05029A, 17CX06007), National Science and Technology Major Project (2017ZX05009001, 2016ZX05060-010).
References
Ambartsumyan, I., Khattatov, E., Yotov, I., Zunino, P., 2017. A Lagrange multiplier method for a Stokes-Biot fluid-poroelastic structure interaction model. Arbogast, T., Brunson, D.S., 2007. A computational method for approximating a Darcy–Stokes system governing a vuggy porous medium. Comput. Geosci. 11, 207–218. Arbogast, T., Lehr, H.L., 2006. Homogenization of a Darcy–Stokes system modeling vuggy porous media. Comput. Geosci. 10, 291–302. Argilaga, A., Papachristos, E., Caillerie, D., Pont, S.D., 2016. Homogenization of a cracked saturated porous medium: Theoretical aspects and numerical implementation. Int. J. Solids Struct. s 94–95, 222–237. Auriault, J., Boutin, C., Geindreau, C., 2009. Homogenization of Coupled Phenomena in Hetrogenous Media. Baber, K., Mosthaf, K., Flemisch, B., Helmig, R., Müthing, S., Wohlmuth, B., 2012. Numerical scheme for coupling two-phase compositional porous-media flow and one-phase compositional free flow. Ima J. Appl. Math. 77, 887–909. Badia, S., Quaini, A., Quarteroni, A., 2009. Coupling Biot and Navier–Stokes equations for modelling fluid–poroelastic media interaction. J. Comput. Phys. 228, 7986–8014. Beavers, G.S., Joseph, D.D., 1967. Boundary conditions at a naturally permeable wall. J. Fluid Mech. 30, 197–207. Biot, M., 1956. General solutions of the equation of elasticity and consolidation for a porous material. Biot, M.A., 1941. General Theory of Three ‐ Dimensional Consolidation. J. Appl. Phys. 12, 155–164. Both, J.W., Borregales, M., Nordbotten, J.M., Kumar, K., Radu, F.A., 2017. Robust fixed stress splitting for Biot’s equations in heterogeneous media. Appl. Math. Lett. 68, 101–108. Bower, A.F., 2009. Applied Mechanics of Solids. CRC Press,. Carraro, T., Goll, C., Marciniak-Czochra, A., Mikelić, A., 2015. Effective interface conditions for the forced infiltration of a viscous fluid into a porous medium using homogenization. Comput. Methods Appl. Mech. Eng. 292, 195–220. Chen, J., Sun, S., Wang, X.P., 2014. A numerical method for a model of two-phase flow in a coupled free flow and porous media system. J. Comput. Phys. 268, 1–16. Chen, Z., Lyons, S.L., Qin, G., 2004. THE MECHANICAL BEHAVIOR OF A POROELASTIC MEDIUM SATURATED WITH A NEWTONIAN FLUID. Int. J. Numer. Anal. Model. 1, 75–97. E, W., 2011. Principles of Multiscale Modeling. Cambridge University Press. Goyeau, B., Lhuillier, D., Gobin, D., Velarde, M.G., 2003. Momentum transport at fluid-porous interface. Int. J. Heat Mass Transf. 46, 4071–4081. Huang, Z.Q., Yao, J., Li, Y., Wang, C., Lv, X., 2011. Numerical Calculation of Equivalent Permeability Tensor for Fractured Vuggy Porous Media Based on Homogenization Theory. Commun. Comput. Phys. 9, 180–204. Huang, Z.Q., Yao, J., Li, Y.J., Wang, C.C., Lü, X.R., 2010. Permeability analysis of fractured vuggy porous media based on homogenization theory. Sci. China Technol. Sci. 53, 839–847. Jaeger, J., Cook, N.G., Zimmerman, R., 2007. Fundamentals of Rock Mechanics, Fourth Revised Edition, 4th ed. Karasuyama, H., Mukai, K., Tsujimura, Y., Obata, K., 2011. On the Importance of the Stokes-Brinkman Equations for Computing Effective Permeability in Karst Reservoirs. Commun. Comput. Phys. 10, 1315–1332. Kouznetsova, V., Brekelmans, W.A.M., Baaijens, F.P.T., 2001. An approach to micro-macro modeling of heterogeneous materials. Comput. Mech. 27, 37–48. Lewandowska, J., Auriault, J.L., 2013. Extension of Biot theory to the problem of saturated microporous elastic media with isolated cracks or/and vugs. Int. J. Numer. Anal. Methods Geomech. 37, 2611–2628. Lucia, F.J., 2007. Carbonate Reservoir Characterization: An Integrated Approach. Springer. Lydzba, D., Shao, J.F., 2000. Study of poroelasticity material coefficients as response of microstructure. Mech. Cohesive ‐ frictional Mater. 52, 149–171. Mikelic, A., Jäger, W., 2000. On The Interface Boundary Condition of Beavers, Joseph, and Saffman. Siam J. Appl. Math. 60, 1111–1127. Mosthaf, K., Baber, K., Flemisch, B., Helmig, R., Leijnse, A., Rybak, I., Wohlmuth, B., 2011. A coupling concept for two ‐ phase compositional porous ‐ medium and single ‐ phase compositional free3