COSINE系統(tǒng)分析程序并行耦合計(jì)算研究

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

(1. 上海交通大學(xué)核科學(xué)與工程學(xué)院，上海 200240；2. 國(guó)家電力投資集團(tuán)科學(xué)技術(shù)研究院有限公司，北京 100029；3. 中國(guó)核動(dòng)力研究設(shè)計(jì)院，四川 成都 610213)

1 引言

全范圍模擬機(jī)可以用于仿真核電站運(yùn)行過(guò)程，并對(duì)核電站的事故工況進(jìn)行安全分析[1-3]。全范圍模擬機(jī)是培訓(xùn)核電站的操縱員的關(guān)鍵設(shè)備，同時(shí)也是核電站的標(biāo)準(zhǔn)配置，準(zhǔn)確性和實(shí)時(shí)性是其關(guān)鍵參數(shù)[4]。當(dāng)模擬機(jī)中的單一仿真模型節(jié)點(diǎn)龐大時(shí)，會(huì)導(dǎo)致程序計(jì)算速度變慢影響實(shí)時(shí)性，開(kāi)展并行計(jì)算研究以提高其計(jì)算速度具有重要的工程意義。

熱工水力仿真程序是開(kāi)發(fā)全范圍模擬機(jī)的基礎(chǔ)，國(guó)外在該方面起步早，現(xiàn)已研發(fā)了多個(gè)知名的熱工水力計(jì)算程序。例如，輕水堆系統(tǒng)的瞬態(tài)分析最佳估算程序RELAP5可以模擬兩相流，用于全范圍模擬機(jī)模型開(kāi)發(fā)，是核電廠分析的重要工具[5]。陳俊杰等針對(duì)換熱問(wèn)題，基于RELAP5程序，使用壁面溫度作為耦合參數(shù)，研究了熱構(gòu)件壁面溫度耦合并行技術(shù)，提高了計(jì)算速度[6]。但是，目前針對(duì)國(guó)產(chǎn)自主化的核電廠熱工水力仿真程序并行計(jì)算的研究較少。

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

2 流體熱工水力基礎(chǔ)方程

在全范圍模擬機(jī)中，熱工水力程序并行耦合計(jì)算通過(guò)子系統(tǒng)間的參數(shù)傳遞實(shí)現(xiàn)。典型的傳遞方法如下，在耦合時(shí)刻和物理空間耦合節(jié)點(diǎn)處，一方面，子系統(tǒng)A通過(guò)流體熱工水力基礎(chǔ)方程將計(jì)算的內(nèi)部節(jié)點(diǎn)參數(shù)傳遞給子系統(tǒng)B用于更新子系統(tǒng)B的邊界參數(shù)；另一方面，子系統(tǒng)B在接收更新的邊界參數(shù)后，同樣調(diào)用流體熱工水力基礎(chǔ)方程完成求解后，也將子系統(tǒng)B的內(nèi)部節(jié)點(diǎn)參數(shù)傳遞給子系統(tǒng)A用于更新子系統(tǒng)A的邊界參數(shù)；依此往復(fù)，實(shí)現(xiàn)熱工水力程序并行耦合計(jì)算。因此，具體耦合的參數(shù)依賴于所選用熱工水力程序的基礎(chǔ)方程。

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

流體在封閉空間內(nèi)流動(dòng)，滿足質(zhì)量守恒，即

(1)

(2)

在計(jì)算過(guò)程中，流體滿足能量守恒，即

(3)

(4)

另外，流體保持動(dòng)量守恒，即

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

(5)

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

(6)

式中，下標(biāo)f是液相，下標(biāo)g是氣相，α是空泡份額，ρ是密度，u是速度，Γ?是由于相變的質(zhì)量生成率，p是壓力，g是重力加速度，θ是控制體與豎直方向之間的夾角，F(xiàn)?I是相間摩擦力，uI是相界面速度，F(xiàn)?W是壁面摩擦力，F(xiàn)?L是局部阻力，h是焓，hI是相界面焓，q?I，f是相界面?zhèn)鞯揭合嗟臒崃?，q?I，g是相界面?zhèn)鞯綒庀嗟臒崃?，q?W，f是壁面?zhèn)鞯揭合嗟臒崃?，q?W，g是壁面?zhèn)鞯綒庀嗟臒崃鳎瑃是時(shí)間，x是距離。

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

