黃鐘民 謝 臻 張易申 彭林欣 ,?,2)

* (廣西大學(xué)土木建筑工程學(xué)院,南寧 530004)

? (廣西大學(xué)廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室,工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室,南寧 530004)

引言

面內(nèi)功能梯度材料薄板結(jié)構(gòu)在土木工程、海洋工程等領(lǐng)域應(yīng)用廣泛,該結(jié)構(gòu)由功能梯度材料所組成(functionally graded material,FGM),其材料特性隨著空間位置的變化而表現(xiàn)為梯度性的變化[1].在現(xiàn)有研究中,于天崇等[2]研究了面內(nèi)變剛度薄板在特定邊界下彎曲問(wèn)題的Levy 解.朱竑禎等[3]研究了周邊固支圓形面內(nèi)變剛度薄板軸對(duì)稱彎曲問(wèn)題的級(jí)數(shù)解答.何建璋等[4]研究了面內(nèi)變剛度矩形薄板自由振動(dòng)問(wèn)題的辛彈性解.以上的理論解答僅針對(duì)特定的功能梯度函數(shù)及特定邊界才成立,一般的情況下難以得出理論解答,而在數(shù)值解法上仍以有限元為主.Santare 和Lambros[5]發(fā)展了一種針對(duì)材料屬性為指數(shù)分布的梯度有限元求解格式.Kim 和Paulino[6]研究了梯度單元以及分層單元在不同荷載下的計(jì)算性能.黃立新等[7]基于分層法思想分析了功能梯度材料的平面應(yīng)力問(wèn)題.田云德和秦世倫[8]采用分層法研究了功能梯度厚板的熱應(yīng)力問(wèn)題.對(duì)于面內(nèi)變剛度功能梯度薄板,采用分層法,薄板求解域采用有限元網(wǎng)格劃分,每個(gè)單元的材料參數(shù)為常數(shù),而其材料參數(shù)則根據(jù)功能梯度函數(shù)由單元內(nèi)特定點(diǎn)進(jìn)行計(jì)算.有限元網(wǎng)格劃分越密其計(jì)算結(jié)果越精確,而在實(shí)際計(jì)算中,越精細(xì)的網(wǎng)格會(huì)導(dǎo)致總體剛度矩陣規(guī)模巨大,需要耗費(fèi)大量的計(jì)算機(jī)內(nèi)存.無(wú)論采用何種數(shù)值方法,其最終目的均是求得面內(nèi)變剛度薄板彎曲控制偏微分方程的近似解答,為進(jìn)一步豐富該類研究,本文擬結(jié)合深度學(xué)習(xí)技術(shù)并發(fā)展求解該類問(wèn)題的新解法.

在早期就有研究[9-10]將人工神經(jīng)網(wǎng)絡(luò)作為一類偏微分方程的求解器用于求解偏微分方程,但由于其對(duì)計(jì)算機(jī)計(jì)算能力的要求過(guò)高以及優(yōu)化算法中存在的問(wèn)題,這一解法在當(dāng)時(shí)并未得到很好的發(fā)展.而如今,自深度學(xué)習(xí)在計(jì)算機(jī)視覺(jué)、語(yǔ)音文字識(shí)別取得成功的應(yīng)用后,深度學(xué)習(xí)技術(shù)也在各個(gè)學(xué)科領(lǐng)域加速發(fā)展.在力學(xué)領(lǐng)域,Weinan 和Yu[11]提出深度Ritz 法,該方法采用變分求解形式對(duì)偏微分方程進(jìn)行求解.Raissi 等[12]提出了用于求解高階非線性偏微分方程的物理驅(qū)動(dòng)的神經(jīng)網(wǎng)絡(luò)(physics-informed neural networks,PINNs).Sirignano 和Spiliopoulos[13]則提出求解高階微分方程的深度伽遼金法(deep galerkin method,DGM).Samaniego 等[14]建立了深度能量法并將其應(yīng)用于求解彈性、超彈性等力學(xué)問(wèn)題.瞿同明等[15]基于深度學(xué)習(xí)技術(shù),研究了細(xì)觀力學(xué)中的顆粒本構(gòu)關(guān)系.謝晨月等[16]發(fā)展了一種模擬湍流大渦的神經(jīng)網(wǎng)絡(luò)方法.劉宇翔等[17]基于卷積神經(jīng)網(wǎng)絡(luò)研究了無(wú)網(wǎng)格方法中影響域的優(yōu)化問(wèn)題.郭宏偉和莊曉瑩[18]采用深度配點(diǎn)法以及深度能量法求解了薄板彎曲問(wèn)題.陳豪龍和柳占立[19]基于數(shù)據(jù)驅(qū)動(dòng)的神經(jīng)網(wǎng)絡(luò)模型求解了熱傳導(dǎo)反問(wèn)題.

在上述研究中,神經(jīng)網(wǎng)絡(luò)解法[11-14]并不像有限元解法一樣可以輕松施加邊界條件,早期的研究采取根據(jù)邊界條件構(gòu)造滿足偏微分方程特解試函數(shù)的形式來(lái)處理邊界條件,但采用該方法會(huì)使得簡(jiǎn)支邊、自由邊試函數(shù)的表達(dá)式變得復(fù)雜,導(dǎo)致程序的實(shí)現(xiàn)較為復(fù)雜.近期的研究則采用罰函數(shù)的方法將邊界處的誤差納入神經(jīng)網(wǎng)絡(luò)的訓(xùn)練誤差中,從而將原問(wèn)題轉(zhuǎn)換為無(wú)約束優(yōu)化問(wèn)題,在實(shí)際計(jì)算中,也會(huì)存在著由于邊界誤差項(xiàng)難以收斂而影響求解精度的情況[20].

同時(shí)由于彎曲剛度函數(shù)是面內(nèi)坐標(biāo)的連續(xù)函數(shù),面內(nèi)變剛度薄板彎曲問(wèn)題的控制方程為一包含了彎曲剛度導(dǎo)數(shù)項(xiàng)的復(fù)雜4 階偏微分方程,在實(shí)際計(jì)算中采用DGM 和PINN 等方法對(duì)其求解時(shí),會(huì)存在由于彎曲剛度偏導(dǎo)數(shù)在域內(nèi)不收斂而導(dǎo)致網(wǎng)絡(luò)擬合不佳的問(wèn)題.

基于上述原因,本文針對(duì)薄板彎曲問(wèn)題求解的特點(diǎn),結(jié)合前面所述的兩種邊界處理方案,建立了一種針對(duì)面內(nèi)變剛度薄板彎曲問(wèn)題的非全連接前饋神經(jīng)網(wǎng)絡(luò)模型,該模型包含撓度網(wǎng)絡(luò)與彎矩網(wǎng)絡(luò):撓度網(wǎng)絡(luò)用于預(yù)測(cè)薄板的撓度,彎矩網(wǎng)絡(luò)用于預(yù)測(cè)薄板的彎矩,進(jìn)而將問(wèn)題轉(zhuǎn)換為求解4 個(gè)二階偏微分方程組.在邊界條件的處理上,本文仍采用罰函數(shù)方法,不同之處在于本文模型的輸出為撓度、彎矩,因而可根據(jù)位移邊界條件對(duì)撓度網(wǎng)絡(luò)構(gòu)造試函數(shù),根據(jù)廣義應(yīng)力邊界條件對(duì)彎矩網(wǎng)絡(luò)構(gòu)造試函數(shù),這使得本文模型對(duì)于常見(jiàn)的邊界條件的施加更為簡(jiǎn)便,進(jìn)而減小邊界誤差項(xiàng)帶來(lái)的影響,同時(shí)計(jì)算效率也得到提高.本文采用Pytorch 深度學(xué)習(xí)框架編寫求解程序,選取不同邊界條件的面內(nèi)變剛度薄板算例,在Ubuntu Kylin 操作系統(tǒng)上進(jìn)行計(jì)算,計(jì)算機(jī)的CPU配置為Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz,8GB 內(nèi)存,并將計(jì)算所得結(jié)果與理論解、有限元解進(jìn)行對(duì)比分析,以驗(yàn)證本文方法的有效性.

1 面內(nèi)變剛度薄板彎曲理論

本文研究變厚度薄板或彈性模量參數(shù)在面內(nèi)變化的薄板的彈性小變形彎曲,設(shè)薄板的厚度函數(shù)為h(x,y),材料的泊松比ν為常數(shù),彈性模量函數(shù)為E(x,y).

根據(jù)Kirchhoff 板理論基本假定,幾何方程為

物理方程為

其中n表示邊界的外法線方向,s表示邊界的切線方向.

將式(4)～ 式(6)代入平衡方程(3)即可得面內(nèi)變剛度薄板彎曲偏微分控制方程

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

2.1 神經(jīng)網(wǎng)絡(luò)模型的構(gòu)建

本文方法并非直接設(shè)計(jì)網(wǎng)絡(luò)來(lái)求解方程(12),而是采用兩個(gè)神經(jīng)網(wǎng)絡(luò)模型來(lái)進(jìn)行求解,如圖1 所示,將待求解的4 階偏微分控制方程轉(zhuǎn)換為求解4 個(gè)二階偏微分方程組,該解法本質(zhì)上仍屬于強(qiáng)形式的求解方案.如果僅以撓度作為預(yù)測(cè)解,在試函數(shù)的構(gòu)造上對(duì)于不同形狀的求解域以及簡(jiǎn)支、自由邊界條件的構(gòu)造會(huì)出現(xiàn)困難.本文采用撓度網(wǎng)絡(luò)預(yù)測(cè)薄板撓度{彎矩網(wǎng)絡(luò)預(yù)測(cè)薄板彎}矩,這樣的做法可以使得位移邊界條件由撓度網(wǎng)絡(luò)施加,廣義應(yīng)力邊界條件由彎矩網(wǎng)絡(luò)施加.

圖1 本文神經(jīng)網(wǎng)絡(luò)模型示意圖Fig.1 The schematic diagram of neural network model in this paper

2.2 誤差函數(shù)

根據(jù)廣義應(yīng)力?應(yīng)變關(guān)系式(6)～ 式(8)可求得撓度二階偏導(dǎo)的彎矩表達(dá)式

誤差函數(shù)的構(gòu)造是神經(jīng)網(wǎng)絡(luò)訓(xùn)練的核心,由于本文方法引入兩個(gè)網(wǎng)絡(luò)進(jìn)行計(jì)算,故在訓(xùn)練中需要考慮兩者之間的耦合誤差.若采用無(wú)約束優(yōu)化方案,本文誤差函數(shù)主要根據(jù)撓度與彎矩網(wǎng)絡(luò)在邊界處的誤差、彎矩網(wǎng)絡(luò)在力平衡方程(3)中的誤差、預(yù)測(cè)的撓度與彎矩通過(guò)式(13)～ 式(15)建立的耦合誤差來(lái)構(gòu)造.

采用均方誤差(mean square error,MSE)來(lái)衡量神經(jīng)網(wǎng)絡(luò)的擬合誤差,記撓度網(wǎng)絡(luò)與彎矩網(wǎng)絡(luò)的內(nèi)部參數(shù)分別為 θw,θm,本文模型的誤差函數(shù)可構(gòu)造為

如果采用構(gòu)造試函數(shù)的形式使得邊界誤差強(qiáng)制滿足,則誤差函數(shù)無(wú)需計(jì)算邊界誤差

其中,向量x=(x,y)表示神經(jīng)網(wǎng)絡(luò)的輸入;n和s分別為邊界的法線、切線方向;Ω 表示求解域;Γ1,Γ2,Γ3分別為固支、簡(jiǎn)支、自由邊界,? Ω=Γ1∪Γ2∪Γ3;為施加于邊界處的彎矩,Mn(x;θm),Qn(x;θm)分別為根據(jù)彎矩網(wǎng)絡(luò)輸出求得的彎矩、剪力;kp為網(wǎng)絡(luò)耦合系數(shù),取值范圍為1～ 1000,該系數(shù)的選取會(huì)影響撓度網(wǎng)絡(luò)與彎矩網(wǎng)絡(luò)之間的耦合效果;k1,k2,k3為邊界處的罰系數(shù),取值范圍為1～ 10 000.

式(22)中

本文的誤差函數(shù)表達(dá)式中包含撓度、彎矩對(duì)自變量的二階偏導(dǎo)數(shù)項(xiàng),對(duì)于這些偏導(dǎo)項(xiàng),一方面可以利用神經(jīng)網(wǎng)絡(luò)的輸出構(gòu)造差分求解格式來(lái)近似求解,但采用該方案需要較大的計(jì)算量才能得到精確的計(jì)算結(jié)果;另一方面,基于計(jì)算圖的自動(dòng)微分技術(shù)(automatic differentiation,AD)可以高效地處理神經(jīng)網(wǎng)絡(luò)對(duì)輸入變量求導(dǎo)過(guò)程,當(dāng)前的深度學(xué)習(xí)框架如Tensoflow,Pytorch,MindsSpore 等均支持自動(dòng)微分.本文基于Pytorch 提供的自動(dòng)微分接口實(shí)現(xiàn)對(duì)上述偏導(dǎo)項(xiàng)及誤差函數(shù)梯度的計(jì)算.

在實(shí)際計(jì)算中,也可靈活采用混合邊界誤差的形式進(jìn)行求解,如對(duì)于部分簡(jiǎn)單的邊界條件構(gòu)造特解,而對(duì)于復(fù)雜的邊界采用相應(yīng)的無(wú)約束優(yōu)化方案.建立本文的誤差函數(shù)后,在每個(gè)訓(xùn)練批次(epoch)中均需計(jì)算其梯度并結(jié)合誤差反向傳播算法更新網(wǎng)絡(luò)的內(nèi)部參數(shù),關(guān)于該過(guò)程,Tang 和Yang[21]對(duì)其進(jìn)行了詳細(xì)的討論.

2.3 學(xué)習(xí)率選取方案

學(xué)習(xí)率的選取可直接影響神經(jīng)網(wǎng)絡(luò)的訓(xùn)練,目前神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)率的選取仍帶有一定的經(jīng)驗(yàn)性,但總體而言,在訓(xùn)練初期選取較大的學(xué)習(xí)率可以加快誤差收斂速度,在訓(xùn)練后期,此時(shí)神經(jīng)網(wǎng)絡(luò)模型已經(jīng)學(xué)習(xí)到相應(yīng)的特征,此時(shí)往往需要降低學(xué)習(xí)率,以便對(duì)神經(jīng)網(wǎng)絡(luò)內(nèi)部參數(shù)進(jìn)行微調(diào),使得誤差波動(dòng)幅度不至于過(guò)大.經(jīng)過(guò)本文的實(shí)踐,本文學(xué)習(xí)率選取方案如下

