Loading [MathJax]/jax/output/SVG/jax.js
Special Issues

Numerical analysis of modular grad-div stability methods for the time-dependent Navier-Stokes/Darcy model

  • In this paper, we construct a modular grad-div stabilization method for the Navier-Stokes/Darcy model, which is based on the first order Backward Euler scheme. This method does not enlarge the accuracy of numerical solution, but also can improve mass conservation and relax the influence of parameters. Herein, we give stability analysis and error estimations. Finally, by some numerical experiment, the scheme our proposed is shown to be valid.

    Citation: Jiangshan Wang, Lingxiong Meng, Hongen Jia. Numerical analysis of modular grad-div stability methods for the time-dependent Navier-Stokes/Darcy model[J]. Electronic Research Archive, 2020, 28(3): 1191-1205. doi: 10.3934/era.2020065

    Related Papers:

    [1] Jiangshan Wang, Lingxiong Meng, Hongen Jia . Numerical analysis of modular grad-div stability methods for the time-dependent Navier-Stokes/Darcy model. Electronic Research Archive, 2020, 28(3): 1191-1205. doi: 10.3934/era.2020065
    [2] Linlin Tan, Bianru Cheng . Global well-posedness of 2D incompressible Navier–Stokes–Darcy flow in a type of generalized time-dependent porosity media. Electronic Research Archive, 2024, 32(10): 5649-5681. doi: 10.3934/era.2024262
    [3] Jiwei Jia, Young-Ju Lee, Yue Feng, Zichan Wang, Zhongshu Zhao . Hybridized weak Galerkin finite element methods for Brinkman equations. Electronic Research Archive, 2021, 29(3): 2489-2516. doi: 10.3934/era.2020126
    [4] Guoliang Ju, Can Chen, Rongliang Chen, Jingzhi Li, Kaitai Li, Shaohui Zhang . Numerical simulation for 3D flow in flow channel of aeroengine turbine fan based on dimension splitting method. Electronic Research Archive, 2020, 28(2): 837-851. doi: 10.3934/era.2020043
    [5] Jie Qi, Weike Wang . Global solutions to the Cauchy problem of BNSP equations in some classes of large data. Electronic Research Archive, 2024, 32(9): 5496-5541. doi: 10.3934/era.2024255
    [6] Jie Zhang, Gaoli Huang, Fan Wu . Energy equality in the isentropic compressible Navier-Stokes-Maxwell equations. Electronic Research Archive, 2023, 31(10): 6412-6424. doi: 10.3934/era.2023324
    [7] Zhonghua Qiao, Xuguang Yang . A multiple-relaxation-time lattice Boltzmann method with Beam-Warming scheme for a coupled chemotaxis-fluid model. Electronic Research Archive, 2020, 28(3): 1207-1225. doi: 10.3934/era.2020066
    [8] Guochun Wu, Han Wang, Yinghui Zhang . Optimal time-decay rates of the compressible Navier–Stokes–Poisson system in $ \mathbb R^3 $. Electronic Research Archive, 2021, 29(6): 3889-3908. doi: 10.3934/era.2021067
    [9] Wei Shi, Xinguang Yang, Xingjie Yan . Determination of the 3D Navier-Stokes equations with damping. Electronic Research Archive, 2022, 30(10): 3872-3886. doi: 10.3934/era.2022197
    [10] Xiaoxia Wang, Jinping Jiang . The uniform asymptotic behavior of solutions for 2D g-Navier-Stokes equations with nonlinear dampness and its dimensions. Electronic Research Archive, 2023, 31(7): 3963-3979. doi: 10.3934/era.2023201
  • In this paper, we construct a modular grad-div stabilization method for the Navier-Stokes/Darcy model, which is based on the first order Backward Euler scheme. This method does not enlarge the accuracy of numerical solution, but also can improve mass conservation and relax the influence of parameters. Herein, we give stability analysis and error estimations. Finally, by some numerical experiment, the scheme our proposed is shown to be valid.



    Numerical methods of Navier-Stokes/Darcy have attracted a lot of attention. So far, a great deal of numerical methods are proposed to solve this model by virtue of different ways, such as finite element methods[12], discontinuous Galerkin finite element methods[5], two-grid methods[1,15,16], modified two-grid methods[6], partitioned time stepping method[7], characteristic stabilized finite element methods[8], mortar finite element methods [2], grad-div stabilized projection finite element method[14], modular grad-div method[11] and so on. The grad-div stabilized method is first introduced in [4], which can penalize mass conservation and improve the solution quality efficiently. Recently, the effectiveness of the grad-div stabilized method has been proved in finite element simulation of Stokes and Navier-Stokes equation[10]. However, this method leads to a singular matrix stemming from grad-div term, and the larger stabilized parameter will cause solver breakdown. As a alternative method for grad-div stabilization, the modular grad-div stabilization method is introduced in [3]. The modular grad-div stabilization method for the Stokes/Darcy model is proposed in [13]. The modular grad-div stabilization method is not only easy to implement, but also avoids the influence of large parameters on the solution as well as preserving the advantages of the grad-div stabilization method.

    In this paper, we extended the modular grad-div stabilization methods from Stokes/Darcy model to Navier-Stokes/Darcy model. Compare to Navier-Stokes model, the convection term causes some difficulties in theoretical analysis and numerical simulation. To deal with convection term, the assumption unf>0 is used[6,9]. Our proposed method not only relax the influence of parameters but also improves the mass conservation.

    The rest of this paper is organized as follows: In section 2, some notations and the time-dependent Navier-Stokes/Darcy model are introduced; In section 3, we give the modular grad-div stabilization method, and stability analysis is provided; In section 4, under some regularity assumptions imposed on the true solution, error estimates are also given; In section 5, Some numerical experiments are given to verify the theoretical result, we compared with the standard scheme, standard grad-div scheme and modular grad-div scheme in the numerical experiment; Finally some conclusions are obtained in section 6.

    The model we considered is confined in a bounded domain Ω Rd(d=2 or 3), which is decomposed into a fluid flow region Ωf and porous media flow region Ωp, see figure 1. Here, ΩfΩp=, ¯Ωf¯Ωp=¯Ω and ¯Ωf¯Ωp=Γ, Γf=ΩfΩ, Γp=ΩpΩ. In the sake of simplicity, we assume that Ωf and Ωp are smooth enough throughout this paper.

    Figure 1.  The global domain Ω.

    The time-dependent Navier-Stokes equation govern the fluid flow in Ωf is expressed as

    utνΔu+(u)u+p=f1(x,t)in Ωf×(0,T),u=0in Ωf×(0,T),u(x,0)=u0(x)in Ωf. (1)

    here u=u(x,t) is the fluid velocity filed, p=p(x,t) is the pressure, function f1 is the external force, and coefficient ν>0 is the kinetic viscosity.

    The Darcy equation govern the porous media flow in Ωp is expressed as :

    S0ϕt+up=f2(x,t)in Ωp×(0,T),up=Kϕin Ωp×(0,T),ϕ(x,0)=ϕ0(x)in Ωp. (2)

    here the first equation is the saturated flow model and the second equation is the Darcy's law. S0 is the specific mass storativity coefficient, up=up(x,t) the velocity, ϕ=ϕ(x,t) the hydraulic head, K is the hydraulic conductivity tensor. We assume that K is a positive symmetric tensor and the function f2 is a source term.

    We know that

    ϕ=z+ppρg

    is the piezometric head, where pp denotes the dynamic pressure, ρ is the height from a reference level.

    Substituting the second formula of (2) into the first formula of (2), the following Darcy equation can be obtained:

    S0ϕt(Kϕ)=f2(x,t)in Ωp×(0,T). (3)

    The interface conditions of the conservation of mass, balance of forces, and the Beavers-Joseph-Saffman condition are imposed on the interface Γ by:

    unf+upnp=0on Γ×(0,T) (4)
    pνnfunf=gϕon Γ×(0,T), (5)
    ντiunf=αgντiKτiuτi,i=1,,d1on  Γ×(0,T), (6)

    where, nf and np are the unit outward normal vectors on Ωf and Ωp, respectively, and τi,i=1,,d1, are the orthonormal tangential unit vectors on the interface Γ, g is the gravitational acceleration, α is a positive parameter depending on the properties of the porous medium and must be determined experimentally.

    For simplicity of analysis, we impose the following boundary conditions on Γf, Γp:

    u=0on Γf×(0,T)ϕ=0on Γp×(0,T) (7)

    In this paper, we assume that:

    unf>0onΓ. (8)

    The assumption is not hold for general case of Navier-Stokes/Darcy Model. But for the gentle river, the water infiltration satisfies the assumption unf>0 on the interface. It is a special case of Navier-Stokes/Darcy Model.

    Next, Hillbert spaces will be introduced:

    Hf={v(H1(Ωf))d:v=0 on Γf},Hp={ψH1(Ωp): ψ=0 on Γp},Q=L2(Ωf).

    where (,)D denotes the L2 inner produce in the domain D, with corresponding norm ||||D. In the rest of this paper, we neglect the subscript. The spaces Hf and Hp are equipped with the following norms:

    ||u||Hf=||u||L2(Ωf)=(u,u)Ωf uHf,||ϕ||Hp=||ϕ||L2(Ωp)=(ϕ,ϕ)Ωf ϕHp.

    Some discrete norms are defined as,

    ||w||2L2(0,T;Hs(Ωf,p))=ΔtNn=0||wn||2Hs(Ωf,p),||w||L(0,T;Hs(Ωf,p))=max0nN||wn||Hs(Ωf,p).

    Due to u=0, we define a trilinear form af,c(;,) as follows

    af,c(u;v,w)=((u)v,w)f+12(u,vw)f=12((u)v,w)f12((u)w,v)f+12vw,unfΓ, (9)

    under the condition (8), we have

    af,c(u;v,v)0. (10)

    Then, we have the following estimates for af,c:

    af,c(u,v,w)C1||u||||v||||w||,af,c(u,v,w)C2||u||||v||2||w||. (11)

    In addition, we recall the Poincarˊe inequality and trace inequality. There exist positive constants cp and ct which only depend on the domain Ωf and exist ˜cp and ˜ct which only depend on the domain Ωp.

    ||v||L2cp||v||f||v||L2(Γ)ct||v||12L2(Ωf)||v||12H1(Ωf),||ψ||L2˜cp||ψ||p||ψ||L2(Γ)~ct||ψ||12L2(Ωp)||ψ||12H1(Ωp),||u||d||u||,d=2,or3.

    Thus, the weak formulation of the time-dependent Naiver-Stokes/Darcy model is to find u:[0,T]Hf,p:[0,T]Q,ϕ:[0,T]Hp, such that

    (ut,v)f+af(u,v)+b(v,p)+aΓ(v,ϕ)+af,c(u;u,v)=(f1,v)fvHf,b(u,q)=0qQ,gS0(ϕt,ψ)p+ap(ϕ,ψ)aΓ(u,ψ)=g(f2,ψ)pψHp. (12)

    Where

    af(u,v)=ν(D(u),D(v))+d1i=1ΓανgτiKτi(uτi)(vτi)ds,ap(ϕ,ψ)=g(Kϕ,ψ)p,aΓ(v,ϕ)=gΓϕvnfds,af,c(u,u,v)=((u)u,v)f,b(v,p)=(p,v)f.

    Bilinear forms af and ap are continuous and coercive:

    af(u,v)C3||u||Hf||v||Hf,af(u,u)ν||u||2Hf,ap(ϕ,ψ)gλmax||ϕ||Hp||ψ||Hp,ap(ϕ,ϕ)gλmin||ϕ||2Hp (13)

    The interface coupling term aΓ satisfies the following estimates:

    |aΓ(u,ϕ)|C4||u||f||ϕ||p|aΓ(u,ϕ)|C5g2h1||u||2f+||ϕ||2p,|aΓ(u,ϕ)|C6g2h1||ϕ||2p+||u||2f (14)

    Lemma 2.1. We assume that

    f1L2(0,T;L2(Ωf)d),f2L2(0,T;L2(Ωp)d),KL(Ωp)d×d,

    and K is uniformly bounded and positive definite in Ωp, there exist two constants kmin>0, kmax>0 such that

    0<kmin|x|2Kxxkmax|x|2 xΩp.

    Furthermore, uL2(Ωp)d,ϕ0L2(Ωp), therefore any solution (u,p,ϕ)(L2(0,T;Hf)H1(0,T;L2(Ωf)d))×L2(0,T;Q)×L2(0,T;Hp) of (1)-(7) is also the solution to the equation (12). The converse of the statement is also true.

    Lemma 2.2. (Discrete Gronwall Lemma). Let Δt, H, an, bn, cn, dn be nonnegative mumbers for n0 such that for N1. If

    aN+ΔtNn=0bnΔtN1n=0dnan+ΔtNn=0cn+H

    then for all Δt>0,

    aN+ΔtNn=0bnexp(ΔtN1n=0dn)(ΔtNn=0cn+H).

    We construct τh is a quasiuniform triangulation of the domain ΩfΩp, depending on a positive parameter h>0, which is made up of triangles if d=2 or tetrahedra if d=3. Then we define the finite element subspace of Hf, Q, Hp as Hfh, Qh, Hph. We assume that the space pairs (Hfh,Qh) satisfies the discrete LBB condition: there exists a positive constant β independent of h, such that qhQh, vhHfh, vfh0

    infqhQhsupvhHfh(qh,vfh)f||qh||Q||vfh||Hfβ. (15)

    Next, we divide time interval [0, T] into: 0=t0<t1<<tN=T with Δt=titi1. This leads to tm=mΔt and T=NΔt as a uniform distribution of discrete time levels. Let (un+1h,pn+1h,ϕn+1h) denote the discrete approximation of (uh(tn+1),ph(tn+1),ϕh(tn+1)). The modular grad-div scheme our proposed can be written as:

    Algorithm 1 (The modular grad-div scheme).

    For vhHfh, qhQh, and ψhHph.

    Step1. Given unhHfh, pnhQh, and ϕnhHph, find ˆun+1h,pn+1h,ϕn+1h satisfying:

    (ˆun+1hunhΔt,vh)+af(ˆun+1h,vh)+b(vh,pn+1h)+aΓ(vh,ϕnh)+af,c(unh;ˆun+1h,vh)=(fn+11,vh)b(ˆun+1h,qh)=0gS0(ϕn+1hϕnhΔt,ψh)+ap(ϕn+1h,ψh)aΓ(unh,ψh)=(fn+12,ψh) (16)

    Step2. Given ˆun+1hHfh, find un+1hHfh satisfying:

    (un+1h,vh)+(β+γΔt)(un+1h,vh)=(ˆun+1h,vh)+β(unh,vh) (17)

    Lemma 3.1. For Algorithm 1, we can obtain the following result,

    ||ˆun+1h||=||un+1h||2+||ˆun+1hun+1h||+2γΔt||un+1h||2+β(||un+1h||2||unh||2+||(un+1hunh)||2).

    Proof. Refer to Lemma 6 of [3] for proof details.

    Theorem 3.2. (Unconditional Stability) For any N>1, the solution of the Algorithm 1, satisfy

    ||uNh||2+β||uNh||2+gs0||ϕNh||2+N1n=0(||ˆun+1hunh||2+||ˆun+1hun+1h||2+β||(un+1hunh)||2+gs0||ϕn+1hϕnh||2)+N1n=02γΔt||un+1h||2+N1n=0νΔt||ˆun+1h||2+N1n=0λmingΔ||ϕn+1h||2C(ΔtN1n=0(2c2pν||fn+11||2+2c2pλming||fn+12||2)+||u0h||2+β||u0h||2+gs0||ϕ0h||2),

    Where the constant C=exp(ΔtN1n=0max2C5ghλmin,2C6g2hν).

    Proof. Taking vh=2Δˆun+1h, qh=2Δpn+1h and ψh=2Δϕn+1h in (16), using the property (10), and we can obtain af,c(unh,ˆun+1h,ˆun+1h)0. Then seeing [13] Theorem 3.1 for a detail proof.

    In this section, we will give some error estimates of our proposed method. Denote un,pn and ϕn be the true solution at time tn=nΔt. Assuming that the true solution have the following regularities,

    uL(0,T;HfHk+1(Ωf)d),utL(0,T;Hk+1(Ωf)d),uttL2(0,T;L2(Ωf)d).pL2(0,T;QHk(Ωf)).ϕL(0,T;HpHk+1(Ωp)),ϕtL(0,T;Hk+1(Ωp)),ϕttL2(0,T;L2(Ωp)) (18)

    Then define a projection operator [12]:

    Ph:(u(t),p(t),ϕ(t))Hf×Q×Hp(Phu(t),Php(t),Phϕ(t))Hfh×Qh×Hph.

    when some regularity conditions on (u(t),p(t),ϕ(t)) are statisfied, (Phu(t),Php(t),Phϕ(t)) is an approximation of (u(t),p(t),ϕ(t)), with the following properties:

    ||Phu(t)u(t)||fChk+1||u(t)||Hk+1(Ωf),||(Phu(t)u(t))||fChk||u(t)||Hk+1(Ωf),||Php(t)p(t)||fChk+1||p(t)||Hk+1(Ωf),||Phϕ(t)ϕ(t)||pChk+1||ϕ(t)||Hk+1(Ωp),||(Phϕ(t)ϕ(t))||pChk||ϕ(t)||Hk+1(Ωp) (19)

    Next, we define the following error equations

    enu=ununh=(unPhun)(unhPhun)=ηnuθnu,enˆu=unˆunh=(unPhun)(ˆunhPhun)=ηnuθnˆuenp=pnpnh=(pnPhpn)(pnhPhpn)=ηnpθnpenϕ=ϕnϕnh=(ϕnPhϕn)(ϕnhPhϕn)=ηnϕθnϕ (20)

    Lemma 4.1. For Algorithm 1, The following inequality holds

    ||θn+1ˆu||2||θn+1u||2+||θn+1ˆuθn+1u||2+β(||θn+1u||2||θnu||2+12||(θn+1uθnu)||2)+γΔt||θn+1u||2βΔt||θnu||2dβ(1+2Δt)||ηu,t||2L2(tn,tn+1;L2(Ωf))dγΔt||ηn+1u||2

    Proof. For the detailed proof process, please refer to Lemma 10 in [3].

    Theorem 4.2. Under the regularity assumption (18). Suppose β>0, then there exists a constant C>0, such that

    ||eNu||2+β||eNu||2+||eNϕ||2+N1n=0(||en+1ˆuen+1u||2+||en+1ˆuenu||2+β2||(en+1ˆuenu)||2+||en+1ϕenϕ||2)+N1n=0γΔt||en+1u||2+N1n=0νΔt||en+1ˆu||2+N1n=0gλminΔt||en+1ϕ||2C(h2k+Δt2+Δth2k)

    Proof. The true solution satisfies the following relations:

    (un+1unΔt,vh)+af(un+1,vh)+b(vh,pn+1)+aΓ(vh,ϕn+1)+af,c(un,un+1,vh)=(fn+11,vh)+(un+1unΔtun+1t,vh)af,c(un+1un,un+1,vh)b(un+1,qh)=0
    gs0(ϕn+1ϕnΔt,ψh)+ap(ϕn+1,ψh)aΓ(un+1,ψh)=(fn+12,ψh)+gs0(ϕn+1ϕnΔtϕn+1t,ψh). (21)

    Subtracting (16) from (21), we arrive that

    (en+1ˆuenuΔt,vh)+af(en+1ˆu,vh)+b(vh,en+1p)+aΓ(vh,ϕn+1ϕnh)+af,c(un,un+1,vh)af,c(unh,ˆun+1h,vh)=(un+1unΔtun+1t,vh)af,c(un+1un,un+1,vh)b(en+1ˆu,qh)=0gso(en+1ϕenϕΔt,ψh)+ap(en+1ϕ,ψh)aΓ(un+1unh,ψh)=gs0(ϕn+1ϕnΔtϕn+1t,ψh) (22)

    Setting vh=2Δtθn+1ˆu,qh=2Δtθn+1p, and ψh=2Δtθn+1ϕ in (22), the resulting equations are added, this yields:

    ||θn+1ˆu||2||θnu||2+||θn+1ϕ||2||θnϕ||2+||θn+1ˆuθnu||2+||θn+1ϕθnϕ||2+2Δtaf(θn+1ˆu,θn+1ˆu)+2Δtap(θn+1ϕ,θn+1ϕ)=2(ηn+1uηnu,θn+1ˆu)+2gs0(ηn+1ϕηnϕ,θn+1ϕ)+2Δtaf(ηn+1u,θn+1ˆu)+2Δtap(ηn+1ϕ,θn+1ϕ)+2Δtaf,c(un,un+1,θn+1ˆu)2Δtaf,c(unh,ˆun+1h,θn+1ˆu)+2ΔtaΓ(θn+1ˆu,ϕn+1ϕnh)2ΔtaΓ(un+1unh,θn+1ϕ)2Δt(un+1unΔtun+1t,θn+1ˆu)+2Δtaf,c(un+1un,un+1,θn+1ˆu)2Δtgs0(ϕn+1ϕnΔtϕn+1t,θn+1ϕ) (23)

    Next, we bound each term on the right hand side of (23), by virtue of Cauchy-Schwarz-Young inequality,

    2(ηn+1uηnu,θn+1ˆu)10c2pε1||ηu,t||2L2(tn,tn+1;L2(Ωf))+ε1Δt10||θn+1ˆu||22gs0(ηn+1ϕηnϕ,θn+1ϕ)6gs20~c2pλmin||ηϕ,t||2L2(tn,tn+1;L2(Ωp))+gλnimΔt6||θn+1ϕ||22Δtaf(ηn+1u,θn+1ˆu)2C3Δt||ηn+1u||||θn+1ˆu||10C23Δtε2||ηn+1u||2+ε2Δt10||θn+1ˆu||22Δtap(ηn+1ϕ,θn+1ϕ)2gλmax||ηn+1ϕ||||θn+1ϕ||6gλ2maxΔtλmin||ηn+1ϕ||2+gλminΔt6||θn+1ϕ||22Δt(un+1unΔtun+1t,θn+1ˆu)2Δtaf,c(un+1un,un+1,θn+1ˆu)CΔt2ε3(||utt||2L2(tn,tn+1;L2(Ωf))+||ut||2L2(tn,tn+1;L2(Ωf)))+ε3Δt10||θn+1ˆu||2
    2Δtgs0(ϕn+1ϕnΔtϕn+1t,θn+1ϕ)6gs02~c2pΔt2λmin||ϕtt||2L2(tn,tn+1;L2(Ωp))+gλminΔt6||θn+1ϕ||2 (24)

    For the interface terms, we treat them as follows:

    2ΔtaΓ(θn+1ˆu,ϕn+1ϕnh)=2ΔtaΓ(θn+1ˆu,ηnϕ)2ΔtaΓ(θn+1ˆu,θnϕ)+2ΔtaΓ(θn+1ˆu,ϕn+1ϕn)2ΔtaΓ(un+1unh,θn+1ϕ)=2ΔtaΓ(ηnu,θn+1ϕ)2ΔtaΓ(θnu,θn+1ϕ)+2ΔtaΓ(un+1un,θn+1ϕ) (25)

    So, we have the following estimates:

    2ΔtaΓ(θn+1ˆu,ηnϕ)10C24Δtε4||ηnϕ||2+ε4Δt10||θn+1ˆu||22ΔtaΓ(θn+1ˆu,θnϕ)2Δt(C6g2h1||θnϕ||2+||θn+1ˆu||2)10C6g2Δthε5||θnϕ||2+ε5Δt10||θn+1ˆu||22ΔtaΓ(θn+1ˆu,ϕn+1ϕn)2Δt(C6g2h1||ϕn+1ϕn||2+||θn+1ˆu||2)10C6g2Δt2hε6||ϕt||2L2(tn,tn+1;L2(Ωp))+ε6Δt10||θn+1ˆu||22ΔtaΓ(ηnu,θn+1ϕ)6C24Δtgλmin||ηnu||2+gλminΔt6||θn+1ϕ||22ΔtaΓ(θnu,θn+1ϕ)2Δ(C5g2h1||θnu||2+||θn+1ϕ||2)6C5gΔthλmin||θnu||2+gλminΔt6||θn+1ϕ||22ΔtaΓ(un+1un,θn+1ϕ)2C4Δt||(un+1un)||||θn+1ϕ||6C24Δt2gλmin||ut||2L2(tn,tn+1;L2(Ωf))+gλminΔt6||θn+1ϕ||2 (26)

    For the trilinear form, we have

    2Δtaf,c(un,un+1,θn+1ˆu)2Δtaf,c(unh,ˆun+1h,θn+1ˆu)=2Δtaf,c(ηnu;un+1,θn+1ˆu)2Δtaf,c(θnu;un+1,θn+1ˆu)+2Δtaf,c(ˆun+1h;ηn+1u,θn+1ˆu)+2Δtaf,c(unhˆun+1h,ηn+1u,θn+1ˆu), (27)

    then

    2Δtaf,c(ηnu;un+1,θn+1ˆu)2C1Δt||ηnu||||un+1||||θn+1ˆu||2C1Δt(||ηnu||2||un+1||22+||θn+1ˆu||22)10C21Δtε7||ηnu||2||un+1||2+ε7Δt10||θn+1ˆu||22Δtaf,c(θnu;un+1,θn+1ˆu)2C2Δt||θnu||||un+1||2||θn+1ˆu||
    10C22Δtε8||un+1||22||θnu||2+ε8Δt10||θn+1ˆu||22Δtaf,c(ˆun+1h;ηn+1u,θn+1ˆu)2C1Δt||ˆun+1h||||ηn+1u||||θn+1ˆu||10C21Δtε9||ˆun+1h||2||ηn+1u||2+ε9Δt10||θn+1ˆu||22Δtaf,c(unhˆun+1h,ηn+1u,θn+1ˆu)2C1Δt||(unhˆun+1h)||||ηn+1u||||θn+1ˆu||10C21Δtε10||(unhˆun+1h)||2||ηn+1u||2+ε10Δt10||θn+1ˆu||220C21Δtε10(||unh||2+||ˆun+1h||2)||ηn+1u||2+ε10Δt10||θn+1ˆu||2 (28)

    Combining Lemma 4.1 with the properties (13), we obtain

    ||θn+1u||2||θnu||2+||θn+1ϕ||2||θnϕ||2+β||θn+1u||2β||θnu||2+||θn+1ˆuθnu||2+||θn+1ˆuθn+1u||2+||θn+1ϕθnϕ||2+β2||(θn+1uθnu)||2+γΔt||θn+1u||2+2νΔt||θn+1ˆu||2+2gλminΔt||θn+1ϕ||22(ηn+1uηnu,θn+1ˆu)+2gs0(ηn+1ϕηnϕ,θn+1ϕ)+2Δtaf(ηn+1u,θn+1ˆu)+2Δtap(ηn+1ϕ,θn+1ϕ)+2Δtaf,c(ηnu;un+1,θn+1ˆu)2Δtaf,c(θnu;un+1,θn+1ˆu)+2Δtaf,c(ˆun+1h;ηn+1u,θn+1ˆu)+2Δtaf,c(unhˆun+1h,ηn+1u,θn+1ˆu)+2ΔtaΓ(θn+1ˆu,ηnϕ)2ΔtaΓ(θn+1ˆu,θnϕ)+2ΔtaΓ(θn+1ˆu,ϕn+1ϕn)2ΔtaΓ(ηnu,θn+1ϕ)+2ΔtaΓ(θnu,θn+1ϕ)2ΔtaΓ(un+1un,θn+1ϕ)2Δt(un+1unΔtun+1t,θn+1ˆu)+2Δtaf,c(un+1un,un+1,θn+1ˆu)2Δtgs0(ϕn+1ϕnΔtϕn+1t,θn+1ϕ)+βΔt||θnu||2+dβ(1+2Δt)||ηu,t||2L2(tn,tn+1;L2(Ωf))+dγΔt||ηn+1u||2. (29)

    Inserting the above results into (29) and let ε1=ε2==ε10=ν, we can get the following estimates:

    ||θn+1u||2||θnu||2+β||θn+1u||2β||θnu||2+||θn+1ϕ||2||θnϕ||2+||θn+1ˆuθnu||2+||θn+1ˆuθn+1u||2+||θn+1ϕθnϕ||2+β2||(θn+1uθnu)||2+γΔt||θn+1u||2+νΔt||θn+1ˆu||2+gλminΔt||θn+1ϕ||210c2pν||ηu,t||2L2(tn,tn+1;L2(Ωf))+(10C23Δtν+dγΔt)||ηn+1u||2+6gS20~c2pλmin||ηϕ,t||2L2(tn,tn+1;L2(Ωp))+6gλ2maxΔtλmin||ηn+1ϕ||2+6C24Δtgλmin||ηnu||2++10C24Δtν||ηnϕ||2+(6C5gΔthλmin+10C22Δtν||un+1||22)||θnu||2+Δtβ||θnu||2+10C6g2Δthν||θnϕ||2+(10C21Δtν||ˆun+1h||2+20C21Δtν(||unh||2+||ˆun+1h||2))||ηn+1u||2
    +10C21Δtν||un+1||2||ηnu||2+CΔt2ν||utt||2L2(tn,tn+1;L2(Ωf))+(CΔt2ν+6C24Δt2gλmin)||ut||2L2(tn,tn+1;L2(Ωf))+6gS20~c2pΔt2λmin||ϕtt||2L2(tn,tn+1;L2(Ωp))+10C6g2Δt2hν||ϕt||2L2(tn,tn+1;L2(Ωp))+dβ(1+2Δt)||ηu,t||2L2(tn,tn+1;L2(Ωf)) (30)

    Denote C=exp(ΔtN1n=0max(6C5ghλmin,10C22ν||u||2,2,β,10C6g2hν)). Sum (30) from n=0 to n=N1, together with Theorem 3.2 and Lemma 2.2, we get the following inequality:

    ||θNu||2+β||θNu||2+||θNϕ||2+N1n=0(||θn+1ˆuθnu||2+||θn+1ˆuθn+1u||2+β2||||θn+1ˆuθnu||2+||θn+1ϕθnϕ||2)+γΔtN1n=0||θn+1u||2+νΔtN1n=0||θn+1ˆu||2+gλminΔtN1n=0||θn+1ϕ||2C[10c2pν||ηu,t||2L2(0,T;L2(Ωf))+(10C23ν+6C24gλmin+dγ+1ν)||ηu||2+6gS20~c2pλmin||ηϕ,t||2L2(0,T;L2(Ωf))+(6λ2maxgλmin+10C24ν)||ηnϕ||2+CΔt2ν||utt||2L2(0,T;L2(Ωf))+(CΔt2ν+6C24Δt2gλmin)||ut||2L2(0,T;L2(Ωf))+6gS20~c2pΔt2λmin||ϕtt||2L2(0,T;L2(Ωp))+10C6g2Δt2hν||ϕt||2L2(0,T;L2(Ωp))+dβ(1+2Δt)||ηu,t||2L2(0,T;L2(Ωf))+||θ0u||2+β||θ0u||2+||θ0ϕ||2] (31)

    Finally, we have

    ||eNu||2+β||eNu||2+||eNϕ||2+N1n=0(||en+1ˆuen+1u||2+||en+1ˆuenu||2+β2||(en+1ˆuenu)||2+||en+1ϕenϕ||2)+N1n=0γΔt||en+1u||2+N1n=0νΔt||en+1ˆu||2+N1n=0gλminΔt||en+1ϕ||2C(h2k+Δt2+Δth2k) (32)

    where C>0 is a constant.

    Remark 1. Here we can propose a second-order backward differentiation formula(BDF2) methods for Stokes/Darcy model and Navier-Stokes/Darcy model, respectively. We will analyze it in the future.

    Algorithm 2 (The BDF2 modular grad-div scheme for Stokes/Darcy model).

    For vhHfh, qhQh, and ψhHph.

    Step1. Given un1h,unhHfh, pnhQh, and ϕn1h,ϕnhHph, find ˆun+1h,pn+1h,ϕn+1h satisfying:

    (3ˆun+1h4unh+un1h2Δt,vh)+af(ˆun+1h,vh)+b(vh,pn+1h)+aΓ(vh,2ϕnhϕn1h)=(fn+11,vh),b(ˆun+1h,qh)=0,gS0(3ϕn+1h4ϕnh+ϕn1h2Δt,ψh)+ap(ϕn+1h,ψh)aΓ(2unhun1h,ψh)=(fn+12,ψh). (33)

    Step2. Given ˆun+1hHfh, find un+1hHfh satisfying:

    (3un+1h3ˆun+1h2Δt,vh)+β(3un+1h4unh+un+1h2Δt,vh)+γ(un+1h,vh)=0. (34)

    Algorithm 3 (The BDF2 modular grad-div scheme for Navier-Stokes/Darcy model).

    For vhHfh, qhQh, and ψhHph.

    Step1. Given un1h,unhHfh, pnhQh, and ϕn1h,ϕnhHph, find ˆun+1h,pn+1h,ϕn+1h satisfying:

    (3ˆun+1h4unh+un1h2Δt,vh)+af(ˆun+1h,vh)+b(vh,pn+1h)+af,c(2unhun1h,ˆun+1h,vh)aΓ(vh,2ϕnϕn1)=(fn+11,vh)b(ˆun+1h,qh)=0gs0(3ϕn+1h4ϕnh+ϕn1h2Δt,ψh)+ap(ϕn+1h,ψh)aΓ(2unun1,ψh)=(fn+12,ψh) (35)

    Step2. Given ˆun+1hHfh, find un+1hHfh satisfying:

    (3un+1h3ˆun+1h2Δt,vh)+β(3un+1h4unh+un1h2Δt,vh)+γ(un+1h,vh)=0 (36)

    In this section, We will compare two grad-div schemes with standard scheme respectively to justify the results of the theoretical analysis. We implement numerical experiments using software Freefem++.

    The domain Ω be decomposed into Ωf=(0,1)×(1,2) and Ωp=(0,1)×(0,1) with the interface Γ=(0,1)×{1}. The exact solution is taken as follows:

    (u1,u2)=([x2(y1)2+y]cos(t), [23x(y1)3]cos(t)+[2πsin(πx)]cos(t)),p=[2πsin(πx)]sin(0.5πy)cos(t),ϕ=[2πsin(πx)][1ycos(πy)]cos(t).

    Here, the parameters n,ρ,g,ν,K,S0 and α are set to 1. The initial conditions, boundary conditions, and the forcing terms follows the exact solution, and we set h=Δt in experiments. The famous Taylor-Hood element(P2-P1) for the Navier-Stokes problem and the piecewise quadratic polynomials(P2) for the Darcy flow are used, respectively. The numerical results are presented in Table 1-Table 4.

    Table 1.  Numerical results at time T = 1 for the standard scheme.
    1h||eu||L2uL2rate||eu||fuHfrate||eu||L2divuL2rate
    40.01589060.03528820.042823
    80.008477820.9064080.01618721.124330.009784382.12983
    160.004367990.9567240.008050771.007650.002301982.08761
    320.002215160.9795590.004043690.9934540.0006099021.91623
    640.001115310.9899660.002031120.9933970.0001314012.2146
    1h||eϕ||L2ϕL2rate||eϕ||pϕHprate||ep||L2pL2rate
    40.04047640.06530310.500655
    80.01824771.149370.01826491.838080.257160.961151
    160.009273250.9765680.007774671.232220.1305690.977854
    320.004686120.9846810.003706711.068640.06583430.987901
    640.002355840.9921520.001838881.011310.03304730.994307

     | Show Table
    DownLoad: CSV
    Table 2.  Numerical results at time T = 1 for the standard grad-div scheme.
    1h||eu||L2uL2rate||eu||fuHfrate||eu||L2divuL2rate
    40.0157960.03499070.0332181
    80.008474770.8983130.01616291.114290.007956382.06179
    160.004367880.9562410.008048561.005880.001925462.04691
    320.002215150.9795290.004043510.9931230.0005518831.80277
    640.001115310.989960.00203110.9933470.0001129742.28837
    1h||eϕ||L2ϕL2rate||eϕ||pϕHprate||ep||L2pL2rate
    40.0403890.06527550.512196
    80.0182391.146940.01826131.837750.2574430.992443
    160.009272650.9759730.007774411.231980.1305820.979297
    320.004686070.9846030.003706691.06860.06583570.988014
    640.002355830.9921430.001838881.01130.03304740.994333

     | Show Table
    DownLoad: CSV
    Table 3.  Numerical results at time T = 1 for the modular grad-div scheme.
    1h||eu||L2uL2rate||eu||fuHfrate||eu||L2divuL2rate
    40.01612270.04222950.00546336
    80.008494780.9244450.01869871.175310.001336662.03116
    160.004368880.9593130.008496081.138070.0002500082.41859
    320.002215260.9797870.00423941.002947.36766e-0051.7627
    640.001115310.9900310.002044151.052366.93936e-0063.40833
    1h||eϕ||L2ϕL2rate||eϕ||pϕHprate||ep||L2pL2rate
    40.04038030.06525230.500113
    80.01823351.147060.01825841.837470.2571520.959633
    160.009272490.9755630.007774361.231760.130570.977798
    320.004686080.9845750.00370671.068590.06583450.987908
    640.002355830.9921460.001838881.011310.03304730.994311

     | Show Table
    DownLoad: CSV
    Table 4.  The ||eu||f for the standard without grad-div scheme, standard scheme and grad-div scheme with vaying hydraulic conductivity tensor K.
    KNon-stabilizedStandard grad-divmodular grad-div
    I0.01691730.01080850.077557
    1e1I0.02029740.01304370.0775562
    1e2I0.04646250.03018790.077554
    1e3I0.1242380.08249880.0775521

     | Show Table
    DownLoad: CSV

    Table 1, Table 2 and Table 3 shows the error and convergence order of velocity u, pressure p and hydraulic head ϕ by using the standard scheme, standard grad-div scheme and modular grad-div scheme, respectively. where γ=1 and β=0.2 for standard grad-div scheme and modular grad-div scheme. By observation and comparison, it can be found that the divergence velocity errors of the standard grad-div and modular grad-div scheme are smaller than that of the standard scheme, and the modular grad-div scheme is more accurate. Moreover both numerical results is consist with the theoretical analysis.

    Next, we set up a numerical experiment by using the different parameter K to show the superiority of modular grad-div scheme. Let's fix Δt=h=140. Table 4 shows the ||eu||f for the standard scheme without stabilization, standard grad-div scheme and modular grad-div scheme with hydraulic conductity K=I, 0.1I, 0.01I, 0.001I. We observe that the divergence of velocity error for all three schemes increase with the decrease of K. Yet their growth is relatively small, in particular there's almost no change in modular grad-div scheme. Such results are also consistent with our theoretical analysis.

    In this paper, we extend the grad-div stabilized method from Stokes/Darcy model to Navier-Stokes/Darcy model. Stability and error estimates are provided. Numerical experiments confirm the theoretical analysis, and show that the modular grad-div scheme is more efficient than that of the standard grad-div scheme.



    [1] Numerical solution to a mixed Navier-Stokes/Darcy model by the two-grid approach. SIMA J. Numer. Anal. (2009) 47: 3325-3338.
    [2] Mortar element method for the time dependent coupling of Stokes and Darcy flows. J. Sci. Comput. (2019) 80: 1310-1329.
    [3] An efficient and modular grad-div stabilization. Comput. Methods Appl. Mech. Engrg. (2018) 355: 327-346.
    [4] Two classes of mixed finite element methods. Comput. Methods Appl. Mech. Engrg. (1988) 69: 89-129.
    [5] DG approximation of coupled Navier-Stokes and Darcy equations by Beaver-Joseph-Saffman interface condition. SIAM J. Numer. Anal. (2009) 47: 2052-2089.
    [6] A modified two-grid decoupling method for the mixed Navier-Stokes/Darcy model. Comput. Math.Appl. (2016) 72: 1142-1152.
    [7] Partitioned time stepping method for fully evolutionary Navier-Stokes/ Darcy flow with BJS interface conditions. Adv. Appl. Math. Mech. (2019) 11: 381-405.
    [8] Decoupled characteristic stabilized finite element method for time-dependent Navier-Stokes/Darcy model. Numer. Meth. Part D. E. (2019) 35: 267-294.
    [9] A decoupling method with different subdomain time steps for the non-stationary Navier-Stokes/Darcy model. J. Comput. Math. (2017) 35: 319-345.
    [10] On the convergence rate of grad-giv stabilized Taylor-Hood to Scott-Vogelius solutions of incompressible flow problems. J. Math. Anal. Appl. (2011) 381: 612-626.
    [11] X. Lu and P. Huang, A modular grad-div stabilization for the 2D/3D nonstationary incompressible magnetohydrodynamic equations, J. Sci. Comput., 82 (2020), Paper No. 3, 24 pp. doi: 10.1007/s10915-019-01114-x
    [12] Decoupled schemes for a non-stationary mixed Stokes-Darcy model. Math. Comput. (2010) 79: 707-731.
    [13] Numerical analysis of two grad-div stabilization methods for the time-dependent Stokes/Darcy model. Comput. Math. Appl. (2020) 79: 817-832.
    [14] Y. Zeng and P. Huang, A grad-div stabilized projection finite element method for a double-diffusive natural convection model, Numer. Heat Tr. B-Fund., (2020). doi: 10.1080/10407790.2020.1747285
    [15] Two-grid finite element methods for the steady Navier-Stokes/Darcy model. East Asian J. Appl. Math. (2016) 6: 60-79.
    [16] A decoupling two-grid algorithm for the mixed Stokes-Darcy model with the Beavers-Joseph interface condition. Numer. Methods Partial Differential Equations (2014) 30: 1066-1082.
  • Reader Comments
  • © 2020 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(2580) PDF downloads(210) Cited by(0)

Figures and Tables

Figures(1)  /  Tables(4)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog