外部圓形裂紋軸對(duì)稱結(jié)構(gòu)彈性分析的多神經(jīng)網(wǎng)絡(luò)聯(lián)合訓(xùn)練方法1)

賀 云 杜 娟 李海濱

*(內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院，呼和浩特 010018)

?(內(nèi)蒙古工業(yè)大學(xué)理學(xué)院，呼和浩特 010051)

**(內(nèi)蒙古財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與數(shù)學(xué)學(xué)院，呼和浩特 010070)

軸對(duì)稱結(jié)構(gòu)裂紋彈性力學(xué)分析是工程實(shí)踐中重要而基礎(chǔ)的問題。傳統(tǒng)的解析求解方法有按位移求解、按應(yīng)力求解和混合求解[1]。這三種方法都要從靜力學(xué)方面、幾何學(xué)方面和物理學(xué)方面考慮。根據(jù)彈性理論可知，彈性力學(xué)問題可以歸結(jié)為具有邊界條件的偏微分方程求解問題。按應(yīng)力求解該問題時(shí)，傳統(tǒng)的求解方法有逆解法和半逆解法[2]。此類方法都需要針對(duì)所要求解的問題，假設(shè)應(yīng)力函數(shù)為參數(shù)未知的某種函數(shù)形式，將其代入應(yīng)力函數(shù)表示的應(yīng)力分量表達(dá)式中。根據(jù)相容方程及邊界條件來考察假設(shè)的應(yīng)力函數(shù)，得到假定應(yīng)力函數(shù)的待求參數(shù)，從而得到應(yīng)力函數(shù)及應(yīng)力分量。但由于彈性體裂紋形狀及邊界條件的復(fù)雜性，一般很難得到彈性力學(xué)問題的解析解。眾所周知，己知描述物理問題的微分方程模型中，只有極少數(shù)部分是可以求出解析解，而絕大多數(shù)只能借助近似數(shù)值方法求解，彈性力學(xué)模型亦不例外，因此尋找有效數(shù)值解法就顯得尤為重要，而有限單元法是迄今為止在科學(xué)與工程計(jì)算領(lǐng)域最活躍、應(yīng)用最廣泛的數(shù)值計(jì)算方法[3-4]。

基于基本的有限單元法，后來一些學(xué)者又提出了一些發(fā)展、改進(jìn)的方法，如穩(wěn)定混合有限單元法[5]、h型混合變階有限單元法[6]、無閉鎖非協(xié)調(diào)有限單元法[7]，自適應(yīng)有限單元法[8]、間斷有限單元法[9]等。但在實(shí)際應(yīng)用中，有許多因素使得彈性力學(xué)問題產(chǎn)生強(qiáng)奇性解，如物理區(qū)域是一個(gè)非光滑的幾何區(qū)域（凹角或裂縫等）、材料常數(shù)不連續(xù)、應(yīng)力集中等。雖然這些特殊彈性力學(xué)問題在進(jìn)行有限單元法求解時(shí)可以通過加密網(wǎng)格來克服求解困難，但是網(wǎng)格加密會(huì)引起計(jì)算量急劇增加。這是因?yàn)橛删W(wǎng)格加密所帶來的自由度成指數(shù)增長，導(dǎo)致計(jì)算量成指數(shù)增加。

由于軸對(duì)稱結(jié)構(gòu)裂紋彈性分析可以歸結(jié)為求解具有邊界條件的偏微分方程，神經(jīng)網(wǎng)絡(luò)方法在此方面也有廣泛的應(yīng)用。Lagaris等[10]研究了通過神經(jīng)網(wǎng)絡(luò)求解含有初始邊界條件的偏微分方程，初始邊界條件的神經(jīng)網(wǎng)絡(luò)不含可調(diào)參數(shù)，偏微分方程為一個(gè)前饋型神經(jīng)網(wǎng)絡(luò)，將此方法的求解結(jié)果與有限單元法進(jìn)行對(duì)比，在得到良好結(jié)果的同時(shí)，執(zhí)行速度有了較大的提升。針對(duì)含初始邊界條件的偏微分方程與神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)及訓(xùn)練集無法有機(jī)結(jié)合在一起的問題，Aarts等[11]提出了一種基于遺傳算法的具有特殊結(jié)構(gòu)神經(jīng)網(wǎng)絡(luò)的求解偏微分方程方法，得到了較好的結(jié)果。張永恒等[12]通過基于Hopflied神經(jīng)網(wǎng)絡(luò)的優(yōu)化算法來實(shí)現(xiàn)打靶法的過程，從而求解了具有初始邊界條件的微分方程，此方法具有對(duì)初值選擇要求較低，而計(jì)算精度較高等優(yōu)點(diǎn)。柴俊霞等[13]設(shè)計(jì)了以多元二次徑向基神經(jīng)網(wǎng)絡(luò)為求解單元的偏微分方程求解方法，并通過算例驗(yàn)證了方法的有效性。Mosleh等[14]運(yùn)用模糊神經(jīng)網(wǎng)絡(luò)求解具有初始值問題的模糊微分方程，并在一個(gè)工程問題中加以運(yùn)用。Guner等[15]運(yùn)用exp函數(shù)方法求解分?jǐn)?shù)階微分方程，同樣得到了較好的結(jié)果。然而，傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)[16-17]求解偏微分方程的難點(diǎn)在于，單個(gè)神經(jīng)網(wǎng)絡(luò)無法同時(shí)滿足偏微分方程及其邊界條件，而訓(xùn)練多個(gè)神經(jīng)網(wǎng)絡(luò)又難以體現(xiàn)神經(jīng)網(wǎng)絡(luò)之間的聯(lián)系，且傳統(tǒng)神經(jīng)網(wǎng)絡(luò)求解偏微分方程問題工作效率低、精度無法保證。

為此，本文在彈性體模型解析求解方法的基礎(chǔ)上，擬提出利用多神經(jīng)網(wǎng)絡(luò)（multi-neural network，MNN）聯(lián)合訓(xùn)練方法實(shí)現(xiàn)軸對(duì)稱結(jié)構(gòu)裂紋彈性力學(xué)分析。首先，對(duì)軸對(duì)稱結(jié)構(gòu)裂紋進(jìn)行彈性力學(xué)描述，由于按應(yīng)力求解彈性力學(xué)應(yīng)力問題可以簡化為應(yīng)力函數(shù)表示的相容方程及其邊界條件的求解問題，所以針對(duì)該問題，隨后將應(yīng)力函數(shù)假設(shè)為統(tǒng)一的神經(jīng)網(wǎng)絡(luò)形式，構(gòu)造應(yīng)力函數(shù)表示的相容方程及應(yīng)力邊界條件的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，從而組成神經(jīng)網(wǎng)絡(luò)系統(tǒng)。運(yùn)用MNN聯(lián)合訓(xùn)練方法訓(xùn)練該神經(jīng)網(wǎng)絡(luò)系統(tǒng)，從而得到應(yīng)力分量的神經(jīng)網(wǎng)絡(luò)數(shù)值解。最后采用兩個(gè)典型算例的計(jì)算和對(duì)比分析驗(yàn)證所提方法的有效性。

1 軸對(duì)稱結(jié)構(gòu)裂紋彈性力學(xué)描述

外部圓形裂紋，裂紋半徑為c，無限彈性體任意半徑為ρ，轉(zhuǎn)過任意角度為φ，如圖1所示。本文將在彈性力學(xué)理論基礎(chǔ)上，用神經(jīng)網(wǎng)絡(luò)方法研究無限彈性體中含有一個(gè)外部圓形裂紋的軸對(duì)稱力學(xué)問題。對(duì)于上述問題，當(dāng)裂紋上下表面作用的外應(yīng)力自相平衡時(shí)，可化歸為具有應(yīng)力邊界的平面軸對(duì)稱結(jié)構(gòu)問題，其精確解可由彈性力學(xué)理論求得。對(duì)于受非對(duì)稱載荷作用的復(fù)雜邊界情況，擬給出裂紋面的神經(jīng)網(wǎng)絡(luò)求解方法。

圖1 外部圓形裂紋Fig.1 External circular crack

根據(jù)軸對(duì)稱結(jié)構(gòu)彈性理論，可以歸結(jié)為求解平衡微分方程、幾何方程和物理方程的方程組。彈性體的平衡微分方程可以表示為

式中，σρ和σφ為極坐標(biāo)系下的正應(yīng)力分量，τρφ為極坐標(biāo)系下的切應(yīng)力分量，fρ和fφ為極坐標(biāo)系下的體力分量。

幾何方程為

