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

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

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

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

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

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

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

1.1 離散單元法

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

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

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

F=Fs+Fn

(1)

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

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

Fn=knUn

(2)

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

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

ΔFs=-ksΔUs

(3)

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

ΔUs=usΔt

(4)

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

總切向接觸力表示為為

(5)

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

Fs,max=μFn

(6)

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

1.2 耦合原理

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

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

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

圖3 耦合邊界力的傳遞

各分量和力矩可以表示為

Fx=fax+fbx

(7)

Fy=fay+fby

(8)

M=-faxya+fayxa-fbxyb+fbyxb

(9)

設(shè)β為差值參數(shù)，則節(jié)點受力可表示為

fax=βFx

(10)

fay=βFy

(11)

fbx=(1-β)Fx

(12)

fby=(1-β)Fy

(13)

其中β由式(14)確定

(14)

2 滾石沖擊砂土墊層數(shù)值模擬

2.1 問題描述

假設(shè)滾石以不同沖擊角度、不同沖擊速度沖擊棚洞砂土墊層，棚洞結(jié)構(gòu)動力響應(yīng)受滾石沖擊角度和沖擊速度的影響，為簡化研究過程，現(xiàn)做出如下假設(shè)。

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

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

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

(4) 不考慮棚洞結(jié)構(gòu)可能發(fā)生的基礎(chǔ)失穩(wěn)。

2.2 物理試驗裝置

為保證數(shù)值模型的正確性，首先開展物理模型試驗對數(shù)值計算結(jié)果進行驗證。圖4為物理模型試驗裝置，鋼架高4.0 m，底部寬2.5 m，頂部寬1.0 m，簡易棚洞結(jié)構(gòu)模型位于鋼架內(nèi)部中心位置，為鋼筋混凝土結(jié)構(gòu)?；炷涟宄叽鐬?.5 m×1.5 m×0.2 m，混凝土板下方4個角為4個混凝土支柱，高0.4 m，長寬均為0.2 m，離板邊界0.2 m，柱與板接觸面安裝壓力傳感器?；炷涟迳细采巴翂|層，厚0.2 m，以直徑為20 cm球形大理石替代下落滾石。試驗時，滾石從鋼架頂部自由落體，沖擊砂土墊層中心位置，力傳感器記錄沖擊過程中支柱所受作用力大小，其相關(guān)參數(shù)如表1所示。

(a) 試驗裝置全貌

表1 傳感器參數(shù)表

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

表2 砂土顆粒級配

2.3 計算模型的構(gòu)建

參照物理模型試驗，按1∶1構(gòu)建數(shù)值計算模型如圖5所示。

(a)

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

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

2.4 參數(shù)取值與校正

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

表3 離散元微觀參數(shù)取值

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

圖6 剪應(yīng)力-剪位移關(guān)系曲線圖

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

棚洞結(jié)構(gòu)采用Drucker-Prager本構(gòu)模型，混凝土材料參數(shù)如表4所示。

表4 混凝土材料參數(shù)取值

以支反力峰值為指標(biāo)，對比試驗與數(shù)值模擬數(shù)據(jù)，擬合校訂所選計算參數(shù)。

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

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

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

圖8 支座反力時程曲線試驗與數(shù)值模擬結(jié)果

3 模擬結(jié)果分析

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

3.1 支座反力響應(yīng)

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

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

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

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

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

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

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

圖10 支座反力峰值與沖擊角度關(guān)系

3.3 峰值沖擊力

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

圖11 滾石沖擊力峰值與角度的關(guān)系

3.4 棚洞頂板中心位移

在沖擊作用下，棚洞頂板發(fā)生變形，研究其變形規(guī)律，分析變形位移值大小可為棚洞結(jié)構(gòu)可靠性提供依據(jù)。

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

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

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

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

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

圖13 頂板位移峰值與沖擊角度的關(guān)系

4 結(jié) 論

基于離散元-有限差分耦合算法開展?jié)L石沖擊棚洞砂土墊層動力響應(yīng)研究，其中采用離散元數(shù)值模擬軟件PFC3D模擬砂土墊層，采用有限差分數(shù)值模擬軟件FLAC3D模擬棚洞結(jié)構(gòu)，開展了不同沖擊角度下棚洞結(jié)構(gòu)的動力力學(xué)響應(yīng)研究，得到結(jié)論如下：

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

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

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