李 琳，姜文全，李婷婷，楊 帆，宿詩(shī)雨，石杰峰

（1.遼寧石油化工大學(xué) 機(jī)械工程學(xué)院，遼寧 撫順113001；2.臨沂市自然資源和規(guī)劃局，山東 臨沂276301；3.遼寧石油化工大學(xué) 石油天然氣工程學(xué)院，遼寧 撫順113001）

液化天然氣（LNG）作為清潔能源，其進(jìn)口量日益增加。由于其低溫特性，應(yīng)用過(guò)程中的傳熱過(guò)程是必要的一環(huán)[1-2]。已有學(xué)者對(duì)不同應(yīng)用背景下的超臨界壓力LNG（或甲烷）的傳熱過(guò)程進(jìn)行了研究。韓昌亮等[3]、王博杰等[4]、李仲珍等[5]、潘杰等[6]對(duì)接收站氣化器管內(nèi)LNG的傳熱過(guò)程進(jìn)行了實(shí)驗(yàn)和數(shù)值研究，分析了壓力、質(zhì)量流量、邊界條件等對(duì)換熱特性的影響，發(fā)現(xiàn)管內(nèi)LNG局部傳熱系數(shù)在擬臨界點(diǎn)附近達(dá)到最大；王亞洲等[7]針對(duì)航空領(lǐng)域低溫推進(jìn)過(guò)程超臨界壓力LNG的傳熱過(guò)程進(jìn)行了研究，發(fā)現(xiàn)臨界區(qū)域附近物性的劇烈變化使傳熱過(guò)程惡化；李瑋哲等[8]、賈丹丹等[9]對(duì)印刷式電路板內(nèi)LNG的傳熱過(guò)程進(jìn)行了數(shù)值研究，討論了通道截面形狀和流道彎曲角對(duì)換熱的影響；H.F.Gu等[10]、郭占魁等[11]對(duì)水平和豎直管道內(nèi)超臨界壓力LNG的傳熱過(guò)程進(jìn)行了研究，均發(fā)現(xiàn)在高熱流、低質(zhì)量流量下會(huì)發(fā)生傳熱惡化現(xiàn)象；D.Huang等[12]對(duì)超臨界壓力水、二氧化碳和碳?xì)淙剂蟼鳠岬膶?shí)驗(yàn)關(guān)聯(lián)式和惡化判據(jù)進(jìn)行了總結(jié)，發(fā)現(xiàn)不同流體關(guān)聯(lián)式和判據(jù)的差異較大；白萬(wàn)金等[13]、Q.Zhang等[14]、B.S.Shiralkar等[15]、曲默豐等[16]發(fā)現(xiàn)，低質(zhì)量流率下會(huì)出現(xiàn)異常傳熱行為；R.N.Xu等[17]及其他學(xué)者[18-28]研究了浮升力和熱加速因子；S.Yildiz等[29]及范辰浩等[30]研究了管徑對(duì)超臨界壓力流體傳熱行為的影響；李志剛等[31]及代寶民等[32]研究了物性變化對(duì)超臨界壓力流體傳熱行為的影響，并進(jìn)行了機(jī)理分析；王彥紅等[33]在超臨界壓力航空煤油的流動(dòng)和傳熱研究中發(fā)現(xiàn)了傳熱不穩(wěn)定性。

綜上所述，超臨界壓力流體因擬臨界區(qū)物性變化劇烈，其傳熱過(guò)程受流動(dòng)工況、流動(dòng)方向、流動(dòng)通道等因素的制約。與超臨界壓力水、二氧化碳和碳?xì)淙剂系膫鳠嵯啾?，超臨界壓力LNG傳熱過(guò)程由于其低溫特性，在軸向和徑向上存在較大溫差，其傳熱規(guī)律和機(jī)理在很大程度上不同于超臨界壓力水和二氧化碳。針對(duì)此問題，本文進(jìn)行了大溫差下超臨界LNG豎直管道內(nèi)的傳熱和不穩(wěn)定性分析，以期為L(zhǎng)NG氣化器、推進(jìn)過(guò)程及以LNG為冷源的傳熱過(guò)程提供理論基礎(chǔ)。

