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

* (中國科學(xué)院力學(xué)研究所流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室,北京 100190)

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

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

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

引言

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

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

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

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

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

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

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

其中,下標(biāo)L表示重流體的物理量,下標(biāo)G表示輕流體的物理量.

2 物理融合神經(jīng)網(wǎng)絡(luò)

物理融合神經(jīng)網(wǎng)絡(luò)主要用于構(gòu)建流體動力學(xué)模型,即建立時(shí)空坐標(biāo)與待求解物理量之間的映射關(guān)系.在原始PF-PINNs 中,該映射關(guān)系可表示為如下形式

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

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

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

為了降低計(jì)算圖的規(guī)模,可在盡可能保證計(jì)算精確的情況下,由軟約束條件限定化學(xué)能的值.改用軟約束條件要求神經(jīng)網(wǎng)絡(luò)將化學(xué)能作為輸出,故采用MIM 思想改進(jìn)后的神經(jīng)網(wǎng)絡(luò)映射關(guān)系如下

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

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

根據(jù)控制方程式(1)～式(4),可給出各方程殘差項(xiàng)的表達(dá)形式

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

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

離散形式的總誤差Ltotal確定以后,神經(jīng)網(wǎng)絡(luò)中權(quán)重和偏差的取值 Θ 可通過隨機(jī)梯度下降算法更新,該算法調(diào)節(jié) Θ 的取值使總誤差下降.本文中采用Adam 方法[28]作為隨機(jī)梯度下降的優(yōu)化算法,并在訓(xùn)練開始時(shí)使用Xavier 初始化方法[29]隨機(jī)初始化神經(jīng)網(wǎng)絡(luò)的權(quán)重和偏差.完成一輪優(yōu)化后,神經(jīng)網(wǎng)絡(luò)將重新在采樣空間中進(jìn)行采樣,并重復(fù)上述流程,直到神經(jīng)網(wǎng)絡(luò)達(dá)到最大的訓(xùn)練次數(shù)或是總誤差降低到容許值內(nèi).整個(gè)算法訓(xùn)練的流程示意圖如圖2 所示.

圖2 基于相場法的物理融合神經(jīng)網(wǎng)絡(luò)(PF-PINNs)訓(xùn)練架構(gòu)圖Fig.2 Framework of physics-informed neural networks for phase-field(PF-PINNs) method

3 針對RT 不穩(wěn)定性問題的神經(jīng)網(wǎng)絡(luò)訓(xùn)練結(jié)果

3.1 算例設(shè)置

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

圖3 RT 不穩(wěn)定性問題的計(jì)算區(qū)域和初始條件Fig.3 Computational domain and initial condition of RT instability

假設(shè)初始時(shí)刻全場的速度為 0 m/s,則所有的初始條件已完備,擾動將在重力的影響下隨時(shí)間發(fā)展.計(jì)算區(qū)域的左側(cè)和右側(cè)設(shè)置為滑移邊界條件,上端和下端則設(shè)置為無滑移邊界條件.RT 不穩(wěn)定性問題的物理參數(shù)以列表的形式展示于表2 中,本問題的關(guān)鍵無量綱數(shù)為雷諾數(shù)Re=ρLDU/μL和阿特伍德數(shù)其中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,本問題中不考慮表面張力的作用,可認(rèn)為 σ0=0.界面相關(guān)的參數(shù)包括界面厚度和界面滲透率,參考之前的文獻(xiàn)研究,界面厚度可設(shè)為 ε=0.01 m,界面滲透率則設(shè)為M0=10 μm·(N·s)-1.

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

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

圖4 t=3 s 時(shí),不同網(wǎng)格密度下PF-PINNs 得到的界面外形.其中紅色實(shí)細(xì)線為粗網(wǎng)格的結(jié)果,藍(lán)色虛細(xì)線為中等網(wǎng)格的結(jié)果,黑色散點(diǎn)為細(xì)網(wǎng)格的結(jié)果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 均勻網(wǎng)格間隔對結(jié)果的影響Table 3 Influence of the interval in x direction and y direction on a uniform grid

本問題的離散采樣點(diǎn)則從均勻網(wǎng)格中隨機(jī)抽取,在每次循環(huán)開始時(shí),PF-PINNs 從采樣區(qū)域內(nèi)部取60000 個(gè)點(diǎn)作為內(nèi)部點(diǎn),從t=0 時(shí)刻抽取30000個(gè)點(diǎn)作為初值條件采樣點(diǎn),從每個(gè)邊界各取2000 個(gè)點(diǎn)作為邊界條件采樣點(diǎn),共計(jì)8000 個(gè)邊界采樣點(diǎn),即=(60000,30000,2000×4).

