計算流體力學應(yīng)用于子通道分析的初步方法研究

臧金光，閆 曉，黃善仿，曾小康，黃彥平

(1.中國核動力研究設(shè)計院 中核核反應(yīng)堆熱工水力技術(shù)重點實驗室，四川 成都 610041；

2.清華大學 工程物理系，北京 100084)

實際反應(yīng)堆系統(tǒng)存在3種主要結(jié)構(gòu)尺度：1) 系統(tǒng)分析尺度，著眼于整個反應(yīng)堆系統(tǒng)的運行特性和響應(yīng)特性；2) 子通道分析尺度，對反應(yīng)堆堆芯采用準二維方法進行建模分析；3) 計算流體力學(CFD)尺度，針對系統(tǒng)中具有明顯三維流場特征的區(qū)域進行計算分析。在開展反應(yīng)堆系統(tǒng)分析時，為尋求精度與效率之間的平衡，同時兼顧不同尺度的物理現(xiàn)象，一個較為普遍的方法是采用多尺度耦合，針對不同尺度的結(jié)構(gòu)選用相應(yīng)尺度的分析程序。一些國家正在積極推進多物理、多尺度、多學科的數(shù)字化耦合平臺研發(fā)，如美國的虛擬反應(yīng)堆計劃，歐盟的NURESIM耦合平臺[1]等，還有一些研究者開發(fā)了一些具體的耦合程序，如CAST3M/ARCTURUS程序[2]，COBRA/CFX耦合程序[3]，RELAP5/CFD耦合程序[4]等。多物理、多尺度的耦合研究是當前研究的熱點之一。

在開展熱工水力多尺度分析時，人為地將各個尺度分解，針對不同尺度發(fā)展了各自的分析程序，然后再建立一耦合平臺將分解后的程序重新組合在一起。實際上，各尺度的分析方法并不存在本質(zhì)差別，它們均源于基本的流體力學控制方程，均基于相同或相近的本構(gòu)關(guān)系式。例如，計算流體力學尺度的控制方程經(jīng)過一定方法的統(tǒng)計平均和適當簡化后可變成子通道控制方程和系統(tǒng)控制方程，因此，計算流體力學方法在理論上已具備求解子通道尺度和系統(tǒng)尺度的條件，從而為計算流體力學程序開展子通道尺度分析和系統(tǒng)尺度分析提供可能。

本工作即是基于計算流體力學方法開展子通道分析的一種嘗試。本工作研究計算流體力學控制方程與子通道控制方程間的聯(lián)系，然后通過一具體組件計算對本工作提出的子通道CFD方法進行比較和分析，為多尺度耦合分析提供一種新思路。

1 控制方程比較

計算流體力學方法所基于的控制方程是局部瞬時Navier-Stokes方程。局部瞬時Navier-Stokes方程是研究各尺度流體流動和傳熱的理論基礎(chǔ)，子通道控制方程的理論形式是基于局部瞬時方程通過一定的數(shù)學變換和簡化假設(shè)推導(dǎo)得到的[5-6]，借助的主要數(shù)學定理包括Leibnitz定理和Gauss定理。Leibnitz定理建立了空間積分與時間微分的交換關(guān)系，可將一體積分的時間變化率變換成一時間變化率的體積分與表面通量之和。Gauss定理揭示了高維積分與低維積分之間的關(guān)系，將散度的體積分轉(zhuǎn)化成它的表面通量。推導(dǎo)子通道控制方程的重要假設(shè)是認為子通道的軸向流動起主要作用，子通道間的橫向流動起次要作用，同時由于軸向尺度遠大于橫向尺度，一些量的軸向梯度可忽略。

進行積分平均不可避免地會帶來信息缺失，引入信息缺失的因素包括：1) 基于類似邊界層理論的假設(shè)，認為軸向尺度遠大于橫向尺度，一些物理量如黏性應(yīng)力、雷諾應(yīng)力等的軸向梯度項可忽略；2) 忽略了各物理場的細致截面分布，認為積分運算與其他代數(shù)運算具有普遍的可交換性，如認為〈ρv〉=C1〈ρ〉〈v〉，〈v2〉=C2〈v〉2，其中，ρ為密度，v為速度，分布系數(shù)C1和C2一般假設(shè)為1；3) 由于認為橫向動量方程在子通道控制方程中不起重要作用，因此在處理橫向動量方程時作了很多簡化，子通道交界面處的湍流交混通過簡單的交混因子進行封閉。大部分情況下，這些假設(shè)均是合理的。

