外延分區(qū)耦合瞬態(tài)計(jì)算流體力學(xué)計(jì)算技術(shù)

張歷軒， 陳廣亮， 田兆斐， 錢(qián)浩， 李瑞

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

堆芯的精細(xì)化瞬態(tài)計(jì)算流體力學(xué)(calculated fluid dynamics，CFD)是了解堆芯瞬態(tài)熱工水力特性，發(fā)現(xiàn)堆芯薄弱環(huán)節(jié)，優(yōu)化堆內(nèi)構(gòu)件設(shè)計(jì)的重要手段[1-2]。典型壓水堆(pressurized water reactor，PWR)堆芯存在數(shù)萬(wàn)個(gè)狹窄流道，且堆芯內(nèi)部包含千余層定位格架，每層格架含數(shù)百個(gè)剛凸、彈簧、攪混翼等具備彎折、扭曲特征的復(fù)雜結(jié)構(gòu)[3]。針對(duì)壓水堆中的5×5典型燃料組件，其需要0.9億的網(wǎng)格量，400 G的內(nèi)存資源。無(wú)法依靠常規(guī)的方法對(duì)其進(jìn)行瞬態(tài)CFD分析[4]。

針對(duì)大規(guī)模堆芯穩(wěn)態(tài)CFD仿真，Navarro等[5]在5×5棒束通道數(shù)值研究中，其將研究流域劃分成7個(gè)子域，并使用壓降作為參考驗(yàn)證分段耦合技術(shù)準(zhǔn)確性，但未對(duì)子域的劃分方法進(jìn)行討論。李小暢[6]針對(duì)軸向計(jì)算域的分段方式及子計(jì)算域出口邊界條件類(lèi)型對(duì)分段求解技術(shù)的數(shù)值模擬結(jié)果影響進(jìn)行了討論。陳廣亮[4]針對(duì)5×5的棒束燃料組件通道，使用分段耦合方法完成了全高度的CFD穩(wěn)態(tài)計(jì)算，并根據(jù)冷卻劑通道不同區(qū)域的特征，采用不同湍流模型耦合來(lái)提高計(jì)算效率。然而，現(xiàn)階段缺少分段耦合技術(shù)在瞬態(tài)CFD仿真中的應(yīng)用研究。

本文基于穩(wěn)態(tài)分段耦合技術(shù)，設(shè)計(jì)了瞬態(tài)分段耦合CFD計(jì)算方案，開(kāi)發(fā)了瞬態(tài)調(diào)控程序，并通過(guò)外延分區(qū)優(yōu)化了子域劃分方案，使其適用于瞬態(tài)耦合，并使用階躍型和周期性的入口條件對(duì)其進(jìn)行驗(yàn)證，討論了其適用性條件。

1 計(jì)算流體力學(xué)計(jì)算模型

1.1 物理模型及數(shù)值計(jì)算方法

選取壓水堆燃料組件帶攪混葉片、3跨高度的AFA-2G 1×2格架模型開(kāi)展研究，格架模型如圖1[4]所示。棒徑為9.5 mm，柵距為12.6 mm，棒壁間距為2.55 mm,定位格架跨距520 mm。

圖1 5×5 AFA2G 格架示意Fig.1 The structure for the 5×5 AFA 2G grid

對(duì)其劃分子域的主要依據(jù)是不同堆內(nèi)構(gòu)件處的流場(chǎng)特性，光棒冷卻劑通道(coolant subchannel, CS)是反應(yīng)堆內(nèi)冷卻劑流動(dòng)和換熱的主要區(qū)域，流動(dòng)阻力較小。定位格架(spacer grid, SG)是對(duì)燃料組件起到支撐和固定作用,其流動(dòng)阻力較大。攪混翼(mixing vane, MV)則是坐落在部分定位格架上，它可以增加冷卻劑的橫向攪混[4]。

