SPH-FEM耦合算法在爆炸螺栓解鎖分離過程中的應(yīng)用

初文華，張阿漫，明付仁，楊文山

(哈爾濱工程大學 船舶工程學院，哈爾濱 150001)

爆炸螺栓的分離解鎖是一個相當復雜的過程，目前尚無十分理想的理論計算方法。由于高能炸藥的爆轟及沖擊波對結(jié)構(gòu)的破壞作用將導致結(jié)構(gòu)大變形等問題，傳統(tǒng)的有限元方法在模擬爆炸螺栓的解鎖分離過程時將遇到極大阻礙，甚至是計算崩潰[1]。相比于傳統(tǒng)的基于網(wǎng)格的有限元法，光滑粒子流體動力學(SPH)方法[1]則在這方面具有很大的優(yōu)勢。它將系統(tǒng)的狀態(tài)用一系列包含著獨自的材料性質(zhì)且按照守恒控制方程規(guī)律運動的粒子來描述，同時具有無網(wǎng)格特性和拉格朗日性質(zhì)，因此能夠克服在計算中與大變形相關(guān)的困難。但是光滑粒子流體動力學(SPH)方法在每個時間步內(nèi)都要對問題域內(nèi)的所有粒子進行搜索，對于三維問題來說，這將消耗大量的計算時間，其計算效率相比于有限元方法大大降低。為了解決上述問題，綜合了FEM方法計算效率高和SPH方法計算精度高的優(yōu)點的耦合SPH-FEM計算方法被提出。

關(guān)于 SPH 與 FEM 的耦合算法，Johnson[2]和 Attaway 等[3]較早開展了相關(guān)研究，此后國內(nèi)外學者[4－8]在其工作的基礎(chǔ)上分別對SPH-FEM耦合算法進行了不同程度的改進，且逐漸將其應(yīng)用到高速沖擊、侵徹等問題的數(shù)值模擬。然而對于該方法在爆炸沖擊問題方面的應(yīng)用，相關(guān)研究開展的則相對較少。因此，本文基于SPH-FEM耦合算法，采用復合鏈表搜索方式，在充分提高計算效率的同時，計及爆炸螺栓與其連接結(jié)構(gòu)間的相互作用，建立爆炸螺栓解鎖分離過程三維數(shù)值模型，模擬分析爆炸螺栓的分離沖擊特性及結(jié)構(gòu)在分離沖擊激勵作用下的動響應(yīng)特性。

1 計算方法

1.1 控制方程

在爆炸沖擊的強烈作用下，固體材料的力學特性將會產(chǎn)生類似流體性質(zhì)的變化，此時運動方程和高壓狀態(tài)方程成為控制材料力學行為特性的關(guān)鍵描述。本文采用計及材料強度的流體動力學控制方程[1]計算在爆炸螺栓沖擊激勵作用下結(jié)構(gòu)的動力學特性，基于SPH方法離散后的控制方程如式(1)所示。

式中:p，ρ，e，vα，xα，σαβ，t分別代表壓力、密度、內(nèi)能、速度分量、空間坐標、總應(yīng)力張量和時間;為光滑核函數(shù)及其導數(shù)，本文采用分段三次樣條函數(shù)作為光滑核函數(shù);同時，式中還考慮了人工粘度∏ij和人工熱量Hi的影響。

1.2 復合鏈表搜索

在SPH方法的實現(xiàn)過程中，鄰近粒子搜索很大程度上決定了計算的效率。在光滑長度為空間常量的情況下，鏈表搜索法[1]是提高搜索效率的一種有效的方法。它通過將胞元看做簿記式的記錄形式來節(jié)省大量的計算時間。將所有的粒子都分布在胞元內(nèi)，并且使用鏈表定義，在NNPS過程中只要使用任一組粒子進行搜索即可，從而在很大程度上節(jié)省計算時間。

在實現(xiàn)鏈表算法時，要在問題域上鋪設(shè)一臨時網(wǎng)格。同時將每個粒子都分布在網(wǎng)格單元內(nèi)，并通過簡單的存儲規(guī)則將每個網(wǎng)格內(nèi)的所有粒子連接起來。若每個單元內(nèi)的平均粒子數(shù)量足夠小，鏈表搜索法將在極大程度上提高計算效率。然而，當問題域的跨度較大，由此導致的臨時網(wǎng)格區(qū)域鋪設(shè)過大時，對于三維問題來說，搜索的網(wǎng)格總數(shù)將會大大增加，從而導致搜索效率在一定程度上降低。本文所計算的模型中，每個螺栓的尺寸并不是很大，然而兩個螺栓之間的距離卻是螺栓尺寸的幾倍，如果將兩個螺栓放在同一個問題域中，采用鏈表搜索法不但不會提高搜索效率，反而會因過多沒有粒子存在的網(wǎng)格參與搜索循環(huán)而使得計算時間大大增加。

