非定常微极流体方程的速度校正投影方法*

阿妮柯孜·奥斯曼，冯新龙，刘德民

(新疆大学数学与系统科学学院，新疆乌鲁木齐 830017)

0 引言

本文考虑了求解非定常微极流体方程的速度校正投影方法, 设区域Ω ⊂Rd(d=2 或3) 是具有光滑边界∂Ω的凸连通区域, 非定常微极流体方程为:

方程中未知函数是速度u, 角速度w, 压力p.ν ＞0 是通常的黏度系数, νr, c0, ca, cd是与非对称应力张量有关的系数, ca, cd和c0满足不等式c0+cd-ca＞0.方程中ν0=ν+νr, c1=ca+cd, c2=c0+cd-ca.当空间维数d=2时, 速度u=(u1,u2,0), 角速度w=(0,0,w3).当d=3 时, u=(u1,u2,u3), w=(w1,w2,w3).

微极流体在工业和工程中有重要的应用价值, 许多学者对其进行了研究.Ortega-Torres 等用全离散加罚有限元法证明了速度、压力和角速度的最优误差估计[1].Nochetto 等提出了微极流体方程的一阶半隐式全离散有限元方法, 使速度和角速度解耦求解, 并给出了算法的无条件稳定性和最优收敛阶[2].Salgado 研究了微极流体方程的全离散分步投影有限元方法, 在连续性方程中加入grad-div 稳定项可以改善质量守恒[3].Galdi 等研究了微极流体方程弱解的存在唯一性[4].非定常微极流体方程数值方法研究的主要困难包括强非线性、无散度限制及多场耦合性等, 使得变分问题成为典型的鞍点问题, 理论分析和数值离散时速度变量和压力变量需满足inf-sup 条件.投影方法在类似于文献[5-6]所提的偏微分方程数值解中有着广泛的应用,就是把原来的鞍点问题转化成椭圆问题进行求解, 从而降低了求解规模并改善了算法稳定性.Jiang 等提出了求解非定常微极流体方程的基于压力校正的投影方法并给出了算法的稳定性和误差估计[7], 第一步先求速度和角速度, 第二步是把速度投影到无散度空间.本文在文献[8] 的基础上给出了微极流体方程的速度校正方法, 算法的第一步是投影, 第二步是校正, 第三步求角速度.与文献[7] 的结果相比, 速度校正方法是一种全解耦的方法, 它不需要压力初值,也不需要速度和压力满足inf-sup 条件[9－10], 并推广到了非齐次Dirichlet 边界条件.

1 预备知识

这里引入一般的标量Sobolev 空间Hm(Ω) (m=0,1,2) 和向量Sobolev 空间Hm(Ω)d, 在不引起混淆的情形下记它们的范数为‖·‖m.特别的, 当m=0 时, 令L2(Ω)=H0(Ω) 和L2(Ω)d=H0(Ω)d, 对应的内积和范数为(·,·), ‖·‖, 记:

2 时间半离散速度校正方法

2.1 一阶时间半离散速度校正方法

下面估计式(10) 右边的项:

将上述估计代入式(10), 并把式(9) 和式(10) 相加, 则有:

当k=0,1,···,m-1 时, 对上面的不等求和并用离散形式的Gronwall 引理, 则结论成立.

为了进行收敛性分析, 引入以下符号:

定理2 令(uk+1,pk+1;～uk+1,wk+1) 是一阶时间半离散速度校正法的解, 0 ≤m ≤N, 则成立:

证明将t=tk+1代入式(1), 则:

将式(13) 的第一个方程减去式(3), 式(13) 的第二个方程减去式(6) 得到:

用u(tk+1)-ν0ΔtΔu(tk+1) 与式(5) 作差, 得到：

将误差方程式(14) 改写为:

2.2 二阶时间半离散速度校正方法

步骤1: 求uk+1, pk+1满足:

步骤2: 求～uk+1满足:

式(25) 也可以改写为:

证明与定理1 类似.

3 全离散速度校正方法

设τh={K} 是Ω 的一致正则三角网格, 且网格大小h=maxK∈τh{diam(K)}, 定义以下离散子空间:

3.1 一阶全离散速度校正方法

证明与定理1 类似.

3.2 二阶全离散速度校正方法

4 数值算例

本节对二维微极流体方程结构化和非结构化网格上协调有限元逼近问题的收敛性进行试验,并对三维微极流体方程进行了结构化网格上的收敛性试验.分别给出了能量稳定性和微极流体在轴承润滑中的应用.

4.1 收敛性试验

设Ω=[0,1]d(d=2, 3), 终止时刻T=0.1 并且ν=νr=ca=cd=1.一阶和二阶速度校正投影方法分别选取Δt=h2, Δt=h.

