田文喜，王明軍，曾春杰，秋穗正，蘇光輝

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

核動(dòng)力系統(tǒng)蒸汽發(fā)生器(SG)作為一、二回路能量轉(zhuǎn)換中心，對(duì)核動(dòng)力系統(tǒng)的安全性和經(jīng)濟(jì)性有著重要影響。國(guó)內(nèi)外針對(duì)SG三維熱工水力數(shù)值模擬開(kāi)展了廣泛研究，開(kāi)發(fā)了系列SG三維熱工水力分析程序，如ATHOS3[1]、TWOPORFLOW[2-3]、PORTHOS[4]、ATHLET[5-6]、GENEPI[7-8]、THERMIT[9-11]、THEDA[12-14]、THIRST[15]、CALIPSOS[16-17]、CUPID-SG[18-20]、FIT-Ⅲ[21]、STAF[22-24]和STAF-CT[25]等，部分已用于核動(dòng)力系統(tǒng)SG設(shè)計(jì)階段的熱工水力計(jì)算及校核。西安交通大學(xué)核反應(yīng)堆熱工水力團(tuán)隊(duì)(NuTHeL)對(duì)標(biāo)EPRI開(kāi)發(fā)的ATHOS程序，并針對(duì)其存在算法和模型不完善等問(wèn)題，開(kāi)發(fā)了更加完善的自主化SG三維熱工水力分析工具STAF1.0程序，通過(guò)與國(guó)際標(biāo)準(zhǔn)題和ATHOS進(jìn)行對(duì)比，證實(shí)了STAF1.0程序的準(zhǔn)確性和可靠性。在STAF1.0程序開(kāi)發(fā)基礎(chǔ)上，為解決一次側(cè)管束區(qū)橫截面的參數(shù)空間分布特性對(duì)二次側(cè)換熱影響及多種瞬態(tài)條件下的蒸汽發(fā)生器性能評(píng)估，團(tuán)隊(duì)開(kāi)發(fā)了STAF-CT程序，將一次側(cè)基于多孔介質(zhì)模型進(jìn)行三維處理，并通過(guò)開(kāi)發(fā)的一、二次側(cè)能量傳遞耦合模塊實(shí)現(xiàn)了兩側(cè)全三維耦合換熱。結(jié)果表明，一次側(cè)三維流場(chǎng)非均勻性對(duì)二次側(cè)的熱工參數(shù)分布具有重要影響。另一方面，無(wú)論是STAF1.0還是STAF-CT程序，兩者均基于ANSYS FLUENT商用CFD軟件平臺(tái)二次開(kāi)發(fā)，為解決所有代碼自主化問(wèn)題，團(tuán)隊(duì)基于開(kāi)源計(jì)算流體力學(xué)軟件OpenFOAM進(jìn)一步開(kāi)發(fā)了STAF3.0版本，在實(shí)現(xiàn)STAF-CT全部計(jì)算功能的前提下，開(kāi)發(fā)實(shí)現(xiàn)了大規(guī)模并行計(jì)算功能，大幅提升了原串行程序的計(jì)算效率，同時(shí)避免了國(guó)外商用軟件平臺(tái)二次開(kāi)發(fā)的限制問(wèn)題。

本文重點(diǎn)介紹西安交通大學(xué)核反應(yīng)堆熱工水力團(tuán)隊(duì)十余年來(lái)從STAF1.0到STAF-CT再到STAF3.0的程序開(kāi)發(fā)歷程和核心數(shù)學(xué)物理模型，并對(duì)STAF系列各種計(jì)算功能模塊及典型應(yīng)用情況進(jìn)行簡(jiǎn)要概述。

1 STAF1.0程序開(kāi)發(fā)

STAF1.0程序是基于多孔介質(zhì)方法和四方程漂移流模型開(kāi)發(fā)的初始STAF版本，其對(duì)標(biāo)EPRI開(kāi)發(fā)的ATHOS程序，團(tuán)隊(duì)旨在開(kāi)發(fā)出一款國(guó)產(chǎn)化SG三維熱工水力分析程序替代ATHOS程序。

1.1 數(shù)學(xué)物理模型

典型核動(dòng)力系統(tǒng)SG內(nèi)有數(shù)千上萬(wàn)根U型傳熱管，開(kāi)展SG全尺寸精細(xì)化三維CFD數(shù)值模擬受到計(jì)算資源和能力的限制目前還無(wú)法實(shí)現(xiàn)?？紤]到SG內(nèi)傳熱管束布置具有規(guī)律性，因此采用各向異性多孔介質(zhì)模型對(duì)管束區(qū)進(jìn)行簡(jiǎn)化，可大幅減少網(wǎng)格數(shù)量。

1.1.1二次側(cè)三維兩相流場(chǎng) SG二次側(cè)流場(chǎng)采用基于多孔介質(zhì)的四方程漂移流模型進(jìn)行描述，控制方程如下。

質(zhì)量守恒方程：

(1)

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

(2)

能量守恒方程：

(3)

空泡輸運(yùn)方程：

(4)

式中：β為孔隙率；ρm為混合物密度，kg·m-3；vm為混合物速度，m·s-1；p為壓力，Pa；μm，eff為混合物有效黏度，Pa·s；g為重力加速度，m·s-2；αg為汽相體積份額；ρg為汽相密度，kg·m-3；αl為液相體積份額；ρl為液相密度，kg·m-3；vgm為汽相相對(duì)于混合物質(zhì)量加權(quán)平均速度的相對(duì)速度，m·s-1；Sv為動(dòng)量源項(xiàng)，N·m-3；Hg為汽相焓，J·kg-1；Hl為液相焓，J·kg-1；vg為汽相速度，m·s-1；vl為液相速度，m·s-1；km為混合物導(dǎo)熱系數(shù)，W·m-1·K-1；SE為能量源項(xiàng)，W·m-3；Sg為汽相質(zhì)量源項(xiàng)，kg·m-3·s-1。

1.1.2一次側(cè)一維單相流場(chǎng) 一次側(cè)簡(jiǎn)化為一維單相流動(dòng)，假設(shè)一次側(cè)每根管束內(nèi)流量均勻分配，則管內(nèi)流速可由流量、密度和管束幾何參數(shù)計(jì)算得出。僅需求解如下一次側(cè)能量方程：

Ac,i(vp,i,inρi,incp,i,inTp,i,in-

vp,i,outρi,outcp,i,outTp,i,out)=qA,iAs,i

(5)

式中：Ac,i為傳熱管管內(nèi)流通橫截面積，m2；vp，i，in為一次側(cè)流體單元進(jìn)口速度，m·s-1；ρi，in為一次側(cè)流體單元進(jìn)口密度，kg·m-3；cp，i，in為一次側(cè)流體單元進(jìn)口比定壓熱容，J·kg-1；Tp，i，in為一次側(cè)流體單元進(jìn)口溫度，K；vp，i，out為一次側(cè)流體單元出口速度，m·s-1；ρi，out為一次側(cè)流體單元出口密度，kg·m-3；cp，i，out為一次側(cè)流體單元出口比定壓熱容，J·kg-1；Tp，i，out為一次側(cè)流體單元出口溫度，K；qA，i為單位面積傳熱管一、二次側(cè)交換的熱流密度，W·m-2；As，i為一次側(cè)流體單元對(duì)應(yīng)傳熱管外壁面面積，m2。

1.1.3一、二次側(cè)耦合換熱模型 SG內(nèi)一、二次側(cè)采用一維-三維耦合換熱方法，單位面積平均熱流密度為：

(6)

式中：qA為基于傳熱管外表面的熱流密度，W·m-2；Tp為二次側(cè)計(jì)算單元對(duì)應(yīng)的一次側(cè)質(zhì)量流量加權(quán)平均流體溫度，K；Ts為二次側(cè)計(jì)算單元溫度，K；hp為一次側(cè)平均換熱系數(shù)，W·m-2·K-1；RM為管壁熱阻，m2·K·W-1；Rf為污垢熱阻，m2·K·W-1；hs為二次側(cè)換熱系數(shù)，W·m-2·K-1。

一次側(cè)管內(nèi)流體為單相強(qiáng)制對(duì)流，采用DB公式[26]計(jì)算對(duì)流換熱系數(shù)hp：

(7)

式中：Re為以傳熱管內(nèi)徑為特征尺寸的一次側(cè)流體雷諾數(shù)；Pr為一次側(cè)流體普朗特?cái)?shù)；kp為一次側(cè)流體導(dǎo)熱系數(shù)，W·m-1·K-1。

二次側(cè)管束區(qū)存在液態(tài)給水加熱產(chǎn)生蒸汽的兩相過(guò)程，即包含單相對(duì)流換熱、過(guò)冷沸騰換熱和飽和沸騰換熱。SG二次側(cè)管束的流動(dòng)存在橫掠、縱掠和斜掠3種工況，換熱系數(shù)計(jì)算需要分別考慮。

橫掠單相對(duì)流換熱系數(shù)hs，c為：

(8)

縱掠單相對(duì)流換熱系數(shù)hs，a為：

(9)

單相對(duì)流斜掠工況采用Singhal等[1]提出的簡(jiǎn)化方法：分別求解斜掠流動(dòng)的橫流、順流速度分量對(duì)應(yīng)的換熱系數(shù)，然后取兩者平方根作為斜掠換熱系數(shù)，即：

(10)

飽和沸騰換熱系數(shù)hb需要迭代計(jì)算，其與飽和沸騰熱流密度qb的關(guān)系[27]為：

hb=Cb(qb)2/3

(11)

式中，Cb為與局部壓力相關(guān)的經(jīng)驗(yàn)常數(shù)：

(12)

式中，ps為兩側(cè)局部壓力,Pa。

飽和沸騰熱流密度qb為：

qb=hb(Tf,o-Tsat)

(13)

式中：Tf,o為污垢外表面溫度，K；Tsat為二次側(cè)流體飽和溫度，K。對(duì)式(11)和式(13)迭代計(jì)算可得飽和沸騰換熱系數(shù)。

