摘要: 介绍了计算流体力学和常用求解方法的特点, 对于定常流动, 通过增加维数的方法把对象划归为同一性质的问题, 避免了控制方程的不确定性给求解造成的困难。以有限体积法为例说明了偏微分方程组向代数方程组转化的方法, 并以燃气轮机涡轮叶片流道内流场的计算为例说明了计算流体力学及求解方法的应用。
关键词: 计算流体力学; 控制方程; 求解方法; 有限体积法; 叶轮机械
中图分类号: O 354 　　文献标志码: A 　　文章编号: 1009 - 444X(2003) 01 - 0016 - 06
Application of Computational Fluid Dynamics
and It s Solution in Turbomachinery
LU J ia2hua
(College of Aeronautic Transportation , Shanghai University of Engineering Science , Shanghai 200336 , China)
Abstract : The significance of computational fluid dynamics (CFD) , which includes parabolic , elliptic and hy2
perbolic problem , is int roduced. Time2march method is adopted to make governing equations of fluid field a u2
niform pattern to steady problem , thus preventing the solving f rom the pre2judgment of the equation charac2
teristics , which simplifies the solution approach. As an example , the fluid field in a blade passage of a gas2tur2
bine is solved by means of Finite Volume Method (FVM) to demonst rate the application of CFD.
Key words : CFD ; Governing equations ; Solution approach ; FVM; Turbomachinery
　　流体的运动, 可以用一组非线性的偏微分方程组来描述。但要用解析法求解这些问题, 只有对极简单的情况方有可能。对于工程上感兴趣的问题, 经典的流体力学就无能为力了。对于这类要么是非线性的, 要么求解域相当复杂的实际工程问题, 只能求助于数值法来求解。
      随着计算机技术的发展, 计算流体力学已发展成为一门与实验流体力学并行的学科。在研究一种新的基本物理现象时, 实验一般说来是领先的, 计算总是力图精确地描述实验得到的物理现象。但是, 偶尔也曾出现过数值计算先于风洞实验指出某种尚未被发现的物理现象的事实。而且数值实验的费用要比实验费用低几个数量级。同时用数值法可以模拟实验难以做到的条件, 如模拟核反应堆的工作过程、爆炸、航空航天、环境保护、气象预报等问题。通过改变某些变量和初始条件, 可方便地得到不同工况下的工作状况。计算流体力学的主要微分方程基于质量守衡、动量守衡和能量守衡的自然规律。另一些方程则通过假定模型对方程组进行封闭, 如时间平均的湍流运动、因描绘湍流黏性引入的湍流模型等等, 这些复杂多变的模型只有通过数值法才能进行近似求解。计算流体力学的基本方法是将微分方程化简为代数方程组(线性或非线性) 的形式, 通过计算机程序求得计算单元(或节点) 上的值。因此可以说, 计算流体力学是近代流体力学, 数值计算方法和计算机应用技术三者有机结合的产物。本文给出描述流体运动的控制方程, 在比较几种常用求解方法的基础上, 以涡轮机械叶片流道流场计算为例, 选用有限体积法对控制方程进行离散和求解, 给出速度和压力分布, 说明计算流体力学和求解方法的具体应用。
1 　控制方程与性质
1. 1 　控制方程
根据Boussinesq 假设和k - ε湍流模型, 相对坐标系下流体运动控制方程通用形式如下:

(1)

等式左边第一项为瞬变项, 第二项为对流项, 等号右边第一项是扩散项, 第二项为源项。通用变量φ,扩散系数Γ和源项S 的表达式见表1 。

表1 　通用变量φ、扩散系数Γ和源项S 的表达式
Tab. 1 　Expression of general variable φ, diffusion coeff icient Γand source S

其中:
(2)

