1、計算流體力學講義2011 第十一講 湍流與轉(zhuǎn)捩 （1） 李新亮 ；力學所主樓219； 82543801,知識點：,1,講義、課件上傳至 (流體中文網(wǎng)） - “流體論壇” -“ CFD基礎理論 ” 講課錄像及講義上傳至網(wǎng)盤 http:/cid-,Copyright by Li Xinliang,線性穩(wěn)定性理論 轉(zhuǎn)捩的預測方法 壁湍流轉(zhuǎn)捩的渦動力學機制,Copyright by Li Xinliang,2, 11.1 線性穩(wěn)定性理論,一、 穩(wěn)定性基本概念,常識：流體中的不穩(wěn)定性,K-H不穩(wěn)定性,A. K-H (Kelvin-Helmholtz）不穩(wěn)定性 自由剪切流的無粘不穩(wěn)定性,混合層 K-H不穩(wěn)

2、定性,K-H不穩(wěn)定性的關鍵： 速度剖面有拐點,Lee-Lin: 速度剖面的拐點是無粘不穩(wěn)定性的必要條件,流體不禁搓，一搓搓出渦,已知某運動狀態(tài)； 在此基礎上施加微小擾動； 如擾動隨時間（或空間）衰減，則稱系統(tǒng)穩(wěn)定，否則為不穩(wěn)定,Copyright by Li Xinliang,3,自然界中 K-H不穩(wěn)定性圖片,智利塞爾扣克島的卡門渦街,澳大利亞Duval山上空的云,KelvinHelmholtz instability clouds in San Francisco,佛蘭格爾島周圍的卡門渦街,高速流,低速流,自由剪切層受到擾動界面變形后的情況 K-H不穩(wěn)定性的產(chǎn)生機理,受阻減速，壓力升高，產(chǎn)生

3、高壓區(qū),高壓導致變形加劇,Copyright by Li Xinliang,4,B. T-S (Tollmien-Schlichting) 不穩(wěn)定性不可壓 壁面剪切流的粘性不穩(wěn)定性,Mack 不穩(wěn)定性 超聲速壁面剪切流的不穩(wěn)定性,不可壓邊界層速度剖面 （Blasius解） 無拐點,可壓縮情況 Mach數(shù)足夠高時會出現(xiàn)廣義拐點 出現(xiàn)無粘不穩(wěn)定性,不可壓縮無粘不穩(wěn)定性 需存在拐點 可壓縮無粘不穩(wěn)定性需存在廣義拐點,Mach 6 鈍錐（1攻角） 不同子午面 的分布,超音速平板邊界層的不穩(wěn)定波,第1模態(tài)（T-S波）,第2模態(tài) （Mack模態(tài)）,Copyright by Li Xinliang,5,激波

4、,密度界面,R-M (Richtmyer-Meshkov)不穩(wěn)定性 激波與密度界面作用的斜壓效應,慣性約束聚變（ICF）示意圖,小知識 渦的產(chǎn)生機制： 粘性、 斜壓、有旋的外力,激波,密度界面,斜壓項,Copyright by Li Xinliang,6,D. R-T (Reyleigh-Taylor)不穩(wěn)定性 重力帶來的不穩(wěn)定性,R-T (Reyleigh-Taylor)不穩(wěn)定性,重介質(zhì),輕介質(zhì),Copyright by Li Xinliang,7,E Barnard熱對流不穩(wěn)定性,其他學科的不穩(wěn)定性：,Euler壓桿的不穩(wěn)定性,Barnard 熱對流的胞格結構,板殼的不穩(wěn)定性,Copyri

5、ght by Li Xinliang,8,二、 穩(wěn)定性問題的常用數(shù)學方法 線性穩(wěn)定性分析,Step 1: 得到線性化的擾動方程,控制方程為：,已知其具有解,最好是精確解，也可用高精度的數(shù)值解,令：,舍棄高階小量，得到線性化的擾動方程,（1）,例如： 平板的Blasius解，槽道的Poiseuille 解,線性方程,Copyright by Li Xinliang,9,Step 2: 求解 的特征值問題,什么條件下具有非零解，非零解如何？,通常假設在某些方向具有周期性，轉(zhuǎn)化為一維問題,數(shù)值方法： 將 （1） 離散代數(shù)方程 何時有非零解， 非零解如何？ 特征值問題,什么條件下有非零解？,特征值問題

6、計算量巨大，目前通常只能處理一維問題,Copyright by Li Xinliang,10,三、 穩(wěn)定性問題示例 不可壓縮槽道流動的線性穩(wěn)定性（LST）理論 (以二維為例),Step 1: 獲得線性化擾動方程,令：,Poiseuille解：,(2),代入方程（2），并舍去高階小量得到線性化的擾動方程,（3）,1） 控制方程及邊界條件,Copyright by Li Xinliang,11,研究擾動發(fā)展的空間模式和時間模式,擾動源,空間模式： 任一點的擾動具有時間周期性 符合物理條件,假設擾動具有如下形式：,沿流向及時間方向具有波動特性 稱為Tolmien-Schlichting（T-S）波,

7、任意擾動可分解為正弦波的疊加 線性系統(tǒng)各成分無相互作用 可獨立研究,為實數(shù),為復數(shù),擾動波的振幅沿流向指數(shù)變化,空間增長率,時間模式： 擾動具有流向的周期性 假設一窗口沿流向運動，研究窗口內(nèi)擾動的演化,為實數(shù),為復數(shù),擾動波的振幅雖時間變化,時間增長率,Copyright by Li Xinliang,12,以時間模式為例：,（4）,（5）,（6）,線性偏微方程（3）轉(zhuǎn)化成為含參數(shù)的線性常微方程組（4）-（6）,譜方法的常規(guī)做法,通過消元法，轉(zhuǎn)化為更高階的常微方程 （不是必須的）,常用做法，通常還可以反向為之： 高階方程轉(zhuǎn)化為低階方程組,消去,Orr-Sommerfeld（O-S) 方程,其中

8、：,最終，控制方程為O-S方程：,Copyright by Li Xinliang,13,邊界條件：,y=1 （固壁）:,y=0 （中心線，對稱）：,可以取計算域-1,1,使用固壁邊界條件； 也可以取計算域-1,0,使用固壁及對稱邊界條件,流函數(shù)形式的O-S方程,引入流函數(shù)，使得：,計算出 后，利用公式,計算其他兩個量,則：,令：,常數(shù)倍,滿足的方程及邊界條件與 完全相同。,如果 恒大于（或恒小于0），則必有,Copyright by Li Xinliang,14,小知識： 關于O-S方程,1） O-S方程適用于不可壓平行流的穩(wěn)定性問題 （不僅槽道流） 2） 準平行流 （流線沿x方向接近平行）

9、也可使用（例如邊界層流動） 3）如果舍去粘性（左端）項，則方程稱為Rayleigh方程,Rayleigh拐點定理： Rayleigh方程存在不穩(wěn)定解的必要條件是速度型存在拐點。即存在某點 使得,若存在無粘不穩(wěn)定性，該項必有0點。,分部積分，并取虛部，得：,不存在非穩(wěn)定解,Copyright by Li Xinliang,15,2） O-S方程的解法,數(shù)學表述 奇性（特征值）問題： 參數(shù) 為何值時，方程有非零解？ 非零解如何？,時間發(fā)展槽道湍流： (通常) 給定Re及a，問 w取何值時，O-S方程有非零解？,增長率,求解步驟： 1） 將O-S方程離散，得到線性代數(shù)方程組 離散方法： 差分法、有限

10、元法、譜方法、打靶法 2） 求w，使得該方程有非零解（奇性或特征值問題）.,求出w,局部法：只求出一個w 全局法：計算出全部的w,試計算 的時間發(fā)展槽流中 （即波長2p）T-S波的頻率及增長率,Copyright by Li Xinliang,16,四: 例題,Step 1: 離散 （差分法）,一維問題,網(wǎng)格： 均勻網(wǎng)格 （簡單，但需要較多網(wǎng)格點 N 300） 非均勻網(wǎng)格,差分離散：,2階格式,4階格式：自行推導 （利用小程序）,非均勻網(wǎng)格：,問題： 會產(chǎn)生大量非物理解 （例： 1000個網(wǎng)格點算出1000個特征值） ； 可通過不同網(wǎng)格的對比，進行篩選,Copyright by Li Xinl

11、iang,17,離散化后得：,Step 2: 求廣義解特征值問題（1） 即w為何值時（1）有非零解,(1),方法1： 全局法 一次計算出全部特征值,常用Q-Z分解法; 其他： 冪法、反冪法、Jacobi法，Householder ),狹義特征值問題,廣義特征值問題,求解特征值是計算數(shù)學的主要研究方向，有大量成熟的方法,可借助軟件包或Lapack庫等 （自行到網(wǎng)上搜索）,通常關心最不穩(wěn)定的擾動波,最大的那個,Copyright by Li Xinliang,18,方法2： 局部法,令：,方法： Newton法， 弦位法，拋物線（Muler）法,Newton 法：,弦位法：,差分化,拋物線法：,已