由于管束區(qū)采用了多孔介質(zhì)模型，無(wú)法從網(wǎng)格層面獲取近壁面流體溫度和主流溫度，因此對(duì)過(guò)冷沸騰采用如下方式判斷：當(dāng)程序計(jì)算壁面溫度高于飽和溫度而主流溫度低于飽和溫度時(shí)，采用單相對(duì)流換熱系數(shù)計(jì)算式和飽和沸騰系數(shù)計(jì)算式分別計(jì)算熱流密度，當(dāng)采用飽和沸騰計(jì)算得到的熱流密度大于單相對(duì)流的熱流密度時(shí)，認(rèn)為過(guò)冷沸騰發(fā)生。

1.1.4阻力模型 SG二次側(cè)流體對(duì)管束沖刷存在橫掠、縱掠和斜掠3種工況，流經(jīng)管束阻力系數(shù)需要分別考慮，對(duì)于豎直管段單相流動(dòng)，橫掠管束的阻力系數(shù)fu[28]為：

fu=0.864Re-0.205

(14)

縱掠管束的阻力系數(shù)fz[29]為：

fz=0.046Re-0.2

(15)

在U型管束彎管區(qū)部分，采用局部坐標(biāo)系變換，在彎管區(qū)流動(dòng)可沿管束方向分解為橫掠、縱掠管束工況。對(duì)于二次側(cè)兩相流動(dòng)阻力，采用兩相摩擦阻力倍增因子方法[29]計(jì)算。

1.1.5湍流模型 基于多孔介質(zhì)方法中湍流模型適用性角度考慮，采用ATHOS程序[30-31]和GENEPI程序[32]所采用的零方程湍流模型計(jì)算兩相流體有效黏度μeff[33]：

μeff=Cρmwmlch

(16)

式中：C為經(jīng)驗(yàn)常數(shù)，C=0.047；lch為特征尺寸，定義如下：

(17)

式中：dwrap為圍板內(nèi)徑，m；dsep為汽水分離器筒體內(nèi)徑，m；δdc為下降段通道厚度，m。

1.1.6漂移流模型 STAF程序采用四方程漂移流模型，引入空泡輸運(yùn)方程克服ATHOS程序的不足。模型忽略x、y方向漂移速度，z方向汽相相對(duì)于混合相質(zhì)量加權(quán)平均速度的漂移速度wgm[30,34]為：

(18)

式中：wgj為汽相相對(duì)于混合相體積加權(quán)平均速度的漂移速度在z方向的分量，m·s-1；co為分布參數(shù)；wm為混合相質(zhì)量加權(quán)速度在z方向的分量，m·s-1；wgj和co定義[35]如下：

(19)

(20)

其中：

κ1=κ0+(1-κ0)(ρg/ρl)0.25

(21)

κ0=min(0.8,(1+e-Re/105)-1)

(22)

(23)

(24)

式中：σ為表面張力，N·m-1；pcr為臨界壓力，Pa。

1.2 STAF1.0程序驗(yàn)證

基于多個(gè)實(shí)驗(yàn)裝置的實(shí)測(cè)數(shù)據(jù)對(duì)STAF1.0程序進(jìn)行了驗(yàn)證，包括利用FRIGG棒束實(shí)驗(yàn)[36-37]和MB-2[38]蒸汽發(fā)生器實(shí)驗(yàn)數(shù)據(jù)對(duì)STAF1.0程序開(kāi)展對(duì)比驗(yàn)證，具體程序驗(yàn)證結(jié)果可參考Cong等[23-24]的文獻(xiàn)。結(jié)果表明：STAF1.0計(jì)算的FRIGG棒束空泡份額分布與實(shí)驗(yàn)值符合較好，且對(duì)于棒束間兩相流場(chǎng)預(yù)測(cè)能力優(yōu)于ATHOS程序；STAF1.0程序?qū)B-2試驗(yàn)段的溫度分布和壓降損失計(jì)算同樣與實(shí)驗(yàn)值符合良好。兩套不同實(shí)驗(yàn)數(shù)據(jù)的支撐，進(jìn)一步驗(yàn)證了STAF1.0程序結(jié)果的可靠性與準(zhǔn)確性。

1.3 基于STAF1.0程序的核電自然循環(huán)SG數(shù)值模型建立

SG一次側(cè)主要包括水室下封頭、管板和U型管束；二次側(cè)主要包括下降段、U型管束區(qū)、管束支承板、提升段以及汽水分離器裝置，如圖1所示。STAF1.0程序計(jì)算域包括下降段、管束區(qū)、管束支承板、提升段以及汽水分離器部分筒體，即圖1中紅框區(qū)域。網(wǎng)格劃分如圖2所示，其中管束支承板部分為面網(wǎng)格，采用多孔階躍模型模擬支承板壓降；二次側(cè)管束區(qū)采用多孔介質(zhì)粗網(wǎng)格劃分方式；一次側(cè)采用一維模擬，在UDF中求解并與二次側(cè)進(jìn)行耦合換熱。

圖1 AP1000典型SG結(jié)構(gòu)[39]Fig.1 Structure of SG for AP1000[39]

2 STAF-CT程序開(kāi)發(fā)

