陳衛(wèi)東， 劉瀾， 路勝卓， 吳世博， 馬敬鑫， 吳培文

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

物質(zhì)點(diǎn)法(material point method，MPM)是一種由粒子算法與網(wǎng)格算法相結(jié)合的數(shù)值計(jì)算方法。該算法最早是由Sulsky等[1]率先提出的，此后國內(nèi)外許多學(xué)者都對該算法進(jìn)行了深入研究，并且在時(shí)間積分[2]、形函數(shù)改進(jìn)[3-5]、接觸算法[6]以及多種算法耦合[7-9]等多個(gè)方面對該算法進(jìn)行了改進(jìn)，從而使該算法不斷趨于成熟。因此在近年來，該算法已經(jīng)在許多工程計(jì)算領(lǐng)域中得到了廣泛應(yīng)用。

但是在爆炸與沖擊類工程計(jì)算領(lǐng)域中，由于計(jì)算對象往往會出現(xiàn)體積劇烈膨脹或結(jié)構(gòu)大尺度拉伸變形等狀況，例如炸藥爆轟產(chǎn)物膨脹擴(kuò)散以及破片飛濺等。所以在應(yīng)用物質(zhì)點(diǎn)法及其現(xiàn)有的改進(jìn)算法求解此類問題時(shí)，當(dāng)模型運(yùn)行至一定階段后，往往會出現(xiàn)背景網(wǎng)格區(qū)域容量不足以及物質(zhì)點(diǎn)陷入拉伸不穩(wěn)定狀態(tài)等不利情況，這些情況會對計(jì)算結(jié)果產(chǎn)生嚴(yán)重影響，從而制約物質(zhì)點(diǎn)法在爆炸與沖擊類工程計(jì)算領(lǐng)域中的應(yīng)用。因此現(xiàn)有必要對物質(zhì)點(diǎn)法作進(jìn)一步改進(jìn)，以使該算法能夠更加適用于求解爆炸與沖擊類工程實(shí)際問題。

1 物質(zhì)點(diǎn)法的顯式時(shí)間積分原理

在物質(zhì)點(diǎn)法中，計(jì)算對象被離散成一系列占據(jù)一定體積的物質(zhì)點(diǎn)，并且被設(shè)置于背景網(wǎng)格之中，物質(zhì)點(diǎn)儲存著其所在區(qū)域的材料的全部物理信息，而背景網(wǎng)格節(jié)點(diǎn)則只用于求解動量方程和空間導(dǎo)數(shù)。因此基于應(yīng)力后更新(update stress last，USL)求解格式的物質(zhì)點(diǎn)法顯式時(shí)間積分原理[10-11]的具體實(shí)現(xiàn)過程可表述如下。

1)將物質(zhì)點(diǎn)所攜帶的質(zhì)量、動量和速度映射到網(wǎng)格節(jié)點(diǎn)上，從而得到網(wǎng)格節(jié)點(diǎn)的質(zhì)量、動量和速度。該步驟的具體表達(dá)式分別為：

(1)

(2)

(3)

2)計(jì)算網(wǎng)格節(jié)點(diǎn)力。網(wǎng)格節(jié)點(diǎn)合力的表達(dá)式為：

(4)

3)計(jì)算下一時(shí)間步的節(jié)點(diǎn)動量。該步驟的表達(dá)式為：

(5)

式中Δtk表示時(shí)間步長。

4)將下一時(shí)間步的節(jié)點(diǎn)動量映射回物質(zhì)點(diǎn)，并更新物質(zhì)點(diǎn)的速度和位置坐標(biāo)。該步驟的表達(dá)式分別為：

(6)

(7)

5)計(jì)算物質(zhì)點(diǎn)的旋率增量和應(yīng)變增量。旋率增量和應(yīng)變增量的表達(dá)式分別為：

(8)

(9)

6)更新物質(zhì)點(diǎn)的密度和幾何尺寸。該步驟的表達(dá)式分別為：

(10)

(11)

7)更新物質(zhì)點(diǎn)的應(yīng)力。該步驟的表達(dá)式為：

(12)

8)重復(fù)上述步驟1)～7)，直至達(dá)到終止計(jì)算時(shí)間。

2 物質(zhì)點(diǎn)法的改進(jìn)

2.1 構(gòu)建物質(zhì)點(diǎn)自由流出邊界條件

2.1.1 物質(zhì)點(diǎn)自由流出邊界條件的作用與構(gòu)建方式

根據(jù)物質(zhì)點(diǎn)法的顯式時(shí)間積分原理可知，背景網(wǎng)格在每次時(shí)間積分計(jì)算結(jié)束后都會自動復(fù)原，并繼續(xù)參與下次時(shí)間積分計(jì)算。因此物質(zhì)點(diǎn)法的背景網(wǎng)格區(qū)域始終處于非嚴(yán)格封閉狀態(tài)，當(dāng)物質(zhì)點(diǎn)擴(kuò)散至背景網(wǎng)格邊界時(shí)，物質(zhì)點(diǎn)及其所攜帶的物理信息都會受到邊界節(jié)點(diǎn)的約束，而無法從背景網(wǎng)格邊界自由流出，導(dǎo)致數(shù)值計(jì)算結(jié)果產(chǎn)生巨大誤差。所以一般應(yīng)設(shè)置容積較大的背景網(wǎng)格區(qū)域，以盡量減小邊界節(jié)點(diǎn)約束對數(shù)值計(jì)算的影響。但是背景網(wǎng)格區(qū)域設(shè)置過大，又勢必會加大數(shù)值計(jì)算的工作量，導(dǎo)致計(jì)算成本大幅度提高，這是物質(zhì)點(diǎn)法所面臨的一種局限。為了破除物質(zhì)點(diǎn)法的這一局限，本文將基于物質(zhì)點(diǎn)法顯式時(shí)間積分原理來構(gòu)建物質(zhì)點(diǎn)自由流出邊界條件，以使物質(zhì)點(diǎn)及其物理信息能夠根據(jù)實(shí)際情況從背景網(wǎng)格邊界自由流出，從而在保證數(shù)值計(jì)算精確度的同時(shí)提高計(jì)算效率。下面對物質(zhì)點(diǎn)自由流出邊界條件的構(gòu)建方式進(jìn)行具體介紹。

在每次時(shí)間積分計(jì)算開始前，都要對模型中的全部物質(zhì)點(diǎn)進(jìn)行位置坐標(biāo)的重新確定，然后按照

(13)

2.1.2 算例分析

為了檢驗(yàn)物質(zhì)點(diǎn)自由流出邊界條件的有效性，本文將該邊界條件應(yīng)用到近物面水下爆炸算例中，并通過算例的對比分析結(jié)果來對該邊界條件的實(shí)際效果以及鋼板前后的正向沖擊波(入射波)、反射波和透射波的形態(tài)進(jìn)行量化描述。在近物面水下爆炸算例中，采用基于Fortran語言自編的物質(zhì)點(diǎn)法程序來構(gòu)建3個(gè)近物面水下爆炸二維模型，這3個(gè)模型在該算例中分別簡稱為模型A、模型B和模型C，這3個(gè)模型的初始形態(tài)如圖1所示，其中，模型A和模型B的矩形背景網(wǎng)格區(qū)域的面積均為100.0 cm×100.0 cm，模型A施加物質(zhì)點(diǎn)自由流出邊界條件，模型B不施加該邊界條件；而模型C的矩形背景網(wǎng)格區(qū)域的面積為200.0 cm×200.0 cm，且不施加物質(zhì)點(diǎn)自由流出邊界條件。除此之外，這3個(gè)模型的其他方面都是相同的，即：直角坐標(biāo)系的坐標(biāo)原點(diǎn)均設(shè)置在背景網(wǎng)格區(qū)域的形心位置，網(wǎng)格單元尺寸均為0.5 cm×0.5 cm，背景網(wǎng)格區(qū)域中均填加了面積為100.0 cm×100.0 cm的矩形物質(zhì)點(diǎn)區(qū)域，物質(zhì)點(diǎn)區(qū)域的中心部位構(gòu)建了尺寸為8.0 cm×8.0 cm的矩形TNT炸藥，炸藥右側(cè)構(gòu)建了厚度為2.0 cm的鋼板，鋼板的前壁面(迎爆面)與炸藥形心的距離為250 cm，物質(zhì)點(diǎn)區(qū)域的其他部位由水填充，離散物質(zhì)點(diǎn)的尺寸為0.25 cm×0.25 cm，在物質(zhì)點(diǎn)區(qū)域中從左至右依次設(shè)置了4個(gè)觀測點(diǎn)，這4個(gè)觀測點(diǎn)的坐標(biāo)見表1。

通過圖1可知，模型C在其物質(zhì)點(diǎn)區(qū)域的四周均預(yù)留出寬度為50 cm的背景網(wǎng)格擴(kuò)充區(qū)域，從而可給物質(zhì)點(diǎn)膨脹擴(kuò)散提供充足的空間，這是應(yīng)用物質(zhì)點(diǎn)法構(gòu)建數(shù)值模型的常規(guī)方式，所以模型C的計(jì)算結(jié)果可以認(rèn)定是正確的，因此本算例采用模型C的計(jì)算結(jié)果來作為模型A和模型B的計(jì)算結(jié)果的正確性檢驗(yàn)標(biāo)準(zhǔn)；而模型A和模型B的背景網(wǎng)格區(qū)域的尺寸與物質(zhì)點(diǎn)區(qū)域的尺寸相同，并沒有給物質(zhì)點(diǎn)膨脹擴(kuò)散提供額外的網(wǎng)格空間，因此可通過將這2個(gè)模型的計(jì)算結(jié)果與模型C的計(jì)算結(jié)果進(jìn)行相互比較，以檢驗(yàn)施加于模型A的物質(zhì)點(diǎn)自由流出邊界條件的有效性。在分別對這3個(gè)模型進(jìn)行數(shù)值計(jì)算后，得到了這3個(gè)模型中4個(gè)觀測點(diǎn)的壓力時(shí)程數(shù)據(jù)以及這3個(gè)模型在各時(shí)刻的壓力云圖。

圖1 3個(gè)模型的初始形態(tài)Fig.1 Initial form of 3 models

現(xiàn)分別提取距離背景網(wǎng)格左側(cè)邊界5 cm的觀測點(diǎn)1和距離背景網(wǎng)格右側(cè)邊界5 cm的觀測點(diǎn)4的壓力時(shí)程數(shù)據(jù)，繪制出圖2所示的這2個(gè)觀測點(diǎn)的壓力時(shí)程曲線。

表1 近物面水下爆炸模型的觀測點(diǎn)坐標(biāo)

通過觀察圖2可知，由于模型C給物質(zhì)點(diǎn)膨脹擴(kuò)散提供了網(wǎng)格空間，所以模型C中水下爆炸形成的正向沖擊波以及鋼板前后形成的反射波和透射波都能夠跟隨物質(zhì)點(diǎn)的膨脹擴(kuò)散而被自由釋放，因此在450 μs的計(jì)算時(shí)間內(nèi)，模型C中這2個(gè)觀測點(diǎn)的壓力時(shí)程曲線都只具有1個(gè)波峰，其中觀測點(diǎn)1的波峰是由正向沖擊波形成的，觀測點(diǎn)4的波峰是由板后透射波形成的；而模型B既沒有給物質(zhì)點(diǎn)膨脹擴(kuò)散提供網(wǎng)格空間，也沒有施加物質(zhì)點(diǎn)自由流出邊界條件，這導(dǎo)致模型B中水下爆炸形成的正向沖擊波以及鋼板后方的透射波在傳播至背景網(wǎng)格邊界后不能被自由釋放，并且在網(wǎng)格邊界產(chǎn)生了反射波，因此模型B在這2個(gè)觀測點(diǎn)位置的壓力時(shí)程曲線都具有2個(gè)波峰，其中第2個(gè)波峰都是由網(wǎng)格邊界產(chǎn)生的反射波所造成的，這與模型C的壓力時(shí)程曲線存在顯著差異；而模型A雖然沒有給物質(zhì)點(diǎn)膨脹擴(kuò)散提供網(wǎng)格空間，但是由于在背景網(wǎng)格邊界施加了物質(zhì)點(diǎn)自由流出邊界條件，所以模型A中水下爆炸形成的正向沖擊波以及鋼板后方的透射波都能夠從網(wǎng)格邊界自由釋放，不會在網(wǎng)格邊界產(chǎn)生反射波，因此模型A能夠得到與模型C十分相近的數(shù)值計(jì)算結(jié)果，從而使模型A和模型C在這2個(gè)觀測點(diǎn)位置的壓力時(shí)程曲線高度重合。

圖2 3個(gè)模型中觀測點(diǎn)1和觀測點(diǎn)4的壓力時(shí)程曲線Fig.2 Pressure history of observation point 1 and 4 in 3 models

再通過觀察圖3所示的3個(gè)模型第300 μs的壓力云圖可知，由于模型C給物質(zhì)點(diǎn)膨脹擴(kuò)散提供了網(wǎng)格空間，因此在模型C的壓力云圖中，能夠清楚的觀察到物質(zhì)點(diǎn)區(qū)域的膨脹變形過程和壓力釋放過程；而模型B既沒有給物質(zhì)點(diǎn)膨脹擴(kuò)散提供網(wǎng)格空間，也沒有施加物質(zhì)點(diǎn)自由流出邊界條件，因此模型B水下爆炸形成的正向沖擊波以及鋼板后方的透射波在傳播至背景網(wǎng)格邊界后會出現(xiàn)明顯的反射現(xiàn)象，反射波壓力峰值會達(dá)到0.3 GPa以上；而模型A雖然也沒有給物質(zhì)點(diǎn)膨脹擴(kuò)散提供網(wǎng)格空間，但是由于在背景網(wǎng)格邊界施加了物質(zhì)點(diǎn)自由流出邊界條件，致使模型A中的部分物質(zhì)點(diǎn)及其物理信息能夠從網(wǎng)格邊界自由流出，所以模型A水下爆炸形成的正向沖擊波以及鋼板后方的透射波在傳播至網(wǎng)格邊界后依然能被成功釋放，因此模型A的物質(zhì)點(diǎn)區(qū)域內(nèi)的壓力分布與變化情況能夠與模型C基本保持一致。

通過以上對比分析結(jié)果可知，由于模型B未施加物質(zhì)點(diǎn)自由流出邊界條件，因此模型B的背景網(wǎng)格邊界是一種非嚴(yán)格封閉邊界，物質(zhì)點(diǎn)及其所攜帶的物理信息不能從網(wǎng)格邊界自由流出。而模型A在背景網(wǎng)格邊界施加了物質(zhì)點(diǎn)自由流出邊界條件，因此物質(zhì)點(diǎn)及其物理信息能夠從網(wǎng)格邊界自由流出，從而使模型A能夠得到與模型C非常相近的數(shù)值計(jì)算結(jié)果，并且由于模型A的網(wǎng)格節(jié)點(diǎn)數(shù)量較少，物質(zhì)點(diǎn)又能夠在計(jì)算過程中從網(wǎng)格邊界大量流出，因此模型A的計(jì)算效率較模型C有明顯提高，計(jì)算結(jié)果顯示，模型A的計(jì)算耗時(shí)約為32 869 s，模型C的計(jì)算耗時(shí)約為50 316 s，模型A的計(jì)算耗時(shí)僅為模型C的0.653倍。由此可知，施加物質(zhì)點(diǎn)自由流出邊界條件的物質(zhì)點(diǎn)模型能夠在保證計(jì)算精確度的同時(shí)大幅度提高計(jì)算效率，因此物質(zhì)點(diǎn)自由流出邊界條件是十分有效的。

此外，為了進(jìn)一步觀察近物面水下爆炸模型中正向沖擊波以及鋼板前后的反射波和透射波的基本形態(tài)，現(xiàn)分別提取模型A中距離鋼板前壁面8 cm的觀測點(diǎn)2和距離鋼板后壁面8 cm的觀測點(diǎn)3的壓力時(shí)程數(shù)據(jù)，繪制出圖4所示的這2個(gè)觀測點(diǎn)的壓力時(shí)程曲線。

圖3 3個(gè)模型第300 μs的壓力云圖Fig.3 Pressure contours regarding 3 models at the 300th μs

通過觀察圖4可知，觀測點(diǎn)2的壓力時(shí)程曲線具有2個(gè)峰值超過0.7 GPa的波峰，這2個(gè)波峰分別是由水下爆炸的正向沖擊波和鋼板前壁面反射波所造成的；而觀測點(diǎn)3的壓力時(shí)程曲線則只具有1個(gè)峰值超過0.25 GPa的波峰，該波峰是由鋼板后壁面透射波所造成的。通過觀察圖4還可知，觀測點(diǎn)2位置的正向沖擊波和板前反射波的峰值較高，但正壓作用時(shí)間較短；而觀測點(diǎn)3位置的板后透射波的峰值較低，但正壓作用時(shí)間較長，這與近物面水下爆炸的真實(shí)情況相符合。

下面分別選取圖5所示的模型A第130 μs和第170 μs的壓力云圖來對該模型中正向沖擊波以及鋼板前后的反射波和透射波的基本形態(tài)做更為直觀的觀察。

通過觀察圖5可知，在模型A中，水下爆炸所形成的正向沖擊波作用到鋼板后，會在鋼板的前壁面和后壁面分別產(chǎn)生反射波和透射波，且正向沖擊波和反射波的壓力峰值要高于透射波的壓力峰值，同時(shí)鋼板會發(fā)生明顯變形，這也符合近物面水下爆炸的真實(shí)情況。

圖5 模型A第130 μs和第170 μs的壓力云圖Fig.5 Pressure contours regarding model A at the 130th μs and the 170th μs

2.2 改進(jìn)物質(zhì)點(diǎn)自適應(yīng)分裂方案

2.2.1 物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案的作用與原理

在標(biāo)準(zhǔn)物質(zhì)點(diǎn)法中，物質(zhì)點(diǎn)數(shù)量在數(shù)值計(jì)算過程中是保持不變的，這會導(dǎo)致在計(jì)算對象沿某一方向膨脹或拉伸變形的過程中，該方向的物質(zhì)點(diǎn)也會隨之變得稀疏，這種情況會對計(jì)算結(jié)果產(chǎn)生一定程度的影響，甚至?xí)箶?shù)值計(jì)算產(chǎn)生嚴(yán)重的拉伸不穩(wěn)定現(xiàn)象，這是物質(zhì)點(diǎn)法所面臨的另一種局限[10,12]。為了解決這一問題，Ma等[12]設(shè)計(jì)并提出了物質(zhì)點(diǎn)自適應(yīng)分裂方案(adaptive particle splitting scheme)，但是該方案只能使物質(zhì)點(diǎn)在每次時(shí)間積分計(jì)算中沿某一個(gè)坐標(biāo)軸方向進(jìn)行自適應(yīng)分裂，而不能使物質(zhì)點(diǎn)在每次時(shí)間積分計(jì)算中沿多個(gè)坐標(biāo)軸方向同時(shí)進(jìn)行自適應(yīng)分裂，所以基于該方案的物質(zhì)點(diǎn)法在多維度物質(zhì)點(diǎn)模型數(shù)值計(jì)算中仍然具有局限性。為了徹底破除物質(zhì)點(diǎn)法的這一局限，本文對物質(zhì)點(diǎn)自適應(yīng)分裂方案進(jìn)行了改進(jìn)，使物質(zhì)點(diǎn)在每次顯式時(shí)間積分計(jì)算中能夠沿多個(gè)坐標(biāo)軸方向同時(shí)進(jìn)行自適應(yīng)分裂，以滿足多維度物質(zhì)點(diǎn)模型的數(shù)值計(jì)算要求，并由此形成物質(zhì)點(diǎn)自適應(yīng)分裂的改進(jìn)方案。下面以二維問題為例，來對物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案的基本原理進(jìn)行詳細(xì)介紹。