1 物理模型

建立豎直圓管內(nèi)超臨界壓力低溫甲烷傳熱三維物理模型，如圖1所示。豎直圓管內(nèi)徑為8 mm，壁厚為1 mm，管長(zhǎng)為3 m，超臨界壓力低溫甲烷沿z軸正方向流動(dòng)。流體豎直向上流動(dòng)時(shí)，重力場(chǎng)方向與流動(dòng)方向相反；流體豎直向下流動(dòng)時(shí)，重力場(chǎng)方向與流動(dòng)方向相同。由于LNG的主要組分為甲烷，因此采用甲烷代替LNG進(jìn)行傳熱過(guò)程分析[4]。超臨界壓力LNG在擬臨界區(qū)物性變化劇烈，采用REFPROP軟件計(jì)算甲烷在5 MPa下的物性參數(shù)，結(jié)果見圖2。采用分段線性方法導(dǎo)入Fluent軟件，以便準(zhǔn)確地模擬沿程的傳熱過(guò)程。

圖1 豎直圓管內(nèi)超臨界壓力低溫甲烷傳熱三維物理模型

圖2 基于REFPROP軟件的甲烷在5 MPa下的物性參數(shù)

2 數(shù)學(xué)模型

2.1 控制方程

將質(zhì)量守恒方程、動(dòng)量守恒方程、能量守恒方程表示為通用控制方程，其具體形式見式（1）[34]。

式中，φ為通用變量；ρφ為廣義密度；Γφ為廣義擴(kuò)散系數(shù)；Sφ為廣義源項(xiàng)；U為速度矢量，m/s；t為時(shí)間，s。各參數(shù)的含義見表1。

表1 各參數(shù)的含義

表1中，ρ為密度，kg/m3，u、v、w分別為x、y、z方向的速度分量，m/s；μ為動(dòng)力黏度，Pa·s；p為流體微元體上的壓力，Pa；λ為導(dǎo)熱系數(shù)，W/(m·K)；ST為體積熱源項(xiàng)，W/m3；T為溫度，K；cp為比定壓熱容，kJ/（kg·K）；fx、fy、fz分別為流體合力在x、y、z方向上的分量，N。

2.2 數(shù)值方法與驗(yàn)證

本文采用RNGk-ε湍流模型。為了保證計(jì)算精度，湍流動(dòng)能方程、耗散率方程均采用二階迎風(fēng)格式離散，采用SIMPLE算法求解壓力速度耦合方程，各方程殘差設(shè)定為10－6，當(dāng)殘差小于設(shè)定值且甲烷出口溫度保持不變時(shí)，認(rèn)為計(jì)算收斂。網(wǎng)格由ICEM生成，徑向加密，其獨(dú)立性驗(yàn)證結(jié)果見表2。為了兼顧模擬的準(zhǔn)確性和計(jì)算速度，選擇網(wǎng)格3為本模型計(jì)算網(wǎng)格。

表2 網(wǎng)格獨(dú)立性驗(yàn)證

為了進(jìn)一步驗(yàn)證模型的準(zhǔn)確性，按照文獻(xiàn)[35]的實(shí)驗(yàn)工況，對(duì)直徑為2 mm的豎直管道中超臨界壓力二氧化碳的豎直向上流動(dòng)過(guò)程進(jìn)行了模擬。沿程壁面溫度的模擬數(shù)據(jù)及實(shí)驗(yàn)數(shù)據(jù)的對(duì)比如圖3所示。圖3中，橫坐標(biāo)x/d為x軸方向的距離與管徑之比。利用圖3中數(shù)據(jù)計(jì)算了熱流密度qw不同時(shí)沿程壁面溫度的最大相對(duì)誤差：當(dāng)qw分別為51 955、39 389、31 882 W/m2時(shí)，沿程最大相對(duì)誤差分別為10.28%、8.63%、11.00%。由此可知，qw＝39 389 W/m2時(shí)吻合最好，充分說(shuō)明了本數(shù)值模型和方法的可靠性。

