Multigrid method for nonlinear poroelasticity equations

In this study, a nonlinear multigrid method is applied for solving the system of incompressible poroelasticity equations considering nonlinear hydraulic conductivity. For the unsteady problem, an additional artiﬁcial term is utilized to stabilize the solutions when the equations are discretized on collocated grids. We employ two nonlinear multigrid methods, i.e. the “full approximation scheme” and “Newton multigrid” for solving the corresponding system of equations arising after discretization. For the steady case, both homogeneous and heterogeneous cases are solved and two different smoothers are examined to search for an efﬁcient multigrid method. Numerical results show a good convergence performance for all the strategies.


Introduction
Shale gas [1,2] is natural gas which is formed by being trapped within shale layer formations. Shale layers have typically low hydraulic conductivity which dramatically reduces the mobility of this so-called unconventional gas. Hydraulic fracturing [3] has been regarded as one of the methods of extracting these gas resources. It is a process in which injection of a highly pressurized fluid creates fractures within rock layers.
The concept of poroelasticity [4][5][6] can be used as a model describing the Earth. It is a well-developed theory that was originally studied by Terzaghi [7] who proposed a model for a one-dimensional consolidation problem. After that, in 1941, Biot [8,9] extended the theory to a general three-dimensional model. The model has been widely used for problems in rock mechanics, and it describes the interaction between the solid (rock) deformation and the fluid motion. It is a coupled model considering the seepage and stress processes together. The 2D poroelasticity problem can be formulated as a system of partial differential equations for the unknowns displacements and pore pressure of the fluid. Existence and uniqueness of the solution of this problem have been studied for example by Showalter [10] and Barucq et al. [11].
Often porous material is assumed to be homogeneous in numerical experiments, however, in fact materials usually have complicated properties composed of different characteristics. Therefore, it is necessary to take heterogeneity into account which can influence the poroelastic behavior in many 123 ways. Heterogeneity means that the coefficients in the equations are not constant and all (or several) characteristics of the main problem follow some distribution. Also the hydraulic conductivity of the material plays a role. There is a significant difference in conductivity once material deformations start to occur. The coefficient of conductivity depends on the stress and fluid pressure, resulting in a nonlinear set of equations. Both heterogeneity and nonlinear conductivity are included in the nonlinear poroelastic model studied here. As an analytic solution is usually not available, we solve the poroelasticity system by means of numerical techniques. We will employ the finite volume method for the nonlinear system of poroelasticity equations. Details about the convergence results of the multi-dimensional finite volume discretization for the nonlinear system of poroelasticity equations is, to our knowledge, not yet known or available in the literature. In one dimension, however, convergence of the discrete solution has been shown in [12] for a staggered arrangement of the poroelasticity unknowns in the nonlinear case. Convergence results of discrete schemes for the (multidimensional) linear system of poroelasticity equations, on a staggered and on a collocated grid, are available, see, for example, [13,14].
We would like to employ the multigrid method [15][16][17][18][19] as the iterative solution method for the discretized partial differential equations. It is known that many classical iterative methods, like Jacobi or Gauss-Seidel cannot efficiently eliminate low-frequency errors appearing in a numerical approximation. This may cause slow convergence, as low frequency error modes tend to disappear extremely slowly from a discrete approximation. Therefore, the multigrid method is chosen as a highly efficient solution method, improving the performance of the basic relaxation schemes. Regarding the nonlinear system, the multigrid iterative method can also be employed as a nonlinear solver. The "full approximation scheme" (FAS) and "Newton multigrid" are both considered for the time-dependent problem. With respect to the smoother in the multigrid algorithm, coupled smoothers called Vanka and point Gauss-Seidel (PGS) are chosen and compared. Up to our knowledge, results about the convergence of the nonlinear multigrid method are unknown.
The organization of our work is as follows: First of all, we present the governing equations of the unsteady and steady poroelastic model, together with the finite volume discretization scheme on a collocated grid in Sect. 2. In Sect. 3, the nonlinear multigrid methods are introduced. Each component of multigrid is clarified. After that, numerical experiments are presented in Sect. 4. All of the results show satisfactory convergence performance of the proposed methods. Finally, a conclusion is given in Sect. 5.

