馮唐思捷 梁 偉

(北京航空航天大學(xué)航空科學(xué)與工程學(xué)院,北京 100191)

引言

隨著計算機硬件的更新迭代,數(shù)值計算資源呈現(xiàn)了爆炸式增長的態(tài)勢.快速降低的計算成本催生了一大批使用大數(shù)據(jù)進行計算的新型數(shù)值方法,例如人工神經(jīng)網(wǎng)絡(luò)和深度學(xué)習(xí)等技術(shù)已經(jīng)在圖像識別[1]和認知科學(xué)[2]等方面取得了突破性進展.但在工程和物理領(lǐng)域,這些新興方法的應(yīng)用推進并不順利.其中一個重要的原因是這些領(lǐng)域的數(shù)據(jù)集的獲取成本高,稀疏性和不確定性大,純粹依靠數(shù)據(jù)驅(qū)動的方法效率低且精度不高.因此,大部分工程和物理問題采用的求解方法都是基于物理模型的正向求解算法.物理模型通常由偏微分方程組的形式刻畫或者描述,例如Navier-Stokes[3]方程、Allen-Cahn[4]方程和Schr?dinger[5]方程等.基于求解偏微分方程的傳統(tǒng)求解方法雖然能夠求解正向問題,但在相應(yīng)的逆問題、高維問題和非線性問題上具有計算成本高,人機交互頻繁,技術(shù)突破難度大的缺點.

近年,同時利用大量計算數(shù)據(jù)和物理規(guī)律來求解偏微分方程組所代表的物理模型成為了一個可能的工程技術(shù)突破口.一方面,在Cybenko[6]和Hornik[7]證明的神經(jīng)網(wǎng)絡(luò)萬能逼近定理的前提下,擁有多于一個隱藏層的神經(jīng)網(wǎng)絡(luò)可以將任意從向量空間Rn到另一個向量空間 Rn進行投影的Borel 函數(shù)的擬合誤差達到足夠小.另一方面,得益于編程語言的發(fā)展,如今神經(jīng)網(wǎng)絡(luò)的主流搭建平臺,例如Tensorflow 和Pytorch 等,都廣泛使用自動微分 (auto-differentiation)[8]來自動計算點與點之間的多階梯度.自動微分使得方程中的偏微分項可以融入到神經(jīng)網(wǎng)絡(luò)的損失函數(shù)設(shè)計當中.這樣,物理模型中包含的物理知識作為一項信息輸入到神經(jīng)網(wǎng)絡(luò)的優(yōu)化過程中,有這項信息輸入的神經(jīng)網(wǎng)絡(luò)就脫離了僅僅將輸入輸出數(shù)據(jù)進行匹配的純數(shù)據(jù)驅(qū)動范疇.這樣綜合了數(shù)據(jù)驅(qū)動和物理信息兩方面的神經(jīng)網(wǎng)絡(luò)被稱之為物理信息神經(jīng)網(wǎng)絡(luò)(PINN).

PINN 最先由Raissi 等[9-10]于2017 年提出,并在2019 年證明對多種偏微分方程有效[11].PINN繼承了深度神經(jīng)網(wǎng)絡(luò)(DNN)的函數(shù)逼近能力,其整體架構(gòu)也與深度神經(jīng)網(wǎng)絡(luò)類似:采用前饋神經(jīng)網(wǎng)絡(luò)作為基本逼近單元,將函數(shù)的自變量作為輸入,在輸出層輸出給定函數(shù)或給定方程的近似解.這種將偏微分方程中的信息內(nèi)嵌到損失函數(shù)的做法使得損失函數(shù)在優(yōu)化的過程中遵循了相關(guān)的物理定律.與傳統(tǒng)算法相比,PINN 不僅不需要任何網(wǎng)格劃分的前處理,也不需要遵循任何人為的先驗假設(shè)[12-13].簡潔的網(wǎng)絡(luò)架構(gòu)使得PINN 能夠?qū)Υ蟛糠制⒎址匠探M進行有效求解,并可以有效攻克一些常規(guī)數(shù)值方法中的難點,例如維數(shù)災(zāi)難和逆問題求解中的參數(shù)識別等[14-19].

