DSPH方法的有效性驗證及應(yīng)用＊

閆 蕊，徐 緋，張岳青

（1.西北工業(yè)大學(xué)航空學(xué)院，陜西 西安 710072；2.北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點實驗室，北京 100081）

光滑粒子流體動力學(xué)（smoothed particle hydrodynamics，SPH）方法是一種典型的無網(wǎng)格拉格朗日粒子法。SPH方法思想簡單，容易實現(xiàn)，一經(jīng)出現(xiàn)便得到了迅速發(fā)展。相比較其他的無網(wǎng)格方法，它最大的優(yōu)點是便于在大變形實際工程中應(yīng)用，這一點是其他無網(wǎng)格方法無法比擬的［1］。在使用SPH進行模擬計算時經(jīng)常會遇到界面處參數(shù)不連續(xù)的問題，通常SPH方法處理界面處的問題時，與連續(xù)體內(nèi)部處理方法相同，會嚴重影響算法的準確性。M.B.Liu等［2］提出了1D不連續(xù)的SPH方法（DSPH），用來模擬不連續(xù)物理量，函數(shù)測試和一維激波管的模擬均顯示了該算法的良好性能。趙燕等［3］推導(dǎo)了2D和3D不連續(xù)公式，并且解決了多維情況下確定不連續(xù)位置困難的問題。本文的主要工作是分析界面處該方法提高計算精度的本質(zhì)原因，探討DSPH方法如何應(yīng)用于實際的物理問題中。

1 DSPH介紹

文獻［3］基于Taylor公式推導(dǎo)了2D和3D的不連續(xù)SPH公式，給出了多維情況下不連續(xù)函數(shù)的粒子近似式和不連續(xù)函數(shù)一階導(dǎo)數(shù)的粒子近似式

式中：W 表示核函數(shù)［1］，mj、ρj分別為粒子j的質(zhì)量及密度，α和β表示空間的維數(shù)。Ω為問題求解域，Ω1和Ω2為Ω的子域（Ω＝Ω1＋Ω2），函數(shù)或物理參數(shù)在Ω1和Ω2上分別連續(xù)，xi為Ω1上靠近不連續(xù)界面的一點，xk是Ω2上的任意一點，具體介紹見文獻［3］，其中

文獻［3］中針對越過材料界面不連續(xù)物理量的計算，給出了大變形計算中確定不連續(xù)位置的方法，即根據(jù)顆粒材料屬性是否相同確定是否為界面處。根據(jù)泰勒展開公式，從理論上給出了確定不連續(xù)公式中點xk的方法，并用數(shù)值方法驗證了此方法的有效性。本文中將從理論上來探究該方法的有效性。

2 DSPH方法函數(shù)值計算的精確度探討

文獻［3］中已經(jīng)證明了DSPH在模擬不連續(xù)問題時相較于SPH和CSPM所表現(xiàn)出來的優(yōu)越性，同時也給出了3種不同粒子分布下DSPH所產(chǎn)生的誤差大小。本文中將進一步針對粒子的均勻性對DSPH的影響進行研究。

粒子的分布情況如圖1所示，粒子間距為1，粒子的自變量m的初始值為粒子下方的編號，核函數(shù)為三次樣條函數(shù)，核函數(shù)的光滑半徑取為h＝1.5。函數(shù)形式為

粒子11移動距離為x，－0.5≤x＜0.5；粒子12移動距離為y，－0.5≤y＜1。其他粒子均保持不動，改變x、y的值可以得到不同粒子分布情況下f（10）的函數(shù)值，不論粒子怎樣移動，用DSPH方法所計算的f（10）的函數(shù)值的誤差均為0。這說明了粒子11、12的值及分布情況對粒子10的計算結(jié)果不產(chǎn)生任何影響。

圖1 粒子分布示意圖Fig.1 The distribution of particles

下面探討DSPH方法能夠精確估算函數(shù)值的本質(zhì)原因，由于式（1）中2部分的分母相同，因此只考慮分子的影響。將式（1）中2部分中的Ω寫成Ω1和Ω2的和，結(jié)合文獻［3］確定xk的方法可將式（1）寫成如式下形式

由式（8）可以看出，DSPH對函數(shù)值的估計的本質(zhì)是：粒子i的函數(shù)值與它不連續(xù)界面異側(cè)的粒子函數(shù)值無關(guān)，它的函數(shù)值取決于i同側(cè)的粒子函數(shù)值，并且加大了粒子i本身在計算中所占的權(quán)重。Ω2中粒子的質(zhì)量、密度及分布情況部分地決定了粒子i在估算點i處函數(shù)值時所占的權(quán)重。當使用DSPH方法處理邊界問題時，Ω2中無粒子存在，DSPH方法與CSPM方法相同［1］。

在式（7）中，左側(cè)粒子的函數(shù)若不是常函數(shù)，計算得到的f（10）的函數(shù)值則可能是不精確的，誤差主要是由于SPH方法本身的特點引起的，因此在此沒有必要討論左側(cè)粒子函數(shù)的變化對粒子10的影響。當搜索域足夠小，搜索域內(nèi)與被估粒子同側(cè)的粒子函數(shù)可以認為是常函數(shù)時，可以保證被估粒子計算的精確性。

3 DSPH方法導(dǎo)數(shù)值計算的精確度探討

與第2節(jié)類似，改變粒子的分布情況，檢測DSPH在計算導(dǎo)數(shù)值時對粒子分布均勻性的敏感度，與圖1相同，改變粒子11和12的位置，函數(shù)形式為

估算點10處的導(dǎo)數(shù)值，經(jīng)過計算我們發(fā)現(xiàn)任意改變粒子11和粒子12的位置都不影響導(dǎo)數(shù)值計算的精確性。

下面通過DSPH方法的導(dǎo)數(shù)公式來分析產(chǎn)生這種現(xiàn)象的原因，與第2節(jié)類似，式（2）可改寫為

在式（10）中，令Fi＝xif′i＋fi，F(xiàn)j＝xjf′i＋fi，則式（10）可寫為

式（11）中的第2項Fi的含義是以點i處導(dǎo)數(shù)值為斜率，以fi為截距虛擬出一條直線來，在新的函數(shù)形式下估算粒子i的導(dǎo)數(shù)值，粒子i導(dǎo)數(shù)值的估算與異側(cè)粒子的原函數(shù)形式?jīng)]有關(guān)聯(lián)。因此當與被估粒子同側(cè)的函數(shù)為線性函數(shù)時，計算得到的被估粒子的導(dǎo)數(shù)值是精確的。由此不難推得，當不連續(xù)函數(shù)為多次函數(shù)時，只要粒子的搜索域足夠小，搜索域內(nèi)與被估粒子同側(cè)的粒子函數(shù)可以認為是線性函數(shù)時，可以保證被估粒子導(dǎo)數(shù)值計算的精確性。

4 物理計算中不連續(xù)參數(shù)分析

SPH方法是一種拉格朗日方法，總的來說拉格朗日描述下的流體控制方程可寫作一系列偏微分方程，即著名的Navier－Stokes方程，用SPH公式解Navier－Stokes方程組就可以得到SPH方法中的密度、應(yīng)變率、加速度和能量的表達式。實施SPH方法的過程是通過應(yīng)用局部區(qū)域內(nèi)的相鄰粒子對應(yīng)的值來疊加求和取代場函數(shù)及其導(dǎo)數(shù)的積分表示形式［1］。從加速度計算的推導(dǎo)過程［4］可知，計算加速度需要考慮作用在該點的力，即使是在界面處，只要有相互作用，就需要考慮作用力，因此僅考慮在計算密度、應(yīng)變率和能量時使用DSPH方法。由Navier－Stokes方程的表達式可知這3個物理量之所以存在不連續(xù)問題，本質(zhì)上都是由于界面處速度的不連續(xù)引起的。根據(jù)式（2）可以得到速度梯度的DSPH表達式

式中：Bi，βγ的含義與式（3）類似。

將式（12）帶入Navier－Stokes方程就可以得到相應(yīng)的密度、應(yīng)變率和能量的DSPH表達式。

4.1 密度

計算模型如圖2所示，鋁塊的初速度為零，粒子間距為5mm，總粒子數(shù)為400。對右側(cè)鋁塊的最右側(cè)一排粒子施加水平方向上的約束。Mie－Grüneisen狀態(tài)方程參數(shù)和Johnson－Cook材料屈服模型參數(shù)均來自文獻［5－6］。各個參數(shù)有如圖3所示的相互影響關(guān)系，在計算中加入了了人工黏性，系數(shù)為0.6、0.6、0.1。

鋼塊的初速度為140m／s，時間步長t＝50ns。圖4中給出了SPH方法和DSPH方法計算得到的1400步時物體的變形情況，SPH的計算結(jié)果中，鋼塊在撞擊鋁塊后粘在了鋁塊上左右震動，如圖4（a）所示，鋁塊的粒子出現(xiàn)斷裂并且鋼鋁界面處粒子雜亂。在DSPH的計算結(jié)果中，鋼塊在撞擊鋁塊后被反彈，兩物塊在界面處變形較小，如圖4（b）所示。圖5中給出了1 400步時左側(cè)鋼塊的密度分布，粒子編號規(guī)律為：先對左側(cè)物體編號，自底部粒子開始，最右側(cè)一排粒子編號為1～10，右起第2排粒子編號為11～20，以此類推?？梢钥闯鼍嚯x界面越近，密度的計算值越不準確。這說明在界面處傳統(tǒng)的SPH方法處理參數(shù)不連續(xù)時所存在的缺陷，同時也說明邊界處理不好，會嚴重影響非邊界處粒子物理量的計算精度，與之相對應(yīng)的DSPH方法所計算出的密度值非常精確，改善了模擬結(jié)果。

圖2 計算模型Fig.2 The model of simulation

圖3 各參數(shù)之間的相互關(guān)系Fig.3 The relationship between the parameters

圖4 1400步時物體的變形情況Fig.4 The deformation of targets in 1 400th step

為了進一步說明DSPH方法計算結(jié)果的準確性，將左側(cè)鋼塊的初速度提高到300m／s，時間步長t＝5ns。圖6是1 200步時SPH和DSPH的計算結(jié)果，圖7是分別使用SPH方法和DSPH方法得到的點A（位置如圖6所示）處的應(yīng)力應(yīng)變曲線，將2種計算結(jié)果與理論值相比較，可以看出DSPH很好地改善了界面處的模擬精度。

圖5 1400步時鋼塊上的密度分布Fig.5 The density in different particles of steel in 1 400th step

圖6 1200步時物體的變形情況Fig.6 The deformation of targets in 1 200th step

連續(xù)介質(zhì)力學(xué)中連續(xù)性方程的推導(dǎo)過程，可反映出對于固體，若界面處存在速度的不連續(xù)，則需要在計算密度時使用DSPH方法［4］，這一結(jié)論同樣適用于同種材料的界面問題。將圖6中的算例改為2個鋼塊對碰，輸出各個粒子的密度值以及點A處的應(yīng)力應(yīng)變曲線，可以發(fā)現(xiàn)，在此類問題中，DSPH方法的計算結(jié)果同樣優(yōu)于SPH方法。當界面兩側(cè)不存在相對速度時，無論是同種材料還是不同材料物體的界面問題，SPH方法都能較為精確地計算密度值，就沒有必要使用DSPH方法了，對此我們也做了相應(yīng)的模擬，2種方法計算得到的結(jié)果相同。

當然以上理論僅針對固體而言，對于等體積內(nèi)質(zhì)量會發(fā)生變化的問題而言，比如不同氣體的滲入、液體混合等就不推薦使用DSPH方法來計算密度。

SPH中計算密度的另一種方法是密度求和法［1］，密度求和法體現(xiàn)了SPH近似法的本質(zhì)，但在實際的程序計算中，該方法計算效果很差。有許多改進該算法的修正方案，比如目前較常用的CSPM方法，但是CSPM方法在求解含間斷算例時卻存在計算不穩(wěn)定的問題［7］，若使用該方法計算密度，則需限定接觸的兩個物體為同種材料。使用DSPH方法可以擺脫這一限制，根據(jù)式（1）可得到密度求和法的DSPH表達式

式（13）右邊的第1項與CSPM方法相同，第2項是通過對密度非連續(xù)性進行處理得到的。使用與圖4相同的計算模型，從圖8中可以看到，CSPM方法的模擬結(jié)果很不理想，鋼塊在界面處收縮，而鋁塊在界面處擴張，最終由于變形過大而導(dǎo)致計算不能繼續(xù)進行。而DSPH方法的計算結(jié)果與圖4（b）相同，密度計算值精確度與連續(xù)性密度法（DSPH）基本相同。這證明了DSPH方法在處理界面不連續(xù)問題時優(yōu)于CSPM方法。

圖7 點A處的應(yīng)力應(yīng)變曲線Fig.7 The stress－strain curve in particleA

圖8 兩物體變形情況Fig.8 The deformation of targets

4.2 應(yīng)變率

由圖3可知，應(yīng)變率的計算主要是影響了偏應(yīng)力和屈服應(yīng)力的計算，為了說明應(yīng)變率在使用DSPH方法后對結(jié)果的影響，左側(cè)的鋼塊以初速20m／s的速度正碰右側(cè)的鋼塊（加速度是由各向同性壓力和偏應(yīng)力共同影響的［1］，速度過大時物體會產(chǎn)生變形，偏應(yīng)力相較于各向同性壓力較小，不利于說明應(yīng)變率改變對結(jié)果的影響），計算模型如圖2所示，時間步長t＝50ns，未使用人工黏性。由同種材料物體的碰撞規(guī)律可知，在撞擊后左側(cè)鋼塊靜止。圖9是計算得到的左側(cè)鋼塊速度隨時間的變化曲線，可以看出應(yīng)變率的DSPH方法對模擬結(jié)果的改善。

并不是在所有界面處都可以使用DSPH方法計算應(yīng)變率，當界面處速度連續(xù)，比如兩種材料組成的組合梁，在梁縱向?qū)ΨQ面內(nèi)作用一對力偶矩，使梁產(chǎn)生彎曲變形，此時界面處應(yīng)力不連續(xù)，但應(yīng)變是連續(xù)的［8］，使用DSPH方法計算應(yīng)變率后會導(dǎo)致界面處變形不一致，兩材料在界面處脫離，傳統(tǒng)的SPH方法保證了變形的連續(xù)性。

圖9 左側(cè)鋼塊速度隨時間步的變化情況Fig.9 The velocity of steel in different time steps

4.3 能量

由圖3可知，能量的計算只影響后續(xù)兩個物理量的計算：各向同性壓力p及屈服強度Y。對于固體來說密度值為常數(shù)，使用 Mie－Grüneisen狀態(tài)方程計算各向同性壓力p，由p的表達式［1］可知若密度計算不準確，將影響p與能量之間的線性關(guān)系。所以考慮在對密度使用DSPH方法的基礎(chǔ)上，討論對能量使用DSPH方法對結(jié)果的影響。計算模型與圖6相同，速度為300m／s，輸出圖6中點A處的能量變化可得圖10，兩曲線分別為對密度使用DSPH方法（d）和對密度、能量同時使用DSPH方法（ed）得到的，兩者有一定的差異。圖11中是兩種方案得到的各向同性壓力p和屈服強度Y的曲線圖，可以看到與p的表達式所顯示的一樣，各向同性壓力p與能量的變化趨勢相同，屈服強度Y的差別較小。提高撞擊速度，兩種算法計算出的能量差異會更大，相應(yīng)的各向同性壓力和屈服強度的差異也較大，但是對加速度并沒有影響，同一時刻兩種方法得到的變形圖基本相同，因此我們認為在固體沖擊問題的計算中沒有必要對能量使用DSPH方法進行計算。

圖10 能量隨時間步的變化Fig.10 The energy in different time steps

圖11 各向同性壓力p和屈服強度Y隨時間步的變化Fig.11 The isotropic pressure pand yield strength Y in different time steps

5 結(jié) 論

通過對DSPH計算公式的推導(dǎo)分析，發(fā)現(xiàn)DSPH方法之所以能夠很好地處理界面問題，本質(zhì)原因在于DSPH方法計算函數(shù)值及導(dǎo)數(shù)值時忽略了在不連續(xù)界面異側(cè)的粒子對被估粒子的影響，按照同等比例加大了被估粒子本身對結(jié)果的影響。

通過對碰模型考察了DSPH方法在計算密度、應(yīng)變率及能量這3個物理量時所產(chǎn)生的影響，三者的實質(zhì)都是考慮了界面處的速度不連續(xù)。密度計算的準確性在整個計算中起著至關(guān)重要的作用；在偏應(yīng)力相較于各向同性壓力不可忽略的問題中使用DSPH方法計算應(yīng)變率會對結(jié)果產(chǎn)生一定的影響；DSPH方法對能量計算的影響并不大。在計算中密度和應(yīng)變率在使用DSPH方法時是相互獨立的，可同時使用，共同提高計算精度。

在本文的模擬中只是使用了簡單的對碰模型，今后需要將此方法應(yīng)用于更為復(fù)雜的物理問題中，以進一步驗證該方法的有效性，更好地發(fā)揮DSPH解決邊界問題的優(yōu)勢。

［1］韓旭，楊剛，強洪夫.光滑粒子流體動力學(xué)—— 一種無網(wǎng)格粒子法［M］.長沙：湖南大學(xué)出版社，2005：35－323.

［2］Liu M B，Liu G R，Lam K Y.A one－dimensional meshfree particle formulation for simulating shock waves［J］.Shock Waves，2003，13（3）：201－211.

［3］趙燕，徐緋，李玉龍.一種可考慮界面不連續(xù)的改進SPH 方法［J］.計算力學(xué)學(xué)報，2009，26（6）：928－934.Zhao Yan，Xu Fei，Li Yu－long.An improved SPH method in the discontinuous interface problem［J］.Chinese Journal of Computational Mechanics，2009，26（6）：928－934.

［4］黃筑平.連續(xù)介質(zhì)力學(xué)基礎(chǔ)［M］.北京：高等教育出版社，2003：125－143.

［5］Morris J P.Analysis of smoothed particle hydrodynamics with application［D］.Australia：Monash University，1996.

［6］Katayama M，Toda S，Kibe S.Numerical simulation of space debris impacts on the Whipple shield［J］.Acta Astronautica，1997，40（12）：859－869.

［7］Zukas J A.High velocity impact［D］.New York：John Wiley and Sons，1990.

［8］茍文選.材料力學(xué)［M］.北京：科學(xué)出版社，2005：215.