將CFD工具運用于子通道分析的前提是二者的控制方程具有統(tǒng)一性，為說明這一問題，以連續(xù)性方程為例來介紹二者之間的聯(lián)系[5-6]。

CFD方法的連續(xù)性方程的形式為：

對上式進行體積平均，認為控制體固定，且積分與微分之間可交換：

其中：〈〈〈 〉〉〉為體積平均算符；t為時間坐標；x、y、z為空間坐標。

再應(yīng)用Gauss定理，轉(zhuǎn)化成子通道的控制方程形式為：

即經(jīng)過積分后的子通道質(zhì)量守恒方程本質(zhì)上相當于以子通道平均量為變量的局部微分方程。對于動量守恒方程和能量守恒方程，可采用類似的方法分析。另外基于CFD工具的子通道分析方法也可考慮子通道間的橫向流動，這種橫向流動主要依靠橫向的壓力梯度驅(qū)動。CFD控制方程與子通道控制方程間的聯(lián)系為將CFD方法應(yīng)用于子通道分析打下了基礎(chǔ)。

2 方法驗證及結(jié)果分析

為進一步驗證將CFD方法應(yīng)用于子通道分析的可行性，選取了中國核動力研究設(shè)計院提出的超臨界水冷堆概念設(shè)計方案中的CSR1000組件作為分析對象，分別用基于CFX軟件的精細CFD方法(CFX-CFD)、基于CFX軟件的子通道方法(CFX-SUB)及基于子通道程序ATHAS的子通道分析方法開展了計算分析，并比較了三者計算結(jié)果的差異。

CSR1000組件采用小組件的設(shè)計方案[7]，1個大組件含有4個小組件，每個小組件由φ9.5 mm燃料元件呈9×9方形布置，中間水棒占有5×5柵元位置，燃料元件分別位于周圍兩排，具體幾何參數(shù)及計算工況列于表1。CSR1000原型組件采用繞肋進行自定位，這里先對光滑組件進行分析，為繞肋的性能評價提供參考。由于組件在幾何上具有1/8對稱性，為節(jié)省計算資源，選取1/8結(jié)構(gòu)進行分析，組件形狀及各子通道編號示于圖1。

表1 幾何參數(shù)及計算工況

圖1 CSR1000組件及子通道編號

CFX-CFD方法即為常規(guī)的CFD方法，對組件通道劃分精細網(wǎng)格，然后開展CFD計算。本次計算選用SST湍流模型。為保證計算結(jié)果的有效性，開展了網(wǎng)格的無關(guān)性分析。圖2示出不同網(wǎng)格參數(shù)下通道摩擦系數(shù)與中心通道平均換熱系數(shù)的對比。由圖2可見，當網(wǎng)格數(shù)為230萬時，就已滿足了網(wǎng)格無關(guān)性要求。

圖2 網(wǎng)格無關(guān)性分析結(jié)果

CFX-SUB方法是利用CFX工具完成子通道分析的功能，以每個冷卻劑通道作為網(wǎng)格的劃分尺度，每個子通道截面只劃分1個網(wǎng)格。這樣進行網(wǎng)格劃分的結(jié)果是CFX軟件無法計算壁面的剪切應(yīng)力分布，需通過附加體積力的方式把摩擦力加入到軸向動量方程中。選用Blausius公式計算摩擦系數(shù)。燃料元件表面的熱流密度也需采用轉(zhuǎn)化成體積熱源的方式添加到軸向能量方程中，記表面熱流密度為qs，子通道流通面積為A，熱周為Pt，熱力直徑為Dt，轉(zhuǎn)化后的等效體積熱流密度qv的表達式為：

在CFX軟件中，將每個子通道作為一流體計算域，子通道與子通道之間通過交界面連接，在每個計算域中設(shè)置對應(yīng)的體積熱源和阻力源。由于只有1層網(wǎng)格，CFX軟件無法識別壁面，也無法對壁面進行求解，因此每個計算域除入口、出口、交界面等邊界條件外，其他邊界均設(shè)置成對稱邊界條件。當子通道數(shù)目較多時，手動為每個子通道設(shè)置條件會很繁瑣，CFX軟件支持CCL命令語言和Perl腳本語言，用于對邊界條件進行批量設(shè)置，可有效簡化設(shè)置的工作量。

