K8型單層球面網(wǎng)殼爆炸動(dòng)力響應(yīng)的簡(jiǎn)化計(jì)算方法研究

蘇倩倩， 翟希梅(.哈爾濱工業(yè)大學(xué) 土木工程學(xué)院，哈爾濱 50090；2.中冶建筑研究總院有限公司，北京 00088)

網(wǎng)殼結(jié)構(gòu)在爆炸荷載下受力復(fù)雜，且多作為大跨空間公共建筑，一旦因爆炸而破壞乃至倒塌，將帶來(lái)重大的生命財(cái)產(chǎn)損失。因此，有必要對(duì)該類(lèi)結(jié)構(gòu)的爆炸動(dòng)力響應(yīng)展開(kāi)相關(guān)研究，以確保網(wǎng)殼結(jié)構(gòu)在重大災(zāi)害中的安全并盡可能減少災(zāi)后損失。目前針對(duì)網(wǎng)殼結(jié)構(gòu)的動(dòng)力響應(yīng)及抗爆性能的研究多局限于有限元數(shù)值模擬，如馬加路等[1]基于網(wǎng)殼在爆炸荷載作用下的破壞及倒塌情況，提出了安全距離的計(jì)算方法。Zhai等[2-3]分析了不同炸藥點(diǎn)位置、支撐條件以及屋面板厚度等參數(shù)下網(wǎng)殼在偏心爆炸荷載下的動(dòng)力響應(yīng)，并得到了結(jié)構(gòu)的四種響應(yīng)模式。高軒能等[4]研究了不同參數(shù)對(duì)單層球面鋼網(wǎng)殼在內(nèi)爆作用下動(dòng)力響應(yīng)的影響規(guī)律，并對(duì)柱殼結(jié)構(gòu)的壓力場(chǎng)進(jìn)行了分析。田力等[5]研究了地下隧道內(nèi)爆炸沖擊作用下雙層柱面網(wǎng)殼的動(dòng)力響應(yīng)。可見(jiàn)，眾多學(xué)者都對(duì)爆炸荷載作用下網(wǎng)殼的動(dòng)力響應(yīng)以及防爆泄爆機(jī)理展開(kāi)較多研究，但相關(guān)網(wǎng)殼的試驗(yàn)和理論研究尚處于空白。

鑒于大跨網(wǎng)殼結(jié)構(gòu)復(fù)雜、體型龐大，其數(shù)值計(jì)算耗時(shí)巨大，且國(guó)內(nèi)外缺少網(wǎng)殼結(jié)構(gòu)的爆炸試驗(yàn)結(jié)果可對(duì)目前的爆炸數(shù)值響應(yīng)模擬結(jié)果進(jìn)行驗(yàn)證。因此，本文希望在網(wǎng)殼結(jié)構(gòu)爆炸響應(yīng)計(jì)算方法上探索新的方向，并為結(jié)構(gòu)爆炸動(dòng)力響應(yīng)分析提供有效理論依據(jù)。根據(jù)哈密頓原理，以拉格朗日方程為基礎(chǔ)[6-7]，本文提出了一種適用于K8型單層球面網(wǎng)殼在爆轟荷載作用下動(dòng)力響應(yīng)求解的理論分析模型，并通過(guò)MATLAB求解微分方程組，進(jìn)而得到結(jié)構(gòu)動(dòng)力響應(yīng)。最后，將分別通過(guò)爆炸荷載下的簡(jiǎn)支梁及K8型單層球面網(wǎng)殼的響應(yīng)結(jié)果與有限元結(jié)果的比較對(duì)本文方法的適應(yīng)性進(jìn)行驗(yàn)證。

1 簡(jiǎn)支梁爆炸動(dòng)力響應(yīng)的理論分析模型

本文采用三角形荷載來(lái)近似考慮爆炸荷載。首先以簡(jiǎn)支梁作為分析對(duì)象，計(jì)算結(jié)構(gòu)在豎向均布爆炸荷載作用下的動(dòng)態(tài)響應(yīng)。簡(jiǎn)支梁模型如圖1所示，跨度為2 m，截面為200 mm×200 mm的矩形。模型選用彈性鋼材，密度為7 850 kg/m3，彈性模量為2.06×1011Pa，在簡(jiǎn)支梁上施加垂直于梁向下的均布線荷載p，p隨時(shí)間的變化如圖2所示。

