李小暢，鄧 堅，楊小磊，郭涌輝，陳廣亮，田瑞峰,*

(1.哈爾濱工程大學(xué) 航天與建筑工程學(xué)院，黑龍江 哈爾濱 150001；2.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點實驗室，四川 成都 610041；3.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院，黑龍江 哈爾濱 150001)

反應(yīng)堆三維精細(xì)化數(shù)值建模與仿真技術(shù)是克服傳統(tǒng)研究方法各種不足及大幅降低研發(fā)成本和周期的有效途徑之一，也是實現(xiàn)核能技術(shù)多學(xué)科耦合協(xié)同研發(fā)與設(shè)計的重要手段，因此受到世界各國重視。各核電大國已投入大量資源進(jìn)行相關(guān)研究并組建了相應(yīng)核能數(shù)值仿真研究中心，并確定了具體研究項目[1-2]，如美國的CASL和NEAMS項目，歐洲地區(qū)的NURESIM系列項目等，中國核動力研究設(shè)計院聯(lián)合哈爾濱工程大學(xué)參考美國CASL項目確定了“數(shù)值反應(yīng)堆”研發(fā)項目。國內(nèi)也加大了推進(jìn)核能一體化數(shù)值建模與仿真技術(shù)開發(fā)力度，以便在新一輪核能研發(fā)領(lǐng)域趕上歐美發(fā)達(dá)國家的步伐。反應(yīng)堆一體化數(shù)值建模與仿真技術(shù)研究對象往往是設(shè)備級的，如全尺寸或1/2尺寸反應(yīng)堆，且需考慮多物理場耦合求解，如熱-流-固-物理耦合。反應(yīng)堆三維精細(xì)化多物理場耦合仿真的前提是解決各單一物理場的三維全尺寸級精細(xì)化求解難題，當(dāng)前仍有諸多基本科學(xué)與技術(shù)問題尚待明確和解決，本文所研究的反應(yīng)堆燃料組件三維精細(xì)化熱工水力數(shù)值仿真求解效率問題則是其中之一。

壓水堆堆芯包含數(shù)十至兩百多個燃料組件不等，典型燃料組件由17×17的燃料棒方形排列而成，同時在豎直方向上設(shè)置約8個起固定和攪混作用的攪混翼定位格架，總高度約4 m。初步估算，即使基于計算量較小的雷諾時均+湍流渦黏模型+多面體網(wǎng)格的計算方法，對單個完整燃料組件進(jìn)行三維精細(xì)化熱工水力數(shù)值求解，所需計算單元總數(shù)約為1億；采用大渦模擬(LES)結(jié)合壁面函數(shù)法則需至少10億網(wǎng)格單元[3]。盡管近20年來，國內(nèi)外針對CFD方法應(yīng)用于反應(yīng)堆燃料組件的相關(guān)研究成果層出不窮，且已證明是一種成熟有效的方法[4-8]，然而這一階段三維精細(xì)化熱工水力分析大多存在一個共同特點：仿真對象主要基于簡化模型或局部模型，如3×3、5×5等，與典型真實燃料組件17×17的規(guī)模還存在較大差距，主要原因一方面是早期研究處于驗證分析階段，采用小規(guī)模棒束可提高研究效率；另一方面則是受到計算能力與計算效率的限制。已有研究表明，真實燃料組件的熱工水力性能無法完全用小規(guī)模棒束通道來表征，不同棒束規(guī)模燃料組件熱工水力特性存在明顯差異[9-10]。可見，目前無論是針對新一輪核能技術(shù)研發(fā)、核安全分析及由此衍生的“虛擬反應(yīng)堆”研發(fā)，還是針對反應(yīng)堆燃料組件熱工水力特性理論研究本身，都對反應(yīng)堆燃料組件大規(guī)模棒束通道的三維精細(xì)化數(shù)值仿真有迫切需求。然而當(dāng)前面臨的一個非?，F(xiàn)實的問題是，如何解決反應(yīng)堆全尺寸級三維精細(xì)化數(shù)值仿真所面臨的網(wǎng)格數(shù)量過多、數(shù)值計算量過大及積累誤差過大等問題。

