Numerical modelling of two-phase flow in fractured rock masses using zero-thickness interface elements

XVI International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS 2021 E. O˜nate, D. Peric, M. Chiumenti and Eduardo de Souza Neto (Eds) NUMERICAL MODELLING OF TWO-PHASE FLOW IN FRACTURED ROCK MASSES USING ZERO-THICKNESS INTERFACE ELEMENTS L. BARANDIAR´ AN∗, J. LIAUDAT†, C. M. L´ OPEZ∗AND I. CAROL∗ ∗ETSECCPB (School of Civil Engineering-Barcelona) Universidad Polit´ecnica de Catalu˜na Campus Norte UPC, 08034 Barcelona, Spain e-mail: [email protected], [email protected], [email protected] †Faculty of Civil Engineering and Geosciences Delft University of Technology Building 23, Stevinweg 1, 2628 CN Delft email: [email protected] Key words: two-phase flow, thermo-hydro-mechanical model, fracture, discontinuity, zero-thickness interface elements, Finite Element Method Abstract. In recent years, the authors and co-workers have developed a 3D finite element model for coupled thermo-hydro-mechanical (THM) problems in fractured rock masses. Zero-thickness interface elements are used for taking into account explicitly the effect of fractures and discontinuities in the fluid flow as well as the effect of fluid pressure in the crack propagation. Furthermore, the use of zero-thickness elements as a discrete modelling approach for fractures and discontinuities makes it possible to account for the heat transport taking place within these elements, even when advection dominates over diffusion (high P´eclet number). The model has been implemented in the finite element code DRAC5, which is equipped with fracture-based interface elements and MPI parallel capabilities. The code was originally developed considering water-saturated porous medium and fractures. The new developments described in the present paper, include the extension of the original formulation to the case of two-phase (liquid and gas) flow within the porous medium and discontinuities. The liquid includes only liquid water species, while the gas phase includes water vapour and gas species. The formulation includes the equilibrium equation, the mass balance of water and gas species and the energy balance equation. The parameters of the retention and relative permeability curves for the interface elements, such as the gas entry value and the residual water saturation, are updated with the variation of the normal aperture. The new capabilities of the model are illustrated with some academic verification examples. 1 INTRODUCTION Within the context of geomechanics, the study of multiphase flow in fractured rock masses considers the simultaneous flow of at least two fluid phases i.e. water, gas and, in some applications, also oil, as well as their interactions with the rock mass. Generally 1 L. Barandiar´an, J. Liaudat, C. M. L´opez and I. Carol speaking, fluid phases are composed by a certain number species that may be transferred from one phase to another due to changes in the multiphase system conditions. Species may be defined as the different types of chemical molecules present in the system whereas phases may be defined as a mixture of the said species isolated by means of a fluid interface [1]. Two-phase flow, a particular case of multiphase flow, takes place in different types of coupled problems which are relevant to the design and performance evaluation of engineering projects such as nuclear waste repositories [2], coalbed methane recovery [3], CO2storage in depleted oil reservoirs [4] and fractured shale gas reservoirs [5]. This paper briefly presents a FE model in which fractures are explicitly represented using zero-thickness interface elements. For this model, an existing THM formulation for a single liquid fluid phase [6] is extended to consider the simultaneous two phase (liquid and gas). The main assumptions of the model are the following: the system is formed by one solid phase and two fluid phases (the water phase and the gas phase), the gas phase is composed by the gas species and eventually also by water species, the water phase is composed only by water species, i.e. gas solubility in the aqueous phase is neglected. The mass transfer between phases is treated with the compositional approach. The mechanical problem is solved assuming small strains. The field variables are: displacements (u), phase pressures (pwand pg) and temperature (T). Ideal gas law is considered and Kelvin-Laplace as well as Clasius-Clapeyron laws are used to calculate the water vapour pressure. van Genuchten [7] characteristic curve is used to compute the water saturation and relative permeabilities (hysteresis is not considered). 2 POROUS MEDIUM FORMULATION 2.1 MECHANICAL PROBLEM The following equilibrium equation of the porous medium is considered: LuTσ0−LuTαm(Swpw+Sgpg) + ρg=0(1) where LuTis the differential operator, σ0is the effective stress vector (Voigt notation), α is the Biot’s coefficient, mis the equivalent in the Voigt notation to the identity tensor, ρis the average density of the porous medium, considering the contributions of the solid and fluid phases, gis the gravity acceleration vector, and Swand Sgare the volumetric saturation degrees of the pore space with respect to the water and gas phases, respectively. The solid phase is deformable and subject to thermal strains. 2.2 FLOW PROBLEM Within the compositional modelling approach, the mass balance equations for the flow problem are not developed in terms of the mass of each phase, but in terms of the mass of each chemical species in each phase. The general mass balance equation of a ψ-species 2 L. Barandiar´an, J. Liaudat, C. M. L´opez and I. Carol within a π-phase is defined as follows: ∂nSπρπψ ∂t +∇T·nSπρπψuπψ+∇T·nSπρπψvπ= ˙mπψ (2) for π= w,g and ψ=w,g where nis the porosity, Sπis the π-phase saturation, ρπψ is the density of the ψ-species within the π-phase, vπis the average velocity of the π-phase in the pore space, uπψ is the diffusive-dispersive velocity of a ψ-species within the π-phase and ˙mπψ is a mass sink/source term accounting for the exchange of the species with other phases. The diffusive-dispersive velocity is assumed to be proportional to the concentration gradient of ψspecies in the π-phase (Fick’s law). The velocity of the π-phase is assumed to be proportional to the gradient of pπ(Darcy’s law). Since dissolution of gas species in the liquid phase is neglected, uwg=uww= 0. The absolute value of the sink/source term in the right hand side of the Equation (2) is the same for the water species in both phases because the mass of evaporated water coming from the liquid phase is equivalent to the mass entering to the gas phase ( ˙mwg=−˙mww). The mass balance equation of a ψ-species in the system is obtained as the sum of the mass balance equations of that ψ-species in each phase, as established by the compositional formulation [8]. 2.3 THERMAL PROBLEM The following assumptions are considered for establishing the energy balance equation: the viscous dissipation related to the fluid flow is neglected, the energy exchange due to the compression-decompression of the fluid phases is neglected and the system is in thermodynamic equilibrium (at a given point the temperature of the solid, liquid and gas phases is the same). Also, the final system equation is obtained as the sum of the thermal balance of each ψ-species within a π-phase. Under this assumptions the following equation is obtained: (ρcp)eff ∂T ∂t +ρwcw pqws +ρgcg pqgs·∇T−∇T·(χeff∇T) = −˙m∆hgw(3) where (ρcp)eff and χeff are the effective heat capacity and the effective thermal conductivity of the porous medium, qws and qgs are the volumetric (Darcy’s) flows of each phase, ∆hgwis the specific latent heat of water evaporation, and ˙mis the mass rate of water evaporation. The first term in the left hand side of Equation (3) is the storage term, the second term represents the heat transport (advection) and the last corresponds to the heat conduction. Finally, the term in the right hand side represents the heat sink/source due to phase change. Heat conduction is calculated by means of a generalized Fourier’s law and heat transport is considered only for a low-advection condition (Pe<1). As a first approach, the time derivatives and the gradients of gas pressure have been neglected given that gas pressure changes are usually very slow and their gradients small in the porous medium [9]. Total volume variation is expected to be negligible, therefore its contribution to the energy balance has not been considered. 3 L. Barandiar´an, J. Liaudat, C. M. L´opez and I. Carol 3 DISCONTINUITY FORMULATION The physical model for discontinuities is based on similar assumptions as for the porous medium. The main differences are the influence of the discontinuity aperture in both the constitutive laws and balance equations, and the absence of the solid phase, given that the discontinuities are modelled as the space between two parallel surfaces where the two-phase flow takes place (Figure 1). Figure 1: Schematic representation of a discontinuity in a continuum porous medium where (n,l) represent the local orthogonal coordinate system, wis the discontinuity width or aperture, MP is the mid-plane of the discontinuity, Qlrepresents the longitudinal local flow and qtop and qbot represent the transversal flow incoming from the porous medium located at the top and bottom of the discontinuity. After Segura [10]. 3.1 MECHANICAL PROBLEM The equilibrium equation of a discontinuity is defined as follows: σmp =σ0 mp −mmpSwpw mp +Sgpg mp=0(4) where σmp and σ0 mp are the total and effective traction vectors at the mid-plane, pw mp and pg mp are the average pressure of water and gas phases in the mid-plane, and mmp = 1 0 0T. The stress traction vector is defined as follows: σmp =σnτl1τl2T(5) where σnis the normal stress and τl1and τl2are the tangential components. The relative displacement at a mid-plane point of the discontinuity is represented by: r=rnrl1rl2T(6) where rnis the normal component (discontinuity aperture) and rl1and rl2are the tangential components. Figure 2 shows the particular case of a one-dimensional discontinuity in a 2D continuum. Given that no solid phase is considered within the discontinuity and that the discontinuity space is assumed to be completely occupied by the two fluid phases, the saturation degrees are defined as the volume fraction of water or gas phase with respect to the total volume of both phases. 4 L. Barandiar´an, J. Liaudat, C. M. L´opez and I. Carol Figure 2: Stresses tractions and relative displacements at the mid-plane of a 1D discontinuity in a 2D case. After Segura [10]. 3.2 FLOW PROBLEM The general mass balance equation of a ψ-species within a π-phase in the discontinuity is defined as follows: ∂rnSπρπψ ∂t +∇T j·rnSπρπψuπψ l+∇T j·rnSπρπψvπ l= ˙mπψ mp (7) for π= w,g and ψ=w,g where rnis the aperture of the discontinuity, ∇j=∂ ∂l1 ∂ ∂l2Tis the partial differential operator for the local in-plane axes, uπψ lis the longitudinal diffusive-dispersive velocity of the water vapour and gas species (only for the gas phase) and vπ lis the longitudinal average velocity of the π-phase. The longitudinal advective flow is calculated by means of a generalized Darcy’s law, considering that the longitudinal transmissivity grows with the cube of the aperture. The saturation degree of the discontinuity is obtained as a function of the normal aperture and the capillary pressure pc=pg−pwat the mid-plane, by means of a modified version of the classic van Genuchten [7] expression for porous media: Sw= (1 −Swr(rn))   1 1 + pc pb(rn)n  m +Swr(rn) (8) where Swr is the residual water saturation, pbis the air-entry pressure value and m,nand λare fitting parameters. In contrast to the original expression, pbis not constant but a monotonic decreasing function of the normal aperture. The effect of the saturation degree on the effective transmissivity of the discontinuity is introduced by means of relative permeability coefficients obtained as linear functions of the the respective saturation degrees. Besides the longitudinal transmissivity, the existence of a discontinuity may also represent an obstacle or resistance to the fluid flow in the transversal direction, for instance due to the transition from a pore system into an open channel and back into a pore system. 5 L. Barandiar´an, J. Liaudat, C. M. L´opez and I. Carol This resistance may reduce the fluid flow in the transversal direction and results in a localized fluid pressure drop across the discontinuity, according to the following equation: Qπψ n=ρπψqπ n=ρπψkπ nˇpπ mp (9) where Qπψ nis the transversal mass flow of a ψ-species within a π-phase and qπ nis the transversal volumetric flow of a π-phase, kπ nthe transversal hydraulic conductivity and ˇpπ mp the transversal pressure drop, which is defined as the difference between the pressure at the top and the bottom of the discontinuity: ˇpπ mp =pπ top−pπ bot. 3.3 THERMAL PROBLEM The general balance of energy along the discontinuity is defined as follows: (ρcp)eff ∂Tmp ∂t +ρwcw pqw l+ρgcg pqg l·∇jTmp −∇T j·(χeff∇jTmp) −rnSg ∂pg mp ∂t −qg l·∇jpg mp +rnpc∂Sw ∂t −ps∂rn ∂t =−˙mmp∆hgw(10) This equation is obtained by adopting similar assumptions as for the porous medium with regard to the viscous flow dissipation, the incompressibility of the water phase, the thermodynamic equilibrium of the system, the substitution of energy by enthalpy and the summation of the energy equation for each species in the fluid phases. The first three terms are analogous to those of Equation (3) and represent the thermal storage term, the advective heat transport and the heat conduction. The subsequent two terms are related to the time derivative and gradients of the gas pressure, which are expected to be important in the discontinuity. It is noteworthy that the negative sign in both terms represents the cooling effect associated to gas expansion. Finally, the remaining terms on the left hand side of the equation are only considered for the discontinuity (and not for the continuum) because the variation of the absolute volume occupied by the fluids may be significant, given that the discontinuity opening is unrestricted. Analogously to the transversal fluid flow, the transversal heat flow is defined as follows: hπψ n=kπψ nˇ Tmp (11) where hπψ nis the transversal heat flow of a ψ-species within a π-phase, kπψ nthe transversal thermal conductivity and ˇ Tmp is the temperature drop, which is defined as the difference between the temperature at the top and the bottom of the discontinuity: ˇ Tmp =Ttop− Tbot. 4 NUMERICAL IMPLEMENTATION The equivalent integral forms of the governing differential equations are obtained by the application of the Galerkin Weighted Residuals method and the Principle of Virtual Work. 3D finite elements are used for the spatial discretization of the porous medium. In the case of the discontinuities, zero-thickness interface elements are used for the spatial 6 L. Barandiar´an, J. Liaudat, C. M. L´opez and I. Carol discretization. The displacements, fluid pressures and temperature at the integration points of the interface elements are approximated by the interpolation of the element node values using mid-plane shape functions matrices, as well as transport matrices, which transform the nodal values into mid-plane values. The generalized trapezoidal rule is considered for the time discretization. Due to the non-linear nature of the system of equations obtained, an iterative procedure based on the Newton-Raphson method is used. The formulation has been implemented in the finite element code DRAC5 following a fully coupled approach. 5 ACADEMIC VERIFICATION EXAMPLES 5.1 CONTINUUM MEDIUM VERIFICATION EXAMPLE The problem consists of a single square continuum element of 1 m side. The boundary conditions are: zero prescribed displacements at the bottom left side of the element, zero prescribed vertical displacements at the bottom right side of the element, no in or out flow of water or gas, and a prescribed temperature increment of 100 K. The initial water and gas pressures in the element are pw= 0.9×106Pa and pg= 1.0×106Pa, respectively, the initial saturation is Sw= 49.5% and the initial temperature is 300 K. The initial porosity is n0= 0.30. The Biot’s coefficient is α= 1.0. The bulk modulus of the solid grains is Ks= 1.0×1012 Pa, and their volumetric coefficient of thermal expansion is βs= 1.0×10−4K−1. The continuum medium is linear elastic with a Young modulus E= 1.0×108Pa and a Poisson coefficient ν= 0.0. Plain strain condition is assumed. The molar mass of the gas species (air) is Mg= 0.016 kg/mol. The water phase density is 997 kg/m3. The water compressibility modulus is Cw= 0 Pa, and its volumetric coefficient of thermal expansion is βw= 2.57 ×10−4K−1. The saturation pressure is pgws 0= 3169.9 Pa and the specific latent heat is ∆hgw= 2441816.06 J/kg. The two-phase characteristic curves are obtained with the closed-form equations proposed by van Genuchten [7] using Swr = 0.20, m=0.3, n=1.43, pb=1.0×104Pa and a=0.5. The left plot in Figure 3 shows the evolution of the water phase pressure (pw), gas phase pressure (pg) and the water vapour pressure (pgw) as temperature is increased from 300 to 400 K. The heating of the porous medium results in the rise of fluid pressures, as normally expected. The porosity remains practically constant because the reduction of porosity due to the thermal expansion of the solid grains is compensated by the increase of porosity due to the volumetric expansion of the porous medium. The final saturation degree is Sw= 48.8%, which is lower than the initial saturation due to water evaporation. In order to remark the contribution of the water vapour pressure to the fluid pressures, additional pressure curves are plottled in the left diagram of Figure 3 with dashed lines corresponding to an additional simulation without taking into account water vapour, i.e. pgws 0= 0.0 Pa. The final saturation degree is Sw= 49.7% which is practically the same as the initial degree (49.5%). The slight difference may be related to gas compressibility. The right plot in Figure 3 shows the evolution of the volumetric strains of the porous medium. As expected, water vapour pressure contributes to a higher volumetric strain. 7 L. Barandiar´an, J. Liaudat, C. M. L´opez and I. Carol 300 320 340 360 380 400 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 106 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Pgv Pwv Pg Pw Pv Porosity 300 320 340 360 380 400 0.0 0.4 0.8 1.2 1.6 2.0 with vapour without vapour Figure 3: Results of the continuum medium verification example. In the left, evolution of gas phase pressure (orange), water phase pressure (blue) and vapour pressure (green). Dashed lines correspond to the evolution of gas and water pressures for a simulation without considering water vapour. In the right, the evolution of volumetric strain. 5.2 DISCONTINUITY VERIFICATION EXAMPLE 1 This example is analogous to the continuum medium verification example. A single horizontal interface element of 1 m of length is considered. The boundary conditions are: zero prescribed displacements of zero at the bottom side of the discontinuity, no in or out flow of water and gas, and a prescribed temperature increment of 100 K. The discontinuity is closed at the beginning of the simulation, therefore the initial retention curve is considered to be the same as the porous medium. The initial water and gas pressures in the element are pw= 0.9×106Pa and pg= 1.0×106Pa, respectively; the initial water saturation is Sw= 49.5% and the initial temperature is 300 K. The discontinuity is elastic and the normal stiffness coefficient is Kn= 1.0×108Pa/m. The molar mass of the gas species is Mg= 0.016 kg/mol and the specific heat capacity is cg p= 2224.72 J/(kg K). The water phase density is 997 kg/m3. The water compressibility modulus is Cw= 0 Pa, its volumetric coefficient of thermal expansion is βw= 2.57 ×10−4K−1and the specific heat capacity is cw p= 4184.5 J/(kg K). The saturation pressure is pgws 0= 3169.9 Pa and the specific latent heat is ∆hgw= 2441816.06 J/kg. The two-phase characteristic curves are obtained with the modified closed-form equations proposed by van Genuchten (8) using Swr = 0.20, m=0.3, n=1.43, pb= 1.0×104(air-entry value pressure) and a=0.5. As a first step, initial pressures, temperature and prescribed displacements are applied in the element resulting in an initial aperture of the discontinuity (rn0= 1.0×10−2m) which is used to update the residual water saturation and the air-entry value pressure. Then, the saturation degree is calculated for the open discontinuity, yielding Sw= 2.27%, which is much lower than the initial degree of the porous medium (Sw= 49.5%). The left plot in Figure 4 shows the evolution of the water phase pressure (pw), gas phase pressure (pg) and the water vapour pressure (pgw) as temperature is increased from 300 to 400 K in the discontinuity. In contrast to the porous medium example, the gas and water phase pressures reach lower final values and exhibit different behaviour: gas pressure grows at higher 8 L. Barandiar´an, J. Liaudat, C. M. L´opez and I. Carol pace than the water pressure. As with the porous medium example, additional pressure curves are plottled in the left diagram of Figure 4 with dashed lines corresponding to an additional simulation without taking into account water vapour, i.e. pgws 0= 0.0 Pa. The effect of water vapour is clearly important in the case of an almost dry discontinuity. The right plot in Figure 4 shows the variation of the discontinuity aperture (∆rn/rn0), which is equivalent to the volumetric strain in the porous medium because Poisson coefficient was considered zero (ν= 0.0). The final apertures of the discontinuity are rn= 1.31×10−2m and rn= 1.15×10−2m for the first and second (without vapour) simulation respectively. As expected, water vapour pressure contributes to a larger final aperture. It is noteworthy that the magnitude of the volumetric strains in the discontinuity are ten times higher than the volumetric strains obtained in the porous medium. 300 320 340 360 380 400 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 106 Pgv Pwv Pg Pw Pv 300 320 340 360 380 400 0 8 16 24 32 with vapour without vapour Figure 4: Results of the discontinuity verification example. In the left, evolution of gas phase pressure (orange), water phase pressure (blue) and vapour pressure (green). Dashed lines correspond to the evolution of gas and water pressures for a simulation without considering water vapour. In the right, the variation of the discontinuity aperture. 5.3 DISCONTINUITY VERIFICATION EXAMPLE 2 The problem consists of a single interface element of 1 m of length. The boundary conditions are: zero prescribed displacements, fixed water pressure, gas null flow, and a prescribed temperature increment of 40 K. The initial water and gas pressures in the element are pw= 0.9×106Pa and pg= 1.0×106Pa, respectively, the initial water saturation is Sw= 49.5% and the initial temperature is 298.15 K. The joint aperture is rn= 1.25 ×10−5 m. The molar mass of the gas species is Mg= 0.016 kg/mol and the specific heat capacity is cg p= 2224.72 J/(kg K). The water phase density is 1000 kg/m3. The water compressibility modulus is Cw= 0 Pa, its volumetric coefficient of thermal expansion is βw= 2.57 ×10−4K−1and the specific heat capacity is cw p= 4184.5 J/(kg K).The saturation pressure is pgws 0= 3169.9 Pa and the specific latent heat is ∆hgw= 2441816.06 J/kg. The two-phase characteristic curves are obtained with the modified closed-form equations proposed by van Genuchten (8) using Swr = 0.20, m=0.3, n=1.43, pb= 1.0×104 9