可壓縮邊界層的拋物化穩定性方程及其在某些問題中的應用張永明2008年10月摘要推導出可壓縮邊界層的拋物化穩定性方程（PSE），驗證了PSE方法的可靠性，再將其應用于非平行效應、二次失穩和轉捩預測三個具體問題中。發現了邊界層非平行性對中性曲線的影響主要集中在臨界雷諾數處；驗證了二次失穩機制的存在，得到了二次失穩問題中三維亞諧波的放大率隨其展向波數和二維基本波幅值的變化關系；得到了準確的轉捩起始位置，并再現了轉捩中breakdown過程的機理。PSE方法是研究有關邊界層內擾動演化的多種問題的有效工具。1引言PSE是研究基本流中擾動演化的比較新的一種方法，此外常用的方法還有線性穩定定性理論（LST）和直接數值模擬（DNS）。與另外二者相比，PSE的優點在于同時具備以下特點：可以考慮基本流的非平行性，可以用于非線性問題，計算量不大。因此，從提出PSE至今不長的時間里，PSE方法得到了充分的發展，并應用到多個問題的研究中。Herbert和Bertolotti（1987）首先提出了拋物化穩定性方程。他們研究的是不可壓邊界層的線性PSE，這是一種比較簡單的情況。Bertolotti等（1992）用線性PSE在大雷諾數處得到的小擾動演化與線性穩定性理論和直接數值模擬的結果一致，驗證了線性PSE計算的可靠性，如圖1-1所示。Bertolotti等（1992）用PSE研究不可壓邊界層的非平行性對其穩定性的影響，得到的中性曲線比用LST得到的更符合Schubauer和Skramstad（1943，1947）、Ross等（1970）、Strazisar等（1977）和Kachanov等（1977）的實驗結果，也與Gaster（1974）、于秀陽和周恒（1986）用另外一種方法得到的計算結果一致，如圖1-2所示。Bertolotti等（1992）提出了不可壓邊界層的非線性PSE，所得擾動演化結果與直接數值模擬所得一致，如圖1-3所示，用計算的方法驗證了非線性PSE的可靠性。而前者計算花費的時間比后者少很多，這顯示出PSE在研究某些非線性問題時的優勢：準確且快速。Joslin等（1992）和Esfahanian等（2001）也做了同樣的驗證。Malik等（1999）將非線性PSE計算所得與實驗測量結果對比，發現二者一致，又從實驗上驗證了其可靠性。在此基礎之上，Malik等（1999）開始用PSE作為可靠的工具來研究不可壓邊界層的穩定性和轉捩問題。Herbert（1997）對不可壓縮流PSE的研究進行了總結。圖1-1不同方法得到的大雷諾數處小擾動增長率。實線為Bertolotti等（1992）的線性PSE結果，虛線為LST結果，圓形標記為DNS結果。圖1-2不同方法得到的中性曲線。實線為Bertolotti等（1992）的線性PSE結果，虛線為LST結果，圓形標記為Gaster（1974）的結果。圖1-3Bertolotti等（1992）的非線性PSE計算與空間模式DNS的對比，包括平均流修正，一階擾動，二階擾動和三階擾動。實線為NPSE結果，方形標記為DNS結果。相對不可壓邊界層，可壓縮邊界層的PSE更為復雜，其結果也較少。但由于航空航天技術發展的需求，這方面的研究逐漸多了起來。Bertolotti和Herbert（1991），Bertolotti（1991）提出了可壓縮平板邊界層的線性PSE，并發現用線性PSE在大雷諾數處得到的小擾動演化與線性穩定性理論的結果一致，驗證了線性PSE計算的可靠性。可壓縮邊界層的非線性PSE研究工作比較少，目前見到的工作主要來自美國國家航空航天局（NASA）Langley研究中心的Chang等人。Chang等（1991）提出了可壓縮邊界層的非線性PSE。為了驗證其可靠性，Chang等（1991，1993）、Pruett等（1995）、Pruett和Chang（1995）、Jiang等（2003，2004）分別將非線性PSE結果與DNS的對比，發現二者一致，但PSE計算花費時間少得多。遺憾的是，上述關于可壓縮邊界層的線性和非線性PSE的驗證工作都不是足夠充分、完整和嚴格的，其原因將在本文第三部分做具體闡述。本文以可壓縮平板邊界層為研究對象。先從完全的Navier－Stokes方程出發，推導出線性和非線性的拋物化穩定性方程（見第二部分），再介紹具體的數值計算方法（見第二部分）。然后對線性PSE和非線性PSE計算方法做充分、完整和嚴格的檢驗，驗證其可靠性（見第三部分）。此后則以PSE作為工具研究具體問題。應用線性PSE研究基本流的非平行性對可壓縮邊界層中性曲線的影響（見第四部分）。應用非線性PSE研究超音速邊界層的二次失穩問題（見第五部分）。應用非線性PSE研究轉捩可壓縮邊界層的轉捩預測問題（見第六部分）。