式中，ερ和εφ為極坐標(biāo)系下的正應(yīng)變分量；γρφ為極坐標(biāo)系下的切應(yīng)變；uρ和uφ為極坐標(biāo)系下的位移分量。

對(duì)于應(yīng)力問題，物理方程為

式中，E為彈性模量，μ為泊松比。

至此，彈性力學(xué)應(yīng)力問題可以轉(zhuǎn)化成求解滿足具有邊界條件的相容方程。該相容方程由應(yīng)力函數(shù)表示，通過應(yīng)力函數(shù)的表達(dá)式，從而可以計(jì)算得到應(yīng)力分量。

對(duì)于彈性體模型應(yīng)力問題，根據(jù)彈性理論，用徑向坐標(biāo)ρ和環(huán)向坐標(biāo)φ表示的極坐標(biāo)系下相容方程為

應(yīng)力邊界條件為

當(dāng)不計(jì)體力時(shí)，通過應(yīng)力函數(shù)φ可以求得極坐標(biāo)系下應(yīng)力分量為

2 MNN聯(lián)合訓(xùn)練方法

對(duì)于包含未知函數(shù)及其偏導(dǎo)數(shù)的等式，一般說來，如果x1,x2,...,xn為變量，以u(píng)? 為未知函數(shù)，對(duì)任意變量xi最高階偏導(dǎo)數(shù)不超過Ai的線性偏微分方程的一般形式可表達(dá)為

在偏微分方程問題中往往還帶有邊界條件，可表示為

由神經(jīng)網(wǎng)絡(luò)相關(guān)理論可知，一個(gè)多層神經(jīng)網(wǎng)絡(luò)可以逼近任意非線性函數(shù)，這對(duì)于無法用解析方法求解及數(shù)值求解較困難的偏微分方程是一種思路。因此，可以將偏微分方程的通解u? 寫為一個(gè)參數(shù)待定的帶權(quán)數(shù)的隱式多項(xiàng)式

根據(jù)如上理論，軸對(duì)稱結(jié)構(gòu)應(yīng)力函數(shù)φ(ρ,φ)表示為變量ρ和φ的帶權(quán)數(shù)的隱式多項(xiàng)式，此式結(jié)構(gòu)和神經(jīng)網(wǎng)絡(luò)表達(dá)式結(jié)構(gòu)類似，即

網(wǎng)絡(luò)隱層單元激活函數(shù)f(x)=ex， 將式（15）代入式（12）～式（14），得到各應(yīng)力分量表達(dá)式為

將式（15）代入式（9），并將式（16）～式（18）代入式（10）和式（11），得到極坐標(biāo)下偏微分方程及其邊界條件的神經(jīng)網(wǎng)絡(luò)表達(dá)式為

式（19）～式（21）就構(gòu)成了具有邊界條件偏微分方程的神經(jīng)網(wǎng)絡(luò)求解表達(dá)式。比較式（19）與式（20）、式（21）對(duì)應(yīng)的神經(jīng)網(wǎng)絡(luò)可知，各網(wǎng)絡(luò)在結(jié)構(gòu)上均為兩輸入、單輸出、n個(gè)隱層單元的三層前饋型神經(jīng)網(wǎng)絡(luò)。網(wǎng)絡(luò)結(jié)構(gòu)如圖2～圖4所示。

圖2 式（19）對(duì)應(yīng)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.2 The neural network structure of Eq.(19)

圖4 式（21）對(duì)應(yīng)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.4 The neural network structure of Eq.(21)

根據(jù)式（19）～式（21）構(gòu)造極坐標(biāo)系下的神經(jīng)網(wǎng)絡(luò)系統(tǒng)。網(wǎng)絡(luò)輸入均為ρ和φ，網(wǎng)絡(luò)輸出分別為h2 ，fˉρ，fˉφ。此外，各網(wǎng)絡(luò)的輸入層到隱層單元連接權(quán)值為kij，其中j=1,2 ，隱層單元閾值為bi。式（19）～式（21）對(duì)應(yīng)的神經(jīng)網(wǎng)絡(luò)與式（15）對(duì)應(yīng)的網(wǎng)絡(luò)隱層到輸出層單元連接權(quán)值之間存在的關(guān)系為