對于壓水堆燃料組件熱工水力數(shù)值求解，造成其數(shù)值計算量過大的根本原因在于基于歐拉法與貼體網(wǎng)格技術(shù)的三維精細(xì)化數(shù)值分析方法不得不采用數(shù)量極為龐大的“非結(jié)構(gòu)化網(wǎng)格”來表征復(fù)雜結(jié)構(gòu)的真實物理邊界，導(dǎo)致計算成本過高。以往針對這一問題的主要解決方法為采用多孔介質(zhì)模型或粗網(wǎng)格技術(shù)等降低計算量[11-14]。直到近年來，學(xué)術(shù)界才陸續(xù)出現(xiàn)專門針對反應(yīng)堆堆芯的三維精細(xì)化數(shù)值仿真計算量與求解效率的相關(guān)研究。Hu等[15-16]首先對鈉冷快堆中的繞絲棒束通道進(jìn)行相關(guān)研究，提出了一個動量源數(shù)理模型來表征繞絲對棒束通道熱工水力的影響。Capone等[17-18]針對壓水堆燃料組件，提出了一個基于高分辨率湍流模型(如大渦模擬)計算結(jié)果的動量源模型，并將此方法應(yīng)用到一5×5的棒束通道中。Mikuz等[19]同樣基于燃料組件完整貼體網(wǎng)格模型的計算結(jié)果，提取了相應(yīng)的動量源項，并在一3×3的裸棒束通道中實現(xiàn)了攪混翼格架對流場與溫度場影響的模擬。上述針對壓水堆的相關(guān)研究，大多基于完整燃料組件的貼體網(wǎng)格計算結(jié)果來提取并擬合出合適的動量源模型，本文則直接基于復(fù)雜結(jié)構(gòu)的影響機(jī)理建立其數(shù)學(xué)模型。

本文主要基于燃料組件攪混翼對流場及溫度場的影響機(jī)理，將復(fù)雜結(jié)構(gòu)對熱工水力特性的綜合影響歸結(jié)為流體質(zhì)點受到的某種虛擬體積力，只需給出這種力的數(shù)學(xué)描述并以動量源項的形式加入動量方程，即可在無翼片通道中獲得攪混通道的熱工水力特性。通過與實驗數(shù)據(jù)及傳統(tǒng)貼體網(wǎng)格建模計算結(jié)果的對比分析，驗證本文所建立虛擬體積力動量源模型的有效性。

1 虛擬體積力動量源建模

本文以Karoutas等[4]的攪混翼棒束通道激光多普勒測速(LDV)實驗為研究對象，建立壓水堆燃料組件攪混翼的虛擬體積力動量源模型。該實驗提供了1個5×5棒束通道攪混翼格架下游不同位置的流速測量結(jié)果，實驗?zāi)Ｐ图俺叽缛鐖D1所示。真實壓水堆燃料組件的攪混翼格架，除攪混翼外，還包含對燃料棒起約束作用的剛凸及彈簧，本文所研究的虛擬體積力動量源模型主要針對攪混翼，暫不考慮剛凸及彈簧。

圖1 5×5棒束實驗?zāi)Ｐ图皟赏ǖ烙嬎隳Ｐ虵ig.1 5×5 rod bundle experiment model and two-subchannel model

