求解不可压缩流动的分步有限元格式

第一篇：求解不可压缩流动的分步有限元格式

求解不可压缩流动的分步有限元格式.txt丶︶￣喜欢的歌，静静的听，喜欢的人，远远的看我笑了当初你不挺傲的吗现在您这是又玩哪出呢？ISSN 1000-0054

CN 11-2223/N 清华大学学报(自然科学版)J Tsinghua Univ(Sci &Tech),2002年第42卷第2期2002, Vol.42, No.236/37278-280

求解不可压缩流动的分步有限元格式

江春波, 徐照明, 李秀丽

(清华大学水利水电工程系,北京100084)

收稿日期: 2000-10-23

基金项目:国家自然科学基金资助项目(59979013);四川大学高速

水力学国家重点实验室访问学者资助项目

作者简介:江春波(1960-),男(汉),吉林,副教授。

摘 要:提出求解不可压缩Navier-Stokes方程的分步有限

元格式,该格式没有高阶微分项产生、程序编制简单,适用于

非线性的多维复杂流动。应用该方法实际模拟了二维圆柱绕

流的旋涡形成与脱落过程,得出了不同Re情况下圆柱绕流的流速分布。计算得到的不同Re下的旋涡脱落频率

(Strouhal数)与前人已有的经典解答符合良好。

关键词:分步有限元法;不可压缩流动;圆柱绕流

中图分类号:TV 131.41

文章编号:1000-0054(2002)02-0278-03文献标识码:A

Fractional step finite element formulation

for solving incompressible flows

JIANG Chunbo,XU Zhaoming,LI Xiuli

(Department of Hydraulic and Hydropower Engineering,Tsinghua University,Beijing 100084,China)

Abstract: A fractional step finite element scheme is proposed for

unsteady incompressible flows.Since no higher order terms are

introduced in the presented computation, this method is suitable for

nonlinear multi-dimensional problems.The unsteady flow around a

two dimensional circular cylinder was analyzed to predict the flow

patterns for different Reynolds numbers.The predicted Strouhal

number for vortex shedding is in good agreement with previous

result.Key words: fractional step finite element method;incompressible

flow;flow around cylinders

求解高Reynold数流动的关键问题之一就是解

决由于对流较强而引起的数值波动问题。为了获得

稳定的数值解,已经设计出了很多迎风格式。目前公

认的精度较高的迎风格式有流线迎风法

(Streamline Upwind/Petrov-Galerkin)[1], Taylor-

Galerkin法[2]和Lax-Wendroff有限元法[3]等,这

些方法都已经应用于求解不可压缩流动和对流扩散

问题[4]。

对于不可压缩流动,由于压力项不在连续方程

中出现,使得数值求解比较困难[5]。本文将压力和流

速分开求解,对压力和流速采用同阶的形函数进行

空间离散,压力通过导出的泊桑方程进行求解,流

速通过显式格式进行求解。

不可压缩粘性流体绕圆柱的流动,是流体力学的经典问题之一,较多的研究者在理论与实验两个

方面进行了大量的工作[6～8]。模拟圆柱绕流流动需

要数值精度高且稳定性好的数值格式,由于本文提

出的基于Taylor展开的分步有限元法数值稳定性

好,并且没有人工调节参数,适用于高Re不可压缩

流动的数值求解。用该方法能够模拟出圆柱绕流尾

流区的非恒定交错排列的Karman涡街,不同Re下

计算得到的涡街脱落频率与前人的经典结果[7,8]符

合良好。由于本文的分步有限元方法具有较大的数

值稳定区域,计算效率较高[2,4]。控制方程和定解条件

不可压缩流动由Navier-Stokes(N-S)方程和连

续方程控制

ui

t+ujui,j=-p,iρ+ν(ui,j+uj,i),j+fi,i= 1,2,3,(1)

ui,i= 0.(2)

式中:ui为流速,p为压力,ρ为流体的密度,ν为

流体的运动粘性系数,fi为外力。设边界由S1,S2

组成,边界条件可表示为

ui=u∧i,在S1上,σij=-pδijρ+ν(ui,j+uj,i)=σ∧ij,在S2上.(3)

这里,σij是应力,u∧i,σ∧ij为已知流速和应力。由于圆柱附近流场比较复杂,含有流动分离和旋涡脱落等

现象,需布置较密的网格;而在远离圆柱的区域,则

可布置较疏的网格,以节省计算工作量。有限元离散格式

流速分量按如下的分步格式求解

un+1/3i=uni+Δt3unit,un+1/2i=uni+Δt2un+1/3it,un+1i=uni+Δtun+1/2it.(4)

从式(4)中可以看出,流速在时间方向上的离散逼近

Taylor展开的三阶精度项。由式(1)和式(4),得到

关于N-S方程的时间离散格式如下

un+1/3i-uni

Δt/3=-unjuni,j-

pn,iρ+ν(uni,j+unj,i),j+fni,un+1/2i-uni

Δt/2=-un+1/3jun+1/3i,j-

pn,iρ+ν(un+1/3i,j+un+1/3j,i),j+fn+1/3i,un+1i-uni

Δt=-un+1/2jun+1/2i,j-

pn+1,jρ+ν(un+1/2i,j+un+1/2j,i),j+fn+1/2i.(5)

采用标准的Galerkin法对式(5)进行空间离散。注

意关于时间导数项的系数矩阵并不是集中质量矩

阵,因此本有限元格式的空间离散精度较高[4]。关于

流速的求解,可以不必形成总体系数矩阵,在单元系

数矩阵的基础上就可以直接求解。根据质量系数矩