欲通過學(xué)習(xí)使上述三個(gè)神經(jīng)網(wǎng)絡(luò)同時(shí)收斂，則三個(gè)神經(jīng)網(wǎng)絡(luò)構(gòu)成網(wǎng)絡(luò)系統(tǒng)的總性能函數(shù)誤差收斂，其中總性能函數(shù)可定義為式（25）。

神經(jīng)網(wǎng)絡(luò)系統(tǒng)中ki1，ki2，bi，wi為獨(dú)立變量，令x=[ki1,ki2,bi,wi]T為網(wǎng)絡(luò)系統(tǒng)待調(diào)節(jié)列向量。系統(tǒng)性能函數(shù)寫為

系統(tǒng)性能函數(shù)E? 為描述神經(jīng)網(wǎng)絡(luò)系統(tǒng)整體誤差的函數(shù)，其值越小，神經(jīng)網(wǎng)絡(luò)計(jì)算精度越高，隱層神經(jīng)元激活函數(shù)為f(x)=ex，文獻(xiàn)[18]驗(yàn)證了隱層單元激活函數(shù)為指數(shù)函數(shù)時(shí)，較logsig和tansig等傳統(tǒng)激活函數(shù)，同樣具有較好的計(jì)算精度和泛化能力。

圖3 式（20）對(duì)應(yīng)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.3 The neural network structure of Eq.(20)

由于運(yùn)用梯度下降法訓(xùn)練網(wǎng)絡(luò)時(shí)，穩(wěn)定性要求學(xué)習(xí)率很小，所以使得訓(xùn)練很慢。動(dòng)量梯度法速度雖然有所提高，但還是無法達(dá)到實(shí)際應(yīng)用的要求。牛頓法收斂速度快，但是在每一次迭代中，需要求出Hessian矩陣，其中包括性能函數(shù)的二次導(dǎo)數(shù)，這就使得計(jì)算量變得很大。還有一些神經(jīng)網(wǎng)絡(luò)的改進(jìn)訓(xùn)練算法[19]，如動(dòng)量梯度下降法、牛頓法、共軛梯度法、 Levenberg-Marquardt（L-M）算法等。

如果神經(jīng)網(wǎng)絡(luò)的性能函數(shù)為平方和的形式，就可以采用L-M訓(xùn)練算法。此算法是神經(jīng)網(wǎng)絡(luò)計(jì)算中較為常用的算法，不需要計(jì)算Hessian陣，Hessian陣可以用矩陣H=JTJ來近似替換。其中，J為網(wǎng)絡(luò)性能函數(shù)E? 對(duì)權(quán)值、閾值一階導(dǎo)數(shù)的雅可比矩陣，它是關(guān)于權(quán)值、閾值的函數(shù)。該矩陣的求解相比Hessian矩陣求解要簡單很多。

L-M訓(xùn)練算法中網(wǎng)絡(luò)參數(shù)通過式xr+1=xr-[JTJ+μI]-1JTE? 進(jìn)行學(xué)習(xí)，式中I為單位矩陣，μ為學(xué)習(xí)率。

L-M訓(xùn)練算法是牛頓法與梯度下降法的結(jié)合，當(dāng)學(xué)習(xí)率μ很小時(shí)，相當(dāng)于牛頓法的迭代步長。當(dāng)μ很大時(shí)，相當(dāng)于梯度下降法的迭代步長。在迭代過程中，如果訓(xùn)練成功，就減小μ的值；如果訓(xùn)練失敗，就增加μ的值，這樣使性能函數(shù)的值逐漸減小。

因此，本文方法結(jié)合L-M訓(xùn)練算法，運(yùn)用上述MNN聯(lián)合訓(xùn)練方法訓(xùn)練神經(jīng)網(wǎng)絡(luò)系統(tǒng)，得到神經(jīng)網(wǎng)絡(luò)參數(shù)。將網(wǎng)絡(luò)參數(shù)代入式（16）～式（18），從而得到極坐標(biāo)系下應(yīng)力分量的數(shù)值解。

3 算例分析

3.1 算例1

圓環(huán)分別受到均布外壓力qb= 2 kN/m，均布內(nèi)壓力qa= 1 kN/m。外徑R= 3 cm，內(nèi)徑r= 2 cm，如圖5，求模型中各點(diǎn)的應(yīng)力分量。

圖5 圓環(huán)結(jié)構(gòu)Fig.5 Ring structure

對(duì)于圓環(huán)軸對(duì)稱問題，相容方程可簡化為

邊界條件為

