童秋實(shí)，鄧康杰，呂 星，曾小康，王艷林，楊 浩，陳 路

（1.蘇州同元軟控信息技術(shù)有限公司，江蘇 蘇州 215123；2.中國(guó)核動(dòng)力研究設(shè)計(jì)院中核核反應(yīng)堆熱工水力技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213）

0 引言

流動(dòng)不穩(wěn)定性（flow instability）指系統(tǒng)受到瞬時(shí)擾動(dòng)后偏離原來(lái)穩(wěn)定狀態(tài)，熱工參數(shù)發(fā)生非周期性漂移或周期性振蕩現(xiàn)象［1］。對(duì)于存在多個(gè)并聯(lián)加熱通道的系統(tǒng)，加熱通道間可能存在并聯(lián)通道之間的流動(dòng)不穩(wěn)定性現(xiàn)象，表現(xiàn)為各加熱通道間的流量發(fā)生周期性脈動(dòng)，此時(shí)將導(dǎo)致臨界熱流密度（critical heat flux）明顯下降或?qū)υO(shè)備產(chǎn)生機(jī)械振蕩，從而危及系統(tǒng)、設(shè)備的安全運(yùn)行。因此，在兩相流研究中并聯(lián)加熱通道流動(dòng)不穩(wěn)定性現(xiàn)象一直是重要的研究課題。

1 相關(guān)研究

針對(duì)并聯(lián)加熱多通道系統(tǒng)的流動(dòng)不穩(wěn)定性現(xiàn)象，許多學(xué)者開(kāi)展了實(shí)驗(yàn)研究。Cheng 等［2］通過(guò)實(shí)驗(yàn)研究自然循環(huán)棒束形并聯(lián)通道流動(dòng)不穩(wěn)定，主要研究穩(wěn)壓器連接在加熱通道上游時(shí)，出現(xiàn)的壓力降流動(dòng)不穩(wěn)定性的實(shí)驗(yàn)現(xiàn)象和產(chǎn)生機(jī)理，并研究入口過(guò)冷度和熱流密度對(duì)流動(dòng)不穩(wěn)定性的影響。唐瑜等［3］在強(qiáng)迫循環(huán)條件下采用矩形并聯(lián)通道開(kāi)展流動(dòng)不穩(wěn)定性的實(shí)驗(yàn)研究，考慮了系統(tǒng)壓力、質(zhì)量流量、入口過(guò)冷度等熱工參數(shù)對(duì)流動(dòng)不穩(wěn)定性界限的影響，通過(guò)實(shí)驗(yàn)數(shù)據(jù)擬合得到了流動(dòng)不穩(wěn)定性無(wú)量綱關(guān)系式。Taleyarhan 等［4］建立平行加熱通道系統(tǒng)模型考慮過(guò)冷沸騰、熱通量分布和管形阻力系數(shù)等因素，描述沸水堆中的熱工水力現(xiàn)象，仿真分析表明在平行加熱通道中功率、流量會(huì)對(duì)流動(dòng)不穩(wěn)定性發(fā)生存在較大影響，實(shí)驗(yàn)數(shù)據(jù)證實(shí)了該模型對(duì)流動(dòng)不穩(wěn)定性邊界的預(yù)測(cè)能力。

目前，主要采用實(shí)驗(yàn)方式研究并聯(lián)加熱通道流動(dòng)不穩(wěn)定性現(xiàn)象［5］，但實(shí)驗(yàn)研究周期長(zhǎng)、實(shí)驗(yàn)工況有限、無(wú)法較好捕捉流動(dòng)不穩(wěn)定性內(nèi)部機(jī)理現(xiàn)象等。因此，有必要采用仿真手段進(jìn)行流動(dòng)不穩(wěn)定性現(xiàn)象模擬，通過(guò)不同工況條件下并聯(lián)加熱通道的快速仿真模擬，研究流動(dòng)不穩(wěn)定機(jī)理與確定流動(dòng)不穩(wěn)定性邊界，為實(shí)驗(yàn)研究提供參考依據(jù)。

