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

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

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

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

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

1 仿真對象介紹

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

圖1 壓力容器及測點(diǎn)分布圖[5-6]

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

2.1 數(shù)值建模

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

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

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

圖2 網(wǎng)格模型

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

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

2.1.2 邊界條件和湍流模型

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

2.2 分段方案研究

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

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

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

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

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

圖6 壓力分布對比

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

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

圖7 DC域中3個(gè)測點(diǎn)的速度

圖8 堆芯入口通道的速度

3 多湍流模型方案研究

3.1 MRANS方案設(shè)計(jì)

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

(1)

(2)

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

(3)

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

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

(4)

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

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

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

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

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

圖9 下降段上部測點(diǎn)標(biāo)量方位角數(shù)據(jù)

圖10 下降段下部測點(diǎn)標(biāo)量方位角數(shù)據(jù)

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

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

3.2 MRANS方案驗(yàn)證

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

圖11 MRANS計(jì)算方案

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

圖12 下降段下部測點(diǎn)的速度與混合標(biāo)量

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

圖13 堆芯入口第16號通道的標(biāo)量隨時(shí)間變化

圖14 堆芯入口第25通道標(biāo)量隨時(shí)間變化

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

表3 MRANS計(jì)算方案與RSM方案計(jì)算效率

4 結(jié) 論

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

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

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