對(duì)于該問題，由于應(yīng)力函數(shù)φ只是半徑ρ的函數(shù)，所以φ的神經(jīng)網(wǎng)絡(luò)形式可簡化為

將式（29）代入相容性方程及邊界條件，構(gòu)造神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)及訓(xùn)練樣本。由于邊界條件式（27）自動(dòng)滿足，所以訓(xùn)練樣本如表1所示。

表1 神經(jīng)網(wǎng)絡(luò)系統(tǒng)訓(xùn)練樣本Table 1 Training samples of neural network system

隱層單元激活函數(shù) ，經(jīng)過調(diào)節(jié)，算例中隱層單元個(gè)數(shù)n= 10。運(yùn)用MNN聯(lián)合訓(xùn)練方法，聯(lián)合訓(xùn)練神經(jīng)網(wǎng)絡(luò)系統(tǒng)20步，網(wǎng)絡(luò)系統(tǒng)收斂，誤差收斂曲線如圖6所示。

圖6 算例1誤差收斂曲線Fig.6 Error convergence curve in example 1

根據(jù)簡化后的應(yīng)力函數(shù)神經(jīng)網(wǎng)絡(luò)形式式（29），得到徑向應(yīng)力分量和環(huán)向應(yīng)力分量的表達(dá)式為

圓環(huán)半徑在[2,3]范圍內(nèi)均勻取11個(gè)點(diǎn)作為測(cè)試樣本，將其代入式（30）和式（31），得到測(cè)試樣本點(diǎn)處應(yīng)力分量的神經(jīng)網(wǎng)絡(luò)數(shù)值解，將徑向應(yīng)力的數(shù)值列于表2中。

根據(jù)彈性力學(xué)解析求解方法得到徑向應(yīng)力的表達(dá)式為

通過式（32）得到該平面問題的彈性力學(xué)解析解作為理論解。為了進(jìn)一步對(duì)比結(jié)果，由于圓環(huán)的對(duì)稱性，取圓環(huán)的1/4進(jìn)行分析，通過Ansys軟件求得該問題的有限單元法解。圖7為計(jì)算精度和時(shí)間隨網(wǎng)格數(shù)量的變化曲線，曲線1表示結(jié)構(gòu)中計(jì)算精度隨網(wǎng)格數(shù)量收斂的一般曲線，曲線2表示計(jì)算時(shí)間隨網(wǎng)格數(shù)量的變化。從圖中可看出網(wǎng)格較少時(shí)，增加網(wǎng)格數(shù)量可以使計(jì)算精度明顯提高，而計(jì)算時(shí)間不會(huì)有大的增加。當(dāng)網(wǎng)格數(shù)量增加到一定程度后，再繼續(xù)增加網(wǎng)格時(shí)，精度提高甚微，而計(jì)算時(shí)間卻有大幅度增加。所以為了增加網(wǎng)格的經(jīng)濟(jì)性，算例中取計(jì)算精度大致收斂時(shí)的網(wǎng)格數(shù)量，P= 160個(gè)。兩種方法的計(jì)算結(jié)果一同列于表2。徑向應(yīng)力的神經(jīng)網(wǎng)絡(luò)方法解與理論解對(duì)比如圖8所示，有限單元法解的徑向應(yīng)力云圖如圖9所示。

圖7 有限單元法計(jì)算精度和時(shí)間曲線Fig.7 Calculation accuracy and time curves of finite element method

圖9 有限單元法的徑向應(yīng)力云圖Fig.9 Radial stress cloud figure of finite element method

表2 徑向應(yīng)力計(jì)算結(jié)果Table 2 Calculation results of radial stress

由表2結(jié)果可以看出，本文方法總體的計(jì)算精度高于有限單元法。并且從圖8中可以看出，在越接近邊界處，本文方法的精度越高，而有限單元法在邊界處的計(jì)算精度卻較差。

圖8 徑向應(yīng)力對(duì)比圖Fig.8 Comparition figure of radial stress

3.2 算例 2

矩形薄板長a= 8 cm，高b= 7 cm，在離開邊界較遠(yuǎn)處有半徑為r=0.5cm 的小圓孔規(guī)則裂紋。上下及左右邊界分別受到均布?jí)毫屠?，大小均為q= 1 kN/m，如圖10所示，求圖中半徑為R的大圓內(nèi)各點(diǎn)的應(yīng)力分量。