其 中,t為訓(xùn)練次數(shù).

2.4 算法流程

算法1.本文算法

3 算例

3.1 雙向面內(nèi)變剛度圓形薄板軸對(duì)稱彎曲

圖2 所示受橫向均布荷載q(x,y)=?q0作用的周邊固支圓形薄板,其半徑R,ν=0.3,彎曲剛度函數(shù)沿半徑變化D(ρ)=D0e?mρ/a,其中m 為梯度系數(shù),q0,D0為常數(shù).該問(wèn)題存在理論解[14].

圖2 圓形面內(nèi)變剛度薄板Fig.2 Circular thin plate with in-plane stiffness gradient

本文選取梯度參數(shù)m分別為0,0.5,1,2 的情況進(jìn)行計(jì)算.本算例在計(jì)算過(guò)程中僅需要施加位移邊界條件,設(shè)w(x,y;θw)為撓度網(wǎng)絡(luò)的輸出,考慮到本算例邊界條件較為簡(jiǎn)單,構(gòu)造撓度試函數(shù)為w?=w(x,y;θw)

本算例的撓度模型以及彎矩模型均采用具有6 層隱藏層,每層隱藏層具備30 個(gè)神經(jīng)元的網(wǎng)絡(luò)結(jié)構(gòu)(記為6 × 30),采用x2作為激活函數(shù),kp=100,NΩ=640 ;采用Adam 優(yōu)化算法,對(duì)各個(gè)工況所采用的學(xué)習(xí)率方案均一致.本文算例在不同的梯度下上述基本參數(shù)不變,僅更改梯度系數(shù).為更加詳細(xì)地顯示誤差函數(shù)的變化情況,如無(wú)特殊說(shuō)明,本文均取誤差的十進(jìn)制對(duì)數(shù)作為等效誤差并繪制相應(yīng)的訓(xùn)練誤差曲線.

每個(gè)梯度參數(shù)下神經(jīng)網(wǎng)絡(luò)模型的計(jì)算誤差如圖3 所示,由此可見(jiàn),m=0 時(shí)為剛度恒定的薄板,此時(shí)誤差函數(shù)收斂較快,相比其他工況其誤差最終的收斂值最小;隨著梯度系數(shù)的增大,訓(xùn)練誤差最終的收斂值出現(xiàn)增大的趨勢(shì).

表1 本文方法計(jì)算與理論解對(duì)比(無(wú)量綱)Table 1 Comparison of dimensionless calculated by neural network method and the theoretical solution

表1 本文方法計(jì)算與理論解對(duì)比(無(wú)量綱)Table 1 Comparison of dimensionless calculated by neural network method and the theoretical solution

注:本文誤差計(jì)算公式為 u 為本文方法計(jì)算結(jié)果,u?為理論解 (Note:The relative error in this paper is calculated by 100%,where u is the calculation results of this paper,u ? is the theoretical solution)

表2 本文方法計(jì)算與理論解對(duì)比(無(wú)量綱)Table 2 Comparison of dimensionless calculated by neural network method and the theoretical solution

表2 本文方法計(jì)算與理論解對(duì)比(無(wú)量綱)Table 2 Comparison of dimensionless calculated by neural network method and the theoretical solution

圖3 神經(jīng)網(wǎng)絡(luò)訓(xùn)練誤差曲線Fig.3 The convergence curve of neural network training error

圖4 PINN 求解圓形面內(nèi)變剛度薄板彎曲問(wèn)題的訓(xùn)練誤差收斂曲線(m=2)Fig.4 Training error convergence curve of PINN (m=2)

為了說(shuō)明本文方法在求解面內(nèi)變剛度薄板彎曲問(wèn)題上的優(yōu)點(diǎn),本算例也利用PINN 來(lái)求解其四階偏微分控制方程(12),采用隱藏層層數(shù)為6,每層隱藏層具有30 個(gè)神經(jīng)元的神經(jīng)網(wǎng)絡(luò)模型對(duì)工況m=2進(jìn)行求解,激活函數(shù)為Tanh 函數(shù),訓(xùn)練的數(shù)據(jù)點(diǎn)由求解域中隨機(jī)生成,數(shù)據(jù)點(diǎn)的產(chǎn)生有兩種方案:

方案(1)為在整個(gè)求解域中隨機(jī)生成訓(xùn)練數(shù)據(jù)點(diǎn),此時(shí)的數(shù)據(jù)點(diǎn)可在原點(diǎn)附近生成;