STAF1.0程序?qū)G兩側(cè)采用一維-三維耦合處理，模型過(guò)于簡(jiǎn)化，無(wú)法滿足SG一、二次側(cè)流場(chǎng)高精度熱工安全分析要求。因此團(tuán)隊(duì)在STAF1.0基礎(chǔ)上開(kāi)發(fā)了具有瞬態(tài)功能且全三維耦合的蒸汽發(fā)生器熱工水力分析程序STAF-CT，其采用多孔介質(zhì)方法簡(jiǎn)化一次側(cè)的管束區(qū)，計(jì)算過(guò)程中一、二次側(cè)管束區(qū)網(wǎng)格分別獨(dú)立計(jì)算流動(dòng)過(guò)程，利用開(kāi)發(fā)的兩側(cè)三維數(shù)據(jù)耦合交換模塊，實(shí)現(xiàn)一次側(cè)和二次側(cè)全三維熱耦合。全三維耦合考慮了一次側(cè)三維流場(chǎng)非均勻性對(duì)二次側(cè)的熱工參數(shù)分布的影響，提高了程序計(jì)算結(jié)果的可靠性。

2.1 數(shù)學(xué)物理模型升級(jí)部分

2.1.1一次側(cè)三維流場(chǎng) 蒸汽發(fā)生器一次側(cè)冷卻劑為高流速的湍流流動(dòng)，在其流動(dòng)過(guò)程中將熱量經(jīng)由傳熱管壁傳遞到二次側(cè)。為考慮STAF程序在瞬態(tài)工況上的應(yīng)用，考慮了時(shí)間項(xiàng)(二次側(cè)同理)，一次側(cè)單相三維流場(chǎng)控制方程如下。

圖2 STAF1.0程序的SG網(wǎng)格模型Fig.2 Mesh model of SG for STAF1.0 code

質(zhì)量守恒方程：

(25)

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

(26)

能量守恒方程：

(27)

本文中所涉及變量的數(shù)據(jù)主要來(lái)源于各?。ㄊ?、區(qū)）發(fā)布的歷年統(tǒng)計(jì)年鑒，部分根據(jù) 《中國(guó)工業(yè)統(tǒng)計(jì)年鑒》和 《中國(guó)能源統(tǒng)計(jì)年鑒》中相關(guān)數(shù)據(jù)直接獲得或計(jì)算得到，其中2017年的數(shù)據(jù)來(lái)源于各省 （市、區(qū)）發(fā)布的2017年國(guó)民經(jīng)濟(jì)和社會(huì)發(fā)展統(tǒng)計(jì)公報(bào)。為了消除各年價(jià)格因素的影響，本文以2000年為基期，利用國(guó)內(nèi)生產(chǎn)總值指數(shù)對(duì)名義GDP進(jìn)行平減處理。

2.1.2一次側(cè)阻力模型 一次側(cè)管內(nèi)流動(dòng)阻力模型采用多孔介質(zhì)模型中分布式阻力實(shí)現(xiàn)，即式(26)中的Si，其表達(dá)式為：

(28)

通過(guò)對(duì)SG中不同曲率半徑U型管開(kāi)展精細(xì)化CFD模擬，獲得流速與壓降的關(guān)系曲線，擬合曲線可得到式(28)中阻力系數(shù)(黏性項(xiàng)和慣性項(xiàng))，實(shí)現(xiàn)對(duì)SG一次側(cè)管束內(nèi)摩擦壓降模擬(具體數(shù)據(jù)參考Zhao等[40]的文獻(xiàn))。同時(shí)考慮多孔介質(zhì)模型各向異性，需要管束軸向布置阻力系數(shù)，徑向則為無(wú)限大阻力以考慮管壁對(duì)流動(dòng)束縛作用。

2.1.3湍流模型 STAF-CT考慮了SG一次側(cè)下封頭、管板至管束流體域，其中下封頭采用精細(xì)化網(wǎng)格建模。結(jié)果發(fā)現(xiàn)冷卻劑在下封頭入口處存在較強(qiáng)射流沖擊，蒸汽發(fā)生器入口和出口腔室存在劇烈的湍流效應(yīng)。因此，STAF-CT程序采用了Transition SST轉(zhuǎn)捩湍流模型。

2.1.4一、二次側(cè)三維耦合換熱 SG一、二次側(cè)流動(dòng)過(guò)程相互獨(dú)立，在耦合模塊中僅需考慮式(27)中能量源項(xiàng)SE在一、二次側(cè)的熱量傳遞。STAF-CT程序采用一、二次側(cè)重疊網(wǎng)格技術(shù)，兩側(cè)的管束區(qū)劃分完全相同的網(wǎng)格，且每個(gè)網(wǎng)格一一對(duì)應(yīng)；開(kāi)發(fā)了三維換熱耦合模塊，獲取管束區(qū)一、二次側(cè)對(duì)應(yīng)網(wǎng)格中溫度、速度等關(guān)鍵熱工參數(shù)，進(jìn)而求得換熱系數(shù)，并利用式(6)中單位面積平均熱流密度qA，通過(guò)轉(zhuǎn)換，獲得控制體內(nèi)的單位體積熱源，即能量源項(xiàng)SE，最終在兩側(cè)間實(shí)現(xiàn)熱量傳遞，如圖3所示。

圖3 兩側(cè)流域?qū)?yīng)網(wǎng)格間能量傳遞示意圖Fig.3 Energy source term transfer between corresponding grid on both sides

