Introduction

Front. Appl. Math. Stat.

Frontiers in Applied Mathematics and Statistics

Front. Appl. Math. Stat.

2297-4687

Frontiers Media S.A.

10.3389/fams.2023.1206500

Applied Mathematics and Statistics

Original Research

A combination of physics-informed neural networks with the fixed-stress splitting iteration for solving Biot's model

Cai

Mingchao

¹ ^* Gu

Huipeng

² Hong

Pengxiang

² Li

Jingzhi

¹Department of Mathematics, Morgan State University, Baltimore, MD, United States ²Department of Mathematics, Southern University of Science and Technology, Shenzhen, Guangdong, China

Edited by: Haizhao Yang, Purdue University, United States

Reviewed by: Kent-Andre Mardal, University of Oslo, Norway; Qifeng Liao, ShanghaiTech University, China

*Correspondence: Mingchao Cai Mingchao.Cai@morgan.edu

03 08 2023

2023

1206500

15 04 2023 14 07 2023

2023

Cai, Gu, Hong and Li

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Introduction

Biot's consolidation model in poroelasticity describes the interaction between the fluid and the deformable porous structure. Based on the fixed-stress splitting iterative method proposed by Mikelic et al. (Computat Geosci, 2013), we present a network approach to solve Biot's consolidation model using physics-informed neural networks (PINNs).

Methods

Two independent and small neural networks are used to solve the displacement and pressure variables separately. Accordingly, separate loss functions are proposed, and the fixed stress splitting iterative algorithm is used to couple these variables. Error analysis is provided to support the capability of the proposed fixed-stress splitting-based PINNs (FS-PINNs).

Results

Several numerical experiments are performed to evaluate the effectiveness and accuracy of our approach, including the pure Dirichlet problem, the mixed partial Neumann and partial Dirichlet problem, and the Barry-Mercer's problem. The performance of FS-PINNs is superior to traditional PINNs, demonstrating the effectiveness of our approach.

Discussion

Our study highlights the successful application of PINNs with the fixed-stress splitting iterative method to tackle Biot's model. The ability to use independent neural networks for displacement and pressure offers computational advantages while maintaining accuracy. The proposed approach shows promising potential for solving other similar geoscientific problems.

physics-informed neural networks the fixed-stress method Biot's model iterative algorithm separated networks

section-at-acceptance

Mathematics of Computation and Data Science

1. Introduction

Biot's consolidation model in poroelasticity describes the interaction between fluid flow and the porous structure it saturates. This model was first proposed by Biot [1], and has a wide range of applications, including biomechanics [2] and petroleum engineering [3]. The partial differential equations (PDEs) for the quasi-static Biot system in a bounded domain Ω ⊂ ℝ^d (where d = 2 or 3) over the time interval (0, T] are as follows:

(1)∂t(1Mp+α∇·u)-∇·K∇p=f in Ω×(0,T],

(2)-∇·σ(u)+α∇p=g in Ω×(0,T].

Here, u is the displacement of solid, p is the fluid pressure, g is the body force, f is a source or sink term, σ(u) = 2με(u) + λdivuI with ε(u)=12[∇u+(∇u)T] being the strain tensor. The Lamé constants λ and μ are expressed in terms of the Young's modulus E and the Poisson ratio ν as

(3)λ=Eν(1+ν)(1-2ν), μ=E2(1+ν).

Other physical parameters are the Biot-Willis constant α > 0, which is close to 1, Biot's modulus M > 0, and hydraulic conductivity K. Equation (1) represents the mass conservation and Equation (2) means forces balance. For ease of presentation, we assume the following pure Dirichlet conditions.

(4)u=ud, p=pd on ∂Ω×(0,T].

The initial conditions are

(5)u=u0, p=p0 in Ω×{0}.

The discussion of the existence and uniqueness of the solution of Biot's system (1)–(5) can be found in [4–6]. In this work, we focus on the algorithm aspect.

Several classical numerical methods have been proposed to solve this problem, including finite volume methods [7], virtual element methods [8], and mixed finite element methods [9]. Biot's model is a multiphysics problem involving both linear elasticity and porous media flow. Numerical difficulties such as elastic locking and pressure oscillations can arise, especially for models based on two-field formulations [10–13]. To overcome these difficulties, various methods have been proposed, such as the discontinuous Galerkin method [14], stabilized finite element methods [11, 15], and three-field or four-field reformulations using inf-sup stable finite element pairs [12, 16–18]. These methods may face challenges in terms of large computational overhead. The fixed-stress splitting iterative method [19, 20] is proposed to address this issue. This method breaks down the original problem into two subproblems and solves them in an iterative manner, rather than solving the entire system at once. The method relies on the contraction mapping principle to prove its convergence [19] and has been shown to be efficient through various studies. Further studies have been carried out based on the fixed-stress splitting iterative method, such as the analysis of the relationship between the convergence rate and the stabilization parameter [21–23] and the implementation of a parallel-in-time strategy to speed up the computations [24].

In recent years, deep neural networks (DNNs) have demonstrated impressive potential in solving partial differential equations (PDEs) with a wide range of applications in various domains. Among DNNs, physics-informed neural networks (PINNs) have become a popular class due to their ability to solve PDEs without meshing. PINNs have proven successful in solving high-dimensional problems and interface problems, and they can tackle inverse problems with slight modifications of the loss function [25–27]. Several studies have used PINNs to tackle the Biot's model [28–30], where the key advantage of PINNs over traditional methods such as finite element methods is their ability to avoid numerical difficulties arising from meshing and not requiring inf-sup stability. Therefore, from a flexibility perspective, deep neural network methods are preferable. However, most of these studies employ a monolithic approach to train the neural network solution, with only one study using a sequential training approach [27], which lacks theoretical analysis and only presents numerical experiments. In this paper, we propose a combination of physics-informed neural networks with fixed-stress splitting method (FS-PINNs) to solve Biot's equations. We employ two PINNs and incorporate them into the fixed-stress splitting iterative method. Our method involves two separate neural networks, one for solid displacement and another for fluid pressure, leading to faster convergence and lower computational cost than classical PINNs [31, 32]. Through a detailed analysis of the monotonic convergence of the fixed-stress splitting method, we present an error estimate for the solution of the proposed FS-PINNs. Future work could include refining the neural network architecture to improve the accuracy and further reduce the computational cost. In addition, the potential application of our approach to other PDE problems and multiphysics problems could be explored. Overall, the proposed FS-PINNs represent a promising approach for solving poroelastic models with faster convergence and lower computational cost.

