可壓縮混合層的渦結構演化與流質(zhì)混合

周 強，何 楓，沈孟育

（清華大學 航天航空學院，北京 100084）

0 引 言

可壓縮平面自由剪切層，是指兩股速度不等的平行氣流在自由空間相互混合，其初始混合面為平面。這種流動是可壓縮湍流研究的典型問題，也是超燃沖壓發(fā)動機內(nèi)基本流動的結構。超聲速平面混合層厚度增長率遠低于亞聲速混合層和不可壓混合層的增長率，所以提高超聲速平面混合層增長率、增強燃料混合是超燃沖壓發(fā)動機需要解決的關鍵技術之一。

Barre［2］等利用紋影技術對超聲速混合層進行了實驗研究，發(fā)現(xiàn)當對流馬赫數(shù)Mc＝0.62時，超聲速混合層的發(fā)展對初始條件十分敏感，但充分發(fā)展段湍流的這種依賴性相對較弱，湍流應力的大小和各向異性特征受到可壓縮性的影響。Clemens ＆ Mungal［3］利用紋影技術、平面激光誘導熒光技術和米散射技術研究了Mc為0.28、0.62和0.79的超聲速平面混合層中大尺度渦結構及其卷吸能力，給出了不同對流馬赫數(shù)下大尺度結構的空間存在形狀。熊紅亮、崔爾杰［4］利用紋影成像技術測出了可壓縮混合層的厚度增長率曲線，紋影照片顯示隨著對流馬赫數(shù)增加，可壓縮混合層厚度增長率顯著降低，渦卷起現(xiàn)象延遲，結果還表明適當增加分隔板厚度有利于平板尾緣可壓縮混合層中的渦提前卷起，有利于促進混合。Zhao［5］等使用高分辨率的NPLS方法對超聲速混合層進行了流動顯示并進行了湍流分形結構研究，結果顯示在轉(zhuǎn)捩區(qū)域分形維數(shù)隨著湍流強度的增強而增長，在自相似區(qū)域分形維數(shù)基本維持不變。

由于超聲速實驗觀測的難度，定量實驗信息仍然非常有限，在分析解讀所獲取的實驗結果上，缺乏足夠的理論依據(jù)，使得現(xiàn)有的實驗結果所能展示的定量流動信息非常有限，數(shù)值模擬便成為重要的研究手段之一。

在早期的研究中，數(shù)值模擬一般都是針對時間發(fā)展混合層［6－10］。在低對流馬赫數(shù)的情況下時間發(fā)展的結果是真實的三維空間發(fā)展混合層結果的良好近似，但在高對流馬赫數(shù)時二者存在有較大差別［11－12］。因此高對流馬赫數(shù)三維空間發(fā)展混合層的精細結構的數(shù)值模擬就顯得特別重要?？蓧嚎s混合層空間發(fā)展模式的數(shù)值模擬，無論是直接數(shù)值模擬（DNS）還是大渦模擬（LES）都得到了混合層發(fā)展早期階段的信息［13－17］，例如都捕捉到初期的Λ渦和馬蹄渦的出現(xiàn)，隨對流馬赫數(shù)的增高動量厚度沿空間增長發(fā)展變緩慢等等物理現(xiàn)象。但對混合層渦結構的中后期湍流發(fā)展的研究不多，關注不同尺度的渦結構對流質(zhì)混合作用的研究也不多，還缺乏足夠全面的信息來加深對可壓縮湍流的認識，這在很大程度上是因為受到計算規(guī)?；蚪?jīng)費的限制，而這方面的研究對認識可壓縮混合層發(fā)展的物理機制，以及提出新的增混舉措或思路有著重要的意義。