圖10 帶圓孔的矩形薄板結(jié)構(gòu)Fig.10 Structure of rectangular thin plate with round hole

該問題為孔口應(yīng)力集中問題，具有局部性，一般孔口的應(yīng)力集中區(qū)域約在距孔邊1.5倍圓孔直徑的范圍內(nèi)。為了不失一般性，本算例取距孔邊為2.5倍圓孔直徑的大圓（虛線所示）范圍內(nèi)進(jìn)行分析。在大圓外部，應(yīng)力情況與無孔時(shí)基本相同，本算例將不作計(jì)算。該問題的相容方程為式（9）。

邊界條件為

運(yùn)用本文方法構(gòu)造網(wǎng)絡(luò)結(jié)構(gòu)及訓(xùn)練樣本，式（9）、式（33）、式（34）分別對(duì)應(yīng)的神經(jīng)網(wǎng)絡(luò)訓(xùn)練樣本如表3～表5所示。

表3 式（9）對(duì)應(yīng)的網(wǎng)絡(luò)訓(xùn)練樣本Table 3 Training samples of Eq.(9) network

表4 式（33）對(duì)應(yīng)的網(wǎng)絡(luò)訓(xùn)練樣本Table 4 Training samples of Eq.(33) network

表5 式（34）對(duì)應(yīng)的網(wǎng)絡(luò)訓(xùn)練樣本Table 5 Training samples of Eq.(34) network

隱層單元激活函數(shù)f=ex，經(jīng)過調(diào)節(jié)，算例中隱層單元個(gè)數(shù)n= 50。運(yùn)用MNN聯(lián)合訓(xùn)練方法，聯(lián)合訓(xùn)練神經(jīng)網(wǎng)絡(luò)系統(tǒng)20步，網(wǎng)絡(luò)系統(tǒng)收斂，誤差收斂曲線如圖11所示。

圖11 算例2誤差收斂曲線Fig.11 Error convergence curve in example 2

將測(cè)試樣本代入式（16）～式（18）， 得到測(cè)試樣本點(diǎn)處的徑向應(yīng)力σρ，環(huán)向應(yīng)力σφ和切應(yīng)力τρφ，如表6～表8所示。

表6 測(cè)試樣本點(diǎn)處徑向應(yīng)力 σρ (單位：kPa)Table 6 Radial stress σ ρ of test sample points (unit: kPa)

表8 測(cè)試樣本點(diǎn)處切應(yīng)力 (單位：kPa)Table 8 Shear stress of test sample points (unit: kPa)

表8 測(cè)試樣本點(diǎn)處切應(yīng)力 (單位：kPa)Table 8 Shear stress of test sample points (unit: kPa)

images/BZ_168_343_2726_422_2759.png images/BZ_168_1338_2693_1430_2727.png0 0.157 1 0.314 2 …… 1.413 7 1.570 8 0.500 0 -0.027 0 -0.021 0 -0.067 6 …… -0.068 6 -0.011 2 0.750 0 -0.005 5 -0.359 4 -0.703 1 …… -0.390 9 -0.118 5 1.000 0 -0.007 2 -0.398 8 -0.711 0 …… -0.417 7 -0.121 7…… …… …… …… …… …… ……2.750 0 -0.001 9 -0.360 2 -0.607 0 …… -0.271 2 -0.113 5 3.000 0 -0.003 9 -0.328 5 -0.605 5 …… -0.310 1 -0.128 2

根據(jù)本文方法得到的徑向應(yīng)力σρ，環(huán)向應(yīng)力σφ和切應(yīng)力τρφ在測(cè)試樣本點(diǎn)處的大小分布如圖12所示。

圖12 算例2本文方法的應(yīng)力分布圖Fig.12 Stress distribution figure of the proposed method in example 2

為了便于進(jìn)行對(duì)比，利用有限單元法求解該問題，根據(jù)圖7的精度收斂曲線，本算例取網(wǎng)格數(shù)量為P= 135個(gè)，并在應(yīng)力集中處進(jìn)行網(wǎng)格加密處理。在Ansys軟件中的應(yīng)力云圖如圖13所示。

表7 測(cè)試樣本點(diǎn)處環(huán)向應(yīng)力 (單位：kPa)Table 7 Hoop stress of test sample points (unit: kPa )

