A Generalized Multiscale Finite Element Method for Poroelasticity Problems II: Nonlinear Coupling ^{†}^{†}thanks: This work was supported by RFBR (project N 153120856)
Abstract
In this paper, we consider the numerical solution of some nonlinear poroelasticity problems that are of Biot type and develop a general algorithm for solving nonlinear coupled systems. We discuss the difficulties associated with flow and mechanics in heterogenous media with nonlinear coupling. The central issue being how to handle the nonlinearities and the multiscale scale nature of the media. To compute an efficient numerical solution we develop and implement a Generalized Multiscale Finite Element Method (GMsFEM) that solves nonlinear problems on a coarse grid by constructing local multiscale basis functions and treating part of the nonlinearity locally as a parametric value. After linearization with a Picard Iteration, the procedure begins with construction of multiscale bases for both displacement and pressure in each coarse block by treating the staggered nonlinearity as a parametric value. Using a snapshot space and local spectral problems, we construct an offline basis of reduced dimension. From here an online, parametric dependent, space is constructed. Finally, after multiplying by a multiscale partitions of unity, the multiscale basis is constructed and the coarse grid problem then can be solved for arbitrary forcing and boundary conditions. We implement this algorithm on a geometry with a linear and nonlinear pressure dependent permeability field and compute error between the multiscale solution with the finescale solutions.
1 Introduction
The applications of mechanics and flow in porous media are wide ranging, as are the challenges involved in simulating some of these problems in nonlinear multiscale contexts. This is particularly true in geomechanical modeling where relevant phenomena may be highly nonlinear, for example in the setting of enhanced production and environmental safety concerns due to overburden subsidence and compaction [19, 20]. Another of the central challenges is the multiscale nature of the media considered in geomechanics problems. Heterogeneity of rock properties should be accurately accounted in the geomechanical model, and this requires a computationally costly a high resolution solve. Moreover, due to the multiphysics nature of the problems, they may involve highly nonlinear relations. This then makes the further requirement of many iterations in a Newton or Picard linearization. Thus, we propose a multiscale method to attempt overcome some of these challenges. The central idea is to linearize in a Picard iteration, and treat the nonlinearities as a parametric value as utilized in [10] and references therein.
As noted in [9], the basic mathematical structure of the poroelasticity models are usually coupled equations for pressure and displacements known as Biot type models [24]. The pressure equations are a parabolic equation coupled to a time derivative of volumetric strain. While the mechanics equations are are given by quasistatic elasticity equations and is coupled by gradients of pressure. In this work however, we focus on the possible nonlinear couplings of the Biot model. There are a myriad of physical and modeling reasons to add nonlinearity to the Biot equations, however, we will primarily focus when the permeability field and elasticity tensors depend nonlinearly on pressure and displacements and their gradients. This is due primarily to us wanting to focus on the nonlinearities effects on our GMsFEM, as nonlinearities in lower order derivative will not interfere with the construction of the local multiscale basis functions.
Nonlinear Poroelastic models of this type have been explored in the literature to incorporate higher order physics considerations. For example, when the viscosity of the fluid heavily depends on the fluid pressure we may obtain relations of permeability of the form
Here is the absolute permeability, and is the pressure and spatially dependent viscosity. This can occur when their are very high pressure gradients [29]. In the setting of complex geomechanical interactions [15] used a relationship between permeability and volumetric strain of the form
where are determined constants and is the volumetric strain. Further in [15], the porosity also depends linearly on however, this is multiplied throughout generating a nonlinearity. In the context of fractured reservoirs, permeability is often computed via the so called ”cubiclaw” through channels and this may be coupled in orientation and magnitude via the displacements in a nonlinear way [30]. With this GMsFEM, we propose a method to efficiently compute solutions to these nonlinear poroelasticity problems with the heterogeneous multiscale properties.
As noted in the prequel [9], there are many very effective multiscale frameworks that have been developed in recent years. There are rigorous approaches based on homogenization of partial differential equations [7, 8]. However, these approaches may have limited computational use. Examples of computational approaches include the Heterogeneous Multiscale Method (HMM) [25, 26], approaches based on the Variational Multiscale Method (see [27]), where coarsegrid quasiinterpolation operators are used to build an orthogonal splitting into a multiscale space and a finescale space [28]. In this paper, we will use the Generalized Multiscale Finite Element Method and its extension to nonlinear poroelasticity problems in the framework of [10]. Specifically, to handle multiscale nonlinear problems, we combine ideas of model reduction, whereby the nonlinearity is replaced locally by a parameter space and offline and online spaces are generated. For a broad presentation of these methods we refer the reader to [1].
The paper is organized similarly to [9], as follows. In Section 2 we provide the mathematical background of the nonlinear poroelasticity problem. We will introduce the Biot type model and highlight where the heterogeneities primarily occur. We again note that the nonlinearities in our model are in the permeability and elasticity tensor as these are second order derivative terms. In Section 3, to outline the difficulties in full direct numerical simulation we introduce the finescale discretizations using coupled timestepping schemes and a Picard iteration technique for linearization. In Section 4, we present our nonlinear GMsFEM algorithm and outline its construction procedure. Finally, numerical implementations are presented in Section 5. Using the GMsFEM, we compare the multiscale solution to finescale solutions and give error estimates. We will present two different examples with permeability being linear and nonlinear with respect to pressure. Additionally, we will implement and discuss different snapshot spaces and coarsegrids choices, and its relation to enrichment and the error.
2 Problem formulation
We denote our computational domain to be a bounded Lipschitz region. We consider a general nonlinear poroelastic system where we wish to find a pressure and displacements satisfying
(1a)  
(1b) 
with initial condition for pressure We write the boundary of the domain into four sections . We suppose the following boundary conditions on each portion
and
Here the symmetric strain is written as and we write
to mean the double contraction of a tensor with a tensor.
As in the linear case, the primary sources of the heterogeneities in the physical properties arise from , the elastic tensor, and , the permeability. In this setting, we suppose these heterogenous parameters can depend in and and their gradients in complicated nonlinear ways. Further, we will denote to be the Biot modulus, is the fluid viscosity, and is the BiotWillis fluidsolid coupling coefficient. Here, is a source term representing injection or production processes and is the unit normal to the boundary. Body forces, such as gravity, are neglected.
Remark: Note that one could also add nonlinearities in the coefficients and , however, these correspond to lower order terms with respect to derivatives. Therefore, these will not contribute to the local problems in the GMsFEM. Hence, we will consider them to be constant throughout.
We recall the setting when these relations become linear. In the case of a linear elastic stressstrain constitutive relation we have that the stress tensor and symmetric strain gradient may be expressed as
where , are Lame coefficients, is the identity tensor. Note here this is not to be confused with what is often used as a parameter. Above we have absorbed into the nonlinear permeability coefficient the fluid viscosity , and in the case of linear permeability, we have
being absolute permeability.
The nonlinear poroelasticity problem (1), can be written in operator matrix form:
(2)  
(3) 
where
and and are gradient and divergence operators and .
3 FineScale Discretization
We will now present finescale approximation and nonlinear solution methods for the above system. We will motivate the need for a multiscale method due to the nonlinearity and the heterogeneity of the poroelasticty problem. To approximate the solution to (1) on finescale grid we will utilize a standard finite element method. The corresponding nonlinear variational form of the continuous problem written as
(4)  
(5) 
for , where
and the test spaces with homogeneous boundary conditions are given by
We define the following nonlinear forms
(6a)  
(6b) 
and bilinear and linear forms
and
Here under the integrand denotes the standard inner product. In Section 5, we will discretize the spaces using a finescale standard FEM and denote them and , being the finegrid size. The FEM using these spaces will serve as a reference solution for our GMsFEM outlined in Section 4.
Nonlinear Solve: We will first consider the time discretizations of the above system, then will discuss resolving the nonlinearity. This discretization leads to several possible couplings between timesteps and the two equations of linear poroelasticity [13, 18]. However in the nonlinear case we will only consider the fully coupled scheme. We proceed by introducing for the nonlinear fully coupled time derivative operators and then the Picard iteration for the linearization of the nonlinear operators.
The standard fully implicit finitedifference scheme, or coupled scheme, can be used for the timediscretization and is given by
(7a)  
(7b) 
with , , where , , and .
We will now consider nonlinear solve in space after time discretization by the fully coupled scheme (7). One could rewrite (7) as a nonlinear system each time step and use a Newton solver, however, for our GMsFEM we prefer to use a linearization based on Picard iteration. Indeed, we may linearize (6) given from a previous iteration step we write
where . We choose this notation in part to emphasize this may be viewed as a parameter in offline phase of the GMsFEM.
Fixing the timestep at , and taking , as data from the previous iteration. For , we wish to find such that
(8a)  
(8b) 
Once the desired convergence criteria is reached, we can set the terminal as previous time data. We then return to the algorithm timestepping and continue the iterative linearization until the terminal time. Note that this process can also be used in an appropriate nonlinear generalization to a fixed stress splitting [13, 18].
4 GMsFEM for nonlinear poroelasticity problem
We will present the offline and online multiscale basis construction in the fluid or pressure solve then its construction in the mechanics or displacement calculation step in this nonlinearly coupled formulation. Similar to the presentation outlined in [9], however, we will focus on the effects of the nonlinearities on the method. Observing the linearized formulation (8), we see that we may consider the nonlinearity as parametric values we are able to successful design a GMsFEM for this nonlinear problem. In this way, we are able to construct an onlineoffline multiscale basis with respect to this nonlinearity. We now outline the general procedure of the GMsFEM algorithm.
We begin briefly with some standard notation. The overall finescale model equations will be solved on a finegrid using spaces and , and will be used for our reference solutions. We now introduce the coarse grid. Let be a standard conforming partition of the computational domain into finite elements. The finegrid, can be taken as a refinement of the coarsegrid. We refer to this partition as the coarsegrid and assume that each coarse element is partitioned into a connected union of fine grid blocks. We use , where is the number of coarse nodes, to denote the vertices of the coarse mesh , and define the neighborhood of the node by
See Figure 1 for an illustration of neighborhoods and elements subordinated to the coarse discretization. We emphasize that the use of is to denote a coarse neighborhood, and we use to denote a coarse element throughout the paper.
For global coupling we use the linearized continuous Galerkin (CG) formulation to find such that
(9a)  
(9b) 
where and denote the online spaces. The online spaces are spanned by multiscale basis functions and for time step and th iteration, each of which is supported in
The indexes represent the numbering of these multiscale basis functions for pressure and displacements, respectively. Here the parameter represents the nonlinear dependence as in (9). Recall that we may take , from the previous timestep and will treat these variables as parametric values on each coarse patch. However, for simplicity we will suppose that the dependence is only on of this nonlinearity.
Remark: Note, the derivative dependent problems, with nonlinear couplings of , may be handled. However, due to the oscillation in these quantities, these terms may not be well approximated by constants on the coarsegrid level. Thus, we would need to have a more enriched parameter space than is utilized here.
We now discuss further how we handle the parametrized nonlinearities. We assume that and are bounded above and below, i.e. and , where and are predefined constants. These may be guessed initially based on initial data or apriori estimates. The intervals and are divided into equal regions:
and
Clearly, if necessary these domains can be partitioned in different number of regions, but for simplicity we suppose they are equal in number. For the parameter we take average values of and in each coarse region . For average of a function we will use the notation
More specifically, we use to represent the dependence of the solution on . The multiscale basis functions will be computed for a selected number of the parameter values , at the offline stage and we will compute multiscale basis functions for each new value of for each at the online stage.
Boadly speaking, the GMsFEM algorithm consist of several steps:

Offline computations:

Generate the coarsegrid, .

Construct the snapshot space, used to compute an offline space, by solving many local problems on the finegrid.

Construct a small dimensional offline space by performing dimension reduction in the space of local snapshots.


Online computations:

In each time step and nonlinear iteration for current value of in each , we compute multiscale basis functions and construct online space by performing dimension reduction in the offline space.

Use small dimensional online space to find the solution of a coarsegrid problem for any force term and/or boundary condition.

We construct multisclate basis functions for pressure and displacements separately. We begin by considering the pressure solve, then, the displacement solve.
4.1 Multiscale basis functions for pressure
In the offline computation, we first construct a snapshot space . Construction of the snapshot space involves solving the local problem for various choices of input parameters and various boundary conditions. These local spatial fields are used then used construct the offline space and the space consists of fields defined on a fine grid. There are a few options available when constructing the snapshot space and we will proceed with the two most natural ways.
Snapshot Space 1: First, we propose a snapshot space generated by harmonic extensions of . For simplicity, we will omit the index when there is no ambiguity. We thus define such that
(10) 
Here are defined by , where denotes the finegrid boundary node on . This is done for each fixed parameter
Snapshot Space 2: Alternatively, we may use local finescale space basis functions within a coarse region and construct local snapshots by solving the following eigenvalue problem with natural boundary conditions
(11) 
Where
are the standard finescale basis functions, and for each fixed parameter values , .
Let be the number of functions in the snapshot space in the region , and define
for each coarse subdomain . We reorder the snapshot functions using a single index to create the matrix
where denotes the total number of functions to keep in the snapshot construction.
To construct the offline space , we perform a dimension reduction of the space of snapshots by using an auxiliary spectral decomposition. More precisely, we solve the eigenvalue problem in the space of snapshots:
(12) 
where
Here
is independent of and are prescribed nonnegative weights. The main objective is to use the offline space to accurately construct a set of multiscale basis functions for each in the online stage. At the offline stage the bilinear forms are chosen to be parameterindependent, such that there is no need to reconstruct the offline space for each .
We then choose the smallest eigenvalues from Eq. (12) and form the corresponding eigenvectors in the space of snapshots by setting , for , where are the coordinates of the vector . We denote the span of this reduced space as .
At the online stage, for a given parameter value , multiscale basis functions are computed based on each local coarse region . The associated online space is the small dimensional subspace of the offline space. In particular, we seek a subspace of the offline space that can approximate any element of the offline space in an appropriate sense. In the the online stage the bilinear forms are chosen to be parameterdependent and we use following eigenvalue problem
(13) 
where
Here and are the fine scale matrices corresponding to the stiffness and mass matrices for given and
Finally, we multiply the partition of unity functions by the eigenfunctions in the online space to construct the resulting basis functions
(14) 
where , is the standard linear partition of unity functions and the denotes the number of online eigenvectors that are chosen for each coarse node . We note that the construction in Eq. (14) yields continuous basis functions due to the multiplication of offline eigenvectors with the initial (continuous) partition of unity. Next, we define the online space as
(15) 
Using a single index notation, we may write , where denotes the total number of basis functions in the spaces and is number of coarse mesh nodes.
Denote the matrix
where are used to denote the nodal values of each basis function defined on the fine grid.
4.2 Multiscale basis functions for displacements
For construction of multiscale basis functions for displacements we use similar algorithm that we used for pressure. We first construct a snapshot space for each parameter . Again, as with pressure we give two possible snapshot space choices.
Snapshot Space 1: As our first possible snapshot space we propose the harmonic extension using . We define as the solution to
(16) 
Again, , and for each fixed parameter values , .
Snapshot Space 2: We could also use the method based on solving an eigenvalue problem with natural boundary conditions given by
(17) 
Where
and, in the case of linear elasticity . In a more complicated relation is related to the lower order operators [4]. Again, are the standard finescale basis functions, and his is done for each fixed parameter values , .
Define
for each coarse subdomain . We denote the corresponding matrix of snapshot functions, again with similar notation, to be
where denotes the total number of functions to keep in the snapshot construction.
Again, we perform a dimension reduction of the space of snapshots by using an auxiliary spectral decomposition. We solve the parameterindependent eigenvalue problem in the space of snapshots
(18) 
where
where and denote fine scale matrices
Here, are finescale basis functions. Further, we have
is independent of and are prescribed nonnegative weights. Recall, the main objective is to use the offline space to accurately construct a set of multiscale basis functions for each in the online stage. As before for the fluids flow module, at the offline stage of the mechanics the bilinear forms are chosen to be parameterindependent, such that there is no need to reconstruct the offline space for each .
We then choose the smallest eigenvalues from Eq. (18) and form the corresponding eigenvectors in the space of snapshots by setting , for , where are the coordinates of the vector . We denote the span of this reduced space as and denote
At the online stage, we use following parameterdependent eigenvalue problem
(19) 
where
Finally, we multiply the linear partition of unity functions by the eigenfunctions in the online space to construct the resulting basis functions
(20) 
where and denotes the number of online eigenvectors that are chosen for each coarse node . Next, we define the online space as
(21) 
Using a single index notation, we may write , where denotes the total number of basis functions in the space .
And after construction we denote the matrix
where are used to denote the nodal values of each basis function defined on the fine grid.
4.3 Global coupling
Now that we have constructed the online spaces for both the fluid and mechanics we now can use this parametrized basis at the global level. Indeed, for global coupling we use system of equations (9) to find , where
Using the matrices
we may write matrix analogue for the variational for (9) that will be used for calculation of multiscale solution Writing (9) in matrix form, using the notation in (2), in the online basis we have
(22)  
(23) 
We also note that matrices and may be analogously used in order to project coarsescale solutions onto the finegrid
5 Numerical Examples
In this section, we present numerical examples to demonstrate the performance of the GMsFEM for computing the solution of the nonlinear poroelasticity problem in heterogenous domains and complex nonlinear dependence on permeability and elastic properties. We use fully coupled scheme for approximation by time with Picard iteration to linearize the nonlinearity. We will implement a single complicated geometry with contrasting parameter values. Indeed, as noted before, there are many possible nonlinear relations, but here we take a an exponential pressure relationship with the permeability. We present the errors with varying number of multiscale basis functions and over time for linear and nonlinear case with parameters.
We proceed as in [9], and we take the computational domain as a unit square , and set the source term in (1). We utilize heterogeneous coefficients that have different values in two subdomains. We denote each region as subdomain 1 and 2, , respectively. We use following coefficients: for the Biot modulus we take in each respective numbered subdomain . For permeability we take a linear and nonlinear relation . More specifically, for the linear regime we have
(24) 
For nonlinear case we use a permeability that depends on pressure
(25) 
For fluidsolid coupling constant we have . For the elastic properties we use following coefficients: elastic modulus is given by in each respective subdomain , the Poisson’s ratio is , and these can be related to the parameters and , for via the relation
in each subdomain. The subdomains for coefficients shown in Figure 2, where the background media in red is the subdomain 1, and isolated particles and strips in blue are the subdomain 2,
As we have chosen we must use boundary conditions to force flow and mechanics. In these tests, as in [9], we use following boundary conditions:
and
and finally,
Here and are left and right boundaries, and are top and bottom boundaries respectively. We set and to drive the flow, and thus, the mechanics.
In Figure 3 we show the two coarse grids and fine grid. The first coarse grid consists of 36 nodes and 50 triangle cells, the second coarse grid contains 121 nodes and 200 triangle cells, and the fine mesh consists of 3721 nodes and 7200 triangle cells. The number of time steps is and the maximal time being set at . As an initial condition for pressure we use . For the nonlinear solve we use Picard iteration for linearization and terminate the iterative loop when , .
The reference solution computed by using a standard FEM (linear basis functions for pressure and displacements) on the fine grid, Picard type linearization, and using a fully coupled timesplitting scheme. The pressure and the displacement fields on the finegrid are presented on the left column of Figure 4 and Figure 5.
The errors will be measured in relative weighted and relative weighted norm for pressure
and for displacements, due to the linearity in our Elasticity in this example, we have
Here and are finescale and coarsescale using GMsFEM solutions, respectively for pressure and displacements.
In our examples, the nonlinearity resides in the pressure solves. Therefore, we will use the nonlinear parameter dependence approach in Section 4.1. For our Elasticity basis construction, we may remain in the linear algorithmic approach to construct the online basis. In general, for simulation using GMsFEM we first generate a snapshot space using first choice ((10), snapshot space 1) or second choice ((11), snapshot space 2), then we use a spectral decomposition to obtain the offline space, and similarly to obtain the online space. For each time step and nonlinear Picard iteration we update the online space for pressure and solve the equation (13) utilizing the previously computed solution . For construction the snapshot space 2 we choose a specified number of eigenfunctions for all . We select the range of solutions and and divide the domain into equally spaced subdomains to obtain discrete points . For simulation we use .
Recall, we will use a few multiscale basis functions per each coarse node , and these number of coarse basis defines the problem size (dimension of online spaces, and ). We suppose that in each patch we take the same number of multiscale basis functions for pressure, , for all . Similarly for displacements we take , for all . Varying the basis functions in both pressure and displacement multiscale spaces we record the errors at the final time. We note that the size of online space and the associated solution accuracy will depend on the number of eigenvectors ( and ) that we keep in the online space construction.
We begin first with the purely linear case with given by (24). In Tables 1 and 2, we present the relative weighted and errors for linear case of the coefficients in geometry Figure 2 using the fully coupled time scheme on two coarse grids. In Table 1 we have a coarsegrid of 36 nodes and in Table 2 we have a refined coarse grid with 121 nodes. We compare these to a finescale solution space with dimension 11163. In these tables, and are number of multiscale basis functions for each neighborhoods, the second column show the dimension of the online space, the next two columns present the relative weighted and errors for pressure and last two columns show the relative weighted and errors for displacements. We note that as the dimension of the online space increases, because we keep more eigenfunctions , in the space construction. We note for the less refined coarsegrid with 36 nodes the relative weighted errors decrease from 36.5% to 0.07% for pressure and from 24.3% to 0.5% for displacements and relative weighted errors decrease from 99.0% to 2.7% for pressure and from 37.7% to 3.4% for displacements. We note for the refined coarsegrid with 121 nodes the relative weighted errors decrease from 14.1% to 0.01% for pressure and from 26.9% to 0.1% for displacements and relative weighted errors decrease from 82.0% to 1.6% for pressure and from 36.1% to 2.5% for displacements. We note that in this example, refining the coarsegrid is not as advantageous to more local basis functions per gridblock. Indeed, with the less refined coarsegrid of 36 nodes and gives a very good percentage error for a space of dimension 1296 when compared to the more refined coarsegrid of 121 nodes and less eigenvectors with space of dimension 1452.
Pressure errors,  Displacements errors,  
dim()  
2  360  0.365  0.990  0.243  0.377 
4  432  0.057  0.435  0.238  0.370 
2  648  0.365  0.990  0.108  0.207 
4  720  0.057  0.435  0.045  0.077 
8  864  0.001  0.059  0.017  0.072 
2  936  0.365  0.990  0.111  0.199 
4  1008  0.057  0.435  0.042  0.045 
8  1152  0.001  0.059  0.007  0.034 
12  1296  0.0007  0.027  0.005  0.034 
Pressure errors,  Displacements errors,  
dim()  
2  1210  0.141  0.827  0.269  0.361 
4  1452  0.007  0.132  0.240  0.352 
2  2178  0.141  0.827  0.069  0.095 
4  2420  0.007  0.132  0.024  0.063 
8  1904  0.001  0.042  0.015  0.062 
2  3148  0.141  0.827  0.059  0.076 
4  3388  0.007  0.132  0.011  0.027 
8  3872  0.001  0.042  0.003  0.025 
12  4356  0.0001  0.016  0.001  0.025 
In a similar setting, we consider the nonlinear case of the coefficient with given by (25). Here we will explore the different snapshot spaces available for us in the nonlinear algorithm. Again as in the linear case we use two coarsegrids and implement this with a fully coupled time scheme and use Picard iterations for the nonlinearity. We present the results in Table 3 for snapshot space 1, the errors are very similar in magnitude when compared to the corresponding linear case. In the left side of Table 3 we present the errors for 36 nodes in the coarsegrid. The relative weighted errors decrease from 8.1% to 0.09% for pressure and from 30.4% to 0.5% for displacements and relative weighted errors decrease from 60.9% to 4.7% for pressure and from 38.0% to 3.4% for displacements. In the right side of Table 3 we present the errors for 121 nodes in the coarsegrid. The relative weighted errors decrease from 4.8% to 0.02% for pressure and from 26.4% to 0.1% for displacements and relative weighted errors decrease from 45.9% to 2.7% for pressure and from 35.7% to 2.5% for displacements. For snapshot space 2 we do precisely the same experiment with two coarsegrids. We present the errors in Table