在彈性力學(xué)中,由于控制方程大多數(shù)為高階偏微分方程,在復(fù)雜幾何區(qū)域下求解困難.近年來,學(xué)界注意到PINN 的良好特性可用于求解固體力學(xué)中的一些實際問題并驗證了其有效性.Tao 等[20]使用PINN 對薄壁圓筒殼在線性屈曲失穩(wěn)下的相對于理論解的塌陷因數(shù)(knockdown factor)進行了研究.結(jié)果證明,PINN 可以在提供更少的標簽數(shù)據(jù)的情況下得到與傳統(tǒng)ANN 表現(xiàn)相似的神經(jīng)網(wǎng)絡(luò),其預(yù)測值不高于標簽數(shù)據(jù)中的實驗數(shù)據(jù),便于提供結(jié)構(gòu)真實屈曲載荷的保守預(yù)測[20].Li 等[21]使用PINN 方法基于偏微分方程形式和Rayleigh-Ritz 能量條件構(gòu)建損失函數(shù)求解了薄板的彎曲問題,結(jié)果和有限元方法的精度相仿.Haghighat 等[22]則在他們的研究中討論了PINN 作為代理模型,利用等幾何分析和有限元解作為標簽數(shù)據(jù)來求解線彈性力學(xué)和彈塑性力學(xué)問題的可能性.Bastek 等[23]則研究了在非歐流型下求解基于Naghdi 板殼理論的變形問題,并討論了待求解方程強弱形式對神經(jīng)網(wǎng)絡(luò)的影響.Yan 等[24]利用PINN 方法和極限學(xué)習(xí)機(extreme learning)研究了復(fù)合材料薄壁結(jié)構(gòu)的靜力和動力問題,并將神經(jīng)網(wǎng)絡(luò)解與這些問題的Navier解進行了對比.雖然PINN 在求解實際問題上仍然處于事實上的起步階段,但其相對于傳統(tǒng)數(shù)值方法擁有不可忽視的潛力.

工程實際中,薄壁結(jié)構(gòu)在工作時主要受軸向壓縮載荷,在這種載荷的作用下,結(jié)構(gòu)容易發(fā)生屈曲失穩(wěn)從而達不到預(yù)期的承載能力[25-26].薄壁結(jié)構(gòu)的屈曲和后屈曲控制方程為復(fù)雜的非線性耦合方程,一般情況下難以得到顯式表達的位移場理論解,而在數(shù)值方法中,也一般采用有限元方法進行求解其弱形式積分方程.之前列舉的一些工作雖然對于固體力學(xué)中的偏微分方程求解提供了新的思路,但求解的力學(xué)問題多局限于線性理論,而對于非線性問題則鮮有涉獵.而屈曲問題是否考慮非線性效應(yīng)將對結(jié)果產(chǎn)生很大的影響.由于無論求解手段如何,屈曲問題和后屈曲問題都可以轉(zhuǎn)化為對其控制方程求擬合解,故而本文將結(jié)合物理信息神經(jīng)網(wǎng)絡(luò)對結(jié)構(gòu)的非線性屈曲問題,基于求強形式控制方程的擬合解來發(fā)展新的解法.本文采用PINN 方法對典型薄壁板殼結(jié)構(gòu)的屈曲和后屈曲問題進行研究,針對不同的邊界條件和載荷形式進行結(jié)構(gòu)屈曲載荷和屈曲平衡位移模態(tài)的研究.本文利用Lu 等[27]建立的開源Python 庫DeepXDE 編寫了3 種典型板殼結(jié)構(gòu)的求解代碼,實搭建的實驗平臺為Python 3.9.16,所有的神經(jīng)網(wǎng)絡(luò)運算均在GPU 上利用NVIDIA 公司開發(fā)的CUDA 運行.PINN 解在最后與有限元解進行了對比以驗證本文方法的有效性.

1 薄板殼非線性屈曲控制方程

長寬均為a,厚度為h的方板,其彈性模量為E,材料泊松比為 ν,受面內(nèi)載荷Px,Py,Pxy作用,引入Airy 應(yīng)力函數(shù)F,將3 種面內(nèi)載荷統(tǒng)一表示為F的偏導(dǎo)數(shù),即

以W?,W分別表示初始缺陷帶來的撓度和因為載荷施加得到的中面撓度,中面x,y方向上的位移分別記為U,V.由von-Kármán Formulation,可知其幾何方程為[28]

式中,上標0 代表中面的應(yīng)變分量和曲率分量.以薄板中面為研究對象,其平衡方程可寫為

假設(shè)材料為各向同性線彈性的,并考慮協(xié)調(diào)方程

以及物理方程

將中面位移全部使用撓度W?,W和Airy 應(yīng)力函數(shù)F表達,則可以得到以下的平板屈曲微分控制方程[29]

邊界條件可分為3 類:簡支、固支與自由,其對應(yīng)的表達式分別為

簡支

固支

自由邊

為了減少神經(jīng)網(wǎng)絡(luò)需要擬合的邊界條件,對于固支邊和簡支邊的位移連續(xù)條件,可以在用對神經(jīng)網(wǎng)絡(luò)的擬合過程中使用硬邊界約束(hard constraint).例如,若某一薄板為四邊簡支,則可以對神經(jīng)網(wǎng)絡(luò)所代表的非線性映射做變換,即令

這樣非線性映射則可自動滿足力學(xué)模型的位移邊界條件.對于邊界條件(8),需要融入到損失函數(shù)的邊界條件減少了一項,其余的邊界條件也可以構(gòu)建類似的試函數(shù)來減少損失函數(shù)項數(shù).