The remaining sections of the paper are organized as follows. Section 2 provides an overview of the fixed-stress splitting iterative method for Biot's model. In Section 3, we introduce the fixed-stress splitting iterative PINNs (FS-PINNs) and present theoretical analyses to demonstrate their approximation properties. Section 4 presents numerical experiments to demonstrate the effectiveness of the proposed method. Finally, conclusions are drawn in Section 5.

2. Fixed-stress splitting method for Biot's equations

In this section, we introduce an iterative scheme proposed in [19], called the fixed-stress splitting method. Instead of solving Biot's model (1)–(5) in a monolithic way, this method decouples the original problem into two subproblems and solves them in an iterative manner. Given a large enough stabilization parameter β_FS and an initial guess (p⁰, u⁰), the standard fixed-stress splitting method computes a sequence of approximations as follows:

The first step: Given pⁿ and uⁿ, we solve for pⁿ⁺¹ satisfying

(6)(1M+βFS)∂tpn+1-∇·K∇pn+1=f-α∇·∂tun+βFS∂tpn in Ω×(0,T],

(7)pn+1=p0 in Ω×{0},

(8)pn+1=pd on ∂Ω×(0,T].

The second step: Using pⁿ⁺¹, we solve for uⁿ⁺¹ satisfying

(9)-∇·σ(un+1)=g-α∇pn+1 in Ω×(0,T],

(10)un+1=u0 in Ω×{0},

(11)un+1=ud on ∂Ω×(0,T].

To derive the variational formulation of the fixed-stress problem (6)–(11), we define the proper functional spaces with H¹(Ω) denoting the Hilbert subspace as follows.

W={p∈H1([0,T];H1(Ω))|p|∂Ω=pd},V={u∈H1([0,T];H1(Ω))|u|∂Ω=ud}.

Multiplying (6), (9) by test functions, and applying the integration by parts, we obtain the variational problems: for a given t ≥ 0, find {(p^{n + 1}, u^{n + 1})} ⊂ W × V such that

(12)(1M+βFS)∫Ω(∂tpn+1)ψdx+K∫Ω(∇pn+1)(∇ψ)dx=∫Ωfψdx−α∫Ω(∇·∂tun)ψdx+βFS∫Ω(∂tpn)ψdx, ∀ψ∈W,

(13)2μ∫Ωε(un+1):ε(v)dx+λ∫Ω(∇·un+1)(∇· v)dx=∫Ωgvdx+α∫Ω(∇· v)pn+1dx, ∀v∈V.

Assume that (p, u) ∈ W × V is the unique solution of Biot's system (1)–(5). Given a large enough stabilization parameter β_FS, the sequence {(pⁿ⁺¹, uⁿ⁺¹)} ⊂ W × V generated by Equations (12, 13) converges to the solution (p, u). More precisely, we can state the following:

THEOREM 2.1. The sequence {(pⁿ⁺¹, uⁿ⁺¹)} ⊂ W × V generated by Equations (6)–(11) converges to (p, u) ∈ W × V for any βFS≥α22(2μd+λ). If one denotes eun+1: =un+1-u and epn+1: =pn+1-p, then there holds

(14)∫0T∫Ω(∂tepn+1)2dxdt≤L∫0T∫Ω(∂tepn)2dxdt,

where L=βFSβFS+2M is a positive constant strictly smaller than 1.

REMARK 2.2. We have included a convergence analysis of the fixed-stress splitting iterative method in the Appendix. In our proof, we demonstrate that the differences epn+1 and eun+1 approach zeros as n tends to infinity. Our proof helps to illustrate the convergence of the fixed-stress splitting method.

3. The classical PINNs and the FS-PINNs for solving Biot's model

In this section, we introduce the training procedure of the classical PINNs for solving Biot's model, and propose an iterative deep learning method combining the fixed-stress splitting method, called FS-PINNs. The proposed FS-PINNs consist of two independent PINNs, one for pressure, and the other for displacement. We then present the training procedure and a theoretical analysis of the proposed FS-PINNs.

3.1. Training procedure of the classical PINNs

The idea of using the classical PINNs to solve PDEs is from [25], which can be easily extended to solve Biot's model. The main idea is to design loss functions for fully connected neural networks. In this context, we provide a brief overview of the training procedure for solving problem (1)–(5), which is illustrated in Figure 1.

Figure 1

The training procedure of the classical PINNs using unified activation functions for a 2D Biot's model.

The functional form of a classical PINN is given as follows.

(15)n^(x;θn)=(p^(x;θn);u^(x;θn))=ΣnLn◦ΣnLn-1◦⋯◦Σn1(x).

Here, n^ is the network approximation of the exact solution (p; u) with L_n hidden layers with the input x=(x1,⋯,xd,t)∈ℝd×[0,T]. The collection of parameter θ_n is given by

(16)θn={Wn1,Wn2,⋯,WnLn,bn1,bn2,⋯,bnLn}.

The mapping function Σni is defined as

(17)xni: =Σni(xni-1)=σni(Wnixni-1+bni), i=1,2,⋯,Ln.

Here, (Wni,bni,σni) represent the weights, biases, and activation functions of the i-th layer for network n^. It is assumed that the number of neurons in each hidden layer is set to be the same N_n and that all activation functions σni are identical.

After setting up the architecture, we can begin training the parameter θ_n using the procedure:

(18)θ*=arg min θnLn(x;θn),

where Ln is the loss function. The loss function for the classical PINNs generally consists of three components: the PDEs loss, the initial condition loss, and the boundary condition loss. The PDEs loss represents the residual of the governing equations, while the initial and boundary condition losses ensure that the predicted solution satisfies the initial and boundary conditions. The loss function is denoted by Ln(x;θn) and has the following representation.

Ln(x;θn)=1N1∑i=1N1(‖p^(x1i;θn)−pi‖L22+‖u^(x1i;θn)−ui‖L22) +1N2∑i=1N2(‖1M∂tp^(x2i;θn)−∇·K∇p^(x2i;θn)+ α∇ · ∂tu^(x2i;θn)−f(x2i)‖L22 +‖−∇·σ(u^(x2i;θn))+ α∇p^(x2i;θn)−g(x2i)‖L22).

The above expression uses the mean square error (MSE) to measure errors. In this expression, {x1i}i=1N1 represents the initial and boundary training data, while {x2i}i=1N2 represents the collocation points in the domain. We then optimize θ_n to obtain θ^*. For cases involving mixed boundary conditions, one can refer to [28] for more details.

3.2. Training procedure of the proposed FS-PINNs