3 耦合方案

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

圖1 整體模型節(jié)點(diǎn)

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

圖2 耦合模型節(jié)點(diǎn)

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

圖3 閥門#3開(kāi)度隨時(shí)間變化

4 數(shù)值計(jì)算分析

為研究耦合方案的適用性，本文選取了耦合處的管道#4和管道#7的質(zhì)量流量、壓力、溫度，將耦合模型與整體模型的上述參數(shù)分別在初始穩(wěn)態(tài)40秒時(shí)刻、中間穩(wěn)態(tài)190秒時(shí)刻、最終穩(wěn)態(tài)300秒時(shí)刻及0-300秒過(guò)程內(nèi)進(jìn)行對(duì)比。

4.1 單液相流體耦合

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

在50秒時(shí)刻，由于閥門開(kāi)度減小，質(zhì)量流量減小，故管道#4和#7溫度開(kāi)始升高。在200秒時(shí)刻，由于閥門開(kāi)度增大，質(zhì)量流量也增大，故管道#4和#7溫度開(kāi)始降低，最終恢復(fù)至初始狀態(tài)。節(jié)點(diǎn)質(zhì)量流量、壓力和溫度如圖4、圖5和圖6所示，參數(shù)偏差如表1所示，其中，F(xiàn)代表質(zhì)量流量，P代表壓力，T代表溫度。從表1可以看出，穩(wěn)態(tài)下單液相流體的整體模型和耦合模型之間的偏差不超過(guò)0.2%，整個(gè)仿真過(guò)程內(nèi)最大偏差不超過(guò)2%。

圖4 節(jié)點(diǎn)質(zhì)量流量

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

圖5 節(jié)點(diǎn)壓力

圖6 節(jié)點(diǎn)溫度

4.2 單氣相流體耦合

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

質(zhì)量流量、壓力和溫度仿真結(jié)果如圖7、圖8和圖9所示，參數(shù)偏差如表2所示。從表2可以看出，穩(wěn)態(tài)下單氣相流體的整體模型和耦合模型之間的參數(shù)偏差不超過(guò)3%，整個(gè)過(guò)程偏差最大不超過(guò)4%。

圖7 節(jié)點(diǎn)質(zhì)量流量

圖8 節(jié)點(diǎn)壓力

圖9 節(jié)點(diǎn)溫度

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

4.3 兩相流體耦合

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

質(zhì)量流量、壓力和溫度仿真結(jié)果如圖10、11和12所示，參數(shù)偏差如表3所示。管道#4和#7壓力降低，故管道#4和#7溫度出現(xiàn)了下降趨勢(shì)。從表3可以看出，穩(wěn)態(tài)下兩相流體的整體模型和耦合模型之間的參數(shù)偏差不超過(guò)3%，整個(gè)過(guò)程偏差最大不超過(guò)6%。

圖10 節(jié)點(diǎn)總質(zhì)量流量

圖11 節(jié)點(diǎn)壓力

圖12 節(jié)點(diǎn)溫度

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

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

表4 CPU計(jì)算消耗時(shí)間

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

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

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

5 結(jié)論

全范圍模擬機(jī)模型龐大時(shí)計(jì)算速度會(huì)下降，難以達(dá)到實(shí)時(shí)仿真。為此，將整體模型進(jìn)行拆分，采用耦合方案進(jìn)行模型耦合，通過(guò)并行計(jì)算仿真了模型在單液相流體、單氣相流體和兩相流體下的工況。通過(guò)整體模型與使用并行計(jì)算的耦合模型之間的節(jié)點(diǎn)參數(shù)對(duì)比，得到如下結(jié)論：

1)整體模型與耦合模型的參數(shù)在穩(wěn)態(tài)下的最大偏差不超過(guò)3%，瞬態(tài)下最大偏差不超過(guò)6%。

2)通過(guò)并行計(jì)算，CPU消耗時(shí)間節(jié)省了18%以上，計(jì)算速度得到了提高。

3)本文所采用的耦合方案適用于全范圍模擬機(jī)的模型耦合，有利于全范圍模擬機(jī)實(shí)現(xiàn)實(shí)時(shí)模擬。