液氮流動沸騰數值模擬研究進展

邵雪鋒 李祥東 汪榮順

1 引言

低溫流體廣泛應用于各種低溫換熱設備中。在最近幾十年里，低溫流體的應用領域不斷擴展，不僅用作火箭發動機的氧化劑和燃料以及用來冷卻超導磁體和電動機，還被用作冷卻介質冷卻計算機硬件以增加計算速度和微電路密度。和其它低溫流體相比，液氮安全，無毒而且廉價，在實驗研究中經常被用來代替其它低溫介質。因此液氮的利用存在著廣闊的前景。由于測試手段有限等原因，液氮流動沸騰的實驗條件較為苛刻。早期研究者進行的實驗工作［1-3］僅僅關注一些宏觀或平均參數來分析低溫氣液兩相流傳熱，也即擬合實驗參數得到換熱系數的經驗關系式。隨著工業技術的快速發展，低溫系統對低溫流體的換熱提出了更高的要求，需要全面了解低溫兩相流動中的相關參數(溫度、壓力、流量、流速、含氣率，空泡份額等)，因為這些參數的變化直接影響流動的換熱與流動狀況。傳統的經驗關系式并不能給出局部參數的變化和分布，需要采用新的方法分析液氮的流動沸騰換熱。

隨著計算機技術的快速發展，計算能力不斷提高，計算流體力學(CFD)技術得以快速發展，并且逐漸發展成為一門比較成熟的學科。作為一種輔助設計和分析的手段，計算流體力學技術已經廣泛地應用于各種行業。特別在核能領域，計算流體力學已經成為一種重要的輔助分析手段，到目前為止已經有諸多代碼應用于水的兩相流動分析中，例如RELAP5、TRAC、CATHARE等。低溫流體和以水為代表的高沸點流體的沸騰在本質上有很多類似之處［4］，用于分析水流動沸騰的數學模型理論上也可用來分析液氮的流動沸騰。但是液氮等低溫流體的物性不同于高沸點流體，例如，液氮的飽和溫度以及接觸角要遠小于水。物性的差別將會對流動沸騰傳熱產生影響，不能簡單將水的分析結果應用到液氮的流動沸騰分析中。要將CFD技術引入到低溫領域，需要針對液氮流動沸騰的特點，建立合適的數學物理模型編制計算程序，并進行大量的實驗來驗證代碼的準確性。只有這樣才可準確的分析和預測流動沸騰過程中局部參數的變化和分布，從而用于液氮系統的優化設計和操作。

本文總結和分析液氮流動沸騰時所采用的數值模擬技術，指出這些技術中存在的問題，同時指出一些迫切需要研究的課題，為將最新的數值模擬技術引入液氮流動沸騰分析打下基礎。

2 液氮流動沸騰的數值模擬研究現狀及存在的問題

由于液氮流動沸騰的實驗系統比較復雜，針對液氮流動沸騰進行的實驗較為有限，液氮流動沸騰的數值模擬研究則更少。鄭正泉［5］采用簡化的漂移通量模型來描述管道中的低溫氣液兩相流的流動過程。和一般漂移通量模型中所用的連續方程和動量方程都是動態方程不同，他采用的模型中只有液體連續方程是動態方程，而氣相的連續方程以及氣相和液相的動量方程都是準平衡態。在該模型中考慮了氣體的可壓縮性，分析了不同流型，最后比較了相分布以及壓力分布的模擬計算結果與實驗結果。陸柳柳［6］等采用一維均相流模型對液氮流動過程進行了穩態數值模擬，得到了循環流動過程中液氮含氣率、壓力和溫度沿管程方向的分布情況，并對由于不同回路條件引起的不同結果進行了分析。Ishimoto［7］等采用基于非穩態漂移通量模型的控制方程描述了垂直通道內液氮流動沸騰的傳熱過程，詳細分析了液氮沸騰兩相流動中的二維結構。孫淑風［8］采用分相流動模型環行通道內的環狀流進行了數值模擬，該模型基于守衡定律，得到了液膜中的壓力梯度，傳熱系數速度以及溫度分布。傳熱系數的實驗值比計算結果偏高30%。上述研究者分別采用不同的數學模型針對液氮流動沸騰進行了數值模擬研究。他們都將數值模擬結果和實驗結果進行了比。通過比較表明這些數學模型在一定程度上可用于分析液氮的流動沸騰，證明了數值模擬方法可以作為實驗方法的一種補充。除Ishimoto外，其它研究者只給出空泡系數沿管程的分布，并未給出局部參數變化情況以及詳細空間分布。表明這些研究者所采用的數學模型不能更深入的研究液氮的流動沸騰以獲的局部參數。李祥東［9-14］將兩流體模型引入到低溫流體的沸騰流動分析中。他采用兩流體模型和CFD技術作為液氮流動沸騰過程研究的框架，根據實驗觀測結果建立適當的封閉方程，然后利用得到的理論模型實現對雙流體模型的修正和完善，同時借助CFD技術實現對垂直通道內液氮流動沸騰過程的預測。他的研究建立了液氮核態沸騰壁面上的傳熱傳質機理模型及流場內的相間傳輸模型，并采用以上模型作為封閉方程完成了對雙流體模型的初步修正，數值模擬的結果與實驗數據在一定范圍內吻合較好。

