龔紅蘭，李 凌

(1.上海理工大學(xué)能源與動力工程學(xué)院，上海 200093；2.上海市動力工程多相流動與傳熱重點實驗室，上海 200093)

1 引言

拉瓦爾噴管廣泛應(yīng)用于工程中，比如超音速噴氣式飛機、火箭等飛行器中的噴管都屬于拉瓦爾噴管[1]。在發(fā)動機工作時常伴隨著高溫高壓燃氣的流動，內(nèi)部熱環(huán)境惡劣，對熱防護層材料性能要求較高，為此，許多文獻對拉瓦爾噴管內(nèi)的流動換熱以及固體壁面導(dǎo)熱進行了研究。噴管內(nèi)的傳熱過程非常復(fù)雜，涉及熱傳導(dǎo)、熱對流等多種傳熱方式，而且噴管的熱狀態(tài)不僅與其工況有關(guān)，還與噴管的結(jié)構(gòu)、材料物性等有關(guān)，所以應(yīng)當將流場與結(jié)構(gòu)溫度場耦合起來進行計算。張曉光等將軟件Fluent與ANSYS結(jié)合起來，完成了噴管喉襯流場及熱結(jié)構(gòu)的模擬分析[2]。孫林、郝雯等建立了二維軸對稱噴管模型，采用流固耦合方法，以噴管內(nèi)壁面溫度分布為邊界條件，研究了噴管在旋轉(zhuǎn)狀態(tài)下的流動及換熱情況[3-4]。白俊華等采用經(jīng)過實驗驗證的流固耦合換熱模型，對影響噴管喉襯熱交換的某些因素進行了研究[5]。吳川、王啟凡等同樣采用流固耦合的方法，借助流體力學(xué)軟件Fluent對噴管的結(jié)構(gòu)溫度場開展了瞬態(tài)仿真[6-8]。

目前，對噴管熱防護層進行數(shù)值計算時都是將噴管內(nèi)的純氣相流場結(jié)果作為熱邊界條件，但是對于某些特殊場合，比如應(yīng)用于固體火箭發(fā)動機這一領(lǐng)域，由于高能推進劑的使用，使得燃燒產(chǎn)物的溫度逐漸增大，并且內(nèi)部伴隨著顆粒的產(chǎn)生，而顆粒對氣相的影響不可忽略，必會導(dǎo)致兩相流場結(jié)果與純氣相不同，此外顆粒還會對噴管壁面造成嚴重的沖刷及燒蝕，給噴管的熱防護帶來了更大的挑戰(zhàn)，所以有必要對兩相流動產(chǎn)生的影響進行更加深入的分析與研究，以便進一步開展噴管的熱分析與結(jié)構(gòu)設(shè)計優(yōu)化工作。

故本文為了提高數(shù)值計算的精度，建立了含熱防護層的噴管兩相流動流固耦合[9]換熱模型，對拉瓦爾噴管內(nèi)部流動換熱進行了研究，并將噴管內(nèi)兩相流場計算結(jié)果作為結(jié)構(gòu)溫度場模擬的熱邊界條件，以獲得噴管熱防護層以及殼體各個時刻的溫度分布情況，為噴管熱分析提供依據(jù)。

2 數(shù)值模型

2.1 物理模型

鑒于本文中噴管的對稱性，建立二維軸對稱數(shù)值模型，如圖1所示。噴管內(nèi)流場部分設(shè)為流體域，其它各部分都為固體域。固體區(qū)域材料性能參數(shù)如表1所示，氣相密度按理想氣體處理，R=320J/(kg·K)，γ=1.2216，導(dǎo)熱率按kinetic-theory給定[10]，粘性按Sutherland定律處理(三系數(shù)形式)[11]；計算區(qū)域都采用四邊形結(jié)構(gòu)網(wǎng)格進行劃分，為了保證結(jié)果的精確性，在流固耦合交界面和噴管喉部進行局部加密，網(wǎng)格劃分情況如圖2所示。

圖1 噴管計算模型示意圖

表1 材料物性參數(shù)

圖2 計算模型網(wǎng)格劃分

為了便于分析，可對模型進行一些簡化，做出如下假設(shè)：

1)忽略內(nèi)壁面的燒蝕與炭化；

2)忽略各層材料之間的接觸熱阻；

3)忽略輻射傳熱和顆粒接觸傳熱；

4)忽略燃氣的化學(xué)反應(yīng)；

5)假設(shè)顆粒內(nèi)部和顆粒之間沒有物理和化學(xué)作用；

2.2 控制方程

本文模擬包含兩個過程，由于流場的發(fā)展速率遠大于噴管固體結(jié)構(gòu)的溫度傳遞速率，因而，在進行噴管結(jié)構(gòu)的熱傳導(dǎo)計算時，噴管流場都處于穩(wěn)定狀態(tài)。所以在計算時先進行噴管兩相流的穩(wěn)態(tài)數(shù)值模擬，然后將其結(jié)果當作熱邊界條件與噴管固體壁面進行瞬態(tài)耦合計算。噴管兩相流計算采用Eulerian方法，對于湍流多相流來說，其控制方程包括流體相與顆粒相的連續(xù)方程、動量方程及能量方程，由于方程較多，動量方程與能量方程見參考文獻[12]，連續(xù)性方程如下所示

1)流體相連續(xù)方程

(1)

顆粒相連續(xù)方程

(2)

式(1)(2)中符號和含義見參考文獻[12]。

為了使方程封閉，還需加上湍流輸運方程，如下所示：

3)流體相湍動能方程k

(3)

顆粒湍能輸運方程kk

(4)

流體湍流耗散方程ε

(5)

式(3)(4)中符號和含義見參見文獻[12]，式(5)中符號和含義見參考文獻[13]。

固體區(qū)域僅考慮導(dǎo)熱，其控制方程如下

(6)

式(6)中符號和含義見參考文獻[5]。

2.3 計算方法與邊界條件

為了更好的體現(xiàn)以兩相流計算結(jié)果和純氣相計算結(jié)果作為熱邊界條件的不同，本文采用擬流體模型(PFM)來模擬噴管內(nèi)的兩相流動，對顆粒相與氣相均采用Eulerian處理，該模型考慮了顆粒相的粘性、導(dǎo)熱及擴散，相比于離散相(DPM)模型，可以全面考慮顆粒相的湍流輸運。湍流模型選用RNG k-ε兩方程模型，采用標準壁面函數(shù)(standard wall function)求解近壁區(qū)域物理量。噴管進口為壓力進口邊界，總壓為7MPa，總溫為3000K，顆粒體積分數(shù)為0.001，氣相按理想可壓流處理；出口取壓力出口條件；噴管內(nèi)流場與固體域接觸表面設(shè)為耦合邊界條件；外壁為對流換熱邊界，h=10 W·m-2·K-1。計算時間5s，時間步長取0.001s。

3 計算結(jié)果分析

3.1 模型驗證

為了驗證本文模型的正確性，對文獻[14]中的JPL(Jet Propulsion Laboratory)噴管進行了純氣相模擬計算，噴管參數(shù)及邊界條件見文獻[15]，將本文結(jié)果與文獻[14]中實驗值進行了對照，如圖3所示，能看出本文結(jié)果和實驗值較吻合。

圖3 模擬結(jié)果與實驗結(jié)果對比

3.2 兩相流場計算結(jié)果分析

