基于稀疏數(shù)據(jù)的流場(chǎng)結(jié)構(gòu)重構(gòu)方法研究

陳曠奇，歐陽晗青，朱志成，郝佳，黃彪

（北京理工大學(xué) 機(jī)械與車輛學(xué)院，北京 100081）

高精度和高可信度的流場(chǎng)數(shù)據(jù)往往伴隨著高昂的時(shí)間成本與經(jīng)濟(jì)成本，以直接數(shù)值模擬(DNS)為例，雖然它對(duì)湍流流場(chǎng)的求解精度得到了學(xué)術(shù)界的廣泛認(rèn)可，并且已經(jīng)有了初步應(yīng)用，但網(wǎng)格數(shù)量至少需要能達(dá)到基本的計(jì)算要求，而急劇增加的網(wǎng)格數(shù)量嚴(yán)重限制了其應(yīng)用范圍.與此同時(shí)，流體實(shí)驗(yàn)所得的有限數(shù)據(jù)也無法提供更為精細(xì)的流場(chǎng)結(jié)構(gòu),在相關(guān)氣/水動(dòng)工程問題中，研究者往往需要根據(jù)有限實(shí)驗(yàn)數(shù)據(jù)反復(fù)迭代，不僅主觀因素強(qiáng)，而且成本極高，極大地制約了研究效率.與此同時(shí)，流動(dòng)條件愈加復(fù)雜，僅通過數(shù)值模擬或?qū)嶒?yàn)測(cè)量獲得高精度數(shù)據(jù)也更加困難，因此，亟需一種數(shù)據(jù)挖掘手段，可以從有限數(shù)據(jù)中獲得所需的研究信息.

機(jī)器學(xué)習(xí)技術(shù)作為數(shù)據(jù)挖掘的重要手段，近年展現(xiàn)出巨大的潛力，而深度學(xué)習(xí)便是其重要分支，能夠在海量數(shù)據(jù)中自動(dòng)提取并學(xué)習(xí)最優(yōu)特征，在圖像識(shí)別[1]、自然語言理解[2]、認(rèn)知科學(xué)[3]，網(wǎng)絡(luò)漏洞分析[4-5]和基因組學(xué)[6]等多領(lǐng)域中已取得突破性進(jìn)展，逐漸在各類基礎(chǔ)研究和工程技術(shù)中發(fā)揮不可替代的作用.目前已有學(xué)者嘗試將機(jī)器學(xué)習(xí)技術(shù)應(yīng)用于流體力學(xué)中，在湍流建模[7]、主動(dòng)流體控制[8]、流場(chǎng)特征識(shí)別[9]以及多相流流態(tài)識(shí)別[10]上都有著突破性成果.傳統(tǒng)機(jī)器學(xué)習(xí)方法主要以數(shù)據(jù)驅(qū)動(dòng)為主，通過神經(jīng)網(wǎng)絡(luò)擬合輸入和輸出間的映射關(guān)系.這種方法實(shí)際上是從已知的信息源中尋找信息，因此對(duì)數(shù)據(jù)的數(shù)量和質(zhì)量有著較高的要求.與此同時(shí)，數(shù)據(jù)驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)將數(shù)據(jù)樣本作為單一信息源，必然導(dǎo)致其對(duì)高質(zhì)量數(shù)據(jù)集的過度依賴，這也是阻礙機(jī)器學(xué)習(xí)應(yīng)用于流體力學(xué)的重要原因.

