用分步有限元法求解三维不可压缩流动

ISSN 1000-0054
CN 11-2223/N　清华大学学报(自然科学版)J Tsinghua Univ (Sci &Tech),2000年第40卷第8期2000, Vol.40, No.829/33110 113
　
用分步有限元法求解三维不可压缩流动*
江春波,　邢秀英,　张庆海
(清华大学水利水电工程系,北京100084)
收稿日期: 1999-08-16
作者简介:江春波(1960-),男(汉),吉林,副教授
　*基金项目:国家自然科学基金项目(59979013)
文　摘:将计算二维不可压缩流动的分步有限元格式扩展
到三维情况,由于该格式没有引入新的高阶空间导数项,适
用于多维空间的非线性问题。实际求解了三维空腔流,在Re
数较低的情况下得到的流态不随时间而变化;在高Re数
(Re=4 000)的情况下,流态十分复杂并且是非定常的,在水
平断面具有多个非定常旋涡,在垂直于流动方向的立面断面
上可以模拟到TGL涡。本格式的迎风效应是Taylor展开式
的高阶精度项,没有人工粘性引入,得到的结果可靠性好。
关键词:分步格式;三维不可压缩流动; TGL涡
中图分类号:TV 131.41文献标识码:A
文章编号: 1000-0054(2000)01-0110-04
求解高Re流动的关键问题之一是解决由于对
流较强而引起的数值波动问题。为了获得稳定的数
值解,已经设计出了很多迎风格式。目前公认的精度
较高的迎风有限元格式有流线迎风格式(streamline
upwind/Petrov-Galerkin, SUPG )[1], Taylor-
Galerkin(T-G)格式[2]等。T-G格式是将计算在时
间方向上进行高阶Taylor展开,在空间用普通的
Galerkin离散得到的格式。已经证明该格式对于线
性对流方程具有三阶精度。但是T-G格式具有高阶
空间导数项,应用于非线性、多维空间问题比较困
难。本文用分步格式近似在时间方向的高阶Taylor
展开,得到的格式没有新的高阶空间导数项出现,适
用于非线性的多维问题。
围绕着如何满足不可压缩条件(连续条件)产生
了各种数值方法[3,4]。文中所介绍的是流速、压力分
离法,该法通过导出压力的泊松方程,使得连续条件
得到良好的满足,求解时对于流速和压力采用同阶
的插值函数,流速分量用显式格式、压力用隐式方法
求解。
应用该法求解二维非定常不可压缩流动的结果
可参考文[5],实际计算表明该法的计算效率较高,
比常规的方法节省计算时间。将该法推广到三维流
动,并分析了在三维情况下的数值精度与稳定性,计
算了三维空腔流,得到了非定常的TGL涡,与前人
的研究结果相符合。
1　控制方程和定解条件
流动控制方程为N-S方程和连续方程:
ui
t+ujui,j=-p,iρ+v(ui,j+uj,i),j+fi,(1)
ui,i= 0. (2)
式中:ui为流速,ui,j=uixj,p为压力,ρ为密度,v
为粘性系数,fi为外力。设边界由S1,S2组成,边界
条件可表示为:
ui=u∧i,
σij=-pδijρ+v(ui,j+uj,i) =σ∧

