基于并行SPH 方法的地震滑坡對橋樁的沖擊作用

王占彬，張衛杰，張健，代登輝，高玉峰?

（1.巖土力學與堤壩工程教育部重點實驗室（河海大學），江蘇南京 210098；2.南京交通職業技術學院建筑工程學院，江蘇南京 211188）

伴隨著我國經濟的高速發展，公路和鐵路等基礎設施的建設和運營需要更安全的自然環境.據交通運輸部統計，截至2019 年，我國公路橋梁總里程達60 634.6 km，鐵路橋梁總里程數超過1 100 km.為振興和發展西部地區經濟，我國積極推進西部大開發和絲綢之路經濟帶戰略，兩大戰略分別涉及山地面積455 萬km2和429 萬km2（按各省份山地面積比例計算得出）.橋梁是溝通山區城市的最主要方式之一.橋梁建設和正常運營對邊坡穩定環境的要求不斷提高.

我國處于環太平洋地震帶和歐亞地震帶之間，地震頻發.截至2019年，我國發生有5人以上死亡或因災損失達到當地GDP 0.1%以上的地震災害近30年平均每年4 次，累計影響人次達0.75 億，直接經濟損失達1 060 億美元.地震導致了大量滑坡事件，橋樁結構因滑坡沖擊發生破壞的現象屢屢出現，如1976 年唐山7.8 級地震中橫跨薊運河的漢沽鐵路橋因河岸滑移嚴重導致橋墩傾斜；2008 年汶川地震滑坡導致10 座橋梁破壞嚴重或損毀，如都汶公路新房子大橋、一碗水中橋、順河大橋等受滑坡沖擊出現傾斜、剪斷、倒塌等不同程度損壞［1］，213 國道百花大橋軟弱場地土層向河心滑移導致橋樁基礎變位、傾斜［2］.橋樁結構對滑坡體沖擊荷載的承受力大小決定著它能否安全運行，因此研究地震荷載下邊坡滑坡沖擊橋樁結構過程具有重要意義.

目前，學者們對地震等誘發滑坡對橋樁等結構形成的沖擊主要以模型試驗、理論解析、數值分析為手段進行研究.在模型試驗方面，王友彪等［3］通過調整滑動體中黏土、砂土、礫土和水的質量分數研究土水混合體密度、滑動形態對橋墩承受沖擊荷載的影響.在理論解析方面，有學者假設沖擊荷載與作用在結構物上靜壓力成正比［4］，但因在理論上沖擊荷載是動力問題而非靜力問題，該假設未能得到后續較深入的研究.同時也有學者［5］認為該沖擊荷載在數值大小上與結構物承受的靜壓力和滑坡體滑動速度的平方成正比.在數值分析方面，基于網格離散的數值方法如有限單元法、有限差分法等在處理流滑大變形問題時將遭遇困難，故越來越多的學者開始采用無網格分析方法對土體的大變形動力學特性進行研究，其中一些研究（如Huang 等［6］、Dai 等［7］、Hu等［8］、王斌等［9］）利用具有無網格近似、追蹤粒子信息和能夠較好模擬不同材料相互作用性質的光滑粒子流體動力學（Smoothed Particle Hydrodynamics，簡寫為SPH）方法與流變模型結合，對地震誘發滑坡的大變形流滑過程進行了數值模擬.Bui等［10］在對巖土體材料滑坡大變形、Chen 等［11］在對巖土顆粒材料滑坡大變形的SPH 研究過程中，巖土體材料的應力應變關系采用非關聯流動法則Drucker-Prager 模型進行描述.Huang 等［12］在對滑坡中液化土體的大變形土體流滑數值模擬研究中采用土-水耦合SPH 方法，土和水分別在不同的計算層，土體采用彈性模型，水體被視為近似不可壓縮流體.唐宇峰等［13］采用SPH 方法對土體滑動大變形規律進行研究，土體應力應變關系采用Drucker-Prager 模型進行描述.Dai 等［14］引入流體-結構相互作用SPH 模型計算研究了文家溝和洪春溝滑坡對節制壩的沖擊影響，滑動土體采用Bingham 流體模型.綜合分析，學者們基于SPH 方法對流態性滑坡的大變形運動特性研究在一定程度上揭示了滑坡的大變形成災機制，但是這些研究缺少地震滑坡沖擊橋樁結構全過程的演化機理分析.