流體力學(xué)研究以數(shù)理模型為主，積累了豐富的理論基礎(chǔ)，但這些信息尚未被數(shù)據(jù)驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)利用.將已有的物理模型作為高精度數(shù)據(jù)的補(bǔ)充融入神經(jīng)網(wǎng)絡(luò)，是深度學(xué)習(xí)技術(shù)與流體力學(xué)領(lǐng)域交叉融合的關(guān)鍵基礎(chǔ).隨著物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)框架(physics-informed neural networks，PINN)[11]被提出，數(shù)據(jù)驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)缺乏可解釋性且極度依賴高質(zhì)量數(shù)據(jù)集的現(xiàn)狀得到改善，立刻受到學(xué)術(shù)界的極大重視，在短時(shí)間內(nèi)出現(xiàn)了大量以PINN 框架為基礎(chǔ)的理論研究，例如非局部型PINN：nPINNs[12]、分?jǐn)?shù)階PINN：fPINNs[13]等.在流體力學(xué)領(lǐng)域，已有不少學(xué)者開始將PINN 框架應(yīng)用到具體流動(dòng)問題中.RAISSI等[14]基于PINN 框架提出了隱式流體力學(xué)（hidden fluid mechanics，HFM），僅利用少于1%的煙霧濃度測(cè)量數(shù)據(jù)，便重構(gòu)出全部流場(chǎng)，并通過與數(shù)值計(jì)算結(jié)果對(duì)比證實(shí)了該方法不受邊界形狀的限制，在保證精度的同時(shí)，有著遠(yuǎn)快于數(shù)值模擬的求解速度.EIVAZI等[15]將PINN 框架應(yīng)用于湍流求解，提出利用數(shù)據(jù)封閉雷諾平均方程，不需要任何湍流模型或假設(shè)即可求解雷諾平均方程.MAO 等[16]證實(shí)了利用PINN 框架近似模擬高速氣動(dòng)力學(xué)中Euler 方程的可行性.同時(shí)建立了PINNs 方法，當(dāng)求解域中存在高梯度區(qū)域時(shí)，將其劃分為若干子域，并對(duì)每一個(gè)子域獨(dú)立使用神經(jīng)網(wǎng)絡(luò)，可以提高求解精度.對(duì)于流體中的正問題，PINNs 的計(jì)算精度不如傳統(tǒng)的數(shù)值方法，但對(duì)于傳統(tǒng)數(shù)值計(jì)算方法無法求解的反問題，PINNs 更具有優(yōu)越性.GAO 等[17]將PINN 的思想引入到卷積神經(jīng)網(wǎng)絡(luò)，提出基于物理信息的DL-SR 方法，利用流體流動(dòng)守恒定律和邊界條件從高維參數(shù)空間的低分辨率數(shù)據(jù)中產(chǎn)生成超分辨流場(chǎng).CHENG 等[18]在PINN 框架的基礎(chǔ)上提出了Res-PINN 模型，通過將Resent 模塊與神經(jīng)網(wǎng)絡(luò)耦合進(jìn)一步增強(qiáng)了深度學(xué)習(xí)的預(yù)測(cè)能力，可以精準(zhǔn)預(yù)測(cè)流場(chǎng)中的速度和壓力信息.同時(shí)，也可以應(yīng)用在流體的反問題中，反演參數(shù)的誤差為0.98%和3.1%；在噪聲數(shù)據(jù)中，反演參數(shù)的誤差僅有0.99%和3.1%.ALMAJID 等[19]利用PINN 框架預(yù)測(cè)多孔介質(zhì)流體流動(dòng).通過融合物理信息和觀測(cè)數(shù)據(jù)來模擬Buckley-Leverett 問題，并通過神經(jīng)網(wǎng)絡(luò)反向求解，推導(dǎo)出最符合實(shí)驗(yàn)結(jié)果的多相流參數(shù).結(jié)果表明，即使沒有觀測(cè)數(shù)據(jù)，PINNs 也能捕捉到解的整體趨勢(shì)，但隨著觀測(cè)數(shù)據(jù)的增加，解的分辨率和精度都有很大提高.

文中首先建立了一種基于邊界數(shù)據(jù)浸入法的數(shù)值模擬方法，并把求解結(jié)果作為神經(jīng)網(wǎng)絡(luò)的高精度數(shù)據(jù)集, 然后基于PINN 框架構(gòu)建了耦合流場(chǎng)物理信息的神經(jīng)網(wǎng)絡(luò)模型，發(fā)展了圓柱繞流流場(chǎng)的數(shù)據(jù)采樣方法和神經(jīng)網(wǎng)絡(luò)模型優(yōu)化策略，基于稀疏流場(chǎng)數(shù)據(jù)建立了物理信息驅(qū)動(dòng)的流場(chǎng)結(jié)構(gòu)重構(gòu)方法.通過分析重構(gòu)流場(chǎng)的水動(dòng)力特性與渦脫落特性揭示了物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)誤差機(jī)理.對(duì)比不同流態(tài)的繞流流場(chǎng)，討論了該神經(jīng)網(wǎng)絡(luò)模型對(duì)不同流場(chǎng)結(jié)構(gòu)的預(yù)測(cè)能力.

1 數(shù)據(jù)集的準(zhǔn)備

高質(zhì)量數(shù)據(jù)集是訓(xùn)練神經(jīng)網(wǎng)絡(luò)模型的必要保證，對(duì)最終預(yù)測(cè)效果的影響甚至遠(yuǎn)遠(yuǎn)超過算法本身.文中基于Fortran 自編程求解了圓柱繞流流場(chǎng)結(jié)構(gòu)，通過兩步投影法結(jié)合泊松方程解耦N-S 方程，同時(shí)通過邊界數(shù)據(jù)浸入法（boundary data immersion method，BDIM）處理固/液邊界，為后續(xù)神經(jīng)網(wǎng)絡(luò)提供了高置信度的數(shù)據(jù)集.

1.1 基本控制方程

本節(jié)的研究對(duì)象剛性圓柱繞流流場(chǎng)，同時(shí)忽視了自由表面對(duì)模型附近流動(dòng)影響的圓柱繞流，計(jì)算中涉及到的介質(zhì)均為不可壓縮流體.因此流動(dòng)控制方程為

①連續(xù)性方程：

②動(dòng)量方程：

式 中： →u和p分 別 為 流 場(chǎng) 的 速 度 和 壓 力； ρ 和 μ分 別 為流 體 的 密 度 和 動(dòng) 力 黏 性 系 數(shù)； →g為 重力加速度； ?為梯度運(yùn)算符.

1.2 求解方法

文中通過Chorin 投影法求解N-S 方程組，同時(shí)采用邊界數(shù)據(jù)浸入法（BDIM）處理固/液邊界，關(guān)于Chorin 投影法和邊界數(shù)據(jù)浸入法的詳細(xì)推導(dǎo)過程，可以見文獻(xiàn)[20]和[21].圖1 展示了利用二步投影法和邊界數(shù)據(jù)浸入法求解不可壓N-S 方程的具體過程.