12、知：,可連成一條拋物線,令：,求出新的值,消元法計算行列式（利用5對角特征）,計算出 后，求解方程（1）就可得到特征向量。,（1）,方程有奇性，可補充一個條件，例如給定某個點的值,Copyright by Li Xinliang,19,效果更好的方法 Malik提出的緊致差分格式求解 參考文獻: 周恒等： 流動穩(wěn)定性 （p. 10-13），國防工業(yè)出版社,非均勻網(wǎng)格的超緊致格式（4階精度）,優(yōu)點： 緊致網(wǎng)格基僅2點 精度高 4階 直接適用于非均勻網(wǎng)格無需坐標變換,原方程組：,形式變換 變?yōu)橐浑A方程組,通常的做法,推導倉促 可能有誤 請仔細推導,與y無關,Copyright by Li Xinl

13、iang,20,令：,帶入差分格式：,得：,Copyright by Li Xinliang,21,令：,邊界處表達式 (C1,D1,CN,DN) 根據(jù)邊界條件而定,內(nèi)點的表達式,特點： 離散形成的代數(shù)方程組呈 塊兩對角特征,求行列式，特征值都非常便利！,其余步驟與前文相同 可用矩陣廣義特征值理論計算全部模態(tài) 也可利用 計算單個模態(tài),計算域可以取為-1,1。 也可取為-1,0， 在中心線給對稱（或反對稱）邊界條件,Copyright by Li Xinliang,22,邊界條件,1） 如果取完整計算域 -1,1,2) 如果取一半計算域-1,0,處可設定對稱或反對稱邊界條件,如設定對稱條件，只能

