(上海核工程研究設(shè)計(jì)院有限公司，上海 200233)

0 引言

蒸汽發(fā)生器作為核電廠一、二回路的樞紐，是核島最重要的設(shè)備之一。其是否安全運(yùn)行直接關(guān)系到整個(gè)核電廠的安全。蒸汽發(fā)生器二次側(cè)流場(chǎng)數(shù)據(jù)是蒸汽發(fā)生器設(shè)計(jì)和傳熱管流致振動(dòng)和磨損分析的重要輸入。由于蒸汽發(fā)生器結(jié)構(gòu)復(fù)雜、體積龐大，采用試驗(yàn)手段獲得二次側(cè)流場(chǎng)數(shù)據(jù)代價(jià)太大，難以實(shí)現(xiàn)；同樣由于結(jié)構(gòu)復(fù)雜，直接針對(duì)管束建模分析很難實(shí)現(xiàn)。因此本文將二次側(cè)管束區(qū)作為多孔介質(zhì)，引入分布阻力概念，對(duì)蒸汽發(fā)生器二次側(cè)三維兩相熱工水力分析程序進(jìn)行開(kāi)發(fā)，并開(kāi)展小規(guī)模管束的流動(dòng)換熱試驗(yàn)，以完成對(duì)程序的試驗(yàn)驗(yàn)證。

最早將多孔介質(zhì)模型應(yīng)用于換熱器數(shù)值模擬分析的是Patankar等[1]，他們于1974年將管殼式換熱器的殼側(cè)作為多孔介質(zhì)，完成了對(duì)管殼式換熱器殼側(cè)流場(chǎng)的數(shù)值模擬。在此基礎(chǔ)上，Sha 等[2]模擬了SG和反應(yīng)堆堆芯中冷卻劑的流動(dòng)；Prithiviraj等[3-4]模擬了三維換熱器中的流動(dòng)。Ferng等[5-7]采用CFX 軟件多孔介質(zhì)模型對(duì)SG管束區(qū)的二次側(cè)流場(chǎng)進(jìn)行了分析。1977年美國(guó)電力研究院(EPRI)開(kāi)始開(kāi)發(fā)蒸汽發(fā)生器熱工水力計(jì)算程序，后來(lái)通過(guò)試驗(yàn)驗(yàn)證和實(shí)堆反饋數(shù)據(jù)的不斷修正，于1984年正式推出了適用于U形管蒸汽發(fā)生器(UTSG)和直流蒸汽發(fā)生器(OTSG)二次側(cè)三維流場(chǎng)分析的CFD程序ATHOS[8]。Chan等[9]于1986年開(kāi)發(fā)了UTSG穩(wěn)態(tài)、瞬態(tài)三維流場(chǎng)分析程序PORTHOS。法國(guó)電力公司(EDF)為研究PWR和LMFBR堆芯和組件內(nèi)的單相、兩相流場(chǎng)而開(kāi)發(fā)了THYC程序[10]，之后推出適用于傳熱器流場(chǎng)分析的版本THYC-EXCHANGER[11]來(lái)分析蒸汽發(fā)生器二次側(cè)管束區(qū)兩相流場(chǎng)。法國(guó)原子能委員會(huì)(CEA)開(kāi)發(fā)并驗(yàn)證了UTSG三維熱工水力分析程序GENEPI[12]?？梢钥闯?，國(guó)際上已有多個(gè)蒸汽發(fā)生器三維熱工水力分析程序，這些軟件均基于多孔介質(zhì)模型，只是其兩相流控制方程、管內(nèi)外傳熱、相變及流動(dòng)阻力計(jì)算方法有所不同。國(guó)內(nèi)也有學(xué)者[13-16]采用CFD軟件對(duì)蒸汽發(fā)生器二次側(cè)進(jìn)行三維熱工水力分析，大多屬于方法性的探索研究，未形成專用分析程序。因此，本文基于已有的理論研究基礎(chǔ)，彌補(bǔ)現(xiàn)有分析程序的缺陷，開(kāi)發(fā)可用于壓水堆蒸汽發(fā)生器三維兩相熱工水力分析程序。

1 程序開(kāi)發(fā)

與常規(guī)流體的兩相流動(dòng)類似，宏觀尺度的多孔介質(zhì)內(nèi)的兩相流動(dòng)研究方法也分為三種：均相流模型、漂移流模型和兩流體模型。在早期的研究工作中，為簡(jiǎn)化計(jì)算，多采用均相流模型，如CALIPSOS[17-18]、THEDA-1[19]、THEDA-2[20-21]等程序以及Ferng等[22-23]的一系列工作。均相流模型將兩相流體等效為單相流體，不能描述兩相間的相對(duì)速度，因此在計(jì)算豎直通道內(nèi)兩相流動(dòng)的時(shí)候，會(huì)出現(xiàn)空泡份額偏大的問(wèn)題。從理論上來(lái)說(shuō)，兩流體模型可以準(zhǔn)確地描述兩相間的熱力學(xué)和水力學(xué)不平衡特性，但由于管束外流動(dòng)的相間作用力及相間傳熱等模型目前尚不完備，采用該模型需引入較多假設(shè)，導(dǎo)致兩流體模型的計(jì)算精度不能保證。漂移流模型同樣將兩相流體視為混合相，但是為描述相間的速度差，在均相方程的基礎(chǔ)上引入相間滑移速度。雖然漂移流模型從理論上相對(duì)于兩流體模型精度較差，但由于該模型需引入的假設(shè)條件較少，計(jì)算精度可滿足SG內(nèi)模擬的要求[8]。因此本程序開(kāi)發(fā)中采用漂移流模型。多孔介質(zhì)內(nèi)的漂移流模型控制方程如下。

質(zhì)量守恒方程：

(1)

式中t——時(shí)間，s；

β——孔隙率；

ρm——混合物密度，kg/m3,ρm=αgρg+αlρl；

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

(2)

式中p——壓力，Pa；

μm,eff——混合物有效黏度，Pa·s,

μm=αgμg+αlμl;

αg——汽相體積份額；

ρg,ρl——汽相、液相密度，kg/m3；

αl——液相體積份額；

能量守恒方程：

(3)

式中Hg——汽相焓，J/kg；

Hl——液相焓，J/kg；

km——混合物導(dǎo)熱系數(shù)，W/(m·K),km=αgkg+αlkl;

SE——能量源項(xiàng),W/m3。

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

(4)

式中Sg——汽相質(zhì)量源項(xiàng),kg/(m3·s)。

程序動(dòng)量方程源項(xiàng)中考慮的阻力主要由以下部件引入，如傳熱管束、汽水分離器、下降通道、支承板等，其中管束、汽水分離器和下降通道的阻力以分布阻力形式添加到這些部件所在網(wǎng)格中，支承板阻力以集中阻力形式添加到支承板所在位置的網(wǎng)格界面上。根據(jù)流動(dòng)方向不同，將流動(dòng)方向沿管束軸向和橫向分解，在軸向和橫向分別引入動(dòng)量源項(xiàng)。對(duì)于順流和橫掠阻力分別采用管束外順流和橫掠[24-25]阻力經(jīng)驗(yàn)關(guān)系式計(jì)算。支承板、下降通道和汽水分離器的阻力系數(shù)，根據(jù)蒸汽發(fā)生器設(shè)計(jì)結(jié)構(gòu)直接給出。

能量方程的源項(xiàng)，即一次側(cè)向二次側(cè)釋熱，被簡(jiǎn)化為一維分布計(jì)算。程序采取一維-三維耦合傳熱計(jì)算方式。該方法一、二次側(cè)網(wǎng)格示意如圖1所示，計(jì)算思路如下：將一次側(cè)沿流體流動(dòng)方向劃分為J個(gè)控制體，局部網(wǎng)格編號(hào)用j表示；將二次側(cè)劃分為K個(gè)控制體，局部網(wǎng)格編號(hào)為k；通過(guò)DEFINE_INIT宏中的幾何處理模塊、根據(jù)網(wǎng)格位置將一、二次側(cè)網(wǎng)格進(jìn)行匹配。

