改進的物理融合神經網絡在瑞利-泰勒不穩定性問題中的應用1)

丘潤荻 王靜竹 黃仁芳 杜特專 王一偉, 黃晨光

* (中國科學院力學研究所流固耦合系統力學重點實驗室,北京 100190)

? (中國科學院大學未來技術學院,北京 100049)

** (中國科學院大學工程科學學院,北京 100049)

?? (中國科學院合肥物質科學研究院,合肥 230031)

引言

兩相流動演化問題普遍存在于海洋工程、生物醫學、工業制造等領域[1-5].數值模擬方法是工程中分析兩相流動問題的最重要手段,兩相界面的捕捉精度直接決定計算與分析的準確程度.因此,提高相界面捕捉精度和降低計算誤差一直是數值模擬兩相流動演化過程的熱點問題[6-9].相場法(phase-field method)由于具有物理含義明確、算法質量守恒和界面捕捉精度高等優勢,被廣泛用于求解黏性與表面張力主導的界面演化、大密度比水汽界面演化、界面強非線性演化、熱毛細對流等問題[10-13].其中,瑞利-泰勒(Rayleigh-Taylor,RT)不穩定性包含非定常、強非線性等因素,對數值方法的計算規模和計算精度要求很高[14-17].經典的RT 不穩定性發生在密度分層流體中,當重密度流體處于輕密度流體之上時,在重力作用下,相界面上的小擾動會逐漸演化發展,重密度流體向下延伸形成一個尖峰形狀的頭部,下方的輕密度流體向上移動形成氣泡的形狀,最終出現不斷擴展的湍流混合層.

針對基于神經網絡的兩相流場求解問題,最近引人關注的是通過物理融合神經網絡(physicsinformed neural networks,PINNs) 求解偏微分方程[18-22],通過添加控制方程和邊界條件作為神經網絡的約束,大大提高流場建模和反演能力.Raissi等[23]基于PINNs 框架,使用已知的濃度分布信息推斷未知的速度場,并將該方法成功應用于圓柱繞流和三維血管造影中.Cai 等[24]在PINNs 框架下,使用實驗獲得的溫度場反推出了低雷諾數范圍下流場的流動情況,并將該方法推廣到傳熱問題的流動分析中.Buhendwad 等[25]探索了PINNs 在正問題和反問題中的應用,計算了包括Couette 流動、Poiseuille 流動、氣泡變形、氣泡震蕩和氣泡上升等兩相流算例,驗證了PINNs 求解兩相流問題的可行性.

Qiu 等[26]使用相場法作為界面捕捉方法,結合納維-斯托克斯(Navier-Stokes,N-S)方程作為控制方程,建立了基于相場法的物理融合神經網絡(PINNs for phase-field method,PF-PINNs),并通過單渦剪切流和氣泡上升算例,均獲得了滿意的效果,證實了PF-PINNs 求解大密度比兩相流的可靠性.然而,在該架構中,PF-PINNs 直接求解相場法方程中最高階導數項,速度仍有待提升.針對PINNs 的效率問題,Lü等[27]提出了深度混合殘差方法(deep mixed residual method,MIM),可以加速神經網絡的訓練過程,該方法通過將一個高階導數降為多個低階導數來提升求解速度,該架構成功用于Poisson 方程、Monge-Ampere方程、Biharmonic 方程和Kortewegde Vries 方程的求解中.

為了提升兩相流PINNs 的計算效率,本文在PF-PINNs 框架基礎上,利用深度混合殘差方法,加入化學能作為神經網絡中的輔助變量,用于求解RT 不穩定性問題,來分析改進的PF-PINNs 的計算精度和計算效率.本文的具體內容包括:第一節中介紹了兩相流問題的基本控制方程;第二節中介紹物理融合神經網絡方法以及計算過程中采用的加速策略;第三節中將改進的PF-PINNs 用于求解RT 不穩定性算例,并分析了改進算法的加速效果;最后一節中對本文進行總結和展望.

1 基于相場模型的二維N-S 方程