由于CFX-SUB方法只能得到子通道的平均質(zhì)量流速和平均溫度，因此首先借助CFX-CFD方法初步分析組件內(nèi)的流場分布特征。圖3示出基于CFX-CFD方法得到的出口質(zhì)量流速和溫度的分布。CFX-CFD方法可得到通道截面各位置的流場信息，截面質(zhì)量流速與溫度的分布在子通道中心與最窄間隙處存在明顯的分布不均。由于子通道中心的流通截面較大，流動阻力較小，流體在各子通道進行流量分配時更傾向于流向?qū)捦ǖ捞?，而窄通道的質(zhì)量流速相對較低，局部位置的最大質(zhì)量流速和最小質(zhì)量流速相差接近兩倍。這種不均衡的質(zhì)量流量分配進一步造成同一截面上的溫度分布存在明顯差異。由于超臨界水在擬臨界點前后物性變化較大，窄間隙內(nèi)的流體會提前跨過擬臨界點，導(dǎo)致密度迅速下降，這又反過來繼續(xù)增加了流動阻力，迫使流體繼續(xù)向?qū)捦ǖ捞幜鲃?。這樣，在元件與元件之間以及元件與組件盒之間的狹窄間隙內(nèi)，很容易出現(xiàn)局部的過熱流體。這說明稠密布置的燃料組件中需要考慮采用一些附加設(shè)計(如繞絲等)以避免局部高溫流體的出現(xiàn)。

子通道分析程序ATHAS[5]可用來開展超臨界條件下的子通道分析。圖4比較了CFX-CFD、CFX-SUB、ATHAS 3種方法計算的出口焓與質(zhì)量流速G的分布。由圖4可見，各子通道中，3種方法預(yù)測的出口焓與質(zhì)量流速分布趨勢整體相同。相同類型子通道的出口焓和質(zhì)量流速相差較小，如中心通道(編號為2，5，8，11)的質(zhì)量流速穩(wěn)定在1 000 kg/(m2·s)左右，邊通道(編號為3，6，7，9，10，12)的質(zhì)量流速在700 kg/(m2·s)左右。盡管相同類型子通道由于與相鄰?fù)ǖ澜换煨阅艿牟煌兴町悾w上保持了相似的熱工水力特性。另外，不同類型子通道由于在幾何結(jié)構(gòu)方面的顯著差異，在吸熱量和質(zhì)量流速方面有明顯不同。CFX-SUB方法預(yù)測的各子通道平均值的整體趨勢與CFX-CFD方法的較為接近。以CFX-CFD方法計算的結(jié)果作為參考，CFX-SUB計算的平均焓的最大相對誤差出現(xiàn)在1號子通道，平均質(zhì)量流速的最大相對誤差出現(xiàn)在8號子通道，二者的最大相對誤差均不超過20%。CFX-SUB方法可反映出角通道過熱、中心通道較為均勻的變化趨勢，可見CFX-SUB方法的預(yù)測結(jié)果是合理的。另一方面，CFX-SUB方法預(yù)測的各子通道之間的差別比CFX-CFD方法和ATHAS程序的計算結(jié)果要大些，CFX-SUB方法預(yù)測的1號子通道的流體溫度過高與其預(yù)測的質(zhì)量流速較低相對應(yīng)。

圖3 CFX-CFD方法的計算結(jié)果

圖4 3種方法計算結(jié)果對比

為進一步分析CFX-SUB與CFX-CFD方法計算結(jié)果的異同，圖5進一步比較了二者在角通道(編號為1)、中心通道(編號為2)和邊通道(編號為3)3種典型子通道的沿程溫度T和質(zhì)量流速的變化。由圖5可見：CFX-SUB方法預(yù)測的整體趨勢與CFX-CFD方法的一致，中心通道沿程質(zhì)量流速上升，邊通道和角通道沿程質(zhì)量流速下降；角通道的溫度分布上升較快，其次是邊通道，中心通道的溫度上升最為平緩。盡管如此，CFX-SUB方法預(yù)測的角通道的流體溫度明顯高于CFX-CFD方法，而預(yù)測的角通道的質(zhì)量流速的差別卻沒有這么大。