ij,
　在S1上,
　
　在S2上.
(3)
其中σij是应力。
2　数值方法
2.1　3步有限元方法
以如下的三维纯对流方程来说明本文的计算格
式:
f
t+αxfx+αyfy+αzfz= 0. (4)
式中:αx,αy,αz分别为x,y,z方向上的流速。在
求解f(t+Δt)时,根据Taylor展开原理将函数f在
时间步上进行展开到三阶精度,得
f(t+Δt) =f(t) +Δtft+Δt222ft2+
Δt3
6
3f
t3+O(Δt4). (5)若将式(4)直接带入(5)中,将出现新的高阶空间导
数项。本文的作法是通过分步形式来实现式(5),具
体步骤是通过如下的3步实现的:
f t+Δt3=f(t) +Δt3f(t)t,
f t+Δt2=f(t) +Δt2f(t+Δt/3)t,
f(t+Δt) =f(t) +Δtf(t+Δt/2)t.
(4)
空间离散仍采用标准的Galerkin有限元法,它保留
了T-G格式的三阶精度,但是不包含任何新的高阶
空间导数项,可用于解决非线性多维流动问题。三维
情况下的数值稳定性分析比较复杂,实际进行分析
的较少。分析过程类似于一维情况[5],可得到三维情
况下本格式放大因子为
G= 1 -12Y2-i-Y-16Y3, (7)
式中,参数Y定义为
Y=Crxf(ξ) +Cryf(η) +Crzf(ζ). (8)
式中:
Crx=Δtaxhx,Cry=Δtayhy,Crz=Δtazhz,
f(ξ)=sinξ
1-23sin2(ξ/2),
f(η)=sinη
1-23sin2(η/2),
f(ζ)=sinζ
1-23sin2(ζ/2).
这里,ξ,η,ζ分别为3个坐标轴方向的波数。
通过式(7)可以得到,在三维情况下分步有限元
格式同样具有三阶精度,并且在库郎数Cr=0～1的
范围内格式均收敛。
2.2　求解不可压缩流动的有限元格式
由于该方法精度高、稳定性好和计算简便的特
点,与T-G法相比,不包含任何新的高阶微分项,可
用于解决非线性多维流动问题。求解N-S方程的分
步有限元公式为
un+1/3i=uni+Δt3unit,
un+1/2i=uni+Δt2un+1/3it,
un+1i=uni+Δtun+1/2it.
(9)
将式(1)的N-S方程代入式(9)中,得到时间离散
格式:
un+1/3i-uni
Δt/3=-unjuni,j-pn,iρ+
v(uni,j+unj,i),j+fni,
un+1/2i-uni
Δt/2=-un+1/3jun+1/3i,j-pn,iρ+
v(un+1/3i,j+un+1/3j,i),j+fn+1/3i,
un+1i-uni
Δt=-un+1/2jun+1/2i,j-pn+1,iρ+
v(un+1/2i,j+un+1/2j,i),j+fn+1/2i.
(10)
式(10)中流速在时间上分3步计算,压力没有采用
分步处理。若将压力也分三步计算,模拟精度没有多
大提高,反而占用较多的计算时间[5]。采用标准的
Galerkin法进行空间离散,可得到最终的有限元方
程如下:
Mαβun+1/3iβ-uniβΔt/3=-Nnαβuniβ+
Liαβpnβρ-vSniα+Mαβfniβ+
∫ΓΦα-pnρ+v(uni,j+unj,i)·nidΓ,
Mαβun+1/2iβ-uniβΔt/2=-Nn+1/3αβun+1/3iβ+
Liαβpnβρ-vSn+1/3iα+Mαβfn+1/3iβ+
∫ΓΦα-pnρ+v(un+1/3i,j+un+1/3j,i)·nidΓ,
Mαβun+1iβ-uniβΔt=-Nn+1/2αβun+1/2iβ+
Liαβpn+1βρ-vSn+1/2iα+Mαβfn+1/2iβ+
∫ΓΦα-pn+1ρ+v(un+1/2i,j+un+1/2j,i)·nidΓ.
(11)
各个矩阵定义为:
Mαβ= ΩeΦαΦβdΩ,
Liαβ= ΩeΦα,iΦβdΩ

,
Nnαβ= ΩeΦαΦγΦβ,junjγdΩ,
Sniα= ΩeΦα,i(Φβ,juniβ+Φβ,iunjβ)dΩ.
(12)
Φα为权函数。由方程(11)可知,流速可由显式求出,
在计算un+1i之前,必须先解出pn+1。在方程(11)的最
后一式两边取散度,并引入不可压条件un+1i,i=0,可
以导出压力泊松方程:
pn+1,ii
ρ=uni,iΔt-un+1/2jun+1/2i,j+
111江春波,等:　用分步有限元法求解三维不可压缩流动v(un+1/2i,j+un+1/2j,i),ji+fn+1/2i,i(13)
应用Galerkin有限元法,可得到式(13)的弱形式:
ΩΦ,ipn+1,iρdΩ=-1Δt ΩΦuni,idΩ-
ΩΦ,iun+1/2jun+1/2i,jdΩ+
ΩvΦ,i(un+1/2i,j+un+1/2j,i),jdΩ+
ΩΦ,ifn+1/2idΩ+∫ΓΦb(u)nidΓ. (14)
式中b(u)为边界积分项。该文将流速和压力分开求
解,并对流速与压力采用相同的形函数,具有计算简
便的特点。在求解流速时,由于格式的稳定性好,可
以使计算步长比普通的有限元格式要长,因此计算
效率较高。上式可进一步整理成
ΩΦ,ipn+1,iρdΩ=-1Δt ΩΦuni,idΩ-
ΩΦ,iun+1/2jun+1/2i,jdΩ+ ΩΦ,ifn+1/2idΩ-
∫ΓΦun+1i-uniΔtnidΓ. (15)
则最终求解压力的有限元方程为:
Gαβpn+1β
ρ=-1ΔtHiαβuniβ-Kjαβun+1/2jβ-Sα, (16)
式中,各个矩阵定义为:
Gαβ= ΩeΦα,iΦβ,idΩ,
Hiαβ= ΩeΦαΦβ,idΩ,
Kjαβ= ΩeΦα,iΦβΦγ,jun+1/2iγdΩ,
Sα=∫ΓΦαun+1i-uniΔtdΓ. (17)
这样,从式(16)解出pn+1之后,就可以从式(11)解出
un+1。
3　计算结果与讨论
三维空间的空腔流是检验数学模型的基本算
例,其实是一个十分复杂的流动,已经观察到在高
Re情况下,流动不再根据几何形状保持对称性,流
场随时间在变化,在垂直流动的立面断面上形成随
时间变化的TGL(Taylor-Gortler-Like)涡[5～7],在
水平断面上也有多个随时间变化的旋涡。常规的数
值格式[8]不可能得到非定常的旋涡,能否模拟出三
维空腔流的TGL涡,是验证求解三维不可压缩流
动数值模型的一个关键。
用文中方法实际模拟了Re分别为100, 1 000
和4 000的流动。在Re分别为100, 1 000的情况
下,流态不随时间而变化;在Re=4 000的情况下,
得到非定常的结果。通过将文中格式的二维模拟结
果与前人的经典结果进行比较,验证了本格式的可
靠性。下面给出Re=4 000情况下的计算结果,计算
网格为30×30×30的均匀六面体单元。
图1所示为顺流方向立面(y=0.5)处不同时刻
的流速分布图,可见主旋涡的大小和形状在随着时
间不断地变化,在主旋涡的左、右下侧有次涡生成。
三维空腔流动的流速分布在水平断面上比图1
的流态要复杂,如图2所示为z=0.5处的水平断面
流速分布,在靠近前后的壁面附近处有旋涡生成,这
些旋涡不断发展变化,没有一个定型不变的状态。但
是可以看出流场关于中心轴对称。
图1　不同时刻顺流x-z断面的流速分布
图2　不同时刻x-y断面的流速分

