陳 林，張 昊，程以炫，郭英冉，馬齊超，周帆帆，楊燕華，趙 萌

(上海交通大學(xué) 核科學(xué)與工程學(xué)院，上海 200240)

壓水堆的堆芯中存在著數(shù)百個燃料組件棒束，冷卻劑存在著沿棒束方向的垂直流動和燃料棒間隙的橫向流動，當(dāng)反應(yīng)堆壓力容器發(fā)生破損事故時，堆芯內(nèi)可能出現(xiàn)劇烈相變并導(dǎo)致冷卻劑流動劇烈波動，繼而影響堆芯冷卻劑的傳熱能力使燃料棒包殼溫度驟升，并引發(fā)潛在的沸騰危機，甚至導(dǎo)致燃料包殼失去完整性[1]。子通道程序通過子通道模型考慮了流體在通道間的流動[2]，并通過一系列經(jīng)驗關(guān)系式，可以較好地預(yù)測壓力容器內(nèi)熱工水力現(xiàn)象，實現(xiàn)對事故下反應(yīng)堆的安全性能準(zhǔn)確評估。

國際上具有代表性的子通道程序包括基于均相流或漂移流模型的COBRA-ⅢC[3]、COBRA-Ⅳ[4]、VIPRE-01[5]、THINC[6]，基于兩流體六方程的VIPRE-02[7]及基于兩流體八方程的COBRA-TF[8]及WCOBRA[9]等。均相流及漂移流模型忽略了相間存在的熱力及流動不平衡，在低壓、高流率等工況下，計算結(jié)果與實際工況具有較大差異。兩流體子通道程序考慮了相間作用，實現(xiàn)程序?qū)峁にW(xué)現(xiàn)象良好的模擬。然而，大多數(shù)子通道程序在模擬變化劇烈的瞬態(tài)過程及強烈擾動的工況時，其計算結(jié)果往往與真實情況存在較大的差距。

國家電投集團(tuán)牽頭研發(fā)的一體化核電廠設(shè)計與安全分析軟件包[10](COSINE)于2011年6月立項。COSINE軟件包中均相流子通道程序cosSUBC[11]已經(jīng)開發(fā)完成并進(jìn)行了一系列測試及應(yīng)用，但該程序采用了均相流模型，計算結(jié)果存在著較大的保守度，導(dǎo)致電廠設(shè)計的經(jīng)濟(jì)性較差。因此，為精確分析多相場的熱工水力學(xué)現(xiàn)象，提高程序的計算精度，本文在COSINE均相流子通道程序的基礎(chǔ)上，開發(fā)新的多相場子通道程序。本文首先介紹COSINE軟件包中的多相場子通道程序的計算模型與求解方法，而后采用經(jīng)典算例及實驗數(shù)據(jù)評估程序計算的準(zhǔn)確性。

1 多相場子通道程序的基本方程及物理模型

多相場子通道程序采用三相八方程子通道模型，將液滴相[12]進(jìn)行單獨考慮。程序采用控制體網(wǎng)格與動量網(wǎng)格的交錯網(wǎng)格方法。程序中包含經(jīng)驗關(guān)系模型以完成方程的封閉，目前，程序采用的經(jīng)驗關(guān)系模型主要基于COSINE軟件包[10]中自主開發(fā)的本構(gòu)關(guān)系式及參考程序[13]中本構(gòu)關(guān)系模型。

1.1 基本方程

1) 質(zhì)量守恒方程

質(zhì)量守恒方程位于控制體網(wǎng)格上，包括氣相、連續(xù)液相、液滴相的3個連續(xù)方程。

(1)

式中：α為相份額；ρ為密度；u為軸向速度；w為間隙橫向速度；x為軸向網(wǎng)格長度；nk為間隙數(shù)量；i代表不同相場。方程右側(cè)1～3項分別為相變及夾帶沉積、壁面換熱、交混導(dǎo)致的各相質(zhì)量變化。

2) 動量守恒方程

動量守恒方程作用于軸向和間隙動量網(wǎng)格中，包括氣相、連續(xù)液相、液滴相共6個動量方程。

