地下洞室地震分析论文

1引言

随着水电事业的发展，在我国西南山区已建、在建和规划了大批地下厂房洞室群。该区域位于板块相交处，为地震高发区，大型地下洞室群的抗震围岩稳定问题受到各方的普遍重视。随着计算机数值计算方法的发展，采用动力时程法模拟地下洞室群的地震过程成为研究洞室群抗震围岩稳定的重要途径之一。其中，岩体的动力损伤本构模型是洞室地震灾变过程模拟的基础。钱七虎等、李海波等在岩石动力学方面，林皋等在坝体混凝土动力特性方面，葛修润等在岩石疲劳破坏方面均做了大量基础试验研究，为研究岩石的动力本构模型积累了丰富的试验资料和理论分析基础。但目前尚无普遍认可和较成熟的适用于地下洞室围岩地震分析的本构模型。鉴此，本文在考虑动荷载作用下岩体强化特性和循环加载情况下岩体疲劳损伤特性的基础上，提出一种考虑岩体动力强化和损伤劣化的动力本构模型，并嵌入到自主开发的显式有限元计算程序中，获得了一些有益的结论。

2岩体弹塑性损伤动力本构模型

2．1弹性本构由已有的材料动力特性试验研究结果，可得“围岩弹性模量随应变率的增加而增加”的定性结论。在静力弹性模量的基础上，围岩动弹性模量E0可表述为：E0＝P（ε）E（1）式中，P（ε）为大于1的应变率函数；E为围岩静弹性模量。P（ε）受材料、围压、试件尺寸等的影响，需根据实际工程试验确定。在工程设计前期无试验确定P（ε）时，从偏安全角度考虑，可取P（ε）＝1．0；对混凝土材料可根据规范［8］建议值，取P（ε）＝1．3。在显式有限元计算中，弹性本构采用剪切模量G和体积模量K描述。根据式（1），动剪切模量G0和动体积模量K0可表示为：G0＝E0／21＋μ［（）］＝P（ε）E／21＋μ［（）］（2）K0＝E0／31－2μ［（）］＝P（ε）E／31－2μ［（）］（3）式中，μ为泊松比。围岩破坏主要以剪切、拉裂破坏为主，而静水压力造成的破坏较少。结合损伤模型，弹性参数受损伤程度的影响可表述为：G＊0＝（1－D）G0（4）K＊0＝K0（5）式中，G＊、K＊分别为考虑损伤的动剪切模量和动体积模量；D为损伤系数。2．2塑性强度准则在带拉伸截止限的Mohr－Coulomb静力屈服准则（图1）基础上，提出考虑应变率和损伤系数对材料粘聚力值影响的动力强度屈服准则。假定压应力为负，拉应力为正。线AB为剪切屈服准则fs＝0，线BC为拉伸屈服准则ft＝0，屈服表达式为：fs＝σ＊1－σ＊3Nφ＋2Tc（ε）槡Nφft＝σ＊3－σt烅烄烆（6）其中Nφ＝1＋sinφ1－sinφ（7）式中，σ＊1、σ＊3分别为损伤后的有效第一、第三主应力；φ为内摩擦角；T为损伤影响系数；c（ε）为与应变率相关的岩体粘聚力；σt为极限抗拉强度，最大值σtmax＝c（ε）／tanφ。二阶张量形式下，有T＝1－D21＋D22槡＋D23（8）式中，D1、D2、D3均为各主应变方向损伤系数，由损伤演化方程求解。文献［3］、［9］的研究表明应变率对材料内摩擦角的影响较小。因此，式（6）、（7）不考虑应变率和损伤系数对围岩内摩擦角的影响。c（ε）是动力屈服准则与静力屈服准则的主要区别。对重大工程，应根据现场试验确定。对于中、小型工程没有条件进行现场试验时，可采用理论推导表达式确定。摩尔－库仑屈服准则采用主应力表示为：σ1－σ32＝σ1＋σ32sinφ＋ccosφ（9）式中，σ1、σ3分别为第一、第三主应力。由此可得，在单轴压缩情况下粘聚力可表示为：c＝σY21－sinφcosφ（10）式中，σY为单轴压缩时的强度极限。欧洲国际混凝土委员会（CEB）在统计众多混凝土动力加载试验结果的基础上，给出混凝土动抗压强度与加载应变率的推荐公式［10］为：fc＝fcsεd／εs（）1．026αεd≤30s－1（11）其中α＝5＋3fcu／（4）－1式中，fc、fcs分别为动态、静态压缩强度；εs为静态应变率，取εs＝3×10－6s－1；εd为动态应变率，介于3×10－6s－1与30s－1之间；fcu为混凝土静态立方体抗压强度。将式（11）应用到洞室围岩中，fcs、fcu取为围岩的抗压强度。考虑到洞室围岩受解理、裂隙、结构面的切割，其强度随应变率的提高比单纯岩块小。将式（11）中的指数1．026α改为1．0α，并与已有的花岗岩、大理岩等试验结果［11］进行对比，见图2。图2抗压强度与加载应变率关系图Fig．2Compressivestrengthandstrainraterelationalgraph由图2可看出，式（11）经过修正后基本能合理描述岩石抗压强度与加载应变率的关系。将式（10）中的σY用式（11）中的fc替代，可得：c＝fcsεd／εs（）1．0α21－sinφcosφεd≤30s－1（12）对于通过直剪试验确定的岩石粘聚力cs，若认为单轴压缩与直剪试验测试的粘聚力相差不大，则式（12）可转换为：c＝csεd／εs（）1．0αεd≤30s－1（13）式（6）、（7）与式（12）或式（13）组成洞室围岩的考虑应变率影响的摩尔－库仑屈服准则。在显式有限元时步计算中，每计算时步内，T、c（ε）为常数，按照围岩非关联性，屈服准则相应的势函数为：gs＝σ＊1－σ＊3Nψgt＝σ＊3｛（14）其中Nψ＝1＋sinψ1－sinψ（15）式中，ψ为围岩的剪胀角。

2．3损伤模型

围岩在地震灾变中的损伤破坏是指随着地震波在围岩介质中的传播，由于张力和剪切作用导致材料颗粒间距增大，围岩中微裂纹产生、扩展，最终导致围岩宏观裂纹形成失稳破坏。为描述这种现象，本文引入二阶损伤张量，并给出损伤判别准则和损伤演化方程。围岩的损伤判据表述如下。在三维情况下，单元应力状态进入塑性，则发生剪损；单元第一主拉应变超过极限拉应变，•83•则发生拉损；单元第一主压应变超过极限压应变，则发生压损。即：D＝0fs≤0；ε3＜εl；ε1＜εyh（εij）fs＞0l（εij）ε3≥εly（εij）ε1≥εy烅烄烆（16）式中，hεij（）为三维剪损演化方程；lεij（）为三维拉损演化方程；yεij（）为三维压损演化方程；fs为剪切屈服准则；ε1、ε3分别为围岩第一、第三主应变；εl、εy分别为围岩极限拉应变和极限压应变；第一主应力为压，第三主应力为拉，第一主应变为压，第三主应变为拉。由于围岩极限剪切应变测定较为困难，根据文献［10］，将围岩应力进入塑性状态作为剪损发生的门槛值。对实际工程中未测试围岩的极限拉、压应变的情况，在本文程序中用岩石抗拉、压强度Rs和弹性模量E的比值代替［12］，即［ε］＝Rs／（kE）（17）式中，Rs为岩石抗压强度；k为安全系数。损伤破坏是围岩微裂隙发展和微观孔隙率增加的过程，是不可逆的。围岩的损伤主要取决于其应变的发展。如式（16）损伤判据，当围岩应力进入塑性状态，即发生剪损破坏，围岩应力和物理参数都随累积的塑性偏应变的增大而减小，累积塑性变形越大，围岩损伤程度越高。同理，围岩进入拉损或压损状态，围岩应力和物理参数都随应变的增大而减小。在围岩材料失稳破坏前的相邻阶段内微裂纹损伤增加是非常迅速的，对于这种损伤破坏的发展特征通常可以用指数函数来描述。三维剪损演化方程为：Dhi＝hεi（）＝1－exp－Rεpiεp槡i（）i＝1，2，3（18）式中，Di、εpi分别为第i个主应变方向上的剪损系数、塑性偏应变；Dhi的上标h表示剪切应变引起的损伤；R为材料的损伤常数。考虑到体积应变不产生剪切破坏，式（18）中采用塑性偏应变计算剪损系数。同理，三维拉、压损演化方程可表示为：Dl3＝lε3（）＝1－exp－Rε3－［εl（］）ε3－［εl（］）（19）Dl1＝lε1（）＝1－exp－Rε1－［εy（］）ε1－［εy（］）槡（）（20）式中，Dl3、Dl1的上标l表示拉、压应变引起的损伤。由于剪切损伤系数、拉裂损伤和压裂损伤系数均采用主应变计算，方向上保持一致，在同时发生情况下，可对损伤系数进行求和计算。为便于在后处理中显示围岩损伤状态，在结果文件输出时，将二阶损伤张量转换为一阶损伤标量D＝D21＋D22槡＋D23／3，表示单元整体损伤状态。Frantziskt等［13］在描述岩石中的混凝土结构材料出现裂隙时，提出由于结构受力后产生微裂隙，在微裂隙区应力将被释放和减小，从而形成应力损伤区。该情况同样适用于洞室围岩应力损伤状态。将应力损伤程度用等效应力的形式表示为：σ＊＝σ1（1－D1h）000σ21－D2（h）000σ3（1－D3h）熿燀燄燅（21）式中，h为修正系数；受拉情况下取h＝1．0，受压情况下取h＝0．2［14］。

