李子亮，于 洋，劉登豐

1 引言

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

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

超音速氣膜冷卻的流動(dòng)包含自由剪切流與附體流，而目前工程中廣泛應(yīng)用的標(biāo)準(zhǔn)k-ε模型[7]只適用于計(jì)算高雷諾數(shù)湍流。Menter[8]提出的SST k-ω模型結(jié)合了k-ε模型與k-ω模型[9]的特點(diǎn)，對(duì)于近壁面低雷諾數(shù)流動(dòng)采用標(biāo)準(zhǔn)k-ω模型計(jì)算，而在遠(yuǎn)壁面完全湍流區(qū)切換為k-ε模型，因此該模型可用于受限空間超聲速氣膜冷卻的數(shù)值模擬。然而，現(xiàn)有SST k-ω湍流模型繼承了k-ε模型的缺陷，在用于預(yù)測(cè)超聲速氣膜冷卻之前，需要進(jìn)行修正。

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

2 模型建立與修正

2.1 控制方程

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

對(duì)于超音速氣體，應(yīng)滿足如式(4)所示理想氣體狀態(tài)方程。

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

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

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

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

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

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

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

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

修正函數(shù)CT的表達(dá)式為式(15)。

式(15)~式(18)中，C1、C2和m為修正函數(shù)CT中的常數(shù)，Mt為湍流馬赫數(shù)，Mt≤0.1時(shí)，f(Mt)=0，即不進(jìn)行可壓縮修正，a為當(dāng)?shù)芈曀?。因此，本模型修正是在可壓縮SST k-ω模型基礎(chǔ)上進(jìn)行溫度修正。

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

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

3 模型驗(yàn)證

3.1 非等溫可壓縮湍流混合

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

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

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

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

圖2 SST k-ω湍流模型修正前后計(jì)算得到的湍動(dòng)能分布(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 為不同湍流模型計(jì)算得到的非等溫可壓縮湍流射流軸向速度分布，從圖中可知溫度與可壓縮修正后的SST k-ω湍流模型計(jì)算結(jié)果與相同條件下的實(shí)驗(yàn)結(jié)果吻合，而可壓縮修正的k-ω模型與SST k-ω模型的計(jì)算結(jié)果與實(shí)驗(yàn)值相差較大，這進(jìn)一步說明了溫度和可壓縮修正的SST k-ω湍流模型對(duì)于大溫度梯度自由剪切流預(yù)測(cè)的準(zhǔn)確性。在此基礎(chǔ)上，將修正模型用于預(yù)測(cè)同時(shí)含自由剪切流和附體流的超音速氣膜冷卻行為的準(zhǔn)確性進(jìn)行驗(yàn)證。

圖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]的絕熱平板超聲速氣膜冷卻實(shí)驗(yàn)數(shù)據(jù)對(duì)修正模型的計(jì)算結(jié)果進(jìn)行驗(yàn)證。并采用與實(shí)驗(yàn)相同幾何尺寸與邊界條件進(jìn)行建模，相關(guān)參數(shù)如表2所示。

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

圖4為平板超聲速氣膜冷卻計(jì)算域，將計(jì)算域劃分為256 089個(gè)網(wǎng)格，對(duì)壁面和冷卻氣體入口附近網(wǎng)格做加密處理，第一層網(wǎng)格節(jié)點(diǎn)距壁面0.005 mm。

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

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

表3 超聲速氣膜冷卻實(shí)驗(yàn)方案與入射參數(shù)Table 3 Experimental scheme and injection parameters of supersonic film cooling

圖5為空氣主流作用下平板壁面熱流軸向分布曲線，圖6為空氣主流作用下凸臺(tái)下游平板壁面熱流軸向分布曲線。從圖中可知，溫度和可壓縮修正的SST k-ω模型計(jì)算得到的無凸臺(tái)平板x＞0 m區(qū)域(圖5)和凸臺(tái)下游(圖6)的壁面熱流與實(shí)驗(yàn)值吻合最好，而采用可壓縮修正的k-ω和SST k-ω模型的計(jì)算值均較實(shí)驗(yàn)值偏高。這是因?yàn)闇囟刃拚腟ST k-ω模型是在可壓縮修正的基礎(chǔ)上考慮了溫度變化對(duì)剪切混合的影響，進(jìn)行了溫度修正，彌補(bǔ)了現(xiàn)有模型過高或過低預(yù)測(cè)剪切層增長(zhǎng)速率的缺陷，模型更加完善，從而使計(jì)算結(jié)果與實(shí)際更為接近。

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

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

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

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

4 結(jié)論

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

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

3)溫度修正的SST k-ω模型有待通過更多的實(shí)驗(yàn)數(shù)據(jù)進(jìn)一步驗(yàn)證和完善，進(jìn)而指導(dǎo)液體火箭發(fā)動(dòng)機(jī)噴管超聲速氣膜冷卻方案的設(shè)計(jì)。