湍流模型概述

大多数飞行器都是在高Re数下飞行，表面的流态是湍流。为了准确地确定湍流流态下的摩阻、热流，湍流成为一个重要而困难的研究课题。（一）DNS目前处理湍流数值计算问题有三种方法， 第一种方法即所谓直接数值模拟方法 （DNS方法），直接求解湍流运动的 N-S方程，得到湍流的瞬时流场，即各种尺度的随机运动，可以获得湍流的全部信息。随着现代计算机的发展和先进的数值方法的研究， DNS方法已经成为解决湍流的一种实际的方法。但由于计算机条件的约束，目前只能限于一些低 Re数的简单流动，不能用于工程应用。目前国际上正在做的湍流直接数值模拟还只限于较低的需诺数(Re～200)和非常简单的流动外形， 如平板边界层、完全发展的槽道流， 以及后台阶流动等。用直接数值模拟方法处理工程中的复杂流动问题，即使是当前最先进的计算机也还差三个量级。（二）LES另一种方法称做大涡模拟方法（LES方法）。这是一种折衷的方法，即对湍流脉动部分直接地模拟，将N-S方程在一个小空间域内进行平均（或称之为滤波），以使从流场中去掉小尺度涡，导出大涡所满足的方程。小涡对大涡的影响会出现在大涡方程中，再通过建立模型（亚格子尺度模型）来模拟小涡的影响。由于湍流的大涡结构强烈地依赖于流场的边界形状和边界条件，难以找出普遍的湍流模型来描述具有不同的边界特征的大涡结构，宜做直接模拟。相反地，小尺度涡对边界条件不存在直接依赖关系，而且一般具有各向同性性质。所以亚格子模型具有更大的普适性，比较容易构造，这是它比雷诺平均方法要优越的地方。自从1970年Deardorff第一次给出具有工程意义的LES计算以来，LES方法已经成为计算湍流的最强有力的工具之一，应用的方向也在逐步扩展，但是仍然受计算机条件等的限制，使之成为解决大量工程问题的成熟方法仍有很长的路要走。（三）RANS目前能够用于工程计算的方法就是模式理论。 所谓湍流模式理论， 就是依据湍流的理论知识、实验数据或直接数值模拟结果，对 Reynolds应力做出各种假设，即假设各种经验的和半经验的本构关系，从而使湍流的平均 Reynolds方程封闭。随着计算流体力学的发展，湍流模式理论也有了很大的进步，有了非常丰硕的成果。从对模式处理的出发点不同，可以将湍流模式理论分类成两大类：一类称为二阶矩封闭模式，另一类称涡粘性封闭模式。（1）雷诺应力模式所谓二阶矩封闭模式，是 从Reynolds应力满足的方程出发，将方程右端未知的项（生成项，扩散项，耗散项等）用平均流动的物理量和湍流的特征尺度表示出来。 典型的平均流动的变量是平均速度和平均温度的空间导数。这种模式理论，由于保留了 Reynolds应力所满足的方程，如果模拟的好，可以较好地反映 Reynolds应力随空间和时间的变化规律，因而可以较好地反映湍流运动规律。因此， 二阶矩模式是一种较高级的模式， 但是，由于保留了Reynolds应力的方程，加上平均运动的方程整个方程组总计 15个方程，是一个庞大的方程组，应用这样一个庞大的方程组来解决实际工程问题， 计算量很大，这就极大地限制了二阶矩模式在工程问题中的应用。（2）涡粘性模式在工程湍流问题中得到广泛应用的模式是涡粘性模式。这是由Boussinesq仿照分子粘性的思路提出的，即设Reynolds应力为，221uiujT(Ui,jUj,i3Uk,kij)3kij（）这里kuiuj是湍动能，T称为涡粘性系数，这是最早提出的基准涡粘性模式，即假设2雷诺应力与平均速度应变率成线性关系，当平均速度应变率确定后，六个雷诺应力只需要通过确定一个涡粘性系数T就可完全确定，且涡粘性系数各向同性，可以通过附加的湍流量来模化，比如湍动能k，耗散率，比耗散率以及其它湍流量k/，lk3/2/，qk，根据引入的湍流量的不同，可以得到不同的涡粘性模式，比如常见的k，k-w模式，以及后来不断得到发展的k，q-w，k-l等模式，涡粘性系数可以分别表示为Ck2/，TkCk，Tq2Ckl.TC，TC，T为了使控制方程封闭，引入多少个附加的湍流量，就要同时求解多少个附加的微分方程，根据求解的附加的微分方程的数目，一般可将涡粘性模式划分为三类：零方程模式，半方程模型，一方程模式，两方程模式。1）零方程模式所谓零方程模式是试图直接用平均流动物理量模化T，而不引入任何湍流量（如k,等）。例如，Prandttl的混合长理论就是一种零方程模式：Tl2U（5.7）y式中l称为混合长。在零方程模式的框架下，得到最为广泛应用的是Baldwin-Lomax模式[22]。该模式是对湍流边界层的内层和外层采用不同的混合长假设。这是因为靠近壁面处，湍流脉动受到很大的抑制，含能涡的尺度减小很多，因此长度尺度减小很多；另一方面，在边界层外缘，湍流呈间歇状，质量、动量和能量的输运能力大大下降，即湍流的扩散能力减小。这样，应用混合长理论来确定涡粘性系数在这两个不同的区域应该有不同的形式。Baldwin-Lomax模式的具体数学描述如下。T

