
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
[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 |
[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
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
The time-dependent Navier-Stokes equation govern the fluid flow in
∂u∂t−νΔ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
The Darcy equation govern the porous media flow in
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.
We know that
ϕ=z+ppρg |
is the piezometric head, where
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
u⋅nf+up⋅np=0on Γ×(0,T) | (4) |
p−νnf∂u∂nf=gϕon Γ×(0,T), | (5) |
−ντi∂u∂nf=α√gν√τi⋅Kτiu⋅τi,i=1,⋅⋅⋅,d−1on Γ×(0,T), | (6) |
where,
For simplicity of analysis, we impose the following boundary conditions on
u=0on Γf×(0,T)ϕ=0on Γp×(0,T) | (7) |
In this paper, we assume that:
u⋅nf>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
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
||u||Hf=||∇u||L2(Ωf)=√(∇u,∇u)Ωf ∀u∈Hf,||ϕ||Hp=||∇ϕ||L2(Ωp)=√(∇ϕ,∇ϕ)Ωf ∀ϕ∈Hp. |
Some discrete norms are defined as,
||w||2L2(0,T;Hs(Ωf,p))=ΔtN∑n=0||wn||2Hs(Ωf,p),||w||L∞(0,T;Hs(Ωf,p))=max0≤n≤N||wn||Hs(Ωf,p). |
Due to
af,c(u;v,w)=((u⋅∇)v,w)f+12(∇⋅u,v⋅w)f=12((u⋅∇)v,w)f−12((u⋅∇)w,v)f+12⟨v⋅w,u⋅nf⟩Γ, | (9) |
under the condition (8), we have
af,c(u;v,v)≥0. | (10) |
Then, we have the following estimates for
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
||v||L2≤cp||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
(ut,v)f+af(u,v)+b(v,p)+aΓ(v,ϕ)+af,c(u;u,v)=(f1,v)f∀v∈Hf,b(u,q)=0∀q∈Q,gS0(ϕt,ψ)p+ap(ϕ,ψ)−aΓ(u,ψ)=g(f2,ψ)p∀ψ∈Hp. | (12) |
Where
af(u,v)=ν(D(u),D(v))+d−1∑i=1∫Γα√νgτi⋅Kτi(u⋅τi)(v⋅τi)ds,ap(ϕ,ψ)=g(K∇ϕ,∇ψ)p,aΓ(v,ϕ)=g∫Γϕv⋅nfds,af,c(u,u,v)=((u⋅∇)u,v)f,b(v,p)=−(p,∇⋅v)f. |
Bilinear forms
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Γ(u,ϕ)|≤C4||∇u||f||∇ϕ||p|aΓ(u,ϕ)|≤C5g2h−1||u||2f+||∇ϕ||2p,|aΓ(u,ϕ)|≤C6g2h−1||ϕ||2p+||∇u||2f | (14) |
Lemma 2.1. We assume that
f1∈L2(0,T;L2(Ωf)d),f2∈L2(0,T;L2(Ωp)d),K∈L∞(Ωp)d×d, |
and
0<kmin|x|2≤Kx⋅x≤kmax|x|2∀ x∈Ωp. |
Furthermore,
Lemma 2.2. (Discrete Gronwall Lemma). Let
aN+ΔtN∑n=0bn≤ΔtN−1∑n=0dnan+ΔtN∑n=0cn+H |
then for all
aN+ΔtN∑n=0bn≤exp(ΔtN−1∑n=0dn)(ΔtN∑n=0cn+H). |
We construct
infqh∈Qhsupvh∈Hfh(qh,∇⋅vfh)f||qh||Q||vfh||Hf≥β. | (15) |
Next, we divide time interval [0, T] into:
Algorithm 1 (The modular grad-div scheme).
For
Step1. Given
(ˆun+1h−unhΔ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+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+1h−un+1h||+2γΔt||∇⋅un+1h||2+β(||∇⋅un+1h||2−||∇⋅unh||2+||∇⋅(un+1h−unh)||2). |
Proof. Refer to Lemma 6 of [3] for proof details.
Theorem 3.2. (Unconditional Stability) For any
||uNh||2+β||∇uNh||2+gs0||ϕNh||2+N−1∑n=0(||ˆun+1h−unh||2+||ˆun+1h−un+1h||2+β||∇⋅(un+1h−unh)||2+gs0||ϕn+1h−ϕnh||2)+N−1∑n=02γΔt||∇⋅un+1h||2+N−1∑n=0νΔt||∇ˆun+1h||2+N−1∑n=0λmingΔ||∇ϕn+1h||2≤C(ΔtN−1∑n=0(2c2pν||fn+11||2+2c2pλming||fn+12||2)+||u0h||2+β||∇⋅u0h||2+gs0||ϕ0h||2), |
Where the constant
Proof. Taking
In this section, we will give some error estimates of our proposed method. Denote
u∈L∞(0,T;Hf∩Hk+1(Ωf)d),ut∈L∞(0,T;Hk+1(Ωf)d),utt∈L2(0,T;L2(Ωf)d).p∈L2(0,T;Q∩Hk(Ωf)).ϕ∈L∞(0,T;Hp∩Hk+1(Ωp)),ϕt∈L∞(0,T;Hk+1(Ωp)),ϕtt∈L2(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
||Phu(t)−u(t)||f≤Chk+1||u(t)||Hk+1(Ωf),||∇(Phu(t)−u(t))||f≤Chk||u(t)||Hk+1(Ωf),||Php(t)−p(t)||f≤Chk+1||p(t)||Hk+1(Ωf),||Phϕ(t)−ϕ(t)||p≤Chk+1||ϕ(t)||Hk+1(Ωp),||∇(Phϕ(t)−ϕ(t))||p≤Chk||ϕ(t)||Hk+1(Ωp) | (19) |
Next, we define the following error equations
enu=un−unh=(un−Phun)−(unh−Phun)=ηnu−θnu,enˆu=un−ˆunh=(un−Phun)−(ˆunh−Phun)=ηnu−θnˆuenp=pn−pnh=(pn−Phpn)−(pnh−Phpn)=ηnp−θnpenϕ=ϕn−ϕnh=(ϕn−Phϕn)−(ϕnh−Phϕ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||2−dβ(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
||eNu||2+β||∇⋅eNu||2+||eNϕ||2+N−1∑n=0(||en+1ˆu−en+1u||2+||en+1ˆu−enu||2+β2||∇⋅(en+1ˆu−enu)||2+||en+1ϕ−enϕ||2)+N−1∑n=0γΔt||∇⋅en+1u||2+N−1∑n=0νΔt||∇en+1ˆu||2+N−1∑n=0gλminΔt||∇en+1ϕ||2≤C(h2k+Δt2+Δth2k) |
Proof. The true solution satisfies the following relations:
(un+1−unΔt,vh)+af(un+1,vh)+b(vh,pn+1)+aΓ(vh,ϕn+1)+af,c(un,un+1,vh)=(fn+11,vh)+(un+1−unΔt−un+1t,vh)−af,c(un+1−un,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ˆu−enuΔ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+1−unΔt−un+1t,vh)−af,c(un+1−un,un+1,vh)b(en+1ˆu,qh)=0gso(en+1ϕ−enϕΔt,ψh)+ap(en+1ϕ,ψh)−aΓ(un+1−unh,ψh)=gs0(ϕn+1−ϕnΔt−ϕn+1t,ψh) | (22) |
Setting
||θ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+1−unh,θn+1ϕ)−2Δt(un+1−unΔt−un+1t,θn+1ˆu)+2Δtaf,c(un+1−un,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+1−unΔt−un+1t,θn+1ˆu)−2Δtaf,c(un+1−un,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+1−unh,θn+1ϕ)=2ΔtaΓ(ηnu,θn+1ϕ)−2ΔtaΓ(θnu,θn+1ϕ)+2ΔtaΓ(un+1−un,θ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(C6g2h−1||θ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(C6g2h−1||ϕ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Δ(C5g2h−1||θnu||2+||∇θn+1ϕ||2)≤6C5gΔthλmin||θnu||2+gλminΔt6||∇θn+1ϕ||22ΔtaΓ(un+1−un,θn+1ϕ)≤2C4Δt||∇(un+1−un)||||∇θ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||2≤20C21Δ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ϕ||2≤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(η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+1−un,θn+1ϕ)−2Δt(un+1−unΔt−un+1t,θn+1ˆu)+2Δtaf,c(un+1−un,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
||θ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ϕ||2≤10c2pν||η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
||θNu||2+β||∇⋅θNu||2+||θNϕ||2+N−1∑n=0(||θn+1ˆu−θnu||2+||θn+1ˆu−θn+1u||2+β2||∇⋅||θn+1ˆu−θnu||2+||θn+1ϕ−θnϕ||2)+γΔtN−1∑n=0||∇⋅θn+1u||2+νΔtN−1∑n=0||∇θn+1ˆu||2+gλminΔtN−1∑n=0||∇θn+1ϕ||2≤C∗[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+N−1∑n=0(||en+1ˆu−en+1u||2+||en+1ˆu−enu||2+β2||∇⋅(en+1ˆu−enu)||2+||en+1ϕ−enϕ||2)+N−1∑n=0γΔt||∇⋅en+1u||2+N−1∑n=0νΔt||∇en+1ˆu||2+N−1∑n=0gλminΔt||∇en+1ϕ||2≤C(h2k+Δt2+Δth2k) | (32) |
where
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
Step1. Given
(3ˆun+1h−4unh+un−1h2Δt,vh)+af(ˆun+1h,vh)+b(vh,pn+1h)+aΓ(vh,2ϕnh−ϕn−1h)=(fn+11,vh),b(ˆun+1h,qh)=0,gS0(3ϕn+1h−4ϕnh+ϕn−1h2Δt,ψh)+ap(ϕn+1h,ψh)−aΓ(2unh−un−1h,ψh)=(fn+12,ψh). | (33) |
Step2. Given
(3un+1h−3ˆun+1h2Δt,vh)+β(∇⋅3un+1h−4unh+un+1h2Δt,∇⋅vh)+γ(∇⋅un+1h,∇⋅vh)=0. | (34) |
Algorithm 3 (The BDF2 modular grad-div scheme for Navier-Stokes/Darcy model).
For
Step1. Given
(3ˆun+1h−4unh+un−1h2Δt,vh)+af(ˆun+1h,vh)+b(vh,pn+1h)+af,c(2unh−un−1h,ˆun+1h,vh)aΓ(vh,2ϕn−ϕn−1)=(fn+11,vh)b(ˆun+1h,qh)=0gs0(3ϕn+1h−4ϕnh+ϕn−1h2Δt,ψh)+ap(ϕn+1h,ψh)−aΓ(2un−un−1,ψh)=(fn+12,ψh) | (35) |
Step2. Given
(3un+1h−3ˆun+1h2Δt,vh)+β(∇⋅3un+1h−4unh+un−1h2Δ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
(u1,u2)=([x2(y−1)2+y]cos(t), [−23x(y−1)3]cos(t)+[2−πsin(πx)]cos(t)),p=[2−πsin(πx)]sin(0.5πy)cos(t),ϕ=[2−πsin(πx)][1−y−cos(πy)]cos(t). |
Here, the parameters
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.0158906 | 0.0352882 | 0.042823 | |||
8 | 0.00847782 | 0.906408 | 0.0161872 | 1.12433 | 0.00978438 | 2.12983 |
16 | 0.00436799 | 0.956724 | 0.00805077 | 1.00765 | 0.00230198 | 2.08761 |
32 | 0.00221516 | 0.979559 | 0.00404369 | 0.993454 | 0.000609902 | 1.91623 |
64 | 0.00111531 | 0.989966 | 0.00203112 | 0.993397 | 0.000131401 | 2.2146 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.0404764 | 0.0653031 | 0.500655 | |||
8 | 0.0182477 | 1.14937 | 0.0182649 | 1.83808 | 0.25716 | 0.961151 |
16 | 0.00927325 | 0.976568 | 0.00777467 | 1.23222 | 0.130569 | 0.977854 |
32 | 0.00468612 | 0.984681 | 0.00370671 | 1.06864 | 0.0658343 | 0.987901 |
64 | 0.00235584 | 0.992152 | 0.00183888 | 1.01131 | 0.0330473 | 0.994307 |
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.015796 | 0.0349907 | 0.0332181 | |||
8 | 0.00847477 | 0.898313 | 0.0161629 | 1.11429 | 0.00795638 | 2.06179 |
16 | 0.00436788 | 0.956241 | 0.00804856 | 1.00588 | 0.00192546 | 2.04691 |
32 | 0.00221515 | 0.979529 | 0.00404351 | 0.993123 | 0.000551883 | 1.80277 |
64 | 0.00111531 | 0.98996 | 0.0020311 | 0.993347 | 0.000112974 | 2.28837 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.040389 | 0.0652755 | 0.512196 | |||
8 | 0.018239 | 1.14694 | 0.0182613 | 1.83775 | 0.257443 | 0.992443 |
16 | 0.00927265 | 0.975973 | 0.00777441 | 1.23198 | 0.130582 | 0.979297 |
32 | 0.00468607 | 0.984603 | 0.00370669 | 1.0686 | 0.0658357 | 0.988014 |
64 | 0.00235583 | 0.992143 | 0.00183888 | 1.0113 | 0.0330474 | 0.994333 |
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.0161227 | 0.0422295 | 0.00546336 | |||
8 | 0.00849478 | 0.924445 | 0.0186987 | 1.17531 | 0.00133666 | 2.03116 |
16 | 0.00436888 | 0.959313 | 0.00849608 | 1.13807 | 0.000250008 | 2.41859 |
32 | 0.00221526 | 0.979787 | 0.0042394 | 1.00294 | 7.36766e-005 | 1.7627 |
64 | 0.00111531 | 0.990031 | 0.00204415 | 1.05236 | 6.93936e-006 | 3.40833 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.0403803 | 0.0652523 | 0.500113 | |||
8 | 0.0182335 | 1.14706 | 0.0182584 | 1.83747 | 0.257152 | 0.959633 |
16 | 0.00927249 | 0.975563 | 0.00777436 | 1.23176 | 0.13057 | 0.977798 |
32 | 0.00468608 | 0.984575 | 0.0037067 | 1.06859 | 0.0658345 | 0.987908 |
64 | 0.00235583 | 0.992146 | 0.00183888 | 1.01131 | 0.0330473 | 0.994311 |
K | Non-stabilized | Standard grad-div | modular grad-div |
I | 0.0169173 | 0.0108085 | 0.077557 |
1e−1I | 0.0202974 | 0.0130437 | 0.0775562 |
1e−2I | 0.0464625 | 0.0301879 | 0.077554 |
1e−3I | 0.124238 | 0.0824988 | 0.0775521 |
Table 1, Table 2 and Table 3 shows the error and convergence order of velocity
Next, we set up a numerical experiment by using the different parameter
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. ![]() |
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.0158906 | 0.0352882 | 0.042823 | |||
8 | 0.00847782 | 0.906408 | 0.0161872 | 1.12433 | 0.00978438 | 2.12983 |
16 | 0.00436799 | 0.956724 | 0.00805077 | 1.00765 | 0.00230198 | 2.08761 |
32 | 0.00221516 | 0.979559 | 0.00404369 | 0.993454 | 0.000609902 | 1.91623 |
64 | 0.00111531 | 0.989966 | 0.00203112 | 0.993397 | 0.000131401 | 2.2146 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.0404764 | 0.0653031 | 0.500655 | |||
8 | 0.0182477 | 1.14937 | 0.0182649 | 1.83808 | 0.25716 | 0.961151 |
16 | 0.00927325 | 0.976568 | 0.00777467 | 1.23222 | 0.130569 | 0.977854 |
32 | 0.00468612 | 0.984681 | 0.00370671 | 1.06864 | 0.0658343 | 0.987901 |
64 | 0.00235584 | 0.992152 | 0.00183888 | 1.01131 | 0.0330473 | 0.994307 |
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.015796 | 0.0349907 | 0.0332181 | |||
8 | 0.00847477 | 0.898313 | 0.0161629 | 1.11429 | 0.00795638 | 2.06179 |
16 | 0.00436788 | 0.956241 | 0.00804856 | 1.00588 | 0.00192546 | 2.04691 |
32 | 0.00221515 | 0.979529 | 0.00404351 | 0.993123 | 0.000551883 | 1.80277 |
64 | 0.00111531 | 0.98996 | 0.0020311 | 0.993347 | 0.000112974 | 2.28837 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.040389 | 0.0652755 | 0.512196 | |||
8 | 0.018239 | 1.14694 | 0.0182613 | 1.83775 | 0.257443 | 0.992443 |
16 | 0.00927265 | 0.975973 | 0.00777441 | 1.23198 | 0.130582 | 0.979297 |
32 | 0.00468607 | 0.984603 | 0.00370669 | 1.0686 | 0.0658357 | 0.988014 |
64 | 0.00235583 | 0.992143 | 0.00183888 | 1.0113 | 0.0330474 | 0.994333 |
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.0161227 | 0.0422295 | 0.00546336 | |||
8 | 0.00849478 | 0.924445 | 0.0186987 | 1.17531 | 0.00133666 | 2.03116 |
16 | 0.00436888 | 0.959313 | 0.00849608 | 1.13807 | 0.000250008 | 2.41859 |
32 | 0.00221526 | 0.979787 | 0.0042394 | 1.00294 | 7.36766e-005 | 1.7627 |
64 | 0.00111531 | 0.990031 | 0.00204415 | 1.05236 | 6.93936e-006 | 3.40833 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.0403803 | 0.0652523 | 0.500113 | |||
8 | 0.0182335 | 1.14706 | 0.0182584 | 1.83747 | 0.257152 | 0.959633 |
16 | 0.00927249 | 0.975563 | 0.00777436 | 1.23176 | 0.13057 | 0.977798 |
32 | 0.00468608 | 0.984575 | 0.0037067 | 1.06859 | 0.0658345 | 0.987908 |
64 | 0.00235583 | 0.992146 | 0.00183888 | 1.01131 | 0.0330473 | 0.994311 |
K | Non-stabilized | Standard grad-div | modular grad-div |
I | 0.0169173 | 0.0108085 | 0.077557 |
1e−1I | 0.0202974 | 0.0130437 | 0.0775562 |
1e−2I | 0.0464625 | 0.0301879 | 0.077554 |
1e−3I | 0.124238 | 0.0824988 | 0.0775521 |
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.0158906 | 0.0352882 | 0.042823 | |||
8 | 0.00847782 | 0.906408 | 0.0161872 | 1.12433 | 0.00978438 | 2.12983 |
16 | 0.00436799 | 0.956724 | 0.00805077 | 1.00765 | 0.00230198 | 2.08761 |
32 | 0.00221516 | 0.979559 | 0.00404369 | 0.993454 | 0.000609902 | 1.91623 |
64 | 0.00111531 | 0.989966 | 0.00203112 | 0.993397 | 0.000131401 | 2.2146 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.0404764 | 0.0653031 | 0.500655 | |||
8 | 0.0182477 | 1.14937 | 0.0182649 | 1.83808 | 0.25716 | 0.961151 |
16 | 0.00927325 | 0.976568 | 0.00777467 | 1.23222 | 0.130569 | 0.977854 |
32 | 0.00468612 | 0.984681 | 0.00370671 | 1.06864 | 0.0658343 | 0.987901 |
64 | 0.00235584 | 0.992152 | 0.00183888 | 1.01131 | 0.0330473 | 0.994307 |
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.015796 | 0.0349907 | 0.0332181 | |||
8 | 0.00847477 | 0.898313 | 0.0161629 | 1.11429 | 0.00795638 | 2.06179 |
16 | 0.00436788 | 0.956241 | 0.00804856 | 1.00588 | 0.00192546 | 2.04691 |
32 | 0.00221515 | 0.979529 | 0.00404351 | 0.993123 | 0.000551883 | 1.80277 |
64 | 0.00111531 | 0.98996 | 0.0020311 | 0.993347 | 0.000112974 | 2.28837 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.040389 | 0.0652755 | 0.512196 | |||
8 | 0.018239 | 1.14694 | 0.0182613 | 1.83775 | 0.257443 | 0.992443 |
16 | 0.00927265 | 0.975973 | 0.00777441 | 1.23198 | 0.130582 | 0.979297 |
32 | 0.00468607 | 0.984603 | 0.00370669 | 1.0686 | 0.0658357 | 0.988014 |
64 | 0.00235583 | 0.992143 | 0.00183888 | 1.0113 | 0.0330474 | 0.994333 |
1h | ||eu||L2 | uL2rate | ||eu||f | uHfrate | ||∇⋅eu||L2 | divuL2rate |
4 | 0.0161227 | 0.0422295 | 0.00546336 | |||
8 | 0.00849478 | 0.924445 | 0.0186987 | 1.17531 | 0.00133666 | 2.03116 |
16 | 0.00436888 | 0.959313 | 0.00849608 | 1.13807 | 0.000250008 | 2.41859 |
32 | 0.00221526 | 0.979787 | 0.0042394 | 1.00294 | 7.36766e-005 | 1.7627 |
64 | 0.00111531 | 0.990031 | 0.00204415 | 1.05236 | 6.93936e-006 | 3.40833 |
1h | ||eϕ||L2 | ϕL2rate | ||eϕ||p | ϕHprate | ||ep||L2 | pL2rate |
4 | 0.0403803 | 0.0652523 | 0.500113 | |||
8 | 0.0182335 | 1.14706 | 0.0182584 | 1.83747 | 0.257152 | 0.959633 |
16 | 0.00927249 | 0.975563 | 0.00777436 | 1.23176 | 0.13057 | 0.977798 |
32 | 0.00468608 | 0.984575 | 0.0037067 | 1.06859 | 0.0658345 | 0.987908 |
64 | 0.00235583 | 0.992146 | 0.00183888 | 1.01131 | 0.0330473 | 0.994311 |
K | Non-stabilized | Standard grad-div | modular grad-div |
I | 0.0169173 | 0.0108085 | 0.077557 |
1e−1I | 0.0202974 | 0.0130437 | 0.0775562 |
1e−2I | 0.0464625 | 0.0301879 | 0.077554 |
1e−3I | 0.124238 | 0.0824988 | 0.0775521 |