典型燃料組件的攪混翼具有周期性重復(fù)排列的特點，最小周期性單元為圖1中藍(lán)色虛線框所示的四通道模型。對于燃料組件的中心區(qū)域棒束通道，已有研究結(jié)果[10]表明，在數(shù)值建模與仿真中，四通道模型最小周期性單元還可采用圖2所示的周期性邊界交叉配對(P1a-P1b、P2a-P2b、P3a-P3b)的方式進(jìn)一步簡化成兩通道模型而不丟失其主要熱工水力特征，進(jìn)一步提高研究過程中的數(shù)值求解效率。對于帶攪混翼的棒束通道模型，由于每個狹長的子通道內(nèi)至少存在2個偏折方向不同的攪混翼，因此在對攪混翼棒束通道進(jìn)行數(shù)值求解時，基于貼體網(wǎng)格理論的CFD方法需將包含攪混翼格架的流體域單獨隔離出來進(jìn)行非結(jié)構(gòu)化網(wǎng)格建模，其他規(guī)則區(qū)域則將面網(wǎng)格拉伸成體網(wǎng)格以減少網(wǎng)格數(shù)量和提高網(wǎng)格質(zhì)量，如圖3所示?？梢姡瑪嚮煲砀窦艿拇嬖谑窃斐少N體網(wǎng)格建模數(shù)值計算量過大的根本原因。反應(yīng)堆堆芯燃料棒的排列本身非常規(guī)則，其數(shù)值分析可采用高質(zhì)量結(jié)構(gòu)化網(wǎng)格建模。如果能通過某種數(shù)學(xué)物理模型來描述攪混翼對棒束通道內(nèi)流場與溫度場的影響，則可將復(fù)雜攪混翼棒束通道的數(shù)值建模簡化為無翼片通道建模，從而大幅降低數(shù)值計算量。本文以圖1中的兩通道模型為研究對象，建立攪混翼的虛擬體積力數(shù)學(xué)模型，以考慮攪混翼對流場與溫度場的影響。

圖2 兩通道模型周期性邊界Fig.2 Periodic boundary conditions of two-subchannel model

圖3 燃料組件攪混翼棒束通道典型幾何與網(wǎng)格模型Fig.3 Typical geometrical and mesh models of fuel assembly rod bundle with mixing vane

通過理論分析可知，攪混翼對棒束通道流場與溫度場的影響主要體現(xiàn)在：攪混翼所產(chǎn)生的二次流及旋渦增強(qiáng)了流體的橫向運動及不同子通道間的交混能力，從而強(qiáng)化了流體與壁面間的對流換熱并降低了溫度不均勻性，同時增大了壓降。進(jìn)一步分析可知，換熱的強(qiáng)化及壓降的增大本質(zhì)上均由攪混翼所產(chǎn)生的二次流及漩渦的交混效應(yīng)所造成，而二次流與漩渦運動又主要由攪混翼的阻流和導(dǎo)流作用引起，棒束通道下游典型橫向流場及渦結(jié)構(gòu)如圖4所示。因此，可認(rèn)為導(dǎo)致上述各物理現(xiàn)象的原始驅(qū)動力是由于攪混翼對流場中的流體質(zhì)點施加了某種作用力，這種作用力改變了流體的運動規(guī)律并導(dǎo)致了復(fù)雜的流動現(xiàn)象(二次流、漩渦、攪混等)及換熱現(xiàn)象(降低溫度不均勻性、強(qiáng)化換熱)。因此，如果能基于攪混翼對流場與溫度場的作用機(jī)理給出這種力的數(shù)學(xué)描述，則在真實燃料組件的數(shù)值建模中就無需再對攪混翼進(jìn)行數(shù)量龐大的非結(jié)構(gòu)化貼體網(wǎng)格建模，從而簡化建模過程及大幅降低計算量、提高計算效率，最終實現(xiàn)在簡單規(guī)則的裸棒束模型通道中即可獲得攪混翼棒束通道中的復(fù)雜流場，而其他所有由流場變化引起的物理現(xiàn)象，如換熱的強(qiáng)化、溫度不均勻性的降低、壓降的增大等均可自然獲得。

圖4 攪混翼格架下游橫向流場及流型示意圖Fig.4 Lateral streamline and flow pattern downstream of mixing vane grid

攪混翼虛擬體積力示意圖如圖5所示。假設(shè)棒束通道中任一攪混翼對流體的虛擬體積力為Fb，結(jié)合翼片形狀可將其分解為3個方向的分力，即與攪混翼垂直的法向分力Fn、與攪混翼平行的分力Fp、與攪混翼相切的分力Ft，則該攪混翼的虛擬體積力合力如式(1)所示。還可進(jìn)一步將攪混翼各方向的虛擬體積力在直角坐標(biāo)系下分解，以便后續(xù)數(shù)值求解中動量源項的建立及加載，如式(2)所示。

Fb=Fn+Fp+Ft

(1)

(2)