兩相流動中的控制方程包括連續性方程、動量方程與界面演化方程.本文的界面演化方程由基于Cahn-Hillard 模型的相場法描述,其考慮不可壓縮不相溶兩相流體組成的混合流體,并假設界面擴散速度和化學能的梯度成正比來精確描述界面附近的擴散效應,確保界面演化方程在滿足質量守恒的同時能準確捕捉界面.關于相場法的推導和本文使用的相分數定義、范圍等可參考Qiu 等[26]的工作.將相場法的界面演化方程與原始的N-S 方程結合,可導出不可壓、不相溶兩相流動的控制方程與表面張力的形式

其中,C為混合物的相分數,本文中相分數的范圍定義在-1 到1 之間,即在重流體中C=1,在輕流體中C=-1,界面附近的區域內-1 ≤C≤1.u=(u,v) 為流場的速度矢量,p表示壓強,為重力加速度矢量;M0為滲透率,φ為化學能,ε為界面寬度;為表面張力,σ0為表面張力系數.ρM和 μM為混合物的密度和黏度,其可分別表示為

其中,下標L表示重流體的物理量,下標G表示輕流體的物理量.

2 物理融合神經網絡

物理融合神經網絡主要用于構建流體動力學模型,即建立時空坐標與待求解物理量之間的映射關系.在原始PF-PINNs 中,該映射關系可表示為如下形式

其中,(x,y,t) 為時空坐標,(u,v,C,p) 為待求物理量,fNN為三個輸入、四個輸出的全連接神經網絡,Θ 為全連接神經網絡的權重和偏差.將式(1)～式(4)通過殘差形式寫入神經網絡,可得到PF-PINNs 的物理約束部分.為了提升PF-PINNs 的計算效率,本文參考MIM[27]的思想,將化學能視為一個輔助變量,也加入到神經網絡的輸出當中,并修改了損失函數中與相分數方程有關的約束形式.本文所述普通形式神經網絡與深度混合殘差形式神經網絡在界面演化方程和相應損失函數的區別如表1 所示.

表1 普通形式神經網絡與深度混合殘差形式神經網絡在界面演化方程中的區別Table 1 Differences between a normal neural network and a deep mixed residual method network in interface evolution equation

由化學能的公式可知,化學能與物理問題之間并不存在直接關聯;它由相場法的定義給出,且在兩相流問題中僅會影響界面附近的區域.但化學能的表達式卻引入了相分數的二階導數,高階導數的存在擴大了神經網絡計算圖的規模,使自動微分的計算成本上升.為了使讀者大致了解MIM 策略生效的原因,本文以自動微分計算二階導uxx為例進行說明.在大部分機器學習庫中,自動微分通過函數ux=tf.gradient(u,x)實現,通過u計算ux需要調用一次計算圖;而對uxx的計算則通過uxx=tf.gradient(ux,x)實現,由于調用ux的同時也需要調用ux的計算圖,故通過ux計算uxx需要調用兩次計算圖,那么通過u計算uxx需要調用三次計算圖.依此類推,通過u計算uxxx需要調用六次計算圖,通過u計算uxxxx需要調用十次計算圖,計算高階導數的開銷可通過降低導數階數得到,若增加一個網絡輸出q,并將uxxxx改寫為qxx,令q=uxx,則計算qxx共需要調用2 × (1+2)=6 次計算圖.雖然將硬約束改為軟約束會使結果損失一定的精度,但是計算所需的開銷大幅下降,這是MIM 策略在PF-PINNs 中生效的主要原因.

為了降低計算圖的規模,可在盡可能保證計算精確的情況下,由軟約束條件限定化學能的值.改用軟約束條件要求神經網絡將化學能作為輸出,故采用MIM 思想改進后的神經網絡映射關系如下

根據改進后的映射關系與物理約束,可建立如圖1 所示的PF-PINNs 框架,圖中的 σ 表示非線性激活函數,本文選用如下形式的Swish 函數:σ(Y)=在每一個迭代步中,程序首先通過神經網絡得到時空坐標到物理量的函數映射關系;隨后使用自動微分計算物理量的各階偏導數;再由物理方程和初邊界條件確定總殘差值,優化器根據總殘差值對神經網絡的參數進行優化;上述過程將被反復執行,直到每次隨機采樣的總殘差值不再下降或迭代步達到最大,此時可認為訓練完成.