對此，本文基于SPH 基本原理提出用于模擬分析地震滑坡沖擊橋樁結構災變過程的并行SPH 計算方法，形成自編程序，采用Arias 烈度增加值方法截取地震波以減少計算空間的浪費，基于該方法研究不同地震加速度幅值、頻譜和持時條件下滑坡對橋樁結構的沖擊規律，并研究不同線程數對并行計算效率的影響.

1 并行SPH地震滑坡沖擊模型的建立

1.1 SPH滑坡沖擊的基本原理

SPH 方法是一種拉格朗日無網格近似方法，將計算區域劃分為一定數目的粒子，每個粒子都包含著自有的場變量，如速度、應力、壓力、質量、密度等.實體運動、變形等所有信息由這些粒子承載，同時追蹤每個粒子的運動，無網格及粒子間相互作用的特點使其更易處理大變形問題.SPH 無網格近似過程（如圖1 所示）包括核近似和粒子近似，假設場函數為f（x），則它的核近似和粒子近似表達式分別如式（1）（2）所示：

圖1 SPH無網格近似示意圖（目標粒子影響域為半徑為κh的圓形區域）Fig.1 SPH approximate schematic diagram（the influence domain of target particle is a circular area with radius κh）

式中：上標α表示坐標方向；特征寬度參數h是定義光滑核函數影響區域的光滑長度；W是光滑核函數；N是粒子i影響域內的粒子數量；mj、ρj分別是粒子j的質量、密度；Wij是粒子i在粒子j位置的值.

SPH 方法的基本控制方程包括連續方程、動量方程和能量方程［15］.本研究中將地震滑坡沖擊橋樁結構的過程假設為等溫過程，因此不考慮能量方程，得到無網格近似后的連續方程、動量方程表達式如式（3）（4）所示：

式中：表示α方向上的體力（如重力）、面力（如外荷載）在i粒子處引起的加速度；分別為速度、應力.

為了提高SPH 計算方法在邊界上的精度，將控制方程表達式用影響域內粒子核函數和的形式進行歸一化［16-18］，得到：

1.2 SPH滑坡沖擊的土-結構模型

研究中結構粒子的控制方程采用的應力形式為Piola-Kirchhoff 第一應力張量，質量守恒和動量守恒以完全Lagrange 方法表示，采用的歸一化表達式如式（7）（8）所示：

式中：含下標“0”的變量表示初始構形的狀態變量；為Piola-Kirchhoff 第一應力張量；Ji為變形梯度的行列式.

由滑動土體對結構體的沖擊動量方程得到：

通過影響域內為土體粒子的無網格粒子近似得到如下形式表達式：

式中：等號右側分子上第一項是土體壓力；第二項是土體粒子與結構粒子相對速度產生的壓力為結構粒子的加速度.

1.3 SPH滑坡沖擊的并行優化

Amdahl 定律給出了算法在運用并行計算之后，運行加速的估算［20］.假設算法按串行計算消耗的時間為t1，按并行計算消耗的時間為tp，tp分為并行計算部分和非并行計算部分，其中非并行計算部分的時間為ts，則可得到

式中：N是并行處理器的數量.得到加速比sp：

式中：假設并行處理器數N趨近無窮大時，加速比為sp=t1/ts.這反映了并行程序的加速比主要受限于非并行計算部分的時間ts，也顯示出了算法的運行中并行計算的部分占比越多，所消耗的時間越少.本研究采用張衛杰等［21］提出的并行計算方法.

1.4 SPH滑坡沖擊的地震加速度施加和驗證

本研究中土體采用的Drucker-Prager 模型是Bui等［10］改進的彈性-完全塑性土體本構模型.地震邊坡滑坡研究中土體阻尼系數和弱化指數借鑒Chen［22］對地震邊坡滑坡的研究.驗證模型中采用下負荷屈服面劍橋模型［23］，參數如表1 所示.計算模型兩側邊界設置為吸能邊界，以此防止地震波傳播到邊界時，邊界再反射回來對試驗模型或計算模型造成影響（如圖2 所示，左側為剛體邊界，右側為吸能邊界），以模擬實際場地的半無限區域特性，地震波由底部輸入，振動方向為水平.

表1 邊坡振動SPH模型參數Tab.1 Parameters of SPH model slope vibration

