過冷沸騰自然對流兩相CFD模擬及應(yīng)用

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

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

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

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

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

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

1.1 側(cè)壁加熱圓柱水箱實驗

為構(gòu)建適用的過冷沸騰自然對流兩相模擬方法，首先采用德國FZDR實驗室的圓柱形水箱過冷沸騰自然對流實驗[10，14]進行兩相CFD模擬研究與驗證。實驗在常壓下進行，實驗裝置為一圓柱形水箱，直徑0.25 m，高0.25 m，水位高0.21 m。壁面均勻加熱。實驗中采用熱電偶及電導(dǎo)探針測量局部溫度及空泡份額，圖1為水箱及溫度、空泡份額測量位置示意圖，通過測量壁面溫度，可計算得到實際由壁面到水的隨時間變化的加熱功率曲線[14]。實驗持續(xù)時間為1 500 s，數(shù)百s后，出現(xiàn)明顯的溫度分層，1 100 s后，空泡探針探測到氣泡產(chǎ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)格進行劃分，并對加熱壁面(右側(cè))進行網(wǎng)格加密。為避免網(wǎng)格相關(guān)性，進行了多種網(wǎng)格劃分(表1)。圖2示出4種網(wǎng)格方案計算所得的速度，選取t=1 000 s、軸向高度0.1 m處進行徑向方向速度對比。圖中x軸0.125 m表示靠近加熱壁面這一側(cè)，而0.0 m表示靠近圓柱水箱中心軸這一側(cè)。由于僅在近壁面處速度較明顯，因此圖中僅展示了x=0.110～0.125 m這一區(qū)域。由圖可知，A方案和B方案計算結(jié)果重合，因此選取網(wǎng)格劃分方案B可滿足網(wǎng)格敏感性要求。圖2中計算采用的網(wǎng)格劃分軸向節(jié)點數(shù)80，徑向節(jié)點數(shù)45。其中第1層網(wǎng)格為2.7×10-4m，y+為3，這是為了能成功捕捉相變過程及獲取更精細的速度場。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2) 不同CFD軟件計算對比

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

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

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

3) 實驗驗證

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

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

圖5 各測量位置計算水溫與實驗水溫的對比Fig.5 Comparison of calculation and experiment water temperatures at different measurement locations

圖6 各測量位置計算空泡份額與實驗空泡份額的對比Fig.6 Comparison of calculation and experiment void fractions at different measurement locations

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

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

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

圖7 實驗裝置流體溫度測量位置Fig.7 Liquid temperature measurement location of experiment facility

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

圖8 測量點計算與實驗水溫隨時間的變化Fig.8 Calculation and experiment water temperatures at measured point vs. time

3 結(jié)論

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

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

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

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