李 琳，姜文全，李婷婷，楊 帆，宿詩雨，石杰峰

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

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

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

1 物理模型

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

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

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

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

2.1 控制方程

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

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

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

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

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

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

表2 網(wǎng)格獨立性驗證

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

圖3 沿程壁面溫度的模擬數(shù)據(jù)及實驗數(shù)據(jù)的對比

3 結(jié)果與討論

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

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

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

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

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

3.2 溫度場與流場

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

圖5 q w＝125 k W/m2時豎直向上/向下流動時的溫度場和流場云圖

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

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

圖6 q w＝125 kW/m2的工況下向上/向下流動時沿程不同位置的徑向流場和湍動能分布

4 結(jié) 論

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

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

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