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

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

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

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

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

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

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

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

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

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

1 理論基礎(chǔ)

2.1SPH基礎(chǔ)[18]

連續(xù)體偏微分方程的求解可轉(zhuǎn)化為粒子的運(yùn)動(dòng)方程而求解, 而每個(gè)粒子的場(chǎng)變量通過積分插值獲得：

(1)

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

(2)

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

(3)

流體動(dòng)力學(xué)的基本控制方程是基于質(zhì)量守恒，動(dòng)量守恒，能量守恒三條基本的物理守恒定律來得到的。

1.2SPH-FEM算法流程

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

圖1 SPH-FEM耦合算法流程

1.3時(shí)間步

SPH采用跳蛙格式求解Navier-Stokes方程，F(xiàn)EM采用中心差分格式求解顯式動(dòng)力學(xué)方程。在Lagrangian框架下要求SPH與FEM的耦合面上積分同步，在同一時(shí)間點(diǎn)進(jìn)行數(shù)據(jù)傳輸[19]。FEM與SPH采用相同時(shí)間步。

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

(4)

(5)

(6)

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

2 沖擊試驗(yàn)

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

圖2 半球殼撞擊試驗(yàn)

試驗(yàn)組次半球殼直徑/m半球殼質(zhì)量/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

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

3 數(shù)值模型

3.1幾何和網(wǎng)格模型

采用LS-Dyna對(duì)半球殼撞擊試驗(yàn)進(jìn)行SPH-FEM耦合方法數(shù)值分析，如圖3所示。半球殼，四周及底面土壤采用有限元單元，內(nèi)部土壤采用SPH粒子。土壤數(shù)值模型尺寸按經(jīng)驗(yàn)確定：沖擊后模型外側(cè)最大應(yīng)變小于0.001%，屬于非常小應(yīng)變[21]。通過模型調(diào)試，土壤數(shù)值模型尺寸1.1 m×1.1 m×0.55 m。大變形部分采用SPH模擬，SPH尺寸0.7 m×0.7 m×0.35 m。有限元單元為實(shí)體單元，按六面體單元?jiǎng)澐志W(wǎng)格，劃分網(wǎng)格后半球殼共有8 908個(gè)單元，四周及底面土壤共有31 616個(gè)單元，內(nèi)部土壤共有171 500個(gè)SPH粒子。

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

3.2材料模型

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

(7)

式中：P為壓力；φ為內(nèi)摩擦角；J2為應(yīng)力張量的第二不變量；k(θ)為張量平面角函數(shù)；C為黏聚力；Ahyp為修正的Mohr_Cloulomb屈服面和標(biāo)準(zhǔn)的Mohr_Cloulomb屈服面相似程度。結(jié)合實(shí)驗(yàn)，土壤模型的主要材料參數(shù)見表2。

3.3模型邊界

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

表2 土壤材料參數(shù)

3.4數(shù)值計(jì)算結(jié)果與試驗(yàn)對(duì)比

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

圖4 試驗(yàn)1和SPH-FEM的加速度曲線對(duì)比

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

圖5 試驗(yàn)3和SPH-FEM的加速度曲線對(duì)比

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

表3 試驗(yàn)結(jié)果和SPH方法最大加速度對(duì)比(g=9.8 kg·m/s2)

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

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

圖6 25 ms時(shí)半球殼沖擊土壤的形狀

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

3.5沖擊能量分析

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

圖7 第一組SPH-FEM數(shù)值模型能量時(shí)程圖

3.6沖擊力曲線

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

圖8 8組數(shù)值模擬沖擊力時(shí)程曲線

4 結(jié) 論

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

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

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

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

[1] 龐勇，劉才山. 低速剛性體沖擊密實(shí)散體介質(zhì)的動(dòng)力學(xué)實(shí)驗(yàn)與理論研究[C]// 中國力學(xué)大會(huì)暨錢學(xué)森誕辰100周年紀(jì)念大會(huì). 哈爾濱市, 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] 黃雨，郝亮，謝攀，等. 土體流動(dòng)大變形的SPH數(shù)值模擬[J]. 巖土工程學(xué)報(bào), 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方法在剪切式碰撞能量吸收器中的應(yīng)用[J]. 振動(dòng)與沖擊, 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] 沈燕鳴，何琨，陳堅(jiān)強(qiáng)，等. SPH同一算法對(duì)自由流體沖擊彈性結(jié)構(gòu)體問題模擬[J]. 振動(dòng)與沖擊, 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 方法對(duì)不同藥型罩線性聚能射流形成及后效侵徹過程的模擬[J]. 振動(dòng)與沖擊, 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方法在巖土工程中的研究應(yīng)用進(jìn)展[J]. 巖土工程學(xué)報(bào), 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] 張志春，強(qiáng)洪夫，高巍然. 一種新型SPH-FEM耦合算法及其在沖擊動(dòng)力學(xué)問題中的應(yīng)用[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] 馬煒. 散體介質(zhì)沖擊荷載作用下力學(xué)行為理論分析與算法實(shí)現(xiàn)[D]. 北京:北京大學(xué), 2008.

[21] 汪中衛(wèi)，王海飆，戚科駿，等. 土體小應(yīng)變?cè)囼?yàn)研究綜述與評(píng)價(jià)[J]. 巖土力學(xué), 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

國家自然科學(xué)基金(50978064)

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

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

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

TD853

： A

10.13465/j.cnki.jvs.2017.17.029