MRANS方案的反應(yīng)堆壓力容器CFD仿真

康慧倫，田兆斐，胡培政，陳廣亮

(核安全與仿真技術(shù)國防重點學(xué)科實驗室(哈爾濱工程大學(xué)) 哈爾濱 150001)

作為堆芯的入口邊界，壓力容器下腔室出口的冷卻劑流量、溫度及硼的分布對反應(yīng)堆的安全性和經(jīng)濟性有著重要影響。因此，對壓力容器中冷卻劑流動特性進行研究成為反應(yīng)堆安全分析的重要組成部分。隨著CFD技術(shù)的快速發(fā)展，利用CFD分析方法對壓力容器進行數(shù)值模擬，逐漸成為與實驗同樣重要的研究手段。中國核動力研究院的眭曦等[1]基于下腔室實驗對CFD計算中不同交混模型進行對比，認為組分輸運模型可以較好的描述由于濃度梯度導(dǎo)致的硼擴散現(xiàn)象。郭超等[2]基于CFX程序進行下腔室全尺寸計算，并對下腔室流場及壓力場進行研究，驗證下腔室流量分配設(shè)計合理性。Cho等[3]基于ROCOM實驗對CFD程序進行驗證，驗證結(jié)果表明CFD程序在適當(dāng)?shù)木W(wǎng)格處理情況下可以合理的預(yù)測壓力容器內(nèi)的硼稀釋現(xiàn)象。

大部分關(guān)于壓力容器的實驗的裝置按比例對真實壓力容器進行縮放，但其復(fù)雜的下腔室結(jié)構(gòu)和依然較大的模擬域?qū)е翪FD模擬所需的網(wǎng)格數(shù)量龐大，對計算使用的硬件提出了較高的要求。同時，冷卻劑在壓力容器下降段與下腔室結(jié)構(gòu)差異很大(下降段為環(huán)段結(jié)構(gòu)，下腔室為球形且內(nèi)含孔板結(jié)構(gòu))，使用同一湍流模型可能無法保證精度與效率的最佳平衡。為解決這些問題，本文將分段技術(shù)和MRANS方案應(yīng)用在壓力容器CFD計算中。Navarro等[4]采用分段模擬技術(shù)對Karoutas等[5]的5×5棒束實驗進行模擬，以減少每次計算的計算負荷；Li等[6]對分段模擬技術(shù)在棒束通道仿真中的應(yīng)用特性進行了研究；Chen等[7-8]針對帶交混翼棒束通道不同區(qū)域流體流動狀態(tài)不同的特點，基于分段模擬技術(shù)提出了多湍流模型(Multiple RANS models, MRANS)方案，在保證計算精度的同時提高計算效率，并提出了DP-CFD方案以解決堆芯大流域CFD仿真問題。值得注意的是，絕大部分的分段CFD分析技術(shù)和MRANS方案的討論都是針對棒束幾何模型，在其他幾何模型中的適用性仍然需要討論。

ROCOM實驗[9-10]被廣泛應(yīng)用在CFD計算方法的驗證[11]中，許多學(xué)者基于ROCOM實驗裝置進行了CFD分析，其中包括網(wǎng)格敏感性、湍流模型和邊界條件效應(yīng)等等[12-14]。本文以ROCOM實驗為對象，針對分段技術(shù)與MRANS方案展開研究，提出合理的分段計算與MRANS方案，減小計算負荷，提高計算效率，并保證CFD計算的準確性。

1 仿真對象介紹

本文以歐洲FLOMIX-R基準實驗中[9-10]的ROCOM實驗裝置為仿真對象，對穩(wěn)態(tài)運行的壓力容器的冷卻劑交混現(xiàn)象進行研究。ROCOM實驗裝置以KONVOI三代壓水堆為原型，按1∶5的比例建造，反應(yīng)堆壓力容器(RPV)為其主體設(shè)備，包括4個冷卻劑入口、出口、下降段(downcomer)及下腔室(lowerplenum)，下腔室布置有空心孔眼滾筒，堆芯入口布置支撐板。裝置具有4個完整的回路，并配有泵、閥門及蒸發(fā)器等設(shè)備，通過對泵的流量控制可以實現(xiàn)不同工況的實驗。在進行混合實驗時，實驗裝置在室溫和環(huán)境壓力下運行，示蹤劑氯化鈉溶液在入口上游注入主冷卻劑系統(tǒng)，并通過裝置進行混合，確保示蹤劑均勻分布在壓力容器入口橫截面處，壓力容器上布置有4處測點(容器入口、下降段上部、下降段下部及堆芯入口)對示蹤劑進行跟蹤，具體參數(shù)如圖1所示。