在工程實(shí)踐中，通常采用一維系統(tǒng)仿真程序?qū)Σ⒙?lián)加熱通道的流動(dòng)不穩(wěn)定性現(xiàn)象進(jìn)行仿真分析，由于強(qiáng)迫循環(huán)下并聯(lián)加熱通道內(nèi)流體狀態(tài)變化劇烈，因此一般將流道分為若干個(gè)控制體，該方法相較于研究流道內(nèi)流場(chǎng)溫度場(chǎng)的三維數(shù)值模擬方法既降低了計(jì)算復(fù)雜性，又表達(dá)了流動(dòng)不穩(wěn)定性熱工水力特性，現(xiàn)已廣泛應(yīng)用于流動(dòng)不穩(wěn)定性的建模仿真分析中［6-7］。常用的仿真程序一般采用Fortran、C語(yǔ)言、Matlab/Simulink 等語(yǔ)言進(jìn)行開(kāi)發(fā)，或直接采用RELAP5、THEATRE 等核工程專(zhuān)業(yè)仿真軟件進(jìn)行計(jì)算［8-9］。然而，上述程序和軟件采用卡片式建模，建模過(guò)程缺乏圖形化界面，用戶難以直觀獲得系統(tǒng)設(shè)備的構(gòu)成、連接關(guān)系、參數(shù)信息等。此外，美國(guó)MMS、芬蘭APROS、中國(guó)STAR-90 等圖形化仿真軟件已廣泛應(yīng)用于動(dòng)力系統(tǒng)仿真領(lǐng)域，通過(guò)拖拽式建模及設(shè)備參數(shù)可視化配置的方式極大降低了建模工作量，但商業(yè)軟件代碼相對(duì)封閉，不便于用戶根據(jù)自身需求進(jìn)行二次開(kāi)發(fā)。

為此，蘇州同元軟控信息技術(shù)有限公司聯(lián)合中國(guó)核動(dòng)力研究設(shè)計(jì)院構(gòu)建了兩相熱工水力特性模型架構(gòu)，采用兩流體六方程數(shù)學(xué)模型開(kāi)發(fā)了基于Modelica 多專(zhuān)業(yè)、多層次、多系統(tǒng)的核反應(yīng)堆統(tǒng)一建模與分析平臺(tái)NUMAP（Nuclear reactor Unified Modeling and Analysis Platform）［10］。具體為，基于Modelica 語(yǔ)言開(kāi)發(fā)的NUMAP 可參考系統(tǒng)拓?fù)鋱D進(jìn)行拖拽式建模，相較于傳統(tǒng)熱工水力軟件具備強(qiáng)大、快捷的可視化交互功能，能極大提升建模效率［11-12］。本文參考某實(shí)驗(yàn)臺(tái)架，在NUMAP 平臺(tái)上建立高溫高壓并聯(lián)加熱雙通道系統(tǒng)模型，通過(guò)分析在均勻、非均勻加熱工況下的流動(dòng)不穩(wěn)定性、實(shí)驗(yàn)臺(tái)架和NUMAP 系統(tǒng)模型的出口質(zhì)量含氣率值，驗(yàn)證了NUMAP 平臺(tái)可用于模擬分析并聯(lián)加熱通道流動(dòng)不穩(wěn)定性瞬態(tài)特性。

2 基于Modelica的兩流體六方程理論模型

NUMAP 熱工水力模型庫(kù)基于Modelica 語(yǔ)言，采用兩流體六方程進(jìn)行描述。兩流體六方程模型是核反應(yīng)堆熱工水力系統(tǒng)分析常用的數(shù)學(xué)方程模型，能較為準(zhǔn)確地模擬核反應(yīng)堆系統(tǒng)在正常運(yùn)行工況、事故工況下汽液兩相的傳質(zhì)傳熱行為，可用于核反應(yīng)堆系統(tǒng)的設(shè)計(jì)驗(yàn)證、運(yùn)行模擬、安全分析等場(chǎng)景。