The primary innovation of our approach is incorporating the fixed-stress strategy into the training procedure of PINNs. The proposed method includes two independent PINNs to solve for displacement and fluid pressure separately. We can then train our FS-PINNs iteratively with respect to their corresponding loss functions. Figure 2 illustrates the training procedure.

Figure 2

The training procedure of the proposed FS-PINNs. Two independent networks are trained iteratively with respect to the corresponding loss functions.

Similar to the classical PINNs, the functional forms of the two independent PINNs are as follows.

(19)p^(x;θp)=ΣpLp◦ΣpLp-1◦⋯◦Σp1(x),

(20)u^(x;θu)=ΣuLu◦ΣuLu-1◦⋯◦Σu1(x).

Here, p^ refers to the network approximation of pressure p with L_p hidden layers, while u^ refers to the network approximation of displacement u with L_u hidden layers. The collection of parameters of p^ and u^ are denoted by θ_p and θ_u, respectively, and are given by:

(21)θp={Wp1,Wp2,⋯,WpLp,bp1,bp2,⋯,bpLp},

(22)θu={Wu1,Wu2,⋯,WuLu,bu1,bu2,⋯,buLu}.

The mapping functions Σpi and Σui are defined as

(23)xpi: =Σpi(xpi-1)=σpi(Wpixpi-1+bpi), i=1,2,⋯,Lp,

(24)xui: =Σui(xui-1)=σui(Wuixui-1+bui), i=1,2,⋯,Lu.

In this context, (Wpi,bpi,σpi) and (Wui,bui,σui) represent the weights, biases, and activation functions of the i-th layer for p^ and u^, respectively. It is also assumed that the number of neurons in each hidden layer of p^ and u^ is set to N_p and N_u, respectively. Furthermore, all activation functions σpi and σui are assumed to be identical.

From now, our attention turns to developing the loss functions. By utilizing the fixed-stress splitting method detailed in Equations (6)–(8) and Equations (9)–(11), we can impose constraints on the networks p^ and u^, respectively. We begin by setting n = 0 and θpn=θun=0. We can then formulate the optimization problem as follows. In the first step, we obtain θpn+1 by solving for the minimum of Lp(x,θun,θpn;θp):

(25)θpn+1=arg min θpLp(x,θun,θpn;θp).

In the second step, we obtain θun+1 by solving for the minimum of Lu(x,θpn+1;θu):

(26)θun+1=arg min θuLu(x,θpn+1;θu).

To be clear, Lp and Lu denote the loss functions for networks p^ and u^, respectively. These loss functions are defined as follows:

(27)Lp(x,θun,θpn;θp)=1N3∑i=1N3‖p^(xpi;θp)−p(xpi)‖2 +1N4∑i=1N4‖(1M+βFS)∂tp^(xfi;θp)−∇·K∇p^(xfi;θp) +α ∇·∂tu^(xfi;θun)−βFS∂tp^(xfi;θpn)−f(xfi)‖2,

(28)Lu(x,θpn+1;θu)=1N5∑i=1N5‖u^(xui;θu)−u(xui)‖2 +1N6∑i=1N6‖−∇·σ(u^(xgi;θu))+α∇p^(xgi;θpn+1)−g(xgi)‖2.

In the above expressions, {xpi}i=1N3 and {xui}i=1N5 denote the initial and boundary training data on p(x) and u(x), respectively. Additionally, {xfi}i=1N4 and {xgi}i=1N6 represent the collocation points for f(x) and g(x) in the domain, respectively. By repeatedly updating θ_p and θ_u using (25) and (26), we can get the network approximations of Biot's model: p^(x;θp) and u^(x;θu).

3.3. Analysis of the proposed FS-PINNs

The classical PINNs have shown great potential in solving PDEs, but quantifying their errors remains an open problem [26]. The accuracy of neural network approximation is influenced by various factors, such as network architecture, training data, etc. In this section, we assume that PINNs can solve the two subproblems within a certain error tolerance. We then employ an iterative strategy for the proposed FS-PINNs. It is crucial to consider the possibility of error accumulation and divergence since network approximations can not guarantee convergence of the fixed-stress iteration. To address this issue, we observe that the fixed-stress splitting method, which is essentially equivalent to the fixed-point iterative scheme, has been shown to mitigate error accumulation. Therefore, we present a mathematical analysis to demonstrate the robustness of our proposed FS-PINNs. Given an initial guess (p^0,u^0), we describe the procedures of (25) and (26) as follows.

Step 1: Given p^n and u^n, we employ (25) to obtain a network solution p^n+1 for the following problem.

(29)(1M+βFS)∂tpn+1-∇·K∇pn+1=f-α∇·∂tu^n+βFS∂tp^n inΩ×(0,T],

(30)pn+1=p0 in Ω×{0},

(31)pn+1=pd on ∂Ω×(0,T].

The corresponding loss function for p^n+1 imposes the initial condition and the boundary condition in the mean square sense.

Step 2: using p^n+1, we employ (26) to obtain a network solution u^n+1 for the following problem.

(32)-∇·σ(un+1)=g-α∇p^n+1 inΩ×(0,T],

(33)un+1=u0 in Ω×{0},

(34)un+1=ud on ∂Ω×(0,T].

The corresponding loss function for u^n+1 imposes the initial condition and the boundary condition in the mean square sense.

To analyze the numerical errors p^n+1-p and u^n+1-u, we can express them as the sum of two terms: (p^n+1-pn+1)+(pn+1-p) and (u^n+1-un+1)+(un+1-u), which can be dealt with separately. We assume that the Physics-Informed Neural Networks (PINNs) used to solve Equations (29) and (32) can achieve a certain error tolerance, denoted by δ, for the approximations p^n+1 and u^n+1 obtained in each step. This is expressed as

(35)‖p^n+1−pn+1‖H1([0,T];H1(Ω))≤δ, ‖u^n+1−un+1‖H1([0,T];H1(Ω))≤δ.

Such assumptions are reasonable if the neural networks are sufficiently deep and wide [33]. For estimating the terms (pⁿ⁺¹−p) and (uⁿ⁺¹−u), we utilize integration by parts and employ the fixed-stress iterative method to analyze convergence. We then build on Theorem 2.1 to perform mathematical analysis.

THEOREM 3.1. Assume that (p, u) ∈ W × V is the unique solution of the problem (1)–(5). Let (p^n+1,u^n+1) be the network approximations of Equations (29)–(34), and (pⁿ⁺¹, uⁿ⁺¹) ∈ W × V be the exact solution of (6)–(11). We define eun+1: =un+1-u and epn+1: =pn+1-p. Under the assumption (35) with a positive δ, if β_FS is large enough, then there exists a positive number L~<1 such that the following estimate holds for all n ≥ 0.