圖1 壓力容器及測點分布圖[5-6]

2 分段CFD技術(shù)研究

2.1 數(shù)值建模

2.1.1 幾何模型和網(wǎng)格劃分

根據(jù)Boumaza等[14]的研究，環(huán)路對發(fā)生在RPV的交混現(xiàn)象的影響很小。因此，在不影響計算準確性的前提下減少計算所需網(wǎng)格，本文使用300 mm的圓管代替環(huán)路作為RPV入口。堆芯采用高度為800 mm的圓柱代替，圓柱以堆芯入口支撐板上直徑30 mm的冷卻劑流量分配孔為基準進行延長。

圖2展示了本文所使用的網(wǎng)格模型。如圖2所示，本文采用多面體網(wǎng)格對RPV模型入口環(huán)腔和下腔室進行網(wǎng)格劃分，這是因為多面體網(wǎng)格與四面體網(wǎng)格相比，在達到相同精度的前提下具有更好的收斂性[6]。多面體網(wǎng)格的生成方式采用將四面體在FLUENT軟件中轉(zhuǎn)化的方式生成。下降段則采用六面體網(wǎng)格。表1展示了在網(wǎng)格敏感性分析中使用的不同劃分方案的網(wǎng)格量。圖3顯示了不同網(wǎng)格方案中堆芯入口處標量最大和平均值隨時間的變化趨勢。由于網(wǎng)格劃分占據(jù)了前處理中的絕大部分時間，同時無止境的增加網(wǎng)格量會增加大量的計算時間。因此，盡管未達到所謂的網(wǎng)格無關(guān)解，考慮精度與計算時間的平衡的情況下，本文選擇了方案2作為本文的計算網(wǎng)格劃分方案。

圖2 網(wǎng)格模型

表1 網(wǎng)格劃分方案

圖3 堆芯入口處標量最大值和平均值計算結(jié)果

2.1.2 邊界條件和湍流模型

本文以ROCOM_STAT_02實驗工況為仿真對象，其是針對低流速下的穩(wěn)態(tài)混合進行的實驗，實驗過程中沒有發(fā)生熱傳導(dǎo)且操作環(huán)境為室溫和環(huán)境壓力，因此選擇25 ℃，1.01×105Pa的水作為計算工質(zhì)，采用0.73 m/s(25%正常流速，即各環(huán)路的質(zhì)量流率為46.25 m3/h)的均勻速度分布作為入口邊界條件，出口采用壓力出口邊界條件。作為硼替代物的混合標量采用均勻分布的方式由入口1中注入，注入時間共35 s。本文采用的時間步長與實驗中測點的時間監(jiān)測間隔相同，均為0.05 s(頻率20 Hz)。本文的分段計算方案研究中,均使用RSM湍流模型。

2.2 分段方案研究

RPV根據(jù)流域特點可以分為：下降段(DC)和下腔室(LP)。在建模和網(wǎng)格劃分過程中，將DC和LP從整體RPV中拆分，得到代表各區(qū)域的獨立模型。本文利用FLUENT軟件帶有的界面數(shù)據(jù)讀寫功能，通過編寫的Python控制程序控制DC和LP之間的數(shù)據(jù)傳輸，完成分段計算的自動計算流程。圖4展示了DC和LP的數(shù)據(jù)傳遞過程。

圖4 分段計算的數(shù)據(jù)傳輸方案

