巢 飛，楊 文，邰 云，邱金榮，單建強(qiáng)

(1.武漢第二船舶設(shè)計(jì)研究所，湖北 武漢 430064；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

0 引 言

自然循環(huán)是指依靠密度差和高度差形成的驅(qū)動(dòng)力驅(qū)動(dòng)流體循環(huán)流動(dòng)的現(xiàn)象。自然循環(huán)在核動(dòng)力系統(tǒng)中非常常見(jiàn)，其對(duì)反應(yīng)堆安全有重要影響，因此準(zhǔn)確模擬自然循環(huán)現(xiàn)象是當(dāng)前反應(yīng)堆熱工水力分析中最重要研究?jī)?nèi)容之一。

目前主流熱工水力系統(tǒng)分析程序如RELAP5，CATHARE，TRACE 等普遍采用低階差分法來(lái)求解基本守恒方程，其在模擬自然循環(huán)不穩(wěn)定性問(wèn)題時(shí)存在預(yù)測(cè)誤差。為提高自然循環(huán)不穩(wěn)定性問(wèn)題的計(jì)算精度，本文采用基于兩流體雙壓力兩相流模型的高精度數(shù)值解法對(duì)Welander 自然循環(huán)進(jìn)行計(jì)算分析。

1 兩流體雙壓力模型

兩流體雙壓力模型基本守恒方程：

體積份額輸運(yùn)方程

式中：αg表示汽相體積份額；Pg和Pf分別代表汽相和液相壓力，Pa；vint為相間界面速度，m·s?1；μ為壓力松弛因子，Pa?1·s?1，表示汽液兩相壓力達(dá)到平衡的松弛速率；A為管道截面積，m2；ρint為相間界面密度，k g·m?3；Γg為單位體積液-汽相間質(zhì)量交換率，kg·m?3·s?1。汽、液兩相體積份額關(guān)系應(yīng)滿足：

其中 τP為壓力松弛時(shí)間，s。

質(zhì)量守恒方程

式中：Γf=?Γg為單位體積汽-液相間質(zhì)量交換率，kg·m?3·s?1；ρ為相密度，kg·m?3；v為相速度，m·s?1；f為液相標(biāo)識(shí)；g為汽相標(biāo)識(shí)。

動(dòng)量守恒方程：

式中：Pint為相間界面壓力，Pa；fwall,k為阻力系數(shù);ρc為連續(xù)相密度，kg·m?3；CiD為相間阻力系數(shù)；Aint為單位體積相間界面面積，m?1；gx為重力加速度在流動(dòng)方向的投影，m·s?2。

能量守恒方程

式中：ek為k相比內(nèi)能/J·kg?1；Qwg為汽相單位體積壁面換熱量，W·m?3；Qwf為液相單位體積壁面換熱量，W·m?3；Hig為單位體積相間界面到汽相的換熱系數(shù)，W·m?3·K?1；Hif為單位體積相間界面到液相換熱系數(shù)，W·m?3·K?1；Ts為相間界面溫度，Tk為k相溫度，K；為相間界面汽相比焓，J·kg?1；為相間界面液相比焓，J·kg?1。

2 高精度數(shù)值算法

2.1 半隱式數(shù)值解法

采用參考文獻(xiàn)[1]作者提出的半隱數(shù)值算法求解兩流體雙壓力模型。為達(dá)到快速計(jì)算的效果，半隱算法只對(duì)時(shí)間常數(shù)小、傳播速度快的量隱式處理，如兩相速度、壓力梯度、相間界面質(zhì)量和動(dòng)量交換項(xiàng)，其他量顯式處理。為了解決兩相和單相之間轉(zhuǎn)變帶來(lái)的數(shù)值不穩(wěn)定性問(wèn)題，動(dòng)量守恒方程和質(zhì)量守恒方程處理成和與差的微分方程，所有守恒方程瞬態(tài)項(xiàng)采用展開(kāi)的形式，再進(jìn)行數(shù)值差分。采用交錯(cuò)網(wǎng)格劃分控制體，如圖1 所示。K，L等表示求解質(zhì)量、能量守恒方程的控制體，存儲(chǔ)流場(chǎng)中的標(biāo)量如壓力、比內(nèi)能、空泡份額等；j代表動(dòng)量控制體（虛線包圍），用于求解動(dòng)量守恒方程，存儲(chǔ)流場(chǎng)中速度矢量。

圖1 交錯(cuò)網(wǎng)格示意圖Fig.1 Schematic drawing of staggered mesh

2.2 瞬態(tài)項(xiàng)二階差分格式