軸向動量守恒方程為：

(2)

橫向動量守恒方程為：

(3)

式(2)右側(cè)1～6項分別為壓力項、重力項、壁面摩擦項、相間拖曳力項、相間傳質(zhì)導(dǎo)致動量變化項、交混導(dǎo)致的動量變化；式(3)右側(cè)1～4項分別為壓力項、壁面阻力項、相間拖曳力項、相間傳質(zhì)導(dǎo)致動量變化項。

3) 能量守恒方程

能量守恒方程作用于控制體網(wǎng)格中，程序中假定連續(xù)液相與液滴相處于熱平衡態(tài)，因此將連續(xù)液相與液滴相的能量方程合并為全液相能量方程，程序中包括氣相與液相2個能量守恒方程。

(4)

方程右側(cè)1～5項分別為壓力項、相間傳質(zhì)引發(fā)的能量變化項、相間換熱項、壁面換熱項、交混導(dǎo)致的能量變化項。

1.2 物理模型

物理模型用于封閉多相場子通道程序守恒方程，程序中包含的物理模型包括：流型圖及相間作用模型、阻力模型、通道交混模型、液滴模型及壁面換熱模型。

1) 流型圖及相間作用模型

多相場子通道程序的相間作用模型依賴于不同的流型，程序流型圖中包括泡狀流、彈狀流、交混流、環(huán)霧流、反環(huán)狀流、反彈狀流、下降液膜流，流型判別方法列于表1。相間作用模型包括相間拖曳力模型與相間換熱模型。程序中認(rèn)為連續(xù)液相與液滴相無直接接觸，因此相間作用模型不含連續(xù)液相與液滴相的相間作用。

相間拖曳力與相間換熱本構(gòu)關(guān)系式分別為：

τij.i=CijρiA?ij|uj-ui|(uj-ui)

(5)

qij,i=htcij,iA?ij(Tsat-Ti)

(6)

式中：τij,i為i與j相的相間摩擦作用在i相上的力；Cij為相間阻力系數(shù)，計算方法見表1，表中后綴GL代表氣相與連續(xù)液相相間作用，GD代表氣相與液滴相相間作用；ρi為各相密度；A?ij為相界面面積；ui為當(dāng)前相速度；uj為相對相的速度；qij,i為i相從相界面獲得的能量；Tsat為飽和溫度；Ti為相溫度；htcij,i為i相相間換熱系數(shù)，相間換熱系數(shù)包括過冷液、過熱液、過冷氣與過熱氣4種狀態(tài)，由于過熱液與過冷氣在物理上并不穩(wěn)定，會迅速達(dá)到飽和態(tài)，因此程序中對過熱液與過冷氣設(shè)置了很大的換熱系數(shù)，過冷液與過熱氣換熱系數(shù)計算方法見表1。

2) 阻力模型

表1 流型圖及相間作用模型Table 1 Flow regime and interphase model

(7)

式中：Cw為壁面摩擦系數(shù)，采用McAdam關(guān)系式[14]；Cf為形阻系數(shù)，由用戶輸入；G為流體的質(zhì)量流率；φ為兩相接觸因子；Dh為通道的水力學(xué)直徑；Δx為控制體長度；液滴相彌散于氣相中且不于壁面直接接觸，不考慮液滴相的阻力。

3) 交混模型

湍流交混是由于流體脈動時自然渦團(tuán)擴散引起的非定向交混，交混模型是子通道程序特有的模型，其基于擴散效應(yīng)，表征了相鄰?fù)ǖ篱g的質(zhì)量、動量以及能量的擴散現(xiàn)象。

質(zhì)量、動量及能量交混關(guān)系式分別為：

(8)

(9)

(10)

式中：下標(biāo)ab、a和ab、b分別代表相鄰?fù)ǖ繿、b作用于通道a和b的量；h為各相焓；u′為脈動速度；Agap為間隙面積；Vc為控制體體積。

脈動速度u′可表示為：

(11)

4) 液滴模型