本文采用約4億規(guī)模的網(wǎng)格以及作者自己構造的七階精度廣義緊致格式離散對流項，用顯式八階精度的中心格式離散粘性項，數(shù)值求解了非定常三維可壓縮Navier－Stokes（NS）方程，直接數(shù)值模擬了對流馬赫數(shù)為0.7的超聲速可壓縮混合層的空間發(fā)展流動，獲得了自初始流動失穩(wěn)直至充分發(fā)展湍流流動結構的精細演化歷程，并著重于研究結構的演化對混合層中流質(zhì)混合的作用。

1 可壓縮混合層直接數(shù)值模擬方法簡述

以高速側(cè)來流參數(shù)作為相應的特征參數(shù)和進口處無擾動速度剖面的渦量厚度為特征長度，對非定常三維可壓縮NS方程及其定解條件和完全氣體狀態(tài)方程進行無量綱化。無量綱化后的NS方程組可參見文獻［6］。

對流項先經(jīng)Steger－Warming矢通量分裂技術［18］分裂成非正和非負兩部分，然后對正負通量均采用作者提出的七階廣義緊致格式進行離散［19］。該緊致格式可表達如下：

其中，u為求導函數(shù)（此處即為分裂后的通量），F(xiàn)為一階導數(shù)，S為二階導數(shù)，h為網(wǎng)格間距。格式系數(shù)am，n和的具體數(shù)值以及非周期邊界的匹配方法可參見文獻［19］。值得注意的是，格式中出現(xiàn)了二階導數(shù)。Fourier分析表明二階導數(shù)的存在使得一階導數(shù)獲得了更好的分辨率性質(zhì)［19］。

在我國西南地區(qū)的深山野林里，生長著一種草本植物，它開著白色的小花，雖然看上去清新可愛，卻也沒什么特別之處。但是，如果下雨之后你再來看，便會發(fā)現(xiàn)原先白色的小花變得晶瑩剔透，就像精心雕琢后的精美藝術品。你知道這種神奇植物的名字嗎？

對NS方程中的粘性項，本文采用普通顯式八階精度的中心格式進行離散。時間積分采用八步四階具有強穩(wěn)定保持性質(zhì)的Runge－Kutta格式［20］。同時求解的還有無量綱被動標量質(zhì)量分數(shù)輸運方程，該方程的具體表達可參見文獻［6］。

初始流向速度場采用上下來流定解條件下的可壓縮邊界層方程的數(shù)值解［21］，橫向速度假設為0；初始溫度剖面由Crocco－Busemann關系求得，普朗特數(shù)設為0.75，混合層上下側(cè)密度相同。

入口邊界在可壓縮邊界層的數(shù)值解基礎上添加擾動。若加入隨機擾動，則流動要經(jīng)歷很長的空間線性穩(wěn)定發(fā)展后，最不穩(wěn)定波才成為主導，隨后流動呈非線性發(fā)展，這將耗費很多計算時間。為減少DNS的計算成本，使流動較快進入轉(zhuǎn)捩階段，本文在入口處直接添加最不穩(wěn)定波形式的擾動［21］。擾動波由一對振幅相同方向相反的最不穩(wěn)定斜波組成，該斜波的具體形式可參見文獻［21－22］。

在流向超聲速出口邊界處無需設定邊界條件，其上流動參數(shù)直接由內(nèi)點外插得到；展向采用周期邊界條件；法向上下兩側(cè)均采用無反射邊界條件［6］。

為顯示轉(zhuǎn)捩后湍流流場的空間演化，直到湍流自相似區(qū)，計算域流向（以x表示）必須足夠長；混合層初始混合平面以y＝0表示，故y為流動的法線方向；計算域的展向（以z表示）的寬度應根據(jù)湍流脈動相關量是否能夠在展向半個寬度內(nèi)衰減為0來設定。

