Oldroyd-B流体的解耦有限元算法

Oldroyd-B流体的解耦有限元算法
周少玲;侯磊
【摘要】建立了求解稳态Oldroyd-B流体模型的有限元算法.为了降低该流体模型的耦合性和有限元方程的求解规模,考虑将Oldroyd-B流体模型解耦为Stokes 方程和本构方程.对于Stokes问题,使用加权最小二乘有限元方法求解;而对于含有对流项的本构方程,则采用具有较好数值稳定性的流线迎风Petrov-Galerkin(SUPG)方法求解.针对本构方程的非线性特点,利用迭代算法将其进行线性化处理.分析了解耦有限元解的先验误差估计.通过粘弹性Oldroy-B流体在管道内的流动问题,验证了算法的有效性和收敛性.
【期刊名称】《中北大学学报（自然科学版）》
【年(卷),期】2016(037)004
【总页数】7页(P323-328,334)
【关键词】Oldroyd-B流体;解耦算法;最小二乘有限元方法;SUPG方法
【作者】周少玲;侯磊
【作者单位】上海大学理学院,上海200444;河北工程大学理学院,河北邯郸056038;上海大学理学院,上海200444
【正文语种】中文
【中图分类】O357.1
近几十年来，随着科学技术的发展，非牛顿流体的相关研究已经取得了巨大的进步，尤其在化工、石油、医药等领域. 非牛顿流体的数值模拟始于20世纪70年
代[1]，虽然有许多算法已经应用于实际问题，但是依然有很多困难需要解决. 例
如随着We数的增大，本构方程的对流占优性和非线性耦合性会逐渐增强. 因此，需要建立有效稳定的算法处理非线性项和数值解的振荡问题[2].
最小二乘有限元算法已被广泛应用于粘弹性流体的数值求解[3-6]. 不同于混合有限元方法，其变分问题是通过极小化由所有方程残量构成的二次泛函得到的. 最小二乘有限元方法具有许多理论和计算上的优势，例如有限元空间无需满足LBB条件，并且离散化后的线性方程组总是对称正定的. Bochev和Gunzburger[7]将最小二
乘法用于求解Stokes方程，证明了数值解的收敛性. Chen等人[8-9]分别使用加
权最小二乘有限元方法求解Carreau流体和Giesekus流体.
Oldroyd-B流体运动方程是一类典型的非牛顿流体运动模型，常用来描述聚合物
流体、生物流体和悬浊液等的流动现象[10-11]. 在形式上， Oldroyd-B 流体运动方程可以看成是Navier-Stokes方程的一个扰动问题. 由于其本构方程的双曲特性，需要用稳定的算法处理对流项，例如可以考虑流线迎风Petrov-Galerkin(SUPG)
方法[12]和Discontinuous Galerkin(DG)方法[13].
本文旨在建立求解Oldroyd-B流体模型的有限元算法. 考虑将整个系统解耦成两
个子系统，分别使用加权最小二乘有限元算法和SUPG方法求解Stokes方程和
本构方程. 利用迭代算法将问题线性化，并分析算法误差.
设Ω是R2内的有界区域，Γ为其Lipschitz连续边界. 记Wm,p(Ω)为通常的Sobolev空间，其范数为‖·‖m,p. 当p=2时，Wm,p(Ω)简记为Hm(Ω)，相应的范数为‖·‖m. 显然，H0(Ω)=L2(Ω)，其范数和内积分别为‖·‖和(·,·). 考虑如下稳态的不可压缩Oldroyd流体
式中：u和p分别为流体的速度和压强；T为额外应力张量；τ为粘弹应力张量；f为体积力； D(u)=(u+uT)/2 为应变率张量；非负常数λ和α分别为We数和粘度比. 这里
且-1≤a≤1. 当a=1时，称为Oldroyd-B流体. 将问题(1)解耦成两个子问题
和
定义未知量u,p,T和τ的解空间分别为
并记X=V×Q×S. 下面将分别采用最小二乘有限元方法和SUPG算法求解上述Stokes问题和本构方程.
本节考虑问题(2)的有限元解法. 假设Ω为多边形区域， Th是其上的正则三角剖分，三角单元记为K. 设hK为K的直径，记hK. 定义有限元空间
其中， Pl(K)表示定义在K上次数不超过l的多项式函数空间，且记
Xh=Vh×Qh×Sh. 定义问题(2)的加权最小二乘泛函为
这里非负权重L用来增强数值解的质量守恒性[14]. Stokes问题的最小二乘有限元解(uh,ph,Th)∈Xh定义为
由变分原理可得，最小二乘有限元解(uh,ph,Th)需满足下面的Euler-Lagrange方程
其中
类似于文献[13]中的证明，可得
引理 1 如果(u,p,T)是(2)的解，则最小二乘有限元解(uh,ph,Th)满足
式中：τh是本构方程(3)的有限元解.
本节将使用SUPG算法处理本构方程中的对流项. 引入下面的线性算子
当u=v时，算子简记为B(u,τ,σ). 利用分步积分公式，可得
令u=v，τ=σ，由式(7)和式(8)，可得B(u,τ,τ)=h‖(u·)τ‖2. 引入检验函数
σu=σ+λh(u·)σ，并定义范数‖‖σ‖2+‖λh1/2(u·)σ‖2. 将本构方程(3)与检验函数
做内积得
定义粘弹应力张量τ的有限元空间为∑h={σ∈Σ;σ|K∈Pl(Ω),∀K∈Th}. 显然有限元解τh满足
假设Oldroy-B流体的准确解(u,p,T,τ)足够光滑，即
且有
则有下面的误差估计成立.
定理 1 若τ为本构方程(3)的准确解，τh为满足方程(9)的有限元解. 假设uh,τh∈L∞(Ω)，且
如果λ足够小，则有
其中，.
证明设为τ在有限元空间Σh上的正交投影，即
则有下面的误差估计成立[15]
显然，准确解τ满足
将式(9)减去式(14)，可得
令得
考虑式(15)等号左边，若h<1，则
另外
故式(15)等号左边满足
下面考虑式(15)等号右边各项. 第一项满足
利用式(8)和(13)，得到
式(15)等号右端的第三项满足
类似地，有
式(15)等号右端的最后一项满足
由式(17)～(21)得，等式(15)的右端满足
综上可得
当λ足够小时，有
证毕.
由式(6)和式(12)，解耦算法的有限元解满足下面的误差估计式.
定理 2 设(u,p,T,τ)∈V×Q×S×Σ是问题(1)的准确解，如果α和λ足够小，并且uh,τh∈L∞(Ω)，则有限元解(uh,ph,Th,τh)满足
本节对粘弹性Oldroyd-B流体在管道内的流动进行数值模拟，计算区域如图1(a)所示(其中x轴为管道的对称轴).
假设速度u=[u1,u2]T、压强p和粘弹应力τ具有以下真解[13]
模型(1)右端由真解确定. 选取如下边值条件[13]
将区域Ω=[0,1]×[0,1]剖分成一致的三角单元(见图1(b))，计算中分别取
m=8,12,18,24. 所有未知量均采用分片连续线性多项式进行插值，使用下面的解耦迭代算法计算：
1) 由前一次的迭代得到，使用加权最小二乘法求解
2) 利用SUPG算法求解
式中：由1)求得. 迭代算法的程序流程图如图 2 所示.
为了验证算法的收敛性，取流体模型(1)中的参数a=1，α=1/9，λ=0.25. 设定迭代初值=0，迭代终止准则为‖)‖<10-3，最大迭代次数为150. 对于参数L，当L取值较大时，会使得线性方程组的条件数变大；当L取值较小时，不能保证数值解的质量守恒性，故取L=10. 图 3，图 4 给出了使用4种不同网格，解耦有限元算法在x=0.5处的计算结果(u1和τ11).
本文研究了流变学中Oldroy-B流体模型的数值求解问题. 考虑到Oldroy-B流体的非线性性，将模型解耦成两个子问题，并利用迭代算法对其进行线性化处理. 对于Stokes问题使用加权最小二乘法求解，并使用SUPG算法消除本构方程中对流项对于数值解的影响. 分析了解耦算法的误差. 通过一算例验证了本文算法的收敛性，且数值解未出现振荡现象，说明算法具有较好的稳定性.。