(36)∫0T∫Ω(∂tepn+1)2dxdt≤L∼∫0T∫Ω(∂tepn)2dxdt+Cδ2,

where C is a positive constant. Moreover, the limit superior exists and satisfies:

(37)lim supn→∞∫0T∫Ω(∂tepn+1)2dxdt≤Cδ21−L∼.

Proof. Firstly, we start with the second step of FS-PINNs. We subtract the original Equation (2) from the n-th iteration of (32). For the resulted equation, taking the derivative with respect to t, we obtain

(38)-∇·σ(∂teun)=-α∇∂tepn+α∇(∂tpn-∂tp^n).

Multiplying (38) by ∂teun, applying the divergence theorem, the Cauchy-Schwarz inequality, and the Poincaré inequality, we shall derive the following inequality with any ϵ₁ > 0.

(39)2μ∫0T∫Ωε(∂teun):ε(∂teun)dxdt+λ∫0T∫Ω(∇·∂teun)2dxdt =α∫0T∫Ω(∇·∂teun)(∂tepn)dxdt+α∫0T∫Ω(∇(∂tp^n−∂tpn))(∂teun)dxdt ≤12(2μd+λ)∫0T∫Ω(∇·∂teun)2dxdt+α22(2μd+λ)∫0T∫Ω(∂tepn)2dxdt +α24ϵ1∫0T∫Ω(∇(∂tp^n−∂tpn))2dxdt+Cpϵ1∫0T∫Ω(∇·∂teun)2dxdt.

On the other hand, there holds [23]

(2μd+λ)∫0T∫Ω(∇·∂teun)2dxdt≤2μ∫0T∫Ωε(∂teun):ε(∂teun)dxdt+λ∫0T∫Ω(∇·∂teun)2dxdt.

Based on the above inequality and (35), we can deduce the following inequality from (39).

(40)2μ∫0T∫Ωε(∂teun):ε(∂teun)dxdt+λ∫0T∫Ω(∇·∂teun)2dxdt≤(1+Cϵ1)[α22μd+λ∫0T∫Ω(∂tepn)2dxdt+α22ϵ1δ2].

Similar to the derivation of (38), subtracting the (n+1)-th iteration Equation (32) from (2) and taking the derivative with respect to t, we obtain

(41)-∇·σ(∂teun+1)=-α∇∂tepn+1+α∇(∂tpn+1-∂tp^n+1).

Subtracting (38) from (41) yields

(42)-∇·σ(∂teun+1-∂teun)=-α∇(∂tepn+1-∂tepn)+α∇[(∂tpn+1-∂tp^n+1)-(∂tpn-∂tp^n)].

We multiply (42) by ∂teun+1-∂teun, and repeat the same argument as that for (40) to derive

(43)2μ∫0T∫Ωε(∂teun+1−∂teun):ε(∂teun+1−∂teun)dxdt +λ∫0T∫Ω(∇·(∂teun+1−∂teun))2dxdt≤(1+Cϵ1)[α22μd+λ∫0T∫Ω(∂tepn+1−∂tepn)2dxdt+2α2ϵ1δ2].

Next, we consider the entire system of the FS-PINNs. Subtracting the original Equation (1) from the (n+1)-th iteration Equation (29) yields

(44)(1M+βFS)∂tepn+1-∇·K∇epn+1=-α∇·∂teun +βFS∂tepn+α∇·(∂tun-∂tu^n)-βFS(∂tpn-∂tp^n).

Multiplying (44) by ∂tepn+1, and (41) by ∂teun, respectively, we then take a summation to derive the following equation.

(45)LHS :=1M∫0T∫Ω(∂tepn+1)2dxdt+K2∫Ω(∇epn+1(T))2dx+ 2μ∫0T∫Ωε(∂teun+1):ε(∂teun)dxdt +λ∫0T∫Ω(∇·∂teun+1)(∇·∂teun)dxdt+ βFS∫0T∫Ω(∂tepn+1)(∂tepn+1−∂tepn)dxdt =−α∫0T∫Ω(∇·(∂tu^n−∂tun))(∂tepn+1)dxdt+ βFS∫0T∫Ω(∂tp^n−∂tpn)(∂tepn+1)dxdt +α∫0T∫Ω(∇(∂tp^n+1−∂tpn+1))(∂teun)dxdt=:RHS.

We recall the following algebraic identities [“(·, ·)” means the l² inner product].

(46)(ξ−η,ξ)=12[(ξ,ξ)−(η,η)+(ξ−η,ξ−η)],(ξ,η)=14[(ξ+η,ξ+η)−(ξ−η,ξ−η)].

Applying (46) to the LHS of (45), and discarding some positive terms, we shall obtain

(47)LHS≥(1M+βFS2)∫0T∫Ω(∂tepn+1)2dxdt− βFS2∫0T∫Ω(∂tepn)2dxdtn −μ2∫0T∫Ωε(∂teun+1−∂teun):ε(∂teun+1−∂teun)dxdt −λ4∫0T∫Ω(∇·(∂teun+1−∂teun))2dxdt+ βFS2∫0T∫Ω(∂tepn+1−∂tepn)2dxdt.

We then deduce the following inequality based on (43).

(48)LHS≥(1M+βFS2)∫0T∫Ω(∂tepn+1)2dxdt− βFS2∫0T∫Ω(∂tepn)2dxdt +(βFS2−1+Cϵ12α22μd+λ)∫0T∫Ω(∂tepn+1−∂tepn)2dxdt −(1+Cϵ1)α2ϵ1δ2.

For estimating an upper bound of the RHS term in (45), we use the Cauchy-Schwarz inequality with any ϵ₂ > 0, and (40). We see that

(49)RHS≤ϵ2∫0T∫Ω(∂tepn+1)2dxdt+ Cpϵ2∫0T∫Ω(∇·∂teun)2dxdt+2α2+βFS4ϵ2δ2 ≤ϵ2∫0T∫Ω(∂tepn+1)2dxdt+Cpϵ2(1+Cϵ1)α22μd+λ∫0T ∫Ω(∂tepn)2dxdt +(2α2+βFS4ϵ2+Cpϵ2(1+Cϵ1)α22ϵ1)δ2.

We let βFS≥(1+Cϵ1)α22μd+λ, by combining (48) and (49), we see that

(50)(1M+βFS2−ϵ2)∫0T∫Ω(∂tepn+1)2dxdt≤(βFS2+C1(ϵ1,ϵ2))∫0T∫Ω(∂tepn)2dxdt+C2(ϵ1,ϵ2)δ2,