本文中空間發(fā)展混合層的DNS計算域選為長方形［0，L1］×［－100，100］×［0，L2］，其中L1、L2分別為流向x和展向z混合層線性穩(wěn)定性分析下各自最不穩(wěn)定波波長的54倍和2倍。參照Pantano ＆Sarkar［23］和 Wu ＆ Moin［24］所做 DNS的經(jīng)驗設定網(wǎng)格密度，網(wǎng)格數(shù)為4160×351×256。所有計算在中科院計算機網(wǎng)絡信息中心進行。

2 可壓縮混合層的渦結構演化與流質(zhì)混合

本文中DNS的對象為：對流馬赫數(shù)為0.7，上下來流馬赫數(shù)分別為2.8和1.4，普朗特數(shù)為0.75，以入口處上下自由來流速度差和動量厚度為特征量定義的雷諾數(shù)為80。Ho ＆Huerre［25］研究表明上述定義的雷諾數(shù)大于50時，不穩(wěn)定波的發(fā)展和渦結構的演化基本不受雷諾數(shù)的影響。

計算保證不同流向x處的平均流向速度剖面沿法向的分布進入自相似區(qū)，流向脈動速度的展向無量綱能譜中高波數(shù)對應的小尺度（Kolmogorov尺度）湍流的粘性子區(qū)的能譜，呈現(xiàn)出Heisenberg［26］通過能量級串理論給出的－7次方律，說明DNS的計算網(wǎng)格密度足以解析粘性子區(qū)的小尺度流動細節(jié)。通過監(jiān)測湍動能輸運方程殘量的大小，發(fā)現(xiàn)數(shù)值殘量相對各輸運項數(shù)值大小是可以忽略的，即數(shù)值耗散遠小于物理耗散，表明物理粘性在計算結果中可得到正確的反映。關于本文計算結果的可靠性驗證的詳細闡述，見文獻［22］。

展現(xiàn)渦結構的形態(tài)有多種方法，本文采用旋轉(zhuǎn)強度（速度梯度的共軛復特征值的虛部的平方）為0.025等值面顯示不同尺度的渦結構形態(tài)。

可壓縮混合層三維結構從早期不穩(wěn)定波到渦結構破碎的演化全過程已在文獻［22］中有詳細描述。為了研究渦結構演化對流質(zhì)混合的作用，本文對渦結構演化過程將做簡單概述。

圖1（a）和圖1（b）顯示了不同視角下渦結構展現(xiàn)的形態(tài)，對應流向坐標為x≈0～250，圖中從藍到紅不同色度反映的是結構在空間法向相對位置的由低到高，由此可以清晰地展示結構在法向方向上的擴張演化。

圖1 流動結構演化（旋轉(zhuǎn)強度為0.025等值面）Fig.1 Vortex structures（isosurfaces of swirl strength with value of 0.025）

在入口處直接加入較強幅度的最不穩(wěn)定波擾動后，混合層離入口僅一個基頻波長就形成了強剪切層（x≈0～11），之后在混合層上下迅速卷起Λ斜波（x≈11～45）。在Λ斜波交接處不連續(xù)，生長著類似帽子的薄殼結構，在剪切作用下這個交接處開始向混合層外側(cè)伸展，并演化為馬蹄渦的形態(tài)（x≈50～80），生長在Λ斜波的兩腿上，而仔細觀察Λ斜波結構上開始附著稍小的結構。馬蹄渦頭部在剪切力的作用下開始失去穩(wěn)定，演變成頸部有半環(huán)形項鏈渦的結構，在位于入口的8倍基頻波長處（x≈88～140）形成復雜的團花形態(tài)，而斜波的形態(tài)已經(jīng)不再清晰了；在12倍流向基頻波長處，團花結構被拉伸并開始破碎，出現(xiàn)大量的各種方位的細管狀結構；約在x≈450之后流動進入湍流自相似區(qū)，充斥著更加細碎的小渦結構，并且這些細碎的小結構管條狀結構整體以一種大的擬序結構向下游遷移。關于渦結構演化的三維視圖和湍流自相似區(qū)的小渦形態(tài)可參看文獻［22］。

