羅 杰， 肖建春， 馬克儉， 毛家意

(1.貴州大學 空間結構研究中心，貴陽 550003； 2.貴陽市人民防空辦公室，貴陽 550003)

半球殼沖擊土壤的SPH-FEM耦合分析方法

羅 杰1， 肖建春1， 馬克儉1， 毛家意2

(1.貴州大學 空間結構研究中心，貴陽 550003； 2.貴陽市人民防空辦公室，貴陽 550003)

SPH方法作為一種具有無網格、拉格朗日性質的動力學求解算法，已經廣泛應用于沖擊動力學的研究。采用SPH-FEM耦合分析和8組試驗進行半球殼沖擊土壤的對比分析，比較加速度時程曲線和局部變形圖，驗證SPH-FEM方法在處理此類問題中的適用性。利用SPH-FEM方法得到的能量和沖擊力時程曲線，歸納出土壤的能量耗散機理。分析表明，數值模擬結果和試驗結果吻合較好，SPH-FEM方法能勝任模擬半球殼沖擊土壤的整個過程。土壤是一種很好的能量緩沖體。

SPH-FEM；沖擊；土壤；半球殼；加速度

土壤是由大量離散固體顆粒組成的集合體，是一種復雜的非平衡態(tài)的能量耗散體系。當受到沖擊時，部分顆粒會呈現流體的現象，并會發(fā)生強烈擠壓和摩擦[1]。散體的這些性質，很難找到一種具體模型來描述。目前散體的動力學行為還沒有成熟的理論，因此研究人員一般采用數值模擬和實驗來進行分析[2-3]。由于試驗復雜且影響因素過多，數值模擬方法是研究土壤受沖擊的重要手段之一。近年來，適合于大變形分析的離散介質力學方法和無網格方法等逐漸發(fā)展起來。其中離散介質力學方法可以真實地表達求解區(qū)域中的幾何狀態(tài)以及大量的不連續(xù)面，易于處理大變形、大位移和和動態(tài)問題；而無網格法則將整個求解區(qū)域離散為獨立的節(jié)點，對網格沒有依賴性，并基于大變形理論建立了不同無網格理論框架下可進行大變形問題分析的計算方法[4]。然而，它們的適用性亦存在一定的缺陷，原因在于，這些方法中非物理參數不易確定，難以直接描述土體的應力-應變關系，限制了其在本研究方向上的應用[5-6]。

光滑粒子流體動力學(SPH)是一種無空間網格的連續(xù)介質動力學計算方法。這種方法最初主要用于研究天體物理現象[7-9]，目前該方法的研究成果已涉及多個領域[10-15]。近年來，SPH方法開始涉及土體材料特性的研究，已被應用在山體滑坡，挖孔成樁等土體大變形問題中[16]。這些研究均證明SPH方法能夠精確描述土體在不同不變形階段的力學性質，并具有較高的計算精度和穩(wěn)定性[17]。與有限元或離散元等數值方法相比，SPH由于其拉格朗日和自適應特性，可以很好地處理大變形和后失穩(wěn)問題，避免了有限元計算中的網格過度扭曲可能導致的計算失敗，或離散單元法計算工作量大和參數標定困難等問題，因此是量化及解釋土體大變形和后失穩(wěn)機理的有力工具。本文采用SPH-FEM方法模擬土壤受沖擊的整個過程，通過與試驗數據對比，并對數值結果進行了詳細的驗證分析，證明該方法在處理土壤受沖擊問題中的適用性。該研究為以后分析散體介質受沖擊的動力學行為提供了一種穩(wěn)定的數值模擬方法，為以后的研究工作奠定了基礎。

1 理論基礎

2.1SPH基礎[18]

連續(xù)體偏微分方程的求解可轉化為粒子的運動方程而求解, 而每個粒子的場變量通過積分插值獲得：

(1)

在SPH方法中，整個系統(tǒng)是由具有獨立的質量，占有獨立空間的有限個粒子表示的。SPH計算時, 通過對支持域內的離散粒子求和來代替式(1)中的連續(xù)積分:

(2)

SPH方法通過使用光滑函數引進積分表示式。光滑函數不僅決定了函數近似式的形式、定義粒子支持區(qū)域的尺寸，而且還決定核近似的一致性和精度。本文采用的光滑核函數是B樣條曲線：

(3)

流體動力學的基本控制方程是基于質量守恒，動量守恒，能量守恒三條基本的物理守恒定律來得到的。

1.2SPH-FEM算法流程

在交界面上SPH粒子與FE單元耦合，通過罰函數將質點的力轉換到FE單元表面。當SPH粒子與FE單元接觸時，更新SPH粒子與FE單元節(jié)點的位移與速度等變量，產生接觸力。流程見圖1。

圖1 SPH-FEM耦合算法流程

1.3時間步

SPH采用跳蛙格式求解Navier-Stokes方程，FEM采用中心差分格式求解顯式動力學方程。在Lagrangian框架下要求SPH與FEM的耦合面上積分同步，在同一時間點進行數據傳輸[19]。FEM與SPH采用相同時間步。

ΔtSPH-FEM=min{ΔtSPH,ΔtFEM}

(4)

(5)

(6)

式中：h為SPH的光滑長度；L為最小單元尺寸；C為材料聲速；α為時間步長比例系數。

2 沖擊試驗

圖2為北京大學所做的撞擊試驗，一個半球殼以一定的初速度豎直撞擊地面。通過不同高度釋放具有不同直徑和質量的半球殼，利用安裝在半球殼中的加速度測量裝置測得了球殼落在土壤中的加速度曲線[20]。根據不同的半球殼直徑，質量以及撞擊速度，實驗共有8組，如表1所示。

圖2 半球殼撞擊試驗

試驗組次半球殼直徑/m半球殼質量/kg撞擊速度/(m·s-1)10.40812.0534.9720.40812.0543.1530.40812.0544.940.40824.531.9450.40824.539.4260.40824.545.3570.66243580.662440

在低速碰撞情況下，土壤顆粒很少跳起，但是球殼的沖擊在土壤中形成沙坑。沙坑的外圍出現空腔，同時球殼周圍一圈的土壤將翻起而略高于土壤表面，如圖1所示。

3 數值模型

3.1幾何和網格模型

采用LS-Dyna對半球殼撞擊試驗進行SPH-FEM耦合方法數值分析，如圖3所示。半球殼，四周及底面土壤采用有限元單元，內部土壤采用SPH粒子。土壤數值模型尺寸按經驗確定：沖擊后模型外側最大應變小于0.001%，屬于非常小應變[21]。通過模型調試，土壤數值模型尺寸1.1 m×1.1 m×0.55 m。大變形部分采用SPH模擬，SPH尺寸0.7 m×0.7 m×0.35 m。有限元單元為實體單元，按六面體單元劃分網格，劃分網格后半球殼共有8 908個單元，四周及底面土壤共有31 616個單元，內部土壤共有171 500個SPH粒子。

圖3 半球殼撞擊有限元模型

3.2材料模型

半球殼采用剛體模型，采用LS-Dyna提供的MAT147材料模型來作為土壤模型。此土壤材料采用修正的Mohr_Cloulomb 屈服準則，是一種適合于實體單元、考慮損傷的各向同性材料模型，其屈服面表達式[22]為

(7)

式中：P為壓力；φ為內摩擦角；J2為應力張量的第二不變量；k(θ)為張量平面角函數；C為黏聚力；Ahyp為修正的Mohr_Cloulomb屈服面和標準的Mohr_Cloulomb屈服面相似程度。結合實驗，土壤模型的主要材料參數見表2。

3.3模型邊界

半球殼與土壤SPH粒子定義為侵蝕接觸，土壤有限元單元與土壤SPH粒子定義為點面接觸，以保證不同算法間的協調一致性。在SPH與FEM耦合處理中，將SPH粒子定義為從節(jié)點，將與SPH粒子接觸界面上的有限元單元表面定義為主面。為了完全消除應力波反射作用，在模型四周及底面有限單元面施加無反射邊界，以描述半無限土體空間。

表2 土壤材料參數

3.4數值計算結果與試驗對比

選取試驗1和試驗3中加速度傳感器得到的半球殼質心的加速度時程曲線與數值模擬得到的半球殼質心的加速度時程曲線進行對比，如圖4，5所示。選取試驗與數值模擬中的最大加速度值進行對比，如表3所示。通過對比發(fā)現，數值模擬得到的最大加速度非常接近試驗得到的最大加速度。數值模擬結果和試驗結果吻合的比較好，SPH-FEM方法能較好的模擬半球殼沖擊土壤的整個過程。

圖4 試驗1和SPH-FEM的加速度曲線對比

Fig.4 Comparison of acceleration curves of first experiment and SPH-FEM

圖5 試驗3和SPH-FEM的加速度曲線對比

Fig.5 Comparison of acceleration curves of third experiment and SPH-FEM

表3 試驗結果和SPH方法最大加速度對比(g=9.8 kg·m/s2)

數值模擬與試驗有一定的誤差，導致該誤差的原因經分析主要在于以下幾點原因：① 加速度傳感器的響應頻率偏小，對高頻的測量不太準確(數值解在不同的時間步長下計算得到的加速度最大值相應的有所變化，時間步長越短得到的加速度值越大，數值解和試驗結果的對比需要在加速度傳感器的最高測量頻率相比較)[20]；② 土壤的緩沖性能與土壤的密實度有密切關系。在一般無網格方法中土壤密實度通過密度來控制[20]，SPH方法均分粒子的質量大小來模擬土壤的密實度。該數值模型中的土壤模型密實度較大，起到的緩沖的作用較小，從而使得加速度下降過程比試驗快(圖3，4)。③ 該SPH算法中粒子與粒子之間忽略摩擦滑移。

半球殼初始位置位于土壤上方，與土壤無接觸；開始瞬間與土壤接觸，土壤由于受到剪切和擠壓作用開始流動；隨著與土壤接觸面積的增大，土壤粒子流動范圍逐漸變大。直到25 ms后觀察土壤變化形狀，如圖6所示。土壤顆粒很少跳起，球殼的沖擊在土壤中形成沙坑。球殼周圍一圈的土壤將翻起而略高于土壤表面，這與試驗現象描述相符。

圖6 25 ms時半球殼沖擊土壤的形狀

Fig.6 Shape of the hemispherical shell impact the soil in 25 ms

3.5沖擊能量分析

圖7為第一組SPH-FEM數值模型的能量時程圖。在沖擊加載過程中，半球殼的初始動能部分轉化為土壤的內能，部分土壤翻起、飛濺，以及單元失效刪除耗散掉，還有一部分能量轉化為接觸能；半球殼與土壤接觸后，土壤內能迅速上升，之后由于土壤單元損傷超過閥值，逐漸出現單元刪除，土壤的內能逐漸下降，轉化為土壤的侵蝕能。半球殼動能大部分被土壤吸收，轉化為土壤內能，僅有很小一部分能量轉化為侵蝕能、接觸能和沙漏能。接觸能、沙漏能為峰值內能的10%以內，是被認為可以接受的范圍。另外SPH方法應用于動力學沖擊時在沖擊波波面應用質量守恒、動量守恒、能量守恒時要將動能轉化為熱能，提供沖擊波波面的耗散，這種能量的轉換往往需要采用人工黏度耗能的形式來表示，因此會產生一部分由人工黏度產生的能量。該耗量對于能量平衡是非常重要的，這使得該算法適合于模擬沖擊波問題。圖中人工黏度能量包含在砂土內能中。

圖7 第一組SPH-FEM數值模型能量時程圖

3.6沖擊力曲線

傳統(tǒng)的SPH方法存在粒子壓力場數值震蕩問題，引入人工黏度參數來考慮粒子壓力場影響。圖8給出了半球殼與砂土之間的沖擊力時程曲線?？梢?，在整個沖擊作用過程中，沖擊力時程曲線包括快速上升，快速下降，穩(wěn)定三個階段，即沖擊力在極短的時間內增加，然后迅速減小到0。土壤是一種很好的緩沖材料，半球殼與顆粒及顆粒間產生激烈的碰撞，使沖擊力得到衰減，從而達到了緩沖的效果。圖中可以得出初始沖擊能量越大，半球殼與土壤之間的沖擊力越大。8組數值模擬沖擊力變成0的時間均一樣，都為10 ms。土壤是一種能量快速耗散體系且沖擊力衰減時間與初始動能和半球殼直徑無關。

圖8 8組數值模擬沖擊力時程曲線

4 結 論

SPH-FEM耦合法能充分利用傳統(tǒng)有限元法的高計算效率和光滑質點流體動力學法處理土體大變形的優(yōu)勢。正是由于這些優(yōu)點，該方法已大量應用于巖土工程的研究中。本文通過采用SPH-FEM方法對土壤受沖擊的過程進行分析，并與試驗對比，得到以下結論：

(1) 利用SPH-FEM耦合算法得到的半球殼與土體碰撞時的局部變形行為以及加速度時程曲線與試驗較為吻合，證明了這種算法是能夠分析土壤受沖擊的整個過程。

(2) 半球殼動能大部分被土壤吸收，轉化為土壤內能，僅有很小一部分能量轉化為侵蝕能、接觸能和沙漏能。

(3) 通過數值模擬得到土壤是一種很好的緩沖材料；初始沖擊能量越大，半球殼與土壤之間的最大沖擊力越大。

[1] 龐勇，劉才山. 低速剛性體沖擊密實散體介質的動力學實驗與理論研究[C]// 中國力學大會暨錢學森誕辰100周年紀念大會. 哈爾濱市, 2011.

[2] GOLDMAN D I, UMBANHOWAR P. Scaling and dynamics of sphere and disk impact into granular media[J]. Physical Review E, 2008,77(1): 255-282.

[3] KATSURAGI H, DURIAN D J. Unified force law for granular impact cratering[J]. Nature Physics, 2007, 3(6): 420-423.

[4] 黃雨，郝亮，謝攀，等. 土體流動大變形的SPH數值模擬[J]. 巖土工程學報, 2009, 31(10): 1520-1524.