2.2 程序驗(yàn)證

STAF-CT程序同樣利用FRIGG棒束實(shí)驗(yàn)[36-37]，將程序計(jì)算值與實(shí)驗(yàn)值進(jìn)行對(duì)比，驗(yàn)證了程序準(zhǔn)確性，詳情參考Zhao等[40]的文獻(xiàn)。

2.3 基于STAF-CT程序的核電自然循環(huán)SG數(shù)值模型建立

相比STAF1.0程序中計(jì)算域僅包括下降段、管束區(qū)、管束支承板、提升段以及汽水分離器的部分筒體，STAF-CT程序拓展了一次側(cè)計(jì)算域，包括進(jìn)出口下腔室、管板以及管束區(qū)，如圖4所示。其中管束區(qū)網(wǎng)格為重疊網(wǎng)格，并采用多孔介質(zhì)模型簡(jiǎn)化，同時(shí)包含一次側(cè)管束區(qū)和二次側(cè)管束區(qū)；進(jìn)出口下腔室流域采用精細(xì)化建模，與之相接的管板采用多孔介質(zhì)模型簡(jiǎn)化。一次側(cè)主體建模采用精細(xì)化建模結(jié)合多孔介質(zhì)方法，在降低計(jì)算量的同時(shí)也保證了計(jì)算精度，實(shí)現(xiàn)了一回路冷卻劑經(jīng)下腔室進(jìn)入管束區(qū)流量分配計(jì)算。相比一次側(cè)一維簡(jiǎn)化，STAF-CT程序可較好地考慮一次側(cè)管束區(qū)流量分配不均對(duì)換熱分布的影響。

圖4 STAF-CT程序的SG網(wǎng)格模型Fig.4 Mesh model of SG for STAF-CT code

3 STAF3.0程序開(kāi)發(fā)

團(tuán)隊(duì)開(kāi)發(fā)的STAF1.0程序和STAF-CT程序均是基于商用CFD平臺(tái)ANSYS FLUENT二次開(kāi)發(fā)，為突破國(guó)外商用軟件對(duì)于模型開(kāi)發(fā)和修改代碼的限制，團(tuán)隊(duì)基于開(kāi)源CFD軟件平臺(tái)OpenFOAM將現(xiàn)有STAF程序進(jìn)行移植和深度開(kāi)發(fā)，并升級(jí)開(kāi)發(fā)了并行計(jì)算能力，顯著提升了程序的計(jì)算效率。下面對(duì)主要模型升級(jí)部分進(jìn)行介紹。

3.1 數(shù)學(xué)物理模型升級(jí)部分

3.1.1支承板阻力模型 STAF1.0和STAF-CT程序中，采用了多孔階躍模型模擬管束區(qū)支承板壓降損失。在STAF3.0程序中，將SG模型支承板處單層面網(wǎng)格復(fù)制成重疊兩層面網(wǎng)格，如圖5所示，兩個(gè)面按流動(dòng)方向分別命名為cyclic B面和cyclic A面，并設(shè)置為cyclic邊界。兩個(gè)cyclic邊界之間可添加自定義壓降模型，為實(shí)現(xiàn)多孔階躍模型，采用如下方程：

(29)

式中，Δm為支承板厚度，m。

圖5 OpenFOAM中多孔階躍模型實(shí)現(xiàn)方法Fig.5 Implementation of porous-jump model in OpenFOAM

3.1.2并行計(jì)算功能 STAF-CT程序開(kāi)發(fā)了一、二次側(cè)全三維耦合換熱模型，涉及重疊網(wǎng)格間數(shù)據(jù)交互。在并行計(jì)算時(shí)，網(wǎng)格會(huì)被自動(dòng)“分割”并分配給不同節(jié)點(diǎn)。因此，實(shí)現(xiàn)并行計(jì)算難點(diǎn)在于管束區(qū)網(wǎng)格不會(huì)被完整分配給1個(gè)節(jié)點(diǎn)，包含管束區(qū)網(wǎng)格節(jié)點(diǎn)需要找到對(duì)應(yīng)的一次側(cè)或二次側(cè)管束區(qū)網(wǎng)格進(jìn)行耦合計(jì)算。因此，團(tuán)隊(duì)提出了一種人工虛擬分配網(wǎng)格至各節(jié)點(diǎn)的數(shù)據(jù)處理方法，利用網(wǎng)格中坐標(biāo)信息分組將所分區(qū)網(wǎng)格中溫度等參數(shù)儲(chǔ)存在數(shù)組中，并通過(guò)事先匹配的一、二次側(cè)網(wǎng)格與所在數(shù)組中的位置信息可實(shí)現(xiàn)兩側(cè)網(wǎng)格一一對(duì)應(yīng)數(shù)據(jù)交換，并將數(shù)組中計(jì)算完成的數(shù)據(jù)返回給網(wǎng)格中，圖6示出一、二次側(cè)并行計(jì)算流程。圖7示出不同核數(shù)迭代求解SG模型(94萬(wàn)網(wǎng)格)100步耗時(shí)對(duì)比，可以看出從單核到6核并行計(jì)算時(shí)間大幅縮短，而繼續(xù)增加核數(shù)帶來(lái)的時(shí)間收益則逐漸減弱(多孔介質(zhì)網(wǎng)格量較少)。整體上并行計(jì)算功能大幅提升了程序計(jì)算效率。

