李松蔚，Riccardo PURAGLIESI，楊 帆，余紅星，沈才芬

(1.中國核動力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213；2.Paul Scherrer Institute, Villigen 5232, Switzerland)

非能動安全系統(tǒng)有助于簡化核電站設(shè)計(jì)及潛在提高其經(jīng)濟(jì)性與安全性。然而，考慮到基于自然循環(huán)的弱驅(qū)動力，必須采用細(xì)致的設(shè)計(jì)和分析來確保系統(tǒng)執(zhí)行其預(yù)期功能。為此，2004—2008年，國際原子能機(jī)構(gòu)開展了關(guān)于“非能動系統(tǒng)的自然循環(huán)現(xiàn)象，建模和可靠性”的合作研究項(xiàng)目(CRP)，OECD/NEA的文件[1-2]中對事故期間非能動系統(tǒng)可能發(fā)生的大量熱工水力現(xiàn)象進(jìn)行了歸類。

大水池中水的行為是關(guān)注重點(diǎn)之一，包括溫度分層、自然/強(qiáng)迫對流及循環(huán)、蒸汽冷凝、上界面?zhèn)鳠醾髻|(zhì)等。實(shí)驗(yàn)重點(diǎn)研究溫度分層過程及其對傳熱能力的影響[3]，包括在瑞士保羅謝爾研究所(PSI)PANDA臺架上開展的以ESBWR(經(jīng)濟(jì)簡化沸水堆)為原型的旨在為系統(tǒng)程序及CFD程序驗(yàn)證提供基準(zhǔn)的ISP-42實(shí)驗(yàn)[4-5]、德國于利希研究中心(FZ Jülich)NOKO臺架上開展的冷凝器帶熱能力綜合實(shí)驗(yàn)[6-7]、德國亥姆霍茲德勒斯登羅森多夫研究中心(FZDR)TOPFLOW臺架上開展的可獲得二次側(cè)三維傳熱特性的傾斜管冷凝實(shí)驗(yàn)[8]、美國普渡大學(xué)PUMA臺架上開展的研究直接接觸冷凝及熱分層實(shí)驗(yàn)[9]等。此外德國FZDR還開展了側(cè)壁加熱圓柱水箱實(shí)驗(yàn)[10]，通過熱電偶及探針測得桶中水的局部溫度及空泡份額。

由于水池體積較大，局部區(qū)域的傳熱并不意味著池中溫度均勻，三維對流影響傳熱過程，導(dǎo)致溫度分層，根據(jù)需要可能需采用計(jì)算流體力學(xué)(CFD)方法。目前計(jì)算仍以系統(tǒng)程序?yàn)橹鳎捎肅FD的計(jì)算較少，Krepper[11]針對NOKO臺架實(shí)驗(yàn)進(jìn)行了單相CFD模擬，針對TOPFLOW的兩相實(shí)驗(yàn)進(jìn)行了兩相CFD模擬[12]，獲得兩維流場和溫場，但未與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較。而ISP-42實(shí)驗(yàn)雖適用于CFD計(jì)算驗(yàn)證，但基準(zhǔn)題計(jì)算提交均為系統(tǒng)程序計(jì)算結(jié)果，并未提交CFD模擬結(jié)果[5]。實(shí)驗(yàn)測量數(shù)據(jù)用于TRACE及GOTHIC程序模擬對比，研究結(jié)果發(fā)現(xiàn)，水池側(cè)的傳熱計(jì)算與實(shí)驗(yàn)有顯著差異[13]，這是因?yàn)門RACE是一維集總參數(shù)系統(tǒng)程序，且并非基于池式沸騰開發(fā)，而GOTHIC程序具有三維計(jì)算能力，但采用非常粗的計(jì)算網(wǎng)格，對局部沸騰特性模擬能力有限。因此，采用兩相CFD進(jìn)行模擬比較顯得尤為必要。本文基于STAR-CCM+軟件，采用兩相泡狀流及壁面沸騰模型，構(gòu)建過冷沸騰自然對流兩相模擬方法，結(jié)合側(cè)壁加熱圓柱水箱實(shí)驗(yàn)數(shù)據(jù)對過冷沸騰自然對流進(jìn)行研究，完成過冷沸騰自然對流兩相CFD模擬，應(yīng)用于ISP-42實(shí)驗(yàn)自然循環(huán)瞬態(tài)模擬研究。