2 神經(jīng)網(wǎng)絡(luò)模型

2.1 網(wǎng)絡(luò)框架

本文采用的PINN 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖1 所示.

圖1 PINN 求解網(wǎng)絡(luò)框架圖Fig.1 A framework of using PINN network as PDE solver

此全連接神經(jīng)網(wǎng)絡(luò)的基本參數(shù)為La×Nu,前者為隱藏層層數(shù),后者為單個隱藏層神經(jīng)元數(shù).神經(jīng)網(wǎng)絡(luò)的輸入為幾何域坐標 (x,y),經(jīng)過神經(jīng)網(wǎng)絡(luò)輸出的擬合解為撓度與應(yīng)力函數(shù)該問題所對應(yīng)的非線性映射可以在數(shù)學(xué)上表示為

式中,參數(shù)Wi,bi分別為每一層神經(jīng)元組成的權(quán)重矩陣和偏置向量,每一層的非線性激活函數(shù)由 σi表示,它們構(gòu)成了需要訓(xùn)練的網(wǎng)絡(luò)參數(shù) Θ.

對于激活函數(shù) σi,根據(jù)Jagtap 等[30]的工作,使用逐層自適應(yīng)的激活函數(shù)(locally adaptive activation function,LAAF)可以進一步提高PINN 的表現(xiàn).在第3 章的算例中,均采用文獻[30]中的自適應(yīng)激活函數(shù),其在DEEPXDE 開源庫中可通過命令LAAF調(diào)取,其基礎(chǔ)的激活函數(shù)為sigmoid 函數(shù),即

2.2 誤差函數(shù)與算法流程

一般來說,PINN 的損失函數(shù) Loss 由3 個部分組成.如果由均方誤差(mean square error,MSE)來衡量神經(jīng)網(wǎng)絡(luò)的擬合誤差,則3 個部分分別是代表偏微分控制方程的誤差MSEP,邊界條件誤差MSEB,以及初值條件誤差MSEI.PINN 方法的待優(yōu)化損失函數(shù)Loss可以帶權(quán)重的方式表示為

MSEP作為整個求解域上的控制方程的損失函數(shù),其首先由神經(jīng)網(wǎng)絡(luò)在整個求解域上預(yù)先布置的采樣點上得到預(yù)測撓度與應(yīng)力函數(shù)再將預(yù)測解代入控制方程求出殘差.設(shè)在域內(nèi)總共布置了NP個采樣點,即有

MSEB作為邊界上的損失函數(shù),其計算方法與MSEP相似.以四邊簡支的邊界條件為例,由于已經(jīng)在神經(jīng)網(wǎng)絡(luò)的輸入端引入了變換層,位移的連續(xù)條件被滿足,故只需要考慮彎矩條件

MSEI作為表示初值的誤差,在屈曲方程中由于不含時,因此初值項為0.

確定PINN 需要優(yōu)化的損失函數(shù)后,采用基于Python 的機器學(xué)習(xí)框架Tensorflow 中的基于計算圖的自動微分(auto differentiation)就可以計算這些函數(shù)中微分項.

損失函數(shù)Loss的表達式中,[wP,wB,wI]分別為3 誤差項對應(yīng)的權(quán)重值向量.在實驗中,作者發(fā)現(xiàn)可以將所有損失函數(shù)項的數(shù)量級調(diào)整至接近一致,這樣做可以提高PINN 算法的訓(xùn)練精度,但具體權(quán)重仍然需要具體問題具體確定.一種辦法是在確定隨機種子后,不使用加權(quán)策略對模型進行迭代數(shù)僅數(shù)百次的預(yù)訓(xùn)練,獲取其損失項大致的數(shù)量級,再確定權(quán)重.流程框圖如下所示.

本文第3 章中的算例均采用此算法確定損失函數(shù)權(quán)重值,確定的損失函數(shù)權(quán)重表可見附錄B.

由于屈曲問題的模態(tài),屈曲路徑和屈曲載荷相對應(yīng),因此必須使用弧長法等迭代方法額外獲得每一步的增量.具體可將當前神經(jīng)網(wǎng)絡(luò)求出的當前載荷的收斂位移Wepoch作為下一步的屈曲幾何缺陷W?代入屈曲控制方程求解屈曲附加撓度,這樣在流程中,PINN 作為偏微分方程求解器使用,而外層使用弧長法進行控制,其迭代求解偏微分控制方程算法流程框圖如圖2 所示.

圖2 以弧長法為外層控制的物理信息神經(jīng)網(wǎng)絡(luò)屈曲/后屈曲求解框圖Fig.2 The framework of PINN as PDE sovler towards buckling/postbuckling problems governed by arcs method