圖1 簡(jiǎn)支梁模型Fig.1 Simplysupportedbeam圖2 p-t曲線Fig.2 p-tcurve

通過(guò)求解式(1)所示的拉格朗日平衡方程，可獲得簡(jiǎn)支梁在三角形荷載作用下的動(dòng)力響應(yīng)

(1)

其中：V=∑U+PE=Um+Ub+PE

式中:Ci(t)為廣義位移;T為動(dòng)能;V為總勢(shì)能;Um為軸向變形能;Ub為彎曲變形能;PE為勢(shì)能損失量。

簡(jiǎn)支梁的變形方程取正弦函數(shù)的疊加，為了不讓計(jì)算過(guò)于復(fù)雜，取兩個(gè)正弦函數(shù)進(jìn)行疊加，如式(2)所示。變形函數(shù)的選取，要滿足邊界條件，w(0,t)=0，w(l,t)=0。

(2)

式中:w(l,t)指節(jié)點(diǎn)的豎向變形量;L指梁長(zhǎng);l指梁長(zhǎng)的變量，從0到L代表從簡(jiǎn)支梁的左端到右端。

根據(jù)變形函數(shù)，可以得到簡(jiǎn)支梁每個(gè)點(diǎn)的豎向變形，根據(jù)式(3)計(jì)算出梁的彎曲變形能

(3)

式中:E為彈性模量;I為桿件截面慣性矩。由于簡(jiǎn)支梁在豎向荷載作用下沒(méi)有軸向變形，所以軸向變形能為零，簡(jiǎn)支梁的內(nèi)能就是彎曲變形能。

根據(jù)變形函數(shù)可以得到簡(jiǎn)支梁上每點(diǎn)的速度如式(4)所示

(4)

(5)

式中:ρ為密度;A為截面面積。

勢(shì)能損失量PE可根據(jù)式(6)進(jìn)行計(jì)算，其中，荷載p為均布線荷載，垂直向下作用于梁

(6)

計(jì)算出簡(jiǎn)支梁的內(nèi)能、動(dòng)能和勢(shì)能損失量，代入式(1)所示的拉格朗日方程，并通過(guò)MATLAB編程求解該微分方程組的數(shù)值解。令l=L/2，可以得到簡(jiǎn)支梁中點(diǎn)的位移響應(yīng)，進(jìn)而計(jì)算得到簡(jiǎn)支梁內(nèi)能和動(dòng)能時(shí)程。

下面，對(duì)上述簡(jiǎn)支梁模型進(jìn)行有限元分析，采用LS-DYNA動(dòng)力分析軟件建立相應(yīng)有限元模型，簡(jiǎn)支梁采用BEAM161單元，均布線荷載p垂直于梁向下。

經(jīng)有限元分析，獲得各個(gè)峰值位移時(shí)刻簡(jiǎn)支梁的豎向變形如圖3所示(變形放大10倍)。將有限元模擬與基于拉格朗日動(dòng)力方程的理論分析模型MATLAB計(jì)算得到的簡(jiǎn)支梁中點(diǎn)位移時(shí)程、簡(jiǎn)支梁系統(tǒng)動(dòng)能和總內(nèi)能(彎曲變形能)結(jié)果對(duì)比如圖4和表1所示。

(a) t=0.041 s

(b) t=0.086 s

(c) t=0.127 s

(d) t=0.172 s

(e) t=0.214 s

(f) t=0.257 s

Fig.3 Deformation of simply supported beam (10 times amplified)

從圖4和表1可以看出，有限元分析與基于拉格朗日方法理論分析模型計(jì)算得到的簡(jiǎn)支梁中點(diǎn)振動(dòng)規(guī)律、峰值位移以及系統(tǒng)的動(dòng)能、總內(nèi)能都極為接近。簡(jiǎn)支梁中點(diǎn)的位移響應(yīng)按照正弦形式振蕩，在0.04 s左右達(dá)到峰值位移，之后節(jié)點(diǎn)位移衰減到負(fù)值再回到正值，峰值位移一直在減小，0.2 s后由于沒(méi)有荷載作用，節(jié)點(diǎn)在平衡位置附近上下振動(dòng)。由于未考慮阻尼，0.2 s后沒(méi)有衰減，以固定的峰值位移一直振動(dòng)下去，且簡(jiǎn)支梁系統(tǒng)的動(dòng)能和內(nèi)能峰值也都不再發(fā)生變化。經(jīng)計(jì)算，簡(jiǎn)支梁中點(diǎn)正向位移峰值誤差僅為0.37%，可見(jiàn)采用該理論分析模型來(lái)近似計(jì)算結(jié)構(gòu)在爆炸荷載作用下的動(dòng)力響應(yīng)是可行的。