2可壓縮邊界層PSE的控制方程和數值求解方法這一部分先介紹可壓縮邊界層的拋物化穩定性方程的推導過程，以及求解方程所用的數值方法。這些是前人做過的工作，但都是PSE方法的基礎，所以有必要做詳細的介紹。首先推導出線性和非線性PSE的控制方程（見2.1）；然后說明數值計算中所用的基本流和入口條件（見2.2）；最后闡述計算網格和差分格式（見2.3）。2.1控制方程推導拋物化穩定性方程需經四個步驟：（1）先由有量綱的完全的N-S方程出發，利用特征量得到無量綱的方程，（2）再推導出完全的擾動方程，（3）然后針對擾動的特點推出穩定性方程，（4）最后利用邊界層厚度增長緩慢的性質將其拋物化而得到拋物化穩定性方程。下面作具體的介紹。用上標表示有量綱的量，在三維笛卡爾坐標系下用x，y，z表示流向、法向和展向，則有量綱的完全的N-S方程為（2.1-1）其中是速度矢量，是時間，是密度，是壓力，是溫度，是定壓比熱，是熱傳導系數，是第一動力粘性系數，是第二動力粘性系數，為普適氣體常數。為粘性耗散函數，表示為，（2.1-2）其中為應變率張量，表達式為。（2.1-3）方程（2.1-1）中第一到第四個方程分別為連續性方程、動量方程、內能方程和狀態方程。本文假設可壓縮流體為完全氣體，使用完全氣體的狀態方程。下面要將N-S方程無量綱化，其中用到的特征量如下。取特征長度為，（2.1-4）其中為平板前緣到計算域入口的距離，為運動粘性系數，運動粘性系數與第一動力粘性系數關系為，為流向速度，下標e表示邊界層外的自由來流中的量。其它的原始特征量為：自由流速度、自由流溫度、自由流密度，自由流第一動力粘性系數，自由流第二動力粘性系數，自由流熱傳導系數。導出的特征量為：特征壓力為，特征時間為。用有量綱量除以特征量即為無量綱量，如。（2.1-5）將上式代入有量綱N－S方程（2.1-1），即得無量綱N－S方程。無量綱N-S方程會用到幾個無量綱參數，即入口出雷諾數、普朗特數和馬赫數，下面分別介紹。平板前緣到入口處的無量綱距離為，（1）而入口處的雷諾數的表達式則為。（2由一直接化簡即可）平板前緣到入口下游某處的距離為，顯見。由此得到各變量對x的一階、二階導數的量級分別為。（2.1-6）這將是后來對方程進行拋物化的依據。普朗特數表達式為。馬赫數表達式為，其中為自由來流中的聲速，其表達式為，其中為比熱比，其表達式為，其中為定容比熱。在整理無量綱N-S方程時，還要注意粘性系數和熱傳導系數的處理。本文使用了Stokes假設（見《流體力學》上冊，68頁）（周光炯2000），第二動力粘性系數可用第一動力粘性系數表示。（2.1-7）而無量綱的第一動力粘性系數則使用Sutherland公式，假設其隨溫度的變化滿足如下關系式，（2.1-8）其中C1＝110.4K為常數。本文中無量綱的熱傳導系數也使用Sutherland公式，形式為。（2.1-9）一般來講可以取C2＝C1，則無量綱的熱傳導系數與第一動力粘性系數相同，（2.1-10）以后的公式中可以都用粘性系數表示。整理之后的無量綱N-S方程形式如下(此種寫法挺好)（a）（b）（c）（d）（e）（f）（2.1-11）其中（a）為連續性方程，（b）、（c）和（d）分別為流向、法向和展向的動量方程，（e）為內能方程，（f）為狀態方程。將（2.1-11）中的狀態方程代入連續性方程、動量方程和內能方程，消去壓力p，則無量綱方程化簡為5個方程，其中的變量也是5個，即：（a）（b）（c）（d）（e）…（2.1-12）方程中有一點需要注意，無量綱長度x、y、z、時間t、入口雷諾數Re0等用到特征長度的量，都與常見的直接數值模擬中的對應量不同。邊界層的直接數值模擬常用入口處位移排移厚度作為特征長度，（常用有兩種，這是其中一種）其表達式為，（2.1-13）其中為馬赫數M的函數，其表達式較復雜，詳見文獻（黃章峰2003第17頁，黃章峰2006第26頁）。對比式（2.1-4）可見，與本文所用特征長度相差倍。那么無量綱量x、y、z、t、Re0等也相差倍。本文后面還會提到一些用特征長度無量綱化的量，也會與DNS中的對應量相差倍，這一點以后要注意。將方程（2.1-12）中的瞬時量設為定常基本流與擾動量之和（2.1-14）其中上標表示基本流，上角標表示擾動量。將（2.1-14）式代入N－S方程（2.1-12），并減去定常基本流滿足的方程，得到擾動方程（方便起見，較復雜的方程用矩陣形式表達）（2.1-15）其中為擾動矢量。系數矩陣Vxx，Vyy，Vzz，Vxy，Vxz，Vyz，Г，A，B，C，D都可分解為線性部分與非線性部分，其中線性部分只包含定常基本流的量，非線性部分則包含擾動量。二者分別用上標l和n表示，則這些矩陣可記為。通過量級分析可以發現，對流項和粘性項的量級分別為1和1/Re0。系數矩陣Г，A，B，C，D同時包含對流項和粘性項，其總的量級為1，系數矩陣Vxx，Vyy，Vzz，Vxy，Vxz，Vyz只包含粘性項，其量級為1/Re0。將線性和非線性部分分別放在方程兩邊，則擾動方程（2.1-15）整理為如下形式（2.1-16）其中右邊的矢量函數表示所有的非線性項，其表達式如下（2.1-17）非線性擾動方程（2.1-16）的系數矩陣和非線性項矢量的元素詳見附錄1。其中需要注意的是μ’的處理。本文利用Sutherland公式（2.1-8）將μ’表示為和的函數，下面介紹詳細的推導過程。由Sutherland公式可知是的函數，（2.1-18）相對來說是很小的量，可用的微分近似，（2.1-19）考慮到，代入上式，則得到的近似表達式，（2.1-20）或簡寫為，（2.1-21）其中，是基本流溫度的函數，即。（2.1-22）（第一部分工作是方程無量綱化，變為關于五個未知量的五個方程，注意特征長度的選擇和以前略有不同，再就是將瞬時量設為定常基本流與擾動量之和減去基本流得到擾動方程）2.1.1線性拋物化穩定性方程對于小擾動情況，非線性擾動方程（2.1-16）中的非線性部分可以略去，從而得到線性擾動方程（注：對于小擾動，可以忽略不同頻率不同展向波長波的之間的非線性作用，對于有限幅值擾動，非線性作用不可忽略，線性拋物化穩定性方程不再適用）（2.1-23）這是一個橢圓型方程（？如何判斷一個方程是何種方程），對其求解需要知道計算域所有邊界上的條件，且計算量很大。考慮到要研究擾動沿空間的演化，假設擾動是時間和展向的周期函數，（2.1-24）其中為復數形式的形狀函數(在哪定義的？)矢量，為復數形式的流向波數，其虛部為i，-i為增長率，為展向波數，ω為頻率，c.c.表示復共軛。此形式與空間模式線性穩定性理論（找是什么）的行波表達式類似，但由于考慮邊界層的增長，形狀函數和波數都是x的緩變函數。將（2.1-24）式代入線性擾動方程（2.1-23）式，得到線性穩定性方程（如何得到？為什么啊），（2.1-25）其中系數矩陣為。至此，線性穩定性方程（2.1-25）仍然是橢圓型的（這又是為什么？），下面介紹如何將其拋物化。用（2.1-6）式估計方程（2.1-25）中各導數，并考慮方程各項系數的量級，得如下的量級估式。11（2.1-26）（，，由矩陣形式判斷量級皆為）（，量級為1？有原因嗎）雷諾數一般是一個大數，則是個小量，是二階小量。去掉方程中二階以上小量（即含的項），就得線性拋物化穩定性方程（有一個特征值為0就是拋物型）。（2.1-27）由于方程（2.1-27）已是拋物型的，計算時只需定出入口條件，即可向下游逐步推進計算，工作量比解橢圓型方程要小得多。線性PSE方程（2.1-27）中系數矩陣的元素詳見附錄2。去掉量級小量的過程即為拋物化（這么定義的？）。它實際上利用了下述事實。由于基本流沿流向只是逐漸改變，擾動形狀函數的變化也是緩慢的，以致其二階導數相對于其一階導數可以忽略。壁面處使用無滑移邊界條件和等溫或絕熱條件。使用等溫條件時為，（2.1-28）使用絕熱條件時為。（2.1-29）邊界層外的自由流條件為。（2.1-30）另外，在上下邊界層都用連續性方程作為邊界條件，使得方程封閉。方程（2.1-27）中有兩個未知函數，即形狀函數和波數，一個方程不足以解這兩個量，這就需要補充一個條件。通常使用的有以下四種補充條件，（如何確定選哪一個），（2.1-31），（2.1-32），（2.1-33），（2.1-34）其中上標c表示復共軛。具體應用某一補充條件時，需要對波數進行迭代求解，直至其收斂于某一固定值。上面四個補充條件分別對應下面四個迭代公式，（2.1-35），（2.1-36），（2.1-37），（2.1-38）其中E為擾動能量。本文采用補充條件（2.1-34），這相當于由形狀函數決定的能量沿流向幾乎是不變的。于是擾動能量的增長就完全體現在波數的虛數部分，而擾動的相位變化則體現在波數的實數部分中。實際應用這一條件（2.1-34）時要用式（2.1-38）算出一新的波數αnew。然后返回去用方程（2.1-27）重算一個新的解。如此反復迭代，直至的變化小于預設的小量（本文在實際計算中取）。具體的講，迭代過程分為以下三個步驟（以流向的第i+1個節點處的計算過程為例）：（1）在xi+1處預估一個流向波數αold，本文取的是前一個流向節點處的波數，即。用這個波數計算得到當地的形狀函數。（2）利用波數迭代公式（2.1-25），用剛得到的形狀函數和已有的流向波數αold解出一個新的波數αnew。（3）判斷新舊流向波數之間的差是否小于預設小量ε。若大于ε，則返回步驟（1），用最新的波數再計算一個解。如此反復迭代，直至小于ε，即認為現有波數和形狀函數已足夠準確，可取為當地真解。2.1.2非線性拋物化穩定性方程在上面的線性PSE中用到了一小擾動假設，即擾動幅值很小，以至于不同頻率、不同展向波長的波之間的非線性作用可以忽略。對于有限幅值擾動的演化，線性PSE則不再適用，而要用非線性PSE進行研究。仍然假設擾動是時間和展向的周期函數。非線性作用將激發高次諧波，它們的波數和頻率與基本波的波數和頻率將是整倍數的關系，因此擾動展成如下的Fourier級數形式，（2.1-39）其中表示階數為（m，n）時的形狀函數，表示階數為（m，n）時的流向波數，其頻率為，展向波數為。當然計算中Fourier級數只能取有限多項，（2.1-40）其中M和N表示級數截斷的階數。為了提高運算速度，本文使用了快速Fourier變換（FFT）。將（2.1-40）擾動展開式代入非線性擾動方程（2.1-16），再去掉方程中含的項，并將波數和頻率相同的項放在一起，則對階數為（m，n）的波，其形狀函數滿足如下的拋物化穩定性方程，（2.1-41）其中為階數為（m，n）時的非線性項，其表達式為，（2.1-42）其中是擾動方程（2.1-16）中非線性項的Fourier級數。（2.1-43）（這是為什么不含波動項？）由于開始時非線性項是外差得到的，不是非常準確，所以也要進行迭代。與前面提到流向波數的迭代過程類似：先用一個不太準的非線性項計算出當地的形狀函數的波數，再用剛得到的形狀函數和波數計算出一個新的非線性項，然后再用這個非線性項算形狀函數和波數，直至結果不再變化為止。非線性PSE中流向波數的處理稍復雜一些。如果入口處擾動波加在階數（m0，n0）上，則對應的波數與線性PSE一樣，要用（2.1-38）進行迭代。其它階數的波數不用迭代，而直接使用“相速度相等”假設，（2.1-44）即假設其它階數的波的相速度（即）與所加擾動波的相同。由于是實數，所以方程（2.1-41）中只需對的階數進行計算，其它的階數為它們的共軛。每一個階數（m，n）下的非線性PSE，都與頻率為，展向波數為的線性PSE基本相同，只是方程右邊多了一個非線性的驅動項。邊界條件、方程求解步驟都與線性PSE的基本相同。但m=n=0的波實際上是對平均流的修正，求解時上邊界的自由流條件應改為。（2.1-45）此邊界條件使得平均流保持質量平衡。2.2基本流和入口條件PSE是研究擾動演化的方法，需要事先知道整個計算域的基本流和入口的擾動條件。本文研究的是可壓縮平板邊界層中三維擾動的演化，而基本流是二維的。它可以取Blasius相似性解，或者取由直接數值模擬得到的定常基本流。Blasius相似性解的解法可文獻（周恒，趙耕夫2004）。這樣計算得到的邊界層厚度沿流向增長，與到平板前緣距離的二分之一次方成正比，基本流顯然不是平行流。拋物化穩定性方程的好處之一是可以考慮基本流的非平行性。由于邊界層方程忽略了某些高階小量，所以Blasius相似性解并不是N-S方程的準確解。將Blasius解作為用N-S方程計算基本流時的入口條件，在計算域入口附近會有一個調整過程，即有一個過渡區。真正的計算域應從這過渡區后開始。入口條件的選取是PSE中一個重要的問題。因為PSE只能計算擾動演化，事先需要給出計算域起始位置的擾動。一般都選取用線性穩定性理論得到的T-S波作為入口擾動。雖然其中用了平行流假設，因而它不是PSE方程的精確解，但在迭代過程中，很快就能調整過來。對于線性拋物化穩定性方程，入口處擾動幅值可以任意給定，本文所取幅值與LST中一樣，為1。而對于非線性PSE，入口處擾動需按實際情況給出具體值。2.3計算網格和差分格式2.3.1計算網格流用等間距的網格。對于線性計算，流向步長可以取1/3擾動波長或更大，而對于非線性計算，流向步長一般可取1/20的擾動基本波波長或更小。法向用變間距網格，靠近壁面處網格較密。為了在法向便于使用高精度等間距差分格式，一般將方程變換到等間距的計算坐標=(y)下。本文所取變換形式如下，（2.3-1）其中為y方向計算域長度，為y方向網格數，，，b為常數，其表達式為，（2.3-2）其中常數k為上邊界網格寬度和下邊界網格寬度之比：。（2.3-3）以=40，=40，k=300為例，圖（2.3-1）畫出了y方向的網格分布。圖2.3-1y方向變間距網格分布2.3.2差分格式（選用何種格式有理由嗎？）流向的一階導數差分化時可用一階或二階精度的單邊差分格式如下：，（2.3-4）。（2.3-5）在法向對坐標的一階和二階導數，用如下的四階精度的中心差分格式，（2.3-6）。（2.3-7）在臨近上下邊界處要適當降低精度。在臨近下邊界的點用如下的差分格式，（2.3-8）。（2.3-9）在臨近上邊界的點用如下的差分格式，（2.3-10）。（2.3-11）在下邊界點，一階導數使用如下的差分格式。（2.3-12）在上邊界點，一階導數使用如下的差分格式。（2.3-14）函數對y的一階和二階導數，則可用下式計算。。（2.3-15）（2.3-16）2.3.3數值不穩定性PSE的方程（2.1-27）和（2.1-41）的求解過程中有一個重要的問題，就是數值不穩定性。PSE方程不是完全的拋物型方程，而存在著殘余橢圓性，殘余橢圓性造成計算中會出現很強的數值不穩定性。殘余橢圓性的存在可以從數值計算和數學原理兩個方面解釋。一方面，前面提到，PSE的系數中存在流向波數，而是事先未知的，需要通過迭代求出。這說明方程的解和方程的系數互相影響，因此方程會存在橢圓性。另一方面，將PSE方程（2.1-27）和（2.1-41）所對應的擾動方程的線性部分單獨拿出來，對其作特征分析可以看出（如何分析？），有一個特征值在靠近壁面處不是正實數，因此方程存在橢圓性。Chang等（1991，1993）具體地介紹了PSE的殘余橢圓性，并提出了消除數值不穩定性的方法。事實上，Chang等只是引用了Vigneron等（1978）在解決拋物化N－S方程（PNS）的殘余橢圓性時所提出的方法，所以Chang等也稱之為Vigneron技術。Haj-Hariri（1994），Li和Malik（1996）則從數學上證明了殘余橢圓性的存在，且證明了Vigneron技術在消除數值不穩定性的同時，仍保持了解的準確性。對于線性PSE來說，當流向步長大于1/3波長時（步長取1/3波長或更大時，計算結果仍可以足夠準確），殘余橢圓性較弱，計算中基本上沒有數值不穩定性，這時可以不用Vigneron技術。當流向步長小于1/3波長時，殘余橢圓性較強，計算中會出現數值不穩定性，這時應該用Vigneron技術。對于非線性PSE來說，只有步長小于1/20基本波長的時候，計算結果才足夠準確，所以必須使用Vigneron技術。具體地說，Vigneron技術就是在流向動量方程中去掉部分或全部的流向壓力梯度，即讓乘以一個大于等于0且小于等于1的實數，即Vigneron系數。因為PSE方程已經用狀態方程消去壓力，所以在PSE線性部分，Vigneron系數都乘在上，而在非線性項中，Vigneron系數則乘在上。Vigneron系數V的取值一般使用下面公式（2.3-17）其中，是當地馬赫數，是當地聲速。如果為了方便，將所有的Vigneron系數都取0也可以，計算結果仍然準確。

