Numerical homogenization method in the modeling of gas hydrate bearing sediments

The study of the mechanical behavior of gas hydrate bearing soils represents a major interest. The hydrate inclusions in sediments change their microstructure and their mechanical properties with it. We developed a numerical homogenization code in order to simulate the macro-mechanical response of a periodic unit-cell using local elastoplastic laws and complex geometries to define the microstructural components of the material. We applied it to a real image of fine-grained sediments containing gas hydrate veins. This technique allowed us to study the effect of an elastoplastic soil matrix and of different volume fractions of gas hydrates with complex shapes on the overall response of the material.


Introduction
Gas hydrates are unstable crystalline compounds consisting of gas molecules trapped in cages of water molecules. They naturally form under high pressure and low temperature conditions, in oceanic sediments around continental margins or in permafrost regions where water and methane gas are present in sufficient quantities. Natural gas hydrates bearing sediments (GHBS) represent a potential energy resource and a geohazard that it is essential to study and control. However, it is difficult and expensive to test natural core samples [1], [2] and most of the experimental knowledge on GBS was obtained by testing synthetic samples in laboratory under in situ conditions [3,4,5,6]. Those results showed that hydrate inclusions can change the microstructure and physical properties of the sediments [7,8]. Several macroscopic multi-physical models have already been applied to GHBS but it is difficult to develop a reliable mechanical constitutive model considering the experimental difficulties.
Using homogenization to overcome inexpensively these difficulties can help investigate the impact of the heterogeneous microstructure on the macroscopic response of GHBS. Homogenization techniques give the macroscopic quantities defining the overall behavior of a composite material by averaging the local or microscopic mechanical fields over a given geometry. A particular numerical homogenization method is based on the use of Fast Fourier Transforms (FFT) to solve the local mechanical problem of an elementary composite cell, under periodic boundary conditions [9], [10]. This original method has led to the development of different FFT-based computational schemes [11][12][13]. These methods do not need any mesh which allows the use of more complex microstructures compared to finite element methods (FEM), and even the use of real images.
In this paper we apply a FFT-based homogenization scheme to study the macroscopic mechanical response of fine-grained host sediments with occurrences of gas hydrate veins. The geometry of the microstructure is based on a real image acquired from micro CT scans of GHBS from Krishna-Godavari Basin [14], and modified to obtain various gas hydrate volume fractions.

Numerical homogenization
Homogenization methods provide the macroscopic behavior of a heterogeneous Representative Elementary Volume (REV). They can help model the mechanical behaviour of GHBS in the absence of reliable experimental properties. In this section we present the theory at the base of the methods we chose to use in our in-house numerical homogenization code.

Periodic homogenization
We remind that the strain and displacement fields are separated as follows in the theory of homogenization, x being the coordinates in real space over a periodic volume Ω: with local fluctuating terms * ( ) and * ( ) having a null average over Ω as they are periodic functions of space.
is then the mean strain over the entire periodic volume Ω: 〈 ( )〉 = . The periodic homogenization consists in solving the local mechanical problem defined for a unit-cell Ω submitted to a macroscopic strain load and periodic boundary conditions. It is governed by the following equations: where the first equation represents the momentum balance; the second equation represents the non-uniform microscopic constitutive law that relates the local stress field to the local strain field ; the third equation represents the compatibility equations and the last two lines represent the periodic and anti-periodic conditions for the local displacements * and traction vectors ⋅ at opposite boundaries. Averaging the local stress field solution to the local problem over the unit-cell gives the macroscopic stress: Σ = 〈 ( )〉.

FFT-based methods
The introduction of a homogeneous elastic reference material, characterized by a stiffness tensor ℂ , and a polarization tensor field leads to the Lippmann-Schwinger (LS) equation as the solution of the local mechanical problem [9]: where ( * ) is a convolution product, is the frequency in Fourier space, and Γ the Green operator associated with ℂ . This equation is more convenient in Fourier space where the convolution product becomes a direct product, hence the use of Fast Fourier Transforms (FFT) introduced by [9] to solve the implicit integral equation (4) and obtain the strain and stress fields solution of the local problem.
Different algorithms have been introduced to solve this equation. All of them are based on the idea of [9] and their 'Basic Scheme' where the unit-cell is discretized into pixels (in 2D) or voxels (in 3D) and the problem equations are also discretized into finite sequences of space in order to compute discrete Fourier transforms (DFT) using FFT algorithms. The LS equation is then solved by means of a fixed point algorithm in the Basic Scheme, but this method showed some limitations such as a slow convergence and a non-convergence when having high stiffness contrasts between phases [10], [11]. New FFT-based homogenization methods were introduced after the original one, using new solving techniques based on Krylov subspace solvers for example [12], [15].

