楊鴻發(fā) 許向寧 楊華陽 蹇代君

摘 要： 為了提高有限元法在解決滑坡較大變形時網格畸變的計算效率和精度，以九寨溝區(qū)域滑坡處地層土壤材料參數(shù)為基礎，采用SPH法建立邊坡模型，通過改變粒子密度和瑞利阻尼參數(shù)，研究了九寨溝區(qū)域滑坡的縮尺斜面模型滑坡面發(fā)生狀況，比較了不同參數(shù)下滑動面產生過程中的差異。結果表明，在應用SPH法分析地震邊坡變形的過程中，平滑距離的調整易使計算精度發(fā)生變化，從而產生一定程度的誤差，而平滑距離也會影響邊坡的最終變形形態(tài);在不同的SPH粒子密度和不同的瑞利阻尼參數(shù)下，滑動面的產生過程也呈現(xiàn)出顯著差異。SPH法有效改善了有限元處理大變形時網格失真等問題，研究結果可為地質災害易發(fā)區(qū)的滑坡災害防治提供一定參考。

關鍵詞： 巖土力學;九寨溝;滑坡;SPH;粒子密度;瑞利阻尼

中圖分類號： P315? ?文獻標識碼：? A

doi：? 10.7535/hbgykj.2020yx06007

Simulation study of Jiuzhaigou landslide based on SPH method

YANG Hongfa 1， XU Xiangning 1， YANG Huayang 1， JIAN Daijun 2

（1. 405 Geological Team， Sichuan Provincial Bureau of Geology and Mineral Exploration and Development， Chengdu， Sichuan? 611830， China; 2. Jiuzhaigou Scenic Spot Administration， Aba Tibetan and Qiang Autonomous Prefecture， Sichuan? 623400， China）

Abstract：

In order to improve the calculation efficiency and accuracy of the finite element method in solving mesh distortion of large landslide deformation， based on the soil materials parameters of the landslides in the Jiuzhaigou area， the SPH method was utilized to establish the slope model. By changing the particle density and Rayleigh damping parameters， the landslide surface occurrence of the scaled slope model of Jiuzhaigou lanolslides was studied and the differences in the sliding surface generation process under different parameters were analyzed. The results show that in the analysis of seismic slope deformation by SPH method， the calculation accuracy is easy to change when adjusting the smoothing distance， which results in a certain? degree? of error， and the final deformation of the slope can be affected by the smoothing distance. With the different SPH particle densities and Rayleigh damping parameters， the generation process of the sliding surface also shows significant differences. SPH method can improve the mesh distortion of large deformation by finite element method， and the results of this study can provide a reference for landslide diaster prevention and control in other geological disaster prone areas.

Keywords：

geotechnical mechanics; Jiuzhaigou; landslide; SPH; particle density; Rayleigh damping

山區(qū)地震發(fā)生時一般伴隨不同程度的崩塌、滑坡和泥石流等災害。2017年8月8日，四川省阿壩藏族羌族自治州九寨溝縣發(fā)生了7.0級地震，震中位于九寨溝縣漳扎鎮(zhèn)世界自然遺產九寨溝風景名勝區(qū)內，此次地震震源深度為20 km，最大烈度為9度，詳細的九寨溝地震烈度圖如圖1所示? [1] 。此次地震共造成四川、甘肅2省8個縣受災，25人死亡，525人受傷，6人失聯(lián)，176 492人受災， 73 671 間房屋不同程度受損。多個自然遺產點（如火花海、樹正群海、五花海、熊貓海）遭到不同程度的破壞，尤其是地震活動引起景區(qū)內大量山體出現(xiàn)不同程度的崩滑，地形，地貌，以及大面積森林植被、耕地和草地遭到嚴重損毀。10年內2次強烈地震破壞的疊加作用，致使九寨溝地質災害隱患更加突出，對九寨溝景區(qū)景觀的生態(tài)安全構成了巨大威脅? [2-3] 。

針對地震誘發(fā)的崩塌、滑坡，常用的研究方法有數(shù)值分析法、人工智能法、監(jiān)測法等? [4-5] 。通過數(shù)值分析法預測并分析地震引起的滑坡崩塌具有重大意義。有限元法是常用于滑坡解析的數(shù)值分析法? [6-8] ，雖然有限元法在預測邊坡穩(wěn)定性方面具有良好的效果，但采用有限元法解析九寨溝邊坡大變形時出現(xiàn)了網格失真、精度大幅降低等問題，對于邊坡的最終變形狀態(tài)和滑動破壞范圍的預測還不夠充分。

