基于SST k-ω的超聲速氣膜冷卻湍流模型溫度修正

李子亮，于 洋，劉登豐

1 引言

超燃沖壓發動機和大推力液體火箭發動機的高熱負荷組件的熱防護面臨嚴峻考驗。長期以來，氣膜冷卻的應用主要集中在以燃氣渦輪冷卻為代表的亞聲速范圍。近年來，隨著先進發動機的發展，超聲速氣膜冷卻以其更高的冷卻效率逐漸應用于超燃沖壓發動機燃燒室和大推力火箭發動機噴管的熱防護。與常規亞聲速氣膜冷卻相比，超聲速氣膜冷卻壓縮性作用顯著，氣體動量與能量緊密耦合，機理十分復雜。目前超聲速氣膜冷卻可壓縮效應尚未完全清楚，溫度梯度和速度梯度耦合作用下的剪切層混合傳熱規律有待進一步闡明。

隨著計算流體力學(Computational Fluid Dynamics，CFD)的發展，CFD模擬成為研究超聲速氣膜冷卻的重要手段。在超聲速氣膜冷卻數值模擬中，從不可壓縮流動發展而來的湍流模型是最具不確定性因素的，在預測高速流動時可壓縮性不可忽略，需要對模型進行可壓縮修正。多年以來，以Sarkar[1-2]和Henze[3]為代表的研究者以標準k-ε模型為基礎進行了大量修正工作。現有湍流模型已包含自由剪切層可壓縮耗散的修正項，但并未將溫度變化對剪切層增長速率的影響考慮在內，導致超聲速氣膜冷卻可壓縮流動傳熱預測結果不理想[4-6]。

超音速氣膜冷卻的流動包含自由剪切流與附體流，而目前工程中廣泛應用的標準k-ε模型[7]只適用于計算高雷諾數湍流。Menter[8]提出的SST k-ω模型結合了k-ε模型與k-ω模型[9]的特點，對于近壁面低雷諾數流動采用標準k-ω模型計算，而在遠壁面完全湍流區切換為k-ε模型，因此該模型可用于受限空間超聲速氣膜冷卻的數值模擬。然而，現有SST k-ω湍流模型繼承了k-ε模型的缺陷，在用于預測超聲速氣膜冷卻之前，需要進行修正。

本文擬在現有SST k-ω湍流模型基礎上，充分考慮剪切層溫度變化對冷卻氣體與主流混合傳熱的影響，對其進行溫度修正，并通過實驗數據對修正后的模型進行驗證。

2 模型建立與修正

2.1 控制方程

超音速氣膜冷卻流動過程滿足粘性可壓縮N-S方程，流體傳輸過程遵循質量守恒、動量守恒和能量守恒，其控制方程如式(1)~(3)所示。

對于超音速氣體，應滿足如式(4)所示理想氣體狀態方程。

式(1)~式(4)中，u為速度，q為熱流量，σij為粘性張量，cp為比熱容，ρ為密度，T為溫度，p為壓力，R為氣體常數。

2.2 SST k-ω湍流模型溫度修正

SST k-ω雙方程湍流模型中湍動能k與比耗散率ω的輸運方程如式(5)~式(10)所示。

SST k-ω模型湍流粘度表達式如式(11)~(13)所示。

式(5)~式(13)中，τij為雷諾應力張量，Ω為渦度張量，μ為層流粘度，μt為湍流粘度，y為到壁面的距離，模型常數 a1=0.31，β?=0.09，κ=0.41，σk1=0.85， σω1=0.5， σk2=1.0， σω2=0.856。

由于超聲速氣膜冷卻剪切層存在較大溫度梯度，本文模型修正的目的是將溫度變化對剪切層發展的影響加入到SST k-ω模型中，方法如下。

通過總溫梯度和湍流長度尺度定義修正變量Tg，其表達式為式(14)。

式中，Tt為氣體總溫，ΔTt為總溫梯度， k1/2/ω為湍流長度尺度。

修正函數CT的表達式為式(15)。

式(15)~式(18)中，C1、C2和m為修正函數CT中的常數，Mt為湍流馬赫數，Mt≤0.1時，f(Mt)=0，即不進行可壓縮修正，a為當地聲速。因此，本模型修正是在可壓縮SST k-ω模型基礎上進行溫度修正。

溫度修正后的SST k-ω湍流粘度為式(19)。

式中，Tj為冷卻氣體溫度，Tm為主流溫度，在超聲速氣膜冷卻中，當冷卻氣體入射溫度小于主流溫度時，取n=-1。該修正模型適用于計算剪切層總溫梯度較大的湍流混合耗散過程。

3 模型驗證

3.1 非等溫可壓縮湍流混合

采用Sumi[5]的超聲速低溫射流與高溫氣體自由剪切混合的實驗數據對模型的修正效果進行初步驗證。建立與實驗相同的二維計算域，并采用相同的邊界條件，如表1和圖1所示。將計算域劃分為173 000個網格，壁面絕熱無滑移，氣體為可壓縮氧氣，采用密度基求解器進行計算。

采用修正前后的SST k-ω湍流模型計算溫度為Tj=190 K超音速氣體射入1002 K高溫氣體環境中時得到的湍動能分布如圖2所示。由圖可知，修正后的模型與修正前相比，剪切層垂直于射流方向，湍動能增長受到抑制，沿射流方向湍動能衰減變慢，這表明溫度修正后的SST k-ω湍流模型有效抑制了大溫度梯度下剪切層的增長速率，修正效果顯著。

表1 非等溫可壓縮自由剪切流邊界條件Table 1 The boundary condition of the non-thermal free shear flow

圖1 非等溫可壓縮湍流混合計算域網格劃分Fig.1 Mesh generation in computational domain for non-isothermal compressible turbulent mixing