弧長法在本文中作為PINN 的外層控制算法,當神經(jīng)網(wǎng)絡(luò)模型在擬合時達到了一定收斂度,便使用弧長法結(jié)束當前的神經(jīng)網(wǎng)絡(luò)迭代并計算下一步載荷增量,再將這一步的位移作為新的缺陷項替換神經(jīng)網(wǎng)絡(luò)中的缺陷項,最后利用神經(jīng)網(wǎng)絡(luò)進行下一步的學(xué)習(xí).

使用弧長法求解PINN 得出的微分方程擬合解的增量并進行下一步求解的方法如算法2 所示.

在本文的流程中,預(yù)訓(xùn)練和正式訓(xùn)練時,默認求解域內(nèi)部訓(xùn)練樣本點NP=200,邊界點NB=100.預(yù)訓(xùn)練時的訓(xùn)練次數(shù)iterpre取500 次,正式訓(xùn)練的最大訓(xùn)練次數(shù)itermax取40000 次.LPF誤差容限 ε=0.01.對于測試集,另取樣本點NT=100,兩種樣本點均采用拉丁超采樣(LHS)方法生成.模型學(xué)習(xí)率初始為lr0=0.0001,設(shè)置學(xué)習(xí)率衰減為逆時衰減策略,即

衰減策略為每T=1000 輪衰減β=0.2.對于網(wǎng)絡(luò)模型的初始化,采用Glorot 正態(tài)分布初始化器對參數(shù)進行初始化.

3 數(shù)值算例

3.1 受面內(nèi)載荷的缺陷方板的屈曲和后屈曲

首先,對于一均質(zhì)薄板的線性與非線性屈曲進行研究以驗證PINN 的有效性.所研究的各向同性矩形板長寬尺寸均為a=1 m,厚度h=2 mm.材質(zhì)為6 系鋁合金,主要的力學(xué)性能為:彈性模量E=70 GPa,泊松比 ν=0.3,屈服應(yīng)力280 MPa.設(shè)立3 種邊界條件.分別為:①四邊簡支(SSSS);②對邊簡支對邊固支(SCSC);③一對對邊簡支,另一對邊一邊固支一邊自由(SCSF).其受載荷形式為在x=0,x=a兩條簡支對邊上受邊載荷Px.

先考慮其不帶缺陷時,在不同的邊界條件下所得到的屈曲模態(tài)和屈曲載荷,對于屈曲模態(tài),由于其是滿足方程邊界條件的非平凡解,因此可以直接使用以載荷比例因子(LPF)為控制的迭代過程.由于PINN 本身的使泛化誤差最小化的特質(zhì),如果其在迭代過程中不手動使載荷增加,而是直接對其設(shè)定估計的屈曲載荷值并讓神經(jīng)網(wǎng)絡(luò)自身逼近,往往會得到平凡解和根據(jù)網(wǎng)絡(luò)自身第一次迭代后的非0 載荷.

由PINN 方法計算的第一階臨界載荷和有限元結(jié)果的對比如表1 所示.

表1 第一階屈曲模態(tài)對應(yīng)的外載對比Table 1 The buckling load of the first buckling mode

可以看出,PINN 對于臨界屈曲載荷能進行較好的預(yù)測,其和有限元相比,相對誤差可以控制在5%以內(nèi).

對于薄板屈曲問題,由于其屈曲撓度分布極其容易受到初缺陷的影響,故將缺陷使用雙三角級數(shù)表示,即

并采用高斯隨機策略以模擬工程中具有初缺陷的薄板,其三角級數(shù)系數(shù)滿足相互獨立的正態(tài)分布

將式(18)作為缺陷代入控制方程,在接下來的算例中,取m=n=4 生成缺陷,并固定隨機種子使得到的初缺陷表達式能夠輸入有限元軟件內(nèi)進行復(fù)現(xiàn),最終得到的結(jié)果和有限元結(jié)果的對比如圖3 所示.

圖3 帶有缺陷的均質(zhì)薄板的非線性屈曲位移云圖對比 (上:PINN 結(jié)果,下:FEM 結(jié)果)Fig.3 The comparison of nonlinear buckling with deficiency (up:PINN results,down:FEM results)

綜上結(jié)果可以看出,PINN 能夠?qū)ν暾“宓木€性屈曲的模態(tài),第一階屈曲載荷和缺陷薄板非線性屈曲的后屈曲做出很好的預(yù)測.這證明采用弧長法進行屈曲控制并通過PINN 求解器計算屈曲問題是可行的.

3.2 受軸壓的缺陷薄壁圓筒殼

在本節(jié)中,使用PINN 方法針對一圓筒殼開展非線性后屈曲分析,圖4 所示的圓筒殼半徑R=1125 mm,高L=600 mm,厚度h=2 mm,材質(zhì)與上節(jié)的薄板相同.邊界條件為上下兩圓周不可滑動簡支,其他區(qū)域不做約束.

圖4 薄壁圓筒殼及坐標系Fig.4 Thin-walled shell and its coordinate system