基于Modelica 語(yǔ)言兩流體六方程模型考慮了實(shí)際兩相流體流動(dòng)過(guò)程中，兩相流體具備的不同物性、流速、溫度，兩相之間的質(zhì)量、能量和動(dòng)量交換，即針對(duì)汽相和液相分別建立質(zhì)量、動(dòng)量和能量守恒方程。為了使控制方程組封閉，還增加了相間界面摩擦、相間傳熱、相間質(zhì)量交換、壁面摩擦、壁面?zhèn)鳠岬缺緲?gòu)方程。

2.1 守恒方程模型

汽液質(zhì)量守恒方程為：

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

汽液能量守恒方程為：

式中：k表示汽相g或液相f；i表示界面；α表示體積百分比；μ表示流速；h表示焓值；Γk表示界面?zhèn)髻|(zhì)；Fwk表示壁面摩擦；Fik表示界面摩擦；Fa表示外力場(chǎng)作用；qik表示界面?zhèn)鳠幔籷wk表示壁面?zhèn)鳠?；Γk、Fwk、Fik、Fa、qik、qwk均為未知量，需要通過(guò)本構(gòu)封閉方程獲得。

2.2 本構(gòu)封閉方程模型

NUMAP 采用的熱工水力守恒方程，在數(shù)學(xué)上封閉求解時(shí)需補(bǔ)充本構(gòu)封閉方程模型，包括流動(dòng)阻力系數(shù)模型、換熱系數(shù)模型、流型圖及流型轉(zhuǎn)換判定條件、水和蒸汽物性計(jì)算等模型［13］。

2.3 方程離散求解

守恒方程與本構(gòu)封閉方程構(gòu)成了求解熱工水力問(wèn)題的非線性方程組。在方程組的求解過(guò)程中，通過(guò)數(shù)值方法進(jìn)行求解模型中壓力—速度耦合關(guān)系。在離散過(guò)程中，采用Patankar、Spalding 在交叉網(wǎng)格上速度和壓力耦合求解［14］的基礎(chǔ)上提出半隱式方法（Semi-Implicit Method for Pressure Linked Equations，SIMPLE），即求解壓力耦合方程的半隱方法。

對(duì)于一維流體，交叉網(wǎng)格指將速度及壓力分別儲(chǔ)存在兩套不同的網(wǎng)格上。其中，壓力儲(chǔ)存在控制體中心（節(jié)點(diǎn)）處，速度儲(chǔ)存在兩控制體間的界面上。在計(jì)算速度修正值時(shí)，無(wú)需考慮相鄰界面上的速度修正值對(duì)該界面上速度修正值的影響，即界面上速度修正值僅由相鄰控制體節(jié)點(diǎn)間的壓力修正值所決定（見(jiàn)圖1）。

Fig.1 Storage of pressure and velocity圖1 壓力與速度存儲(chǔ)

3 并聯(lián)通道系統(tǒng)模型建立與工程驗(yàn)證

3.1 并聯(lián)雙通道系統(tǒng)模型搭建與驗(yàn)證

并聯(lián)加熱雙通道實(shí)驗(yàn)在高溫高壓汽水兩相熱工水力實(shí)驗(yàn)裝置上開(kāi)展，裝置設(shè)計(jì)壓力17.2 MPa，溫度350 ℃，實(shí)驗(yàn)質(zhì)量流量范圍100～5 000 kg/（m2·s）［15］。參考實(shí)驗(yàn)臺(tái)架設(shè)備，利用NUMAP 平臺(tái)熱工水力模型庫(kù)中的管道、分支及進(jìn)出口邊界等模型，采用拖拽式建模方式快速建立并聯(lián)加熱雙通道系統(tǒng)模型，如圖2所示。