在傳統(tǒng)的中學(xué)地理教學(xué)活動中，我們老師可能或多或少的忽視了中學(xué)生地理學(xué)習(xí)中的知識情感、想象以及知識點的領(lǐng)悟等多方面的發(fā)展，。在地理課堂教學(xué)活動中老師過多地強調(diào)知識的記憶、模仿，抑制了學(xué)生的頭腦意識，表達(dá)動手能力，導(dǎo)致學(xué)生地理思維能力欠缺、地理課堂活力不足，從根本上制約了學(xué)生的主動性和創(chuàng)造性,從而讓我的地理課堂教學(xué)變得固執(zhí)、沉悶，缺少了靈性，缺乏了活力。在新課改的要求下，教師必須改變教學(xué)方法使地理課堂教學(xué)活力從根本上得到提升。

液滴存在于臨界后流型及臨界前的交混流、環(huán)霧流中，這些流型的計算中需要考慮液滴模型。液滴模型包括夾帶沉積模型及液滴相界面面積模型。其中，夾帶與沉積模型表征連續(xù)液相與液滴相間的質(zhì)量傳遞；液滴面積輸運方程用于液滴相界面面積的計算，用于相間本構(gòu)關(guān)系式。

夾帶與沉積模型為：

Γd,ld=-Γl,ld=Γentr-Γdepo

(12)

式中：Γd,ld為液滴生成率；Γl,ld為液滴向連續(xù)液相中的凈傳質(zhì)；Γentr為液滴夾帶率；Γdepo為液滴沉積率。程序中液滴夾帶沉積模型采用參考程序[13]的方法，夾帶包括在環(huán)霧狀流夾帶、反環(huán)狀流夾帶、下降液膜流夾帶，沉積主要包括環(huán)霧狀流沉積、橫向流沉積以及面積變化帶來的沉積。

程序中采用的液滴面積模型包括液滴相界面輸運方程與臨界韋伯?dāng)?shù)模型兩種方法。液滴相界面輸運模型為：

(13)

根據(jù)韋伯?dāng)?shù)(We)的定義，可得Wecrit與液滴直徑Dd的關(guān)系：

(14)

程序中假定液滴為圓球形，因此，液滴相界面面積濃度A?gd與Dd的關(guān)系為：

A?gd=6αd/Dd

(15)

式(14)、(15)為臨界韋伯?dāng)?shù)液滴面積模型，式中：σ為表面張力；ugd為氣相與液滴相對速度；ρd為液滴密度；αd為液滴份額。

5) 壁面換熱模型

壁面換熱模型表征壁面與流體間的傳熱關(guān)系。壁面換熱模型依賴于不同的壁面換熱模式，程序中包括8種壁面換熱模式，換熱模式及對應(yīng)的傳熱系數(shù)模型列于表2。壁面換熱會導(dǎo)致流體發(fā)生相變，程序中僅考慮壁面換熱所導(dǎo)致的液體蒸發(fā)效應(yīng)，而不考慮壁面過冷導(dǎo)致的蒸汽冷凝現(xiàn)象。壁面換熱關(guān)系式為：

qw,i=htcw,i(Tw-Ti)

(16)

(17)

式中：htcw,i為壁面向各相的傳熱系數(shù)，計算方法列于表2；Tw為壁面溫度；qw,i為壁面向流體的傳熱量；式(17)為壁面換熱導(dǎo)致質(zhì)量傳遞，htcw,fg為壁面相變換熱系數(shù)；hfg為汽化潛熱；Msliq為網(wǎng)格中的液相總質(zhì)量。為防止壁面換熱導(dǎo)致液相蒸發(fā)過大，設(shè)定為單個時間步內(nèi)的網(wǎng)格蒸發(fā)量不超過該網(wǎng)格連續(xù)液相總質(zhì)量的80%。

表2 流型圖及傳熱系數(shù)模型Table 2 Flow regime and wall heat transfer coefficient model

6) 多相場子通道程序的計算流程

程序采用半隱式計算方法及Newton-Raphson求解器，并通過牛頓下山法提高計算收斂性，采用殘差限值和庫朗數(shù)限值判斷計算是否收斂，求解流程如圖1所示。

