收稿日期：2002-02-06；修订日期：2002-04-18
作者简介：陈广业（1976-），男，江苏张家港人，北京航空航天大学动力系硕士生.
第18卷第1期2003年2月
航空动力学报
Journal of Aerospace Power
VoI.18No.1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Feb.2003
文章编号：1000-8055（2003）01-0002-06
使用GAO-YONG 方程组
对不可压转捩/湍流平板边界层的计算
陈广业，高
歌，尹幸愉
（北京航空航天大学动力系，北京100083）
摘要：与以往的湍流模型不同的是GAO-YONG 不可压湍流控制方程组不需要任何经验系数及壁面函数。本文运用SIMPLE 方法OUIC 格式求解GAO-YONG 方程组，
对二维零压力梯度下平板转捩/湍流边界层进行了数值模拟。结果表明，方程不仅对于边界层流动的各项细节（如表面摩擦系数!f ，对数律，亏损律等）能作出良好的预测，成功地解决了以往一般模型不能同时计算近壁区和远壁区的难题，而且能够预报层流-湍流转捩过程。本文还对机械能方程如何影响边界层近壁区特性进行了数值研究。关
键
词：航空、航天推进系统；湍流；GAO-YONG 理性方程；二维不可压平板边界层；转捩
中图分类号：O357.5
文献标识码：A
The Calculation of Incompressible Wall Boundary Layer Using GAO-YONG Turbulence Eguations
CHEN Guang-ye ，GAO Ge ，YIN Xing-yu
（Beijing University of Aeronautics and Astronautics ，Beijing 100083，China ）Abstract ：UnIike wideIy ud turbuIence modeIs ，the GAO-YONG eguations u no empiricaI coef-ficients or waII functions.In this paper ，two-dimensionaI ，stationary ，and incompressibIe transitionaI /turbu-Ient boundary Iayer with zero pressure gradient is numericaIIy simuIated with GAO-YONG turbuIence mod-eI.EmpIoying staggered grid arrangement ，an incompressibIe SIMPLE code with a cond-order OUIC scheme is ud.The caIcuIation successfuIIy simuIates the fIow and predicts the fundamentaI properties of the boundary Iayer such as !f ，the Iog Iaw ，and the defect Iaw ，etc ，which are cIo to the experimentaI da-ta.AII the resuIts further verify the adaptabiIity of the eguations to reaI turbuIent fIow.
Key words ：aerospace propuIsion ；turbuIence ；GAO-YONG rationaI modeI ；
incompressibIe waII boundary Iayer ；transition
1前言
流体力学中的湍流问题是一个重要而又古老的问题，一直是流体力学中心问题之一，也被认为是经典
物理中留下的最大难题。这曾经吸引了不少最著名的力学和物理学家参与研究，但至今仍是流体力学中有待进一步研究的最重要的问题之
一。20世纪由于计算机和计算方法的发展，曾经有人认为力学问题几乎都可以用计算的方法解决。但实际上由于湍流理论模式没有得到妥善解决，使得很多问题无法精确计算。在流体工程中，人们遇见的也大多为湍流。湍流对于摩擦阻力、热交换、流动的分离及边界层的增长都有重要影响。因此，对湍流的研究，不仅有理论意义，更有
重大的应用价值。
20世纪90年代，
高歌和熊焰在侧偏平均的基础上，发展出了一套完整的湍流模式理论。由于方程组是建立在对湍流物理的深入理解之上，现有的计算表明，它具有很强的生命力和极好的发展前景。比如该方程在l999年完成了验证三个代表性算例：平板边界层流动、平面射流及后台
阶流动［l ~3］；在2000年又对圆射流进行了模拟，成功的解决了对于一般模型而言的圆射流/平面
射流异常现象（现有的微分方程湍流模型在对二维平面射流进行数值模拟时，计算得到的射流扩散率（spreading rate ）与实验值吻合的较好；而使用
相同的模型计算轴对称圆射流时得到的射流扩散率要超出实验值约40%，且大于二维平面射流的
计算值，这与实验结果相反）［4］。本文的工作主要
阳痿原因是对转捩/湍流平板边界层进行更加深入细致的研究。
平板边界层是形式简单的流动，但其中的物理机制极为复杂，至今未被完全理解，其结果就是即便当今最好的湍流模型，也难于完满的予以解决。转捩一般分为自然转捩和by-pass 转捩，自然转捩发生在来流湍流强度很低的情况下，当超过一定的雷诺数，边界层内小扰动形成二维T-S 波，
再发展成三维波并快速增长形成湍流斑，最后发展为湍流。当来流湍流强度增大，湍流斑在边界层中直接形成，称为by-pass 转捩。本文主要研究自然转捩。湍流边界层中，远离壁面的区域的上层区域，速度梯度很小，而靠近壁面的区域则有极大的速度梯度，造成边界层内层（粘性底层，过渡层，对数律层）和外层（亏损律层）有很大的不同，比如湍流粘性，涡团尺度，平均流动能和湍动能之间的能量转换等等，而描述方式也不一样，内层由壁面律描述：U +=f （y +），U +=U T ，y +=P U T y /I ，外层可用亏损律描述：（U max -U ）/U T =g
（y /S ）［5］
。对于湍流边界层的计算，
涡团粘性的湍流模型，若要把近壁远壁区域同时算好，不仅需要联
合两种模型分别计算不同区域，而且还要调整经
验系数［6］
。而对于多方程模型，或者要进行壁面
修正，或者使用壁面函数，另外还要使用更多的经验系数。本文尝试不使用任何经验系数及壁面函数对转捩/湍流平板边界层进行计算，对Gao -Yong 湍流模型进行验证。
!控制方程
［!，"］
N-S 方程经系综平均而得的平均流连续方
程、平均动量方程、以及运用GAO-YONG 模型经加权侧偏平均而得的漂移流连续方程、漂移流动量方程和利用平均流与漂移流之间相互作用推导出来的机械能方程，分别如下：
1982年什么命·
U =0（l ）
U t +（
U · ） U =-l P P +U L 2 U +l P
· T T （2） ·
U =0（3）
U t +（
U · ） U =-l P
P +U L 2（ U - U ）+l P
（ · T T - · T T ）-（ U · ） U （4）
L
J =l
（R · ）J U · U J ！=R ·-l P P +U L
2
（
U - U ）[+l P
如何喷香水·（ T T - T T ）-（ U · ） U -I U ]
It （5）式中： U ， P 分别为平均流的速度和压力； U ， P 分别为漂移流的速度和压力；P ，
U L 为平均流密度及层流运动粘性系数； T T ， T T 为平均流及漂移流的
湍流粘性应力，其计算方法详见文献［2，3］
，其中严格推导得出的实度系数在二维边界层时为l /8。R 为脉动场的平均漂移位矢，
其模量由（5）式决定，方向为平均拉伸主轴方向。
"计算方法
使用SIMPLE 方法［7］
及交错网格对流场进行
计算。SIMPLE 算法的主要思想是，采用分离式的迭代求法（gregated method ），即先求解U 速度场，
再求解1速度场，
最后利用连续方程使假定的压力场随迭代不断改进，直至收敛。对于本例，则方
程组的求解步骤为：（l ）求解漂移流 U 速度场，再求解漂移流 1速度场，
最后利用漂移流连续方程修正漂移压力场和速度场；（2）求解机械能方程，得出漂移位矢R ，进而得出湍流粘性系数；（3）求解U 速度场，再求解1速度场，最后利用漂移流
连续方程修正压力场和平均速度场；（4）若未收敛，回到（l ）
再行迭代求解。计算域取在长高2m X l m 的矩形区域内，为
了验证网格对模型的影响，本文采用了三种网格即稀疏网格（40X 40）、中等密度网格（70X 50）和
8
9l
较密网格（l20>60），水平x 方向平均划分，y 方向网格使用双曲函数分布，第一层网格的量级大约y +=0.l ~l 。计算使用无量纲化方程组，所取雷诺数为Re =（U >l 0）/!=6.3>l06，
其中U 是来流速度，l 0取单位长度。来流条件为：7u =l ，
7u =0，7p =l ，~u =0，~u =0，~p =0，
边界条件的处理是：平板处使用无滑移条件，出口及上边界速度均外推。动量方程中对流项的差分格式使用二阶
OUICK
［8~ll ］
格式，扩散项及源项使用具有二阶精度的中心差分，离散后生成五对角代数方程组，用交替方向的逐线迭代的TDMA 算法［7］求解。
!
计算结果
本文的主要目的是在不使用经验系数和壁面
函数的情况下，对二维零压力梯度下平板自然转捩和湍流进行模拟。为了验证程序及所用方法和精度，首先计算了层流情况，图l 是计算的层流速度型与Blasius 解的比较，沿流线方向取了7个站点的剖面，由图可以看到，其相似性（图中"=y Ue /（
ux \），u /Ue =g （"））及与理论解的符合性都很好。图2中，将层流摩擦系数C f 的理论值
小规模纳税人企业所得税优惠政策C f =0.664/Re \()x 与计算结果作比较，可以看到
符合的较好。程序在此基础上加入湍流模型，并且计算过程中分别对指数格式（一阶）及
OUICK 格式（二阶）进行计算，发现结果差别不明显，这里均
图l
层流计算速度型
Fig.l
Comparison of the calculated velocity profiles
with the Blasus solution.
采用二阶精度的OUICK 格式下的结果。在计算过程中发现，入口条件的确定是影响计算结果的重要因素。本文采用的方法是：首先
进行层流计算，使用层流计算的结果，营求入口的速度型及摩擦系数都要满足Blasius 解，又因为本
文计算的是自然转捩，其典型的转捩雷诺数［5］
是Re x tr =3>l06，
本文所取入口Re x in =l >l06，出口Re x out =2
>l07。图3所示为由层流到湍流速度型的发展。
图2层流计算摩擦系数分布Fig.2
Comparison of calculated
laminar friction coefficient with Blasus solution
图3转捩过程中速度型的变化Fig.3
Mean velocity profiles in the transitional
zone of flat-plat flow
图4所示为完全湍流发展区内近壁速度型。计算的结果与对数律吻合的相当好。图5列出了
在3套不同网格下得到的靠近壁面的速度型之间的比较。在40>40的网格下计算值略低于对数律，加密网格之后，可以看到速度型与对数律吻合。图6是湍流充分发展区计算所得亏损律与实验值（Claur ，l956年）
的比较。图7画出在三套网格下湍流充分发展区整个边界层内的无量纲速
9
9l
度型并与实验值（Klebanoff ，1955年）
比较。图4充分发展湍流区内层速度型Fig.4
Predictions of the velocity profiles for
the wall region in the full developed
turbulence
zone
图5不同网格下计算所得速度型的比较Fig.5
Predictions of the velocity profiles for the wall region using different grids
对湍流进行数值模拟，最重要的参数是表面摩擦系数C f ，因为它直接关系到流动的阻力，是工程上最关心的内容。图8描述的即是沿流向表面摩擦系数C f =
!w 0.5"U 2
e
的变化情况（三套网格），与经验估计（层流部分C f =0.664Re !x ，湍流部分C f =
0.455
ln
（0.06Re x )
踢毽子技巧）比较可见吻合的很好。图9为沿流向形状因子H 的变化，一般层流的形状因
子为2.6
，而湍流为1.4左右，可以看到计算结果也符合的很好。
图6计算所得亏损律与实验值的比较
Fig.6
Comparison of predictions with the experiment data for the
outer region of flat-plate
boundary layer
图7计算所得外层速度型与实验值的比较Fig.7
Mean velocity profiles in the fully developed turbulence zone of flat-plat flow
!讨论
通常，平均流动能方程及扰动流动能方程都是动量方程乘以速度后得出的，方程本身不独立，而机械能方程（5）描述了平均流和漂移流之间的能量转换关系，是一个完全独立的能量方程。方程的左端是一个泰勒级数展开项，表示平均流因漂移流阻力作用在#距离上丧失的动能，右侧则是漂移流阻力在#距离上所作的功。本文的计算都取三阶，即j =3，以达到较好的求解精度。数值计算表明，取二或三阶时差别并不明显，以壁面的速
山西的旅游景点02
度型为例如图10。在求解该方程的过程中发现，方程右侧最后一项，在近壁区的计算精度方面，有至关重要的作用。如果去掉该项，依然可以算好外层速度型，但内层的摩擦速度u !="!u !}w 的计算值偏
大，使得壁面速度型的对数律的斜率偏小。方程右边加入该项后，就能得到如上图4的好结果。机械能方程中有无该项对近壁计算结果的比较见图11。在靠近壁面附近，剪切非常剧烈提供了边界层70%湍流能量的猝发现象（burst-ing ）
也产生于该区。该区内平均流的特性在很大程度上受扰动流挟制。从能量方面考虑，在平均能和湍动能之间的能量转换传递，在大部分范围内，
是平均流动
图8转捩过程中壁面摩擦系数C f 的变化Fig.8
Computed variation of friction coefficient
for the transitionaI zone of fIat-pIat
fIow
图9
形状因子的变化Fig.9芭蕾舞者
Downstream variation of computed
shape factors
能降低，能量传给湍动能，但实验也发现在近壁区，大量生成的湍动能竟有一部分又回传给了平
均流，使平均流动能增加［5］
。
D "u
D t ·#
可以表示为
返
图10机械能方程阶数对壁面律影响的比较Fig.10
Comparison of the different jth orders of the mechanicaI energy
eguation
图11
机械能方程特性比较
夸奖老师的词语Fig.11
Comparison of the characteristics of the mechanicaI energy eguation
还给平均流的能量，该项在近壁面边界层的精确计算中发挥发挥关键作用。
!结论
对零压力梯度下光滑平板的模拟结果表明，GAO-YONG 模型对于固体边界层强湍流区域和弱湍流区域都有很强的适应能力，对近壁区和远壁
区可以使用统一的模型，而且不需要传统的壁面修正和经验系数。对于自然转捩及充分发展湍流区内的各项参数如对数律、壁面C f 、形状因子H 、内外层速度型等都能做良好的预测。机械能方程
1
02
