基于SPH-FEM耦合方法的泥石流沖擊柱形結構物動力響應分析

韓俊輝， 姚昌榮， 余勁松， 匡 睿， 強 斌， 李亞東

(西南交通大學土木工程學院，四川成都 610031)

[通信作者]姚昌榮(1974～)，男，博士，副教授,主要從事橋梁防災減災、結構損傷識別及健康監測、橋梁施工控制研究工作。

泥石流是世界范圍內常見的自然災害，多發于山地、丘陵和高原地區[1]。我國正在修建的川藏鐵路是繼青藏鐵路之后又一進藏鐵路干線，將西藏納入國家快速鐵路網中，可加速西藏經濟發展。川藏鐵路地質條件復雜，沿線泥石流災害頻發，而川藏鐵路橋梁占比高達13.2 %，橋梁面臨的泥石流危害也不容小覷。

常用的泥石流研究方法有野外實測調查法、數值模擬法和水槽試驗法等，其中數值模擬法以其較低的成本和研究上的便捷，逐漸成為研究泥石流廣泛采用的方法。傳統的泥石流數值模擬方法一般為基于網格的方法，但基于網格的方法在處理自由液面、大變形及多相流問題時有著難以捕捉自由液面、網格纏繞及網格需重新劃分等不可避免的缺陷，近些年來，無網格方法越來越多的被用來模擬泥石流[2][4]。SPH方法作為無網格方法中較早提出且目前較為成熟的一種方法，越來越受到研究泥石流學者們的重視。

為分析不同橫截面形式對泥石流沖擊柱體動力響應的影響，本文以川藏鐵路某泥石流溝的泥石流流變參數為例，采用非線性動力學數值模擬軟件LS-DYNA，基于SPH方法建立泥石流漿體，FEM方法建立柱體，研究泥石流發生概率P分別為10 %、5 %、2 %和1 %條件下各不同截面柱體受泥石流沖擊的動力響應。

1 SPH基礎理論

1.1 SPH方法簡介

SPH方法是一種基于粒子的純拉格朗日方法，主要用來模擬連續介質。該方法克服了基于網格方法的難以捕捉自由液面、難以處理自由邊界及物體發生大變形時網格纏繞等缺點，在模擬流體時具有基于網格方法不可比擬的優勢[3]。

SPH方法最早于1977年由Gingold、Monaghan[5]和Lucy[6]提出用來解決天體物理問題，1992年Monaghan[7]首先將SPH方法拓展到不可壓縮性自由表面流的模擬中。在計算流體領域中經過二十多年的發展，SPH方法在泥石流的模擬中已較為成熟。趙宏亮[8]采用Python語言基于SPH方法自編了模擬黏性泥石流運動過程的程序；Chen等[9]采用GIS建立的實際地形，基于SPH方法分析了在泥石流溝不同位置設置不同數量的攔擋壩對泥石流運動過程的影響；柳春等[10]采用SPH-FEM耦合的方法模擬了含大塊石泥石流沖擊攔擋壩的動力響應。

1.2 SPH方法的計算理論

SPH方法將計算域離散為攜帶物理屬性的粒子，通過支持域內粒子核函數近似的形式求得某一粒子的物理屬性[11](圖1)。

圖1 粒子近似示意

(1)

式中：i和j分別代表粒子i和j；x為待求物理屬性；N為粒子i支持域內的粒子總數；m為粒子的質量；ρ為粒子的密度；h為支持域半徑，又稱光滑半徑；W(xi-xj,h)為核函數，又稱光滑函數、權函數，形式如式(2)。

(2)

式中:d為所求問題的維數。

核函數形式有Bell-Shaped核函數、Gaussian核函數以及三次B-Spline核函數等，本文采用LS-DYNA中默認的三次B-Spline核函數，以θ(x)表示的形式為式(3)。

(3)

在離散的粒子框架下求解控制方程，流體控制方程Navier-Stoke方程的粒子近似形式為式(4)、式(5)[9]。

(4)

(5)

式中：t為時間；v為速度；α、β為坐標方向；p為壓力；F為外力，如重力、摩擦力等。

2 建立模型

2.1 流變參數

根據中科院山地所的勘察報告，川藏鐵路某泥石流溝泥石流為典型的暴雨型泥石流，其細顆粒黏粒含量較少，高頻率條件下泥石流以稀性泥石流為主，低頻率條件下則以過渡性泥石流或高容重低粘度的水石流為主。流變參數見表1。

表1 泥石流流變參數

2.2 幾何及網格模型