3數值計算的檢驗PSE作為一種新的計算方法，其數值計算的結果需要經受檢驗。下面分別介紹對線性和非線性PSE計算的檢驗。3.1線性PSE計算的檢驗用線性PSE得到的可壓縮邊界層中小擾動的演化結果與傳統的LST的結果進行比較。因為在大雷諾數處邊界層近似平行，LST結果足夠準確，可以作為檢驗的標準。為了循序漸進，先檢驗線性PSE在平行流中的計算結果（3.1.1），再檢驗線性PSE在非平行流中的計算結果（3.1.2）。3.1.1線性PSE在平行流中計算的檢驗把線性PSE方法用于平行流中擾動演化的計算，并將其結果與線性穩定性理論的結果相比較。PSE沿流向要推進計算，LST只計算一個剖面下的結果即可。對亞音速（馬赫數M＝0.3）邊界層和超音速（馬赫數M＝4.5）邊界層中的擾動演化結果都進行了檢驗。對于亞音速邊界層來說，不穩定擾動只有一個模態。對于馬赫數為4.5的超音速邊界層來說，不穩定模態有兩個，通常稱作第一和第二模態。3.1.1.1線性PSE在亞音速平行流中計算的檢驗基本流取為平板二維Blasius相似性解，但將其局部平行化，馬赫數M＝0.3，自由來流溫度255.7K，普朗特數Pr＝0.707，雷諾數Re0＝1706；入口擾動為平行流假設下線性穩定性理論空間模式的增長的三維T－S波，其流向波數α=0.07873589-0.0008705923i，頻率ω＝0.02844068，展向波數β＝0.1706441，壁面用無滑移條件和定溫條件，溫度Tw=259.6K。圖3.1-1和圖3.1-2分別畫出了線性PSE和線性穩定性理論得到的流向波數的實部和虛部沿流向的變化，其中虛部的負值為擾動增長率。可見，在經過入口處附近的調整之后，線性PSE到x＝300處后已不再變化，線性穩定性理論所給之值的相對偏差為10－4量級。圖3.1-3至圖3.1-7畫出了x＝400處線性PSE得到的擾動形狀函數（下游結果與此處完全相同），以及線性穩定性理論得到的結果。二者是一致的。圖3.1-1流向波數實部αr圖3.1-2流向波數虛部αi圖3.1-3流向速度擾動剖面圖3.1-4法向速度擾動剖面圖3.1-5展向速度擾動剖面圖3.1-6密度擾動剖面圖3.1-7溫度擾動剖面3.1.1.2線性PSE在超音速平行流第一模態擾動演化中計算的檢驗基本流為平板二維Blasius相似性解，但將其局部平行化。馬赫數M＝4.5，自由來流溫度255.7K，普朗特數Pr＝0.707，入口雷諾數Re0＝10494；入口擾動為線性穩定性理論空間模式第一模態增長的三維T－S波，其流向波數α=0.06479422-0.0004128761i，頻率ω＝0.05351234，展向波數β＝0.02331896，壁面用無滑移條件和定溫條件，溫度Tw=1116.8K。圖3.1-8和圖3.1-9分別畫出了線性PSE和線性穩定性理論得到的流向波數的實部和虛部。線性PSE所給出的值在經過入口處附近的調整之后，到x＝100后不再變化，與線性穩定性理論所給值的相對偏差為10－4量級（如何看出來的偏差？）。圖3.1-10至圖3.1-14畫出了x＝100處線性PSE得到的擾動形狀函數（下游結果與此處完全相同），以及線性穩定性理論得到的結果。二者是一致的。圖3.1-8流向波數實部αr圖3.1-9流向波數虛部αi圖3.1-10流向速度擾動剖面圖3.1-11法向速度擾動剖面圖3.1-12展向速度擾動剖面圖3.1-13密度擾動剖面圖3.1-14溫度擾動剖面3.1.1.3線性PSE在超音速平行流第二模態擾動演化中計算的檢驗基本流為平板二維平行的Blasius相似性解，但將其局部平行化。馬赫數M＝4.5，自由來流溫度255.7K，普朗特數Pr＝0.707，入口雷諾數Re0＝11251；入口擾動為線性穩定性理論給出的空間模式第二模態增長的三維T－S波，其流向波數α=0.2740305-0.004373686i，頻率ω＝0.2518447，展向波數β＝0.04663792，壁面用無滑移條件和定溫條件，溫度Tw=1116.8K。圖3.1-15和圖3.1-16分別畫出了線性PSE和線性穩定性理論得到的流向波數的實部和虛部。可見，線性PSE所給流向波數的實部和虛部，在經過入口處附近的調整之后，到x＝1000后已不再變化，與線性穩定性所給之值的相對偏差為10－4量級。圖3.1-17至圖3.1-21畫出了x＝1000處線性PSE得到的擾動形狀函數（下游結果與此處完全相同），以及線性穩定性理論得到的結果。二者是一致的。圖3.1-15流向波數實部αr圖3.1-16流向波數虛部αi圖3.1-17流向速度擾動剖面圖3.1-18法向速度擾動剖面圖3.1-19展向速度擾動剖面圖3.1-20密度擾動剖面圖3.1-21溫度擾動剖面通過上面三個算例的檢驗可見，使用本文所用的入口擾動，迭代條件和差分格式，在平行的基本流下計算擾動演化，線性PSE得到的結果和由線性穩定性理論所得結果是可比的。線性PSE在平行流中的計算是可靠的。3.1.2線性PSE在非平行流中計算的檢驗把線性PSE方法用于非平行流中擾動演化的計算，并將其結果與線性穩定性理論的結果相比較。PSE沿流向要推進計算，LST在流向每個位置都要單獨計算一次，并將當地的結果與PSE的對比。對亞音速（馬赫數M＝0.3）邊界層和超音速（馬赫數M＝4.5）邊界層中的擾動演化結果都進行了檢驗。3.1.1.2線性PSE在亞音速非平行流中計算的檢驗基本流為平板二維非平行的Blasius相似性解，馬赫數M＝0.3，自由來流溫度255.7K，普朗特數Pr＝0.707，入口雷諾數Re0＝5119；入口擾動為空間模式線性穩定性理論給出的三維T－S波，其流向波數α=0.03518606+0.0004978373i，頻率ω＝0.006768882，展向波數β＝0.01137627，壁面用無滑移條件和定溫條件，溫度Tw=259.6K。圖3.1-22和圖3.1-23分別畫出了線性PSE和線性穩定性理論得到的流向波數的實部和虛部沿流向的變化。該三維T－S波沿流向經歷了由衰減到增長，再由增長到衰減的過程。兩種方法所得擾動流向波數實部αr和增長率-αi相符。圖3.1-24到圖3.1-31畫出了入口下游x＝10548，21097，31645，42193四個位置處的擾動剖面，兩種方法所得結果也相符。圖3.2-22波數實部αr沿流向變化圖3.1-23波數虛部αi沿流向變化圖3.1-24x＝10548處速度擾動剖面圖3.1-25x＝10548處溫度擾動剖面圖3.1-26x＝21097處速度擾動剖面圖3.1-27x＝21097處溫度擾動剖面圖3.1-28x＝31645處速度擾動剖面圖3.1-29x＝31645處溫度擾動剖面圖3.1-30x＝42193處速度擾動剖面圖3.1-31x＝42193處溫度擾動剖面3.1.2.2線性PSE在超音速非平行流第一模態擾動演化中計算的檢驗基本流為平板二維非平行的Blasius相似性解，馬赫數M＝4.5，自由來流溫度255.7K，普朗特數Pr＝0.707，入口雷諾數Re0＝11251；入口擾動為空間模式線性穩定性理論給出的第一模態三維T－S波，其流向波數α=0.01564881-0.00003612885i，頻率ω＝0.01223779，展向頻率β＝0.002331896，壁面用無滑移條件和定溫條件，溫度Tw=1116.8K。圖3.1-32和圖3.1-33分別畫出了線性PSE和線性穩定性理論得到的流向波數的實部和虛部沿流向的變化。該三維T－S波沿流向經過了增長率逐漸變大，再逐漸減小，直到衰減的過程。兩種方法所得結果非常接近。圖3.1-34到圖3.1-41畫出了入口下游x＝171534，343069，514603，686137四個位置處的擾動剖面。兩種方法所得結果吻合。圖3.1-32波數實部αr沿流向變化圖3.1-33波數虛部αi沿流向變化圖3.1-34x＝171534處速度擾動剖面圖3.1-35x＝171534處溫度擾動剖面圖3.1-36x＝343069處速度擾動剖面圖3.1-37x＝343069處溫度擾動剖面圖3.1-38x＝514603處速度擾動剖面圖3.1-39x＝514603處溫度擾動剖面圖3.1-40x＝686137處速度擾動剖面圖3.1-41x＝686137處溫度擾動剖面3.1.2.3線性PSE在超音速非平行流第二模態擾動演化中計算的檢驗基本流為平板二維非平行的Blasius相似性解，馬赫數M＝4.5，自由來流溫度255.7K，普朗特數Pr＝0.707，入口雷諾數Re0＝11251；入口擾動為空間模式線性穩定性理論給出的第二模態三維T－S波，其流向波數α=0.2189046+0.001522132i，頻率ω＝0.2170412，展向波數β＝0.04663792，壁面用無滑移條件和定溫條件，溫度Tw=1116.8K。圖3.1-42和圖3.1-43分別畫出了線性PSE和線性穩定性理論得到的流向波數的實部和虛部沿流向的變化。該三維T－S波沿流向經歷了由衰減到增長，再到衰減的過程。兩種方法所得結果基本吻合。圖3.1-44到圖3.2-51畫出了入口下游x＝2745，5489，8234，10978四個位置處的擾動剖面。兩種方法所得結果吻合。圖3.1-42波數實部αr沿流向變化圖3.1-43波數虛部αi沿流向變化圖3.1-44x＝2745處速度擾動剖面圖3.1-45x＝2745處溫度擾動剖面圖3.1-46x＝5489處速度擾動剖面圖3.1-47x＝5489處溫度擾動剖面圖3.1-48x＝8234處速度擾動剖面圖3.1-49x＝8234處溫度擾動剖面圖3.1-50x＝10978處速度擾動剖面圖3.1-51x＝10978處溫度擾動剖面通過上面三個算例的檢驗可見，無論是亞音速還是超音速邊界層，用線性PSE計算得到的擾動的波數，增長率和剖面都是可信的。線性PSE在非平行流中的計算是可靠的。3.2非線性PSE計算的檢驗非線性拋物化穩定性方程在不可壓縮邊界層中應用的檢驗，已經很充分了。Herbert和Bertolotti（1992），Esfahanian等（2001）將非線性PSE的計算結果與空間模式的直接數值模擬的結果做對比。他們發現二者計算結果幾乎完全一致，包括平均流修正、基本波、高次諧波的剖面形狀和幅值大小，如圖1-3所示。Malik等（1999）則將非線性PSE計算所得各階波的幅值與實驗結果對比，發現二者一致。可壓縮的非線性拋物化穩定性方程更為復雜，其在邊界層中應用的檢驗尚不充分。目前見到的工作主要來自NASA下屬Langley研究中心（NASALangleyResearchCenter）的Chang等人（1991）。他們用非線性的PSE計算了馬赫數4.5的邊界層中第二模態T—S波的演化，并將擾動的二次諧波的密度剖面與時間模式DNS的結果對比。結果表明兩種計算結果基本一致。Pruett和Chang（1995）、Pruett等人（1995）改用空間模式DNS與PSE比較，但是只比較了平均流修正，沒有比較二次諧波。要充分驗證非線性PSE在可壓縮邊界層中擾動演化問題中的應用，上述工作是不夠的。充分、完整的檢驗工作應該同時滿足下面4個要求：（1）應將非線性PSE的結果與空間模式，而不是時間模式DNS結果進行比較。（2）應該比較多種典型流動情況，包括不同的馬赫數和不同的擾動模態，而不僅限于第二模態擾動。（3）應該比較各階擾動，包括平均流修正、基本擾動和高次諧波，而不僅限于二次諧波。（4）對各階擾動，應該都比較其幅值大小和剖面形狀，而不僅是形狀。本文將選取三種典型流動情況，將非線性PSE計算結果與空間模式DNS結果進行對比，對比內容包括平均流修正、基本擾動和二次諧波的幅值大小和剖面形狀。另外，前面提到基本流可以取DNS計算得到的定常基本流，也可以取Blasius相似性解。這兩種基本流差別很小，但是用DNS計算擾動演化時為了避免數值振蕩，只能使用前者。本文在兩種基本流下都做了PSE計算，觀察結果是否有所差別。如果沒有差別，則以后用PSE計算時都可使用Blasius相似性解，從而避免了用DNS計算定常基本流所花的時間。3.2.1非線性PSE（即為NPSE）在亞音速邊界層中計算的檢驗3.2.1.1基本流為直接數值模擬所得的情況基本流為由直接數值模擬所得的定常層流，馬赫數M＝0.3，自由來流溫度300K，普朗特數Pr＝0.707，計算域入口處的雷諾數Re0＝1138；入口擾動為空間模式線性穩定性理論給出的增長的二維T－S波，其流向波數α=0.09498115-0.00424745i，頻率為ω＝0.02844742，擾動幅值為0.0001，壁面使用無滑移條件和等溫條件，溫度Tw=304.5K。非線性PSE和DNS計算時取同樣的網格，法向計算域長度為70.3，約合20倍邊界層名義厚度，法向共200個網格點，其中邊界層內不少于80個網格點，流向的計算步長為3.515，約為基本擾動波長的1/20。在用非線性PSE進行計算時，考慮到入口處加入的是二維擾動，所以只需考慮二維情況，式（2.1-40）和（2.1-41）中m取-7到7，n取0。T－S波對應的階數是（1，0），為基本擾動。由Fourier級數的知識可知，擾動在物理空間的幅值是0.0001的話，在Fourier展開譜空間的階數為（1，0）的波上的幅值就是0.0001/2，同時在階數為（-1，-0）的波上的幅值也是0.0001/2。階數為（2，0）的擾動就是二次諧波，階數為（3，0）的擾動就是三次諧波，以此類推，直至七次諧波，而階數為（0，0）的擾動就是平均流修正。非線性PSE計算可直接得到各階數的擾動，而DNS得到的則是各個空間點和各個時刻的物理量。因此，要將PSE的結果和DNS的結果相比，就要對DNS結果作譜分析。所用公式如下，，（3.2-1）其中u’為擾動速度，即瞬時速度減去基本流的速度，為階數為（m，n）的擾動的形狀函數，β為展向基本波數，ω為基頻。當然計算中譜分析也只需取有限多項，（3.2-2）其中M和N表示譜分析截斷的階數。圖3.2-1畫出了非線性PSE和DNS計算得到的入口下游x＝598，x＝773，x＝949，x＝1125四處的（0，0）階擾動剖面，也就是平均流修正，其中圖（a）畫的是平均速度修正剖面，圖（b）畫的是平均溫度修正剖面。非線性PSE和DNS計算得到的結果完全一致。圖3.2-2畫出了同樣位置處的（1，0）階，即基本擾動幅值剖面。為了便于與DNS的結果比較，PSE得到基本擾動的幅值都已經乘上了放大倍數N，即放大率的積分。由圖可見，非線性PSE和DNS計算得到的基本擾動完全一致。圖3.2-3畫出了同樣位置處的（2，0）階擾動，即二次諧波的幅值剖面。由圖可見，非線性PSE和DNS計算得到的二次諧波完全一致。這里，基本擾動的高次諧波只比較到了二次諧波，更高次的未作比較，這是因為直接數值模擬的分辨率不夠，無法準確分辨三次及三次以上的諧波。非線性拋物化穩定性方程只能計算有限次的諧波，但對每一個的諧波來說，其數值計算的分辨率都不成問題。問題是如果擾動幅值更大，則計算時需要將更高階的諧波包括進來。而對直接數值模擬而言，看起來似乎直接包含了很高次的諧波，但實際上，對高次諧波來說，其分辨率受到網格密度的制約。例如，本節中直接數值模擬在一個基本波長里約有20個網格點，在2次諧波的一個波長里就只有約10個網格點，做譜分析時，還勉強可以。但是對3次諧波來說，一個波長里只有約6個網格點，做譜分析（對空間數據進行格局、尺度分析的一種數學方法,運用傅里葉變換求取觀測數據分解產生正弦波及擬合最優波函數。）時，其精度就不夠了，更高次的諧波更是如此。如果想用直接數值模擬準確分辨更高次的諧波，就需要進一步加密計算網格，而這就要大大增加計算量。對于本算例來說，非線性PSE計算擾動演化結果與DNS所得完全一致，而前者所用時間比后者的小了約兩個數量級。（a）擾動速度剖面（b）擾動溫度剖面圖3.2-1入口下游x＝598，773，949，1125處（0，0）階擾動剖面（圖中四種線型表示PSE計算結果，四種標記符號表示DNS結果。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-2入口下游x＝598，773，949，1125處（1，0）階擾動幅值剖面（圖中四種線型表示PSE計算結果，四種標記符號表示DNS結果。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-3入口下游x＝598，773，949，1125處（2，0）階擾動幅值剖面（圖中四種線型表示PSE計算結果，四種標記符號表示DNS結果。）3.2.1.2基本流為Blasius相似性解的情況直接數值模擬計算的是完全的N－S方程，在入口加入擾動之前基本流需要計算至定常。而PSE則不同，它計算的是擾動量，選擇其它的基本流在計算上沒有什么困難，所以如果能有足夠近似的基本流，也能滿足精度要求。為此以下我們試用Blasius相似性解作為基本流，看其對非線性PSE計算結果的精度有何影響。圖3.2-4－3.2-7分別畫出了使用兩種基本流，通過非線性PSE計算得到的入口下游x＝598，x＝773，x＝949，x＝1125處的（0，0）、（1，0）、（2，0）和（7，0）階擾動的剖面。由圖可見，二者所得得各階擾動的幅值和剖面形狀完全一致。這就表明對于平板二維亞音速邊界層來說，用Blasius相似性解作為基本流，同樣可保證PSE計算結果的精度。（a）擾動速度剖面（b）擾動溫度剖面圖3.2-4入口下游x＝598，773，949，1125處PSE得到的（0，0）階擾動剖面（圖中四種線型表示用的是DNS計算的定常基本流，四種標記符號表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-5入口下游x＝598，773，949，1125處PSE得到的（1，0）階擾動幅值剖面（圖中四種線型表示用的是DNS計算的定常基本流，四種標記符號表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-6入口下游x＝598，773，949，1125處PSE得到的（2，0）階擾動幅值剖面（圖中四種線型表示用的是DNS計算的定常基本流，四種標記符號表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-7入口下游x＝598，773，949，1125處PSE得到的（7，0）階擾動幅值剖面（圖中四種線型表示用的是DNS計算的定常基本流，四種標記符號表示用的是Blasius相似性解。）3.2.2非線性PSE在超音速邊界層第一模態擾動演化中計算的檢驗3.2.2.1基本流為直接數值模擬所得的情況基本流為由直接數值模擬所得的定常層流，馬赫數M＝4.5，自由來流溫度255.7K，普朗特數Pr＝0.72，計算域入口處的雷諾數Re0＝10471；入口擾動為空間模式線性穩定性理論給出的第一模態二維T－S波，其流向波數α=0.0642482-0.0003461351i，頻率ω＝0.05333153，擾動幅值為0.01，壁面用無滑移條件和絕熱條件。非線性PSE和DNS計算時取同樣的網格，法向計算域長度為171.9，約合15倍邊界層名義厚度，共300個網格點，其中邊界層內不少于200個網格點，流向的計算步長為2.15，約為基本擾動波長的1/40。在用非線性PSE進行計算時，考慮到入口處加入的是二維擾動，所以只需考慮二維情況，式（2.1-40）和（2.1-41）中m取-7到7，n取0。T－S波對應的階數是（1，0）。圖3.2-8畫出了非線性PSE和DNS計算得到的入口下游x＝516處的（0，0）階擾動剖面，也就是平均流修正，其中圖（a）畫的是平均速度修正剖面，圖（b）畫的是平均溫度修正剖面。由圖可見，二者所給結果基本符合，其中溫度比速度符合得更好。圖3.2-9畫出了非線性PSE和DNS計算得到的入口下游x＝516處的（1，0）階，也就是基本擾動的幅值剖面。PSE得到基本擾動的幅值都已經乘上了放大倍數N。由圖可見，二者所給結果完全符合。圖3.2-10畫出了非線性PSE和DNS計算得到的入口下游x＝516處的（2，0）階擾動幅值剖面。由圖可見，二者所給結果基本符合，其中溫度比速度符合得更好。（a）擾動速度剖面（b）擾動溫度剖面圖3.2-8入口下游x＝516處（0，0）階擾動剖面（圖中實線表示PSE計算結果，帶有圓形標記的虛線表示DNS結果。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-9入口下游x＝516處（1，0）階擾動幅值剖面（圖中實線表示PSE計算結果，帶有圓形標記的虛線表示DNS結果。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-10入口下游x＝516處（2，0）（第一個表示頻率ω，第二個表示波數β二維模式）階擾動幅值剖面（圖中實線表示PSE計算結果，帶有圓形標記的虛線表示DNS結果。）黃章峰等（2005）和曹偉等（2006）的研究發現，在層流-湍流轉捩的breakdown過程中，層流剖面得以快速轉變為湍流剖面的機理在于擾動增長后，對平均剖面進行的修正導致了其穩定性特征的顯著變化，表現為其對應的中性曲線和不穩定波的增長率急劇增大。為了以后能夠將非線性PSE方法用于轉捩預測，有必要對PSE和DNS所得修正后的平均流剖面進行穩定性分析，并作比較。本文在入口下游x＝516將得到的平均流修正，即（0，0）階擾動，與當地的基本流相加，就得到了此處修正后的平均流修正剖面。對其進行線性穩定性分析，尋找其中性曲線。圖3.2-11（a）在（β，ω）平面上畫出了第一模態和第二模態擾動的中性曲線（只畫了β≥0的一半），圖3.2-11（b）畫出了第一模態和第二模態擾動的最大增長率隨展向波數β的變化（只畫了β≥0的一半）。可以看出，PSE和DNS計算得到的修正后的平均流剖面，其穩定性分析結果完全相符。（a）（b）圖3.2-11入口下游x＝516用PSE和DNS得到的修正后的平均流剖面的穩定性分析結果：（a）（β，ω）平面上三維擾動的中性曲線；（b）三維擾動的最大增長率隨展向波數β的變化。對于本算例來說，非線性PSE計算擾動演化結果與DNS所得基本符合，而前者所用的計算量遠遠少于后者的。3.2.2.2基本流為Blasius相似性解的情況上面用的基本流是直接數值模擬計算得到的定常基本流，本文還用Blasius相似性解作為基本流，使用非線性PSE計算擾動演化。圖3.2-12－3.2-15分別畫出了使用兩種基本流，通過非線性PSE計算得到的入口下游x＝516處的（0，0）、（1，0）、（2，0）和（7，0）階擾動幅值剖面。由圖可見，基本流改用Blasius相似性解之后，非線性PSE所得結果不變。而用Blasius相似性解作為基本流，比用直接數值方法求基本流，可以節省不少時間。（a）擾動速度剖面（b）擾動溫度剖面圖3.2-12入口下游x＝516處PSE得到的（0，0）階擾動剖面（圖中實線表示用的是DNS計算的定常基本流，圓形標記表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-13入口下游x＝516處PSE得到的（1，0）階擾動幅值剖面（圖中實線表示用的是DNS計算的定常基本流，圓形標記表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-14入口下游x＝516處PSE得到的（2，0）階擾動幅值剖面（圖中實線表示用的是DNS計算的定常基本流，圓形標記表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-15入口下游x＝516處PSE得到的（7，0）階擾動幅值剖面（圖中實線表示用的是DNS計算的定常基本流，圓形標記表示用的是Blasius相似性解。）3.2.3非線性PSE在超音速邊界層第二模態擾動演化中計算的檢驗3.2.3.1基本流為直接數值模擬所得的情況基本流為由直接數值模擬所得的定常層流，馬赫數M＝4.5，自由來流溫度255.7K，普朗特數Pr＝0.72，計算域入口處的雷諾數Re0＝11635；入口擾動為空間模式線性穩定性理論給出的第二模態二維T－S波，其流向波數α=0.2826772-0.005469101i，頻率ω＝0.2580685，擾動幅值為0.001，壁面用無滑移條件和絕熱條件。非線性PSE和DNS計算時取同樣的網格，法向計算域長度為171.9，約合15倍邊界層名義厚度，法向共256個網格點，其中邊界層內不少于130個網格點，流向的計算步長為0.8595，約為基本擾動波長的1/25。在用非線性PSE進行計算時，式（2.1-27）和（2.1-28）中m取-7到7，n取0，T－S波對應的階數為（1，0）。圖3.2-16畫出了非線性PSE和DNS計算得到的入口下游x＝138處的（0，0）階擾動剖面，也就是平均流修正，其中圖（a）畫的是平均速度修正剖面，圖（b）畫的是平均溫度修正剖面。由圖可見，非線性PSE和DNS計算得到的平均流修正基本符合，其中溫度比速度符合得更好。圖3.2-17畫出了非線性PSE和DNS計算得到的入口下游x＝138處的（1，0）階，也就是基本擾動幅值剖面。PSE得到基本擾動的幅值都已經乘上了放大倍數N。由圖可見，二者所給結果完全吻合。圖3.2-18畫出了非線性PSE和DNS計算得到的入口下游x＝138處的（2，0）階擾動幅值剖面。由圖可見，二者所給結果基本吻合，其中溫度比速度符合得更好。（a）擾動速度剖面（b）擾動溫度剖面圖3.2-16入口下游x＝138處（0，0）階擾動剖面（圖中實線表示PSE計算結果，帶有圓形標記的虛線表示DNS結果。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-17入口下游x＝138處（1，0）階擾動幅值剖面（圖中實線表示PSE計算結果，帶有圓形標記的虛線表示DNS結果。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-18入口下游x＝138處（2，0）階擾動幅值剖面（圖中實線表示PSE計算結果，帶有圓形標記的虛線表示DNS結果。）在入口下游x＝138將得到的平均流修正，即（0，0）階擾動，與當地的基本流相加，就得到了此處修正后的平均流修正剖面，然后對其進行線性穩定性分析。圖3.2-19（a）在（β，ω）平面上畫出了第一模態和第二模態擾動的中性曲線，圖3.2-19（b）畫出了第一模態和第二模態擾動的最大增長率隨展向波數β的變化。可以看出，PSE和DNS計算得到的修正后的平均流剖面，其穩定性分析結果完全相符。（a）（b）圖3.2-19入口下游x＝138用PSE和DNS得到的修正后的平均流剖面的穩定性分析結果：（a）（β，ω）平面上三維擾動的中性曲線；（b）三維擾動的最大增長率隨展向波數β的變化。對于本算例來說，非線性PSE計算擾動演化所用時間比直接數值模擬的小約兩個數量級。且本算例中DNS目前無法用單機計算，所以只能使用多個機器并行計算，而非線性PSE則使用單機計算。對于本算例來說，非線性PSE計算擾動演化結果與DNS所得基本符合，而前者所用的計算量遠遠少于后者的。3.2.3.2基本流為Blasius相似性解的情況下面再用Blasius相似性解作為基本流，使用非線性PSE計算擾動演化。圖3.2-20－3.2-23畫出了使用兩種基本流所得的入口下游x＝138處的（0，0）、（1，0）、（2，0）和（7，0）階擾動剖面。由圖可見，基本流改用Blasius相似性解之后，非線性PSE所得結果——包括各個階數的擾動幅值大小和剖面形狀——完全不變。這就表明對于平板二維超音速邊界層來說，研究其中的第二模態擾動演化時也可以選用Blasius相似性解作為基本流，這樣可以節省計算定常基本流的時間。（a）擾動速度剖面（b）擾動溫度剖面圖3.2-20入口下游x＝138處PSE得到的（0，0）階擾動剖面（圖中實線表示用的是DNS計算的定常基本流，圓形標記表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-21入口下游x＝138處PSE得到的（1，0）階擾動幅值剖面（圖中實線表示用的是DNS計算的定常基本流，圓形標記表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-22入口下游x＝138處PSE得到的（2，0）階擾動幅值剖面（圖中實線表示用的是DNS計算的定常基本流，圓形標記表示用的是Blasius相似性解。）（a）擾動速度幅值剖面（b）擾動溫度幅值剖面圖3.2-23入口下游x＝138處PSE得到的（7，0）階擾動幅值剖面（圖中實線表示用的是DNS計算的定常基本流，圓形標記表示用的是Blasius相似性解。）通過上面三個算例的檢驗可見，無論是亞音速還是超音速邊界層，用非線性PSE計算得到的擾動結果可信的。非線性PSE的計算是可靠的。另外，基本流可以使用Blasius相似性解，所得結果足夠準確。綜上所述，線性和非線性PSE的計算都得到了充分檢驗，其結果可靠。下面可以用PSE作為工具研究某些具體問題。