在文獻[4,6-7,15]的分段計算中，將計算區(qū)域沿流動方向分為若干段，并將前一段的出口的速度場結(jié)果作為下一段的入口邊界條件，如圖4中黑線的傳輸方案。這種計算方式相當(dāng)于把流場下游對上游的影響切斷，意味著假設(shè)了下游的幾何與流動狀態(tài)對上游沒有影響。這種數(shù)據(jù)傳輸方式可以保證在沒有回流模型中速度場的準確傳播，但由于缺乏下游模型，上游的出口壓力場將無法準確得到。綜上所述，由于缺少數(shù)據(jù)通信，會導(dǎo)致上游出口面與下游入口面的壓力分布不一致，無疑會引入一定的誤差。

針對上、下游交界面壓力不一致的問題，本文提出了一種基于區(qū)域重疊的雙向數(shù)據(jù)傳遞分段計算方案。區(qū)域重疊方案如圖5所示，DC域向下拉伸了150 mm的圓環(huán)柱與LP域重疊。在計算過程中，DC將LP的入口所在面的速度和其他標量數(shù)據(jù)作為入口邊界條件傳遞至LP，而LP將DC的出口所在面的壓力分布數(shù)據(jù)作為出口邊界條件傳遞至DC，其與文獻[4,6-7]的單向數(shù)據(jù)傳輸?shù)牟煌幵谟跀?shù)據(jù)傳輸包括了圖4中紅線的傳輸內(nèi)容(壓力場)。圖6展示了不同的分段方案交界面處的壓力對比圖。

圖5 區(qū)域劃分方案

圖6 壓力分布對比

如圖6所示，相較于pressure-outlet出口邊界的均勻壓力場，outflow出口邊界的壓力分布在趨勢上與整體計算相似，這與Li等[6]所得到結(jié)果相同。然而，由于單向數(shù)據(jù)傳輸?shù)淖陨砣毕?，使用outflow邊界在交界面處的壓差與整體計算依然有一定差距(整體計算中界面壓力最高點與壓力最低點的壓差為100 Pa，而outflow計算為40 Pa)。而雙向數(shù)據(jù)傳遞的分段方案考慮了上下游的相互反饋，從圖6可以看出，雙向數(shù)據(jù)傳遞的分段方案與整體計算最接近。

圖7、8分別為不同的分段方案與整體計算的結(jié)果對比。圖7為DC下降段下部的3個測點的速度。考慮到計算幾何具有對稱性,3個測點選取在同一個象限內(nèi)。圖8為LP中16號通道與42號通道堆芯入口平面的速度。選取這兩個通道是因為兩個通道分別在下腔室空心孔眼滾筒的內(nèi)外兩側(cè)，具有一定的代表性。從圖中可以看出，采用雙向數(shù)據(jù)傳遞的分段方案與整體計算結(jié)果對比具有較好的一致性。

圖7 DC域中3個測點的速度

圖8 堆芯入口通道的速度

3 多湍流模型方案研究

3.1 MRANS方案設(shè)計

受到計算資源的限制，DNS和LES模型在工業(yè)CFD仿真領(lǐng)域使用較少，通常采用基于雷諾時均方法的RANS模型進行計算。雷諾時均的N-S(Navier-Stokes)方程可以表示如下：

(1)

(2)

標準k-ε模型(SKE)、realizablek-ε模型(RKE)、renormalization groupk-ε模型(RNG)及shear stress transportk-ε模型(SST)等2方程模型基于Boussinesq假設(shè)和湍流黏度的概念，通過下式求解雷諾應(yīng)力。

(3)

式中：k為湍流動能,μt為湍流黏度。湍流黏度是湍流動能和湍流耗散率或比耗散率的函數(shù)，因此求解兩方程模型的雷諾應(yīng)力需要建立湍流動能輸運方程和湍流耗散率或比耗散率輸運方程。

雷諾應(yīng)力模型(RSM)摒棄了各向同性的湍流黏度概念，直接構(gòu)造了關(guān)于雷諾應(yīng)力的輸運模型，一般如下式所示。

(4)

式中：等號左側(cè)第1項為應(yīng)力時間項,第2項為應(yīng)力對流項,右側(cè)Dij為湍流擴散項,φij為壓力應(yīng)變項,Gij為應(yīng)力產(chǎn)生項,εij為耗散項。

