液氮重力自循環(huán)系統(tǒng)的仿真及設(shè)計

高 天，王安良

(北京航空航天大學(xué)宇航學(xué)院， 北京 100191)

simulation

液氮重力自循環(huán)系統(tǒng)作為一種低溫?zé)岢料到y(tǒng)，其工作原理如圖1所示，其主要物理機制是熱虹吸流和自然對流沸騰。流體在傳熱的作用下流動，熱沉管內(nèi)液氮因壁面受熱而產(chǎn)生自然對流或部分汽化，致使上升段與下降段產(chǎn)生重力壓差，且流動同時受摩擦阻力的限制，流動達到穩(wěn)定時重力壓差與系統(tǒng)總阻力達到動態(tài)平衡。

國內(nèi)，重力自循環(huán)系統(tǒng)直接應(yīng)用于工程領(lǐng)域的研究，工作較少且集中在結(jié)構(gòu)設(shè)計與實驗方面[1]。國外，Podkorytov等[2]實驗研究了毛細力對液氮自然對流循環(huán)過程中傳熱的影響。韓國Choi等[3]、Kim和Chang[4]以及Chang[5]課題組把液氮重力自循環(huán)制冷方法用于高溫超導(dǎo)研究領(lǐng)域，他們的研究結(jié)果表明，液氮重力自循環(huán)系統(tǒng)可用于高溫超導(dǎo)(high temperature superconductor,HTS)領(lǐng)域，如果液氮有一定的過冷度，效果更好；隨著表面熱流密度的升高，被冷卻壁面溫度也在升高。日本中部大學(xué)聯(lián)合美國MIT等多家單位[6-9]針對過冷液氮重力自循環(huán)系統(tǒng)在HTS領(lǐng)域的應(yīng)用開展了實驗與理論研究工作。

法國的Baudouy[10]開展了低熱流密度下液氮自然循環(huán)系統(tǒng)的實驗研究，并采用Gr數(shù)，Nu數(shù)分析了兩相流區(qū)段的換熱特性，并給出了一個修正的計算自然對流關(guān)聯(lián)式的公式

(1)

式中，NuFC為強迫對流無量綱數(shù)。

由于重力自循環(huán)系統(tǒng)涉及到氣液兩相流動，物理機制較復(fù)雜，因而國內(nèi)外的研究大多用CFD軟件單獨模擬熱沉部分的流動，計算耗時較長，且不能得到系統(tǒng)的整體特性。

本文提出了針對重力自循環(huán)系統(tǒng)的一維仿真模型，研究系統(tǒng)在不同加熱功率下的總體流量及壓降特性。作者研究了不同系統(tǒng)結(jié)構(gòu)參數(shù)和邊界條件下系統(tǒng)的性能，可為一般重力自循環(huán)系統(tǒng)的設(shè)計提供參考。

1 數(shù)學(xué)模型及邊界條件

由于計算工況中干度較小[見圖3(c)]，因而采取一維均相流動模型，將氣液混合物看作擁有平均物理性質(zhì)的贗流體。這里進行簡化，認為沸騰起始點(ONB)點之前，流體干度x等于0，ONB點后流體處于飽和狀態(tài)。

ONB點位置與壁溫密切相關(guān)，因而需要選取合適的換熱系數(shù)得到壁溫，采用式(2)計算[11]

(2)

式中，Pr為流體的普朗特數(shù)。則壁溫可計算如下

(3)

其中，l為管徑d。得到壁溫后，即可根據(jù)式(4)確定ONB點位置

(4)

流動滿足質(zhì)量守恒方程

(5)

(6)

其中：ρm、ρl、ρg分別為混合密度、液氮密度、氣氮密度；A為管截面積；v為流體速度。

能量守恒方程

(7)

其中：q為熱流密度；D為管徑；G為系統(tǒng)流量；H為流體的比焓；v為流體的速度；g為重力加速度；z為豎直高度；θ為管傾斜角

H=cpl·Tf·(1-x)+r·x

(8)

式中：cpl為流體比熱容；Tf為流體溫度；r為汽化潛熱；x為干度。起泡前，基本可以認為x等于0，起泡后，由于流體已達飽和狀態(tài)，可以認為液體部分溫度Tf=Tsat。根據(jù)能量守恒，干度計算如下

(9)

其中，zFDB表示起泡點所在高度。

流體溫度可由式(7)推出

(10)

