基于Modelica的并聯多通道系統流動不穩定性仿真研究

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

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

0 引言

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

1 相關研究

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

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

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

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

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

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

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

2.1 守恒方程模型

汽液質量守恒方程為：

汽液動量守恒方程為：

汽液能量守恒方程為：

式中：k表示汽相g或液相f；i表示界面；α表示體積百分比；μ表示流速；h表示焓值；Γk表示界面傳質；Fwk表示壁面摩擦；Fik表示界面摩擦；Fa表示外力場作用；qik表示界面傳熱；qwk表示壁面傳熱；Γk、Fwk、Fik、Fa、qik、qwk均為未知量，需要通過本構封閉方程獲得。

2.2 本構封閉方程模型

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

2.3 方程離散求解

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

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

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

3 并聯通道系統模型建立與工程驗證

3.1 并聯雙通道系統模型搭建與驗證

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

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

搭建完畢后，參考實驗臺架通過參數面板設置結構參數與初始化參數，建立系統模型，如圖3所示。

Fig.3 Model of parallel heating dual-channel system圖3 并聯加熱雙通道系統模型

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

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

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

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

3.2 不同加熱通道數對并聯多通道流動不穩定邊界的影響

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

3.2.1 流動不穩定性邊界

Ishii［16］研究發現，流動不穩定性邊界可由無量綱相變數Npch和過冷度數Nsub構成。目前許多學者采用Npch-Nsub圖來繪制流動不穩定邊界和不穩定區域的準確性，效果較好，因此本文也沿用該方法進行實驗。無量綱相變數Npch和無量綱過冷度數Nsub為：

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

3.2.2 多通道流動不穩定性邊界

根據上述實驗工況進行并聯雙、三、四通道系統進行流動不穩定性邊界的研究，仿真曲線如圖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 四通道對稱均勻加熱流量曲線

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

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

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

4 結語

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

同時，本文使用NUMAP 軟件進行并聯多通道加熱模型的仿真研究，在加熱通道數分別為2、3、4 時，計算流動不穩定發生點過冷度數Nsub和相變數Npch的值，發現其流動不穩定性邊界相差均在±5%以內，證實了在多個并聯加熱通道時，采用并聯雙通道結構可獲得其流動不穩定邊界。