赵宝友,马震岳,梁 冰,徐 伟,3,许新勇

（1.大连理工大学土木水利学院,辽宁 大连 116085;

2.辽宁工程技术大学力学与工程学院,辽宁 阜新 123000;

3.沈阳农业大学水利学院,辽宁 沈阳 110161）

1 引 言

地下结构抗爆性能的研究目前大体有4种方法:现场爆炸原型观测试验、小比例物理模型试验、理论解析和数值模拟。现场爆炸原型观测试验无疑是最直接地研究地下结构抗爆性能的方法,观测数据也最有价值,但该方法不仅耗资巨大,且影响生态环境,为此,迄今为止这方面可借鉴的资料较少。小比例物理模型试验在一定程度上可再现原型试验的效果,相对原型试验,耗资小,在这方面已开展了一定的研究工作[1-6]。顾金才等[1]通过物理模型试验,研究爆炸荷载作用下不同加固方案的对比抗爆效果,研究结果对改进地下抗爆结构加固措施和提高其抗爆能力具有参考价值;杨苏杭等[3]较系统地研究了深埋岩石洞室在爆炸应力波作用下的破坏效应。但地下结构自身的稳定性受很多因素影响,如埋深、地质构造应力、岩体特性等,在考虑上述因素的基础上,要获得较合理的地下结构抗爆性能,需要重复多次小比例模型试验。理论解析方法[7-8]主要是以波动理论为基础,按拟静力法求解具有一定边界条件和简单规则几何形状断面洞室的抗爆性能,对于一些洞室断面形状略微复杂的地下结构,理论解析方法就显得无能为力,不得不求助于近几十年飞速发展起来的有限元理论。以有限元为手段的数值模拟方法不仅可以同时考虑上述多种影响因素,而且模拟精度高和可重复性好[9-12]。赵以贤等[10]采用非线性有限元,从加速度、速度和应力角度研究了爆炸荷载作用下土体中地下拱形结构的动力响应;孙钧[11]系统地阐述了地下抗爆结构有限元计算理论,并指出有限元计算中应注意的若干问题。

已建或将建的水电站地下厂房洞室大都建于岩体较坚硬的岩体内,这类洞室边墙高、跨度大、轴线长。在正常运行期间该洞室可为水利发电服务,在战备时期,是储备物资和人员隐蔽的有利场所。但如此大型的地下岩体洞室群结构,一旦受到恐怖袭击,在强大的常规武器或核武器的地表爆炸荷载作用下,洞室表面岩体或混凝土衬砌结构就会发生损伤,产生裂缝,甚至形成大小不同的岩石或混凝土碎块。这种碎块够携带较大的动能,速度可高达几十米每秒甚至数百米每秒,因而具有很大的杀伤和破坏威力,对洞室结构内部人员和设备的安全构成威胁。因此有必要从力学的角度研究这类大型洞室的抗爆能力,找出洞室最易破坏的部位,进行合理支护,以保证洞室的稳定性,使洞室无论在正常运行期还是在战时都能更好地发挥经济、军事效益。

本文中利用非线性有限元软件ABAQUS[13],建立3维水电站地下主厂房洞室动力分析模型,以洞室围岩和混凝土衬砌结构的损伤为指标,研究在地表强爆炸荷载作用下,洞室结构随埋深、地应力侧压因数和岩体强度的变化的抗爆能力。数值计算中对洞室岩体和混凝土衬砌结构均采用弹塑性软化损伤本构模型。另外,大量数值计算表明[14],对于建于高、中强度岩体内的水电站地下厂房洞室群来说,当相邻洞室间围岩厚度超过相邻洞室中最大洞室宽度的1～2倍时,由开挖引起的洞室群应力效应对相邻洞室动力响应的影响不大,可不考虑这种洞室群效应,而大多数水电站地下岩体洞室群间围岩的厚度大都大于1～2倍的最大洞室宽度;为此,本文中只对单一洞室进行抗爆动力分析。

2 岩体洞室动力数值计算模型和弹塑性损伤本构理论

2.1 岩体洞室动力数值计算模型

