高 天，王安良

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

simulation

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

國(guó)內(nèi)，重力自循環(huán)系統(tǒng)直接應(yīng)用于工程領(lǐng)域的研究，工作較少且集中在結(jié)構(gòu)設(shè)計(jì)與實(shí)驗(yàn)方面[1]。國(guó)外，Podkorytov等[2]實(shí)驗(yàn)研究了毛細(xì)力對(duì)液氮自然對(duì)流循環(huán)過(guò)程中傳熱的影響。韓國(guó)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)域，如果液氮有一定的過(guò)冷度，效果更好；隨著表面熱流密度的升高，被冷卻壁面溫度也在升高。日本中部大學(xué)聯(lián)合美國(guó)MIT等多家單位[6-9]針對(duì)過(guò)冷液氮重力自循環(huán)系統(tǒng)在HTS領(lǐng)域的應(yīng)用開展了實(shí)驗(yàn)與理論研究工作。

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

(1)

式中，NuFC為強(qiáng)迫對(duì)流無(wú)量綱數(shù)。

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

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

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

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

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

(2)

式中，Pr為流體的普朗特?cái)?shù)。則壁溫可計(jì)算如下

(3)

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

(4)

流動(dòng)滿足質(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為干度。起泡前，基本可以認(rèn)為x等于0，起泡后，由于流體已達(dá)飽和狀態(tài)，可以認(rèn)為液體部分溫度Tf=Tsat。根據(jù)能量守恒，干度計(jì)算如下

(9)

其中，zFDB表示起泡點(diǎn)所在高度。

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

(10)

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

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

(11)

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

(12)

(13)

(14)

(15)

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

重力壓降與加速壓降計(jì)算如下：

(16)

(17)

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

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

(18)

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

pout=pa

(19)

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

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

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

2 實(shí)驗(yàn)裝置介紹與算法說(shuō)明

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

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

3 結(jié)果分析

3.1 模型正確性驗(yàn)證

以Baudouy的實(shí)驗(yàn)系統(tǒng)為對(duì)象進(jìn)行仿真計(jì)算，計(jì)算值與實(shí)驗(yàn)值的對(duì)比如圖3所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

參考文獻(xiàn)：

[1] 單巍巍，劉波濤，丁文靜,等.重力式自循環(huán)系統(tǒng)中熱沉結(jié)構(gòu)設(shè)計(jì)方法研究[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] 徐濟(jì)鋆.沸騰傳熱和氣液兩相流[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.