圖1（a）至圖1（b）還清晰地展現(xiàn)了混合層上下側(cè)不對稱的發(fā)展演化。下側(cè)渦結構演化略早于上側(cè)，下側(cè)的渦花結構產(chǎn)生較早，在較強的剪切下很快破碎，而上側(cè)的渦花結構維持過程相對長些。不對稱發(fā)展的根本原因在于渦結構的抬起使得渦結構的行進速度受到了主流的影響。關于該點的詳細解釋見參考文獻［22］的圖9以及與圖9相關的論述。結構這種上下不對稱的發(fā)展使得流質(zhì)的混合也向下側(cè)偏移發(fā)展。

式（2）為沿流向動量厚度計算公式［27］，其中ˉρ為展向平均的雷諾時均統(tǒng)計值；ρ1為初始來流高速側(cè)密度。采用式（2），可以得動量厚度沿流向的變化曲線（文獻［22］中的圖14b）。

從動量厚度曲線可以看出在流向初始階段（x＜200），動量厚度呈現(xiàn)快速增長階段，如果對應圖1可知，這一階段是從Λ渦發(fā)展到團花渦結構豐富的演化過程階段；x＞200團花失穩(wěn)破碎，流動進入轉(zhuǎn)捩。在進入充分發(fā)展湍流階段，動量厚度的曲線趨于線性增長時期，轉(zhuǎn)捩后的動量厚度增長率基本不變。x＜200是大尺度結構演化階段，是在主流的強剪切作用下、混合層快速吸收主流動量和能量的發(fā)展階段，這個過程越早發(fā)生或者延長這個階段的生存期，會促使混合層動量厚度曲線在x＞200的后階段向上平移，雖然在充分發(fā)展階段的沿流向的增長率未必提高，但最終動量厚度的絕對值會提高。

在混合層發(fā)展的不同空間階段，對應不同尺度結構，這些形態(tài)各異的結構對流質(zhì)混合的作用也是不同的。圖2～圖7展示了在不同渦結構的旋轉(zhuǎn)卷吸作用下，由被動標量質(zhì)量分數(shù)顯示的混合效果，白色為0是下側(cè)來流初始質(zhì)量分數(shù)，黑色為1是上側(cè)來流初始質(zhì)量分數(shù)，同時也給出了相應的渦量云圖。

圖2（a）為混合層初始發(fā)展起來的Λ渦結構，如果對流場信息做平行流動方向的投影z＝8.83鉛垂切面A，該投影面恰好穿過Λ渦頭部，如圖2（a）所示。圖2（b）和圖2（c）為該切面上的展向渦量分布云圖和質(zhì)量分數(shù)分布云圖。Λ渦頭部幾乎是展向方向的強旋轉(zhuǎn)促進了流體混合，并使混合層界面沿流向產(chǎn)生很大的波動，由于混合層法向上下各發(fā)展的Λ渦略有所不同，在y＜0低速側(cè)的Λ渦頭部具有更集中和強度更大的負展向渦量。

圖2 Λ渦在展向z＝8.83切面上的展向渦量云圖和被動標量云圖Fig.2 Vorticity contours and passive scalar contours corresponding to vortexΛat the plane z＝8.83

圖3 Λ渦在展向z＝4.43切面上的展向渦量云圖和被動標量云圖Fig.3 Vorticity contours and passive scalar contours corresponding toΛvortex at the plane z＝4.43

圖4 Λ渦在流向x＝25.3切面上的流向渦量云圖和被動標量云圖Fig.4 Vorticity contours and passive scalar contours corresponding toΛvortex at the plane x＝25.3

圖5 馬蹄渦作用下的混合被動標量云圖Fig.5 Passive scalar contours corresponding to horseshoe shaped vortex at the plane z＝8.83、17.0respectively

圖6 團花結構對應的混合被動標量云圖Fig.6 Passive scalar contours corresponding to flower－ball shaped vortex at the plane z＝4.43、13respectively