式中，θ為攪混翼相對燃料棒軸向的偏折角，壓水堆燃料組件約為25°。結(jié)合實際攪混翼對流體的阻流、導(dǎo)流及攪混現(xiàn)象，從圖5可看出，垂直于攪混翼平面的力Fn一方面對流體流動產(chǎn)生阻礙作用，同時在x正方向提供一定程度的橫向驅(qū)動力；平行于攪混翼方向的力Fp主要提供額外的摩擦阻力，同時對流體在x軸負(fù)方向提供一定驅(qū)動力；與攪混翼相切的力Ft則主要影響流體在y方向的流動。這些力既能把流體導(dǎo)向某個方向并加速，也能對流體產(chǎn)生一定的阻礙作用，其綜合效果體現(xiàn)了攪混翼對流場的影響。

圖5 攪混翼虛擬體積力示意圖Fig.5 Schematic of virtual body force for mixing vane

根據(jù)式(3)所示的不可壓縮黏性流體動量守恒方程，攪混翼的虛擬體積力將以動量源項的形式加入到動量方程的體積力項f(N/m3)中，其物理意義為單位體積流體所受到的體積力。

(3)

式中，ρ、V、p、t分別為密度、流速、壓強(qiáng)及時間。

在攪混翼法向方向，流體質(zhì)點受到攪混翼的阻礙及導(dǎo)向作用后的加速度可由流速對時間的質(zhì)點導(dǎo)數(shù)(全導(dǎo)數(shù))來表示：

(4)

式中，nn、nt、np分別表示攪混翼的法向、平行及切向方向。穩(wěn)態(tài)工況下，式(4)右側(cè)的當(dāng)?shù)丶铀俣软?Vn/?t=0，只剩遷移加速度項。考慮翼片法向的體積力Fn是引起流體在該方向上動量變化的主要原因，則法向虛擬體積力導(dǎo)致的控制體單位體積內(nèi)流體動量的變化率可表示為：

(5)

式中：Vn、Vt、Vp分別表示翼片法向、平行及切向方向的流速；nn為法向單位矢量。對于式(5)中法向流速在各方向的速度梯度，本文結(jié)合攪混翼幾何特征，采用式(6)所示形式。

(6)

式中：k為乘數(shù)因子，可基于具體燃料組件的模型實驗數(shù)據(jù)對計算結(jié)果進(jìn)行修正；ln為攪混翼法向厚度；at、ap分別為切向厚度修正因子及平行厚度修正因子?？紤]到攪混翼實際導(dǎo)流效果，法向速度在翼片平行方向的梯度必定大于在翼片切向方向的梯度，因此需保證at>ap。

對于與攪混翼平面平行和相切的兩個方向的虛擬體積力Fp及Ft，考慮到這兩個方向的力主要由翼片與流體之間的切向力引起，因此基于達(dá)西公式建立其動量源模型：

(7)

式中：np、nt分別為攪混翼平行向及切向單位矢量；ζp、ζt分別為各自方向上的阻力系數(shù)，由McAdams關(guān)系式[20]確定。

ζ=0.184/Re0.23×104

(8)

由式(5)～(8)可看出，攪混翼對流體各方向的虛擬體積力均可表示成當(dāng)?shù)厮俣鹊暮瘮?shù)。數(shù)值求解過程如圖6所示，首先基于無攪混翼的棒束通道數(shù)值模型，根據(jù)初始流場計算原攪混翼所在空間網(wǎng)格節(jié)點各方向的虛擬體積力動量源項。將所計算的動量源項加載到動量守恒方程中并迭代，獲得新的流場后判斷結(jié)果是否收斂，若未收斂，則基于新的速度場重新計算動量源項，重復(fù)上述過程直至計算收斂。虛擬體積力動量源項與當(dāng)?shù)厮俣葓鲋g為自動耦合關(guān)系，基于當(dāng)?shù)厮俣鹊玫絼恿吭错棧瑒恿吭错椩僮饔糜诹鲌?，直到所計算的動量源項恰好能維持該速度場的穩(wěn)定?？梢?，整個動量源計算及求解過程只需基于當(dāng)?shù)厮俣?，無需引入其他變量。虛擬體積力模型在獲得攪混翼下游流場的同時，翼片對溫度場的影響可自然獲得。這是因為，所形成的橫向流場將一個通道的冷卻劑輸送至另一個通道降低了溫度分布不均勻性；各通道內(nèi)的漩渦增加了流體的攪混效應(yīng)，并迫使邊界層內(nèi)的高溫流體向主流低溫區(qū)運動，提高了加熱壁面與流體間的溫差，實現(xiàn)了換熱的強(qiáng)化。

