周少玲， 侯　磊

(1. 上海大学 理学院， 上海 200444； 2. 河北工程大学 理学院， 河北 邯郸 056038)



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

周少玲1,2， 侯磊1

(1. 上海大学 理学院， 上海 200444； 2. 河北工程大学 理学院， 河北 邯郸 056038)

摘要:建立了求解稳态Oldroyd-B流体模型的有限元算法. 为了降低该流体模型的耦合性和有限元方程的求解规模， 考虑将Oldroyd-B流体模型解耦为Stokes方程和本构方程. 对于Stokes问题， 使用加权最小二乘有限元方法求解； 而对于含有对流项的本构方程， 则采用具有较好数值稳定性的流线迎风Petrov-Galerkin(SUPG)方法求解. 针对本构方程的非线性特点， 利用迭代算法将其进行线性化处理. 分析了解耦有限元解的先验误差估计. 通过粘弹性Oldroyd-B流体在管道内的流动问题， 验证了算法的有效性和收敛性.

关键词:Oldroyd-B流体； 解耦算法； 最小二乘有限元方法； SUPG方法

0引言

近几十年来， 随着科学技术的发展， 非牛顿流体的相关研究已经取得了巨大的进步， 尤其在化工、 石油、 医药等领域. 非牛顿流体的数值模拟始于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方程和本构方程. 利用迭代算法将问题线性化， 并分析算法误差.

1预备知识和解耦算法

设Ω是R2内的有界区域，Γ为其Lipschitz连续边界. 记Wm,p(Ω)为通常的Sobolev空间， 其范数为‖·‖m,p. 当p=2时，Wm,p(Ω)简记为Hm(Ω)， 相应的范数为‖·‖m. 显然，H0(Ω)=L2(Ω)， 其范数和内积分别为‖·‖和(·,·). 考虑如下稳态的不可压缩Oldroyd流体

(1)

式中：u和p分别为流体的速度和压强； T为额外应力张量； τ为粘弹应力张量； f为体积力；D(u)=(u+uT)/2 为应变率张量； 非负常数λ和α分别为We数和粘度比. 这里

且-1≤a≤1. 当a=1时， 称为Oldroyd-B流体. 将问题(1)解耦成两个子问题

(2)

和

(3)

定义未知量u,p,T和τ的解空间分别为

并记X=V×Q×S. 下面将分别采用最小二乘有限元方法和SUPG算法求解上述Stokes问题和本构方程.

2加权最小二乘有限元方法

其中， Pl(K)表示定义在K上次数不超过l的多项式函数空间， 且记Xh=Vh×Qh×Sh. 定义问题(2)的加权最小二乘泛函为

这里非负权重L用来增强数值解的质量守恒性[14].Stokes问题的最小二乘有限元解(uh,ph,Th)∈Xh定义为

(4)

由变分原理可得， 最小二乘有限元解(uh,ph,Th)需满足下面的Euler-Lagrange方程

(5)

其中

类似于文献[13]中的证明， 可得

引理 1如果(u,p,T)是(2)的解， 则最小二乘有限元解(uh,ph,Th)满足

(6)

式中：τh是本构方程(3)的有限元解.

3本构方程求解

本节将使用SUPG算法处理本构方程中的对流项. 引入下面的线性算子

(7)

当u=v时， 算子简记为B(u,τ,σ). 利用分步积分公式， 可得

(8)

令u=v，τ=σ， 由式(7)和式(8)， 可得B(u,τ,τ)=h‖(u·)τ‖2. 引入检验函数σu=σ+λh(u·)σ， 并定义范数‖‖σ‖2+‖λh1/2(u·)σ‖2. 将本构方程(3)与检验函数做内积得

定义粘弹应力张量τ的有限元空间为∑h={σ∈Σ;σ|K∈Pl(Ω),∀K∈Th}. 显然有限元解τh满足

(9)

假设Oldroy-B流体的准确解(u,p,T,τ)足够光滑， 即

且有

(10)

则有下面的误差估计成立.

定理 1若τ为本构方程(3)的准确解， τh为满足方程(9)的有限元解. 假设uh,τh∈L∞(Ω)， 且

(11)

如果λ足够小， 则有

(12)

则有下面的误差估计成立[15]

(13)

显然， 准确解τ满足

(14)

将式(9)减去式(14)， 可得

(15)

考虑式(15)等号左边， 若h<1， 则

另外

故式(15)等号左边满足

(16)

下面考虑式(15)等号右边各项. 第一项满足

(17)

利用式(8)和(13)， 得到

(18)

式(15)等号右端的第三项满足

(19)

类似地， 有

(20)

式(15)等号右端的最后一项满足

(21)

由式(17)～(21)得， 等式(15)的右端满足

综上可得

当λ足够小时， 有

证毕.

由式(6)和式(12)， 解耦算法的有限元解满足下面的误差估计式.

定理 2设(u,p,T,τ)∈V×Q×S×Σ是问题(1)的准确解， 如果α和λ足够小， 并且uh,τh∈L∞(Ω)， 则有限元解(uh,ph,Th,τh)满足

4数值算例

本节对粘弹性Oldroyd-B流体在管道内的流动进行数值模拟， 计算区域如图1(a)所示(其中x轴为管道的对称轴).

图 1　计算区域及其三角剖分Fig.1　Computational domain and triangular mesh

