非定常Navier-Stokes方程的一种非线性局部投影稳定化有限元方法

李 西， 罗加福， 冯民富

(1.四川大学数学学院, 成都 610064； 2.成都体育学院， 成都 610041)

1 引 言

粘性不可压缩Navier-Stokes方程(NS方程)的混合有限元离散通常会出现两个难点:需满足经典的inf-sup条件和需要克服由高雷诺数(Reynolds number,Re)导致的数值震荡.

为了解决上述两个问题, 大量专家学者提出了不同的稳定化方法.比如，为了克服对流-扩散方程中对流占优时所引起的数值震荡,Brooks和Hughes[1]提出了Streamline Upwind Petrov-Galerkin方法(SUPG),并将其推广到不可压缩非定常NS方程.随后, Hughes 和Franca等[2]发现SUPG法不仅可以保持格式的相容性,不损失解的逼近性,还能增强离散解的稳定性——可利用SUPG 法避免inf-sup条件从而得到了速度和压力的等阶有限元插值.这种用来稳定离散压力的Petrov-Galerkin变分又被称为Pressure Stabilized Petrov-Galerkin(PSPG). 虽然SUPG/PSPG既可以做到避免inf-sup条件,又可以得到对高Re依然有效的先验误差估计，但需在经典的Galerkin变分格式中引入最小二乘项,导致稳定项中含有非物理耦合项,从而在用高阶元近似时稳定化格式中需要计算二阶导数[3-4].随后,Becker和Braack提出的Local Projection Stabilization方法(LPS)[5]很好地解决了上述不足.与SUPG相比, LPS中添加的稳定项中不含有速度-压力的耦合项,并且不需要计算二阶导数，也可以避免inf-sup条件和克服对流占优.正因如此,LPS从提出到现在得到了充分地发展.从文献[5]中首次采用LPS来处理Stokes问题之后，LPS被推广到了运输方程[6],Oseen方程的低阶离散[7],NS方程[8-9]以及其它一些相关工作[10-12].

在处理由高Re带来的数值震荡时,Burman和Fernandez等[13]采用了Bertoluzza在文献[14]中证明了的离散commutator性质,得到的误差估计随着流体Re的增大依然有效,即该数值格式可以克服由高Re带来的数值震荡.之后, Chen和Feng等[9]以及Frutos和Garcia等人[12,15]将上述技巧用在NS方程的LPS方法中,同样得到了对高Re依然有效的误差估计.

本文在文献[9,11-12,16]基础上,提出了一种基于非线性项的局部投影稳定格式，其中时间半隐格式中的非线性局部投影项等价于文献[16]中的SUPG-type稳定项和文献[11]中的高阶term-by-term稳定项(因而本文提出的非线性局部投影可看做是该两类稳定方法的非线性推广)，对速度-压力我们采用等阶元(Pk,Pk)进行逼近,并添加局部压力梯度投影稳定化项克服inf-sup条件.在分析离散解的先验误差估计时,我们采用Bertoluzza[14]提出的离散commutator性质，得到的误差估计右端项的常系数中不含有粘性系数的倒数1/ν,使得当流体的Re增大(即流体的粘性系数ν减小)时该误差估计依然有效.

2 半离散稳定化有限元格式

Hm(Ω)=Wm,2(Ω),

‖·‖m,M=‖·‖Hm(M),

‖·‖k,∞,M=‖·‖Wk,∞(M).

当M=Ω时，我们省略下标M,即‖·‖m=‖·‖Hm(Ω),‖·‖k,∞=‖·‖Wk,∞(Ω), 其中k=1,2是一个整数. 我们分别用Lp,Hs和Wm,p来表示相应于Lp,Hs和Wm,p的向量值Sobolev 空间,以及

设X为Sobolev空间, 定义映射φ(x,t):[0,T]→X,

并简记‖φ‖Lp(X)=‖φ‖Lp(0,T;X),p=2或∞.

取时间区间I=[0,T], 其中T为一个固定的正常数. 非定常不可压缩粘性流体的流动由以下非定常Navier-Stokes方程表示:

(1)

其中u=u(x,t)∈Rd,p=p(x,t)∈R,f=f(x,t)∈Rd分别表示不可压缩粘性流体流动的速度场, 压力场和外力场. 设特征长度L和特征速度U均为单位1. 于是该系统的运动粘度系数ν=Re-1, 其中Re为雷诺数.