(a) 中點(diǎn)位移時(shí)程曲線

(b) 動(dòng)能時(shí)程曲線

(c) 內(nèi)能時(shí)程曲線

動(dòng)力響應(yīng)正向峰值/m與有限元誤差值/%負(fù)向峰值/m與有限元誤差值/%中點(diǎn)位移有限元模擬0.0136400.007210理論計(jì)算0.01359-0.370.00717-0.55動(dòng)能有限元模擬0.5563000理論計(jì)算0.5563000內(nèi)能有限元模擬1.5603000理論計(jì)算1.56120.0600

2 K8型網(wǎng)殼基于拉格朗日方法動(dòng)力響應(yīng)計(jì)算

下面，我們將基于拉格朗日方法的理論分析模型應(yīng)用到K8型單層球面網(wǎng)殼結(jié)構(gòu)，由于網(wǎng)殼結(jié)構(gòu)比較復(fù)雜，桿件較多，本文暫時(shí)僅考慮網(wǎng)殼桿件，并選取分頻數(shù)為2的K8型單層球面網(wǎng)殼進(jìn)行計(jì)算分析。根據(jù)變形函數(shù)中廣義位移參數(shù)的個(gè)數(shù)，本節(jié)分為兩參數(shù)分析、三參數(shù)分析以及改進(jìn)后的三參數(shù)分析。

2.1 爆炸荷載作用下K8型網(wǎng)殼動(dòng)力響應(yīng)兩參數(shù)分析模型

模型為分頻數(shù)為2的K8型單層球面網(wǎng)殼桿件，跨度為5 m、矢跨比為1/5，主桿采用Φ14×0.5 mm圓鋼管，網(wǎng)殼底部節(jié)點(diǎn)采用固定約束，模型如圖5所示。模型選用彈性鋼材材料，參數(shù)設(shè)置同第1節(jié)。

在網(wǎng)殼的每個(gè)頂點(diǎn)上施加集中荷載P(t)，P(t)為簡(jiǎn)化爆炸荷載——三角形荷載，P(t)隨時(shí)間的變化如圖6所示，在t=0時(shí)刻P0為荷載峰值，荷載P(t)的方向沿徑向向外，荷載施加方式見(jiàn)圖5。

(a) 軸側(cè)圖

(b) 正視圖

由于網(wǎng)殼模型以及三角形荷載的對(duì)稱性，本文取1/8模型進(jìn)行理論分析計(jì)算，這樣可以很大程度的減少計(jì)算量和編程量，理論分析模型如圖7所示。圖中，1～6表示網(wǎng)殼節(jié)點(diǎn)，L1～L6表示網(wǎng)殼桿件。

從網(wǎng)殼結(jié)構(gòu)整體變形方面假設(shè)網(wǎng)殼的變形函數(shù)，如式(7)所示。通過(guò)變形函數(shù)可以確定6個(gè)節(jié)點(diǎn)的變形。由于該變形函數(shù)中有兩個(gè)待求解廣義位移參數(shù)，因此本節(jié)為兩參數(shù)分析。變形函數(shù)滿足邊界條件，w(0,t)=0

(7)

式中:w(l,t)指節(jié)點(diǎn)的徑向變形量;L指弧長(zhǎng);l指弧長(zhǎng)的變量，從0到L，代表從網(wǎng)殼邊緣到頂點(diǎn)。

圖6 P-t曲線Fig.6 P-tcurve圖7 理論分析模型Fig.7 Theoreticalanalysismodel

根據(jù)變形函數(shù)，求得網(wǎng)殼6個(gè)節(jié)點(diǎn)的徑向變形量，然后計(jì)算出9根桿件變形后的長(zhǎng)度以及變形量Δi。繼而根據(jù)式(8)可以求得9根桿件的軸向變形能之和，Um是關(guān)于廣義位移Ci(t)的函數(shù)

(8)

式中:Lbi指每根桿件變形前的長(zhǎng)度。網(wǎng)殼6個(gè)節(jié)點(diǎn)相應(yīng)的速度可由式(9)求得

(9)

(10)

求解過(guò)程中，對(duì)稱截面位置處桿件的截面取Φ7×0.5 mm圓鋼管，其他桿件截面為Φ14×0.5 mm圓鋼管。

勢(shì)能損失量PE可根據(jù)式(11)進(jìn)行計(jì)算，由于模型的對(duì)稱性，作用在1點(diǎn)的荷載為P(t)/8，作用在2點(diǎn)和3點(diǎn)的荷載為P(t)/2。

PE=-∑P(t)×w(l,t)

(11)

計(jì)算出網(wǎng)殼的內(nèi)能、動(dòng)能和勢(shì)能損失量，代入式(1)所示的拉格朗日方程，通過(guò)MATLAB編程求解，即可得到結(jié)構(gòu)6個(gè)節(jié)點(diǎn)的位移響應(yīng)，進(jìn)而計(jì)算得到網(wǎng)殼結(jié)構(gòu)的內(nèi)能和動(dòng)能時(shí)程。

2.2 爆炸荷載作用下K8型網(wǎng)殼動(dòng)力響應(yīng)三參數(shù)分析模型

本節(jié)在2.1節(jié)中網(wǎng)殼整體變形函數(shù)(式(7))的基礎(chǔ)上，再假設(shè)每根桿件的變形如式(12)所示。這樣，本節(jié)分析中，變形函數(shù)中共有三個(gè)待求解廣義位移參數(shù)，稱為三參數(shù)分析。

(12)

每根桿件變形后的長(zhǎng)度由兩部分組成，第一部分是由于節(jié)點(diǎn)徑向位移w(l,t)導(dǎo)致的桿件變形量，該部分就是2.1節(jié)中桿件變形后的長(zhǎng)度，第二部分是由于桿件彎曲導(dǎo)致的變形量，由變形函數(shù)wb(x,t)以及以直代曲的近似計(jì)算方法，可以得到每根桿件的變形量如式(13)所示。繼而根據(jù)式(8)可以得到9根桿件的軸向變形能之和。桿件的彎曲變形能由式(14)計(jì)算。

(13)

(14)

網(wǎng)殼6個(gè)節(jié)點(diǎn)相應(yīng)的速度如式(9)所示，網(wǎng)殼桿件上任一點(diǎn)的速度由兩部分組成，第一部分是由于網(wǎng)殼節(jié)點(diǎn)運(yùn)動(dòng)導(dǎo)致的桿件上線性分布的速度，這一部分同2.1節(jié)，第二部分是由于桿件發(fā)生彎曲變形導(dǎo)致的速度，可以由式(15)計(jì)算得到，將這兩部分速度疊加，即可得到網(wǎng)殼桿件上任意一點(diǎn)的速度。然后根據(jù)式(10)可以計(jì)算9根桿件的動(dòng)能之和。

(15)

勢(shì)能損失量PE同2.1節(jié)一樣，將網(wǎng)殼的內(nèi)能、動(dòng)能和勢(shì)能損失量代入式(1)所示的拉格朗日方程，并通過(guò)MATLAB編程求解，即可得到結(jié)構(gòu)6個(gè)節(jié)點(diǎn)的位移響應(yīng)，進(jìn)而計(jì)算得到網(wǎng)殼結(jié)構(gòu)的內(nèi)能和動(dòng)能時(shí)程。

2.3 K8型網(wǎng)殼理論計(jì)算與有限元結(jié)果的對(duì)比分析

為了與基于拉格朗日動(dòng)力方程的理論分析結(jié)果對(duì)比，本節(jié)采用LS-DYNA軟件建立了一個(gè)與理論計(jì)算模型相對(duì)應(yīng)的K8型單層球面網(wǎng)殼1/4模型，主桿采用Φ14×0.5 mm圓鋼管，對(duì)稱截面處采用Φ7×0.5 mm圓鋼管，采用BEAM161單元。有限元模型如圖8所示。

(a) 俯視圖

(b) 正視圖

鋼材選用彈性材料，參數(shù)設(shè)置同第1節(jié)。在網(wǎng)殼頂點(diǎn)及第一環(huán)的三個(gè)頂點(diǎn)上施加集中荷載，由于模型的對(duì)稱性，作用在1點(diǎn)的荷載為P(t)/4，作用在2點(diǎn)的荷載為P(t)，作用在22點(diǎn)和42點(diǎn)的荷載為P(t)/2。P(t)隨時(shí)間的變化同2.1節(jié)。

約束模型底部5個(gè)節(jié)點(diǎn)的三向位移、對(duì)節(jié)點(diǎn)1和22施加Y軸對(duì)稱約束、對(duì)節(jié)點(diǎn)1和42施加X(jué)軸對(duì)稱約束，不同荷載下的響應(yīng)結(jié)果皆顯示：X向和Y向位移較小且數(shù)值接近，兩者的位移峰值響應(yīng)僅為Z向位移響應(yīng)的28.36%，表明位移響應(yīng)是以Z向?yàn)橹鳎虼吮疚暮罄m(xù)分析中僅給出Z向位移響應(yīng)結(jié)果。當(dāng)荷載峰值P0為50 kN時(shí)，提取K8型球面網(wǎng)殼結(jié)構(gòu)在不同時(shí)刻的Z向位移如圖9所示。

從圖9可以看出，三角形爆炸荷載作用下網(wǎng)殼結(jié)構(gòu)的位移響應(yīng)受高階振型影響非常大，這是由網(wǎng)殼自身結(jié)構(gòu)特性決定的。同時(shí)，網(wǎng)殼桿件上節(jié)點(diǎn)相較于桿件與桿件的交點(diǎn)來(lái)說(shuō)，位移響應(yīng)有滯后的效應(yīng)。網(wǎng)殼底部桿件發(fā)生響應(yīng)的時(shí)間最遲，且位移響應(yīng)較小。

荷載峰值P0分別為10 kN、20 kN、30 kN、40 kN、50 kN時(shí)，將有限元數(shù)值分析、基于拉格朗日方程的理論兩參數(shù)和三參數(shù)分析得到的網(wǎng)殼最大節(jié)點(diǎn)位移(Z向)時(shí)程對(duì)比如圖10所示。

(a) t=0.000 25 s

(b) t=0.000 5 s

(c) t=0.000 75 s

(d) t=0.001 s

(e) t=0.001 5 s

(f) t=0.002 s

Fig.9Z-displacement of K8 single-layer reticulated shell(P0=50 kN)

(a) P0=10 kN

(b) P0=20 kN

(c) P0=30 kN

(d) P0=40 kN

(e) P0=50 kN

從圖10可以看出，兩參數(shù)和三參數(shù)理論分析模型得到的最大位移節(jié)點(diǎn)按照正弦形式振蕩，節(jié)點(diǎn)都在0.001 5 s左右達(dá)到正向峰值位移，然后開(kāi)始衰減，經(jīng)過(guò)平衡位置后負(fù)向位移逐漸增大，并且三參數(shù)分析中節(jié)點(diǎn)再次回到平衡位置的時(shí)間要比兩參數(shù)分析滯后，荷載峰值P0越大，這種規(guī)律越明顯。當(dāng)荷載峰值P0從10 kN增大到50 kN時(shí)，有限元數(shù)值模擬與基于拉格朗日動(dòng)力方程計(jì)算得到的網(wǎng)殼最大位移節(jié)點(diǎn)時(shí)程(Z向)的正向峰值位移對(duì)應(yīng)較好，但理論分析達(dá)到正向峰值位移的時(shí)刻比有限元分析滯后一倍左右。同時(shí)，基于拉格朗日動(dòng)力方程兩參數(shù)計(jì)算得到的位移響應(yīng)比三參數(shù)計(jì)算結(jié)果偏大，這是因?yàn)閮蓞?shù)計(jì)算過(guò)程中沒(méi)有考慮桿件的彎曲，外力功僅轉(zhuǎn)變?yōu)檩S向變形能和動(dòng)能，而三參數(shù)計(jì)算過(guò)程中外力功轉(zhuǎn)變?yōu)檩S向變形能、彎曲變形能和動(dòng)能。

對(duì)比前1-2個(gè)振動(dòng)周期，基于拉格朗日動(dòng)力方程的理論分析模型和有限元產(chǎn)生誤差的原因有：①由于網(wǎng)殼結(jié)構(gòu)自身特性，高階振型影響非常大，而理論分析模型中假設(shè)的變形函數(shù)比較簡(jiǎn)單，最多考慮前三階振型，且該變形函數(shù)和有限元的變形模式有偏差，不能完全相符；②有限元分析中，三角形爆炸荷載作用下網(wǎng)殼桿件上節(jié)點(diǎn)相較于桿件與桿件的交點(diǎn)來(lái)說(shuō)，位移響應(yīng)有滯后的效應(yīng)，而理論分析模型中假設(shè)變形函數(shù)時(shí)，是以桿件交點(diǎn)與桿件同時(shí)達(dá)到峰值位移進(jìn)行計(jì)算的；③MATLAB求解拉格朗日方程過(guò)程中，采用數(shù)值解，有累計(jì)的數(shù)值求解誤差。

其中，第一個(gè)和第二個(gè)原因最關(guān)鍵，第二個(gè)原因也是導(dǎo)致理論分析模型兩參數(shù)和三參數(shù)分析結(jié)果相差較小的主要原因?？梢酝ㄟ^(guò)增加理論分析模型變形函數(shù)中的參數(shù)個(gè)數(shù)以及單獨(dú)假設(shè)每根桿件的變形函數(shù)來(lái)改善，但這樣將導(dǎo)致MATLAB編程量大幅度增加。

2.4 爆炸荷載下K8型網(wǎng)殼動(dòng)力響應(yīng)計(jì)算方法的修正

根據(jù)2.3節(jié)中基于拉格朗日動(dòng)力方程理論分析與有限元分析產(chǎn)生誤差的原因，對(duì)2.2節(jié)中的理論分析模型進(jìn)行改進(jìn)，改進(jìn)方法介紹如下。

有限元模型中，網(wǎng)殼結(jié)構(gòu)的最大位移節(jié)點(diǎn)(2點(diǎn))幾乎都在0.000 7 s左右達(dá)到峰值，當(dāng)荷載峰值P0為30 kN時(shí)，在0.000 738 s達(dá)到最大位移0.020 3 m，針對(duì)該算例，提取t=0.000 738 s時(shí)刻節(jié)點(diǎn)1和節(jié)點(diǎn)2之間桿件節(jié)點(diǎn)、節(jié)點(diǎn)2和節(jié)點(diǎn)101之間桿件節(jié)點(diǎn)的Z向位移如圖11所示。為了更直觀的將變形曲線表達(dá)清晰，繪制時(shí)將位移放大10倍。

圖11 桿件變形圖

從圖11可以看出，上部桿件的變形接近1/2周期的三角函數(shù)，下部桿件的變形接近1/4周期的三角函數(shù)，于是，在網(wǎng)殼桿件交點(diǎn)變形函數(shù)(式(7))的基礎(chǔ)上，對(duì)2.2節(jié)中桿件的變形函數(shù)進(jìn)行修正。其中，理論分析模型中桿件編號(hào)見(jiàn)圖7所示。

桿件1、桿件2和桿件5的變形函數(shù)取式(16)，桿件3、桿件6、桿件7和桿件4的變形函數(shù)取式(17)。

(16)

(17)

根據(jù)變形函數(shù)，先由2.1節(jié)中整體變形函數(shù)得到網(wǎng)殼節(jié)點(diǎn)變形后的坐標(biāo)以及桿件變形后長(zhǎng)度，再以直代曲的近似計(jì)算方法，桿件1、桿件2和桿件5的最終變形量由式(18)計(jì)算，桿件3、桿件6、桿件7和桿件4的最終變形量由式(19)計(jì)算，進(jìn)而根據(jù)式(8)計(jì)算所有桿件軸向變形能之和。

(18)

(19)

桿件1、桿件2和桿件5的彎曲變形能由式(20)計(jì)算，桿件3、桿件6、桿件7和桿件4的彎曲變形能由式(21)計(jì)算，桿件8和桿件9彎曲變形能為零。

(20)

(21)

同2.2節(jié)一樣，網(wǎng)殼桿件上任一點(diǎn)的速度由兩部分組成，其中，桿件1、桿件2和桿件5任意一點(diǎn)由彎曲導(dǎo)致的速度由式(22)計(jì)算，桿件3、桿件6、桿件7和桿件4任意一點(diǎn)由彎曲導(dǎo)致的速度由式(23)計(jì)算。然后根據(jù)式(10)可以計(jì)算9根桿件的動(dòng)能之和。

(22)

(23)