4用線性PSE研究非平行性對可壓縮邊界層中性曲線的影響線性穩定性理論假設流動是平行的。嚴格的平行流動，在自然界極少見，在工程技術中也只在有固壁限制的情況下才能見到。非平行流動則是常見的，例如邊界層流、自由剪切層等。只有在雷諾數很大時，流線才近似于平行。而對不可壓平板邊界層，臨界雷諾數不是很大，非平行性不可忽略，因而由線性穩定性理論在平行流假定下給出的臨界雷諾數大于實驗結果。為此，Gaster（1974）、于秀陽和周恒（1986）等研究了在線性穩定性理論的框架內如何改進結果的問題。后來PSE方法出來后，用以研究這一問題，得到了滿意的結果。對可壓縮邊界層，顯然有同樣問題。本文將用線性拋物化穩定性方程和線性穩定性理論分別求臨界雷諾數附近的中性曲線，并進行對比，以考察非平行性對中性曲線的影響。4.1非平行性對亞音速邊界層中性曲線的影響基本流為平板二維Blasius相似性解，馬赫數M＝0.3，自由來流溫度255.7K，普朗特數Pr＝0.72，壁面用無滑移條件和定溫條件（看粘流無滑移條件就是指流體和固壁之間沒有相對滑動,比如水流過管道的時候,無滑移條件就是指在靠近管壁面處水的流動速度為零.有滑移條件就是他們之間有相對速度網上查的），溫度Tw=259.6K，從雷諾數Re0＝170.6處開始找中性曲線。在入口加入衰減的二維T－S波，其流向波數為α=0.1159899+0.009840206i，頻率為ω＝0.04550159。圖4.1-1顯示出了線性PSE計算得到的該頻率擾動的放大率（）沿流向的變化情況。放大率在x<x1和x>x2處小于0，在x1<x<x2處大于0，x＝x1和x＝x2處等于0，其中，x1＝252.344，x2＝412.909，兩處對應的當地雷諾數（）分別為93729.88和471189.3。它們在頻率、雷諾數平面上代表的就是中性曲線上的兩個點，其中第一個點在中性曲線第一分支上，第二個點在中性曲線第二分支上。在入口加入衰減的二維T－S波，其流向波數為α=0.1604609+0.006466334i，頻率為ω＝0.06825238。圖4.1-2顯示出了線性PSE計算得到的該頻率擾動的放大率（）沿流向的變化情況。可見，放大率沿流向始終小于0，這表明該頻率的波在中性曲線以外。將不同頻率擾動加在入口，并將得到的中性點都在頻率、雷諾數平面上標出，就得到中性曲線，如圖4.1-3中實線所示（其中無量綱頻率，為有量綱頻率）。圖4.1-3中虛線為用線性穩定性理論（是什么理論）得到的中性曲線。對比兩種方法得到的中性曲線可以看出，二者在小雷諾數處存在差別，隨著雷諾數的增大，二者逐漸趨于一致，在大雷諾數處則完全重合。用線性拋物化穩定性方程得到（如何得到）的臨界雷諾數要小于線性穩定性理論的結果，而在臨界雷諾數處，用線性PSE所得的不穩定擾動的頻率范圍更大。大雷諾數處兩條曲線完全一致，這表明此處非平行性對邊界層中性曲線的影響很小。這是因為平板邊界層厚度與到前緣距離的1/2次方成正比（δ～x1/2），大雷諾數處邊界層厚度沿流向的變化就很小（dδ/dx～x-1/2，當x很大時dδ/dx趨近于0），因而流動接近于平行。上述結論與人們用實驗、理論和計算的方法（Schubauer和Skramstad19431947，Gaster1974，于秀陽和周恒1986，Herbert1997，Bertolotti等1992，Ross等1970，Strazisar1977，Kachanov等1977）研究非平行性對不可壓縮邊界層中性曲線的影響時所得結論相同（見圖1-2）。圖4.1-1某一頻率（ω＝0.04550159）擾動的放大率沿流向的變化圖4.1-2某一頻率（ω＝0.06825238）擾動的放大率沿流向的變化圖4.1-3馬赫數為0.3的邊界層中二維擾動的中性曲線（，）4.2非平行性對超音速邊界層中性曲線的影響基本流為平板二維Blasius相似性解，馬赫數M＝4.5，自由來流溫度255.7K，普朗特數Pr＝0.72，壁面用無滑移條件和定溫條件，溫度Tw=1125.2K，開始計算處的雷諾數Re0＝349。求第二模態T－S波的中性曲線（如何求？）。在入口加入衰減的二維T－S波（為什么這么加？衰減是由于e的指數決定的指數大于零時是增長的，小于時是衰減的，即流向波數的虛部），其流向波數α=0.216862+0.003333925i，頻率ω＝0.1977932。圖4.2-1顯示出了線性PSE計算得到的該頻率擾動的放大率（）沿流向的變化情況。和上一節相似，可見，放大率在x<x1和x>x2處小于0，在x1<x<x2處大于0，而x1＝111.7，x2＝199.6，對應的當地雷諾數（）分別為160888和191594。因而可得中性曲線上得兩個點，分別位于中性曲線的第一分支和第二分支上。若加入另一個二維T－S波，其流向波數α=0.24118969+0.001540257i，頻率ω＝0.221063，則圖4.2-2顯示出該擾動波的放大率（）沿流向始終小于0，這表明該頻率的擾動波在頻率、雷諾數平面上處在中性曲線以外。將不同頻率擾動加在入口，并將得到中性點都在頻率、雷諾數平面上表示出來，就得到了馬赫數4.5第二模態T－S波的中性曲線，如圖4.2-3中實線所示。圖4.2-3中虛線為用線性穩定性理論得到的中性曲線。對比兩種方法得到的中性曲線可以看出，結果與亞音速情況類似，兩條中性曲線在小雷諾數處存在差別，隨著雷諾數的增大，二者逐漸趨于一致，在大雷諾數處則完全重合。線性拋物化穩定性方程得到的臨界雷諾數要小于線性穩定性理論的結果，在臨界雷諾數處，線性PSE得到的增長擾動的頻率范圍更大。圖4.2-1某一頻率（ω＝0.1977932）擾動的放大率沿流向的變化圖4.2-2某一頻率（ω＝0.221063）擾動的放大率沿流向的變化圖4.2-3馬赫數為4.5的邊界層中第二模態二維T－S波的中性曲線（，）4.3小結本文使用線性拋物化穩定性方程尋找可壓縮邊界層中二維小擾動的中性曲線，并與空間模式線性穩定性理論所得進行比較。結果顯示，無論是在亞音速還是超音速邊界層中，非平行性的影響在臨界雷諾數處比較明顯，具體有兩個方面：一個是使得臨界雷諾數變小；另一個是使得臨界雷諾數處增長擾動的頻率范圍變得更大。在大雷諾數處，由于流動近似平行，非平行性對邊界層中性曲線的影響很弱。上述結論與人們研究非平行性對不可壓縮邊界層中性曲線的影響時所得結論相同，但目前尚未見到可壓縮流的實驗結果。

