基于SPH-FEM耦合方法的落石沖擊拱形鋼筋混凝土棚洞數(shù)值模擬

柳 春， 余志祥,2， 郭立平， 駱麗茹， 趙世春,2

(1. 西南交通大學(xué) 土木工程學(xué)院， 成都 610031；2. 陸地交通地質(zhì)災(zāi)害防治技術(shù)國家工程實(shí)驗(yàn)室， 成都 610031)

落石多發(fā)于高山峽谷，對鐵路和公路沿線造成了巨大潛在威脅[1]。由于落石具有廣布偶發(fā)和集中群發(fā)的特點(diǎn)，主動(dòng)防護(hù)較為困難[2]。棚洞作為一種被廣泛應(yīng)用的落石被動(dòng)防護(hù)結(jié)構(gòu)，主要適用于落石多、不宜采用集中清理的中小型崩塌災(zāi)害區(qū)域[3]。傳統(tǒng)框架鋼筋混凝土棚洞一般上覆土墊層，通過土墊層耗散沖擊能量(圖1(a))，其安全性得到了實(shí)際工程的驗(yàn)證，但其結(jié)構(gòu)自重大、跨度小的問題仍然未被解決[4]。目前，利用拱形棚洞能解決上述問題(圖1(b))，其結(jié)構(gòu)形式為拱形，不存在大跨度橫梁，解決了結(jié)構(gòu)跨度小的問題，另外由于結(jié)構(gòu)形式的改善和輕質(zhì)緩沖墊層的使用，也解決了自重過大的問題[5]。

目前對鋼筋混凝土棚洞動(dòng)力學(xué)行為研究方法有試驗(yàn)方法、理論方法、數(shù)值方法。Mougin等[6]等采用試驗(yàn)方法研究了滾石對棚洞混凝土板的沖擊。Kishi等[7]對預(yù)應(yīng)力鋼筋混凝土棚洞進(jìn)行了足尺落石沖擊測試，以評估棚洞的極限承載能力。Delhomm等[8]通過在混凝土棚洞支座處增設(shè)耗能減震器(SDR)以替代砂石墊層吸收落石沖擊能量，并進(jìn)行了縮尺試驗(yàn)研究其動(dòng)力學(xué)行為。Kishi等[9]為了優(yōu)化沖擊力的荷載分配，在鋼筋混凝土棚洞上覆三層緩沖層，并進(jìn)行了足尺試驗(yàn)，最大沖擊能量達(dá)到了1 500 kJ。Bhatti等[10]進(jìn)行了拱形棚洞在不同沖擊能量下的足尺試驗(yàn)研究。王爽等通過室內(nèi)模型試驗(yàn)研究了拱形棚洞的抗落石沖擊的力學(xué)性能，并進(jìn)行了不同沖擊能量、不同沖擊位置、不同緩沖層厚度的一些參數(shù)分析。但有限的試驗(yàn)并不能代表復(fù)雜多變的現(xiàn)場落石沖擊條件，若要通過試驗(yàn)全面研究拱棚洞動(dòng)力學(xué)性能，則會(huì)造成試驗(yàn)費(fèi)用較高、周期長。理論研究多集中于落石沖擊力研究。劉茂[11]基于彈塑性赫茲接觸理論提出了落石沖擊力計(jì)算公式。易偉等[12]基于量綱分析法提出了改進(jìn)的落石沖擊力計(jì)算公式。王林峰等[13]根據(jù)消能棚洞的結(jié)構(gòu)特性和受力條件,建立消能棚洞的振動(dòng)模型,推導(dǎo)消能棚洞剛度系數(shù)的計(jì)算公式,進(jìn)而得到消能棚洞的落石沖擊力計(jì)算公式。國外常用沖擊力計(jì)算公式基于半經(jīng)驗(yàn)半理論算法。Montani等[14]以緩沖層的彈性模量、落石速度、質(zhì)量等為主要考慮因素計(jì)算沖擊力，但沒有考慮沖擊體與被沖擊體的材料特性，其計(jì)算局限性大。日本道路協(xié)會(huì)根據(jù)Kawahara等[15]的研究成果制定了相關(guān)落石沖擊力規(guī)范，但其計(jì)算公式的物理量量綱不和諧，其計(jì)算結(jié)果有一定的偏差。目前，通過理論方法獲得物理參量的解析解都存在一定的偏差或地域性。因此數(shù)值模擬方法以其經(jīng)濟(jì)性與高效性逐漸成為此類沖擊問題重要的研究手段，目前主流算法是有限元方法[16-17]，但由于拱棚洞的砂土墊層是由大量離散顆粒體組成的集合體，當(dāng)受到?jīng)_擊時(shí)，會(huì)出現(xiàn)類似流體的大變形性質(zhì)，若用有限元方法計(jì)算，很難處理網(wǎng)格大變形而導(dǎo)致計(jì)算失敗的問題[18]。而離散元方法存在計(jì)算量大和參數(shù)標(biāo)定困難的問題[19]。拱棚洞受沖擊過程中存在材料大變形、高速度、高應(yīng)變率等問題，為研究其動(dòng)力學(xué)行為，則需要一個(gè)合理的數(shù)值計(jì)算方法。