目前現(xiàn)有的系統(tǒng)程序普遍采用一階時(shí)間差分（BDF1）離散守恒方程瞬態(tài)項(xiàng)。BDF1 離散誤差大，給守恒方程的差分方程引入了“人工”擴(kuò)散，宜采用精度更高的二階時(shí)間差分（BDF2）來(lái)降低瞬態(tài)項(xiàng)的離散誤差。二階時(shí)間差分根據(jù)Taylor 級(jí)數(shù)展開(kāi)法得到[2]：

其中：?φ/?t表示守恒方程中的瞬態(tài)項(xiàng)，Δt=tn+1?tn，Δtold=tn?tn?1。該式考慮了時(shí)間步長(zhǎng)的變化，具有一般性，在系統(tǒng)程序中應(yīng)用更為方便。BDF1 和BDF2 可處理成統(tǒng)一形式：

2.3 對(duì)流項(xiàng)三階TVD 差分格式

目前主流熱工水力程序采用一階施主元法離散對(duì)流項(xiàng)。一階施主元法本質(zhì)上是一階迎風(fēng)，計(jì)算精度低，因此宜采用穩(wěn)定高階施主元法，以提高對(duì)流項(xiàng)的計(jì)算精度。雙壓力模型守恒方程中與速度相關(guān)的項(xiàng)均采用施主元法計(jì)算。以能量方程對(duì)流項(xiàng)為例，其離散形式為：

由于速度本身定義在接管上，因此對(duì)流項(xiàng)的計(jì)算精度很大程度上取決于接管上的標(biāo)量（如）的計(jì)算方案。接管上標(biāo)量數(shù)值計(jì)算的統(tǒng)一形式可以寫成：

其中：Δx是控制體長(zhǎng)度；?φ/?x是控制體邊界上的梯度；ψ(r)為通量限制函數(shù)；r為梯度比，由下式給出：

空間離散方案統(tǒng)一形式（13）由兩部分組成，第一部分為穩(wěn)定的一階迎風(fēng)方案，第二部分為反擴(kuò)散項(xiàng)。因此，物理量 φ˙計(jì)算精度由反擴(kuò)散項(xiàng)中的通量限制函數(shù)ψ(r)決定。

穩(wěn)定三階TVD 差分格式（命名為TOU_TVD_CFL）為[2]：

式 中C為Courant 數(shù)。

3 計(jì)算分析

3.1 Welander 自然循環(huán)問(wèn)題

Welander 提出了豎直回路中單相自然循環(huán)問(wèn)題來(lái)研究流體的不穩(wěn)定性行為[3]。Welander 自然循環(huán)回路由2 個(gè)平行的豎直絕熱管，頂部冷源和底部熱源組成，如圖2 所示，回路總長(zhǎng)為L(zhǎng)，截面積為A，熱阱和熱源長(zhǎng)度為 Δs。通過(guò)頂部熱阱的冷卻和底部熱源加熱，流體在回路頂部和底部產(chǎn)生密度差，從而驅(qū)動(dòng)流體發(fā)生自然循環(huán)流動(dòng)。通過(guò)選擇合適的浮力與摩擦阻力的比值，流動(dòng)可能會(huì)出現(xiàn)穩(wěn)定、弱不穩(wěn)定和強(qiáng)不穩(wěn)定行為。在文獻(xiàn)[3]中，Welander 推導(dǎo)了實(shí)驗(yàn)回路層流下流動(dòng)不穩(wěn)定性理論邊界，Ambrosini 和Ferreri[4]將不穩(wěn)定性邊界擴(kuò)展湍流情況，如圖3 所示。Welander 自然循環(huán)穩(wěn)定性邊界繪制成無(wú)量綱重力系數(shù) α和無(wú)量綱摩擦系數(shù) ε的關(guān)系式圖，無(wú)量綱量具體表達(dá)式為：

圖2 Welander 自然循環(huán)問(wèn)題示意圖Fig.2 Schematic diagram of Welander natural circulation

圖3 Welander 自然循環(huán)問(wèn)題流動(dòng)不穩(wěn)定性邊界Fig.3 Flow instability map of Welander natural circulation

式中：Welander 定義的層流摩擦系數(shù)為R=64/Ref·Welander 定義的換熱系數(shù)為溫差ΔT=當(dāng)a=0.316 4/4,b=0.25時(shí)，表示湍流條件；當(dāng)a=16,b=1時(shí)，表示層流條件。