本文先對噴管內(nèi)的兩相流場進行了模擬，圖4是噴管內(nèi)兩相流場與純氣相流場計算結(jié)果的對照情況，其中圖4(a)是兩種工況下噴管內(nèi)氣相壓強沿軸向變化曲線，從圖中可以看出，顆粒相的存在對于噴管內(nèi)流場壓強分布的影響并不明顯，這是由于顆粒不會像氣體一樣膨脹做功導(dǎo)致的。圖4(b)為兩種工況下噴管內(nèi)氣相溫度沿軸向變化曲線，從圖中可以看出，加入顆粒相之后，同一位置處燃氣溫度要大于純氣相情況下的溫度，而且在噴管喉部和擴張段位置處顆粒相對氣相溫度的影響比較明顯，這是因為顆粒的比熱容較大，容易保持原來的溫度，而且在噴管喉部之后燃氣由于膨脹其溫度減小，與此同時，顆粒相與氣相發(fā)生熱交換而使得氣相溫度增加。但在噴管喉部之前顆粒相對氣相溫度的影響并不明顯，是因為這些位置燃氣溫度較高，顆粒相與氣相的熱交換較小。圖4(c)(d)分別是噴管內(nèi)氣相速度與馬赫數(shù)沿軸向的變化曲線，由這兩個圖可以看出，兩相流情況下相同位置處燃氣速度與馬赫數(shù)都比純氣相時要小，這是因為顆粒具有較大的慣性，存在速度滯后對氣相產(chǎn)生了阻力導(dǎo)致氣流流速減小。從以上可以看出，兩相流情況下的計算結(jié)果與純氣相時有較大不同，所以在進行噴管結(jié)構(gòu)溫度場數(shù)值計算時應(yīng)該以兩相流流場計算結(jié)果作為瞬態(tài)熱傳導(dǎo)模擬的初始值，而本文正是如此。

圖4 兩相流計算結(jié)果與純氣相計算結(jié)果對比

3.3 結(jié)構(gòu)溫度場結(jié)果分析

3.3.1 瞬時溫度分布

圖5是工作初期噴管喉部附近的溫度云圖，從圖中可以看出，由于喉襯材料導(dǎo)熱系數(shù)較大，炭/炭喉襯的溫度傳導(dǎo)速率明顯高于燒蝕層和絕熱層，喉襯與燒蝕層和絕熱層的溫差較大。而且時間越往后，可以發(fā)現(xiàn)炭/炭喉襯處熱量向兩側(cè)燒蝕層傳遞，它們的溫差逐步減小。圖6是t=2，3，4，5 s時噴管的溫度場分布情況，如圖6所示，隨著時間的推移，噴管壁面溫度逐漸增加，熱量在不斷地往固體域內(nèi)部傳遞，氣體與內(nèi)壁的對流換熱和固體域內(nèi)的熱傳導(dǎo)，使得噴管熱防護層在徑向呈現(xiàn)出清晰的溫度梯度，由于燒蝕層、絕熱層等的材料不同，其熱力學(xué)性能有較大差異，使得各固體域的溫度分布不同，炭/炭材料喉襯導(dǎo)熱性能好，而絕熱層和燒蝕層材料導(dǎo)熱系數(shù)較小，溫度上升緩慢，由于絕熱層對熱傳遞的阻隔，殼體溫度也上升緩慢，靠近噴管內(nèi)壁面的地方溫度較高，噴管壁面的熱量傳遞呈二維特性。

圖5 工作初期噴管喉部附近溫度云圖