圖1 改進的基于相場法的物理融合神經網絡(PF-PINNs)簡圖Fig.1 Sketch of modified physics-informed neural networks for phase-field (PF-PINNs) method

根據控制方程式(1)～式(4),可給出各方程殘差項的表達形式

其中,?a表示而 ?ab表示LM,LU,LV,LC,Lφ為各控制方程的殘差項,其分別對應式(1)～式(4)所對應的五個標量形式方程式,包括連續性方程、兩個方向的動量方程、界面演化方程與化學能方程.神經網絡的訓練過程主要是在離散采樣點上完成的,故式(10)～式(14)均需表示為離散形式的誤差,以離散采樣點取值表示的總誤差、方程誤差、初始條件誤差和邊界條件誤差可寫為

其中,Ltotal,LEqn,LICs,LBCs分別為總誤差、方程誤差、初始條件誤差和邊界條件誤差.下標“Eqn”、“ICs”和“BCs”分別代表物理量或坐標在內部區域、初始條件和邊界條件上的取值,上標“i”代表的是第i個采樣點上的信息.NEqn,NICs,NBCs分別為內場區域、初值條件和邊值條件上采樣點的個數,物理量在第i個初邊界采樣點上的取值則根據物理問題的初邊值條件確定.

離散形式的總誤差Ltotal確定以后,神經網絡中權重和偏差的取值 Θ 可通過隨機梯度下降算法更新,該算法調節 Θ 的取值使總誤差下降.本文中采用Adam 方法[28]作為隨機梯度下降的優化算法,并在訓練開始時使用Xavier 初始化方法[29]隨機初始化神經網絡的權重和偏差.完成一輪優化后,神經網絡將重新在采樣空間中進行采樣,并重復上述流程,直到神經網絡達到最大的訓練次數或是總誤差降低到容許值內.整個算法訓練的流程示意圖如圖2 所示.

圖2 基于相場法的物理融合神經網絡(PF-PINNs)訓練架構圖Fig.2 Framework of physics-informed neural networks for phase-field(PF-PINNs) method

3 針對RT 不穩定性問題的神經網絡訓練結果

3.1 算例設置

RT 不穩定性問題的計算域和邊界條件如圖3 所示,其中整個區域的空間大小為Ω ∈[0,0.5] m×[0,4] m,計算總時長為t∈[0,3] s.計算區域上方為重流體L,下方為輕流體G,兩流體間的初始分界線可定義為一正弦函數形式的擾動,其擾動函數為h=-0.1cos(2πx).根據擾動函數的物理含義,可定義初始時刻相分數場滿足的關系式為

圖3 RT 不穩定性問題的計算區域和初始條件Fig.3 Computational domain and initial condition of RT instability

假設初始時刻全場的速度為 0 m/s,則所有的初始條件已完備,擾動將在重力的影響下隨時間發展.計算區域的左側和右側設置為滑移邊界條件,上端和下端則設置為無滑移邊界條件.RT 不穩定性問題的物理參數以列表的形式展示于表2 中,本問題的關鍵無量綱數為雷諾數Re=ρLDU/μL和阿特伍德數其中D為特征長度,為本問題的特征速度;本問題中特征長度為1 m,特征速度為1 m/s.重流體密度為 ρL=3 kg/m3,輕流體密 度則 為 ρG=1 kg/m3,兩種流體的黏度均為μL=μG=0.001 N·s/m2,重力加速度為g=(0,1) m/s2,本問題中不考慮表面張力的作用,可認為 σ0=0.界面相關的參數包括界面厚度和界面滲透率,參考之前的文獻研究,界面厚度可設為 ε=0.01 m,界面滲透率則設為M0=10 μm·(N·s)-1.

表2 RT 不穩定性算例中的物理參數和無量綱值Table 2 Physical parameters and non-dimensional value in RT instability