因此，基于模型距離跨度大的特點，本文采用復合鏈表搜索方法來提高搜索效率。如圖1所示，在每個螺栓所在的區(qū)域分別鋪設(shè)一局部臨時網(wǎng)格，由于兩個螺栓距離很遠，一個螺栓的粒子不會存在于另一個螺栓的支持域內(nèi)。由此，每個局部臨時網(wǎng)格內(nèi)的粒子僅參與該螺栓所在區(qū)域的粒子搜索，而不會對另一個螺栓區(qū)域的粒子搜索產(chǎn)生影響。此時，兩個螺栓之間的區(qū)域由于沒有粒子的存在則不需要鋪設(shè)網(wǎng)格，網(wǎng)格總數(shù)大大降低，同時每個螺栓所在區(qū)域分別采用各自的鏈表搜索系統(tǒng)，則每個搜索區(qū)域的網(wǎng)格數(shù)僅為整個問題域網(wǎng)格總數(shù)的一半，因此搜索效率在極大程度上得到了提高。

圖1 復合鏈表搜索方式示意圖Fig.1 Compound linked list search mode

1.3 耦合SPH-FEM技術(shù)

在本文中，光滑粒子流體動力學(SPH)方法與有限元(FEM)方法耦合的基本思想是在初始時刻爆炸螺栓及其連接結(jié)構(gòu)分別采用SPH粒子和有限元網(wǎng)格進行建模，在計算過程中，使交界面上的粒子與有限元邊界保持接觸但不允許相互侵入。

如圖2所示，當接觸面處的光滑粒子侵入有限元界面時，按下式對該邊界粒子進行位置和速度的調(diào)整:

圖2 SPH-FEM耦合計算方法原理圖Fig.2 Schematic diagram of coupling SPH-FEM method

與此同時，為了實現(xiàn)耦合計算，接觸面處的有限元界面上的單元壓力由各單元中心附近的光滑粒子壓力值取平均值得到:

其中:Pb_FEM為接觸面處有限元單元上的壓力，Pjb_SPH為該單元中心附近的光滑粒子壓力值;K為該有限元單元中心點“支持域”內(nèi)的光滑粒子總數(shù)，本文中該“支持域”尺寸定義為與光滑粒子的支持域尺寸相同。

2 計算模型

本文針對爆炸螺栓基本結(jié)構(gòu)及工作原理，采用耦合SPH-FEM算法建立爆炸螺栓爆炸解鎖分離數(shù)值模型。爆炸螺栓及其內(nèi)部裝藥采用SPH粒子進行離散，而其連接結(jié)構(gòu)則采用有限元網(wǎng)格單元進行模擬計算。

2.1 材料模型

本文的計算模型中，爆炸螺栓及其連接結(jié)構(gòu)均采用鋼材料進行模擬。在采用SPH方法對爆炸螺栓進行建模時，將Jaumann應(yīng)力率應(yīng)用到材料本構(gòu)模型的基本方程中，即在式(1)中，剪切應(yīng)力ταβi與應(yīng)變率εαβi存在如下關(guān)系:

高校應(yīng)高度重視政府會計制度落實，高校主要領(lǐng)導要高度重視，確保這項工作能夠落到實處。高校要成立由單位主要領(lǐng)導牽頭，財務(wù)、資產(chǎn)等相關(guān)業(yè)務(wù)部門組成的領(lǐng)導小組，規(guī)劃方案、明確目標、落實責任，加強部門協(xié)同，按責任分工。

式中:Rβγ為扭轉(zhuǎn)率應(yīng)變率 εαβ=

在應(yīng)力更新的過程中，采用Von Mises屈服準則來判斷材料的屈服狀態(tài)，材料的屈服應(yīng)力Y采用Johnsoncook材料屈服模型計算，其表達式為:

若SPH方法的計算過程中，粒子i的等效應(yīng)力σeff超過材料的屈服強度Y，則該粒子進入塑性屈服狀態(tài)，其剪切應(yīng)力ταβ按下式進行修正:

2.2 爆炸螺栓三維SPH 模型

本文基于常用爆炸螺栓的基本結(jié)構(gòu)，建立簡化的爆炸螺栓三維SPH數(shù)值模型如圖3所示，對于爆炸螺栓基本結(jié)構(gòu)中的螺栓頭、螺母、裝藥及用以減少爆炸螺栓沖擊環(huán)境的防護盒，該SPH模型中均進行了模擬。

在SPH-FEM耦合計算方法中，一個關(guān)鍵問題就是SPH粒子與有限元網(wǎng)格在接觸面上的耦合，而解決這一問題的首要因素就是保證二者初始時刻在接觸面處吻合。本文在自主編寫計算程序進行爆炸螺栓三維SPH建模的過程中，嚴格按照有限元模型中兩個螺栓孔的尺寸和位置，分別建立了上下兩個爆炸螺栓的數(shù)值模型。

圖3 爆炸螺栓三維簡化SPH模型Fig.3 Three-dimensional simplified SPH model of explosion bolt

2.3 分離結(jié)構(gòu)有限元模型

在數(shù)值建模的過程中，將爆炸螺栓的連接結(jié)構(gòu)簡化為如圖4(a)所示的對稱半圓框結(jié)構(gòu)。圖中的圓孔部分代表四個相同尺寸爆炸螺栓的安裝位置。

圖4 分離結(jié)構(gòu)及爆炸螺栓位置示意圖Fig.4 Separating structure and position of the four explosion bolts

3 計算結(jié)果

3.1 爆炸螺栓爆炸沖擊特性分析

為分析爆炸螺栓解鎖分離過程沖擊激勵的特點，圖5給出了爆炸螺栓解鎖分離過程中典型時刻的壓力云圖。

從圖中可以看到，兩爆炸螺栓從各自的裝藥中心同時起爆后，沖擊波從炸藥中心向螺栓頭、螺母等周圍結(jié)構(gòu)中傳播，并隨時間迅速衰減;炸藥周圍的結(jié)構(gòu)逐漸膨脹變形，直至最終產(chǎn)生破壞斷裂。由于模型尺寸較小，且為三維模型，為了更清楚地看到?jīng)_擊壓力在結(jié)構(gòu)中的分布及其隨時間的變化，圖中給出了放大后的爆炸螺栓結(jié)構(gòu)壓力分布圖，并將圖例中壓力最大值的顯示范圍調(diào)小。此外還可以看到，由于上下兩螺栓結(jié)構(gòu)尺寸完全相同，裝藥量也完全相同，兩爆炸螺栓在解鎖分離過程中的沖擊特性也幾乎完全相同。

圖5 爆炸螺栓解鎖分離過程沖擊壓力分布圖Fig.5 Impact pressure at different moments in the unlocking process of explosion bolt

圖6 沖擊激勵作用下爆炸螺栓連接結(jié)構(gòu)不同時刻應(yīng)力分布圖Fig.6 Stress distribution of the structure connected with explosion bolts under the impact excitation at different moments

3.2 爆炸螺栓沖擊激勵作用下的結(jié)構(gòu)響應(yīng)分析

圖6 給出了爆炸螺栓分離沖擊激勵作用下不同時刻連接結(jié)構(gòu)的MISES應(yīng)力分布圖。從圖6中可以看出，當爆炸螺栓內(nèi)裝藥起爆后，沖擊波通過接觸面向連接結(jié)構(gòu)中傳播;隨著時間增加，沖擊波在傳播過程中逐漸衰減，結(jié)構(gòu)中的MISES應(yīng)力響應(yīng)也逐漸減小，最終趨于穩(wěn)定;而爆炸螺栓連接的前后兩部分結(jié)構(gòu)也在爆炸分離沖擊激勵的作用下最終分離。

從圖6中還可以看到，應(yīng)力較大的位置始終出現(xiàn)在上下兩螺栓孔附近區(qū)域。圖7進一步給出了沖擊激勵作用下連接結(jié)構(gòu)上下兩螺栓孔附近位置處的MISES應(yīng)力時程曲線。

從圖7中可以看到，上下兩螺栓孔附近測點的MISES應(yīng)力響應(yīng)峰值均達近240 MPa，因此在爆炸螺栓解鎖分離過程的前期，強烈的沖擊激勵極有可能對爆炸螺栓附近連接結(jié)構(gòu)可能產(chǎn)生較高的破壞效應(yīng)。

為進一步分析爆炸螺栓解鎖分離過程中連接結(jié)構(gòu)在沖擊激勵作用下的響應(yīng)特性，圖8給出了連接結(jié)構(gòu)在上下兩螺栓孔附近區(qū)域的三方向沖擊加速度時程曲線。