圖6 一、二次側(cè)并行求解流程 Fig.6 Parallel computing process between primary side and secondary side

圖7 迭代時(shí)間隨分配核數(shù)的變化Fig.7 Iteration time change with core number

3.2 模型驗(yàn)證

STAF3.0程序基于開(kāi)源CFD軟件OpenFOAM開(kāi)發(fā)，其作為新平臺(tái)，利用FRIGG棒束實(shí)驗(yàn)[36-37]和MB-2[38]蒸汽發(fā)生器實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了更加詳細(xì)的驗(yàn)證對(duì)比，程序計(jì)算值與實(shí)驗(yàn)值進(jìn)行了對(duì)比驗(yàn)證，詳細(xì)數(shù)據(jù)參考He等[41]的文獻(xiàn)，最終驗(yàn)證了STAF3.0程序的準(zhǔn)確性與可靠性。

4 典型應(yīng)用

4.1 SG三維熱工水力特性分析

STAF系列程序基礎(chǔ)功能是求解SG三維熱工水力參數(shù)分布，圖8示出AP1000蒸汽發(fā)生器對(duì)稱截面上速度、溫度、壓力及空泡份額等關(guān)鍵熱工參數(shù)分布情況。以速度場(chǎng)為例，在SG熱側(cè)(右側(cè)為熱側(cè))，由于冷、熱側(cè)換熱量差異(圖9示出一、二次側(cè)換熱分布)，熱側(cè)濕蒸汽流速顯著大于冷側(cè)，且管束外圍間隙和管廊區(qū)流速遠(yuǎn)大于管束區(qū)域，詳細(xì)熱工水力分析可參考Zhao等[40]的成果。高精度熱工水力場(chǎng)可用于流致振動(dòng)以及腐蝕沉積分析計(jì)算，對(duì)SG長(zhǎng)期運(yùn)行過(guò)程安全分析至關(guān)重要。

圖8 SG對(duì)稱截面熱工參數(shù)分布 Fig.8 Thermal parameter distribution on SG symmetrical cross-section

4.2 帶軸流式預(yù)熱器SG模擬

與常規(guī)SG相比，帶軸流式預(yù)熱器可在保持相同換熱能力前提下，擁有更緊湊體積。其原因在于軸流式預(yù)熱器結(jié)構(gòu)將SG下降段分為冷側(cè)和熱側(cè)兩部分，并將冷態(tài)給水和再循環(huán)水重新按冷側(cè)100%給水+10%再循環(huán)水和熱側(cè)剩余90%再循環(huán)水分配，然后在管束區(qū)管廊處設(shè)置一定高度中央隔板，使來(lái)自下降段流體流動(dòng)換熱過(guò)程在管束中央隔板區(qū)域內(nèi)獨(dú)立進(jìn)行。在STAF程序中可通過(guò)網(wǎng)格標(biāo)記法實(shí)現(xiàn)對(duì)帶軸流式預(yù)熱器結(jié)構(gòu)模擬。如圖10所示，可以看出在熱側(cè)因供給為再循環(huán)水，一直處于飽和狀態(tài)；冷側(cè)主要為冷態(tài)給水，溫度更低，需要更長(zhǎng)流動(dòng)長(zhǎng)度達(dá)到飽和；兩種SG一、二次側(cè)溫差的差異如圖11所示，與常規(guī)SG相比，在熱側(cè)(橫坐標(biāo)0～0.5)帶軸流式預(yù)熱器SG一、二次側(cè)溫差小幅減??；但是在冷側(cè)(橫坐標(biāo)0.5～1)帶軸流式預(yù)熱器SG一、二次側(cè)溫差被大幅拉大，有效提升了換熱效率。針對(duì)帶軸流式預(yù)熱器SG三維熱工水力研究可為新一代高效率緊湊型SG設(shè)計(jì)提供設(shè)計(jì)參考。詳細(xì)研究可參考Zeng等[42]的工作。

圖9 SG兩側(cè)換熱源項(xiàng)分布Fig.9 Energy source term distribution between SG two sides

a——常規(guī)SG；b——帶軸流式預(yù)熱器SG圖10 溫度分布對(duì)比 Fig.10 Comparison of temperature distribution

4.3 堵管工況

核電廠蒸汽發(fā)生器U型管束在長(zhǎng)期運(yùn)行過(guò)程中難免會(huì)出現(xiàn)破損失效情況。少數(shù)傳熱管失效可通過(guò)堵管方式進(jìn)行處理。但堵管會(huì)對(duì)一次側(cè)管束內(nèi)流量分配造成影響，且堵管位置處損失了換熱能力，會(huì)降低SG整體換熱性能。STAF程序可通過(guò)采用網(wǎng)格標(biāo)記法實(shí)現(xiàn)對(duì)一次側(cè)不同堵管情況的模擬。圖12示出一次側(cè)在管束區(qū)內(nèi)環(huán)發(fā)生10%堵管時(shí)的速度分布，可發(fā)現(xiàn)在堵管區(qū)域右側(cè)管束流量分配變大，即最大流量管束位置因堵管發(fā)生了向右側(cè)移動(dòng)。圖13示出二次側(cè)在管束區(qū)內(nèi)環(huán)發(fā)生10%堵管時(shí)二次側(cè)空泡份額分布特性，結(jié)果表明：堵管區(qū)域?qū)Χ蝹?cè)冷側(cè)換熱區(qū)影響明顯，冷側(cè)空泡份額顯著降低，對(duì)冷側(cè)管束上方蒸汽區(qū)蒸汽品質(zhì)參數(shù)產(chǎn)生影響，具體分析可參考趙曉晗等[43]的研究成果。