一次側(cè)為單相液體對(duì)流換熱，采用Dittus-Boelter 公式[26]計(jì)算；二次側(cè)換熱包括單相對(duì)流換熱、過(guò)冷沸騰換熱及飽和沸騰換熱，對(duì)于單相對(duì)流換熱，由于目前的經(jīng)驗(yàn)關(guān)系式多為管外順流和橫掠關(guān)系式，無(wú)傾斜沖刷管束的換熱關(guān)系式，因此將流速分解為順流和橫掠流速，分別計(jì)算順流換熱系數(shù)[24-25]和橫掠換熱系數(shù)[27]，并將兩者的加權(quán)和[8]作為實(shí)際的斜掠換熱系數(shù)。在計(jì)算中，給定一次側(cè)入口的質(zhì)量流速和溫度，便可計(jì)算一、二次側(cè)間的流動(dòng)換熱。

圖1 程序的一、二次側(cè)網(wǎng)格示意

考慮到多孔介質(zhì)模型對(duì)局部參數(shù)的平均效應(yīng)，以及漂移流模型對(duì)相界面密度的依賴性較低，因此，程序中對(duì)沸騰模型進(jìn)行簡(jiǎn)化，僅根據(jù)能量平衡計(jì)算兩相流溫度及蒸發(fā)率，不計(jì)算相界面輸運(yùn)過(guò)程。引入熱量密度二分方法，將總熱流密度分為蒸發(fā)熱流密度和對(duì)流熱流密度，分別用于蒸發(fā)液相和加熱混合相。

空泡份額方程中的源項(xiàng)表示氣相質(zhì)量源項(xiàng)，其值等于蒸發(fā)率與冷凝率之差，對(duì)于穩(wěn)態(tài)計(jì)算忽略SG內(nèi)的非穩(wěn)態(tài)冷凝現(xiàn)象，氣相質(zhì)量源項(xiàng)等于液相蒸發(fā)率。

對(duì)于湍流模型，已公開(kāi)發(fā)表的文獻(xiàn)多是針對(duì)純流體區(qū)域的精細(xì)網(wǎng)格提出[28]，Teruel等[29-31]從宏觀尺度研究多孔介質(zhì)內(nèi)的湍流流動(dòng)，他們?cè)诤暧^網(wǎng)格內(nèi)對(duì)k-e方程進(jìn)行積分，得到宏觀k-e二方程湍流模型，并且分別采用微觀湍流模型和宏觀湍流模型計(jì)算幾何規(guī)則的多孔介質(zhì)內(nèi)的流動(dòng)，從而驗(yàn)證這些宏觀k-e湍流模型的正確性。但是這些模型均針對(duì)多孔介質(zhì)內(nèi)的單相流體進(jìn)行開(kāi)發(fā)，并不適用于多孔介質(zhì)內(nèi)的兩相流動(dòng)。因此，本程序開(kāi)發(fā)中采用了商用蒸汽發(fā)生器熱工水力分析程序ATHOS[32]和GENEPI[33]中采用的零方程湍流模型計(jì)算兩相混合物的有效黏度。

針對(duì)壓水堆蒸汽發(fā)生器結(jié)構(gòu)特點(diǎn)，以CFD軟件多孔介質(zhì)模型為基礎(chǔ)，采用UDF將以上傳熱、阻力模型添加到求解器中，并通過(guò)UDF求解一次側(cè)流場(chǎng)及一、二次側(cè)耦合傳熱特性，完成蒸汽發(fā)生器二次側(cè)三維兩相熱工水力分析程序的開(kāi)發(fā)。

2 程序驗(yàn)證