Governing equations
We deal with a deformable fluid-saturated porous medium, whose solid matrix is elastic and the fluid is viscous. Both solid matrix and pore network are considered to be continuous, and thus fully connected. Biot's poroelastic theory [8,9] is based on the coupling between the coherent solid skeleton and the pore fluid flow. A change in the applied stress of the skeleton will affect the pressure or mass of the fluid, and a change in fluid pressure will lead to a change in the volume of the porous material.
A poroelasticity system is constructed on the description of the fluid pressure, stress, displacement and strain in the medium, and mass and momentum conservation principles. Supposing that the porous matrix is fluid-filled, the total Cauchy stress tensor σ i j can be divided into two parts, pore pressure (fluid) p and effective stress of the soil skeleton (solid) σ i j . The effective stress is defined as a subtraction of pore pressure from the total stress. Pore pressure only influences the normal stress. The total momentum balance reads: where F i is the body force in the i-th direction. The summation convention is used when repeated subscripts occur. The strain quantity, ε i j , is a measure of the solid deformation with respect to an initial state. Variables u i and ε i j are related according to the compatibility condition: Moreover, ε v = ε 11 + ε 22 + ε 33 , is the volume strain. The constitutive equation in Biot's model is based on the assumptions of linearity between stress and strain: where λ and G are the effective Lamé constant and the effective shear modulus of the porous matrix, respectively; δ i j is the Kronecker delta; α denotes the coefficient of pore pressure, which is also called effective stress coefficient.
Darcy's law [20] describes the rate at which a fluid flows through a permeable medium. The seepage equation in the poroelasticity model is obtained when substituting Darcy's law into the continuity equation for the fluid, i.e., with k the hydraulic conductivity, and γ is the product of the porosity and the compressibility of the fluid. When an element of rock undergoes elastic deformation, its hydraulic properties will change. Hydraulic conductivity regarding as a transport property describes how easily fluid flows through the rock material which will be influenced by variation of stress. Supposing that the conductivity and stress follow a negative exponential function, the conductivity reads where k 0 is the initial conductivity of the rock element; β is a coupling coefficient which reflects the influence of stress on the coefficient of conductivity; σ ii /3 is the average stress; and ξ is a mutation coefficient to account for the increase in conductivity of the material during fracture formation. When ξ > 1, the above expression gives a higher conductivity caused by damage. In other words, the conductivity in Eq. (5) will dramatically increase when there is a "failure" in the element.
For the unsteady-state case in 2D, the governing equations have a time dependent term and a Navier-type equation for displacement vector u resulting by applying the constitutive Eq. (3) and the geometric relation (2) where ε i j is expressed in terms of the displacement gradient to the balance Eq. (1), The source terms g = (g 1 , g 2 ) and f are supposed to be in (L 2 ( )) 2 and L 2 ( ), respectively. They are used to represent a density of applied body forces and a forced fluid extraction or injection process for each case [13]. The term (α grad p) presents the stress due to fluid pressure within the medium, and div u is the volume change rate. Here μ is the Lamé coefficient that equals to the shear modulus G.
Coefficients λ and μ are related to Young's modulus and Poisson's ratio by