對于此算例,圓柱殼還特別地有周期性條件

對于此算例,考慮由初缺陷帶來的幾何非線性,使用3.1 中的初缺陷生成辦法生成波形幾何缺陷,進行非線性屈曲的研究.加載方式采用弧長法進行迭代,將無缺陷的圓筒殼的臨界屈曲載荷作為基準,使用比例載荷因子進行加載控制.此外,為了規(guī)避數(shù)值計算中不同量級數(shù)據(jù)的舍入誤差,對于此方程,引入無量綱化參數(shù)

式中,W與W?分別是在屈曲作用下的撓度和因圓筒殼初始缺陷帶來的幾何初撓度,F為Airy 應(yīng)力函數(shù).將上述變換代入方程(6)中,得到對于薄壁圓筒殼的全局屈曲控制方程為

其中非線性項為

工程上通常關(guān)心前屈曲階段的最大載荷因子,而后屈曲階段的研究則主要以端部的最大縮短距離的變化趨勢為主.

使用算法2,將迭代跳出判斷條件改為

并將誤差容限 ε 設(shè)置為足夠小的正數(shù),即可捕捉到屈曲載荷達到下降段前的最大比例載荷因子.對于后屈曲路徑中最大縮短距離,取作為迭代終止條件,其表達式可用撓度表達為

兩種方法生成的位移云圖對比和不同網(wǎng)絡(luò)參數(shù)得到的最大比例載荷因子見圖5 和表2.

表2 帶缺陷的圓筒殼屈曲行為最大比例載荷因子(LPFmax)與無量綱撓度 (wmax) 誤差對比Table 2 The L PFmax and wmax of shell buckling with initial deficiency and their relative error

圖5 帶有缺陷的均質(zhì)圓筒殼在最大 L PFmax 處非線性屈曲位移云圖對比Fig.5 The comparison of nonlinear buckling displacement atLPFmax

圖5 帶有缺陷的均質(zhì)圓筒殼在最大 L PFmax 處非線性屈曲位移云圖對比 (續(xù))Fig.5 The comparison of nonlinear buckling displacement atLPFmax(continued)

從結(jié)果可以看出,神經(jīng)網(wǎng)絡(luò)方法,尤其是采用基于一階導(dǎo)數(shù)優(yōu)化的Adam 優(yōu)化器進行神經(jīng)網(wǎng)絡(luò)的訓(xùn)練時,對最大承載時的位移模式求解較好,其位移分布和最大位移大小均和有限元結(jié)果接近.由PINN方法求出的軸壓作用下后屈曲載荷-端部縮短曲線如圖6 所示.

圖6 后屈曲載荷-端部縮短曲線Fig.6 The end shrinkage plot of the cylindrical shell against LPF

3.3 圓柱筒殼的缺陷敏感度

以上兩節(jié)的研究表明,初始幾何缺陷的存在會對結(jié)構(gòu)在屈曲過程中能達到的最大屈曲因子產(chǎn)生較大的影響,且根據(jù)文獻[28],初始幾何缺陷對屈曲載荷的影響對后屈曲最小載荷的影響更大.

圖7 給出了在PINN 算法:有限元算法和Koiter理論解中對于此參數(shù)的圓筒殼缺陷敏感度的影響曲線.圖中y軸為缺陷圓筒殼最大比例因子相對完整圓筒殼最大比例因子的比值,x軸為缺陷大小相對于筒殼厚度的比值.為了和已有結(jié)論做對比,調(diào)整厚度使半徑與厚度比R/h=200.可以看出,PINN 計算出的缺陷敏感度要大于有限元的缺陷敏感度.值得注意的是,Koiter 理論由于使用同一個缺陷敏感度來描述缺陷敏感度,而沒有考慮到筒殼本身的參數(shù),對于此算例,其在經(jīng)典圓筒殼理論中的Batdorf 殼體參數(shù)

圖7 圓筒殼的缺陷敏感度影響曲線Fig.7 Comparison of deficiency sensitivity of cylindrical shells

約為60,屬中長殼范疇,結(jié)構(gòu)后屈曲對缺陷敏感,PINN 和有限元算法都能較好地體現(xiàn)這一性質(zhì).

4 神經(jīng)網(wǎng)絡(luò)參數(shù)對收斂性影響

4.1 損失函數(shù)權(quán)重