在本文開(kāi)發(fā)的程序用于工程之前，需對(duì)其進(jìn)行試驗(yàn)驗(yàn)證。由于蒸汽發(fā)生器結(jié)構(gòu)龐大，對(duì)原型蒸汽發(fā)生器進(jìn)行試驗(yàn)研究需要極大的蒸汽流量和加熱功率，成本很高。因此，基于?；治?，采用小規(guī)模管束進(jìn)行流動(dòng)傳熱試驗(yàn)，模擬蒸汽發(fā)生器一、二次側(cè)流場(chǎng)，獲得試驗(yàn)件一、二次側(cè)流場(chǎng)信息，同時(shí)采用該程序?qū)υ囼?yàn)件建模分析，預(yù)測(cè)試驗(yàn)工況下的試驗(yàn)件的一、二次流場(chǎng)數(shù)據(jù)，通過(guò)與試驗(yàn)測(cè)量結(jié)果的比較，達(dá)到對(duì)程序驗(yàn)證的目的。

2.1 蒸汽發(fā)生器二次側(cè)流動(dòng)傳熱試驗(yàn)

試驗(yàn)中采用40根傳熱管模擬二次側(cè)管束的流動(dòng)換熱狀態(tài)，管束橫截面如圖2所示。為準(zhǔn)確體現(xiàn)蒸汽發(fā)生器原型的熱工水力參數(shù)，試驗(yàn)件管束高度上與原型保持相等，試驗(yàn)?zāi)M體中采用的傳熱管材料、管徑、管間距、布置方式、管束長(zhǎng)度和支承板結(jié)構(gòu)與原型相同。試驗(yàn)中一、二次側(cè)進(jìn)口溫度、系統(tǒng)運(yùn)行壓力和循環(huán)倍率等參數(shù)與原型相同，分別為：傳熱管規(guī)格為?17.48 mm×1.01 mm，直段高度9.2 m，相鄰傳熱管的中心距24.89 mm，呈三角形排列，一次側(cè)壓力15.5 MPa，進(jìn)口溫度322 ℃；二次側(cè)壓力5.6 MPa，給水溫度226.7 ℃。試驗(yàn)?zāi)M體如圖3所示，包括上下筒體、上下封頭、管板、管束、下降管、管束支承板、汽水分離裝置等部件。為了避免測(cè)量探頭同時(shí)穿過(guò)筒體和套筒兩層壁面，采用兩根圓形下降管代替蒸汽發(fā)生器環(huán)形下降通道，同時(shí)在下降管上安裝了調(diào)節(jié)球閥，調(diào)節(jié)下降管阻力，以便更好地控制試驗(yàn)?zāi)M體運(yùn)動(dòng)的循環(huán)倍率，下降管與管束區(qū)入口的連接如圖4所示。試驗(yàn)體運(yùn)行及設(shè)計(jì)參數(shù)如表1,2所示。

圖2 管束橫截面示意

圖3 試驗(yàn)?zāi)M體示意

圖4 下降管與管束區(qū)入口連接示意

表1 試驗(yàn)體額定運(yùn)行參數(shù)

表2 試驗(yàn)體設(shè)計(jì)參數(shù)

試驗(yàn)中測(cè)量的參數(shù)包括流量、溫度、壓力、差壓、空泡份額、流體速度等。一次側(cè)進(jìn)、出口溫度通過(guò)在一次側(cè)進(jìn)、出口管道上安裝鎧裝熱電偶測(cè)量；一次側(cè)流量通過(guò)在一次側(cè)進(jìn)口管道上安裝孔板流量計(jì)測(cè)量。二次側(cè)出口蒸汽壓力在蒸發(fā)器出口管上利用壓力傳感器測(cè)量，采用在蒸汽出口管道上安裝的孔板流量計(jì)測(cè)量蒸汽出口流量，給水流量測(cè)量與蒸汽流量測(cè)量方法一致，在給水管道上安裝鎧裝熱電偶測(cè)量給水溫度。采用光纖探針和伽馬射線法測(cè)量空泡份額。通過(guò)測(cè)量得到如下宏觀參數(shù)：熱功率、一回路壓力、一回路進(jìn)、出口溫度、給水流量、給水溫度、給水壓力、蒸汽流量、蒸汽溫度、蒸汽壓力、二次側(cè)水位等。下降通道上的流量采用超聲波流量計(jì)進(jìn)行測(cè)量。在下降段底部布置2個(gè)熱電偶用于測(cè)量進(jìn)入管束區(qū)的流體溫度。

