王占彬，張衛(wèi)杰，張健，代登輝，高玉峰?

（1.巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室（河海大學(xué)），江蘇南京 210098；2.南京交通職業(yè)技術(shù)學(xué)院建筑工程學(xué)院，江蘇南京 211188）

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

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

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

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

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

1.1 SPH滑坡沖擊的基本原理

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

圖1 SPH無網(wǎng)格近似示意圖（目標(biāo)粒子影響域?yàn)榘霃綖棣蔴的圓形區(qū)域）Fig.1 SPH approximate schematic diagram（the influence domain of target particle is a circular area with radius κh）

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

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

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

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

1.2 SPH滑坡沖擊的土-結(jié)構(gòu)模型

研究中結(jié)構(gòu)粒子的控制方程采用的應(yīng)力形式為Piola-Kirchhoff 第一應(yīng)力張量，質(zhì)量守恒和動(dòng)量守恒以完全Lagrange 方法表示，采用的歸一化表達(dá)式如式（7）（8）所示：

式中：含下標(biāo)“0”的變量表示初始構(gòu)形的狀態(tài)變量；為Piola-Kirchhoff 第一應(yīng)力張量；Ji為變形梯度的行列式.

由滑動(dòng)土體對(duì)結(jié)構(gòu)體的沖擊動(dòng)量方程得到：

通過影響域內(nèi)為土體粒子的無網(wǎng)格粒子近似得到如下形式表達(dá)式：

式中：等號(hào)右側(cè)分子上第一項(xiàng)是土體壓力；第二項(xiàng)是土體粒子與結(jié)構(gòu)粒子相對(duì)速度產(chǎn)生的壓力為結(jié)構(gòu)粒子的加速度.

1.3 SPH滑坡沖擊的并行優(yōu)化

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

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

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

1.4 SPH滑坡沖擊的地震加速度施加和驗(yàn)證

本研究中土體采用的Drucker-Prager 模型是Bui等［10］改進(jìn)的彈性-完全塑性土體本構(gòu)模型.地震邊坡滑坡研究中土體阻尼系數(shù)和弱化指數(shù)借鑒Chen［22］對(duì)地震邊坡滑坡的研究.驗(yàn)證模型中采用下負(fù)荷屈服面劍橋模型［23］，參數(shù)如表1 所示.計(jì)算模型兩側(cè)邊界設(shè)置為吸能邊界，以此防止地震波傳播到邊界時(shí)，邊界再反射回來對(duì)試驗(yàn)?zāi)Ｐ突蛴?jì)算模型造成影響（如圖2 所示，左側(cè)為剛體邊界，右側(cè)為吸能邊界），以模擬實(shí)際場(chǎng)地的半無限區(qū)域特性，地震波由底部輸入，振動(dòng)方向?yàn)樗?

表1 邊坡振動(dòng)SPH模型參數(shù)Tab.1 Parameters of SPH model slope vibration

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

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

圖3 邊坡振動(dòng)SPH模型Fig.3 SPH model for slope vibration

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

圖4 振動(dòng)后土體粒子分布圖Fig.4 Soil particles distribution after vibration

圖5 振動(dòng)后土體粒子最大剪應(yīng)變?cè)茍DFig.5 Maximum shear strain cloud figure of soil particles after vibration

1.5 SPH滑坡沖擊橋樁結(jié)構(gòu)的數(shù)值模擬工況

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

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

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

2 地震滑坡對(duì)橋樁結(jié)構(gòu)的沖擊規(guī)律

2.1 原波與截取波荷載下滑坡結(jié)果對(duì)比

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

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

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

圖8 42 s Kobe波邊界測(cè)點(diǎn)B1輸出位移時(shí)程Fig.8 Displacement time history of boundary point B1 under 42 s Kobe wave

表2 邊坡振動(dòng)SPH模型參數(shù)Tab.2 Parameters of SPH model slope vibration

圖9 15 s Kobe截取波邊界測(cè)點(diǎn)B1輸出位移時(shí)程Fig.9 Displacement time history of boundary point B1 under 15 s Kobe intercepted wave

圖10～圖13 分別給出了4 個(gè)測(cè)點(diǎn)S1、S2、S3、S4在42 s Kobe 波和15 s Kobe 截取波作用下的水平向和豎向位移對(duì)比.從圖中可以看出S1 測(cè)點(diǎn)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 測(cè)點(diǎn)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 測(cè)點(diǎn)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 測(cè)點(diǎn)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截取波邊坡測(cè)點(diǎn)S1位移對(duì)比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截取波邊坡測(cè)點(diǎn)S2位移對(duì)比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截取波邊坡測(cè)點(diǎn)S3位移對(duì)比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截取波邊坡測(cè)點(diǎn)S4位移對(duì)比Fig.13 Displacement comparison at S4 in slope under 42 s Kobe wave and 15 s Kobe intercepted wave

圖14、圖15 分別給出了兩種地震波作用下的總位移云圖.從圖中可以看出滑動(dòng)土體粒子的位移規(guī)律是相似的，即滑動(dòng)機(jī)理沒有因?yàn)榻厝〉卣鸩ㄝ斎攵霈F(xiàn)與原波差距較大的變化.由兩者的計(jì)算結(jié)果比較可知橋樁結(jié)構(gòu)所受荷載、滑動(dòng)土體位移在兩種地震動(dòng)輸入情況下相差很小，為提高研究效率，選取此15 s 時(shí)段地震加速度作為本小節(jié)數(shù)值模型計(jì)算研究的地震輸入.

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

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

2.2 地震加速度幅值的影響

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

