Abstract
Understanding the gas migration in highly water saturated sedimentary rock formations is of great importance for safety of radioactive waste repositories which may use these host rocks as barrier. Recent experiments on drainage in argillite samples have demonstrated that they cannot be represented in terms of standard two-phase flow Darcy model. It has been suggested that gas flows along highly localized dilatant pathways. Due to very small pore size and the opacity of the material, it is not possible to observe this two-phase flow directly. In order to better understand the gas transport, a numerical coupled hydraulic-mechanical model at the pore scale is proposed. The model is formulated in terms of Smoothed Particle Hydrodynamics (SPH) and is applied to simulate drainage within a sample reconstructed from the Focused Ion Beam (FIB) images of Callovo-Oxfordian claystone. A damage model is incorporated to take into account the degradation of elastic solid properties due to local conditions, which may lead to formation of new pathways and thus to modifications of fluid transport. The influence of the damage model as well as the possible importance of rigid inclusions is demonstrated and discussed.
1. Introduction
The principal objective of this study is to improve the representation and understanding of gas migration in highly water saturated clays. This objective is closely related to several projects that use a deep sedimentary rock formation, such as Callovo-Oxfordian (COx) clay [1], as a host rock and a natural barrier for radioactive waste repository [2]. These sedimentary formations are characterized by a low permeability and are water saturated in their initial state. However, during the postclosure phase, the gases, mainly hydrogen, are being produced mostly by anaerobic corrosion of iron components of the repository. This important gas production may lead to appearance of a free gas phase and to a local increase of pressure, which could influence confining properties of the repository because of desaturation of clay pores and because of the mechanical disturbance caused by a possible microfracturing in the vicinity of the installation.
In classical approaches, several gas transport mechanisms can be active, separately or simultaneously [3, 4]. With increasing rates of gas production, gas should be transported first in dissolved state by advection and diffusion [5] and then by the flow of small gas bubbles inside the water phase and by viscocapillary flow of gas inside the water phase, or by the dilatation of existing and creation of new localized pathways [3, 6, 7]. The ultimate stage appearing at very high gas production rates would be macroscopic fracturing of the host rock; however, it is rather not expected in the case of repositories, where the gas source terms are not sufficient. Experimental evidences accumulating in recent years indicate that in case of indurated clays the gas migration may directly switch from dissolution-diffusion stage to the dilatant pathways mechanism. Indeed, when an increasing gas pressure is applied to a saturated argillite sample, at a certain gas pressure (“entry pressure”), a rapid but intermittent increase of the gas outflow is usually observed [3, 6, 7]. This gas flow is not associated with a significant sample desaturation, as it would be the case in viscocapillary flows. On the other hand, the gas flow causes a macroscopic deformation of the clay sample and the associated permeability tends to increase when repeating gas pulse tests. The prevalence of such a dilatant flow can be justified by the fact that the major part of the porosity of indurated clays exhibits pore size of a few nanometers [4]. Thus, due to associated high capillary pressure in such water saturated nanosize pores, the viscocapillary and bubbly flows would be restricted to largest pores which do not form a connected macroscopic paths and are not suffecient for the macroscopic gas transfer. On the other hand, the low tensile strength of the claystone facilitates the solid matrix cracking and subsequently the gas percolation [8]. Therefore, in order to get a better understanding of dilatant transport mode, it becomes necessary to take into account the coupling of the two phase fluid/gas flow inside the pore space with solid matrix mechanics, which constitute the scope of our study together with modelling of modification of argillite mechanical properties due to damage induced by gas high pressures. The phenomena of gas dissolution and diffusion are not considered in our model, since they are out of scope of this paper. Gas slippage is also excluded from considerations, since at high gas pressures these effects are of limited importance.
The coupling between the fluid flow hydrodynamics and solid matrix deformation has been tackled recently by several authors. In [9–11], a coupled continuous hydromechanical model with embedded fractures is presented where simple planar fractures are characterized by aperture which varies depending on local conditions. This allows taking into account the dependence of medium permeability on pore fluid pressure and improving to a certain extent the predictions of the fluid flow. In [12], the author studied dilatant flows by means of a single, nonpropagating elliptical crack containing a wetting and a nonwetting fluids and embedded in a linearly elastic solid. For this simple geometry it was shown that the onset of dilatant two-phase flow depends on a dimensionless parameter, and analytic solutions were obtained for some physical properties like crack aperture, capillary pressure, and permeabilities. However, the above-mentioned approaches do not take into account neither the fracture propagation process, nor the complicated underlying pore sizes distribution and connectivity. In [13–16], a more sophisticated continuous hydromechanical model is introduced which not only takes into account variation of permeability due to local stress and fluid pressure state, but also considers variation of solid matrix mechanical properties according to an ad hoc damage model which corresponds to emergence of new microfractures. However, in these works, the details of processes taking place at the pore scale are neglected; that is, each elementary volume of the model is characterized by its permeability; thus the actual geometry of the pore space is replaced by a coarsened property together with governing equations (in particular the fluid flow is described in terms of Darcy equation).
In recent years, due to X-ray computer tomography CT-scan and Focused Ion Beam Scanning Electron Microscope (FIB-SEM) imaging, it become possible to acquire models of porous opaque materials with submicronic resolution (from 40 μm down to 5–10 nm). Thus it is possible to conduct flow simulations inside realistic pore geometries. It is to note that, for indurated argillites, the finest possible linear resolution (5–10 nm) is not yet sufficient to capture the connectivity of pores, and that most often the percolating components of pore space are formed by desiccation cracks, which can be assimilated to rough fractures. In [17], the two-phase flow was studied by Lattice Boltzmann Method inside a rough fracture extracted from a X-ray microtomography of Opalinus clay. It has been shown that, under realistic physical conditions and parameters (surface tension, pressure gradient), the nonwetting phase permeability is zero whenever this phase saturation is lower than about , while the nonwetting phase splits into bubbles which become rapidly stuck on fracture roughness. One of the limitations of this work was the lack of deformable character of the pore space, which hindered the possibility of transferring the conclusions to real argillites fissures.
In order to overcome this difficulty, we propose to implement a model taking into account both two-phase flow inside pore space and its interaction with deformable rock matrix. We aim as well at modelling of the appearance of new pore connections, and for this reason we have given up the Lattice Boltzmann Model, which is solved on underlying regular mesh and thus is not well suited for capturing of fracture initiation. It is important to stress here that these “new” connections are seen not only as true microcracks, but also as a way to take into account the underlying pores connectivity that is not visible on the CT-scans. This will allow us to work with samples that are apparently nonpercolating and to ultimately restore visible connectivity during gas migration.
The model presented in this paper is based on the framework of the Smoothed Particle Hydrodynamics (SPH) which is a Lagrangian mesh-free method [18, 19]. Recently it was successfully applied to simulations of miscible, immiscible [20], and unstable [21] fluid flows as well as to drainage [22] at pore scale. In [23], the SPH model was applied to direct simulations of the fluid flow through the samples obtained from the sandstone micro-CT images. A good agreement between the SPH simulations of two-phase flow and micromodel experiments was reported in [24]. It was also demonstrated [25, 26] that SPH is capable of simulating elastic media and that solid-fluid interactions can be relatively easily coupled in the framework of the same approach [27–30]. An interesting feature of the SPH method is its almost local character, which makes it an ideal choice for massively parallel computing, for example, with Compute Unified Device Architecture (CUDA) on Graphics Processing Units (GPU). For all these reasons and combined with the fact that different physical phenomena can be relatively easily expressed in terms of the SPH model, this method can be considered as an attractive choice for our application.
The article is organized as follows. In Section 2, the physical phenomena taken into account in our model are described. The numerical aspects concerning the discretization of the model equations are provided in Section 3. Some of the validation tests are presented in Section 4. The simulations of drainage in the sample obtained from the FIB images of the COx claystone are discussed in Section 5. The paper ends with conclusions in Section 6.
2. Physical Model
2.1. Fluid Phase
Fluid flow inside a porous medium can be described by different equations depending mainly on flow characteristics and the scale chosen for the flow description. At macroscopic scales the Darcy equation and its modifications are widely used [11, 13]; at microscopic scale, down to several nanometers of confinement, fluid flow can be described by the Navier-Stokes equations; when the confinement is limited to a few nanometers, fluid flow can only be described by molecular dynamics and stochastic methods which directly take into account intermolecular interactions. As it was mentioned before, most of the pores inside of a clayey material are of the size of a few nanometers [4]; however, these pores cannot be identified on images of real clay sample considered here with voxel size of 10 nm. Thus, we can employ the Navier-Stokes equations with a good accuracy. In the scope of the Navier-Stokes equations, the fluid is considered as a continuum matter. In the absence of sources and sinks, the conservation of mass is expressed in Lagrangian frame aswhere is the density, is the velocity, is the time, and denotes the material derivative.
Momentum conservation equation in Lagrangian frame can be written aswhere is the pressure, is the gravitational acceleration, and is the dynamic viscosity. The momentum equation for incompressible Newtonian fluid in Lagrangian frame can be stated asIn order to complete the description of fluid state at one point in terms of density , pressure , and velocity , one has to provide the equation of state which usually relates density to pressure. A common choice, when thermodynamics effects are neglected, is the ideal gas equation of statewhere is the speed of sound in the fluid and is the initial fluid density.
2.1.1. Multiphase Flow
The key problem in the modelling of immiscible multiphase flow is the description of interface dynamics. When both of the fluid phases are described by means of the Navier-Stokes equations (1) and (3), one has to assure the continuity of viscous stress tensor at the fluid-fluid interface, where the difference in interactions and/or density between both fluids causes surface tension effects. The resulting force acting on molecules at the fluid interface tends to minimize the surface area. Surface tension generates a pressure step for a curved interface between two fluids, which can be written in terms of the Laplace lawwhere is the surface tension coefficient, and are the principal radii of curvature, and is the unit interface normal.
A widespread approach consists in introducing the continuum surface force [31–35]. We follow the version presented in [36].
In the framework of this model, a volumetric surface tension force is introduced at the fluid-fluid interface. The direction of this force is normal to the interface and the magnitude is proportional to the local curvature of the interface. For the implementation of this method, the main challenge consists in calculating the local curvature of the interface. This can be done by introduction of the color field defined aswhere subscripts and denote liquid and gas, respectively. The values of the color field can be used to calculate the local curvature of the interface as follows:where the unit normal is given byThe volumetric surface tension force is expressed asTo complete the description of the interface phenomena, the contact angle can be introduced by assigning the color function value to the solid phase so thatThese values are then used in calculations of the unit normal near the solid surface.
2.2. Solid Phase
Clay material is usually composed of a mixture of various minerals. At the macroscopic scale this mixture can be considered as a continuum medium with homogeneous mechanical properties. Obviously, this is no longer true at the microscopic scale where the clayey matrix and the mineral inclusions of various sizes can be clearly distinguished. Due to their mechanical properties, the mineral inclusions can be considered as incompressible rigid bodies, while the clayey matrix exhibits mainly elastic response to the applied stress [2, 37] up to a certain point where it starts to fail. Therefore, in our model we are going to implement the possibility to take into account both rigid inclusions and elastic clayey matrix. Generally, clay matrix exhibits elastic properties of transversely isotropic material [1], which, for sake of simplicity, is replaced with isotropic elastic behaviour in our model.
2.2.1. Elastic Material
The deformations inside an elastic material are given by the strain or deformation tensorwhere is the local displacement. The Hooke law linearly relates the strain tensor to the stress tensor aswhere is the stiffness tensor. For an isotropic medium, the components of the stiffness tensor are given bywhere and are the Lamé coefficients and is the Kronecker delta function. Then, (12) can be simplified aswhere is the unit tensor. This expression can be also written aswhere is the bulk modulus and can be associated with pressure inside the elastic medium. The dynamic behaviour of an elastic medium is described byThe volumetric elastic force density components of an isotropic medium can be expressed in terms of partial derivatives of the local displacement field as
2.2.2. Damage Model
If pressure build-up or mechanical stress applied to the clay is high enough, its transport properties may be modified because of the elastic deformation of clayey matrix, as a result of, for example, dilatation of existing pores creating preferential gas paths. However, in the case of nonpercolating pore space the transport of gas phase can only take place while creating new microcracks.
Emergence of such microcracks can be modelled as a modification of mechanical properties of the clayey material under certain conditions. Since the exact models for this kind of features are not available, we propose here to reuse the ideas related to damage and fracture propagation in elastic materials. In the simplest isotropic scalar model, the process of degradation can be quantified by introduction of the damage variable [13, 14, 38, 39] which varies from 0 for intact to 1 for completely damaged material. Young’s modulus of the damaged elastic material is changed as follows [13]:where is Young’s modulus of the intact material and is the total damage variable.
In order to relate the damage variable to mechanical load inside the clayey matrix, we adopt the criterion proposed by [14] with modifications introduced by [13]. According to this criterion, two damage modes are accounted for by introduction of damage variables and responsible for tensile and compressive or shear damage, respectively.
To verify if the elastic material fails under the tensile load, the Rankine criterion is applied [14]:where is the tensile strength and is the smallest principal stress. If the criterion is satisfied, the damage variable associated with it is defined as [14]where is the residual tensile strength, is the strain at the elastic limit, and is the ultimate tensile strain at which the material can be considered as completely damaged.
To verify if the elastic phase fails under the shear or compressive load, the Mohr-Coulomb criterion is applied [14]where is the internal friction angle, is the uniaxial compressive strength, and is the largest principal stress. The damage variable associated with this criterion is given by [14]where is the residual compressive strength and is the compressive strain at the elastic limit. Following [13, 14], in (20) and (22) in case can be replaced with principal strains and , respectively. A schematic shape of the constitutive damage law given by (20) and (22) is presented in Figure 1.
It can be assumed that, for a certain load range, the material follows linear elasticity law, and outside this range a nonlinear damage behaviour is observed, which coincides with observations made for clayey materials [2, 37]. It is well known that the Mohr-Coulomb criterion works fine for materials whose tensile strength is significantly smaller than their uniaxial compressive strength which is the case for the COx clay. Although, based on the experimental data [2], other specific damage criteria are available for the COx clay, they can be fitted with the Mohr-Coulomb criterion for the sake of generality, which is done in [2]. Based on these data, the parameters adopted for COx clay are summarized in Table 1.
The total damage is a combination of and and can be calculated as [13]where the coefficients and are calculated from the local stress-strain state as described in [13, 38]. It must be noted that the scalar damage variable can be replaced with a tensor damage variable [39] which may influence the orientation of microcracks.
In order to avoid sharp restoration of the damage variable which may cause oscillations of elastic properties, the relaxation is applied to obtain a restored damage value based on value calculated by (23)where is the dimensionless relaxation constant which is set to 10 in our simulations.
2.2.3. Rigid Bodies
In their natural state, clayey rocks are highly heterogeneous, especially at the microscale where various rigid inclusions, such as quartz ( GPa) or calcite, are observed in a softer clayey matrix. Contact lines between these rigid inclusions and the clayey matrix are considered as probable sites for emergence of the new gas passages. Therefore, this component should be included into our model for further investigation.
Rigid body dynamics is fully described by the Newton-Euler equations,where is the position of the center of mass and its velocity, is the total external force acting on the body, is the mass of the body, is the angular velocity, is the tensor of the moment of inertia, is the total torque acting on the body, and are the angles between the body axis and principal coordinate axes. The first two equations describe movement of the center of mass, while the last two describe the rigid body rotation. The position of the center of mass can be found asThe tensor of the moment of inertia is given byThe total torque acting on the body can be calculated aswhere is the external force density.
The Newton-Euler equations as described by (25a), (25b), (25c), and (25d) can be integrated by any appropriate numerical method.
2.3. Boundary Conditions
In order to provide a correct interaction between phases at the solid-fluid interface, the boundary conditions must be defined. In most cases, the velocity of both phases at the contact point is supposed to be the same:This condition defines mutual nonpenetration of the phases and the no-slip boundary condition and must be completed with the condition of continuity of normal stresses at the solid-fluid interfacewhere and are fluid and solid phase stress tensors, respectively, is the unit normal pointing out of the fluid phase, and .
At the interface between the solid phases, the continuity of normal stresses is imposed:For a more realistic model, one should also add the friction force between solid components. In this work the friction force is omitted, but it may be added in the future when major phenomena will be sufficiently well understood.
3. Numerical Modeling
3.1. Smoothed Particle Hydrodynamics
The model presented in previous section must be solved numerically. In this paper, this is done by means of the Smoothed Particle Hydrodynamics (SPH) which is a Lagrangian mesh-free method originally applied for solution of astrophysical problems [18, 19]. Nowadays the method is widely used in various engineering and scientific areas [40, 41]. For our purposes it is important that the method can be successfully implemented to simulate fluid hydrodynamics, elastic, plastic, and rigid bodies as well as coupled solid-fluid systems. All these simulations can be performed within the same conceptual framework which decreases significantly the amount and the complexity of work to be done. Another considerable advantage is that the method is meshless which simplifies significantly simulations of multiphase flow with moving interfaces and also allows the description of some physical processes (such as fracture propagation) without mesh size dependency. Also the method can be easily parallelized and may profit of CUDA technology for massively parallel computations on GPU.
The SPH model follows the idea that, for a given function , its value can be found as [40]where is the Dirac delta function and the integral is taken over the entire space. The Dirac delta function has nonzero value only in one point. However, from (32) one can expect that the value of a function at a point can be approximated by this function values in close proximity of the point by means of a similar expression [40]where is a smoothing function, is a smoothing length, and is the approximation error. The integral in (33) represents the smoothed value of the function [40]The approximation error can be estimated as by taking Taylor expansion of the right part of (34). By using the divergence theorem and integration by parts, the smoothed approximation of the spatial derivatives can be obtained as [40]The higher order derivatives can be obtained in the similar way.
3.1.1. Smoothing Function
The smoothing function must satisfy certain conditions. It must have the compact support, be normalized, and converge to the Dirac delta function when tends to zero. Most of the smoothing functions are positively defined inside the support, which has certain physical meaning. In order to get the same contribution from the points located at the same distance from the center, must be an even function. It is also observed that smoothing functions with smooth first and higher order derivatives usually provide more stable numerical results, because of smaller sensitivity to particle disorder [40].
A variety of smoothing functions satisfying these and some additional conditions [40] can be constructed, so one can try to find a smoothing function which best fits the precise application needs in terms of accuracy, numerical stability, and performance. One of the most popular smoothing functions, especially in hydrodynamics simulations, is the cubic spline which is adopted in this paper, and in space it can be expressed as [40]where . It is clear that has units of inverse of volume.
Smoothing functions proposed by [42] are designed to have some useful properties. Desbrun’s spiky smoothing function [43] is appropriate for pressure forces calculation since the spiky shape of its gradient prevents the formation of particle clusters when they get close to each other. In [40] the quadratic smoothing function was proposed which is similar to the cubic spline function but is not a piecewise-defined function which may be advantageous for numerical calculation. Smoothing functions are mostly represented by polynomials for computational convenience; however, one may consider a Gaussian smoothing function [18] when looking for stability and accuracy.
3.1.2. SPH Model
In SPH model, the medium is approximated by a set of particles which are initially uniformly or randomly distributed in space. Each SPH particle is a Lagrangian particle and is characterized by its position and velocity. It also has certain physical variables such as mass and density associated with it. The particle volume can be defined asThe continuous integrals (34) and (35) can be approximated by their discrete analogues,where the summation in function term (38a) is carried out for all neighbouring particles within the distance from the th particle including ; note that in the derivative term (38b). According to [44], the error introduced by the finite sum approximation of continuous integrals can be estimated as , where is the length scale characterizing the average distance between the particles.
When is a density, (38a) is transformed toUsually, the masses of SPH particles representing the same material are equal. Another useful quantity is the number density, which can be calculated asand also represents the inverse of the particle volume. When applied directly, (39) has a drawback to underestimate the density near the material interface, which can cause consistency issues. This can be corrected by using the normalization [40]If the material is in contact with another material represented by SPH particles, the density deficiency near the interface can be corrected by considering the virtual contribution to density from another material particles, which means that summation in (39) is simply extended to neighbouring material particles located within the smoothing function support.
The direct use of (38b) for calculation of derivatives is generally avoided because of low accuracy. For example, the derivative of a constant function calculated by (38b) does not vanish. This drawback can be corrected by using the following identity [41, 44]:where is a differentiable function. Choosing various functions, one obtains different formulas for derivative term. The most frequently used functions are , , and , substitution of which in (42) yields the following expressions of SPH derivative term:The calculation of a Laplacian of some functions is often required for the transport equation. It can be obtained asThis formula also can be symmetrized just as SPH derivatives terms; however, even symmetrized formulations may suffer from changing a sign for some smoothing functions. This may cause nonphysical behaviour when (44) is applied to calculate forces. This formulation is also sensible to particle disorder [41].
Another way to calculate the divergence of a scalar function can be resumed as follows [41, 44]:where and are some scalar functions. It is useful to note that [44]where .
Using these formulations for scalars and derivative terms, it is possible to discretize various physical equations in order to obtain their SPH analogues which can be solved numerically. For example, the movement of the Lagrangian particles can be described by the following system of equations:where is the particle velocity and is the acceleration which results from the sum of the forces acting on the particle. This system can be integrated in time by using one of the explicit time integration schemes. The Leapfrog integration scheme can be stated as [44]It is similar to the velocity Verlet algorithm, which can be expressed as [44]Both methods are second-order integration schemes. For numerical stability the time step must obey the Courant-Friedrich-Levy condition [35]and particle acceleration related conditionWhen viscous dissipation term is present, the following condition must be satisfied:
The SPH is usually used to solve systems of partial differential equations which relate the time derivative of variables of interest to spatial derivatives. Since the equations of interest here describing fluid flow, elastic body deformation, and solid body movement satisfy these conditions, let us consider their discretization in detail.
3.1.3. SPH Discretization of Navier-Stokes Equations
In the framework of SPH, the continuity equation (1) can be discretized in various ways depending on the formulation adopted for the spatial derivatives (43a), (43b), and (43c). The most widespread variants are [41]The first one is preferable when simulating a multiphase flow with large density contrasts since it provides better accuracy [41]. It is worth mentioning that these expressions conserve the mass only approximately (error of the order of ) [44]. Another possibility to follow the evolution of fluid density consists in using (39) which conserves the mass exactly.
The right part of the momentum conservation equation (3) can be discretized by using (43a), (43b), and (43c). The pressure gradient term can be represented in SPH form by one of the following expressions:Both expressions provide similar results in terms of accuracy; however, the second expression is more stable with respect to particle disorder [41]. The use of (43a) and (43b) must be avoided for discretization of the pressure gradient term, since it becomes unstable for negative pressure values. The viscous dissipation term can be discretized by using (44) or (45). The use of the (45) should be preferred for better accuracy and numerical stability, which leads us to the following expression:At the interface between two fluids (of when viscosity term is a spatial variable) one can use the original expression for viscous dissipation term based on (45) which corresponds to the arithmetic mean approximation; however, for better precision it is preferable to use the expression based on the harmonic meanThus, the discretized momentum equation can be written as follows:where is the surface tension force.
In weakly compressible formulation adopted in our model, the incompressibility of the fluid is imposed by the proper choice of the speed of sound in the equation of state (4). It must be prescribed based on the desired density variation, keeping in mind that high values of the speed of sound lead to small time steps and noisy pressure distributions [44]. An alternative approach consists in using the projection method [45].
Two principal approaches exist to model the surface tension by using the SPH model. The first one consists in using the pairwise interparticle force which mimics the intermolecular attraction forces [44]. This allows obtaining surface tension effects naturally by adjusting the interaction force constants for different types of particles. However, the introduction of such a force can modify the viscosity of fluids and the equations of state [44].
In our model, the surface tension effects are introduced by discretizing the continuum force model (Section 2.1.1) in terms of SPH. First, a certain value of color function must be attributed to the fluid and gas SPH particles. The simplest way to do this is to use constant values and for fluid and gas SPH particles, respectively. Based on these values one can define a smoothed value of the color function associated with SPH particleAlso, the color function can be calculated based on (6). The color function normal is calculated asThe interface local curvature given by (7) is calculated asIn order to improve the accuracy of the curvature calculation, we employ the concept of reliable normal introduced in [35]. Then, the continuous surface tension force is calculated by (9). Finally, the contact angle is prescribed by assigning a color function value to the solid and elastic SPH particles based on (10). Integration is performed using the Leapfrog integration scheme (48a) and (48b) while the fluid densities are updated using (39).
3.1.4. SPH Discretization of Elastic Bodies Dynamics
The most frequently used method to model elastic bodies by means of SPH consists in integration of rate of change of the deviatoric stress tensor [25, 26, 28–30] in the framework of the hypoelastic constitutive model. Another approach consists in calculation of elastic forces from expansion of the elastic strain energy [46, 47]. In our model, for the sake of simplicity we discretize the equation governing the dynamic behaviour of the elastic medium (16) by means of SPH.
In order to approximate the force components by SPH, second derivatives of the local displacement can be calculated using the same approach as for the first derivatives. However, this formulation has several disadvantages [41], such as high sensitivity to particle disorder. A better alternative is to use the following approximation of second derivatives [41]:where the Greek indexes and indicate spatial dimensions and is the Kronecker delta function.
When elastic forces acting on elastic particles are calculated, the particle trajectory can be easily obtained by numerical integration of corresponding equations (47a) and (47b) of motion. In order to limit vibrations and dissipate elastic energy, a viscous damping term similar to the one used in the Navier-Stokes equation can be added to (17); another possibility consists in introducing the viscosity-like term in the Verlet algorithm [48]. Both options are implemented in our model.
3.2. Discretization of Rigid Body Dynamics
A rigid body can be discretized with a set of Lagrangian points distributed inside the body volume and at the body surface. This spatial description allows a natural coupling of discretized equations of rigid body dynamics with the other equations in the frame of SPH [27].
In the context of our model, the movement of rigid inclusions results from forces applied on their surfaces by other phases as well as from contact forces between inclusions themselves. In order to decrease the memory use and to save calculation time, one can describe the rigid body by taking into account only the points lying at the body surface and in its close proximity (i.e., the points which may experience the external forces from other phases) without any loss of accuracy since the expected rate of rotation is low, and the final position of the rigid grain is most likely defined by its shape and distribution of external forces at the grain surface rather than by distribution of mass inside this body.
The center of mass of the rigid body represented by a set of points of mass can be found fromThe total external force acting on the rigid body is given byThe tensor of the moment of inertia is obtained fromThe total torque acting on the body can be calculated aswhere is the acceleration of the point due to external forces.
The Newton-Euler equations (25a), (25b), (25c), and (25d) can be integrated using the Leapfrog algorithm as follows [27]. First, the acceleration of the center of mass and the angular acceleration of the rigid body are calculated as [27]Then, the velocity of the center of mass and the angular velocity of the rigid body are updated [27]Finally, the position of the center of mass is updated by [27]The velocities and the positions of the rigid body points can be obtained as follows [27]where is the rotation matrix. It can be calculated as [27]where , , and are the matrices which yield the rotation about principal coordinate axes
3.3. Boundary Conditions
Interaction of the fluid phase described by SPH with solid objects can be obtained in various ways. Boundary forces taking into account viscous stress and fluid pressure are described in [49, 50]; coupling of SPH fluid with rigid solid bodies can be found in [27], while coupling with deformable solids is considered in [51]; coupled models where both the fluid and elastic solid are described in terms of SPH can be found in [28–30].
SPH particles belonging to two different phases or two different objects of the same phase are considered as boundary particles if the distance between them is smaller than smoothing length . Interaction between these particles is defined by the boundary conditions described in Section 2.3.
In our model, we adopted the simplest way to satisfy the conditions (29) and (30) at the solid-viscous fluid interface by extending the summation in momentum conservation equations for all types of interacting particles [29]. For nonviscous fluid (the viscosity is equal to zero) a more complicated approach is necessary [27, 29]. Pressure values for solid phase particles can be obtained similarly to (4) by introduction of the speed of sound:where is the bulk modulus of the solid phase and the initial solid phase density. for the elastic phase is defined by its elastic properties, and for the rigid phase is set 8 times larger than that for elastic phase. Using this approach, the no-slip boundary condition is imposed at the solid SPH particles.
Interaction between boundary particles belonging to different solid phases or objects can be modelled by introduction of pairwise interparticle forces based on Lennard-Jones potential [52], or repulsive forces modelling conservative full elastic interaction [27, 41]:where is the force acting on particle , is the symmetric function which defines the interaction strength, and is the local boundary normal. We have also successfully tested the interparticle interaction term based on the repulsive part of the pressure term from the momentum conservation equation.
3.4. XSPH Correction
In our model, we implement the XSPH correction which is a frequently used technique [25, 28, 30, 41] allowing the significant increase in the simulations stability and decrease in nonphysical noisy fluctuations. It consists in using the corrected velocity for update of the SPH particles positions in the adopted numerical integration scheme instead of the usual particle velocity which is used in the momentum equation. The corrected velocity is calculated by [41]where is a constant parameter which varies from 0 to 1 and usually is set to 0.5. The summation in (76) is taken over the particles of the same phase. This modification makes particles move with velocities which are not too different for close particles, while the linear and angular momenta are not modified [41].
4. Validation Tests
4.1. Poiseuille-Couette Flow
The capacity of the model to represent correctly important physical phenomena was tested in various validation tests. A simple test to verify the capacity of the model to simulate fluid flow is the transient Poiseuille-Couette flow between two plane walls and . The fluid velocities imposed at these walls are and , and the analytical solution for the velocity profile is given by [53]
For SPH fluid simulations, the cubic spline smoothing function is selected, which is most often used for fluid flow since it provides most accurate results. Fluid particles are initially placed at regular cubic lattice nodes with spacing; the solid walls are represented in the same way by immobile SPH particles. The no-slip boundary conditions are applied at these solid phase particles. Spatially periodic boundary conditions are imposed along - and -axis. The smoothing length is m, the distance between the walls is , the body force is , and the imposed velocity is . The fluid density and viscosity are and Pa·s, respectively. In all cases shown in this paper, small values of body force are imposed following [53].
The obtained velocity profiles at different times are compared with analytical solutions in Figure 2. The agreement between the numerical and analytical solution is very good. The maximal velocity error was smaller than 0.4%. It must be noted that the accuracy of the simulation depends on the smoothing function and the way chosen to discretize the components of the Navier-Stokes equation, especially the viscous dissipation term. The best precision was obtained for the cubic spline and the viscous dissipation term discretization as given in (45).
4.2. Periodic Arrays of Spheres
The capacity of the SPH to simulate fluid flow in porous media can be tested for a nontrivial geometry represented by periodic array of spheres for which an analytical expression of permeability is available [54, 55].
The permeability of a porous medium can be found from Darcy’s lawwhere is the viscosity, is the mean fluid seepage velocity, is the pressure gradient, and is the body force. For a spatially periodic simple cubic array of spheres, the permeability can be found analytically [54, 55] aswhere is the volume of the unit cell, the radius of the sphere, and the inverse of the dimensionless force exerted by the fluid flow on the array sphere [55]:Values of can be found in [54, 55] for various values of solid volume fraction
The ability of SPH to simulate fluid flow for periodic array of spheres with a good accuracy was demonstrated in [53] for accurately discretized geometries; about 30 SPH particles across the pore throat are necessary in order to obtain the results within about of analytical solution. In our case, we are interested in estimation of possible differences between the analytical solution and the simulation result for relatively coarse geometry and particle (i.e., the number of particles along the smoothing length ) discretizations, since for the clay samples images obtained by the imaging methods, the overall sample size is too big to discretize everywhere the pore space with a large number of SPH particles.
For numerical simulation the unit cell size is , where m is the smoothing length. The center of the solid sphere is located at the center of the unit cell. Fluid and solid SPH particles are initially placed at regular cubic lattice nodes with interparticle spacing. The no-slip boundary conditions are applied at the middle distance between the solid and fluid particles, while spatially periodic boundary conditions are imposed at the sides of the unit cell. In order to mimic a macroscopic pressure gradient, the flow is driven by the imposed body force , the fluid density is , and viscosity is Pa·s. The speed of sound is . High values of speed of sound result in very small time steps which slow down the simulation. Therefore, the common approach when using SPH is to adopt smaller values of speed of sound [44] keeping in mind that they should be large enough to ensure the relative fluid density variation inferior to 1% during the simulation.
In order to test the sensitivity of the method to various values, three values of solid concentration are considered , , and . Coarse solid spheres are discretized with elementary cubes of size ; that is, each elementary cube contains 8 solid SPH particles. The elementary cube is considered solid if its center satisfies the condition . In order to test the influence of the sphere surface discretization, spheres are also discretized with solid SPH particles which satisfy the condition , where is the location of the particle. This way provides a refined discretization and equivalent to discretizing the sphere with elementary solid cubes of size , each of which contains a single SPH solid particle located in its center. The difference between the coarse and refined discretizations is displayed in Figure 3 for . The sensitivity of the method to particle discretization is tested by performing the same series of simulations for the interparticle spacing.
Numerical errors for the calculated permeability values are summarized in Table 2. For the interparticle spacing, errors of order of are observed with little influence of the sphere surface discretization, which means that precision is mainly affected by the coarse discretization of flow equations in terms of SPH. For the interparticle spacing, the errors due to particle discretization are less significant, and clear influence of the sphere surface discretization can be observed. The errors of order of 20% are observed for the coarse sphere surface discretization, while for the refined sphere surface discretization the error is less than 10%. As mentioned at the beginning of this section, the pore throat must be very well discretized in order to obtain the permeability value within few percentages of analytical solution when using SPH.
It must be noted (see Table 3) that the precision of the simulation results depends strongly on the number of SPH particles within the smoothing length, which is closely related to the initial interparticle spacing. The larger number of particles provides a better discretization of equations in terms of SPH and consequently more precise results at the expense of significantly increased computational cost. The influence of the number of SPH particles within the smoothing length on the simulation precision can be stronger than accuracy of geometry discretization. Thus, the permeability calculated for , , setup is only 7.02% accurate compared to 3.97% accuracy for , , setup, in spite of larger total number of SPH particles and better sphere surface discretization.
Based on these simulation results, we can expect for our drainage simulations the precision of about 30%, which is satisfactory for simulations performed for coarse geometries with interparticle spacing.
4.3. Two-Phase Poiseuille Flow
In order to test the capacity of the model to simulate multiphase flow, let us consider a two-phase (liquid-gas) Poiseuille flow between parallel plane walls located at and . The liquid-gas interface is also a plane located at . The gas is located in and the liquid in region, respectively. The analytical solution for fluid velocities is the following:where is a body force or a pressure gradient.
The simulation domain is the same as for the considered above transient Poiseuille-Couette flow simulation. The no-slip boundary conditions are applied at the solid phase SPH particles, while spatially periodic boundary conditions are imposed along - and -axis. The cubic spline smoothing function is used with smoothing length , the gas/liquid region width is , and the body force is equal to . The liquid density is and viscosity Pa·s; the gas density is and viscosity Pa·s.
The numerical velocity profile is compared to the analytical solution in Figure 4. A very good agreement is obtained with maximal velocity deviation smaller than 1%. Several simulations performed for even larger ratios of water/gas physical properties also demonstrate a very good accuracy, which proves the capability of the SPH model to simulate accurately multiphase mixtures with highly contrasting physical properties.
4.4. Laplace Law
In order to test the continuum force algorithm described in Section 2.1.1, the compliance with Laplace’s law for spherical bubbles is tested. Simulation is initiated with a cube of size consisting of gas SPH particles, placed at the center of the cubic unit cell of size and surrounded by fluid SPH particles. Under the influence of surface tension forces, the gas cube is transformed into a sphere of radius . After the equilibrium state has been obtained, the pressure difference across the fluid-gas interface is calculated and plotted for different sphere radii in Figure 5. The difference between the numerical and analytical pressure drop across the interface does not exceed 2% even for the smallest bubbles.
4.5. Reaction of Elastic Phase to Fluid Pressure
The elastic-fluid coupling can be tested by calculation of the relative volume variation of an elastic cube due to the pressure step applied by the surrounding fluid:
The volume variations for several pressure steps are compared with analytical solution given by (83) in Figure 6. Relative error varies from 0.33% for small volume variations to 11.4% for relatively large volume variations. This error is mostly due to the fact that the adopted approximation for second-order derivative (61) in elastic force calculation generates errors in the range of 5–10% for large deformations . However, from Table 1 it follows that clayey matrix elastic material fails at much lower strains which limits the influence of second-order derivative errors on the simulation results when applied to an argillite sample.
5. Application to Clay Sample
In order to test the ability of our code to simulate gas-water flow inside a realistic pore space in repository conditions, we make a very first application to a small sample of Callovo-Oxfordian (COx) clay.
5.1. Characterization of the Sample
The COx claystone sample EST27405-1 is represented by a series of 103 images of pixels size corresponding to consecutive scans orthogonal to -axis obtained by FIB method [56]. The voxel size is . The images are treated based on gray level in order to recognize the pore-space and clayey matrix. The porosity of the sample is . A more complex morphological analysis was also applied to the sample in order to distinguish various solid inclusions which are most likely small grains of calcite and quartz [56]. However, for this first approach, in the following we will consider the nonporous part of the sample as a continuous elastic phase without rigid inclusions. This simplified approach allows us to use macroscopic mechanical properties of COx rock for solid matrix parameterization which are summarized in Table 1.
In order to decrease the computational costs, a coarsening procedure was applied to the original sample consisting in replacing of 8 neighbouring voxels with a single two times larger voxel whose nature (pore space or clayey matrix) was defined based on the majority rule. Two samples of and size were obtained by applying this procedure once and twice, respectively. The pore space of the sample consists of several isolated pores. The percolating pores can be assimilated to fissures which connect opposite sides of the sample along the -axis. They are surrounded by multiple small disconnected pores. The percolating pore components are displayed in Figure 7. As it was mentioned before, most of the pores inside the claystone are of the nanometer size; therefore, some connectivity information can be missing because of insufficient resolution of FIB imaging.
Due to significant computational costs (one simulation for size unit cell requires more than 3 million time steps and lasts about 30 days on our basic GPU card), only one percolating pore component was considered (see Figure 7). Since the individual shape of the selected percolating component influences the details of the drainage, the obtained results cannot be generalized for all the percolating pore components; however, we expect that they will provide an insight into the phenomena which may occur in these conditions. For simulations we considered a percolating pore (the blue pore in Figure 7) extracted from once and twice coarsened samples. The chosen pore has a fracture like shape with larger entry at the top side and smaller outlet at the bottom side. The results obtained for the percolating pore (see Figure 8(a)) extracted from the sample coarsened twice are very noisy due to insufficient particle resolution; therefore only the results for less coarsened sample will be discussed here.
(a)
(b)
For the present study, the OpenMP [57] and CUDA GPU parallelization techniques were applied with a significant gain in numerical performance. One drainage simulation as presented further in the chapter took about 30 days on NVidia Quadro K5200 standard graphic card, which is about 52 times faster than a single thread processor version. However, this performance could be greatly improved (speed up of about 12 times) by using dedicated GPU units like Nvidia Maxwell P100.
5.2. Numerical Simulations
Several drainage tests are performed in order to study the model behavior under various conditions. A unit cell for simulations is obtained by cutting a part of the sample containing the selected percolating pore component as well as small surrounding disconnected pores. In order to avoid fluid movement across the cell borders, an elastic layer is added to each lateral side of the sample. Finally, the sample is placed between horizontal solid walls in order to prevent its spatial movement during the simulation as displayed in Figure 8(b). Pores inside the sample are initially filled with water phase, while the gas phase is located in the reservoir outside the sample. The height of the reservoir is . The volume of gas is sufficient to completely drain the sample. Spatially periodic boundary conditions are applied to the unit cell along the -axis. A body force is applied to the gas in the reservoir outside the sample, which creates a pressure difference between upper () and lower () sides of the sample (see Figure 8(b)) and causes drainage. Finally, the confining stress is applied to the lateral sides of the sample to study its impact on gas migration and damage of the clay. A cross-section of the unit cell of size is given in Figure 8(b). One can see that the gas phase (green) penetrates the sample and liquid phase (yellow) is leaving the sample from through the bottom. The entry and exit pore throats are not visible since they are located in different cross-sections than the one displayed in Figure 8(b).
Usually, it may be possible to prescribe real values of physical parameters of all the phases considered in the model and to solve numerically the corresponding equations. However, in this case the time step would be prohibitively small ( s) due to the explicit integration scheme adopted in the model. Also, the time to completely drain the sample can be very long. For this reason, the physical parameters values are selected based on the dimensionless numbers which define the dominant physical phenomena.
Two-phase flow in porous media can be characterized by several dimensionless numbers:The Reynolds number which represents the ratio between the inertial and viscous forces is very small for gas flow in clay due to small flow velocity and small characteristic length which is typically the pore size. The mobility number gives the ratio between water and gas viscosities. The capillary number Ca measures the ratio of the viscous forces to the surface tension forces and is also usually very small for the considered situation. The Bond number characterizes the ratio between external forces (pressure gradient, gravitation, and body forces) and surface tension and describes the shape of bubbles. The dimensionless numbers values expected for the gas-water flow in real clay microfractures [17] and the values applied in our simulations are presented in Table 4. It must be noted that the and numbers given in [17] were calculated for micrometer pore sizes. However, the direct extrapolation to nanometer size pores would lead to still smaller numbers, but it does not change the flow regime.
It can be observed that the Reynolds and capillary numbers observed in our simulations (for mean gas velocity before the percolation) are several orders of magnitude higher than the expected real ones as a consequence of significantly reduced viscosity ( Pa·s, Pa·s) and surface tension ( N·m−1) values; however, they are still sufficiently small and should describe adequately laminar two-phase flow dominated by capillary forces.
As it was previously mentioned, in weakly compressible formulation the speed of sound in the fluid must be high enough to effectively penalize density variations and at the same time small enough to allow reasonably large time steps. The adopted value m·s−1 for both fluid phases satisfies these conditions with relative fluid density variation of order of 0.12% observed in our simulations for the pressure difference = 0.3 Pa imposed between sample sides. In order to accelerate the simulations and to more easily observe the effects of damage, the gas pressure at the inlet is set to = 15 MPa which is slightly higher than the maximal pressure which could be expected within the repository in COx (the vertical confinement is of about 12 MPa and the gas entry pressure to the sound argillite may be as high as 10 MPa [58]). Simulations with smaller pressure values were conducted but will not be presented in this paper. So, the ratio between and simulation = 0.3 Pa pressure values at the sample entry is . Since in these conditions the fluid pressure values calculated by (4) are much smaller than those expected in real conditions, the elastic properties from Table 1 are scaled so that the ratio between the gas pressure at the sample entry and Young’s modulus of the elastic phase used in the simulationstays the same as in repository conditions where . For the sake of better understanding with respect to the aimed application, we are going to comment all simulations by scaling the numerical results back to the repository conditions values (using the factor ). These rescaled values will be preceded by .
5.2.1. Drainage with and without Damage
The pressure difference is imposed between upper and lower sample sides, while the confining stress of is applied to the lateral sides. The drainage tests are performed with and without damage model in the elastic phase. Another test is performed with a completely incompressible rigid porous matrix. Due to the pressure difference between the opposite sides, gas phase enters the sample and pushes the water phase out. Since the gas phase is nonwetting, it propagates mostly inside the wetting water phase.
Since the ratio is very small, only small local deformations will be observed in the intact elastic matrix; however, when accumulated through the sample, these deformations will result in faster drainage compared to the completely incompressible rigid porous matrix as it can be observed in Figure 9(a) where the evolution of saturation with time is shown. It can be observed in (Figure 9(a)) that at the very beginning the gas enters the sample with a larger velocity regardless the application of the damage model. This is caused by the establishment of the pressure gradient through the sample, which causes deformation of the elastic matrix increasing the available pore space. The effect is more pronounced for the model with damage due to the softening of the elastic matrix. The compressibility of the fluid phases can also cause the similar effect if the speed of sound is not high enough.
(a)
(b)
After the initial compression phase, the saturation evolves almost linearly until the percolation, where the curve changes the slope and the drainage continues with a decreasing velocity due to the fact that water phase is now mainly located near the walls. It can be expected that the sample will be completely drained after a longer time due to viscous drag force. Incorporation of the damage model (Figure 9(a)) increases the drainage rate at all the stages of simulation and the nonwetting fluid percolation time is accelerated (from s to s). However, for this particular geometry of pores, which are percolating already in the initial sample, this is the most important change in behaviour, which stays qualitatively similar with or without damage.
5.2.2. Drainage under Different Confining Pressures
The next series of drainage tests is performed for various confining stresses , 8 MPa, 6 MPa, and 4 MPa applied to lateral sides and for the same pressure difference .
Same as for (Figure 9(a)), the initial compression phase can be observed in Figure 9(b). It is more pronounced for where the largest deformation of the elastic matrix by the entering gas is expected due to the smallest applied confining stress. What is more, in this case we note the largest number of damaged particles (Figure 10) causing the elastic matrix softening. As can be seen in Figure 9(b), these effects are largely attenuated with increasing confining stress.
(a)
(b)
The percolation times and gas saturation at percolation are summarized in Table 5. As expected, the percolation occurs earlier for smaller confining stresses due to the larger deformations and as a consequence due to increased permeability that facilitates drainage. The largest corresponds to the incompressible porous matrix. The gas saturation at percolation is close to 0.5 for all the cases and increases slightly with confining stress.
Evolution of the number of damaged points and of the total damage is presented in Figure 10. After an initial rapid increase, number of damaged points continues to grow with a reduced rate for all considered confining stresses. A slightly faster increase can be observed just before percolation which is due to the increase of the water pressure (see Figure 13). After percolation the number of damaged points stays almost constant with a very slow grow rate. The number of damaged points is larger for small confining stresses than for high confining stresses since it is more difficult to fracture compressed elastic matrix which results in slower drainage in the last case. Initially damage occurs at the edgy parts of the elastic-fluid interface where the most of the traction stress is accumulated; then it propagates through the elastic phase in the direction of the sample sides (see Figure 11). Also a lot of damaged elastic points were observed in a region where the pore space approaches the sample lateral sides. The total accumulated damage presented in Figure 10(b) follows the same pattern as the number of damaged particles before the percolation, after which it decreases rapidly to a new relatively stable level as a consequence of fast relaxation of gas and water pressures after the percolation (see Figure 13(a)). As it can be observed in Figures 11(c) and 11(d), this reduction of pressure causes also a decrease of number of damaged elastic particles: first particles to restore their original elastic properties are those located the furthest from the damage departure location.
(a)
(b)
(c)
(d)
The seepage velocities of fluid phases during drainage are plotted in Figure 12. At the initial stage, the velocities are high due to multiple factors mentioned above. During the linear stage of drainage, the velocities are almost constant and close to each other for all the considered cases. Water velocity gradually decreases before the percolation, since the water which is left in the sample is located close to the pore walls where its velocity tends to zero. After the percolation it stabilizes at a significantly lower level. As expected, gas velocity increases several times after the percolation.
(a)
(b)
(a)
(b)
The water and gas pressures averaged over the percolating pore volume and the capillary pressure are presented in Figure 13. The water phase pressure increases rapidly until the establishment of the pressure gradient through the sample; then it grows gradually until the percolation after which it drops and stabilizes at a new level. The gas phase pressure decreases at the initial stage of drainage; then it stays relatively constant till the percolation after which it drops and stabilizes at a new level similarly to the water phase. In some cases, the drop of the pressure can lead to reclosure of the percolating pore passes opened during the gas passage. Capillary pressure decrease with growing gas saturation is related to the geometry of the percolating pore, in which aperture increases from entry to the central part where it reaches the maximum value. Capillary pressure curves are almost identical for the considered cases starting from . Just before the percolation the capillary pressure decreases faster than at the earlier stage and stabilizes after the percolation. This decrease is related to the growth of the water pressure which is caused by the gas pushing the water under low pressure from the bottom part of the sample. The results for the smallest confining stress are quite noisy due to important number of damaged elastic points which can cause significant pore surface displacements influencing the stability of the numerical method.
It must be noted that the numerically obtained capillary pressure values displayed in Figure 13 are about 3.6 times smaller than those which can be expected for this pore size in the case of water/air interface with , since water/air surface tension value was decreased 3.6 times in simulations, in order to speed them up.
5.2.3. Damage Model versus Numerical Fracturing
From the presented simulations it can be concluded that the mechanical interaction between the solid matrix and the fluid phases can influence the drainage process in clays in repository conditions. It has been seen as well that the resulting damaged area was relatively wide and did not localize into new fractures. We have tested several approaches in order to better localize damaged particles. We observed that the details of damage model greatly influenced simulation outcome. The results presented in previous chapters were obtained supposing the reversibility of damage variable; that is, elastic properties of damaged particles can be restored if this is indicated by the damage criterion. Also, in order to preserve the volume of the elastic phase, damaged elastic particles continue to interact by pressure gradient term (second term in (15)). Under these conditions, the number of damaged elastic particles was bigger for smaller confining stresses, and these particles formed connected regions as it can be seen in Figure 11. These regions can be considered as containing microfractures, which present apertures too small to appear at the model scale. For these reasons, the fluids particles cannot penetrate in these microfractures and create new connections between pores. However, these passages can be created if the fluid pressure force acting on the pore space surface exceeds the confining stress as illustrated in Figure 14(b) where a formation of new connections visible at the model scale is shown. This crack is big enough to be penetrated by gas particles. It can be observed that the pore space is noticeably larger than at the initial time (see Figure 14(a)).
(a)
(b)
(c)
(d)
When using SPH model for the elastic phase, a fracturing can be observed even without applying a damage model. For the simulations presented in Sections 5.2.1 and 5.2.2 the ratio is very small, which means that very small local deformations will be observed in the intact elastic matrix. Therefore, without application of damage model, there would be no initiation of new fractures. However, for larger values, for example, , the deformations may be big enough to separate the SPH particles at distance larger than , which can be interpreted as a creation of a new fracture at the model scale. An example of such a simulation, without application of a damage model can be seen in Figure 14(c). This phenomenon is usually called the numerical fracturing. It allows in some conditions to localize nicely new fractures; however, it is impossible to parameterize this process in order to fit it to a particular porous material.
A still different behaviour can be observed if damaged elastic particles are considered as ghost particles, which no longer interact with any other particles. In this case, the release of the stress accumulated inside the elastic matrix may result in formation of narrow fractures as displayed in Figure 14(d). However, the total volume of the elastic phase is no longer conserved in these conditions.
5.2.4. Solid Inclusions
In order to consider the possible influence of the rigid inclusions on the flow conditions, two rigid spheres were incorporated into the elastic matrix as displayed in Figure 15. The only implemented interaction between solid and elastic particles was a pressure gradient repulsion. All other materials parameters were the same as in Section 5.2.2. It can be observed (see Figure 15(a)) that if the sphere is not deeply incorporated into the elastic matrix, its position can be altered by the flow damaging the surrounding elastic particles, and creating new narrow passages between the rigid inclusion and the elastic matrix which can be penetrated by fluid phase particles. For deeply incorporated rigid inclusion (see Figure 15(b)), very little influence on the surrounding elastic solid is observed. It can be concluded that the presence of rigid inclusions in the model provides, as expected, a different way of creating new paths for gas migration. In order to obtain a better description of the interaction between rigid inclusions and elastic matrix, friction forces must be included at the interface which will decrease to a certain extent the tangential movement of rigid inclusions in contact with elastic matrix.
(a)
(b)
6. Conclusions
The aim of this paper was to propose a modelling approach for the dilatant two-phase flow which is assumed to occur during gas migration within initially water saturated indurated clay rock. This type of flow, unlike classical viscoelastic flows, require an interaction of flow with a deformable solid matrix.
We developed a numerical model which allows simulating a two-phase fluid flow inside a porous medium and taking into account the mechanics of the solid matrix, possibly composed of two different media: a linearly elastic continuous medium and isolated rigid inclusions. Then a damage model is introduced in order to take into account in a controlled manner the degradation of elastic material properties due to local stress-strain state. The model can be applied to the FIB-SEM images of the real clay samples in order not only to study various phenomena at pore scale but also to calculate effective transport properties, that is, permeability for varying pressure gradients and confining stresses taking into account deformation and degradation of the elastic matrix.
The model was solved in the framework of Smoothed Particle Hydrodynamics, which is a Lagrangian mesh-free method. Classical methods to numerically simulate fluid flow are mostly based on Finite Volume Methods (FVM), while in solid mechanics the Finite Element Methods (FEM) are widely applied. These methods allow solving several problems more efficiently than SPH; however, in some complicated situations they may suffer from drawbacks related to their mesh based nature. This includes the cases with deformable and moving interfaces, dealing with large deformations and propagation of fractures. Use of SPH model is very appropriate in this context since it provides the opportunity to describe various types of interactions between different materials within the same framework. The presented SPH model was coded using C++ programming language implementing all the phenomena mentioned above. The SPH is computationally expensive. The main load comes from calculation of neighbouring particles contributions when evaluating densities and forces. Its cost grows significantly with the number of particles located within the smoothing length. Therefore, the original code was parallelized using OpenMP [57] and CUDA GPU techniques. Parallelization of SPH can be performed relatively easily, since densities and forces are calculated by summing contributions from neighbouring particles; therefore, each particle can be treated separately by a computational thread.
For the present study, the OpenMP and CUDA GPU parallelization techniques were applied with a significant gain in numerical performance when run on an ordinary graphical card.
With several elementary tests we have shown that the model is sufficiently stable and robust for use in conditions representative for application to gas-water flow in argillites and that a limited resolution was sufficient to obtain permeability values with a reasonable precision (about 30%).
The model was then applied to simulate drainage in a small sample obtained from FIB reconstruction of COx claystone. The simulation parameters were set to approximate as close as possible the experimental conditions while keeping in mind numerical constraints. The , , , and number values of our simulations describe qualitatively two-phase flow regime in clay microfractures. The elastic properties are scaled so that proportions between them and are the same as in repository conditions. The values of are slightly larger than expected in repository conditions in order to speed up the simulations.
For the particular geometry of our sample containing a percolating pore, the overall drainage process stayed similar, though the softening of clay though damaging resulted in accelerated drainage velocity. Also, the obtained relation for the suction curve of an isolated pore is decreasing with gas saturation, which is the opposite from the experimental relationships obtained for macroscopic samples. One remarkable feature of capillary pressure curves in Figure 13(b) is that they are almost identical for saturation of gas larger than 0.1.
The capability of the model to create new pathways due to localized damage was also investigated. We have demonstrated that for a small sample under applied conditions it is difficult to create macroscopic fractures and damage zone stays diffuse. However, taking into account the rigid inclusions resulted, as expected, in creation of new pores between them and the elastic matrix.
We conclude that the developed SPH code is capable of representing all required phenomena for FIB-SEM based geometries in conditions close to the real ones. To be able to obtain more quantitative information about drainage in argillites some further work needs to be conducted. In particular, it is not completely clear whether the macroscopic damage model is appropriate for use if the solid matrix is divided into clayey matrix and rigid inclusions (for COx they constitute about 40% of volume). At the nanometric scale also the mass exchange phenomena can be very important, especially the possibility for gas to dissolve and for water to evaporate, since in tortuous pore space the continuous phase flow may be slowed down or even stopped. Finally, our model should be applied to larger samples in order to minimize the boundary effects and to provide data independent of individual pore space geometry.
Disclosure
A part of this work was presented during JEMP 2016 conference, held on 12–14 October 2016 in Anglet, France.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The financial support from NEEDS MIPOR is gratefully acknowledged. The help of Alain Genty (CEA) who provided a GPU cluster for simulations is greatly appreciated.