對(duì)于Welander 單相振蕩自然循環(huán)問(wèn)題，選取Ambrosinia 和Ferreri[5]設(shè)置的幾何參數(shù)和邊界條件。管內(nèi)徑為0.1 m，熱源和冷源長(zhǎng)度為0.1 m。管道充滿壓力為0.1 MPa 的過(guò)冷水。熱源溫度為30℃，與流體換熱系數(shù)為20 000 W·m?2·K?1，冷源溫度維持20℃，與流體換熱系數(shù)為20 000 W·m?2·K?1。本文研究Welander 單相振蕩自然循環(huán)問(wèn)題3 個(gè)不同的實(shí)驗(yàn)工況：強(qiáng)不穩(wěn)定模式、弱不穩(wěn)定模式和穩(wěn)定模式，見(jiàn)表1。圖4 給出了Welander 自然循環(huán)回路數(shù)值模擬節(jié)點(diǎn)圖。在進(jìn)行模擬計(jì)算時(shí)，時(shí)間步長(zhǎng)為0.5 s，熱源和冷源分別采用一個(gè)節(jié)點(diǎn)模擬，兩垂直絕熱管都均分為N個(gè)控制體，因此環(huán)路控制體總數(shù)為2N+2。

表1 Welander 自然循環(huán)問(wèn)題實(shí)驗(yàn)工況Tab.1 Test conditions of Welander natural circulation

圖4 Welander 自然循環(huán)節(jié)點(diǎn)圖Fig.4 Nodalization of Welander natural circulation

3.2 數(shù)值結(jié)果

對(duì)于強(qiáng)不穩(wěn)定模式工況，環(huán)路長(zhǎng)度取L=20.2 m，此時(shí) (α,ε)為 (343.79,2.346 84)。強(qiáng)不穩(wěn)定模式時(shí)，實(shí)驗(yàn)回路流量出現(xiàn)流動(dòng)反轉(zhuǎn)現(xiàn)象，并作平均速度為零的等幅振蕩。為便于觀測(cè)到流動(dòng)不穩(wěn)定性行為，實(shí)驗(yàn)之初，給定回路一個(gè)極小流量的擾動(dòng)。

圖5 給出了使用低階差分方案BDF1+FOU 在不同控制體數(shù)下的數(shù)值結(jié)果。由圖可知控制體數(shù)低于62 時(shí)，數(shù)值誤差過(guò)大，無(wú)法預(yù)測(cè)出不穩(wěn)定行為；控制體數(shù)達(dá)到102 時(shí)才預(yù)測(cè)出不穩(wěn)定性行為，控制體數(shù)越多，數(shù)值誤差越小，預(yù)測(cè)的不穩(wěn)定行為振蕩周期越短。由于數(shù)值誤差的阻尼作用，粗糙網(wǎng)格可以穩(wěn)定物理上不穩(wěn)定的系統(tǒng)。采用不同計(jì)算節(jié)點(diǎn)，低階差分?jǐn)?shù)值結(jié)果會(huì)出現(xiàn)系統(tǒng)穩(wěn)定和不穩(wěn)定2 種相矛盾的情況。因此數(shù)值誤差對(duì)不穩(wěn)定性行為準(zhǔn)確預(yù)測(cè)有很不利的影響，低階差分方案BDF1+FOU 精度低，不利于研究流動(dòng)不穩(wěn)定行為。

圖5 強(qiáng)不穩(wěn)定模式質(zhì)量流量：低階BDF1+FOU方案的計(jì)算結(jié)果Fig.5 Mass flow of strong instability mode:calculation results of BDF1+FOU scheme

圖6 給出了高階差分方案與低階差分方案結(jié)果的對(duì)比。高階差分方案模擬流動(dòng)不穩(wěn)定性行為具有相當(dāng)大的優(yōu)勢(shì)。高階差分方案TOU_TVD_CFL 只需22 個(gè)控制體就可以預(yù)測(cè)出不穩(wěn)定行為，因此高階算法有效減小了數(shù)值擴(kuò)散，提高了計(jì)算精度。高階差分算法采用22 個(gè)控制體的數(shù)值精度可達(dá)到低階差分算法控制體數(shù)為182 的計(jì)算精度，振蕩周期相近。

圖6 強(qiáng)不穩(wěn)定模式質(zhì)量流量：高階和低階方案的計(jì)算結(jié)果對(duì)比Fig.6 Mass flow of strong instability mode:comparison between low order and high order difference scheme