本文用于計算RT 不穩定性的神經網絡包含15 個隱藏層,每層由100 個神經元組成.該神經網絡的初始學習率設置為10-3,并隨著循環次數的增加而下降至10-5,總循環次數為20 萬次.本問題的均勻網格在時空采樣區域 (x,y,t) 中等間距生成,為了保證該問題中采樣尺寸足夠合適,在正式開始計算之前本文對不同尺度的均勻網格工況進行了驗證.本文在驗證時分別選取了細、中、粗三套均勻網格,每套網格在x和y方向的間隔如表3 所示.本文以高精度譜方法在t=2.75 s 的速度場和相場作為初始值,隨后用神經網絡計算t=2.75～ 3 s 的流場結果,并于圖4 中給出了三種方法在t=3 s 時的界面位置,以證明三種方法在復雜流場下最終計算結果基本一致.結果表明,三個網格數量下,神經網絡對相分數場的計算結果相同.由于三個網格數量下,計算結果沒有顯著的差別,為了提升收斂速度,本文采用細網格作為所選網格,選擇的均勻網格在空間兩個方向上的數量分別為nx,ny=(65,512),即網格尺寸分別為 Δx=7.7 mm,Δy=7.7 mm;均勻網格在時間上的時間步長為Δt=2 ms.

圖4 t=3 s 時,不同網格密度下PF-PINNs 得到的界面外形.其中紅色實細線為粗網格的結果,藍色虛細線為中等網格的結果,黑色散點為細網格的結果Fig.4 The shape of interface in various intervals when t=3 s.The red solid line denotes the result from coarse mesh,the blue dashed line denotes the result from medium mesh,and the black scattered line denotes the result from fine mesh

表3 均勻網格間隔對結果的影響Table 3 Influence of the interval in x direction and y direction on a uniform grid

本問題的離散采樣點則從均勻網格中隨機抽取,在每次循環開始時,PF-PINNs 從采樣區域內部取60000 個點作為內部點,從t=0 時刻抽取30000個點作為初值條件采樣點,從每個邊界各取2000 個點作為邊界條件采樣點,共計8000 個邊界采樣點,即=(60000,30000,2000×4).

在神經網絡訓練過程中,本文使用了時間分割策略.使用多個神經網絡依次訓練的方法借鑒了Wight 等[30]提出的“time marching strategy”的思想,第一個網絡訓練時間區間t∈[0,0.5] s 的結果,待第一個網絡訓練收斂后,將第一個網絡在t=0.5 s 時的預測值傳出至第二個網絡,作為第二個網絡的初值條件.第二個網絡負責訓練時間區間t∈[0.5,1] s 的結果,訓練完成后將結果傳給第三個網絡,依此類推,直到所有時間區間上的網絡均被訓練.根據RT 不穩定性問題的特性,本文采用了如下的分割策略:當t∈[0,1.5] s 時,以0.5 s 作為時間間隔分割計算區域;當t∈[1.5,3] s 時,以0.25 s 作為時間間隔分割計算區域.綜上,本問題需要訓練9 個神經網絡以順利求解,其占據的時間區間分別為 [0,0.5] s,[0.5,1] s,[1,1.5]s 以及[1.5,1.75]s,[1.75,2]s,[2,2.25]s,[2.25,2.5]s,[2.5,2.75]s和 [2.75,3] s.這樣的設置是由于RT 不穩定性問題的物理特性所導致的,當t∈[0,1.5] s 時界面形態較為穩定,流場特征并不豐富,采用較大的時間間隔可更高效地進行計算.但當t∈[1.5,3] s 時界面迅速變形,速度場變化較大,為保證神經網絡結果收斂,需要縮小每個網絡的時間區間范圍.時間區間的長度需要根據物理問題的特征合理控制,以平衡神經網絡訓練的效率與精度.

3.2 訓練結果