CFX-SUB方法未正確考慮子通道交界面上的湍流交混，導(dǎo)致子通道間流體熱交換能力減弱。子通道交界上的湍流交混是子通道間能量交換的重要方式，由于交界面上的渦團脈動，使得兩個子通道內(nèi)的流體微團產(chǎn)生凈能量交換。在推導(dǎo)子通道控制方程中，引入了湍流熱通量的二階相關(guān)統(tǒng)計項，使得方程組不封閉。為解決這種不封閉問題，子通道程序中通常定義湍流交混系數(shù)——無量綱數(shù)β來表征相鄰子通道之間的能量交換量。

圖5 CFX-CFD與CFX-SUB方法計算結(jié)果的比較

為了分析湍流交混系數(shù)對子通道分析計算結(jié)果的影響，使用ATHAS程序計算了不同湍流交混系數(shù)下各子通道的溫度t的分布(圖6)。隨著湍流交混系數(shù)的降低，角通道的流體出口溫度逐步上升，當湍流交混系數(shù)降至0.001時，ATHAS計算的角通道流體出口溫度與CFX-SUB方法計算的結(jié)果相近。可見，CFX-SUB方法計算結(jié)果偏高的原因很有可能是未考慮子通道間的湍流交混，需要在后續(xù)的工作中改進。

圖6 湍流交混因子影響分析

3 結(jié)論

多尺度熱工水力耦合是目前開展反應(yīng)堆安全分析的一種發(fā)展趨勢，普通的做法是通過建立不同程序的數(shù)據(jù)交換平臺來實現(xiàn)。本工作提出另外一種開展多尺度耦合的思路，以計算流體力學方法為手段，開展宏觀的多尺度分析，并以子通道分析為例驗證這種方法，得到的主要結(jié)論如下。

1) 計算流體力學方法既可用于滿足連續(xù)介質(zhì)假設(shè)的流體微團分析，又可用于子通道尺度或更大尺度的流體問題分析，區(qū)別不同尺度的模擬在于網(wǎng)格尺度的大小及對于幾何結(jié)構(gòu)的描述能力。通過對子通道控制方程和局部瞬時控制方程的對比，認為子通道守恒方程本質(zhì)上相當于以子通道平均量為變量的局部微分方程，從理論上說明了計算流體力學應(yīng)用于子通道分析的可能性。

2) 以一典型組件為例，分別采用計算流體力學方法、子通道方法及計算流體力學應(yīng)用于子通道這3種方法進行計算，并對計算結(jié)果進行比較分析，說明了計算流體力學應(yīng)用于子通道方法的合理性以及具體操作的可實現(xiàn)性。另外，由于新方法中尚未考慮子通道間的湍流交混，使得這種方法仍存在進一步的改進空間。

參考文獻：

[1] CHAULIAC C, ARAGONéS J M, BESTION D, et al. NURESIM: A European simulation platform for nuclear reactor safety: Multi-scale and multi-physics calculations, sensitivity and uncertainty analysis[J]. Nuclear Engineering and Design, 2011, 241(9): 3 416-3 426.

[2] STUDER E, BECCANTINI A, GOUNANDA S, et al. CAST3M/ARCTURUS: A coupled heat transfer CFD code for thermal-hydraulic analyzes of gas cooled reactors[J]. Nuclear Engineering and Design, 2007, 237(15): 1 814-1 828.

[3] 劉余，張虹，賈寶山. COBRA-Ⅳ與CFX程序耦合研究[J]. 核動力工程，2010，31(2)：43-46.

LIU Yu, ZHANG Hong, JIA Baoshan. Research on coupling of COBRA-Ⅳ and CFX[J]. Nuclear Power Engineering, 2010, 31(2): 43-46(in Chinese).

[4] ANDERSON N, HASSAN Y, SCHULTZ R. Analysis of the hot gas flow in the outlet plenum of the very high temperature reactor using coupled RELAP5-3D system code and a CFD code[J]. Nuclear Engineering and Design, 2008, 238(1): 274-279.

[5] 單建強，李昌瑩，陳偉. 子通道分析程序ATHAS理論說明書[R]. 西安：西安交通大學，2008.

[6] 郝老迷. 沸騰傳熱與氣(汽)液兩相流動[M]. 北京：核工業(yè)研究生部，2007.

[7] 夏榜樣，楊平，王連杰，等. 超臨界水冷堆CSR1000初步概念設(shè)計[J]. 核動力工程，2013，34(1)：9-14.

XIA Bangyang, YANG Ping, WANG Lianjie, et al. Core preliminary conceptual design of supercritical water-cooled CSR1000[J]. Nuclear Power Engineering, 2013, 34(1): 9-14(in Chinese).