圖2 吸能邊界與剛體邊界模型示意圖Fig.2 Schematic diagram of energy absorption boundary and rigid body boundary

為驗證所提出方法在地震作用下土體大變形滑坡方面的應用，建立邊坡振動SPH 模型如圖3 所示.輸入的地震波采用修正Kobe 波，地震加速度為水平向.計算模型滑動土體部分由砂土粒子和軟黏土粒子間隔等間距分布、硬土層由砂土粒子和硬黏土粒子間隔等間距分布，軟黏土和砂土混合而成的邊坡土體厚度為20.1 cm，坡度為27.9°，硬黏土和砂土混合而成的邊坡土體左側厚4.8 cm、右側厚28 cm.模型箱寬1.6 m，高0.384 m.初始粒子間距0.8 cm，軟黏土粒子數為1 267，硬黏土粒子數為1 211，砂土粒子數為2 479，水體粒子數為4 957，地震邊界粒子數為897.左右兩側各設置3 層吸能邊界粒子，其數值為270.砂土密度為2 130 kg/m3，采用Drucker-Prager 模型，剪切波速79.5 m/s，泊松比0.3，內摩擦角30°，滲透系數1×10-4m/s，軟黏土和硬黏土采用能夠描述軟化和硬化行為的下負荷面劍橋模型，材料參數如表1所示，混合土體密度均為1 760 kg/m3，混合土體為飽和狀態，計算結果與振動臺試驗數據（Wartman等［24］）進行對比.

圖3 邊坡振動SPH模型Fig.3 SPH model for slope vibration

圖4 和圖5 分別給出了振動后土體粒子分布圖和振動后土體粒子最大剪應變云圖.從圖4 中可以看出，計算結果振后軟黏土粒子外輪廓線與物理模型試驗振后土體輪廓線較接近，表明計算模型的振后土體滑動位移與物理模型的土體滑動位移大致相近.試驗土體變形較計算值大，原因是試驗中軟化土層的下部和硬化土層的上部分別被嵌附一層0.5 mm厚光滑的高密度聚乙烯土工膜作為預設滑動面，從云圖中可以看出，土體最大剪應變位置沿軟化土層和硬化土層接觸面出現，由此可得出結論：計算值和試驗值的土坡滑動規律大致相同.

圖4 振動后土體粒子分布圖Fig.4 Soil particles distribution after vibration

圖5 振動后土體粒子最大剪應變云圖Fig.5 Maximum shear strain cloud figure of soil particles after vibration

1.5 SPH滑坡沖擊橋樁結構的數值模擬工況

Arias 烈度（Ia）［25］是反映觀測點地震加速度產生能量大小的重要指標，Ia綜合考慮了地震動幅值、地震動頻率和地震動持續時間的影響，已被很多學者［26-28］用于地震邊坡滑坡潛在危險分析.0.8g峰值加速度的Kobe 地震波加速度持續時間為42 s，時程曲線如圖6 所示，其Arias 烈度達0.583 m/s，在5.29～20.29 s 時段內Arias 烈度增長幅度為0.575 m/s，占總增長的99.02%，說明該地震在此測點15 s 時間段內能夠釋放大部分的能量.選取5.98～20.98 s時段Kobe地震波和全時程地震波作為計算模型的地震動輸入.后文將42 s Kobe 地 震波和5.29～20.29 s 時段Kobe地震波分別簡稱為Kobe波和Kobe截取波.

圖6 0.8g峰值Kobe地震波加速度時程及Arias烈度曲線Fig.6 Acceleration time history and Arias intensity curves of 0.8g Kobe seismic wave

本研究計算采用的地震峰值加速度為0.2g、0.4g、0.6g、0.8g；地震反應譜以時間變化系數的形式進行改變：0.5、1.0、2.0；地震持續時間選取1 次和疊加2次相同的典型時段截取波.

2 地震滑坡對橋樁結構的沖擊規律

2.1 原波與截取波荷載下滑坡結果對比

為比較Kobe 截取波和Kobe 波在本研究算例的差別，選取5.98～20.98 s時段Kobe地震波和全時程地震波作為計算模型的地震動輸入，建立模型如圖7所示.布置4 個邊坡測點和1 個邊界測點，編號分別為S1、S2、S3、S4、B1，其中S1、S2 分別是42 s Kobe 地震波和15 s Kobe截取地震波計算模型總位移最大值點，S3、S4 分別是邊坡坡腳和坡頂的測點，B1 用于測量邊界在施加42 s Kobe 地震波和15 s Kobe 截取地震波加速度之后得到該點的位移，驗證地震加速度施加情況.

