1. 引言
非定常渗透对流模型指的是由流体力学中的非定常不可压缩Navier-Stokes方程和热传导方程所耦合的非线性偏微分方程组,在自然世界和工业制造中有广泛的应用[1]-[3]。对于有界的凸多边形区域
和正常数
,未知函数速度u、温度
和压力p在
内满足如下的非定常渗透对流模型:
.(1)
.(2)
.(3)
其中
是单位基向量,正常数
和
分别是粘性系数和导热系数。
对该模型,我们考虑如下的初始条件和边界条件:
.(4)
.(5)
关于非定常渗透对流模型数值方法及其收敛性已有一些研究工作。Ravindran在[4]中使用二阶BDF格式分析了该模型的时间和空间离散化的外推两步后向差分,给出了全离散有限元格式的最优误差估计。Cao等在文章[5]中采用Mini元和分片线性有限元分别逼近速度场、压力场和温度场,构造了线性化的Crank-Nicolson有限元格式,得到了L2范数下的最优误差估计。Wan等在文献[6]中基于grad-div稳定化方法,研究了非定常渗透对流模型的一阶后向Euler和二阶BDF2有限元全离散格式,所得到的误差估计界不依赖于粘性系数和导热系数。
由于非定常渗透对流模型是由Navier-Stokes方程和热传导方程所耦合,使得构造求解该模型的数值格式时,需考虑Navier-Stokes方程的不可压缩性和非线性性,这也是构造求解Navier-Stokes方程数值方法的难点之一。常用的方法之一为投影方法[7] [8],然而投影方法中过渡的中间速度场不满足原方法的边界条件,导致数值边界层的出现。Blasco[9]等人构造了求解Navier-Stokes方程的一阶分数步长算法,在每个时间步上,该算法的基本思想是:首先通过求解一非线性椭圆方程得到过渡的中间速度场,然后通过求解一广义Stokes方程得到该时间步上的速度场和压力,从而实现方程的不可压缩性和非线性的分离,这是求解Navier-Stokes方程的一个高效稳定数值算法。基于该构造思想,Wu等在[10]中研究了求解由非定常Navier-Stokes方程和定常Maxwell方程所耦合的混合MHD系统的一阶分数步长有限元算法。An在[11]中研究了三维变密度不可压缩Navier-Stokes方程的一阶分数步长时间离散格式。目前,还没有关于求解非定常渗透对流模型分数步长算法的研究。
本文将基于Blasco等所构造的一阶分数步长算法,构造求解非定常渗透对流模型的分数步长算法,使得所构造的数值算法是无条件稳定的。同时,在解的正则性假设下,我们给出了速度场和温度场的一阶时间误差估计。
2. 预备知识
首先引入下列符号。当
,
时,记
为Sobolev空间[12]。当
时,
定义为
。
表示
内积。用
,
和
分别表示
,
和
的范数。
记:
,
,
,
,
,
其中n表示边界
上的单位向外法向量。
定义如下三线性项:
,
。
当
时,通过分部积分,容易验证三线性项
和
满足如下反对称性:
,
且有:
.
记
是时间区间
的一个均匀化分割,对于每个时间步长
和
,
。下面引入离散的Gronwall不等式[13]:
引理2.1对整数
,记
和B都是非负数。若成立:
.
则当
时有:
,
其中
。
3. 数值格式
下面给出求解非定常渗透对流模型的一阶分数步长时间离散算法。对
:
第一步:求
使得:
.(6)
第二步:求中间速度
使得:
,(7)
.(8)
第三步:求
和
使得:
,(9)
.(10)
下面给出上述分数步长算法的稳定性分析。
定理3.1假设正则性(A)成立,对所有的
和
有:
,(11)
和
(12)
证明:将(6)与
作内积,利用公式
,则会有:
(13)
将式(13)从
加到
(
),由于解的正则性会有:
(14)
另一方面,将(9)与
作内积可得:
(15)
将式(7)与
作内积可得:
(16)
将式(15)和(16)分别从
加到
,可以得到:
(17)
通过上述推导和分析,定理3.1成立。证毕。
4. 时间误差估计
记速度、温度以及压力的误差函数为:
(18)
当
时,容易验证真解满足如下的离散抛物系统:
,(19)
,(20)
其中截断函数为:
下面我们开始证明所构造的一阶分数步长算法具有一阶时间收敛精度。为此,我们假设真解满足如下正则性:
(A)
定理4.1假设正则性(A)成立,则当
和对足够小的
,我们有:
(21)
证明:利用(19)减去(6),可以获得:
.(22)
对于非线性项,我们做如下处理:
将(22)与
作内积,用Hölder不等式和Young不等式估计其右端项,对每个分量的处理如下:
1) 非线性项
,
2) 截断误差项
,
根据上述推导,我们有:
(23)
另一方面,利用(20)减去(7),可得:
(24)
将(24)与
作内积,对(24)的右侧非线性项的处理如下:
(25)
类似地,对于式子(24)的最后两项可以分裂成:
(26)
这样整理到(24)就会有:
(27)
下面分别对(27)的右端每一项进行估计:
1) 截断误差项
,(28)
2) 压力梯度项
,(29)
这里用到了
,
3) 非线性项
,(30)
,(31)
4) 与温度相关的项
,(32)
,(33)
,(34)
,(35)
,(36)
根据上述推导,将(28)~(36)代入(27)中得到:
(37)
将式(9)重新写为:
,(38)
将(38)与
作内积,再由于
,可得:
(39)
将(23),(37)和(39)从
加到
,可以得到:
(40)
对(40)利用引理2.1和连续解的正则性质(A),使得定理4.1成立。证毕。
定理4.1得到了在
中一致稳定的速度和温度,换句话来说就是,存在一个与时间步长
不相关的常数C (
),使得对所有
,可得:
(41)
定理4.2 假设正则性(A)成立,则当
和对足够小的
有:
(42)
证明:将(22)与
作内积,可以得到:
(43)
下面分别对(43)的右端的每一项进行估计:
1) 截断误差项
,(44)
2) 非线性项
,(45)
.(46)
将式(7)和(9)相加,可得:
.(47)
将式(20)减去(47),可得:
(48)
将(48)与
作内积,又由于在原先的假设下有
,可以得到:
(49)
类似式(25)的分裂,对(49)的右端项进行估计:
1) 截断误差项
,(50)
2) 非线性项
(51)
下面分别估计
和
,由于
的反对称性,
,可得:
(52)
.(53)
下面估计其余的非线性项:
,(54)
(55)
3) 与温度相关的项
(56)
(57)
将(43)和(49)从
加到
,再由(44)~(57)可以得到:
(58)
对足够小的
,可以对式子(58)应用引理2.1,最后一项被左边吸收,定理4.2成立。证毕。
定理4.3 假设正则性(A)成立,则当
和对足够小的
有:
.(59)
证明:重写(48)为:
(60)
利用离散的LBB条件:
,(61)
对任意的
,需要估计(60)的右边:
,(62)
,(63)
.(64)
对于非线性项,有如下的分裂:
(65)
下面分别估计G1,G2和G3:
(66)
,(67)
.(68)
对于式(60)右端剩下的两项,可以进行如下分裂:
(69)
下面分别估计(69)的右端项,可得:
(70)
,(71)
,(72)
,(73)
.(74)
因此,整理(62)~(74),可以得到:
(75)
从而也会有:
(76)
不等式(59)可以由定理4.2和与连续性的解u相关的(A)给出。证毕。
5. 数值结果
下面将给出数值实验来验证理论推导的最优误差估计。为了简便,考虑在单位正方形[0, 1] × [0, 1]上的问题(1)~(5)。在计算过程中,通过取适当的f和g,使得真解有如下形式:
,
,
.
Table 1. Numerical error and convergence order under L2
表1. L2范数下的数值误差和收敛阶数
M |
|
阶数 |
|
阶数 |
|
阶数 |
10 |
1.11E−03 |
− |
5.62E−03 |
− |
1.10E−03 |
- |
20 |
5.49E−04 |
1.01 |
2.82E−03 |
0.99 |
4.00E−04 |
1.46 |
30 |
3.57E−04 |
1.07 |
1.87E−03 |
1.02 |
2.22E−04 |
1.45 |
40 |
2.61E−04 |
1.09 |
1.39E−03 |
1.03 |
1.46E−04 |
1.45 |
60 |
1.67E−04 |
1.10 |
9.06E−04 |
1.05 |
8.15E−05 |
1.44 |
80 |
1.22E−04 |
1.10 |
6.64E−04 |
1.08 |
5.43E−05 |
1.41 |
Table 2. Numerical error and convergence order under L2
表2. L2范数下的数值误差和收敛阶数
M |
|
阶数 |
|
阶数 |
10 |
4.43E−03 |
− |
9.47E−03 |
- |
20 |
1.62E−03 |
1.46 |
4.16E−03 |
1.19 |
30 |
9.46E−04 |
1.32 |
2.86E−03 |
0.92 |
40 |
6.76E−04 |
1.17 |
2.30E−03 |
0.77 |
60 |
4.55E−04 |
0.97 |
1.76E−03 |
0.65 |
80 |
3.62E−04 |
0.79 |
1.49E−03 |
0.59 |
为了验证收敛阶,固定
,取时间步长
,其中M为正方形边长上的节点数。对于定理4.2中给出的时间误差估计,分别取
,得到的数值误差及收敛阶如表1和表2所示。从表1和表2可以看出,当网格逐渐变小时,得到的一阶收敛与定理4.2中的理论相一致。