1 過冷沸騰自然對流兩相CFD模擬

1.1 側(cè)壁加熱圓柱水箱實(shí)驗(yàn)

為構(gòu)建適用的過冷沸騰自然對流兩相模擬方法，首先采用德國FZDR實(shí)驗(yàn)室的圓柱形水箱過冷沸騰自然對流實(shí)驗(yàn)[10，14]進(jìn)行兩相CFD模擬研究與驗(yàn)證。實(shí)驗(yàn)在常壓下進(jìn)行，實(shí)驗(yàn)裝置為一圓柱形水箱，直徑0.25 m，高0.25 m，水位高0.21 m。壁面均勻加熱。實(shí)驗(yàn)中采用熱電偶及電導(dǎo)探針測量局部溫度及空泡份額，圖1為水箱及溫度、空泡份額測量位置示意圖，通過測量壁面溫度，可計(jì)算得到實(shí)際由壁面到水的隨時(shí)間變化的加熱功率曲線[14]。實(shí)驗(yàn)持續(xù)時(shí)間為1 500 s，數(shù)百s后，出現(xiàn)明顯的溫度分層，1 100 s后，空泡探針探測到氣泡產(chǎn)生，在整個(gè)實(shí)驗(yàn)過程中，水箱中心未監(jiān)測到氣泡產(chǎn)生。

圖1 溫度及空泡份額測量位置Fig.1 Locations of temperature and volume fraction measurements

1.2 幾何建模與網(wǎng)格劃分

模擬采用兩維軸對稱幾何建模。幾何建模高0.21 m，半徑0.125 m。采用結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分，并對加熱壁面(右側(cè))進(jìn)行網(wǎng)格加密。為避免網(wǎng)格相關(guān)性，進(jìn)行了多種網(wǎng)格劃分(表1)。圖2示出4種網(wǎng)格方案計(jì)算所得的速度，選取t=1 000 s、軸向高度0.1 m處進(jìn)行徑向方向速度對比。圖中x軸0.125 m表示靠近加熱壁面這一側(cè)，而0.0 m表示靠近圓柱水箱中心軸這一側(cè)。由于僅在近壁面處速度較明顯，因此圖中僅展示了x=0.110～0.125 m這一區(qū)域。由圖可知，A方案和B方案計(jì)算結(jié)果重合，因此選取網(wǎng)格劃分方案B可滿足網(wǎng)格敏感性要求。圖2中計(jì)算采用的網(wǎng)格劃分軸向節(jié)點(diǎn)數(shù)80，徑向節(jié)點(diǎn)數(shù)45。其中第1層網(wǎng)格為2.7×10-4m，y+為3，這是為了能成功捕捉相變過程及獲取更精細(xì)的速度場。

表1 網(wǎng)格劃分方案Table 1 Mesh setting case

圖2 不同網(wǎng)格方案近壁面液相速度計(jì)算結(jié)果Fig.2 Calculated liquid velocity in near wall region for different mesh setting cases

1.3 數(shù)值模型與不確定性分析

兩相CFD過冷沸騰自然循環(huán)是瞬態(tài)計(jì)算，模擬求解基于STAR-CCM+的歐拉-歐拉法[15]。壁面沸騰求解基于Kurul與Podowski[16]提出的壁面熱流分配(RPI)模型，在該模型中，來自壁面的給定熱通量被分為單相傳熱、基于成核點(diǎn)和淬火建立的氣泡生長蒸發(fā)傳熱、氣泡成長到一定大小離開壁面時(shí)由水代替原來氣泡占據(jù)位置所帶來的傳熱影響。模型計(jì)算了每部分的熱通量和壁溫。