對(duì)于一次側(cè)沿程流體溫度，采用K型鎧裝熱電偶，熱電偶采用管接頭方式從傳熱管內(nèi)部引出，再采用管接頭將其引出筒體。安裝方式如圖5所示。管接頭作為鎧裝熱電偶穿過(guò)壓力邊界的密封元件，在與傳熱管焊接安裝前進(jìn)行焊接工藝試驗(yàn)，驗(yàn)證了密封可靠性。

圖5 一次側(cè)測(cè)溫?zé)犭娕脊芙宇^安裝方式

空泡份額的測(cè)量分別采用光纖探針?lè)ā①ゑR射線法和壓差法獲得局部點(diǎn)、截面和體積含汽率。

光纖探針?lè)ㄊ抢霉鈱W(xué)的全反射原理，傳輸光在光纖探針敏感頭界面與媒質(zhì)界面處符合全反射條件時(shí)，這一部分光將返回敏感頭內(nèi)部，如果敏感頭接觸的媒質(zhì)折射率產(chǎn)生變化時(shí)，由于可能將不滿足全反射條件，此時(shí)不能產(chǎn)生全反射現(xiàn)象。通過(guò)檢測(cè)經(jīng)敏感頭反射回接收端的光強(qiáng)的強(qiáng)弱變化而產(chǎn)生的高低不同的電信號(hào)來(lái)分辨光纖探針敏感頭測(cè)量到的是氣相還是液相，從而根據(jù)產(chǎn)生的連續(xù)變化的電壓信號(hào)測(cè)量局部截面持氣率，光纖探針敏感頭測(cè)量原理圖如圖6所示。試驗(yàn)中測(cè)點(diǎn)位置見(jiàn)圖7。具體測(cè)點(diǎn)位置和探針數(shù)量見(jiàn)表3。

射線法是一種非接觸式的測(cè)量方法，其原理是根據(jù)汽、液兩種物質(zhì)對(duì)射線的吸收程度不同推算出射線照射區(qū)域的含汽率。該方法不會(huì)破壞管道中的流場(chǎng)和溫度場(chǎng)的自然分布，適用于高溫高壓條件下測(cè)量，原理如圖8所示。試驗(yàn)測(cè)點(diǎn)位置如圖9所示。

圖6 光纖探針?lè)y(cè)量含汽率的原理示意

圖7 空泡份額測(cè)點(diǎn)布置

表3 光纖探針空泡份額測(cè)點(diǎn)位置

圖8 射線法原理示意

圖9 射線法測(cè)點(diǎn)位置示意

壓差法則是通過(guò)測(cè)量?jī)蓚€(gè)取壓截面的壓力差，計(jì)算兩個(gè)截面之間的流體體積含汽率。

為了獲得較多的試驗(yàn)數(shù)據(jù)，試驗(yàn)工況圍繞額度工況，分別改變一次側(cè)進(jìn)出口溫度、給水溫度、加熱功率和循環(huán)倍率等，完成了11個(gè)試驗(yàn)工況。

試驗(yàn)測(cè)量的一、二次側(cè)溫度分布如圖10所示。三種方法空泡份額測(cè)量結(jié)果見(jiàn)表4。

2.2 試驗(yàn)件模擬分析

采用本文開(kāi)發(fā)的程序?qū)υ囼?yàn)件進(jìn)行數(shù)值模擬分析，預(yù)測(cè)二次流場(chǎng)數(shù)據(jù)。計(jì)算中選取試驗(yàn)件管板以上至分離器入口為分析對(duì)象。對(duì)計(jì)算區(qū)域進(jìn)行建模，管束區(qū)采用多孔介質(zhì)模型進(jìn)行模擬，為了方便設(shè)置支承板阻力系數(shù)，支承板阻力采用多孔階躍邊界描述，在邊界上設(shè)定支承板的阻力系數(shù)參數(shù)和厚度。為了保證數(shù)值計(jì)算穩(wěn)定性與準(zhǔn)確性，采用一體化網(wǎng)格，對(duì)整個(gè)計(jì)算區(qū)域進(jìn)行建模，然后劃分網(wǎng)格，所劃分網(wǎng)格全為六面體網(wǎng)格，網(wǎng)格總數(shù)為196 107，計(jì)算模型網(wǎng)格如圖11所示。