求(u,p)∈X×Q, 满足

(∂tu,v)+b(u;u,v)+ν(∇u,∇v)-(p,∇·v)=

(2)

Xh(M)={vh|M:vh∈Vh,vh=0 on ΩM}.

由文献[17], LPS方法误差分析的关键在于存在一个插值算子j满足如下最优近似性质:

假设2.1令以下局部inf-sup条件成立：

β>0

(3)

由该局部inf-sup条件可以推导出插值算子ju和jp具有如下的正交性和近似性[17]:

引理2.2令假设2.1成立.则存在两个插值算子ju:X→Xh,jp:Q→Qh满足以下的正交性和近似性[18]:

(u-juu,vh)=0,∀u∈X,∀vh∈Dh,

(p-jpp,qh)=0,∀p∈Q,∀ph∈Dh,

‖u-juu‖0+h|u-juu|1≤

Chr+1‖u‖r+1,∀u∈X∩Hr+1(Ω)，

‖p-jpp‖0+h|p-jpp|1≤

Chr+1‖p‖r+1,∀p∈Q∩Hr+1(Ω)，

‖u-juu‖0,∞+h|u-juu|1,∞≤

Ch‖u‖1,∞,∀u∈X∩W1,∞(Ω)

(4)

下面我们给出NS方程的非线性局部投影稳定的空间半离散格式.∀t∈(0,T), 求(uh(t),ph(t))∈Xh×Qh, 满足,∀(vh,qh),

(∂tuh,vh)+ν(∇uh,∇vh)+b(uh;uh,vh)-

(ph,∇·vh)+Sconv(uh;uh,vh)=f,vh，

(qh,∇·uh)+Spres(ph,qh)=0

(5)

其中

Sconv(uh;uh,vh)=

其中α1,M和α2,M分别为相应于对流项和压力梯度项的局部投影稳定参数, 并且满足以下假设.

α1,M=α1hm1,α2,M=α2hm2，

其中m1,m2为稳定参数的阶,其具体数值在后文的误差分析中给出.

3 含稳定项的全离散格式

接下来我们对空间半离散格式(5)做时间上的有限差分.一般说来,时间差分上的全隐格式为无条件稳定,但每层时间上需解一个非线性方程组，全显格式在计算模拟时具有优势,但为了格式的稳定性需对时间步长Δt有诸多限制.一种通常的做法是线性项采用隐式格式,非线性项采用显示格式. 下面我们给出(5)的不同时间差分格式.

(6)

其中

注意到该格式中的非线性局部投影项实为线性型:

此时从程序实现时做数值积分的角度可知, 该项等价于以下两种线性化对流项的局部投影：

(1) 高阶term-by-term稳定项[11]

(2) SUPG-type的局部投影稳定项[16]

(κM((uM·∇)uh),κM((uM·∇)vh))M.

其中uM为流场速度u的局部投影平均

(7)

(8)

其中

下文我们将给出格式2的稳定性和收敛性分析.在前面的讨论中我们知道,格式1中的非线性局部投影稳定项可等价于文献[16]中提出的SUPG-type的稳定项和文献[11]中的高阶term-by-term稳定项.为简洁起见,我们略去格式3中离散解的稳定性和收敛性分析.

4 稳定性与收敛性分析

本节我们给出全离散格式(7)的稳定性与收敛性分析.为此我们需要一些必要的假设和不等式, 并且为了记法的简洁,我们定义

∀(vh,qh)∈Xh×Qh.

在算子Π的定义中，曾假设Π具有局部稳定性及近似性[18]，即

‖Π(v)‖0,M≤C‖v‖0,M.

注1有好几类算子满足上式给出的最优近似性和局部稳定性,如文献[19]中的Scott-Zhang-like算子以及文献[20]中的插值算子.

如文献[12,21]中所述, 误差分析要做到与1/ν一致的关键在于Bertoluzza所证明的离散commutator 性质, 即

引理4.2对于任意的u∈W1,∞(Ω),vh∈Vh，以下不等式成立:

‖(I-Π)(u·vh)‖0≤Ch‖u‖1,∞‖vh‖0.

(9)

再次利用Cauchy-Schwarz不等式得

利用上两式,再将(7)式乘上2Δt,并从0加到n得

对任意的1+Δt

利用

C0(N+1)≥Δt⟺