從上面可以看出，到目前為止只有少量的研究者采用有限的數學模型分析了液氮的流動沸騰。只有Ishimoto和李祥東分別采用漂移通量模型和兩流體模型得到液氮流場中局部參數的詳細分布。漂移通量模型沒有完全考慮到氣液兩相間的相互作用，和兩流體模型相比，它忽略了氣液兩相間傳輸的局部特征，一般被看作是一種簡化的兩流體模型。和均相模型以及分相模型相比，兩流體模型可以獲得更為詳細和準確的預測結果。將兩流體模型引入低溫流體的分析是一個很大的進步。但到目前為止兩流體模型僅用于分析泡狀流。并沒有研究者采用兩流體模型分析液氮流動沸騰中出現的其它流型或流型過渡區。在實際的低溫系統中，例如汽化器，液氮的流動沸騰過程中包括多種流型。如果采用傳統的設計方式(先根據流型判斷準則判斷并確定各流型，之后再針對各流型采用經驗關系式)，將帶來兩步誤差。

計算流體力學大師Patankar［15］指出，理論的預測出自一個數學模型的結果，而不是出自于一個實際的物理模型。準確的模擬液氮流動沸騰中的流型轉變區需要考慮各影響因素，建立準確的數學模型。已經有不少研究者采用數值模擬的方法研究水流動沸騰流型轉變區，但國內外幾乎沒有關于液氮沸騰流型轉變區的數值模擬研究，因此需要借鑒水的流動沸騰數學模型，根據液氮的特點建立準確的封閉方程。在液氮泡彈狀流型過渡區，明顯的存在著氣泡破裂和聚合等現象，在兩流體模型中考慮氣泡的破裂和聚合現象，將大大拓展兩流體模型的使用范圍。

3 可用于液氮流動沸騰數值模擬研究的數學模型

在兩相流中，最重要的基本幾何參數是空泡系數，相界面面積濃度，當地氣泡尺寸。空泡系數表征氣液相的分布，它是換熱系統中熱流體設計的必需參數。兩流體模型中的很多參數，比如氣相在過冷液相中的冷凝速度，相間曳力等是相界面面積濃度或當地氣泡尺寸的函數，因此相界面面積濃度和當地氣泡尺寸強烈的影響相間的動量、質量和能量傳輸，是兩流體模型中的必需參數。在目前的兩流體模型中，當地氣泡尺寸通常由經驗關系式確定，相界面面積濃度則由空泡系數和當地氣泡尺寸的關系得到。因此，準確的確定當地氣泡尺寸顯得格外重要。Tu等［16］采用兩流體模型分析低壓水流動沸騰，發現數值模擬得到的徑向當量氣泡直徑以及氣液相速度分布與實驗結果存在較大偏差。這是因為在兩流體模型中，過冷液體中關于氣泡尺寸的經驗關系式只能用來預測宏觀沸騰現象，并不能用它們準確的預測實驗中可以觀察到的氣泡聚合和破裂等微觀現象。為了更準確的模擬液氮流動沸騰，需要考慮氣泡之間的相互作用以建立更為準確的氣泡直徑分布函數。也即需要在兩流體模型的基礎上建立一個關于氣泡直徑或者相界面面積濃度的補充方程，在這個補充方程中考慮氣泡的破裂和聚合等因素，從而完善兩流體模型。群體平衡模型的成熟為補充方程的建立奠定了基礎。