Discretization
When discretizing the poroelasticity equations, we employ a collocated grid [13,18,21] on which all variablesdisplacements vector u = (u, v) and fluid pressure p, are placed at the same grid points. Collocated grid arrangements are convenient for numerical iteration methods like multigrid. Whereas for unsteady poroelasticity simulations an artificial stabilization has to be used, this is not the case for the steady poroelasticity case, as oscillations of pressure do not occur there.
The finite volume method is employed as the discretization scheme. We take the first equation in (6) as an example. It can be written as where σ denotes the two-dimensional stress matrix σ x x σ xy σ yx σ yy . In Fig. 1, the square with dotted lines represents a control volume V i, j used for collocated grids, and i, j denotes the boundary of V i, j . Now, we consider the first equation in (8) for example which can be expressed as, By integration of (9) over an arbitrary control volume V i, j , one obtains, The Gauss divergence theorem converts volume integrals to surface integrals, resulting here in where n is the unit normal vector to the volume, and τ i (i = 1, 2, 3, 4) is the right, top, left and bottom boundary of the control volume V i, j , respectively. By applying to (12) the following relations, the resulting discrete equation is given by where h is the grid size of the uniform grid used for the space discretization. Central differences are applied to the firstorder derivatives u x , u y , v x , and v y . Similarly, the discrete equation for the second equation in (8) can also be obtained. It is convenient to apply stress boundary conditions in finite volume Eq. (12), by substituting the given stress at boundaries, of course, adapting the control volume.
Regarding the seepage equation, we discretize it on a collocated grid as well. However, such discretization for unsteady problems may be unstable, because oscillations may appear in the first time steps of a numerical solution. After this phase, the solution becomes smoother and these oscillations tend to disappear. Some special care is needed to construct a stable discretization for the whole process. The stabilized discretization can be achieved by adding an artificial elliptic pressure term [13,21] to the seepage equation in (6), as follows The artificial term is ∂/∂t (∇ · ( ∇ p)), with ε = h 2 4(λ+2μ) , see [13]. When the grid size h → 0, the artificial term tends to 0. Now the system of equations with a stabilization term added, can be discretized on a collocated grid with the finite volume method. Since this term is proportional to h 2 , second order accuracy can be maintained.
Since Eq. (15) is a time-dependent problem, a time discretization is required. The discrete problems have to be solved in each time step. In each of the finite volumes V i, j , the integral formulation of (15) has the form, Using the Gaussian theorem, we have where t is the step-size in time direction and the approximation of conductivity k i+ 1 2 , j can be expressed as For other values of permeabilities, a similar formulation can be obtained without difficulties. Obviously, the values of unknowns {u n+1 , p n+1 } at a new time t + t can be calculated immediately from the values of the previous time. From discrete formula (17), the backward Euler and Crank-Nicolson schemes can both be used by choosing with θ = 1 and θ = 0.5, respectively. We prefer to use the Crank-Nicolson scheme so that second order accuracy in time can be obtained.

Governing equations
Under steady-state conditions, (4) becomes a Poisson-like equation as From Eqs. (6) and (19), the governing equations for the 2D steady poroelasticity model can be obtained as Here we simply consider the case that the source terms are all zeros.

Discretization
Since the seepage equation in (20) is a Poisson type equation, second-order accuracy can be achieved again by applying finite volumes on a uniform rectangular grid. Note that the artificial stabilization term is not needed for the steady case. The discretization form of the seepage equation reads The left-hand side of (21) can be reformulated to an integral over the boundary of the volume V i, j , Replacing p x and p y by a centered approximation, one obtains It can be seen that Lamé coefficients and permeabilities are required at the four mid points (i + 1 2 , j), (i − 1 2 , j), (i, j + 1 2 ) and (i, j − 1 2 ) that are at the control volume boundaries. If the material is heterogeneous, these parameters are determined from a stochastic distribution, which means they are different at each vertex of the collocated grid. Due to this, the averaged values of two adjacent vertices are used for those middle points in the discrete formulas. Boundary conditions will be specified in Sect. 4.

Numerical method
The multigrid method-an efficient numerical technique for solving systems of linear and nonlinear equations-is employed for the solution of the discretized poroelastic partial differential equations, based on earlier multigrid work on poroelasticity model problems [13,21].

Nonlinear multigrid method
To deal with the nonlinearity, there are basically two approaches.
Newton multigrid there is no doubt that Newton's method is the most important method for solving nonlinear equations. Newton multigrid is based on global linearization. Newton's method is applied to linearize the equations, and multigrid solves the resulting linear Jacobian system.
Full approximation scheme the other multigrid technique suitable for nonlinear problems is the FAS [15,18] which treats directly the nonlinear equations on fine and coarse grids. In FAS, a nonlinear iteration, such as the nonlinear Gauss-Seidel method is applied to smooth the error. Differently from linear multigrid, the full-scale equation is solved on the coarse grid instead of the residual equation, because of the nonlinearity.