圖1 不可壓縮N-S 方程求解流程圖Fig.1 Flow chart of solving incompressible N-S equation

1.3 計(jì)算域及邊界條件設(shè)置

本算例是基于Fortran 自編程求解，流場(chǎng)計(jì)算域及邊界條件設(shè)置如圖2 所示，計(jì)算域建立在笛卡爾坐標(biāo)系中，其中x代表順流方向.以圓柱體的直徑D=0.064 m 為特征長(zhǎng)度，計(jì)算域選取在長(zhǎng)24D、寬5D、高5D的矩形流域內(nèi).計(jì)算采用入口速度條件，出口壓力條件設(shè)置，圓柱體的中軸線距離入口邊界和出口邊界分別為2.5D和12D，與上下邊界的距離為2.5D，入口邊界的距離為2.5D.

圖2 計(jì)算域和邊界條件設(shè)置Fig.2 Calculation domain and boundary condition setting

在文中的數(shù)值計(jì)算中,液體的密度和動(dòng)力黏度分別取ρg=997 kg/m3，μg=10-3Pa·s.25 °C 下，飽和蒸汽壓取pv=3 169 Pa，出口處壓力設(shè)為p0=43 169 Pa.

在本算例中采用高分辨率的精細(xì)化均布網(wǎng)格，在節(jié)省計(jì)算資源的同時(shí)，盡可能多地捕捉流場(chǎng)流動(dòng)細(xì)節(jié)，提高計(jì)算精度.文中針對(duì)工況Re=100 000，參考其他學(xué)者的實(shí)驗(yàn)結(jié)果(Cd=1.25，St=0.19)[20-22]，采用表1所示的四組不同尺度網(wǎng)格對(duì)圓柱的平均阻力系數(shù)和斯特勞哈爾書進(jìn)行了對(duì)比計(jì)算.

表1 網(wǎng)格參數(shù)設(shè)置及獨(dú)立性檢驗(yàn)Tab.1 Mesh parameter setting and independence verification.

從結(jié)果可以看出，隨著總節(jié)點(diǎn)數(shù)的增加，平均阻力系數(shù)逐漸收斂于試驗(yàn)結(jié)果，并且網(wǎng)格3 和網(wǎng)格4的結(jié)果非常相似.綜合考慮計(jì)算精度和計(jì)算資源，網(wǎng)格3 較為合理，將在本研究中使用.

為了探討PINN 框架對(duì)不同流態(tài)流場(chǎng)的重構(gòu)能力，文中分別計(jì)算了雷諾數(shù)Re=100，3 900，20 000 和100 000 四種工況，同時(shí)與其他學(xué)者實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，結(jié)果如表2 和表3 所示.

表2 平均阻力系數(shù)對(duì)比Tab.2 Comparison of average resistance coefficient

表3 斯特勞哈爾數(shù)對(duì)比Tab.3 Comparison of Strouhal numbers.

文中對(duì)比了不同雷諾數(shù)下平均阻力系數(shù)對(duì)比與斯特勞哈爾數(shù)，并且對(duì)比了Re=3 900 工況下時(shí)均流速的實(shí)驗(yàn)數(shù)據(jù)[22],可以看出，文中的數(shù)值計(jì)算方法在中低雷諾數(shù)下具有較高的精度.

如圖3 所示，平均阻力系數(shù)，斯特勞哈爾數(shù)以及時(shí)均流速均與實(shí)驗(yàn)吻合較好，因此可以認(rèn)為文中數(shù)值計(jì)算方法提供的數(shù)據(jù)集真實(shí).

圖3 軸向中心線上沿流動(dòng)方向的時(shí)均流速分布對(duì)比（Re=3 900）Fig.3 Comparison of time-averaged velocity distribution along the flow direction on the axial centerline (Re=3 900)

2 神經(jīng)網(wǎng)絡(luò)模型搭建與訓(xùn)練

PINN 框架將物理信息與神經(jīng)網(wǎng)絡(luò)結(jié)合，通過補(bǔ)充物理機(jī)理便可使基于稀疏數(shù)據(jù)的神經(jīng)網(wǎng)絡(luò)收斂.該框架耦合物理信息的基本原理是：基于物理方程構(gòu)造方程項(xiàng)殘差并與數(shù)據(jù)項(xiàng)損失函數(shù)加權(quán)，最終統(tǒng)一為神經(jīng)網(wǎng)絡(luò)的損失函數(shù).隨著網(wǎng)絡(luò)參數(shù)的迭代優(yōu)化，損失函數(shù)值逐漸減小至預(yù)定精度，本質(zhì)上是對(duì)控制方程的求解.該技術(shù)與求解偏微分方程的變分方法相似，不同之處在于使用深度神經(jīng)網(wǎng)絡(luò)模型作為試函數(shù)，利用機(jī)器學(xué)習(xí)方法進(jìn)行權(quán)重的迭代更新.

2.1 稀疏數(shù)據(jù)采樣