圖3 沿程壁面溫度的模擬數(shù)據(jù)及實(shí)驗(yàn)數(shù)據(jù)的對(duì)比

3 結(jié)果與討論

3.1 溫度與對(duì)流傳熱系數(shù)

不同熱流密度下流體向上/向下流動(dòng)時(shí)的溫度及對(duì)流換熱系數(shù)如圖4所示。

圖4 不同熱流密度下流體向上/向下流動(dòng)時(shí)的溫度及對(duì)流換熱系數(shù)

由圖4可以看出，在不同的熱流密度下，當(dāng)主流溫度跨越擬臨界溫度（Tcp＝193 K）時(shí)，壁面溫度和平均溫度突然增大；熱流密度越大，主流溫度跨越擬臨界溫度越迅速，壁面溫度和平均溫度突變極值越大，突變區(qū)域越小，向上流動(dòng)和向下流動(dòng)的突變溫差越大；在相同的熱流密度下，向上流動(dòng)的壁面溫度和平均溫度均高于向下流動(dòng)的壁面溫度，而主流溫度則基本一致；在主流溫度超過(guò)擬臨界溫度后，向上流動(dòng)和向下流動(dòng)的壁面溫差消失，壁面溫度呈平緩上升趨勢(shì)，符合第二類邊界條件下的管內(nèi)強(qiáng)制對(duì)流換熱的一般規(guī)律。在高熱流密度（125 k W/m2和100 k W/m2）下，壁面溫度和平均溫度在惡化區(qū)域均出現(xiàn)了震蕩現(xiàn)象，同一熱流密度下向上流動(dòng)的震蕩程度強(qiáng)于向下流動(dòng)，熱流密度越大震蕩越劇烈，溫度的不穩(wěn)定性越強(qiáng)，熱流密度為125 kW/m2的工況下震蕩現(xiàn)象強(qiáng)于熱流密度為100 kW/m2的工況；在熱流密度為125 kW/m2時(shí)，向上流動(dòng)的主流區(qū)溫度出現(xiàn)了明顯的上升滯緩現(xiàn)象。

由圖4（d）可以看出，主流區(qū)流體接近擬臨界溫度時(shí)即為沿程對(duì)流傳熱系數(shù)下降區(qū)域，在此區(qū)域傳熱被惡化；隨著熱流密度的增大，對(duì)流傳熱系數(shù)平均值增大，惡化區(qū)間被縮短；相同熱流密度下的豎直向上和向下流動(dòng)的換熱系數(shù)沿程均有一定差異，向上流動(dòng)的換熱系數(shù)的變化先于向下流動(dòng)；隨著熱流密度的減小，差異區(qū)間被拉長(zhǎng)，這是因?yàn)闊崃髅芏仍叫?，甲烷跨越擬臨界溫度區(qū)間較長(zhǎng)；在低熱流密度（75 kW/m2和50 kW/m2）下，局部對(duì)流傳熱系數(shù)具有不穩(wěn)定性，向上流動(dòng)的不穩(wěn)定性強(qiáng)于向下流動(dòng)。

3.2 溫度場(chǎng)與流場(chǎng)

為了解釋相同熱流密度下向上流動(dòng)和向下流動(dòng)傳熱特性和不穩(wěn)定性的差異，研究了qw＝125 k W/m2時(shí)豎直向上/向下流動(dòng)的溫度場(chǎng)和流場(chǎng)，結(jié)果如圖5所示。由圖5（a）可以看出，在相同的熱流密度下，流動(dòng)方向不同時(shí)主流區(qū)溫度場(chǎng)分布基本一致，但向上流動(dòng)時(shí)的壁面溫度與向下流動(dòng)時(shí)的壁面溫度差別較大，向上流動(dòng)的熱擾動(dòng)較深，熱影響層較厚，在壁面上形成的類氣膜較厚，由此產(chǎn)生的不穩(wěn)定性較強(qiáng)，而向下流動(dòng)的熱擾動(dòng)較淺，僅停留在壁面上，熱影響層和類氣膜較薄，不穩(wěn)定性較弱。由圖5（b）可以看出，在軸線位置甲烷的主流區(qū)流速分布基本一致，而在擬臨界區(qū)間，軸向上的溫差使其密度也產(chǎn)生差異，在重力場(chǎng)的作用下向上流動(dòng)的流場(chǎng)呈“M”型分布，向下流動(dòng)的流場(chǎng)則呈“D”型分布。