圖7 給出了低階差分方案BDF1+FOU 計(jì)算的流體溫度與質(zhì)量流量關(guān)系圖，隨著控制體數(shù)增多，精度越高，混沌行為（Chaotic Behavior）越明顯。圖8 給出了高階BDF2+TOU_TVD_CFL 方案計(jì)算的流體溫度與質(zhì)量流量關(guān)系圖。同樣可以看出，隨著控制體數(shù)增多，可以預(yù)測(cè)到流動(dòng)反轉(zhuǎn)的混沌演化過(guò)程，流體溫度與質(zhì)量流量關(guān)系圖類似蝴蝶形狀。

圖7 低階方案BDF1+FOU 計(jì)算的流體溫度與質(zhì)量流量關(guān)系圖Fig.7 Fluid temperature versus mass flow calculated by low order scheme BDF1+FOU

圖8 高階BDF2+TOU_TVD_CFL 方案計(jì)算的流體溫度與質(zhì)量流量關(guān)系圖Fig.8 Fluid temperature versus mass flow calculated by high order scheme BDF2+TOU_TVD_CFL

針對(duì)弱不穩(wěn)定工況，回路總長(zhǎng)為80.2 m。弱不穩(wěn)定模式時(shí)，回路流量不會(huì)出現(xiàn)流動(dòng)反轉(zhuǎn)現(xiàn)象，其維持一個(gè)流動(dòng)方向上的近似等幅振蕩。圖9 給出了采用低階方案計(jì)算的質(zhì)量流量隨時(shí)間變化圖。控制體數(shù)為722 時(shí)質(zhì)量流量的數(shù)值結(jié)果依舊出現(xiàn)緩慢的衰減，低階差分結(jié)果誤差大。圖10 給出了弱不穩(wěn)定工況下采用高階差分方案計(jì)算的質(zhì)量流量隨時(shí)間變化圖。高階算法的結(jié)果有效減小人工粘性，抑制振蕩衰減，控制體數(shù)為102 時(shí)就可以模擬出等幅振蕩的流動(dòng)不穩(wěn)定性行為。

圖9 弱不穩(wěn)定模式質(zhì)量流量低階方案的計(jì)算結(jié)果Fig.9 Calculation results of low order scheme by mass flow of weak instability mode

圖10 弱不穩(wěn)定模式質(zhì)量流量高階方案的計(jì)算結(jié)果Fig.10 Calculation results of high order scheme by mass flow of weak instability mode

針對(duì)L=5.2 m 的穩(wěn)定模式工況，表現(xiàn)為：質(zhì)量流量初始增加到最大值，然后經(jīng)過(guò)一段時(shí)間的衰減型振蕩后維持穩(wěn)定運(yùn)動(dòng)。圖11 和圖12 給出了穩(wěn)定模式工況下質(zhì)量流量的計(jì)算結(jié)果，高階方案與低階方案計(jì)算結(jié)果對(duì)比見(jiàn)圖13。因?yàn)槿斯ふ承缘挠绊?，低階差分預(yù)測(cè)更小的振蕩振幅，更快達(dá)到穩(wěn)定狀態(tài)；控制體數(shù)為102 的低階結(jié)果精度比控制體數(shù)為22 的高階結(jié)果精度還低，振蕩衰減更快。由于高階差分對(duì)穩(wěn)態(tài)沒(méi)影響，最終高階差分預(yù)測(cè)的穩(wěn)態(tài)流量與低階結(jié)果一致。

圖11 穩(wěn)定模式質(zhì)量流量低階方案的計(jì)算結(jié)果Fig.11 Calculation results of low order scheme by mass flow of stable mode

圖12 穩(wěn)定模式質(zhì)量流量高階方案的計(jì)算結(jié)果Fig.12 Calculation results of high order scheme by mass flow of stable mode

圖13 穩(wěn)定模式質(zhì)量流量高階和低階方案的計(jì)算結(jié)果對(duì)比Fig.13 Comparison between low order and high order scheme by mass flow of stable mode

4 結(jié) 語(yǔ)

本文采用兩流體雙壓力兩相流模型高精度數(shù)值算法對(duì)Welander 自然循環(huán)的流動(dòng)不穩(wěn)定性問(wèn)題進(jìn)行計(jì)算分析，數(shù)值結(jié)果表明：高階差分格式可以以較少的網(wǎng)格精確地預(yù)測(cè)不穩(wěn)定性邊界和混沌行為，而采用低階差分在網(wǎng)格不夠精細(xì)的情況下有可能將自然循環(huán)不穩(wěn)定模式預(yù)測(cè)為穩(wěn)定模式，得到錯(cuò)誤的結(jié)果。因此高精度數(shù)值解法能有效減小數(shù)值擴(kuò)散，提高自然循環(huán)問(wèn)題的預(yù)測(cè)精度。