光滑粒子流體動(dòng)力學(xué)方法是一種無網(wǎng)格的連續(xù)介質(zhì)力學(xué)計(jì)算方法，具有很好的自適應(yīng)特點(diǎn)，可以很好地處理大變形和后失穩(wěn)問題，避免了FEM中網(wǎng)格大變形導(dǎo)致計(jì)算失敗，也避免了DEM中計(jì)算量大和參數(shù)標(biāo)定困難等問題。該方法能夠描述砂土材料的大變形特征，并具有較高的計(jì)算精度和穩(wěn)定性，但計(jì)算時(shí)間會(huì)有所增加[20-21]。目前，對于鋼筋和混凝土的模擬，其材料模型及其有限元算法已比較成熟，使用有限元模擬已能保證足夠的精度[22-23]。

(a) 框架棚洞(b) 拱棚洞

圖1 鋼筋混凝土棚洞

Fig.1 Reinforced concrete tunnel

本文提出了有限元法與光滑粒子流體動(dòng)力學(xué)相耦合的方法研究鋼筋混凝土拱棚洞的動(dòng)力學(xué)行為。為提高計(jì)算效率和精度，利用SPH粒子模擬落石沖擊區(qū)域的砂土，非沖擊區(qū)域砂土用有限元單元模擬；將混凝土、鋼筋、巖石、沖擊錘等劃分為有限元網(wǎng)格。然后將SPH與FEM進(jìn)行耦合計(jì)算。將該方法得到的模擬結(jié)果與試驗(yàn)結(jié)果進(jìn)行對比分析，進(jìn)而為拱形棚洞的研究設(shè)計(jì)和應(yīng)用提供一個(gè)合理的計(jì)算方法。

1 計(jì)算理論

1.1 SPH算法基礎(chǔ)

SPH方法是一種無網(wǎng)格粒子方法，通過使用一系列任意分布的粒子來求解具有各種邊界條件的積分方程或偏微分方程，得到精確穩(wěn)定的數(shù)值解。SPH方程的構(gòu)造按兩個(gè)關(guān)鍵步驟進(jìn)行，第一步為核函數(shù)插值，實(shí)現(xiàn)場變量或場變量梯度的插值；第二步為粒子近似，實(shí)現(xiàn)對核函數(shù)估計(jì)積分表達(dá)式的粒子離散[24]。

在SPH方法中,任意宏觀變量f(r)在空間某一點(diǎn)r上的核估計(jì)可以通過函數(shù)f(r)在定義域中的積分獲得

(1)