where

C1(ϵ1,ϵ2)=Cpϵ2(1+Cϵ1)α22μd+λ,C2(ϵ1,ϵ2)=2α2+βFS4ϵ2+(Cpϵ22+1)(1+Cϵ1)α2ϵ1.

We choose ϵ₁ and ϵ₂ small enough, such that the contraction coefficient L~, defined as follows, is strictly smaller than 1.

L~: =βFS+2C1(ϵ1,ϵ2)2M+βFS-2ϵ2<1.

Then, we can rewrite (50) in the following form.

(51)∫0T∫Ω(∂tepn+1)2dxdt≤L∼∫0T∫Ω(∂tepn)2dxdt+Cδ2.

This leads to the limit supremum.

(52)lim supn→∞∫0T∫Ω(∂tepn+1)2dxdt≤lim supn→∞[Ln+1∼∫0T∫Ω(∂tep0)2dxdt+(Ln∼+⋯+1)Cδ2]≤Cδ21−L∼.

This completes our proof.

4. Benchmark tests

In this section, the efficacy of the proposed FS-PINNs is presented through three case studies. The first study involves a pure Dirichlet boundary value problem, while the second study considers a problem with mixed Dirichlet-Neumann boundary conditions. Finally, the third study investigates Barry-Mercer's problem. All experiments were implemented in the open-source machine learning tool PyTorch [34] on a single NVIDIA GeForce RTX 3060 GPU. Our code is available at https://github.com/newpolarbear/FS-PINNs. We evaluate the computational efficiency of our proposed method by comparing its true errors and runtime overheads with those of classical PINNs.

4.1. Models generation

In practice, the training procedure of classical PINNs is heavily influenced by several factors. One of these factors is the design of the neural network architecture, which does not have a universal standard. It's well-known that the network's representational capability is affected by its depth and width. While deeper and wider networks have greater representation capacity, they may suffer from slower training, compromises can be made to reduce runtime overhead. Another factor is the choice of nonlinear activation functions and optimization methods. In this work, we used tanh(x)=ex-e-xex+e-x as the activation function and Adam optimizer with a fixed learning rate of 0.0001. In addition, some factors are specific to the problem being solved, such as the number of epochs, collocation points, and initial and boundary training data. In this study, we randomly select 10, 000 points as collocation points within the domain. The initial training data and each boundary training data are set to 500 and 200, respectively.

In addition to the standard configurations of classical PINNs, there are other aspects that come into play when training FS-PINNs. One FS-PINN is composed of two separate networks, which directly impact the performance. Moreover, training FS-PINNs involves networks optimization for each fixed-stress iteration, as shown in Figure 2. In the n-th fixed-stress iteration, we carry out the optimization procedure of network parameters θ_p and θ_u. For simplicity, we assume that the number of epochs for optimizing θ_p equals that for optimizing θ_u. Since the fixed-stress splitting method is an iterative algorithm, achieving accurate solutions at each step is not so necessary. For simplicity, we set the fixed-stress stabilization parameter to βFS=α2/(2μd+λ).

The training procedure of FS-PINNs is dynamic, loss functions Lp and Lu can only show the performance in the current fixed-stress iteration. To evaluate the proposed method, we consider the true errors between network solutions and exact solutions as follows.

MSE for pressure p:MSEp=1N∑i=1N||p^(xi)-p(xi)||2,MSE forx-displacement u1:MSEu1=1N∑i=1N||û1(xi)-u1(xi)||2,MSE fory-displacement u2:MSEu2=1N∑i=1N||û2(xi)-u2(xi)||2.

Here, p^, û₁, and û₂ represent the neural network solutions of pressure, x-displacement, and y-displacement, respectively, while p, u₁, and u₂ denote their corresponding exact solutions. The total mean squared error can be computed by taking the summation of these errors:

(53)Total MSE: MSEtotal=MSEp+MSEu1+MSEu2.

4.2. Example 1

The first numerical experiment is based on the example in [16]. The domain is Ω = [0, 1]² and the final time is T = 0.5, where pure Dirichlet boundary conditions are applied. In this experiment, the source or sink term f, body force g, initial conditions, and boundary conditions are chosen such that the exact solutions are provided as follows:

(54)u(x,y,t)=t2(x2,y2)T, p(x,y,t)=sin(x+y)et.

The model parameters are set to E = 1.0, ν = 0.3, α = 1.0, M = 1.0, and K = 1.0.

4.2.1. Robustness of the FS-PINNs architecture

The performance of FS-PINNs is highly dependent on its architecture. We fixed the total width to N_p+N_u = 30 and depths to L_p = L_u = 3, where N_p, N_u are the number of neurons, and L_p, L_u are the number of layers. We then designed an experiment with five sets of parameter selection for N_p:N_u. To ensure that the networks are well-trained, we applied 20 fixed-stress iterations and 3,000 epochs for each fixed-stress step. The experiment is tested repeatedly to ensure the reliability of the results.

Table 1 shows the runtime overheads for various parameter selections, revealing that a more uniform separation exhibits lower runtime overhead. Notably, the case N_p = 5:N_u = 25 performs poorly compared to other cases, because wider networks experience slower training, which leads to a significant increase in the runtime overhead of network û. Conversely, the other four cases exhibit similar performance, suggesting that the proposed FS-PINNs with reasonable architecture are resilient to runtime overhead.

Table 1

Runtime overheads of five different parameter selections for FS-PINNs architecture of Example 1.

Parameter selection	Runtime for network p^	Runtime for network u^	Total runtime
N_p = 5:N_u = 25	514.9 s (35.4%)	939.1 s (64.6%)	1454.0 s
N_p = 10:N_u = 20	471.4 s (40.8%)	685.4 s (59.2%)	1156.8 s
N_p = 15:N_u = 15	497.9 s (42.0%)	688.0 s (58.0%)	1185.9 s
N_p = 20:N_u = 10	483.0 s (43.1%)	637.0 s (56.9%)	1120.0 s
N_p = 25:N_u = 5	569.9 s (47.4%)	632.5 s (52.6%)	1202.4 s

Figure 3 tracks the trends of all MSE terms during different parameter selections training procedures, revealing that the case N_p = 25:N_u = 5 performs worst, and the case N_p = 5:N_u = 25 is less stable. The inadequate representation ability of tiny networks leads to inaccurate solutions. Network approximations p^ and u^ depend on the interdependent loss functions for FS-PINNs. On the other hand, the trends of the remaining three cases are quite similar, which implies that FS-PINNs with reasonable architecture can withstand accuracy issues.

Figure 3