本文流動(dòng)介質(zhì)為水，入口溫度為292.8 ℃，速度為3.76 m/s。設(shè)置周期性邊界條件，壁面使用自定義函數(shù)加載熱流密度，無(wú)滑移壁面，出口設(shè)置為自由流出。使用Fluent進(jìn)行計(jì)算，壓力-速度耦合選擇SIMPLEC算法，湍流動(dòng)能、耗散、能量項(xiàng)采用二階迎風(fēng)格式進(jìn)行離散求解。能量方程收斂殘差設(shè)置為1×10-6，其他項(xiàng)設(shè)置為1×10-4,采用標(biāo)準(zhǔn)壁面函數(shù)。

湍流模型選擇雷諾應(yīng)力模型(Reynolds stress model,RSM)。雖然RSM對(duì)計(jì)算資源需求較大[7]，但由于冷卻劑在攪混翼區(qū)存在較強(qiáng)的橫向攪混，具有較強(qiáng)的各項(xiàng)異性，采用RSM模型計(jì)算結(jié)果較為準(zhǔn)確[8]。因此各子域間依次傳遞的信息不僅包括溫度、壓力、速度、湍動(dòng)能、耗散率等，還應(yīng)包括雷諾應(yīng)力張量。RSM的控制方程推導(dǎo)如下所示：

湍流場(chǎng)中，各物理量均可以被分解成平均值與脈動(dòng)值之和：

(1)

雷諾通過(guò)對(duì)納維-斯托克斯方程(N-S)進(jìn)行時(shí)間平均處理可得：

(2)

引入不封閉項(xiàng)-雷諾應(yīng)力[9-10]。針對(duì)雷諾應(yīng)力直接構(gòu)造封閉方程，即為雷諾應(yīng)力模型[11]。針對(duì)其所構(gòu)造的輸運(yùn)方程為[12]:

(3)

針對(duì)右式各項(xiàng)分別建模，從而完成N-S方程的封閉。

將3跨高度1×2燃料組件從入口端至出口端，采用外延分區(qū)方案，按照光棒冷卻劑通道、定位格架、攪混翼，依次劃分為10個(gè)計(jì)算子域，具體子域劃分方案見(jiàn)圖2。針對(duì)3跨高度的1×2燃料組件和各個(gè)子域進(jìn)行網(wǎng)格繪制。整體網(wǎng)格與各個(gè)子域網(wǎng)格均在光棒區(qū)使用六面體網(wǎng)格，在格架區(qū)使用四面體網(wǎng)格，設(shè)置3層邊界層，y+值大約為100。

圖2 外延分區(qū)子域劃分方案Fig.2 The design of the simulation scheme using overlapping domain division

所有網(wǎng)格均使用壓降進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證，如表1、2所示。對(duì)于整體計(jì)算網(wǎng)格，選擇網(wǎng)格數(shù)量為175萬(wàn)的方案。ΔP為不同數(shù)量網(wǎng)格壓降變化率，各子域也均采用壓降進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證，以子域1為例，即冷卻劑入口段，子域1選擇網(wǎng)格數(shù)量為21萬(wàn)的方案。

表1 整體網(wǎng)格獨(dú)立性驗(yàn)證Table 1 Overall grid independence verification

表2 子域1(CS)網(wǎng)格獨(dú)立性驗(yàn)證

1.2 模型驗(yàn)證

