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

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

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

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

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

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

1.1 離散單元法

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

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

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

F=Fs+Fn

(1)

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

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

Fn=knUn

(2)

式中,kn為接觸點(diǎn)的法向剛度。

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

ΔFs=-ksΔUs

(3)

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

ΔUs=usΔt

(4)

式中：us為相對(duì)切向速度;Δt為時(shí)間步長(zhǎng)增量。

總切向接觸力表示為為

(5)

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

Fs,max=μFn

(6)

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

1.2 耦合原理

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

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

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

圖3 耦合邊界力的傳遞

各分量和力矩可以表示為

Fx=fax+fbx

(7)

Fy=fay+fby

(8)

M=-faxya+fayxa-fbxyb+fbyxb

(9)

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

fax=βFx

(10)

fay=βFy

(11)

fbx=(1-β)Fx

(12)

fby=(1-β)Fy

(13)

其中β由式(14)確定

(14)

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

2.1 問題描述

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

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

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

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

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

2.2 物理試驗(yàn)裝置

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

(a) 試驗(yàn)裝置全貌

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

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

表2 砂土顆粒級(jí)配

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

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

(a)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 模擬結(jié)果分析

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

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

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

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

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

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

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

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

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

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

3.3 峰值沖擊力

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

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

3.4 棚洞頂板中心位移

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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