在每個時間步內(nèi)，程序首先進(jìn)行預(yù)處理并計算出求解器所需的各項，而后使用Newton-Raphson求解器求解新時刻的各主變量，其中，多相場子通道程序主變量包括11個主變量，分別為體積分?jǐn)?shù)(αg、αd)，軸向與橫向速度(ug、ul、ud、wg、wl、wd)，焓(hv、hf)和壓力(p)；隨后使用牛頓下山法判斷是否需要松弛迭代，最后使用收斂準(zhǔn)則判斷程序是否計算收斂，并調(diào)整時間步長，進(jìn)入下一次計算。

圖1 程序計算流程 Fig.1 Calculation process of code

2 多相場子通道程序的評估

選取典型測試算例及實驗工況對程序進(jìn)行評估，驗證程序計算的合理性與正確性。

2.1 液滴模型測試

采用參考程序COBRA進(jìn)行對照，評估多相場子通道程序?qū)σ旱涡袨榈挠嬎隳芰?，并進(jìn)一步分析液滴尺寸對流場分布的影響。

1) 液滴模型評估

選取長度為1 m，直徑為2 cm的單根圓管算例進(jìn)行測試。如圖2所示，通道劃分為20個控制體，底部入口為氣相份額為0.925的流動邊界，入口處氣相與連續(xù)液相速度分別為5.18 m/s及4.1 m/s，各相焓均采用飽和值，頂部出口為壓力邊界，壓力為10 MPa。觀察計算穩(wěn)定后通道中流場分布情況。

液滴相彌散于氣相中且具有較大的相界面面積，因此液滴相與氣相間存在著較強的相間作用，圖3展示了液滴測試中本程序與參考程序沿軸向通道的流場分布對照，本程序與參考程序的計算趨勢符合良好。其中圖3a展示了液滴模型開閉情況下氣相速度的分布情況，開啟液滴模型后，液滴份額沿通道逐漸增加，液滴相重力較大，因此速度較慢，在相間拖曳力的作用下，氣相損失動能，氣相速度會逐漸下降，而液滴相速度逐漸增加；關(guān)閉液滴模型后，氣相僅存在與連續(xù)液相的相間拖曳力，該相間拖曳力較小，因此氣相速度沿通道逐漸增加并最終達(dá)到穩(wěn)定。圖4為本程序與參考程序在流場參數(shù)計算中的誤差分布，可以發(fā)現(xiàn)，氣相份額、氣相速度及液滴速度的計算相對誤差小于5%，液滴份額相對誤差小于10%，程序總體計算誤差較小。因此，多相場子通道程序可以合理計算液滴相運動行為。

圖2 液滴測試建模Fig.2 Model of droplet behavior test

圖3 液滴測試計算結(jié)果比較Fig.3 Comparison of calculation in droplet test

圖4 計算誤差分布Fig.4 Distribution of error in droplet test

2) 液滴尺寸對流場影響

多相場子通道程序中假定液滴為圓球形，液滴尺寸與液滴相界面面積呈反比。本測試通道長度為4 m，劃分為40個節(jié)點，邊界條件與算例1相同。本節(jié)采用臨界韋伯?dāng)?shù)(Wecrit)方法控制液體的尺寸，設(shè)定Wecrit分別為2.0、3.0、6.0、12.0，測試不同尺寸的液滴對流場分布的影響(圖5)。

由圖5a、b可知，Wecrit越小，則液滴直徑越小，液滴相界面面積越大，從而液滴的相間作用越強烈。由圖5c、d可知，Wecrit越小，相間面積越大，氣相與液滴相的相間拖曳力越強，因此計算穩(wěn)定后的氣相/液滴相的相對速度越小，對比不同Wecrit下的通道出口的液滴直徑與氣相/液滴相相對速度，可以發(fā)現(xiàn)，兩者基本呈線性分布，即液滴直徑越大，則流場中的氣相/液滴相相對速度越大。

圖5 不同液滴直徑計算結(jié)果Fig.5 Calculation of droplet test with different diameters

2.2 湍流交混測試

