基于神經(jīng)網(wǎng)絡(luò)的CFD粗網(wǎng)格模擬優(yōu)化方法研究

靳 爽 劉曉晶 程 旭

（上海交通大學(xué) 核科學(xué)與工程學(xué)院 上海200240）

在反應(yīng)堆計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）數(shù)值模擬中，計(jì)算精度和計(jì)算時(shí)間是一對突出的矛盾。為了能夠較好地兼顧計(jì)算效率和計(jì)算精度，實(shí)現(xiàn)在較高計(jì)算效率下取得滿意的計(jì)算結(jié)果，研究人員嘗試了多種新技術(shù)方法，比如并行計(jì)算、對網(wǎng)格劃分方式進(jìn)行優(yōu)化［1-2］等。隨著大數(shù)據(jù)和人工智能等領(lǐng)域的發(fā)展，機(jī)器學(xué)習(xí)也成為解決上述問題的一種可行路徑。一方面，大量先前已有的數(shù)值模擬結(jié)果、理論研究成果和實(shí)驗(yàn)結(jié)果積累為數(shù)據(jù)訓(xùn)練提供了廣泛的訓(xùn)練資源；另一方面，隨著反應(yīng)堆熱工水力研究和應(yīng)用的不斷發(fā)展，對“物理過程相似、工況參數(shù)不同的大批量多組工況”的CFD數(shù)值模擬需求日益凸顯，例如研究多種不同入口速度下的棒束流動傳熱CFD結(jié)果、多種搖擺條件（搖擺俯角、搖擺周期的各類組合）下的反應(yīng)堆CFD計(jì)算結(jié)果等，在該類問題中，基于神經(jīng)網(wǎng)絡(luò)的CFD粗網(wǎng)格模擬優(yōu)化方法對于“實(shí)現(xiàn)在較高計(jì)算效率下實(shí)現(xiàn)較高精度的數(shù)值模擬”具有顯著效果。

根據(jù)美國愛達(dá)荷國家實(shí)驗(yàn)室于2018年發(fā)布的報(bào)告［3］，在數(shù)據(jù)驅(qū)動的機(jī)器學(xué)習(xí)架構(gòu)中可以使低分辨率的模型計(jì)算出高分辨率模型下的模擬結(jié)果，并對絕熱腔室中射入空氣湍流模擬結(jié)果進(jìn)行了優(yōu)化。目前，通過機(jī)器學(xué)習(xí)針對反應(yīng)堆熱工水力數(shù)值模擬進(jìn)行優(yōu)化的案例研究還較少，本文擬應(yīng)用機(jī)器學(xué)習(xí)方法，針對粗網(wǎng)格下5×5棒束出口截面冷卻劑子通道溫度的CFD結(jié)果進(jìn)行優(yōu)化，以期為機(jī)器學(xué)習(xí)方法在反應(yīng)堆熱工領(lǐng)域?qū)崿F(xiàn)較高計(jì)算效率和較高計(jì)算精度的數(shù)值模擬應(yīng)用提供參考。

1 數(shù)值模擬優(yōu)化方法

1.1 棒束模型

本文計(jì)算的棒束模型是截短型組件的5×5棒束柵元，其主要參數(shù)參考一體化模塊小型反應(yīng)堆常用的截短型組件參數(shù)而設(shè)定［4-7］，如表1所示。

針對本文計(jì)算的棒束模型，可以劃分為36個(gè)子通道，從左上到右下依次編號為子通道1（sub 1）、子通道2（sub 2）、……、子通道36（sub 36），如圖1所示。

表1 5×5棒束柵元主要參數(shù)Table 1 Main parameters of 5×5 rod bundle

圖1 5×5棒束柵元子通道劃分方式Fig.1 Subchannel division method of 5×5 rod bundle

1.2 問題描述