式中：W(r-r′,h)為核函數(shù)，滿足正則化、緊支性、非負(fù)性、對稱性等條件；h為光滑長度，本文取為1.2d0，d0為相鄰粒子間距離。

為了能夠最終獲得粒子的離散控制方程，通過對支持域內(nèi)一系列粒子的離散化求和來實(shí)現(xiàn)，即粒子i處的場函數(shù)值可以通過核函數(shù)對該粒子緊支域內(nèi)所有粒子的函數(shù)值加權(quán)平均得到

(2)

式中：ρj、mj分別為粒子j的密度和質(zhì)量。

通過式(2)，質(zhì)量、動(dòng)量、能量守恒的離散化方程可相應(yīng)得到。

SPH方法通過使用光滑核函數(shù)進(jìn)行積分表示，它不僅決定了函數(shù)近似式的形式、定義粒子支持域的尺寸，還決定了核近似和粒子近似的一致性和精度。本文采用三次樣條型核函數(shù)來處理三維問題

(3)

1.2 SPH-FEM耦合

在交界面處，有限單元與SPH粒子進(jìn)行耦合。當(dāng)FEM單元與SPH粒子接觸時(shí)，其接觸條件為龍格-庫塔條件

(4)

式中：g為間隙函數(shù)；t為接觸力。

SPH和FEM耦合首先是進(jìn)行接觸搜索，尋找臨近單元，LS-DYNA中的自動(dòng)接觸算法是基于段的接觸搜索方式，在計(jì)算過程中，SPH被定義為從節(jié)點(diǎn)，F(xiàn)EM被定義為主段。采用點(diǎn)面接觸來實(shí)現(xiàn)兩者的耦合，

它是通過罰函數(shù)算法實(shí)現(xiàn)將從節(jié)點(diǎn)的力作用到有限單元上的，罰函數(shù)的基本原理相當(dāng)于在SPH和FEM之間加上法向接觸彈簧限制質(zhì)點(diǎn)穿透主面。在每一個(gè)計(jì)算時(shí)間步檢查各節(jié)點(diǎn)是否穿透主面，若沒有穿透不做處理，若穿透則引入接觸力，其大小與主面剛度和穿透深度成正比，具體計(jì)算方法[25]不在此贅述。有限元單元與光滑粒子界面接觸模型[26]如圖2所示。耦合計(jì)算流程如圖3所示。張志春等[27]借鑒了無網(wǎng)格接觸算法的思想，提出了背景粒子搜索方式，將有限元節(jié)點(diǎn)視為背景粒子，被動(dòng)地被SPH粒子搜索，任何位于有限元節(jié)點(diǎn)支持域內(nèi)的SPH粒子都會(huì)對節(jié)點(diǎn)產(chǎn)生接觸力。本文的接觸搜索方式與背景粒子搜索方式不同，且不需要根據(jù)耦合界面附近有限元單元的尺寸調(diào)整耦合界面附近SPH粒子的光滑長度。

在求解耦合方程時(shí)，SPH采用條件穩(wěn)定的跳蛙顯式積分方法，F(xiàn)EM也采用條件穩(wěn)定的中心差分法，兩者耦合要求兩者積分必須同步，這就要求兩者在同一計(jì)算框架下每一步計(jì)算采用相同的計(jì)算步長，時(shí)間步長ΔtSPH-FEM取兩者的較小值(式(5))。

ΔtSPH-FEM=min(ΔtSPH,ΔtFEM)

(5)

式中：SPH時(shí)間步長ΔtSPH=βh/c; FEM時(shí)間步長ΔtFEM≤Lmin/c,其中c為材料聲速，β為時(shí)間步長比例系數(shù)，Lmin為最小單元尺寸。

圖2 SPH與FEM接觸

圖3 SPH-FEM接觸算法流程

2 拱棚洞沖擊試驗(yàn)

2.1 試驗(yàn)?zāi)Ｐ兔枋?/h3>

