离散元-边界元动力耦合模型在地下结构动力分析中的应用

摘要
利用离散元-边界元耦合模型对溪洛渡工程地下厂房洞室群静、动力响应进行了分析，对离散元模拟地下结构的一些问题进行了探讨，得到了一些结论：离散元模型能够很好地模拟地下洞室群的变形，通过与边界元的耦合，可以模拟辐射阻尼的影响，由于辐射阻尼的影响，地下洞室群的地震响应比地面结构小很多。
关键词：离散、辐射阻尼、地下结构、动力分析
离散元-边界元动力耦合模型[1]将二维可变形体离散元[2]与二维动力时域边界元[3]两种计算方法耦合起来，充分发挥离散元与边界元的优点，将非连续体的模拟与无限介质辐射阻尼的模拟统一到一个模型中，为地下结构的抗震稳定分析提供了全新的手段。本文利用这一耦合模型对溪洛渡工程地下厂房洞室群进行了静、动力分析，对离散元模拟地下结构的一些问题进行了探讨，得到了一些结论。
1. 工程简况[4]
溪洛渡水电站在长江上游干流金沙江上，是一座以发电为主，兼有拦沙、防洪、漂木、航运等综合效益的巨型水电站。坝址区位于溪洛渡峡谷中段，山体浑厚，谷坡陡峻，地形完整。经国家地震局烈度评定委员会审定，地震烈度为VIII度。厂房发电系统拟定为库区内地下厂房，左、右岸对称布置，且地下主厂房、主变室和尾水调压室均呈"品" 字形分布；由上游至下游依次是地下厂房、主变室和尾水调压室。两岸地下厂房尺寸规模相同，跨度32m，长414m，高74.7m，底板高程335.5m。厂房区岩体由含斑玄武岩和玄武质角砾熔岩组成，岩体一般新鲜较完整，无大的断层切割，层间错动带一般不发育，产状平缓，挤压紧密，一般干燥无水。层内错动带在局部层内相对较发育，错动带一般宽5～10cm，挤压紧密，为岩块岩屑型。裂隙以陡裂和缓裂为主，中裂一般不发育。厂房区地下水活动微弱，厂房区地应力属中等地应力区。
2. 计算条件和简化假定
2.1 计算断面的选取
限于时间和条件，目前开发的离散元-边界元动力耦合分析软件只能进行二维模拟，因此必须选取适当的二维计算断面。以右岸厂房为例，初步选定通过12号机组中心线的厂房断面，得到一个向两侧无限延伸的断面，如图1所示。然后分别再划定边界元和离散元的离散范围。离散元的离散范围为以洞室群为基准，向下延伸长度大于最大洞室高度，向两侧延伸约两倍最大洞室高度，即图中折线ABCDEA所包围的区域。折线FABCDG以下的半无限平面则采用边界元方法来模拟，并假定为各向同性的连续介质，AF和DG均为水平线，根据计算经验，各用15个70m的边界单元离散，并在F和G处截断。折线ABCD用24个边界单元离散，每个单元长度大约为70m，整个地表FAEDG的边界条件均为自由表面。计算中对洞室断面作了相应的简化，对于采用离散元模拟的ABCDEA区域，必须进行块体的划分和材料基本参数的确定。根据成勘院[4]提供的工程地质条件，首先确定选取断面内的主要裂隙的参数，根据这些参数可以有计算机自动生成裂隙并以此来对计算区域进行块体划分。岩体中实际存在的裂隙往往是间断性的，存在岩桥，而目前离散元的研究对象必须是被裂隙完全切割的各种形状的块体。因此必须对实际的裂隙进行必要的处理，将较短的裂隙进行延长，使它们相交，将整个岩块比较彻底地切割成各种块体，以便进行计算。为保证如此处理后仍能模拟原有的岩体，将根据裂隙延长的长度对计算的力学参数如摩擦系数f、粘滞力c、接触弹性系数Kn、Ks等进行加权修正，以保证等效的抗滑力、弹性力相同。最后形成的离散块体见图2，图2中所示阴影部分表示的块体为耦合模型分析所需的界面块[1]。为清晰起见，图2中未示出各离散元块体内部的差分网格，在洞室群附近差分网格的大小均小于4m，离洞室群稍远处网格的大小均小于10m，上部山体的网格大小均小于40m，整个计算模拟的块体数大于500块。整个完成一次静、动力计算在Pentium III/450上约需80~100小时。
2.2 地应力及其模拟方法
地应力是影响围岩稳定的一个重要因素。溪洛渡地下厂房区存在着较高的初始地应力，在计算模型中应该充分予以考虑。在采用耦合模型进行数值分析中，针对溪洛渡工程的具体情况，可以在竖直方向，通过增大重力加速度，即施加一个大于g的重力加速度，相当于给岩体施加附加的竖向荷载，有效地模拟竖向地应力，再通过施加水平力来挤紧岩块，来模拟水平方向的应力，这样，通过调整施加的重力加速度及水平力的大小，可以使垂直地应力和水平地应力均达到预期要求的数值。
本文主要考虑地下厂房洞室群的稳定性，因此计算中主要考虑洞室群及其附近区域的地应力场。实际的地下厂房洞室群地应力场数据已经由国家电力公司成都勘测设计研究院提供，它是根据实测的地应力加以反演推得的。经过试算，计算模型中在水平方向共施加4.5×106kN压力，在竖直方向施加1.15g的重力加速度。
2.3 渗流场的模拟
水库蓄水后，水位抬高，会在两岸山体内形成渗流场。地下洞室开挖以后，水会由洞壁渗漏到洞室中，地下水将在洞室周围形成局部汇流区。地下渗流场的存在对洞室的稳定会产生一定的影响，因此计算模型中必须考虑渗流场。
耦合模型中主要考虑静水头的影响，即根据各个节点的静水头确定作用在节点上的渗透压力，目前，静水头暂时按渗流面考虑，未能充分模拟洞室本身的排水效应。当然，将来也可按渗流等压面考虑静水头。地震动水压力的影响由于时间所限，未能引入模型，也未予模拟。本文分析中的渗透水压是根据河海大学渗流实验室提供的《溪洛渡电站地下厂房渗流控制计算工作汇报》的渗流场施加的。由于在水库蓄水前，地下水位很低，水库蓄水后，地下水位上升，洞室才开始受渗透水压力的影响。因此，渗透压力的施加是在开挖完成，且离散元系统计算达到新的稳定以后。
2.4 地下厂房系统开挖顺序
由于耦合模型分析是仿真分析，施工开挖顺序会对分析结果产生较大影响，计算中必须模拟施工开挖的顺序。在施工过程中，溪洛渡地下厂房系统的开挖有一定的先后顺序，须按开挖施工进度表的安排进行。在主厂房开挖的过程中，同时还要进行主变室和调压室的开挖及支护。但实际的开挖施工进度安排比较复杂，且本研究的重点又是动力稳定，故在计算中进行了一定的简化，将洞室群的开挖简化为三步进行，如图3所示，其中I、II、III分别代表第一、二、三步开挖的部分。计算时每一步开挖完成后给予系统一定的时间，使之达到稳定状态，以便进行下一步计算。
3. 施工期与蓄水期分析
3.1 模拟过程
　　离散元计算的初始状态假定块体间没有相互作用力。计算开始以后，在自重和地应力（以及可能的地下水渗透压力）作用下，块体之间会产生相互作用力并引起变形，直至达到平衡状态，可称为"固结过程"。因此，采用耦合模型进行计算时，首先要经过一个固结过程，直到系统到达稳定后，才开始施工模拟计算。施工过程分三步开挖来模拟，每步开挖均要在达到稳定后再进行下一步开挖。施工开挖模拟结束后，再施加渗透压力荷载，模拟蓄水过程，为简化起见，蓄水过程假定为一步。
3.2 变形分析
　　取一典型a点（见图3），耦合模型计算的该点在施工过程中的水平位移曲线如图4a所示（图中的时间没有实际意义），可以分为三个阶段。ab段为从计算初始状态开始的稳定阶段（固结阶段），即达到施工前初始地应力的状态；bc阶段为开挖阶段，共分三次完成开挖；cd阶段为蓄水阶段。图4b为图4a的bd段局部放大图，表示了每一步洞室开挖以后，a点发生的回弹变形。由于该点位于洞室侧壁的中部，第一步只开挖洞室的上部，因此第一步a点的回弹位移较小，为0.3mm；第二步开挖到洞室的中部，这时a点产生了较大的回弹位移，为10.4mm；第三步开挖洞室的下部，这一步开挖量较大，a点产生了很大的回弹位移，大小为29.5mm。水库蓄水以后，由于上游侧水位高于下游侧，岩体受到渗压的作用，也会产生新的位移，运动趋势为向下游方向。
　　其他各点也有类似规律，值得注意的是主变室的位移。它的开挖在第二步一次完成，它的回弹位移也一次完成，规律同主厂房和调压室。但在第三步开挖主厂房和调压室的下部时，受调压室开挖的影响，它会产生向下游的水平位移。取其洞壁上的b点（见图3）进行分析，在第二步开挖中水平回弹位移大小为4.7mm，由于调压室所影响，它在开挖第三步产生了向下游方向的位移，大小为20.0mm，影响是很大的。因此洞室群的开挖要十分注意洞室开挖的相互影响。
三个洞室因开挖产生的总的回弹位移大小如图5所示，最大位移发生在调压室的侧壁上，最大回弹位移值为53mm。另外主厂房下游侧底部也存在向下游的变形，这也是主变室和调压室开挖的影响造成的。
4. 地震响应分析
　　在开挖与蓄水过程静力分析的基础上，利用离散元-边界元耦合模型分析了地震波波动输入条件下，各洞室的地震响应。地震荷载采用从无限远处波动输入，峰值加速度是采用溪洛渡工程的地震危险性分析提供的峰值加速度0.18g，地震加速度反应历程由抗震规范规定的标准反应谱反演得到，地震历时24.46s，见图6，通过数值积分并进行归零调整得到地震位移历程，波动输入的幅值按该地震波在半平面表面的响应为0.18g反演，在SV波竖向输入时，在无限远处输入的地震波动幅值为0.09g，此时的地震输入与通常的水平地震输入等价。
4.1 地表地震响应
　　首先分析在地震波输入时，山体表面的地震响应。计算结果显示山体地表各点的峰值加速度略有不同，图1中山顶E点附近的加速度较大，峰值加速度为1.94 m/s2，为0.18g的1.1倍；河谷底部A点附近的加速度较小，峰值加速度为1.36 m/s2，为0.18g的0.77倍。这一分布与大量地震实际观测的现象吻合。地表峰值加速度平均值为1.63 m/s2，为0.18g的0.925倍，略小于0.18g，这与离散块体的阻尼和地震分布有关，总的来说，地表的地震加速度与设计地震相当。
4.1 洞室周边各点的地震响应
　　图7为耦合模型分析得到的a点（见图3）地震加速度响应历程，它是由计算得到的位移历程通过数值微分得到的，可以注意到在25s附近有一较大峰值，这是由于地震时程只有24.46s记录，突然终止造成的虚假响应。其他各点的地震加速度与之相似。图8和图9分别是采用离散元模型与耦合模型计算的各洞室周边各点的加速度峰值与设计地震峰值的比值分布图，其中，离散元的计算网格与图2中耦合模型的离散元网格系统，仅是没有与边界元耦合，界面块体为刚性块，地震输入是在这个刚性块上给定设计地震来施加的，因此它们之间的差别显示了辐射阻尼和地震波动输入的影响。从图中可以分析出以下结论：
(1) 常规的离散元分析中，各洞室的地震响应与输入地震相近，最大响应的放大倍数为1.64，发生在调压室中下部。
(2) 由于辐射阻尼和地震波动输入，耦合模型计算的各洞室的地震响应与地表响应及常规离散元分析的地震响应相比有非常明显的下降，最大响应仅有设计地震的0.357倍，这与实际地震中地下结构地震反应较小的实际观测结果相吻合，说明在溪洛渡地下厂房洞室群这样埋深较大的地下结构，考虑辐射阻尼和地震波动输入能够更好地符合实际情况。
(3) 耦合模型计算中，由于地震波从下向上入射，洞室底部是迎波面，加速度响应较大，而顶拱的加速度响应较小。残余位移的分布呈上大下小的分布，最大残余位移发生在调压室顶部，最大残余位移值为10mm。地震残余位移较小，对洞室群的围岩稳定威胁较小。
4.2 地震波不同入射方向的影响
考虑到地震波的输入方向受到多种因素的影响，难以准确确定，从工程安全的角度出发，对地震波的不同入射方向进行了分析，给出了SV波与竖直方向夹角为30°和-30°入射的结果。可以分析出以下结论：
(1) 不同入射角工况条件下，洞室周边的地震响应略有不同。
(2) 在各种入射角条件下，一般迎波面的响应比背波面大，在地震从下向上的入射条件下，洞室底部的加速度响应较大，而顶拱的加速度响应较小。
(3) 残余位移的分布仍呈上大下小的分布。各种输入条件下的残余位移略有不同，其中30°输入条件下，不仅加速度响应较大，残余位移也大，最大值为29mm，也发生在调压室顶部。对地下洞室群的整体围岩稳定威胁不大。
5. 主要结论
　　采用离散元-边界元动力耦合模型对溪洛渡工程地下厂房洞室群的静、动力响应进行了分析，得出以下主要结论：
　　1． 开挖过程中，由于应力回弹作用，地下厂房会发生一定的变形。变形的最大值发生在调压室侧壁的中段，最大值为53mm。
　　2． 开挖过程中，三个洞室之间会产生相互的影响，因此在实际施工时要合理安排开挖顺序。
　　3． 利用上述模型分析溪洛渡工程地下厂房洞室群，得出了地下洞室群地震响应明显小于地表响应的结论，与地震实际观测的一般结论吻合。
　　4． 不同的地震入射方向，对计算结果有一定的影响，在目前的条件下，建议在工程分析中应该计算可能的不同入射方向，选取最危险的情况进行设计。
　　5． 溪洛渡工程地下厂房洞室群的围岩完整性较好，具有较高的抗震稳定能力，在设计地震作用下发生围岩整体失稳破坏的可能性不大。