表中: x i 为笛卡儿坐标, i = 1 , 2 , 3 分别代表x , y , z 三个坐标方向; ui 是x i 坐标方向的时间平均速度分量; p 为压强;νef f , ν,νt 分别代表有效运动粘性系数、运动粘性系数、湍流运动粘性系数; ρ代表密度; k 为湍流动能; ε为湍动能耗散率; Cμ , C1 , C2 , σk , σε为湍流常数, 根据Spalding 等的建议值, 取Cμ = 0. 09 , C1 = 1. 43 , C2 = 1. 92 , σk = 1. 0 , σε= 1. 3 。
f i 为相对坐标系上的单位质量力:
(3)
式中Ω为旋转角速度。
1. 2 　控制方程的性质
      控制方程是考虑各方面的影响因素而得出的, 因而方程的各项都很复杂, 为使求解方法更有针对性, 对方程的性质必须了解。在流体流动中, 如果在坐标上的一个给定位置上处的条件, 要受该位置两侧条件变化的影响, 那么, 这个坐标就是一个双向坐标。反之, 一个坐标上的给定位置处的条件, 受该位置一侧条件变化的影响, 这个坐标就是单向坐标。
若在一个坐标方向上有限位置的单向流动, 各种主要的影响从上游传播到下游, 那么在某点的流体状态主要受其上游条件的影响, 而受其下游条件的影响很小。空间坐标的单向特性是一种近似, 当对流扩散的强度与Peclet 数相比很大(即高Reynolds 数问题) 时, 空间坐标就近似单向坐标了。
      计算方法上的单向与双向坐标这两个概念在微分方程的分类中, 分别为抛物形和椭圆形,即单向坐标对应的微分方程为抛物形, 双向坐标对应与椭圆形。双曲形不能准确地归入计算方法上的分类。双曲形问题具有一种单向的特性, 这种单向特性不是沿坐标的方向, 而是沿特征线的方向。对于特征线的计算方法只能用双曲形问题解决, 也可以把双曲形问题看成是椭圆形问题的一个组成部分(即全部是双向坐标) 。如由于射流泵有限射流具有一个明显的主流方向, 可以认为是抛物形流动; 有旋流场中, 用流函数作为求解函数求解亚音速流场的问题属于椭圆形问题; 超音速流动中, 始终存在着特征线和特征面,其控制方程是双曲形的。
      然而, 在定常的跨音速流场中, 在求出其解之前, 并不可能确切地知道场中任何一点上的流动是亚音速还是超音速的, 从而并不知道控制方程的类型, 给控制方程求解带来很大困难, 为此, 通过增维(时间相关法[1 ]) 的方法加以解决, 就是将定常问题转化为非定常的双曲形问题求解。
2 　计算流体力学的求解方法
2. 1 　主要求解方法的特点
      主要求解方法有:有限差分法、有限元法、边界元法、有限体积法和有限分析法。
      有限差分法是一个较为成熟的方法, 解的唯一性、收敛性、误差估计、稳定性等数学基础比较严密,发展比较完善。在处理高雷诺数问题时比有限元、边界元优越, 而程序较有限元和边界元简单。但边界必须规则, 否则必须采用贴体坐标系进行变换、计算较为复杂。
      有限元尺寸可根据需要来划分, 单元分析可以规格化、统一化, 接替能力强, 而且边界条件容易处理。
      但同差分法比较, 其统一性、收敛性、误差估计、稳定性等数学基础尚不严密, 计算结果易受人为因素的影响。
边界元法可将问题的维数减少一维, 计算工作量大为减少, 而且易处理无限区域问题。这种方法要求采用解析函数基本解, 目前还只使用于线性问题, 此外方程组的系数矩阵不对称, 并为满秩。

E , W , N , S , T , B 为主节点
图1 　三维有限体积示意图
Fig. 1 　Illustration of 3OD f inite volume

      上世纪70 年代发展起来的有限体积法(FVM) 将微分方程积分, 变为积分方程, 用单元中心点函数表示其平均值, 适用于复杂边界条件及不稳定流动情况。但要求单元要小, 否则误差较大。
2. 2 　有限体积法的应用
      大多情况下, 流体流动是不稳定的, 因此用有限体积法求解流场是比较合适的。采用有限体积法对控制方程进行离散化, 主要基于以下考虑:所得