a——常規(guī)SG；b——帶軸流式預(yù)熱器SG 圖11 兩側(cè)沿管束方向溫度分布Fig.11 Temperature distribution along tube bundle on two sides

圖12 堵管情況下一次側(cè)速度分布Fig.12 Velocity distribution of primary side under tube plugging

圖13 堵管情況下二次側(cè)空泡份額分布Fig.13 VOF distribution of secondary side under tube plugging

4.4 二次側(cè)16N遷移

16N核素是一回路16O的活化產(chǎn)物，其β衰變時(shí)會(huì)釋放兩種能量的γ射線，可通過(guò)探測(cè)上述衰變用于SG傳熱管裂紋完整性在線監(jiān)測(cè)，根據(jù)泄漏量制定SG維護(hù)方案。在STAF程序中植入薄壁微裂紋閃蒸相變模型，可實(shí)現(xiàn)傳熱管微裂紋泄漏情況下16N遷移擴(kuò)散計(jì)算能力，程序結(jié)果可為SG傳熱管微裂紋擴(kuò)大破損前在線監(jiān)測(cè)提供關(guān)鍵數(shù)據(jù)。圖14示出在一假定泄漏位置處發(fā)生泄漏時(shí)，16N在二次側(cè)分布情況，可看出從泄漏點(diǎn)沿管束向上流動(dòng)方向，16N質(zhì)量份額不斷降低，同時(shí)橫向擴(kuò)散在增強(qiáng)，特別是在橫流較強(qiáng)的彎管區(qū)域。圖15示出不同截面高度下出口處16N泄漏質(zhì)量份額最大值隨遷移時(shí)間的變化。結(jié)果表明：在泄漏位置小于9.0 m時(shí)(位于直管段)，其出口最大16N質(zhì)量份額相近；泄漏高度越高，其達(dá)到最大泄漏值所需時(shí)間越短，且位置在彎管區(qū)時(shí)，由于該位置有較強(qiáng)橫向湍流擴(kuò)散作用，在出口處有著遠(yuǎn)高于直管段位置的16N質(zhì)量份額。

圖14 16N泄漏的質(zhì)量份額和泄漏位置Fig.14 16N leaking quality share and leaking position

圖15 16N泄漏質(zhì)量份額最大值隨遷移時(shí)間的變化Fig.15 Change of maximum 16N leakage mass share with migration time

4.5 傳熱管腐蝕沉積

SG二次側(cè)由于腐蝕、泄漏等問(wèn)題產(chǎn)生雜質(zhì)會(huì)引起積污，嚴(yán)重影響SG運(yùn)行安全。因此，針對(duì)雜質(zhì)在SG二次側(cè)遷移沉積過(guò)程及其對(duì)熱工水力性能影響研究很有必要?；赟TAF程序，考慮雜質(zhì)輸運(yùn)過(guò)程和壁面沉積過(guò)程開(kāi)發(fā)了蒸汽發(fā)生器二次側(cè)雜質(zhì)遷移沉積模型，可用于計(jì)算沉積物在SG二次側(cè)的分布特性。圖16示出不含污垢SG經(jīng)過(guò)300 d沉積所得到沉積物質(zhì)量分布。結(jié)果表明：由于沸騰作用，熱側(cè)比冷側(cè)沉積物質(zhì)量更大，且沿著管束向上沉積量也在增加，彎管區(qū)沉積更加明顯，這是由于管束上部液相汽化較多，雜質(zhì)被濃縮沉積。管束表面沉積污垢會(huì)對(duì)換熱造成影響，圖17示出沉積前后SG一、二次側(cè)換熱量，受污垢影響，熱側(cè)換熱量降低，雖然冷側(cè)換熱量增大補(bǔ)償了一部分熱側(cè)換熱量損失，但SG整體換熱量下降。詳細(xì)分析可參考Mu等[44]的成果。

圖16 沉積300 d后沉積物的分布Fig.16 Distribution of deposit after 300 d of deposition

4.6 海洋條件下SG熱工水力響應(yīng)特性