Fig.2 Thermal hydraulic model library of NUMAP圖2 NUMAP熱工水力模型庫(kù)

搭建完畢后，參考實(shí)驗(yàn)臺(tái)架通過(guò)參數(shù)面板設(shè)置結(jié)構(gòu)參數(shù)與初始化參數(shù)，建立系統(tǒng)模型，如圖3所示。

Fig.3 Model of parallel heating dual-channel system圖3 并聯(lián)加熱雙通道系統(tǒng)模型

仿真計(jì)算工況的開(kāi)展方式是在系統(tǒng)壓力及入口參數(shù)保持不變的前提下，階梯式緩慢提升加熱功率，每提升一次功率后保持2 min，直至流動(dòng)不穩(wěn)定現(xiàn)象發(fā)生時(shí)加熱功率不變，并記錄通道出口的質(zhì)量含氣率，仿真計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果如圖4、圖5所示。

Fig.4 Calculation error of uniform heating圖4 均勻加熱計(jì)算誤差

Fig.5 Calculation error of non-uniform heating圖5 非均勻加熱計(jì)算誤差

參考實(shí)驗(yàn)工況，對(duì)均勻加熱和非均勻并聯(lián)加熱雙通道流動(dòng)不穩(wěn)定性分別進(jìn)行仿真模擬，選取流動(dòng)不穩(wěn)定性發(fā)生點(diǎn)通道出口質(zhì)量含氣率的仿真計(jì)算結(jié)果與實(shí)驗(yàn)值進(jìn)行比較，均勻加熱和非均勻加熱工況下仿真計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果誤差最大值分別為18.87%、15.7%，滿足該實(shí)驗(yàn)工況下20%的誤差要求，表明NUMAP 軟件可用模擬分析并聯(lián)通道流動(dòng)不穩(wěn)定性瞬態(tài)特性。

3.2 不同加熱通道數(shù)對(duì)并聯(lián)多通道流動(dòng)不穩(wěn)定邊界的影響

由于并聯(lián)多通道系統(tǒng)存在多種核電設(shè)備中，其流動(dòng)不穩(wěn)定性會(huì)導(dǎo)致臨界熱流密度（CHF）明顯下降或設(shè)備發(fā)生機(jī)械振蕩，影響核電站正常運(yùn)行。因此，在核電站的運(yùn)行過(guò)程中，如何避免并聯(lián)多通道系統(tǒng)流動(dòng)不穩(wěn)定現(xiàn)象發(fā)生一直是熱議的研究課題。本文在驗(yàn)證NUMAP 軟件對(duì)并聯(lián)加熱雙通道流動(dòng)不穩(wěn)定性瞬態(tài)計(jì)算能力的基礎(chǔ)上，對(duì)并聯(lián)多通道流動(dòng)不穩(wěn)定性邊界進(jìn)行仿真計(jì)算。

3.2.1 流動(dòng)不穩(wěn)定性邊界

Ishii［16］研究發(fā)現(xiàn)，流動(dòng)不穩(wěn)定性邊界可由無(wú)量綱相變數(shù)Npch和過(guò)冷度數(shù)Nsub構(gòu)成。目前許多學(xué)者采用Npch-Nsub圖來(lái)繪制流動(dòng)不穩(wěn)定邊界和不穩(wěn)定區(qū)域的準(zhǔn)確性，效果較好，因此本文也沿用該方法進(jìn)行實(shí)驗(yàn)。無(wú)量綱相變數(shù)Npch和無(wú)量綱過(guò)冷度數(shù)Nsub為：

式中：Q為總加熱功率，單位為W；W為通道進(jìn)口總流量，單位為kg/s；hfg為汽化潛熱，單位為J/kg；Δhsub為進(jìn)口欠焓，單位為J/kg；vfg為飽和氣液兩相比容差，單位為m3/kg；vf為飽和液相比容，單位為m3/kg。