在由直角坐標(biāo)系所表征的二維問題中，計(jì)算對象具有沿x軸方向和y軸方向進(jìn)行膨脹或拉伸變形的可能性，因此在物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案中，應(yīng)將物質(zhì)點(diǎn)自適應(yīng)分裂分為以下3種情況來分別進(jìn)行討論。

1)只沿x軸方向進(jìn)行分裂。

當(dāng)物質(zhì)點(diǎn)沿x軸方向發(fā)生膨脹或拉伸變形時(shí)，如果物質(zhì)點(diǎn)在x軸方向上的特征尺寸達(dá)到分裂判定條件：

(14)

(15)

2)只沿y軸方向進(jìn)行分裂。

當(dāng)物質(zhì)點(diǎn)沿y軸方向發(fā)生膨脹或拉伸變形時(shí)，如果物質(zhì)點(diǎn)在y軸方向上的特征尺寸達(dá)到分裂判定條件：

(16)

3)同時(shí)沿x軸方向和y軸方向進(jìn)行分裂。

當(dāng)物質(zhì)點(diǎn)沿x軸方向和y軸方向發(fā)生膨脹或拉伸變形時(shí)，如果物質(zhì)點(diǎn)在這2個(gè)方向上的特征尺寸同時(shí)達(dá)到式(14)和式(16)所示的分裂判定條件，則使物質(zhì)點(diǎn)沿這2個(gè)方向同時(shí)進(jìn)行均等分裂，分裂以后將產(chǎn)生4個(gè)新物質(zhì)點(diǎn)，這4個(gè)新物質(zhì)點(diǎn)的各項(xiàng)廣延量的數(shù)值為原物質(zhì)點(diǎn)的0.25倍，而各項(xiàng)強(qiáng)度量的數(shù)值則保持不變，因此新物質(zhì)點(diǎn)的各項(xiàng)物理量的取值方式為：

(17)

以上即為二維條件下的物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案的基本原理。而三維條件下的物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案的基本原理也與之類似，只是需要考慮物質(zhì)點(diǎn)沿z軸方向進(jìn)行自適應(yīng)分裂的情況，因此存在7種分裂方式，此處不再逐一贅述。

2.2.2 算例分析

為了檢驗(yàn)物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案的有效性，本文將該改進(jìn)方案應(yīng)用到聚能射流算例中，并通過算例的對比分析結(jié)果來對該改進(jìn)方案的實(shí)際效果進(jìn)行量化描述。在聚能射流算例中，同樣采用自編的物質(zhì)點(diǎn)法程序構(gòu)建了圖6所示的聚能射流模型(本算例中簡稱物質(zhì)點(diǎn)聚能射流模型)，該模型的背景網(wǎng)格區(qū)域的面積為25.0 cm×10.0 cm，背景網(wǎng)格區(qū)域中構(gòu)建了藥型罩錐角為60°的射流彈，射流彈的藥型罩的材料為銅，殼體材料為鋼，內(nèi)裝藥為PBX9501炸藥，背景網(wǎng)格區(qū)域的其他部位用空氣填充，射流彈各結(jié)構(gòu)的物質(zhì)點(diǎn)的尺寸為0.05 cm×0.05 cm，空氣物質(zhì)點(diǎn)的尺寸為0.1 cm×0.1 cm，網(wǎng)格單元的尺寸為0.1 cm×0.1 cm，整個(gè)模型共離散出29 440個(gè)物質(zhì)點(diǎn)和25 000個(gè)網(wǎng)格單元，背景網(wǎng)格邊界施加了物質(zhì)點(diǎn)自由流出邊界條件，在背景網(wǎng)格區(qū)域內(nèi)沿著射流路徑依次設(shè)置了6個(gè)觀測點(diǎn)，這6個(gè)觀測點(diǎn)的坐標(biāo)見表2。

圖6 聚能射流模型Fig.6 Shaped charge jet mode