為模擬主流區(qū)的兩相流，液相采用了可實(shí)現(xiàn)的k-ε雙層湍流模型，雙層全y+壁面處理。曳力采用標(biāo)準(zhǔn)阻力系數(shù)的Schiller-Naumann關(guān)系式，該模型適用于球形固體顆粒、液滴和小直徑(球形)氣泡。未來將考慮升力和湍流分散力。

采用Ranz和Marshall模型考慮相間傳熱傳質(zhì)。該模型認(rèn)為液體溫度低于飽和溫度則冷凝，液體溫度高于飽和溫度則蒸發(fā)。蒸汽氣泡直徑設(shè)置為單分散尺寸分布，它們在流動區(qū)域中的某個(gè)位置處的直徑，線性依賴于該位置處的液體溫度。這種簡化處理可能會影響相間界面面積，并因此影響冷凝速率。

由于初始水為靜止?fàn)顟B(tài)，而此后過冷沸騰產(chǎn)生氣泡會逸出或冷凝，因此幾何建模左側(cè)設(shè)置為軸向?qū)ΨQ邊界條件，右側(cè)設(shè)置為加熱壁面，底部為壁面無滑移。水箱上表面設(shè)置為氣相滲透，該邊界條件是使用無滑移壁面及氣相滲透壁面對氣體逸出表面進(jìn)行建模，可與自由出口或壓力出口邊界條件共同使用。STAR-CCM+自動計(jì)算通過氣相滲透邊界的氣體質(zhì)量流量(逸氣率)。采用該邊界條件，允許氣相自由逸出，但不允許液相逸出，考慮到此時(shí)處于過冷沸騰狀態(tài)，水箱中心液體溫度仍低于飽和溫度，逸出的蒸汽較少，不足以降低自由液面，計(jì)算未考慮自由液面高度降低。

自然對流由流體密度不同所驅(qū)動，液相須考慮隨溫度變化導(dǎo)致的密度變化以及浮力作用，因此水物性(密度、比定壓熱容、動力黏度、熱傳導(dǎo)率)采用隨溫度變化的函數(shù)，氣相密度采用隨當(dāng)?shù)貕毫ψ兓暮瘮?shù)，浮力采用Boussinesq浮力模型。

由于在低壓情況下由壓降引起的飽和溫度下降較明顯，水箱內(nèi)不同高度壓力不同，因此飽和溫度被設(shè)置為與壓力相關(guān)的函數(shù)。參考壓力設(shè)置為0.1 MPa。

計(jì)算主要關(guān)注溫度變化，針對與此相關(guān)的過冷沸騰模型中的氣泡脫離直徑進(jìn)行了敏感性分析，計(jì)算對比了模型默認(rèn)氣泡脫離直徑、2倍模型默認(rèn)氣泡脫離直徑及1/2倍模型默認(rèn)氣泡脫離直徑，溫度計(jì)算結(jié)果無明顯改變，后繼計(jì)算模型參數(shù)均采用默認(rèn)值。

從模型研發(fā)上來說，采用CFD計(jì)算自然循環(huán)需考慮對高雷諾數(shù)流動可能超出Boussinesq近似的有效范圍；雷諾應(yīng)力和雷諾通量均基于各向同性開發(fā)，而自然循環(huán)重力起重要作用，需考慮各向異性[3]。計(jì)算未出現(xiàn)高雷諾數(shù)流動情形，同時(shí)已考慮重力影響，模型適用于當(dāng)前研究情形。

1.4 計(jì)算結(jié)果與分析

1) 兩相CFD模擬結(jié)果