群體平衡模型描述的是系統中某些實體的數量的平衡關系，對于不同的系統，實體可以是固體顆粒、液滴、氣泡或者細胞等。對于沸騰兩相流系統，這些實體就是流場中的氣泡。對于某一給定的氣液兩相流空間，除了不斷有氣泡進入或離開該空間外，空間內的氣泡還會由于各種原因(如聚合、破裂以及相變等)不斷產生或消亡。大量的研究者基于群體平衡模型建立了兩流體模型的補充方程，這些補充方程可分為3類:

第一類補充方程將相界面面積濃度作為一種傳輸量，根據氣泡數量密度和相界面面積濃度的關系建立相界面面積傳輸方程，在相界面面積方程的源項中考慮氣泡破裂和聚合等因素。Ishii等［17］認為由于兩相流動中存在較多的可變形的移動界面，相間傳輸機理比較復雜，因此，兩流體模型中最大的不足是很難準確的封閉相間傳輸項。相間傳輸項和相界面面積濃度以及驅動力成正比。驅動力則取決于局部傳輸機理，包括相界面附近的湍流等。因此，相界面面積濃度是一個必需的量。

第二類補充方程根據氣泡直徑的大小將氣泡分為若干組(N組)，通過氣泡數量密度和氣泡尺寸的關系將群體平衡模型與兩流體模型結合起來，可在較大范圍預測兩相流場內的氣泡直徑分布。這類方程組(稱為多尺寸組方程)中需要計算每組氣泡的連續方程，液相的連續方程以及氣液相的動量和能量方程與兩流體模型并無區別。在氣泡的連續方程中可以考慮氣泡的破裂和聚合。Tu等采用此類方程模擬了低壓水的過冷流動沸騰，數值模擬結果和實驗結果吻合的較好。

第三類補充方程將平均氣泡體積作為傳輸量，基于群體平衡模型建立了平均氣泡體積的傳輸方程。Lehr等［18］采用該類補充方程分析了鼓泡床中的氣泡特征。

3類補充方程都可以考慮氣泡之間的相互作用，可被用于分析流型轉變區域。加入此類補充方程將使得兩流體模型的使用范圍大大拓展。前兩種方法都被認為是研究兩相流動有前途的方法，第三類方法也引起了一些研究者的注意。國內外學者針對水的流動沸騰進行了廣泛而深入的研究，相界面傳輸方程和多尺寸方程得到了廣泛的關注。但是國內外很少有研究者將這些模型引入到低溫流體的流動沸騰中。準確數學模型的缺乏限制了數值模擬技術在液氮流動沸騰中的應用，因此，在低溫領域中引入這些數學模型是必需的。李祥東成功的將兩流體模型引入到液氮的流動沸騰分析表明準確的數學模型是數值模擬成功的基礎。關于水的實驗證明兩流體模型并非完全的機理模型，需要針對液氮流動沸騰的特點建立準確的補充方程從而完善兩流體模型。

4 液氮流動沸騰數值模擬研究展望

由于液氮的物性和水存在著較大的區別，液氮的流動沸騰過程中存在著更加復雜的問題。這些問題包括液氮兩相流動中存在較多可變形的移動界面，以及液氮物性的不連續性和氣液界面處存在著復雜流場等。使得準確建立液氮流動沸騰的數學模型并不容易。

Ishii指出在兩相流動中，至少有4種尺度是非常重要的。這4種尺度包括系統尺度、連續介質假設所要求的宏觀尺度、局部結構所要求的中等尺度和分子結構有關的微觀尺度。其中系統尺度主要關注系統的瞬態效應以及組件的相互作用，宏觀尺度包括相界面結構以及質量、動量和能量的傳輸，中等尺度強調湍流對動量和能量的影響以及相界面處質量、動量和能量的傳輸。要準確的采用數值模擬方法預測液氮的流動沸騰，需要在這4種尺度上進行廣泛而深入的研究。相界面傳輸方程和多尺寸組方程所考慮的氣泡聚合和破裂現象屬于微觀度，需要深入和集中的研究。和多尺寸組方程相比，相界面傳輸方程求解的只是相同控制容積中基于平均分布的相界面面積和空泡系數得到的平均氣泡尺寸，因此相界面面積方程可看作是多尺寸組方程的簡化形式。完善兩流體模型的第一步工作是將多尺寸組方程引入到液氮的流動沸騰分析中，因為準確分析流場中氣泡尺寸的分布是建立準確相間傳輸模型的基礎。在多尺寸組方程中，壁面傳熱傳質模型依然是求解成功與否的關鍵模型之一。李祥東的研究表明，豎直通道內液氮流動沸騰中的活化核心密度，氣泡脫離頻率可分別采用Kirichenko［19］，Cole［20］提出的公式計算。氣泡脫離直徑則可采用Kirichenko［21］提出的公式。數值模擬中采用這類經驗關系式得到的數據可在一定范圍內與實驗數據吻合較好。但為了拓展兩流體模型的使用范圍和增加兩流體模型的預測精度，需要采用更加機理的方法減少經驗關系式帶來的不確定性，這需要詳細的研究氣泡的核化，生長，滑移和脫離。因此，完善兩流體模型的第二步工作是準確的建立液氮流動沸騰中氣泡的受力分析模型，得到準確的氣泡脫離直徑公式和氣泡脫離頻率公式。另外目前兩流體模型中的壁面熱量拆分模型只適合于氣泡迅速離開加熱壁面的過冷沸騰，該模型中并未考慮氣泡的滑移，第三步工作需要修正目前兩流體模型中的壁面熱量拆分模型。