RT 不穩定性問題的界面變化情況如圖5 所示.其中圖5(a)為改進的PF-PINNs 計算結果界面變化,界面處的初始擾動隨著時間的增加而擴大,重流體從左側邊界的凹陷處下沉,輕流體則從右側的凸起處上浮,分別形成了尖釘和氣泡.在t∈[1.5,2] s 時,界面處的流體發生劇烈卷曲,在尖釘附近形成了復雜的渦結構,隨輕重流體的反向運動被拉伸并向外擴張,重流體旋流結構非常復雜.圖5(b)給出了高精度譜元法[32]的計算結果用作對比,其與圖5(a)形態非常接近,表明改進的PF-PINNs 方法可接近傳統算法的計算精度.圖6 還給出了神經網絡結果與高精度譜元法的計算結果在t=3.0 s 時 界面附近((x,y)∈[0,0.5]×[-1,0])速度場矢量的對比以及t=3.0 s 時兩種算法的界面位置對比,結果表明,改進的PF-PINNs 計算的渦結構與傳統算法得到的渦結構形狀吻合,且流場中速度矢量方向與大小基本一致,矢量場各處疏密程度無明顯差別.改進的PF-PINNs 的計算結果與高精度譜元法相比,在卷曲處的界面演化稍有差別,但在流場中大部分區域內界面形狀保持一致,這說明了改進的PF-PINNs 的計算結果符合物理規律.

圖5 RT 不穩定性問題的相分數演化和文獻[33]結果的對比Fig.5 Evolution of volume fraction of RT instability and comparison with Ref.[33]

圖6 在 t=3.0 s 時RT 不穩定性問題中:(a) 改進的PF-PINNs,(b) 高精度譜元法的速度矢量對比和(c) 改進的PF-PINNs 與高精度譜元法的界面形態對比.輕重流體交界面由 C=-0.8 時的等值線定義,速度矢量的大小在圖中由顏色條表示.圖6(c)中紅線為改進的PF-PINNs 的界面,藍線為高精度譜元法的界面Fig.6 Comparison of the velocity vector when t=3.0 s between (a) modified PF-PINNs,(b) high-order spectral element method and (c) interface evolution in RT instability.The interface between heavier fluid and lighter fluid are defined using isoline C=-0.8,and the magnitude of velocity vector is shown with color bar.The red line is the interface from the modified PF-PINNs,while the blue line is the interface from the high-order spectral element method in Fig.6(c)

為了定量分析神經網絡的計算結果,本文以高精度譜元法的結果作為精確解,計算了神經網絡結果相較于高精度譜元法在各個典型時刻的相對誤差,以及各時刻計算結束時的總誤差Ltotal,并列于表4 中.本文的相對誤差定義為其中qnn為神經網絡中物理量q的結果,qaccurate為精確解中物理量q的結果.通過表4 的結果可以看出,神經網絡訓練結果在0～ 2.5 s 內均維持著較好的精度,在2.5～ 3 s 間由于速度場的非線性演化導致誤差快速上升.總體而言,隨著訓練時間的推移,神經網絡的總誤差值在不斷上升,這也導致了速度相對誤差的不斷上升,當誤差過大時計算結果將偏離真實解.本文中訓練次數、學習率等參數是一直不變的,當界面流動趨于復雜時,采用原有的超參數應對當前問題可能無法保證結果有較好的收斂性.適當調節每個子網絡在不同時刻的超參數可以緩解這種現象,提升神經網絡的精度.圖7 展示了尖釘和氣泡在y方向上的位置隨時間的變化曲線,并與文獻[31-33]給出的結果進行了對比.其中紅色線條為本文神經網絡計算所得結果,其他形狀的散點為參考文獻給出的結果.對比結果表明,改進的PF-PINNs的模擬結果與參考文獻相吻合,體現了改進的PFPINNs 在宏觀結果上的計算精度.

表4 各個典型時刻神經網絡與高精度譜元法的相對誤差和訓練結束時的總誤差Table 4 Relative errors in various typical times between neural network and accurate spectral element method and the total loss when training ends

圖7 氣泡與尖釘位置隨時間演化過程與參考文獻結果對比Fig.7 Comparison of the position of bubble and spike between references and modified PF-PINNs

3.3 改進的PF-PINNs 加速效果分析