圖3示出了兩相CFD計(jì)算不同時(shí)刻的水溫云圖與液相速度場。兩相CFD模擬可成功再現(xiàn)實(shí)驗(yàn)中由于氣泡產(chǎn)生流動加速所導(dǎo)致的第1次溫度階躍(800～1 200 s)，以及由于水箱整體水溫達(dá)到飽和出現(xiàn)的第2次階躍(約1 300 s)。模擬展示出隨時(shí)間的推移，沸騰在加熱面水箱上表面出現(xiàn)，由于水箱上部蒸汽產(chǎn)生的影響，自然循環(huán)顯著，溫度分布相同；而下部由于還是單相，仍保持穩(wěn)定的熱分層狀態(tài)。隨時(shí)間的推移，溫度分布相同的區(qū)域和熱分層的區(qū)域水平界線將逐漸向下推移。兩相CFD模擬展示了沸騰出現(xiàn)后水箱中顯著的自然循環(huán)，如采用單相CFD模擬，由于不考慮相變，無法正確模擬過冷沸騰出現(xiàn)后的溫度特性，水箱中不會出現(xiàn)類似自然循環(huán)結(jié)果，僅在加熱面水箱上表面流動較顯著。

2) 不同CFD軟件計(jì)算對比

將本文STAR-CCM+軟件計(jì)算結(jié)果與文獻(xiàn)[14]中采用CFX軟件計(jì)算結(jié)果進(jìn)行對比。圖4示出了CFX軟件和STAR-CCM+軟件典型位置計(jì)算溫度與實(shí)驗(yàn)測量溫度隨時(shí)間的變化對比。從圖4可見，CFX對水箱上表面的計(jì)算更準(zhǔn)確，但對水箱下部的計(jì)算誤差較大，出現(xiàn)溫度第1次階躍后，CFX的預(yù)測值較實(shí)驗(yàn)值偏低。STAR-CCM+對水箱上表面溫度計(jì)算過高的原因是在STAR-CCM+中沸騰模型與壁面滑移無法同時(shí)使用，因此計(jì)算的水箱上表面速度非常小，這使得溫度交混很弱，導(dǎo)致溫度被高估。而在實(shí)驗(yàn)中，由于水箱上表面為自由表面，設(shè)置為可滑移邊界條件是更符合實(shí)際的。該問題可通過未來STAR-CCM+中在邊界構(gòu)造特殊的速度函數(shù)解決。兩種軟件均具有自然對流過冷沸騰兩相模擬能力，STAR-CCM+除水箱上表面溫度計(jì)算過高外，其余部分溫度模擬均比較準(zhǔn)確，因此可采用STAR-CCM+進(jìn)行下一步研究。

圖3 兩相CFD計(jì)算水溫云圖與液相速度場Fig.3 Two-phase CFD calculation temperature and liquid velocity fields

圖4 不同位置計(jì)算及測量溫度隨時(shí)間的變化Fig.4 Calculation and experiment temperatures of different locations vs. time

3) 實(shí)驗(yàn)驗(yàn)證

用STAR-CCM+對實(shí)驗(yàn)進(jìn)行模擬驗(yàn)證。圖5示出了水箱近壁面測量位置計(jì)算的溫度與實(shí)驗(yàn)測量溫度隨時(shí)間變化的對比。從圖5可知，計(jì)算和實(shí)驗(yàn)結(jié)果符合良好，合理再現(xiàn)了實(shí)驗(yàn)中的兩次溫度階躍。

圖6示出了計(jì)算空泡份額與實(shí)驗(yàn)空泡份額隨時(shí)間變化的比較結(jié)果?？张莘蓊~上升趨勢及出現(xiàn)位置、出現(xiàn)時(shí)間均與實(shí)驗(yàn)符合較好。剛出現(xiàn)兩相時(shí)，僅在接近上表面的加熱壁面附近出現(xiàn)氣泡，隨時(shí)間的推移，蒸汽存在區(qū)域由壁面向中心遷移。由于實(shí)驗(yàn)中空泡探針探頭直徑為1 mm，無法探測1.5 mm直徑以下的氣泡，因此計(jì)算空泡份額較實(shí)驗(yàn)空泡份額高是合理的。

圖5 各測量位置計(jì)算水溫與實(shí)驗(yàn)水溫的對比Fig.5 Comparison of calculation and experiment water temperatures at different measurement locations

圖6 各測量位置計(jì)算空泡份額與實(shí)驗(yàn)空泡份額的對比Fig.6 Comparison of calculation and experiment void fractions at different measurement locations

