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

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

引 言

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

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

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

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

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

1.1 基礎(chǔ)框架

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

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

如圖1 所示，求解T(0)(x,y,θ)可按以下三個(gè)步驟進(jìn)行：（1）首先將輸入{x,y,θ}轉(zhuǎn)換成神經(jīng)元的形式，生成初始的神經(jīng)網(wǎng)絡(luò)t0(x,y,θ)；（2）借助自動(dòng)微分算法[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)，同時(shí)也使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

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

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

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

1.2 軟邊界PINN

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

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

1.3 硬邊界PINN

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

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

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

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

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

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

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

2.1 案例模型

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

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

假設(shè)熱源保持穩(wěn)定（q=1），材料熱導(dǎo)率在{κ ∈R:0.5 ≤κ ≤1}內(nèi)變化，相應(yīng)的參數(shù)θ 變化范圍為{θ = q/κ ∈R:1 ≤θ ≤2}。以空間坐標(biāo)和控制方程的參數(shù){x, y, θ}作三維輸入，建立PINN 代理模型。訓(xùn)練集中，樣本坐標(biāo)值（x 和y）使用均勻采樣生成51×51 的網(wǎng)格，參數(shù)θ 在區(qū)間內(nèi)參數(shù)化，樣本量為10。對(duì)于軟邊界PINN，邊界上使用Monte Carlo采樣方法生成1000 個(gè)樣本。此外，設(shè)置PINN 含有3 個(gè)隱藏層，每個(gè)隱藏層含有50 個(gè)神經(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 溫度場(chǎng)預(yù)測(cè)結(jié)果

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

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

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

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

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

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

3.1 案例模型

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

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

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

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

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

其中，L 代表平板長(zhǎng)度，x∈[0, L]。基于這一定義，可構(gòu)造方程式(7)所示的擴(kuò)展函數(shù)：

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

3.2 溫度場(chǎng)預(yù)測(cè)結(jié)果

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

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

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

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

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

4 結(jié) 論

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