圖7 Kobe截取波SPH地震計算模型Fig.7 SPH seismic calculation model with Kobe intercepted wave

如圖7 所示，所建模型采用的粒子間距為0.25 m，粒子總數為13 623，其中滑動土體粒子數為8 592，橋樁結構粒子數為1 350，邊界粒子數為3 681.土體采用DP 模型，參數見表2.邊界測點B1 施加兩種地震波后的計算結果如圖8、圖9 所示，從圖中可看出SPH 測得位移時程與實際加速度積分后得到的理論值基本一致，隨著時間的推移，理論位移偏移終值分別為0.05 mm、-20.56 mm，SPH 計算值分別為1.23 mm、-41.15 mm，位移累積偏移差值分別為1.18 mm、-20.59 mm.

圖8 42 s Kobe波邊界測點B1輸出位移時程Fig.8 Displacement time history of boundary point B1 under 42 s Kobe wave

表2 邊坡振動SPH模型參數Tab.2 Parameters of SPH model slope vibration

圖9 15 s Kobe截取波邊界測點B1輸出位移時程Fig.9 Displacement time history of boundary point B1 under 15 s Kobe intercepted wave

圖10～圖13 分別給出了4 個測點S1、S2、S3、S4在42 s Kobe 波和15 s Kobe 截取波作用下的水平向和豎向位移對比.從圖中可以看出S1 測點15 s 地震波輸入工況相較于42 s 地震波輸入工況在水平向和豎向存在最大3.4 s 的延遲位移，42 s 地震波輸入工況最大水平位移為40.60 m、豎向位移為-15.29 m，15 s 地震波輸入工況最大水平位移為40.42 m、豎向位移為-14.58 m，水平向和豎向最大位移差分別為0.18 m、0.71 m，差值占比分別為0.44%、4.65%；S2 測點15 s 地震波輸入工況相較于42 s 地震波輸入工況在水平向和豎向存在最大3.4 s 的延遲位移，42 s 地震波輸入工況最大水平位移為40.54 m、豎向位移為-15.36 m，15s 地震波輸入工況最大水平位移為40.53 m、豎向位移為-14.72 m，水平向和豎向最大位移差分別為0.01 m、0.64 m，差值占比分別為0.02%、4.14%；S3 測點15 s 地震波輸入工況相較于42 s 地震波輸入工況在水平向和豎向存在最大3.3 s的延遲位移，42 s地震波輸入工況最大水平位移為34.30 m、豎向位移為-13.76 m，15 s地震波輸入工況最大水平位移為34.10 m、豎向位移為-13.97 m，水平向和豎向最大位移差分別為0.20 m、0.21 m，差值占比分別為0.58%、1.53%；S4 測點15 s 地震波輸入工況相較于42 s 地震波輸入工況在水平向和豎向存在最大3.5 s的延遲位移，42 s 地震波輸入工況最大水平位移為24.84 m、豎向位移為-19.81 m，15 s 地震波輸入工況最大水平位移為23.57 m、豎向位移為-19.73 m，水平向和豎向最大位移差分別為1.27 m、0.08 m，差值占比分別為5.12%、0.40%.

圖10 42 s Kobe波和15 s Kobe截取波邊坡測點S1位移對比Fig.10 Displacement comparison at S1 in slope under 42 s Kobe wave and 15 s Kobe intercepted wave

圖11 42 s Kobe波和15 s Kobe截取波邊坡測點S2位移對比Fig.11 Displacement comparison at S2 in slope under 42 s Kobe wave and 15 s Kobe intercepted wave

圖12 42 s Kobe波和15 s Kobe截取波邊坡測點S3位移對比Fig.12 Displacement comparison at S3 in slope under 42 s Kobe wave and 15 s Kobe intercepted wave

圖13 42 s Kobe波和15 s Kobe截取波邊坡測點S4位移對比Fig.13 Displacement comparison at S4 in slope under 42 s Kobe wave and 15 s Kobe intercepted wave