圖5 q w＝125 k W/m2時(shí)豎直向上/向下流動(dòng)時(shí)的溫度場(chǎng)和流場(chǎng)云圖

qw＝125 k W/m2的工況下向上/向下流動(dòng)時(shí)沿程不同位置的徑向流場(chǎng)分布和湍動(dòng)能分布如圖6所示。由圖6（a）可以看出，在熱影響區(qū)z＝0.90 m的截面出現(xiàn)了明顯的“M”型流場(chǎng)分布，而在z＝0.75 m的截面上湍動(dòng)能大幅度減小，并且影響區(qū)域向管內(nèi)軸向方向擴(kuò)展，湍流的發(fā)展受到抑制；在z＝1.25 m的截面上“M”型流場(chǎng)開始消失，且湍動(dòng)能開始恢復(fù)，且其最大值向管壁貼近，湍流得以發(fā)展，傳熱進(jìn)入正常階段。

由圖6（b）可以看出，與向上流動(dòng)的情況不同，豎直向下流動(dòng)的流場(chǎng)轉(zhuǎn)變晚于向上流動(dòng)，在z＝1.25 m的截面開始出現(xiàn)軸線上的“D”型流場(chǎng)，使此截面上湍動(dòng)能發(fā)展受到抑制，在z＝1.50 m截面流場(chǎng)恢復(fù)，湍動(dòng)能也得以恢復(fù)。由此可見，“M”型流動(dòng)和“D”型流動(dòng)均會(huì)降低湍動(dòng)能，使湍流發(fā)展受阻，傳熱被惡化。

圖6 q w＝125 kW/m2的工況下向上/向下流動(dòng)時(shí)沿程不同位置的徑向流場(chǎng)和湍動(dòng)能分布

4 結(jié) 論

（1）主流區(qū)甲烷溫度可用于預(yù)測(cè)傳熱惡化的發(fā)生與否，當(dāng)甲烷溫度接近擬臨界溫度時(shí)，管壁溫度迅速上升，沿程傳熱系數(shù)減小，傳熱被惡化；甲烷溫度超過(guò)擬臨界溫度后，壁面溫度恢復(fù)正常值，傳熱恢復(fù)正常。

（2）傳熱的不穩(wěn)定性均出現(xiàn)在傳熱惡化區(qū)間，隨著熱流密度的減小，傳熱惡化區(qū)被拉長(zhǎng)，但壁面溫度突升的極大值減?。桓邿崃髅芏认鹿軆?nèi)壁溫度和平均溫度具有不穩(wěn)定性，低熱流密度下沿程傳熱系數(shù)具有不穩(wěn)定性。

（3）在相同的熱流密度下，向上流動(dòng)時(shí)在惡化區(qū)間流場(chǎng)呈“M”型，而向下流動(dòng)時(shí)在惡化區(qū)間流場(chǎng)呈“D”型分布；“M”型和“D”型流場(chǎng)均使湍動(dòng)能減小，湍流的發(fā)展受到抑制，是傳熱惡化的主要原因；向上流動(dòng)的熱影響層大于向下流動(dòng)，傳熱惡化向正常傳熱的轉(zhuǎn)變先于向下流動(dòng)，而向上流動(dòng)的不穩(wěn)定性均強(qiáng)于向下流動(dòng)，熱影響層內(nèi)的類氣膜是傳熱不穩(wěn)定的主要因素。