圖6 數(shù)值求解過程Fig.6 Numerical solution process

圖7 攪混翼類型Fig.7 Type of mixing vane

從上述建模過程可看出，各方向的虛擬體積力動量源數(shù)學(xué)模型的具體形式及其在直角坐標(biāo)系下的分解與翼片位置及偏折方向有關(guān)。整個反應(yīng)堆堆芯燃料組件的攪混翼是4種類型翼片(圖7)周期性重復(fù)排列的結(jié)果，與坐標(biāo)軸相對位置及偏折方向相同的攪混翼具有相同形式的虛擬體積力動量源項。因此，只需對4種類型的攪混翼分別建立其虛擬體積力動量源模型，并加載到同類型的翼片上，即可實現(xiàn)完整燃料組件的數(shù)值建模。另外3種類型的攪混翼虛擬體積力動量源建模方法與翼片類型1類似。

2 數(shù)值驗證

基于Karoutas等[4]的攪混翼棒束通道激光多普勒測速(LDV)實驗，以圖1中兩通道模型為研究對象，分別建立其虛擬體積力動量源數(shù)值模型及貼體網(wǎng)格數(shù)值模型，其中側(cè)邊界采用圖2所示的交叉配對周期性邊界條件，以考慮該計算通道與周圍通道之間的交混。該實驗來流速度為6.79 m/s，溫度為26.67 ℃，系統(tǒng)壓力為48.3 kPa。為研究虛擬體積力動量源模型對溫度場的模擬效果，計算中對燃料棒施加707 kW/m2的定熱流密度。數(shù)值求解基于ANSYS FLUENTL，分別建立兩通道模型中4個攪混翼的虛擬體積力動量源數(shù)學(xué)模型，通過用戶自定義函數(shù)(UDF)以動量源項的形式加入動量方程中。通過對比分析虛擬體積力動量源模型計算結(jié)果、常規(guī)貼體網(wǎng)格計算結(jié)果及實驗數(shù)據(jù)來驗證模型的有效性。

對于攪混翼格架棒束通道的網(wǎng)格建模，當(dāng)攪混翼部分的網(wǎng)格基準(zhǔn)尺寸約為0.25 mm時，可避免網(wǎng)格尺寸對計算結(jié)果的影響[10]，當(dāng)采用SSTk-ω湍流模型結(jié)合y+為30左右的邊界層網(wǎng)格時，對流場、壓降及溫度場有較好的計算效果[8]?；趦赏ǖ滥Ｐ停謩e建立帶攪混翼的完整貼體網(wǎng)格模型及無攪混翼的混合網(wǎng)格模型，前者用于傳統(tǒng)貼體網(wǎng)格求解，后者用于虛擬體積力動量源求解。從虛擬體積力動量源數(shù)學(xué)建模過程可知，所建立的動量源模型須加載到原攪混翼所在區(qū)域。鑒于本文主要是驗證虛擬體積力動量源數(shù)學(xué)模型的有效性，在驗證計算時首先將無攪混翼模型的流體域分為兩部分，一部分為攪混翼所在區(qū)域的流體網(wǎng)格，另一部分為其他區(qū)域流體網(wǎng)格，以方便動量源模型加載到原翼片所在區(qū)域。經(jīng)虛擬體積力動量源計算結(jié)果與實驗數(shù)據(jù)及貼體網(wǎng)格建模結(jié)果的綜合對比分析后，確定虛擬體積力動量源模型參數(shù)如下：k=2.2，at=1.3，ap=1.7。