Nonlinear algorithm
In our code we chose to adopt the resolution scheme introduced by [16] instead of the basic scheme from [9]. [16] chose to formulate the LS equation as a residual: (5) and to solve it with a general Newton-Raphson (NR) algorithm. The time dependency is introduced here as local non-linear constitutive laws can be used and local stresses and internal variables are function of the previous states. The linearized problem corresponding to the residual (5) in the NR algorithm is the following: where ℂ * is the 'tangent stiffness matrix': The expression of the Jacobian matrix / ( ) is not straightforward as it involves a convolution product in real space, but again the equation (7) is easily computed in Fourier space with the use of FFT. As it is not straightforward, one cannot use a simple matrix inversion algorithm and the linearized problem of the NR algorithm is solved by means of a Conjugate Gradient (CG) method. The solution given by the CG is then used to update the strains of a NR iteration step.
The stress equilibrium condition is chosen as the convergence criterion of the NR algorithm [10], [16]: We used the CG-based solver from the PETSc library [17] in our homogenization code.
Different factors like the gas hydrate volume fraction, or their complex geometries have to be incorporated in the modeling of the mechanical behavior of GHBS. The nonlinear behavior of each microscopic constituent seems also to play a part in the overall mechanical response of the GHBS. Modeling the behavior of GHBS through macroscopic constitutive laws require a significant number of material parameters. Instead, FFT-based homogenization techniques provide a simpler way to integrate these nonlinearities through realistic microstructures and simple elastoplastic laws at the microscopic level.
We will illustrate it in the following section through the application of our code to a real image of GHBS with finegrained sediments and gas hydrate veins using non-linear constitutive laws at the microscale.

Image processing
A 256 pixels resolution unit-cell was taken from micro-CT scans of a depressurized and cooled core sample at the Krishna-Godavari Basin [14]. Figure 1 represents the vertical slice of the scans from which we extracted our image, and it corresponds to a sample area where the clayey matrix is crossed over by multiple veins of gas hydrates.
Since the extracted image is a non-periodic unit-cell, performing a numerical periodic homogenization over it would result in disturbed mechanical fields at the edges of this geometry. Instead we chose to modify the image prior to calculations to make it periodic, using Brisard's Python implementation [18] of the periodic-plus-smooth decomposition algorithm [19].
This image processing is illustrated in Fig. 1 and, even though it gives us a periodic unit-cell to work with, it creates a new image which is also disturbed compared to the original. The choice between working with either a nonperiodic image and disturbed mechanical fields or a modified periodic image with undisturbed solution fields depends on the author. Here, we chose to work on an example of application of our code, so we needed a real image, but we could allow ourselves a slight modification of the image.
Here, the periodic microstructure consists of two constituting phases: the fine-grained soil matrix and the gas hydrates forming veins (Fig. 1). We had to segment the image and for this we decided to test two threshold values while looking at the image histogram: (A) the minimum value of the valley between the two principal peaks of the histogram and (B) the average value of these two peaks (Fig. 2).
The two segmented images corresponding to the (A) and (B) threshold values are represented in Fig. 3. They are characterized by gas hydrate volume fractions of 16.59 % and 23.63 % respectively so we can see that the segmentation process has also an impact on the geometry of the unit-cell, which certainly will have an impact on the homogenization simulations. Finally, we applied morphological operations, erosion and dilation, to these two images just to obtain new segmented images with artificial Fig. 1. Vertical slice of the micro-CT scans of natural core section 46-66~cm published by [14], and the periodic-plus-smooth decomposition of the extracted image we used for our application. gas hydrate volume fractions of 14 %, 16 %, 23 % and 34 % (Fig. 3).

Calculations
We describe here the main parameters used for the homogenization calculations we ran with our code. The unit-cell was subjected to a mixed macroscopic loading, where the horizontal stresses were kept constant and a vertical strain was imposed: Δ = 5% and ΔΣ = ΔΣ = ΔΣ = 0 ∀ ≠ (10) The calculations started with an initial mean stress of 1.25 MPa over the unit-cell, and the tolerance for the NR convergence criterion was taken equal to 10 -5 .
Concerning the material parameters of each constituting phase, we chose an elastic perfectly plastic model with a von Mises' criterion for the gas hydrate phases, and a modified Cam Clay model for the soil matrix. All the parameters required to define the constitutive law of these two material phases are summarized in Tab. 1.

Results and discussion
The macroscopic stress-strain responses of each unit-cells for different and from different segmentation calculations (A or B) are represented in Fig. 4. The elastoplastic behavior of the sedimentary matrix has clearly an effect on the macroscopic response of the unit-cell. This macroscopic response is non-linear and can be divided in two parts during the loading: before and after all the pixels belonging to the soil matrix phase have reached the plastic regime.
There is a clear distinction between the mechanical responses of the two geometries directly obtained from the image segmentation process using either the threshold values (A) or (B) (the ones represented with a dashed frame around it in Fig. 3). It appears though that the choice of the threshold value for the image segmentation has little effect on the response when comparing results for the same gas hydrate volume fraction . Figure 5 gives the fraction of yielding pixels for each constituting phase of the unit-cells obtained with the   threshold value (A), as a function of the deviator stress q during shearing. We can see that the soil matrix yields more rapidly than the gas hydrate phase. However, when 100 % of the pixels of the gas hydrate phase have reached the constant deviator stress characterizing the perfectly plastic regime of this material component, it is not visible on the macroscopic response shown by Fig. 4. The macroscopic hardening response of the soil matrix prevails over the gas hydrate phase as they represent 76 % of the total number of pixels of the unit-cells.

Conclusions
We used a FFT-based homogenization method to simulate the macroscopic response of unit-cells representing finegrained sediments with complex gas hydrate morphologies. We based the geometries of these unit-cells on a real image extracted from micro-CT scans from the Krishna-Godavari Basin. This allowed us to study the impact of the image segmentation when having two constituting phases for the microstructure. This example showed the potential of the FFT-based homogenization methods applied to GHBS, in particular to fine-grained sediments. This could lead to the development of macroscopic constitutive laws for GHBS using complex microstructures without the experimental difficulties.