馬赫數6柱-裙激波/邊界層干擾直接模擬

孫東, 劉朋欣，沈鵬飛, 童福林, 郭啟龍

空氣動力學國家重點實驗室，綿陽 621000

激波/邊界層干擾是超聲速以及高超聲速中常見的流動現象，是高速飛行器設計中必須考慮的重要問題之一。20年前，Dolling[1]在對當時激波/邊界層干擾研究進行總結后指出，激波/邊界層干擾低頻非定常效應和三維效應將是未來激波/邊界層干擾研究的重點方向。

與激波的非定常運動相比，激波/邊界層干擾的三維效應受到的重視存在明顯不足。實際上，即使是二維構型上(如壓縮拐角、斜激波反射)的激波/邊界層干擾也存在顯著的三維效應。Bookey等[10]在激波反射的實驗中觀察到分離區附近的拓撲結構存在類似“Owl eye”的流動結構，這預示著分離區在展向存在大尺度的流動結構。這種展向三維結構在直接數值模擬結果[11-12]中也曾被發現。Loginov等[13]通過對二維壓縮拐角的大渦模擬，發現在干擾區后存在流線的匯聚與分散。這與Zheltovodov[14]的實驗結果一致。他們認為這種流向結構可能是G?rtler渦。Grilli等[15]在對壓縮膨脹角的激波干擾研究中對這種G?rtler渦的展向差異進行了研究。Tong等[16]對拐角倒圓的壓縮拐角的激波干擾進行了研究，并得到了G?rtler渦的精細結構。Helm和Martin[17]對激波干擾中的G?rtler渦結構開展了系統的研究，指出G?rtler渦與分離區和上游流動結構的相互干擾密切相關，有可能是分離激波運動的原因之一。

總體而言，目前對激波/邊界層干擾的研究仍集中在超聲速范圍，對于高超聲速情況下的激波干擾還鮮有報道。與超聲速情況相比，更高的馬赫數將帶來顯著的可壓縮效應，并且激波與邊界層發生強烈相互干擾，由此引起的壓力和熱流載荷的增加要遠高于較低馬赫數下的流動。為了研究高超聲速情況下的激波/邊界層干擾問題，本文采用高超聲速柱-裙構型開展直接數值模擬，對激波/邊界層干擾后的G?rtler結構及其對分離泡和表面壓力、熱流的影響進行詳細研究。

內容安排如下：第1節給出計算采用的數值方法和算例設置；第2節給出結果分析，包含基本流場信息、G?rtler渦對分離泡和表面壓力、熱流的影響以及對不同展向截面進行正交特征分解；第3節給出結論。

1 數值方法與計算設置

本文采用的控制方程為曲線坐標系下的Navier-Stokes方程，其具體表達式可參考文獻[18]。無量綱參考長度為毫米，參考速度和溫度等采用無窮遠處的流場變量值。若無特別說明，本文中長度均采用毫米進行無量綱化。采用高精度有限差分計算軟件OpenCFD-SC進行計算，該軟件已經成功應用于超聲速/高超聲速激波/邊界層干擾的數值模擬中，計算的可靠性和準確性得到了充分的驗證[19-21]。對流項計算格式為優化的六階單調保持格式OMP6[22]，并且在強激波附近，添加自適應濾波來進一步穩定計算。黏性項計算格式為六階中心格式，時間推進為三階TVD性質的龍格-庫塔方法。來流條件設置如下：來流馬赫數Ma=6，毫米雷諾數Remm=12 000，來流溫度T∞= 65 K，壁面溫度Tw=305.5 K。

計算構型為柱-裙構型，如圖1所示。圖中r為徑向，前體圓柱半徑R0為137，計算域法向高度為40；z為流向，總長約為450，其中圓柱面長度L為375。s表示展向坐標(在圖中未顯示)，弧度為0.102，圓柱表面的展向弧長為14。裙體與柱面夾角θ為33°。在圖1中給出了柱-裙構型二維示意圖以及網格分布。進口為層流剖面，上邊界為無反射邊界[23]，出口為超聲速出口條件，并對出口附近網格進行拉伸，以進一步抑制反射波的前傳。展向為周期邊界條件。壁面條件為等溫無滑移壁面，并在進口下游設置吹吸擾動(圖中紅色標注)[24]以促發轉捩。吹吸擾動vbs通過在壁面處引入法向速度實現，具體表達式為