方案(2)為在求解域中隨機(jī)生成的數(shù)據(jù)點(diǎn)但離原點(diǎn)較遠(yuǎn).此時(shí)兩方案在訓(xùn)練過(guò)程中的誤差收斂情況如圖4 所示,可見(jiàn)采用相同的模型,而生成的數(shù)據(jù)點(diǎn)不同則會(huì)導(dǎo)致模型的訓(xùn)練出現(xiàn)不同的結(jié)果,雖然采用數(shù)據(jù)點(diǎn)生成方案(1)的模型訓(xùn)練也收斂,但由于其誤差此時(shí)收斂于一個(gè)較大的值,得不到正確解.經(jīng)過(guò)本文分析,這主要是由于本算例的彎曲剛度函數(shù)D的二階導(dǎo)數(shù)在靠近原點(diǎn)區(qū)域出現(xiàn)“爆炸”式變化的原因,即剛度函數(shù)的二階偏導(dǎo)在原點(diǎn)處不收斂,在靠近原點(diǎn)處等的解答急劇增大,這會(huì)導(dǎo)致PINN 采用方案(1)訓(xùn)練時(shí),遇到靠近原點(diǎn)處的點(diǎn),計(jì)算所得域內(nèi)誤差突然增大,進(jìn)而導(dǎo)致誤差訓(xùn)練難以收斂.PINN 最初提出時(shí)并未考慮求解域內(nèi)存在奇異點(diǎn)的情況,對(duì)于該情況,一般情況下可在生成的數(shù)據(jù)點(diǎn)中排除掉奇異點(diǎn),但對(duì)于本算例中奇異點(diǎn)處被施予荷載的情況,如果不能很好地處理則會(huì)影響求解的精度.對(duì)此,本文認(rèn)為可以弱化相應(yīng)的偏微分控制方程再利用PINN 求解,也可參考本文思路,結(jié)合神經(jīng)網(wǎng)絡(luò)解法的特點(diǎn),根據(jù)具體問(wèn)題對(duì)原偏微分控制方程等效化處理.本文方法在求解時(shí)并非直接從方程(12)入手,而是通過(guò)求解一系列偏微分方程組來(lái)逼近真實(shí)解,避開(kāi)了對(duì)彎曲剛度函數(shù)求偏導(dǎo)數(shù),故其求解僅與域內(nèi)的剛度值有關(guān),其適應(yīng)性更強(qiáng),對(duì)薄板彎曲問(wèn)題的求解更具“魯棒性”.

3.2 單向面內(nèi)變剛度方形薄板受非線性荷載作用

如圖5 所示邊長(zhǎng)為a的方形薄板,厚度h,ν=0.3,1,2,3 邊固支,4 邊簡(jiǎn)支,其彎曲剛度函數(shù)為為常數(shù),m為梯度系數(shù),受橫向非線性荷載作用,利用本文方法求解m分別為0,1,2 情況下的撓度、內(nèi)力.

圖5 方形面內(nèi)變剛度薄板Fig.5 Thin square plates with in-plane stiffness gradient

本算例選取6 × 30 的撓度網(wǎng)絡(luò)模型,5 × 50 的彎矩網(wǎng)絡(luò)模型,激活函數(shù)選擇Tanh 函數(shù),kp=100,NΩ=450,采用Adam 優(yōu)化算法.訓(xùn)練誤差曲線如圖6 所示.設(shè)w(x,y;θw)為撓度網(wǎng)絡(luò)的輸出,彎矩網(wǎng)絡(luò)的輸出為M(y},考慮本算例?的位移及廣義應(yīng)力邊界條件,撓度及彎矩試函數(shù)構(gòu)造為

圖6 神經(jīng)網(wǎng)絡(luò)訓(xùn)練誤差曲線Fig.6 Neural network training error-curve

如果將本算例的應(yīng)力邊界條件也通過(guò)撓度網(wǎng)絡(luò)進(jìn)行構(gòu)造,則撓度表達(dá)式將復(fù)雜許多.

將本算例的計(jì)算結(jié)果與有限元解答對(duì)比,有限元計(jì)算中每個(gè)單元的彎曲剛度根據(jù)單元的形心坐標(biāo)計(jì)算,采用50× 50 的矩形薄板非協(xié)調(diào)單元來(lái)對(duì)求解域進(jìn)行離散,離散方案通過(guò)小片測(cè)試,對(duì)本文的3 種工況該網(wǎng)格離散方案均收斂.將本文與有限元計(jì)算結(jié)果的無(wú)量綱撓度、彎矩進(jìn)行對(duì)比分析,無(wú)量綱計(jì)算公式為w?和M?為實(shí)際計(jì)算所得的撓度、彎矩值.

由撓度計(jì)算圖7 可知,在撓度的求解上,本文解法與有限元解法一致.由圖8 的彎矩對(duì)比圖可發(fā)現(xiàn),本文彎矩解與有限元解答基本吻合,而當(dāng)梯度系數(shù)m=2 時(shí),雖然本文解與有限元解在部分點(diǎn)上彎矩的相對(duì)誤差增大,但整體上解答吻合.

圖7 不同梯度參數(shù)下本文撓度計(jì)算結(jié)果與有限元對(duì)比(y=0.5 m)Fig.7 Comparison of dimensionless deflection calculation results in this paper with FEM when different gradient parameters (y=0.5 m)