SPH（smoothed particle hydrodynmics）法將連續(xù)體看作多數(shù)粒子的集合，屬于粒子分析法，是一種適合處理大變形問題的解析方法，主要應用于流體分析領域。作為對壓縮性流體的分析方法，SPH法最初被用于宇宙物理學領域? [9] 。由于不需要網格并具有適合大變形解析等優(yōu)點，學者們開始將其運用于普通流體的解析，后來在一般彈塑性體的大變形解析中也加入了SPH法? [10-12] 。BUI等? [13] 將Drucker-Prager準則導入SPH法解析中，解析結果與砂土模型實驗結果高度相似。NAILI等? [14] 將液化后的地基視為黏性流體并使用SPH法解析，成功分析出流動地基對建筑物的地下構造物的荷載作用。

因此，在本研究中利用SPH法代替有限元法對地震時的邊坡狀態(tài)進行了分析驗證，并提出了瑞利阻尼的導入方法。瑞利阻尼是一種廣泛采用的正交阻尼模型，與頻率相關。在進行土層地震反應時有較好的解耦性能。

通過對九寨溝滑坡發(fā)生處的地層土壤材料建立邊坡模型并進行解析，確認了粒子密度（SPH法中特有的參數(shù)）以及瑞利阻尼在邊坡崩塌過程中的影響。

1 SPH法的理論基礎

1.1 基本理論

SPH法是以連續(xù)體作為研究對象，使用多數(shù)粒子集合的近似計算方法。一般進行核近似和粒子近似2種近似計算方式。核近似如式（1）所示。

SPH粒子和核函數(shù)的詳細關系可參考圖2，求解物理量 f時相當于求出離中央的黑色粒子一定距離內的粒子（圖中的深色粒子）所具有的值的加權平均值，在淺色粒子的位置 ，核函數(shù)的值為零則不含在計算中。

平滑距離 是常用參數(shù)，取值過小易導致SPH近似計算的粒子數(shù)減少，從而導致計算精度的降低。為了解決計算誤差問題，使用CHEN等? [15] 提出的CSPM（corrective smoothed particle method）法，將式（2）和式（3）用式（7）和式（8）代替：

1.3 彈塑性本構方程及應力增量

在 SPH法彈塑性分析中，可以使用有限元法中求應力增量的彈塑性矩陣法。形變速度張量和速度梯度由式（13）—式（15）求得：

1.4 瑞利阻尼的使用

在地震分析中經常會用到瑞利阻尼，有限元法中的瑞利阻尼矩陣[ C ]利用系數(shù)αR及βR表示 如下：

此時，式（17）中的剛性矩陣[ K ]在分析對象時相當于切線剛度矩陣。

2 模型和解析條件

2.1 模型建立

圖3為典型地震滑坡高分影像特征，九寨溝滑坡發(fā)生處地層土壤材料的物理性質見表1。圖4中的圓弧線為使用表1數(shù)據進行的經典斜面模型實驗所得到的滑動面，九寨溝景區(qū)熊貓海沿岸發(fā)生的滑坡按等比例縮小的模型實驗如圖5所示，實曲線即為模型實驗中產生的滑動面。筆者采用4種不同粒子排布密度的模型進行解析，如圖6所示。SPH粒子用規(guī)則格狀排布，粒子 間隔定義為dp，粒子總數(shù)定義為np。邊界條件為底面固定，側面只在水平向固定，粒子可以上下方向自由運動，平滑距離h為2.6dp。形變速度和加速度由上述各式分別計算得到。瑞利阻尼在一次模式和二次模式中使用相同的衰減常數(shù)。為了考 察不同衰減對分析結果產生的影響，使用瑞利阻尼分別為2%，5%，10%時進行對比分析。

2.2 邊坡自重作用分析

為重現(xiàn)邊坡的崩塌過程，對模型進行了自重作用分析。解析中因為彈塑性本構方程的使用，在震動時有可能發(fā)生粒子塑性化。圖7為重力加速度逐漸增加時的分析。為了減少自重作用分析所需的時間，采用BUI等? [16] 提出的衰減項代替運動方程：

式中 ξ為衰減強度參數(shù)，一般為0.002～0.005，本文采用ξ=0.005。

圖8表示了坡頂角垂直方向速度隨時間發(fā)生的變化。在粒子密度不同的模型中，衰減項均發(fā)揮了充分的效果，可以確認在分析結束時粒子不會產生震動。而圖9表示了自重作用分析結束時沿垂直方向的軸應力 σ? yy 的分布。所有模型的結果都很相似，對于初始應力狀態(tài)，粒子密度沒有太大的差異。