(a)一次側(cè)溫度分布

(b)二次側(cè)溫度分布

表4 空泡份額測(cè)量結(jié)果對(duì)比

根據(jù)試驗(yàn)件和試驗(yàn)工況，計(jì)算中入口為下部開(kāi)口，出口為分離器入口，入口邊界設(shè)為速度入口，不同工況給出不同入口流速，出口邊界設(shè)為壓力出口。計(jì)算中采用COUPLE算法求解控制方程組。本分析采用UDF輸入自定義模塊，如附加源項(xiàng)、各附加阻力項(xiàng)，因此計(jì)算過(guò)程中，除求解控制方程，還需調(diào)用UDF，在User-Defined Profile定義計(jì)算區(qū)域孔隙率；在User-Defined Init中定義熱源分布、網(wǎng)格潤(rùn)濕面積，以及其他需要在計(jì)算進(jìn)行之前計(jì)算的變量；在User-Defined Adjust中計(jì)算求解阻力源項(xiàng)所需的各個(gè)變量，如各控制體的阻力系數(shù)。

計(jì)算結(jié)果如圖12所示。為了圖片更加清晰，圖片在高度方向縮短了15倍。

圖11 計(jì)算模型網(wǎng)格

圖12 工況1的二次側(cè)溫度云圖及氣相體積份額云圖

2.3 程序?qū)Ρ闰?yàn)證

為了驗(yàn)證程序分析的準(zhǔn)確性，將相同工況下程序預(yù)測(cè)結(jié)果與試驗(yàn)測(cè)量結(jié)果進(jìn)行比較。二次側(cè)溫度進(jìn)行了歸一化處理，用測(cè)量點(diǎn)的溫度除以該工況下的飽和溫度，熱側(cè)和冷側(cè)歸一化比較結(jié)果如圖13所示。空泡份額比較如圖14所示。

(a)熱側(cè)

(b)冷側(cè)

(a)射線法測(cè)量的空泡份額

(b)探針?lè)y(cè)量的二次側(cè)空泡份額

3 結(jié)論

本文針對(duì)壓水堆核電廠蒸汽發(fā)生器結(jié)構(gòu)特點(diǎn)，以CFD軟件中多孔介質(zhì)模型為基礎(chǔ)，采用UDF對(duì)CFD軟件進(jìn)行了二次開(kāi)發(fā)，完成了蒸汽發(fā)生器二次側(cè)三維兩相熱工水力分析程序的開(kāi)發(fā)。并基于模化分析，采用小規(guī)模管束進(jìn)行流動(dòng)傳熱試驗(yàn)，模擬蒸汽發(fā)生器一、二次側(cè)流場(chǎng)，獲得試驗(yàn)體試驗(yàn)工況下的一、二次側(cè)流場(chǎng)數(shù)據(jù)，采用該程序?qū)υ囼?yàn)件建模分析。結(jié)果表明：程序能夠?qū)崿F(xiàn)對(duì)蒸汽發(fā)生器試驗(yàn)件的模擬，預(yù)測(cè)結(jié)果與試驗(yàn)測(cè)量結(jié)果符合良好，證明了該程序計(jì)算結(jié)果的可靠性。

未來(lái)蒸汽發(fā)生器二次側(cè)三維兩相熱工水力分析程序還需開(kāi)展進(jìn)一步工作。經(jīng)后續(xù)不斷改進(jìn)完善，程序?qū)?shí)現(xiàn)在工程上的應(yīng)用，該程序可成為我國(guó)核電廠蒸汽發(fā)生器設(shè)計(jì)研發(fā)的重要工具。