圖8 不同梯度參數(shù)下本文彎矩 計(jì)算結(jié)果與有限元對(duì)比(y=0.5 m)Fig.8 Comparison of dimensionless bending moment calculation results of this paper with FEM when different gradient parameters(y=0.5 m)

3.3 三角形面內(nèi)變剛度薄板受線性分布荷載作用

圖9 三角形面內(nèi)變剛度薄板Fig.9 Thin triangular plate with in-plane stiffness gradient

圖10 三角形面內(nèi)變剛度薄板沿軸線x=0 上的撓度 分布Fig.10 Dimensionless deflection variation of thin triangular plate with in-plane stiffness gradient along axis x=0

圖11 m=0.2 時(shí)三角形面內(nèi)變剛度薄板彎矩的有限元計(jì)算結(jié)果Fig.11 Finite element calculation of bending moment of thin plate with triangular in-plane variable stiffness (m=0.2)

圖12 m=0.2 時(shí)三角形面內(nèi)變剛度薄板彎矩的本文計(jì)算結(jié)果Fig.12 Neural network method calculation of bending moment of thin plate with triangular in-plane variable stiffness (m=0.2)

4 本文方法計(jì)算效率及其有限元對(duì)比分析

以T表示算例各個(gè)工況下的平均用時(shí),T′表示各個(gè)工況下誤差函數(shù)的構(gòu)建及其梯度計(jì)算的平均耗時(shí),本文各算例求解時(shí)迭代所需的數(shù)據(jù)點(diǎn)數(shù)及所需的平均時(shí)間、內(nèi)存如表3 所示,可看出本文誤差函數(shù)的構(gòu)建及其梯度的求解占據(jù)了神經(jīng)網(wǎng)絡(luò)訓(xùn)練總時(shí)長(zhǎng)的70%左右.

表3 本文各算例求解所需的數(shù)據(jù)點(diǎn)數(shù)、內(nèi)存、時(shí)間Table 3 The number of training data points,computational memory and computing time of numerical examples in this paper

本文的有限元求解程序采用python 語(yǔ)言進(jìn)行編寫,選用每個(gè)節(jié)點(diǎn)有3 個(gè)自由度的薄板彎曲單元計(jì)算,剛度矩陣以 compressed sparse column (CSC)格式的稀疏矩陣存儲(chǔ),利用科學(xué)計(jì)算庫(kù)scipy 中的線性求解器求解剛度方程.根據(jù)有限元解答的最小撓度判斷有限元解答是否收斂,算例2 和算例3 中各個(gè)工況下有限元計(jì)算收斂時(shí)的節(jié)點(diǎn)數(shù)目如表4 所示,本文有限元計(jì)算所需的節(jié)點(diǎn)數(shù)與計(jì)算所需內(nèi)存的關(guān)系如圖13 所示.

圖13 有限元計(jì)算所需內(nèi)存與節(jié)點(diǎn)數(shù)的關(guān)系(薄板單元)Fig.13 The relationship between the number of nodes and the memory needed in finite element calculation using thin plate bending element

表4 算例2、算例3 的有限元求解收斂所需節(jié)點(diǎn)數(shù)目Table 4 The number of nodes needed for the convergence of the finite element solution of numerical example 2 and 3

可以發(fā)現(xiàn)有限元解答收斂時(shí)所需的節(jié)點(diǎn)(網(wǎng)格)數(shù)目隨著梯度系數(shù)的增大而增大,節(jié)點(diǎn)數(shù)與所需內(nèi)存并非呈現(xiàn)線性關(guān)系,節(jié)點(diǎn)數(shù)的增多會(huì)導(dǎo)致所需計(jì)算內(nèi)存的急劇增大,而神經(jīng)網(wǎng)絡(luò)方法在梯度系數(shù)增大時(shí),在單次迭代時(shí)仍可以較少的數(shù)據(jù)點(diǎn)進(jìn)行迭代,這使得本文方法在求解時(shí)所需內(nèi)存較小.在時(shí)間上,以算例3 中梯度系數(shù)m=0.5 為例,此時(shí)有限元求解收斂時(shí),所需內(nèi)存為291.5 MiB,剛度方程從組裝到求解用時(shí)1.798 s,可以看出本文解法的求解速度明顯慢于有限元,一方面是由于本文方法求解的是強(qiáng)形式的偏微分控制方程,與求解弱形式的方程相比,往往需要更多的迭代次數(shù),另一方面,神經(jīng)網(wǎng)絡(luò)方法是一類以數(shù)據(jù)為驅(qū)動(dòng)的解法,在求解過(guò)程中需要往復(fù)迭代,而對(duì)于本文研究的線性問(wèn)題,有限元僅需求解一次剛度方程.

5 網(wǎng)絡(luò)耦合系數(shù) kp 、隱藏層層數(shù)與神經(jīng)元個(gè)數(shù)、激活函數(shù)對(duì)網(wǎng)絡(luò)收斂性的影響分析

5.1 網(wǎng)絡(luò)耦合系數(shù) kp