洞室几何尺寸为:边墙高65.5 m,拱高9.0 m,跨度32 m,轴线（水平y向）方向长200 m,混凝土衬砌厚0.5 m;洞室四周向外均取150 m岩体,洞室底部取100 m岩体,洞室拱顶一直建至地表,将以上范围内的岩体做为有限域岩体;有限域岩体四周和底部均建立一层50 m厚的岩体做为无限域岩体。图1给出埋深为100 m的洞室的有限域模型。另外,数值计算结果表明,超过一定的埋深后,由于岩体辐射阻尼的影响,地表凸凹不平的山体场地形成的反射波对岩体洞室动力响应的影响较小,为此,假定计算模型地表面为水平,不考虑山体场地效应。

图1 开挖后地下有限域岩体洞室数值模型Fig.1 Numerical model of underground rock cavern in finite field after excavation

2.2 弹塑性损伤本构理论

在实际地下岩体洞室的服务期间,洞室岩体可能会经历如地震等周期性循环动荷载的作用,即岩石经受多次加、卸荷载的作用。然而与岩石初始线性阶段的刚度相比,岩石的后继屈服阶段时的卸荷刚度表现明显的弱化现象。ABAQUS中的弹塑性损伤本构模型[15-16]通过引入损伤因子,来反映混凝土和岩石等准脆性材料在周期荷动荷载作用下后继屈服的损伤机理。本文中采用该模型来模拟混凝土衬砌和洞室周围有限域岩体的动力特性。

根据塑性增量理论,总应变张量ε可分解为弹性部分εe和等效塑性部分εp

混凝土未发生损伤时,在传统塑性力学中混凝土的弹塑性应力应变关系可表示为

式中:σ为总应力,De为弹性刚度矩阵。

混凝土发生损伤时,引入损伤因子描述材料的刚度退化,则损伤后混凝土的应力应变关系可表达为

式中:σef为有效应力;d为刚度退化因子,描述混凝土不同损伤状态的刚度恢复。

周期循环荷载作用下,混凝土力学损伤机理十分复杂,尤其当混凝土由受拉状态转为受压状态时,混凝土的弹性刚度会得到部分恢复,即所谓的“单边效应”。为很好地考虑此效应,损伤变量d假定符合如下关系式

式中:dc和dt分别是压缩和拉伸状态下对应的刚度恢复因子;st和sc是与应力变向有关的刚度恢复下的应力状态的函数,可定义为

式中:wt和wc为刚度恢复权重因子,与材料属性有关,控制着反向荷载下拉、压刚度的恢复。σefi（i=1,2,3）为主应力分量。r（σef）为多轴应力权重因子,是主应力的函数。〈·〉定义为〈x〉=（x+|x|）/2。

此模型的屈服函数中考虑了拉、压荷载作用下混凝土材料的强度演化,以有效应力形式表达如下

该损伤塑性模型应用Drucker-Prager非关联流动准则,流动势G为双曲线函数

式中:φ为高围压下的剪胀角,σt0为单轴拉伸极限强度;Π为双曲线离心率的变量,定义了函数趋近于渐近线的速率。

该模型同样适合于准脆性材料,如岩石等。计算采用的3种强度的岩体和混凝土衬砌的损伤变量dt与开裂位移的曲线和峰后软化段应力与开裂位移曲线分别见图2～3,曲线的确定参照ABAQ US[13]帮助中的方法。其他岩体和混凝土衬砌的物理力学参数见表1,其中E为弹性模量,ν为泊松比,ρ为密度,B为刚度阻尼,φ为剪胀角,σc0为初始压缩强度,σu为极限压缩强度,σt0为极限拉伸强度。

图2 岩体和混凝土衬砌拉伸损伤-开裂位移曲线Fig.2 Tensile damage against cracking displacement for rock and concrete liner

图3 岩体和混凝土衬砌峰后软化段拉应力-开裂位移曲线Fig.3 Post-failure tensile stress against cracking displacement for rock and concrete liner

表1 动力计算中岩体和混凝土衬砌的物理力学参数Table 1 Mechanics characteristic parameters of rock mass and concrete liner for dynamic simulation

2.3 动力边界条件和爆炸荷载

2.3.1 动力边界条件

采用无限元单元作为静、动力边界条件。在静力分析中无限元能提供静力边界条件,动力分析中它能吸收射向边界外侧的地震波,因此无限元单元可以近似模拟无限域地基辐射阻尼效应。

2.3.2 爆炸荷载

选取图4所示的爆炸压应力荷载时程,垂直施加在洞室正上方25 m×26 m（x×y）的水平地表面。荷载最大幅值为60 MPa,升压 0.015 s,动荷载作用0.15 s,有限元动力计算总时间为0.51 s。

3 岩体洞室抗爆性能分析

图4 爆炸压应力荷载时程Fig.4 Pressure history of explosion

数值实现过程:首先在重力作用下,进行一次开挖与衬砌支护的静力分析,形成洞室开挖稳定后的初始应力场,并将模型的位移清零;然后在此基础上,在地表输入爆炸压应力荷载时程,进行动力计算。静、动力计算均考虑模型的几何大变形和材料非线性效应。

从洞室岩体和衬砌的损伤角度来分析埋深、岩体强度和地应力对洞室抗爆性能的影响。

3.1 埋深的影响

建立 50、100、200、400、800和 1 200 m 等 6 种不同埋深的洞室计算模型,每个模型平行轴线方向（y向）的水平侧压因数λ均为1.2,垂直洞室轴线方向（x向）的水平侧压因数λ均为1.0,对有限域岩体采用岩体3的力学参数。以图4的压力时程作为单位参考动荷载,改变其时程的荷载幅值,对每个模型进行多次试算,以洞室混凝土衬砌发生初始损伤为参考指标,得到不同埋深下的洞室进入初始损伤的临界爆炸荷载幅值,见图5。

数值模拟计算得到的临界爆炸荷载在很大程度上体现了不同埋深下洞室的抗爆性能。按曲线斜率的不同,图5的曲线以横坐标100、400 m为分界点可划分为3段,埋深为50～100 m这一段曲线的斜率最小,埋深为100～400 m这一段曲线的斜率居中,埋深为400～1 200 m这段曲线斜率最大。这说明随洞室埋深的增加,使洞室混凝土衬砌发生损伤的外在爆炸荷载的幅值逐渐增加;浅埋深洞室的抗爆性能较差,如埋深为50 m的洞室受到30 MPa（0.50×60 MPa）的地表爆炸荷载就会发生损伤;埋深较深的洞室的抗爆性能较强,如使埋深为400 m的洞室发生损伤,作用在地表的爆炸荷载的幅值要达到1 500 MPa,其抗爆能力是埋深50 m洞室的50倍。这主要与洞室覆岩的强滤波耗能性能有关,洞室覆岩越厚即洞室埋深越深,岩体这种滤波耗能作用越强,进而使深埋深洞室发生损坏的外在地表爆炸荷载幅值越大。

图5 爆洞室发生初始损伤的临界爆炸荷载幅值随埋深的变化Fig.5 Critical explosive pressure amplitude versus embedded depth

3.2 围岩岩体强度的影响

建立3种不同围岩岩体强度下洞室结构的抗爆分析模型,每个模型的洞室的埋深均为400 m,平行轴线方向的水平侧压力因数λ均为1.0,垂直轴线方向的水平侧压力因数λ均为0.8,在模型正上方的地表面均施加45×60 MPa幅值的爆炸荷载。

由于模拟结果以过洞室中轴线的竖直面呈对称分布,为此,只给出对称面一侧的计算结果的云图。图6和图7均为计算结束后的模拟结果的云图分布。分析图6中不同围岩条件下洞室混凝土衬砌的损伤分布和损伤程度可知,随着围岩强度的提高,不仅混凝土衬砌结构的损伤程度减小,而且损伤的面积明显减小,损伤区域也发生显著变化。如围岩材料为岩体1时,大部分拱顶和底板的衬砌均产生不同程度的损伤,最大损伤发生在拱顶,损伤因数为0.95;当围岩材料为岩体2时,洞室底板处的衬砌的损伤消失,只有拱顶处的衬砌发生损伤,且最大损伤因数降至0.6;当围岩材料为岩体3时,只有拱顶沿轴线方向的前后段处的衬砌发生损伤,且最大损伤因数为0.14。