3实例

我国西南某水利枢纽工程位于金沙江下游，场区地震基本烈度为Ⅷ度。该枢纽左岸地下厂房洞室群由主厂房、主变洞、尾调室、引水洞、尾水洞、母线洞及交通洞等附属洞室组成，规模庞大。其中，主厂房尺寸310．00m×31．80m×84．20m（长×宽×高）。岩体以厚层灰岩为主，岩体类别为Ⅱ、Ⅲ类。初始地应力场第一主应力为11～20MPa，第三主应力为6～10MPa，属中高等地应力场。左岸地下厂房平行布置6台水轮发电机组，单机功率850MW。采用本文的计算理论对该地下厂房洞室群的抗震稳定进行分析，以指导该洞室的抗震支护设计。

3．1计算条件

（1）计算模型。计算模型涵盖左岸地下厂房洞室群主要建筑物，共划分188906个六面体单元。模型坐标原点位于＃6机组中心线，X轴与厂房轴线方向垂直，指向下游为正，计算范围382．30m；Y轴与厂房轴线方向重合，从＃6机组指向＃1机组为正，计算范围471．20m；Z轴与大地坐标重合，指向上为正，计算范围345．07m。模型网格最大单元特征尺寸10．00m，以最大网格单元特征尺寸不大于1／10波长计，按照最小波速1484．78m／s计算，适用计算地震波最高频率为15Hz。模型边界采用自由场人工边界。计算模型见图3。

（2）计算参数。动力时程计算的洞周应力场采用有支护静力开挖完成后的洞周围岩扰动二次应力场。岩体材料参数见表1。主洞室的支护参数主要为喷层＋锚杆＋锚索。锚杆Φ32＠1．5m×1．5m，L＝6m／9m，锚索2000kN＠4．5m×4．5m，L＝25m。