圖1 柱-裙外形以及網格分布

(1)

式中：A=0.2，為吹吸擾動幅值；ωs=128，用于控制吹吸擾動周向周期數；θ為周向角度；φx為控制吹吸的流向強度；xb和xe分別為壁面擾動開始和結束的流向位置。

為考察網格分布對計算結果的影響，共采用3套網格進行網格收斂性分析。網格點數以及在湍流充分發展區域內(z=-240)網格間距等信息在表1中給出。Δx+、Δr+和Δz+為采用黏性長度無量綱化的流向、徑向和展向的網格間距。從表中可以看出，從Grid-C到Grid-F，網格逐步加密。

表1 網格總數和網格間距信息

圖2給出了不同網格的表面壓力Pw/P∞、摩阻系數Cf和熱流系數Ch的計算結果，橫坐標為采用湍流邊界層約化后的流向距離。邊界層厚度δ≈6.5，為z=-240處的邊界層厚度，z0為拐角的流向位置。從圖中可以看出，在不同網格密度下，上游邊界層上的壓力、摩阻和熱流的分布差別很小，而在拐角附近，Grid-M和Grid-F的結果相差較小，Grid-C預測的分離區偏小。在干擾之后，不同網格密度的結果存在一定的差別。但總體而言，Grid-F和Grid-M預測的壓力、摩阻和熱流的分布趨勢都較為一致。以上結果體現出了較好的網格收斂性。考慮到計算精度和分辨率，下面采用Grid-F的結果進行分析。

圖2 不同網格間距下表面壓力、摩阻和熱流系數沿流向分布

2 結果分析

2.1 上游湍流邊界層

圖3 瞬時密度云圖

圖4 干擾區前湍流邊界層剖面和雷諾應力分布

由于計算構型半徑遠大于弧長，可以將局部的圓柱表面近似為平板，并與其他超聲速以及高超聲速的平板DNS結果[25-27]比較。從湍流邊界層速度剖面可以看出，平均速度剖面具有典型的線性區、對數區以及尾跡區特征。在對數區截距為5.5，較一般超聲速平板的截距偏高，這一現象在其他高超聲速平板的文獻中也有所體現[25-27]。從雷諾應力分布上看，目前雷諾應力分布在r/δ>0.2的邊界層內，與其他超聲速湍流平板數據的雷諾應力分布較為一致。雷諾應力分布的主要差異體現在r/δ<0.2的范圍內，考慮到圓柱表面與平板流動的差異，這一差異在可以接受的范圍內。這表明，本文算例中干擾區前的湍流邊界層在平均速度剖面和雷諾應力分布上都具有典型合理的湍流特征。

2.2 基本流場結構

圖5中給出了Q準則顯示的瞬時旋渦結構，并用法向距離進行著色，從圖中可以看出來流湍流結構在經過干擾之后流動結構的尺度得到顯著增強。為顯示干擾區附近的流動結構，進一步考察溫度和密度的分布情況。

圖5 壁面法向距離著色的瞬時旋渦結構

從圖6的溫度分布上看，溫度的極高值集中出現在分離之后的干擾區內，之后迅速降低，這種分布規律與二維壓縮拐角有顯著不同，主要原因可能與柱-裙結構裙體的膨脹過程有關。在平均溫度云圖中，黑色實線給出了流向速度為0.99倍來流值以顯示邊界層厚度，黑色虛線顯示聲速線，可以看出高溫度區域集中在聲速線以下的分離泡內部。

圖6 干擾區附近瞬時(上)和平均(下)溫度云圖

圖7給出了瞬時數值紋影和平均密度分布以顯示柱-裙激波干擾中的波系結構。數值紋影的計算公式[28]為

NS=c1exp[-c2(x-xmin)/(xmax-xmin)]

(2)

式中：x為密度梯度幅值；系數c1=0.8，c2=20。

從數值紋影中可以看到分離激波與湍流邊界層之間存在強烈的相互作用，再附位置邊界層被嚴重壓縮在激波之下，形成較為明顯的反射激波，之后邊界層結構出現快速增長，出現了少量大尺度的流動結構。在平均密度圖上可以看出，經過再附之后出現明顯的膨脹過程，并且在壁面附近的密度值出現了升高，這與圖6中溫度降低密切相關。

