Processing math: 34%
Research article Special Issues

Regularity of all minimizers of a class of spectral partition problems

  • We study a rather broad class of optimal partition problems with respect to monotone and coercive functional costs that involve the Dirichlet eigenvalues of the partitions. We show a sharp regularity result for the entire set of minimizers for a natural relaxed version of the original problem, together with the regularity of eigenfunctions and a universal free boundary condition. Among others, our result covers the cases of the following functional costs (ω1,,ωm)mi=1(kij=1λj(ωi)pi)1/pi,mi=1(kij=1λj(ωi)),mi=1(kij=1λj(ωi)) where (ω1,,ωm) are the sets of the partition and λj(ωi) is the j-th Laplace eigenvalue of the set ωi with zero Dirichlet boundary conditions.

    Citation: Hugo Tavares, Alessandro Zilio. Regularity of all minimizers of a class of spectral partition problems[J]. Mathematics in Engineering, 2021, 3(1): 1-31. doi: 10.3934/mine.2021002

    Related Papers:

    [1] Heesung Shin, Jiang Zeng . More bijections for Entringer and Arnold families. Electronic Research Archive, 2021, 29(2): 2167-2185. doi: 10.3934/era.2020111
    [2] Bin Han . Some multivariate polynomials for doubled permutations. Electronic Research Archive, 2021, 29(2): 1925-1944. doi: 10.3934/era.2020098
    [3] Shishuo Fu, Zhicong Lin, Yaling Wang . Refined Wilf-equivalences by Comtet statistics. Electronic Research Archive, 2021, 29(5): 2877-2913. doi: 10.3934/era.2021018
    [4] Jiafan Zhang . On the distribution of primitive roots and Lehmer numbers. Electronic Research Archive, 2023, 31(11): 6913-6927. doi: 10.3934/era.2023350
    [5] Dmitry Krachun, Zhi-Wei Sun . On sums of four pentagonal numbers with coefficients. Electronic Research Archive, 2020, 28(1): 559-566. doi: 10.3934/era.2020029
    [6] Massimo Grossi . On the number of critical points of solutions of semilinear elliptic equations. Electronic Research Archive, 2021, 29(6): 4215-4228. doi: 10.3934/era.2021080
    [7] Ji-Cai Liu . Proof of Sun's conjectural supercongruence involving Catalan numbers. Electronic Research Archive, 2020, 28(2): 1023-1030. doi: 10.3934/era.2020054
    [8] Hongjian Li, Kaili Yang, Pingzhi Yuan . The asymptotic behavior of the reciprocal sum of generalized Fibonacci numbers. Electronic Research Archive, 2025, 33(1): 409-432. doi: 10.3934/era.2025020
    [9] Hanpeng Gao, Yunlong Zhou, Yuanfeng Zhang . Sincere wide τ-tilting modules. Electronic Research Archive, 2025, 33(4): 2275-2284. doi: 10.3934/era.2025099
    [10] Taboka Prince Chalebgwa, Sidney A. Morris . Number theoretic subsets of the real line of full or null measure. Electronic Research Archive, 2025, 33(2): 1037-1044. doi: 10.3934/era.2025046
  • We study a rather broad class of optimal partition problems with respect to monotone and coercive functional costs that involve the Dirichlet eigenvalues of the partitions. We show a sharp regularity result for the entire set of minimizers for a natural relaxed version of the original problem, together with the regularity of eigenfunctions and a universal free boundary condition. Among others, our result covers the cases of the following functional costs (ω1,,ωm)mi=1(kij=1λj(ωi)pi)1/pi,mi=1(kij=1λj(ωi)),mi=1(kij=1λj(ωi)) where (ω1,,ωm) are the sets of the partition and λj(ωi) is the j-th Laplace eigenvalue of the set ωi with zero Dirichlet boundary conditions.


    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 f,f1, ϕ, and ϕ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 η and problem (2.9) and (2.10) is a coupled problem associated with the diffusion problem, we can solve a generalized Stokes problem for u and ξ by using the solution of η at the previous time-step. After that, we solve the diffusion problem for η and p.

    To derive the error estimates of the MFEM, we define L2(Ω)-projection Qh:L2(Ω)Xkh by

    (Qhφ,ψh)=(φ,ψh),ψhXkh and φL2(Ω), (3.7)

    where Xkh:={ψhC0;ψh|KPk(K)KTh}, k is the degree of piecewise polynomial on K.

    Next, for any φH1(Ω), we define the elliptic projection Sh:H1(Ω)Xkh as

    (KShφ,φh)=(Kφ,φh),φhXkh, (3.8)
    (Shφ,1)=(φ,1), (3.9)

    and for any vH1(Ω), we define its elliptic projection Rh:H1(Ω)Vkh as

    (ε(Rhv),ε(wh))=(ε(v),ε(wh))whVkh, (3.10)

    where Vkh:={vhC0;vh|KPk(K),(vh,r)=0rRM}.

    It follows from [28] that for any 0sk, Qh,Sh and Rh satisfy

    QhφφL2(Ω)+h(Qhφφ)L2(Ω)Chs+1φHs+1(Ω),φHs+1(Ω), (3.11)
    ShφφL2(Ω)+h(Shφφ)L2(Ω)Chs+1φHs+1(Ω),φHs+1(Ω), (3.12)
    RhvvL2(Ω)+h(Rhvv)L2(Ω)Chs+1vHs+1(Ω),vHs+1(Ω). (3.13)

    Furthermore, we introduce the following error notations:

    Enu=u(tn)unh,Enξ=ξ(tn)ξnh,Enη=η(tn)ηnh,Enp=p(tn)pnh.

    Also, we denote

    Enu=u(tn)Rh(u(tn))+Rh(u(tn))unh:=Ynu+Znu,Eni=i(tn)Sh(i(tn))+Sh(i(tn))inh:=Yni+Zni,i=ξ,η,p,Eni=i(tn)Qh(i(tn))+Qh(i(tn))inh:=Fni+Gni,i=ξ,η,p.

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

    Theorem 4. Let {(unh,ξnh,ηnh,pnh)}n0 be solutions of the coupled algorithm (i.e., Algorithm 1). Then it holds

    max0nN[μ(u(tn)unh)L2(Ω)+κ2η(tn)ηnhL2(Ω)+κ3ξ(tn)ξnhL2(Ω)]+(τNn=01μfK12(p(tn)pnh)2L2(Ω))12C1τ+C2h2. (3.14)

    Theorem 5. Let {(unh,ξnh,ηnh,pnh)}n0 be solutions of the decoupled algorithm (i.e., Algorithm 2). Then for the case c0>0 and λ<, it holds

    max0nN[μ(u(tn)unh)L2(Ω)+κ2η(tn)ηnhL2(Ω)+κ3ξ(tn)ξnhL2(Ω)]+(τNn=01μfK12(p(tn)pnh)2L2(Ω))12ˉC1τ+ˉC2h2. (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 τ=O(h2) in the error estimates for the algorithm with θ=0 [17]. This constraint arises because the choice of the test function ψh results in the variables p,ξ, and η 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 ψh that ensures the variables p,ξ, and η reside in the same time layer. To do that, we first examine the following error equations:

    2κ1(Gn+1ξ,ωh)+κ2(Gn+1η,ωh)(Gn+1p,ωh)=0,ωhWh, (3.16)
    (dtGn+1η,ψh)+1μf(KZn+1p,ψh)=(Rn+11,η,ψh),ψhWh, (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 η1h by

    κ1(η1h,φh)=κ3(ξ0h,φh)+(divu0h,φh).

    Setting ωh=dtGnη,ψh=Znp in (3.16) and (3.17) separately after lowering the super index n+1 to n, applying the summation operator τln=0 to both sides of (3.17) and utilizing the fact that G1η=0, there holds

    κ22Glη2L2(Ω)+τln=0[κ1(dtGnη,Gn+1ξ)+κ2τ2dtGnη2L2(Ω)+1μfK12Znp2L2(Ω)]=τln=0[(dtGnη,YnpFnp)+(dtη(tn+1)ηt(tn+1),Znp)+κ1τ(dtGn+1ξ,dtGnη)]. (3.18)

    Note that the choice of ψh=Znp is critical. Unlike the setting in [17], where ψh=ˆZnp:=κ1Zn+1ξ+κ2Znη is employed, selecting ψh=Znp 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

    κ1τ2ln=0(dtGn+1ξ,dtGnη)κ1τ2ln=0[12ϵdtGn+1ξ2L2(Ω)+ϵ2dtGnη)2L2(Ω)], (3.19)

    where ϵ(2κ2κ3,+) is a positive constant from Young inequality.

    Additionally, we point out that:

    (i) The aim of ϵ(2κ2κ3,+) is to ensure the signs of terms κ2κ1ϵ2τ2ln=0dtGnη2L2(Ω) and (κ34κ12ϵ)τ2ln=0dtGnη2L2(Ω) to be positive. In addition, for the case of λ+, there does not exist the term κ1τ2ln=0(dtGn+1ξ,dtGnη) in error estimates. For the case of c0=0, it has a similar discussion for the term κ1τ2ln=0(dtGn+1ξ,dtGnη).

    (ii) For the decoupled algorithm, it does not exist the constraint on the time step τ 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

    vn+12:=(vn+1+vn)/2andDtvn=3vn4vn1+vn22τ. (4.1)

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

    Algorithm 3 BDF2 algorithm
    Input: Evaluate u0hVh,p0hMh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      Step 1: Solve for (u1h,ξ1h,η1h)Vh×Mh×Wh such that
            2μ(ε(u12h),ε(vh))(ξ12h,divvh)=(f,vh)+f1,vh,vhVh,(4.2)
            κ3(ξ12h,φh)+(divu12h,φh)=κ1(η12h,φh),φhMh,(4.3)
            (dtη1h,ψh)+1μf(K(κ1ξ12h+κ2η12h),ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh.(4.4)
      Step 2:
      for n=1,2,, do
        (i) Solve for (un+1h,ξn+1h,ηn+1h)Vh×Mh×Wh such that
            2μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,(4.5)
            κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηn+1h,φh),φhMh,(4.6)
            (Dtηn+1h,ψh)+1μf(K(κ1ξn+1h+κ2ηn+1h),ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh.(4.7)
        (ii) Update pn+1h and qn+1h by
            pn+1h=κ1ξn+1h+κ2ηn+1h,qn+1h=κ1ηn+1hκ3ξn+1h,(4.8)
      end for

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

    Lemma 4.1. Let {(unh,ξnh,ηnh)}n0 be defined by the BDF2 algorithm. Then it holds

    μ8ε(ul+1h)2L2(Ω)+κ34ξl+1h2L2(Ω)+κ24ηl+1h2L2(Ω)+τ2μfln=1K12(κ1ξn+1h+κ2ηn+1h)2L2(Ω)C3(f2L2(Ω)+f12L2(Ω)+ε(u0h)2L2(Ω)+ξ0h2L2(Ω)+η0h2L2(Ω)+ϕ2L2(Ω)+ϕ12L2(Ω)),l1, (4.9)

    where C3=C3(ck,ca,cb,cp,μ,κ3,κ2,T,K,μf).

    Proof. Setting vh=dtu1h in (4.2), φh=dtξ1h in (4.3) and ψh=p12h in (4.4), we get

    2μ(ε(u12h),ε(dtu1h))(ξ12h,divdtu1h)=(f,dtu1h)+f1,dtu1h, (4.10)
    κ3(ξ12h,dtξ1h)+(divu12h,dtξ1h)=κ1(η12h,dtξ1h), (4.11)
    (dtη1h,p12h)+1μf(K(κ1ξ12h+κ2η12h),p12h)=(ϕ,p12h)+ϕ1,p12h. (4.12)

    Adding (4.10)–(4.12), there holds

    μ2τ(ε(u1h)2L2(Ω)ε(u0h)2L2(Ω))+κ32τ(ξ1h2L2(Ω)ξ0h2L2(Ω))+κ22τ(η1h2L2(Ω)η0h2L2(Ω))+1μfK12(κ1ξ12h+κ2η12h)2L2(Ω)=(f,dtu1h)+f1,dtu1h+(ϕ,p12h)+ϕ1,p12h+1τ((ξ0h,divu1h)(ξ1h,divu0h))+κ1τ((η0h,ξ1h)(η1h,ξ0h)). (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:

    (f,dtu1h)1τ(5c2kμf2L2(Ω)+μ16ε(u1h)2L2(Ω)+μ4ε(u0h)2L2(Ω)), (4.14)
    f1,dtu1h1τ(5c2ac2kμf12L2(Ω)+μ16ε(u1h)2L2(Ω)+μ4ε(u0h)2L2(Ω)), (4.15)
    (ϕ,p12h)c2pμfK12ϕ2L2(Ω)+14μfK12p12h2L2(Ω), (4.16)
    ϕ1,p12hc2ac2pμfK12ϕ12L2(Ω)+14μfK12p12h2L2(Ω), (4.17)
    (ξ0h,divu1h)(ξ1h,divu0h)2c2bμξ0h2L2(Ω)+μ8ε(u1h)2L2(Ω)+κ38ξ1h2L2(Ω) (4.18)
    +2c2bκ3ε(u0h)2L2(Ω),κ1(η0h,ξ1h)κ1(η1h,ξ0h)2κ21κ3η0h2L2(Ω)+κ38ξ1h2L2(Ω)+κ24η1h2L2(Ω)+κ21κ2ξ0h2L2(Ω), (4.19)

    where ca is a positive constant related to Trace inequality.

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

    μ4ε(u1h)2L2(Ω)+κ34ξ1h2L2(Ω)+κ24η1h2L2(Ω)+τ2μfK12(κ1ξ12h+κ2η12h)2L2(Ω)5c2kμf2L2(Ω)+5c2ac2kμf12L2(Ω)+c2pμfτK12ϕ2L2(Ω)+c2ac2pμfτK12ϕ12L2(Ω)+(μ+2c2bκ3)ε(u0h)2L2(Ω)+(κ32+2c2bμ+κ21κ2)ξ0h2L2(Ω)+(κ22+2κ21κ3)η0h2L2(Ω). (4.20)

    Next, setting vh=Dtun+1h,φh=ξn+1h and ψh=pn+1h in (4.5)–(4.7), we obtain

    2μ(ε(un+1h),ε(Dtun+1h))(ξn+1h,divDtun+1h)=(f,Dtun+1h)+f1,Dtun+1h, (4.21)
    κ3(Dtξn+1h,ξn+1h)+(divDtun+1h,ξn+1h)=κ1(Dtηn+1h,ξn+1h), (4.22)
    (Dtηn+1h,pn+1h)+1μf(K(κ1ξn+1h+κ2ηn+1h),pn+1h)=(ϕ,pn+1h)+ϕ1,pn+1h. (4.23)

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

    (Dtvn+1h,vn+1h)=14τ(vn+1h2L2(Ω)vnh2L2(Ω)+2vn+1hvnh2L2(Ω)2vnhvn1h2L2(Ω)+vn+1h2vnh+vn1h2L2(Ω)), (4.24)

    we have

    μ4τ(ε(un+1h)2L2(Ω)ε(unh)2L2(Ω)+2ε(un+1h)ε(unh)2L2(Ω)2ε(unh)ε(un1h)2L2(Ω)+ε(un+1h)2ε(unh)+ε(un1h)2L2(Ω))+κ34τ(ξn+1h2L2(Ω)ξnh2L2(Ω)+2ξn+1hξnh2L2(Ω)2ξnhξn1h2L2(Ω)+ξn+1h2ξnh+ξn1h2L2(Ω))+κ24τ(ηn+1h2L2(Ω)ηnh2L2(Ω)+2ηn+1hηnh2L2(Ω)2ηnhηn1h2L2(Ω)+ηn+1h2ηnh+ηn1h2L2(Ω))+1μfK12(κ1ξn+1h+κ2ηn+1h)2L2(Ω)=(f,Dtun+1h)+f1,Dtun+1h+(ϕ,pn+1h)+ϕ1,pn+1h. (4.25)

    Applying the summation operator τln=1 to both sides of (4.25), there holds

    μ4(ε(ul+1h)2L2(Ω)+2ε(ul+1h)ε(ulh)2L2(Ω))+κ34(ξl+1h2L2(Ω)+2ξl+1hξlh2L2(Ω))+κ24(ηl+1h2L2(Ω)+2ηl+1hηlh2L2(Ω))+ln=1(μ2ε(un+1h)2ε(unh)+ε(un1h)2L2(Ω)+κ34ξn+1h2ξnh+ξn1h2L2(Ω)+κ24ηn+1h2ηnh+ηn1h2L2(Ω)+τμfK12(κ1ξn+1h+κ2ηn+1h)2L2(Ω))=(f,3ul+1hulh3u1h+u0h2)+f1,3ul+1hulh3u1h+u0h2+τln=1[(ϕ,pn+1h)+ϕ1,pn+1h]+μ4(ε(u1h)2L2(Ω)+2ε(u1h)ε(u0h)2L2(Ω))+κ34(ξ1h2L2(Ω)+2ξ1hξ0h2L2(Ω))+κ24(η1h2L2(Ω)+2η1hη0h2L2(Ω)). (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] Alper O (2020) On the singular set of free interfacein an optimal partition problem. Commun Pure Appl Math 73: 855-915.
    [2] Alt H, Caffarelli L, Friedman A (1984) Variational problems with two phases and their free boundaries. T Am Math Soc 282: 431-461.
    [3] Band R, Berkolaiko G, Raz H, et al.(2012) The number of nodal domains on quantum graphs as a stability index of graph partitions. Commun Math Phys 311: 815-838.
    [4] Berkolaiko G, Kuchment P, Smilansky U (2012) Critical partitions and nodal deficiency of billiard eigenfunctions. Geom Funct Anal 22: 1517-1540.
    [5] Bonnaillie-Noël V, Helffer B, Vial G (2010) Numerical simulations for nodal domains and spectral minimal partitions. ESAIM Contr Optim Ca 16: 221-246.
    [6] Bourdin B, Bucur D, Oudet E (2009) Optimal partitions for eigenvalues. SIAM J Sci Comput 31: 4100-4114.
    [7] Bucur D, Buttazzo G (2005) Variational Methods in Shape Optimization Problems, Boston: Birkhäuser Boston Inc.
    [8] Bucur D, Buttazzo G, Henrot A (1998) Existence results for some optimal partition problems. Adv Math Sci Appl 8: 571-579.
    [9] Caffarelli L, Lin FH (2007) An optimal partition problem for eigenvalues. J Sci Comput 31: 5-18.
    [10] Caffarelli L, Lin FH (2010) Analysis on the junctions of domain walls. Discrete Contin Dyn Syst 28: 915-929.
    [11] Chang SM, Lin CS, Lin TC, et al (2004) Segregated nodal domains of two-dimensional multispecies Bose-Einstein condensates. Phys D 196: 341-361.
    [12] Conti M, Terracini S, Verzini G (2002) Nehari's problem and competing species systems. Ann I H Poincaré Anal Non Linéaire 19: 871-888.
    [13] Conti M, Terracini S, Verzini G (2003) An optimal partition problem related to nonlinear eigenvalues. J Funct Anal 198: 160-196.
    [14] Conti M, Terracini S, Verzini G (2005) On a class of optimal partition problems related to the Fučík spectrum and to the monotonicity formulae. Calc Var Partial Dif 22: 45-72.
    [15] De Philippis G, Lamboley G, Pierre M, et al.(2018) Regularity of minimizers of shape optimization problems involving perimeter. J Math Pure Appl 109: 147-181.
    [16] Hecht F (2012) New development in freefem++. J Numer Math 20: 251-265.
    [17] Helffer B, Hoffmann-Ostenhof T, Terracini S (2009) Nodal domains and spectral minimal partitions. Ann I H Poincaré Anal Non Linéaire 26: 101-138.
    [18] Helffer B, Hoffmann-Ostenhof T, Terracini S (2010) Nodal minimal partitions in dimension 3. Discrete Contin Dyn Syst 28: 617-635.
    [19] Noris B, Tavares H, Terracini S, et al. (2010) Uniform Hölder bounds for nonlinear Schrödinger systems with strong competition. Commun Pure Appl Math 63: 267-302.
    [20] Ramos M, Tavares H, Terracini S (2016) Extremality conditions and regularity of solutions to optimal partition problems involving Laplacian eigenvalues. Arch Ration Mech Anal 220: 363-443.
    [21] Soave N, Tavares H, Terracini S, et al. (2016) Hölder bounds and regularity of emerging free boundaries for strongly competing Schrödinger equations with nontrivial grouping. Nonlinear Anal 138: 388-427.
    [22] Soave N, Terracini S (2015) Liouville theorems and 1-dimensional symmetry for solutions of an elliptic system modelling phase separation. Adv Math 279: 29-66.
    [23] Tavares H, Terracini S (2012) Regularity of the nodal set of segregated critical configurations under a weak reflection law. Calc Var Partial Dif 45: 273-317.
    [24] Tavares H, Terracini S (2012) Sign-changing solutions of competition-diffusion elliptic systems and optimal partition problems. Ann I H Poincaré Anal Non Linéaire 29: 279-300.
    [25] Terracini S, Verzini G, Zilio A (2016) Uniform Hölder bounds for strongly competing systems involving the square root of the laplacian. J Eur Math Soc 18: 2865-2924.
  • This article has been cited by:

    1. Shishuo Fu, Jiaxi Lu, Yuanzhe Ding, A skeleton model to enumerate standard puzzle sequences, 2021, 30, 2688-1594, 179, 10.3934/era.2022010
    2. Kanasottu Anil Naik, Rayappa David Amar Raj, Chepuri Venkateswara Rao, Thanikanti Sudhakar Babu, Generalized cryptographic image processing approaches using integer-series transformation for solar power optimization under partial shading, 2022, 272, 01968904, 116376, 10.1016/j.enconman.2022.116376
    3. Sen-Peng Eu, Tung-Shan Fu, Springer Numbers and Arnold Families Revisited, 2024, 10, 2199-6792, 125, 10.1007/s40598-023-00230-9
  • Reader Comments
  • © 2021 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(4574) PDF downloads(648) Cited by(0)

Figures and Tables

Figures(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog