球?qū)ΨQ形式物質(zhì)點(diǎn)法及其在水下爆炸中的應(yīng)用＊

陳衛(wèi)東，楊文淼，嚴(yán) 涵，張豐超

（哈爾濱工程大學(xué)航天與建筑工程學(xué)院，黑龍江 哈爾濱 150001）

對(duì)水下爆炸問(wèn)題的研究一直是爆炸沖擊動(dòng)力學(xué)領(lǐng)域研究的熱點(diǎn)，在軍事工業(yè)以及民用建設(shè)領(lǐng)域都具有重要的意義。水下爆炸涉及高能炸藥的爆轟、沖擊波在水介質(zhì)中的傳播和氣泡的脈沖等一系列復(fù)雜過(guò)程，理論與實(shí)驗(yàn)研究都具有一定的局限性。數(shù)值研究具有低成本、可重復(fù)和過(guò)程細(xì)節(jié)易于捕捉等優(yōu)勢(shì)，在水下爆炸研究領(lǐng)域中起著越來(lái)越重要的作用。水下爆炸過(guò)程存在大變形問(wèn)題，使得基于隨體網(wǎng)格的拉格朗日方法產(chǎn)生嚴(yán)重的網(wǎng)格畸變，影響計(jì)算精度，甚至使計(jì)算無(wú)法進(jìn)行，所以大都采用具有空間固定坐標(biāo)的歐拉方法對(duì)水下爆炸進(jìn)行數(shù)值分析［1－4］。但歐拉方法不易觀察材料的變形歷史，并且物質(zhì)界面的處理也是歐拉方法面臨的一個(gè)難題，這些都限制了歐拉方法的應(yīng)用。近年來(lái)無(wú)網(wǎng)格方法如SPH（smoothed particle hydrodynamics）、MPM（material－point method）等因?yàn)椴淮嬖诰W(wǎng)格畸變，可以很自然地處理大變形問(wèn)題，在爆炸沖擊問(wèn)題中得到了廣泛應(yīng)用［5－12］。MPM由流體動(dòng)力學(xué)領(lǐng)域的PIC（particle－in－cell）法發(fā)展而來(lái)，兼具拉格朗日法和歐拉法的優(yōu)點(diǎn)，在處理大變形以及多介質(zhì)耦合問(wèn)題時(shí)具有突出優(yōu)勢(shì)［13］。MPM將問(wèn)題域離散成一系列攜帶材料性質(zhì)的物質(zhì)點(diǎn)集，對(duì)于三維水下爆炸問(wèn)題，尤其是三維遠(yuǎn)場(chǎng)爆炸問(wèn)題，問(wèn)題域的空間尺度是巨大的，并且隨著問(wèn)題維度的增加，物質(zhì)點(diǎn)數(shù)量也會(huì)成指數(shù)增長(zhǎng)。對(duì)于遠(yuǎn)場(chǎng)水下爆炸，人們關(guān)心的是結(jié)構(gòu)在爆炸沖擊波作用下的沖擊響應(yīng)，而對(duì)爆炸沖擊波在水介質(zhì)中傳播過(guò)程的計(jì)算會(huì)占用大量的計(jì)算資源，嚴(yán)重影響求解效率。當(dāng)忽略界面或結(jié)構(gòu)的影響時(shí)，水下爆炸產(chǎn)生的沖擊波在水介質(zhì)中具有球?qū)ΨQ性。利用水下爆炸這一特點(diǎn)，將沖擊波在自由場(chǎng)中的傳播過(guò)程簡(jiǎn)化為一維球?qū)ΨQ模型進(jìn)行計(jì)算，這樣可以在保證求解精度的前提下提高計(jì)算效率。基于AUTODYN中的映射求解思想［14］，建立了基于物質(zhì)點(diǎn)法的水下爆炸重映射求解模型，當(dāng)爆炸沖擊波到達(dá)界面或者結(jié)構(gòu)之前，采用一維球?qū)ΨQ模型進(jìn)行計(jì)算，然后將計(jì)算結(jié)果映射到三維模型中，并將各計(jì)算量作為初值在三維模型中繼續(xù)完成計(jì)算，有效避免了計(jì)算資源在沖擊波傳播過(guò)程上的消耗。

本文中，首先推導(dǎo)一維球?qū)ΨQ形式的物質(zhì)點(diǎn)法，然后對(duì)球形炸藥水下爆炸過(guò)程進(jìn)行數(shù)值計(jì)算，并與Cole經(jīng)驗(yàn)公式［15］以及DYNA計(jì)算結(jié)果進(jìn)行比較，結(jié)果吻合較好，驗(yàn)證了一維球?qū)ΨQ形式物質(zhì)點(diǎn)法的準(zhǔn)確性。在此基礎(chǔ)上，提出了基于物質(zhì)點(diǎn)法的重映射求解算法，將一維球?qū)ΨQ結(jié)果作為初始條件映射到笛卡爾坐標(biāo)系下的三維模型中，從而可以在三維模型中繼續(xù)求解沖擊波的傳播以及結(jié)構(gòu)的沖擊響應(yīng)過(guò)程。

1 球?qū)ΨQ形式的物質(zhì)點(diǎn)法

圖1 球坐標(biāo)系下的離散單元Fig.1 Discrete elements in spherical coordinate system

在一維球?qū)ΨQ坐標(biāo)下離散的質(zhì)點(diǎn)質(zhì)量是不均勻的，其質(zhì)量大小與質(zhì)點(diǎn)的密度分布以及徑向坐標(biāo)都有關(guān)系。如圖1所示，在球坐標(biāo)系下沿徑向離散的單元坐標(biāo)為r1和r2，單位立體角內(nèi)離散單個(gè)質(zhì)點(diǎn)的質(zhì)量為：

式中：ρ為密度，a為加速度，σ為柯西應(yīng)力張量，b為體力，τ為邊界外力，w為試函數(shù)。

在MPM中，將連續(xù)體離散成一系列物質(zhì)點(diǎn)集，其密度可以寫成：

物質(zhì)點(diǎn)法中需要布置歐拉網(wǎng)格覆蓋求解域，同有限元法類似，背景網(wǎng)格單元內(nèi)的場(chǎng)變量可表示為：

式中：g（x，t）為任意場(chǎng)變量，可代表速度、加速度和試函數(shù)等；Ni（x）為單元形函數(shù)。

利用式（3）和式（4），式（2）可寫成如下形式：

在球?qū)ΨQ形式下，各物理量只是r的函數(shù)，與θ和φ無(wú)關(guān)，應(yīng)力張量有如下形式：

在球?qū)ΨQ形式下，其內(nèi)力向量只有沿徑向的分量，其徑向分量為：

采用顯式差分格式求解式（6），從而得到節(jié)點(diǎn)的更新速度：

式中：上標(biāo)k和k＋1分別代表上一時(shí)間步和更新后的時(shí)間步。

利用更新后的節(jié)點(diǎn)速度更新物質(zhì)點(diǎn)的位置和速度：

在球?qū)ΨQ形式下，物質(zhì)點(diǎn)的應(yīng)變?cè)隽恐慌cr有關(guān)，可以寫成：

得到物質(zhì)點(diǎn)的應(yīng)變?cè)隽亢?，通過(guò)材料的本構(gòu)方程便可完成物質(zhì)點(diǎn)的應(yīng)力以及其他內(nèi)部變量的更新。在球?qū)ΨQ形式下的應(yīng)力更新過(guò)程與笛卡爾坐標(biāo)系下的計(jì)算過(guò)程沒(méi)有太大的區(qū)別。

2 水下爆炸的數(shù)值模擬

2.1 數(shù)學(xué)模型描述

對(duì)球形TNT裝藥在水下自由場(chǎng)中的爆炸過(guò)程進(jìn)行了數(shù)值模擬。為了簡(jiǎn)化理論分析，忽略水下爆炸過(guò)程中熱傳導(dǎo)和黏性的影響。水下爆炸過(guò)程可由下列一組歐拉方程和恰當(dāng)?shù)臓顟B(tài)方程來(lái)模擬：

式中：ρ、e、p、v和t分別為密度、內(nèi)能、壓力、速度矢量和時(shí)間。

式（13）中前3個(gè)方程分別描述了質(zhì)量守恒、動(dòng)量守恒和能量守恒，最后一個(gè)方程代表不同材料的狀態(tài)方程。狀態(tài)方程用來(lái)定義固體或流體在不同狀態(tài)下的壓力和密度以及比內(nèi)能之間的函數(shù)關(guān)系。正確選取狀態(tài)方程以及對(duì)狀態(tài)方程中的參數(shù)設(shè)定合適的數(shù)值對(duì)于計(jì)算結(jié)果至關(guān)重要。

本文中，對(duì)水采用多項(xiàng)式狀態(tài)方程：

對(duì)TNT炸藥采用JWL狀態(tài)方程：

在求解水下爆炸沖擊問(wèn)題時(shí)，沖擊波將導(dǎo)致壓力、密度、質(zhì)點(diǎn)速度和能量等在波陣面前后發(fā)生突變（跳躍），這種不連續(xù)性給運(yùn)動(dòng)微分方程的求解帶來(lái)了較大的困難。為了解決此問(wèn)題，引入了人工黏性的概念。人工黏性將應(yīng)力波的強(qiáng)間斷面模糊成一個(gè)在相當(dāng)狹窄的過(guò)渡區(qū)域內(nèi)連續(xù)變化的波陣面，從而給工程應(yīng)用帶來(lái)極大的方便［16］。

本文中采用線性形式和二次形式相結(jié)合的人工黏性形式：

式中：c0和c1為人工黏性系數(shù)；le為特征長(zhǎng)度，在物質(zhì)點(diǎn)法中，取背景網(wǎng)格最小尺寸；ρ為現(xiàn)時(shí)構(gòu)形密度為體積應(yīng)變率，c為材料聲速。

2.2 計(jì)算結(jié)果與分析

圖2 一維球?qū)ΨQ計(jì)算模型Fig.2 One－dimensional spherically－symmetric model

計(jì)算模型如圖2所示，球形TNT裝藥半徑為0.1m，水域半徑為5m，離散網(wǎng)格尺寸為0.01m，每個(gè)網(wǎng)格放置2個(gè)物質(zhì)點(diǎn)，總粒子數(shù)為1000。

圖3給出了分別采用Cole經(jīng)驗(yàn)公式、MPM和DYNA程序得到的水中不同位置處沖擊波壓力峰值，圖中a為炸藥半徑，R為觀察點(diǎn)到起爆中心的距離。從圖3可以看出壓力峰值都隨水中沖擊波的傳播而逐漸減小，而且都呈指數(shù)形式衰減，這符合水下爆炸沖擊波的基本規(guī)律。在離起爆中心較近的地方，采用MPM和DYNA計(jì)算的結(jié)果與采用Cole經(jīng)驗(yàn)公式計(jì)算結(jié)果相對(duì)誤差較小，但隨著離起爆點(diǎn)距離的增加，相對(duì)誤差有增大的趨勢(shì)。與采用Cole經(jīng)驗(yàn)公式計(jì)算的結(jié)果相比，采用MPM和DYNA計(jì)算的壓力峰值衰減更快。這主要是因?yàn)镸PM和DYNA均采用理想化的計(jì)算模型，忽略了水的黏性、密度不均勻性以及重力等影響因素，而實(shí)際的水下爆炸過(guò)程很復(fù)雜，Cole經(jīng)驗(yàn)公式是由大量的實(shí)驗(yàn)數(shù)據(jù)擬合得到的，所以采用MPM和DYNA計(jì)算的結(jié)果與采用Cole經(jīng)驗(yàn)公式計(jì)算的結(jié)果之間的誤差是合理的。不同爆距處，采用MPM和DYNA計(jì)算結(jié)果與采用Cole經(jīng)驗(yàn)公式計(jì)算結(jié)果之間的相對(duì)誤差εr，MPM和εr，DYNA見(jiàn)表1。

表1 不同爆距處，采用MPM和DYNA計(jì)算的結(jié)果與采用Cole經(jīng)驗(yàn)公式計(jì)算的結(jié)果之間的相對(duì)誤差Table1 Relative errors between calculated results by MPM and DYNA and ones by Cole formula at different locations

圖4給出了幾個(gè)典型時(shí)刻水下沖擊波壓力沿徑向的分布曲線。從圖4可以看出，采用MPM計(jì)算得到的沖擊波前沿位置與沖擊波峰值大小與采用DYNA的計(jì)算結(jié)果吻合較好，沖擊波沿徑向以近似指數(shù)衰減的階梯式增長(zhǎng)。在炸藥爆轟產(chǎn)物與水介質(zhì)的交界面附近，有一定的數(shù)值振蕩產(chǎn)生，這是由爆轟產(chǎn)物與水介質(zhì)的不均一性導(dǎo)致的，但數(shù)值振蕩與沖擊波大小相比是可以接受的。隨著沖擊波背離起爆中心向外傳播，其波形被不斷拉寬，且峰值具有明顯的衰減現(xiàn)象。

本算例如果采用三維模型進(jìn)行計(jì)算，并假設(shè)采用與一維球?qū)ΨQ模型相同的網(wǎng)格離散精度，即在0.01m×0.01m×0.01m的離散網(wǎng)格內(nèi)放置8個(gè)物質(zhì)點(diǎn)，即使采用1／8模型，離散物質(zhì)點(diǎn)總數(shù)也將達(dá)到1000×1000×1000。三維模型計(jì)算規(guī)模遠(yuǎn)大于離散物質(zhì)點(diǎn)總數(shù)為1000的一維球?qū)ΨQ計(jì)算模型，使數(shù)據(jù)的計(jì)算量和存儲(chǔ)量變得很龐大，很難處理，需在大型并行機(jī)上完成計(jì)算。雖然通過(guò)降低網(wǎng)格離散精度，可以在一定程度上減小三維模型的計(jì)算規(guī)模，但即使采用0.1m×0.1m×0.1m網(wǎng)格尺寸進(jìn)行離散，其物質(zhì)點(diǎn)總數(shù)為100×100×100與一維球?qū)ΨQ模型相比，其計(jì)算規(guī)模仍很巨大。與三維計(jì)算模型相比，一維球?qū)ΨQ模型在保證網(wǎng)格離散精度的同時(shí)可以有效減小計(jì)算規(guī)模。

圖3 壓力峰值沿徑向的分布曲線Fig.3 Pressure peak distribution along radial direction

圖4 典型時(shí)刻的壓力分布曲線Fig.4 Pressure distribution curves at typical times

2.3 基于物質(zhì)點(diǎn)法的重映射算法

重映射計(jì)算需要用到球坐標(biāo)系與笛卡爾坐標(biāo)系的轉(zhuǎn)換矩陣，2個(gè)坐標(biāo)系下基矢存在如下關(guān)系：

式中：er、eθ和eφ為球坐標(biāo)系下單位基矢量，ex、ey和ez為笛卡爾坐標(biāo)系下單位基矢量，R為轉(zhuǎn)換矩陣。

對(duì)于矢量a，如速度、加速度和力等，在笛卡爾坐標(biāo)系與球坐標(biāo)系間有如下的轉(zhuǎn)換關(guān)系：

對(duì)于張量B，如應(yīng)力張量，在笛卡爾坐標(biāo)系與球坐標(biāo)系間有如下的轉(zhuǎn)換關(guān)系：

基于物質(zhì)點(diǎn)法的重映射計(jì)算，采用如下的求解步驟：

（1）將球?qū)ΨQ形式下物質(zhì)點(diǎn)攜帶的內(nèi)部變量，如密度、壓力和能量等映射到球?qū)ΨQ的網(wǎng)格節(jié)點(diǎn)上；

（2）通過(guò)式（4）得到球坐標(biāo)系下場(chǎng)變量的空間分布；

（3）通過(guò)式（18）和（19）將場(chǎng)變量轉(zhuǎn)換到笛卡爾坐標(biāo)系下；

（4）將笛卡爾坐標(biāo)系下場(chǎng)變量的分布空間作為一個(gè)連續(xù)體，采用物質(zhì)點(diǎn)進(jìn)行離散，其中在一維球?qū)ΨQ形式下的物質(zhì)分界點(diǎn)變成三維空間中的物質(zhì)分界面。

以上4步便完成了將一維球?qū)ΨQ計(jì)算結(jié)果到三維笛卡爾坐標(biāo)系下的映射計(jì)算過(guò)程。

圖5 映射后的三維計(jì)算模型Fig.5 Mapping to the three－dimensional computational model

圖5給出了算例中t＝1ms和t＝3ms時(shí)刻計(jì)算結(jié)果映射后的三維計(jì)算模型。一維球?qū)ΨQ形式下的計(jì)算結(jié)果將作為初值賦值給三維模型中離散的物質(zhì)點(diǎn)，以便在三維空間中繼續(xù)完成結(jié)構(gòu)沖擊響應(yīng)的計(jì)算。為了便于觀察，圖中顯示的是1／8模型。

3 結(jié) 論

針對(duì)求解三維遠(yuǎn)場(chǎng)水下爆炸沖擊響應(yīng)問(wèn)題，提出了一維球?qū)ΨQ形式物質(zhì)點(diǎn)法，并進(jìn)行了數(shù)值模擬驗(yàn)證，取得了如下結(jié)論：

（1）對(duì)比Cole公式、DYNA和一維球?qū)ΨQ形式物質(zhì)點(diǎn)法計(jì)算結(jié)果，在起爆中心附近MPM、DYNA和Cole公式計(jì)算結(jié)果間相對(duì)誤差較小，隨起爆點(diǎn)距離增加，相對(duì)誤差有增大趨勢(shì)，但仍在誤差允許范圍之內(nèi)。因此，采用一維球?qū)ΨQ形式物質(zhì)點(diǎn)法能夠較好模擬水下爆炸過(guò)程；

（2）在爆炸模擬過(guò)程中，通過(guò)重映射算法，將一維球?qū)ΨQ計(jì)算結(jié)果映射到三維計(jì)算模型中，避免了計(jì)算資源在沖擊波傳播過(guò)程計(jì)算上的消耗，顯著降低了三維遠(yuǎn)場(chǎng)水下爆炸沖擊響應(yīng)求解復(fù)雜度，具有較高的工程實(shí)用價(jià)值。

［1］師華強(qiáng)，汪玉，宗智，等.近水面水下爆炸二維Level＿set數(shù)值模擬［J］.計(jì)算力學(xué)學(xué)報(bào)，2010，27（3）：409－414.Shi Hua－qiang，Wang Yu，Zong Zhi，et al.2Dnumerical simulation of underwater explosion near free surface based on Level＿set method［J］.Chinese Journal of Computational Mechanics，2010，27（3）：409－414.

［2］李健，榮吉利，項(xiàng)大林.裝藥量及水深對(duì)水下爆炸氣泡動(dòng)態(tài)特性的影響［J］.爆炸與沖擊，2010，30（4）：324－348.Li Jian，Rong Ji－li，Xiang Da－lin.Effects of charge mass and water depth on dynamic behaviors of an underwater explosion bubble［J］.Explosion and Shock Waves，2010，30（4）：324－348.

［3］梁龍河，曹菊珍，袁仙春.水下爆炸特性的二維數(shù)值模擬研究［J］，高壓物理學(xué)報(bào)，2004，18（3）：203－208.Liang Long－h(huán)e，Cao Ju－zhen，Yuan Xian－chun.2DNumerical simulation of characteristics of underwater explosions［J］.Chinese Journal of High Pressure Physics，2004，18（3）：203－208.