為了進一步分析噴管徑向溫度分布，分別在圖7所示位置取截面，5s時截面上的熱防護層和殼體溫度分布如圖8所示，橫坐標均已無量綱化，橫坐標為0的點是流體域與固體域的交界面位置處。由圖8可以看出，每個截面的溫度均沿徑向逐漸降低，在炭/炭喉襯與絕熱層以及燒蝕層與絕熱層交界位置處曲線發(fā)生了轉(zhuǎn)折，曲線斜率不同，但由于燒蝕層與絕熱層材料導(dǎo)熱系數(shù)相差較小導(dǎo)致圖8(a)中效果不是太明顯，而炭、炭喉襯與絕熱層材料導(dǎo)熱系數(shù)相差比較大，溫度曲線斜率相差很大。殼體材質(zhì)為鋼，從表1可以知道其導(dǎo)熱系數(shù)比噴管其它部分大得多，所以殼體溫度保持在大于300K的小范圍溫度區(qū)間，圖8中殼體位置的溫度曲線斜率接近為零。另外，從圖中可以看出，噴管內(nèi)壁面處收斂段截面溫度最高，然后是喉部的三個截面，左邊截面溫度最大，右邊截面最小，且擴張段截面溫度小于喉部截面溫度，入口段截面溫度小于收斂段截面溫度。說明噴管內(nèi)壁面最高溫度位于收斂段中部至喉部前半段這一范圍內(nèi)，從圖6可以看出，這是因為該區(qū)域容易積聚大量顆粒，受熱嚴重，且顆粒的沖刷會對壁面造成不好的影響，加劇了該區(qū)域的燒蝕。

圖6 不同時刻噴管溫度場分布云圖

圖7 噴管各部位截面位置

圖8 5s時噴管各部位截面溫度曲線

3.3.2 對比分析

喉部是噴管的關(guān)鍵部位，且所處熱環(huán)境惡劣，所以有必要考慮喉部溫度分布及變化情況。圖9是t=1，3，5 s時以兩相流和純氣相流場計算結(jié)果為初始值的噴管喉部截面2上熱防護層和殼體的溫度分布曲線，圖中橫坐標原點為噴管內(nèi)流場與固體域的交界面位置，橫坐標為1的點為噴管外壁面位置，橫坐標都已無量綱化。從圖中可以看出，無論兩相流情況下還是純氣相情況下，隨著時間的推進，熱防護層溫度逐漸升高，在各個時刻，兩相流情況下，噴管喉部截面2上沿徑向溫度分布比純氣相條件下要高，但殼體位置效果還不太明顯。出現(xiàn)這種現(xiàn)象是因為顆粒比熱容較大，容易保持原來的溫度，而氣相在喉部時已開始膨脹降溫，由于溫差，噴管內(nèi)顆粒與燃氣之間必會發(fā)生熱交換，使得同一位置處氣相的溫度提高，同時粒子、燃氣與壁面之間也存在著大量的熱傳遞，進而熱量以導(dǎo)熱方式向熱防護層傳遞。顆粒的存在對噴管熱防護層、殼體的溫度分布產(chǎn)生了較大影響，使得其溫度值明顯提高，這對噴管的安全是不利的。因此，在對噴管進行熱分析時，不能忽略顆粒相的影響。

圖9 不同時刻噴管喉部溫度曲線對比

4 結(jié)論

相較于以往文獻采用純氣相流場結(jié)果進行后續(xù)數(shù)值計算，本文不僅采用了流固耦合的方法，并考慮了顆粒的影響，以噴管兩相流場結(jié)果作為固體壁面熱傳導(dǎo)初始條件對拉瓦爾噴管內(nèi)的流動換熱情況進行了模擬研究，得到了以下結(jié)論：

1)發(fā)動機工作過程中，燃氣熱量由噴管內(nèi)壁向外壁傳遞，固體壁面溫度逐漸升高，形成明顯的徑向溫度梯度，同時喉襯導(dǎo)熱系數(shù)大，不斷往兩側(cè)傳遞熱量，噴管壁面的熱量傳遞呈二維特性。

2)噴管固體壁面最高溫度位于內(nèi)壁面收斂段中部至喉部前半部分這一區(qū)間內(nèi)，此處容易積聚大量顆粒，在噴管中也最為脆弱。

3)顆粒的存在使得噴管內(nèi)流場結(jié)構(gòu)發(fā)生顯著變化，導(dǎo)致噴管熱防護層溫度升高，對熱防護層材料性能要求更高，因此，在進行噴管熱分析時應(yīng)考慮顆粒的作用，從而更加精確地對噴管熱防護層進行設(shè)計。