表2 聚能射流模型的觀測點(diǎn)坐標(biāo)

此外，由于Autodyn軟件是目前公認(rèn)的數(shù)值模擬爆炸與沖擊動力學(xué)問題效果最好的商業(yè)軟件之一，因此本算例運(yùn)用Autodyn軟件的Euler-2D Multi-material求解器來構(gòu)建與物質(zhì)點(diǎn)聚能射流模型相對應(yīng)的Autodyn二維軸對稱聚能射流模型(本算例中簡稱Autodyn聚能射流模型)，并且采用Autodyn聚能射流模型的計(jì)算結(jié)果來作為物質(zhì)點(diǎn)聚能射流模型的計(jì)算結(jié)果的正確性檢驗(yàn)標(biāo)準(zhǔn)。Autodyn聚能射流模型的Euler網(wǎng)格單元尺寸為0.025 cm×0.025 cm，網(wǎng)格邊界添加Autodyn自帶的流出邊界條件(Flow_Out邊界條件)。

然后分別對物質(zhì)點(diǎn)聚能射流模型和Autodyn聚能射流模型進(jìn)行數(shù)值計(jì)算，其中，物質(zhì)點(diǎn)聚能射流模型的數(shù)值計(jì)算分別在啟動物質(zhì)點(diǎn)自適應(yīng)分裂子程序和不啟動該子程序這2種條件下進(jìn)行，以檢驗(yàn)物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案的有效性。數(shù)值計(jì)算完成后，得到圖7所示的第35 μs模型外觀形態(tài)對比結(jié)果以及表3所示的聚能射流頭部到達(dá)6個(gè)觀測點(diǎn)位置時(shí)頭部速度的對比結(jié)果，且為了便于觀察聚能射流的外觀形態(tài)，空氣介質(zhì)在圖7的3張圖片中均不予顯示。

圖7 3個(gè)模型第35 μs的外觀形態(tài)對比結(jié)果Fig.7 Comparative results of the visual forms regarding 3 models at the 35th μs

首先通過觀察圖7可知，在未啟動物質(zhì)點(diǎn)自適應(yīng)分裂子程序的物質(zhì)點(diǎn)聚能射流模型數(shù)值計(jì)算中，隨著內(nèi)裝藥爆轟產(chǎn)物的持續(xù)擴(kuò)散和銅質(zhì)藥型罩的持續(xù)擠壓拉伸變形，內(nèi)裝藥爆轟產(chǎn)物物質(zhì)點(diǎn)和藥型罩物質(zhì)點(diǎn)在背景網(wǎng)格區(qū)域內(nèi)的空間分布密度逐漸降低，數(shù)值計(jì)算過程產(chǎn)生了嚴(yán)重的拉伸不穩(wěn)定現(xiàn)象，以至于當(dāng)物質(zhì)點(diǎn)聚能射流模型運(yùn)行至35 μs時(shí)，內(nèi)裝藥爆轟產(chǎn)物區(qū)顯露出大面積空白區(qū)域，并且藥型罩物質(zhì)點(diǎn)之間也發(fā)生了明顯的非物理性分離，因此所形成的聚能射流并不連續(xù)，這與Autodyn聚能射流模型的聚能射流外觀形態(tài)存在本質(zhì)性差別；而在啟動物質(zhì)點(diǎn)自適應(yīng)分裂子程序的物質(zhì)點(diǎn)聚能射流模型數(shù)值計(jì)算中，隨著內(nèi)裝藥爆轟產(chǎn)物的持續(xù)擴(kuò)散和藥型罩的持續(xù)擠壓拉伸變形，內(nèi)裝藥爆轟產(chǎn)物物質(zhì)點(diǎn)和藥型罩物質(zhì)點(diǎn)在x軸方向和y軸方向上不斷進(jìn)行自適應(yīng)分裂，從而可保證物質(zhì)點(diǎn)在背景網(wǎng)格區(qū)域內(nèi)的空間分布密度不會明顯降低，有效避免了拉伸不穩(wěn)定現(xiàn)象的產(chǎn)生，因此當(dāng)物質(zhì)點(diǎn)聚能射流模型運(yùn)行至35 μs時(shí)，共產(chǎn)生出84 600個(gè)物質(zhì)點(diǎn)(其中包括29 440個(gè)初始物質(zhì)點(diǎn))，此時(shí)的物質(zhì)點(diǎn)總數(shù)約是初始物質(zhì)點(diǎn)總數(shù)的2.87倍，內(nèi)裝藥爆轟產(chǎn)物區(qū)沒有出現(xiàn)大面積空白區(qū)域，并且藥型罩物質(zhì)點(diǎn)之間也沒有發(fā)生非物理性分離，因此所形成的聚能射流的連續(xù)性非常好，這與Autodyn聚能射流模型的聚能射流外觀形態(tài)具有較高的吻合度。