［4］張振華，朱錫，白雪飛.水下爆炸沖擊波的數(shù)值模擬研究［J］.爆炸與沖擊，2004，24（2）：182－188.Zhang Zhen－h(huán)ua，Zhu Xi，Bai Xue－fei.The study on numerical simulation of underwater blast wave［J］.Explosion and Shock Waves，2004，24（2）：182－188.

［5］強(qiáng)洪夫，王坤鵬，高巍然.基于完全變光滑長(zhǎng)度SPH方法的高能炸藥爆轟過(guò)程數(shù)值試驗(yàn)［J］.含能材料，2009，17（1）：27－31.Qiang Hong－fu，Wang Kun－peng，Gao Wei－ran.Numerical simulation of high explosive detonation process using SPH method with fully variable smoothing lengths［J］.Chinese Journal of Energetic Materials，2009，17（1）：27－31.

［6］楊剛，韓旭，龍述堯.應(yīng)用SPH 方法模擬近水面爆炸［J］.工程力學(xué)，2008，25（4）：204－208.Yang Gang，Han Xu，Long Shu－yao.Simulation of underwater explosion near air－water surface by SPH method［J］.Engineering Mechanics，2008，25（4）：204－208.

［7］Liu G R，Liu M B，Zong Z，et al.Numerical simulation of underwater explosion by SPH［C］∥Atluri S N，Brust F W.Proceedings of the Advances in Computational Engineering ＆ Sciences.US：Tech Science Press，2000：1475－1480.

［8］Liu M B，Liu G R，Lam K Y，et al.Smoothed particle hydrodynamics for numerical simulation of underwater explosions［J］.Computational Mechanics，2003，30（2）：106－118.

［9］馬上，張雄.聚能裝藥射流形成的自適應(yīng)物質(zhì)點(diǎn)法模擬［J］.固體力學(xué)學(xué)報(bào)，2009，30（5）：504－508.Ma Shang，Zhang Xiong.Adaptive material point method for shaped charge jet formation［J］.Chinese Journal of Solid Mechanics，2009，30（5）：504－508.

［10］張忠，陳衛(wèi)東，楊文淼.非均質(zhì)固體炸藥沖擊起爆的物質(zhì)點(diǎn)法［J］.爆炸與沖擊，2011，31（1）：25－30.Zhang Zhong，Chen Wei－dong，Yang Wen－miao.The material point method for shock－to－detonation transition of heterogeneous solid explosive［J］.Explosion and Shock Waves，2011，31（1）：25－30.

［11］張忠，陳衛(wèi)東.基于物質(zhì)點(diǎn)法的超高速碰撞問(wèn)題研究［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2010，31（10）：1312－1316.Zhang Zhong，Chen Wei－dong.Solving hypervelocity impact problems with the material point method［J］.Journal of Harbin Engineering University，2010，31（10）：1312－1316.

［12］王宇新，陳震，張洪武，等.多層抗爆結(jié)構(gòu)沖擊響應(yīng)無(wú)網(wǎng)格 MPM 法分析［J］.工程力學(xué)，2007，24（12）：186－192.Wang Yu－xin，Chen Zhen，Zhang Hong－wu，et al.Response of multi－layered structure due to impact load using material point method［J］.Engineering Mechanics，2007，24（12）：186－192.

［13］Sulsky D，Chen Z，Schreyer H L.A particle method for history－dependent materials［J］.Computer Methods in Applied Mechanics and Engineering，1994，118（1／2）：179－196.

［14］Tham C Y.Numerical simulation on the interaction of blast waves with a series of aluminum cylinders at near－field［J］.International Journal of Impact Engineering，2009，36（1）：122－131.

［15］Cole R H.Underwater explosion［M］.New Jersey：LISA，Princeton University Press，1948：31－70.

［16］張雄，王天舒.計(jì)算動(dòng)力學(xué)［M］.北京：清華大學(xué)出版社，2007：273－274.