沈艷榮，楊 冬，陸啟亮，姜 勇，王少飛

（西安交通大學(xué)動力工程多相流國家重點實驗室，陜西西安 710049）

非均勻加熱蒸發(fā)管內(nèi)兩相流動不穩(wěn)定性數(shù)值分析

沈艷榮，楊 冬，陸啟亮，姜 勇，王少飛

（西安交通大學(xué)動力工程多相流國家重點實驗室，陜西西安 710049）

針對蒸發(fā)管非均勻加熱工況，采用包含水平、傾斜上升、垂直上升和垂直下降流動的優(yōu)化內(nèi)螺紋管，建立了基于均相流模型的蒸發(fā)管數(shù)學(xué)模型，采用時域法模擬了蒸發(fā)管內(nèi)兩相流動不穩(wěn)定性。計算了不穩(wěn)定發(fā)生時，進口流量在各時間點的數(shù)值，并與Siemens公司的計算結(jié)果進行了對比，結(jié)果符合良好，表明本文采用的數(shù)學(xué)模型和數(shù)值方法在模擬兩相流動不穩(wěn)定性上具有一定的可靠性。分析了系統(tǒng)參數(shù)對兩相流動不穩(wěn)定性的影響，并與均勻加熱工況進行了比較。結(jié)果表明：非均勻加熱工況系統(tǒng)參數(shù)影響規(guī)律與均勻加熱工況具有相似性，增加進口壓力或進口流量系統(tǒng)的穩(wěn)定性提高；減小進口阻力系數(shù)或增大出口阻力系數(shù)系統(tǒng)的穩(wěn)定性降低。

非均勻加熱；兩相流動不穩(wěn)定性；時域法；數(shù)值分析

兩相流動不穩(wěn)定性廣泛存在于核電、火電、石油化工等許多領(lǐng)域的系統(tǒng)和設(shè)備中，它不僅會引起部件的強制機械振動，還會引起管壁溫度脈動，使管材疲勞失效，甚至導(dǎo)致傳熱惡化，嚴(yán)重威脅系統(tǒng)的安全運行［1］。因此，兩相流動不穩(wěn)定性的研究對核反應(yīng)堆蒸汽發(fā)生器和大型直流鍋爐等工業(yè)設(shè)備的設(shè)計及安全運行有著十分重要的意義。

在過去的幾十年中，氣液兩相流動不穩(wěn)定性的研究大多采用熱負(fù)荷均勻分布的假設(shè)［2］，而實際上熱負(fù)荷分布通常是不均勻的，且研究對象大多為流動形式簡單的光管。因此，本文針對蒸發(fā)管非均勻加熱工況，充分考慮蒸發(fā)管幾何結(jié)構(gòu)的多樣性，采用包含水平、傾斜上升、垂直上升、垂直下降多種流動形式的優(yōu)化內(nèi)螺紋管，建立基于均相流模型的蒸發(fā)管數(shù)學(xué)模型，采用時域法對非均勻加熱蒸發(fā)管內(nèi)兩相流動不穩(wěn)定性進行數(shù)值模擬，并與Siemens公司的計算結(jié)果進行對比，以及對非均勻加熱工況和均勻加熱工況下兩相流動不穩(wěn)定性的影響因素進行比較。

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

通常采用的兩相流動不穩(wěn)定性數(shù)值分析方法有頻域法和時域法［3－4］兩種。頻域法［5－6］是通過將控制方程無量綱化、線性化及Laplace變換得到系統(tǒng)的特征方程，只能獲得脈動發(fā)生時的閾值。時域法［7－8］是通過對控制方程進行離散，采用數(shù)值分析方法獲得脈動發(fā)生前、發(fā)生時和發(fā)生后各參數(shù)隨時間的變化，重現(xiàn)兩相流動不穩(wěn)定現(xiàn)象，有助于探索其發(fā)生機理。本文采用時域法。

1.1 基本假設(shè)

為簡化模型，作了如下假設(shè)：1）采用一維近似，考慮壓縮性和熱膨脹性；2）氣液兩相間處于熱力平衡，即不考慮欠熱沸騰和相間熱力弛豫，兩相區(qū)采用均相流模型來描述；3）介質(zhì)與金屬壁只在徑向方向進行換熱，而不考慮軸向換熱；4）在同一截面內(nèi)工質(zhì)溫度和速度分布是均勻的，且管內(nèi)介質(zhì)只沿軸向流動，無內(nèi)部環(huán)流；5）在能量方程中忽略黏性耗散、動能和勢能的影響。

