K8型單層球面網殼爆炸動力響應的簡化計算方法研究

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

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

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

1 簡支梁爆炸動力響應的理論分析模型

本文采用三角形荷載來近似考慮爆炸荷載。首先以簡支梁作為分析對象，計算結構在豎向均布爆炸荷載作用下的動態響應。簡支梁模型如圖1所示，跨度為2 m，截面為200 mm×200 mm的矩形。模型選用彈性鋼材，密度為7 850 kg/m3，彈性模量為2.06×1011Pa，在簡支梁上施加垂直于梁向下的均布線荷載p，p隨時間的變化如圖2所示。

圖1 簡支梁模型Fig.1 Simplysupportedbeam圖2 p-t曲線Fig.2 p-tcurve

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

(1)

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

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

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

(2)

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

根據變形函數，可以得到簡支梁每個點的豎向變形，根據式(3)計算出梁的彎曲變形能

(3)

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

根據變形函數可以得到簡支梁上每點的速度如式(4)所示

(4)

(5)

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

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

(6)

計算出簡支梁的內能、動能和勢能損失量，代入式(1)所示的拉格朗日方程，并通過MATLAB編程求解該微分方程組的數值解。令l=L/2，可以得到簡支梁中點的位移響應，進而計算得到簡支梁內能和動能時程。

下面，對上述簡支梁模型進行有限元分析，采用LS-DYNA動力分析軟件建立相應有限元模型，簡支梁采用BEAM161單元，均布線荷載p垂直于梁向下。

經有限元分析，獲得各個峰值位移時刻簡支梁的豎向變形如圖3所示(變形放大10倍)。將有限元模擬與基于拉格朗日動力方程的理論分析模型MATLAB計算得到的簡支梁中點位移時程、簡支梁系統動能和總內能(彎曲變形能)結果對比如圖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可以看出，有限元分析與基于拉格朗日方法理論分析模型計算得到的簡支梁中點振動規律、峰值位移以及系統的動能、總內能都極為接近。簡支梁中點的位移響應按照正弦形式振蕩，在0.04 s左右達到峰值位移，之后節點位移衰減到負值再回到正值，峰值位移一直在減小，0.2 s后由于沒有荷載作用，節點在平衡位置附近上下振動。由于未考慮阻尼，0.2 s后沒有衰減，以固定的峰值位移一直振動下去，且簡支梁系統的動能和內能峰值也都不再發生變化。經計算，簡支梁中點正向位移峰值誤差僅為0.37%，可見采用該理論分析模型來近似計算結構在爆炸荷載作用下的動力響應是可行的。

(a) 中點位移時程曲線

(b) 動能時程曲線

(c) 內能時程曲線

動力響應正向峰值/m與有限元誤差值/%負向峰值/m與有限元誤差值/%中點位移有限元模擬0.0136400.007210理論計算0.01359-0.370.00717-0.55動能有限元模擬0.5563000理論計算0.5563000內能有限元模擬1.5603000理論計算1.56120.0600

2 K8型網殼基于拉格朗日方法動力響應計算

下面，我們將基于拉格朗日方法的理論分析模型應用到K8型單層球面網殼結構，由于網殼結構比較復雜，桿件較多，本文暫時僅考慮網殼桿件，并選取分頻數為2的K8型單層球面網殼進行計算分析。根據變形函數中廣義位移參數的個數，本節分為兩參數分析、三參數分析以及改進后的三參數分析。

2.1 爆炸荷載作用下K8型網殼動力響應兩參數分析模型

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

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

(a) 軸側圖

(b) 正視圖

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

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

(7)

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

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

根據變形函數，求得網殼6個節點的徑向變形量，然后計算出9根桿件變形后的長度以及變形量Δi。繼而根據式(8)可以求得9根桿件的軸向變形能之和，Um是關于廣義位移Ci(t)的函數

(8)

式中:Lbi指每根桿件變形前的長度。網殼6個節點相應的速度可由式(9)求得

(9)

(10)

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

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

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

(11)

計算出網殼的內能、動能和勢能損失量，代入式(1)所示的拉格朗日方程，通過MATLAB編程求解，即可得到結構6個節點的位移響應，進而計算得到網殼結構的內能和動能時程。

2.2 爆炸荷載作用下K8型網殼動力響應三參數分析模型

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

(12)

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

(13)

(14)

網殼6個節點相應的速度如式(9)所示，網殼桿件上任一點的速度由兩部分組成，第一部分是由于網殼節點運動導致的桿件上線性分布的速度，這一部分同2.1節，第二部分是由于桿件發生彎曲變形導致的速度，可以由式(15)計算得到，將這兩部分速度疊加，即可得到網殼桿件上任意一點的速度。然后根據式(10)可以計算9根桿件的動能之和。

(15)

勢能損失量PE同2.1節一樣，將網殼的內能、動能和勢能損失量代入式(1)所示的拉格朗日方程，并通過MATLAB編程求解，即可得到結構6個節點的位移響應，進而計算得到網殼結構的內能和動能時程。

2.3 K8型網殼理論計算與有限元結果的對比分析

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

(a) 俯視圖

(b) 正視圖

鋼材選用彈性材料，參數設置同第1節。在網殼頂點及第一環的三個頂點上施加集中荷載，由于模型的對稱性，作用在1點的荷載為P(t)/4，作用在2點的荷載為P(t)，作用在22點和42點的荷載為P(t)/2。P(t)隨時間的變化同2.1節。

約束模型底部5個節點的三向位移、對節點1和22施加Y軸對稱約束、對節點1和42施加X軸對稱約束，不同荷載下的響應結果皆顯示：X向和Y向位移較小且數值接近，兩者的位移峰值響應僅為Z向位移響應的28.36%，表明位移響應是以Z向為主，因此本文后續分析中僅給出Z向位移響應結果。當荷載峰值P0為50 kN時，提取K8型球面網殼結構在不同時刻的Z向位移如圖9所示。