图6 不同岩性下混凝土衬砌结构的拉伸损伤云图Fig.6 Tensile damage contour for liner located in surrounding rocks with different strengths

由图7可知,3种岩体强度下的洞室在爆炸荷载作用后,混凝土衬砌的大部分最大主应力均由静力开挖支护稳定后的压应力变为拉应力,这主要与爆炸地震波在洞室表面发波动生反射所引起的拉伸应力波有关,洞室表面的岩体和混凝土衬砌在强爆炸荷载反射拉伸应力波作用下产生了过大的弹性变形,甚至是不可恢复的非弹性变形,最终导致洞室表面岩体和混凝土衬砌的损坏。图7还表明,混凝土衬砌结构的最大主应力分布规律与损伤相似,即随岩体强度的增大,衬砌最大主应力不仅数值上明显减小,而且最大主应力为拉应力的分布区域也显著减小。围岩材料为岩体1的洞室的拱顶、和底板处的混凝土衬砌均处于拉应力状态,拱顶及侧墙与底板的连接处为较大的拉应力区。随岩体强度的增加,洞室拱顶和底板处的混凝土衬砌拉应力逐渐减小,侧墙衬砌的压应力区域逐渐扩大,可见,高强度的围岩有利于洞室结构自身的稳定性,同时也提高其抗爆的能力。

图7 不同岩性下混凝土衬砌结构的最大主应力云图Fig.7 Maximum principal stress contour for liner located in surrounding rocks with different strengths

另外,数值计算结果还表明,由于地下厂房洞室拱顶效应的影响,再加之爆破冲击波的影响,动力计算过程中,厂房洞室拱肩和底拱角（尤其是洞室的拱肩）附近的剪应力最大,最大剪应力达到12 MPa,严重影响厂房洞室的静动力稳定性。如图8所示埋深400 m洞室的围岩剪应力分布云图和图9的洞室中截面右拱肩某单元剪应力时程曲线更清晰地说明了洞室周边围岩的剪应力分布规律,因此,建议实际的地下厂房洞室应加强拱肩、拱顶和拱角的支护工作。

图8 洞室围岩剪应力云图Fig.8 Shear stress contour for the rock cavern

图9 洞室中截面拱肩剪应力时程曲线Fig.9 Time-history curve of shear stress for the spandrel

3.3 地应力侧压因数的影响

仍以400 m埋深的洞室为例,考虑4种水平两向均相同的侧压因数（λ=0.5,0.8,2.0,3.0）对洞室抗爆性能的影响,洞室围岩材料取岩体3,施加在地表面的爆炸荷载幅值为45×60 MPa。图10～12同样给出计算结束后的1/2对称模型的模拟结果。

如图10所示,小侧压因数λ下,洞室混凝土衬砌的拉伸损伤均发生在洞室拱顶的衬砌之上,拉伸损伤因数均很小,且损伤面积变化不大;随侧压因数λ的增大,尤其是λ＞1后,洞室混凝土结构的损伤程度和损伤区域显著增大。当λ＞1,如λ=2时,由于洞室承受较大的水平侧向压力,进而抵消一部分或者抵消掉爆炸冲击波所产生的拉应力荷载,使洞室拱顶不再产生拉伸损伤,而洞室底板与侧墙连接处附近的衬砌首先产生拉伸损伤;当λ=3时,拉伸损伤扩大到洞室的大部分的底板和侧墙的衬砌之上,这仍与洞室受水平向的侧压力的增大有关。对比分析图10和图11可知,与洞室混凝土衬砌的损伤变化规律类似,当λ＜1时,随λ的增大,洞室周围岩体的损伤程度和损伤区域变化不大,洞室岩体的损伤只发生在洞室底板和侧墙连接处附近;当λ＞1后,随λ的增大,洞室周围岩体的损伤程度和损伤面积显著增加,如λ=3时,整个侧墙均发生拉伸损伤,最大拉伸损伤因数为0.9。

