基于離散元-有限差分耦合的滾石沖擊棚洞墊層動力響應研究

王東坡， 劉 浩， 裴向軍， 孫新坡， 周良坤， 劉彥輝

(1.成都理工大學 地質災害防治與地質環境保護國家重點實驗室，成都 610059；2.四川輕化工大學 土木工程學院，四川 自貢 643000;3.中鐵第一勘察設計院集團有限公司，西安 710043)

滾石沖擊棚洞過程是一個瞬態動力學過程，往往涉及到結構大變形和復雜的能量轉換，目前的理論模型難以定量描述沖擊的動力響應，各個學者對沖擊過程進行了一些列研究[1-3]。盡管結構原型沖擊試驗可以得到真實的沖擊過程數據，且正在被越來越多的研究者所重視，然而，由于物理模型試驗投資巨大，周期長，相關的研究成果還較少。近年來，伴隨數值模擬方法的發展，越來越多的學者采用數值模擬手段研究滾石沖擊棚洞的動力響應。何思明等[4]采用有限元程序ABAQUS研究了滾石沖擊荷載下棚洞結構動力響應問題，包括沖擊接觸壓力與接觸位移的時間關系、不同沖擊角度時的彈坑形狀、墊層材料對沖擊壓力分布的影響等。楊璐等[5]利用有限元軟件對滾石沖擊過程進行數值模擬,研究了不同速度和入射角下滾石對棚洞的沖擊作用。王東坡等[6]建立有限元模型分析了滾石沖擊荷載下棚洞鋼筋混凝土板動力響應特征，并研究了EPS(expanded Polystyrene)墊層及雙層泡沫鋁夾心板等新型墊層材料的耗能減震作用[7-8]。裴向軍等[9]采用動力有限元對滾石沖擊砂土墊層進行數值仿真計算，分析不同沖擊能量下多組砂土墊層厚度組合的動力響應及耗能緩沖機理。王玉鎖等[10]基于模型試驗，建立數值模型開展落石沖擊下無回填土拱形明洞破壞特征及失效模式、極限承載力、落石沖擊荷載及極限狀態表達式等的研究。張群利等[11]借助ANSYS/LS-DYNA有限元軟件模擬棚洞結構在落石沖擊荷載下的動力響應過程，通過分析比較類棚洞結構的受力與變形的特征，研究不同結構類型的抗沖擊性能。

上述的數值模擬研究多采用有限元方法進行，該方法可較好的描述棚洞頂板、梁、柱等鋼筋混凝土結構，然而針對棚洞頂板的緩沖墊層，則難以有效的模擬砂土顆粒的離散性，且無法考慮顆粒間復雜的力的傳導及顆粒的位移，不能準確的描述其沖擊荷載下的響應特征。為此，針對該問題，采用離散元-有限差分耦合算法對墊層和棚洞兩部分分別進行模擬，充分發揮兩種模擬方法各自優勢，開展滾石沖擊棚洞墊層的研究。PFC-FALC是常用的離散-連續介質耦合軟件，Cai等[12]采用PFC-FLAC耦合的方法研究了日本某地下硐室巖體聲發射特征；Song等[13]利用耦合手段模擬了隧道的開挖變形；Jia等[14]采用PFC-FLAC耦合方法，揭示了粒狀土在動態壓實過程中的宏觀和微觀機制。該耦合算法發展較為成熟，但針對滾石沖擊棚洞頂板問題還鮮有使用。

棚洞結構的設計多以滾石自由落體垂直沖擊棚洞結構為依據。然而，滾石往往斜向沖擊棚洞頂板[15]，其棚洞結構的動力響應與垂直沖擊相比會有顯著不同， 為此，文章擬采用離散元-有限差分耦合算法開展棚洞結構在滾石不同沖擊角度下的動力力學響應研究，為棚洞工程設計提供參考依據。

1 離散元-有限差分耦合算法

1.1 離散單元法

假設砂土墊層顆粒間無黏聚力，視為散體結構，采用基于離散單元法的顆粒流程序(particle flow code，PFC)來模擬砂土墊層及滾石的動力力學行為。

顆粒間的接觸剛度模型采用線性接觸模型。在接觸模型中，顆粒接觸點的接觸力和相對位移可分解為沿法向和切向的分量，由法向剛度和切向剛度通過力-位移定律分別將法向力和法向位移、切向力和切向相對位移相聯系[16]，如圖1所示。接觸點合力為F，表示為

圖1 接觸模型示意圖(據Su等[17]的研究)

F=Fs+Fn

(1)

式中,Fs,Fn分別為接觸力的切向和法向分量。

顆粒之間的法向位移由法向重疊量Un表示，則接觸力法向分量Fn可表示為

Fn=knUn

(2)

式中,kn為接觸點的法向剛度。