目前，第一部分工作已初步完成。本文針對幾何原型為一內徑6 mm，外徑20 mm，長1.5 m的豎直環行管道進行了數值計算。過冷度為5 K的液氮從環形通道下部垂直向上流動，入口的液相速度為0.144 m/s;內壁均勻加熱，加熱量為18 000 W/m2;出口壓力為0.7 MPa。由于壁面熱流量均一，為了充分利用環形幾何形狀，只將其1/4作為模擬區域。采用適體網格系統在環形通道中產生三維網格，最后得到總14(徑向)×30(高度)×3(沿周長方向)個控制容積。活化核心密度，氣泡脫離頻率分別采Kirichenko，Cole提出的公式計算，氣泡脫離直徑則可采用Kirichenko提出的公式。

圖1給出了3個軸向位置的平均當量氣泡直徑的徑向分布，從數值模擬的結果可以看出，當量氣泡直徑的最大值出現在管道中而不是加熱壁面附近，數值分析的結果與液氮流動沸騰的可視化實驗相符。沿流動方向，當量平均氣泡直徑的最大位置慢慢趨近于管道中央。在出口處最大氣泡直徑出現在管道正中間，表明此時在管道中間形成了大的彈狀氣泡。

圖1 平均當量氣泡直徑的徑向分布Fig.1 Local mean radial profiles of Sauter diameter at different height

圖2 給出了空泡系數徑向分布，在高度為0.925 m處，空泡系數的最大值出現在加熱壁面，沿流動方向空泡系數的最大值漸漸遠離加熱壁面，接近出口處時最大空泡系數出現在管道中間，同樣表明此時流動處于彈狀流狀態。

圖2 空泡系數的徑向分布Fig.2 Local mean radial profiles of void fraction at different height

圖3 和圖4給出了沿流動方向分布的平均當量氣泡直徑和空泡系數的分布。從圖3可以看出，在某一軸向位置平均當量氣泡直徑開始超過氣泡起飛直徑，此后氣泡的聚合作用越來越明顯，形成較大的氣泡，直到接近出口處產生彈狀氣泡。圖4的空泡系數分布可以看出，沿著流動方向，剛開始空泡系數的增長速度較平緩，到某一位置時，空泡系數的增長速度急劇增加。該結果符合過冷流動沸騰的特點。

圖3 平均當量氣泡直徑沿流動方向的分布Fig.3 Axial distribution of mean Sauter diameter

圖4 空泡系數沿流動方向的變化Fig.4 Axial distribution of mean void fraction

這4個圖表明多尺寸方程可以通過氣泡尺寸和空泡系數的變化和分布捕捉到氣泡破裂和聚合的相關信息，還可以判斷流動沸騰處于何種流型，因此多尺寸組方程可以擴展兩流體模型的使用范圍。

5 結語

只有少量的研究者采用數值模擬的方法研究了液氮的流動沸騰，只有Ishimoto和李祥東得到了液氮流場中局部參數的詳細分布。Ishimoto的方法中忽略了氣液兩相間傳輸的局部特征，兩流體模型中的氣泡直徑分布采用經驗關系式，這些因素限制了漂移通量模型和兩流體模型的使用范圍。拓展兩流體模型的范圍需要考慮流場中氣泡的破裂和聚合因素。數量密度方程的成熟為兩流體模型的完善奠定了基礎。建立在數量密度基礎上的相界面傳輸方程和多尺寸組方程被認為是很有前途的研究兩相流的方法，受到了廣泛的關注。將這些模型引入到低溫流體的沸騰分析中，可以更為準確的分析流場中的局部數據，判斷流型轉變。從而可以輔助低溫換熱設備的設計與分析。