在PINN 的訓(xùn)練過程中,損失權(quán)重的選取會極大地影響到神經(jīng)網(wǎng)絡(luò)的訓(xùn)練成功與否.對于構(gòu)成復(fù)雜、組成部分數(shù)量級相差較大的損失函數(shù),若不對其進行基于某種策略的平衡,勢必會造成神經(jīng)網(wǎng)絡(luò)更偏重某一項損失函數(shù)的學(xué)習(xí),而忽視了其他損失函數(shù)項.此處采用算例1 中帶有缺陷的四邊簡支薄板來分析網(wǎng)絡(luò)參數(shù)的變化對計算過程中誤差函數(shù)的影響.圖8 展示了未使用基于第2 章的權(quán)重調(diào)整策略(即權(quán)重分量默認為1)和使用了基于預(yù)訓(xùn)練權(quán)重調(diào)整策略對四邊簡支薄板損失函數(shù)的訓(xùn)練影響.其中橫軸為訓(xùn)練周期數(shù),縱軸為當前訓(xùn)練周期(epoch)的損失函數(shù)Lossepoch與第0 次(未訓(xùn)練)損失函數(shù)值Loss0比值的對數(shù).

可以看出,更改權(quán)重對于此問題影響非常關(guān)鍵,采用了權(quán)重調(diào)整策略的PINN 模型對于損失函數(shù)的優(yōu)化更好.同時為了排除PINN 撒點策略對結(jié)果產(chǎn)生的隨機性影響,本章的算例均采用隨機種子進行了多次實驗.本章中曲線圖由其平均值繪制,方差表可見附錄C.

4.2 隱藏層層數(shù) La 與神經(jīng)元數(shù)目Nu

在實際求解過程中,隱藏層層數(shù)和每層的神經(jīng)元個數(shù)選取會極大地影響到神經(jīng)網(wǎng)絡(luò)的訓(xùn)練效率和收斂性.在神經(jīng)元總數(shù)過小或過大時,神經(jīng)網(wǎng)絡(luò)訓(xùn)練過程中會產(chǎn)生欠擬合和過擬合現(xiàn)象[31].而隱藏層層數(shù)在增加到一定程度時,神經(jīng)網(wǎng)絡(luò)損失函數(shù)隨著訓(xùn)練次數(shù)的增加震蕩會加劇[32].隱藏層層數(shù)和神經(jīng)元數(shù)目對模型訓(xùn)練的結(jié)果影響如圖9 所示.

圖9 不同網(wǎng)絡(luò)參數(shù)對損失函數(shù)收斂特性的影響Fig.9 The influence of different network parameters to loss functions

從圖9 中可以看出,對于算例1,隱藏層層數(shù)為3,每層神經(jīng)元個數(shù)為60 左右可以達到最佳精度.相比之下,更多的隱藏層數(shù)目和神經(jīng)元數(shù)目會加重計算硬件的負擔(dān),且不能保證更好的精度.總的來說,兩者數(shù)量的選擇仍然具有很強的經(jīng)驗性.另外,從圖上可以看出,過大的隱藏層數(shù)目會加劇迭代過程中的震蕩,而過多的每層神經(jīng)元個數(shù)并沒有對計算精度帶來提升,更多神經(jīng)網(wǎng)絡(luò)參數(shù)對于計算精度的理論影響可參照文獻[33-36].

4.3 求解器

本次實驗所使用的深度學(xué)習(xí)平臺為Tensorflow,Tensorflow 中自帶了Adadelta,AdaGrad,Adam 和L-BFGS 等多種求解器,其基于的數(shù)學(xué)原理各不相同.以Adam 為基準,各優(yōu)化器在其他神經(jīng)網(wǎng)絡(luò)參數(shù)相同的情況下對神經(jīng)網(wǎng)絡(luò)結(jié)果的影響如圖10 和圖11所示.

圖11 使用各優(yōu)化器時對損失函數(shù)的影響Fig.11 The influence of optimizers to loss functions

從圖10 和圖11 可以看出,L-BFGS 優(yōu)化器由于是基于二階梯度的方法,在所有優(yōu)化器中有著最高的效率.在實際求解中,L-BFGS 優(yōu)化器經(jīng)常在達到最大訓(xùn)練次數(shù)之前就能使迭代步長小于誤差容限ε從而停止繼續(xù)訓(xùn)練.但根據(jù)重復(fù)實驗的結(jié)果,L-BFGS 求解器更容易陷入局部最優(yōu)解,從而不能得到真實的結(jié)果.一次典型的L-BFGS 陷入的局部最優(yōu)解如圖12 所示,可見其最大撓度過小.

圖12 L-BFGS 陷入的局部最優(yōu)解Fig.12 The local optimum result that L-BFGS stuck with

4.4 激活函數(shù)