Multigrid components
To apply multigrid efficiently, each component of multigrid needs to be selected with special attention. Those components are related to the collocated grid arrangement of the poroelasticity discretization chosen here. The multigrid transfer operators, restrictions and prolongations, are then well-known in geometric multigrid. Regarding the most important component in this setting, i.e., the smoothing scheme, we consider two methods for the time-dependent system of equations-Gauss-Seidel relaxation and boxrelaxation [13,21].
Smoothers the choice of smoother is the most significant part which can crucially affect the performance of multigrid. Several branches of robust smoothers have been developed for the poroelasticity equations. They all fall into two major categories: decoupled and coupled smoothers. With respect to the decoupled smoother, which is also denoted as equation-wise relaxation, distributive Gauss-Seidel (DGS) is the original method, introduced in [22]. Here we focus on the other type of smoothers-coupled relaxation. A coupled smoother is a state-of-the-art relaxation scheme for saddle point problems. The technique is based on processing each grid cell in some order and to relax all unknowns associated with that cell simultaneously.
A point-wise collective Gauss-Seidel (PGS) relaxation is chosen, which processes three unknowns u i, j , v i, j and p i, j at grid point (i, j) simultaneously. A small 3 × 3 system is solved for each grid point. We consider the correction equation during smoothing for convenience, i.e., ⎛ ⎝ a 1,1 a 1,2 a 1,3 a 2,1 a 2,2 a 2,3 a 3,1 a 3,2 a 3,3 where eu m+1 i, j = u m+1 i, j − u m i, j is an increment to the solution u corresponding to node (i, j). It is the same for v and p. ru i, j , r v i, j and r p i, j represent the corresponding equation of the system. After solving the residual equations, the computed increments will be added to the current solution (take u for example), with an underrelaxation parameter ω ranging from 0 to 1. The use of a collective smoother in the case of nonlinear systems of equations requires some more explanation. In principle, we have a variety of choices of local linearization and grid point processing at our disposal, like Newton-Gauss-Seidel, Newton-Jacobi, Picard-Gauss-Seidel, Picard-Jacobi relaxation, and many others, see, for example [17,18].
We here employ a "straightforward" basic linearization variant to deal with the nonlinearity locally.
When we process a certain grid point, in the collocated grid, we assume that we work in a lexicographical Gauss-Seidel fashion, and that we have already updated the poroelasticity unknowns on previous grid points. Because of the form of the nonlinearity in (20), we will use the most recent updated unknowns when we set up the matrix element related to the nonlinearity, in the system (24).
Due to the nonlinear term, k(σ, p) = ξ k 0 e −β((λ+μ)div u−αp) , we need to define div u at a certain grid point (i, j), and we will always take the latest updates of the displacement unknowns to define this matrix element. This means, the just updated neighbouring u-values are combined with the "old" displacement values for the points that have not yet been processed. For the pressure in the nonlinear term, we always take the "old" value, in the linearization process. This is a pragmatic way of linearization, related to Gauss-Seidel-Picard linearization. We will perform the same local linearization when dealing with the Vanka smoother in (26), and carefully consider each equation, take the latest values for all unknowns. These are sometimes recently updated values, and sometimes the values from a previous iteration.
The second smoothing scheme considered is the Vanka smoother which is also called box relaxation. It was originally proposed by Vanka [23] for solving the Navier-Stokes equations discretized by the finite difference method. This approach can be applied in a wide range of problems, which we extend to the poroelasticity equations, see also [13]. Instead of the three unknowns at each grid point as  Fig. 2). At each grid-point, a (5 × 5)-system is solved.
Coarse grid correction with respect to the coarse grid correction, we choose geometric grid coarsening as we will deal with regular Cartesian grids in these experiments. The sequence of coarse grids is obtained by doubling the mesh size in each spatial direction and we use a (2 × 2)-grid as the coarsest grid. As in the scalar case, any suitable solver can be applied on the coarsest level.
Transfer operators supposing we have performed several smoothing steps to obtain an updated solution and a sequence of coarse grids is ready, the next step is to restrict the approximation to the next coarser grid. Since we are trying to solve a system of equations, there are more than one equation at each grid point. The restriction is done separately for each of these equations in the system. As a standard choice for the scalar problem, the full weighting (FW) operator is chosen for the restriction in our case. With respect to the prolongation, a typical choice is bilinear interpolation for each unknown grid function from coarse grid to the next finer grid.