在爆炸沖擊激勵作用下，應(yīng)力波在結(jié)構(gòu)中以三維空間波的方式傳播，結(jié)構(gòu)中的加速度響應(yīng)在極短時間內(nèi)迅速突躍到最大值，然后隨著時間增加幅值逐漸減小為零，整個持續(xù)時間不過幾十毫秒，然而加速度脈寬僅為零點幾毫秒，可見應(yīng)力波作用期間波形變化極快。此外，即使是同一點，加速度峰值在三個方向的大小也是不同的。

圖7 上下兩螺栓孔附近MISES應(yīng)力時程曲線Fig.7 MISES stress-time history curves of points near bolt holes

工程應(yīng)用中，主要把加速度及其計算譜作為爆炸沖擊環(huán)境評估的主要參數(shù)，設(shè)計規(guī)范和試驗規(guī)范都是以沖擊響應(yīng)譜為參數(shù)給出的。最大譜是爆炸沖擊試驗中最常用的一種，它是一組理想單自由度振子對基礎(chǔ)運動的最大響應(yīng)隨振子頻率ω變化的圖譜。圖9給出了根據(jù)圖8中六組沖擊加速度響應(yīng)時程曲線得到的六組最大沖擊響應(yīng)譜。

圖8 上下兩螺栓三個方向的加速度響應(yīng)時程曲線Fig.8 Acceleration-time history curves of the two explosion bolts in the three direction

圖9 上下兩螺栓三個方向的加速度響應(yīng)譜Fig.9 Acceleration response spectrums of the two explosion bolts in the three direction

4 結(jié)論

本文針對爆炸螺栓解鎖分離過程的工作原理及載荷特點，采用光滑粒子流體動力學(SPH)方法與有限元(FEM)方法相結(jié)合的技術(shù)手段，對爆炸螺栓的解鎖分離過程進行了數(shù)值模擬研究，得出了以下主要結(jié)論:

(1)兩爆炸螺栓分別從各自的裝藥中心同時起爆后，沖擊波分別向兩爆炸螺栓結(jié)構(gòu)中傳播，并隨時間迅速衰減;同時，炸藥周圍的結(jié)構(gòu)逐漸膨脹變形直至最終產(chǎn)生破壞斷裂;

(2)在爆炸螺栓分離沖擊激勵作用下，沖擊波從螺栓孔附近開始向結(jié)構(gòu)四周傳播，并隨時間增加逐漸衰減，結(jié)構(gòu)中的MISES應(yīng)力響應(yīng)也隨之逐漸減小，最終趨于穩(wěn)定;爆炸螺栓連接的前后兩部分結(jié)構(gòu)在爆炸分離沖擊激勵的作用下最終分離;

(3)在爆炸沖擊激勵作用下，應(yīng)力波在結(jié)構(gòu)中以三維空間波的方式傳播，結(jié)構(gòu)中的加速度響應(yīng)在極短時間內(nèi)迅速突躍到最大值，然后隨著時間增加幅值逐漸減小為零，整個持續(xù)時間不過幾十毫秒。

[1]韓 旭，楊 剛，強洪夫譯.光滑粒子流體動力學——一種無網(wǎng)格粒子法[M].長沙:湖南大學出版社，2005.

[2]Johnson G R，Stryk R A，Beissel S R，et al.An algorithm to automatically convert distorted finite elements into meshless particle during dynamics deformation[J].Int J Impact Eng.，2002，27(10):997－1013.

[3]Attaway S，Heinstein M，Swegle J.Coupling of smoothed particle hydrodynamics with the finite element method [J].Nucl.Eng.Des.，1994，150(2 －3):199 －205.

[4]王 吉，王肖鈞，卞 梁.光滑粒子法與有限元的耦合算法及其在沖擊動力學中的應(yīng)用[J].爆炸與沖擊，2007，27(6):522－528.

[5]卞 梁，王肖鈞，章 杰.SPH/FEM耦合算法在陶瓷復合靶抗侵徹數(shù)值模擬中的應(yīng)用[J].高壓物理學報，2010，24(3):161－167.

[6]Jun-xiang X U，LIU X L.Analysis of structural response under blast loads using the coupled SPH-FEM approach[J].Journal of Zhejiang University，2008，9(9):1184 －1192.

[7]Zang M Y，Chen H，Lei Z.Simulation on high velocity impact process of windshield by SPH/FEM coupling method [J].2010 WASE International Conference on Information Engineering，2010:381 －384.

[8]Fourey G，Oger G，Touze D Le，et al.Violent fluid-structure Interaction simulations using a coupled SPH/FEM method[J].Materials Science and Engineering，2010，10:1－8.