柱體高10m，采用實體單元建模，截面形狀按TB10002-2017《鐵路橋涵設計規范》[12]和JTGD60-2015《公路橋涵設計通用規范》[13]提供的流水壓力計算的5種橋墩截面，分別為圓形、正方形、菱形、矩形和圓端形。柱體迎流寬度均設為1.5m，如圖2所示。設置矩形和圓端形截面的目的是研究長寬比對沖擊動力響應的影響。

圖2 柱體截面形狀(單位：m)

在柱體周圍建立一個14.2m×12m×10m的空間作為泥石流的流動空間，左側和下側用殼單元模擬，在左側殼壁面右側0.2m處建立SPH初始位置點，以產生SPH粒子流，如圖3所示。

圖3 幾何模型(單位：m)

2.3 材料模型

因該泥石流溝細顆粒黏粒含量較少，高頻率條件下泥石流以稀性泥石流為主，低頻率條件下以過渡性泥石流或高容重低粘度的水石流為主，故不考慮泥石流粘度，泥石流材料模型采用*MAT_NULL模型，流變參數見表1。

泥石流的內部壓強根據狀態方程計算，采用Murnaghan狀態方程，計算式為式(6)。

(6)

式中：ρ0為粒子的初始密度；γ為系數, 對流體常取γ=7；B為系數，計算式為式(7)。

(7)

式中：c為人工聲速，c≥10vmax，vmax為流體運動過程中的最大速度。

柱體采用混凝土連續蓋帽模型[14](*MAT_CSCM_CONCRETE)，該模型是美國聯邦公路管理局為研究公路防撞護欄而開發，適用于低速撞擊等情形。材料密度取2 400kg/m3，其余參數采用默認值及參考文獻[14]，具體參數見表2。

表2 柱體材料參數

左側和下側殼單元邊界不考慮變形，采用剛體材料模型。

2.4 邊界條件

SPH粒子兩側采用LS-DYNA中SPH對稱平面邊界(*BOUNDARY_SPH_SYMMETRY_PLANE)，如圖4所示。該方法為SPH邊界處理方法中的虛粒子法，在對稱平面處生成鏡像粒子，以避免算法造成的SPH粒子向邊界聚集問題。此外，為節約計算時間，對SPH粒子設置激活區域，運動出此區域的粒子將失活，不再參與計算，如圖5所示。

圖4 SPH對稱邊界條件

圖5 SPH粒子激活區域

柱體邊界條件為約束底面節點全部自由度。殼單元在材料模型里約束全部自由度。

泥石流和柱體以及下側河道采用點面自動接觸的方式。泥石流和柱體的摩擦系數取0.12[10]。采用動力松弛方式初始化柱體模型，以消除重力導致的柱體初始振蕩；對柱體進行剛性類型沙漏控制，對泥石流進行粘性類型沙漏控制。

3 計算結果

3.1 沖擊過程

沖擊過程以發生概率P=10%的泥石流沖擊圓形柱體為例，如圖6所示。t=0.95s時泥石流前端到達柱體，開始對柱體產生沖擊，由于泥石流前端上層泥石流在重力作用下重力勢能轉化為動能，所以泥石流前端速度較大；t=2.60s時泥石流沖擊爬高，持續沖擊柱體，在柱體的背流側形成無流區；t=7.65s時，隨著泥石流沖擊的繼續進行，柱體背流區會逐漸淤積，沖擊已基本穩定。

圖6 沖擊過程

3.2 沖擊力

3.2.1 不同截面形狀對沖擊力的影響

泥石流發生概率P為10 %時，五種截面柱體所受來流方向沖擊力大小呈增大趨勢，如圖7所示。初始沖擊時，正方形柱體和矩形柱體所受沖擊力大小基本相同，圓形柱體、圓端形柱體和菱形柱體所受沖擊力基本相同；隨著沖擊的進行，矩形柱體所受沖擊力略大于正方形柱體，均大于其它三種柱體；圓形柱體、菱形柱體和圓端形柱體所受沖擊力有微小的差值，所受沖擊力從小到大順序沖擊過程中首先依次為圓形柱體、菱形柱體和圓端形柱體，三種截面柱體中，菱形柱體沖擊力最后穩定，最終穩定時菱形柱體沖擊力略大于圓端形柱體。為了消除隨機噪聲，對柱體沖擊力采用移動平均法處理[15]，時間窗口取1s。

圖7 泥石流發生概率P=10%條件下沖擊力

對移動平均處理后的沖擊力最大值進行分析，見表3。從圓端形柱體所受沖擊力大于圓形柱體，矩形柱體所受沖擊力大于正方形柱體可以得出：相同迎流面下，長寬比對柱體所受沖擊力有一定的影響，截面長寬比越大，柱體所受沖擊力越大。這是因為隨著長寬比的增大，柱體側面受到的泥石流摩擦力增大，并且與背流側無流區淤積等多種因素有關。從正方形柱體和矩形柱體所受沖擊力大于圓形柱體、圓端形柱體和菱形柱體可以得出：圓形迎流面和菱形迎流面有助于減小沖擊力。這與規范給出的圓形柱體所受泥石流沖擊力系數較小相一致。

以圓形柱體為基準，泥石流發生概率P=10%條件下，各截面柱體所受來流方向沖擊力的移動平均最大值與圓形柱體的對比見表3。由計算可知，矩形柱體所受沖擊力相較圓形柱體增加最大，為44.0 %，菱形柱體增加最小，為8.8 %。本文分析所用菱形柱體和圓端形柱體的泥石流沖擊力移動平均最大值始終相近，但圓端形柱體所受沖擊力受截面長寬比的影響，調整圓端形柱體橫截面的長寬比可使圓端形柱體所受沖擊力發生變化。

表3 沖擊力最大值(移動平均)對比

3.2.2 不同發生概率泥石流對沖擊力的影響

分析各截面柱體在不同發生概率泥石流條件下的沖擊力，并采用移動平均法計算得到其最大沖擊力，如圖8所示?？傮w來看，不同發生概率泥石流沖擊力，均是矩形、正方形截面柱最大，其余三種截面最大沖擊力接近。以泥石流發生概率P=10%為基準，各截面柱體所受來流方向沖擊力的移動平均最大值與P=10%條件下的對比見表4。由表4可知，P=5%條件下，各截面柱體所受沖擊力增加40 %～50 %左右；P=2 %條件下，各截面柱體所受沖擊力增加100 %～110 %左右；P=1%條件下，各截面柱體所受沖擊力增加190 %～210 %左右。

圖8 不同概率條件下最大沖擊力大小對比

表4 不同概率條件下最大沖擊力比較

3.3 柱頂水平位移

泥石流發生概率P=10%條件下菱形柱體水平位移最大，其次依次為正方形柱體、圓形柱體、圓端形柱體和矩形柱體，見圖9。柱頂水平位移與其截面的慣性矩及所受水平力的大小和分布有關。本文所選擇的5種截面形式，圓形、正方形、菱形、矩形、圓端形的截面慣性矩比值為1.0∶1.69∶0.42∶3.38∶2.70。菱形的截面慣性矩最小，且其所受水平力比圓形還大，因此其柱頂的水平位移最大。隨著泥石流發生概率的減小，由于泥石流的容重、速度均增大，因此各截面柱的沖擊力增大，導致柱頂水平位移增大，如表5所示。

圖9 P=10%條件下柱頂關鍵點位移時程

表5 柱頂最大水平位移mm

但是當P=2%、P=1%時，菱形截面柱體的位移急劇增大，且在新的位置來回振動，如圖10所示。主要是由于菱形截面柱的剛度最小，在較大泥石流的沖擊下，菱形截面柱體混凝土出現了大面積的損傷，柱體受力性能變差，使得柱體不能維持原來的平衡。

圖10 菱形柱體柱頂關鍵點水平位移時程

4 結論

SPH-FEM耦合的方法能夠很好地模擬泥石流沖擊柱體時迎流側沖擊爬高、背流側無流區域以及背流區泥石流淤積等沖擊過程。

(1)相同發生概率泥石流沖擊下，沖擊穩定后矩形柱體所受沖擊力最大；柱體橫截面長寬比對沖擊力有一定的影響，長寬比越大，所受沖擊力也越大。

(2)不同發生概率泥石流條件下，由于泥石流的容重、速度均增加，導致沖擊力均有不同程度的增加。

(3)同等條件下，菱形截面柱由于其剛度最小，導致其柱頂水平位移大于其它幾種截面柱，其受力性能最差，也最容易發生損毀，在結構設計中應避免選擇這種截面形式。

在本文的研究中仍存在一些局限：在采用SPH方法模擬泥石流漿體時，未考慮泥石流漿體的粘性；對泥石流建模時只考慮了泥石流漿體，未考慮泥石流中的粗顆粒和大塊石；對柱體底面的約束過于理想化等。這些需要繼續研究。