通過對圓柱形水箱兩相自然對流實(shí)驗(yàn)進(jìn)行兩相CFD數(shù)值模擬計(jì)算，完成了網(wǎng)格敏感性分析、時(shí)間步長選擇、模型選擇，并利用實(shí)驗(yàn)數(shù)據(jù)對數(shù)值計(jì)算結(jié)果進(jìn)行了驗(yàn)證和評價(jià)，數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)在溫度、空泡出現(xiàn)時(shí)間、出現(xiàn)位置上均符合較好，該套計(jì)算方法可用于非能動安全殼冷卻系統(tǒng)(PCCS)計(jì)算。

2 ISP-42實(shí)驗(yàn)?zāi)M應(yīng)用

在PANDA臺架上開展的ISP-42實(shí)驗(yàn)原型為ESBWR非能動安全殼冷卻系統(tǒng)。實(shí)驗(yàn)包括6個(gè)階段(A～F)，代表一系列操作模式或過程，每個(gè)階段實(shí)際上是1個(gè)獨(dú)立的實(shí)驗(yàn)，具有自己的初始條件和邊界條件。模擬所關(guān)注的是B階段，由于蒸汽進(jìn)入管道而導(dǎo)致的PCCS水箱傳熱沸騰。裝置(圖7)為高5 m的水箱，中間有20根通蒸汽的管道，水箱通過自然循環(huán)帶走管道的熱量。兩相模擬主要關(guān)注3 000～3 500 s時(shí)間段。在這個(gè)時(shí)間段，流體溫度上升，壁面已有氣泡產(chǎn)生，但主流仍無氣泡存在[5]。

圖7 實(shí)驗(yàn)裝置流體溫度測量位置Fig.7 Liquid temperature measurement location of experiment facility

計(jì)算采用二維幾何建模，建模選擇中間剖面。圖8示出了實(shí)驗(yàn)中不同監(jiān)測點(diǎn)水溫隨時(shí)間的變化與計(jì)算結(jié)果的對比，計(jì)算與實(shí)驗(yàn)結(jié)果符合較好，可較好跟隨流體溫度變化。再次證明了該兩相模擬方法的適用性。蒸汽的產(chǎn)生加速了水的流動，單相時(shí)自然對流不強(qiáng)，而兩相的出現(xiàn)加強(qiáng)了自然對流。

圖8 測量點(diǎn)計(jì)算與實(shí)驗(yàn)水溫隨時(shí)間的變化Fig.8 Calculation and experiment water temperatures at measured point vs. time

3 結(jié)論

通過模擬側(cè)壁加熱圓柱形水箱過冷沸騰自然循環(huán)實(shí)驗(yàn)，成功構(gòu)建了兩相CFD模擬方法，實(shí)現(xiàn)了過冷沸騰自然對流瞬態(tài)模擬。

1) 兩相CFD模擬可獲得與實(shí)驗(yàn)較一致的溫度分布、氣泡產(chǎn)生時(shí)間與產(chǎn)生位置?？烧_模擬沸騰出現(xiàn)后的溫度階躍，以及水箱內(nèi)顯著的自然循環(huán)，單相CFD無法模擬該現(xiàn)象。兩相CFD計(jì)算結(jié)果能合理解釋實(shí)驗(yàn)現(xiàn)象的發(fā)生，證實(shí)了在過冷沸騰自然對流上進(jìn)行數(shù)值計(jì)算的合理性與準(zhǔn)確性。

2) 采用相同的兩相CFD模擬對PANDA臺架上ISP-42實(shí)驗(yàn)進(jìn)行模擬(基于歐洲ESBWR PCCS原型)，可獲得與實(shí)驗(yàn)較一致的溫度分布，再次驗(yàn)證了該模擬方法的正確性。

該兩相CFD模擬方法可用于構(gòu)建自然對流傳熱特性模型，為非能動安全殼冷卻系統(tǒng)及非能動余熱排出系統(tǒng)設(shè)計(jì)提供新的設(shè)計(jì)手段。