在流場(chǎng)數(shù)值模擬中，每一個(gè)網(wǎng)格節(jié)點(diǎn)均代表一個(gè)數(shù)據(jù)點(diǎn)，考慮到圓柱繞流流場(chǎng)的周期性，僅在計(jì)算域中選取圖4 所示虛線區(qū)域作為重構(gòu)區(qū)域，三個(gè)方向均是等比例均勻采點(diǎn)，數(shù)據(jù)點(diǎn)個(gè)數(shù)為240×100×100=2 400 000.以數(shù)值模擬一個(gè)時(shí)間步為一個(gè)單位，間隔選取流動(dòng)穩(wěn)定后的50 個(gè)時(shí)間步用于訓(xùn)練.

圖4 神經(jīng)網(wǎng)絡(luò)重構(gòu)區(qū)域Fig.4 Areas reconstructed by neural network.

物理機(jī)理地引入使得神經(jīng)網(wǎng)絡(luò)僅利用稀疏數(shù)據(jù)便可實(shí)現(xiàn)網(wǎng)絡(luò)收斂.考慮到實(shí)驗(yàn)中數(shù)據(jù)獲取的合理性，文中將數(shù)值計(jì)算的部分速度數(shù)據(jù)和圓柱表面壓力數(shù)據(jù)作為神經(jīng)網(wǎng)絡(luò)的訓(xùn)練數(shù)據(jù).流場(chǎng)速度數(shù)據(jù)如圖5 所示.

圖5 速度數(shù)據(jù)采樣方式Fig.5 Sampling method of velocity data

圖中灰色區(qū)域代表采樣點(diǎn)，x和y方向均采集數(shù)據(jù)點(diǎn)的1/4，最終流場(chǎng)采樣數(shù)據(jù)點(diǎn)為60×25×25=37 500；壓力數(shù)據(jù)采樣也采用等間隔采樣方式，由于液體的黏性，固體表面的速度理論上為0，但由于邊界數(shù)據(jù)浸入法會(huì)使圓柱表面存在一個(gè)圖6 所示的過渡區(qū)域.為了保證預(yù)測(cè)的準(zhǔn)確性，文中只將圖6 所示的邊界壓力數(shù)據(jù)代入神經(jīng)網(wǎng)絡(luò).

圓柱表面節(jié)點(diǎn)的采樣方式如圖7 所示.在周向和軸向方向上均采用均勻間隔采樣，最終界面采點(diǎn)數(shù)量為30×15=450.

圖7 壓力數(shù)據(jù)采樣方式Fig.7 Sampling method of pressure data

綜上，在固定時(shí)間步下，文中從流場(chǎng)中采樣數(shù)據(jù)如表4 所示.

表4 固定時(shí)間步下采樣數(shù)據(jù)匯總Tab.4 Sampling data summary under fixed time step

在PINN 框架中各物理變量需要代入方程求解方程殘差，由于min-max 標(biāo)準(zhǔn)化或Z-score 標(biāo)準(zhǔn)化會(huì)破壞變量之間的量綱關(guān)系，因此文中僅對(duì)數(shù)據(jù)僅進(jìn)行無量綱化處理，處理式為

式中：u，p和x分別為有量綱的速度，壓力和長(zhǎng)度；u*，p*和x*分別為無量綱的速度，壓力和長(zhǎng)度（下標(biāo)i分別為x,y,z三個(gè)方向），分別被來流速度u0、圓柱直徑D和流體密度ρ無量綱化.

2.2 神經(jīng)網(wǎng)絡(luò)構(gòu)建與訓(xùn)練方法

基于PINN 框架，文中通過耦合流動(dòng)機(jī)理得到如圖8 的物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)模型.

圖8 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.8 Architectures of neural network

該模型將時(shí)空信息（即流場(chǎng)位置坐標(biāo)與時(shí)間）作為神經(jīng)網(wǎng)絡(luò)的輸入，通過神經(jīng)網(wǎng)絡(luò)計(jì)算流場(chǎng)各物理量，包括三個(gè)方向的速度和壓力.首先，將神經(jīng)網(wǎng)絡(luò)的輸出流場(chǎng)結(jié)構(gòu)信息{u,v,w,p}與其標(biāo)簽信息{x,y,z,t}共同帶入偏微分方程，通過自動(dòng)微分技術(shù)求得方程誤差.方程誤差包括質(zhì)量方程誤差L1和動(dòng)量方程誤差L2、L3、L4，分別為

式中：x,y,z,t為樣本的標(biāo)簽物理量，文中指樣本點(diǎn)的時(shí)空坐標(biāo)(x,y,z,t)；u,v,w,p為神經(jīng)網(wǎng)絡(luò)的輸入數(shù)據(jù)集，分別為流向速度u、法向速度v、垂向速度w、流場(chǎng)壓力p；u',v',w',p':神經(jīng)網(wǎng)絡(luò)的輸出數(shù)據(jù)集，分別代表預(yù)測(cè)流向速度u'、法向速度v'、垂向速度w'、流場(chǎng)壓力p';N: 神經(jīng)網(wǎng)絡(luò)的輸入樣本點(diǎn)總量；i為第i個(gè)樣本數(shù)據(jù)的對(duì)應(yīng)編號(hào).之后求解神經(jīng)網(wǎng)絡(luò)的輸出流場(chǎng)信息{u',v',w',p'}與該時(shí)空坐標(biāo)下已知信息{u,v,w,p}間的誤差，基于本算例的采樣情況，該誤差分為流域速度誤差L5、L6、L7與邊界壓力誤差L8，其計(jì)算式為

