中科院计算流体力学讲义CFD第11讲湍流及转捩1

1、1计算流体力学讲义计算流体力学讲义2011 第十一讲第十一讲 湍流与转捩湍流与转捩 （1） ；力学所主楼；力学所主楼219； 82543801 知识点：知识点： 讲义、课件上传至讲义、课件上传至 (流体中文网）流体中文网） - “流体论坛流体论坛” -“ CFD基础理论基础理论 ”讲课录像及讲义上传至网盘讲课录像及讲义上传至网盘 线性稳定性理论线性稳定性理论 转捩的预测方法转捩的预测方法 壁湍流转捩的涡动力学机制壁湍流转捩的涡动力学机制2 11.1 线性稳定性理论线性稳定性理论一、一、 稳定性基本概念稳定性基本概念常识：流体中的不稳定性常识：流体中的不稳定性K-H不稳定性不稳定性A. K-H

2、(Kelvin-Helmholtz）不稳定性）不稳定性 自由剪切流自由剪切流的的无粘无粘不稳定性不稳定性混合层混合层 K-H不稳定性不稳定性K-H不稳定性的关键：不稳定性的关键：速度剖面有拐点速度剖面有拐点 Lee-Lin: 速度剖面的拐点是无粘不稳定性的必要条件速度剖面的拐点是无粘不稳定性的必要条件流体不禁搓，流体不禁搓，一搓搓出涡一搓搓出涡已知某运动状态；已知某运动状态；在此基础上施加微小扰动；在此基础上施加微小扰动；如扰动随时间（或空间）衰减，则称系统稳定，否则为不稳定如扰动随时间（或空间）衰减，则称系统稳定，否则为不稳定3自然界中自然界中 K-H不稳定性图片不稳定性图片智利塞尔扣克岛智

3、利塞尔扣克岛的卡门涡街的卡门涡街澳大利亚澳大利亚Duval山上空的云山上空的云 KelvinHelmholtz instability clouds in San Francisco佛兰格尔岛周围的卡佛兰格尔岛周围的卡门涡街门涡街高速流低速流自由剪切层受到扰动界面变形后的情况自由剪切层受到扰动界面变形后的情况 K-H不稳定性的产生机理不稳定性的产生机理受阻减速，压力升高，受阻减速，压力升高，产生高压区产生高压区高压导致变高压导致变形加剧形加剧4B. T-S (Tollmien-Schlichting) 不稳定性不稳定性不可压不可压 壁面剪切流壁面剪切流的的粘性粘性不稳定性不稳定性Mack 不稳

4、定性不稳定性 超声速壁面剪切流的不稳定性超声速壁面剪切流的不稳定性不可压边界层速度剖面不可压边界层速度剖面 （Blasius解）解） 无拐点无拐点可压缩情况可压缩情况 Mach数足够高时会出现广义拐数足够高时会出现广义拐点点 出现无粘不稳定性出现无粘不稳定性 y/00.511.500.5z=400du/dydegreedegreedegreeGIPs0)(dydudyd不可压缩无粘不可压缩无粘不稳定性不稳定性 需存在拐点需存在拐点可压缩无粘可压缩无粘不稳定性不稳定性需存在广义拐点需存在广义拐点022dyudMach 6 钝锥（钝锥（1攻角）攻角）不同子午面不同子午面 的分布的分布dydu超音速

5、平板边界超音速平板边界层的不稳定波层的不稳定波第第1模态（模态（T-S波）波）第第2模态模态 （Mack模态）模态）5激波激波密度界面密度界面R-M (Richtmyer-Meshkov)不稳定性不稳定性 激波与密度界面作用的激波与密度界面作用的斜压斜压效应效应惯性约束聚变（惯性约束聚变（ICF）示意图示意图小知识小知识 涡的产生机制：涡的产生机制： 粘性、粘性、 斜压、有旋的外力斜压、有旋的外力0p激波密度界面p)(31)(1)()(2vvFvvpdtd斜压项斜压项6D. R-T (Reyleigh-Taylor)不稳定性不稳定性 重力带来的不稳定性重力带来的不稳定性R-T (Reyleig

6、h-Taylor)不稳定性不稳定性重重介介质质轻轻介介质质7E Barnard热对流不稳定性热对流不稳定性其他学科的不稳定性：其他学科的不稳定性：Euler压杆压杆的不稳定性的不稳定性Barnard 热对流热对流的胞格结构的胞格结构板壳的不板壳的不稳定性稳定性8二、二、 稳定性问题的常用数学方法稳定性问题的常用数学方法 线性稳定性分析线性稳定性分析Step 1: 得到线性化的扰动方程得到线性化的扰动方程0PU控制方程为：控制方程为：已知其具有解已知其具有解 0U0P0U最好是精确解，也可最好是精确解，也可用高精度的数值解用高精度的数值解令令：UUU00)(UUU0PP舍弃高阶小量，得到线性化的

7、舍弃高阶小量，得到线性化的扰动方程扰动方程0UL（1）xuuxuuxuuxuuuuxuuxuu) () (000000例如：例如： 平板的平板的Blasius解，槽道的解，槽道的Poiseuille 解解线性方程线性方程9Step 2: 求解求解 的特征值问题的特征值问题什么条件下具有非零解，非零解如何？什么条件下具有非零解，非零解如何？通常假设在某些方向具有周期性，转化为一维通常假设在某些方向具有周期性，转化为一维问题问题)()(tzxieyUU0UL数值方法：数值方法： 将将 （1） 离散离散代数方程代数方程何时有非零解，何时有非零解， 非零解如何？非零解如何？ 特征值问题特征值问题xAx

8、0)(xIA0Bx什么条件下有什么条件下有非零解？非零解？特征值问题计算量巨大，目特征值问题计算量巨大，目前通常只能处理一维问题前通常只能处理一维问题10三、三、 稳定性问题示例稳定性问题示例 不可压缩槽道流动不可压缩槽道流动的线性稳定性（的线性稳定性（LST）理论）理论 (以二维为例以二维为例)uuuuu2Re10pt uuuStep 1: 获得线性化扰动方程获得线性化扰动方程 令令：Re2, 0,12xpvyuPoiseuille解：解： ppp(2)代入方程（代入方程（2），并舍去高阶小量得到线性化的），并舍去高阶小量得到线性化的扰动方程扰动方程xuuxuuxuuxuuuuxuuxuu)

9、 () (uuuuuuu 2Re10pt（3）1） 控制方程及边界条件控制方程及边界条件11研究扰动发展的空间模式和时间模式扰动源)()( )( )( txieypyvyupvu空间模式：空间模式： 任一点的扰动具有时间周期性任一点的扰动具有时间周期性 符合物理条件符合物理条件 假设扰动具有如下形式：假设扰动具有如下形式：沿流向及时间方向具有波动特性沿流向及时间方向具有波动特性称为称为Tolmien-Schlichting（T-S）波）波任意扰动可分解为正弦波的叠加任意扰动可分解为正弦波的叠加 线性系统各成分无线性系统各成分无相互作用相互作用 可独立研究可独立研究为实数为实数iri为复数为复数

10、 扰动波的振幅沿流向指数变化扰动波的振幅沿流向指数变化xieAxA)0(/ )(空间增长率中性0扰动衰减0扰动增长0i时间模式：时间模式： 扰动具有流向的周期性扰动具有流向的周期性 假设一窗口沿流向运动，研究窗口内扰动的演化假设一窗口沿流向运动，研究窗口内扰动的演化为实数为实数为复数为复数iri 扰动波的振幅虽时间变化扰动波的振幅虽时间变化tieAtA)0(/ )(时间增长率中性扰动衰减扰动增长000i12以时间模式为例：以时间模式为例：)(Re1)(Re1022222222yvxvypxvutvyuxuxpyuvxuutuyvxu)()( )( )( txieypyvyupvu0yvuiuy

11、pivyuuuii)(Re1 222)()()(,txitxitxieuitueyuyueuixuvyypvuii )(Re1 222（4）（5）（6）线性偏微方程（线性偏微方程（3）转化成为含参数的线性常微方程组（）转化成为含参数的线性常微方程组（4）-（6）谱方法的常谱方法的常规做法规做法通过消元法，转化为更高阶的常微方程通过消元法，转化为更高阶的常微方程 （不是必须的）（不是必须的）常用做法，通常还可以常用做法，通常还可以反向为之：反向为之： 高阶方程转高阶方程转化为低阶方程组化为低阶方程组)4()5(iy消去p yviu1vyuyuivyRe222222222Orr-Sommerfel

12、d（O-S) 方程方程22224422222yyy其中：13最终，控制方程为最终，控制方程为O-S方程：方程：vyuyuivyRe222222222边界条件：边界条件：1y0y1y0 v0yvui0 u0yvy=1 （固壁）:y=0 （中心线，对称）（中心线，对称）：033yvyv可以取计算域可以取计算域-1,1,使用固壁边界条件；使用固壁边界条件；也可以取计算域也可以取计算域-1,0,使用固壁及对称使用固壁及对称边界条件边界条件流函数形式的流函数形式的O-S方程方程xvyu,引入流函数，使得：引入流函数，使得：计算出计算出 后，利用公式后，利用公式0yvuiuypivyuuuii)(Re1

13、222计算其他两个量计算其他两个量v 则：则：)()( txiey令：iv 常数倍满足的方程及边界条件与满足的方程及边界条件与 完全相同完全相同。v 14如果 恒大于（或恒小于0），则必有 小知识：小知识： 关于关于O-S方程方程vyuyuivyRe2222222221） O-S方程适用于不可压方程适用于不可压平行流平行流的稳定性问题的稳定性问题 （不仅槽道流）（不仅槽道流）2） 准平行流准平行流 （流线沿（流线沿x方向接近平行）也可使用（例如边界层流动）方向接近平行）也可使用（例如边界层流动）3）如果舍去粘性（左端）项，则方程称为）如果舍去粘性（左端）项，则方程称为Rayleigh方程方程0

14、Re1Rayleigh拐点定理：拐点定理： Rayleigh方程存在不稳定解的必要条件是速度型存在拐点。即存在某点方程存在不稳定解的必要条件是速度型存在拐点。即存在某点 使得使得sy022syyu若存在无粘不稳定性，该项必有0点。022222vyuyu01122222*dyvyuyuv分部积分，并取虚部，得：01122 dyucuvci/cu 0/iic01122 dyucuv不存在非稳定解152） O-S方程的解法方程的解法 数学表述数学表述 奇性（特征值）问题：奇性（特征值）问题： 参数参数 为何值时，方程有非零解？为何值时，方程有非零解？非零解如何？非零解如何？vyuyuivyRe222

15、222222,Re,0),(Re,F0i时间发展槽道湍流：时间发展槽道湍流： (通常通常) 给定给定Re及及 ，问，问 取何值时，取何值时，O-SO-S方方程有非零解？程有非零解？iri增长率求解步骤：求解步骤： 1） 将将O-S方程离散，得到线性代数方程组方程离散，得到线性代数方程组 离散方法：离散方法： 差分法、有限元法、谱方法、差分法、有限元法、谱方法、打靶法打靶法 2） 求求 ，使得该方程有非零解（奇性或特征，使得该方程有非零解（奇性或特征值问题）值问题）. 0)(xATnvvvx),. ,(210)(A求出求出局部法：只求出一个局部法：只求出一个 全局法：计算出全部的全局法：计算出全

16、部的 16试计算试计算 的时间发展槽流中的时间发展槽流中 （即波长（即波长2p p）T-S波的频率及增长率波的频率及增长率四四: 例题例题7500Re 1p2xL2yL21yuvyuyuivyRe222222222Step 1: 离散离散 （差分法）（差分法）一维问题网格：网格： 均匀网格均匀网格 （简单，但需要较多网格点（简单，但需要较多网格点 N 300） 非均匀网格非均匀网格)(yy 差分离散：差分离散：3211244321123321122112/ )464 (2/ )22 (/ )2 (/ ) (hvvvvvyvhvvvvyvhvvvyvhvvyvjjjjjjjjjjjjjjjjjj

17、2阶格式4阶格式：自行推导阶格式：自行推导 （利用小程序）（利用小程序）yvyv非均匀网格：非均匀网格：22224422222yyy17问题：问题： 会产生大量非物理解会产生大量非物理解 （例：（例： 1000个网格点算个网格点算出出1000个特征值）个特征值） ； 可通过不同网格的对比，可通过不同网格的对比，进行筛选进行筛选离散化后得离散化后得： 0ReRe222222222222vyivyuyuiy0)(xBAvyuyuivyRe222222222Tnvvv),. ,(21xStep 2: 求广义解特征值问题（求广义解特征值问题（1） 即即 为何值时（为何值时（1）有非零解）有非零解 (1

18、)方法方法1： 全局法全局法 一次计算出全部特征值一次计算出全部特征值常用常用Q-Z分解法分解法;其他：其他： 幂法、反幂法、幂法、反幂法、Jacobi法，法，Householder )0)(xIA狭义特征值问题0)(xBA广义特征值问题求解特征值是计求解特征值是计算数学的主要研算数学的主要研究方向，有大量究方向，有大量成熟的方法成熟的方法可借助软件包或可借助软件包或Lapack库等库等 （自行到网上搜索）（自行到网上搜索）通常关心最不稳定的扰通常关心最不稳定的扰动波动波i最大的那个最大的那个18方法方法2： 局部法局部法0)(xBABF A)(令：0)(F方法：方法： Newton法，法，

19、弦位法，抛物线（弦位法，抛物线（Muler）法）法Newton 法：n1n)(F)( / )(1nnnnFF弦位法：弦位法：n1n)(F1n)(nF)(1nF)(nF1n)()()(111nnnnnnnFFF差分化抛物线法：抛物线法：n1n)(F1n)(nF)(1nF已知：已知：nnn,12可连成一条抛物线可连成一条抛物线CBAf202CBA令：求出新的值1n消元法计算行列式（利用5对角特征）计算出计算出 后，求解方程（后，求解方程（1）就可得到特征向量）就可得到特征向量。（1）方程有奇性，可补充一个条件，例如给定某个点的值方程有奇性，可补充一个条件，例如给定某个点的值19 效果更好的方法效果

20、更好的方法 Malik提出的紧致差分格式求解提出的紧致差分格式求解 参考文献参考文献: 周恒等：周恒等： 流动稳定性流动稳定性 （p. 10-13），国防工业出版社），国防工业出版社)(12)(21211jjjjjjyy 22,dyddyd 非均匀网格的超紧致格非均匀网格的超紧致格式（式（4阶精度）阶精度）优点：优点： 紧致紧致网格基仅网格基仅2点点 精度高精度高 4阶阶 直接适用于非均匀网格直接适用于非均匀网格无需坐标变换无需坐标变换0dyvduiudydpivdyuduuii)(Re1 222vdyddypdvuii )(Re1 222原方程组：原方程组：形式变换形式变换变为一阶方程组变为

21、一阶方程组 upvuupvudydAdyudu dyuddyud22通常的通常的做法做法uidyudidyvd22000Re000000000000ReReRe0Re000010000ReReReRe0Re000010002222iiidyuduiiuiiidyuduiiiuiiiA推导仓促推导仓促可能有误可能有误请仔细推导请仔细推导20 与与y无关无关令：AdydBAAAA2)(22dyddyddyddyd)()(12)()(21211jjjjjjyyBBAA upvu21yu,ReReRe0ReRe0Re00000ReReRe2232222uiiiudyudidyudiiidyuduiiA

22、00ReRe00000000000222dyuddyudidyudidydA01jjjjDC带入差分格式：)(12)(21211jjjjjjyy 得：jjjyyBAIC1222121122jjjyyBAID21NX.21令：0)(XGNNDCDCDCDCG.332211边界处表达式边界处表达式 (C1,D1,CN,DN)根据边界条件根据边界条件而定而定内点的表达式内点的表达式特点：特点： 离散形成的代数方程组呈离散形成的代数方程组呈 块两对角块两对角特征特征求行列式，特征值都非常便利！求行列式，特征值都非常便利！其余步骤与前文相同其余步骤与前文相同 可用矩阵广义特征值理论计算全部模态可用矩阵广

23、义特征值理论计算全部模态 也可利用也可利用 计算单个模态计算单个模态0)(G计算域可以取为计算域可以取为-1,1。也可取为也可取为-1,0， 在中心在中心线给对称（或反对称）边线给对称（或反对称）边界条件界条件221jNj 边界条件边界条件1） 如果取完整计算域如果取完整计算域 -1,110yatvu2) 如果取一半计算域如果取一半计算域-1,00y处可设定对称或反对称边界条件处可设定对称或反对称边界条件0反对对称0对称：2233dyvdvdyvddyvd如设定对称条件，只能计算出对称扰动模态如设定对称条件，只能计算出对称扰动模态如设定反对称条件，只能计算出反对称模态如设定反对称条件，只能计算