14、計算出對稱擾動模態(tài) 如設定反對稱條件，只能計算出反對稱模態(tài),對于槽道流，最不穩(wěn)定模態(tài)是對稱模態(tài),但有些情況下，穩(wěn)定模態(tài)在轉(zhuǎn)捩過程中也發(fā)揮作用（例如感受性過程，見Zhong et al. JFM 556,55-103,2006）,按照內(nèi)點方法計算,先根據(jù)處理內(nèi)點的方法，求出所有點上的C,D值，再對邊界點進行特殊處理（D的前兩行設為0， C的前兩行設為(1,0,0,0)及（0，1，0，0）,Copyright by Li Xinliang,23,具體解法（局部法）,求解,顯然,因此只需計算每個4*4矩陣的行列式即可 （可直接寫出表達式，也可用消元法計算）,編制好計算 的子程序后，可利用Newton

15、法，弦位法，拋物線法等求解 ，得到復增長率,Copyright by Li Xinliang,24,計算出 后，利用消元法求解方程 ，即可得到特征向量,消元過程中，充分利用兩對角塊矩陣的性質(zhì)，可減少計算量,獨立消元,注： 由于方程有奇性，消元后，最后一個方程為0=0, 舍棄最后一個方程，并令最后一個未知數(shù)為1 （ ），即可解出 顯然 （a為任意常數(shù)）也是原方程的解，因此 可用某一值（例如 ）歸一化。,最終，得到 條件下，波數(shù)為 的最不穩(wěn)定擾動波：,擾動波型函數(shù)的分布,Copyright by Li Xinliang,25, 11.2 邊界層轉(zhuǎn)捩的預測方法,1. 經(jīng)驗公式法,轉(zhuǎn)捩位置,Mach