(T)innyyc）(T)onty（5.8yc这里yc是(T)inn (T)ont的离壁面最小距离 y值。对于内层，即 y yc，有(T)inn l2 （5.9）是涡量，ijkUk,j,l是长度尺度lky(1exy(y/A))其中k=0.4是Karman常数，A+是模化常数，y是无量纲法向距离：y Uy/ w而u是摩擦速度，其含义为Uu ,yw此处下标w表示壁面。对于外层，即yyc，有(T)out(FwakeFkleb(y)其中Fwakemin(ymaxFmax,CwkymaxUdif2/Fmax)Fmax是下列函数的最大值：

5.10）5.11）F(y)y(1exp(y/A))而ymax是F(y)达到最大值的位置。Fkleb是所谓的Klebanoff间歇函数：5.5(Ckleby)1Fkleb(y)16ymaxUdif是平均速度分布中最大值和最小值之差。几个模化常数的值如下：A 26.0; C 0.02668; Ckleb 0.3; Cwk 1.0,; K 0.4.由上述模化关系中可以看出， Reynolds应力完全地由当时当地的平均流参数用代数关系式所决定。平均流场的任何变化立刻为当地的湍流所感知， 这表明零方程模式是一个平衡态模式，假定湍流运动永远处于和平均运动的平衡之中。 实际上对大多数湍流运动而言，并非如此，特别是对平均流空间和时间有剧烈变化的情形， 再有因为坐标 y显式地出现在湍流模式中，零方程模式不具有张量不变性， 当将它应用到复杂几何外形的流动的数值模拟会带来困难。当流动发生分离时， Baldwin-Lomax 模式会遇到困难，这是因为在分离点和再附点附近，摩擦速度 u为零，此时要引入一些人为的干涉来消除这些困难。计算实践表明，只要流动是附体的，零方程模式一般都可以较好地确定压强分布，但是摩阻和传热率的估算不够准确，特别是当流动有分离和再附时。这是因为附体流压强分布对湍流应力不敏感。总之，对附体流动，如果只关心压强分布，应用零方程模式通常可以给出满意的结果，而且模式应用起来十分简便。但是对于我们计算摩阻的需求，零方程模式是不能满足要求。对于有分离、再附等复杂流动，零方程模式是不适用的。2）半方程模式为了能计算具有较强压强梯度，特别是较强逆压梯度的非平衡湍流边界层,Johnson-King 于1985年提出了一个非平衡代数模型， 该模型仍采用涡粘性假设， 把涡粘性的分布与最大剪切应力联系在一起， 内层涡粘性与外层涡粘性分布用一个指数函数作光滑拟合，外层涡粘性系数作为一个自由参数，由描述最大剪切应力沿流向变化的常微分方程来确定，此常微分方程是由湍流动能方程导出的，故此模型又称为半方程模型。JK模型虽然仍采用涡粘性假设， 却包含有雷诺应力模型的特点。 由于求解常微分方程比一方程，二方程模型中求解偏微分方程要简单，省时的多，故用 JK模型的工作量只略高于通常平衡状态的零方程代数模型的工作量JK模型后又经不断修正，发展了 JK1990A，JK1990J以及JK1992等改进型3）一方程模式Baldwin-Barth(BB) 模型是在二方程模型中， 将某一导出的应变量作为基本物理量而得到的，应用此一方程模型可避免求解两方程时会遇到的某些数值困难。 BB一方程模型所选择的导出应变量为“湍流雷诺数” Rt。BB模型对计算网格的要求低，壁面的网格可以与采用 BL代数模型的相当，而不象两方程 k-e模型那样要求壁面网格很细，这样就避免了在k-e模型中流场求解的刚性问题。Spalart-Allmaras(SA) 模型与BB模型不同，不是直接利用 k-e模型两方程模型加于简化而得，而是从经验和量纲分析出发， 由针对简单流动在逐渐补充发展而适用于带有层流流动的固壁湍流流动的一方程模型，模型中选用的应变量是与涡粘性

T相关的量

~，除在粘性次层外，

~与

T是相等的。上述两种一方程模型具有相似的特点， 它们不象代数模型那样需要分为内层模型，层模型或壁面模型，尾流模型，同时亦不需要沿法向网格寻找最大值，因此易于用到非结构网格中去；但由于在每个时间步长内，需要对整个流场求解一组偏微分方程，故比和JK模型更费机时4）两方程模式

外BL2.1k- 两方程模式 标准k- 两方程模式k- 模式是最为人所知和应用最广泛的两方程涡粘性模式，为积分到壁面的不可压缩可压缩湍流的两方程涡粘性模式，各种不同版本的 k- 模式常见于各种文献中，选择

/Jones-Launderk-度尺度

模式最初的发展是为了改善混合长 (mixing-length)(turbulent length scale)的代数表示(algebraic

模式和避免复杂流动中湍流长prescription) 。它求解两个湍流标量

k和

的输运方程。

k方程表示湍动能输运方程，

方程表示湍动能的耗散率。该模式对较小压力梯度 (relativelysmallpressuregradients) 下的自由剪切流(free-shear-layer flows) 具有较好的结果。对于壁面流动 (wall boundedflows) ，在零或者小平均压力梯度下， 模式结果和实验结果符合得较为一致， 但是对大的逆压梯度 (adversepressuregradients) ，其结果就不太正确了。另外，在壁面附近，该模式需要壁面衰减函数(wall-dampingfunctions) 和较好的网格分布。模式方程雷诺应力的涡粘性模型为tijuuj2t(SijSnnij3)2kij3这里t为涡粘性(eddyviscosity)，Sij为平均速度应变率张量(mean-velocitystrain-ratetensor)，为流体密度，k为湍动能，ij为克罗内克算子(Kroneckerdelta)。涡粘性定义为湍动能k和湍流耗散率的函数tcfk2基于量纲分析，涡粘性由流体密度，湍流速度尺度(turbulentvelocityscale)k2和长度尺度(length-scale)k32来标度，衰减函数f由湍流雷诺数Retk2来模化。湍流输运方程可表示成以下形式湍流能量输运方程kujkktijSijφktxjxjxjk能量耗散输运方程2tujxjc1tijSijc2f2φxjkk这里右端项分别表示生成项(productionterm)耗散项(dissipationterm)和壁面项(wallterm)。模式常数和参数模式中各常数的定义为c0.09c11.45c21.92k1.01.3Prt0.9近壁衰减函数fexp(3.4(10.02Ret)2)k2f210.3exp(Ret2)和Ret壁面项22φkk和φ2t2us2y2y这里us为平行于壁面的流动速度。边界条件积分到壁面的无滑移边界条件为k002.1.2可实现性k模式上述标准k模式，对于高平均切变率流动会出现非物理的结果（例如当Sk/3.7时，其中S2SijSij）。为了保证模式的可实现性，模式函数C不应该是常数，而应当是平均庆变率的函数。实验表明，对边界层流动和均匀切变流，C的值是非常不同的。为此人们根据可实现性对模式的约束条件，建议采用以下形式的C（Reynolds,1987,Shih,1994）1（5.19）CA0AsU*k式中U*SijSij**ijij*2ijkijijkijijijkk而ij是在以角速度k旋转的旋转坐标系中得到的平均旋转速率。As6cos,1cos1(6w)3）SijSjkSki~（5.20W,SijSij.~3SS上述关系式中唯一未确定的系数是。为简单起见，可以设其为常数。对边界层流动。A可以取0=4.0。对其他流动，0的数值可以调节。AA2.1.3低Reynolds数k上述几种k模式适用于高Reynolds数情形。但是对近壁区，湍流需诺数很低，对湍流动力学而言粘性效应非常重要，此时湍流Reynolds数的效应必须加以考虑。我们研究摩阻的计算关注的恰恰是近壁区，因此低Reynolds数k模式的研究是十分重要的。现将有关结果整理如下：低Reynolds数下的涡粘性和k模式方程为cfk(k)（5.22）T(k),t(Uik),i(T)k,juiujUi,j(5.23)k,jT),j2TS,jS,j(),t(Ui),i(c1f1SC2f2C3,jk式中S2SijSij,Sij1(Ui,jUj,i)2uiuj2kijT(Ui,jUj,i2Uk,kij)33所有模化常数如下：1C4 AsU*kC1C2kU*Sij*ff1f2

max0.43,51.9,C31.01.0,1.2Sij*Sij*ijij,Sk1Sij Uk,kij31expa1Ra2R2a3R3a4R4a5R51expa1Ra2R2a3R3a4R4a5R510.22expRt236其中13k2(k)2k2R,Rt此处f和f1,f2称为阻尼函数，是一个经验公式用来反映近壁区低雷诺数效应，系数ai和ai列表如下：I12345ai-3-6.0×10-5-7-3.6-9-123.3×106.6×10×108.4×10ai-3-5.7×10-5-7-3.6-9-122.53×106.55×10×108.3×10 常见

k-

两方程模式在文献中有许多种

k

涡粘性模式。为了便于比较，我们将几种常见的

k

模式作一归类。它们的主要区别一是在

k和 的方程及其边界条件，另一方面是阻尼函数

f

的取法。模 式代号

作 者ChLBNTLSMKYSS&LCMOTT

Chien,1982LamandBremhorst,1981NaganoandTagawa,1990LaunderandSharma,1974MyongandKasagi,1988YangandShih,1991ShihandLumley,1993ZhuandShih,1995所有上述八种模式都可以用一个统一的方程组表示：tCfk(5)dktkUiD(6)dtxitijxjkxidtC1f11UiC2f2E(7)dtxitijxjxi有关的项，D，E，T列表如下：ModelTDEChk2k2exp(0.5y)y2y2LBk0NTk02k k

002U

2LS2yMKk0YSk0

2 T y202U 2y2S&Lk0

2U

2y2CMOTTk0

2U

2阻尼函数 f,f1和f2对不同的模式有不同的表示式。

y2ModelCh

ff1f21exp(0.115y)11.22exp(Rt2)3620.5).053LB(1e0.165Rk)2(111exp(Rt2)Rtf[1.3exp((Rt)2)]24.16.5NTy11exp1[1exp(y)]226Rt3/46LSexp3.411.3exp(Rt2)(1Rt/50)23.45y[12exp(Rt2)]MK(1exp(1936)(1))Rt70[1exp(y)]251exp(.004y5e5y2YS12e6y38e4)8y1exp(.004Rk5e5Rk2S&L12e6Rk38e4)8Rk1exp(.004Rk5e5Rk2CMOTT12e6Rk348e8Rk)式中Rk,y和Rt定义为Rkky,yuy,Rtk2.模式中出现的模化常数分别为ModelCC1C2Ch.091.351.80LB.091.441.92NT.091.451.90LS.091.441.92MK.091.401.80YS.091.441.92S&L.091.441.92CMOTTEq(5.19)1.441.92对不同的模式有不同的处理连界条件的方法：

1.22exp(Rt2)361.22exp(Rt2)361.22exp(Rt2)36k1.01.31.01.31.41.31.01.31.41.31.01.31.01.31.01.3Model B.C.forkw B.C.for wCh00LB2k02yNT0LS0MK0YS0S&L0.25u2CMOTT0.25u22.2 其它双方程模式涡粘性系数的量纲为速度×长度，当用T Ck2/ 。我们注意到，对标准 k

2ky202ky22k2yu40.251u40.251k, 来模化时，它们之间的关系为模式的 方程，在固壁上有奇点问题（壁面上湍动能k 0），这是因为模式不尽合理带来的非物理的奇点。此外在计算中由于 k, 在壁面附近变化剧烈， 必须在物面附近将网格划分得非常小， 才能得到合理的结果。 为了克服这些困难，人们试图寻找其它的湍流量来代替 k, 。可能的选择有/k,k/,lk3/2/,qk,相应地，涡粘性系数可表示成：TCk,TCk,TCq2,TCkl.现在就来介绍几种典型的模式：2.2.1k-w两方程模式(Wilcox)k-模式是最为人所知和应用最广泛的两方程涡粘性模式，为积分到壁面的不可压缩/可压缩湍流的两方程涡粘性模式，最主要文献来自Wilcox。求解湍动能k和它的/k,k/,lk3/2/,qk,(specificdissipationrate) 的对流输运方程已经证明 Wilcoxk- 模式在粘性子层比 k- 具有更好的数值稳定性。由于壁面附近，值较大，模式不象 k- 模式或者其它两方程模式，它不需要显式的壁面衰减函数。对于比较缓的逆压梯度流动，该模式在对数区域给出的结果和实验数据较为符合。模式方程雷诺应力的涡粘性模型为tij 2t(Sij Snn ij 3) 2k ij 3这里 t为涡粘性(eddy viscosity) ，Sij为平均速度应变率张量(mean-velocitystrain-ratetensor), 为流体密度 , k为湍动能 , ij 为克罗内克算子 (Kroneckerdelta) 。涡粘性定义为湍动能 k和比耗散率 的函数t

kk和 的输运方程为kujk(*ktijSij-*ktt)xjxjuj(t)tijSij-2txjxjk模式常数和参数模式中各常数的定义为53*99401000.5*0.5Prt0.9边界条件对边界层流动，壁面无滑移边界条件为k0和6102(y1)这里y1为离开壁面第一个点的距离，且+。y1<1对称边界条件采用零梯度条件，各种附加的边界条件将在具体流动中讨论。 两方程模式(Menter)k- SST剪切应力输运(shear-stress-transport) 模式在近壁处采用 Wilcoxk- 模式，在边界层边缘 boundarylayeredges) 和自由剪切层 free-shearlayers) 采用k-模式 k- 形式 ，其间通过一个混合函数 blending function 来过渡，属于积分到壁面的不可压缩/可压缩湍流的两方程涡粘性模式。为了有效结合 k- 和k- 模式，统一写成 k- 形式模式方程涡粘性定义为a1kTmaxa1;F2这里是涡量的绝对值，a0.31，F是混合函数。122F2tanhmax2ky,5000.99y2T的形式解决了湍流剪切应力在逆压梯度边界层的输运。 k和 由相应的模式输运方程得到。湍动能输运方程kujk(kt)kkttijSijxjxj湍流比耗散率方程uj(t)P22(1F1)2ktxjxjxjxj上式中最后一项代表交错扩散项(cross-diffusionterm)，生成项P2SijSnnij3Sij2模式常数和参数2F1tanhminmaxk,50042ky,y20.99y2CDk这里CDk22k,1020maxxjxj这里CDk代表k-模式中的交叉扩散(cross-diffusion)。SST模式常数a10.31*0.090.41模式参数,,k,由来表示，用1,2分别表示原始k-模式系数和转化的k-模式系数F11112这里k,,,Innermodel 系数:k10.8510.51*2*111Outermodel 系数:k2 1.0 2 0.856 2

0.0750.5530.0828*2*0.4402222.2.3k模式方程为DkTk,iuiujUi,jk(5.29)Dtk,iDTT2C1kuiujUi,jC2(5.30)Dt,i,i,i,i模化常数为C1C2k22e20.92,k0.4,T1.5.C1,T对低雷诺数流动有TCfk(5.31)kf1exp(a1Ra2R2a3R3)(5.32)其中Rk(uup)C21.92[10.22exp(Rt2/36)]1在k模式的框架下，Speziale(1990)[33]提出了下列的模式：kTCfk,(5.33)DkTk,iuiujUi,jk(5.34)Dtk,iDT2TDt,i,i,i2,i2(5.35)2Tk,i,i(C11)uiujUi,j(C2f21)k1k模化常数及参数分列如下：Mod