本文旨在對不同入口速度工況中出口截面上36個(gè)子通道各自的平均溫度值進(jìn)行修正，粗網(wǎng)格和細(xì)網(wǎng)格均采用STAR-CCM中的結(jié)構(gòu)化定向網(wǎng)格［8］，基礎(chǔ)網(wǎng)格為多邊形網(wǎng)格，對于邊界層處的網(wǎng)格，均使用棱柱層網(wǎng)格生成，粗網(wǎng)格和細(xì)網(wǎng)格下均設(shè)置棱柱層層數(shù)為三層，棱柱層增長率取默認(rèn)值1.5，棱柱層總厚度取默認(rèn)值即網(wǎng)格尺寸基數(shù)的33.3%。其網(wǎng)格參數(shù)的具體設(shè)置如表2所示。

細(xì)網(wǎng)格參數(shù)的選取達(dá)到網(wǎng)格無關(guān)性的要求，基礎(chǔ)尺寸為0.5 mm、0.6 mm、0.7 mm網(wǎng)格下的子通道1、子通道2、子通道16（分別為角通道、邊通道、中心通道）的溫度計(jì)算結(jié)果如表3所示。

本文根據(jù)入口速度的不同，將20組工況作為訓(xùn)練組用于機(jī)器學(xué)習(xí)，另外5組工況作為驗(yàn)證組以檢驗(yàn)修正的效果。其中，5組驗(yàn)證組工況中既包括相對于訓(xùn)練組的內(nèi)含參數(shù)工況，也包括外延參數(shù)工況。訓(xùn)練組工況及驗(yàn)證組工況的入口速度設(shè)定情況如表4所示。鑒于本文計(jì)算的物理過程并不復(fù)雜，在本文涉及的計(jì)算工況中，所選用的湍流模型均為可實(shí)現(xiàn)的k-ε（Realizablek-ε）模型。

表2 粗網(wǎng)格參數(shù)和細(xì)網(wǎng)格參數(shù)Table 2 Parameters of coarse grid and fine grid

表3 網(wǎng)格無關(guān)性驗(yàn)證結(jié)果Table 3 Results of grid independence verification

表4 訓(xùn)練組和驗(yàn)證組工況入口速度設(shè)定情況Table 4 Inlet speed of training group and verification group working conditions

1.3 基于機(jī)器學(xué)習(xí)的優(yōu)化方法

數(shù)值模擬的誤差來源主要包括以下三類：第一類是由于物理簡化和（或）數(shù)學(xué)近似導(dǎo)致的模型誤差；第二類是由于網(wǎng)格的存在，導(dǎo)致守恒方程或源項(xiàng)在時(shí)空平均方法下產(chǎn)生信息丟失所引起的網(wǎng)格誤差；第三類是由于迭代收斂、算法選擇等原因引起的其他誤差。

對于第一類誤差，針對流體和壁面間傳熱、質(zhì)量交換、動量交換等現(xiàn)象存在多種關(guān)系模型進(jìn)行模擬，與網(wǎng)格尺寸相關(guān)的特征長度是這些模型中的一個(gè)重要參數(shù)。對于第二類誤差，同樣由網(wǎng)格尺寸決定。而第三類誤差相較前兩類誤差對計(jì)算結(jié)果的影響要小得多。因此，數(shù)值模擬的誤差與網(wǎng)格尺寸密切相關(guān)，且實(shí)際分析時(shí)往往無法區(qū)分模型誤差和網(wǎng)格誤差，應(yīng)將其作為一個(gè)整體來考慮。

我們可以運(yùn)用模式識別和統(tǒng)計(jì)分析等機(jī)器學(xué)習(xí)方法，從數(shù)據(jù)中直接獲得容易被人忽視的信息，試圖尋找到包括網(wǎng)格信息、模型信息、物理信息在內(nèi)的參數(shù)與誤差之間的定量關(guān)系。在此基礎(chǔ)上，運(yùn)用該定量關(guān)系對粗網(wǎng)格模擬結(jié)果進(jìn)行修正，實(shí)現(xiàn)滿足計(jì)算要求的快速數(shù)值模擬。

根據(jù)問題的需要，我們可以將模擬對象的幾何參數(shù)、模擬工況的物理參數(shù)及數(shù)值計(jì)算的網(wǎng)格參數(shù)納入機(jī)器學(xué)習(xí)的輸入變量中。針對本文研究的工況，每組工況均采用同一套粗網(wǎng)格和同一套細(xì)網(wǎng)格，各組工況在進(jìn)行CFD計(jì)算時(shí)選用的湍流模型、傳熱模型一致，僅有冷卻劑入口速度設(shè)定不同。因此，本文通過機(jī)器學(xué)習(xí)需要得出的是入口速度、粗網(wǎng)格下出口截面各個(gè)子通道溫度、細(xì)網(wǎng)格下出口截面各個(gè)子通道溫度之間的隱含關(guān)系。

根據(jù)Chang等［9］的文章，目前數(shù)據(jù)驅(qū)動的熱工水力模擬機(jī)器學(xué)習(xí)架構(gòu)可以依據(jù)數(shù)值模擬是否包含偏微分方程、偏微分方程的形式是否已知、偏微分方程本身是否作為機(jī)器學(xué)習(xí)對象、模型的建立是否需要在不同尺度上分別進(jìn)行等四項(xiàng)標(biāo)準(zhǔn)分成5種類型，本文的研究屬于上述五種架構(gòu)類型中的第二類，其主要步驟如圖2所示。

圖2 基于機(jī)器學(xué)習(xí)的數(shù)值模擬優(yōu)化框架[9]Fig.2 Optimization framework for numerical simulation based on machine learning[9]

本文中，定義粗網(wǎng)格下子通道i（subi）的平均溫度為TLi，作為待修正值；細(xì)網(wǎng)格下子通道i的溫度記為THi，作為精確值（修正目標(biāo)值）。定義誤差函數(shù)如下：

根據(jù)前述的優(yōu)化原理，將誤差函數(shù)E i設(shè)置為如下的形式：

式中：V為冷卻劑入口速度，子通道編號i=1，2，3，...，36，誤差函數(shù)中的具體函數(shù)形式由機(jī)器學(xué)習(xí)方法得到。將粗網(wǎng)格計(jì)算結(jié)果按照誤差函數(shù)修正后的子通道i的溫度記為TMi，即有：

在機(jī)器學(xué)習(xí)方法的選擇上，目前主要使用深度學(xué)習(xí)的方法。原則上，任何兩層以上的神經(jīng)網(wǎng)絡(luò)都可以被認(rèn)為是深度學(xué)習(xí)。Hornik等［10］在文章中論證了多層神經(jīng)網(wǎng)絡(luò)是一種通用方法，可以捕捉任何可度量信息的特性。神經(jīng)網(wǎng)絡(luò)包括卷積神經(jīng)網(wǎng)絡(luò)和前饋神經(jīng)網(wǎng)絡(luò)兩種方式。其中，卷積神經(jīng)網(wǎng)絡(luò)通常具有更好的預(yù)測能力，但所需的數(shù)據(jù)訓(xùn)練量也更多。鑒于本文研究的工況相對簡單，本文選用所需數(shù)據(jù)量較少的前饋神經(jīng)網(wǎng)絡(luò)。對于神經(jīng)網(wǎng)絡(luò)層數(shù)和神經(jīng)元數(shù)量的選取，海軍工程大學(xué)曹植珺等［11］的綜述文章中指出目前缺少很好的參數(shù)調(diào)節(jié)的手段，包括訓(xùn)練神經(jīng)網(wǎng)絡(luò)時(shí)隱含層數(shù)、隱含節(jié)點(diǎn)數(shù)、輸入節(jié)點(diǎn)數(shù)等基本依靠經(jīng)驗(yàn)。對于不復(fù)雜的關(guān)系結(jié)構(gòu)，神經(jīng)網(wǎng)絡(luò)的層數(shù)通常為2～3層，神經(jīng)元數(shù)量不超過10個(gè)，本文選擇的神經(jīng)網(wǎng)絡(luò)層數(shù)為兩層，神經(jīng)元數(shù)量為4個(gè)，其結(jié)構(gòu)如圖3所示。

第一層為輸入變量層，包括工況冷卻劑入口速度V以及粗網(wǎng)格下出口截面子通道溫度TLi(i=1，2，3，...，36)共計(jì)37個(gè)輸入變量；第二層為神經(jīng)元層，共包括4個(gè)神經(jīng)元N1，…，N4；第三層為輸出變量層，包括粗網(wǎng)格與細(xì)網(wǎng)格下出口截面各個(gè)子通道溫度的誤差預(yù)估值E1，E2，…，E36，共計(jì)36個(gè)輸出變量。根據(jù)神經(jīng)網(wǎng)絡(luò)的基本原理，各層元素之間的關(guān)系如下：

圖3 前饋神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)架構(gòu)Fig.3 Framework of feedforward neural network(FNN)

式中：f為神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)的函數(shù)形式，為表達(dá)方便，記：

本文采用均方根誤差來評估優(yōu)化效果，分別定義未經(jīng)優(yōu)化的粗網(wǎng)格下出口截面36個(gè)子通道溫度與精確值（細(xì)網(wǎng)格下相應(yīng)計(jì)算結(jié)果）之間的均方根誤差為RMSEL，經(jīng)過優(yōu)化的粗網(wǎng)格下出口截面36個(gè)子通道溫度與精確值（細(xì)網(wǎng)格計(jì)算結(jié)果）之間的均方根誤差為RMSEM，其形式如下：

通過比較RMSEM與RMSEL的大小差異來評估優(yōu)化效果，如果0～RMSEM?RMSEL，則可認(rèn)為達(dá)到預(yù)期優(yōu)化效果。

2 優(yōu)化結(jié)果

2.1 計(jì)算精度

使用機(jī)器學(xué)習(xí)得到的誤差函數(shù)對訓(xùn)練組和驗(yàn)證組的粗網(wǎng)格結(jié)果進(jìn)行優(yōu)化前后，出口截面36個(gè)子通道溫度與細(xì)網(wǎng)格下模擬結(jié)果的均方根誤差如圖4所示。

其中，5組驗(yàn)證組工況的粗網(wǎng)格結(jié)果優(yōu)化前后的均方根誤差如表5所示。

圖4 粗網(wǎng)格結(jié)果優(yōu)化前后的均方根誤差Fig.4 Root mean square error(RMSE)before and after optimization

表5 驗(yàn)證組工況粗網(wǎng)格結(jié)果優(yōu)化前后的均方根誤差Table 5 Root mean square error before and after optimization of the verification group results

通過圖4和表5可以看出，訓(xùn)練組和驗(yàn)證組粗網(wǎng)格下修正優(yōu)化后的結(jié)果相比于優(yōu)化前，與細(xì)網(wǎng)格下結(jié)果的均方根誤差大大降低，5組檢驗(yàn)組工況出口截面36個(gè)子通道溫度均方根誤差平均值從1.16 K降低為9.24×10-3K，滿足0～RMSEM?RMSEL，達(dá)到預(yù)期的優(yōu)化效果。

以入口速度為2.25 m?s-1的工況為例，其細(xì)網(wǎng)格下、粗網(wǎng)格下結(jié)果修正前后36個(gè)子通道的溫度如圖5所示。通過圖5可以看出，修正優(yōu)化后的各個(gè)子通道溫度值相比于優(yōu)化前，明顯更加接近細(xì)網(wǎng)格下的結(jié)果；且對于邊通道和角通道，優(yōu)化提升的效果更加明顯。

2.2 計(jì)算效率

使用本文的方法可以對粗網(wǎng)格下的計(jì)算結(jié)果進(jìn)行修正，取得了良好的優(yōu)化效果，這也使得CFD計(jì)算效率大幅提升。

對于粗網(wǎng)格、細(xì)網(wǎng)格工況，在運(yùn)算時(shí)均以連續(xù)20個(gè)迭代步內(nèi)出口平均溫度變化不超過0.01℃為判斷計(jì)算達(dá)到收斂的條件。本文涉及的表4中共計(jì)25個(gè)工況在粗網(wǎng)格、細(xì)網(wǎng)格下進(jìn)行CFD運(yùn)算所需要的總CPU系統(tǒng)時(shí)間如表6所示。

圖5 V=2.25 m?s-1工況下粗網(wǎng)格結(jié)果優(yōu)化前后子通道溫度值與細(xì)網(wǎng)格下結(jié)果的對比Fig.5 Comparison of subchannel temperature before and after optimization under the condition of V=2.25 m?s-1

從表6可以看出，不同工況下粗網(wǎng)格計(jì)算所需的總CPU系統(tǒng)時(shí)間遠(yuǎn)遠(yuǎn)小于細(xì)網(wǎng)格下所需時(shí)間，上述25個(gè)工況粗網(wǎng)格下總CPU系統(tǒng)時(shí)間的平均值約為細(xì)網(wǎng)格下平均值的2.54%。就本文研究的問題而言，若需要得到入口速度從1.65～2.85 m·s-1間等距100個(gè)工況的計(jì)算結(jié)果，傳統(tǒng)方法需運(yùn)行100個(gè)細(xì)網(wǎng)格工況，而基于神經(jīng)網(wǎng)絡(luò)的粗網(wǎng)格優(yōu)化方法下，僅需計(jì)算25個(gè)細(xì)網(wǎng)格工況和100個(gè)粗網(wǎng)格工況（以表6的CPU系統(tǒng)時(shí)間平均值計(jì)算，約相當(dāng)于三個(gè)細(xì)網(wǎng)格工況的計(jì)算量），總體的計(jì)算效率得到大幅提高。

3 結(jié)語

本研究利用前饋神經(jīng)網(wǎng)絡(luò)的機(jī)器學(xué)習(xí)方法，對不同入口速度工況下出口截面子通道冷卻劑溫度粗網(wǎng)格下的CFD數(shù)值模擬結(jié)果進(jìn)行了優(yōu)化，得到如下結(jié)論：

1）采用神經(jīng)網(wǎng)絡(luò)方法對本文工況粗網(wǎng)格下的結(jié)果進(jìn)行修正后，5組檢驗(yàn)組工況出口截面36個(gè)子通道溫度的均方根誤差平均值從1.16 K降低為9.24×10-3K，誤差大大減小，優(yōu)化效果顯著。

2）在本文研究的問題中，粗網(wǎng)格下總CPU系統(tǒng)時(shí)間的平均值約為細(xì)網(wǎng)格下平均值的2.54%。通過神經(jīng)網(wǎng)絡(luò)對粗網(wǎng)格結(jié)果優(yōu)化的方式，對于批量化的CFD計(jì)算或者已有較為充分的訓(xùn)練數(shù)據(jù)的CFD計(jì)算，其計(jì)算效率將大幅提高。

3）就本文算例來看，基于前饋神經(jīng)網(wǎng)絡(luò)的機(jī)器學(xué)習(xí)框架對于內(nèi)部預(yù)測（即驗(yàn)證工況的速度值在訓(xùn)練集速度值范圍內(nèi)）的優(yōu)化效果好于外部預(yù)測。但由于本文算例較為局限，這一結(jié)論還需要更多算例結(jié)果的支持驗(yàn)證。

本文的研究為利用機(jī)器學(xué)習(xí)框架對粗網(wǎng)格數(shù)值模擬進(jìn)行優(yōu)化，從而兼顧計(jì)算效率與計(jì)算精度提供了可行參照。未來，可以在反應(yīng)堆熱工水力領(lǐng)域的更多問題特別是多耦合或非穩(wěn)態(tài)問題上進(jìn)行更廣泛的驗(yàn)證。

表6 25個(gè)工況在粗網(wǎng)格、細(xì)網(wǎng)格下運(yùn)行所需的總CPU系統(tǒng)時(shí)間Table 6 The total CPU system time required for 25 operating conditions to run under coarse and fine grids