Unsteady case
We present some comparison results for the usteady case for a model problem.
Initial settings for the unsteady problem, we first consider the case of a homogeneous, isotropic and incompressible medium. The computational domain is (0, 1) × (0, 1). We numerically solve the poroelasticity problem with a simple analytic solution given by u = cos(π x)sin(π y)sin(π t), v = sin(π x)cos(π y)sin(π t), p = −2(λ + 2μ)π sin(π x)sin(π y)sin(π t). (27) Source term g and f are consequently determined from (6). Dirichlet boundary conditions are considered. Before fluid starts to flow, the following initial condition is considered, The nonlinear conductivity k(σ, p) = e (λ+μ)·div u− p is used in the computations. We discretize the incompressible and time-dependent poroelasticity equations with an artificial pressure term by finite volume methods and Crank-Nicolson scheme in time. The PGS smoother and box-relaxation, i.e. Vanka smoother are both taken into consideration. The Lamé coefficients are taken as λ = μ = 0.1 and the final time is set as t = 0.5.
FAS firstly, FAS is employed to solve the time-dependent problem. It is known that in FAS the full equation is solved on the coarse grids. It is required to transfer the current approximation from the fine grid to the coarse grid. After this process, the coarse grid error is subtracted from the solution. The correction is then interpolated and added to the fine grid approximation [24].
The underrelaxation parameters for the two smoothers, i.e. PGS and Vanka, are 1.0 and 0.7, respectively. The results are generated by an F-cycle with two pre-and two postsmoothing steps. The stopping criterion is that the maximum norm of the residual r u ∞ + r v ∞ + r p ∞ should be less than 10 −7 . Table 1 presents the relative errors between analytic and numerical (with subscript h) solutions in the maximum norm at final time t = 0.5 with the artificial pressure term for different target grids. With the decrease of the mesh size, the relative error is one quarter of the previous one. The CPU time and multigrid convergence factors are also shown in Table 1. The convergence factor represents an average residual reduction factor over previous iterations. From Table 1, we can conclude that second order accuracy is maintained and FAS performs very well.
When the nonlinear conductivity is an extremely small value, i.e. k(σ, p) = 10 −13 · e (λ+μ)·div u− p , the convergence performance of Vanka is much better than PGS, see Table 2. With very fine grids and a small time step, we will have a saddle point problem. PGS doesn't work, which results in a divergent solution. Vanka is still efficient for this kind of problem.
Newton multigrid different from the idea of FAS, Newton multigrid is employed too. A standard Newton method is applied to linearize the equations, then multigrid follows for the solution of the (linear) Jacobian system in each iteration. It is a combination of Newton's method for the outer iteration and multigrid for the inner iteration.
In this test, only one F-cycle is used to solve the Jacobian system. The convergence factors in Table 3 are the corresponding to the Newton method. A comparison between Table 1 and Table 3 indicates that convergence of FAS is 123  a bit faster than that of Newton multigrid. In general, the convergence performance of both methods is very satisfactory.
Regarding the small value of the nonlinear conductivity, we can get the same conclusion. Vanka still works fine even with fine grids and a small time step. However, PGS fails to get a convergent solution, see Table 4. Here we still applied one multigrid cycle for the Jacobian system. If more multigrid cycles are used, PGS is available for some easy cases that the grid and time step are not very small.

Steady case
Due to the efficiency and robustness of FAS for the timedependent nonlinear problems, now we only consider FAS as the solver for the steady numerical tests. In this section, we will also consider poroelasticity systems with heterogeneous coefficients, which is closer to the real engineering applications.

Homogeneous test
First, a numerical test with homogeneous material is chosen for the steady problem (20). The mechanical parameters are constants at each grid point.
System (20) is solved iteratively by multigrid with both the PGS smoother and the Vanka smoother on different grid sizes. Table 5 shows the multigrid convergence results by using an F-cycle. F(3, 3) means three pre-and three postsmoothing steps. The numbers in the table denote multigrid convergence factors, and the corresponding CPU time in seconds is presented in brackets. In general, the performance of both smoothing schemes is satisfactory, resulting in hindependent convergence of the multigrid method. PGS takes less CPU time than the box relaxation method, as expected.