流動達到穩(wěn)定時，系統(tǒng)的重力壓差與流動阻力達到動態(tài)平衡，因而正確計算系統(tǒng)各部分的壓降至關(guān)重要。系統(tǒng)的阻力包括兩部分：沿途阻力損失與局部阻力損失。

沿途壓降可根據(jù)動量守恒方程計算

(11)

兩相流的壓降可分為3部分：重力壓降、摩擦壓降與加速壓降。壓降根據(jù)Martinelli-Chisholm關(guān)聯(lián)式[12]計算，鄧東等證實此關(guān)聯(lián)式適用于液氮[13]。即：

(12)

(13)

(14)

(15)

其中：μl、μg分別為液氮與氣氮的動力黏度；ρl、ρg分別為液氮與氣氮的密度；x為干度。

重力壓降與加速壓降計算如下：

(16)

(17)

其中：G為系統(tǒng)流量；A為管道截面積。

局部壓力損失系數(shù)為：

(18)

系統(tǒng)在流動達到穩(wěn)定時，整體的動力壓頭與總阻力是相等的，因而出口處流體的靜壓強應(yīng)等于外界環(huán)境壓強，否則系統(tǒng)流量仍會波動，如圖1所示，即有：

pout=pa

(19)

計算中，出口處壓強pout根據(jù)入口壓強與沿途壓降相加得到，環(huán)境壓強pa是已知的。

式(3)、式(5)、式(9)、式(10)、式(11)、式(19)六式組成的方程組中，在zzFDB時均只有5個獨立方程起作用，其中共涉及G、x、p、Tw、Tf五個獨立變量，因此方程組的解是唯一的。在外界環(huán)境壓力、系統(tǒng)參數(shù)等條件給定的情況下，可以聯(lián)立求解得到系統(tǒng)穩(wěn)態(tài)時的各變量值。若壁溫不符合設(shè)計(比如小于100 K)要求，則進一步調(diào)整系統(tǒng)結(jié)構(gòu)參數(shù)以滿足設(shè)計目標(biāo)。

本研究采用Baudouy的實驗結(jié)果[10]驗證模型的正確性，相關(guān)邊界條件如下：①pa取1個大氣壓；② 熱流密度q變化范圍為2～12 kW/m2；③ 管長、管徑等結(jié)構(gòu)參數(shù)見第2節(jié)。

2 實驗裝置介紹與算法說明

本研究首先采用Baudouy的實驗數(shù)據(jù)驗證模型的正確性，其實驗裝置工作狀態(tài)如圖1所示，實驗裝置由循環(huán)回路與外部的絕熱罩兩部分組成。循環(huán)回路的上升段與下降段通過U型與液氮儲箱連接，下降段管徑為40 mm，上升段管徑為10 mm。液氮儲箱直徑為0.45 m，高0.3 m。上升段長1.45 m，其中受熱部分位于上升段的下部，長0.95 m，非受熱段長0.5 m。實驗通過安置文氏流量計測量系統(tǒng)流量，文氏流量計長0.4 m，入口直徑為40 mm，喉部直徑為10 mm。實驗裝置的具體結(jié)構(gòu)參數(shù)見文獻[10,14]。

計算時，初設(shè)流量后進行迭代對數(shù)學(xué)模型進行求解，其流程如圖2所示，算法在Matlab中實現(xiàn)。

3 結(jié)果分析

3.1 模型正確性驗證

以Baudouy的實驗系統(tǒng)為對象進行仿真計算，計算值與實驗值的對比如圖3所示。

由圖3可見，隨著熱流密度增大，計算流量先增大，達到極大值后開始減小，流量的相對誤差絕大部分均在±10%以內(nèi)；受熱段總壓降隨著熱流密度的升高不斷減小，壓降的相對誤差絕大部分在±10%以內(nèi)。流量與壓降隨熱流密度的變化趨勢均與實驗結(jié)果相符，且相對誤差較小，驗證了仿真模型的正確性。

3.2 管徑對系統(tǒng)工作特性的影響

熱沉管徑D對流體的流動與換熱特性有著重要影響。在其他參數(shù)不變時，研究不同管徑對系統(tǒng)工作特性的影響。