表7 測(cè)試樣本點(diǎn)處環(huán)向應(yīng)力 (單位：kPa)Table 7 Hoop stress of test sample points (unit: kPa )

images/BZ_168_346_2032_425_2066.png images/BZ_168_1341_2000_1433_2033.png0 0.157 1 0.314 2 …… 1.413 7 1.570 8 0.500 0 -4.031 2 -3.744 3 -3.313 7 …… 3.758 5 4.023 0 0.750 0 -1.621 4 -1.441 0 -1.252 2 …… 1.468 7 1.646 9 1.000 0 -1.166 9 -1.188 0 -1.007 6 …… 1.053 7 1.157 5…… …… …… …… …… …… ……2.750 0 -1.061 4 -0.903 2 -0.822 8 …… 0.904 8 0.990 9 3.000 0 -0.945 4 -0.954 8 -0.760 3 …… 0.920 8 0.930 2

根據(jù)彈性力學(xué)半逆解法，得到應(yīng)力分量的求解表達(dá)式為

通過式（35）～式（37）可以求得該平面問題應(yīng)力分量的彈性力學(xué)解析解作為理論解。有限單元法、本文方法求得的數(shù)值解與理論解相比的平均絕對(duì)誤差列于表9。

表9 測(cè)試樣本點(diǎn)處應(yīng)力的平均絕對(duì)誤差Table 9 Average absolute error of test sample points stress

由表9結(jié)果可以看出，對(duì)于徑向應(yīng)力和切應(yīng)力，有限單元法和本文方法的計(jì)算誤差相差不大。而對(duì)于應(yīng)力集中較為嚴(yán)重的環(huán)向應(yīng)力，本文方法的計(jì)算精度高于有限單元法。與傳統(tǒng)有限單元法求解時(shí)，需要對(duì)物理模型進(jìn)行網(wǎng)格劃分相比，本文方法從彈性力學(xué)的基本方程出發(fā)，引入邊界條件方程可以直接求解。

3.3 算例 3

無限大彈性體，在上邊界處有半徑為a=0.2cm的半圓孔裂紋。上邊界受到相反的兩個(gè)力，大小均為F= 10 kN，如圖14所示，分析圖中半徑為b=0.4cm的半圓內(nèi)各點(diǎn)的應(yīng)力分量。

圖14 帶半圓孔裂紋的彈性體結(jié)構(gòu)Fig.14 Structure of elastomer with semi-circular hole crack

根據(jù)相容方程式（9）和式（38），式（39）的邊界條件，構(gòu)造MNN結(jié)構(gòu)，取對(duì)稱結(jié)構(gòu)的1/2進(jìn)行分析，運(yùn)用本文方法計(jì)算得到應(yīng)力分量，如圖15所示。

圖15 算例3本文方法的應(yīng)力分布圖Fig.15 Stress distribution figure of the proposed method in example 3

利用有限單元法求解該問題，根據(jù)精度收斂曲線，本算例取網(wǎng)格數(shù)量為P= 1 400個(gè)，網(wǎng)格劃分如圖16所示。

圖16 有限元網(wǎng)格劃分圖Fig.16 The figure of finite element mesh division

應(yīng)力分量的理論解求解表達(dá)式為

有限單元法、本文方法求得的數(shù)值解與理論解相比的平均絕對(duì)誤差列于表10。

由表10結(jié)果可以看出：（1）相比算例2，網(wǎng)格劃分加密后，有限單元法計(jì)算精度略有提高；（2）對(duì)于理論解表達(dá)式非線性較強(qiáng)時(shí)，本文方法精度略有降低，但是對(duì)于徑向應(yīng)力和切應(yīng)力，兩種方法的計(jì)算誤差相差不大；（3）對(duì)于應(yīng)力集中較為嚴(yán)重的環(huán)向應(yīng)力，本文方法在計(jì)算精度方面的優(yōu)越性明顯高于有限單元法。

表10 測(cè)試樣本點(diǎn)處應(yīng)力的平均絕對(duì)誤差Table 10 Average absolute error of test sample points stress

3.4 算例4

載人潛水器的耐壓球殼直徑設(shè)計(jì)為2.1 m，殼體上有一個(gè)小端直徑為200 mm的主觀察窗，兩個(gè)小端直徑為120 mm的側(cè)觀察窗對(duì)稱分布于主觀察窗兩側(cè)。主側(cè)觀察窗的厚度為234 mm、均采用45o角錐臺(tái)形結(jié)構(gòu)，結(jié)構(gòu)尺寸如圖17所示，外部邊界處的壓強(qiáng)為162 MPa。

