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

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

(1.中國(guó)核動(dòng)力研究設(shè)計(jì)院 中核核反應(yīng)堆熱工水力技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610041；

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

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

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

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

1 控制方程比較

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

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

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

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

對(duì)上式進(jìn)行體積平均，認(rèn)為控制體固定，且積分與微分之間可交換：

其中：〈〈〈 〉〉〉為體積平均算符；t為時(shí)間坐標(biāo)；x、y、z為空間坐標(biāo)。

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

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

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

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

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

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

圖1 CSR1000組件及子通道編號(hào)

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

參考文獻(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]. 核動(dòng)力工程，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] 單建強(qiáng)，李昌瑩，陳偉. 子通道分析程序ATHAS理論說(shuō)明書(shū)[R]. 西安：西安交通大學(xué)，2008.

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

[7] 夏榜樣，楊平，王連杰，等. 超臨界水冷堆CSR1000初步概念設(shè)計(jì)[J]. 核動(dòng)力工程，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).