從圖9可以看出，三角形爆炸荷載作用下網殼結構的位移響應受高階振型影響非常大，這是由網殼自身結構特性決定的。同時，網殼桿件上節點相較于桿件與桿件的交點來說，位移響應有滯后的效應。網殼底部桿件發生響應的時間最遲，且位移響應較小。

荷載峰值P0分別為10 kN、20 kN、30 kN、40 kN、50 kN時，將有限元數值分析、基于拉格朗日方程的理論兩參數和三參數分析得到的網殼最大節點位移(Z向)時程對比如圖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可以看出，兩參數和三參數理論分析模型得到的最大位移節點按照正弦形式振蕩，節點都在0.001 5 s左右達到正向峰值位移，然后開始衰減，經過平衡位置后負向位移逐漸增大，并且三參數分析中節點再次回到平衡位置的時間要比兩參數分析滯后，荷載峰值P0越大，這種規律越明顯。當荷載峰值P0從10 kN增大到50 kN時，有限元數值模擬與基于拉格朗日動力方程計算得到的網殼最大位移節點時程(Z向)的正向峰值位移對應較好，但理論分析達到正向峰值位移的時刻比有限元分析滯后一倍左右。同時，基于拉格朗日動力方程兩參數計算得到的位移響應比三參數計算結果偏大，這是因為兩參數計算過程中沒有考慮桿件的彎曲，外力功僅轉變為軸向變形能和動能，而三參數計算過程中外力功轉變為軸向變形能、彎曲變形能和動能。

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

其中，第一個和第二個原因最關鍵，第二個原因也是導致理論分析模型兩參數和三參數分析結果相差較小的主要原因。可以通過增加理論分析模型變形函數中的參數個數以及單獨假設每根桿件的變形函數來改善，但這樣將導致MATLAB編程量大幅度增加。

2.4 爆炸荷載下K8型網殼動力響應計算方法的修正

根據2.3節中基于拉格朗日動力方程理論分析與有限元分析產生誤差的原因，對2.2節中的理論分析模型進行改進，改進方法介紹如下。

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

圖11 桿件變形圖

從圖11可以看出，上部桿件的變形接近1/2周期的三角函數，下部桿件的變形接近1/4周期的三角函數，于是，在網殼桿件交點變形函數(式(7))的基礎上，對2.2節中桿件的變形函數進行修正。其中，理論分析模型中桿件編號見圖7所示。

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

(16)

(17)

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

(18)

(19)

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

(20)

(21)

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

(22)

(23)

勢能損失量PE同2.1節一致，計算出網殼的內能、動能和勢能損失量，代入式(1)所示的拉格朗日方程，并通過MATLAB編程求解，即可得到結構6個節點的位移響應，進而計算得到結構的內能和動能時程。

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

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

(a) P0=10 kN

(b) P0=20 kN

(c) P0=30 kN

(d) P0=40 kN

(e) P0=50 kN

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

綜合圖12和表2，本文認為改進后的三參數理論分析模型得到的網殼最大節點位移(Z向)振動規律比原三參數分析更貼合有限元分析結果。

將有限元分析、基于拉格朗日方程的理論三參數分析以及改進后的三參數分析得到的網殼系統動能、內能時程對比分別如圖13和圖14所示。

從圖13可以看出，改進后的三參數分析得到的網殼結構的動能時程以及達到動能峰值的時刻與有限元結果的振動規律更為接近，理論分析得到的總動能最大值均出現在時程的第二個峰值，有限元的總動能最大值出現在時程的第二個或第三個峰值。

表2 網殼位移響應對比

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

(e) P0=50 kN

從圖14可以看出，改進后的三參數理論分析模型得到的內能峰值、達到內能峰值所用的時間均比原三參數分析結果更接近有限元分析，當荷載峰值P0為30 kN和40 kN時，改進后的三參數理論分析得到的內能最大值與有限元分析結果誤差最小。

(a) P0=10 kN

(b) P0=20 kN

(c) P0=30 kN

(d) P0=40 kN

(e) P0=50 kN

綜合考慮最大節點位移響應、網殼結構的動能和內能情況，建議采用基于拉格朗日動力方程的理論分析模型計算網殼在爆炸荷載作用下的動力響應時，荷載峰值P0不小于20 kN。

3 結 論

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

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

(2) 當荷載峰值P0從10 kN增大到50 kN時，有限元數值模擬與基于拉格朗日動力方程計算得到的網殼最大位移節點時程(Z向)的正向峰值位移對應較好，但理論分析達到正向峰值位移的時刻比有限元分析滯后一倍左右。

(3) 改進后三參數分析模型得到的最大節點位移(Z向)時程、網殼結構的總動能和總內能與有限元分析結果更為接近。建議如采用該理論分析模型計算網殼動力響應時，荷載峰值P0應不小于20 kN，且假設的變形函數應接近網殼實際變形。

[1] 馬加路，支旭東，范峰,等. 外爆荷載下Kiewitt8型單層球面網殼的動力響應[J]. 振動與沖擊, 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] 翟希梅，黃明. 偏心爆炸荷載下網殼結構的動力響應分析[J]. 振動與沖擊. 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] 高軒能，李超，江媛. 單層球面鋼網殼結構在內爆炸作用下的動力響應[J]. 天津大學學報(自然科學與工程技術版), 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] 田力，李靈聚. 雙層柱面網殼結構在地下隧道內爆炸沖擊下的動力響應分析[J]. 振動工程學報, 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.