馮唐思捷 梁 偉

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

引言

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

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

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

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

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

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

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

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

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

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

以及物理方程

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

邊界條件可分為3 類(lèi):簡(jiǎn)支、固支與自由,其對(duì)應(yīng)的表達(dá)式分別為

簡(jiǎn)支

固支

自由邊

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

這樣非線性映射則可自動(dòng)滿足力學(xué)模型的位移邊界條件.對(duì)于邊界條件(8),需要融入到損失函數(shù)的邊界條件減少了一項(xiàng),其余的邊界條件也可以構(gòu)建類(lèi)似的試函數(shù)來(lái)減少損失函數(shù)項(xiàng)數(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ù)為L(zhǎng)a×Nu,前者為隱藏層層數(shù),后者為單個(gè)隱藏層神經(jīng)元數(shù).神經(jīng)網(wǎng)絡(luò)的輸入為幾何域坐標(biāo) (x,y),經(jīng)過(guò)神經(jīng)網(wǎng)絡(luò)輸出的擬合解為撓度與應(yīng)力函數(shù)該問(wèn)題所對(duì)應(yīng)的非線性映射可以在數(shù)學(xué)上表示為

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

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

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

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

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

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

MSEI作為表示初值的誤差,在屈曲方程中由于不含時(shí),因此初值項(xiàng)為0.

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

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

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

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

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

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

使用弧長(zhǎng)法求解PINN 得出的微分方程擬合解的增量并進(jìn)行下一步求解的方法如算法2 所示.

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

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

3 數(shù)值算例

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

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

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

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

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

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

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

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

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

圖3 帶有缺陷的均質(zhì)薄板的非線性屈曲位移云圖對(duì)比 (上: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ù)測(cè).這證明采用弧長(zhǎng)法進(jìn)行屈曲控制并通過(guò)PINN 求解器計(jì)算屈曲問(wèn)題是可行的.

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

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

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

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

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

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

其中非線性項(xiàng)為

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

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

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

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

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

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

圖5 帶有缺陷的均質(zhì)圓筒殼在最大 L PFmax 處非線性屈曲位移云圖對(duì)比 (續(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ìn)行神經(jīng)網(wǎng)絡(luò)的訓(xùn)練時(shí),對(duì)最大承載時(shí)的位移模式求解較好,其位移分布和最大位移大小均和有限元結(jié)果接近.由PINN方法求出的軸壓作用下后屈曲載荷-端部縮短曲線如圖6 所示.

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

3.3 圓柱筒殼的缺陷敏感度

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

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

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

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

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

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

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

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

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

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

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

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

4.3 求解器

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

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

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

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

圖13 均質(zhì)薄板算例中不同激活函數(shù)對(duì)損失函數(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 額外標(biāo)簽數(shù)據(jù)對(duì)第一階屈曲載荷的影響

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

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

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

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

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

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

5 PINN 與有限元效率對(duì)比

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

圖16 本文算例所需時(shí)間及其組成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 方法相對(duì)成熟度高的商業(yè)有限元軟件計(jì)算效率偏低,其原因主要有以下幾點(diǎn):(1) PINN 求解的是強(qiáng)形式的控制方程,其損失函數(shù)雖然能夠表征物理規(guī)律,但對(duì)應(yīng)的優(yōu)化問(wèn)題有可能高度非凸,總體需要的計(jì)算量大.(2) 商業(yè)有限元軟件具有成熟高效的求解流程,在每一步的計(jì)算效率較高.(3) PINN在求解時(shí)需要更多的單輪迭代步驟.即使有了額外的數(shù)據(jù)標(biāo)簽,其求解效率也僅有部分提高.有限元軟件是直接求解弱形式的物理方程,單輪迭代一次僅需要求解一次方程.因此計(jì)算效率更高.

6 結(jié)論與展望

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

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

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

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

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

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

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

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

附錄A:本文主要符號(hào)變量表