①流域速度誤差

②邊界壓力誤差

式中：Nl為輸入神經(jīng)網(wǎng)絡(luò)的流域樣本點(diǎn)總量；Nb為輸入神經(jīng)網(wǎng)絡(luò)的邊界樣本點(diǎn)總量；u,v,w,p為采樣點(diǎn)的流場(chǎng)信息.

將兩種誤差相加，并通過平衡權(quán)重α調(diào)整數(shù)據(jù)誤差在總誤差中所占的比例，即可得到神經(jīng)網(wǎng)絡(luò)的損失函數(shù)

從表2 中可以看出，速度v量級(jí)遠(yuǎn)小于速度u,w，這是由于本算例的研究對(duì)象的是在軸線方向上沒有曲率變化的光滑圓柱.從表2 中可以看出，增大速度v的量級(jí)相比于其他變量過小，因此為了增大速度v對(duì)神經(jīng)網(wǎng)絡(luò)權(quán)重更新的影響，文中對(duì)速度v的數(shù)據(jù)誤差乘以修正權(quán)重100，使其與速度u，速度w處于同一量級(jí)，因此神經(jīng)網(wǎng)絡(luò)損失函數(shù)為

通過神經(jīng)網(wǎng)絡(luò)的反向傳播與梯度下降算法最小化損失函數(shù)，更新神經(jīng)網(wǎng)絡(luò)中的權(quán)重，直至神經(jīng)網(wǎng)絡(luò)收斂.綜合上述內(nèi)容，神經(jīng)網(wǎng)絡(luò)的訓(xùn)練流程如圖9所示.

圖9 神經(jīng)網(wǎng)絡(luò)訓(xùn)練流程圖Fig.9 Flow chart of neural network training

文中神經(jīng)網(wǎng)絡(luò)基于全連接神經(jīng)網(wǎng)絡(luò)搭建，包含10 個(gè)隱藏層，每層包含500 個(gè)神經(jīng)元.本次訓(xùn)練采用了Adam 優(yōu)化器，學(xué)習(xí)率設(shè)為10-3，預(yù)設(shè)訓(xùn)練時(shí)間為40h，微分的計(jì)算由tensorflow 中的自動(dòng)微分完成.激活函數(shù)選取sigmoid 函數(shù).

2.3 神經(jīng)網(wǎng)絡(luò)權(quán)重選擇

文中首先以雷諾數(shù)Re=3 900 工況為例，為了便于誤差分析，文中將神經(jīng)網(wǎng)絡(luò)損失函數(shù)分4 部分輸出，分別為基于流場(chǎng)數(shù)據(jù)的流場(chǎng)誤差loss1；基于表面壓力數(shù)據(jù)的界面誤差loss2；基于物理方程的方程誤差loss3以及三項(xiàng)誤差之和的總誤差loss_all，其計(jì)算式為

流域誤差loss1:

界面誤差loss2:

方程誤差loss3:

其中，L1至L7的計(jì)算方式見式(4)～(5).

為了確定平衡權(quán)重α的最佳取值，文中首先將α取1，即數(shù)據(jù)項(xiàng)和方程項(xiàng)對(duì)神經(jīng)網(wǎng)絡(luò)權(quán)重更新的影響比例相同，此時(shí)流域誤差為loss1=L5+L6+L7.圖10 為α=1 時(shí)神經(jīng)網(wǎng)絡(luò)四類誤差隨訓(xùn)練輪數(shù)（epoch）的輸出曲線可以看出，四項(xiàng)誤差均隨著迭代次數(shù)的增加而明顯減小.方程誤差loss3的最終收斂數(shù)值在3 類誤差（不包括總誤差）中最大，穩(wěn)定在10-4.2左右，此時(shí)方程誤差對(duì)總誤差的影響較大，導(dǎo)致在迭代過程中，神經(jīng)網(wǎng)絡(luò)很難從數(shù)據(jù)項(xiàng)中獲得訓(xùn)練信息.因此通過增大損失函數(shù)中平衡權(quán)重系數(shù)α增大流域誤差對(duì)神經(jīng)網(wǎng)絡(luò)訓(xùn)練的影響.對(duì)于權(quán)重的選擇參考文獻(xiàn)[25]結(jié)果，文中重點(diǎn)關(guān)注α取1，100，300，500 這幾種情況，其中方程誤差曲線如圖11 所示.

圖10 神經(jīng)網(wǎng)絡(luò)誤差輸出曲線（α=1）Fig.10 Error output curve of neural network (α=1)

圖11 不同權(quán)重系數(shù)下神經(jīng)網(wǎng)絡(luò)誤差曲線對(duì)比圖Fig.11 Comparison of neural network error curves under different weight coefficients

可以看出，隨著權(quán)重系數(shù)的不斷增大，方程誤差的收斂值在逐漸降低，當(dāng)權(quán)重等于300 和500 時(shí)，方程誤差收斂值均為10-4.9.不同點(diǎn)在于α=500 時(shí)誤差曲線波動(dòng)較大穩(wěn)定性較差.