1.2 控制方程

對于均相流模型，管中的兩相流動和單相流動均用以下守恒方程描述［9］。

質(zhì)量守恒方程：

動量守恒方程：

能量守恒方程：

式（2）中未計入形阻。狀態(tài)方程：ρ＝f（p，h）（4）式中：ρ、h、M、p分別為流體的密度、比焓、質(zhì)量流量和壓力；ql為線密度；A為管的橫截面積；θ為流動方向與水平面的夾角；dn為蒸發(fā)管內(nèi)徑；f為摩擦阻力系數(shù)，兩相時需進行修正［10］。

當(dāng)為單相時，上述方程的物性參數(shù)對應(yīng)該相的物性參數(shù)；當(dāng)為兩相時，采用折和參數(shù)。

2 數(shù)值計算方法

2.1 離散方法

采用內(nèi)節(jié)點法進行網(wǎng)格劃分。控制體劃分示于圖1。采用控制容積法對方程進行離散。界面上的物性參數(shù)采用一階迎風(fēng)差分［11］，即物性參數(shù)取來流方向上最后1個節(jié)點值，在時間上采用全隱式。

根據(jù)圖1的控制體劃分方式對控制體離散如下。

圖1 控制體示意圖Fig.1 Schematic diagram of control volume

對于節(jié)點i，質(zhì)量守恒方程為：

動量守恒方程為：

能量守恒方程為：

狀態(tài)方程為：

2.2 邊界條件

邊界條件為：1）進口焓為常數(shù)；2）進口壓力為常數(shù)；3）進、出口壓降為常數(shù)。

2.3 數(shù)值計算步驟

2.4 系統(tǒng)結(jié)構(gòu)及分析方法

蒸發(fā)管結(jié)構(gòu)示于圖2，包括水平流動段AB、FG，垂直上升段CD、EF，垂直下降段GH和傾斜上升段BC、DE。蒸發(fā)管從進口到出口總長度為48.44 m，根據(jù)熱負(fù)荷分布將蒸發(fā)管分為16段，每段長度如圖2所示。每段內(nèi)的線密度相等，其中AB、FG、GH為不受熱段，其余受熱段的線密度如圖3所示。

圖2 蒸發(fā)管結(jié)構(gòu)示意圖Fig.2 Schematic diagram of evaporative tube

圖3 線密度沿蒸發(fā)管高度分布Fig.3 Distribution of linear density along height of evaporative tube

本文在流體穩(wěn)定流動時對熱負(fù)荷施加一定擾動（增加5%），計算進、出口流量隨時間的變化情況，從而判斷系統(tǒng)的穩(wěn)定性。當(dāng)進、出口流量衰減振蕩時，系統(tǒng)穩(wěn)定；當(dāng)進、出口流量發(fā)散振蕩時，系統(tǒng)不穩(wěn)定；當(dāng)進、出口流量等幅振蕩時，系統(tǒng)臨界穩(wěn)定。

2.5 網(wǎng)格無關(guān)性分析

為確定數(shù)值計算結(jié)果的準(zhǔn)確性，需進行網(wǎng)格無關(guān)性分析。取進口壓力為12.25 MPa，進口焓為1 256.6 kJ／kg，穩(wěn)態(tài)出口含氣率為1.199，初始穩(wěn)態(tài)進口流量為0.14 kg／s。將16段分別均勻劃分為若干個管段，管段的長度即為空間網(wǎng)格的長度。

首先進行時間網(wǎng)格無關(guān)性分析，為了保證時間步長與空間步長的獨立性，取空間步長為0.02 m。時間步長從1.6 s逐漸減小到0.8 s，進口流量先增大后減小，從1.5 s減小到0.8 s時，進口流量變化不大，說明取時間步長為1.5 s對計算結(jié)果無影響，如圖4a所示。與時間網(wǎng)格無關(guān)性分析類似，取時間步長為1.5 s，空間步長從0.3 m逐漸減小到0.02 m，進口流量逐漸減小，從0.05 m逐漸減小到0.02 m時，進口流量變化不大，說明取空間步長為0.05 m對計算結(jié)果無影響，如圖4b所示。

圖4 時間步長和空間步長對進口流量的影響Fig.4 Effect of temporal and spatial grid sizes on inlet flow rate

3 計算結(jié)果及驗證