24、出反对称模态对于槽道流，最不稳定模态是对称模态对于槽道流，最不稳定模态是对称模态但有些情况下，稳定模态在转捩过程中也发挥作用（例但有些情况下，稳定模态在转捩过程中也发挥作用（例如感受性过程，见如感受性过程，见Zhong et al. JFM 556,55-103,2006）01jjjjDC02111DC upvu*001000011C按照内点按照内点方法计算方法计算*000000001D先根据处理内点的方法，求出所有点上的先根据处理内点的方法，求出所有点上的C,D值，再对边界点进行特殊处理值，再对边界点进行特殊处理（D的前两行设为的前两行设为0， C的前两行设为的前两行设为(1,0,0,0)及

25、（及（0，1，0，0）23NNDCDCDCDCG.332211具体解法（局部法）具体解法（局部法）0)()(GF求解求解显然显然NCCCF21)( G因此只需计算每个因此只需计算每个4*4矩阵的行列式即可矩阵的行列式即可（可直接写出表达式，也可用消元法计算）（可直接写出表达式，也可用消元法计算）编制好计算编制好计算 的子程序后，可利用的子程序后，可利用Newton法，弦位法，抛物线法法，弦位法，抛物线法等求解等求解 ，得到复增长率，得到复增长率 G)(F0)(F)( / )(1nnnnFF)()()(111nnnnnnnFFF02CBA1n24计算出计算出 后，利用消元法求解方程后，利用消元法

26、求解方程 ，即，即可得到特征向量可得到特征向量0)(XGXNNDCDCDCDCG.332211消元过程中，充分利用消元过程中，充分利用两对角块矩阵的性质，两对角块矩阵的性质，可减少计算量可减少计算量独立消元注：注： 由于方程有奇性，消元后，最后一个方程为由于方程有奇性，消元后，最后一个方程为0=0, 舍弃最后一个方程，并令最后一个未知数为舍弃最后一个方程，并令最后一个未知数为1 （ ），即可解出），即可解出 显然显然 （a为任意常数）也是原方程的解，因为任意常数）也是原方程的解，因此此 可用某一值（例如可用某一值（例如 ）归一化。）归一化。NX.21X1Nu upvu最终，得到最终，得到 条件

27、下，波数为条件下，波数为 的最不稳定扰动波：的最不稳定扰动波：7500Re 1)()( )( )( txieypyvyupvuXa)max(juy-1-0.500.51-1-0.500.51uruiy-1-0.500.51-1-0.500.51vrviy-1-0.500.51-1-0.500.51prpi扰动波型函数的分布扰动波型函数的分布25 11.2 边界层转捩的预测方法边界层转捩的预测方法1. 经验公式法经验公式法),(ReReMafxx转捩位置z200400600800100000.0010.002transition flowlaminar flowBlasiusKarmanblow

28、 and suctionperturbation regionTransition onset ofHorvaths experment(a)CfMach 6 钝锥边界层表面的摩擦系数分布钝锥边界层表面的摩擦系数分布(Li et al. Phys. Fluid. 22， 025105， 2010. Li et al. AIAA J. 46（11），），2899-2913，2008)x 摩阻或摩阻或热流热流转捩起始点转捩起始点（transition onset)转捩峰（转捩峰（transition peak)充分发展湍流区充分发展湍流区2004015. 02325. 037. 510ReeeMM