图12中给出了不同侧压因数λ下混凝土衬砌的最大主应力云图,进一步说明了不同侧压因数对地下厂房洞室这类大跨度、高边墙的洞室结构动力稳定性的影响。在小侧压因数下,如λ=0.5,0.8,由于洞室承受竖向自重应力为主,在地表强爆炸荷载的作用下,洞室拱顶和边墙与拱顶、底板连接处的衬砌均出现拉应力,最大拉应力的数值分别为0.64和0.72 MPa,增大幅度不大,且边墙基本处于受压状态;大侧压因数下,如λ=2.0,3.0,加之地表的爆炸荷载联合作用下,洞室边墙的混凝土衬砌出现大面积的拉应力区,最大拉应力的数值分别为1.95和2.12 MPa。

由此可见,在外地表强爆炸荷载的扰动下,较大的侧压力因数虽然有利于洞室拱顶的稳定,但洞室的高边墙却出现了大面积的损伤和拉应力区域,这并不利于水电站地下厂房这类高边墙的静、动力稳定;这也与采用数值计算与现场监测方法研究初始地应力场对钻爆开挖过程中围岩振动的影响所得出的结论[17]一致。在水平侧压力为主的情况下,应考虑适当改变洞室的高跨比,以适应水平为主的地应力条件,进而提高洞室的静、动力稳定性[18]。

图10 不同侧压因数下混凝土衬砌结构的拉伸损伤云图Fig.10 Tensile damage contour for liner under different lateral pressure coefficients

图11 不同侧压因数下洞室围岩的拉伸损伤云图Fig.11 Tensile damage contour for surrounding rocks under different lateral pressure coefficients

图12 不同侧压因数下混凝土衬砌结构的最大主应力云图Fig.12 Maximum principal stress contour for liner under different lateral pressure coefficients

4 结 论

通过数值模拟,研究了洞室埋深、围岩岩体强度和地应力对水电站地下厂房这类高边墙、大跨度、长轴线的大型岩体洞室结构的抗爆性能的影响:

（1）浅埋深的洞室抗爆性能较差,深埋深的洞室抗爆性能较好;在相同地质条件下,埋深增加,洞室抗爆能力增大。

（2）围岩的岩体强度越高,洞室的抗爆性能越强,随着围岩强度的提高,不仅洞室结构的损伤程度减小,而且损伤的面积明显减小,损伤区域也发生显著变化。

（3）当侧压因数小于1时,随侧压因数的增加,洞室的抗爆性能变化并不明显;当地应力侧压因数大于1后,洞室稳定性受水平地应力影响较大,随侧压因数的增加,洞室的抗爆能力显著降低;在水平侧压力为主的情况下,应考虑适当改变洞室的高跨比,进而提高洞室的静、动力稳定性。

（4）由于洞室拱顶效应影响,实际工程应加强对洞室拱肩、拱顶、底脚的加固支护,特别是洞室的拱肩部位。

[1] 顾金才,陈安敏,徐景茂,等.在爆炸荷载条件下锚固洞室破坏形态对比试验研究[J].岩石力学与工程学报,2008,27（7）:1315-1320.

GU Jin-cai,CHEN An-min,XU Jing-mao,et al.Model test study of failure patterns of anchored tunnel subjected to explosion load[J].Chinese Journal of Rock Mechanics and Engineering,2008,27（7）:1315-1320.

[2] 杨苏杭,梁斌,顾金才,等.锚固洞室抗爆模型试验锚索预应力变化特性研究[J].岩石力学与工程学报,2006,25（2）:3749-3756.

YANG Su-hang,LIANG Bin,GU Jin-cai,et al.Research on characteristics of prestress change of anchorage cable in anti-explosion model test of anchored cavern[J].Chinese Journal of Rock Mechanics and Engineering,2006,25（2）:3749-3756.

[3] 陈剑杰.深埋岩石硐室在爆炸应力波荷载作用下的破坏效应[J].岩石力学与工程学报,2001,20（4）:589-589.

CHEN Jian-jie.Failure effect of underground rock openings under stress wave induced by blast[J].Chinese Journal of Rock Mechanics and Engineering,2001,20（4）:589-589.

[4] 唐德高,毕佳,孙乃光.常规武器爆炸作用下地震动的试验研究[J].爆炸与冲击,1997,17（3）:207-213.