在神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過程中,隱藏層的自適應(yīng)非線性激活函數(shù) σi可以基于多種基礎(chǔ)的非線性函數(shù)進行構(gòu)建.激活函數(shù)賦予神經(jīng)網(wǎng)絡(luò)的非線性性質(zhì)能夠讓神經(jīng)網(wǎng)絡(luò)擬合各種非線性函數(shù),常用的用于構(gòu)建自適應(yīng)激活函數(shù)的非線性函數(shù)有線性整流函數(shù)(relu),邏輯函數(shù)(sigmoid)和雙曲正切函數(shù)(tanh)3 種,它們對均質(zhì)薄板(SCSF)算例的影響如圖13所示.其余兩種函數(shù)的表現(xiàn)上,對于此學(xué)習(xí)率,雙曲正切函數(shù)的收斂速度更好,但精度不如邏輯函數(shù).值得注意的是,自適應(yīng)relu 函數(shù)雖然取得了最小的損失函數(shù)數(shù)值,但在實驗中作者觀察到,由于其在x=0附近較為嚴重的導(dǎo)數(shù)不連續(xù)性,在算例中容易造成神經(jīng)元死亡現(xiàn)象,因此不能正確地計算控制方程的損失函數(shù)項.在損失函數(shù)的變化曲線上,神經(jīng)元死亡表現(xiàn)為在整個過程中曲線震蕩較其他激活函數(shù)的曲線嚴重.一次典型的由于relu 函數(shù)特性導(dǎo)致其神經(jīng)網(wǎng)絡(luò)求出不符合物理規(guī)律的解的示意圖如圖14所示.

圖13 均質(zhì)薄板算例中不同激活函數(shù)對損失函數(shù)的影響Fig.13 The influence of the activation function on the loss function

圖14 relu 激活函數(shù)導(dǎo)致的求出不符合物理規(guī)律的解Fig.14 A result obtained by relu activation function that does not meet the governing equation

4.5 額外標簽數(shù)據(jù)對第一階屈曲載荷的影響

對于完整的板殼,由于其受到的殼邊載荷在達到其第一階屈曲載荷之前,其屈曲方程總是有平凡解,而非平凡解的變形模式必須在載荷大于臨界載荷時才能出現(xiàn).據(jù)此,可以將利用Rayleigh-Ritz 等能量方法假設(shè)出的變形形式作為標簽數(shù)據(jù)輸入到網(wǎng)絡(luò)里.假設(shè)3 種薄板算例的第一階屈曲模態(tài)分別為

系數(shù)A理論上在無缺陷板上可以取任意小的值,考慮到數(shù)值計算的特性,這里取A=0.1h,根據(jù)文獻[29]中的圖表,此時對應(yīng)的屈曲載荷大約為線性屈曲理論的 101%～105% 左右,本文取103% 作為對照值.將上述表達式在求解域上以均勻布點法選取40 個點,其標簽記為 (xT,yT),T∈[1,40],獲取的位移數(shù)據(jù)記為WT(xT,yT),并將這些位移數(shù)據(jù)作為位移得到標簽數(shù)據(jù)輸入到神經(jīng)網(wǎng)絡(luò).這樣PINN 就能額外在這些點上建立損失函數(shù),其損失函數(shù)項為

之后進行物理信息神經(jīng)網(wǎng)絡(luò)訓(xùn)練,值得注意的是,由于已經(jīng)提前讓網(wǎng)絡(luò)學(xué)習(xí)了屈曲的位移模式,網(wǎng)絡(luò)可以直接求解控制方程而無需進行外層迭代.訓(xùn)練中可以發(fā)現(xiàn),由于無需迭代和新增的標簽數(shù)據(jù)損失函數(shù),在擁有同樣網(wǎng)絡(luò)參數(shù)的物理信息神經(jīng)網(wǎng)絡(luò)下,擁有標簽數(shù)據(jù)的神經(jīng)網(wǎng)絡(luò)訓(xùn)練較原有的神經(jīng)網(wǎng)絡(luò)收斂更快.但由于提前跳出,最終損失函數(shù)較迭代法更大,且相對誤差更大.相應(yīng)地,訓(xùn)練周期可減少到原有的一半左右.訓(xùn)練后的結(jié)果具體可見圖15 與表3.

表3 擁有標簽數(shù)據(jù)訓(xùn)練得到的第一階屈曲載荷Table 3 The buckling load of the first buckling mode trained with labeled data

圖15 標簽數(shù)據(jù)對損失函數(shù)收斂特性的影響Fig.15 The influence of labeled data to loss functions

需要指出的是,在實際工程中,由于真實結(jié)構(gòu)的邊界條件十分復(fù)雜,要得到其屈曲模態(tài)幾乎是不可能的.因此在工程中,要獲得這樣的標簽數(shù)據(jù)難度較大,對于看重整個過程的屈曲問題更是如此.不可否認的是,標簽數(shù)據(jù)可以提升神經(jīng)網(wǎng)絡(luò)求解效率,減少所需的計算資源,因此從工程實際中獲得部分既容易測量又可供神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)的高質(zhì)量標簽數(shù)據(jù)將是作者今后工作的重點之一.

5 PINN 與有限元效率對比