在神經(jīng)網(wǎng)絡(luò)訓(xùn)練過程中,本文使用了時(shí)間分割策略.使用多個(gè)神經(jīng)網(wǎng)絡(luò)依次訓(xùn)練的方法借鑒了Wight 等[30]提出的“time marching strategy”的思想,第一個(gè)網(wǎng)絡(luò)訓(xùn)練時(shí)間區(qū)間t∈[0,0.5] s 的結(jié)果,待第一個(gè)網(wǎng)絡(luò)訓(xùn)練收斂后,將第一個(gè)網(wǎng)絡(luò)在t=0.5 s 時(shí)的預(yù)測值傳出至第二個(gè)網(wǎng)絡(luò),作為第二個(gè)網(wǎng)絡(luò)的初值條件.第二個(gè)網(wǎng)絡(luò)負(fù)責(zé)訓(xùn)練時(shí)間區(qū)間t∈[0.5,1] s 的結(jié)果,訓(xùn)練完成后將結(jié)果傳給第三個(gè)網(wǎng)絡(luò),依此類推,直到所有時(shí)間區(qū)間上的網(wǎng)絡(luò)均被訓(xùn)練.根據(jù)RT 不穩(wěn)定性問題的特性,本文采用了如下的分割策略:當(dāng)t∈[0,1.5] s 時(shí),以0.5 s 作為時(shí)間間隔分割計(jì)算區(qū)域;當(dāng)t∈[1.5,3] s 時(shí),以0.25 s 作為時(shí)間間隔分割計(jì)算區(qū)域.綜上,本問題需要訓(xùn)練9 個(gè)神經(jīng)網(wǎng)絡(luò)以順利求解,其占據(jù)的時(shí)間區(qū)間分別為 [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.這樣的設(shè)置是由于RT 不穩(wěn)定性問題的物理特性所導(dǎo)致的,當(dāng)t∈[0,1.5] s 時(shí)界面形態(tài)較為穩(wěn)定,流場特征并不豐富,采用較大的時(shí)間間隔可更高效地進(jìn)行計(jì)算.但當(dāng)t∈[1.5,3] s 時(shí)界面迅速變形,速度場變化較大,為保證神經(jīng)網(wǎng)絡(luò)結(jié)果收斂,需要縮小每個(gè)網(wǎng)絡(luò)的時(shí)間區(qū)間范圍.時(shí)間區(qū)間的長度需要根據(jù)物理問題的特征合理控制,以平衡神經(jīng)網(wǎng)絡(luò)訓(xùn)練的效率與精度.

3.2 訓(xùn)練結(jié)果

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

圖5 RT 不穩(wěn)定性問題的相分?jǐn)?shù)演化和文獻(xiàn)[33]結(jié)果的對比Fig.5 Evolution of volume fraction of RT instability and comparison with Ref.[33]

圖6 在 t=3.0 s 時(shí)RT 不穩(wěn)定性問題中:(a) 改進(jìn)的PF-PINNs,(b) 高精度譜元法的速度矢量對比和(c) 改進(jìn)的PF-PINNs 與高精度譜元法的界面形態(tài)對比.輕重流體交界面由 C=-0.8 時(shí)的等值線定義,速度矢量的大小在圖中由顏色條表示.圖6(c)中紅線為改進(jìn)的PF-PINNs 的界面,藍(lán)線為高精度譜元法的界面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)

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

表4 各個(gè)典型時(shí)刻神經(jīng)網(wǎng)絡(luò)與高精度譜元法的相對誤差和訓(xùn)練結(jié)束時(shí)的總誤差Table 4 Relative errors in various typical times between neural network and accurate spectral element method and the total loss when training ends

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

3.3 改進(jìn)的PF-PINNs 加速效果分析

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

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

圖8 t ∈[0,0.5] s 時(shí)采用了混合殘差方法(紅線)與未采用混合殘差方法(藍(lán)線)的PF-PINNs:(a) 總殘差隨迭代次數(shù)的變化,(b) 總殘差隨時(shí)間的變化.其中圖8 (b)中的時(shí)間指的是GPU 計(jì)算用時(shí)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 時(shí)采用了混合殘差方法(紅線)與未采用混合殘差方法(藍(lán)線)的PF-PINNs:(a) 總殘差隨迭代次數(shù)的變化,(b) 總殘差隨時(shí)間的變化.其中(b)中的時(shí)間指的是GPU 計(jì)算用時(shí)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 結(jié)論

為了實(shí)現(xiàn)兩相流場界面的高精度和高效捕捉,本文進(jìn)一步改進(jìn)了基于相場法的物理融合神經(jīng)網(wǎng)絡(luò)(PF-PINNs),參考深度混合殘差方法,將輔助變量化學(xué)能作為神經(jīng)網(wǎng)絡(luò)的獨(dú)立輸出,從而在相分?jǐn)?shù)輸運(yùn)方程中拆分出關(guān)于化學(xué)能的額外控制方程,同時(shí)考慮兩個(gè)方程的損失函數(shù)..

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

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