29、trL球锥的转捩Reynolds数eM边界层外缘的Mach数3048. 0/197. 0240ReeMtrLe动量厚度定义的转捩Reynolds数262. eN 方法方法)()()(txixtxirie fee ffxLST理论：理论：积分起积分起始点始点x),(xi0),(xi0 xxixdxexAxA0),(0)(/ )( 扰动波进入中性曲线后，开始增长，扰动波进入中性曲线后，开始增长，局部增长率为局部增长率为),(xi0 xxxixdxexAxA0),(0)(/ )(eN理论：扰动波增长到理论：扰动波增长到eN倍，即发生转捩倍，即发生转捩N值需要由实验（或经验）确定，通常为值需要由实验（

30、或经验）确定，通常为811, 即扰动增长即扰动增长4个量级（个量级（10000倍）左右倍）左右。 在不可压缩及航空领域（亚、跨、超）较在不可压缩及航空领域（亚、跨、超）较为成熟。为成熟。 在航天领域（高超声速），还有待检验。在航天领域（高超声速），还有待检验。不足之处：不足之处： 未考虑扰动波进入中性曲线前的衰减过程，没考虑感受性过程。未考虑扰动波进入中性曲线前的衰减过程，没考虑感受性过程。他人的改进：他人的改进： 苏彩虹，周恒等苏彩虹，周恒等 考虑衰减过程考虑衰减过程 C. H. Su, and H. Zhou, Science in China G, 52 (1):115-123 (200

31、9). 273. PSE （抛物化扰动方程）法（抛物化扰动方程）法 优点：优点： 1） 无需平行流假设；无需平行流假设； 2） 可处理非可处理非线性线性(N-PSE)xyStep 1: 得到扰动的控制方程得到扰动的控制方程qqq已知解已知解0 qN线性化线性化0 qLL-PSEN-PSEStep 2: 假设扰动具有波动形式假设扰动具有波动形式),(),( ),(txieyxqtyxq振幅，沿振幅，沿x方向是个方向是个缓变量缓变量 （相对（相对y方向而言）方向而言）Step 3: 带入扰动方程，得到振幅带入扰动方程，得到振幅 的控制方程的控制方程“缓变量缓变量”是个是个很有用的概念，很有用的概念

32、，可用来简化方程可用来简化方程Plantdl边界层理论就边界层理论就是利用是利用“缓变量缓变量”的的概念进行简化的。概念进行简化的。pvuqq LST的的 方程是一维的方程是一维的PSE的的 方程是二维的方程是二维的q q Step 4: 利用缓变量性质，舍弃方程中的椭圆项利用缓变量性质，舍弃方程中的椭圆项（为高阶小量），得到抛物化的扰动方程（为高阶小量），得到抛物化的扰动方程22xqqNxq沿沿x方向推进求解方向推进求解 （类似时间方向的处理），计算量（类似时间方向的处理），计算量相当于一维问题。相当于一维问题。 （“抛物化抛物化”的优势）的优势）非线性项的处理方非线性项的处理方法与谱方法相

33、似法与谱方法相似284. 转捩模型法转捩模型法 （间歇因子模型）（间歇因子模型）实际粘性系数tlc)1 ( 层流粘性系数 湍流粘性系数（由湍流模型给定）湍流间歇因子湍流间歇因子 （0表示纯层表示纯层流，流，1表示纯湍流）表示纯湍流）方法方法 1） 根据经验公式，给定根据经验公式，给定 沿流向的分布沿流向的分布方法方法2） 给出给出 的发展方程，进行求解的发展方程，进行求解 xy29常识：常识： 湍流的间歇性湍流的间歇性外间歇性：外间歇性： 层流及湍流交替出现的现象层流及湍流交替出现的现象Mach 6 钝锥边界层的密度分布钝锥边界层的密度分布 （Li et al. PoF 22, 2010)边界层有清晰锐利的界边界层有清晰锐利的界面（层流区、湍流区面（层流区、湍流区“泾渭分明泾渭分明”）t湍流信号湍流信号层流信号层流信号内间歇性：内间歇性： 湍流脉动的概率密度分布偏湍流脉动的概率密度分布偏离离Gauss分布（随机分布）分布（随机分布）)log(PDF u2xe概率论（中心极限定理）概率论（中心极限定理）独立随机事独立随机事件满足件满足Gauss分布分布偏离偏离Gauss分布分布 湍流不是完全随机的湍