王东坡， 刘 浩， 裴向军， 孙新坡， 周良坤， 刘彦辉

(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) 滚石冲击在顶板中心位置时，同一冲击角度下不同位置的支柱竖向支座反力大小近似相等。