圖4(a)中，不同管徑下，流量隨熱流密度升高的變化趨勢是一致的，都是先增大后減小。在同一熱流密度下，大管徑對應(yīng)著較大的流量。而其他參數(shù)相同時，大管徑所對應(yīng)的單位長度上的換熱量較大，因而對應(yīng)著較大的干度，所以重力壓差較大；另一方面，摩擦壓降與管徑成反比關(guān)系，兩種作用使得大管徑系統(tǒng)流量更大。

隨著管徑增大，使系統(tǒng)達到最大流量的臨界熱流密度也在逐漸增加大。在管徑為0.005 m時，臨界熱流密度在 2 000 W/m2左右，而在管徑為0.02 m時，臨界熱流密度已經(jīng)超過了10 000 W/m2。圖4(c)展示了不同管徑下系統(tǒng)出口處的干度值，在管徑較小時，受熱管內(nèi)干度較大，摩擦阻力較大，使得流量在較小的熱流密度下達到極大值。

3.3 受熱段長度對系統(tǒng)工作特性的影響

液氮在受熱段吸熱汽化，受熱段長度決定著液氮總吸熱量，進而影響系統(tǒng)內(nèi)氣相組分大小。在Baudouy的實驗裝置基礎(chǔ)上，保持上升段長度為1.45 m不變，取受熱段起始點高度H0=0，改變受熱部分的長度，計算結(jié)果如圖5所示。

對比圖5發(fā)現(xiàn)，受熱長度對系統(tǒng)流量的影響與分受熱長度大小有關(guān)。在受熱長度較小(Lh/L<0.4)時，受熱長度增加使系統(tǒng)流量增大；在受熱長度滿足0.6>Lh/L>0.4時，受熱長度增加在低熱流密度下使系統(tǒng)流量增大，在高熱流密度下使系統(tǒng)流量減小。在受熱長度較大(Lh/L>0.4)時，受熱長度增加反而使系統(tǒng)流量減小。

上升段內(nèi)的平均干度受熱流密度與加熱長度的共同影響。存在一個臨界的平均干度xc，在平均干度小于xc時，上升段內(nèi)重力壓降起主要作用，在平均干度大于xc時，摩擦壓降替代重力壓降起主要作用。在受熱長度較小時，上升段內(nèi)干度較小，因而受熱長度增加使平均干度增加，上升段與下降段的重力壓差隨之增大，導(dǎo)致系統(tǒng)流量增大。在受熱長度較大時，上升段內(nèi)平均干度較大，而摩擦壓降隨干度的增加而增加，因而系統(tǒng)流量反而隨受熱長度的增大而減小。在受熱長度處于兩者之間時，受熱長度對流量的影響在高熱流密度區(qū)與低熱流密度區(qū)截然相反。

3.4 受熱段位置對系統(tǒng)工作特性的影響

加熱段長度一定時，加熱位置不同會改變流動的阻力特性。在Baudouy的實驗裝置基礎(chǔ)上，保持上升段長度為1.45 m不變，取Lh/L=0.2，改變受熱段起始點高度H0，結(jié)果如圖6所示。

總體上，受熱位置上升會使得系統(tǒng)流量減小。分析認為，系統(tǒng)循環(huán)流動的動力在于流體受熱汽化后在下降段與上升段之間產(chǎn)生了密度差，密度差導(dǎo)致重力壓差。而兩相區(qū)的長度越大，重力壓差越大。因而，加熱段位于下部時，液氮汽化后流經(jīng)的長度更大，使得上升段與下降段的壓強差更大，流量更大。圖6(b)證實了這一點。

干度與受熱長度近似成正比關(guān)系，只有受熱長度達到一定值時，干度才能產(chǎn)生較大影響。因此可以推斷：受熱長度越大，受熱位置上升使系統(tǒng)流量減小的趨勢越明顯。

3.5 上升段總管長對系統(tǒng)工作特性的影響

在Baudouy實驗裝置中，非受熱段處于受熱段的上方，因而上升段與下降段的重力壓差由受熱段與非受熱段共同提供，因此上升段的總管長對系統(tǒng)工作特性有著重要影響。Baudouy的實驗裝置中，系統(tǒng)上升段總長度為1.45 m，在實際工程中管長有時候會達到數(shù)十米乃至百米[7]。為研究管長對系統(tǒng)工作特性的影響，這里保持系統(tǒng)的受熱段長度為0.95 m不變，將上升段的管長Lup增大，同時相應(yīng)地增大下降段的管長Ldown，計算結(jié)果如圖7所示。