TANG De-gao,BI Jia,SUN Nai-guang.Test and study on ground shock caused by conventional weapon explosion[J].Explosion and Shock Waves,1997,17（3）:207-213.

[5] Charles E J,George S.Brick model tests of shallow underground magazines[R].Department of the Army Waterways Experiment Station Corps of Engineers,1992.

[6] 王祥林,郭靖华.冲击荷载下孔洞附近的动力集中[J].地震工程与工程振动,1984,4（4）:50-60.

WANG Xiang-lin,GUO Jing-hua.On the dynamic stress concentration in the neighborhood of a cavity under impulse loading[J].Earthquake Engineering and Engineering Vibration,1984,4（4）:50-60.

[7] 任云燕,张莉,韩峰.地下结构物在爆炸冲击波作用下的动力分析[J].应用数学和力学,2006,27（9）:1122-1128.

REN Yun-yan,ZHANG Li,HAN Feng.Dynamic load analysis of underground structure under the effect of blast wave[J].Applied Mathematics and Mechanics,2006,27（9）:1122-1128.

[8] 盖秉政.弹性波在平面多连通域中的绕射与动应力集中[J].应用数学和力学,1986,7（1）:25-36.

GAI Bing-zheng.Diffraction of elastic waves in the plane multiply connected region and dynamic stress concentration[J].Applied Mathematics and Mechanics,1986,7（1）:25-36.

[9] 赵晓兵,方秦.地下结构变形速度在抗爆结构分析中的作用[J].应用数学和力学,2002,23（4）:436-440.

ZHAO Xiao-bing,FANG Qin.Role of underground structure deformation velocity in the analysis of blast resistant structures[J].Applied Mathematics and Mechanics,2002,23（4）:436-440.

[10] 赵以贤,王良国.爆炸荷载作用下地下拱形结构动态分析[J].爆炸与冲击,1995,15（3）:201-211.

ZHAO Yi-xian,WANG Liang-guo.Dynamic analysis of responses of an underground arc structure to explosive loading[J].Explosion and Shock Waves,1995,15（3）:201-211.

[11] 孙钧.地下抗爆结构有限元数值分析的若干课题[J].岩石力学与工程学报,1983,2（1）:121-129.

SUN Jun.Several problems in the finite element numerical analysis of underground anti-blast structures[J].Chinese Journal of Rock Mechanics and Engineering,1983,2（1）:121-129.

[12] 鞠杨,环小丰,宋振铎,等.损伤围岩中爆炸应力波动的数值模拟[J].爆炸与冲击,2007,27（2）:136-142.

JU Yang,HUAN Xiao-feng,SONG Zhen-duo,et al.Numerical analyses of blast wave stress propagation and damage evolution in rock masses[J].Explosion and Shock Waves,2007,27（2）:136-142.

[13] ABAQUS,Inc.ABAQUS/Standard User's Manual[M].Version 6.6.Providence,RI,2006.

[14] ZHAO Bao-you,M A Zhen-yue.Influence of cavern spacing on the stability of large cavern groups in a hydraulic power station[J].International Journal of Rock Mechanics and Mining Sciences,2009,46（3）:506-513.

[15] Lubliner J,Oliver J,Oller S,et al.A plastic damage model for concrete[J].International Journal of Solids and Structures,1989,25（3）:299-326.

[16] Lee J,Fenves G L.Plastic damage model for cyclic loading of concrete structures[J].Journal of Engineering Mechanics,1998,124（8）:892-900.

[17] 严鹏,卢文波,陈明,等.初始地应力场对钻爆开挖过程中围岩振动的影响研究[J].岩石力学与工程学报,2008,27（5）:1036-1045.

YAN Peng,LU Wen-bo,CHEN Ming,et al.Study on impaction of initial geostress feild on vibration of surrounding rock during excavation with drilling and blasting method[J].Chinese Journal of Rock Mechanics and Engineering,2008,27（5）:1036-1045.

[18] Broch E,Myrvang A M,Stjern G.Support of large rock caverns in Norway[J].Tunneling and Underground Space Technology,1996,11（1）:11-19.