HUANG Yu, HAO Liang, XIE Pan, et al. Numerical simulation of large deformation of soil flow based on SPH method[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(10): 1520-1524.

[5] SASAKI T, OHNISHI Y, YOHINAKA R. Discontinues deformation analysis and its application to rock mechanics problems[J]. Journal of Geotechnical Engineering, JSCE, 1994,Ⅲ-27: 11-20.

[6] BAXTER G W, BEHRINGER R P. Cellular automata models of granular flow[J]. Physical Review A, 1990, 42(2): 1017-1020.

[7] LUCY L B. A numerical approach to testing of the fission hypothesis[J]. The Astronomical Journal, 1977, 82: 1013-1024.

[8] GINGOLD R A, MONAGHAN J J. Smoothed particle hydrodynamics theory and application to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society, 1977, 181: 375-389.

[9] MONAGHAN J J, LATTANZIO J C. A refined particle method for astrophysical problem[J]. Astronomy and Astrophysics, 1985, 149(1): 135-143.

[10] MICHEL Y, CHEVALIER J M, DURIN C, et al. Hypervelocity impacts on thin brittle targets: experimental data and SPH simulations[J]. International Journal of Impact Engineering, 2006, 33: 441-451.

[11] CLEARY P W, SINNOTT M, MORRISON R. Prediction of slurry transport in SAG mills using SPH fluid flow in a dynamic DEM based porous media[J]. Minerals Engineering, 2006, 19(15): 1517-1527.

[12] FANG J N, OWENS R G, TACHER L, et al. A numerical study of the SPH method for simulating transient viscoelastic free surface flows[J]. Journal of Non-Newtonian Fluid Mechanics, 2006, 139(1/2): 68-84.

[13] 許慶新，沈榮瀛，周海亭. SPH方法在剪切式碰撞能量吸收器中的應用[J]. 振動與沖擊, 2007, 26(9): 108-111.

XU Qingxin, SHEN Rongying, ZHOU Haiting. Applying SPH method to shearing energy absorbers[J]. Journal of Vibration and Shock, 2007, 26(9): 108-111.

[14] 沈燕鳴，何琨，陳堅強，等. SPH同一算法對自由流體沖擊彈性結構體問題模擬[J]. 振動與沖擊, 2015, 34(16): 60-65.

SHEN Yanming, HE Kun, CHEN Jianqiang, et al. Numerical simulation of free surface flow impacting elastic structure with SPH uniform method[J]. Journal of Vibration and Shock, 2015, 34(16): 60-65.

[15] 楊剛，傅奕軻，鄭建民，等. 基于SPH 方法對不同藥型罩線性聚能射流形成及后效侵徹過程的模擬[J]. 振動與沖擊, 2016, 35(4): 56-61.

YANG Gang, FU Yike, ZHENG Jianmin, et al. Simulation of formation and subsequent penetration process of linear shaped charge jets with different liners based on SPH method[J]. Journal of Vibration and Shock, 2016, 35(4): 56-61.

[16] BUI H H, FUKAGAWA R, SAKO K, et al. Lagrangian mesh-free particle method (SPH) for large deformation and post-failure of geomaterial using elastic-plastic soil constitutive model[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32: 1537-1570.

[17] 黃雨，郝亮，野々山澩人. SPH方法在巖土工程中的研究應用進展[J]. 巖土工程學報, 2008, 30(2): 256-262.

HUANG Yu, HAO Liang, NONOYAMA Hideto. The state of the art of SPH method applied in geotechnical engineering[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(2): 256-262.

[18] LIU G R, LIU M B. Smoothed particle hydrodynamics-a meshfree particle method[M]. New Jersey: World Scientific Publishing Company, 2003.

[19] 張志春，強洪夫，高巍然. 一種新型SPH-FEM耦合算法及其在沖擊動力學問題中的應用[J]. 爆炸與沖擊, 2011, 31(3): 243-249.

ZHANG Zhichun, QIANG Hongfu, GAO Weiran. A new coupled SPH-FEM algorithm and its application to impact dynamics[J]. Explosion and Shock Waves, 2011, 31(3): 243-249.

[20] 馬煒. 散體介質沖擊荷載作用下力學行為理論分析與算法實現[D]. 北京:北京大學, 2008.

[21] 汪中衛(wèi)，王海飆，戚科駿，等. 土體小應變試驗研究綜述與評價[J]. 巖土力學, 2007, 28(7): 1518-1523.

WANG Zhongwei, WANG Haibiao, QI Kejun, et al. Summary and evaluation of experimental investigation on small strain of soil[J]. Rock and Soil Mechanics, 2007, 28(7): 1518-1523.

[22] LEWIS B A. Manual for LS-DYNA soil material model 147[R]. USA: Federal Highway Administration Research and Development Turner Fairbank Highway Research Center, 2004.

SPH-FEMcoupledmethodforanalyzingahemisphericalshellimpactsoil

LUO Jie1, XIAO Jianchun1, MA Kejian1, MAO Jiayi2

(1. Space Structure Research Center, Guizhou University, Guiyang 550003, China;2. Guiyang City People’s Air Defense Office, Guiyang 550003, China)

SPH method as a meshless Lagrange dynamic solving method is widely applied in studying impact dynamics. Here, the SPH-FEM coupled method and 8 tests were adopted to analyze a hemispherical shell impact soil. The acceleration time history curves and local deformations obtained with SPH-FEM coupled method were compared with those gained in tests to validate the applicability of SPH-FEM method. The energy dissipation mechanism of soil was generalized with the energy and impact force time history curves obtained using SPH-FEM method. The results showed that the numerical simulation results agree well with the test ones; SPH-FEM coupled method can be used to correctly simulate the whole process of a hemispherical shell impact soil, soil is a good energy buffer.

SPH-FEM; impact; soil; hemispherical shell; acceleration

國家自然科學基金(50978064)

2016-08-12 修改稿收到日期：2016-12-11

羅杰 男，博士生，1989年8月

肖建春，男，博士，教授，1967年9月

TD853

： A

10.13465/j.cnki.jvs.2017.17.029