在海洋條件下核動(dòng)力系統(tǒng)會(huì)產(chǎn)生六自由度的運(yùn)動(dòng)，其引入的附加力會(huì)對(duì)SG一、二次側(cè)流動(dòng)換熱特性產(chǎn)生影響，進(jìn)而對(duì)SG二次側(cè)空泡份額、溫度、流速等關(guān)鍵熱工參數(shù)分布造成影響，因此有必要對(duì)運(yùn)動(dòng)條件下核動(dòng)力系統(tǒng)SG三維熱工水力特性開(kāi)展研究?；赟TAF程序，在動(dòng)量方程中引入不同方向附加力模型可實(shí)現(xiàn)對(duì)海洋條件下的模擬。圖18為5 s搖擺周期下，管束區(qū)局部位置(冷側(cè)和熱側(cè))不同搖擺幅度下橫向流速的變化，結(jié)果表明：橫向流速隨搖擺幅值增加而增加，且發(fā)生方向的周期性改變，導(dǎo)致傳熱管流致振動(dòng)風(fēng)險(xiǎn)的增加，詳細(xì)分析可參考文獻(xiàn)[43]。

a——沉積前；b——沉積300 d后圖17 沉積前后換熱量分布Fig.17 Heat transfer distribution before and after deposition

a——SG熱側(cè)；b——SG冷側(cè)圖18 不同搖擺角度下橫向流速的變化Fig.18 Velocity variation under different swing angles

5 總結(jié)與展望

本文對(duì)西安交通大學(xué)核反應(yīng)堆熱工水力團(tuán)隊(duì)開(kāi)發(fā)的SG三維熱工水力分析系列程序STAF發(fā)展歷程和主要數(shù)學(xué)物理模型進(jìn)行了闡述。STAF系列程序發(fā)展歷程如圖19所示。STAF1.0基于多孔介質(zhì)方法和漂移流模型開(kāi)發(fā)，使國(guó)內(nèi)擁有了可與ATHOS程序?qū)?biāo)的自主化蒸汽發(fā)生器三維熱工水力程序；STAF-CT程序進(jìn)一步拓展了計(jì)算域和功能，實(shí)現(xiàn)了SG一、二次側(cè)全三維耦合及瞬態(tài)分析能力；STAF3.0程序?qū)TAF-CT程序移植到OpenFOAM開(kāi)源平臺(tái)，升級(jí)并行計(jì)算功能，實(shí)現(xiàn)了代碼自主化，避免了國(guó)外商用軟件二次開(kāi)發(fā)限制。另一方面，對(duì)STAF系列程序典型應(yīng)用場(chǎng)景進(jìn)行了介紹，如SG三維熱工水力穩(wěn)態(tài)及瞬態(tài)分析、16N遷移、腐蝕沉積和海洋條件等。STAF系列程序功能全面、應(yīng)用廣泛，可為我國(guó)自主化核動(dòng)力系統(tǒng)SG三維熱工水力分析程序開(kāi)發(fā)提供重要借鑒，也為蒸汽發(fā)生器設(shè)計(jì)優(yōu)化以及安全分析提供工具支撐。

圖19 核動(dòng)力蒸汽發(fā)生器三維熱工水力分析程序STFA系列研發(fā)歷程Fig.19 History of STAF series at XJTU NuTHeL

當(dāng)前STAF程序采用漂移流模型，簡(jiǎn)化了兩相間復(fù)雜作用力，且在瞬態(tài)計(jì)算中會(huì)忽略相間熱力學(xué)不平衡效應(yīng)，同時(shí)多孔介質(zhì)模型也限制了對(duì)于局部流場(chǎng)的精細(xì)化捕捉能力。團(tuán)隊(duì)持續(xù)對(duì)STAF系列程序進(jìn)行功能升級(jí)和模型完善，正在開(kāi)展全SG流域、多功能、多物理場(chǎng)耦合等領(lǐng)域的研究工作，對(duì)STAF程序幾點(diǎn)展望如下。

1) 先進(jìn)兩相流模型。目前團(tuán)隊(duì)正在開(kāi)發(fā)以兩流體模型為基礎(chǔ)的蒸汽發(fā)生器三維熱工水力程序，將持續(xù)提升STAF程序在核動(dòng)力系統(tǒng)蒸汽發(fā)生器設(shè)計(jì)及研發(fā)過(guò)程中的適用性。

2) 多物理場(chǎng)耦合能力?；赟TAF平臺(tái)，結(jié)合局部精細(xì)化分析方法，進(jìn)一步引入流致振動(dòng)、傳熱管磨損退化等物理場(chǎng)，最終實(shí)現(xiàn)流場(chǎng)、振動(dòng)和磨損多物理場(chǎng)耦合計(jì)算。

3) 跨維度耦合。綜合系統(tǒng)程序和CFD程序的優(yōu)點(diǎn)，開(kāi)發(fā)系統(tǒng)程序與STAF程序的耦合模塊，實(shí)現(xiàn)SG與反應(yīng)堆一回路系統(tǒng)耦合分析，開(kāi)展反應(yīng)堆系統(tǒng)變工況對(duì)SG二次側(cè)熱工參數(shù)的影響研究。

4) 數(shù)字孿生技術(shù)?；谏疃葘W(xué)習(xí)方法和模型降階，團(tuán)隊(duì)正在開(kāi)展核動(dòng)力系統(tǒng)蒸汽發(fā)生器數(shù)字孿生技術(shù)研發(fā)，實(shí)現(xiàn)蒸汽發(fā)生器設(shè)計(jì)與安全分析的快速迭代和全壽命周期狀態(tài)預(yù)測(cè)。