阵具有对称及对角占优之特点,用简单雅可比迭代

法进行2至3次迭代就可以获得收敛解。在计算流

速un+1i之前,必须先解出压力pn+1。在方程(5)的最

后一式两边取散度,并引入不可压条件,可以导出压

力泊松方程

pn+1,ii

ρ=uni,iΔt-(un+1/2jun+1/2i,j),i+

ν(un+1/2i,j+un+1/2j,i),ji+fn+1/2i,i.(6)

应用Galerkin有限元法可对式(6)进行空间离散。

本文将流速和压力分开求解,并对流速与压力采用

相同的形函数,具有计算简便的特点。在求解流速

时,由于格式的稳定性好[4],可以使计算空间步长比

普通的有限元格式要长,因此计算效率较高。同时每时间步内显式求解三次流速值,用隐式法求解一次

压力,在求解压力时可以用较大的时间步长,因此本

有限元法比较节省计算时间。计算结果与讨论

3.1 计算区域与计算网格

为了使两侧边界不影响圆柱附近的流动,计算

区域在横向(垂直来流方向)取为10d, d为圆柱的直径;在纵向取为27d。由于圆柱附近流场比较复

杂,含有流动分离和旋涡脱落等现象,须布置较密的网格;而在远离圆柱的区域,则可布置较疏的网格,以节省整个计算的工作量,又考虑到圆与矩形的衔

接,故这里采用三角形网格,图1为圆柱绕流网格划

分示意图,总结点数为2 815,划分单元数为5 424。

图1 计算区域与网格划分

3.2 不同Re下流动情况的比较

定义Strouhal数为Sr=fd/u0,其中f为旋涡的生成脱落频率, d为圆柱的直径, u0为均匀来流

流速。图2给出了不同Re下计算得到的流线图。在Re较小(Re=10, t=200s)时,流动在圆柱背后形

成对称的旋涡。随着Re的增高,圆柱后面的旋涡逐

渐趋于非对称,在圆柱后面形成交错排列的Karman涡街。用本文的有限元格式计算,不必给人

工扰动就可以模拟出非恒定的Karman涡街。图2

中分别给出了Re=10(t=200s), 100(t=180s),Re=1 000(t=80s)的流线分布情况。从计算结果中

可见,旋涡的生成脱落半周期约为6s,相应的Sr为

0.166,与Williamson[7]的结果符合良好;在Re=1

000的情况下,对流作用强,流态变化复杂,用本文的分步有限元格式可以得到稳定的数值解。计算过

程没有人工粘性系数的导入,也没有引入需要调整的人工参数。

3.3 固定Re不同时刻流线图

为了形象地考察圆柱后面的旋涡生成与脱落过

程,计算了Re=200的流动,图3给出了时间t=

154～159s每秒时的流线图。注意t=154s时在圆

柱的右下方有一个小旋涡生成,随着时间的增长这

个小旋涡在不断的扩大并向圆柱的右上方移动,在t=159s时这个旋涡破裂,同时一个新的小旋涡在279江春波,等: 求解不可压缩流动的分步有限元格式图2 不同Re下圆柱附近流线 圆柱的右上方生成。这一过程记为旋涡生成脱落的半周期,时间间隔约为5.5s,相应的Sr为0.188,与前人的经典结果(图4)[7]符合良好,也符合半周期与Re的关系曲线。

图3 不同时刻圆柱附近流线分布(Re=200)

图4 文[7]的流线图(Re=200)结 论

基于Taylor展开的分步有限元格式由于没有

新的高阶项,故此法适用于有复杂边界形状的非线

性多维问题。本方法具有算法简单,计算效率高的优

点。本文将该有限元法应用于求解不可压缩流动,得

到了圆柱绕流交错排列的Karman涡街。算例表明

此法计算效率较高,计算结果与已有文献结果符合良好。

参考文献(References)

[1] Tezduyar T E, Ganjoo D K.Petrov-Galerin formulations

with weighting functions dependent upon spatial and temporal

discretion: application to transient convection-diffusion

problems [J].Comput Methods Appl Mech Eng, 1985,59:

71.[2] Donea J, Quartapelle L, Selmin V.An analysis of time

discretization in the finite element solution of hyperbolic problems [J].J Comput Phys, 1978,70: 463 499.[3] Kawahara M.Convergence of finite element Lax-Wendroff method for linear hyperbolic differential equation [A].Proc JSCE [C].Tokyo, 1976, 253: 95 107.[4] Jiang C B, Kawahara M, Kashiyama K.A Taylor-Galerkin based finite element method for turbulent flows [J].Fluid Dynamics Research, 1992,9: 165 178.[5] Gresho P M, Chan S T, Lee R L, et al.A modified finite element method for solving the time-dependent,incompressible Navier-Stokes equations, Part 1: Theory [J].Int J Num Methods in Fluids, 1984,4: 557 598.[6] Johnson A A, Tezduyar T E.Parallel computation of incompressible flows with complex geometries [J].Int J Numerical Methods in Fluids, 1997,24: 1321 1340.[7] Williamson C H.Three-dimensional wake transition [J].J Fluid Mech, 1996,328: 345 407.[8]方正,李炜.粘性流体圆柱绕流的混合有限分析解[J].水动 力学研究与进展, A辑, 1998, 13(2): 133 139.FANG Zheng, LI Wei.Simulations of Viscous flow around

circular cylinder by finite analysis method [J].J Hydraulic Research,Series A, 1998,13(2): 133 139.(in Chinese)280清华大学学报(自然科学版)2002,42(2)