總管長的增加對系統(tǒng)流量的影響需要分低熱流密度區(qū)與高熱流密度區(qū)討論。熱流密度較低時，管長增加對系統(tǒng)流量的影響不大，熱流密度較高時，管長增加會使得系統(tǒng)流量有一定程度的增大。分析認為，在熱流密度較低時，上升段內(nèi)主要是單相流，這使得上升段與下降段間壓差減小，使得系統(tǒng)在較小的流量下即可實現(xiàn)換熱和動量平衡。在熱流密度達到3 000 W/m2以上時，由于受熱管內(nèi)干度與熱流密度近似成正比關(guān)系，受熱段內(nèi)干度足夠大，使得上升段內(nèi)重力壓降起主要作用，因而管長增加會使得系統(tǒng)流量增大。

4 結(jié)論

1) 仿真模型的仿真結(jié)果與Baudouy的實驗數(shù)據(jù)相當(dāng)一致，說明該模型能對一般重力自循環(huán)系統(tǒng)的穩(wěn)態(tài)特性進行正確計算。

2) 仿真研究發(fā)現(xiàn)，系統(tǒng)的結(jié)構(gòu)參數(shù)如管徑等對系統(tǒng)的流量與壓降特性有重要影響。

3) 本文研究工作可用于估算已有的重力自循環(huán)熱沉系統(tǒng)的性能，也可用于設(shè)計新系統(tǒng)。

參考文獻：

[1] 單巍巍，劉波濤，丁文靜,等.重力式自循環(huán)系統(tǒng)中熱沉結(jié)構(gòu)設(shè)計方法研究[J].航天器環(huán)境工程,2010,27：489-492.

[2] DMITRIY PODKORYTOV,LEONID TIMKIN,VLADIMIR CHECHOVICH.The investigation of flow characteristics and heat transfer at the nitrogen natural circulation in capillary[J].Cryogenics,1992,32(1):316-319.

[3] YEON SUK CHOI,HO-MYUNG CHANG,STEVEN W.Van Sciver.Performance of extended surface from a cryocooler for subcooling liquid nitrogen by natural convection[J].Cryogenics,2006,46(5):396-402.

[4] KIM M J,CHANG H.M.Natural circulation loop of subcooled liquid nitrogen[C]//Advances in Cryogenic Engineering.Chattanooga,2008.

[5] CHANG HO-MYUNG,KIM MIN-JEE,SIM JUNG WOOK.A compact cryocooling system with subcooled liquid nitrogen for small HTS magnets[C]//20th International Conference on Magnet Technology.Philadelphia,2007.

[6] YU.IVANOV,A.RADOVINSKY,A.ZHUKOVSKY,et al.Compact counter-flow cooling system with subcooled gravity-fed circulating liquid nitrogen[J].Physica C:Superconductivity,2010,470 (20):1895-18989.

[7] IVANOV Y,RADOVINSKY A,ZHUKOVSKY A,et al.A compact cooling system for HTS power cable based on thermal siphon for circulation of LN2[C]//Transactions of the Cryogenic Engineering Conference—cec:Advances in Cryogenic Engineering.American Institute of Physics,2010:865-870.

[8] YAMAGUCHI S,IVANOV Y,SUN J,et al.Experiment of the 200-Meter Superconducting DC Transmission Power Cable in Chubu University[J].Physics Procedia,2012,36:1131-1136.

[9] YURY IVANOV,HIROFUMI WATANABE,MAKOTO HAMABE,et al.Observation of the thermosiphon effect in the circulation of liquid nitrogen in HTS cable cooling system[J].Physics Procedia,2012,27:368-371.

[10] BAUDOUY.Experimental study of a nitrogen natural circulation loop at low heat flux[C]//Joint Cryogenic Engineering Conference.Tucson,2009.

[11] 李友榮.高等傳熱學(xué)[M].北京：科學(xué)出版社，2013.

[12] 徐濟鋆.沸騰傳熱和氣液兩相流[M].北京:原子能出版社，2001.

[13] DENG D,XIE S W,WANG R S.Two-phase flow pressure drop of liquid nitrogen boiling in the straight section downstream of U-bend[J].Journal of Shanghai Jiaotong University (Science),2014,19(4):495-50187.

[14] BENKHEIRA L,BAUDOUY B,SOUHAR M.Heat transfer characteristics of two-phase He I (4.2 K) thermosiphon flow[J].International Journal of Heat and Mass Transfer,2007,50(17):3534-3544.