圖14、圖15 分別給出了兩種地震波作用下的總位移云圖.從圖中可以看出滑動土體粒子的位移規律是相似的，即滑動機理沒有因為截取地震波輸入而出現與原波差距較大的變化.由兩者的計算結果比較可知橋樁結構所受荷載、滑動土體位移在兩種地震動輸入情況下相差很小，為提高研究效率，選取此15 s 時段地震加速度作為本小節數值模型計算研究的地震輸入.

圖14 42 s Kobe波地震荷載下計算模型總位移云圖（單位：m）Fig.14 Total displacement cloud figure of calculation model under 42 s Kobe wave（unit：m）

圖15 15 s Kobe截取波地震荷載下計算模型總位移云圖（單位：m）Fig.15 Total displacement cloud figure of calculation model under 15 s Kobe intercepted wave（unit：m）

2.2 地震加速度幅值的影響

建立計算模型如圖16 所示.布置測點S1，坐標（52.75，16.25），用于監測位移變化；布置A-A斷面，橫坐標為54.25，用于監測滑坡土體通過此斷面的體積量.所研究地震加速度幅值分別為0.2g、0.4g、0.6g、0.8g.土體參數如表2所示.

圖16 地震加速度幅值影響計算模型及測點布置Fig.16 Model of seismic acceleration amplitude influence and layout of measuring point

圖17～圖20 分別給出了0.2g、0.4g、0.6g、0.8g峰值加速度作用下的位移云圖.從圖17、圖18 中可以看出，0.2g峰值加速度和0.4g峰值加速度作用下土體最終位移和滑動的土體體積量很接近（0.2g工況S1 測點最大水平向位移為34.71 m，最大豎向位移為-13.89 m；0.4g工況S1 測點最大水平向位移為34.76 m，最大豎向位移為-13.83 m），4種峰值加速度作用下的S1 測點最終位移均很接近的原因是該測點與橋樁結構的水平距離為35.28 m，0.6g工況S1測點最大水平向位移為34.71 m，0.8g工況S1 測點最大水平向位移為35.46 m，與該測點和橋樁結構水平距離相差很小，該測點經過地震荷載作用已滑動到接近橋樁結構的位置.

圖17 0.2g峰值加速度下計算模型總位移云圖（單位：m）Fig.17 Total displacement cloud figure（PGA=0.2g）（unit：m）

圖18 0.4g峰值加速度下計算模型總位移云圖（單位：m）Fig.18 Total displacement cloud figure（PGA=0.4g）（unit：m）

圖19 0.6g峰值加速度下計算模型總位移云圖（單位：m）Fig.19 Total displacement cloud figure（PGA=0.6g）（unit：m）

圖20 0.8g峰值加速度下計算模型總位移云圖（單位：m）Fig.20 Total displacement cloud figure（PGA=0.8g）（unit：m）

圖21 給出了4 種峰值加速度作用下S1 測點水平向和豎向位移-時間關系對比曲線，從圖中可以看出，0.4g峰值加速度作用下土體位移快速增大的時間要早于0.2g峰值加速度作用下的土體位移快速增大的時間4.2 s.0.6g峰值加速度作用下土體位移快速增大的時間要早于0.2g峰值加速度作用下的土體位移快速增大的時間7.6 s.0.8g峰值加速度作用下土體位移快速增大的時間要早于0.2g峰值加速度作用下的土體位移快速增大的時間9.5 s.這說明峰值加速度的提高加速了邊坡的破壞和加大了邊坡的滑坡體積量.監測斷面A-A，0.2g通過此斷面體積量為51.81 m3，占比9.65%；0.4g為54.94 m3，占比10.23%；0.6g為90.19 m3，占比16.79%；0.8g為266.125 m3，占比49.56%.監測斷面數據也說明0.2g峰值加速度工況和0.4g峰值加速度工況滑坡體積量很相近，0.6g峰值加速度工況滑坡體積量稍大于0.2g和0.4g工況滑坡體積量，0.8g峰值加速度工況滑坡體積量最大，且超過0.2g工況滑坡5 倍之多.地震加速度幅值的提高在增加土體邊坡滑坡體積量的同時，也加快了邊坡土體的滑動破壞.

圖21 不同峰值加速度下S1測點位移-時間對比曲線Fig.21 Displacement-time comparison curves at S1 under different PGAs