4.1.1 二维微极流体方程在结构化网格上的收敛性试验

考虑具有光滑解析解的流动问题, 方程(1) 的右端由精确解决定.

表1 和表2 分别给出了二维微极流体方程的一阶和二阶速度校正法的速度、角速度和压力变量在结构化网格上的数值结果.由表1 和表2 可知,一阶和二阶格式得到的速度、角速度、压力都能达到最优收敛阶.由此可知,速度校正方法是有效和可靠的.

表1 一阶速度校正法在Δt=h2 时的误差和收敛阶

表2 二阶速度校正法在Δt=h 时的误差和收敛阶

4.1.2 二维微极流体方程在非结构化网格上的收敛性试验

考虑非结构化网格上一阶和二阶速度校正方法的误差和收敛阶.其中, 网格特征尺寸分别为h=(h1,h2,h3,h4,h5)=(0.270 282, 0.146 33, 0.073 319, 0.035 511, 0.018 824 5).选取精确解为:

表3 和表4 分别给出了微极流体方程一阶和二阶速度校正法的速度、角速度和压力变量在非结构化网格上的数值结果.由表3 和表4 可知, 一阶和二阶格式得到的速度、角速度、压力都能达到最优收敛阶.由此可知, 速度校正方法是有效和可靠的.

表3 一阶速度校正方法在Δt=0.000 1 时的误差和收敛阶

表4 二阶速度校正方法在Δt=0.000 1 时的误差和收敛阶

4.1.3 三维微极流体方程在结构化网格上的收敛性试验

为了更好地说明速度校正投影法的有效性和可靠性, 给出了三维微极流体方程对应的收敛性检验.精确解如下:

如表5 和表6 所示, 给出了三维微极流体方程一阶和二阶速度校正法在非结构化网格上的数值结果.由表5 和表6 可知, 一阶和二阶格式得到的速度、角速度、压力都能达到最优收敛阶.

表6 二阶速度校正法在Δt=h 时的误差和收敛阶

4.2 能量稳定性

考虑含有周期边界条件的求解区域Ω=[0,1]2, 分析在给定初值条件下, 外力为f =g=0 时系统的能量耗散过程.选择初始条件

方程(1) 中, 当ν=νr=ca=cd=c0=1, T=2, h=1/60 时.给出了一阶和二阶速度校正投影方法的能量耗散图.

考虑系统能量En:=‖un‖2+‖vn‖2+‖wn‖2, 图1 给出了系统在不同时间步长下的能量耗散过程.四个能量曲线在所有时间步长上都呈现单调衰减.其中Δt=0.1,0.05,0.025,0.012 5.

图1 在不同时间步长下一阶(左)和二阶(右)格式的系统能量

4.3 微极性流体在轴承润滑中的应用

微极流体方程在轴承润滑问题中有着极其重要的应用.在本节中, 将具有非线性项的微极流体方程(u·∇)u和j(u·∇)w 加到方程(1) 的动量方程和角动量方程中, 并将对应的离散形式(uk·∇)uk+1和j(uk·∇)wk+1加到所有算法的第一步和第三步.显然, 非线性项的引入不会给数值分析带来额外的困难.其中d=2, ν=νr=ca=cd=c0=1.选取P1b-P1-P1b 有限元对.给出了不同转速Ωr下水平速度u1,垂直速度u2,压力和角速度用一阶速度校正法求解微极流体方程, 并给出微极流体方程在轴承润滑中应用的演化图.设外圆半径r1=1, 内圆半径r2=0.3, Γ1为外边界, Γ2为内边界.Γ1满足齐次边界条件u=0, Γ2满足u=(-Ωrr2sinθ,-Ωrr2cosθ), w=0,这里θ=arctan(y/x).分别取Ωr=100,500,1 000, Δt=0.01, T =1.

如图2 所示,随着转速Ωr的增大,流体的水平速度u1、垂直速度u2、角速度、压力也增大,同时流体的非对称性也逐渐明显.将线速度分量的大小与角速度的大小进行比较, 可知流体的极性效应仍然很小.当Ωr=1 000时我们的算法依然有效.

图2 Ωr=100, 500, 1 000时的水平速度、垂直速度、角速度和压力的演化

5 结论

本文提出了微极流体方程的速度校正投影法,针对线性的微极流体方程给出了稳定性分析和有限元逼近的误差估计结果,数值算例中, 将格式推广到非线性微极流体方程.数值算例中分别给出了解析解问题、周期边界的初值问题以及非齐次边界的同心轴承润滑问题.速度校正法可以有效降低求解规模, 数值算例验证了相应的结论.在实际问题的模拟中还需要根据具体的物理工况进行参数设置和格式调整, 这将在后期进一步研究.