為驗(yàn)證模型與計(jì)算方法的準(zhǔn)確性，仿真結(jié)果與中國(guó)核動(dòng)力設(shè)計(jì)研究院的實(shí)驗(yàn)數(shù)據(jù)[13-15]進(jìn)行對(duì)比。選取距格架入口端50 mm和320 mm的下游處進(jìn)行橫流狀態(tài)驗(yàn)證。在對(duì)堆芯進(jìn)行CFD分析時(shí)，橫流的狀態(tài)[16]會(huì)影響冷卻劑與壁面的換熱以及流動(dòng)阻力。常用的反映橫流狀態(tài)方法：1) 在關(guān)注截面上做一取樣路徑，對(duì)路徑上檢測(cè)點(diǎn)的橫向速度進(jìn)行速度分解；2) 將檢測(cè)點(diǎn)的與取樣路徑垂直的橫向速度與該截面的平均軸流速度做商；3) 以取樣路徑上各監(jiān)測(cè)點(diǎn)的相對(duì)位置作為x坐標(biāo)，以各點(diǎn)的商作為y坐標(biāo)，繪制橫向速度規(guī)律曲線(xiàn)。橫向取樣路徑如圖3所示。穩(wěn)態(tài)工況下橫流對(duì)比情況如圖4所示，仿真結(jié)果與實(shí)驗(yàn)數(shù)據(jù)在趨勢(shì)和數(shù)值上均比較接近，從而驗(yàn)證了模型的適用性。

圖3 取樣路徑示意Fig.3 Sampling pathline for the analysis on lateral flow

圖4 AFA2G格架軸向不同截面流速分布分析Fig.4 Analysis on the lateral flow status at different vertical cross-sections for AFA2G grid

1.3 瞬態(tài)分段耦合計(jì)算方案

穩(wěn)態(tài)分段耦合方案是通過(guò)將研究區(qū)域劃分為多個(gè)子域，依次傳遞邊界條件從而實(shí)現(xiàn)全流域的CFD計(jì)算，其結(jié)果已得到驗(yàn)證[4-6]。

本文開(kāi)發(fā)了調(diào)控程序，調(diào)用Fluent，使各子域在每個(gè)時(shí)間步上交換邊界條件，順次完成全流域的瞬態(tài)CFD計(jì)算。每個(gè)子域完成一個(gè)時(shí)間步長(zhǎng)的計(jì)算后將出口參數(shù)傳遞給下一個(gè)子域，同時(shí)該子域繼續(xù)進(jìn)行瞬態(tài)計(jì)算，后續(xù)子域依次獲取前一子域的邊界條件，完成瞬態(tài)計(jì)算，直至整個(gè)流域瞬態(tài)仿真結(jié)束。

同時(shí)，該程序還具有保存各子域瞬態(tài)出口邊界信息，調(diào)控特定子域參與計(jì)算和暫停等待的功能，既可以使全流域各個(gè)子域同步完成瞬態(tài)計(jì)算，也可以先完成前序部分子域的瞬態(tài)計(jì)算，保存各時(shí)間步的出口邊界條件，以備后續(xù)子域調(diào)用，從而合理的規(guī)劃CPU和內(nèi)存的使用，實(shí)現(xiàn)多子域瞬態(tài)計(jì)算的分布并行。具體計(jì)算流程如圖5所示。

圖5 瞬態(tài)分段耦合計(jì)算方案Fig.5 Domain-divided solving technique in transient CFD

1.4 外延分區(qū)方案

在穩(wěn)態(tài)分段耦合時(shí)增加了攪混翼格架下游計(jì)算域的長(zhǎng)度保證了復(fù)雜流場(chǎng)計(jì)算結(jié)果的準(zhǔn)確性[6]。然而將適用于穩(wěn)態(tài)的子域劃分方案應(yīng)用至瞬態(tài)分段耦合中，其計(jì)算結(jié)果存在很大的誤差。例如，在子域1的入口段引入流量階躍后，監(jiān)測(cè)子域2的出口截面平均流速，并將其與整體計(jì)算方案進(jìn)行對(duì)比，如圖6所示，其計(jì)算結(jié)果誤差較大。

圖6 原子域劃分方案計(jì)算結(jié)果Fig.6 Results of the original strategy

