Processing math: 20%
Research article

Convergence and parameter-robust analysis of several locking-free algorithms for a linear poroelasticity model

  • Received: 26 February 2025 Revised: 26 April 2025 Accepted: 07 May 2025 Published: 21 May 2025
  • MSC : 65N30, 65M12

  • The quasi-static linear poroelasticity model is widely applied in science, geophysics, biomechanics, and engineering. The simulations of this model exhibits locking phenomena and may depend on some model parameters. In this work, we follow the idea proposed in [Feng, Ge, and Li, IMA J. Numer. Anal., 2018,330–359] to transform the linear poroelasticity model into a four-field problem. The new four-field problem has a built-in mechanism to circumvent the locking phenomena existing in the original problem. We first prove that the inf-sup condition holds uniformly independent of the model parameters for the four-field problem and design several coupled, decoupled, and BDF2 algorithms in time. After that, we establish the analysis of the unconditionally optimal convergence and parameter robustness for the proposed algorithms. Finally, numerical examples are provided to investigate the convergence and parameter robustness of the proposed algorithms, which have no locking phenomena and are consistent with our theory. In addition, we also apply the algorithms to simulate brain edema, which is aligned with the experiment results.

    Citation: Peizhen Wang, Wenlong He, Yajuan Di. Convergence and parameter-robust analysis of several locking-free algorithms for a linear poroelasticity model[J]. AIMS Mathematics, 2025, 10(5): 11627-11655. doi: 10.3934/math.2025527

    Related Papers:

    [1] Abdul Qadeer Khan, Fakhra Bibi, Saud Fahad Aldosary . Bifurcation analysis and chaos in a discrete Hepatitis B virus model. AIMS Mathematics, 2024, 9(7): 19597-19625. doi: 10.3934/math.2024956
    [2] Huiling Wang, Zhaolu Tian, Yufeng Nie . The HSS splitting hierarchical identification algorithms for solving the Sylvester matrix equation. AIMS Mathematics, 2025, 10(6): 13476-13497. doi: 10.3934/math.2025605
    [3] Fatimah Alshahrani, Wahiba Bouabsa, Ibrahim M. Almanjahie, Mohammed Kadi Attouch . Robust kernel regression function with uncertain scale parameter for high dimensional ergodic data using k-nearest neighbor estimation. AIMS Mathematics, 2023, 8(6): 13000-13023. doi: 10.3934/math.2023655
    [4] Xue Wang, Weihu Cheng, Clécio S. Ferreira . Estimation and diagnostic for a skewed generalized normal partially linear models. AIMS Mathematics, 2025, 10(7): 15698-15719. doi: 10.3934/math.2025703
    [5] Lidiya Melnikova, Valeriy Rozenberg . One dynamical input reconstruction problem: tuning of solving algorithm via numerical experiments. AIMS Mathematics, 2019, 4(3): 699-713. doi: 10.3934/math.2019.3.699
    [6] Liangkun Xu, Hai Bi . A multigrid discretization scheme of discontinuous Galerkin method for the Steklov-Lamé eigenproblem. AIMS Mathematics, 2023, 8(6): 14207-14231. doi: 10.3934/math.2023727
    [7] Shuhua Wang . Convergence of online learning algorithm with a parameterized loss. AIMS Mathematics, 2022, 7(11): 20066-20084. doi: 10.3934/math.20221098
    [8] Jagbir Kaur, Vivek Sangwan . Stability estimates for singularly perturbed Fisher's equation using element-free Galerkin algorithm. AIMS Mathematics, 2022, 7(10): 19105-19125. doi: 10.3934/math.20221049
    [9] Giuseppe Maria Coclite, Lorenzo di Ruvo . Discontinuous solutions for the short-pulse master mode-locking equation. AIMS Mathematics, 2019, 4(3): 437-462. doi: 10.3934/math.2019.3.437
    [10] Adisorn Kittisopaporn, Pattrawut Chansangiam . Approximate solutions of the 2D space-time fractional diffusion equation via a gradient-descent iterative algorithm with Grünwald-Letnikov approximation. AIMS Mathematics, 2022, 7(5): 8471-8490. doi: 10.3934/math.2022472
  • The quasi-static linear poroelasticity model is widely applied in science, geophysics, biomechanics, and engineering. The simulations of this model exhibits locking phenomena and may depend on some model parameters. In this work, we follow the idea proposed in [Feng, Ge, and Li, IMA J. Numer. Anal., 2018,330–359] to transform the linear poroelasticity model into a four-field problem. The new four-field problem has a built-in mechanism to circumvent the locking phenomena existing in the original problem. We first prove that the inf-sup condition holds uniformly independent of the model parameters for the four-field problem and design several coupled, decoupled, and BDF2 algorithms in time. After that, we establish the analysis of the unconditionally optimal convergence and parameter robustness for the proposed algorithms. Finally, numerical examples are provided to investigate the convergence and parameter robustness of the proposed algorithms, which have no locking phenomena and are consistent with our theory. In addition, we also apply the algorithms to simulate brain edema, which is aligned with the experiment results.



    Poroelastic materials play an important role in various fields such as geophysics [1,2,3], biomechanics, chemical engineering [4,5], and materials science [6,7]. The physical process in porous media is usually described by a poroelasticity model. In homogeneous isotropic linear elastic porous media, the governing equations for the quasi-static poroelasticity model are given as

    2divσ(u)+αp=fin ΩT:=Ω×(0,T)Rd×(0,T), (1.1)
    (c0p+αdivu)t+divvf=ϕin ΩT, (1.2)

    where

    σ(u):=με(u)+λtr(ε(u))I,ε(u):=12(u+Tu), (1.3)
    vf:=Kμfp. (1.4)

    Here u and p represent the displacement vector of solid and the pressure of solvent, σ(u) is the (effective) stress tensor, vf is the volumetric solvent flux, and Eq (1.4) is so-called Darcy's law. I denotes the d×d identity matrix and ε(u) is the Green strain tensor. f is the body force and ϕ is the source term. The permeability tensor K=K(x) is assumed to be symmetric and uniformly positive definite in the sense that there exist positive constants K1 and K2 such that K1|ζ|2K(x)ζζK2|ζ|2 for a.e. xΩ and ζRd. The parameters μf,α and c0 are the solvent viscosity, Biot-Willis constant and constrained specific storage coefficient, respectively. λ and μ (Note that we use μ instead of 2μ in (1.3), as often seen in the literature, to simplify the notation) are Lamé constants, which can be computed from the Young's modulus E and Poisson ratio ν as

    λ=Eν(1+ν)(12ν),μ=E2(1+ν).

    Much effort has been carried out for the above poroelasticity model. Biot [1] proposed the three-dimensional consolidation theory, providing important theoretical support for the development of the poroelasticity model. Showalter [8] showed the well-posedness of a weak solution of a linear poroelasticity model. Phillips and Wheeler proposed and analyzed a continuous-in-time linear poroelasticity model [9] and a discrete-in-time linear poroelasticity model [10] with the continuous Galerkin mixed finite element method (MFEM) in space. After that, a discontinuous Galerkin MFEM for the model is also developed [11]. The nonconforming FEM was studied for the linear poroelasticity model [12,13], and the MINI and stabilized P1P1 elements were employed to solve the linear poroelasticity model [14]. Ju et al. [15] developed parameter-robust numerical algorithms for the Biot model and applied the algorithms to simulations of brain edema. Yi [16] studied two kinds of locking phenomena (Poisson locking and pressure oscillation) in linear poroelasticity and proposed a new family of MFEM to overcome them. To reveal the underlying deformation and diffusion processes of the original model, Feng et al. reformulated the linear poroelasticity model to a fluid-fluid coupled system and proposed a multiphysics MFEM to overcome the locking phenomena [17]. One can refer to more numerical methods for solving the poroelasticity model, such as the Galerkin least square methods [18,19], two-field finite element method [20,21], four-field finite element method[22,23], finite volume method [24], and Weak Galerkin finite element method [25,26].

    In practical problems, the model parameters have different ranges, such as the Poisson ratio ν=0.10.3 and the permeability K=10141016m2 in geophysics [6]. This requires us to construct numerical methods to be robust for the selection of model parameters. Following the idea of [17], we first reformulate the original problem into a four-field by introducing some auxiliary variables. After that, we discretize the four-field problem by the backward Euler method in time and design the coupled and decoupled algorithms by using MFEM with P2P1P1P1 element pairs in space. The inf-sup condition is proved to be uniformly independent of the physical parameters, and the optimal convergence rates both in space and time are also obtained for the proposed algorithms. It is pointed out that, for the decoupled algorithm [17], there exists the constraint τ=O(h2) on the mesh size h and time step τ while performing the error estimates. But there is no constraint on h and τ for the decoupled algorithm arising from our four-field problem. The reason is that the decoupled algorithm ensures that the variables p,ξ, and η are in the same time layer, thereby avoiding the use of inverse inequality. To improve computational efficiency and achieve high-accuracy in time, we further design a two-step backward differentiation formula (BDF2) Galerkin algorithm, which combines the MFEM with P2P1P1 element pairs for the space variables. The corresponding stability and error estimates of the BDF2 algorithm are performed, and the optimal convergence rate in time is also obtained in Theorem 7. Numerical tests are also provided to verify that our proposed algorithms are robust for selection of ν and K, overcome the Poisson locking in stress field and nonphysical oscillation in the pressure field, and have the optimal convergence rates. In addition, our algorithms are also used to simulate the brain edema. The numerical results on pressure are consistent with the experimental results [27]. The influence of the Poisson rate ν on the displacement of brain tissue and maximum pressure is investigated, which shows that the values of ν have little influence on the maximum pressure, and the maximum tissue deformation becomes smaller as ν gets bigger.

    The remainder of this paper is organized as follows: In Section 2, the original problem is equivalently transformed into a four-field problem, and the four-field problem is also proved to satisfy the inf-sup condition. In Section 3, the coupled and decoupled algorithms with the Euler method are proposed for the four-field problem, and the error estimates are also presented. In Section 4, the BDF2 algorithm for the poroelasticity problem is also designed, and the stability analysis and error estimates are established. In Section 5, numerical examples are provided to verify that the algorithms are robust with respect to the key physical parameters, reach the optimal convergence rates, and overcome the locking phenomenon. In Section 6, the algorithms are applied to simulate the brain edema. Finally, we draw a conclusion to summarize the main results of this paper.

    The standard function space notations are adopted; their precise definitions can be found in [28,29,30]. In particular, (,) and , denote respectively the standard L2(Ω) and L2(Ω) inner products, where ΩRd(d=2,3) is a bounded polygonal domain with the boundary Ω. For any Banach space B, we let B=[B]d, and use B to denote its dual space. We also introduce the function spaces

    L20(Ω):={qL2(Ω);(q,1)=0},X:=H1(Ω).

    Define the space of infinitesimal rigid motions by RM={r:=a+b×x;a,b,xRd}. As in [28,30,31], RM is the kernel of the operator ε, i.e., rRM if and only if ε(r)=0. Hence, we have

    ε(r)=0,divr=0,rRM. (2.1)

    Denote L2(Ω) and H1(Ω) respectively by the subspaces of L2(Ω) and H1(Ω), which are orthogonal to RM, and are given as,

    H1(Ω):={vH1(Ω);(v,r)=0rRM},L2(Ω):={gL2(Ω);g,r=0rRM}.

    It follows from [17,32] that there exists a constant ck>0 such that

    vL2(Ω)=infrRMv+rL2(Ω)rL2(Ω)ckε(v)L2(Ω),vH1(Ω). (2.2)

    To close the above system (1.1)–(1.4), here we set the following boundary and initial conditions:

    2σ(u)nαpn=f1on ΩT:=Ω×(0,T), (2.3)
    1μfKpn=ϕ1on ΩT, (2.4)
    u=u0,p=p0in Ω×{t=0}. (2.5)

    By introducing the variables q=divu,η=c0p+αq,ξ=αpλq, one has

    p=κ1ξ+κ2η,q=κ1ηκ3ξ, (2.6)

    with

    κ1=αα2+λc0,κ2=λα2+λc0,κ3=c0α2+λc0.

    Thus, the problem (1.1)–(1.4) can be equivalently written as

    2μdivε(u)+ξ=fin ΩT, (2.7)
    κ3ξ+divu=κ1ηin ΩT, (2.8)
    κ1ξ+κ2ηp=0in ΩT, (2.9)
    ηt1μfdiv(Kp)=ϕin ΩT. (2.10)

    Remark 1. We point out that the following results with respect to the pure Neumann condition can be extended to the case of the pure Dirichlet condition or mixed Dirichlet-Neumann condition. If the boundary condition (2.3) is replaced by a pure Dirichlet condition or the mixed Dirichlet-Neumann condition, there is no need to introduce the space H1(Ω) for the variable u.

    We here consider an equidistant partition of the time interval [0,T] with a time step τ=T/N and tn=nτ (n=1,2,,N), where NN. The focus of this section is on the proof of the inf-sup condition, which will hold uniformly independently of physical parameters. To do that, we first discretize the time-dependent problem (2.7)–(2.10) with (2.3)–(2.5) by the backward Euler method and BDF2 scheme in time to make it a static problem in each time step. Note that the discussion on inf-sup condition is similar for the backward Euler method and BDF2 scheme. Below, we mainly consider the backward Euler method. By the backward Euler method, the weak formulation of problem (2.7)–(2.10) with (2.3)–(2.5) reads: (v,φ,ω,ψ)H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω), (u,ξ,η,p)H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω) such that

    2μ(ε(u),ε(v))(ξ,divv)=(f,v)+f1,v, (2.11)
    κ3(ξ,φ)+(divu,φ)=κ1(η,φ), (2.12)
    (κ1ξ,ω)+(κ2η,ω)(p,ω)=0, (2.13)
    (η,ψ)+τμf(Kp,ψ)=(˜ϕ,ψ)+˜ϕ1,ψ, (2.14)

    where ˜ϕ(x,tn)=τϕ(x,tn)+η(x,tn1) and ˜ϕ1=τϕ1.

    Define

    A(u,ξ,η,p)=(μdivε00divκ3κ100κ1κ21001τμfdivK)(uξηp).

    Taking Λ=H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω) and ξ0 as the mean-value zero part of ξ, we define

    (u,ξ,η,p)Λ=μ12ε(u)L2(Ω)+κ123ξL2(Ω)+δ120ξ0L2(Ω)+κ122ηL2(Ω)+(τμf)12K12pL2(Ω), (2.15)

    where δ0 is a positive constant, uH1(Ω), ξL2(Ω), ηL2(Ω) and pH1(Ω).

    The following theorem shows that A is invertible and the inverse is a map from Λ to Λ.

    Theorem 1. There exists a constant β>0 for the problem (2.11)(2.14), independent of μ,κ2,κ3 and Kτ/μf, such that the following inf-sup condition holds:

    inf(u,ξ,η,p)Λsup(v,φ,ω,ψ)Λ(A(u,ξ,η,p),(v,φ,ω,ψ))(u,ξ,η,p)Λ(v,ϕ,ω,ψ)Λβ. (2.16)

    Proof. To prove (2.16), we need to prove the following inequalities: (0,0,0,0)(u,ξ,η,p)Λ, (v,φ,ω,ψ)Λ such that

    (v,φ,ω,ψ)Λβ1(u,ξ,η,p)Λ,(A(u,ξ,η,p),(v,φ,ω,ψ))β2(u,ξ,η,p)2Λ,

    where β1 and β2 are positive constants and independent of μ,κ1,κ2,κ3 and Kτ/μf.

    Suppose that (0,0,0,0)(u,ξ,η,p)Λ is given. It follows from the theory of the Stokes equation that there exists a constant β0>0, depending only on the domain Ω, and wH1(Ω) such that

    (divw,ξ)=ξ02L2(Ω),(ε(w),ε(w))β20ξ02L2(Ω). (2.17)

    Taking v=uδ0w, ϕ=ξ, ω=η, ψ=p in (A(u,ξ,η,p),(v,φ,ω,ψ)), we have

    (A(u,ξ,η,p),(v,φ,ω,ψ))=μ(ε(u),ε(u))δ0μ(ε(u),ε(w))(ξ,divu)+δ0(ξ,divw)+κ3(ξ,ξ)+(divu,ξ)κ1(η,ξ)+(κ1ξ+κ2η,η)(p,η)+(η,p)+τμfKp2L2(Ω)=με(u)2L2(Ω)δ0μ(ε(u),ε(w))+δ0(ξ,divw)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω). (2.18)

    Using Young's inequality and (2.17), we have

    (ε(u),ε(w))12θ0(ε(u),ε(u))+θ02(ε(w),ε(w))12θ0(ε(u),ε(u))+θ0β202(ξ0,ξ0),θ0>0. (2.19)

    Inserting (2.19) into (2.18), using (2.17) and choosing θ0=δ0=β20μ1, we obtain

    (A(u,ξ,η,p),(v,φ,ω,ψ))μ2ε(u)2L2(Ω)+δ02ξ02L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω). (2.20)

    Noting that

    (u,ξ,η,p)2Λ=(μ12ε(u)L2(Ω)+κ123ξL2(Ω)+κ122ηL2(Ω)+(τμf)12KpL2(Ω))24(με(u)2L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω))4(με(u)2L2(Ω)+δ0ξ02L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω)), (2.21)

    and using (2.21) in (2.20), we have

    (A(u,ξ,η,p),(v,φ,ω,ψ))18(u,ξ,η,p)2Λ. (2.22)

    It is easy to check that

    (v,φ,ω,ψ)Λ(1+μ12δ120β0)(u,ξ,η,p)Λ. (2.23)

    Combining (2.22) and (2.23), the (2.16) holds and β=1/(8(1+μ12δ120β0)). The proof is complete.

    Theorem 1 shows that the bilinear form (A(),) is coercive. Moreover, it is easy to check that (A(),) is bounded. According to the Lax-Milgram theorem, there exists a unique solution (u,ξ,η,p)Λ for the linear system (2.11)–(2.14). Thus, there holds the following theorem:

    Theorem 2. There exists a unique weak solution to problem (2.7)(2.10) with (2.3)(2.5). Likewise, there exists a unique weak solution to problem (1.1)(1.4) and (2.3)(2.5).

    Let Th be a quasi-uniform triangulation or rectangular partition of Ω with maximum mesh size h, and ˉΩ=KThˉK. In this part, we use the backward Euler method in time discretization. For a smooth function v defined on [0,T], we denote dtvn=(vnvn1)/τ.

    A number of stable mixed finite element spaces (Xh,Mh) have been presented in the literature [33]. A typical example is the following so-called Taylor-Hood element (cf. [33,34]):

    Xh={vhC0(¯Ω);vh|KP2(K)KTh},Mh={φhC0(¯Ω);φh|KP1(K)KTh}.

    Finite element approximation space Wh for p and η variables can be chosen independently; any piecewise polynomial space is acceptable provided that WhMh. A convenient choice is Wh=Mh.

    Define

    Vh:={vhXh;(vh,r)=0,rRM}. (3.1)

    If the problem (2.11)–(2.14) is discretized by the finite element spaces Vh,Mh, and Wh, then the discrete counterpart is to find (uh,ξh,ηh,ph)Vh×Mh×Wh×Wh such that

    2μ(ε(uh),ε(vh))(ξh,divvh)=(f,vh)+f1,vh,vhVh, (3.2)
    κ3(ξh,φh)+(divuh,φh)=κ1(ηh,φh),φhMh, (3.3)
    κ1(ξh,ωh)+κ2(ηh,ωh)(p,ωh)=0,ωhWh, (3.4)
    (ηh,ψh)+τμf(Kph,ψh)=(˜ϕ,ψh)+˜ϕ1,ψh,ψhWh. (3.5)

    Set Λh=Vh×Mh×Wh×Wh, and define the operator Ah:ΛhΛh arisen from (3.2)–(3.5). Then we have the following theorem:

    Theorem 3. For the problem (3.2)(3.5), there exists a constant β>0, independent of μ,κ2,κ3, and Kτ/μf, such that the following inf-sup condition holds:

    inf(uh,ξh,ηh,ph)Λhsup(vh,φh,ωh,ψh)Λh(Ah(uh,ξh,ηh,ph),(vh,φh,ωh,ψh))(uh,ξh,ηh,ph)Λh(vh,φh,ωh,ψh)Λhβ. (3.6)

    We omit the proof of Theorem 3, since it is completely analogous to the proof of Theorem 1. Moreover, Theorem 3 implies that Ah is invertible and the inverse is a map from Λh to Λh with operator norm independent of the parameters.

    Algorithm 1 A coupled algorithm
    Input: Evaluate u0hVh,p0hWh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      for n=0,1,2,, do
        solve for (un+1h,ξn+1h,ηn+1h,pn+1h)Vh×Mh×Wh×Wh such that
                μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηn+1h,φh),φhMh,κ1(ξn+1h,ωh)+κ2(ηn+1h,ωh)(pn+1h,ωh)=0,ωhWh,(dtηn+1h,ψh)+1μf(Kpn+1h,ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh,
      end for

    Algorithm 2 A decoupled algorithm
    Input: Evaluate u0hVh,p0hWh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      for n=0,1,2,, do
        (i) solve for (un+1h,ξn+1h)Vh×Mh such that
                μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηnh,φh),φhMh.
        (ii) Using (un+1h,ξn+1h) obtained in (i), solve for (ηn+1h,pn+1h)Wh×Wh by
                κ1(ξn+1h,ωh)+κ2(ηn+1h,ωh)(pn+1h,ωh)=0,ωhWh,(dtηn+1h,ψh)+1μf(Kpn+1h,ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh,
      end for

    Below, we assume \mathbf{f}, \; \mathbf{f}_{1} , \phi , and \phi_{1} all are independent of t . Note that all the results can be easily extended to the case of time-dependent source functions.

    Now, we consider the fully discrete MFEM for the problem (2.7)–(2.10) with (2.3)–(2.5) based on the backward Euler method, which involves the coupled (i.e., Algorithm 1) and decoupled (i.e., Algorithm 2) algorithms. The coupled algorithm solves four unknowns simultaneously based on problem (2.7)–(2.10). The decoupled algorithm only solves two subproblems. Noting that problem (2.7) and (2.8) can be regarded as a generalized Stokes problem for a given \eta and problem (2.9) and (2.10) is a coupled problem associated with the diffusion problem, we can solve a generalized Stokes problem for \mathbf{u} and \xi by using the solution of \eta at the previous time-step. After that, we solve the diffusion problem for \eta and p .

    To derive the error estimates of the MFEM, we define L^{2}(\Omega) -projection \mathcal{Q}_{h}: L^{2}(\Omega)\rightarrow X^{k}_{h} by

    \begin{equation} (\mathcal{Q}_{h}\varphi,\psi_{h}) = (\varphi,\psi_{h}), \quad \forall \; \psi_{h}\in X^{k}_{h}\text{ and } \varphi \in L^{2}(\Omega), \end{equation} (3.7)

    where X^{k}_{h}: = \{\psi_{h}\in C^{0};\; \psi_{h}|_{\mathcal{K}}\in P_{k}(\mathcal{K})\; \forall \mathcal{K}\in\mathcal{T}_{h}\} , k is the degree of piecewise polynomial on \mathcal{K} .

    Next, for any \varphi\in H^{1}(\Omega) , we define the elliptic projection \mathcal{S}_{h}: H^{1}(\Omega)\rightarrow X^{k}_{h} as

    \begin{align} (K\nabla \mathcal{S}_{h}\varphi,&\nabla\varphi_{h}) = (K\nabla\varphi,\nabla\varphi_{h}), \quad \forall \; \varphi_{h}\in X^{k}_{h}, \end{align} (3.8)
    \begin{align} &(\mathcal{S}_{h}\varphi,1) = (\varphi,1), \end{align} (3.9)

    and for any \textbf{v}\in\textbf{H}^{1}(\Omega) , we define its elliptic projection \mathcal{R}_{h}: \textbf{H}^{1}(\Omega)\rightarrow \textbf{V}^{k}_{h} as

    \begin{eqnarray} (\varepsilon(\mathcal{R}_{h}\textbf{v}),\varepsilon(\textbf{w}_{h})) = (\varepsilon(\textbf{v}),\varepsilon(\textbf{w}_{h}))\; \; \; \; \; \forall\textbf{w}_{h}\in \textbf{V}^{k}_{h}, \end{eqnarray} (3.10)

    where \textbf{V}^k_{h}: = \{\textbf{v}_{h} \in \textbf{C}^{0};\; \textbf{v}_{h}|_{\mathcal{K}}\in \textbf{P}_{k}(\mathcal{K}), (\textbf{v}_{h}, \textbf{r}) = 0 \; \forall \textbf{r} \in \textbf{RM}\} .

    It follows from [28] that for any 0\leq s\leq k , \mathcal{Q}_{h}, \mathcal{S}_{h} and \mathcal{R}_{h} satisfy

    \begin{align} \|\mathcal{Q}_{h}\varphi-\varphi\|_{L^{2}(\Omega)}+&h\|\nabla(\mathcal{Q}_{h}\varphi-\varphi)\|_{L^{2}(\Omega)}\leq Ch^{s+1}\|\varphi\|_{H^{s+1}(\Omega)},\; \; \forall \varphi\in H^{s+1}(\Omega), \end{align} (3.11)
    \begin{align} \|\mathcal{S}_{h}\varphi-\varphi\|_{L^{2}(\Omega)}+&h\|\nabla(\mathcal{S}_{h}\varphi-\varphi)\|_{L^{2}(\Omega)}\leq Ch^{s+1}\|\varphi\|_{H^{s+1}(\Omega)},\; \; \forall \varphi\in H^{s+1}(\Omega), \end{align} (3.12)
    \begin{align} \|\mathcal{R}_{h}\textbf{v}-\textbf{v}\|_{L^{2}(\Omega)}+&h\|\nabla(\mathcal{R}_{h}\textbf{v}-\textbf{v})\|_{L^{2}(\Omega)}\leq Ch^{s+1}\|\textbf{v}\|_{H^{s+1}(\Omega)},\; \; \forall \textbf{v}\in \textbf{H}^{s+1}(\Omega). \end{align} (3.13)

    Furthermore, we introduce the following error notations:

    \begin{align*} E_{{\mathbf{u}}}^{n} = {\mathbf{u}}(t_{n})-{\mathbf{u}}_{h}^{n},\; \; \; \; E_{\xi}^{n} = \xi(t_{n})-\xi_{h}^{n},\; \; \; \; E_{\eta}^{n} = \eta(t_{n})-\eta_{h}^{n},\; \; \; \; E_{p}^{n} = p(t_{n})-p_{h}^{n}. \end{align*}

    Also, we denote

    \begin{align*} &E_{{\mathbf{u}}}^{n} = {\mathbf{u}}(t_{n})-\mathcal{R}_{h}({\mathbf{u}}(t_{n}))+\mathcal{R}_{h}({\mathbf{u}}(t_{n}))-{\mathbf{u}}_{h}^{n}: = Y_{{\mathbf{u}}}^{n}+Z_{{\mathbf{u}}}^{n},\\ &E_{i}^{n} = i(t_{n})-\mathcal{S}_{h}(i(t_{n}))+\mathcal{S}_{h}(i(t_{n}))-i_{h}^{n}: = Y_{i}^{n}+Z_{i}^{n},\; i = \xi,\; \eta,\; p,\\ &E_{i}^{n} = i(t_{n})-\mathcal{Q}_{h}(i(t_{n}))+\mathcal{Q}_{h}(i(t_{n}))-i_{h}^{n}: = F_{i}^{n}+G_{i}^{n},\; i = \xi,\; \eta,\; p. \end{align*}

    Then we present the error estimates for the coupled and decoupled algorithms by the following theorems.

    Theorem 4. Let \left\lbrace ({\mathbf{u}}_{h}^{n}, \xi_{h}^{n}, \eta_{h}^{n}, p_{h}^{n})\right\rbrace_{n\geq0} be solutions of the coupled algorithm (i.e., Algorithm 1). Then it holds

    \begin{eqnarray} &&\max\limits_{0\leq n\leq N}\left[\sqrt{\mu}\|\nabla({\mathbf{u}}(t_{n})-{\mathbf{u}}_{h}^{n})\|_{L^{2}(\Omega )}+\sqrt\kappa_{2}\|\eta(t_{n})-\eta_{h}^{n}\|_{L^{2}(\Omega )}+\sqrt{\kappa_{3}}\|\xi(t_{n})-\xi_{h}^{n}\|_{L^{2}(\Omega )} \right]\\ &&+\Big(\tau\sum\limits_{n = 0}^{N}\dfrac{1}{\mu_{f}}\|K^{\frac{1}{2}}( p(t_{n})-p_{h}^{n})\|_{L^{2}(\Omega )}^{2} \Big)^{\frac{1}{2}}\leq C_{1}\tau+C_{2}h^{2}. \end{eqnarray} (3.14)

    Theorem 5. Let \left\lbrace ({\mathbf{u}}_{h}^{n}, \xi_{h}^{n}, \eta_{h}^{n}, p_{h}^{n})\right\rbrace_{n\geq0} be solutions of the decoupled algorithm (i.e., Algorithm 2). Then for the case c_{0} > 0 and \lambda < \infty , it holds

    \begin{eqnarray} &&\max\limits_{0\leq n\leq N}\left[\sqrt{\mu}\|\nabla({\mathbf{u}}(t_{n})-{\mathbf{u}}_{h}^{n})\|_{L^{2}(\Omega )}+\sqrt{\kappa_{2}}\|\eta(t_{n})-\eta_{h}^{n}\|_{L^{2}(\Omega )}+\sqrt{\kappa_{3}}\|\xi(t_{n})-\xi_{h}^{n}\|_{L^{2}(\Omega )} \right]\\ &&+\Big(\tau\sum\limits_{n = 0}^{N}\dfrac{1}{\mu_{f}}\|K^{\frac{1}{2}}( p(t_{n})-p_{h}^{n})\|_{L^{2}(\Omega )}^{2} \Big)^{\frac{1}{2}}\leq \bar{C}_{1}\tau+\bar{C}_{2}h^{2}. \end{eqnarray} (3.15)

    The proof of the aforementioned theorems follows a similar approach to that of [17, Theorem 3.7]; thus, further details are omitted here. Our primary focus is on the constraint condition \tau = O(h^2) in the error estimates for the algorithm with \theta = 0 [17]. This constraint arises because the choice of the test function \psi_{h} results in the variables p, \; \xi , and \eta being in different time layers, which in turn leads to the application of an inverse inequality in the estimation of the right-hand side term. To address this issue, a straightforward solution is to select an appropriate test function \psi_{h} that ensures the variables p, \; \xi , and \eta reside in the same time layer. To do that, we first examine the following error equations:

    \begin{align} {2} &\kappa_{1}\left(G_\xi^{n+1},\omega_{h}\right)+\kappa_{2}\left(G_\eta^{n+1},\omega_{h}\right) -\left(G_p^{n+1},\omega_{h}\right) = 0, \quad &&\forall \omega_{h}\in W_{h}, \end{align} (3.16)
    \begin{align} &\bigl(d_t G_\eta^{n+1}, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K{{\nabla}} Z_p^{n+1},{{\nabla}}\psi_h \bigr) = (R_{1,\eta}^{n+1},\psi_{h}), \quad && \forall \psi_h\in W_h, \end{align} (3.17)

    which are derived from the fourth equation of Algorithm 2 and the continuous problem (2.10) the fourth equation of Algorithm 2 and the continuous problem (2.10). Based on Algorithm 2, we can define \eta_{h}^{-1} by

    \begin{align} \kappa_{1}(\eta_{h}^{-1},\varphi_{h}) = \kappa_{3}(\xi_{h}^{0},\varphi_{h})+({\rm div}\mathbf{u}_{h}^{0},\varphi_{h}). \end{align}

    Setting \omega_{h} = d_{t}G_{\eta}^{n}, \; \psi_{h} = Z_{p}^{n} in (3.16) and (3.17) separately after lowering the super index n+1 to n , applying the summation operator \tau\sum_{n = 0}^{l} to both sides of (3.17) and utilizing the fact that G_{\eta}^{-1} = 0 , there holds

    \begin{eqnarray} &&\dfrac{\kappa_{2}}{2}\|G_{\eta}^{l}\|_{L^{2}(\Omega)}^{2}+\tau\sum\limits_{n = 0}^{l}\left[\kappa_{1}(d_{t}G_{\eta}^{n},G_{\xi}^{n+1})+\dfrac{\kappa_{2}\tau}{2}\|d_{t}G_{\eta}^{n}\|_{L^{2}(\Omega)}^{2}+\frac{1}{\mu_f}\|K^{\frac{1}{2}}{{\nabla}} Z_p^{n}\|_{L^{2}(\Omega)}^{2}\right] \\ && = \tau\sum\limits_{n = 0}^{l}\left[\bigl(d_t G_\eta^{n},Y_{p}^{n}-F_{p}^{n}\bigr)+(d_{t}\eta(t_{n+1})-\eta_{t}(t_{n+1}),Z_{p}^{n})+\kappa_{1}\tau(d_{t}G_{\xi}^{n+1},d_{t}G_{\eta}^{n})\right]. \end{eqnarray} (3.18)

    Note that the choice of \psi_{h} = Z_{p}^{n} is critical. Unlike the setting in [17], where \psi_{h} = \hat{Z}_{p}^{n}: = \kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n} is employed, selecting \psi_{h} = Z_{p}^{n} eliminates the need to apply the inverse inequality to bound a gradient term (as seen in [17, estimate (3.52)]). Instead, we only need to apply the Cauchy-Schwarz and Young inequalities to the last term on the right-hand side of (3.18) to derive

    \begin{eqnarray} \kappa_{1}\tau^{2}\sum\limits_{n = 0}^{l}(d_{t}G_{\xi}^{n+1},d_{t}G_{\eta}^{n})\leq \kappa_{1}\tau^{2}\sum\limits_{n = 0}^{l}\left[\dfrac{1}{2\epsilon}\|d_{t}G_{\xi}^{n+1}\|_{L^{2}(\Omega)}^{2}+\dfrac{\epsilon}{2}\|d_{t}G_{\eta}^{n})\|_{L^{2}(\Omega)}^{2}\right], \end{eqnarray} (3.19)

    where \epsilon\in\big(\sqrt{\frac{2\kappa_{2}}{\kappa_{3}}}, +\infty\big) is a positive constant from Young inequality.

    Additionally, we point out that:

    (i) The aim of \epsilon\in\big(\sqrt{\frac{2\kappa_{2}}{\kappa_{3}}}, +\infty\big) is to ensure the signs of terms \frac{\kappa_{2}-\kappa_{1}\epsilon}{2}\tau^{2} \sum_{n = 0}^{l}\|d_{t}G_{\eta}^{n}\|_{L^{2}(\Omega)}^{2} and \big(\frac{\kappa_{3}}{4}-\frac{\kappa_{1}}{2\epsilon}\big)\tau^{2}\sum_{n = 0}^{l}\|d_{t}G_{\eta}^{n}\|_{L^{2}(\Omega)}^{2} to be positive. In addition, for the case of \lambda\rightarrow +\infty , there does not exist the term \kappa_{1}\tau^{2}\sum_{n = 0}^{l}(d_{t}G_{\xi}^{n+1}, d_{t}G_{\eta}^{n}) in error estimates. For the case of c_{0} = 0 , it has a similar discussion for the term \kappa_{1}\tau^{2}\sum_{n = 0}^{l}(d_{t}G_{\xi}^{n+1}, d_{t}G_{\eta}^{n}) .

    (ii) For the decoupled algorithm, it does not exist the constraint on the time step \tau and mesh size h .

    We now consider the BDF2 scheme in time discretization for the problem (2.7)–(2.10) with (2.3)–(2.5). For a smooth function v defined on [0, T] , we denote

    \begin{equation} v^{n+\frac{1}{2}}: = (v^{n+1}+v^{n})/2 \quad \text{and} \quad D_{t}v^{n} = \frac{3v^{n}-4v^{n-1}+v^{n-2}}{2\tau}. \end{equation} (4.1)

    Based on (4.1), we introduce the BDF2 fully discrete scheme in Algorithm 3.

    Algorithm 3 BDF2 algorithm
    Input: Evaluate \mathbf{u}^0_h\in \mathbf{V}_h, \; p_{h}^{0}\in M_{h} , \xi^0_h\in M_h by \xi^0_h =\alpha p_{h}^{0}-\lambda{\rm div}\mathbf{u}_{h}^{0} and \eta_{h}^{0}\in W_{h} by \eta_{h}^{0}=c_{0}p_{h}^{0}+\alpha{\rm div}\mathbf{u}_{h}^{0} .
      Step 1: Solve for (\mathbf{u}^{1}_h, \xi^{1}_h, \eta^{1}_h)\in \mathbf{V}_h\times M_h \times W_h such that
             \begin{align} {2} & \; \; \mu\bigl(\varepsilon(\mathbf{u}_{h}^{\frac{1}{2}}), \varepsilon(\mathbf{v}_h) \bigr)-\bigl(\xi^{\frac{1}{2}}_h, {\rm div} \mathbf{v}_h \bigr) = (\mathbf{f}, \mathbf{v}_h)+\langle \mathbf{f}_1, \mathbf{v}_h\rangle, & & \quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align}\quad\quad\quad\quad\quad (4.2)
             \begin{align} & \; \; \kappa_3\bigl(\xi^{\frac{1}{2}}_h, \varphi_h \bigr) +\bigl({\rm div}\mathbf{u}^{\frac{1}{2}}_h, \varphi_h \bigr) =\kappa_1\bigl(\eta^{\frac{1}{2}}_h, \varphi_h \bigr), & & \quad \forall \varphi_h \in M_h, \end{align} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad(4.3)
             \begin{align} & \; \; \bigl(d_t\eta^{1}_h, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}\xi_{h}^{\frac{1}{2}}+\kappa_{2}\eta_{h}^{\frac{1}{2}}), \nabla\psi_h \bigr)=(\phi, \psi_h)+\langle \phi_1, \psi_h\rangle, & & \quad \forall \psi_h\in W_h. \end{align}\; (4.4)
      Step 2:
      for n=1, 2, \cdots , do
        (i) Solve for (\mathbf{u}^{n+1}_h, \xi^{n+1}_h, \eta^{n+1}_h)\in \mathbf{V}_h\times M_h \times W_h such that
             \begin{align} {2} & \; \; \mu\bigl(\varepsilon(\mathbf{u}_{h}^{n+1}), \varepsilon(\mathbf{v}_h) \bigr)-\bigl(\xi^{n+1}_h, {\rm div} \mathbf{v}_h \bigr) = (\mathbf{f}, \mathbf{v}_h)+\langle \mathbf{f}_1, \mathbf{v}_h\rangle, & & \quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align}\quad\quad\quad\quad\quad\quad (4.5)
             \begin{align} & \; \; \kappa_3\bigl(\xi^{n+1}_h, \varphi_h \bigr) +\bigl({\rm div}\mathbf{u}^{n+1}_h, \varphi_h \bigr) =\kappa_1\bigl(\eta^{n+1}_h, \varphi_h \bigr), & & \quad \forall \varphi_h \in M_h, \end{align} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad(4.6)
             \begin{align} & \; \; \bigl(D_t\eta^{n+1}_h, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}\xi_{h}^{n+1}+\kappa_{2}\eta_{h}^{n+1}), \nabla\psi_h \bigr)=(\phi, \psi_h)+\langle \phi_1, \psi_h\rangle, & & \quad \forall \psi_h\in W_h. \end{align} (4.7)
        (ii) Update p^{n+1}_h and q^{n+1}_h by
             \begin{align} &p^{n+1}_h=\kappa_1\xi^{n+1}_h +\kappa_2\eta^{n+1}_h, \quad q^{n+1}_h=\kappa_1\eta^{n+1}_h-\kappa_3\xi^{n+1}_h, \end{align}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (4.8)
      end for

    The next lemma establishes a discrete energy inequality for the BDF2 algorithm.

    Lemma 4.1. Let \{(\mathbf{u}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0} be defined by the BDF2 algorithm. Then it holds

    \begin{align} &\frac{\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|\xi_{h}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|\eta_{h}^{l+1}\|_{L^{2}(\Omega)}^{2} \\ &\quad+\frac{\tau}{2\mu_{f}}\sum\limits_{n = 1}^{l}\|K^{\frac{1}{2}}\nabla (\kappa_{1}\xi_{h}^{n+1}+\kappa_{2}\eta_{h}^{n+1})\|_{L^{2}(\Omega)}^{2}\\ &\leq C_{3}\big(\|\mathbf{f}\|_{L^{2}(\Omega)}^{2}+\|\mathbf{f}_{1}\|_{L^{2}(\partial\Omega)}^{2}+\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2}+\|\xi_{h}^{0}\|_{L^{2}(\Omega)}^{2}+\|\eta_{h}^{0}\|_{L^{2}(\Omega)}^{2}\\ &\quad+\|\phi\|_{L^{2}(\Omega)}^{2}+\|\phi_{1}\|_{L^{2}(\partial\Omega)}^{2}\big),\quad \forall\; l\geq1, \end{align} (4.9)

    where C_{3} = C_{3}(c_{k}, c_{a}, c_{b}, c_{p}, \mu, \kappa_{3}, \kappa_{2}, T, K, \mu_{f}) .

    Proof. Setting \mathbf{v}_{h} = d_{t}\mathbf{u}_{h}^{1} in (4.2), \varphi_{h} = d_{t}\xi_{h}^{1} in (4.3) and \psi_{h} = p_{h}^{\frac{1}{2}} in (4.4), we get

    \begin{align} {2} &\; \; \mu\bigl(\varepsilon(\mathbf{u}_{h}^{\frac{1}{2}}), \varepsilon(d_{t}\mathbf{u}_{h}^{1}) \bigr)-\bigl( \xi^{\frac{1}{2}}_h, {\rm div} d_{t}\mathbf{u}_{h}^{1} \bigr) = (\mathbf{f}, d_{t}\mathbf{u}_{h}^{1})+\langle \mathbf{f}_1,d_{t}\mathbf{u}_{h}^{1}\rangle, && \end{align} (4.10)
    \begin{align} &\; \; \kappa_3\bigl(\xi^{\frac{1}{2}}_h, d_{t}\xi_{h}^{1} \bigr) +\bigl({\rm div}\mathbf{u}^{\frac{1}{2}}_h, d_{t}\xi_{h}^{1}\bigr) = \kappa_1\bigl(\eta^{\frac{1}{2}}_h, d_{t}\xi_{h}^{1} \bigr), && \end{align} (4.11)
    \begin{align} &\; \; \bigl(d_t\eta^{1}_h, p_{h}^{\frac{1}{2}} \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}\xi_{h}^{\frac{1}{2}}+\kappa_{2}\eta_{h}^{\frac{1}{2}}),\nabla p_{h}^{\frac{1}{2}} \bigr) = (\phi, p_{h}^{\frac{1}{2}})+\langle \phi_1,p_{h}^{\frac{1}{2}} \rangle. && \end{align} (4.12)

    Adding (4.10)–(4.12), there holds

    \begin{align} &\frac{\mu}{2\tau}\big(\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}-\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{2\tau}\big(\|\xi_{h}^{1}\|_{L^{2}(\Omega)}^{2}-\|\xi_{h}^{0}\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{2}}{2\tau}\big(\|\eta_{h}^{1}\|_{L^{2}(\Omega)}^{2}-\|\eta_{h}^{0}\|_{L^{2}(\Omega)}^{2}\big) \\ &+\frac{1}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}\xi_{h}^{\frac{1}{2}}+\kappa_{2}\eta_{h}^{\frac{1}{2}})\|_{L^{2}(\Omega)}^{2} = (\mathbf{f}, d_{t}\mathbf{u}_{h}^{1})+\langle \mathbf{f}_1,d_{t}\mathbf{u}_{h}^{1}\rangle+(\phi, p_{h}^{\frac{1}{2}})+\langle \phi_1,p_{h}^{\frac{1}{2}}\rangle+\frac{1}{\tau}\Big((\xi_{h}^{0},{\rm div}\mathbf{u}_{h}^{1})\\ &-(\xi_{h}^{1},{\rm div}\mathbf{u}_{h}^{0})\Big)+\frac{\kappa_{1}}{\tau}\Big((\eta_{h}^{0},\xi_{h}^{1})-(\eta_{h}^{1},\xi_{h}^{0})\Big). \end{align} (4.13)

    Using the Cauchy-Schwarz inequality, triangle inequality, (2.2), (4.32), trace inequality, Poincaré inequality, and Young inequality, we can bound the right-side terms of (4.13) as follows:

    \begin{align} (\mathbf{f}, d_{t}\mathbf{u}_{h}^{1}) &\leq\frac{1}{\tau}\big(\frac{5c_{k}^2}{\mu}\|\mathbf{f}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{16}\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{4}\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2}\big), \end{align} (4.14)
    \begin{align} \langle \mathbf{f}_1,d_{t}\mathbf{u}_{h}^{1}\rangle &\leq\frac{1}{\tau}\big(\frac{5c_{a}^{2}c_{k}^2}{\mu}\|\mathbf{f}_{1}\|_{L^{2}(\partial\Omega)}^{2}+\frac{\mu}{16}\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{4}\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2}\big), \end{align} (4.15)
    \begin{align} (\phi, p_{h}^{\frac{1}{2}}) &\leq c_{p}^{2}\mu_{f}\|K^{-\frac{1}{2}}\phi\|_{L^{2}(\Omega)}^{2}+ \frac{1}{4\mu_{f}}\|K^{\frac{1}{2}}\nabla p_{h}^{\frac{1}{2}}\|_{L^{2}(\Omega)}^{2}, \end{align} (4.16)
    \begin{align} \langle \phi_1,p_{h}^{\frac{1}{2}}\rangle &\leq c_{a}^{2}c_{p}^{2}\mu_{f}\|K^{-\frac{1}{2}}\phi_{1}\|_{L^{2}(\partial\Omega)}^{2}+\frac{1}{4\mu_{f}} \|K^{\frac{1}{2}}\nabla p_{h}^{\frac{1}{2}}\|_{L^{2}(\Omega)}^{2}, \end{align} (4.17)
    \begin{align} (\xi_{h}^{0},{\rm div}\mathbf{u}_{h}^{1})-(\xi_{h}^{1},{\rm div}\mathbf{u}_{h}^{0})&\leq\frac{2c_{b}^{2}}{\mu}\|\xi_{h}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{8}\|\xi_{h}^{1}\|_{L^{2}(\Omega)}^{2} \end{align} (4.18)
    \begin{align} &\quad+\frac{2c_{b}^{2}}{\kappa_{3}}\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2},\\ \kappa_{1}(\eta_{h}^{0},\xi_{h}^{1})-\kappa_{1}(\eta_{h}^{1},\xi_{h}^{0})&\leq\frac{2\kappa_{1}^{2}}{\kappa_{3}}\|\eta_{h}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{8}\|\xi_{h}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|\eta_{h}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{1}^{2}}{\kappa_{2}}\|\xi_{h}^{0}\|_{L^{2}(\Omega)}^{2}, \end{align} (4.19)

    where c_{a} is a positive constant related to Trace inequality.

    Substituting (4.14)–(4.19) into (4.13), we have

    \begin{align} &\frac{\mu}{4}\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|\xi_{h}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|\eta_{h}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\tau}{2\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}\xi_{h}^{\frac{1}{2}}+\kappa_{2}\eta_{h}^{\frac{1}{2}})\|_{L^{2}(\Omega)}^{2} \\ &\leq\frac{5c_{k}^2}{\mu}\|\mathbf{f}\|_{L^{2}(\Omega)}^{2}+\frac{5c_{a}^{2}c_{k}^2}{\mu}\|\mathbf{f}_{1}\|_{L^{2}(\partial\Omega)}^{2}+c_{p}^{2}\mu_{f}\tau\|K^{-\frac{1}{2}}\phi\|_{L^{2}(\Omega)}^{2}+c_{a}^{2}c_{p}^{2}\mu_{f}\tau\|K^{-\frac{1}{2}}\phi_{1}\|_{L^{2}(\partial\Omega)}^{2}\\ &\quad+\big(\mu+\frac{2c_{b}^{2}}{\kappa_{3}}\big)\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2}+\big(\frac{\kappa_{3}}{2}+\frac{2c_{b}^{2}}{\mu}+\frac{\kappa_{1}^{2}}{\kappa_{2}}\big)\|\xi_{h}^{0}\|_{L^{2}(\Omega)}^{2}+\big(\frac{\kappa_{2}}{2}+\frac{2\kappa_{1}^{2}}{\kappa_{3}}\big)\|\eta_{h}^{0}\|_{L^{2}(\Omega)}^{2}. \end{align} (4.20)

    Next, setting \mathbf{v}_{h} = D_{t}\mathbf{u}_{h}^{n+1}, \; \varphi_{h} = \xi_{h}^{n+1} and \psi_{h} = p_{h}^{n+1} in (4.5)–(4.7), we obtain

    \begin{align} {2} &\; \; \mu\bigl(\varepsilon(\mathbf{u}_{h}^{n+1}), \varepsilon(D_{t}\mathbf{u}_{h}^{n+1}) \bigr)-\bigl( \xi^{n+1}_h, {\rm div} D_{t}\mathbf{u}_{h}^{n+1} \bigr) = (\mathbf{f}, D_{t}\mathbf{u}_{h}^{n+1})+\langle \mathbf{f}_1,D_{t}\mathbf{u}_{h}^{n+1}\rangle, && \end{align} (4.21)
    \begin{align} &\; \; \kappa_3\bigl(D_{t}\xi^{n+1}_h, \xi_{h}^{n+1} \bigr) +\bigl({\rm div}D_{t}\mathbf{u}^{n+1}_h, \xi_{h}^{n+1} \bigr) = \kappa_1\bigl( D_{t}\eta^{n+1}_h, \xi_{h}^{n+1} \bigr), && \end{align} (4.22)
    \begin{align} &\; \; \bigl(D_t\eta^{n+1}_h, p_{h}^{n+1} \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}\xi_{h}^{n+1}+\kappa_{2}\eta_{h}^{n+1}),\nabla p_{h}^{n+1} \bigr) = (\phi, p_{h}^{n+1})+\langle \phi_1,p_{h}^{n+1}\rangle. && \end{align} (4.23)

    Adding (4.21)–(4.23) and using the fact of

    \begin{align} (D_{t}v_{h}^{n+1},v_{h}^{n+1}) = &\frac{1}{4\tau}\big(\|v_{h}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|v_{h}^{n}\|_{L^{2}(\Omega)}^{2}+\|2v_{h}^{n+1}-v_{h}^{n}\|_{L^{2}(\Omega)}^{2}-\|2v_{h}^{n}-v_{h}^{n-1}\|_{L^{2}(\Omega)}^{2} \\ &+\|v_{h}^{n+1}-2v_{h}^{n}+v_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}\big), \end{align} (4.24)

    we have

    \begin{align} &\frac{\mu}{4\tau}\big(\|\varepsilon(\mathbf{u}_{h}^{n+1})\|_{L^{2}(\Omega)}^{2}-\|\varepsilon(\mathbf{u}_{h}^{n})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(\mathbf{u}_{h}^{n+1})-\varepsilon(\mathbf{u}_{h}^{n})\|_{L^{2}(\Omega)}^{2}-\|2\varepsilon(\mathbf{u}_{h}^{n})-\varepsilon(\mathbf{u}_{h}^{n-1})\|_{L^{2}(\Omega)}^{2} \\ &+\|\varepsilon(\mathbf{u}_{h}^{n+1})-2\varepsilon(\mathbf{u}_{h}^{n})+\varepsilon(\mathbf{u}_{h}^{n-1})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4\tau}\big(\|\xi_{h}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|\xi_{h}^{n}\|_{L^{2}(\Omega)}^{2}+\|2\xi_{h}^{n+1}-\xi_{h}^{n}\|_{L^{2}(\Omega)}^{2}\\ &-\|2\xi_{h}^{n}-\xi_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}+\|\xi_{h}^{n+1}-2\xi_{h}^{n}+\xi_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{2}}{4\tau}\big(\|\eta_{h}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|\eta_{h}^{n}\|_{L^{2}(\Omega)}^{2}+\|2\eta_{h}^{n+1}-\eta_{h}^{n}\|_{L^{2}(\Omega)}^{2}\\ &-\|2\eta_{h}^{n}-\eta_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}+\|\eta_{h}^{n+1}-2\eta_{h}^{n}+\eta_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}\big)+\frac{1}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}\xi_{h}^{n+1}+\kappa_{2}\eta_{h}^{n+1})\|_{L^{2}(\Omega)}^{2}\\ & = (\mathbf{f}, D_{t}\mathbf{u}_{h}^{n+1})+\langle \mathbf{f}_1,D_{t}\mathbf{u}_{h}^{n+1}\rangle+(\phi, p_{h}^{n+1})+\langle \phi_1,p_{h}^{n+1}\rangle. \end{align} (4.25)

    Applying the summation operator \tau\sum_{n = 1}^{l} to both sides of (4.25), there holds

    \begin{align} &\frac{\mu}{4}\big(\|\varepsilon(\mathbf{u}_{h}^{l+1})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(\mathbf{u}_{h}^{l+1})-\varepsilon(\mathbf{u}_{h}^{l})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4}\big(\|\xi_{h}^{l+1}\|_{L^{2}(\Omega)}^{2}+\|2\xi_{h}^{l+1}-\xi_{h}^{l}\|_{L^{2}(\Omega)}^{2}\big) \\ &+\frac{\kappa_{2}}{4}\big(\|\eta_{h}^{l+1}\|_{L^{2}(\Omega)}^{2}+\|2\eta_{h}^{l+1}-\eta_{h}^{l}\|_{L^{2}(\Omega)}^{2}\big)+\sum\limits_{n = 1}^{l}\big(\frac{\mu}{2}\|\varepsilon(\mathbf{u}_{h}^{n+1})-2\varepsilon(\mathbf{u}_{h}^{n})+\varepsilon(\mathbf{u}_{h}^{n-1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\\ &\cdot\|\xi_{h}^{n+1}-2\xi_{h}^{n}+\xi_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|\eta_{h}^{n+1}-2\eta_{h}^{n}+\eta_{h}^{n-1}\|_{L^{2}(\Omega)}^{2}+\frac{\tau}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}\xi_{h}^{n+1}+\kappa_{2}\eta_{h}^{n+1})\|_{L^{2}(\Omega)}^{2}\big)\\ & = \big(\mathbf{f}, \frac{3\mathbf{u}_{h}^{l+1}-\mathbf{u}_{h}^{l}-3\mathbf{u}_{h}^{1}+\mathbf{u}_{h}^{0}}{2}\big)+\langle \mathbf{f}_1,\frac{3\mathbf{u}_{h}^{l+1}-\mathbf{u}_{h}^{l}-3\mathbf{u}_{h}^{1}+\mathbf{u}_{h}^{0}}{2}\rangle+\tau\sum\limits_{n = 1}^{l}\left[(\phi, p_{h}^{n+1})+\langle \phi_1,p_{h}^{n+1}\rangle\right]\\ &+\frac{\mu}{4}\big(\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(\mathbf{u}_{h}^{1})-\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4}\big(\|\xi_{h}^{1}\|_{L^{2}(\Omega)}^{2}+\|2\xi_{h}^{1}-\xi_{h}^{0}\|_{L^{2}(\Omega)}^{2}\big)\\ &+\frac{\kappa_{2}}{4}\big(\|\eta_{h}^{1}\|_{L^{2}(\Omega)}^{2}+\|2\eta_{h}^{1}-\eta_{h}^{0}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.26)

    Similar to (4.13), we can bound the first four terms of the right side of (4.26) as follows:

    \begin{align} \big(\mathbf{f}, \frac{3\mathbf{u}_{h}^{l+1}-\mathbf{u}_{h}^{l}-3\mathbf{u}_{h}^{1}+\mathbf{u}_{h}^{0}}{2}\big)\leq&\frac{4c_{k}^{2}}{\mu}\|\mathbf{f}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{16}\|2\varepsilon(\mathbf{u}_{h}^{l+1})-\varepsilon(\mathbf{u}_{h}^{l})\|_{L^{2}(\Omega)}^{2} \end{align} (4.27)
    \begin{align} &+\frac{\mu}{16}\|\varepsilon(\mathbf{u}_{h}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{3\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2},\\ \langle \mathbf{f}_1,\frac{3\mathbf{u}_{h}^{l+1}-\mathbf{u}_{h}^{l}-3\mathbf{u}_{h}^{1}+\mathbf{u}_{h}^{0}}{2}\rangle\leq&\frac{4c_{a}^{2}c_{k}^{2}}{\mu}\|\mathbf{f}_{1}\|_{L^{2}(\partial\Omega)}^{2}+\frac{\mu}{16}\|2\varepsilon(\mathbf{u}_{h}^{l+1})-\varepsilon(\mathbf{u}_{h}^{l})\|_{L^{2}(\Omega)}^{2} \end{align} (4.28)
    \begin{align} &+\frac{\mu}{16}\|\varepsilon(\mathbf{u}_{h}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{3\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon(\mathbf{u}_{h}^{0})\|_{L^{2}(\Omega)}^{2},\\ \tau\sum\limits_{n = 1}^{l}(\phi, p_{h}^{n+1})\leq&\tau\sum\limits_{n = 1}^{l}\big( c_{p}^{2}\mu_{f}\|K^{-\frac{1}{2}}\phi\|_{L^{2}(\Omega)}^{2}+ \frac{1}{4\mu_{f}}\|K^{\frac{1}{2}}\nabla p_{h}^{n+1}\|_{L^{2}(\Omega)}^{2}\big), \end{align} (4.29)
    \begin{align} \tau\sum\limits_{n = 1}^{l}\langle \phi_1,p_{h}^{\frac{1}{2}}\rangle\leq& \tau\sum\limits_{n = 1}^{l}\big(c_{a}^{2}c_{p}^{2}\mu_{f}\|K^{-\frac{1}{2}}\phi_{1}\|_{L^{2}(\partial\Omega)}^{2}+\frac{1}{4\mu_{f}} \|K^{\frac{1}{2}}\nabla p_{h}^{\frac{1}{2}}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.30)

    Substituting (4.27)–(4.30) into (4.26) and using (4.20) for the resulting inequality, we can deduce that (4.9) holds. The proof is complete.

    In this subsection, we consider the optimal convergence analysis of the BDF2 algorithm. Before that, we have the following estimate for the BDF2 algorithm.

    Theorem 6. Let \{(\mathbf{u}_h^{n}, \xi_h^{n}, \eta_h^{n})\}_{n\geq 0} be defined by the BDF2 algorithm. Then there exists a positive constant C_{4} independent of mesh size h and time step \tau such that

    \begin{align} &\frac{\mu}{8}\|\varepsilon(Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2} \\ &+\frac{\tau}{2\mu_{f}}\sum\limits_{n = 1}^{l}\|K^{\frac{1}{2}} (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1})\|_{L^{2}(\Omega)}^{2}\leq C_{4}(h^{4}+\tau^{4}). \end{align} (4.31)

    Proof. First, it is easy to check that there exists a positive constant c_{b} such that

    \begin{align} \left\|\mathrm{div}\mathbf{w} \right\|_{L^{2}(\Omega)}\leq c_{b}\left\|\varepsilon (\mathbf{w}) \right\|_{L^{2}(\Omega)},\; \forall\; \mathbf{w}\in \mathbf{H}^{1}(\Omega). \end{align} (4.32)

    For n\geq 1 , at t = t_{n+1} , we have for (2.7)–(2.10) that

    \begin{align} {2} &\mu\bigl(\varepsilon(\mathbf{u}(t_{n+1})), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( \xi(t_{n+1}), {\rm div} \mathbf{v}_h \bigr) = (\mathbf{f}, \mathbf{v}_h)+\langle \mathbf{f}_1,\mathbf{v}_h\rangle, &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.33)
    \begin{align} &\kappa_3\bigl(\xi(t_{n+1}), \varphi_h \bigr) +\bigl({\rm div}\mathbf{u}(t_{n+1}), \varphi_h \bigr) = \kappa_1\bigl( \eta(t_{n+1}), \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.34)
    \begin{align} &\bigl(\eta_{t}(t_{n+1}), \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}\xi(t_{n+1})+\kappa_{2}\eta(t_{n+1})),\nabla\psi_h \bigr) = (\phi, \psi_h)+\langle \phi_1,\psi_h\rangle, &&\quad \forall \psi_h\in W_h. \end{align} (4.35)

    Moreover, we have from (2.7)–(2.10) at t = t_{0} and t = t_{1} that

    \begin{align} {2} &\mu\bigl(\varepsilon(\mathbf{u}(t_{\frac{1}{2}})), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( \xi(t_{\frac{1}{2}}), {\rm div} \mathbf{v}_h \bigr) = (\mathbf{f}, \mathbf{v}_h)+\langle \mathbf{f}_1,\mathbf{v}_h\rangle, &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.36)
    \begin{align} &\kappa_3\bigl(\xi(t_{\frac{1}{2}}), \varphi_h \bigr) +\bigl({\rm div}\mathbf{u}(t_{\frac{1}{2}}), \varphi_h \bigr) = \kappa_1\bigl( \eta(t_{\frac{1}{2}}), \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.37)
    \begin{align} &\bigl(\eta_{t}(t_{\frac{1}{2}}), \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}\xi(t_{\frac{1}{2}})+\kappa_{2}\eta(t_{\frac{1}{2}})),\nabla\psi_h \bigr) = (\phi, \psi_h)+\langle \phi_1,\psi_h\rangle, &&\quad \forall \psi_h\in W_h, \end{align} (4.38)

    where v(t_{\frac{1}{2}}) = (v(t_{0})+v(t_{1}))/2 .

    Subtracting (4.2)–(4.4) from (4.36)–(4.38), there holds

    \begin{align} {2} &\mu\bigl(\varepsilon(E_{\mathbf{u}}^{\frac{1}{2}}), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( E_{\xi}^{\frac{1}{2}}, {\rm div} \mathbf{v}_h \bigr) = 0, &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.39)
    \begin{align} &\kappa_3\bigl(E_{\xi}^{\frac{1}{2}}, \varphi_h \bigr) +\bigl({\rm div}E_{\mathbf{u}}^{\frac{1}{2}}, \varphi_h \bigr) = \kappa_1\bigl( E_{\eta}^{\frac{1}{2}}, \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.40)
    \begin{align} &\bigl(d_{t}E_{\eta}^{1}, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}E_{\xi}^{\frac{1}{2}}+\kappa_{2}E_{\eta}^{\frac{1}{2}}),\nabla\psi_h \bigr) = (R_{1,\eta},\psi_{h})+(R_{2,\eta},\psi_{h}), &&\quad \forall \psi_h\in W_h, \end{align} (4.41)

    where

    \begin{align*} &R_{1,\eta}: = \frac{1}{2\tau}\int_{t_{0}}^{t_{1}}\Gamma^{2}(t)\eta_{ttt}dt,\qquad R_{2,\eta}: = \frac{1}{2}\int_{t_{0}}^{t_{1}}\Gamma(t)\eta_{ttt}dt,\\ &\Gamma(t): = \left\lbrace \begin{aligned} &t_{1}-t &&{\rm when}\; t\geq \frac{t_{0}+t_{1}}{2},\\ &t-t_{0}&&{\rm when}\; t < \frac{t_{0}+t_{1}}{2}. \end{aligned}\right. \end{align*}

    Using the definitions of the projection operators \mathcal{Q}_{h}, \mathcal{S}_{h}, \mathcal{R}_{h} for (4.39)–(4.41), we have

    \begin{align} {2} &\mu\bigl(\varepsilon(Z_{\mathbf{u}}^{\frac{1}{2}}), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( G_{\xi}^{\frac{1}{2}}, {\rm div} \mathbf{v}_h \bigr) = \bigl( F_{\xi}^{\frac{1}{2}}, {\rm div} \mathbf{v}_h \bigr), &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.42)
    \begin{align} &\kappa_3\bigl(G_{\xi}^{\frac{1}{2}}, \varphi_h \bigr) +\bigl({\rm div}Z_{\mathbf{u}}^{\frac{1}{2}}, \varphi_h \bigr) = \kappa_1\bigl( G_{\eta}^{\frac{1}{2}}, \varphi_h \bigr)-\bigl({\rm div}Y_{\mathbf{u}}^{\frac{1}{2}}, \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.43)
    \begin{align} &\bigl(d_{t}G_{\eta}^{1}, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}Z_{\xi}^{\frac{1}{2}}+\kappa_{2}Z_{\eta}^{\frac{1}{2}}),\nabla\psi_h \bigr) = (R_{1,\eta},\psi_{h})+(R_{2,\eta},\psi_{h}), &&\quad \forall \psi_h\in W_h. \end{align} (4.44)

    Setting \mathbf{v}_{h} = d_{t}Z_{\mathbf{u}}^{1} in (4.42), \varphi_{h} = d_{t}G_{\xi}^{1} in (4.43) and \psi_{h} = Z_{p}^{\frac{1}{2}} in (4.44), we get

    \begin{align} {2} &\frac{\mu}{2\tau}\big(\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}-\|\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}\big)-\bigl( G_{\xi}^{\frac{1}{2}}, {\rm div} d_{t}Z_{\mathbf{u}}^{1} \bigr) = \bigl( F_{\xi}^{\frac{1}{2}}, {\rm div} d_{t}Z_{\mathbf{u}}^{1} \bigr), && \end{align} (4.45)
    \begin{align} &\frac{\kappa_{3}}{2\tau}\big(\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}-\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}\big) +\bigl({\rm div}Z_{\mathbf{u}}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr) = \kappa_1\bigl( G_{\eta}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr)-\bigl({\rm div}Y_{\mathbf{u}}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr), && \end{align} (4.46)
    \begin{align} &\frac{\kappa_{2}}{2\tau}\big(\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}-\|G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}\big)+\kappa_{1}\bigl(d_{t}G_{\eta}^{\frac{1}{2}}, G_{\xi}^{\frac{1}{2}} \bigr) +\frac{1}{\mu_f} \|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{\frac{1}{2}}+\kappa_{2}Z_{\eta}^{\frac{1}{2}})\|_{L^{2}(\Omega)}^{2}\\ = &(R_{1,\eta},Z_{p}^{\frac{1}{2}})+(R_{2,\eta},Z_{p}^{\frac{1}{2}})+(d_{t}G_{\eta}^{1},Y_{p}^{\frac{1}{2}}-F_{p}^{\frac{1}{2}}). \end{align} (4.47)

    Adding (4.45)–(4.47), there holds

    \begin{align} &\frac{\mu}{2\tau}\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{2\tau}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{2\tau}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{1}{\mu_f} \|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{\frac{1}{2}}+\kappa_{2}Z_{\eta}^{\frac{1}{2}})\|_{L^{2}(\Omega)}^{2} \\ & = \frac{\mu}{2\tau}\|\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{2\tau}\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{2\tau}\|G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}+ \bigl( F_{\xi}^{\frac{1}{2}}, {\rm div} d_{t}Z_{\mathbf{u}}^{1} \bigr)-\bigl({\rm div}Y_{\mathbf{u}}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr)\\ &\quad+(R_{1,\eta},Z_{p}^{\frac{1}{2}})+(R_{2,\eta},Z_{p}^{\frac{1}{2}})+(d_{t}G_{\eta}^{1},Y_{p}^{\frac{1}{2}}-F_{p}^{\frac{1}{2}})+\frac{1}{\tau}\Big((G_{\xi}^{0},{\rm div}Z_{\mathbf{u}}^{1})-(G_{\xi}^{1},{\rm div}Z_{\mathbf{u}}^{0})\Big)\\ &\quad+\frac{\kappa_{1}}{\tau}\Big((G_{\eta}^{0},G_{\xi}^{1})-(G_{\eta}^{1},G_{\xi}^{0})\Big). \end{align} (4.48)

    Then, we bound the fourth to tenth terms on the right side of (4.48). Using the Cauchy-Schwarz inequality, triangle equality, and (4.32), Poincaré and Young inequalities, we have

    \begin{align} \bigl( F_{\xi}^{\frac{1}{2}}, {\rm div} d_{t}Z_{\mathbf{u}}^{1} \bigr) &\leq\frac{1}{2\tau}\big(\frac{2c_{b}^{2}}{\mu}\|F_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{2c_{b}^{2}}{\mu}\|F_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{2}\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2} \end{align} (4.49)
    \begin{align} &\quad+\frac{\mu}{2}\|\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}\big),\\ -\bigl({\rm div}Y_{\mathbf{u}}^{\frac{1}{2}}, d_{t}G_{\xi}^{1} \bigr) &\leq\frac{1}{2\tau}\big(\frac{4c_{b}^{2}}{\kappa_{3}}\|\varepsilon(Y_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{4c_{b}^{2}}{\kappa_{3}}\|\varepsilon(Y_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2} \end{align} (4.50)
    \begin{align} &\quad+\frac{\kappa_{3}}{4}\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}\big),\\ (R_{1,\eta},Z_{p}^{\frac{1}{2}}) &\leq \frac{c_{p}^{2}\mu_{f}\tau^{3}}{640K_{1}}\|\eta_{ttt}\|_{L^{2}(t_{0},t_{1};L^{2}(\Omega))}^{2}+\frac{1}{4\mu_{f}}\|K^{\frac{1}{2}}\nabla Z_{p}^{\frac{1}{2}}\|_{L^{2}(\Omega)}^{2}, \end{align} (4.51)
    \begin{align} (R_{2,\eta},Z_{p}^{\frac{1}{2}}) &\leq \frac{c_{p}^{2}\mu_{f}\tau^{3}}{96K_1}\|\eta_{ttt}\|_{L^{2}(t_{0},t_{1};L^{2}(\Omega))}^{2}+\frac{1}{4\mu_{f}}\|K^{\frac{1}{2}}\nabla Z_{p}^{\frac{1}{2}}\|_{L^{2}(\Omega)}^{2}, \end{align} (4.52)
    \begin{align} (d_{t}G_{\eta}^{1},Y_{p}^{\frac{1}{2}}-F_{p}^{\frac{1}{2}}) &\leq\frac{1}{2\tau}\big(\frac{\kappa_{2}}{2}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{2}\|G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{4}{\kappa_{2}}\|Y_{p}^{1}\|_{L^{2}(\Omega)}^{2} \end{align} (4.53)
    \begin{align} &\quad+\frac{4}{\kappa_{2}}\|Y_{p}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{4}{\kappa_{2}}\|F_{p}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{4}{\kappa_{2}}\|F_{p}^{0}\|_{L^{2}(\Omega)}^{2}\big),\\ (G_{\xi}^{0},{\rm div}Z_{\mathbf{u}}^{1})-(G_{\xi}^{1},{\rm div}Z_{\mathbf{u}}^{0})&\leq\frac{2c_{b}^{2}}{\mu}\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{8}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2} \end{align} (4.54)
    \begin{align} &\quad+\frac{2c_{b}^{2}}{\kappa_{3}}\|\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2},\\ \kappa_{1}(G_{\eta}^{0},G_{\xi}^{1})-\kappa_{1}(G_{\eta}^{1},G_{\xi}^{0})&\leq\frac{2\kappa_{1}^{2}}{\kappa_{3}}\|G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{8}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{1}^{2}}{\kappa_{2}}\|G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}. \end{align} (4.55)

    Substituting (4.49)–(4.55) into (4.48) and using (3.11)–(3.13), the facts of Z_{\mathbf{u}}^{0} = \mathbf{0}, \; G_{\xi}^{0} = 0 , and G_{\eta}^{0} = 0 , we find there exists a positive constant \widetilde{C} = \widetilde{C}\big(\|\mathbf{u}\|_{L^{\infty}(t_{0}, t_{1};H^{2}(\Omega))}, \|\xi\|_{L^{\infty}(t_{0}, t_{1};H^{2}(\Omega))}, \|\eta\|_{L^{\infty}(t_{0}, t_{1};H^{2}(\Omega))}, \\ \|\eta_{ttt}\|_{L^{2}(t_{0}, t_{1};L^{2}(\Omega))}, c_{b}, \mu, \kappa_{3}, c_{p}, \mu_{f}, K_{1}, \kappa_{2}, \kappa_{1}\big) such that

    \begin{align} &\mu\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\kappa_{3}\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\kappa_{2}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\frac{\tau}{\mu_f} \|K\nabla (\kappa_{1}Z_{\xi}^{\frac{1}{2}}+\kappa_{2}Z_{\eta}^{\frac{1}{2}})\|_{L^{2}(\Omega)}^{2}\leq \widetilde{C}(h^{4}+\tau^{4}). \end{align} (4.56)

    Next, we show the error estimate at t = t_{n+1} . Subtracting (4.5)–(4.7) from (4.33)–(4.35), it holds

    \begin{align} {2} &\mu\bigl(\varepsilon(Z_{\mathbf{u}}^{n+1}), \varepsilon(\mathbf{v}_h) \bigr)-\bigl( G_{\xi}^{n+1}, {\rm div} \mathbf{v}_h \bigr) = \bigl( F_{\xi}^{n+1}, {\rm div} \mathbf{v}_h \bigr), &&\quad \forall \mathbf{v}_h\in \mathbf{V}_h, \end{align} (4.57)
    \begin{align} &\kappa_3\bigl(G_{\xi}^{n+1}, \varphi_h \bigr) +\bigl({\rm div}Z_{\mathbf{u}}^{n+1}, \varphi_h \bigr) = \kappa_1\bigl( G_{\eta}^{n+1}, \varphi_h \bigr)-\bigl({\rm div}Y_{\mathbf{u}}^{n+1}, \varphi_h \bigr), &&\quad \forall \varphi_h \in M_h, \end{align} (4.58)
    \begin{align} &\bigl(D_{t}G_{\eta}^{n+1}, \psi_h \bigr) +\frac{1}{\mu_f} \bigl(K\nabla (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1}),\nabla\psi_h \bigr) = \bigl(R_{3,\eta}, \psi_h \bigr), &&\quad \forall \psi_h\in W_h, \end{align} (4.59)

    where

    \begin{align*} R_{3,\eta} = \frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}(t-t_{n})^{2}\eta_{ttt}dt-\frac{1}{4\tau}\int_{t_{n-1}}^{t_{n+1}}(t-t_{n+1})^{2}\eta_{ttt}dt. \end{align*}

    Setting \mathbf{v}_{h} = D_{t}Z_{\mathbf{u}}^{n+1} in (4.57), \varphi_{h} = G_{\xi}^{n+1} in (4.58) after using the operator D_{t} to both sides of (4.58), \psi_{h} = Z_{p}^{n+1} in (4.59), using (4.24) and adding the result equation, we get

    \begin{align} &\frac{\mu}{4\tau}\big(\|\varepsilon(Z_{\mathbf{u}}^{n+1})\|_{L^{2}(\Omega)}^{2}-\|\varepsilon(Z_{\mathbf{u}}^{n})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(Z_{\mathbf{u}}^{n+1})-\varepsilon(Z_{\mathbf{u}}^{n})\|_{L^{2}(\Omega)}^{2}-\|2\varepsilon(Z_{\mathbf{u}}^{n})-\varepsilon(Z_{\mathbf{u}}^{n-1})\|_{L^{2}(\Omega)}^{2} \\ &+\|\varepsilon(Z_{\mathbf{u}}^{n+1})-2\varepsilon(Z_{\mathbf{u}}^{n})+\varepsilon(Z_{\mathbf{u}}^{n-1})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4\tau}\big(\|G_{\xi}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|G_{\xi}^{n}\|_{L^{2}(\Omega)}^{2}+\|2G_{\xi}^{n+1}-G_{\xi}^{n}\|_{L^{2}(\Omega)}^{2}\\ &-\|2G_{\xi}^{n}-G_{\xi}^{n-1}\|_{L^{2}(\Omega)}^{2}+\|G_{\xi}^{n+1}-2G_{\xi}^{n}+G_{\xi}^{n-1}\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{2}}{4\tau}\big(\|G_{\eta}^{n+1}\|_{L^{2}(\Omega)}^{2}-\|G_{\eta}^{n}\|_{L^{2}(\Omega)}^{2}\\ &+\|2G_{\eta}^{n+1}-G_{\eta}^{n}\|_{L^{2}(\Omega)}^{2}-\|2G_{\eta}^{n}-G_{\eta}^{n-1}\|_{L^{2}(\Omega)}^{2}+\|G_{\eta}^{n+1}-2G_{\eta}^{n}+G_{\eta}^{n-1}\|_{L^{2}(\Omega)}^{2}\big)\\ &+\frac{1}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1})\|_{L^{2}(\Omega)}^{2} = \bigl( F_{\xi}^{n+1}, {\rm div}D_{t}Z_{\mathbf{u}}^{n+1} \bigr)-\bigl({\rm div}D_{t}Y_{\mathbf{u}}^{n+1}, G_{\xi}^{n+1}\bigr)+\bigl(R_{3,\eta}, Z_{p}^{n+1} \bigr)\\ &+\big(D_{t}G_{\eta}^{n+1},Y_{p}^{n+1}-F_{p}^{n+1}\big). \end{align} (4.60)

    Applying the summation operator \tau\sum_{n = 1}^{l} to both sides of (4.60), there holds

    \begin{align} &\frac{\mu}{4}\big(\|\varepsilon(Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(Z_{\mathbf{u}}^{l+1})-\varepsilon(Z_{\mathbf{u}}^{l})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4}\big(\|G_{\xi}^{l+1}\|_{L^{2}(\Omega)}^{2}+\|2G_{\xi}^{l+1}-G_{\xi}^{l}\|_{L^{2}(\Omega)}^{2}\big) \\ &+\frac{\kappa_{2}}{4}\big(\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2}+\|2G_{\eta}^{l+1}-G_{\eta}^{l}\|_{L^{2}(\Omega)}^{2}\big)+\sum\limits_{n = 1}^{l}\big(\frac{\mu}{2}\|\varepsilon(Z_{\mathbf{u}}^{n+1})-2\varepsilon(Z_{\mathbf{u}}^{n})+\varepsilon(Z_{\mathbf{u}}^{n-1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\\ &\cdot\|G_{\xi}^{n+1}-2G_{\xi}^{n}+G_{\xi}^{n-1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{4}\|G_{\eta}^{n+1}-2G_{\eta}^{n}+G_{\eta}^{n-1}\|_{L^{2}(\Omega)}^{2}+\frac{\tau}{\mu_{f}}\|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1})\|_{L^{2}(\Omega)}^{2}\big)\\ & = \tau\sum\limits_{n = 1}^{l}\left[\bigl( F_{\xi}^{n+1}, {\rm div}D_{t}Z_{\mathbf{u}}^{n+1} \bigr)-\bigl({\rm div}D_{t}Y_{\mathbf{u}}^{n+1}, G_{\xi}^{n+1}\bigr)+\bigl(R_{3,\eta}, Z_{p}^{n+1} \bigr)+\big(D_{t}G_{\eta}^{n+1},Y_{p}^{n+1}-F_{p}^{n+1}\big)\right]\\ &+\frac{\mu}{4}\big(\|\varepsilon(Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}+\|2\varepsilon(Z_{\mathbf{u}}^{1})-\varepsilon(Z_{\mathbf{u}}^{0})\|_{L^{2}(\Omega)}^{2}\big)+\frac{\kappa_{3}}{4}\big(\|G_{\xi}^{1}\|_{L^{2}(\Omega)}^{2}+\|2G_{\xi}^{1}-G_{\xi}^{0}\|_{L^{2}(\Omega)}^{2}\big)\\ &+\frac{\kappa_{2}}{4}\big(\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}+\|2G_{\eta}^{1}-G_{\eta}^{0}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.61)

    Now, we begin to estimate the first four terms on the right side of (4.61). For the first term, it is easy to check that

    \begin{align} \tau\sum\limits_{n = 1}^{l}\bigl( F_{\xi}^{n+1}, {\rm div}D_{t}Z_{\mathbf{u}}^{n+1} \bigr) = &\frac{1}{2}\Big(\big(F_{\xi}^{l+1},{\rm div}(3Z_{\mathbf{u}}^{l+1}-Z_{\mathbf{u}}^{l})\big)-\big(F_{\xi}^{l+1}-F_{\xi}^{l},{\rm div}Z_{\mathbf{u}}^{l}\big)-\big(3F_{\xi}^{2},{\rm div}Z_{\mathbf{u}}^{1}\big) \\ &+\big(F_{\xi}^{3},{\rm div}Z_{\mathbf{u}}^{1}\big)\Big)+\tau\sum\limits_{n = 3}^{l}\big(\frac{F_{\xi}^{n+1}-4F_{\xi}^{n}+3F_{\xi}^{n-1}}{2\tau},{\rm div}Z_{\mathbf{u}}^{n-1}\big). \end{align} (4.62)

    Applying the Cauchy-Schwarz inequality, (4.32), triangle and Young inequalities to (4.62), we have

    \begin{align} \tau\sum\limits_{n = 1}^{l}\bigl( F_{\xi}^{n+1}, {\rm div}D_{t}Z_{\mathbf{u}}^{n+1} \bigr)\leq&\frac{1}{2}\Big(\frac{4c_{b}^{2}}{\mu}\|F_{\xi}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|2\varepsilon (Z_{\mathbf{u}}^{l+1})-\varepsilon (Z_{\mathbf{u}}^{l})\|_{L^{2}(\Omega)}^{2} \\ &+\frac{\mu}{8}\|\varepsilon (Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{10c_{b}^{2}}{\mu}\|F_{\xi}^{l+1}-F_{\xi}^{l}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{8}\|\varepsilon (Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}\\ &+\frac{\mu}{8}\|2\varepsilon (Z_{\mathbf{u}}^{l+1})-\varepsilon (Z_{\mathbf{u}}^{l})\|_{L^{2}(\Omega)}^{2}+\frac{c_{b}^{2}}{\mu}\|F_{\xi}^{3}-3F_{\xi}^{2}\|_{L^{2}(\Omega)}^{2}+\frac{\mu}{4}\|\varepsilon (Z_{\mathbf{u}}^{1})\|_{L^{2}(\Omega)}^{2}\Big)\\ &+\tau\sum\limits_{n = 3}^{l}\Big(\frac{\mu}{8}\|\varepsilon (Z_{\mathbf{u}}^{n-1})\|_{L^{2}(\Omega)}^{2}+\frac{2c_{b}^{2}}{\mu}\|\frac{F_{\xi}^{n+1}-4F_{\xi}^{n}+3F_{\xi}^{n-1}}{2\tau}\|_{L^{2}(\Omega)}^{2}\Big). \end{align} (4.63)

    For the second term, using the Canchy-Schwarz and Young inequalities, we obtain

    \begin{align} -\tau\sum\limits_{n = 1}^{l}\bigl({\rm div}D_{t}Y_{\mathbf{u}}^{n+1}, G_{\xi}^{n+1}\bigr)\leq\tau\sum\limits_{n = 1}^{l}\big(\frac{1}{\kappa_{3}}\|{\rm div}D_{t}Y_{\mathbf{u}}^{n+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{n+1}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.64)

    Applying the Canchy-Schwarz inequality, Poincaré inequality, Young inequality, and the fact of

    \begin{align*} \|R_{3,\eta}\|_{L^{2}(\Omega)}^{2}\leq\frac{\tau^{3}}{5}\|\eta_{ttt}\|_{L^{2}(t_{n},t_{n+1};L^{2}(\Omega))}^{2}+\frac{2\tau^{3}}{5}\|\eta_{ttt}\|_{L^{2}(t_{n-1},t_{n+1};L^{2}(\Omega))}^{2} \end{align*}

    to the third term, there holds

    \begin{align} \tau\sum\limits_{n = 1}^{l}\bigl(R_{3,\eta}, Z_{p}^{n+1} \bigr)&\leq\tau\sum\limits_{n = 1}^{l} \big(\frac{c_{p}^{2}\mu_{f}}{2}\|R_{3,\eta}\|_{L^{2}(\Omega)}^{2}+\frac{1}{2\mu_{f}}\|\nabla Z_{p}^{n+1}\|_{L^{2}(\Omega)}^{2}\big)\\ &\leq\tau\sum\limits_{n = 1}^{l} \big(\frac{3\tau^{3}c_{p}^{2}\mu_{f}}{10K_{1}}\|\eta_{ttt}\|_{L^{2}(t_{n-1},t_{n+1};L^{2}(\Omega))}^{2}+\frac{1}{2\mu_{f}}\|K^{\frac{1}{2}}\nabla Z_{p}^{n+1}\|_{L^{2}(\Omega)}^{2}\big). \end{align} (4.65)

    Similar to the first term, we can bound the fourth term as follows:

    \begin{align} \tau\sum\limits_{n = 1}^{l}\big(D_{t}G_{\eta}^{n+1},Y_{p}^{n+1}-F_{p}^{n+1}\big)\leq&\frac{1}{2}\Big(\frac{4}{\kappa_{2}}\|Y_{p}^{l+1}-F_{p}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|2G_{\eta}^{l+1}-G_{\eta}^{l}\|_{L^{2}(\Omega)}^{2} \\ &+\frac{\kappa_{2}}{8}\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{10}{\kappa_{2}}\|Y_{p}^{l+1}-Y_{p}^{l}-F_{p}^{l+1}+F_{p}^{l}\|_{L^{2}(\Omega)}^{2}\\ &+\frac{\kappa_{2}}{8}\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|2G_{\eta}^{l+1}-G_{\eta}^{l}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|G_{\eta}^{1}\|_{L^{2}(\Omega)}^{2}\\ &+\frac{2}{\kappa_{2}}\|Y_{p}^{3}-3Y_{p}^{2}-F_{p}^{3}+3F_{p}^{2}\|_{L^{2}(\Omega)}^{2}\Big)+\tau\sum\limits_{n = 3}^{l}\Big(\frac{\kappa_{2}}{8}\|G_{\eta}^{n-1}\|_{L^{2}(\Omega)}^{2}\\ &+\frac{2}{\kappa_{2}}\|\frac{Y_{p}^{n+1}-4Y_{p}^{n}+3Y_{p}^{n-1}}{2\tau}-\frac{F_{p}^{n+1}-4F_{p}^{n}+3F_{p}^{n-1}}{2\tau}\|_{L^{2}(\Omega)}^{2}\Big). \end{align} (4.66)

    Substituting (4.63)–(4.66) into (4.61), using the facts of Z_{\mathbf{u}}^{0} = \mathbf{0}, \; G_{\xi}^{0} = 0, \; G_{\eta}^{0} = 0 , (4.56), (3.11)–(3.13), there exists a positive constant \widehat{C} = \widehat{C}(c_{b}, \mu, \kappa_{3}, \mu_{f}, K_{1}, c_{p}, \kappa_{2}) such that

    \begin{align} &\frac{\mu}{8}\|\varepsilon(Z_{\mathbf{u}}^{l+1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{l+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|G_{\eta}^{l+1}\|_{L^{2}(\Omega)}^{2} \\ &\quad+\frac{\tau}{2\mu_{f}}\sum\limits_{n = 1}^{l}\|K^{\frac{1}{2}}\nabla (\kappa_{1}Z_{\xi}^{n+1}+\kappa_{2}Z_{\eta}^{n+1})\|_{L^{2}(\Omega)}^{2}\\ &\leq \tau\sum\limits_{n = 1}^{l}\big(\frac{\mu}{8}\|\varepsilon (Z_{\mathbf{u}}^{n+1})\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{3}}{4}\|G_{\xi}^{n+1}\|_{L^{2}(\Omega)}^{2}+\frac{\kappa_{2}}{8}\|G_{\eta}^{n+1}\|_{L^{2}(\Omega)}^{2}\big)+ \widehat{C}\big(h^{4}\|\xi\|_{L^{\infty}(0,T;H^{2}(\Omega))}\\ &\quad+h^{4}\|\xi_{t}\|_{L^{2}(0,T;H^{2}(\Omega))}+h^{4}\|\nabla\mathbf{u}_{t}\|_{L^{2}(0,T;H^{2}(\Omega))}+\tau^{4}\|\eta_{ttt}\|_{L^{2}(0,T;L^{2}(\Omega))}+h^{4}\|\eta\|_{L^{\infty}(0,T;H^{2}(\Omega))}\\ &\quad+h^{4}\|\eta_{t}\|_{L^{2}(0,T;H^{2}(\Omega))}+\widetilde{C}(h^{4}+\tau^{4})\big). \end{align} (4.67)

    Applying the Gronwall inequality to (4.67), one has that (4.31) holds. The proof is complete.

    Furthermore, we state the main result of this part by a theorem.

    Theorem 7. There exists a positive constant C_{5} independent of h and \tau such that the solution of the BDF2 algorithm satisfies the following error estimates:

    \begin{align} &\max\limits_{0\leq n\leq N}\left[\|\nabla(\mathbf{u}(t_{n})-\mathbf{u}_{h}^{n})\|_{L^{2}(\Omega)}+\sqrt\kappa_{2}\|\eta(t_{n})-\eta_{h}^{n}\|_{L^{2}(\Omega )}\right. \end{align} (4.68)
    \begin{align} &\left.+\sqrt{\kappa_{3}}\|\xi(t_{n})-\xi_{h}^{n}\|_{L^{2}(\Omega )} \right]\leq C_{5}(h^{2}+\tau^{2}),\\ &\max\limits_{0\leq n\leq N}\|p(t_{n})-p_{h}^{n}\|_{L^{2}(\Omega )}\leq C_{5}(h^{2}+\tau^{2}). \end{align} (4.69)

    Proof. Using (3.11)–(3.13), the relation p = \kappa_{1}\xi+\kappa_{2}\eta , Theorem 6, (4.56), and the triangle inequality on

    \begin{align*} &\mathbf{u}(t_{n})-\mathbf{u}_{h}^{n} = Y_{\mathbf{u}}^{n}+Z_{\mathbf{u}}^{n},\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \xi(t_{n})-\xi_{h}^{n} = Y_{\xi}^{n}+Z_{\xi}^{n} = F_{\xi}^{n}+G_{\xi}^{n},\\ &\eta(t_{n})-\eta_{h}^{n} = Y_{\eta}^{n}+Z_{\eta}^{n} = F_{\eta}^{n}+G_{\eta}^{n},\; \; \; \; \; \; \; p(t_{n})-p_{h}^{n} = Y_{p}^{n}+Z_{p}^{n} = F_{p}^{n}+G_{p}^{n}, \end{align*}

    one can see that (4.68) and (4.69) hold. The proof is complete.

    Here the \mathbf{P}_{2}-P_{1}-P_{1}-P_{1} element pairs are employed for the coupled algorithm and decoupled algorithms, and \mathbf{P}_{2}-P_{1}-P_{1} element pairs are employed for the BDF2 algorithm. To verify the robustness of our algorithms with respect to parameters, we compute the spatial errors and convergence rates of the variables \mathbf{u} and p for different values of Poisson rates and permeability tensor separately. In the meantime, we show the effectiveness of our algorithms in overcoming the Poisson locking as \nu\rightarrow 0.5 . After that, to illustrate that there is no nonphysical oscillation in the pressure field, we use the decoupled algorithm to test a Barry-Mercer problem. Finally, we use the following examples to compute the convergence rates for the above proposed algorithms in subsections 5.1, 5.2, and 5.4.

    Example 1. Let \Omega = (0, 1)^2 , \Gamma_1 = \{(1, x_2); 0\leq x_2\leq 1\} , \Gamma_2 = \{(x_1, 0); 0\leq x_1\leq 1\} , \Gamma_3 = \{(0, x_2); 0\leq x_2\leq 1\} , \Gamma_4 = \{(x_1, 1); 0\leq x_1\leq 1\} and T = 1 . The boundary and initial conditions are

    \begin{alignat*} {2} p = t\sin(\pi x_{1})\cos(\pi x_{2}) &\qquad\mbox{on }\partial\Omega_T,\\ u_1 = e^{-t}\sin(2\pi x_{2})(-1+\cos(2\pi x_{1}))+\dfrac{1}{\mu+\lambda}\sin(\pi x_{1})\sin(\pi x_{2}) &\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 1,3,\\ u_2 = e^{-t}\sin(2\pi x_{1})(1-\cos(2\pi x_{2}))+\dfrac{1}{\mu+\lambda}\sin(\pi x_{1})\sin(\pi x_{2}) &\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 2,4,\\ \sigma({\mathbf{u}})\mathbf{n}-\alpha p \mathbf{n} = \mathbf{f}_1 &\qquad \mbox{on } \partial\Omega_T,\\ \mathbf{u}(\mathbf{x},0) = \mathbf{0}, \quad p(\mathbf{x},0) = 0.& \end{alignat*}

    The exact solutions of (1.1)–(1.4) and (2.3)–(2.5) are p = t\sin(\pi x_{1})\cos(\pi x_{2}) and

    \begin{alignat*} {2} u_1 & = e^{-t}\sin(2\pi x_{2})(-1+\cos(2\pi x_{1}))+\dfrac{1}{\mu+\lambda}\sin(\pi x_{1})\sin(\pi x_{2}),\\ u_2 & = e^{-t}\sin(2\pi x_{1})(1-\cos(2\pi x_{2}))+\dfrac{1}{\mu+\lambda}\sin(\pi x_{1})\sin(\pi x_{2}). \end{alignat*}

    We first test the robustness of the coupled and decoupled algorithms on the Poisson ratio \nu by fixing \alpha = 1 , \mu_{f} = 1 , E = 2.8\text{e}3 , K = 1 e-7 \mathbf{I} , and c_{0} = 0.2 for \nu = 0.49 and \nu = 0.4999999 , respectively.

    Tables 1 and 2 show the spatial errors and convergence rates of the coupled algorithm for \nu = 0.49 and \nu = 0.4999999 , respectively. Numerical results show that the H^{1} convergence rates of \mathbf{u} and the L^{2} convergence rates of p are all approaching 2 . The L^{2} convergence rates of \mathbf{u} are around 3 . The H^{1} convergence rates of p are around 1 . It can be concluded that the convergence rates of \mathbf{u} and p are optimal, which is consistent with Theorem 4.

    Table 1.  Spatial errors and convergence rates of \mathbf{u} and p of the coupled algorithm \nu = 0.49 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 8.5755e-2 4.2637e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 1.7262e-2 2.31 2.0452e-1 1.06
    h=1/16, \tau=1/256 6.5936e-4 3.10 1.1349e-2 1.96 3.8013e-3 2.18 9.9924e-2 1.03
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 9.8977e-4 1.94 4.9462e-2 1.01

     | Show Table
    DownLoad: CSV
    Table 2.  Spatial errors and convergence rates of \mathbf{u} and p of the coupled algorithm for \nu = 0.4999999 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7651e-2 1.6055e-1 8.5700e-2 4.2574e-1
    h=1/8, \tau=1/64 5.6669e-3 3.07 4.4078e-2 1.86 1.7255e-2 2.31 2.0459e-1 1.06
    h=1/16, \tau=1/256 6.5948e-4 3.10 1.1349e-2 1.96 3.7745e-3 2.19 9.9988e-2 1.03
    h=1/32, \tau=1/1024 7.9662e-5 3.05 2.8614e-3 1.99 8.7440e-4 2.11 4.9501e-2 1.01

     | Show Table
    DownLoad: CSV

    For the decoupled algorithm, we also take \nu = 0.49 and \nu = 0.4999999 separately to observe the convergence rates of \mathbf{u} and p . Tables 3 and 4 show that the convergence rates of all variables are optimal.

    Table 3.  Spatial errors and convergence rates of \mathbf{u} and p of the decoupled algorithm for \nu = 0.49 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 8.5755e-2 4.2637e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 1.7262e-2 2.31 2.0452e-1 1.06
    h=1/16, \tau=1/256 6.5935e-4 3.10 1.1349e-2 1.96 3.8013e-3 2.18 9.9924e-2 1.03
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 9.8977e-4 1.94 4.9462e-2 1.01

     | Show Table
    DownLoad: CSV
    Table 4.  Spatial errors and convergence rates of \mathbf{u} and p of the decoupled algorithm for \nu = 0.4999999 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7651e-2 1.6055e-1 8.5700e-2 4.2574e-1
    h=1/8, \tau=1/64 5.6669e-3 3.07 4.4078e-2 1.86 1.7255e-2 2.31 2.0459e-1 1.06
    h=1/16, \tau=1/256 6.5948e-4 3.10 1.1349e-2 1.96 3.7745e-3 2.19 9.9988e-2 1.03
    h=1/32, \tau=1/1024 7.9662e-5 3.05 2.8614e-3 1.99 8.7440e-4 2.11 4.9501e-2 1.01

     | Show Table
    DownLoad: CSV

    Thus, Tables 14 imply the coupled and decoupled algorithms are robust on the Poisson rate \nu .

    Table 5 reports the CPU time for both the coupled and decoupled algorithms. The results indicate that the decoupled algorithm significantly outperforms the coupled one in terms of runtime, demonstrating its superior efficiency.

    Table 5.  The CPU time of coupled and decoupled algorithms for \nu = 0.49 .
    h, \tau the CPU time of coupled algorithm the CPU time of decoupled algorithm
    h=1/4, \tau=1/16 0.19 s 0.13 s
    h=1/8, \tau=1/64 2.61 s 1.63 s
    h=1/16, \tau=1/256 51.82 s 30.32 s
    h=1/32, \tau=1/1024 915.23 s 503.59 s

     | Show Table
    DownLoad: CSV

    We now test the robustness of the algorithms with respect to the permeability tensor K by fixing \alpha = 1 , \mu_{f} = 1 , E = 2.8\text{e}3 , \nu = 0.49 , and c_{0} = 0.2 for K = 1\mathbf{I} and K = 1 e- 4\mathbf{I} , respectively.

    For the cases of K = 1\mathbf{I} and K = 1 e- 4\mathbf{I} , Tables 6 and 7 exhibit the spatial errors and convergence rates of \mathbf{u} and p for the coupled algorithm. The convergence rates of \mathbf{u} with respect to H^{1} norm and L^{2} norm are all around 3 and 2 separately. The convergence rates of p with respect to H^{1} norm and L^{2} norm are all around 2 and 1, respectively. This implies that the coupled algorithm is robust for K .

    Table 6.  Spatial errors and convergence rates of \mathbf{u} and p of the coupled algorithm for K = 1\mathbf{I} .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 1.3174e-1 3.7941e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 3.5460e-2 1.89 1.9461e-1 0.96
    h=1/16, \tau=1/256 6.5936e-4 3.10 1.1349e-2 1.96 9.0399e-3 1.97 9.7954e-2 0.99
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 2.2712e-3 1.99 4.9060e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 7.  Spatial errors and convergence rates of \mathbf{u} and p of the coupled algorithm for K = 1 e- 4\mathbf{I} .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 8.5753e-2 4.2449e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 1.7287e-2 2.31 2.0376e-1 1.06
    h=1/16, \tau=1/256 6.5936e-4 3.10 1.1349e-2 1.96 3.8309e-3 2.17 9.9511e-2 1.03
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 1.0018e-3 1.93 4.9279e-2 1.01

     | Show Table
    DownLoad: CSV

    Again, for the cases of K = 1\mathbf{I} and K = 1 e- 4\mathbf{I} , we compute the convergence rates of \mathbf{u} and p by the decoupled algorithm. Tables 8 and 9 show that the convergence rates of \mathbf{u} and p are optimal, which is consistent with the results in Theorem 5 and shows that the decoupled algorithm is robust for K .

    Table 8.  Spatial errors and convergence rates of \mathbf{u} and p of the decoupled algorithm for K = 1\mathbf{I} .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 1.3174e-1 3.7941e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 3.5460e-2 1.89 1.9461e-1 0.96
    h=1/16, \tau=1/256 6.5935e-4 3.10 1.1349e-2 1.96 9.0399e-3 1.97 9.7954e-2 0.99
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 2.2712e-3 1.99 4.9060e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 9.  Spatial errors and convergence rates of \mathbf{u} and p of the decoupled algorithm for K = 1 e- 4\mathbf{I} .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/16 4.7630e-2 1.6047e-1 8.5753e-2 4.2449e-1
    h=1/8, \tau=1/64 5.6649e-3 3.07 4.4075e-2 1.86 1.7287e-2 2.31 2.0376e-1 1.06
    h=1/16, \tau=1/256 6.5935e-4 3.10 1.1349e-2 1.96 3.8309e-3 2.17 9.9511e-2 1.03
    h=1/32, \tau=1/1024 7.9656e-5 3.05 2.8614e-3 1.99 1.0018e-3 1.94 4.9279e-2 1.01

     | Show Table
    DownLoad: CSV

    We now verify that there is no nonphysical oscillation in the pressure field for the proposed algorithms. For brevity, we only consider the decoupled algorithm. The benchmark problem for the poroelasticity model is given as follows:

    Example 2. (Barry-Mercer problem) Again we take \Omega = (0, 1)^2 and T = 1 . The boundary segments \Gamma_{j}(j = 1, 2, 3, 4) are defined in Example 1. The source functions \mathbf{f} = 0 and \phi = 0 , and the boundary and initial conditions are

    \begin{alignat*} {2} p & = 0 &&\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 1,3,4,\\ p& = p_{1} &&\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 2,\\ u_1 & = 0 &&\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 2,4,\\ u_2 & = 0 &&\qquad\mbox{on }\Gamma_j\times (0,T),\, j = 1,3,\\ \sigma(\mathbf{u})\mathbf{n}-\alpha p \mathbf{n} & = \mathbf{f}_1: = (0,\alpha p)^{T} &&\qquad \mbox{on } \partial\Omega_T,\\ \mathbf{u}(\mathbf{x},0) = \mathbf{0}, \quad p(\mathbf{x},0) & = 0 &&\qquad\mbox{in } \Omega, \end{alignat*}

    where

    \begin{align*} &p_{1} = \left\lbrace \begin{aligned} &\sin t &&{\rm when}\; ( x_{1},t)\in[0.2,0.8)\times(0,T),\\ &0&&{\rm others}. \end{aligned}\right. \end{align*}

    Set \alpha = 1 , \mu_{f} = 1 , E = 2.8 e 3 , \nu = 0.4 , K = 1 e- 8\mathbf{I} , and c_{0} = 0 . For the case of c_{0} = 0 , low permeability K , and small time step \tau , it will lead to \mathrm{div}\mathbf{u}^{1}\approx 0 for the discrete form of Eq (1.2). If using \mathbf{P}_2-P_1 element pairs to solve the original problem, there exists the nonphysical oscillation in the pressure field at a very short time see Figure 1. But Figure 2 shows that there is no nonphysical oscillation in the pressure field for the decoupled algorithm with \mathbf{P}_2-P_1-P_1-P_1 element pairs. In addition, the arrows near the boundary match very well with the inner arrows in Figure 3. This implies that the decoupled algorithm is locking-free and stable.

    Figure 1.  The numerical solution p_{h}^{n+1} by using the \mathbf{P}_2-P_1 element pairs for the original problem at t = 0.00001 .
    Figure 2.  The numerical solution p_{h}^{n+1} by using the decoupled algorithm with \mathbf{P}_2-P_1-P_1-P_1 element pairs at t = 0.00001 .
    Figure 3.  Arrow plot of numerical solution \mathbf{u}_h by using the decoupled algorithm with \mathbf{P}_2-P_1-P_1-P_1 element pairs at t = 0.00001 .

    We here test the convergence rates for the BDF2 algorithm by fixing \alpha = 1 , \mu_{f} = 1 , E = 2.8 e 3 , K = 1 e- 7\mathbf{I} and c_{0} = 0.2 for \nu = 0.49 and \nu = 0.4999999 , respectively.

    Tables 10 and 11 show the spatial errors and convergence rates of the BDF2 algorithm for \nu = 0.49 and \nu = 0.4999999 , respectively. The numerical results show that the H^{1} convergence rates of \mathbf{u} and the L^{2} convergence rates of p are all approaching 2 . The L^{2} convergence rates of \mathbf{u} are around 3 . The H^{1} convergence rates of p are around 1 . This implies that the BDF2 algorithm has the optimal convergence rates for the different values of \nu .

    Table 10.  Spatial errors and convergence rates of \mathbf{u} and p of the BDF2 algorithm for \nu = 0.49 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/4 4.8482e-2 1.6013e-1 5.7104e-2 4.0384e-1
    h=1/8, \tau=1/8 5.7075e-3 3.09 4.4041e-2 1.86 1.3284e-2 2.10 1.9872e-1 1.02
    h=1/16, \tau=1/16 6.6098e-4 3.11 1.1349e-2 1.96 3.2745e-3 2.02 9.8508e-2 1.01
    h=1/32, \tau=1/32 7.9759e-5 3.05 2.8615e-3 1.99 9.3353e-4 1.81 4.9128e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 11.  Spatial errors and convergence rates of \mathbf{u} and p of the BDF2 algorithm for \nu = 0.4999999 .
    h, \tau \dfrac{\|e_{\mathbf{u}}\|_{L^2(\Omega)}}{\|\mathbf{u}\|_{L^2(\Omega)}} Rate \dfrac{\|e_{\mathbf{u}}\|_{H^1(\Omega)}}{\|\mathbf{u}\|_{H^1(\Omega)}} Rate \dfrac{\|e_{p}\|_{L^2(\Omega)}}{\|p\|_{L^2(\Omega)}} Rate \dfrac{\|e_{p}\|_{H^1(\Omega)}}{\|p\|_{H^1(\Omega)}} Rate
    h=1/4, \tau=1/4 4.8493e-2 1.6012e-1 5.7040e-2 4.0369e-1
    h=1/8, \tau=1/8 5.7096e-3 3.09 4.4040e-2 1.86 1.3273e-2 2.10 1.9874e-1 1.02
    h=1/16, \tau=1/16 6.6116e-4 3.11 1.1349e-2 1.96 3.2409e-3 2.03 9.8517e-2 1.01
    h=1/32, \tau=1/32 7.9771e-5 3.05 2.8615e-3 1.99 8.0517e-4 2.01 4.9130e-2 1.00

     | Show Table
    DownLoad: CSV

    In this section, the coupled and decoupled algorithms with \mathbf{P}_2-P_1-P_1-P_1 element pairs and BDF2 algorithm with \mathbf{P}_2-P_1-P_1 element pairs are used to simulate brain swelling caused by brain injury. Since the coupled and decoupled algorithms present almost the same result, we only report the simulation result based on the coupled algorithm. For more details about \Omega, \; \Gamma = \Gamma_{1}+\Gamma_{2} , source functions, boundary conditions, and the values of parameters, one can see [15]. Here we list the required information, which is obtained from [15]. The \Omega is a 2D cross-section of a 3D model of the brain (see Figure 4), and the values of parameters are shown in Table 12. The injured area is marked with a yellow line in the right-hand figure of Figure 4.

    Figure 4.  A slice of the magnetic resonance imaging for a human brain (left) and the finite-element mesh (right).
    Table 12.  Values of parameters.
    Parameters Values Parameters Values
    c_{0} 4.5e-7 {\rm Pa}^{-1} K 1.4e-9 {\rm mm}^{2}
    c_{b} 3e-5 mm/min/Pa \alpha 1
    p_{\text{ SAS}} 1070 Pa \nu 0.35
    \mu_{f} 1.48e-5 {\rm Pa\cdot min} E 9010 Pa

     | Show Table
    DownLoad: CSV

    Moreover, the source functions are \mathbf{f} = \mathbf{0}, \; \phi = 9 e- 3\; \text{ml}/{\rm min} , and the boundary conditions are

    \begin{alignat*} {2} \mathbf{u} & = \mathbf{0} &&\qquad\mbox{on }\Gamma_1,\\ (K\nabla p)\cdot\mathbf{n}& = c_{b}(p_{\text{ SAS}}-p) &&\qquad\mbox{on }\Gamma_1,\\ p & = 1100\; {\rm Pa} &&\qquad\mbox{on }\Gamma_{2},\\ (\sigma-\alpha p\mathbf{I})\mathbf{n} & = -p\mathbf{n} &&\qquad \mbox{on } \Gamma_{2}. \end{alignat*}

    Figures 5 and 6 show the simulation results of the coupled and BDF2 algorithms for an injured brain with respect to \nu = 0.35, \; 0.49 . From Figure 5, one can see that the maximum pressures of the coupled algorithm are 3025.17 Pa and 3025.69 Pa separately, which are consistent with that in [27]. And Figure 6 illustrates the similar facts for the BDF2 algorithm by the values of maximum pressures ( 3028.13 Pa and 3025.87 Pa). In addition, the brain tissue can deform due to the influence of the total stress. In our simulation, the values of maximum tissue deformation are 0.99342 mm and 0.048082 mm for the coupled algorithm and 1.0018 mm and 0.048174 mm for the BDF2 algorithm, respectively. It is easy to find that the change of \nu has little influence on the maximum pressure, and the maximum tissue deformation becomes smaller as \nu gets bigger. When the Poisson ratio is approaching 0.5 , the brain tissue is nearly incompressible such that the tissue deformation is very small. Moreover, from Table 13, the execution time of the CPU of the coupled algorithm is longer than the BDF2 algorithm; thus the latter has higher computational efficiency.

    Figure 5.  Simulation results of the coupled algorithm with \mathbf{P}_2-P_1-P_1-P_1 element pairs.
    Figure 6.  Simulation results of the BDF2 algorithm with \mathbf{P}_2-P_1-P_1 element pairs.
    Table 13.  Execution time of CPU.
    \nu Coupled algorithm BDF2 algorithm
    0.35 1547.65s 1239.13s
    0.496 126.97s 87.21s

     | Show Table
    DownLoad: CSV

    In this paper, we transform the poroelasticity model into a four-field problem by introducing the new variables and discuss the inf-sup condition about the four-field problem. Besides, we design the coupled and decoupled algorithms for the four-field problem and perform the error estimates. Optimal convergence is obtained in space and time. It can be observed that there is no constraint on mesh size and time step in error estimates. Then we introduce the BDF2 algorithm and perform stability analysis and error estimates. Furthermore, some numerical examples are presented to verify the theoretical results. More importantly, we successfully apply the algorithms to simulate the brain edema and observe the influence on the maximum tissue deformation and maximum pressure.

    Peizhen Wang: Conceptualization, investigation, writing-original draft; Wenlong He: Validation; Yajuan Di: Writing-review, editing. All authors have read and approved the final version of the manuscript for publication.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This work was supported by the National Natural Science Foundation of China (Grant No.11901197) and the Young Teacher Foundation of Henan Province (Grant No.2021GGJS080). The numerical simulations have been done on the supercomputing system in the Supercomputing Center of Wuhan University.

    The authors declare that they have no competing interests.



    [1] M. Biot, Theory of elasticity and consolidation for a porous anisotropic media, J. Appl. Phys., 26 (1955), 182–185. https://doi.org/10.1063/1.1721956 doi: 10.1063/1.1721956
    [2] W. Pao, R. Lewis, I. Masters, A fully coupled hydro-thermo-poro-mechanical model for black oil reservoir simulation, Int. J. Numer. Anal. Met., 25 (2001), 1229–1256. https://doi.org/10.1002/nag.174 doi: 10.1002/nag.174
    [3] D. Gawin, P. Baggio, B. Schrefler, Coupled heat, water and gas flow in deformable porous media, Int. J. Numer. Meth. Fl., 20 (1995), 969–978. https://doi.org/10.1002/fld.1650200817 doi: 10.1002/fld.1650200817
    [4] J. Hudson, O. Stephansson, J. Andersson, C. Tsang, L. Ling, Coupled T–H–M issues related to radioactive waste repository design and performance, Int. J. Rock Mech. Min., 38 (2001), 143–161. https://doi.org/10.1016/S1365-1609(00)00070-8 doi: 10.1016/S1365-1609(00)00070-8
    [5] D. Nemec, J. Levec, Flow through packed bed reactors: 1. single-phase flow, Int. J. Chem. Eng., 60 (2005), 6947–6957. https://doi.org/10.1016/j.ces.2005.05.068 doi: 10.1016/j.ces.2005.05.068
    [6] O. Coussy, Poromechanics, Wiley & Sons, 2004.
    [7] M. Doi, S. Edwards, The theory of polymer dynamics, Oxford: Clarendon Press, 1986.
    [8] R. Showalter, Diffusion in poro-elastic media, J. Math. Anal. Appl., 251 (2000), 310–340. https://doi.org/10.1006/jmaa.2000.7048 doi: 10.1006/jmaa.2000.7048
    [9] P. Phillips, M. Wheeler, A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous in time case, Computat. Geosci., 11 (2007), 131–144. https://doi.org/10.1007/s10596-007-9045-y doi: 10.1007/s10596-007-9045-y
    [10] P. Phillips, M. Wheeler, A coupling of mixed and continuous Galerkin fnite element methods for poroelasticity II: The discrete-in-time case, Computat. Geosci., 11 (2007), 145–158. https://doi.org/10.1007/s10596-007-9044-z doi: 10.1007/s10596-007-9044-z
    [11] P. Phillips, M. Wheeler, A coupling of mixed and discontinuous Galerkin finite-element methods for poroelasticity, Computat. Geosci., 12 (2008), 417–435. https://doi.org/10.1007/s10596-008-9082-1 doi: 10.1007/s10596-008-9082-1
    [12] S. Yi, A coupling of nonconforming and mixed finite element methods for Biot's consolidation model, Numer. Meth. Part. D. E., 9 (2013), 1949–1777. https://doi.org/10.1002/num.21775 doi: 10.1002/num.21775
    [13] X. Hu, C. Rodrigo, F. Gaspar, L. Zikatanov, A nonconforming finite element method for the Biot's considation model in poroelasticity, J. Comput. Appl. Math., 310 (2017), 143–154. https://doi.org/10.1016/j.cam.2016.06.003 doi: 10.1016/j.cam.2016.06.003
    [14] C. Rodrigo, F. Gaspar, X. Hu, L. Zikatanov, Stability and montonicity for some discretizations of the Biot's considation model, Comput. Method. Appl. M., 298 (2016), 183–204. https://doi.org/10.1016/j.cma.2015.09.019 doi: 10.1016/j.cma.2015.09.019
    [15] G. Ju, M. Cai, J. Li, J. Tian, Parameter-robust multiphysics algorithms for Biot model with application in brain edema simulation, Math. Comput. Simulat., 177 (2020), 385–403. https://doi.org/10.1016/j.matcom.2020.04.027 doi: 10.1016/j.matcom.2020.04.027
    [16] S. Yi, A study of two models of locking in poroelasticity, SIAM J. Numer. Anal., 55 (2017), 1915–1936. https://doi.org/10.1137/16M1056109 doi: 10.1137/16M1056109
    [17] X. Feng, Z. Ge, Y. Li, Analysis of a multiphysics finite element method for a poroelasticity model, IMA J. Numer. Anal., 38 (2018), 330–359. https://doi.org/10.1093/imanum/drx003 doi: 10.1093/imanum/drx003
    [18] J. Korsawe, G. Starke, A least-squares mixed finite element method for Biot's consolidation problem in porous media, SIAM J. Numer. Anal., 43 (2005), 318–339. https://doi.org/10.1137/S0036142903432929 doi: 10.1137/S0036142903432929
    [19] J. Korsawe, G. Starke, W. Wang, O. Kolditz, Finite element analysis of poro-elastic consolidation in porous media: Standard and mixed approaches, Comput. Method. Appl. M., 195 (2006), 1096–1115. https://doi.org/10.1016/j.cma.2005.04.011 doi: 10.1016/j.cma.2005.04.011
    [20] M. Murad, A. Loula, On stability and convergence of finite element approximations of Biot's consolidation problem, Int. J. Numer. Meth. Eng., 37 (1994), 645–667. https://doi.org/10.1002/nme.1620370407 doi: 10.1002/nme.1620370407
    [21] M. Murad, V. Thomée, A. Loula, Asymptotic behavior of semidiscrete finite-element approximations of Biot's consolidation problem, SIAM J. Numer. Anal., 33 (1996), 1065–1083. https://doi.org/ 10.1137/0733052 doi: 10.1137/0733052
    [22] S. Yi, M. Bean, Iteratively coupled solution strategies for a four-field mixed finite element method for poroelasticity, Int. J. Numer. Anal. Met., 41 (2017), 159–179. https://doi.org/10.1002/nag.2538 doi: 10.1002/nag.2538
    [23] S. Yi, Convergence analysis of a new mixed finite element method for Biot's consolidation model, Numer. Meth. Part. D. E., 30 (2014), 1189–1210. https://doi.org/10.1002/num.21865 doi: 10.1002/num.21865
    [24] A. Naumovich, On finite volume discretization of the three-dimensional biot poroelasticity system in multilayer domains, Comput. Method. Appl. M., 6 (2006), 306–325. https://doi.org/10.2478/cmam-2006-0017 doi: 10.2478/cmam-2006-0017
    [25] H. Peng, W. Qi, Weak Galerkin finite element method with the total pressure variable for Biot's consolidation model, Appl. Numer. Math., 207 (2025), 450–469. https://doi.org/10.1016/j.apnum.2024.09.017 doi: 10.1016/j.apnum.2024.09.017
    [26] S. Gu, S. Chai, C. Zhou, J. Zhou, Weak Galerkin finite element method for linear poroelasticity problems, Appl. Numer. Math., 190 (2023), 200–219. https://doi.org/10.1016/j.apnum.2023.04.015 doi: 10.1016/j.apnum.2023.04.015
    [27] X. Li, H. Holst, S. Kleiven, Influence of gravity for optimal head positions in the treatment of head injury patients, Acta Neurochir., 153 (2011), 2057–2064. https://doi.org/10.1007/s00701-011-1078-2 doi: 10.1007/s00701-011-1078-2
    [28] S. Brenner, L. R. Scott, The mathematical theory of finite element methods, 3 Eds., Springer, 2008. http://dx.doi.org/10.1016/B978-0-12-775850-3.50017-0
    [29] P. Ciarlet, The finite element method for elliptic problems, North-Holland, Amsterdam, 1978. http://doi.org/10.1016/0378-4754(80)90078-6
    [30] R. Temam, Navier-Stokes equations, Studies in Mathematics and its Applications, North-Holland, 33 (1977).
    [31] V. Girault, P. Raviart, Finite element method for Navier-Stokes equations: Theory and algorithms, Springer Science & Business Media, 5 (2012).
    [32] R. Dautray, J. Lions, Mathematical analysis and numerical methods for science and technology, Springer Verlag, 1 (1990).
    [33] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, New York: Springer, 1992.
    [34] M. Bercovier, O. Pironneau, Error estimates for finite element solution of the Stokes problem in the primitive variables, Numer. Math., 33 (1979), 211–224. https://doi.org/10.1007/bf01399555 doi: 10.1007/bf01399555
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(238) PDF downloads(25) Cited by(0)

Figures and Tables

Figures(6)  /  Tables(13)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog