基于物理信息神經(jīng)網(wǎng)絡(luò)的傳熱過程物理場代理模型的構(gòu)建

陸至彬，瞿景輝，劉樺，何暢，3，張冰劍，3，陳清林，3

（1 中山大學(xué)材料科學(xué)與工程學(xué)院，廣東廣州510275； 2 中山大學(xué)化學(xué)工程與技術(shù)學(xué)院，廣東珠海519082； 3 廣東省石化過程節(jié)能工程技術(shù)研究中心，廣東廣州510275）

引 言

復(fù)雜流體在自然界和工業(yè)過程中普遍存在，精確地模擬流體基本傳熱行為在化工、熱力和航空航天等諸多學(xué)科領(lǐng)域的研究中不可或缺。通常，研究人員需要測量大量的實驗數(shù)據(jù)（例如溫度、速度、壓力等）來構(gòu)建模型[1]，以便進行參數(shù)分析、性能表征和模型優(yōu)化等工作。傳統(tǒng)的機器學(xué)習(xí)方法（深度/卷積/遞歸神經(jīng)網(wǎng)絡(luò)）[2-4]可以最大化利用上述數(shù)據(jù)蘊含的信息[5]，構(gòu)建高保真的傳熱模型以準確描述“輸入—輸出”間關(guān)聯(lián)關(guān)系，進而可用于狀態(tài)監(jiān)視、產(chǎn)品質(zhì)量控制和風險預(yù)測。但是，對于大部分具有微觀、封閉特征的復(fù)雜換熱器件（如：微尺度換熱器和高溫輻射換熱器），其內(nèi)部物理場的關(guān)鍵信息數(shù)據(jù)的測量成本高昂甚至難以測量，研究人員往往需要在數(shù)據(jù)稀疏的小數(shù)據(jù)條件下構(gòu)建模型或者做出決策，這導(dǎo)致絕大多數(shù)傳統(tǒng)機器學(xué)習(xí)方法缺乏魯棒性，無法提供收斂的保證[6-7]。另外，需要注意的是換熱器件內(nèi)部流體系統(tǒng)通常受高度非線性的Navier-Stokes 方程組控制，物理信息間受到基本物理規(guī)律的約束，盲目地以數(shù)據(jù)驅(qū)動方法學(xué)習(xí)物理場信息，可能會得到難以解釋的甚至違背物理規(guī)律的結(jié)果。

近年來，隨著深度殘差神經(jīng)網(wǎng)絡(luò)（ResNet）[8-9]以及物理信息神經(jīng)網(wǎng)絡(luò)（PINN）[7]的相繼提出，深度學(xué)習(xí)算法在求偏微分方程和認識物理規(guī)律方面取得了重大進展，突破了傳統(tǒng)機器學(xué)習(xí)受限于數(shù)據(jù)驅(qū)動和黑箱建模的局面。特別是，PINN將神經(jīng)網(wǎng)絡(luò)的預(yù)測結(jié)果約束在物理規(guī)律之內(nèi)，使機器學(xué)習(xí)方法擺脫了對實驗或模型數(shù)據(jù)的根本性依賴，在提高精度同時也大大提高了模型的可理解性。例如：Raissi等[7,10-11]將PINN 用于求解確定性的一維Burger 偏微分方程，以及求解具有中等數(shù)量訓(xùn)練數(shù)據(jù)的二維/三維偏微分方程約束反問題，可反演推導(dǎo)Navier-Stokes 方程中的不確定參數(shù)。 在此基礎(chǔ)上Karniadakis等[12-14]成功將PINN應(yīng)用于反應(yīng)擴散方程和歐拉方程的求解。Wang 等[15]利用PINN 來預(yù)測基于雷諾時均方程模型的流體雷諾應(yīng)力的差異，并取得了良好的預(yù)測效果。類似地，PINN還被應(yīng)用于學(xué)習(xí)達西流中的本構(gòu)關(guān)系[16]、材料的原子模擬[17]、多物理場亞表面?zhèn)鬏攩栴}[18]，以及不確定性量化問題[19-21]。

PINN耦合了換熱器件內(nèi)部的物理規(guī)律，因而僅僅需要實際測得的邊界條件即可準確預(yù)測整個物理場的信息。但仍需注意到，在訓(xùn)練PINN 時使用不同的邊界條件設(shè)置方法，對訓(xùn)練結(jié)果有著顯著的影響。PINN 的損失函數(shù)包含域損失項和邊界損失項，通常損失函數(shù)通過“軟邊界”的方式設(shè)置邊界損失項，由懲罰系數(shù)來控制邊界損失項和域損失項的權(quán)重，以加速優(yōu)化問題的收斂。但是，懲罰系數(shù)的選擇往往依靠經(jīng)驗來調(diào)整，不當?shù)膽土P系數(shù)容易導(dǎo)致求解結(jié)果出現(xiàn)畸形解。針對上述問題，Sun等[22]引入距離函數(shù)，使用“硬邊界”的方式設(shè)置邊界條件，將邊界損失和域損失統(tǒng)一在同一項中，避免了懲罰系數(shù)的使用，并提升了預(yù)測精度。

PINN提供了求解偏微分方程的新途徑，同時其耦合物理定律的特性呈現(xiàn)出其直接作為代理模型的潛力[22-25]。本工作基于PINN 理論，以二維空間坐標和控制方程的參數(shù)構(gòu)成三維輸入，預(yù)測了物理場的云圖分布，并與商業(yè)軟件模擬結(jié)果進行對比驗證。在這一過程中，首先通過簡單的具有內(nèi)熱源的二維穩(wěn)態(tài)導(dǎo)熱案例分析對比軟邊界PINN 與硬邊界PINN 的性能，接著將硬邊界PINN 應(yīng)用到板式換熱器的流動傳熱案例當中，討論其在構(gòu)建復(fù)雜物理場代理模型方面的性能。

1 基于物理信息的神經(jīng)網(wǎng)絡(luò)

1.1 基礎(chǔ)框架

物理過程均可以用多個偏微分方程進行描述，求解得到的原函數(shù)，無論其形式是解析解還是數(shù)值解，均包含了大量的物理場信息。對原函數(shù)T(0)(x,y,θ)，其二維偏微分方程g(T)如下：

其中，T(n)(x,y,θ)表示T(0)(x,y,θ)的n階導(dǎo)數(shù)項，n={0,1,2}，C為常數(shù)，x、y為坐標，θ為參數(shù)。

如圖1 所示，求解T(0)(x,y,θ)可按以下三個步驟進行：（1）首先將輸入{x,y,θ}轉(zhuǎn)換成神經(jīng)元的形式，生成初始的神經(jīng)網(wǎng)絡(luò)t0(x,y,θ)；（2）借助自動微分算法[26]求t0(x, y, θ)的梯度，構(gòu)成偏微分方程g(t0)，并以g(t0)與g(T)之間的偏差模量作為損失函數(shù)；（3）最后利用梯度下降算法[27]，調(diào)整神經(jīng)網(wǎng)絡(luò)的權(quán)重w 和b，使g(t0)逼近g(T)，同時也使t(x, y, θ)逼近T(0)(x, y, θ)。多次迭代訓(xùn)練后得到的t(x, y, θ)即是T(0)(x, y, θ)的PINN代理模型。

圖1 PINN結(jié)構(gòu)示意圖Fig.1 PINN structure diagram

其中，上角標l 代表神經(jīng)網(wǎng)絡(luò)的層序號，下角標i，j代表神經(jīng)元索引，例如：和分別為第l-1 層中第i 個神經(jīng)元到第l 層中第j 個神經(jīng)元的權(quán)重和偏置；a=σ(z)為激活函數(shù)[28-29]，可選擇雙曲正切函數(shù)或sigmoid函數(shù)：

其中，z 和a 分別為神經(jīng)元的輸入值和輸出值。最后一層第L層為輸出層，輸出擬合函數(shù)t：

由此，基于式(2)給定任意的輸入值即可快速計算出t(x,y,θ)。

1.2 軟邊界PINN

神經(jīng)網(wǎng)絡(luò)各層之間的權(quán)重、偏置可以通過最小化均方誤差損失來學(xué)習(xí)，定義損失函數(shù)為

對于給定的參數(shù)，{xb,yb,T(xb,yb)}為用于訓(xùn)練的邊界點數(shù)據(jù)，{xf,yf}為區(qū)域內(nèi)部用于訓(xùn)練的配置點數(shù)據(jù)，λ 是懲罰系數(shù)。損失函數(shù)Loss 第一項對應(yīng)于邊界數(shù)據(jù)；而第二項為充當懲罰作用的正則化項，保證了神經(jīng)網(wǎng)絡(luò)在有限的配置點集合上強制執(zhí)行偏微分方程式(1)所施加的物理定律約束。例如，在不可壓縮流體動力學(xué)問題中，任何違反質(zhì)量守恒原理的流動解都會被丟棄，從而將解空間約束在物理框架下[7,30]。

1.3 硬邊界PINN

軟邊界PINN 對邊界的處理采用的是有約束優(yōu)化，其局限性在于損失函數(shù)的優(yōu)化性能取決于每個項的相對大小。因此，如何確定λ 值使得損失函數(shù)每一項在優(yōu)化過程中處于相同的數(shù)量級是很有挑戰(zhàn)性的工作。Sun 等[22,31]提出一種靈活處理邊界條件的方法，即構(gòu)造一個滿足邊界條件特解的函數(shù)，具體的表達式如下：

其中，G(x,y,θ)為邊界函數(shù)，是邊界數(shù)據(jù)的光滑擴展函數(shù)，D(x,y,θ)為距離函數(shù)，是與內(nèi)部點到邊界最小距離有關(guān)的光滑函數(shù)。

對于Dirichlet 邊界條件，邊界上各點的邊界值為T(xb, yb, θ)，坐標為(xb, yb)。此時，由于邊界上的D(xb, yb, θ)=0，因此有h(xb, yb, θ)=G(xb, yb, θ) = T(xb, yb,θ)。而對于不在邊界的內(nèi)部點(xf,yf)，其D(xf,yf,θ)取值不為0，則h(x,y,θ)受到邊界函數(shù)和距離函數(shù)共同控制。通過將Dirichlet邊界條件耦合在所構(gòu)造的擴展函數(shù)h(x, y, θ)當中，損失函數(shù)被進一步簡化為只包括偏導(dǎo)數(shù)的形式：

其中，{xt,yt}為用于訓(xùn)練的所有點數(shù)據(jù)（包括邊界）。在上述情況下，距離函數(shù)與DNN 輸出的乘積D(x,y,θ)aL(x,y,θ;w,b)在數(shù)學(xué)上相當于此偏微分方程的通解，而邊界函數(shù)則包含了偏微分方程的特解。這種“通解+特解”的數(shù)學(xué)結(jié)構(gòu)，賦予求解結(jié)果自動滿足特解的先行條件，加速PINN 向合理梯度方向收斂。

值得注意的是，上述推導(dǎo)是基于Dirichlet 邊界條件，對于梯度形式的Neumann 邊界（h(1)(x, y, θ)），則將邊界梯度直接引入損失函數(shù)當中，以硬邊界PINN損失函數(shù)為例：

式中，{xB, yB}為用于訓(xùn)練的Neumann 邊界點數(shù)據(jù)，C0為Neumann邊界處的梯度值。

2 案例一：具有內(nèi)熱源的二維穩(wěn)態(tài)導(dǎo)熱方程

2.1 案例模型

本節(jié)首先以一個簡單的案例探究軟、硬兩種邊界條件設(shè)置方式下，PINN預(yù)測結(jié)果的差異。以具有內(nèi)熱源的二維穩(wěn)態(tài)導(dǎo)熱方程為例，其表示式如下

其中，q為熱源，κ為熱導(dǎo)率。案例選取邊長為1的正方形區(qū)域作為研究對象，以方形中心為原點建立笛卡兒坐標系。

假設(shè)熱源保持穩(wěn)定（q=1），材料熱導(dǎo)率在{κ ∈R:0.5 ≤κ ≤1}內(nèi)變化，相應(yīng)的參數(shù)θ 變化范圍為{θ = q/κ ∈R:1 ≤θ ≤2}。以空間坐標和控制方程的參數(shù){x, y, θ}作三維輸入，建立PINN 代理模型。訓(xùn)練集中，樣本坐標值（x 和y）使用均勻采樣生成51×51 的網(wǎng)格，參數(shù)θ 在區(qū)間內(nèi)參數(shù)化，樣本量為10。對于軟邊界PINN，邊界上使用Monte Carlo采樣方法生成1000 個樣本。此外，設(shè)置PINN 含有3 個隱藏層，每個隱藏層含有50 個神經(jīng)元，激活函數(shù)為sigmoid函數(shù)，見式(4)。

軟邊界PINN 的損失函數(shù)設(shè)置為式(6)的形式，且方形四周的邊界溫度恒為1。而硬邊界PINN 的損失函數(shù)設(shè)置為式(8)的形式，其中邊界函數(shù)G(x,y,θ)和距離函數(shù)D(x,y,θ)分別如下：

2.2 溫度場預(yù)測結(jié)果

本工作所涉及的采樣方案和神經(jīng)網(wǎng)絡(luò)訓(xùn)練均基于Python3.7.7 和Tensorflow1.14 環(huán)境，在Spyder 中編譯完成。有限元模擬均使用COMSOL 5.5。計算機配置為Intel i7-8550 CPU。訓(xùn)練后使用θ=1.5 時的溫度場與COMSOL 模擬結(jié)果對比，以驗證PINN的預(yù)測效果。如圖2所示，兩種PINN 方法構(gòu)建的代理模型預(yù)測結(jié)果均與模擬結(jié)果非常接近，均準確地捕捉到了溫度的變化趨勢。結(jié)合相對誤差場圖（軟邊界PINN 見圖3，硬邊界PINN 見圖4），軟邊界PINN 和硬邊界PINN 最大相對誤差分別為0.60%和0.015%，再次表明兩種PINN 方法訓(xùn)練得到的代理模型都具有很高的精度，但硬邊界PINN 的相對誤差要更低。

圖2 PINN代理模型和COMSOL模擬對溫度場的預(yù)測結(jié)果對比Fig.2 Comparsion of temperature fields predicted by soft PINN,hard PINN,and COMSOL simulation

圖3 軟邊界PINN溫度等值線的預(yù)測結(jié)果對比Fig.3 Comparsion of temperature contours predicted by soft PINN

圖4 硬邊界PINN溫度等值線的預(yù)測結(jié)果對比Fig.4 Comparsion of temperature contours predicted by hard PINN

在分別對比軟邊界和硬邊界與COMSOL 模擬結(jié)果的等值線圖可見，中心區(qū)域附近，使用這兩種邊界設(shè)置方法的預(yù)測結(jié)果均與模擬結(jié)果基本重合。但仍需注意到，軟邊界PINN 在四個角落區(qū)域的預(yù)測效果較差，而硬邊界PINN 在四個角落區(qū)域的預(yù)測效果依然良好。造成這種現(xiàn)象的原因可能是軟邊界PINN 用于訓(xùn)練的邊界點數(shù)據(jù)量較少。為此，將邊界點的樣本量增加到4000后重新訓(xùn)練，預(yù)測精度提升13.3%，但是角落區(qū)域預(yù)測誤差較大的問題依然存在。此外，硬邊界PINN 的計算時間更小。硬邊界PINN 所用的訓(xùn)練時間為7 min，遠低于軟邊界PINN 所用的訓(xùn)練時間39 min（1000 個邊界樣本點）和55 min（4000 個邊界樣本點）。因此可以認為，硬邊界PINN要優(yōu)于軟邊界PINN。

3 案例二：平板間的二維穩(wěn)態(tài)對流傳熱方程

3.1 案例模型

在換熱器件內(nèi)部，動量傳遞、熱量傳遞和質(zhì)量傳遞往往聯(lián)系密切，因此研究多物理場對指導(dǎo)實際換熱過程具有重要的意義。例如在板式換熱器中，平板間的溫度傳遞不僅由熱傳導(dǎo)引起，還受到流體流動影響，本節(jié)以此為例利用PINN 來預(yù)測二維穩(wěn)態(tài)對流傳熱方程。所使用模型幾何如圖5 所示，換熱工質(zhì)在兩塊平板間流動，并存在熱傳導(dǎo)和熱對流。假設(shè)壁面無滑移，且兩平板間的不可壓縮穩(wěn)態(tài)流動已經(jīng)充分發(fā)展，則垂直于平板方向的速度可由Navier-Stokes方程推導(dǎo)得出拋物線形的解析解：

其中，u 表示速度場的x 分量，umax為x 分量的最大速度，R為平板壁面到中軸線的垂直距離。

對于有內(nèi)熱源的穩(wěn)態(tài)層流傳熱，忽略黏性耗散后，其能量守恒式可簡化為

其中，ρ 和Cp分別表示流體的密度和比定壓熱容，顯然該方程是在方程式(10)的基礎(chǔ)上引入了第一項對流項。假設(shè)平板隔熱性能良好，與外界環(huán)境絕熱，則在平板上下壁面需設(shè)置Neumann 邊界條件：

本案例采用硬邊界PINN 的方法建立代理模型，假設(shè)進口溫度的Dirichlet邊界Tin已知，則可構(gòu)造以下滿足條件的邊界函數(shù)和距離函數(shù)：

其中，L 代表平板長度，x∈[0, L]?；谶@一定義，可構(gòu)造方程式(7)所示的擴展函數(shù)：

令空氣為板間流體，物性參數(shù)為：ρ=1.239 kg/m3，Cp=1,005 J/(kg·K)，κ=0.024 W/(m2·K)。進口溫度Tin=30℃，平板間的通道尺寸為L=1 m，R=0.05 m。以{x,y,umax}作為輸入變量的三個維度。其中，x，y為兩平板間均勻生成的50×100 網(wǎng)格的坐標。最大速度umax在區(qū)間{0.01～0.03 m/s}內(nèi)參數(shù)化，取10 個樣本。設(shè)置PINN 隱藏層4 層，每層含有50 個神經(jīng)元，使用sigmoid函數(shù)[式(4)]作為激活函數(shù)。

3.2 溫度場預(yù)測結(jié)果

訓(xùn)練代理模型所用的時間為54 min，訓(xùn)練后損失函數(shù)值為2.6×10-4。使用umax=0.02 m/s時的溫度場驗證預(yù)測效果，如圖6 所示。由圖可知硬邊界PINN構(gòu)建的代理模型所預(yù)測的流體流動過程中的溫度變化趨勢與實際模擬結(jié)果是一致的。值得注意的是，在相同的橫坐標下，流體在中間的溫度要低于兩側(cè)，而中間的流速高于兩側(cè)，由此說明流動加速了流體的熱量交換，流速越大熱交換越明顯。此外，溫度場表現(xiàn)出與速度場類似的對稱性，再次說明穩(wěn)態(tài)熱量傳遞與動量傳遞密切相關(guān)。圖6的溫度場絕對誤差云圖表明，平板區(qū)域內(nèi)的溫度預(yù)測誤差基本接近于0，僅在四個角落處誤差略微變大，但依然在10-3數(shù)量級內(nèi)，最大絕對誤差為0.0038℃，再次凸顯了硬邊界PINN 構(gòu)建的代理模型精度高的優(yōu)勢。

圖5 平板的幾何結(jié)構(gòu)示意圖Fig.5 Geometric diagram of the flat plate

圖6 COMSOL模擬溫度場與硬邊界PINN預(yù)測溫度場對比Fig.6 Comparsion of temperature fields predicted by hard PINN and COMSOL simulation

圖7 硬邊界PINN代理模型（紅）與COMSOL模擬（黑）預(yù)測溫度等值線結(jié)果對比Fig.7 Comparsion of temperature contours predicted by hard PINN(red line)and COMSOL simulation(black line)

硬邊界PINN 預(yù)測的溫度等值線與COMSOL 模擬的溫度等值線在圖7 中給出。顯然，等值線在平板間中軸線處的延伸趨勢比壁面處的更大，這一變化直觀地呈現(xiàn)了流體流動對溫度分布的影響。需要指出的是，PINN預(yù)測與COMSOL模擬結(jié)果的差距沿流體流動方向擴大。在進口處兩者的溫度等值線基本重合，而在遠離進口處兩者的差異逐漸增加，并且在中軸線附近的差距最大。這是由于速度場的非線性造成溫度場非線性的增加，這會使代理模型中產(chǎn)生更大的偏差，這也是未來PINN 需要面對的主要挑戰(zhàn)。

4 結(jié) 論

本工作將PINN 模型引入到換熱物理場的代理模型的構(gòu)建中，預(yù)測給定參數(shù)下的溫度場分布，并與商業(yè)軟件的模擬結(jié)果進行驗證。通過簡單的具有內(nèi)熱源的二維穩(wěn)態(tài)導(dǎo)熱案例分析，軟邊界PINN和硬邊界PINN 最大相對誤差分別為0.60% 和0.015%，結(jié)果表明對邊界條件進行耦合的硬邊界PINN 在建立代理模型時更勝一籌。而將硬邊界PINN應(yīng)用到板式換熱器的流動傳熱案例當中，所建立的代理模型的預(yù)測結(jié)果與模擬結(jié)果吻合良好，溫度場最大絕對誤差為0.0038℃，再次證明了硬邊界PINN構(gòu)建高精度代理模型的能力，并為換熱器的開發(fā)提供了良好的工具。