到的结果将意味着在任何一控制容积内, 亦即在整个计算区域内, 诸如质量、动量及能量这样一些物理量的积分守衡都可以精确地得到满足, 对于任意数目的网格节点, 这一特征都存在(不只是限于网格节点变得很大时的极限意义上) , 因而即使是粗网格节点解也能保证严格的积分平衡, 数值处理也比较简便。
      用有限体积法离散控制方程, 就是将计算区域划分为一系列的组, 使其在不重合的有限体积内包含有一个计算节点P , 将控制方程在每个有限体积内积分, 得出一组离散方程, 其中的未知数是网格节点上的因变量φ的值。图1 为一个三维有限体积的示意图。
将方程(1) 的三维通用方程写为
(4)

式中J x , Jy , Jz 为总的对流加上扩散的流量, 定义为
(5)
式中u , v , w 表示在x , y , z 方向上的速度分量。对方程(2) 在图1 中所示的有限体积上积分, 得
(6)
其中源项S 可按Patankar[2 ] 推荐的方式进行了线性化, 并且对于非稳定项, 假定和在整个有限体内占控制地位。时间步长开始时的“旧值”用和 表示。按照全隐式的做法, 所有其它的没有角标的值认为是“新”值。量Je , J w , J n , Js , Jt , Jb 都是在有限体积各面上积分的总流量, 即 等等。
令式(1) 中φ= 1 , Γ= 0 , S = 0 , 可得连续方程:
(7)
以类似的方法, 在有限体积面上积分连续方程(7) , 得:
(8)
式中Fe , Fw , Fn , Fs , Ft , Fb 表示有限体积的各个面上流动强度, 分别为:

　如果用φp 乘以方程(8) , 并用方程(6) 减去它, 可得

式中, 在有限体积面上的所有量可用网格节点上的值(角标为大写字母表示) 表示, 即

整理后, 最终的离散化方程写为
(9)
式中:


其中P 为Peclet 数, 取为F 和D 之比, 如Pe = Fe/ De , 以此类推。各种格式的函数A (| P| ) 与选用的差分格式密切有关[3 ] 。扩散的传导性定义为

至此, 式(9) 将控制方程转化为了一组代数方程, 可用三对角阵算法( TDMA) 求解。
3 　计算实例
　　某双级燃气轮机透平叶片流道如图2 所示。

共两级, 按等内径设计。静叶为直叶片, 动叶为扭叶片, 轮毂直径700 mm , 汽缸壁的扩张角为15°, 透平初参数和主要热力参数参见表2 。
叶片流道采用贴体网格, 每叶片排按周向、径向和轴向划分为6 ×10 ×34 计算单元, 图3 给出的是第一级动叶网格划分情况。流场数值模拟
中, 各叶片排的出口静压沿径向分布按径向平衡方程得到, 其它边界条件和相关处理可参考文献[3 ] 。
图4 和图5 分别为叶片平均叶高处的静压分布和速度矢量。




4 　结　论
计算流体力学及求解方法的研究与计算机技术的发展使复杂流场的求解成为可能。以叶轮机械流场的计算为例, 在上世纪70 年代还处在二维和准三维的计算水平, 而今全三维、超跨音速黏性流场的求解也不再是难事。这对于优化叶轮机械设计, 提高机组热效率提供了理论依据, 其经济意义是明显的。随着理论模型的不断完善, 借助实验对模型的不断修正, 计算流体力学将会有更广的应用前景。
参考文献:
[1 ] 　DENTON J D. An improved time2marching method for turbomachinery flow calculation[J ] . ASME J of Engineering for power , 1983 , 105
(4) :514～524.
[ 2 ] 　帕坦卡S V. 传热和流体流动的数值方法[M] . 合肥:安徽科学技术出版社, 1984.
[ 3 ] 　鲁嘉华, 凌志光. 某双级跨音速燃气透平全三元黏性流场的数值分析[J ] . 能源技术, 2001 , 22 (4) :141～145.
　第1 期鲁嘉华:计算流体力学与求解方法在叶轮机械中的应用·21 ·　　
© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.

作者简介: 鲁嘉华(1960 - ) , 男, 浙江绍兴人, 副教授, 工学博士. 主要研究方向:计算流体动力学试验研究和数值计算.