為了進(jìn)一步驗(yàn)證權(quán)重系數(shù)對(duì)預(yù)測(cè)精度的影響，在圖4 所示的流場(chǎng)區(qū)域中隨機(jī)選取100 個(gè)點(diǎn)，計(jì)算在各個(gè)權(quán)重系數(shù)下流場(chǎng)各變量與真實(shí)值之間的平均相對(duì)誤差，最終結(jié)果見表5.

表5 各權(quán)重系數(shù)下流場(chǎng)變量預(yù)測(cè)平均相對(duì)誤差Tab.5 The prediction average relative error of flow field variables under each weight coefficient

可以看出，隨著α從1 增加到500，速度v相對(duì)誤差下降明顯，考慮α=500 時(shí)誤差曲線有明顯波動(dòng)，因此平衡權(quán)重α的取300 最合適.

3 預(yù)測(cè)結(jié)果分析

圖12 展示了圓柱軸向中心平面上神經(jīng)網(wǎng)絡(luò)重構(gòu)結(jié)果.從訓(xùn)練數(shù)據(jù)中選取第500 步，t=12.35s 的數(shù)值模擬結(jié)果（此時(shí)間步下的流場(chǎng)數(shù)據(jù)并不參與神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過程過程），結(jié)果均已無量綱化處理.

圖12 神經(jīng)網(wǎng)絡(luò)重構(gòu)結(jié)果對(duì)比（t=12.35 s）Fig.12 Comparison of prediction results of neural network (t=12.35 s)

為了定量表示神經(jīng)網(wǎng)絡(luò)的重構(gòu)誤差，圖12(c)給出了相對(duì)誤差云圖.其計(jì)算方法為

式中：α為流場(chǎng)中的變量；i為第i個(gè)樣本點(diǎn)對(duì)應(yīng)編號(hào).可以看出，神經(jīng)網(wǎng)絡(luò)重構(gòu)流場(chǎng)與數(shù)值模擬結(jié)果吻合度較高.

圖13 給出了數(shù)值模擬和神經(jīng)網(wǎng)絡(luò)重構(gòu)流場(chǎng)的速度u局部流線圖，可以看出神經(jīng)網(wǎng)絡(luò)可以精準(zhǔn)捕捉流場(chǎng)流動(dòng)方向.

圖13 uNS 和uPINN 局部流線圖對(duì)比Fig.13 Comparison of local streamlines of uNS and uPINN

圖14 為神經(jīng)網(wǎng)絡(luò)重構(gòu)流場(chǎng)的升阻力系數(shù)變化對(duì)比圖.可以看出，神經(jīng)網(wǎng)絡(luò)重構(gòu)結(jié)果僅能較好地吻合升阻力系數(shù)變化周期，這與神經(jīng)網(wǎng)絡(luò)訓(xùn)練時(shí)參與的數(shù)據(jù)有關(guān).文中選取50%的稀疏時(shí)間參與訓(xùn)練，可以充分體現(xiàn)出各變量的周期變化規(guī)律.但由于文中采用邊界數(shù)據(jù)浸入法處理固/液表面，難以獲得準(zhǔn)確數(shù)據(jù)，導(dǎo)致神經(jīng)網(wǎng)絡(luò)難準(zhǔn)捕捉渦的初生與脫落位置.

圖14 升阻力系數(shù)對(duì)比Fig.14 Comparison of lift and drag coefficients

表6 為神經(jīng)網(wǎng)絡(luò)重構(gòu)流場(chǎng)兩個(gè)周期內(nèi)的渦量變圖，其中tref=10 s.可以看出，神經(jīng)網(wǎng)絡(luò)能精準(zhǔn)捕捉渦的初生與脫落位置.在圖中用兩條虛線標(biāo)注出一個(gè)正漩渦從初生到脫落的完整階段，可以看出在二者結(jié)果中，渦核脫落頻率和運(yùn)動(dòng)速度大致相等.

表6 尾渦演化預(yù)測(cè)對(duì)比（Re=100）Tab.6 Comparison of tail vortex evolution prediction (Re = 100)

為了精準(zhǔn)對(duì)比神經(jīng)網(wǎng)絡(luò)對(duì)渦量的預(yù)測(cè)精度，圖15為t=0.5tref時(shí)刻渦量相對(duì)誤差云圖.可以看出渦量誤差集中的區(qū)域也呈現(xiàn)出一定規(guī)律性，集中在速度變化劇烈的旋渦結(jié)構(gòu)附近，分別為紅色區(qū)域的高誤差區(qū)①，相對(duì)誤差最高可以達(dá)到30%左右.灰色區(qū)域的低誤差區(qū)②，③，④，相對(duì)誤差最高在15%左右.高誤差區(qū)位于正漩渦的渦脫落位置和反漩渦的初生位置之間，并且越靠近中心位置，誤差越大.低誤差區(qū)位于正漩渦和反漩渦中間.結(jié)合圖16 的渦量等值線圖可以看出，在高誤差區(qū)，等值線密集，渦量梯度下降最快，在低誤差區(qū)，渦量趨近于0.高誤差區(qū)和低誤差區(qū)的存在，和物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)的擬合方式有關(guān)，一方面，神經(jīng)網(wǎng)絡(luò)的數(shù)據(jù)項(xiàng)與正則化項(xiàng)在訓(xùn)練過程中會(huì)相互干擾，當(dāng)正則化項(xiàng)梯度過大時(shí)，神經(jīng)網(wǎng)絡(luò)會(huì)忽略數(shù)據(jù)項(xiàng)的信息，從而難以擬合微分方程的特解.另一方面，在低誤差區(qū)中，渦量的變化梯度趨近于0 且數(shù)值較小，越靠近兩渦中心位置渦量的數(shù)值就越小，給神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)帶來困難，因此誤差較大.但總體來看，除了大梯度變化區(qū)域，神經(jīng)網(wǎng)絡(luò)可以較為準(zhǔn)確地預(yù)測(cè)流場(chǎng)渦量變化.