Trends of all the MSE terms at T = 0.5 for five different FS-PINNs architecture during training with respect to epoch.

4.2.2. Acceleration effect of the FS-PINNs training

This section demonstrates the acceleration effect of the FS-PINNs over the classical PINNs. Biot's equations are solved using both methods with the same setup as in Example 1. The number of neurons for the two networks in FS-PINNs is fixed to N_p = N_u = 15 (as indicated by Table 1), while the total number of neurons in classical PINNs is set to N_n = 30 with the number of layers L_n = 3. To examine the effect of the number of epochs optimizing networks and the number of fixed-stress iterations, we then conduct the test with a constant number of the multiplication of the above two number, which is set to 20, 000.

Figure 4 shows the trends of MSE terms for solving the problem in Example 1. We present the results generated by one classical PINN and four FS-PINNs under identical conditions, which demonstrates the acceleration effect of FS-PINNs. By decoupling the original optimization problem into two small sub-problems, together with the fixed-stress coefficient L < 1, FS-PINNs increase efficiency and maintain certain accuracy. From this perspective, FS-PINNs offer more flexibility in solving the problem through different setting of optimizing epoch and fixed-stress iteration. In Figure 4, we can see that a smaller number of epochs lead to faster convergence initially but lower accuracy eventually. Since the fixed-stress splitting method is a numerical approach, high accuracy is not essential for each fixed-stress step. It implies that a large number of optimizing epochs is not so helpful.

Figure 4

Trends of all the MSE terms at T = 0.5 for one classical PINN and four FS-PINNs with [the number of epochs training networks, the number of fixed-stress iterations] during training with respect to runtime overhead.

Among all the results generated by FS-PINNs, here we consider in-depth observations for the case of [400, 50]. Here, 400 means the number of epochs, and 50 means the number of fixed-stress iterations. To test the cases with a larger parameter λ, we also include the testing result when the Poisson ratio ν = 0.49 for Example 1. Figure 5 provides a comparison between the analytical solution and the numerical solution given by FS-PINNs with [400, 50]. Theses results show that the proposed FS-PINNs with ν = 0.3 and ν = 0.49 accurately predict the pressure p and displacement u with errors within an acceptable tolerance range. However, it is true that when λ is large, to obtain a more accurate solution, we need to put a relatively large weight for the first equation in the total loss function. And also, we need to increase the number of epochs and set a finer learning rate for the training of the network for u.

Figure 5

A comparison of the analytical solution (left), the FS-PINN solution with ν = 0.3 (middle), and the FS-PINN solution ν = 0.49 (right) at T = 0.5 for Example 1.

4.3. Example 2

The Biot's problem with mixed boundary conditions presented in [2] is tested. We consider the square domain Ω = [0, 1]², with Dirichlet boundary Γ₁∪Γ₃ and Neumann boundary Γ₂∪Γ₄. Specifically, Γ₁ = (1, y);0 ≤ y ≤ 1, Γ₂ = (x, 0);0 ≤ x ≤ 1, Γ₃ = (0, y);0 ≤ y ≤ 1, and Γ₄ = (x, 1);0 ≤ x ≤ 1. Setting the final time to T = 0.2, we choose the body force g, the source or sink term f, as well as the initial and boundary conditions such that exact solutions are as follows.

u1=e-t(sin(2πy)(cos(2πx)-1)+1μ+λsin(πx)sin(πy)),u2=e-t(sin(2πx)(1-cos(2πy))+1μ+λsin(πx)sin(πy)),p=e-tsin(πx)sin(πy).

The parameters of the model are E = 1.0, ν = 0.3, α = 1.0, M = 1.0, K = 1.

In dealing with mixed boundary conditions, it is necessary to select suitable loss functions capable of handling both Dirichlet and Neumann boundary conditions. To this end, we follow the approach taken in [28, 32] and employ the corresponding loss functions Lp and Lu for pressure and displacement, respectively. In this test, following the result in Section 4.2.2., we employ the FS-PINN with [400, 50] using the same architecture to solve Example 2.

In Figure 6, we compare the analytical solution with our FS-PINN solution and observe that our approach can approximate solutions for problems with mixed boundary conditions. However, its performance is not as good as that of solving the pure Dirichlet problem in Example 1. The network approximations in the boundary region are sub-optimal because all the parameters used were selected to solve Example 1. Hence, the complexity of Example 2 presents a challenge to our network, which could potentially be overcome by optimizing the parameters, such as by using additional collocation points in the boundary region.

Figure 6

A comparison of the analytical solution (left) and the FS-PINN solution with [400, 50] (right) at T = 0.2 for Example 2.

4.4. Example 3

The Barry-Mercer's problem is a classical benchmark problem of porous media involving injection of fluids [14, 31, 35, 36]. It is a two-dimensional problem in the square domain Ω = [0, 1]² that has a point source f located at (x₀, y₀) = (0.25, 0.25), which makes it challenging to solve numerically. Figure 7 depicts the boundary conditions under consideration.

Figure 7

Boundary conditions for the Barry-Mercer's problem.

As for the initial conditions, we set

(55)u=0, p=0 in Ω×{0}.

To define the body force term g and the source/sink term f, we use the following expressions:

(56)g=0, f=2δ(x-x0)δ(y-y0)sin(t),

where δ denotes the Dirac function. We approximate the Dirac function using the Gaussian distribution [31], that is, δ(x)≈1cgπe-(x/cg)2 with c_g = 0.04. The physical parameters are specified as follows:

(57) 1/M=0, α=1.0, λ=0.2, μ=0.4, K=1.0.

Denote λ_n = nπ, λ_q = qπ, and λnq=λn2+λq2, then the analytical solution is as follows:

p(x,y,t)=-4∑n=1∞∑q=1∞p^(n,q,t)sin(λnx)sin(λqy),u1(x,y,t)=4∑n=1∞∑q=1∞û1(n,q,t)cos(λnx)sin(λqy),u2(x,y,t)=4∑n=1∞∑q=1∞û2(n,q,t)sin(λnx)cos(λqy),

where

p^(n,q,t)=-2sin(λnx0)sin(λqy0)λnq2+1(λnqsint-cost+e-λnqt),û1(n,q,t)=λnλnqp^(n,q,t), û2(n,q,t)=λqλnqp^(n,q,t).