布
如图3所示为不同时刻在垂直于流动方向的立
面断面上(x=0.5)的流态,数值模拟得到了TGL
旋涡,这是三维空腔流的一个十分重要的特征,该涡
的存在已经被实验所证实[4,6,7]。模拟得到旋涡主要
在边界底部附近生成,旋涡的大小与形状随时间在
不断地变化。文中的模拟结果与前人的实验结论符
合良好。事实上,分步法近似高阶导数项相当于迎风
有限元格式的迎风项,在这里迎风效应不为人工粘
性,因此得到的结果比较真实。且在程序的编制上比
较简便,适用于多维的非线性问题。
112清华大学学报(自然科学版) 2000,40(8)图3　不同时刻垂y-z断面的流速分布
4　结　论
基于T-G格式的分步有限元法保留了原T-G
格式的三阶精度和Courant数Cr等于1的稳定条
件。由于没有新的高阶项,故此法适用于有复杂边界
条件的非线性多维问题。由于算法简单,不必经过任
何修改就可用于三维问题。该法用于求解三维空腔
流,数值结果表明当Re=4 000时,流动为非定常,
在垂直流动的立面上出现TGL涡,算例表明此法
计算效率较高,计算结果与已有文献结果符合良好。
[参考文献]　References
[1] Johnson A A, Tezduyar T E. Parallel computation
of incompressible flows with complex geometries
[J]. Int J for Numerical Methods in Fluids, 1997,
24: 1321 1340.
[2] Donea J A. A Taylor-Galerkin method for convective
transport problems [J]. Int J Numerical Method
Engineering, 1984, 20: 101 119.
[3] Cao S L, Goulas A, Yakinthis K,et al. Numerical
simulation of three dimensional turbulent flow in a
contrifugal pump impeller [A]. In: Proc of the 3rd
Internatioanl Conference on Pumps and Fans [C].
Beijing, 1998. 411 418.
[4] Koseff J R, Street R L. Visualization studies of a
shear driven three-dimensional recirculation flow
[J]. J Fluids Eng, 1984, 106: 21-29.
[5] Jiang C B, Kawahara M. The analysis of
incompressible flows by a three-step finite element
method [J]. Int J Numerical Method In Fluids,
1993, 16: 793 811.
[6] Koseff J R, Street R L.On end wall effects in a
lid-driven cavity flows [J]. J Fluids Eng, 1984,
106: 385 389.
[7] Koseff J R, Street R L.The lid-driven cavity flow:
synthesis of qualitative and quantitative observations
[J]. J Fluids Eng, 1984, 106: 390 398.
[8] Ku Hu, Hirsh R S, Taylor T D. A pseudospectral
method for solution of the three-dimensional cavity
flows [J]. J Computational Physics, 1987, 70: 439
462.
Fractional step finite element
method for three-dimensional
incompressible flows
JIANG Chunbo,XING Xiuying,ZHANG Qinghai
(Department of Hydraulic and Hydropower Engineering,
Tsinghua University, Beijing 100084, China)
Abstract: The two-dimensional fractional step finite
element scheme was extended to three-dimensional cases.
The stability and numerical accuracy of the 3D case was
analyzed showing that the present fractional

step finite
element scheme has third order accuracy and an extended
stability domain. Three-dimensional cavity flows were
simulated. For small Reynolds numbers, the results are
steady. For large Reynolds numbers (i.e.,Re= 4 000),
the results are unsteady and the flows are very complicated.
The TGL vortices obtained in the present computation are
in good agreement with observations.
Key words: fractional step scheme; 3D incompressible
flow; TGL vortices
113江春波,等:　用分步有限元法求解三维不可压缩流动