（3）地震波荷载。对相关部门批复的地震波加速度时程进行滤波和基线校正后，计算采用的地震波时程见图4、5。岩体阻尼比取5％。主要考虑顺水流向和竖直向地震作用，其中竖直向取顺水流向地震波的2／3。图4地震波加速度时程Fig．4Seismicwaveaccelerationtimehistory图5地震波位移时程Fig．5Seismicwavedisplacementtimehistory

3．2洞周围岩破坏区发展规律

地震过程中，洞室群围岩受地震波附加应力场的影响，岩体应力和应变处于波动状态。洞周围岩塑性区、开裂区体积时程见图6、7。由图6可看出，地震过程中洞周围岩塑性区体积处于波动状态，塑性体积较大时段主要集中于2～8s，这与地震波位移时程第一个峰谷波段基本吻合（图5）；后续虽有同样的峰谷波段（8～14s），塑性区体积却维持在较低量值波动。说明地震过程中洞周围岩塑性破坏主要发生在第一个峰谷波段。图6洞周围岩塑性区体积时程Fig．6Surroundingrockplasticvolumetimehistory图7洞周围岩开裂区体积时程Fig．7Surroundingrockcrackingvolumetimehistory由图7可看出，随着地震过程的进行，围岩开裂区体积处于逐渐增加的状态。这主要是因为地震过程中围岩在地震波的作用下表现为循环加卸载过程，弹性变形恢复，而塑性变形则累积增加，致使围岩总体应变逐渐增加，并逐步超出围岩的极限应变，从而使得洞周围岩开裂体积逐渐增大。由图6、7的对比也可看出这一点，图6中在13．7s之前洞周围岩已经历过多个加卸载循环，洞周围岩塑性变形得到充分的累积，致使13．7s后洞周围岩总应变逐渐达到拉裂破坏标准，围岩开裂体积开始急剧增加。说明洞周围岩的开裂破坏主要受地震波的持时影响，即使峰值较低的长持时地震，也可能造成洞室围岩失稳。因此，对于低峰值长持时的地震，采用拟静力法计算未考虑地震波持时，可能未反映其主要地震特性。随着地震波的传播，洞周围岩破坏状态处于动态调整过程，但任何破坏形式产生的岩体损伤会逐步累积增加。洞室开挖完毕损伤较严重的区域主要分布在尾调室多洞交口处，其中上游岩柱底部、尾调室间下部隔岩等部位较严重，损伤系数达到0．045；地震作用完毕后，损伤较严重区域分布规律与洞室开挖完毕时基本相同，损伤系数最大达到0．170，增加278％。说明洞室开挖中损伤较严重的部位在地震作用下会更加严重。因此，在抗震设计中开挖损伤破坏较大部位依然是抗震支护的重点区域。地震后，主厂房最大开裂区深度由6．6m增加到8．5m，主变洞由3．7m增加到5．7m，尾闸室由5．7m增加到8．7m。说明地震对洞室群围岩稳定有一定程度影响，但开裂区深度均在已有支护措施的控制范围内，洞室群围岩基本稳定。

4结论

a．在已有大量岩体动力试验研究的基础上，考虑地震动作用下岩体的强化特性和循环加载情况下岩体的疲劳劣化特性，提出了一种岩体弹塑性损伤动力本构模型，并将该模型应用到我国西南某在建大型地下厂房洞室群的动力时程计算中，得到一些有益结论。

b．地震过程中洞周围岩塑性破坏主要发生在第一个峰谷波段，洞周围岩的塑性破坏受地震波峰值影响较大。

c．地震过程中洞室周围岩体的开裂区范围呈逐渐增加的趋势。洞周围岩体的开裂破坏受地震波的持时影响较大，即使峰值较低的长持时地震，也可能造成洞室围岩失稳。

d．围岩在开挖中损伤较严重的区域，在地震过程中其损伤程度可能会成倍增加。在抗震设计中开挖损伤破坏较大部位是抗震支护的重点区域

作者：张志国 牟春来 邵年 单位：长江勘测规划设计研究有限责任公司