COSINE系統分析程序并行耦合計算研究

胡碩文，林 萌*，劉 京，陳冠宇

(1. 上海交通大學核科學與工程學院，上海 200240；2. 國家電力投資集團科學技術研究院有限公司，北京 100029；3. 中國核動力研究設計院，四川 成都 610213)

1 引言

全范圍模擬機可以用于仿真核電站運行過程，并對核電站的事故工況進行安全分析[1-3]。全范圍模擬機是培訓核電站的操縱員的關鍵設備，同時也是核電站的標準配置，準確性和實時性是其關鍵參數[4]。當模擬機中的單一仿真模型節點龐大時，會導致程序計算速度變慢影響實時性，開展并行計算研究以提高其計算速度具有重要的工程意義。

熱工水力仿真程序是開發全范圍模擬機的基礎，國外在該方面起步早，現已研發了多個知名的熱工水力計算程序。例如，輕水堆系統的瞬態分析最佳估算程序RELAP5可以模擬兩相流，用于全范圍模擬機模型開發，是核電廠分析的重要工具[5]。陳俊杰等針對換熱問題，基于RELAP5程序，使用壁面溫度作為耦合參數，研究了熱構件壁面溫度耦合并行技術，提高了計算速度[6]。但是，目前針對國產自主化的核電廠熱工水力仿真程序并行計算的研究較少。

為提高全范圍模擬機的計算速度，并促進我國核電熱工水力仿真程序的發展，本文使用具有我國自主知識產權的COSINE(COre and System INtegrated Engine for design and analysis)軟件包里的cosSyst(系統分析程序)[7]，進行了熱工模型搭建，通過并行計算，縮短了計算時間，全范圍模擬機更容易滿足實時性的要求。

2 流體熱工水力基礎方程

在全范圍模擬機中，熱工水力程序并行耦合計算通過子系統間的參數傳遞實現。典型的傳遞方法如下，在耦合時刻和物理空間耦合節點處，一方面，子系統A通過流體熱工水力基礎方程將計算的內部節點參數傳遞給子系統B用于更新子系統B的邊界參數；另一方面，子系統B在接收更新的邊界參數后，同樣調用流體熱工水力基礎方程完成求解后，也將子系統B的內部節點參數傳遞給子系統A用于更新子系統A的邊界參數；依此往復，實現熱工水力程序并行耦合計算。因此，具體耦合的參數依賴于所選用熱工水力程序的基礎方程。

本文使用的是cosSyst程序，其采用一維瞬態非熱平衡非力平衡的兩流體模型，基于半隱式差分離散格式，求解基礎變量壓力、液相焓值、氣相焓值、液相速度、氣相速度和空泡份額等。

流體在封閉空間內流動，滿足質量守恒，即

(1)

(2)

在計算過程中，流體滿足能量守恒，即

(3)

(4)

另外，流體保持動量守恒，即

+Γ?fuI+F?I，f+F?W，f+F?L，f

(5)

+Γ?guI+F?I，g+F?W，g+F?L，g

(6)

式中，下標f是液相，下標g是氣相，α是空泡份額，ρ是密度，u是速度，Γ?是由于相變的質量生成率，p是壓力，g是重力加速度，θ是控制體與豎直方向之間的夾角，F?I是相間摩擦力，uI是相界面速度，F?W是壁面摩擦力，F?L是局部阻力，h是焓，hI是相界面焓，q?I，f是相界面傳到液相的熱流，q?I，g是相界面傳到氣相的熱流，q?W，f是壁面傳到液相的熱流，q?W，g是壁面傳到氣相的熱流，t是時間，x是距離。

壓力、兩相焓值、兩相流速與空泡份額是上述方程重要的因變量，因此選擇壓力、焓值、質量流量與空泡份額作為模型耦合參數。

3 耦合方案

模擬機熱工水力建模使用的基本模塊包含管道、閥門和加熱器等模塊，由多種模塊組成復雜的熱工水力系統。cosSyst在求解模塊的熱工參數時，都是使用上述方程進行求解，求解原理是相同的。因此，本文以管道、閥門和加熱器為基礎建立簡單模型，驗證耦合方案在cosSyst的適用性。模型拆分前為整體模型，將模型進行拆分得到的兩個模型為耦合模型。整體模型如圖1所示，耦合模型如圖2所示。部件#1、#5、#6與#10為邊界，部件#2與#9模擬由5個控制體組成的管道，部件#3與#8模擬閥門，部件#4與#7模擬單控制體管道，部件#11模擬加熱器。

圖1 整體模型節點

耦合模型之間通過傳遞參數進行耦合，計算時間每過100毫秒(設置的同步周期)上游模型與下游模型的壓力P、焓值h、空泡份額α和質量流量F會通過數據庫進行數據交互，耦合方案如圖2所示。

圖2 耦合模型節點

在核電廠流體系統中，流體存在單液相、單氣相與兩相這三種典型狀態，故本文對整體模型與耦合模型在上述三種流體狀態下仿真，對比整體模型與耦合模型的參數偏差。同時，并行耦合計算的熱工水力程序需要在穩態工況下和瞬態工況下均能得到準確的計算結果，測試將覆蓋穩態和瞬態工況。針對圖1和圖2所示的熱工模型，在達到穩定計算的基礎上，調節閥門#3的開度。在50至60秒期間，閥門開度由0.5線性降至0.25；在200至210秒期間，閥門開度由0.25線性升至0.5，閥門#3開度變化如圖3所示。通過調節閥門開度改變局部阻力造成流體流速與壓力場分布變化，以制造典型的瞬態工況。