圖5示出蒸發(fā)管采用φ31.8 mm×6.21 mm的優(yōu)化內(nèi)螺紋管、進口壓力為12.25 MPa、進口焓為1 256.6 kJ／kg、穩(wěn)態(tài)出口含氣率為1.199、初始穩(wěn)態(tài)進口流量為0.14 kg／s條件下不穩(wěn)定發(fā)生時進口流量脈動曲線。作為對本文計算模型可靠性的驗證，與相同計算參數(shù)下Siemens公司的計算結(jié)果進行對比，由圖5可看出，振幅相對誤差約為8%，周期相差約7～10 s，表明本文采用的數(shù)學(xué)模型和數(shù)值方法在模擬兩相流動不穩(wěn)定性上具有一定的可靠性。

圖5 進口流量脈動曲線Fig.5 Oscillation curve for inlet flow rate

4 系統(tǒng)參數(shù)對穩(wěn)定性的影響

圖6 進口壓力對系統(tǒng)穩(wěn)定性的影響Fig.6 Effect of inlet pressure on system stability

圖6示出其他參數(shù)不變時，進口壓力對系統(tǒng)穩(wěn)定性的影響。由圖6可見，隨著壓力的增大，流量脈動的振幅逐漸減小，回復(fù)到穩(wěn)定狀態(tài)的時間也逐漸減少，系統(tǒng)的穩(wěn)定性提高，與文獻(xiàn)［12］中均勻加熱工況的結(jié)論相似。

圖7a示出其他參數(shù)不變時，進口流量對系統(tǒng)穩(wěn)定性的影響。隨進口流量的增大，流量脈動的振幅逐漸減小，恢復(fù)到穩(wěn)定狀態(tài)的時間也逐漸減少，系統(tǒng)的穩(wěn)定性提高，與文獻(xiàn)［11］中均勻加熱工況的結(jié)論相似。圖7b示出其他參數(shù)不變時，進口阻力系數(shù)kin、出口阻力系數(shù)kout對系統(tǒng)穩(wěn)定性的影響。圖7b表明，減小kin或增大kout會降低系統(tǒng)的穩(wěn)定性，與文獻(xiàn)［13］中均勻加熱工況的結(jié)論相似。

圖7 進口流量和進出口阻力系數(shù)對穩(wěn)定性的影響Fig.7 Effect of inlet flow rate and inlet／outlet pressure drop coefficients on system stability

5 結(jié)論

本文針對蒸發(fā)管非均勻加熱工況，充分考慮蒸發(fā)管幾何結(jié)構(gòu)的多樣性，采用包含水平、傾斜上升、垂直上升、垂直下降多種流動形式的優(yōu)化內(nèi)螺紋管，建立了基于均相流模型的蒸發(fā)管數(shù)學(xué)模型，采用隱式迎風(fēng)差分法和迭代法對非均勻加熱蒸發(fā)管內(nèi)兩相流動不穩(wěn)定進行了模擬，結(jié)論如下。

1）兩相流動不穩(wěn)定發(fā)生時，采用本文模型計算的進口流量脈動曲線與Siemens公司的計算結(jié)果符合良好，表明本文采用的數(shù)學(xué)模型和數(shù)值方法在模擬兩相流動不穩(wěn)定性上具有一定的可靠性。

2）非均勻加熱工況系統(tǒng)參數(shù)影響規(guī)律與均勻加熱工況具有相似性，增加蒸發(fā)管進口壓力或進口流量系統(tǒng)的穩(wěn)定性提高，減小進口阻力系數(shù)或增大出口阻力系數(shù)系統(tǒng)的穩(wěn)定性降低。

［1］ ZHANG T J，WEN J T，PELES Y，et al.Twophase refrigerant flow instability analysis and active control in transient electronics cooling systems［J］.International Journal of Multiphase Flow，2011，37（1）：84-97.

［2］ 周源，閆曉，王艷林.加熱雙通道密度波流動不穩(wěn)定性數(shù)值研究［J］.原子能科學(xué)技術(shù)，2013，47（4）：552-556.

ZHOU Yuan，YAN Xiao，WANG Yanlin.Numerical investigation on density wave oscillation of two heated parallel channels［J］.Atomic Energy Science and Technology，2013，47（4）：552-556（in Chinese）.

［3］ CHATOORGOON V.Stability of supercritical fluid flow in a single-channel natural-convection loop［J］.International Journal of Heat and Mass Transfer，2001，44（10）：1 963-1 972.

［4］ ZHOU J，PODOWSKI M Z.Modeling and analysis of hydrodynamic instabilities in two-phase flow using two-fluid model［J］.Nuclear Engineering and Design，2001，204（1）：129-142.