圖6 GE3×3實驗建模Fig.6 Model of GE3×3 experiment

選取GE3×3-1C棒束實驗[29]進(jìn)行湍流交混測試，實驗中棒束長度為1.828 8 m，單根棒直徑為0.014 5 m，棒中心距為0.014 7 m，棒束出口壓力為6.895 MPa，通道中流入過冷水，過冷度為264 K，各通道中的入口質(zhì)量流率均為1 342.67 kg/(m2·s)。3×3通道具有幾何對稱性，因此在程序中對棒束進(jìn)行1/8建模。如圖6所示，模型中包括3個通道，每個通道中設(shè)置有72個控制體，其中CH1為角通道，CH2為邊通道，CH3為中心通道，各通道入口處的質(zhì)量流率為1 342.67 kg/(m2·s)，觀察計算穩(wěn)定后各通道中質(zhì)量流率沿通道的分布情況。

CH1為角通道，水力學(xué)直徑最小，相同質(zhì)量流率下壁面阻力較大，因而通道壓降最大；CH3為中心通道，水力學(xué)直徑最大，相同質(zhì)量流率下壁面阻力最小，通道壓降最小。在通道間壓差作用下，CH1中流體通過間隙向外流出導(dǎo)致通道質(zhì)量流量降低，而CH3中有來自間隙的流體流入導(dǎo)致通道質(zhì)量流量增加。湍流交混體現(xiàn)了通道間的擴散效應(yīng)，用于降低不同通道中的質(zhì)量流率差異。圖7為湍流交混測試計算結(jié)果。由圖7可知，與關(guān)閉模型的結(jié)果相比，開啟模型后CH1～CH3通道間的質(zhì)量流率差值明顯減小，CH1出口質(zhì)量流率升高，而CH3出口質(zhì)量流率降低，計算結(jié)果與理論分析相符。此外，計算結(jié)果與實驗值符合良好，CH1～CH3通道出口流量的計算誤差分別為0.62%、1.70%、0.098%，因此，多相場子通道程序可以準(zhǔn)確計算通道間的湍流交混現(xiàn)象。

2.3 壁面?zhèn)鳠釡y試

Bennet-Hewitt實驗[30]是Bennet臨界后實驗，實驗壓力為6.89 MPa，入口處流入過冷度為12 K過冷水，入口質(zhì)量流量為393 kg/(m2·s)，流體沿通道向上流動并被加熱，壁面熱流密度為0.393 MW/m2。流體沿通道流動被加熱蒸發(fā)，壁面?zhèn)鳠崮Ｊ接膳R界前進(jìn)入臨界后，壁面溫度在臨界熱流密度處發(fā)生突跳。裝置建模如圖8所示，通道劃分為40個節(jié)點，入口和出口處分別為流動邊界與壓力邊界。觀察計算穩(wěn)定后，壁面溫度沿通道的分布情況。

圖7 湍流交混測試計算結(jié)果Fig.7 Result of turbulent mixing test

圖8 Bennet-Hewitt實驗建模Fig.8 Model of Bennet-Hewitt experiment

計算結(jié)果表明，由入口到出口，壁面?zhèn)鳠崮Ｊ接膳R界前過冷沸騰、飽和沸騰發(fā)展至臨界后彌散液滴模態(tài)沸騰、單相氣換熱。壁面溫度分布如圖9所示，在臨界前后傳熱模式轉(zhuǎn)換點，壁面熱流密度達(dá)到了臨界熱流密度(CHF)，壁面溫度發(fā)生突跳，實驗得到溫度起跳位置在0.426處，計算所得起跳位置在0.45處，程序計算出的溫度起跳位置較實驗中延后。溫度起跳點左側(cè)的臨界前傳熱模式，程序計算得到的壁面溫度較實驗值偏高，臨界前壁面溫度的計算誤差小于+15%，說明程序采用的臨界前壁面換熱模型在計算中低估了壁面換熱能力；臨界后的計算溫度與實驗數(shù)值符合良好，實驗值與計算值誤差在±5%以內(nèi)，說明由于具有液滴相換熱并采用了合理的臨界后壁面換熱模型，程序可以很好地模擬臨界后壁面換熱工況。綜上，程序可以很好地模擬跨臨界的壁面換熱工況。