圖7 充分發(fā)展小尺度湍流結構作用下的混合被動標量云圖，z＝8.83和z＝13切面Fig.7 Passive scalar contours corresponding to full development flow at the plane z＝8.83、13respectively

如果對流場信息做平行流動方向的投影z＝4.43鉛垂切面B，如圖3（a）所示，切面錯過了Λ渦頭部，由圖3（b）和圖3（c）可以看到該切面上的展向渦量分布云圖和質(zhì)量分數(shù)分布云圖與圖2（b）、圖2（c）是不同的。圖3（b）中Λ渦腿之間的渦量，使得切面上呈現(xiàn)出類似不可壓混合層的渦卷序列，而實際是不同空間位置的Λ渦腿誘導出的展向渦量分布。從圖3（c）被動標量質(zhì)量分數(shù)分布可以看出，其渦腿對上下流質(zhì)混合是有一定的卷吸作用，但相對Λ渦交接的頭部而言（圖2c），圖3（c）對應的混合層直觀可視厚度還是相對小的。結合圖2（c）和圖3（c）表明混合層沿展向z的混合效果是不相同不均勻的，混合接觸界面是三維波動形式。

如果做圖4（a）所示的x＝25.3切面C，該切面是切在了同一個渦的雙腿間。可以看出，在主流強剪切作用下，渦的雙腿被拉伸并變成扁橢圓形，特別是這個渦渦腿之間的相對旋轉(zhuǎn)造成渦腿間的射流，主要使得混合層低速側(cè)的流質(zhì)向高速側(cè)噴注（或高速側(cè)向低速側(cè)噴注），相對旋轉(zhuǎn)的兩渦腿內(nèi)混合均勻。如圖4（c）所示混合層界面產(chǎn)生間歇性的大波動，為后續(xù)混合發(fā)展奠定了基礎。

隨著Λ渦頭部向混合層外側(cè)伸展，渦的形態(tài)演化為馬蹄渦。圖5顯示了主要演化出的馬蹄渦和展向z＝8.83和z＝17.0切面上的被動標量的云圖。馬蹄渦的寬頭部的強旋轉(zhuǎn)使得對應的混合層界面的局部區(qū)域可視混合厚度相對較大。這時混合層的界面開始破碎形成更多的接觸面，上下流體開始在渦的作用下進行混合。

圖6是馬蹄渦結構頭部演化成團花結構時期在不同的展向切面上對應的混合狀態(tài)被動標量云圖，以及如圖7所示進入充分發(fā)展湍流階段不同展向切面上的混合狀態(tài)?？梢钥吹皆谕牧鞯男〕叨冉Y構作用下混合趨于均勻并沿展向混合狀態(tài)基本相同，這個時期對應x＞200，結合動量厚度曲線（文獻［22］的圖14（b））可知混合層厚度沿流向的增長變得很緩慢。

如果將被動標量0＜Sc＜1顯示的質(zhì)量分數(shù)的混合，看成是定性的可視混合層厚度，由圖8不同展向切面上的可視混合可以看出，混合層向低速側(cè)偏移發(fā)展，與圖1（b）顯示的低速側(cè)的大結構發(fā)展先于高速側(cè)相呼應。最終x＞450充分發(fā)展湍流階段的混合基本上是在x＝100～200間的可視混合厚度量級基礎上進行的，這個階段可視混合在法向上的擴展并不大，小尺度的湍流結構混合耗散的使得流質(zhì)混合趨于均勻，由此也可見先期大尺度結構的發(fā)展使得混合層從主流獲取動量和能量，主宰著混合層可視厚度的發(fā)展量級，而充分發(fā)展階段對混合的主要貢獻是使這個厚度量級內(nèi)的流質(zhì)混合均勻化。

圖8 不同z展向切面上的可視混合效果Fig.8 Visual mixing layer shown by passive scalar contours at the z－plane