［5］ GóMEZ T O，CLASS A，LAHEY R T，et al.Stability analysis of a uniformly heated channel with supercritical water［J］.Nuclear Engineering and Design，2008，238（8）：1 930-1 939.

［6］ 鄢炳火，于雷.改進的并聯(lián)管內(nèi)兩相流密度波脈動線性均相模型［J］.原子能科學(xué)技術(shù)，2008，42（增刊）：146-149.

YAN Binghuo，YU Lei.Advanced homogeneous flow model for density wave oscillation of twophase flow in parallel tubes［J］.Atomic EnergyScience and Technology，2008，42（Suppl.）：146-149（in Chinese）.

［7］ 田文喜，田曉艷，馮健，等.超臨界水冷堆CSR1000流動不穩(wěn)定性研究［J］.核動力工程，2013，34（1）：45-51.

TIAN Wenxi，TIAN Xiaoyan，F(xiàn)ENG Jian，et al.Flow instability analysis of supercritical water-cooled reactor CSR1000 based on frequency domain［J］.Nuclear Power Engineering，2013，34（1）：45-51（in Chinese）.

［8］ 姚偉，匡波，楊燕華，等.沸騰兩相自然循環(huán)系統(tǒng)動態(tài)不穩(wěn)定性的數(shù)值分析［J］.原子能科學(xué)技術(shù)，2002，36（1）：62-66.

YAO Wei，KUANG Bo，YANG Yanhua，et al.Numerical analysis of dynamic instability of boiling two-phase natural circulation system［J］.Atomic Energy Science and Technology，2002，36（1）：62-66（in Chinese）.

［9］ XIONG Ting，YAN Xiao，HUANG Hanfang，et al.Modeling and analysis of supercritical flow instability in parallel channels［J］.International Journal of Heat and Mass Transfer，2013，57（2）：549-557.

［10］車得福.鍋爐［M］.西安：西安交通大學(xué)出版社，2008：406-407.

［11］陶文銓.數(shù)值傳熱學(xué)［M］.西安：西安交通大學(xué)出版社，2011：140-141.

［12］黃可欣，高峰，畢勤成，等.超臨界鍋爐并聯(lián)傾斜內(nèi)螺紋管兩相流密度波型脈動數(shù)值研究［J］.核動力工程，2008，29（4）：63-68，73.

HUANG Kexin，GAO Feng，BI Qincheng，et al.Numerical research on density wave oscillation of two-phase flow in parallel inclined internally ribbed pipes for supercritical pressure boiler［J］.Nuclear Power Engineering，2008，29（4）：63-68，73（in Chinese）.

［13］熊挺，閆曉，俞冀陽，等.并行雙通道內(nèi)超臨界水流動不穩(wěn)定性數(shù)值分析［J］.核動力工程，2012，33（2）：62-65.

XIONG Ting，YAN Xiao，YU Jiyang，et al.Numerical analysis of supercritical flow instability in parallel dual channel［J］.Nuclear Power Engineering，2012，33（2）：62-65（in Chinese）.

Numerical Analysis of Two-phase Flow Instability in Evaporating Tube under Non-uniform Heating Condition

SHEN Yan-rong，YANG Dong，LU Qi-liang，JIANG Yong，WANG Shao-fei
（State Key Laboratory of Multiphase Flow in Power Engineering，Xi’an Jiaotong University，Xi’an 710049，China）

Based on the homogeneous flow model，the mathematic model of optimization of inside rifled tube under non-uniform heating condition which included horizontal flow，inclined upward flow，vertical upward flow and vertical downward flow was established.The two-phase flow instability was simulated by the time-domain approach.With the same parameters as the Siemens，the variation of inlet flow rate with time was obtained.It is shown that the results have a good agreement with those of the Siemens.Therefore，the mathematic model and numerical method were reliable to simulate the two-phase flow instability.The parametric effects on the two-phase flow instability under non-uniform heating condition were analyzed and compared with the effects under uniform heating condition.The results show that they are similar to those underuniform heating condition，which indicates that increasing inlet pressure or inlet flow rate will improve the system stability，and decreasing the inlet pressure drop coefficient or increasing the outlet pressure drop coefficient will reduce the system stability.

non-uniform heating；two-phase flow instability；time-domain approach；numerical analysis

TL333

A

1000-6931（2015）09-1573-06

10.7538／yzk.2015.49.09.1573

2014-05-15；

2014-07-17

中國科學(xué)院戰(zhàn)略性先導(dǎo)科技專項資助（XDA07030100）

沈艷榮（1988—），女（滿族），河北保定人，碩士研究生，熱能工程專業(yè)