In this case, we use the same architecture as in Section 4.2.2, which is the FS-PINN with [400, 50], to solve Barry-Mercer's problem. As mentioned in the literature [37], networks training relies on the propagation of information from initial and boundary points to the interior points. To improve the performance of our FS-PINN, we apply weighted loss functions proposed in [32], and increase the weights of the initial and boundary loss terms. Specifically, we replace 1N3, 1N5 with 1000N3, 100N5 in the loss functions (27), (28), respectively. The solution obtained from our FS-PINN closely matches the analytic solution at T = π/4 for Example 3, as demonstrated in Figure 8. This indicates the effectiveness of our approach in approximating the analytic solution.

Figure 8

A comparison of the analytical solution (left) and the FS-PINN solution (right) at T = π/4 for Example 3.

5. Conclusion

This paper presents a new approach to solving Biot's consolidation model by combining Physics-informed neural networks (PINNs) with the fixed-stress splitting method. We demonstrate that our method achieves faster convergence and lower computational costs than the neural network methods which solve Biot's problem in a monolithic way. We provide a theoretical analysis to show the convergence of our proposed FS-PINNs. The numerical experiments illustrate the effectiveness of our approach on various test cases with different boundary conditions and physical parameters.

While our approach has several advantages over some traditional methods for solving poroelastic models, such as being meshless, not requiring inf-sup stability, there are also some limitations and challenges that need to be addressed in future work. These include choosing optimal network architectures and parameters for different problems, improving accuracy and reducing computational costs even further, and extending the approach to other PDE problems and multi-physics problems. We hope that our work can inspire more research on combining deep learning techniques with iterative methods for solving PDEs.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

MC: conceptualization, investigation, writing—original draft, validation, methodology, writing—review and editing, supervision, formal analysis, and funding acquisition. HG: investigation, writing—original draft, validation, visualization, writing—review and editing, formal analysis, and data curation. PH: investigation, writing—original draft, validation, visualization, writing—review and editing, software, and data curation. JL: conceptualization, investigation, funding acquisition, validation, methodology, project administration, and supervision. All authors contributed to the article and approved the submitted version.

Funding

The work of MC was partially supported by the NIH-RCMI grant through 347 U54MD013376, the affiliated project award from the Center for Equitable Artificial Intelligence and Machine Learning Systems (CEAMLS) at Morgan State University (project ID 02232301), and the National Science Foundation awards (1831950 and 2228010). The work of HG, PH, and JL were partially supported by the NSF of China No. 11971221, Guangdong NSF Major Fund No. 2021ZDZX1001, the Shenzhen Sci-Tech Fund Nos. RCJC20200714114556020, JCYJ20200109115422828, and JCYJ20190809150413261, National Center for Applied Mathematics Shenzhen, and SUSTech International Center for Mathematics.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References 1. Biot

. General theory of three-dimensional consolidation. J Appl Phys. (1941) 12:155–64. 10.1063/1.1712886 2. Ju

Cai

Tian

. Parameter-robust multiphysics algorithms for Biot model with application in brain edema simulation. Math Comput Simul. (2020) 177:385–403. 10.1016/j.matcom.2020.04.027

32546894

3. Kim

Tchelepi

Juanes

. Stability and convergence of sequential methods for coupled flow and geomechanics: fixed-stress and fixed-strain splits. Comput Methods Appl Mech Eng. (2011) 200:1591–606. 10.1016/j.cma.2010.12.022 4. Ženíšek

. The existence and uniqueness theorem in Biot's consolidation theory. Aplik Matem. (1984) 29:194–211. 10.21136/AM.1984.104085 5. Showalter

. Diffusion in poro-elastic media. J Math Anal Appl. (2000) 251:310–40. 10.1006/jmaa.2000.7048

36971445

6. Phillips

Wheeler

. A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous in time case. Comput Geosci. (2007) 11:131–44. 10.1007/s10596-007-9045-y 7. Nordbotten

. Stable cell-centered finite volume discretization for Biot equations. SIAM J Numer Anal. (2016) 54:942–68. 10.1137/15M1014280 8. Coulet

Faille

Girault

Guy

Nataf

. A fully coupled scheme using virtual element method and finite volume for poroelasticity. Comput Geosci. (2020) 24:381–403. 10.1007/s10596-019-09831-w 9. Yi

. A coupling of nonconforming and mixed finite element methods for Biot's consolidation model. Numer Methods Part Diff Equat. (2013) 29:1749–77. 10.1002/num.21775 10. Yi

. A study of two modes of locking in poroelasticity. SIAM J Numer Anal. (2017) 55:1915–36. 10.1137/16M1056109 11. Rodrigo

Gaspar

Zikatanov

. Stability and monotonicity for some discretizations of the Biot's consolidation model. Comput Methods Appl Mech Eng. (2016) 298:183–204. 10.1016/j.cma.2015.09.019 12. Lee

Mardal

Winther

. Parameter-robust discretization and preconditioning of Biot's consolidation model. SIAM J Sci Comput. (2017) 39:A1–24. 10.1137/15M1029473 13. Cai

Zhang

. Comparisons of some iterative algorithms for Biot equations. Int J Evol Equat. (2017) 10:267–82.30364836 14. Phillips

Wheeler

. Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach. Comput Geosci. (2009) 13:5–12. 10.1007/s10596-008-9114-x 15. Rodrigo

Ohm

Adler

Gaspar

Zikatanov

. New stabilized discretizations for poroelasticity and the Stoke's equations. Comput Methods Appl Mech Eng. (2018) 341:467–84. 10.1016/j.cma.2018.07.003 16. Feng

. Analysis of a multiphysics finite element method for a poroelasticity model. IMA J Numer Anal. (2018) 38:330–59. 10.1093/imanum/drx003 17. Oyarzúa

Ruiz-Baier

. Locking-free finite element methods for poroelasticity. SIAM Journal on Numerical Anal. (2016) 54:2951–73. 10.1137/15M1050082 18. Yi

Bean

. Iteratively coupled solution strategies for a four-field mixed finite element method for poroelasticity. Int J Numer Analyt Methods Geomech. (2017) 41:159–79. 10.1002/nag.2538 19. Mikelić

Wheeler

. Convergence of iterative coupling for coupled flow and geomechanics. Comput Geosci. (2013) 17:455–61. 10.1007/s10596-012-9318-y 20. Mikelić

Wang

Wheeler

. Numerical convergence study of iterative coupling for coupled flow and geomechanics. Comput Geosci. (2014) 18:325–41. 10.1007/s10596-013-9393-8 21. Both

Borregales

Nordbotten

Kumar

Radu

. Robust fixed stress splitting for Biot's equations in heterogeneous media. Appl Math Lett. (2017) 68:101–8. 10.1016/j.aml.2016.12.019 22. Bause

Radu

Köcher