圖8給出了平均馬赫數和湍流馬赫數云圖。湍流馬赫數為速度脈動的均方根與平均聲速的比值用以顯示流場中的可壓縮性。從馬赫數云圖上看，在經過分離激波之后馬赫數迅速降低至2左右，這與斜激波關系式是符合的。從湍流馬赫數云圖上看，在干擾區附近湍流馬赫數急劇升高，體現了極強的可壓縮效應，這可能對流動結構造成顯著的影響。

圖8 干擾區附近平均馬赫數(上)和湍流馬赫數(下)云圖

圖9給出了流向速度等于零的瞬時等值面以顯示分離泡的瞬時結構，并采用壁面距離著色。從圖中可以看出，來流邊界層在分離泡前形成了多個展向的條帶結構，分離泡內部在展向呈現出強烈的三維效應，并且在分離泡下游出現了3條長條形的流向結構，推測為G?rtler渦結構。

圖9 壁面法向距離著色的流向速度為零等值面

在壓縮拐角和柱-裙的干擾區下游，由于分離泡引起邊界層抬升，引起當地等效物面曲率變化，在離心不穩定性的作用下形成G?rtler渦。為進一步證實此推測，圖10給出了表面極限流線。可以看出，分離附近流動拓撲呈現出結點和鞍點交替分布的特征，在分離區內流動在展向存在大尺度結構，在再附位置出現兩個結點，結點之間為鞍點。在結點之后出現了流線的匯聚和分散，可以說明在展向存在兩個大尺度的流向渦結構。進一步對再附線的間距和分離線的間距進行了估計，發現間距均為6.7左右，略大于來流邊界層厚度。在離心不穩定性的分析中，當地邊界層厚度和流線曲率的比值δ/R常用來作為G?rtler不穩定性的判別準則。圖11中給出了δ/R的比值。沿著法向從上之下選取3條流線，可以看到δ/R的峰值位置逐步后移。Smits和Dussauge[29]在Ma=3的湍流給出δ/R的閾值為0.03，而本算例的3條流線在拐角附近大于此閾值，因此認為此處發生的為G?rtler不穩定性。

圖10 表面極限流線(N表示結點，S表示鞍點)

圖11 當地邊界層厚度與流線曲率的比值

在G?rtler渦的影響下，干擾區內的展向結構呈現除了顯著的三維效應，并進而影響分離泡結構以及表面壓力、熱流的展向分布，接下來對這些方面進行詳細分析。

2.3 分離泡結構

首先對不同展向位置上的分離泡結構進行分析。根據上文的分析，選取如下展向位置：展向位置1位于展向中線上，同時也是一條G?rtler渦的分離線；展向位置2位于展向1/4位置，同時也是一條G?rtler渦的再附線。兩個展向位置分別對應了G?rtler渦結構在展向的分離和再附流動過程。為方便敘述，定義展向分離位置為位置1，展向再附位置為位置2。圖12給出了瞬時流向速度w/w∞云圖，并且以粉色實線表示平均流向速度為零的等值線。從此等值線可以看出，位置2的回流區域比位置1明顯大很多。但從流向速度云圖上看，在不同展向位置上，出現回流(流向速度小于零)的區域并不規則，不是比較完整的整塊區域。比如，位置2處的回流區域在z=-180位置分成了兩個前后相對獨立的區域。為了對兩個展向位置處的回流區域進行定量的比較，對回流出現的概率密度函數(PDF)進行了計算，結果如圖13所示，粉色實線仍為平均流向速度為零的等值線，背景云圖為平均流向速度云圖，以顯示分離泡結構和邊界層信息，概率密度以等值線的形式體現。根據Simpson[30]的定義，概率等于0.01為初始分離(Incipient Detachment)、概率等于0.2為間歇性瞬變分離(Intermittent Transitory Detachement)、概率等于0.5為瞬變分離(Transitory Detachment)，最后一種分離狀態為平均意義上的分離。從圖13中可以看出位置2的瞬變分離區域要遠大于位置1的瞬變分離區域。

圖12 不同展向位置瞬時流向速度云圖