在DEEPXDE 中,PINN 的求解需要經(jīng)歷3 個步驟,分別為建立微分方程及數(shù)據(jù)(build),編譯損失函數(shù)(compile)和訓(xùn)練神經(jīng)網(wǎng)絡(luò)(train).本文部分算例求解時迭代所需的平均時間及時間的組成如圖16所示.從圖中可以看出,PINN 方法耗時最長的部分是訓(xùn)練神經(jīng)網(wǎng)絡(luò)部分,不論是迭代到誤差容限以下結(jié)束訓(xùn)練,還是迭代至最大訓(xùn)練數(shù),訓(xùn)練損失函數(shù)部分通常需要占據(jù)神經(jīng)網(wǎng)絡(luò)訓(xùn)練總時長的80%～90%左右,并且由于非線性屈曲問題的迭代特性,使得時間復(fù)雜度進一步提升.

圖16 本文算例所需時間及其組成Fig.16 The computation time of the numerical examples and its composition

本文的有限元結(jié)果由商業(yè)有限元軟件Abaqus 2021 得到,其中線性屈曲載荷由求解器“Linear perturbation,Buckle”求得,而非線性屈曲模態(tài)和最大載荷因子則由求解器“Static,Riks”求得.從圖16 可以看出PINN 方法相對成熟度高的商業(yè)有限元軟件計算效率偏低,其原因主要有以下幾點:(1) PINN 求解的是強形式的控制方程,其損失函數(shù)雖然能夠表征物理規(guī)律,但對應(yīng)的優(yōu)化問題有可能高度非凸,總體需要的計算量大.(2) 商業(yè)有限元軟件具有成熟高效的求解流程,在每一步的計算效率較高.(3) PINN在求解時需要更多的單輪迭代步驟.即使有了額外的數(shù)據(jù)標簽,其求解效率也僅有部分提高.有限元軟件是直接求解弱形式的物理方程,單輪迭代一次僅需要求解一次方程.因此計算效率更高.

6 結(jié)論與展望

本文基于深度學(xué)習(xí)技術(shù)、屈曲控制方程的強形式求解與迭代法建立了一種求解薄壁板殼結(jié)構(gòu)屈曲問題的神經(jīng)網(wǎng)絡(luò)方法.通過兩種典型結(jié)構(gòu)的算例,可以得出以下幾個結(jié)論.

(1) 本方法的結(jié)果與有限元結(jié)果的誤差大多時候能控制在10%以內(nèi),且在數(shù)據(jù)驅(qū)動方面不需要額外的撓度和力學(xué)量數(shù)據(jù),這一特性展現(xiàn)了此方法在工程中的實用性和合理性.

(2) 本方法的結(jié)果是對強形式偏微分控制方程進行求解的結(jié)果,因此數(shù)值解具備相對于控制方程的連續(xù)性,滿足物理定律,這一點有別于一般數(shù)值驅(qū)動神經(jīng)網(wǎng)絡(luò)所帶來的結(jié)果,但由于屈曲問題的數(shù)據(jù)集收集較為困難,本文未能將PINN 結(jié)果和一般數(shù)據(jù)驅(qū)動的神經(jīng)網(wǎng)絡(luò)結(jié)果作對比.

(3) 由于屈曲臨界載荷和物理信息神經(jīng)網(wǎng)絡(luò)本身都需要進行迭代,因此本文方法在計算部分的時間復(fù)雜度上劣于傳統(tǒng)的有限元方法,但其不需要網(wǎng)格前處理的特性有望減少解決復(fù)雜問題時的整體工作量.

根據(jù)目前所得出的結(jié)論,進一步工作可以從以下幾個方面展開.

(1) 本文所建立的神經(jīng)網(wǎng)絡(luò)模型采用了硬邊界約束(hard constraint),其本質(zhì)是構(gòu)造滿足邊界條件的試函數(shù)并在輸入層進行變換使神經(jīng)網(wǎng)絡(luò)的擬合解能夠滿足邊界條件,這一部分仍然具有較強的經(jīng)驗性,有待實現(xiàn)自動化,有興趣的讀者可參考文獻[36].

(2) 由4.5 節(jié)可見,標簽數(shù)據(jù)對增加PINN 訓(xùn)練效率、減少計算資源使用有著較為顯著的作用,但本文算例中給出的標簽數(shù)據(jù)可以較為簡單地獲取,而實際工程中則不然.因此,如何從工程實際結(jié)構(gòu)中獲得容易測量又可作為額外標簽的物理量提供給PINN 學(xué)習(xí)將會是PINN 算法在實際工程運用中一個十分有前景的問題.

(3) 本文所建立的神經(jīng)網(wǎng)絡(luò)模型仍然是建立在人為經(jīng)驗的基礎(chǔ)上,因此具有較大的優(yōu)化空間,其中一個可改進的方面是優(yōu)化隱藏層的布置,例如卷積層等其他神經(jīng)函數(shù)隱藏層用于神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)屈曲模式.此有望用于在復(fù)雜結(jié)構(gòu)的屈曲過程中判斷屈曲模式,為實際工程優(yōu)化問題提供進一步解決方案.

附錄A:本文主要符號變量表