3 結 論

采用DNS的方法研究對流馬赫數(shù)為0.7的可壓縮混合層問題，得到了自混合層失穩(wěn)渦開始卷起直至充分發(fā)展湍流的整個流動結構演化的歷程?；旌蠈釉谧畈环€(wěn)定波主導下開始失穩(wěn)，旋轉(zhuǎn)強度等值面顯示下的流動結構開始具有初始的渦卷起結構，隨后形成Λ渦。Λ渦渦腿的交接處會在剪切作用下發(fā)展為馬蹄渦，馬蹄渦是生長在Λ渦渦腿上的，并且法向方向生長沿流向傾斜演化。馬蹄渦又在剪切作用下其頭部演化更加復雜的團花形態(tài)，而Λ渦渦腿形態(tài)上開始包含著稍小尺度的結構。團花結構最終在剪切的作用下失穩(wěn)破碎，混合層流場中的結構逐漸趨于管狀的渦結構形態(tài)。在充分發(fā)展階段，混合層中充斥著更加細碎的小尺度三維管狀渦結構。這些位于不同階段的各種渦結構形態(tài)，對混合層的接觸面和流質(zhì)混合產(chǎn)生了不同的作用，特別是初期演化的Λ渦結構、馬蹄渦、團花等大尺度結構從主流獲得的能量和質(zhì)量，使得混合層的界面在法向和展向產(chǎn)生很大波動，使混合層的動量厚度沿流向急劇增長，而在充分發(fā)展湍流階段，細碎渦結構作用下的均勻混合主要是在大尺度營造的可視混合層厚度量級的基礎上進行的。

［1］LELE S K.Compressibility effects on turbulence［J］.Annual Review of Fluid Mechanics，1994，26：211－254.

［2］BARRE S，QUINE C，DUSSAUGE J P.Compressibility effects on the structure of supersonic mixing layers：experimental results［J］.Journal of Fluid Mechanics，1994，259：47－78.

［3］CLEMENS N T，MUNGAL M G.Large－scale structure and entrainment in the supersonic mixing layer［J］.Journal of Fluid Mechanics，1995，284：171－216.

［4］XIONG H L，CUI E J.Studies on flow dynamics of compressible mixing layer and mixing enhancement［A］.Proceedings of frontier research symposium on aerodynamics［C］，2003：491－497.（in Chinese）熊紅亮，崔爾杰.可壓縮混合層流動特性與增混措施研究［A］.空氣動力學前沿研究學術研討會論文集［C］，2003：491－497.

［5］ZHAO Y X，YI S H，TIAN L F，et al.The fractal measurement of experimental images of supersonic turbulent mixing layer［J］.Science in China Series G：Physics Mechanics and Astronomy，2008，51（N8）：1134－1143.

［6］SANDHAM N D，REYNOLDS W C.Three dimensional simulations of large eddies in the compressible mixing layer［J］.Journal of Fluid Mechanics，1991，224：133－58.

［7］VREMAN B，KUERTEN H，GEURTS B.Shocks in direct numerical simulation of the confined three－dimensional mixing layer［J］.Physics of Fluids，1995，7（9）：2105－2107.

［8］FU D X，MA Y W.Direct numerical simulation of three－dimensional temporally evolving compressible mixing layer［J］.ACTA Mechanica Sinica，1998，30（2）：129－137.（in Chinese）傅德薰，馬延文.時間發(fā)展平面混合流的三維演化［J］.力學學報，1998，30（2）：129－137.

［9］FU D X，MA Y W，ZHANG L B.Direct numerical simulation of transition and turbulence in compressible mixing layer［J］.Science in China（Series A），2000，43（4）：421－429.

［10］KOURTA A，SAUVAGE R.Computation of supersonic mixing layers［J］.Physics of Fluids，2002，14（11）：3790－3797.