經(jīng)過?；蟮睦字Z應(yīng)力模型中的雷諾應(yīng)力為二階張量，需要求解6個雷諾應(yīng)力分量方程和1個湍流耗散輸運方程。RSM模型在求解雷諾應(yīng)力時需要建立更多的方程，導(dǎo)致其計算效率低于2方程模型，但由于RSM模型直接求解雷諾應(yīng)力張量，可以更加準確的反應(yīng)雷諾應(yīng)力的各向異性，一般來說在復(fù)雜湍流應(yīng)用中RSM模型的適用性和精確性優(yōu)于兩方程模型[15]。

由于在RPV整個流域中，平均雷諾數(shù)達到了104量級，所以本文的多湍流模型方案選取SKE、RKE以及RSM這3個高雷諾數(shù)模型進行研究。上述3種模型有相同的壁面處理方式，可使用同樣類型與密度的網(wǎng)格對多湍流模型混合計算精度與效率進行對比分析。MRANS計算時需要在區(qū)域邊界面上對求解動量方程與湍流方程的重要參數(shù)數(shù)據(jù)進行傳遞，表2顯示了分段界面上各模型計算所需的參數(shù)。

表2 湍流模型參數(shù)

圖9、10為應(yīng)用不同的湍流模型的計算結(jié)果與實驗值[11]的對比，分別對應(yīng)著DC不同高度(upper和lower)測點的混合標量值。

圖9 下降段上部測點標量方位角數(shù)據(jù)

圖10 下降段下部測點標量方位角數(shù)據(jù)

如圖9、10所示，在下降段上部2方程模型與RSM模型的計算結(jié)果差異很小且均與實驗值符合較好。在下降段下部，2方程模型預(yù)測的標量場空間分布的坡度比實驗值更加陡峭，而RSM明顯與實驗值的分布更加接近?？梢钥闯?，2方程模型明顯的低估了下降段湍流交混的能力，這可能是將各向同性的2方程湍流模型應(yīng)用在具有高的各向異性的結(jié)構(gòu)上所導(dǎo)致的。同時，這種交混能力的低估在很多的有關(guān)FLOMIX的結(jié)果中被發(fā)現(xiàn)[10]。

根據(jù)下降段的計算結(jié)果對比，本文基于SKE，RKE，RSM 3個模型進行互相匹配，設(shè)計多湍流模型耦合計算方案。因為多湍流模型計算方案的意義在于解決模擬大尺度部件時計算時間與計算精度難以平衡的問題，本文將遵循保證計算精度同時減小計算所需的資源的原則進行MRANS方案的設(shè)計。圖11為本文設(shè)計的MRANS計算方案，其中垂線連接的模型屬于單一模型計算方案，斜線連接屬于MRANS計算方案。通過各計算方案對整個流域的求解計算，可以從計算精度與效率兩個方面綜合評價MRANS方案與傳統(tǒng)方案的區(qū)別。從表2中可以看出，7方程模型包含有2方程模型所需要的全部參數(shù)，因此，在將7方程模型應(yīng)用在2方程模型的上游時，無需補充參數(shù)。

3.2 MRANS方案驗證

根據(jù)圖11的MRANS計算方案，分別對5種SKE,RKE,RSM的單一/混合模型進行計算。作為下腔室的入口，下降段出口處的流場及標量場結(jié)果的準確性會對下腔室流場的計算造成極大的影響。圖12分別展示了MRANS方案和單一湍流模型方案在下降段下部測點預(yù)測的速度和標量。圖中黑色實線為實驗值，虛線為實驗值±30%范圍。圖12(a)中速度為z軸方向速度。如圖12所示，由于MRANS模型在下降段使用RSM模型進行計算，因此在下降段下部測點的速度分布與RSM模型計算結(jié)果高度一致。SKE模型和RKE模型的計算結(jié)果也表現(xiàn)出一致性。在速度分布的計算中，RSM模型和MRANS模型在大部分區(qū)域與2方程模型精度相差不大，但在峰/谷值處與實驗值更一致。

圖11 MRANS計算方案

