一種考慮界面不連續(xù)的改進的有限粒子法

王 璐，楊 揚，徐 緋

(西北工業(yè)大學(xué)航空學(xué)院，陜西 西安 710072)

光滑粒子流體動力學(xué)(smoothed particle hydrodynamics，SPH)是一種典型的拉格朗日型無網(wǎng)格粒子法，具有概念簡單、自適應(yīng)性強等諸多優(yōu)點，工程應(yīng)用十分廣泛[1-3]。但其發(fā)展至今仍然存在一定的數(shù)值缺陷，其中在粒子缺失較明顯的模型邊界或物理量變化較大的非連續(xù)界面附近表現(xiàn)出的核估計精度不足尤為顯著[4-5]。研究發(fā)現(xiàn)，可以通過提高核函數(shù)的一致性提高邊界區(qū)域的計算精度。為此，Liu等[6]在SPH的基礎(chǔ)上通過重構(gòu)核函數(shù)提出了再生核質(zhì)點法(reproducing kernel particle method，RKPM)，而后Chen等[7]以泰勒級數(shù)展開為基礎(chǔ)提出了一種對SPH估計式進行正則化的方法，即修正光滑粒子法(corrective smoothed particle method，CSPM)，并得到了成功應(yīng)用[8]；隨后， Liu等[9]基于泰勒級數(shù)展開提出了有限粒子法(finite particle method，F(xiàn)PM)，相比傳統(tǒng)的SPH方法，F(xiàn)PM在核函數(shù)的選取上較寬松，邊界區(qū)域計算精度高，高階擴展性強。相比CSPM，F(xiàn)PM可以同時估計核函數(shù)及其導(dǎo)數(shù)值，一定程度上降低了CSPM中低階導(dǎo)數(shù)用于高階導(dǎo)數(shù)計算而產(chǎn)生的累積誤差。然而FPM也存在一些固有缺陷：在處理不連續(xù)問題時，若與連續(xù)體內(nèi)部處理方法一致，算法的準(zhǔn)確性會大大降低；求解線性方程組一定程度上增加了計算復(fù)雜度，計算時間較長；系數(shù)矩陣是否奇異不可控，影響計算穩(wěn)定性。因而使用FPM方法時需對界面進行更精細的處理[10-11]。基于以上缺陷，本文中針對界面不連續(xù)問題對FPM方法進行改進，提高界面處的精度和計算效率。

針對傳統(tǒng)SPH方法，目前已發(fā)展了不同的界面處理方法[12-13]，其中較經(jīng)典的是Liu等[14]基于CSPM提出的非連續(xù)SPH(discontinuous SPH，DSPH)方法，DSPH方法改善了傳統(tǒng)SPH方法在非連續(xù)區(qū)域上的計算精度。本文中基于DSPH方法處理非連續(xù)區(qū)域的思想，從FPM出發(fā)，推導(dǎo)DSFPM(discontinuous special FPM)方程，在理論層面改進FPM算法在界面處的計算精度。

1 DSFPM

1.1 FPM基本理論

根據(jù)Liu等[9]基于泰勒展開和矩陣理論給出的FPM基本方程，經(jīng)過具體推導(dǎo)，給出FPM中一維以及二維情形下函數(shù)及其一階導(dǎo)數(shù)的粒子近似式：

(1)

(2)

式中:f(xi)為點xi處的函數(shù)值，fx(xi)為點xi處的函數(shù)導(dǎo)數(shù)值；類似地，f(xi,yi)、fx(xi,yi)和fy(xi,yi)分別為點(xi,yi)處的函數(shù)值及偏x和偏y導(dǎo)數(shù)；j為粒子i支持域內(nèi)的粒子，Wij為基函數(shù)，Wij,x=?Wij/?x，Wij,y=?Wij/?y分別為基函數(shù)的偏導(dǎo)數(shù)，mj、ρj分別為粒子j的質(zhì)量和密度，N為粒子i完整支持域內(nèi)的粒子數(shù)。

從式(1)～(2)可以看出，式(1)～(2)等號右邊第一項矩陣只和粒子i支持域內(nèi)粒子的位置有關(guān)，與其函數(shù)值無關(guān)，粒子的分布決定了矩陣是否奇異和FPM計算的穩(wěn)定性。

1.2 DSPH基本理論

不同于FPM在粒子i完整支持域內(nèi)進行泰勒展開，DSPH基于界面處物理量的不連續(xù)性，在非連續(xù)區(qū)域兩邊分別進行泰勒展開。Xu等[15]和閆蕊等[16]在Liu等[14]提出的一維不連續(xù)SPH方法的基礎(chǔ)上將DSPH擴展到多維情形，并且避免了多維情形下確定不連續(xù)位置，降低了計算的復(fù)雜性。式(3)～(4)為DSPH粒子估計式，如圖1 所示，Ω為粒子的整個支持域，Ω1和Ω2為其子域，Ω1和Ω2之間不連續(xù)界面為Γ(Ω1∪Ω2∪Γ=Ω)。

圖1 界面附近粒子支持域Fig.1 The supported domain of the particle near the interface

(3)

(4)

式中:fα(xi)=?f(xi)/?xα，α、β表示坐標(biāo)方向(一維為x,二維為x、y)。方程(3)和(4)等號右邊第一項與CSPM保持一致，第二項是考慮了不連續(xù)界面后產(chǎn)生的附加項，附加項中只提取了異側(cè)粒子的位置、質(zhì)量和密度信息，方程式很明顯為隱式格式，保證了計算的穩(wěn)定性但較耗時。

1.3 DSFPM推導(dǎo)

從式(1)可以看到，當(dāng)遇到非連續(xù)物理時，F(xiàn)PM在完整支持域內(nèi)求和的做法對于求解非連續(xù)問題的精度是遠遠不夠的，本節(jié)將基于DSPH方法處理非連續(xù)問題的思想對FPM進行理論改進。

從方程式(3)出發(fā)，對其進行變形推導(dǎo)得到:

(5)

(6)

從式(5)～(6)可以看出DSPH公式經(jīng)過變形后，本質(zhì)是只在子域Ω1內(nèi)對被估粒子進行一致化修正，這樣被估粒子的物理信息就完全與異側(cè)粒子分隔開。基于該想法，對FPM基本方程即式(1)～(2)進行下列修正[17]，得到DFPM(discontinuous FPM)的基本方程:

(7)

(8)

其中，DFPM只考慮了鄰域內(nèi)與粒子i同側(cè)的粒子信息。該修正式既保證了高階導(dǎo)數(shù)的計算精確性，也實現(xiàn)了不連續(xù)界面處粒子信息的準(zhǔn)確估計。但是式(7)～(8)等號右邊第一項矩陣仍不能準(zhǔn)確保證計算的穩(wěn)定性，并且矩陣的各項求解也使得計算過程變得繁瑣費時，因此下面進一步對式(7)～(8)進行變形整理。

從式(7)出發(fā)，將式(7)整理為：

(9)

一維情形下，粒子所占體積為粒子間距dj，區(qū)域Ω1內(nèi)的粒子數(shù)為N1，類似文獻[18]的處理方法，在式(9)中，記：

(10)

(12)

式(9)變形為：

KDCf=KDF

(13)

圖2 一維情形下DSFPM粒子選取方式Fig.2 Particle selection mode for 1-D DSFPM

圖3 二維情形下DSFPM粒子選取方式Fig.3 Particle selection mode for 2-D DSFPM

可以看出，K只與核函數(shù)有關(guān)，D為粒子間距信息，C為區(qū)域Ω1內(nèi)的粒子坐標(biāo)信息，F(xiàn)反映了區(qū)域Ω1內(nèi)的粒子函數(shù)值信息。此時，f為待求未知量。由于FPM中核函數(shù)選取自由，因此我們只考慮K為滿秩的情形，如果在區(qū)域Ω1內(nèi)只選取距離粒子i最近的2個粒子，如圖2所示，紅虛線為界面，此時，K為可逆方陣，D也為可逆方陣。在式(13)兩邊同時乘以D-1K-1，式(13)簡化為：

Cf=F

(14)

(15)

類似地，二維情形下，如果在粒子(xi,yi)的支持域子域Ω1內(nèi)選取距離粒子(xi,yi)最近的3個粒子，如圖3所示，紅虛線為界面，此時K和D亦均為可逆矩陣，同樣可得到式(14)。具體表達式為：

(16)

式(15)～(16)即為DSFPM基本方程，相比于式(7)～(8)，DSFPM消除了核函數(shù)對于矩陣奇異性的影響，提高了計算的穩(wěn)定性，計算過程大大簡化。

2 DSFPM粒子近似精度分析

一般情況下，如果一種近似格式能夠再生k階多項式，那么則稱該格式具有Ck階連續(xù)性。類似地，我們定義如果近似格式能夠再生k階非連續(xù)多項式，那么該格式具有Ck階非連續(xù)性。接下來，通過一系列數(shù)值算例對傳統(tǒng)FPM、DSPH 和DSFPM等3種方法的粒子近似精度進行對比分析。在數(shù)值計算中，選用常用的三次B-樣條函數(shù)作為光滑函數(shù)[19]。

2.1 C0階非連續(xù)性分析

C0階表征近似格式對于非連續(xù)常值函數(shù)的再生能力，考慮區(qū)間[0,1]上的非連續(xù)常值函數(shù)：

(17)

自x=0至x=1均勻分布21個粒子，粒子間距d=0.05，光滑長度h=1.2d。圖4給出了FPM、DSPH和DSFPM等3種方法估計的函數(shù)值誤差φ及其一階導(dǎo)數(shù)誤差ψ,紅虛線為界面。從圖4可以看出，DSPH和DSFPM均可以準(zhǔn)確計算界面附近粒子函數(shù)及其導(dǎo)數(shù)值，F(xiàn)PM在界面處計算精度低。這說明FPM處理不連續(xù)問題有明顯缺陷，而其余2種方法改善了界面處的精度缺陷。

圖4 3種方法估計非連續(xù)常值函數(shù)誤差及其導(dǎo)數(shù)誤差Fig.4 Estimation errors for the discontinuous 0-order function and its derivative by three methods

2.2 C1階非連續(xù)性分析

C1階表征近似格式對于非連續(xù)線性函數(shù)的再生能力，考慮一維不連續(xù)函數(shù)：

(18)

自x=0至x=1均勻分布21個粒子，粒子間距及光滑長度不變。圖5給出FPM、DSPH和DSFPM等3種方法近似的函數(shù)值誤差φ及其一階導(dǎo)數(shù)誤差ψ。

從圖5可以看出：(1)對線性函數(shù)值，F(xiàn)PM和DSPH均不能準(zhǔn)確估計，F(xiàn)PM在界面(x=0.525)附近的誤差值遠大于DSPH相應(yīng)的誤差值。這說明DSPH在一定程度上改善了界面處的計算精度。而DSPH在邊界處(x=0，x=1)不能準(zhǔn)確估計函數(shù)值，這是由于邊界粒子支持域內(nèi)粒子信息缺失導(dǎo)致核函數(shù)估計精度不足。而DSFPM能夠很好地估計界面和邊界處線性函數(shù)值。(2)對線性函數(shù)一階導(dǎo)數(shù)值，F(xiàn)PM在界面附近的計算誤差較大，DSPH對邊界和界面附近的函數(shù)值估計有誤差導(dǎo)致DSPH估計一階導(dǎo)數(shù)時仍不精確，這體現(xiàn)了DSPH利用低階導(dǎo)數(shù)計算高階導(dǎo)數(shù)帶來的誤差累積效應(yīng)，DSFPM仍可以準(zhǔn)確再生一階導(dǎo)數(shù)值，大大改善了界面以及邊界附近的計算精度。DSFPM不僅延續(xù)了FPM在邊界處核估計的高精度性，也繼承了DSPH對于非連續(xù)區(qū)域計算的修正。

圖5 3種方法估計非連續(xù)一次函數(shù)誤差及其導(dǎo)數(shù)誤差Fig.5 Estimation errors for the discontinuous 1-order function and its derivative by three methods

2.3 C2階非連續(xù)性分析

這一部分將驗證DSFPM對二次函數(shù)的再生能力。考慮二維情形下，區(qū)域 [0,1]×[0,1]上的非連續(xù)二次函數(shù)：

(19)

在二維區(qū)域內(nèi)均勻分布441個粒子，光滑長度h=1.2d。圖6為估計二次函數(shù)及其一階偏導(dǎo)(由于界面垂直于y軸,因此只關(guān)注偏y導(dǎo)數(shù)誤差)的誤差。從圖6可以看出，DSFPM不能對二次函數(shù)進行準(zhǔn)確再生。這說明DSFPM只能對一次非連續(xù)函數(shù)進行準(zhǔn)確估計，具有C1階非連續(xù)性。這是因為DSFPM推導(dǎo)過程中利用的FPM基本方程只進行了一階泰勒展開，并未考慮更高階項。如果將f(x,y)泰勒展開到更高階項，那么可以得到更高階精度的DSFPM，這也體現(xiàn)了DSFPM的高階可擴展性。

圖6 DSFPM計算二次非連續(xù)函數(shù)本身和其偏y導(dǎo)數(shù)誤差分布Fig.6 Estimation error distribution for the discontinuous 2-order function and its partial y derivative by DSFPM

3 碰撞算例模擬

為了進一步驗證DSFPM的工程實用我們對碰撞算例進行模擬。由于DSFPM所選取參與計算的粒子數(shù)有限，且參與計算粒子的構(gòu)型不能退化為直線，因此對于大變形碰撞問題，DSFPM可能會產(chǎn)生精度不足或不穩(wěn)定現(xiàn)象。為了避免DSFPM處理大變形問題的缺陷，采用如圖7所示的算法流程圖。一方面，在變形初始階段仍使用DSFPM，若粒子應(yīng)變達到一定閾值時，DSFPM將自動轉(zhuǎn)化為DFPM，自適應(yīng)地實現(xiàn)大變形碰撞問題；另一方面，當(dāng)粒子的應(yīng)變未達到設(shè)定的閾值，而此時DSFPM中所選取的有限數(shù)目粒子的構(gòu)型不能滿足穩(wěn)定性，算法也轉(zhuǎn)化為DFPM算法。這節(jié)將對小變形和大變形問題分別進行模擬。

圖7 處理碰撞沖擊問題的算法框圖Fig.7 Algorithm diagram of the collision problem

3.1 離散公式

彈性體運動控制方程為：

(20)

式中：ρ、p、x和u分別為彈性體的密度、各向同性壓力、位移和速度。σ、g、c和ρ0分別為應(yīng)力、重力加速度、聲速和參考密度。利用SPH方法離散后，用α、β表示坐標(biāo)方向，控制方程(20)用指標(biāo)法可寫為：

(21)

σαβ=Sαβ-δαβp

(22)

(23)

(24)

(25)

3.2 小變形碰撞

這部分以鋁塊碰撞為例對小變形問題進行模擬，模型如圖8所示，兩鋁塊長20 mm,高10 mm,初始時左側(cè)鋁塊以20 m/s的初始速度向右運動，右側(cè)鋁塊靜止，M點為右側(cè)鋁塊碰撞界面處的中點。鋁參考密度為2 785 kg/m3，聲速取為5 328 m/s，楊氏模量為72 GPa，泊松比為0.3。粒子間距d=0.5 mm，光滑長度h=1.2d。

圖8 鋁塊碰撞模型Fig.8 The model for the aluminum block collision

分別將式(4)、(6)、(16)中的f(xi,yi)替換為vα(xi,yi)，實現(xiàn)FPM、DSPH和DSFPM等3種方法對于速度梯度的估計。為了保證粒子相互作用的對稱性，動量方程仍都采用式(21)進行計算。式(6)、(16)中只考慮了與粒子i同側(cè)的粒子信息，但是由于2個鋁塊之間存在能量傳遞，因此在DSPH和DSFPM計算中加入接觸力作為2個鋁塊相互作用的體現(xiàn)，當(dāng)粒子i支持域內(nèi)搜索到異側(cè)粒子s時，粒子s對粒子i施加一定大小的作用力，接觸力的具體表達式[20]為：

(26)

式中：ci為粒子i的聲速，mi和ms分別為粒子i和粒子s的質(zhì)量，W(ris)為核函數(shù)，ris=ri-rs。該接觸力直接加到動量方程中，且滿足mif(|ris|)=-msf(|rsi|)。

3.2.1速度分布

圖9 3種方法模擬得到的鋁塊速度隨時間的變化曲線Fig.9 Velocity-time cuvers of aluminum blocks by three methods

圖10 有限元模擬得到的鋁塊速度隨時間的變化曲線Fig.10 Velocity-time cuvers of aluminum blocks by finite element method

圖11給出了FPM、DSPH和DSFPM模擬得到的在鋁塊碰撞過程中某一時刻的速度分布圖。從圖11可以看出，該時刻兩側(cè)鋁塊正在進行速度交換，DSFPM和DSPH模擬得到的速度連續(xù)且分布均勻，但FPM模擬得到的速度在2個鋁塊界面附近出現(xiàn)很明顯的不連續(xù)，震蕩較明顯。這由于鋁塊速度相差較大，F(xiàn)PM的全支持域搜索會導(dǎo)致粒子速度梯度估計過大，造成FPM在界面附近速度梯度計算精度不足。

圖11 3種方法模擬得到的鋁塊碰撞過程中的速度分布Fig.11 Velocity distribution of the aluminum blocks by three methods

3.2.2應(yīng)力分布

2個鋁塊發(fā)生碰撞時，在鋁塊內(nèi)部會產(chǎn)生沖擊波，壓縮(拉伸)波在彈性體內(nèi)傳遞，圖12給出碰撞過程中某一時刻3種方法模擬的x方向應(yīng)力分布，該時刻2鋁塊應(yīng)力為負，受到擠壓。從圖12可以看出:(1)FPM模擬的界面附近出現(xiàn)很嚴重的正負應(yīng)力交錯現(xiàn)象，說明FPM在界面處的精度最低。(2)DSPH雖然大大改善了界面處應(yīng)力交錯現(xiàn)象，但其在模型四周邊界處出現(xiàn)了應(yīng)力不均勻的缺陷，體現(xiàn)出DSPH在邊界處核近似精度不足。(3)DSFPM模擬得到的應(yīng)力不論是在碰撞界面還是模型邊界都很均勻，說明了DSFPM處理不連續(xù)界面和邊界的優(yōu)勢。

圖12 3種方法模擬得到鋁塊碰撞過程中的x方向應(yīng)力分布Fig.12 x-direction total stress distribution of the aluminum blocks simulated by three methods

圖13 3種方法模擬得到點M處x方向應(yīng)力隨時間的變化曲線Fig.13 x-directional total stress changes at point M simulated by three methods