顆粒流模型中，以增量形式表示剪切力的計算，在顆粒集合體模型生成時，接觸力切向分量為0，后續計算過程中，每一計算時步內顆粒位移的變化引起接觸力切向分量的變化逐一累加到前一計算時步的數值上，切向力計算時步內的增量表示為

ΔFs=-ksΔUs

(3)

式中：ks為顆粒接觸點的切向剛度;ΔUs為每個計算時步內切向位移增量。由式(4)計算

ΔUs=usΔt

(4)

式中：us為相對切向速度;Δt為時間步長增量。

總切向接觸力表示為為

(5)

砂土顆粒不存在黏聚力，無法向抗拉強度，顆粒在其抗剪強度范圍內發生滑動，顆粒間發生滑動的辨別條件為

Fs,max=μFn

(6)

式中,μ為顆粒間的摩擦因數，當兩顆粒的摩擦因數不相等時取小值。若Fs>Fs,max，則顆粒間發生滑動，滑動后的切向接觸力取Fs,max。

1.2 耦合原理

PFC-FLAC耦合模擬計算采用基于邊界控制墻體的方法[18]，PFC,FLAC分別從宏觀、細觀上模擬連續域、離散域內介質的力學行為，在連續域和離散域的接觸邊界二者相互耦合，并借助Scoket O/I接口進行不同域間計算數據的傳輸與轉換,如圖2所示。

圖2 PFC-FLAC耦合計算原理(據石崇等研究)

FLAC連續域與PFC離散域的接觸面指定為PFC的墻單元，墻上的接觸力和彎矩采用等效方法分配到墻面頂點，而這些頂點附著在FLAC單元網格點上，因此墻面頂點與實體單元網格點同步運動，這些力參與連續域中的分析，連續域網格點的變動也帶動墻面頂點的變動，進而將位置和速度信息傳遞到離散域中的顆粒。PFC-FLAC耦合邊界的力傳導,如圖3所示。Fx,Fy分別為作用在PFC墻體上力的分量，M為合力矩;fax,fay和fbx,fby分別為對應的FLAC單元節點a、節點b上力的分量。

圖3 耦合邊界力的傳遞

各分量和力矩可以表示為

Fx=fax+fbx

(7)

Fy=fay+fby

(8)

M=-faxya+fayxa-fbxyb+fbyxb

(9)

設β為差值參數，則節點受力可表示為

fax=βFx

(10)

fay=βFy

(11)

fbx=(1-β)Fx

(12)

fby=(1-β)Fy

(13)

其中β由式(14)確定

(14)

2 滾石沖擊砂土墊層數值模擬

2.1 問題描述

假設滾石以不同沖擊角度、不同沖擊速度沖擊棚洞砂土墊層，棚洞結構動力響應受滾石沖擊角度和沖擊速度的影響，為簡化研究過程，現做出如下假設。

(1) 將滾石視為質量分布均勻的剛性球體，密度為2 000 kg/m3。

(2) 砂土墊層視為粒徑不一的剛性球形顆粒。

(3) 滾石沖擊位置為墊層頂部中心處。

(4) 不考慮棚洞結構可能發生的基礎失穩。

2.2 物理試驗裝置

為保證數值模型的正確性，首先開展物理模型試驗對數值計算結果進行驗證。圖4為物理模型試驗裝置，鋼架高4.0 m，底部寬2.5 m，頂部寬1.0 m，簡易棚洞結構模型位于鋼架內部中心位置，為鋼筋混凝土結構。混凝土板尺寸為1.5 m×1.5 m×0.2 m，混凝土板下方4個角為4個混凝土支柱，高0.4 m，長寬均為0.2 m，離板邊界0.2 m，柱與板接觸面安裝壓力傳感器。混凝土板上覆砂土墊層，厚0.2 m，以直徑為20 cm球形大理石替代下落滾石。試驗時，滾石從鋼架頂部自由落體，沖擊砂土墊層中心位置，力傳感器記錄沖擊過程中支柱所受作用力大小，其相關參數如表1所示。

(a) 試驗裝置全貌

表1 傳感器參數表

所選砂土孔隙率0.38，顆粒級配如表2所示。

表2 砂土顆粒級配

2.3 計算模型的構建

參照物理模型試驗，按1∶1構建數值計算模型如圖5所示。

(a)

考慮到計算機性能，設定砂土墊層顆粒粒徑8～13 mm，在棚洞頂板中間1.0 m×1.0 m×0.2 m內共生成25 680個顆粒，四周以墻單元進行圍欄，模擬物理模型試驗砂箱，防止顆粒逸散，對棚洞頂板上部中心處豎向位移和支柱頂端壓力進行監測。

棚洞頂板與支柱采用共節點連接，最小劃分單元2 cm×2 cm×2 cm，共計58 476個單元。