2.3 地震加速度頻譜的影響

為研究不同地震加速度頻譜對邊坡滑坡沖擊橋樁荷載的影響，將15 s Kobe 截取地震波壓縮至7.5 s和延長至30 s，定義時間變化系數為現研究地震波時長與原時長的比值.本節所選取的時間變化系數tvc為0.5、1.0、2.0.3種時間變化系數的地震加速度時程如圖22 所示.它們經過傅里葉變換以頻域的形式展示如圖23 所示，它們的波峰頻率分別為0.27、0.54、1.08，隨著時間變化系數的提高，波峰頻率增大.計算模型如圖16 所示，布置測點S1，坐標（52.75，16.25），用于監測位移變化；布置A-A斷面，橫坐標為54.25，用于監測滑坡土體通過此斷面的體積量.土體參數如表2所示.

圖22 不同時間變化系數下地震加速度時程Fig.22 Seismic acceleration time history under different tvc

圖23 不同時間變化系數下地震反應譜頻域分布Fig.23 Frequency domain distribution of seismic response spectrum under different tvc

圖24、圖25 分別給出了時間變化系數為0.5 和2.0 的邊坡地震位移云圖，圖26 給出了3 種時間變化系數下地震邊坡模型S1 測點水平向位移、豎向位移與時間的對比關系曲線.從位移云圖中可以看出，時間變化系數為0.5 時最大位移和最小位移分別為37.4 6m、2.05 m，通過斷面A-A的體積量為59.31 m3，占滑坡體總體積量的11.05%；時間變化系數為1.0時最大位移和最小位移分別為43.13 m、20.53 m，通過斷面A-A的體積量為266.13 m3，占滑坡體總體積量的49.56%；時間變化系數為2.0時最大位移和最小位移分別為43.65 m、17.71 m，通過斷面A-A的體積量為273.44 m3，占滑坡體總體積量的46.08%.時間變化系數由0.5 增大至1.0，土體滑坡位移和體積量會明顯增多，時間變化系數由1.0 增大到2.0，土體滑坡最終位移和滑動體積量無明顯變化.圖26 中時間變化系數2.0 相較于1.0 時在15～30 s 出現的水平位移波動由邊界位移波動引起，邊界位移波動如圖9 所示.從中也可以看出，時間變化系數由1.0 減小至0.5，延緩了滑坡體的位移變化（延緩約2.6 s），也減小了滑坡體的滑動速度，時間變化系數由1.0增大至2.0，滑坡位移未出現較大變化，提前約0.9 s，較小程度地提高滑坡體的滑動速度，原因可能是邊坡土體在時間變化系數為1.0 的情況下已經達到足夠大的滑坡位移.

圖24 時間變化系數為0.5計算模型總位移云圖（單位：m）Fig.24 Total displacement cloud figure（tvc=0.5）（unit：m）

圖25 時間變化系數為2.0計算模型總位移云圖（單位：m）Fig.25 Total displacement cloud figure（tvc=2.0）（unit：m）

圖26 不同時間變化系數下S1測點位移-時間對比曲線Fig.26 Displacement-time comparison curves of measuring point S1 under different tvc

2.4 地震加速度持時的影響

建立計算模型如圖27 所示，布置測點S1、S2，坐標分別為（52.75，16.25）和（38.75，28），用于監測位移變化；布置A-A斷面，橫坐標為54.25，用于監測滑坡土體通過此斷面的體積量.所考慮的地震加速度持時影響是選取1 次的15 s Kobe 截取波和2 次截取波時間上的疊加，持續時間分別為15 s、30 s.土體參數如表2所示.

圖27 地震加速度持時影響計算模型及測點布置Fig.27 Calculation model of seismic acceleration duration influence and layout of measuring points