本文第2 節指出,深度混合殘差方法[27]能夠顯著提升高階微分方程的計算速度,同時保證神經網絡的精度.為了驗證改進方法在簡單流場與復雜流場中的效果,本文對比了原始PF-PINNs 與改進的PF-PINNs 在t∈[0,0.5] s 和t∈[2.75,3.0] s 時總誤差隨計算時間的變化曲線以及總誤差隨迭代次數的變化曲線,相關數據列于圖8 與圖9 中.用于對比的神經網絡大小為10 層100 節點,學習率為10-3,迭代次數為80000 次,其余設置與3.1 節中給出的設置完全一致.

由圖8(a)與圖9(a)的結果可以看出,在其他所有條件相同的情況下,使用MIM 與不使用MIM 得到的總誤差沒有太大的區別,兩者之間的總誤差均在同一量級且差別不大.但在圖8(b)的結果中,在其他所有條件相同的情況下,采用MIM 策略訓練用時為40720 s,平均單步用時為0.509 s;而未采用MIM 策略訓練用時為115652 s,平均單步用時為1.446 s;采用MIM 策略使平均單步用時減少了64.80%.在圖9(b)的結果中,在其他所有條件相同的情況下,采用MIM 策略訓練用時為34281 s,平均單步用時為0.429 s;而未采用MIM 策略訓練用時為113640 s,平均單步用時為1.420 s;采用MIM 策略使平均單步用時減少了69.79%.以上結果證實了深度混合殘差方法在計算初始階段和流場復雜演化階段的有效性與可靠性.

圖8 t ∈[0,0.5] s 時采用了混合殘差方法(紅線)與未采用混合殘差方法(藍線)的PF-PINNs:(a) 總殘差隨迭代次數的變化,(b) 總殘差隨時間的變化.其中圖8 (b)中的時間指的是GPU 計算用時Fig.8 Variations of total loss with (a) iterations and (b) times using MIM (red line) and no MIM (blue line) in t ∈[0,0.5] s.The label“time”in Fig.8(b) means the GPU computational time of PF-PINNs

圖9 t ∈[2.75,3.0] s 時采用了混合殘差方法(紅線)與未采用混合殘差方法(藍線)的PF-PINNs:(a) 總殘差隨迭代次數的變化,(b) 總殘差隨時間的變化.其中(b)中的時間指的是GPU 計算用時Fig.9 Variations of total loss with (a) iterations and (b) times using MIM (red line) and no MIM (blue line) in t ∈[2.75,3.0] s.The label“time”in (b) means the GPU computational time of PF-PINNs

4 結論

為了實現兩相流場界面的高精度和高效捕捉,本文進一步改進了基于相場法的物理融合神經網絡(PF-PINNs),參考深度混合殘差方法,將輔助變量化學能作為神經網絡的獨立輸出,從而在相分數輸運方程中拆分出關于化學能的額外控制方程,同時考慮兩個方程的損失函數..

本文基于改進的PF-PINNs,求解了經典的Rayleigh-Taylor 不穩定問題,從而考核當存在強非線性因素條件下這一方法的求解能力、精度和效率.與高精度譜元法結果對比表明,本文的改進的PF-PINNs 有能力捕捉到以RT 不穩定為代表的兩相界面的強非線性演化過程,計算結果符合物理規律.尖釘和氣泡位移結果與參考文獻結果定量吻合良好.與此同時,所采用的深度混合殘差方法能夠有效地提升計算速度,RT 不穩定算例中,在其他所有條件相同的情況下,采用MIM 策略可使訓練用時減少60%以上,在保證精度的前提下帶來了較明顯的效率提升.

上述結果表明,PF-PINNs 方法以及深度混合殘差改進方法可用于精細求解兩相流界面問題.考慮到相比CFD 方法,PINNs 方法以及神經網絡方法在并行效率、軟硬件結合等方面具有較強的優勢,在反問題求解和數據融合等方面均具有較好的應用靈活性,因此PINNs 相關方法在未來應用的潛力很大.本文工作可為進一步提升神經網絡的訓練速度、探索高精度數據同化方法等提供參考.