相界面傳輸方程只是多尺寸組方程的一種簡化形式，因此，完善兩流體模型需要首先將多尺寸組方程引入到低溫流體的沸騰分析中。兩流體模型中，核態沸騰壁面傳熱傳質模型需要繼續加以完善，建立更為機理的氣泡起飛直徑公式和氣泡脫離頻率公式。同時需要考慮氣泡的滑移進一步改善壁面熱量拆分模型。本文已經初步完成第一步工作，即將多尺寸組方程引入到液氮的流動沸騰分析中。結果表明多尺寸組方程能考慮實驗觀察到的氣泡破裂和聚合現象，并且能對流型轉變區域進行模擬。低溫流體流動沸騰的實驗數據相當有限，需要進行更多的實驗獲的準確的局部數據，進一步修正氣泡起飛直徑，氣泡脫離頻率等公式，從而驗證代碼的可靠性。

1 Klimenko V V.Heat transfer intensity at forced flow boiling of cryogenic liquids in tubes［J］.Cryogenics，1982，22(6):569-576.

2 Steiner D，Schlunder E U.Heat transfer and pressure drop for boiling nitrogen flowing in a horizontal tube1saturated flow boiling［J］.Cryogenics，1976，16(3):387-399.

3 Steiner D，Schlunder E U.Heat transfer and pressure drop for boiling nitrogen flowing in a horizontal tube2.saturated flow boiling［J］.Cryogenics，1976，16(3):457-764.

4 Barron R F.Cryogenic Heat Transfer［M］.Taylor＆Francis，Philadelphia，1999.

5 鄭正泉.低溫汽液兩相流動態模擬及實驗研究［J］.低溫與超導，1999，27(4):21-26.

6 陸柳柳，匡 波，陳 宏.低溫氣液兩相流數值計算［J］.低溫與超導，2005，33(1):27-31.

7 Ishimoto J，Oike M，Kamijo K.Two-dimensional numerical simulation of boiling two-phase flow of liquid nitrogen［J］.Transactions of the Japan Society for Aeronautical and Space Sciences，2000，43(141):114-121.

8 Chen L F，Wu Y Y，Liu Y Z，et al..The study on performance of heat transfer in a quasi-annular flow condenser-evaporator［J］.Cryogenics，1999，39(3):209-216.

9 李祥東，汪榮順，石玉美.垂直通道內低溫液體過冷流動沸騰傳熱的數值預測模型［J］.低溫工程，2006，34(1):6-11.

10 李祥東.豎直通道內液氮流動沸騰的雙流體模型及沸騰兩相流不穩定性研究［D］.上海:上海交通大學大學，2007.

11 李祥東，汪榮順，石玉美.低溫液體核態流動沸騰傳熱的機理模型［J］.低溫與超導，2006，34(3)，168-171.

12 Li X D，Wang R S，Gu A Z.Numerical simulation of flow boiling of cryogenic liquids in tubes［C］.In:20th International Cryogenic Engineering Conference，Beijing，2004.

13 李祥東，汪榮順，顧安忠.低溫液體流動沸騰數值計算中的動量模擬［J］.低溫與超導，2004，32(4):53-58.

14 李祥東，汪榮順，顧安忠.低溫液體流動沸騰數值計算中的相間傳熱模型［J］.低溫與超導，2005，33(2):26-29.

15 Patankar SV.Numerical heat transfer and fluid flow［M］.McGraw:McGraw-Hill，1980.

16 Yeoh G Y，Tu JY.Population balance modeling for bubbly flows with heat and mass transfer［J］.Chemical Engineering Science，2004，59(15):3125-3139.

17 Ishiii M，Hibiki T.Thermo-fluid dynamics of two-phase flow［M］.Springer，New York，2006.

18 Lehr F，Mewes D.A transport equation for the interfacial area density applied to bubble columns［J］.Chemical Engineering Science，2001，56(3):1159-1166.

19 Kirichenko Y A，Dolgoy M L，Levchenko N M.A study of the boiling of cryogenic liquids［J］.Heat Transfer-Soviet Research，1976，8(4):63-72.

20 Cole R.A photographic study of pool boiling in the region of the critical heat flux［J］.AIChE Journal，1960，6(3):533-542.

21 Kirichenko Y A，Slobozhanin L A，Shcherbakova N S.Analysis of quasi-static conditions of boiling onset and bubble departure［J］.Cryogenics，1983，23(2):110-112.