因?yàn)橄鄬?duì)于穩(wěn)態(tài)計(jì)算，瞬態(tài)計(jì)算中攪混翼出口處不能保證為充分發(fā)展湍流，無(wú)法滿(mǎn)足自由出口的應(yīng)用條件[17-19]，即零擴(kuò)散通量假設(shè)和質(zhì)量平衡法則。因此在建模的過(guò)程中需要對(duì)其進(jìn)行修正，將數(shù)值出口盡量設(shè)置到遠(yuǎn)離物理出口的位置，從而盡可能大的減小誤差。

由此提出了外延分區(qū)的子域劃分方案，重新進(jìn)行了幾何建模以及網(wǎng)格劃分，將原子域的出口部分進(jìn)行了延長(zhǎng)，并重新繪制了網(wǎng)格。子域間數(shù)據(jù)傳遞時(shí)，將物理出口的邊界條件傳遞給下一子域的進(jìn)口。

嘗試將原子域延長(zhǎng)自身軸向長(zhǎng)度的5%、10%、20%。對(duì)比其均方根誤差和計(jì)算所需時(shí)間來(lái)選擇子域延長(zhǎng)方案。本文將整體計(jì)算方案仿真的數(shù)值作為真實(shí)值，并以此來(lái)計(jì)算外延分區(qū)耦合方案的誤差率和均方根誤差(root mean square error, RMSE)：

(4)

(5)

以格架攪混翼區(qū)為例，對(duì)比在周期性入口條件下，不同外延方案計(jì)算1 000個(gè)時(shí)間步出口截面軸流速度的RMSE和計(jì)算耗時(shí)。

經(jīng)對(duì)比，延長(zhǎng)子域軸向長(zhǎng)度的20%可以降低與整體計(jì)算方案相較的RMSE，同時(shí)其計(jì)算耗時(shí)也不會(huì)大幅增加。因此本文的計(jì)算模型均采用外延子域自身20%軸向長(zhǎng)度的計(jì)算方案。

表3 不同外延方案效率對(duì)比Table 3 Analysis of the performance of different strategies

2 瞬態(tài)計(jì)算適用性分析

2.1 流量階躍響應(yīng)特性分析

本文使用階躍型和周期性入口邊界條件對(duì)外延分區(qū)耦合瞬態(tài)CFD計(jì)算方案的有效性進(jìn)行了驗(yàn)證。穩(wěn)態(tài)時(shí)，將子域1的入口速度設(shè)置為3.76 m/s,所有子域都計(jì)算收斂后，將其入口的流量降低20%，調(diào)控所有子域進(jìn)行瞬態(tài)計(jì)算。由圖7可知，通過(guò)外延分區(qū)后，整體計(jì)算方案與分段耦合方案在瞬態(tài)仿真時(shí)，子域2出口的平均軸向流速的誤差率在整個(gè)仿真時(shí)間內(nèi)都很小，不到千分之一。相較于圖6計(jì)算結(jié)果有了很大改善。

圖7 外延分區(qū)耦合方案Fig.7 Results of the overlapping domain division

同樣采用橫向速度與軸向速度比值來(lái)分析橫流特性，取樣位置如圖2所示。當(dāng)引入流量下降的階躍入口邊界條件時(shí)，子域6(MV)出口截面的橫向速度各時(shí)刻分布規(guī)律變化如圖8所示，可以看到降低流量時(shí)，攪混翼出口橫流速度與軸流速度的比值快速增大。2種方案均能反應(yīng)這一特性，并且2種方案所反映的橫流特性也基本一致，存在的誤差從初始穩(wěn)態(tài)至瞬態(tài)計(jì)算過(guò)程中未增加，從而驗(yàn)證了外延分區(qū)耦合瞬態(tài)CFD計(jì)算方案在階躍型邊界條件下仿真橫流狀態(tài)的準(zhǔn)確性。

圖8 不同時(shí)刻攪混翼出口橫流特性分析Fig.8 Analysis on the lateral flow status at different times

2.2 橫向邊界適用性分析