5非線性PSE在超音速邊界層二次失穩問題中的應用層流到湍流的轉捩一直是人們感興趣的問題。對于不可壓縮流來說，轉捩可以由二維小擾動的線性放大開始。當二維擾動幅值大到一定程度，非線性作用就會激發出三維擾動，它開始對轉捩起重要作用并最終導致流動變為湍流。Herbert（1983a，1983b，1984）所引入的二次失穩機制，就是三維擾動得以被激發的可能原因之一。用流動顯示技術發現的Λ結構，從實驗（Saric等1984）上證實了該機制的存在。此外，實驗上還得到了三維亞諧波的增長率隨其展向波數的變化關系（Thomas1987），理論（Herbert1988）和計算（Spalart和Yang1986）的結果都與之一致。我們關心的問題是這一機制對超音速邊界層是否起作用。目前，這方面的工作還很少。Maslov所做的實驗是其中之一，他在轉捩過程中發現了亞諧波，但這只是用熱線風速儀測量到了亞諧信號，很難確定它就是由二次失穩產生的。另一個就是董亞妮和周恒（2006）所作的數值研究。他們通過直接數值模擬證實了當基本擾動是第二模態T-S波時二次失穩機制存在，并且找到了亞諧波的增長率隨其展向波數和二維基本波幅值的變化關系。但是DNS耗時太多，研究的算例很有限，所以董亞妮和周恒沒有給出能產生不穩定亞諧波的基本波的門限值，也沒有研