基于物理信息神經網絡的傳熱過程物理場代理模型的構建

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

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

引 言

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

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

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

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

1 基于物理信息的神經網絡

1.1 基礎框架

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

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

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

圖1 PINN結構示意圖Fig.1 PINN structure diagram

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

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

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

1.2 軟邊界PINN

神經網絡各層之間的權重、偏置可以通過最小化均方誤差損失來學習，定義損失函數為

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

1.3 硬邊界PINN

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

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

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

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

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

式中，{xB, yB}為用于訓練的Neumann 邊界點數據，C0為Neumann邊界處的梯度值。

2 案例一：具有內熱源的二維穩態導熱方程

2.1 案例模型

本節首先以一個簡單的案例探究軟、硬兩種邊界條件設置方式下，PINN預測結果的差異。以具有內熱源的二維穩態導熱方程為例，其表示式如下

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

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

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

2.2 溫度場預測結果

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

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

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

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

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

3 案例二：平板間的二維穩態對流傳熱方程

3.1 案例模型

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

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

對于有內熱源的穩態層流傳熱，忽略黏性耗散后，其能量守恒式可簡化為

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

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

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

令空氣為板間流體，物性參數為：ρ=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 網格的坐標。最大速度umax在區間{0.01～0.03 m/s}內參數化，取10 個樣本。設置PINN 隱藏層4 層，每層含有50 個神經元，使用sigmoid函數[式(4)]作為激活函數。

3.2 溫度場預測結果

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

圖5 平板的幾何結構示意圖Fig.5 Geometric diagram of the flat plate

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

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

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

4 結 論

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