勢(shì)能損失量PE同2.1節(jié)一致，計(jì)算出網(wǎng)殼的內(nèi)能、動(dòng)能和勢(shì)能損失量，代入式(1)所示的拉格朗日方程，并通過(guò)MATLAB編程求解，即可得到結(jié)構(gòu)6個(gè)節(jié)點(diǎn)的位移響應(yīng)，進(jìn)而計(jì)算得到結(jié)構(gòu)的內(nèi)能和動(dòng)能時(shí)程。

當(dāng)荷載峰值P0分別為10 kN、20 kN、30 kN、40 kN、50 kN時(shí)，將有限元分析、基于拉格朗日方程的理論三參數(shù)分析以及改進(jìn)后的三參數(shù)分析得到的網(wǎng)殼最大節(jié)點(diǎn)位移(Z向)時(shí)程對(duì)比如圖12所示，數(shù)值見(jiàn)表2。

從圖12可以看出，改進(jìn)后的三參數(shù)分析得到的網(wǎng)殼的最大節(jié)點(diǎn)位移(Z向)時(shí)程與有限元分析更為接近。并且，荷載峰值P0越大，誤差越小，尤其是負(fù)向峰值位移和達(dá)到負(fù)向峰值位移的時(shí)間與有限元結(jié)果更為貼近。

(a) P0=10 kN

(b) P0=20 kN

(c) P0=30 kN

(d) P0=40 kN

(e) P0=50 kN

從表2可以看出，當(dāng)荷載峰值P0從10 kN增大到30 kN時(shí)，原三參數(shù)分析得到的正向峰值位移值均比有限元分析結(jié)果要小，當(dāng)P0為40 kN和50 kN時(shí)，則反之。而改進(jìn)后的三參數(shù)分析結(jié)果均比有限元分析結(jié)果偏大，最大誤差為26.8%。除了當(dāng)P0為10 kN時(shí)，負(fù)向峰值位移誤差較大，其他4組荷載下，負(fù)向峰值位移都吻合較好，原三參數(shù)分析最大誤差為35.9%，改進(jìn)后的三參數(shù)分析最大誤差為16.1%。

綜合圖12和表2，本文認(rèn)為改進(jìn)后的三參數(shù)理論分析模型得到的網(wǎng)殼最大節(jié)點(diǎn)位移(Z向)振動(dòng)規(guī)律比原三參數(shù)分析更貼合有限元分析結(jié)果。

將有限元分析、基于拉格朗日方程的理論三參數(shù)分析以及改進(jìn)后的三參數(shù)分析得到的網(wǎng)殼系統(tǒng)動(dòng)能、內(nèi)能時(shí)程對(duì)比分別如圖13和圖14所示。

從圖13可以看出，改進(jìn)后的三參數(shù)分析得到的網(wǎng)殼結(jié)構(gòu)的動(dòng)能時(shí)程以及達(dá)到動(dòng)能峰值的時(shí)刻與有限元結(jié)果的振動(dòng)規(guī)律更為接近，理論分析得到的總動(dòng)能最大值均出現(xiàn)在時(shí)程的第二個(gè)峰值，有限元的總動(dòng)能最大值出現(xiàn)在時(shí)程的第二個(gè)或第三個(gè)峰值。

表2 網(wǎng)殼位移響應(yīng)對(duì)比

(a)P0=10kN(b)P0=20kN(c)P0=30kN(d)P0=40kN

(e) P0=50 kN

從圖14可以看出，改進(jìn)后的三參數(shù)理論分析模型得到的內(nèi)能峰值、達(dá)到內(nèi)能峰值所用的時(shí)間均比原三參數(shù)分析結(jié)果更接近有限元分析，當(dāng)荷載峰值P0為30 kN和40 kN時(shí)，改進(jìn)后的三參數(shù)理論分析得到的內(nèi)能最大值與有限元分析結(jié)果誤差最小。

(a) P0=10 kN

(b) P0=20 kN

(c) P0=30 kN

(d) P0=40 kN

(e) P0=50 kN

綜合考慮最大節(jié)點(diǎn)位移響應(yīng)、網(wǎng)殼結(jié)構(gòu)的動(dòng)能和內(nèi)能情況，建議采用基于拉格朗日動(dòng)力方程的理論分析模型計(jì)算網(wǎng)殼在爆炸荷載作用下的動(dòng)力響應(yīng)時(shí)，荷載峰值P0不小于20 kN。