采用傳統(tǒng)貼體網(wǎng)格建模及虛擬體積力動量源建模所計算的橫向流速與實驗數(shù)據(jù)的對比示于圖8，數(shù)據(jù)均取自圖1中的實驗測量路徑。圖中，Vy為y方向的流速，Vavg為取值路徑上的平均流速，相對位置為x坐標(biāo)值與棒間距P的比值，z為所在路徑與攪混翼格架末端的軸向距離。從圖8可看出，基于貼體網(wǎng)格建模與虛擬體積力動量源建模所計算的橫向流速在分布趨勢及大小上均與實驗結(jié)果基本吻合，僅在攪混翼格架近下游區(qū)域的2個側(cè)端面處與實驗結(jié)果存在一定偏差。從動量源法與傳統(tǒng)貼體網(wǎng)格法的對比結(jié)果看，在攪混翼格架的近下游區(qū)域，二者非常接近，遠(yuǎn)離格架則出現(xiàn)一定偏差，動量源法計算的橫向流速衰減略偏大。

格架下游的橫向流場本質(zhì)上是攪混翼所引起的二次流現(xiàn)象，體現(xiàn)了攪混翼對下游流場的攪混與導(dǎo)流效果，也是換熱強(qiáng)化、壓降增大的根本原因。為更直觀地展示動量源模型對攪混翼下游橫向流場的模擬效果，采用兩種方法計算了格架下游不同位置橫向流場，結(jié)果示于圖9。從圖9可看出，兩種方法所計算的攪混翼格架下游橫向流場與圖4所示的典型橫向流型一致，通道中心存在1個大尺度漩渦，流體沿著燃料棒邊緣從一個通道流向另一個通道，且在通道出口處存在小尺度漩渦。動量源法所計算的通道中心大尺度漩渦位置及形狀與傳統(tǒng)貼體網(wǎng)格計算結(jié)果基本一致，但二者在z=12.7 mm截面燃料棒近壁面區(qū)域的小尺度漩渦計算結(jié)果存在一定差別，貼體網(wǎng)格計算結(jié)果為單個小尺度漩渦，動量源法結(jié)果則分裂為2個，且位置略有偏差。圖中各截面云圖為軸向流速分布，可見動量源法計算的軸向流速的大小及分布與貼體網(wǎng)格結(jié)果也基本吻合。

圖8 攪混翼格架下游橫向流速分布Fig.8 Lateral velocity profile downstream of mixing vane grid

圖9 攪混翼格架下游橫向流場對比Fig.9 Comparison of lateral flow fields downstream of mixing vane grid

攪混翼格架及其下游壓降計算結(jié)果示于圖10，圖中實驗關(guān)聯(lián)式結(jié)果基于Chun等[21]的半經(jīng)驗公式，該壓降關(guān)聯(lián)式由格架形狀阻力、格架摩擦阻力、攪混翼形狀阻力及燃料棒摩擦阻力4部分組成。從圖10可看出，動量源法與傳統(tǒng)貼體網(wǎng)格法所計算的各段壓降與實驗關(guān)聯(lián)式均較為接近，局部最大偏差在10%以內(nèi)。動量源法與貼體網(wǎng)格法計算的壓降與實驗關(guān)聯(lián)式的偏差主要出現(xiàn)在攪混翼格架段，格架下游區(qū)域的摩擦阻力段則十分接近，但與實驗關(guān)聯(lián)式相比均略偏高。對比兩種計算方法的結(jié)果可知，二者壓降曲線在攪混翼格架下游部分重合度較高，說明動量源法很好地捕捉到了格架下游的阻力特性。

圖10 攪混翼棒束通道壓降Fig.10 Pressure loss along rod bundle with mixing vane grid

貼體網(wǎng)格法及虛擬體積力動量源法所計算的Nusselt數(shù)(Nu)沿燃料棒軸向及周向的分布規(guī)律示于圖11。Nu的物理意義為壁面附近流體的無量綱溫度梯度，其大小表征了流體與固體壁面間對流換熱的強(qiáng)弱。從圖11可看出，動量源法所計算的Nu沿軸向的分布與傳統(tǒng)貼體網(wǎng)格計算結(jié)果吻合良好，最大偏差在5%以內(nèi)。周向分布上，二者在攪混翼格架的近下游區(qū)域(12.7 mm)吻合較好，在遠(yuǎn)下游區(qū)域(50.4 mm)則存在一定偏差。遠(yuǎn)下游區(qū)域的偏差主要體現(xiàn)在Nu的局部極值在方位角上有所偏移，且該偏移與圖8中的橫向速度偏差存在一定關(guān)聯(lián)性，原因在于Nu的局部大小主要受當(dāng)?shù)亓魉俚挠绊??？傮w上，虛擬體積力動量源法所計算的Nu無論在數(shù)值還是分布規(guī)律上，均與傳統(tǒng)貼體網(wǎng)格法計算結(jié)果吻合較好。從本文動量源數(shù)理建模及計算過程來看，在未對溫度場求解進(jìn)行任何特殊處理的情況下較好地捕捉到了攪混翼格架對棒束下游傳熱特性的影響，說明本文所給出的攪混翼對棒束通道流場與溫度場的影響機(jī)理分析是合理的，也驗證了虛擬體積力動量源建模思想的可行性。