2.4 參數取值與校正

砂土顆粒間、砂土與墻體間采用線性接觸模型，模型的微觀參數取值如表3所示。

表3 離散元微觀參數取值

砂土直剪物理試驗與直剪的數值模擬剪應力-剪位移曲線,如圖6所示。

圖6 剪應力-剪位移關系曲線圖

圖6反映了不同法向應力下物理試驗和數值計算的剪應力-剪位移關系，圖中S為試驗值，M為模擬計算值。從圖6可看出100 kPa和200 kPa法向應力下試驗曲線與數值計算得到的曲線能夠很好的擬合，在300 kPa和400 kPa下雖部分區段兩條曲線有明顯的差距，但整體的趨勢和數值均相差較小，因此，總體擬合結果較好，所選參數能夠反映試驗砂土性質。

棚洞結構采用Drucker-Prager本構模型，混凝土材料參數如表4所示。

表4 混凝土材料參數取值

以支反力峰值為指標，對比試驗與數值模擬數據，擬合校訂所選計算參數。

圖7(a)、圖7(b)分別表示支反力峰值隨墊層厚度和下落高度的變化。隨著下落高度的增大，支反力峰值明顯加大，試驗條件下，下落高度從2 m增加到4 m，支反力峰值從3.9 kN增加到6.4 kN，每米的增加幅度分別為18%和39%，總的峰值沖擊力增加了64%。數值計算的結果每米的增加幅度分別為21%,33%，總的峰值沖擊力增加幅度為72%，與試驗數據差值分別為4%，6%和8%，在下落高度3.5 m時曲線出現拐點，這是由于下落高度較大時，高度的增大導致速度增大的幅度更大，使得支反力峰值更大。不同墊層厚度下的支反力峰值，其數值計算結果與物理試驗結果差值最大不超過0.3 kN，誤差均小于10%，表明數值計算的結果能較好的反應支反力峰值的大小及變化趨勢，所取計算參數適宜。

(a) 下落高度2 m,墊層厚度擬合

圖8為墊層厚度30 cm、滾石下落高度2 m條件下支座反力時程曲線試驗與數值模擬結果，二者在峰值力、響應周期及各時刻的支座反力上能夠較好擬合，證明離散-有限差分的方法適用于該問題的分析。

圖8 支座反力時程曲線試驗與數值模擬結果

3 模擬結果分析

以沖擊方向與水平面夾角為沖擊角度，研究不同沖擊角度下棚洞結構的動力響應特征，選取角度分別為30°,45°,60°,75°和90°，每個沖擊角度下沖擊速度分別取20 m/s,30 m/s和40 m/s，共計15組研究工況。

3.1 支座反力響應

圖9為不同沖擊速度及不同沖擊角度下支座反力與沖擊時間的關系曲線。隨著滾石接觸棚洞墊層，開始對棚洞結構產生沖擊作用，支座反力短時間內急劇增大，達到最大值后迅速減小為0并轉為負值，隨后快速增大至第二個峰值后開始減小，直到達到下一時刻正的峰值，如此往復。整個過程中，支座反力交替表現為壓和拉兩種作用方式的力，混凝土板在這個過程中不斷振蕩耗能至支座反力逐漸減小近似為0。

由圖9可知伴隨沖擊速度的增大，不同沖擊角度下的支座反力峰值均呈現出明顯增大的趨勢，且出現峰值支座反力的時間隨之提前。在沖擊速度較小時，5個不同沖擊角度的支座反力峰值相差較小，并且小角度沖擊時，峰值支座反力的差值較大，而當沖擊角度較大時差值較小。如當沖擊角度為75°和90°時，峰值支座反力差值均在3 kN以內，這是由于75°時速度在豎直方向的分量較大，接近總速度，因而所造成的峰值支座反力相差不大。

(a) 沖擊速度為20 m/s

相同沖擊速度時，較小沖擊角度下曲線震蕩時間及幅度較小，支座反力更快的趨于0。特別是沖擊角度為30°時，曲線波動小并且震蕩時間短，各沖擊速度下平衡時間分別為95 ms，103 ms和112 ms，這是由于當沖擊角度為30°時，速度的水平分量較大，豎直分量小，對墊層的沖擊作用弱。

但沖擊速度較大時，小沖擊角度比較大沖擊角度平衡時間更長。沖擊速度為40 m/s時，90°沖擊角度下約160 ms曲線震蕩結束，達到平衡狀態；而當沖擊角度為75°時，曲線的震蕩時間約為180 ms，所需平衡時間更長。

3.2 支座反力峰值隨角度變化