圖16 地震加速度幅值影響計(jì)算模型及測(cè)點(diǎn)布置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峰值加速度作用下土體最終位移和滑動(dòng)的土體體積量很接近（0.2g工況S1 測(cè)點(diǎn)最大水平向位移為34.71 m，最大豎向位移為-13.89 m；0.4g工況S1 測(cè)點(diǎn)最大水平向位移為34.76 m，最大豎向位移為-13.83 m），4種峰值加速度作用下的S1 測(cè)點(diǎn)最終位移均很接近的原因是該測(cè)點(diǎn)與橋樁結(jié)構(gòu)的水平距離為35.28 m，0.6g工況S1測(cè)點(diǎn)最大水平向位移為34.71 m，0.8g工況S1 測(cè)點(diǎn)最大水平向位移為35.46 m，與該測(cè)點(diǎn)和橋樁結(jié)構(gòu)水平距離相差很小，該測(cè)點(diǎn)經(jīng)過地震荷載作用已滑動(dòng)到接近橋樁結(jié)構(gòu)的位置.

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

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

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

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

圖21 給出了4 種峰值加速度作用下S1 測(cè)點(diǎn)水平向和豎向位移-時(shí)間關(guān)系對(duì)比曲線，從圖中可以看出，0.4g峰值加速度作用下土體位移快速增大的時(shí)間要早于0.2g峰值加速度作用下的土體位移快速增大的時(shí)間4.2 s.0.6g峰值加速度作用下土體位移快速增大的時(shí)間要早于0.2g峰值加速度作用下的土體位移快速增大的時(shí)間7.6 s.0.8g峰值加速度作用下土體位移快速增大的時(shí)間要早于0.2g峰值加速度作用下的土體位移快速增大的時(shí)間9.5 s.這說明峰值加速度的提高加速了邊坡的破壞和加大了邊坡的滑坡體積量.監(jiān)測(cè)斷面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%.監(jiān)測(cè)斷面數(shù)據(jù)也說明0.2g峰值加速度工況和0.4g峰值加速度工況滑坡體積量很相近，0.6g峰值加速度工況滑坡體積量稍大于0.2g和0.4g工況滑坡體積量，0.8g峰值加速度工況滑坡體積量最大，且超過0.2g工況滑坡5 倍之多.地震加速度幅值的提高在增加土體邊坡滑坡體積量的同時(shí)，也加快了邊坡土體的滑動(dòng)破壞.

圖21 不同峰值加速度下S1測(cè)點(diǎn)位移-時(shí)間對(duì)比曲線Fig.21 Displacement-time comparison curves at S1 under different PGAs

2.3 地震加速度頻譜的影響

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

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

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

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

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

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

圖26 不同時(shí)間變化系數(shù)下S1測(cè)點(diǎn)位移-時(shí)間對(duì)比曲線Fig.26 Displacement-time comparison curves of measuring point S1 under different tvc

2.4 地震加速度持時(shí)的影響

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

圖27 地震加速度持時(shí)影響計(jì)算模型及測(cè)點(diǎn)布置Fig.27 Calculation model of seismic acceleration duration influence and layout of measuring points

圖20、圖28 分別給出了1 次15 s Kobe 截取地震波和2次疊加截取地震波的位移云圖，圖29、圖30給出了兩種持續(xù)時(shí)間下地震邊坡模型S1、S2 測(cè)點(diǎn)水平向位移、豎向位移與時(shí)間的關(guān)系對(duì)比曲線.從云圖中可以看出，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%，滑坡體積量均很大.從位移-時(shí)間關(guān)系曲線中可以看出，2 次疊加截取地震波與1 次截取地震波作用的邊坡，對(duì)S1 測(cè)點(diǎn)位移的影響差別很小，但對(duì)S2 測(cè)點(diǎn)在15 s（即第1 次截取地震波作用結(jié)束后）出現(xiàn)增大，S2 測(cè)點(diǎn)最終水平位移、豎向位移在1次和2次疊加截取地震波作用下分別為40.53 m、14.72 m 和43.03 m、15.06 m，2 次多于1 次2.5 m、0.34 m.這說明2 次疊加截取地震波對(duì)滑坡體的中部和后緣部分位移影響較大，對(duì)前緣部分無較大影響.

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

圖29 1次和2次截取地震波持時(shí)疊加測(cè)點(diǎn)S1位移-時(shí)間對(duì)比曲線Fig.29 Displacement-time comparison curve of S1 under single and superpositioned wave

圖30 1次和2次截取地震波持時(shí)疊加測(cè)點(diǎn)S2位移-時(shí)間對(duì)比曲線Fig.30 Displacement-time comparison curve of S2 under single and superpositioned wave

2.5 SPH滑坡沖擊模型的并行優(yōu)化效率

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

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

3 結(jié)論

1）基于SPH 方法的基本原理，建立了考慮靜壓力和相對(duì)速度產(chǎn)生的壓力的土-結(jié)構(gòu)粒子相互作用模型，通過地震滑坡模型驗(yàn)證了SPH 方法模擬振動(dòng)臺(tái)試驗(yàn)的準(zhǔn)確性.

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

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

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

5）地震加速度持時(shí)對(duì)滑坡體后緣部分影響較大，對(duì)后緣測(cè)點(diǎn)的位移增加明顯高于前緣測(cè)點(diǎn)，滑動(dòng)到橋墩結(jié)構(gòu)的累積土體體積量隨持時(shí)的增加而增加.

6）對(duì)多種工況的SPH 程序進(jìn)行不同線程數(shù)的并行優(yōu)化，總結(jié)得出線程數(shù)從1 增加到8，加速比呈快速增長趨勢(shì)；線程數(shù)從8 增加到20，加速比增長速度稍緩；線程數(shù)從20 增加到28，加速度增長很緩慢，表明計(jì)算中線程數(shù)為20～28 時(shí)，計(jì)算機(jī)處理效能的利用已接近或達(dá)到最大狀態(tài).