再通過表3可知，物質(zhì)點(diǎn)聚能射流模型在未啟動物質(zhì)點(diǎn)自適應(yīng)分裂子程序的條件下，其與Autodyn聚能射流模型只是在觀測點(diǎn)1位置的射流頭部速度相對差值的絕對值低于4%，而在其余5個(gè)觀測點(diǎn)位置的射流頭部速度相對差值的絕對值均高于4%；而物質(zhì)點(diǎn)聚能射流模型在啟動物質(zhì)點(diǎn)自適應(yīng)分裂子程序的條件下，其與Autodyn聚能射流模型在6個(gè)觀測點(diǎn)位置的射流頭部速度相對差值的絕對值均低于2.5%。由此可知，啟動物質(zhì)點(diǎn)自適應(yīng)分裂子程序條件下的物質(zhì)點(diǎn)聚能射流模型與Autodyn聚能射流模型的射流頭部速度相差較小，即該條件下的物質(zhì)點(diǎn)聚能射流模型的計(jì)算精確度更高。

表3 3個(gè)模型在6個(gè)觀測點(diǎn)位置的射流頭部速度對比結(jié)果Table 3 Comparative results of shaped charge jet tip velocity by the location of 6 observation points regarding 3 models

通過本算例的上述對比分析結(jié)果可知，物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案能夠使物質(zhì)點(diǎn)在每次顯式時(shí)間積分計(jì)算中沿多個(gè)坐標(biāo)軸方向同時(shí)進(jìn)行自適應(yīng)分裂，從而可避免物質(zhì)點(diǎn)模型的數(shù)值計(jì)算產(chǎn)生拉伸不穩(wěn)定現(xiàn)象，進(jìn)而使數(shù)值計(jì)算的精確度得以顯著提高。由此可知，物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案對物質(zhì)點(diǎn)法的改進(jìn)效果是十分顯著的。

3 結(jié)論

1) 物質(zhì)點(diǎn)自由流出邊界條件能夠使物質(zhì)點(diǎn)及其所攜帶的物理信息在每次顯式時(shí)間積分計(jì)算中根據(jù)實(shí)際情況從背景網(wǎng)格邊界自由流出，所以構(gòu)建背景網(wǎng)格時(shí)并不需要給物質(zhì)點(diǎn)膨脹擴(kuò)散提供額外的網(wǎng)格空間，因此可使物質(zhì)點(diǎn)模型在保證計(jì)算精確度的同時(shí)大幅度提高計(jì)算效率；

2) 物質(zhì)點(diǎn)自適應(yīng)分裂改進(jìn)方案能夠使物質(zhì)點(diǎn)在每次顯式時(shí)間積分計(jì)算中根據(jù)分裂判定條件沿多個(gè)坐標(biāo)軸方向同時(shí)進(jìn)行自適應(yīng)分裂，從而有效避免了數(shù)值計(jì)算過程中拉伸不穩(wěn)定現(xiàn)象的發(fā)生，使多維度物質(zhì)點(diǎn)模型的計(jì)算精確度得以顯著提高；

3) 經(jīng)過這2項(xiàng)改進(jìn)的物質(zhì)點(diǎn)法能夠更加適用于求解爆炸與沖擊類工程實(shí)際問題。