為進(jìn)一步綜合評估虛擬體積力動量源模型對攪混翼棒束通道下游流動與傳熱特性的數(shù)值模擬效果，對攪混翼格架下游燃料棒近壁面對數(shù)律區(qū)流速沿軸向分布進(jìn)行對比，結(jié)果示于圖12。其中Dh為棒束通道水力直徑。近壁面對數(shù)律區(qū)是邊界層中的重要區(qū)域，該區(qū)域流速對固體壁面的對流換熱起決定性作用，其大小與分布體現(xiàn)了攪混翼對下游流場與溫度場的影響。從圖12可看出，基于虛擬體積力動量源法所計算的燃料棒表面對數(shù)律區(qū)流速，無論是分布趨勢還是大小均與完整貼體網(wǎng)格建模結(jié)果基本吻合。虛擬體積力動量源所計算的近壁面對數(shù)律區(qū)流速雖總體上略高于傳統(tǒng)貼體網(wǎng)格計算結(jié)果，但最大偏差小于5%。此外，對比圖11可知，Nu與對數(shù)律區(qū)流速的軸向分布存在一定相似性，符合流動與換熱關(guān)系的理論預(yù)期。

圖11 燃料棒軸向及周向Nu分布對比Fig.11 Comparison of axial and azimuthal variations in Nusselt number

為驗證虛擬體積力動量源模型對不同來流工況的適應(yīng)性，對不同流速(V0)工況下棒束通道內(nèi)橫向流速及對流換熱系數(shù)與傳統(tǒng)貼體網(wǎng)格建模法進(jìn)行對比，結(jié)果示于圖13、14。從圖13可看出，不同來流工況下，兩種方法所計算的格架下游12.7 mm處的橫向流速均十分接近。可見從流場的角度，虛擬體積力動量源模型在各流速工況下均具有較好的適應(yīng)性。從圖14可看出，兩種方法的計算結(jié)果在小流速工況下吻合更好，隨著流速的增大，二者偏差呈現(xiàn)略增大的趨勢?？傮w上，不同來流工況下虛擬體積力動量源法與傳統(tǒng)貼體網(wǎng)格建模法所計算的對流換熱系數(shù)在大小及分布趨勢上均有較高的吻合度。

圖12 燃料棒近壁面對數(shù)律區(qū)速度分布Fig.12 Profile of velocity in log-law region near rod wall

為明確虛擬體積力動量源法對棒束通道全結(jié)構(gòu)化網(wǎng)格模型的敏感性，以評估計算效率的提升效果，本文進(jìn)一步對兩通道模型進(jìn)行了六面體全結(jié)構(gòu)化網(wǎng)格建模。求解時先標(biāo)記原攪混翼所在空間的網(wǎng)格單元，用于動量源模型的加載。計算表明，當(dāng)采用全六面體結(jié)構(gòu)化網(wǎng)格建模，兩通道模型網(wǎng)格單元數(shù)達(dá)到10.8萬時，動量源計算結(jié)果不受網(wǎng)格影響，與貼體網(wǎng)格的對比如圖15所示。