在流域橫向使用周期性條件對(duì)外延分區(qū)耦合方案進(jìn)行測(cè)試，其初始穩(wěn)態(tài)與流量階躍案例一致。通過(guò)對(duì)子域1的入口流量引入周期為0.025 s的20%流量波動(dòng)，在瞬態(tài)計(jì)算中，監(jiān)測(cè)各個(gè)子域的出口截面平均流速，并與相同邊界條件的整體計(jì)算方案進(jìn)行對(duì)比，各橫向截面對(duì)比結(jié)果如圖9所示。各子域的劃分方案如圖5所示。

圖9 不同子域出口軸向流速Fig.9 Axial velocity at the outlet of different subdomains

以整體計(jì)算方案參數(shù)為真實(shí)值進(jìn)行對(duì)比，在子域1的出口，其平均流速的誤差率僅有1×10-4，效果良好，如圖9(a)所示。隨著耦合子域數(shù)量的增加，后序子域的誤差逐漸增大，最大誤差率不超過(guò)7×10-3，最大RMSE不超過(guò)2×10-2。

圖10比較了經(jīng)過(guò)一個(gè)周期后，t=2×10-5s時(shí)刻外延分區(qū)耦合方案與整體計(jì)算方案的壓降與溫度在軸向上的分布，其溫度分布在軸向上吻合良好，壓降的數(shù)值和趨勢(shì)也基本一致。

通過(guò)對(duì)比橫向截面與軸向的參數(shù)，驗(yàn)證了外延分區(qū)耦合瞬態(tài)CFD計(jì)算技術(shù)在階躍性和周期性邊界條件下的適用性與準(zhǔn)確性。在穩(wěn)態(tài)分段耦合計(jì)算中，出口邊界參數(shù)法向梯度值較小，以光棒區(qū)出口為例，其穩(wěn)態(tài)計(jì)算中，如表4所示，出口截面參數(shù)的法向梯度相較于瞬態(tài)過(guò)程小很多。因此基于物理量的局部單相化假設(shè)條件[20]，以及出口零法向梯度的充分發(fā)展邊界條件[21]，將出口設(shè)置為自由流出滿(mǎn)足數(shù)值計(jì)算條件要求，計(jì)算結(jié)果準(zhǔn)確。將出口的邊界條件條件依次傳遞給下游子域，仍然能夠保證計(jì)算結(jié)果的準(zhǔn)確性。而瞬態(tài)變化過(guò)程中，熱工參數(shù)波動(dòng)非常劇烈的，出口邊界參數(shù)法向梯度很大，因此各子域出口處不能保證為充分發(fā)展的湍流，無(wú)法滿(mǎn)足自由流出條件要求，導(dǎo)致仿真結(jié)果失真。

圖10 t=2×10-5 s時(shí)刻整體計(jì)算方案與外延分區(qū)耦合方案軸向壓降與溫度對(duì)比Fig.10 Comparison of axial parameters between overlapping domain division and overall calculation scheme at 2×10-5 seconds

表4 子域1出口處參數(shù)的法向梯度

而通過(guò)外延分區(qū)后，將數(shù)值出口設(shè)置在遠(yuǎn)離物理出口的地方，減少了不合理出口邊界條件給計(jì)算引入的誤差。

3 結(jié)論

1)本文研究了高適用且高效的瞬態(tài)CFD計(jì)算分析技術(shù)，提出的外延分區(qū)耦合瞬態(tài)CFD計(jì)算分析技術(shù)可以實(shí)現(xiàn)調(diào)控CFD程序?qū)崟r(shí)切換邊界條件，優(yōu)化子域劃分方案。

2)經(jīng)對(duì)比分析驗(yàn)證該方案在瞬態(tài)計(jì)算中橫截面速度、橫流狀態(tài)以及軸向溫度、壓降分布等均與整體計(jì)算方案吻合較好，具備工程適用性，并且降低了單次計(jì)算負(fù)荷，解決了堆芯大流域瞬態(tài)CFD計(jì)算對(duì)CPU和內(nèi)存要求過(guò)高的問(wèn)題。