Heterogeneous test
Weibull distribution usually, a material has complicated properties and boundaries. Rock mass is assumed to be a heterogeneous deformable body composed of different material properties. Heterogeneity is a concept relating to nonuniformity (composition or character) in a substance. Rock specimen in numerical tests is subdivided into square elements with randomly distributed mechanical properties in each element. Due to the size and shape consistency, there is no geometric preference orientation in the specimen. In order to simulate random heterogeneity in a rock, a statistical approach is used. In [25] the material heterogeneity is defined by a Weibull distribution, with the probability density function given by where s denotes a given mechanical property, such as Young's modulus, the coefficient of conductivity or the strength;s is the mean value; and m represents the homogeneity index which defines the shape of the distribution function. The corresponding Weibull distribution function is given by Figure 4 displays, for different values of the homogeneity index m, the probability density function in terms of the ratio between s ands. It is obvious that a higher value of m represents more homogeneous material and vice versa, as for higher m, the values of s are concentrated arounds. As an example, we consider initial Young's modulus E 0 with a mean valueĒ 0 of 50 GPa and homogeneity index m = 3. In this case, Fig. 5 shows a possible randomly distributed Young's modulus E 0 in each element.
In our model, the Young's modulus and initial conductivity are modeled in this way. So, these parameters differ for each element. Figure 6 presents histograms of the observed Young's modulus in three simulations with different homogeneity indices, m = 1, 3, 5. The chosen average is the same as in Fig. 5. Notice that the randomly distributed parameters on coarser grids are generated from finer grids. In this way, the coarse system is related to the fine system, thus the same characteristic of the material is guaranteed. The standard restriction operator "Full Weighting" is used to transfer these parameters to each level of the grids. Then we perform a regular finite volume method discretization with these averaged coefficients on the coarse grids.
Some results in this test we solve the steady equations for a heterogeneous material. Two different homogeneity indices m are chosen to investigate the convergence results of multigrid. Figure 7 shows the randomly distributed Young's modulusĒ 0 = 50 GPa with different homogeneity indices m = 1 and m = 5. Obviously, m = 5 denotes a more homo- geneous material compared to m = 1. Conductivity also follows the stochastic distribution and the average value of k 0 is taken as 0.01 m 2 . All the other parameters, i.e. Poisson ration ν and coupling coefficient β, are the same as for the homogeneous test.
We consider a challenge here to deal with these heterogeneous cases using the standard multigrid components described above. So, we do not use the Galerkin coarse matrices or operator-dependent prolongation and restriction here. We merely wish to study the impact of the heterogeneity on the FAS multigrid convergence using the standard geometric multigrid components, like the direct discretization of the operator on the coarse grids, and the two described smoothers. Tables 6 and 7 show the multigrid convergence results for the heterogeneous tests with m = 1 and m = 5, respectively. It can be seen that the convergence factors are larger than the results from the homogeneous test case, as expected. Comparing the results in Table 6 with those in Table 7, the multigrid convergence is better when the distributions of the mechanical parameters are more homogeneous. PGS is still faster than the Vanka method for this steady problem in CPU time. Overall, the multigrid convergence is still highly satisfactory for the heterogeneous test cases used here with the standard multigrid components. Convergence can be further improved in this framework by using multigrid as a preconditioner for a Krylov subspace iteration, like in the setting of the recombination technique in [26,27].

Conclusions
In this study, we have used the nonlinear multigrid method to solve the poroelasticity system considering material heterogeneity and a nonlinear conductivity. First, we have solved an unsteady problem. Since oscillations may occur in the first time steps of the solution when the equations are discretized on a collocated grid by the finite volume method, some special care is needed in discretization. The stabilization can be achieved by adding an artificial term in one of the continuous equations. A simple analytic test is employed to verify our method. This is done by FAS and Newton multigrid with a coupled smoother "Vanka" which solves the discrete equations locally cell by cell, and the PGS method. Numerical results show that both nonlinear strategies are stable and second order accurate. Vanka is more efficient and robust for some difficult problems where the values of the coefficient are extremely small or the system is discretized on a very fine grid.
Next to these tests, we consider steady problems for homogeneous and heterogeneous cases. Heterogeneity means that some characteristics of the main problem follow a random distribution. The Full Approximation Scheme with collective point-wise Gauss-Seidel smoother, that is, updating all unknowns at each grid point together, shows highly satisfactory multigrid convergence. For the tests performed, we do not need the commonly used box relaxation scheme as smoother.