［11］ROGERS M M，MOSER R D.The three－dimensional evolution of a plane mixing layer Part 1.the Kelvin－Helmholtz rollup［R］.NASA TM－103856，N92－11303，1992.

［12］LI Q，ZHANG H X，GAO S C.Numerical simulations on supersonic shear layer flow［J］.ACTA Aerodynamica Sinica，2000，18（S1）：67－77.（in Chinese）李沁，張涵信，高樹椿.關于超聲速剪切流動的數(shù)值模擬［J］.空氣動力學學報，2000，18（S1）：67－77.

［13］LI Q B.Numerical study of compressible mixing layer with BGK scheme［D］.Beijing：Tsinghua University，2002.（in Chinese）李啟兵.應用BGK格式對可壓縮混合層的數(shù)值研究［D］.［博士學位論文］.北京：清華大學，2002.

［14］FU S，LI Q B.Numerical simulation of compressible mixing layers［J］.Int.J.Heat Fluid Flow，2006，27：895－901.

［15］PAN H L，MA H D，WANG Q.Transition coherent structures and shocklets in 3－D spatial developing mixing layers at high convective Mach numbers［J］.ACTA Aerodynamica Sinica，2008，26（3）：275－281.（in Chinese）潘宏祿，馬漢東，王強.高對流Mach數(shù)三維混合層轉(zhuǎn)捩特性分析及小激波結構模擬［J］.空氣動力學學報，2008，26（3）：275－281.

［16］LI Q，LUO J R，ZHANG H X.Studies on several issues of three－dimensional temporally and spatially developing mixing layer［J］.Physics of Gases－Theory and Applications，2009，4（1）：45－57.（in Chinese）李沁，羅俊榮，張涵信.三維時間和空間發(fā)展混合層若干問題研究［J］.氣體物理－理論與應用，2009，4（1）：45－57.

［17］ZHOU Q，HE F，SHEN M Y.The evolution of three－dimensional temporally evolving plane mixing layers under strong vortex disturbances［J］.Tsinghua Science and Technology，2009，14（S2）：17－21.

［18］STEGER J L，WARMING R F.Flux vector splitting of the inviscid gasdynamic equations with application to finite difference methods［J］.Journal of Computational Physics，1981，40：263－293.

［19］ZHOU Q，YAO Z H，HE F，et al.A new family of high－order compact upwind difference schemes with good spectral resolu－tion［J］.Journal of Computational Physics，2007，227（2）：1306－1339.

［20］SPITERI R J，RUUTH S J.Non－linear evolution using optimal fourth－order strong－stability－preserving Runge－Kutta methods［J］.Math.Comput.Simul.，2003，62：125－135.

［21］SANDHAM N D，REYNOLDS W C.Compressible mixing layer：linear theory and direct simulation［J］.AIAA J.，1990，28：618－624.

［22］ZHOU Q，HE F，SHEN M Y.Direct numerical simulation of a spatially developing compressible plane mixing layer：flow structures and mean flow properties［J］.Journal of Fluid Mechanics，2012，711：437－468.

［23］PANTANO C，SARKAR S.A study of compressibility effects in the high－speed turbulent shear layer using direct simulation［J］.Journal of Fluid Mechanics，2002，451：329－371.

［24］WU X H，MOIN P.Direct numerical simulation of turbulence in a nominally zero－pressure－gradient flat－plate boundary layer［J］.Journal of Fluid Mechanics，2009，630：5－41.

［25］HO C M，HUERRE P.Perturbed free shear layers［J］.Annual Review of Fluid Mechanics，1984，16：365－424.

［26］HEISENBERG W.Zur statistischen theorie der turbulenz［J］.Z.Phys.，1948，124：628－657.

［27］BROWAND F K，LATIGO B O.Growth of the two－dimensional mixing layer from a turbulent and nonturbulent boundary layer［J］.Physics of Fluids，1979，22（6）：1011－1019.