. Space-time finite element approximation of the Biot poroelasticity system with iterative coupling. Comput Methods Appl Mech Eng. (2017) 320:745–68. 10.1016/j.cma.2017.03.017 23. Storvik

Both

Kumar

Nordbotten

Radu

. On the optimization of the fixed-stress splitting for Biot's equations. Int J Numer Methods Eng. (2019) 120:179–94. 10.1002/nme.6130 24. Borregales

Kumar

Radu

Rodrigo

Gaspar

. A partially parallel-in-time fixed-stress splitting method for Biot's consolidation model. Comput Math Appl. (2019) 77:1466–78. 10.1016/j.camwa.2018.09.005 25. Raissi

Perdikaris

Karniadakis

. Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys. (2019) 378:686–707. 10.1016/j.jcp.2018.10.045 26. Lu

Meng

Mao

Karniadakis

. DeepXDE: a deep learning library for solving differential equations. SIAM Rev. (2021) 63:208–28. 10.1137/19M1274067 27. Amini

Haghighat

Juanes

. Inverse modeling of nonisothermal multiphase poromechanics using physics-informed neural networks. arXiv preprint arXiv:220903276. (2022). 10.1016/j.jcp.2023.112323 28. Kadeethum

Jørgensen

Nick

. Physics-informed neural networks for solving nonlinear diffusivity and Biot's equations. PLoS ONE. (2020) 15:e0232683. 10.1371/journal.pone.0232683

32374751

29. Millevoi

Spiezia

Ferronato

. On Physics-Informed Neural Networks Architecture for Coupled Hydro-Poromechanical Problems. (2022). 10.2139/ssrn.4074416 30. Bekele

. Physics-informed deep learning for flow and deformation in poroelastic media. arXiv preprint arXiv:201015426. (2020). 31. Haghighat

Amini

Juanes

. Physics-informed neural network simulation of multiphase poroelasticity using stress-split sequential training. Comput Methods Appl Mech Eng. (2022) 397:115141. 10.1016/j.cma.2022.115141 32. Amini

Haghighat

Juanes

. Physics-informed neural network solution of thermo-hydro-mechanical (THM) processes in porous media. arXiv preprint arXiv:220301514. (2022). 10.1061/(ASCE)EM.1943-7889.0002156

29515898

33. De Ryck

Jagtap

Mishra

. Error estimates for physics informed neural networks approximating the Navier-Stokes equations. arXiv preprint arXiv:220309346. (2022). 10.1093/imanum/drac085 34. Paszke

Gross

Massa

Lerer

Bradbury

Chanan

. PyTorch: An imperative style, high-performance deep learning library. In:Wallach

Larochelle

Beygelzimer

Alche-Buc

Fox

Garnett

, editors. Advances in Neural Information Processing Systems, Vol. 32. Curran Associates (2019). Available online at: https://proceedings.neurips.cc/paper_files/paper/2019/file/bdbca288fee7f92f2bfa9f7012727740-Paper.pdf 35. Barry

Mercer

. Exact solutions for two-dimensional time-dependent flow and deformation within a poroelastic medium. J Appl Mech. (1999) 66:536–40. 10.1115/1.2791080 36. Phillips

. Finite Element Methods in Linear Poroelasticity: Theoretical and Computational Results. Austin, TX: The University of Texas at Austin (2005).16880152 37. Daw

Wang

Perdikaris

Karpatne

. Mitigating propagation failures in pinns using evolutionary sampling, conference paper at ICLR 2023?. (Under review).

Appendix Proof of Theorem 2.1

Proof. Subtracting (6), (9) from (1), (2), respectively, we shall derive

(58)(1M+βFS)∂tepn+1-∇·K∇epn+1=-α∇·∂teun+βFS∂tepn,

(59)-∇·σ(eun+1)+α∇epn+1=0.

After differentiating (2) with respect to t, one has

(60)-∇·σ(∂teun+1)+α∇∂tepn+1=0.

We multiply (1), (3) by ∂tepn+1, ∂teun, respectively, and use integration by parts to obtain

(61)(1M+βFS)∫0T∫Ω(∂tepn+1)2dxdt+K2∫Ω(∇epn+1(T))2dx =−α∫0T∫Ω(∇·∂teun)(∂tepn+1)dxdt+βFS∫0T∫Ω(∂tepn)(∂tepn+1)dxdt,

(62)2μ∫0T∫Ωε(∂teun+1):ε(∂teun)dxdt+λ∫0T∫Ω(∇·∂teun+1)(∇·∂teun)dxdt=α∫0T∫Ω(∇·∂teun)(∂tepn+1)dxdt.

Combining Equations (4) and (5) results in the following equation.

(63)1M∫0T∫Ω(∂tepn+1)2dxdt+K2∫Ω(∇epn+1(T))2dx+2μ∫0T∫Ωε(∂teun+1):ε(∂teun)dxdt + λ∫0T∫Ω(∇·∂teun+1)(∇·∂teun)dxdt βFS∫0T∫Ω(∂tepn+1)(∂tepn+1−∂tepn)dxdt=0.

Applying the algebraic identities (46) to (6), we derive

(64)(1M+βFS2)∫0T∫Ω(∂tepn+1)2dxdt+K2∫Ω(∇epn+1(T))2dx +μ2∫0T∫Ωε(∂teun+1+∂teun):ε(∂teun+1+∂teun)dxdt +λ4∫0T∫Ω(∇·(∂teun+1+∂teun))2dxdt++βFS2∫0T∫Ω(∂tepn+1−∂tepn)2dxdt=μ2∫0T∫Ωε(∂teun+1−∂teun):ε(∂teun+1−∂teun)dxdt +λ4∫0T∫Ω(∇·(∂teun+1−∂teun))2dxdt+βFS2∫0T∫Ω(∂tepn)2dxdt.

The difference of two successive iterations based on (2) will yield

(65)-∇·σ(∂teun+1-∂teun)+α∇(∂tepn+1-∂tepn)=0.

Multiplying (8) by (∂teun+1-∂teun), integration by parts, and then applying Cauchy-Schwarz inequality, we derive that

(66)2μ∫0T∫Ωε(∂teun+1−∂teun):ε(∂teun+1−∂teun)dxdt + λ∫0T∫Ω(∇·(∂teun+1−∂teun))2dxdt+≤α22μd+λ∫0T∫Ω(∂tepn+1−∂tepn)2dxdt.

Discarding the second, third, and fourth terms in (7), we choose βFS≥α22(2μd+λ) to obtain

(67)∫0T∫Ω(∂tepn+1)2dxdt≤βFSβFS+2M∫0T∫Ω(∂tepn)2dxdt,

which directly implies (14).