圖2 SST k-ω湍流模型修正前后計算得到的湍動能分布(T j=190 K,T a=1002 K)Fig.2 The calculated turbulent kinetic energy distribution before and after the SST k-ωmodel modification(T j=190 K,T a=1002 K)

圖3 為不同湍流模型計算得到的非等溫可壓縮湍流射流軸向速度分布，從圖中可知溫度與可壓縮修正后的SST k-ω湍流模型計算結果與相同條件下的實驗結果吻合，而可壓縮修正的k-ω模型與SST k-ω模型的計算結果與實驗值相差較大，這進一步說明了溫度和可壓縮修正的SST k-ω湍流模型對于大溫度梯度自由剪切流預測的準確性。在此基礎上，將修正模型用于預測同時含自由剪切流和附體流的超音速氣膜冷卻行為的準確性進行驗證。

圖3 非等溫可壓縮湍流射流軸向速度分布(T j=190 K,T a=1002 K)Fig.3 Axial velocity distribution of the non-isothermal compressible turbulent(T j=190 K,T a=1002 K)

3.2 超聲速氣膜冷卻傳熱

采用Holden[10]的絕熱平板超聲速氣膜冷卻實驗數據對修正模型的計算結果進行驗證。并采用與實驗相同幾何尺寸與邊界條件進行建模，相關參數如表2所示。

表2 平板氣膜冷卻計算域尺寸參數Table 2 The calculating domain size of flat plate film cooling /m

圖4為平板超聲速氣膜冷卻計算域，將計算域劃分為256 089個網格，對壁面和冷卻氣體入口附近網格做加密處理，第一層網格節點距壁面0.005 mm。

圖4 平板超聲速氣膜冷卻計算域網格劃分Fig.4 Mesh generation in computational domain for supersonic film cooling of flat pate

超聲速氣膜冷卻主流為空氣，冷流為氦氣，氦氣通過小噴管加速后在凸臺下游與空氣進行摻混，將小噴管出口參數作為冷卻氣體在流場中的入射邊界參數。具體的實驗方案如表3所示。方案1為無氣膜冷卻下的超聲速空氣主流在無凸臺平板上流動，方案2為無氣膜冷卻的超聲速空氣主流在下游壁面流動，方案3為氣膜冷卻下超聲速氦氣與空氣主流在下游混合流動。壁面條件設置為光滑絕熱無滑移，壁面溫度為294.26 K，選擇密度基求解器進行計算求解。

表3 超聲速氣膜冷卻實驗方案與入射參數Table 3 Experimental scheme and injection parameters of supersonic film cooling

圖5為空氣主流作用下平板壁面熱流軸向分布曲線，圖6為空氣主流作用下凸臺下游平板壁面熱流軸向分布曲線。從圖中可知，溫度和可壓縮修正的SST k-ω模型計算得到的無凸臺平板x＞0 m區域(圖5)和凸臺下游(圖6)的壁面熱流與實驗值吻合最好，而采用可壓縮修正的k-ω和SST k-ω模型的計算值均較實驗值偏高。這是因為溫度修正的SST k-ω模型是在可壓縮修正的基礎上考慮了溫度變化對剪切混合的影響，進行了溫度修正，彌補了現有模型過高或過低預測剪切層增長速率的缺陷，模型更加完善，從而使計算結果與實際更為接近。

圖5 空氣主流作用下平板壁面熱流分布Fig.5 Heat flow distribution on flat plate wall under the action of air mainstream

圖6 空氣主流作用下凸臺下游平板壁面熱流分布Fig.6 Heat flow distribution on upstream flat plate wall of the boss under the action of air mainstream

圖7 為平板凸臺下游超聲速冷卻氦氣與空氣相互作用下壁面熱流的軸向分布曲線。從圖中可知，在初始混合區(x＜0.1 m)計算結果與實驗結果存在加大差距，這是由于本文計算中將冷卻氣體入射邊界y方向設置為均一參數，使初始區壁面邊界層氣體溫度過低，熱流偏高。比較圖中3個湍流修正模型在充分混合區的結果可知，溫度和可壓縮修正的SST k-ω模型與實驗結果最為接近。這是因為經過溫度修正后的模型可以更加準確預測冷卻氣體與主流剪切層的混合發展速率；此外，比較與溫度和可壓縮修正的SST k-ω模型的計算結果可知，在前半段與溫度修正模型計算結果接近，而在充分發展區二者計算結果差距變大，可壓縮修正的k-ω模型對于湍流充分發展區的預測表現出一定的局限性。由此可知，在可壓縮修正的SST k-ω模型基礎上構建的溫度修正模型適用于預測超聲速氣膜冷卻湍流剪切層的發展變化，確保了冷卻氣體與高溫主流剪切混合傳熱計算結果的準確性。

圖7 氣膜冷卻作用下凸臺下游平板壁面熱流分布Fig.7 Heat flow distribution on upstream flat plate wall of the boss under air film cooling

4 結論

1)溫度變化對超聲速氣膜冷卻剪切層湍流混合耗散影響顯著，基于總溫梯度變化的SST k-ω修正模型彌補了現有湍流模型的理論缺陷，大幅度提高了非等溫超聲速自由剪切流的計算精度。

2)溫度修正的SST k-ω模型對絕熱平板超聲速氣膜冷卻的計算結果與實驗結果吻合，對于剪切層總溫大梯度變化的超聲速氣膜冷卻數值模擬，該模型與現有可壓縮修正模型相比具有顯著優越性。

3)溫度修正的SST k-ω模型有待通過更多的實驗數據進一步驗證和完善，進而指導液體火箭發動機噴管超聲速氣膜冷卻方案的設計。