圖10為不同沖擊速度下左右兩側(即圖5中1、2監測位置)支柱峰值支座反力與沖擊角度關系曲線。相同沖擊角度和沖擊速度下左右兩側支柱峰值支座反力相差較小。較小沖擊角度下，支座反力峰值隨速度變化較小，隨著沖擊角度增大，支座反力峰值越來越大，沖擊速度的影響也愈加凸顯。沖擊角度為30°時速度每增加10 m/s，支反力峰值增加量分別為6.3%和13.0%，在90°時則變為56.0%和16.3%，可見在小沖擊角度時沖擊速度的增大對支座反力影響較小，而在較大沖擊角度，這種影響相對變強。

圖10 支座反力峰值與沖擊角度關系

3.3 峰值沖擊力

圖11為各沖擊速度下滾石沖擊力峰值隨沖擊角度的變化關系。沖擊力峰值與沖擊角度呈正相關，隨著沖擊角度的增大，滾石峰值沖擊力增大的趨勢愈發明顯。在較小沖擊速度下，如當沖擊速度20 m/s時沖擊角度的增大對峰值沖擊力的提升不大，當沖擊角為30°時峰值沖擊力為19.6 kN，而當沖擊角度為90°時峰值沖擊力為36.94 kN，增大約1.8倍；而在沖擊速度為40 m/s時，增大了2.1倍，可見在較大沖擊角度下，角度對沖擊力的影響更加顯著。圖中曲線斜率隨沖擊速度增大而增大，因此相對于沖擊角度，沖擊速度對沖擊力峰值的影響更加顯著。

圖11 滾石沖擊力峰值與角度的關系

3.4 棚洞頂板中心位移

在沖擊作用下，棚洞頂板發生變形，研究其變形規律，分析變形位移值大小可為棚洞結構可靠性提供依據。

圖12為不同沖擊速度和不同沖擊角度下棚洞頂板中心處豎向位移值與時間關系曲線，同支座反力曲線相似，頂板中心處位移值隨沖擊角度和沖擊速度增大而增大，且沖擊速度、沖擊角度越大，出現峰值的時間越早。

(a) 沖擊速度為20 m/s

當沖擊角度為30°時曲線波動幅度很小，且沖擊發生后很快便達到平衡狀態，證明在整個沖擊過程中頂板中心處只發生小幅度的位移形變，并且很快便不再發生振蕩，達到平衡。小沖擊角度下，沖擊角度的提升對位移影響顯著，隨著沖擊角度的增大，沖擊角度的變化對位移值的影響逐漸減小。

值得注意的是當沖擊速度不大時，不同沖擊角度下的位移曲線變化較一致，特別90°和75°沖擊角度下的位移-時間曲線變化有較好的同步性，這是由于二者速度的豎直分量接近，且75°的水平分量比較小。在大沖擊速度下，沖擊角度在75°,60°,45°和30°的都有較大的橫向速度分量，90°的橫向速度為0，導致位移-時間曲線有較大不同。

圖13給出了不同沖擊速度下混凝土頂板豎向峰值位移隨沖擊角度的變化。隨著沖擊角度和沖擊速度增大，頂板峰值位移也逐漸增大，當沖擊角度為30°、當沖擊速度20 m/s時最小，為0.52 mm；當沖擊角度為90°、沖擊速度40 m/s時最大，達到3.2 mm。各點連線在較小角度段呈小幅上凹，在較大角度段上凸，表明較小沖擊角度下，角度對頂板峰值位移影響較大，而沖擊角度較大時，其影響減弱。

圖13 頂板位移峰值與沖擊角度的關系

4 結 論

基于離散元-有限差分耦合算法開展滾石沖擊棚洞砂土墊層動力響應研究，其中采用離散元數值模擬軟件PFC3D模擬砂土墊層，采用有限差分數值模擬軟件FLAC3D模擬棚洞結構，開展了不同沖擊角度下棚洞結構的動力力學響應研究，得到結論如下：

(1) 對比數值模擬與物理試驗在不同下落高度和不同墊層厚度下的支座反力峰值，二者差值小于0.5 kN，誤差小于10%，數值計算與物理模型試驗結果較為吻合，說明采用基于離散元-有限差分耦合算法的數值計算手段研究滾石沖擊棚洞墊層是可行的。

(2) 沖擊角度和沖擊速度對支座反力以及頂板中心位置豎向位移有顯著影響。隨著沖擊角度和沖擊速度的增加，支座反力峰值和頂板位移峰值逐漸增大；且沖擊速度越大，棚洞的動力響應越明顯，達到峰值支座反力和峰值頂板位移的時間越短；小沖擊角度下沖擊角度對支反力峰值和頂板位移峰值影響較大，而沖擊角度較大時，沖擊速度的影響更加顯著。

(3) 滾石沖擊在頂板中心位置時，同一沖擊角度下不同位置的支柱豎向支座反力大小近似相等。