3 結(jié) 論

本文以拉格朗日動(dòng)力方程為基礎(chǔ)，提出一種適用于簡(jiǎn)支梁和K8型單層球面網(wǎng)殼在簡(jiǎn)化爆炸荷載作用下動(dòng)力響應(yīng)的計(jì)算方法，并與LS-DYNA有限元分析結(jié)果進(jìn)行對(duì)比和誤差分析，得到以下結(jié)論：

(1) 利用該方法獲得的簡(jiǎn)支梁在三角形爆炸荷載作用下的跨中振動(dòng)規(guī)律、峰值位移、速度以及系統(tǒng)的動(dòng)能和內(nèi)能等與有限元結(jié)果十分吻合，其中最大節(jié)點(diǎn)位移峰值誤差僅為0.37%。

(2) 當(dāng)荷載峰值P0從10 kN增大到50 kN時(shí)，有限元數(shù)值模擬與基于拉格朗日動(dòng)力方程計(jì)算得到的網(wǎng)殼最大位移節(jié)點(diǎn)時(shí)程(Z向)的正向峰值位移對(duì)應(yīng)較好，但理論分析達(dá)到正向峰值位移的時(shí)刻比有限元分析滯后一倍左右。

(3) 改進(jìn)后三參數(shù)分析模型得到的最大節(jié)點(diǎn)位移(Z向)時(shí)程、網(wǎng)殼結(jié)構(gòu)的總動(dòng)能和總內(nèi)能與有限元分析結(jié)果更為接近。建議如采用該理論分析模型計(jì)算網(wǎng)殼動(dòng)力響應(yīng)時(shí)，荷載峰值P0應(yīng)不小于20 kN，且假設(shè)的變形函數(shù)應(yīng)接近網(wǎng)殼實(shí)際變形。

[1] 馬加路，支旭東，范峰,等. 外爆荷載下Kiewitt8型單層球面網(wǎng)殼的動(dòng)力響應(yīng)[J]. 振動(dòng)與沖擊, 2015,34(21):65-70.

MA Jialu, ZHI Xudong, FAN Feng, et al. Dynamic responses of Kiewitt 8 single-layer reticulated domes subjected to outside explosion loads[J]. Journal of Vibration and Shock, 2015,34(21):65-70.

[2] ZHAI Ximei, WANG Yonghui, HUANG Ming. Performance and protection approach of single-layer reticulated dome subjected to blast loading[J]. Thin-Walled Structures, 2013,73:57-67.

[3] 翟希梅，黃明. 偏心爆炸荷載下網(wǎng)殼結(jié)構(gòu)的動(dòng)力響應(yīng)分析[J]. 振動(dòng)與沖擊. 2014,33(4):178-184.

ZHAI Ximei, HUANG Ming. Dynamic response of a reticulated shell under eccentric blast loading[J]. Journal of Vibration and Shock, 2014,33(4): 178-184.

[4] 高軒能，李超，江媛. 單層球面鋼網(wǎng)殼結(jié)構(gòu)在內(nèi)爆炸作用下的動(dòng)力響應(yīng)[J]. 天津大學(xué)學(xué)報(bào)(自然科學(xué)與工程技術(shù)版), 2015, 48(增刊1):102-109.

GAO Xuanneng, LI Chao, JIANG Yuan. Dynamic responses of single-layer spherical steel reticulated shell under internal explosions[J]. Journal of Tianjin University (Science and Technology), 2015,48 (Suppl):102-109.

[5] 田力，李靈聚. 雙層柱面網(wǎng)殼結(jié)構(gòu)在地下隧道內(nèi)爆炸沖擊下的動(dòng)力響應(yīng)分析[J]. 振動(dòng)工程學(xué)報(bào), 2011, 24(5):482-490.

TIAN Li, LI Lingju. Dynamic response analysis of the cylindrical reticulated shell structure due to an underground in-tunnel explosion[J]. Journal of Vibration Engineering, 2011, 24(5):482-490.

[6] SCHLEYER G K, HSU S S. A modelling scheme for predicting the response of elastic-plastic structures to pulse pressure loading[J]. International Journal of Impact Engineering, 2000(24):759-777.

[7] DONALDSON B K. Introduction to structural dynamics[M]. Cambridge: Cambridge University Press, 2012: 46-70.