圖15 t=0.5 tref 渦量相對(duì)誤差云圖Fig.15 Relative error cloud chart of vorticity (t=0.5tref)

圖16 t=0.5 tref 數(shù)值模擬二維渦量等值線圖Fig.16 Numerical simulation of two-dimensional vorticity contour map(t=0.5 tref)

為了進(jìn)一步觀察重構(gòu)流場(chǎng)中的渦結(jié)構(gòu)，圖17 給出了基于準(zhǔn)則下的渦識(shí)別圖.可以看出，數(shù)值方法和神經(jīng)網(wǎng)絡(luò)對(duì)渦管中心位置的捕捉幾乎相同，但神經(jīng)網(wǎng)絡(luò)存在部分特征丟失.

圖17 渦結(jié)構(gòu)對(duì)比圖（Q=0.2）Fig.17 Comparison of Vortex Structure (Q=0.2)

為探討神經(jīng)網(wǎng)絡(luò)對(duì)不同流態(tài)流場(chǎng)的重構(gòu)能力，文中分別取雷諾數(shù)為100，3 900，20 000，100 000 四種工況，尾渦對(duì)應(yīng)4 種不同的脫落狀態(tài)，分別為初始脫落，穩(wěn)定脫落，向湍流過渡的亞穩(wěn)定脫落以及發(fā)展成湍流的不穩(wěn)定脫落，下圖為不同雷諾數(shù)下神經(jīng)網(wǎng)絡(luò)重構(gòu)流場(chǎng)的渦量對(duì)比結(jié)果如表7 所示.

表7 不同雷諾數(shù)下二維渦量對(duì)比圖Tab.7 Comparison of two-dimensional vorticity at different Reynolds numbers

當(dāng)雷諾數(shù)Re=100 時(shí)，雷諾數(shù)很小，首次產(chǎn)生渦脫落現(xiàn)象，但是漩渦范圍較小.可以看出，神經(jīng)網(wǎng)絡(luò)對(duì)渦量的重構(gòu)結(jié)果與數(shù)值模擬基本一致，其中包括漩渦的中心位置，脫落位置以及漩渦的長(zhǎng)度.當(dāng)雷諾數(shù)增加到Re=20 000 時(shí)，流場(chǎng)處于臨界流態(tài)，尾渦開始逐漸表現(xiàn)出非定常特性，流動(dòng)較為紊亂，脫落頻率和渦的長(zhǎng)度不再固定，同時(shí)在流場(chǎng)中開始出現(xiàn)小尺度的離散渦.神經(jīng)網(wǎng)絡(luò)的重構(gòu)結(jié)果顯示，物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)在此流態(tài)下已無法精準(zhǔn)重構(gòu)流場(chǎng)渦的大小以及渦核的位置，高誤差區(qū)的預(yù)測(cè)相對(duì)誤差甚至高達(dá)40%左右.當(dāng)雷諾數(shù)增加至Re=100 000 時(shí)，此時(shí)圓柱尾流完全演變?yōu)橥牧?，流?chǎng)處于跨臨界流態(tài).尾渦以大量的小尺度和離散渦為主，繞流極為凌亂、不規(guī)則，此時(shí)流場(chǎng)表現(xiàn)出強(qiáng)非線性性和隨機(jī)性，使得神經(jīng)網(wǎng)絡(luò)重構(gòu)精度進(jìn)一步下降.在前文中已經(jīng)論述，神經(jīng)網(wǎng)絡(luò)無法捕捉到流場(chǎng)中的小尺度渦結(jié)構(gòu)，而在此流態(tài)下流場(chǎng)以小尺度渦為主，因此當(dāng)雷諾數(shù)Re=100 000 時(shí)，神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)結(jié)果中沒有太多明顯的渦結(jié)構(gòu).對(duì)比四種流態(tài)下的結(jié)果可以看出，當(dāng)雷諾數(shù)增大時(shí)，由于渦結(jié)構(gòu)逐漸細(xì)小化和離散化，導(dǎo)致神經(jīng)網(wǎng)絡(luò)重構(gòu)精度逐漸降低.為了展示神經(jīng)網(wǎng)絡(luò)在雷諾數(shù)Re=20 000 和Re=100 000 兩種工況下對(duì)渦結(jié)構(gòu)的精細(xì)捕捉，采用Q準(zhǔn)則對(duì)兩種工況下的渦結(jié)構(gòu)進(jìn)行識(shí)別，如圖18 所示.

圖18 不同流態(tài)下神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)流場(chǎng)渦結(jié)構(gòu)對(duì)比（Q 準(zhǔn)則）Fig.18 Comparison of vortex structure predicted by neural network under different flow regimes(Q criterion)