16、6 鈍錐邊界層表面的摩擦系數(shù)分布 (Li et al. Phys. Fluid. 22， 025105， 2010. Li et al. AIAA J. 46（11），2899-2913，2008),x,摩阻或熱流,轉(zhuǎn)捩起始點（transition onset),轉(zhuǎn)捩峰（transition peak),充分發(fā)展湍流區(qū),球錐的轉(zhuǎn)捩Reynolds數(shù),邊界層外緣的Mach數(shù),動量厚度定義的轉(zhuǎn)捩Reynolds數(shù),Copyright by Li Xinliang,26,2. eN 方法,LST理論：,積分起始點,擾動波進入中性曲線后，開始增長，局部增長率為,eN理論：擾動波增長到eN倍，即發(fā)生轉(zhuǎn)捩

17、,N值需要由實驗（或經(jīng)驗）確定，通常為811, 即擾動增長4個量級（10000倍）左右。 在不可壓縮及航空領域（亞、跨、超）較為成熟。 在航天領域（高超聲速），還有待檢驗。,不足之處： 未考慮擾動波進入中性曲線前的衰減過程，沒考慮感受性過程。,他人的改進： 蘇彩虹，周恒等 考慮衰減過程 C. H. Su, and H. Zhou, Science in China G, 52 (1):115-123 (2009).,Copyright by Li Xinliang,27,3. PSE （拋物化擾動方程）法,優(yōu)點： 1） 無需平行流假設； 2） 可處理非線性(N-PSE),Step 1: 得到擾

18、動的控制方程,已知解,線性化,L-PSE,N-PSE,Step 2: 假設擾動具有波動形式,振幅，沿x方向是個緩變量 （相對y方向而言）,Step 3: 帶入擾動方程，得到振幅 的控制方程,“緩變量”是個很有用的概念，可用來簡化方程,Plantdl邊界層理論就是利用“緩變量”的概念進行簡化的。,LST的 方程是一維的 PSE的 方程是二維的,Step 4: 利用緩變量性質(zhì)，舍棄方程中的橢圓項（為高階小量），得到拋物化的擾動方程,沿x方向推進求解 （類似時間方向的處理），計算量相當于一維問題。 （“拋物化”的優(yōu)勢）,非線性項的處理方法與譜方法相似,Copyright by Li Xinliang

19、,28,4. 轉(zhuǎn)捩模型法 （間歇因子模型）,實際粘性系數(shù),層流粘性系數(shù),湍流粘性系數(shù) （由湍流模型給定）,湍流間歇因子 （0表示純層流，1表示純湍流）,方法 1） 根據(jù)經(jīng)驗公式，給定 沿流向的分布 方法2） 給出 的發(fā)展方程，進行求解,Copyright by Li Xinliang,29,常識： 湍流的間歇性,外間歇性： 層流及湍流交替出現(xiàn)的現(xiàn)象,Mach 6 鈍錐邊界層的密度分布 （Li et al. PoF 22, 2010),邊界層有清晰銳利的界面（層流區(qū)、湍流區(qū)“涇渭分明”）,湍流信號,層流信號,內(nèi)間歇性： 湍流脈動的概率密度分布偏離Gauss分布（隨機分布）,概率論（中心極限定理）獨立隨機事件滿足Gauss分布