基于Modelica的并聯(lián)多通道系統(tǒng)流動不穩(wěn)定性仿真研究

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

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

0 引言

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

1 相關(guān)研究

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

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

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

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

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

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

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

2.1 守恒方程模型

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

汽液動量守恒方程為：

汽液能量守恒方程為：

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

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

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

2.3 方程離散求解

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)語

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

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