對回流區面積進行統計，具體做法為遍歷這兩個展向位置的干擾區域，統計流向速度不大于零的網格單元面積之和，作為瞬時回流面積。統計結果在圖14中給出。采用藍色點劃線和紅色虛線分別表示展向位置1和2的平均回流面積。位置1的回流面積約為16，而位置2的面積為約為28。這與圖12和圖13的分析結果是一致的。進一步對圖14中的數據進行功率譜密度分析，具體結果在圖15中給出。可以看出，兩個展向位置上回流面積上的能量主要以中低頻為主，高頻能量占比很小。中頻的能量集中在fδ/U∞=0.07附近，低頻的能量集中在fδ/U∞=0.02附近，這個能量與激波振蕩的低頻能量十分接近。在圖15中還可以看出，位置1的中頻部分附近能量要高于位置2，相反在低頻部分附近能量低于位置2。這種變化可能與它們各自的回流面積相關，位置1處的回流面積要明顯小于位置2。

圖13 不同展向位置處干擾區內回流的概率密度函數

圖14 不同展向位置處干擾區內逆流面積隨時間變化歷程

圖15 不同展向位置處干擾區逆流的功率譜密度

進一步分析兩個展向位置處回流面積運動過程的相互關系，對其變化歷程進行了相關性分析。互相關系數的定義為

(3)

式中：τ為延遲時間。

圖16 兩個展向站位回流面積的時間相關系數

2.4 物面壓力和熱流

對G?rtler渦引起的物面壓力和熱流的展向差異展開分析。圖17給出了時間平均壓力P/P∞和熱流系數Ch云圖，并在圖中標出了平均分離S、拐角O和再附位置R。從圖中可以看出在拐角之前，壓力和熱流在展向上分布較為均勻。在再附之后，壓力出現壓力極高值點，然后緩慢降低；熱流的極高值點出現在平均再附線之前，并在經過短暫升高之后迅速下降。

圖17 物面平均壓力(上)和熱流系數(下)云圖

圖18給出了兩個展向位置上壓力和熱流沿流線的分布。展向位置的選取與2.3節相同。從壓力分布上看，位置2的壓力峰值比位置1更為靠前，并且壓力峰值比位置1的壓力峰值高13%，之后開始緩慢降低。從熱流分布上看，兩個展向位置的峰值位置十分接近，但熱流峰值上，位置2比位置1要高出16.2%。一般而言，熱流和壓力的升高主要來自于激波干擾中激波的壓縮效應，而G?rtler渦進一步對其展向分布進行了調節。位置2為再附線附近，此時由于流向渦旋轉的作用，流體直接撞擊壁面，因此壓力和熱流會明顯高于位置1。可見，G?rtler渦對壓力和熱流在展向分布造成了顯著的影響。

圖18 物面平均壓力(上)和熱流系數(下)沿流向分布

圖19中壓力脈動和熱流脈動的云圖呈現出了相似的特征。在再附線附近出現兩個壓力和熱流的極高值區域，但很快就出現了下降。壓力脈動在拐角之后出現了明顯增長，而熱流的增長則始于分離線之后。壓力脈動的極高值約為平均壓力極高值的1/3以上，而熱流脈動的極高值與熱流的平均值相當。因此在實際設計中，應該充分考慮熱流脈動的影響。在考慮脈動熱流后，實際熱流值會顯著高于平均熱流。為定量考慮壓力和熱流脈動的展向差異，圖20給出了不同展向位置脈動沿流向的分布。結果顯示，位置1的脈動值要明顯低于位置2的脈動值。位置2的壓力脈動值比位置1高28%，熱流脈動值比位置1高20%。

圖19 壓力(上)和熱流系數(下)脈動的均方根云圖

圖20 不同展向位置平均壓力(上)和熱流系數(下)均方根的流向分布曲線

對位置1和位置2上的壓力和熱流進行了功率譜密度(PSD)分析，以研究不同位置上頻譜特性的差異，結果見圖21和圖22。從圖21的壓力頻譜中可以看出，不同展向位置上在平均分離位置均捕捉到了低頻信號，并且信號的強度和位置相差不大。這是由于不同展向位置的平均分離位置相差不大(如圖12和圖13所示)，而低頻信號主要出現在平均分離位置。從圖22的熱流功率譜密度中可以看出，兩個位置的熱流信號都沒有捕捉到低頻信號，熱流的峰值頻率干擾之前約為fδ/U∞=1以上，到干擾區下游移動到fδ/U∞=0.5附近。圖21和圖22顯示，在不同的展向位置壓力和熱流信號的頻譜基本相同。

圖21 不同展向位置物面壓力功率譜密度分析

圖22 不同展向位置物面熱流功率譜密度分析

