賀 云,李海濱,杜 娟
(1.內(nèi)蒙古農(nóng)業(yè)大學(xué) 水利與土木建筑工程學(xué)院,呼和浩特 010018; 2.內(nèi)蒙古工業(yè)大學(xué) 理學(xué)院,呼和浩特 010051;3.內(nèi)蒙古財經(jīng)大學(xué) 統(tǒng)計與數(shù)學(xué)學(xué)院,呼和浩特 010070)
對于固體火箭發(fā)動機(jī)藥柱粘彈性理論,國內(nèi)外學(xué)者都進(jìn)行了大量研究。在國外,CHARLES系統(tǒng)描述了固體推進(jìn)劑粘彈性表征[1];JUNG等提出了一種考慮大變形及循環(huán)載荷對發(fā)動機(jī)裝藥影響的非線性本構(gòu)方程[2];HIMANSHU等提出固體推進(jìn)劑Maxwell模型,可以準(zhǔn)確預(yù)測推進(jìn)劑材料系數(shù)[3]。在國內(nèi),張海聯(lián)等建立Lagrangian法隨機(jī)虛功方程,研究了藥柱的具有隨機(jī)參數(shù)結(jié)構(gòu)響應(yīng)問題[4];張書俊等對固體發(fā)動機(jī)藥柱進(jìn)行粘彈性動態(tài)可靠度分析與不可壓縮粘彈性隨機(jī)有限元分析[5];張建偉等通過局部平均方法對隨機(jī)場進(jìn)行離散,得到隨機(jī)粘彈性有限元方程[6];田四朋根據(jù)Herrmann變分原理和對應(yīng)原理出發(fā),推導(dǎo)得到了不可壓和近似不可壓粘彈性材料積分型本構(gòu)關(guān)系[7];職世君根據(jù)藥柱在固化降溫情況下的結(jié)構(gòu)完整性分析,在藥柱施加點火沖擊的壓力載荷時,分別進(jìn)行動態(tài)和準(zhǔn)靜態(tài)計算,得出較好的結(jié)果[8]。
綜上所述,對于固體火箭發(fā)動機(jī)藥柱粘彈性結(jié)構(gòu)的力學(xué)性能分析,在線性粘彈性理論中對應(yīng)原理分析方法較為成熟,是一種行之有效的方法,可以大大簡化粘彈性問題的求解過程[9-14],但是對于非線性問題往往更貼近藥柱實際情況。關(guān)于非線性粘彈性力學(xué)問題的研究,工程中應(yīng)用較早的半經(jīng)驗單積分方法主要有四類,分別為修正的疊加方法[15]、Beinstin彈性流體理論[16]、Schapery熱力學(xué)方法[17]及Valanis內(nèi)時理論[18]。但這些方法模型中涉及較多需要待定的參數(shù),本構(gòu)關(guān)系復(fù)雜,給工程實際應(yīng)用帶來困難。SCHAPERY提出了一個基于偽應(yīng)力和偽應(yīng)變的對應(yīng)原理,通過遺傳性積分化,將現(xiàn)時應(yīng)力和現(xiàn)時應(yīng)變轉(zhuǎn)換為偽應(yīng)力和偽應(yīng)變,給出了服從非線性彈性關(guān)系的偽應(yīng)力和偽應(yīng)變,但偽應(yīng)力和偽應(yīng)變的所有加載曲線不能完全重合,且卸載過程不能同時回到原點,不符合彈性力學(xué)性能要求[19-20]。胡強等提出基于多重積分的粘彈性本構(gòu)關(guān)系及多重Laplace變換對應(yīng)原理[21],但由于本構(gòu)關(guān)系包含太多材料參數(shù)函數(shù),很難解決實際問題。
神經(jīng)網(wǎng)絡(luò)作為在結(jié)構(gòu)分析等領(lǐng)域被廣泛應(yīng)用的力學(xué)分析方法,同樣在粘彈性研究方面也有涉及。ZHANG等進(jìn)行了基于神經(jīng)網(wǎng)絡(luò)的短碳纖維增強復(fù)合材料聚四氟乙烯(PTFE)動態(tài)力學(xué)性能預(yù)測研究[22]。陳國榮在三元件粘彈性力學(xué)模型基礎(chǔ)上,利用BP神經(jīng)網(wǎng)路方法對地下洞室圍巖粘彈性材料參數(shù)進(jìn)行反分析[23]。伍振志等提出了利用遺傳的神經(jīng)網(wǎng)絡(luò)算法進(jìn)行粘彈性巖體力學(xué)參數(shù)反演[24]。TAO等基于分?jǐn)?shù)階導(dǎo)數(shù)模型分析了粘彈性材料的粘彈性特性,通過建立BP神經(jīng)網(wǎng)絡(luò),實現(xiàn)了粘彈性參數(shù)的準(zhǔn)確識別,從而分析了粘彈性材料的動態(tài)特性[25]。ZOPF等通過引入淺層神經(jīng)網(wǎng)絡(luò),實現(xiàn)了彈性本構(gòu)關(guān)系的準(zhǔn)確描述,對于無法用模型表示的粘彈性材料,通過循環(huán)神經(jīng)網(wǎng)絡(luò)的深度學(xué)習(xí),獲得了粘彈性本構(gòu)關(guān)系的準(zhǔn)確表達(dá),以天然橡膠為例,取得了良好的結(jié)果[26]。DEVRIES等通過一類深層神經(jīng)網(wǎng)絡(luò)的深度學(xué)習(xí),得到了在任何時間及位置的流變學(xué)結(jié)構(gòu)中粘彈性問題的有效解決方案,使得神經(jīng)網(wǎng)絡(luò)計算快速、可靠及具有高空間和時間分辨率,從而實現(xiàn)了大型粘彈性地震周期模型的計算[27]。
基于此,本文通過對應(yīng)原理把粘彈性問題對應(yīng)為彈性問題,在彈性力學(xué)問題解析求解方法的基礎(chǔ)上,針對固體火箭發(fā)動機(jī)點火工作狀態(tài)下,擬提出利用對偶神經(jīng)網(wǎng)絡(luò)方法求解管形藥柱粘彈性力學(xué)問題新方法。
非線性粘彈性本構(gòu)關(guān)系在一維情況下的Volterra-Frechet展開關(guān)系式可寫為
(1)
參照Rabotnov[28],設(shè)Ji為i個不同變量的類函數(shù)乘積:
(2)
ε(t)=a1σe(t)+a2[σe(t)]2+a3[σe(t)]3+…
(3)
進(jìn)而式(3)的反函數(shù),即瞬時彈性應(yīng)力可以表示為現(xiàn)時應(yīng)變的函數(shù)為
σe(t)=φ[ε(t)]
(4)
在小變形范圍內(nèi),基于非線性本構(gòu)關(guān)系找到一種非線性粘彈性對應(yīng)原理。
平衡方程:
σij,j+Fi=0
(5)
幾何方程:
(6)
本構(gòu)關(guān)系:
(7)
或
(8)
當(dāng)邊界條件為給定位移邊界時,在邊界s1上的面力及在邊界s2上的位移分別為
(σijnj)s1=Ti=0
(9)
(ui)s2=Ui
(10)
式中nj為邊界面外法線方向余弦。
物體內(nèi)部體力Ti=0時,該粘彈性問題的解為
(11)
(12)
(13)
當(dāng)邊界條件為給定應(yīng)力邊界時,在邊界s1上的面力及在邊界s2上的位移分別為
(σijnj)s1=Ti
(14)
(ui)s2=Ui=0
(15)
則該粘彈性問題的解為
(16)
(17)
(18)
根據(jù)彈性理論可知,彈性力學(xué)平面問題可以歸結(jié)為含有邊界條件的偏微分方程求解[31]。
對于彈性力學(xué)平面問題,按應(yīng)力求解時,得到直角坐標(biāo)系下的相容方程:
(19)
式中φ為應(yīng)力函數(shù)。
假設(shè)邊界條件全部為應(yīng)力邊界條件,則在邊界s上有如下表達(dá)式:
(20)
(21)
同理,根據(jù)彈性力學(xué)理論,極坐標(biāo)系下相容方程可表示為
(22)
應(yīng)力邊界條件為
(23)
(24)
將坐標(biāo)變量x、y和應(yīng)力函數(shù)φ(x,y)的關(guān)系式寫為如下神經(jīng)網(wǎng)絡(luò)形式:
(25)
式中ki1、ki2為輸入層到隱層單元連接權(quán)值;bi為隱層單元閾值;wi為隱層到輸出層單元連接權(quán)值;c為輸出層單元閾值。
針對式(19)、式(25)對坐標(biāo)變量求偏導(dǎo)可得下列式(26)~式(28)。
(26)
(27)
(28)
文獻(xiàn)[32]指出隱層單元激活函數(shù)為exp時,較logsig、tansig等傳統(tǒng)激活函數(shù),同樣具有很好的泛化能力。為了方便構(gòu)造神經(jīng)網(wǎng)絡(luò)之間的求導(dǎo)關(guān)系,本文方法網(wǎng)絡(luò)隱層單元激活函數(shù)選取為f(x)=ex。將式(26)~式(28)代入式(19)等式左端,可得神經(jīng)網(wǎng)絡(luò)形式為
f(ki1x+ki2y+bi)=h1(x,y)
(29)
在彈性理論中,結(jié)構(gòu)除必須滿足式(29)控制方程之外,對于具體問題,還需滿足相應(yīng)的邊界條件,然后根據(jù)方程給出其神經(jīng)網(wǎng)絡(luò)形式。
不計體力,根據(jù)應(yīng)力函數(shù)φ可以得到應(yīng)力分量:
(30)
(31)
將式(30)~式(32)帶入邊界條件式(20)、式(21)中,得到對應(yīng)的神經(jīng)網(wǎng)絡(luò)形式:
式(29)、式(33)、式(34)就構(gòu)成了彈性力學(xué)問題的神經(jīng)網(wǎng)絡(luò)求解表達(dá)式。
2.2.1 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)
比較式(29)、式(33)與式(34),各網(wǎng)絡(luò)具有相同的結(jié)構(gòu),為雙輸入、單輸出、n個隱層單元的前饋型神經(jīng)網(wǎng)絡(luò),如圖1所示。
(a) Eq.(29) (b)Eq.(33) (c)Eq.(34)圖1 式(29)、式(33)、式(34)對應(yīng)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 Neural network structure of Eq.(29), Eq.(33) and Eq.(34)
Wi=wi[(ki1)4+2(ki1)2(ki2)2+(ki2)4]
W1i=wiki2(lki2+mki1)
W2i=wiki1(mki1+lki2)
運用多神經(jīng)網(wǎng)絡(luò)聯(lián)合訓(xùn)練算法對上述神經(jīng)網(wǎng)絡(luò)系統(tǒng)進(jìn)行聯(lián)合訓(xùn)練,訓(xùn)練后的網(wǎng)絡(luò)參數(shù)ki1、ki2、bi、wi帶入式(30)~式(32),得到應(yīng)力分量的神經(jīng)網(wǎng)絡(luò)數(shù)值解。
2.2.2 多神經(jīng)網(wǎng)絡(luò)聯(lián)合訓(xùn)練算法
(35)
其中,
f(ki1x+ki2y+bi)-h1(x,y)
因此,結(jié)合L-M訓(xùn)練算法,多神經(jīng)網(wǎng)絡(luò)聯(lián)合訓(xùn)練算法的計算步驟如下:
(1)初始化網(wǎng)絡(luò)參數(shù),得到訓(xùn)練誤差e、常數(shù)μ0和學(xué)習(xí)率μ的增加率β(0<β<1),令迭代步數(shù)r=1,μ=μ0。
(3)計算雅可比矩陣J(xr)、梯度矩陣g(xr)、擬海瑟矩陣H(xr)。
(4)根據(jù)L-M訓(xùn)練算法中的迭代算式Δx=-[H(xr)+μI]-1g(xr),計算權(quán)值、閾值的變化矩陣Δx。
固體火箭發(fā)動機(jī)點火工作狀態(tài)下,其藥柱的粘彈性特性對振動情況有影響,所以不能忽視。由于固體火箭發(fā)動機(jī)藥柱呈現(xiàn)出粘彈非線性,所以對其進(jìn)行非線性粘彈性力學(xué)分析是必要的。
固體火箭發(fā)動機(jī)截面可簡化為圖2的形式[33]。管形藥柱為粘彈特性的推進(jìn)劑,內(nèi)徑a=0.12 m,外徑b=0.2 m。殼體為線彈性鋼殼,厚度為δ1=0.003 m。如考慮殼體為線彈性材料,藥柱及殼體的相關(guān)參數(shù)為藥柱瞬態(tài)彈性模量E(0)=5.5766 MPa,藥柱泊松比μ=0.489,藥柱密度ρ=1635.2 kg/m3;殼體彈性模量Ec=196.5 GPa,殼體泊松比μc=0.29,殼體密度ρc=7850 kg/m3。設(shè)發(fā)動機(jī)圓筒段長L=3 mm,粘彈性材料的材料阻尼c=0.01,總體結(jié)構(gòu)瑞利阻尼阻尼系數(shù)為0.000 5。
圖2 發(fā)動機(jī)截面簡圖Fig.2 Simplified cross section of solid rocket motor
在點火狀態(tài)下,發(fā)動機(jī)工作時其圓筒段處于無約束狀態(tài)且無預(yù)應(yīng)力等影響。由于藥柱在燃燒時的強大氣流對內(nèi)壁產(chǎn)生壓力,所以分析內(nèi)壓作用下藥柱的力學(xué)行為較為貼合實際情況。根據(jù)文獻(xiàn)[34],載荷q(t)是關(guān)于時間的漸變函數(shù),關(guān)系式為
q(t)=mp0[1-exp(-nt)]
(36)
式中p0為漸變載荷幅值,p0=2.0 MPa;m、n為系數(shù),m=1,n=0.4;t為發(fā)動機(jī)工作時刻,t=1,2,3,…,30 s。
當(dāng)邊界條件為應(yīng)力邊界時,得到粘彈性應(yīng)力分量的解為式(17),現(xiàn)時應(yīng)力和瞬時應(yīng)力相等。根據(jù)文獻(xiàn)[29]中關(guān)于粘彈性材料應(yīng)力、應(yīng)變關(guān)系的闡述,粘彈性材料的瞬時應(yīng)力、應(yīng)變分量服從彈性力學(xué)的規(guī)律,上述粘彈性力學(xué)問題可以轉(zhuǎn)換為彈性力學(xué)問題。這樣就可以通過對偶神經(jīng)網(wǎng)絡(luò)對彈性力學(xué)的求解方法進(jìn)行數(shù)值求解。由于考慮殼體的影響,藥柱和殼體的材料性質(zhì)不同,不符合均勻性假設(shè),必須考慮交界面的接觸條件,由多連體的位移單值條件,取藥柱的應(yīng)力表達(dá)式為
殼體的應(yīng)力表達(dá)式為
并考慮邊界條件和接觸條件,注意到這里是平面應(yīng)變問題,藥柱和殼體的徑向位移表達(dá)式為
在接觸面上,藥柱和殼體具有相同的位移:
φ取任何數(shù)值時,上式均成立。
通過以上條件,得到待定常數(shù)A、C、A′、C′,得到ρ=b處的徑向應(yīng)力邊界條件為
(37)
其中,
根據(jù)彈性力學(xué)應(yīng)力求解相容方程可寫為
(38)
補充的邊界條件有
(τρφ)ρ=a=0,(τρφ)ρ=b=0
(39)
(σρ)ρ=a=-q(t)
(40)
將極坐標(biāo)系下應(yīng)力函數(shù)的神經(jīng)網(wǎng)絡(luò)形式帶入相容性方程及邊界條件,由于邊界條件式(39)自動滿足,所以只需構(gòu)造式(37)、式(38)及式(40)對應(yīng)的神經(jīng)網(wǎng)絡(luò)。兩個網(wǎng)絡(luò)都為單輸入、單輸出、三層前饋型神經(jīng)網(wǎng)絡(luò)。如表1所示,藥柱半徑為[0.12,0.2]范圍內(nèi)均勻取9個點作為神經(jīng)網(wǎng)絡(luò)訓(xùn)練輸入樣本,表1中y1和y2分別為相容方程式(38)和邊界條件式(37)、式(40)所對應(yīng)的網(wǎng)絡(luò)輸出。
表1 神經(jīng)網(wǎng)絡(luò)系統(tǒng)訓(xùn)練樣本Table 1 Neural network system training samples
算例中隱層單元個數(shù)n=10,運用多神經(jīng)網(wǎng)絡(luò)聯(lián)合訓(xùn)練算法訓(xùn)練網(wǎng)絡(luò)系統(tǒng)100步,t=1 s時誤差收斂曲線如圖3所示。
圖3 網(wǎng)絡(luò)系統(tǒng)誤差收斂曲線Fig.3 Error convergence curve of network systems
針對不同的時間t對應(yīng)的漸變載荷q(t),神經(jīng)網(wǎng)絡(luò)方法循環(huán)計算30次,得到t=1,2,3,…,30 s時的危險點藥柱內(nèi)表面瞬時應(yīng)力分量,如圖4所示。
(a) Instantaneous radial stress (b) Instantaneous hoop stress圖4 瞬時應(yīng)力隨時間的變化曲線Fig.4 Curves of instantaneous stress with time varies
根據(jù)Stieltjes卷積的交換律及Gurtin和Sternberg的處理方法,式(17)可以寫為
(40)
根據(jù)粘彈性理論相對松弛模量計算方法,得到相對蠕變?nèi)崃康谋磉_(dá)式為
(41)
進(jìn)而
(42)
(43)
由于粘彈性材料的瞬時應(yīng)力、應(yīng)變分量服從彈性力學(xué)規(guī)律,所以瞬時應(yīng)變分量可以寫為
(44)
(45)
根據(jù)文獻(xiàn)[35],固體火箭發(fā)動機(jī)藥柱粘彈性材料泊松比可以用自定義神經(jīng)網(wǎng)絡(luò)通過如下擬合表達(dá)式進(jìn)行擬合,且擬合效果較好。
(46)
擬合表達(dá)式中各參數(shù)如表2所示,擬合曲線如圖5所示。
表2 泊松比擬合表達(dá)式中各參數(shù)Table 2 The parameters in Poisson 's fitting expression
圖5 泊松比擬合曲線Fig.5 Poisson 's fitting curve
這樣,式(44)中的所有參量都可以得到,進(jìn)而瞬時徑向應(yīng)變就可以寫為
(47)
根據(jù)式(43)、式(47)、式(48),得到現(xiàn)時徑向應(yīng)變的表達(dá)式為
(49)
同理,也可以得到現(xiàn)時環(huán)向應(yīng)變的表達(dá)式。
根據(jù)對偶神經(jīng)網(wǎng)絡(luò)構(gòu)造及積分方法,計算現(xiàn)時應(yīng)變分量表達(dá)式中的積分項,從而得到現(xiàn)時應(yīng)變分量隨時間變化的曲線如圖6所示,對應(yīng)的現(xiàn)時應(yīng)力-應(yīng)變曲線如圖7所示。
(a) Current radial strain (b)Current hoop strain圖6 現(xiàn)時應(yīng)變隨時間的變化曲線Fig.6 Curves of current strain with time varies
(a) Current radial stress-strain (b) Current hoop stress-strain圖7 現(xiàn)時應(yīng)力-應(yīng)變曲線Fig.7 Current stress-strain curves
為了驗證方法的有效性,作為對比對圓筒進(jìn)行有限元仿真, 根據(jù)對稱性選取結(jié)構(gòu)的1/4,劃分1400個單元格,如圖8所示。
圖8 有限元網(wǎng)格劃分及觀測點Fig.8 Finite element mesh division and observation points
同樣選取式(36)漸變載荷幅值p0=2.0 MPa,m=1,n=0.4。首先,漸變載荷歷時30 s,觀察藥柱響應(yīng)隨時間的變化規(guī)律;然后,卸掉載荷再觀察藥柱的30 s回復(fù)過程。選取4項Prony級數(shù)擬合模量函數(shù),有限元仿真模擬得到1~30 s時刻位移,圖9為30 s時的徑向位移云圖。
圖9 30 s時的徑向位移云圖Fig.9 Radial displacement contour of the grain at 30 s
從圖9看出,圓筒內(nèi)徑處位移最大。由于藥柱內(nèi)表面不能產(chǎn)生過大形變是壓力載荷下藥柱的失效判據(jù),即藥柱危險點處徑向最大位移應(yīng)當(dāng)小于其極限位移,所以取內(nèi)徑處的徑向應(yīng)變作為對比。表3為本文方法與有限元方法結(jié)果對比。從表3中可看出,以劃分1400個單元格的有限元方法作為理論解,本文方法體現(xiàn)出很好的精度,最大相對誤差為 2.941 2%,符合工程實際<5%的要求。而當(dāng)有限元方法劃分20個單元格時,最大相對誤差達(dá)到7.202 2%。在計算效率方面,本文方法優(yōu)于有限元方法,有限元方法劃分1400個單元格,而本文對偶神經(jīng)網(wǎng)絡(luò)僅訓(xùn)練100步,體現(xiàn)出較好的計算效率。
表3 1~30 s危險點處徑向應(yīng)變結(jié)果對比Table 3 Comparison of radial strain results at the danger point 1~30 s
本文給出了一種基于對偶神經(jīng)網(wǎng)絡(luò)的固體火箭發(fā)動機(jī)藥柱非線性粘彈性力學(xué)分析方法,通過對發(fā)動機(jī)藥柱的簡化模型進(jìn)行粘彈性力學(xué)分析,根據(jù)對應(yīng)原理及對偶神經(jīng)網(wǎng)絡(luò)方法,得到了該算例模型各點處的現(xiàn)時應(yīng)力、應(yīng)變分量,由于粘彈性材料內(nèi)部恢復(fù)力的存在,現(xiàn)時徑向應(yīng)變分量隨時間的增加逐漸減小,現(xiàn)時環(huán)向應(yīng)變分量隨時間的增加逐漸增大,這符合粘彈性材料的規(guī)律。
由于藥柱內(nèi)表面不能產(chǎn)生過大形變是藥柱的失效判據(jù),所以作為對比,采用本文方法和有限元方法同時分析圓筒內(nèi)徑處的徑向位移,得到1~30 s時的計算結(jié)果,以劃分1400個單元格的有限元方法作為理論解,本文方法的最大相對誤差為2.941 2%,符合工程實際要求。
有限元收斂的標(biāo)準(zhǔn)為:單元越多,精度越高,越接近精確解,而計算過程中不能達(dá)到無數(shù)多個單元,所以只能越來越接近精確解。有限元軟件一般采用位移法有限元,需要進(jìn)行應(yīng)力修勻,在網(wǎng)格較少的情況下,傳統(tǒng)的應(yīng)力修勻方法在邊界處的應(yīng)力誤差很大,本文固體火箭發(fā)動機(jī)藥柱模型危險點恰處于內(nèi)邊緣處,從而使有限元方法計算受到局限。神經(jīng)網(wǎng)絡(luò)方法計算體現(xiàn)出較好的優(yōu)勢,網(wǎng)絡(luò)訓(xùn)練步數(shù)較少,而在邊界處計算精度可以保障,不受應(yīng)力修勻的影響,從而在保證計算精度的情況下,本文方法的計算效率優(yōu)于有限元方法。
本文貢獻(xiàn)在于:(1)在粘彈性力學(xué)分析中,將應(yīng)力函數(shù)寫成神經(jīng)網(wǎng)絡(luò)表達(dá)式,無需根據(jù)問題的不同假設(shè)不同的應(yīng)力函數(shù),克服了應(yīng)力函數(shù)的分析困難;(2)由于本文神經(jīng)網(wǎng)絡(luò)方法直接從理論求解粘彈性力學(xué)問題的偏微分方程出發(fā),所以相比傳統(tǒng)有限元方法的求解過程,本文方法無需對粘彈性體物理模型進(jìn)行網(wǎng)格劃分;(3)通過本文方法的求解,使粘彈性力學(xué)問題的數(shù)值解具有良好的計算精度和效率,對固體火箭發(fā)動機(jī)粘彈性體力學(xué)分析提供了一種新的思路。
同時,可以通過本文方法對模型中所有點進(jìn)行力學(xué)分析,從中找出危險點,為分析固體火箭發(fā)動機(jī)藥柱可靠性提供條件。