本文模型由于引入兩個(gè)神經(jīng)網(wǎng)絡(luò)模型對(duì)面內(nèi)變剛度薄板彎曲問(wèn)題進(jìn)行計(jì)算,與采用單網(wǎng)絡(luò)的神經(jīng)網(wǎng)絡(luò)方法相比,在訓(xùn)練過(guò)程中需考慮網(wǎng)絡(luò)之間的耦合誤差,本文引入網(wǎng)絡(luò)耦合系數(shù)kp以加強(qiáng)網(wǎng)絡(luò)之間的耦合,kp的取值過(guò)小會(huì)導(dǎo)致網(wǎng)絡(luò)耦合不佳,過(guò)大則可能會(huì)導(dǎo)致訓(xùn)練不收斂,因此需要對(duì)kp的取值進(jìn)行討論,選取kp=1,10,100,1000 四種情況對(duì)算例2 進(jìn)行計(jì)算分析(梯度系數(shù)m=2),以網(wǎng)絡(luò)訓(xùn)練時(shí)薄板的變形能、外力功變化情況衡量不同kp下本文模型的訓(xùn)練效果.

薄板變形能計(jì)算公式為

k為根據(jù)彎矩網(wǎng)絡(luò)的輸出由式(8)和式(9)計(jì)算所得.

外力所作實(shí)功為

其則根據(jù)撓度網(wǎng)絡(luò)的輸出值計(jì)算.

采用高斯積分,由上述公式計(jì)算神經(jīng)網(wǎng)絡(luò)訓(xùn)練時(shí)薄板的變形能、應(yīng)變能的變化情況,根據(jù)撓度網(wǎng)絡(luò)輸出計(jì)算外力實(shí)功Wext,根據(jù)彎矩網(wǎng)絡(luò)輸出計(jì)算薄板變形能Wint.當(dāng)神經(jīng)網(wǎng)絡(luò)解答收斂時(shí),根據(jù)能量原理,其外力實(shí)功應(yīng)與變形能相等.此時(shí)各工況下的變形能與外力實(shí)功在神經(jīng)網(wǎng)絡(luò)訓(xùn)練過(guò)程中的變化情況如圖14 所示.可以看出,隨著訓(xùn)練的進(jìn)行,外力實(shí)功與薄板的變形能逐漸收斂,在各個(gè)工況下,變形能的收斂速度相比于外力功要慢,進(jìn)而說(shuō)明撓度網(wǎng)絡(luò)收斂速度較快.當(dāng)kp=1 時(shí),此時(shí)薄板變形能與外力功收斂時(shí),兩者相差較大;而隨著kp的增大,神經(jīng)網(wǎng)絡(luò)計(jì)算的變形能收斂加快,此時(shí)當(dāng)變形能與外力功收斂時(shí),兩者之間的差別減小,而當(dāng)kp=1000 時(shí),其對(duì)訓(xùn)練的影響與kp=100 時(shí)差別不大,在實(shí)際計(jì)算中,該系數(shù)不能過(guò)大,過(guò)大會(huì)導(dǎo)致神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)不到正確的特征,進(jìn)而導(dǎo)致訓(xùn)練不收斂.

圖14 神經(jīng)網(wǎng)絡(luò)訓(xùn)練時(shí)不同 kp 下薄板變形能與外力功的變化情況Fig.14 Changes in the deformation energy and the external force work during the training process with different kp

5.2 隱藏層層數(shù)與神經(jīng)元個(gè)數(shù)

隱藏層層數(shù)與每層的神經(jīng)元個(gè)數(shù)的選取均會(huì)影響到神經(jīng)網(wǎng)絡(luò)的訓(xùn)練效率,不失一般性,在其他計(jì)算參數(shù)不變的情況下,本節(jié)選取算例2 中m=2 的情況分別討論隱藏層層數(shù)、神經(jīng)元個(gè)數(shù)的改變對(duì)計(jì)算過(guò)程中誤差函數(shù)收斂情況的影響,計(jì)算結(jié)果如圖15所示.

圖15 隱藏層數(shù)、每層神經(jīng)元個(gè)數(shù)對(duì)計(jì)算誤差的收斂影響對(duì)比Fig.15 Comparison of the effects of different number of hidden layers and neurons on the convergence of loss function

可以發(fā)現(xiàn),在一定程度內(nèi)誤差函數(shù)的收斂速度隨著網(wǎng)絡(luò)層數(shù)、隱藏層每層神經(jīng)元數(shù)的增大而加快,當(dāng)兩者都增大到一定程度時(shí),此時(shí)誤差的收斂速度會(huì)趨向于一個(gè)“飽和”狀態(tài).適當(dāng)?shù)卦黾訉訑?shù)或神經(jīng)元的個(gè)數(shù)有利于誤差函數(shù)的收斂,目前在利用神經(jīng)網(wǎng)絡(luò)方法求解偏微分方程的研究中,在兩者的選取上仍帶有一定的經(jīng)驗(yàn)性.結(jié)合上述的計(jì)算結(jié)果,考慮到計(jì)算機(jī)硬件能力的限制,本文算例的隱藏層層數(shù)在4～ 6 層間選取,每層神經(jīng)元數(shù)目在30～ 50 之間選取.