假设速度u=[u1,u2]T、 压强p和粘弹应力τ具有以下真解[13]

模型(1)右端由真解确定. 选取如下边值条件[13]

将区域Ω=[0,1]×[0,1]剖分成一致的三角单元(见图1(b))， 计算中分别取m=8,12,18,24. 所有未知量均采用分片连续线性多项式进行插值， 使用下面的解耦迭代算法计算：

(22)

2) 利用SUPG算法求解

(23)

图 2　程序流程图Fig.2　Program flow chart

图 3　x=0.5处的速度分量u1Fig.3　u1 along the line x=0.5

图 4　x=0.5处的速度分量τ11Fig.4　τ11 along the line x=0.5

5结论

本文研究了流变学中Oldroy-B流体模型的数值求解问题. 考虑到Oldroy-B流体的非线性性， 将模型解耦成两个子问题， 并利用迭代算法对其进行线性化处理. 对于Stokes问题使用加权最小二乘法求解， 并使用SUPG算法消除本构方程中对流项对于数值解的影响. 分析了解耦算法的误差. 通过一算例验证了本文算法的收敛性， 且数值解未出现振荡现象， 说明算法具有较好的稳定性.

参考文献:

[1]Baaijens F P T. Mixed finite element methods for viscoelastic flow analysis: a review[J]. Journal of Non-Newtonian Fluid Mechanics, 1998, 79: 361-385.

[2]Cai Z, Westphal C R. An adaptive mixed least-squares finite element method for viscoelastic fluids of Oldroyd type[J]. Journal of Non-Newtonian Fluid Mechanics, 2009, 159: 72-80.

[3]Frey S, Fonseca C, Zinani F, et al. Galerkin least-squares approximations for upper-convected Maxwell fluid flows[J]. Mechanics Research Communications, 2010, 37: 666-671.

[4]Lee H C. A nonlinear weighted least-squares finite element method for the Oldroyd-B viscoelastic flow[J]. Applied Mathematics and Computation, 2012, 219: 421-434.

[5]Lee H C, Chen T F. Adaptive least-squares finite element approximations to Stokes equations [J]. Journal of Computational and Applied Mathematics, 2015, 280: 396-412.

[6]Chaudhry J H, Bond S D, Olson L N. A weighted adaptive least-squares finite element method for the Poisson-Boltzmann equation[J]. Applied Mathematics and Computation, 2012, 218(9): 4892-4902.

[7]Bochev P B, Gunzburger M D. Least-squares methods for the velocity-pressure-stress formulation of the Stokes equations[J]. Computer Methods in Applied Mechanics and Engineering, 1995, 126: 267-287.

[8]Chen T F, Cox C L, Lee H C, Tung K L. Least-squares finite element methods for generalized Newtonian and viscoelastic flows[J]. Applied Numerical Mathematics, 2010, 60: 1024-1040.

[9]Lee H C. An adaptively refined least-squares finite element method for generalized Newtonian fluid flows using the Carreau model[J]. SIAM J. Sci. Comput., 2014, 36(1): 193-218.

[10]Fortin A, Guenette R, Pierre R. On the discrete EVSS method[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 189: 121-139.[11] Wang Kun. A new discrete EVSS method for the viscoelastic flows[J]. Computers and Mathematics with Applications, 2013, 65: 609-615.

[12] Brooks A N. Streamline upwind/petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[J]. Computer Methods in Applied Mechanics and Engineering, 1982, 32: 199-259.

[13] Chen T F, Lee H, Liu C C. Numerical approximation of the Oldroyd-B model by the weighted least-squares/discontinuous Galerkin method[J]. Numerical Methods for Partial Differential Equations, 2013, 29: 531-548.

[14] Bolton P, Thatcher R W. On mass conservation in least-squares methods[J]. Journal of Computational Physics, 2005, 203: 287-304.

[15] Najib K, Sandri D. On a decoupled algorithm for sloving a finite element problem for the approximation of viscoelastic fluid flow[J]. Numerische Mathematik, 1995, 72: 223-238.

文章编号：1673-3193(2016)04-0323-06

收稿日期：2015-09-29

基金项目：国家自然科学基金资助项目(11271247)； 河北省自然科学基金资助项目(G2013402063)

作者简介：周少玲(1979-)， 女， 讲师， 博士， 主要从事微分方程数值解研究.

中图分类号:O357.1

文献标识码:A

doi:10.3969/j.issn.1673-3193.2016.04.001

Decoupled Finite Element Algorithm for the Oldroyd-B Model

ZHOU Shao-ling1,2, HOU Lei1

(1. School of Science, Shanghai University, Shanghai 200444, China;2. School of Science, Hebei University of Engineering, Handan 056038, China)

Abstract:A new finite element algorithm is presented for the steady Oldroyd-B fluid. In order to reduce the coupling and the unknowns in the linear algebraic system, the fluid model is decoupled into a Stokes-like problem and a constitutive equation. The weighted least-squares finite element method is used to solve the Stokes-like problem. Due to the convective term, Streamline Upwind/Petrov-Galerkin method is applied to deal with the constitutive equation. The nonlinear terms in constitutive equation are linearized by an iterative algorithm. A prior error estimate for the decoupled algorithm is analyzed. For numerical test, the flow in a planer channel is considered. Numerical results show that the method is efficient and convergent.

Key words:Oldroyd-B fluid; decoupled algorithm; least-squares finite element method; SUPG