2.5 特征正交分解

為了研究不同展向截面上流動結構的差異，采用特征正交分解(POD)對截面上的流向速度進行模態分解。特征正交分解的具體表達式可見文獻[31]。本文共選擇了1 154個樣本，采樣間隔為0.256δ/U∞，時間跨度大約為兩個低頻周期。

圖23給出了兩個展向位置上POD分解的特征值分布情況。可以看到位置1上低階模態包含更多的能量，而位置2則在第10個模態后的高階模態能量高于展向位置1。

圖23 不同展向位置POD分解特征值能量百分比分布

圖24和圖25分別給出了兩個展向位置上POD模態的空間分布情況，并采用黑色虛線標出平均聲速線的位置。從圖24和圖25中可以看出，展向截面上剪切層結構在流動中具有主導作用。從Mode-1～Mode-6的空間分布上看，POD的能量主要集中在聲速線附近，并且模態結構從Mode-1中在聲速線附近集中的單一結構，逐步演化為沿著聲速線出現正負交替的流動結構。在位置1(圖24)中Mode-3的主要結構分別位于聲速線上下兩側，這與其他低階模態有著明顯不同，推測該模態可能與分離泡結構相關。該模態能量占比為5%。在位置2(圖25)上同樣存在類似結構的模態，為Mode-4，能量占比為3%，低于位置1的相應模態的能量占比。

圖24 展向位置1上的POD模態空間分布

圖25 展向位置2上的POD模態空間分布

在得了流向速度的模態分布之后，可以采用不同數量的POD模態重構流場，以定量分析不同展向位置上POD模態的能量收斂速度。在得到重構的流向速度后，按照2.3節中的方式得到干擾區中回流面積隨時間的變化過程。圖26(a)給出了位置1上采用不同數量POD模態(模態數量分別為10，20，40和80)重構流場后，回流面積隨時間的變化過程，并且給出了tU∞/δ=100～200范圍內的結果，黑色虛線為平均分離泡大小。從圖中可以看出采用10個模態已經能夠較為準確地捕捉到分離泡的變化過程；隨著模態數的增加，一些局部的變化細節也可以被捕捉到。

圖26(b)給出了位置2上重構的回流面積隨時間的變化過程。可以看出，在采用10個模態的時候，回流面積的時間變化歷程已經能夠得到較為準確的刻畫。

圖26 不同展向位置POD模態重構的分離泡面積

為了定量顯示重構流場與原流場的差異，圖27給出了重構誤差的均方根分布。可以看出，位置2的重構誤差大于位置1的重構誤差。同時可以得到誤差與模態數量的對數函數關系，Error=klg(Δm)，其中Error為重構誤差，Δm=m-m0為當前模態數m減去起始模態數m0(本文中m0為10)，k為斜率。從圖27中可以看出位置1上斜率的絕對值為2.8，而位置2上斜率的絕對值為4.5。位置2上的重構流場收斂速度要高于位置1。

圖27 不同展向位置處POD模態的重構誤差

3 結 論

采用直接數值模擬對Ma=6條件下的柱-裙激波/湍流邊界層干擾進行了詳細研究。計算結果表明，流動分離會顯著影響分離泡附近的流線曲率，并繼而導致流動在離心不穩定性驅動下形成G?rtler渦結構。本文通過選取G?rtler渦分離和再附的兩個特征展向位置研究G?rtler渦對分離泡結構、物面壓力和熱流的影響，并采用特征正交分解研究了這兩個展向位置的模態結構特征，得到以下結論：

1)展向分離位置上分離泡的平均大小要明顯小于展向再附位置，并且兩個展向位置上分離泡的運動基本是同步的，分離泡面積的功率譜密度分析顯示在低頻和中頻范圍內都存在峰值。

2)壓力和熱流在展向再附線附近出現顯著的不均勻性，展向再附位置的平均壓力和熱流要比展向分離位置分別高13%和16.2%；其脈動壓力和熱流比展向分離位置高28%和20%，并且熱流的脈動峰值與平均值相當。因此在熱流計算中，應充分考慮熱流的非定常特征。

3)特征正交分解結果顯示兩個展向位置上的POD能量主要集中在剪切層附近，并且展向再附位置上低頻能量占比更高。采用POD模態重構流場并分析分離區面積發現，展向分離位置的重構誤差更小，而展向再附位置上的重構誤差收斂更快。