圖17 觀察窗結(jié)構(gòu)圖（單位：mm）Fig.17 Structure of observation window (unit: mm)

此軸對(duì)稱結(jié)構(gòu)觀察窗和窗座的接觸面可看作為一個(gè)外部圓形裂紋，由于此處容易出現(xiàn)應(yīng)力集中，從而導(dǎo)致觀察窗的結(jié)構(gòu)破壞，該算例將對(duì)此處進(jìn)行力學(xué)分析。根據(jù)相容方程和邊界條件，構(gòu)造MNN結(jié)構(gòu)，取對(duì)稱結(jié)構(gòu)的1/2進(jìn)行分析，運(yùn)用本文方法計(jì)算得到應(yīng)力分量。作為對(duì)比，利用有限單元法求解該問題，根據(jù)精度收斂曲線，本算例取網(wǎng)格數(shù)量為P= 400個(gè)，應(yīng)力云圖如圖18所示，本文方法與有限單元法在測(cè)試點(diǎn)處的結(jié)果對(duì)比如圖19所示。

圖18 觀察窗有限單元法應(yīng)力云圖Fig.18 Stress cloud of finite element method in observation window

測(cè)試樣本點(diǎn)處有限單元法、本文方法求得的數(shù)值解及最大相對(duì)誤差列于表11。

由表11結(jié)果可以看出：本文方法相比有限單元法，在樣本點(diǎn)處徑向應(yīng)力、環(huán)向應(yīng)力的最大相對(duì)誤差均小于4%，從而說明對(duì)于復(fù)雜結(jié)構(gòu)，本文方法有較好的適用性。從圖19中可以看出：結(jié)構(gòu)邊界處兩種方法相差較大，這是由于有限單元法在邊界處的精度較差，而本文方法更接近于理論值。

表11 測(cè)試樣本點(diǎn)處的應(yīng)力Table 11 The stress of test sample points

圖19 本文方法和有限元方法應(yīng)力對(duì)比圖Fig.19 Stress comparison between the proposed method and finite element method

4 結(jié)論

本文提出了軸對(duì)稱結(jié)構(gòu)外部圓形裂紋彈性力學(xué)分析的神經(jīng)網(wǎng)絡(luò)方法。假設(shè)應(yīng)力函數(shù)為神經(jīng)網(wǎng)絡(luò)形式，通過相容方程及邊界條件與應(yīng)力函數(shù)的導(dǎo)數(shù)關(guān)系，分別構(gòu)造參數(shù)相關(guān)聯(lián)的多個(gè)神經(jīng)網(wǎng)絡(luò)，通過聯(lián)合訓(xùn)練方法訓(xùn)練網(wǎng)絡(luò)系統(tǒng)，實(shí)現(xiàn)神經(jīng)網(wǎng)絡(luò)的高精度數(shù)值求解。數(shù)值算例表明，本文方法對(duì)包含應(yīng)力集中在內(nèi)的彈性力學(xué)分析有較好的適用性，在計(jì)算精度方面，相比有限單元法的計(jì)算結(jié)果有所提高。由于本文介紹的神經(jīng)網(wǎng)絡(luò)方法直接從理論求解彈性力學(xué)問題的偏微分方程出發(fā)，所以相比傳統(tǒng)有限單元法的求解過程，本文方法在數(shù)值求解時(shí)無需對(duì)彈性體物理模型進(jìn)行網(wǎng)格劃分。

本文貢獻(xiàn)在于：（1）將應(yīng)力函數(shù)寫成統(tǒng)一的神經(jīng)網(wǎng)絡(luò)形式，無需根據(jù)具體問題假設(shè)應(yīng)力函數(shù)，從而克服了應(yīng)力函數(shù)的求解困難；（2）通過神經(jīng)網(wǎng)絡(luò)方法的求解，使軸對(duì)稱結(jié)構(gòu)外部圓形裂紋彈性力學(xué)分析的數(shù)值解精度更高，對(duì)復(fù)雜形狀及復(fù)雜邊界的彈性體力學(xué)分析提供了一種新的思路。另外，本文方法也可以拓展到按位移求解、混合求解的彈性力學(xué)分析中。