2.3 粒子密度對地震解析結果的影響

水平向輸入的地震加速度如圖10所示，地震作用下解析得到的邊坡變形情況如圖11所示，在震動開始后經過10 s，?dp =0.5 m時，最大剪切變形出現(xiàn)的地方和經典模型試驗結果非常相似。 dp =1.0 m和2.0 m時，最大剪切變形出現(xiàn)的地方與經典斜面模型試驗相比相對靠下。 dp =5.0 m時沒有出現(xiàn)明顯的滑動線。如圖11 b）所示，開始后經過30 s震動結束時， dp 為0.5，1.0，2.0 m時，震后得到的結果形狀幾乎相同。 dp 為10， 2.0 m時，剪切變形超過60%的區(qū)域基本一致。但在 dp= 0.5 m時，在較深的位置也出現(xiàn)了剪切形變超過60%的滑動線并延伸到了頂部。另一方面， dp =5.0 m時， 坡趾處出現(xiàn)剪切變形超過60%的區(qū)域但沒有延伸到頂部，與其他情況有明顯差別。

從解析結果可以看出，利用SPH法對邊坡崩塌的分析中，根據設定的不同粒子密度，其結果也會受到一定影響。

dp 為1.0 m和2.0 m時，得到了相似的結果。當粒子密度更大（ dp =0.5 m）時，出現(xiàn)了更多的滑動線，由于分析條件的差異，可以得出完全不同的崩塌狀態(tài)。經典斜面模型實驗結果是邊坡在圖4和圖5所示圓弧線處滑落，而解析結果更傾向于整個坡面崩塌。

2.4 瑞利阻尼的作用對地震結果的影響

圖12為瑞利阻尼分別為2%，5%，10%時的解析結果。震后10 s和震后30 s最大剪切變形出現(xiàn)的地方相對一致。瑞利阻尼為2%，5%，10%時，剪切形變超過60%的區(qū)域全部延伸到了頂部。隨著瑞利阻尼從10%下降到2%，剪切形變超過60%的部分在逐漸增多，可以清晰地看到至少產生2個滑動面。從圖12可以觀察到，施加的衰減越小，邊坡的剪切變形越大，并且區(qū)域分布也不相同。從解析結果看，滑動面的下部情況與經典斜面模型實驗相對吻合，滑動面上部最大剪切變形部分仍有待分析。綜上所述，由于瑞利阻尼可導致不同斜面的崩潰，需要進一步對使用的多個參數(shù)進行研究。

3 結 語

針對九寨溝區(qū)域分布的滑坡進行了分析，采用相對于傳統(tǒng)有限元法更為精確的SPH法建模并與室內斜面模型實驗結果進行了對比，研究了不同粒子密度以及瑞利阻尼作用對于邊坡滑動的影響。

1） 模型解析時，通過采用不同的SPH粒子密度和不同的瑞利阻尼，滑動面的產生過程隨之呈現(xiàn)出顯著區(qū)別。粒子密度為1.0 m和2.0 m時，計算結果與模型實驗較為相似。解析時施加的瑞利阻尼越小，邊坡的剪切變形越大。

2）采用SPH法對地震邊坡分析時，調整平滑距離易使計算精度發(fā)生變化，從而產生一定程度的誤差，而平滑距離會顯著影響邊坡的最終變形，因此計算時需要充分考慮使用合理的平滑距離。

3）利用SPH法對九寨溝發(fā)生的真實滑坡進行了縮小模型的再現(xiàn)分析，使用單純的粒子解析避免大變形失真問題，結果可為九寨溝等地質災害易發(fā)區(qū)域的滑坡災害防治提供重要參考。

雖然本研究對粒子密度和瑞利阻尼進行了一定的參數(shù)分析，但是SPH法模擬結果與斜面模型實驗產生的滑動面上部未能較好地吻合，原因可能在于各類參數(shù)的適用性范圍還存在一定的差別。在今后的研究中還需要對粒子密度、瑞利阻尼等各種影響參數(shù)的適用范圍作進一步分析。

參考文獻/References：

[1]? 徐錫偉，陳桂華，王啟欣，等.九寨溝地震發(fā)震斷層屬性及青藏高原東南緣現(xiàn)今應變狀態(tài)討論[J].地球物理學報，2017，60（10）：4018-4026.

XU Xiwei， CHEN Guihua， WANG Qixin，et al. Discussion on seismogenic structure of Jiuzhaigou earthquake and its implication for current strain state in the southeastern Qinghai-Tibet Plateau[J]. Chinese Journal of Geophysics， 2017，60（10）：4018-4026.

[2]? 房立華， 吳建平， 蘇金蓉， 等. 四川九寨溝Ms7.0地震主震及其余震序列精定位[J]. 科學通報，2018，63（7）：649-662.

FANG Lihua， WU Jianping， SU Jinrong，et al. Relocation of mainshock and aftershock sequence of the Ms7.0 Sichuan? Jiuzhaigou? earthquake[J]. Chinese Science Bulletin，2018，63（7）：649-662.

[3]? 李勇， 邵崇建， 李芃宇， 等. 九寨溝Ms 7.0級地震的左旋走滑作用與動力機制[J]. 成都理工大學學報（自然科學版）， 2017， 44（6）：641-648.

LI Yong，SHAO Chongjian，LI Pengyu，et al. Left-lateral strike-slip effect and dynamic mechanism of the Jiuzhaigou? Ms 7.0? earthquake in the eastern margin of Tibetan Plateau， China [J]. Journal of Chengdu University of Technology （Science and Technology Edition），2017， 44（6）：641-648.

[4]? 黎璽克.遺傳算法優(yōu)化BP神經網絡的巖質邊坡穩(wěn)定性預測[J].河北工業(yè)科技，2020，37（3）：164-169.

LI Xike. Prediction of rock slope stability based on BP neural network optimized by genetic algorithm[J]. Hebei Journal of Industrial Science and Technology， 2020，37（3）：164-169.

[5]? 于遠亮，靳靜.公路高邊坡穩(wěn)定性遠程監(jiān)控技術研究[J].河北工業(yè)科技，2014，31（3）：239-242.

YU Yuanliang， JIN Jing. Study on remote monitoring technology about highway slope stability[J]. Hebei Journal of Industrial Science and Technology， 2014，31（3）：239-242.

[6]? WATKINSON I M， HALL R. Impact of communal irrigation on the 2018 Palu earthquake-triggered landslides[J]. Nature Geoscience， 2019， 12（11）： 940-945.

[7]? FRANCI A， DE-POUPLANA I， CELIGUETA M A， et al. Application of the particle finite element method （PFEM） to the numerical simulation of landslides[C]// International Conference of Natural Hazards and Infrastructure （ICONHIC）. Athens： National Technical University of Athens， 2019：? 1-10.

[8]? WANG? Liang， ZHANG Xue， ZANIBONI F， et al. Mathematical optimization problems for particle finite element analysis applied to 2D landslide modeling[J]. Mathematical Geosciences， 2019.

doi： 10.1007/s11004-019-09837-1.

[9]? MONAGHAN J J. Simulating free surface flows with SPH[J]. Journal of Computational Physics， 1994， 110（2）： 399-406.

[10]? HU Wei， GUO Guannan， HU Xiaozhe， et al. A consistent spatially adaptive smoothed particle hydrodynamics method for fluid-structure interactions[J]. Computer Methods in? Applied? Mechanics and Engineering， 2019， 347： 402-424.

[11]? BAYAREH? M， NOURBAKHSH A， ROUZBAHANI F，? et al . Simulation of sand particles flow using a weakly compressible smoothed particle hydrodynamics method （WCSPH）[J]. Annales de Chimie： Science des Materiaux， 2019， 43（1）： 43-51.

[12]? ONO? Y， OKAMOTO R， KOHNO M， et al.3D simulation of the aratozawa landslides base on smoothed particle hydrodynamics method[J]. Journal of Japan Society of Civil Engineers， 2017， 73（4）： 346-356.

[13]? BUI? H H， FUKAGAWA R， SAKO K， et al. Lagrangian meshfree particles method （SPH） for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[J]. International Journal for Numerical and Analytical Methods in Geomechanics， 2008， 32（12）： 1537-1570.

[14]? NAILI M， MTSUSHIMA T， YAMADA Y.A 2D smoothed particle hydrodynamics method for liquefaction induced lateral spreading analysis [J]. Journal of Applied Mechanics， 2015， 8：591-599.

[15]? CHEN J K， BERAUN J E， JIH C J. Completeness of corrective smoothed particle method for linear elastodynamics[J]. Computational Mechanics， 1999， 24（4）： 273-285.

[16]? BUI H， FUKAGAWA R， SAKO K. A study of the matter of SPH application to saturated soil problems[C]// Proceedings of the 5th International SPHERIC Workshop.Manchester：? University? of Manchester， 2015： 354-361.