圖3 閥門#3開度隨時間變化

4 數值計算分析

為研究耦合方案的適用性，本文選取了耦合處的管道#4和管道#7的質量流量、壓力、溫度，將耦合模型與整體模型的上述參數分別在初始穩態40秒時刻、中間穩態190秒時刻、最終穩態300秒時刻及0-300秒過程內進行對比。

4.1 單液相流體耦合

單液相流體為過冷水，設置入口邊界#1壓力為7.0 MPa，液相比焓為1214 kJ/kg，設置出口邊界#10壓力為6.8 MPa。

在50秒時刻，由于閥門開度減小，質量流量減小，故管道#4和#7溫度開始升高。在200秒時刻，由于閥門開度增大，質量流量也增大，故管道#4和#7溫度開始降低，最終恢復至初始狀態。節點質量流量、壓力和溫度如圖4、圖5和圖6所示，參數偏差如表1所示，其中，F代表質量流量，P代表壓力，T代表溫度。從表1可以看出，穩態下單液相流體的整體模型和耦合模型之間的偏差不超過0.2%，整個仿真過程內最大偏差不超過2%。

圖4 節點質量流量

表1 單液相流體整體模型與耦合模型參數偏差

圖5 節點壓力

圖6 節點溫度

4.2 單氣相流體耦合

單氣相流體為過熱蒸汽，設置入口邊界#1壓力為6.0MPa，氣相比焓為2975kJ/kg，設置出口邊界#10壓力為5.6MPa。

質量流量、壓力和溫度仿真結果如圖7、圖8和圖9所示，參數偏差如表2所示。從表2可以看出，穩態下單氣相流體的整體模型和耦合模型之間的參數偏差不超過3%，整個過程偏差最大不超過4%。

圖7 節點質量流量

圖8 節點壓力

圖9 節點溫度

表2 單氣相流體整體模型與耦合模型參數偏差

4.3 兩相流體耦合

流體為飽和水與飽和蒸汽的混合物，設置入口邊界#1壓力為6.0MPa，液相體積份額為0.5，氣相體積份額為0.5，設置出口邊界#10壓力為5.6MPa。

質量流量、壓力和溫度仿真結果如圖10、11和12所示，參數偏差如表3所示。管道#4和#7壓力降低，故管道#4和#7溫度出現了下降趨勢。從表3可以看出，穩態下兩相流體的整體模型和耦合模型之間的參數偏差不超過3%，整個過程偏差最大不超過6%。

圖10 節點總質量流量

圖11 節點壓力

圖12 節點溫度

表3 兩相流體整體模型與耦合模型參數偏差

單液相、單氣相與兩相模型仿真300秒的CPU計算消耗時間如表4所示。由表4可以看出，將模型拆分為耦合模型，通過并行計算可以縮短CPU計算消耗時間，提高計算速度。

表4 CPU計算消耗時間

4.4 耦合方法適用性理論分析

本文所采用的耦合方案實質是基于熱工水力程序基礎守恒方程的一種邊界值動態顯式耦合方法，通過不斷傳遞和修改不同熱工水力系統邊界上的參數值實現。即，被計算的熱工水力系統的邊界值在本計算時間步長內保持不變，該邊界值由另一熱工水力系統在前一耦合時刻根據計算動態確定。根據原理分析，此種耦合方法對于系統內部的節點數或系統運行狀態的改變并無依賴。因此，理論上該方法既適用于小節點規模，也適用于大節點規模的熱工水力模型，既適用于穩態或瞬態計算，也適用于事故工況計算。值得注意的是，此種耦合方法由于采用顯式耦合方式，耦合頻率對耦合計算結果有影響。對于參數變化劇烈的過程(參數變化響應時間與耦合周期相當，一般為0.1s)，如果耦合頻率過低，可能會造成模型參數計算震蕩。

在上述耦合方案基礎上，針對熱工水力模型，本文以簡單對象為例進行了三種典型算例計算，包括單相液、單相氣和氣液兩相。對于兩流體的熱工水力程序而言，雖然本文采用的對象簡單，但從計算節點守恒方程的角度進行了全覆蓋，同時也測試了方程中所涉及的熱工水力參數的主要狀態。相較于全范圍模擬機更復雜的熱工水力系統模型、更多的計算節點和狀態，其計算方法和原理與本文的簡單對象模型是完全一致的。因此，理論上本文所采用的耦合方法可推廣至全范圍模擬機的熱工水力模型耦合中，此適用性在目前已開展的全范圍模擬機熱工水力建模與耦合中得到了部分證實。

5 結論

全范圍模擬機模型龐大時計算速度會下降，難以達到實時仿真。為此，將整體模型進行拆分，采用耦合方案進行模型耦合，通過并行計算仿真了模型在單液相流體、單氣相流體和兩相流體下的工況。通過整體模型與使用并行計算的耦合模型之間的節點參數對比，得到如下結論：

1)整體模型與耦合模型的參數在穩態下的最大偏差不超過3%，瞬態下最大偏差不超過6%。

2)通過并行計算，CPU消耗時間節省了18%以上，計算速度得到了提高。

3)本文所采用的耦合方案適用于全范圍模擬機的模型耦合，有利于全范圍模擬機實現實時模擬。