圖4為鋼筋混凝土拱棚洞足尺沖擊試驗(yàn)?zāi)Ｐ?通過吊車提升質(zhì)量為10 000 kg沖擊器到2.5 m，5 m，10 m，20 m來分別獲得250 kJ, 500 kJ, 1 000 kJ, 2 000 kJ的沖擊能量。沖擊器為一個(gè)組合體，上半部分為圓柱體，高度0.95 m，直徑為1.25 m，下半部分為一個(gè)被切割后的球體，高度為0.3 m，直徑為2 m(圖5(a))。其中砂墊層被裝在一個(gè)木箱中，高度0.9 m，跨度10 m，沿道路軸向5 m。拱形棚洞斷面的幾何尺寸和鋼筋分布位置如圖5所示，其中，D19表示直徑為19 mm的鋼筋。其他具體鋼筋連接設(shè)計(jì)細(xì)節(jié)參考日本規(guī)范[28]?？v向方向沿道路長為6 m。

圖4 拱棚洞沖擊試驗(yàn)?zāi)Ｐ蚚10]

(a) 斷面尺寸圖

(b) 鋼筋分布圖

2.2 測試方法

沖擊力時(shí)程曲線用一個(gè)量程為500g，采集頻率為5 000 Hz，安裝在沖擊器內(nèi)部的加速度計(jì)采集。位移計(jì)安裝在拱形棚洞的拱圈內(nèi)表面，測量其豎直方向位移，量程和采集頻率分別為500 mm和1 000 Hz。

3 計(jì)算模型

3.1 計(jì)算單元

LS-DYNA在分析結(jié)構(gòu)受沖擊時(shí)的幾何非線性、邊界非線性、材料非線性計(jì)算上有著獨(dú)特的優(yōu)勢。本文采用其進(jìn)行耦合計(jì)算，計(jì)算模型如圖6所示。被沖擊砂土采用SPH粒子，被沖擊區(qū)域尺寸按經(jīng)驗(yàn)調(diào)試確定，取為2 m(x)×2 m(y)×0.6 m(z)，相鄰粒子間距0.05 m；非沖擊區(qū)域的砂土采用8節(jié)點(diǎn)完全積分六面體單元，單元尺寸0.1 m，在這里將砂土劃分為兩個(gè)區(qū)域主要考慮以下兩點(diǎn)：一個(gè)是SPH粒子計(jì)算開銷大，非沖擊區(qū)域沖擊變形小，可用有限元網(wǎng)格模擬；另一個(gè)是由于SPH的光滑核函數(shù)不具有Kronecher delta 張量的性質(zhì)，且邊界粒子本身存在缺陷，不適合作為邊界?；炷痢r土、沖擊器采用6節(jié)點(diǎn)/8節(jié)點(diǎn)完全積分的實(shí)體單元?；炷羻卧W(wǎng)格及其臨近的巖土網(wǎng)格尺寸為0.1 m，為提高計(jì)算效率，遠(yuǎn)端巖土網(wǎng)格尺寸較大。鋼筋采用桁架單元模擬[29]。劃分網(wǎng)格后，共有22 919個(gè)SPH粒子，219 848個(gè)實(shí)體單元，13 249個(gè)桁架單元。

圖6 拱棚洞計(jì)算模型(mm)

3.2 材料模型

表1所示為各材料的基本材料特性。各材料的本構(gòu)模型如下所述。

砂：砂土材料是一種多孔介質(zhì)材料，在強(qiáng)沖擊作用下，會(huì)表現(xiàn)出可壓縮泡沫材料的力學(xué)特性，為了準(zhǔn)確描述砂土在沖擊作用下的力學(xué)行為，利用了式(6)的拋物線描述其應(yīng)力-應(yīng)變曲線(圖7(a))?？捎肕AT_CRUSHABLE_FOAM進(jìn)行描述。