圖9 壁面?zhèn)鳠釡y試計算結(jié)果Fig.9 Result of wall heat transfer test

2.4 再淹沒實驗驗證

圖10 FLECHT-TSP實驗建模Fig.10 Model of FLECHT-TSP

FLECHT-TSP實驗[31]為臨界后低淹沒速率實驗，實驗段長度為3.66 m，實驗工況壓力為275.79 kPa，加熱棒峰值功率為2.3 kW/m，其中加熱棒軸向功率分布沿頂部傾斜，最大峰均功率比為1.35。為驗證程序預(yù)測頂部傾斜功率分布的再淹沒工況的能力，對實驗裝置進(jìn)行建模，如圖10所示，包括20個軸向節(jié)點，流道內(nèi)設(shè)置有加熱棒，實驗中過冷水以0.02 m/s的速度由入口沿著流道向上流動，計算時間為500 s，觀測各監(jiān)測點的壁面溫度變化及驟冷時間。其中，實驗段nd1～nd5觀測點距離加熱棒底端距離分別為1.016、1.676 4、2.286、2.895 6、3.200 4 m。圖11為FLECHT-TSP計算結(jié)果。圖12為計算誤差分布。

圖11 FLECHT-TSP計算結(jié)果Fig.11 Result of FLECHT-TSP

圖12 計算誤差分布Fig.12 Distribution of error

由圖11可知，在計算開始后，由于加熱棒功率較高且通道中僅含換熱能力較差的單相氣，加熱棒表面溫度迅速升高并進(jìn)入臨界后狀態(tài)，隨后底部灌水浸潤了加熱棒表面，加熱棒由下至上，表面溫度依次急劇下降并進(jìn)入臨界前狀態(tài)，程序能較好地模擬出各觀測點的壁面溫度變化情況；但程序計算出的壁面溫度下降速度較實驗更快，這是由于程序中采用的壁面換熱模型在臨界前后的換熱能力差異過大，當(dāng)加熱棒表面被浸潤時，壁面?zhèn)鳠崮Ｊ接膳R界后進(jìn)入臨界前，換熱系數(shù)迅速增大，流體從壁面處帶走了大量熱量，導(dǎo)致壁面溫度迅速下降。由圖12可知，計算得到的各監(jiān)測點峰值溫度較實驗值稍高，各監(jiān)測點峰值溫度相對誤差均小于7%，表明程序計算出的峰值溫度是準(zhǔn)確的，且結(jié)果具有一定的保守度。綜上，程序計算結(jié)果與實驗值整體符合良好，程序能較好地模擬堆芯再淹沒工況。

3 結(jié)論

COSINE系列軟件包經(jīng)過多年的發(fā)展，已具備了壓水堆核電站的熱工水力分析等各種功能。為滿足先進(jìn)反應(yīng)堆“國和一號”的工程設(shè)計與安全分析要求，需開發(fā)具有更高計算精度的多相場子通道分析程序。本文介紹了多相場子通道程序的物理模型和數(shù)值方法，并通過程序?qū)φ辗椒ㄔu估了液滴模型計算的合理性；隨后，通過GE交混實驗，評估了程序?qū)ψ油ǖ澜换飕F(xiàn)象的計算能力；最后通過Bennet傳熱實驗及FLECHT實驗，驗證程序在計算臨界前后傳熱問題及再淹沒問題上的模擬能力。結(jié)果表明，程序在含液滴工況的計算是合理的，程序?qū)换鞂嶒灥挠嬎憬Y(jié)果與實驗數(shù)據(jù)符合良好；在臨界前后壁面換熱及再淹沒問題的計算中，程序可很好地模擬壁面溫度的分布及再淹沒工況下熱工水力學(xué)參數(shù)的變化。因此，多相場子通道程序?qū)Χ研緹峁にW(xué)現(xiàn)象的計算是合理且正確的，程序已具備了對堆芯子通道熱工水力現(xiàn)象的計算模擬能力。