定理4.4(离散压力的稳定性) 对任意的ph∈Qh, 存在一个与h,Δt和ν无关的正实数β, 使得

β‖ph‖0

(10)

证明 由连续inf-sup条件可知, 对任意ph∈Qh⊂Q, 存在v∈X使得

∇·v=ph,|v|1≤C‖ph‖0.

由Green公式,(4)式中给出的算子ju的正交性和近似性以及Poincaré不等式可得

(∇·(v-juv),ph)+(∇·juv,ph)=

-(v-juv,κh∇ph)+(∇·juv,ph)≤

Ch|v|1·‖κh∇ph‖0+(∇·juv,ph).

又由ju:X→Xh有

定理得证.

定理4.5(离散速度的误差估计) 令(u,p)∈X×Q为问题(1)的解.我们假设∂tf∈L∞(L2)及(u,p)有额外的正则性,即(u,p)∈C0(Hs+1)×C0(Hs+1).又令存在常数C, 使得

‖∂tu‖L∞(Hs+1)+‖∂ttu‖L∞(L2)+

‖u‖L∞(W1,∞)≤C.

(11)

证明 记

E1+E2+E3+E4+E5+E6.

由Young不等式可得

由

Cauchy-Schwarz不等式, Young不等式及(4)式可得

关于对流项之差，我们将其分成如下形式：

E21+E22+E23+E24+E25,

这里我们利用了b(u;v,v)=0. 由Cauchy-Schwarz不等式,u的正则性,(4)式,逆不等式以及Young不等式有

类似地,有

及

对于第二项, 我们采用文献[9]中的技巧,有

令

在方程式(2)中取时间t=tn+1后减去(7)式, 取测试函数为

(v,q)=(vh,qh)=

后有

对于上式(4)右端第一项, 我们有

将(4)式及逆不等式应用于上式第一项可得

对于最后一项, 我们再一次利用不等式

‖juun‖Wk,∞≤‖un‖Wk,∞,k=0,1

以及(4)式, 逆不等式和算子Π的稳定性,得

由逆不等式和(4)式,我们得到了最后一项的估计如下：

由引理4.2有

由以上估计, 我们有

对于第三项E3, 运用Cauchy-Schwarz不等式, Young不等式以及(4)式有

然后, 由分部积分, 算子ju的正交性, Cauchy-Schwarz不等式及Young不等式有

对于E5, 我们将其分成两部分

由Young不等式可得

于是

C+ChM‖un+1‖1,M,∞+‖un+1‖1,M,∞.

对于剩下的部分, 由三角不等式及(4)式

从而由u∈L∞(W1,∞(Ω))可得

C(h2+1)h2s+m1‖u‖L∞(Hs+1).

于是

C(h2+1)h2s+m1‖u‖L∞(Hs+1).

类似可得

最后, 我们得到E6的估计

结合以上估计式,分别取m1,m2=0,2,再由逆不等式和范数‖·‖LPS的定义可得

在上式不等号左右同时乘上Δt,并且从0加到n,再由引理4.2可推出

最终, 由三角不等式和算子ju,jp的近似性, 定理得证.

定理4.6(离散压力的误差估计) 令速度误差估计中的条件成立,且时间和空间的划分具有相同的阶, 即存在与h和Δt无关的常数C, 使得Δt≤Ch.于是我们可得如下的压力误差估计：

C(1+(1+ν2+Eu(1+h-d))/β2)Eu+

C/β2(1+h2+Eu)Su

(12)

其中我们令离散速度的稳定性和误差估计的右端项分别为

和

(Δt)2.

证明 我们采用速度误差方程中相同的误差分解, 在(2)式中令时间t=tn+1后减去(7)式, 并令测试函数(v,q)=(vh,qh)=(vh,0)可得压力的误差方程

采用估计E1时的类似技巧和Poincaré不等式可得

不同于E2的估计, 这里我们采用文献[9]中给出的方法来估计对流项之差

|T2|=|b(un+1;un+1,vh)-

对上式第二项,我们采用分部积分得

T21+T22.

由逆不等式知

于是

由Cauchy-Schwarz不等式,‖·‖LPS的定义和速度误差估计的运用可以得到

至于T5,由离散速度的稳定性和离散Cauchy-Schwarz不等式可得

即

最后, 我们有

利用Δt≤Ch,结合以上的估计式可得

最终,由三角不等式,(4)式以及离散压力的稳定性，定理得证.