5.3 激活函數(shù)

在神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過(guò)程中,有多種非線性激活函數(shù)可以選擇,非線性的激活函數(shù)是使神經(jīng)網(wǎng)絡(luò)具備擬合非線性函數(shù)能力的重要原因,常用的激活函數(shù)有Tanh,ReLU,Sigmoid,Swish 函數(shù)等.為討論激活函數(shù)對(duì)神經(jīng)網(wǎng)絡(luò)訓(xùn)練的影響,其余計(jì)算參數(shù)不變,本節(jié)選取x2作為激活函數(shù)并與Tanh,Swish 函數(shù)進(jìn)行對(duì)比分析,對(duì)算例1、算例3 進(jìn)行計(jì)算,訓(xùn)練過(guò)程中的誤差走向如圖16 所示.在誤差函數(shù)的收斂速度上,x2優(yōu)于Tanh,Swish 函數(shù),同時(shí)由于其函數(shù)形式較為簡(jiǎn)單,在自動(dòng)微分計(jì)算中其所需計(jì)算量較小.本文經(jīng)驗(yàn)表明,在薄板彎曲問(wèn)題的求解上,采用多項(xiàng)式函 數(shù)x2作為激活函數(shù)可加快神經(jīng)網(wǎng)絡(luò)的收斂.

圖16 不同的激活函數(shù)對(duì)計(jì)算誤差的收斂影響對(duì)比Fig.16 Comparison of the effects of different activation functions on the convergence of loss function

6 結(jié)論與展望

本文基于深度學(xué)習(xí)技術(shù)與強(qiáng)形式的求解方案建立了一種直角坐標(biāo)下求解面內(nèi)變剛度薄板彎曲問(wèn)題的神經(jīng)網(wǎng)絡(luò)方法,通過(guò)幾個(gè)算例分析,得出以下結(jié)論:

(1)本文解答與理論解、有限元解吻合,證明了本文方法在求解面內(nèi)變剛度薄板彎曲問(wèn)題上的正確性,本文的神經(jīng)網(wǎng)絡(luò)模型不需要對(duì)彎曲剛度函數(shù)求偏導(dǎo),其適應(yīng)性更強(qiáng),同時(shí)在薄板的位移邊界條件、應(yīng)力邊界條件的施加上較為方便.

(2)本文方法屬于強(qiáng)形式的數(shù)值解法,其計(jì)算所得結(jié)果具備連續(xù)性與可導(dǎo)性.理論上,本文方法可以求解彈性模量以及厚度在面內(nèi)連續(xù)變化的薄板彎曲問(wèn)題.彎矩網(wǎng)絡(luò)的求解受到梯度系數(shù)的影響,在梯度變化較大處彎矩網(wǎng)絡(luò)的求解精度受到一定的影響,但對(duì)撓度網(wǎng)絡(luò)的求解精度影響不大.

(3)由于神經(jīng)網(wǎng)絡(luò)方法為迭代類解法,本文方法在薄板線性彎曲問(wèn)題求解上的收斂速度較有限元慢,但其計(jì)算所需內(nèi)存較小.通過(guò)本文的模型結(jié)構(gòu)可看出,神經(jīng)網(wǎng)絡(luò)方法具備相當(dāng)大的靈活性,根據(jù)這一特點(diǎn),可進(jìn)一步發(fā)展求解面內(nèi)變剛度功能梯度薄板非線性彎曲問(wèn)題的神經(jīng)網(wǎng)絡(luò)方法,神經(jīng)網(wǎng)絡(luò)方法在非線性問(wèn)題的求解中具備潛在優(yōu)勢(shì).

(4)本文模型仍存在優(yōu)化空間:一方面在本文模型的訓(xùn)練過(guò)程中,誤差函數(shù)及其梯度的計(jì)算在整個(gè)訓(xùn)練過(guò)程中占據(jù)大部分的時(shí)間,可以考慮優(yōu)化誤差函數(shù)的構(gòu)建過(guò)程,如引入有限元中形函數(shù)的思想對(duì)算法進(jìn)行優(yōu)化;另一方面為了使得撓度、彎矩網(wǎng)絡(luò)具備較強(qiáng)的表達(dá)能力,本文模型采用了兩個(gè)具有獨(dú)立參數(shù)的網(wǎng)絡(luò)進(jìn)行計(jì)算,這導(dǎo)致了本文模型的訓(xùn)練參數(shù)較多,為此后續(xù)優(yōu)化中可將本文的兩個(gè)網(wǎng)絡(luò)合并為一個(gè)網(wǎng)絡(luò)(2 個(gè)輸入,4 個(gè)輸出),對(duì)網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行改良,以減少訓(xùn)練參數(shù).