基于結(jié)構(gòu)化和非結(jié)構(gòu)化網(wǎng)格建模時動量源法所計算的當(dāng)?shù)貦M向速度及Nu軸向分布與帶攪混翼貼體網(wǎng)格結(jié)果的對比示于圖16。從圖16可看出，當(dāng)采用全結(jié)構(gòu)化網(wǎng)格時，動量源法所計算的橫向速度及Nu均與非結(jié)構(gòu)化網(wǎng)格動量源模型及貼體網(wǎng)格模型整體吻合良好，均捕捉到了各參數(shù)的主要變化特征。仔細(xì)觀察橫向速度對比結(jié)果可知，結(jié)構(gòu)化網(wǎng)格動量源模型所計算的攪混翼格架下游12.7 mm處的橫向速度在左側(cè)通道中較另外兩種方法偏低，速度變化更趨平緩?？赏茰y，該偏差產(chǎn)生的原因在于全結(jié)構(gòu)化網(wǎng)格模型中施加動量源的網(wǎng)格單元無法與真實結(jié)構(gòu)的邊界完全貼合，丟失了部分結(jié)構(gòu)特征對流場的影響。上述分析表明，基于全結(jié)構(gòu)化網(wǎng)格的動量源法在較好地預(yù)測到攪混翼格架下游流動與傳熱特性的同時，可顯著提高棒束通道的數(shù)值求解效率。以本文單跨兩通道模型為例，采用傳統(tǒng)貼體網(wǎng)格建模所需單元數(shù)為43.2萬，基于結(jié)構(gòu)化網(wǎng)格的動量源法則只需10.8萬，網(wǎng)格數(shù)量減少了76%?？深A(yù)測，將動量源法應(yīng)用于單個完整燃料組件或堆芯，計算效率的提升將更可觀。

圖13 z=12.7 mm處不同流速工況下橫向速度計算結(jié)果對比Fig.13 Comparison of calculated lateral velocities under different flow rates at z=12.7 mm

圖14 不同流速工況下?lián)Q熱系數(shù)計算結(jié)果對比Fig.14 Comparison of calculated heat transfer coefficients under different flow rates

圖15 用于動量源模型的結(jié)構(gòu)化網(wǎng)格與用于完整結(jié)構(gòu)的貼體網(wǎng)格Fig.15 Structured mesh for momentum source model and body-fitted mesh for full structure

圖16 基于不同網(wǎng)格的動量源模型與貼體網(wǎng)格模型結(jié)果對比Fig.16 Comparison between results of momentun source model based on different meshes and body-fitting mesh model

3 結(jié)論

1) 壓水堆燃料組件棒束通道中的攪混翼對流動與傳熱特性的綜合影響，可由復(fù)雜結(jié)構(gòu)對流體質(zhì)點產(chǎn)生的某種虛擬體積力來表征，只需對這種體積力建立合適的數(shù)學(xué)模型并以動量源形式加入無攪混翼棒束通道模型中，即可同時獲得攪混翼對棒束通道內(nèi)流場與溫度場的影響。

2) 基于攪混翼的阻流、導(dǎo)流、攪混及換熱強(qiáng)化等機(jī)理所構(gòu)建的虛擬體積力動量源模型可表示成當(dāng)?shù)厮俣葓龅暮瘮?shù)，無需引入其他變量，且數(shù)值求解為自動耦合過程，當(dāng)計算收斂時，合適的動量源大小即被確定。

3) 虛擬體積力動量源模型所計算的橫向流場與實驗數(shù)據(jù)及體貼網(wǎng)格模型結(jié)果相比，在攪混翼格架近下游區(qū)域與二者均吻合較好，遠(yuǎn)下游區(qū)域則出現(xiàn)流速衰減偏大的現(xiàn)象。針對攪混翼棒束通道的壓降、Nu、流速等關(guān)鍵熱工水力參數(shù)，虛擬體積力動量源模型計算結(jié)果在數(shù)值與分布趨勢上均與貼體網(wǎng)格模型結(jié)果基本吻合，局部最大偏差在5%以內(nèi)。此外，虛擬體積力動量源模型對棒束通道來流速度不敏感，能適應(yīng)不同流速工況。驗證結(jié)果證明了虛擬體積力動量源法的合理性與可行性。

本文研究主要針對攪混翼，后續(xù)將研究其他復(fù)雜結(jié)構(gòu)(如格架帶條、剛凸與彈簧等)的類似建模方法，有望將燃料組件簡化為裸棒束通道而不丟失攪混翼格架的主要熱工水力特性，進(jìn)一步提高精細(xì)化建模與仿真效率。