圖13給出3種方法模擬得到右側(cè)鋁塊點M的x方向應(yīng)力隨時間變化曲線。從圖13可以看到，該點應(yīng)力約在0.04 ms達到最低值，壓縮應(yīng)力達到最大，隨后，應(yīng)力迅速上升。最后DSFPM和DSPH模擬的應(yīng)力便保持穩(wěn)定，而FPM計算的應(yīng)力迅速上升后達到正值，出現(xiàn)應(yīng)力震蕩現(xiàn)象，和圖12的現(xiàn)象一致。

3.2.3計算效率

DSFPM計算耗時最短，為99.49 s；DSPH和FPM計算所用時間相當(dāng)，分別為171.57、173.89 s。這是由于在粒子法計算過程中，粒子鄰域搜索這一環(huán)節(jié)占計算總時間比重大，而根據(jù)框圖7，DSFPM是基于初始構(gòu)型對有限個粒子進行選取的，在計算過程中不需要進行鄰域搜索，這大大降低了計算時間。DSPH和FPM都需進行鄰域搜索，兩者僅在速度梯度計算方法上有所差異，因此計算時間相差很小。對該算例，DSFPM可以節(jié)省約40%的計算時間。

3.3 大變形碰撞

基于框圖7算法，對如下大變形問題進行模擬，如圖14所示，左側(cè)物塊以50 m/s撞擊右側(cè)板，板兩端固定，左側(cè)物塊長20 mm，高5 mm，密度為7 200 kg/m3，彈性模量E1=270 GPa，泊松比為0.3。右側(cè)板長2.5 mm，高20 mm，密度為2 785 kg/m3，彈性模量E2為72 GPa，泊松比為0.3。粒子間距d=0.2 mm，光滑長度h=1.2d。

圖14 物塊沖擊模型Fig.14 The model for block collosion

該算例中選取的應(yīng)變閾值為0.01，當(dāng)粒子應(yīng)變超過該值時，DSFPM算法將被轉(zhuǎn)化為DFPM算法。圖15給出了不同時刻框圖7算法模擬的速度分布。從圖15可以看出，右側(cè)板從中部開始變形，兩端固支附近彎矩最大，最后板兩端和中間附近出現(xiàn)裂口。該算法模擬的速度分布較均勻，體現(xiàn)了本文提出的算法在界面處計算精度的優(yōu)勢。

若其他條件不變，減小右側(cè)板彈性模量E2，便可以得到更大變形的結(jié)果。圖16為右側(cè)板彈性模量為72 MPa時不同時刻的碰撞速度分布圖。從圖16可以看出，當(dāng)變形增大時，框圖7算法仍然能夠模擬出連續(xù)均勻的速度。由此得知，DSFPM結(jié)合DFPM后，它不僅避免了DSFPM在處理大變形問題中參與粒子數(shù)較少的問題，也規(guī)避了DSFPM中有限粒子構(gòu)型可能帶來的穩(wěn)定性問題，是處理大變形問題的有效方法。

圖16 E2=72 MPa時，不同時刻框圖7算法模擬的速度云圖Fig.16 Velocity distribution of the blocks simulted by the algorithm in Fig.7 at different times when E2=72 MPa

4 結(jié) 論

針對FPM處理非連續(xù)物理場的缺陷，基于DSPH方法提出了DSFPM。通過對DSFPM進行近似精度分析發(fā)現(xiàn)DSFPM具有二階精度，可以準(zhǔn)確再生一次非連續(xù)函數(shù)。在鋁塊沖擊碰撞小變形問題中，通過對比分析DSFPM、DSPH和FPM等3種方法模擬的鋁塊速度、應(yīng)力以及計算時間，驗證了DSFPM在不連續(xù)界面處的計算優(yōu)勢。DSFPM大大改善了界面區(qū)域的計算精度，提高了計算效率，是FPM方法的有效改進。由于DSFPM算法本身選取粒子的特殊性，DSFPM在大變形問題中具有一定局限性，為了使改進算法仍然適用于大變形，我們提出了DSFPM自適應(yīng)轉(zhuǎn)化為DFPM的方法，并通過沖擊大變形算例驗證了該算法的可行性。