圖20、圖28 分別給出了1 次15 s Kobe 截取地震波和2次疊加截取地震波的位移云圖，圖29、圖30給出了兩種持續時間下地震邊坡模型S1、S2 測點水平向位移、豎向位移與時間的關系對比曲線.從云圖中可以看出，2次疊加截取地震波最終最大和最小位移分別為45.59 m、21.26 m，高于1次截取地震波邊坡最大位移和最小位移2.46 m、0.73 m，1 次截取地震波邊坡最大和最小位移分別為43.13 m、20.53 m.1次截取地震波和2 次疊加截取地震波滑坡土體通過斷面A-A的體積量相同，均為266.125 m3，占滑坡體總體積量的49.56%，滑坡體積量均很大.從位移-時間關系曲線中可以看出，2 次疊加截取地震波與1 次截取地震波作用的邊坡，對S1 測點位移的影響差別很小，但對S2 測點在15 s（即第1 次截取地震波作用結束后）出現增大，S2 測點最終水平位移、豎向位移在1次和2次疊加截取地震波作用下分別為40.53 m、14.72 m 和43.03 m、15.06 m，2 次多于1 次2.5 m、0.34 m.這說明2 次疊加截取地震波對滑坡體的中部和后緣部分位移影響較大，對前緣部分無較大影響.

圖28 2次截取地震波持時疊加計算模型總位移云圖（單位：m）Fig.28 Total displacement cloud figure under wave superposition（unit：m）

圖29 1次和2次截取地震波持時疊加測點S1位移-時間對比曲線Fig.29 Displacement-time comparison curve of S1 under single and superpositioned wave

圖30 1次和2次截取地震波持時疊加測點S2位移-時間對比曲線Fig.30 Displacement-time comparison curve of S2 under single and superpositioned wave

2.5 SPH滑坡沖擊模型的并行優化效率

為驗證并行優化的效率，計算研究了線程數分別為1、2、4、6、8、10、12、16、20、24、28，峰值加速度分別為0.2g、0.4g、0.6g、0.8g，時間變化系數Rt分別為0.5、1.0、2.0，施加1 次截取波、循環施加2 次截取波地震作用下邊坡土體滑動沖擊橋樁結構工況（其中峰值加速度為0.8g，時間變化系數Rt為1.0，1 次15 s截取Kobe 波為同一工況，在圖中簡寫為0.8g），按式（13）得到各線程數下的加速比sp（多線程和單線程計算消耗時間的比值），如圖31 所示.圖中虛線為加速比隨線程數增加的均值趨勢線，從圖31 中可以看出，在該計算平臺下，線程數從1 增加到8，加速比呈快速增長趨勢；線程數從8 增加到20，加速比增長速度稍緩；線程數從20 增加到28，加速度增長很緩慢，表明計算中線程數為20～28 時，計算機處理效能的利用已接近或達到最大狀態.線程數為28 時最大達6.1，平均5.7，其計算時間近似是單線程計算時間的1/6.該結果說明本研究土-結構耦合SPH 并行方法極大地提高了地震作用下邊坡土體滑移計算效率，減少了計算時間.

圖31 地震滑坡沖擊橋樁結構分析的加速比隨線程數的變化曲線Fig.31 Acceleration ratio curve with the number of threads in landslide impact analysis of bridge pile under earthquake

3 結論

1）基于SPH 方法的基本原理，建立了考慮靜壓力和相對速度產生的壓力的土-結構粒子相互作用模型，通過地震滑坡模型驗證了SPH 方法模擬振動臺試驗的準確性.

2）根據Arias 烈度增長值確定15 s 截取Kobe 地震波代替42 s 原Kobe 地震波，減小了對計算空間的浪費.

3）高地震加速度幅值會大幅增加可滑動土體的體積量，0.4g相對于0.2g滑坡體積量增加不明顯，0.8g滑坡體積量超過0.2g滑坡5倍之多，低地震加速度幅值會大幅加快滑動速度，0.4g下土體位移快速增大的時間要早于0.2g下的時間4.2 s，0.8g下土體位移快速增大的時間要早于0.6g下的時間1.9 s.

4）地震加速度頻率由0.27 增大至0.54，土體滑坡體積量增大近4 倍，滑動位移加快2.6 s，由0.54 增大到1.08，土體滑坡體積量無明顯變化，滑動位移加快0.9 s.

5）地震加速度持時對滑坡體后緣部分影響較大，對后緣測點的位移增加明顯高于前緣測點，滑動到橋墩結構的累積土體體積量隨持時的增加而增加.

6）對多種工況的SPH 程序進行不同線程數的并行優化，總結得出線程數從1 增加到8，加速比呈快速增長趨勢；線程數從8 增加到20，加速比增長速度稍緩；線程數從20 增加到28，加速度增長很緩慢，表明計算中線程數為20～28 時，計算機處理效能的利用已接近或達到最大狀態.