考慮到下降段標量場的影響范圍，圖12(b)將標量未擴散的區(qū)域的數(shù)據(jù)予以剔除。相較于7方程模型，單一的兩方程模型低估了低濃度區(qū)域的混合標量濃度，并且高估了高濃度區(qū)域的混和標量濃度。這與圖10的計算結(jié)果相符合。RSM的混合標量預(yù)測結(jié)果與實驗值符合最好，而MRANS模型也僅在極少的位置與RSM模型有差異。無論是速度場還是標量場，MRANS模型都展現(xiàn)出了比SKE或RKE更優(yōu)秀的精度。

圖12 下降段下部測點的速度與混合標量

圖13、14展示了不同RANS方案的RPV堆芯入口第16號和第25號通道的標量濃度的計算結(jié)果。在靠近下降段的第16號通道，單一的SKE和RKE模型預(yù)測的混合標量濃度上升非常迅速，在10 s內(nèi)達到了充分混合的平穩(wěn)階段(在實驗中上升過程約為15 s)，且過高的預(yù)測了平穩(wěn)階段的混合標量濃度。相較于2方程模型，RSM模型在上升階段和平穩(wěn)階段的預(yù)測結(jié)果均與實驗符合較好，兩種MRANS方案中RSM-SKE模型在16通道均與RSM模型符合較好，且能較準確的描述標量隨時間的變化過程。單一的SKE和RKE模型計算結(jié)果相類似，在上升階段和平穩(wěn)階段預(yù)測結(jié)果與實驗值的平均偏差分別為32%和17%，MRANS模型為11%和9%，RSM模型為5%和7%。而在相對更靠近中心的25號通道，單一的兩方程模型預(yù)測的上升過程與實驗符合較好，但仍低估了平穩(wěn)階段的標量濃度。所有的模型在上升階段的預(yù)測與實驗的平均偏差都在10%左右。在平穩(wěn)階段，單一的SKE和RKE模型的平均偏差為19%，MRANS模型為11%，RSM模型為5%。兩種MRANS方案在上升階段和平穩(wěn)階段的前半段預(yù)測符合較好，但在平穩(wěn)階段的后半段的混合標量濃度預(yù)測較RSM仍有較大偏差。之所以產(chǎn)生這種情況是因為，16通道為邊緣通道，靠近下降段，流體從下降段進入下腔室后迅速進入通道，因此MRANS方案對該區(qū)域的準確性的提升較大。而25通道為中心通道，流體進入下腔室后需要流經(jīng)分配孔板等幾何結(jié)構(gòu)后才能進入通道，因此MRANS方案對該區(qū)域的提升較小。

圖13 堆芯入口第16號通道的標量隨時間變化

圖14 堆芯入口第25通道標量隨時間變化

通過上述分析可以看出，基于RSM模型的MRANS方案可以得到與單一RSM模型符合較好的計算結(jié)果。由于不是全域使用更耗時的RSM模型，采用MRANS方案減少了大量的計算時間。表3為本文所設(shè)計的MRANS方案與單一的RSM模型在Inter(R) Xeon(R) CPU E3-1280 v5，3.7 GHz，64 G的單機8核處理器條件下的計算效率。如表3所示，MRANS方案相較于單一的RSM模型更加高效，在本文的網(wǎng)格規(guī)模下可以縮短約30%的CPU時間，且隨著網(wǎng)格數(shù)量增加或LP網(wǎng)格數(shù)量在整體RPV模型的比重增加，可以增大這種MRANS方案的優(yōu)勢。

表3 MRANS計算方案與RSM方案計算效率

4 結(jié) 論

1)本文將棒束CFD分析中提出的分段方案運用到壓力容器CFD分析中，設(shè)計了基于區(qū)域重疊方法的雙向數(shù)據(jù)傳遞分段計算方案。

2)分析了不同的湍流模型在各區(qū)域的計算結(jié)果，為分區(qū)計算選擇合適的湍流模型。

3)通過運用MRANS方案和所設(shè)計的雙向數(shù)據(jù)傳遞分段計算方案構(gòu)建了SKE-RSM與RKE-RSM兩種MRANS方案，相比于傳統(tǒng)的單一RSM湍流模型計算方案，本文提出的MRANS方案具有更高的效率。研究成果對提高RPV CFD仿真效率提供了一種新思路和途徑。