表1 各種材料特性

(6)

式中：σsand為應(yīng)力，單位為MPa，εsand為體積應(yīng)變。

鋼筋：采用各向同性的彈塑性本構(gòu)進(jìn)行模擬，塑性硬化模量取為彈性模量的0.01倍(圖7(c))，初始屈服強(qiáng)度為338 MPa，采用Von Mises屈服準(zhǔn)則?？捎肕AT_PLASTIC_KINEMATIC進(jìn)行描述。

巖體：巖體采用彈性體來模擬，可用MAT_ELASTIC描述。

沖擊器：由于沖擊器幾乎不變形，為加快計(jì)算效率，沖擊器采用剛體模擬，可用MAT_RIGID描述。

(a) 砂

(b) 混凝土

(c) 鋼筋

圖7 主要材料本構(gòu)模型

Fig.7 Main material constitutive models

3.3 模型邊界

沖擊器與砂土SPH粒子間定義為侵蝕接觸；砂土有限單元與砂土SPH粒子間定義為點(diǎn)面接觸；砂土有限單元與混凝土單元采用面面接觸；砂土有限單元與巖土單元采用面面接觸；混凝土單元與巖土采用共節(jié)點(diǎn)連接；鋼筋單元嵌固在混凝土單元內(nèi)由于鋼筋為小應(yīng)變，可用CONSTRAINED_LAGRANGE_IN_SOLID來定義鋼筋與混凝土的約束。

巖土底部采用全約束邊界，巖土側(cè)壁可采用無反射邊界條件。沖擊器與砂土接觸時(shí)開始計(jì)時(shí)總的沖擊時(shí)間定義為200 ms。

4 結(jié)果分析

4.1 數(shù)值計(jì)算與試驗(yàn)對比

選取沖擊能量為1 000 kJ的沖擊錘沖擊力時(shí)程曲線、拱中點(diǎn)位移時(shí)程曲線與數(shù)值模擬的結(jié)果對比，如圖8所示，可看出，沖擊力時(shí)程曲線包括快速上升、快速下降、穩(wěn)定三個(gè)階段，可視其為一個(gè)脈沖荷載，兩者的時(shí)程曲線具有較好地吻合性。選取試驗(yàn)與數(shù)值模擬的沖擊器沖擊力峰值、拱中點(diǎn)位移峰值進(jìn)行對比，如表2所示。

通過對比發(fā)現(xiàn)，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果吻合較好，沖擊力峰值和拱中點(diǎn)位移峰值最大誤差均沒有超過10%。SPH-FEM方法能很好地模擬沖擊器沖擊拱棚洞的整個(gè)過程。隨著沖擊能量的增大，沖擊力峰值和拱中點(diǎn)位移峰值也逐漸增大。數(shù)值模擬結(jié)果與試驗(yàn)有一定的誤差，數(shù)值模擬的沖擊力峰值整體稍偏大，拱中點(diǎn)位移峰值稍偏小，主要由以下幾點(diǎn)原因造成：SPH方法忽略了粒子間的摩擦耗能；砂土的緩沖性功能與密實(shí)度有關(guān)，所用材料本構(gòu)無法考慮到砂土的密實(shí)度影響。雖然雙線性混凝土材料本構(gòu)簡便，但混凝土本構(gòu)非常復(fù)雜，需用試驗(yàn)確定，造成拱中點(diǎn)位移時(shí)程后半段有一定出入，但總體趨勢一致。

(a) 沖擊力時(shí)程

(b) 拱中點(diǎn)位移時(shí)程

表2 試驗(yàn)結(jié)果與仿真結(jié)果對比

4.2 沖擊過程現(xiàn)象分析

圖9所示，為沖擊錘沖擊砂土的過程，沖擊能量為1 000 kJ，當(dāng)開始沖擊錘與砂土接觸時(shí)，砂土由于受到擠壓和剪切作用，開始向四周流動(dòng)(T=30 ms),此時(shí)砂土應(yīng)力較??；隨著沖擊錘與砂土接觸面積的增大，砂土被沖擊錘擠向四周，在砂土表面形成一個(gè)微隆起區(qū)，內(nèi)部區(qū)域形成了一個(gè)砂坑，在T=125 ms時(shí)，沖擊錘位移達(dá)到最大值0.825 m，此時(shí)砂土應(yīng)力最大值為12.74 MPa，對應(yīng)的體積應(yīng)變達(dá)到0.505，說明砂土體在沖擊過程中存在著大變形；在T=200 ms時(shí)，沖擊過程趨于靜止，靜止后沖擊位移為0.804 m，卸載恢復(fù)的彈性變形很小，此時(shí)砂土體卸載導(dǎo)致其應(yīng)力也變小。整個(gè)沖擊過程與試驗(yàn)現(xiàn)象有較好的一致性。

T=0 ms

T=30 ms

T=125 ms

T=200 ms

圖9 沖擊器沖擊砂土過程

Fig.9 The process of penetration by impactor

4.3 耗能分析

圖10為沖擊能量1 000 kJ的沖擊工況的能量時(shí)程圖。在沖擊過程中，沖擊器的動(dòng)能和勢能轉(zhuǎn)化為砂墊層的內(nèi)能、摩擦耗能、棚洞內(nèi)能等，砂墊層耗散了很大一部分能量?；颇Σ聊苤饕蓻_擊過程中沖擊器與砂墊層的摩擦所致。棚洞和巖石幾乎沒有耗散沖擊耗能，說明塑性變形很小，砂墊層起到了很好的保護(hù)作用。從表3可看出，隨著沖擊能量的增大，土墊層耗能逐漸減小。這是由于沖擊能量較小時(shí)，鋼筋混凝土和巖石處于彈性階段，不耗散沖擊能量，隨著沖擊能量增大，鋼筋混凝土開始出現(xiàn)塑性變形，耗散了一部分能量。但砂土墊層耗能仍然占總初始沖擊動(dòng)能的85%以上，它是一種很好的緩沖耗能材料。砂土墊層使得棚洞實(shí)現(xiàn)了分級(jí)耗能抗沖擊，小沖擊能量，棚洞不損傷，大沖擊能量，棚洞不倒。

圖10 能量時(shí)程圖

表3 砂土墊層耗能比例

5 結(jié) 論

SPH-FEM方法充分結(jié)合有限元與光滑粒子流體力學(xué)的優(yōu)勢，能夠模擬大變形材料。本文通過SPH-FEM方法對落石沖擊鋼筋混凝土拱棚洞進(jìn)行了全過程數(shù)值模擬，并與試驗(yàn)結(jié)果對比，得到以下結(jié)論：

利用SPH-FEM耦合方法得到的沖擊錘沖擊力時(shí)程、拱中點(diǎn)位移時(shí)程趨勢與試驗(yàn)結(jié)果較為吻合。隨著沖擊能量的增大，沖擊力峰值和拱中點(diǎn)位移峰值也逐漸增大。數(shù)值模擬的沖擊力峰值和拱中點(diǎn)位移峰值最大誤差均沒有超過10%，證明了該算法的準(zhǔn)確性。

利用SPH-FEM耦合方法形象地再現(xiàn)了砂土成坑的現(xiàn)象。砂墊層幾乎耗散了所有的沖擊能量，砂土墊層是一種很好的緩沖耗能材料。砂土墊層使得棚洞實(shí)現(xiàn)了分級(jí)耗能抗沖擊，小沖擊能量，棚洞不損傷，大沖擊能量，棚洞不倒。

SPH-FEM耦合方法在模擬大變形緩沖層材料方面具有可行性和有效性，為鋼筋混凝土拱棚洞的設(shè)計(jì)提供了一種數(shù)值手段。