可以看出，Re=20 000 和Re=100 000 兩種工況下，數(shù)值模擬結(jié)果中流場(chǎng)的多尺度效應(yīng)明顯， 神經(jīng)網(wǎng)絡(luò)雖然無法捕捉小尺度渦結(jié)構(gòu)，但是卻可以預(yù)測(cè)出圖中已由黑色虛線框標(biāo)出的大尺度渦的位置和形狀.

圖19 給出了不同雷諾數(shù)下的一個(gè)周期tref=10 s內(nèi)的圓柱表面時(shí)均升阻力系數(shù)的相對(duì)誤差，從圖中可以看出，當(dāng)雷諾數(shù)Re=100 和3 900 時(shí)，圓柱的尾流區(qū)域處于穩(wěn)定渦渦脫落狀態(tài)，對(duì)升阻力系數(shù)的重構(gòu)誤差均在25%以內(nèi)，雷諾數(shù)Re=3 900 時(shí)均升阻力系數(shù)會(huì)略有增大，但變化較小，均在5%以內(nèi).隨著雷諾數(shù)的增大，流動(dòng)結(jié)構(gòu)逐漸復(fù)雜，重構(gòu)誤差逐漸增大，在雷諾數(shù)Re=100 000 時(shí)，時(shí)均升阻力系數(shù)的預(yù)測(cè)誤差達(dá)到55%左右，這與神經(jīng)網(wǎng)絡(luò)對(duì)大梯度變化不敏感有關(guān)，在上文已詳細(xì)闡述.可以看出，隨著流場(chǎng)結(jié)構(gòu)逐漸復(fù)雜，物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)對(duì)流場(chǎng)水動(dòng)力特性的重構(gòu)能力在逐漸降低.

圖19 不同雷諾數(shù)下升阻力系數(shù)時(shí)均相對(duì)誤差Fig.19 Average relative error of lift-drag coefficient at different Reynolds numbers

4 結(jié)束語

本章首先介紹了物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)在流體力學(xué)中的應(yīng)用情況.接著發(fā)展了基于Fortran 自編程序的圓柱繞流數(shù)值求解方法，并將其作為神經(jīng)網(wǎng)絡(luò)的輸入數(shù)據(jù)集.以PINN 框架為基礎(chǔ)構(gòu)建了耦合流場(chǎng)物理信息的神經(jīng)網(wǎng)絡(luò)模型，建立了物理信息驅(qū)動(dòng)的流場(chǎng)結(jié)構(gòu)重構(gòu)方法，同時(shí)根據(jù)具體流場(chǎng)結(jié)構(gòu)針對(duì)性提出神經(jīng)網(wǎng)絡(luò)的數(shù)據(jù)采樣方式和優(yōu)化策略.最后探討了該方法對(duì)不同流態(tài)的流場(chǎng)重構(gòu)能力，揭示了神經(jīng)網(wǎng)絡(luò)的重構(gòu)誤差機(jī)理，具體內(nèi)容如下：

①通過輸出誤差的對(duì)比將神經(jīng)網(wǎng)絡(luò)確定為雙權(quán)重修正模型.確定數(shù)據(jù)項(xiàng)平衡權(quán)重α取300 最為合適，垂向速度v誤差修正權(quán)重λ取100 最為合適.

②對(duì)比分析物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)的重構(gòu)精度.結(jié)果證明，神經(jīng)網(wǎng)絡(luò)利用稀疏數(shù)據(jù)和流動(dòng)控制方程即可實(shí)現(xiàn)對(duì)流場(chǎng)的精準(zhǔn)重構(gòu).通過對(duì)比渦量云圖，證實(shí)了神經(jīng)網(wǎng)絡(luò)可以精準(zhǔn)重構(gòu)流場(chǎng)渦結(jié)構(gòu)，但在渦量大梯度變化且數(shù)值較小處預(yù)測(cè)能力欠佳.

③對(duì)比分析了物理信息驅(qū)動(dòng)型神經(jīng)網(wǎng)絡(luò)對(duì)不同流態(tài)流場(chǎng)的重構(gòu)能力.針對(duì)尾渦初始脫落（Re=100），尾渦穩(wěn)定脫落（Re=3 900），尾渦亞穩(wěn)定脫落（Re=20 000）以及尾渦不穩(wěn)定脫落（Re=100 000）的四種流動(dòng)狀態(tài)，探討了神經(jīng)網(wǎng)絡(luò)在不同雷諾數(shù)下對(duì)單相不可壓縮流場(chǎng)的重構(gòu)情況.在Re=100 時(shí)，數(shù)值模擬和神經(jīng)網(wǎng)絡(luò)模型對(duì)流場(chǎng)結(jié)構(gòu)的結(jié)果基本一致.在Re=20 000 和Re=100 000 時(shí)，神經(jīng)網(wǎng)絡(luò)受限于流場(chǎng)多尺度效應(yīng)，僅能較為精準(zhǔn)的預(yù)測(cè)剪切層破裂前的大尺度渦結(jié)構(gòu).同時(shí)隨著雷諾數(shù)增大，神經(jīng)網(wǎng)絡(luò)對(duì)流場(chǎng)升阻力系數(shù)預(yù)測(cè)能力逐漸降低.

（責(zé)任編輯：孫竹鳳）