3.2.2 多通道流動(dòng)不穩(wěn)定性邊界

根據(jù)上述實(shí)驗(yàn)工況進(jìn)行并聯(lián)雙、三、四通道系統(tǒng)進(jìn)行流動(dòng)不穩(wěn)定性邊界的研究，仿真曲線如圖6—圖8所示。

Fig.6 Flow rate curve for parallel dual-channel system with symmetric uniform heating圖6 雙通道對(duì)稱(chēng)均勻加熱流量曲線

Fig.7 Flow rate curve for parallel triple-channel system with symmetric uniform heating圖7 三通道對(duì)稱(chēng)均勻加熱流量曲線

Fig.8 Flow rate curve for parallel quadruple-channel system with symmetric uniform heating圖8 四通道對(duì)稱(chēng)均勻加熱流量曲線

在并聯(lián)加熱多通道系統(tǒng)仿真中，隨著加熱功率上升，加熱通道出口達(dá)到飽和后開(kāi)始出現(xiàn)氣液兩相流動(dòng)。隨著氣液兩相流動(dòng)的含氣率逐漸增大，當(dāng)通道出口含氣率增大到一定程度后，并聯(lián)多通道系統(tǒng)中各通道的流量出現(xiàn)周期性流量脈動(dòng)，即發(fā)生了流動(dòng)不穩(wěn)定性，如圖6—圖8 中局部放大部分所示。

在不同工況下設(shè)置加熱通道數(shù)目分別為2、3、4 時(shí)，計(jì)算得到流動(dòng)不穩(wěn)定發(fā)生點(diǎn)過(guò)冷度數(shù)Nsub和相變數(shù)Npch，并在此基礎(chǔ)上進(jìn)行并聯(lián)多通道系統(tǒng)流動(dòng)不穩(wěn)定性邊界的比較，如圖9 所示。由此可見(jiàn)，流動(dòng)不穩(wěn)定性邊界相差在±5%以內(nèi)，仿真計(jì)算的結(jié)果與王艷林等［15］實(shí)驗(yàn)結(jié)果一致，即存在多個(gè)并聯(lián)加熱通道時(shí)，采用并聯(lián)雙通道結(jié)構(gòu)即可獲得其流動(dòng)不穩(wěn)定性邊界。

Fig.9 Flow instability boundaries under different heating channel numbers圖9 不同加熱通道數(shù)目下流動(dòng)不穩(wěn)定性邊界

4 結(jié)語(yǔ)

流動(dòng)不穩(wěn)定性是核反應(yīng)堆熱工安全的重要研究課題之一，仿真分析可作為開(kāi)展流動(dòng)不穩(wěn)定性研究的重要手段。本文基于Modelica 語(yǔ)言開(kāi)發(fā)的NUMAP 軟件搭建了并聯(lián)雙通道加熱模型，仿真模擬結(jié)果表明，實(shí)驗(yàn)誤差滿足實(shí)驗(yàn)工況下的誤差要求，驗(yàn)證了模型對(duì)流動(dòng)不穩(wěn)定現(xiàn)象的瞬態(tài)計(jì)算分析能力。

同時(shí)，本文使用NUMAP 軟件進(jìn)行并聯(lián)多通道加熱模型的仿真研究，在加熱通道數(shù)分別為2、3、4 時(shí)，計(jì)算流動(dòng)不穩(wěn)定發(fā)生點(diǎn)過(guò)冷度數(shù)Nsub和相變數(shù)Npch的值，發(fā)現(xiàn)其流動(dòng)不穩(wěn)定性邊界相差均在±5%以內(nèi)，證實(shí)了在多個(gè)并聯(lián)加熱通道時(shí)，采用并聯(lián)雙通道結(jié)構(gòu)可獲得其流動(dòng)不穩(wěn)定邊界。