高水平放射性废物(简称高放废物)具有放射性强、毒性大、半衰期长的特点,对其进行最终安全处置难度极大,面临一系列的科学、技术和工程挑战^[1]。目前,深地质处置是国际上普遍认可的长期处置高放废物的可行方法,即将高放废物埋置在500~1 000 m深度范围内稳定的地质体中^[2]。高放废物地质处置库工程具有建设条件复杂、安全等级高、服务期限长(数万年计)等特点,研发难度很大。为此,为了进行场址特征评价、开发和测试处置概念和技术,获取处置库围岩热—水—力—化学—生物—放射性的耦合特征,并完成工程屏障,评估其长期性能和安全性,世界范围内陆续建设了多个地下实验室,如瑞典的Äspö、加拿大的Whiteshell、芬兰的ONKALO和法国的Meuse/Haute Marne地下实验等^[3⇓⇓-6]。高放废物地下实验室是高放废物最终处置不可或缺的一环。

北山地下实验室南侧的河西走廊和青藏高原北缘构造环境复杂,强震多发。如公元180年甘肃高台西7¹/₂级地震, 1932年甘肃玉门昌马7.6级地震。近几年在青藏高原北缘也时有发生,如2021青海玛多M 4.7级地震、2022年门源M 6.6级地震。这些区域潜在强震的发生可能对地下实验室场址的稳定性产生影响。然而,由于地下岩石硐室的抗震性能比地面结构要好,导致针对大型洞室群地震反应及其抗震性能的研究不多,自1995年日本阪神地震中出现地下结构严重的震害后,才开始重视地下结构的抗震性能研究^[7]。如杜修力等^[8]、刘晶波等^[9]、Yuan等^[10]和Zhang等^[11]分别对土质地基中地下结构的地震反应特征、地震破坏机理和抗震性能水平等进行了全面且系统地研究,何川等^[12]、崔臻等^[13]、Hashash等^[14]和Corigliano等^[15]分别对岩质地基中隧道的抗震性能和动力反应特性进行了研究,主要考察了近断层的影响。然而,已有对岩质地基地下洞室群的地震响应影响的研究较为有限,尤其是对考虑复杂岩体构造的花岗岩硬岩场区大型地下洞室群整体抗震性能的研究较为缺乏。

根据我国“十四五”的核工业布局和规划,我国已于2021年7月在甘肃北山开工建设高放废物地下实验室,而相关地下工程在设计、施工等过程中均需考虑抗震问题^[16]。同时,高放废物地下实验室由多个非规则地下结构所组成,具有明显的大尺度空间分布特征,因此研究地下实验室地下结构群的地震反应特征已迫在眉睫。本文基于这一背景,结合现阶段近场强震非均匀地震动场模拟及地下工程抗震的研究成果,考虑岩体断层的分布,建立了岩体—地下结构体系精细化三维大尺度整体有限元计算模型,开展关键断裂对高放废物地质处置地下实验室近场地震安全性影响研究,研究成果探明了花岗岩硬岩场地的地震动传播规律及其特征,同时能够为地下实验室场址地壳稳定性评价提供技术支撑。

数值分析模型针对位于甘肃北山的高放废物处置地下实验室洞室群,综合震源—传播路劲—局部场地特性提供计算地震动输入,采用有限元—无限元耦合方法建立综合模型,开展体系响应分析。分析方法涉及到地震动模拟、地下工程高效数值建模、无限域有限化等关键技术。

采用缩减积分单元进行模拟,设置沙漏刚度为0.01,单元数约605万。综合考虑网格控制,优化不同部位网格尺度,以适应计算与分析的要求。在几何变化较大的部位(尤其是突变处),采用比较密集的网格;而在几何变化较小的部位,为减小计算规模,则采用相对稀疏的网格。数值模型网格跨尺度方法既可以保持相当的计算精度,又可使网格数量减小。其次,网格划分时考虑了截止频率所对应的波长沿深度的变化,竖向网格最大尺寸h_max取截止频率对应波长的1/8~1/10,按下式计算确定:

其中:V_s为剪切波速;f_max为截止频率。

表1~表2给出了新场岩体代表性钻孔岩心基本物理性质参数、吸水性参数、弹性波速参数及其随深度的变化规律。截止频率f_max取为20 Hz时,据现场钻孔提供的波速,竖向单元尺寸25~31.5 m。横向单元尺寸13.5~17.0 m。综合考虑计算精度及稳定性,本项目在竖井及断层位置局部加密网格,采用1~5 m不等的单元尺寸。

1973年Ungless等^[17]提出了无限元传输边界的基本思想,随后由Bettess^[18]、Beer^[19]等对其进行了发展,目前已被广泛应用于半无限域的研究^[20-21]。本文建立模型时,采用有限元网格模拟感兴趣的近场区域,使用无限元模拟远离结构的远场区域。在无限元分析理论中,通过边界上引入阻尼吸收波的辐射能量,使波的反射对近场区域的影响可以忽略。以二维无限元传输边界为例,为了抑制波的反射,设置一个阻尼边界条件,其阻尼应力为:

式中:u·x、u·y为振动速度;d_p、d_s为介质密度,与波速有关,可表述为:

式中:ρ为密度;λ和G为各向同性线弹性材料拉梅常数;c_p、c_s分别为P波和S波波速。

基于前期地质调查和工程设计资料并结合上述建模技术,建立了高放废物地质处置地下实验室精细化数值模型。由于BS32孔未来是地下实验室人员竖井,位于地下实验室场址的中心位移,因此在建立区域地质模型时,将该孔获得的岩石力学参数作为控制参数。数值模型为3 km×3 km×1.5 km的地质空间岩体,选入的区域范围见图1a。选入的断层为F31、F34、F29-1、F33、F32,如图1b所示。本文重点对F29-1是否穿越地下实验室开展讨论。根据提供的地质图确定了断层宽度及走向,模型边界采用无限元设置模拟半无限空间的能量开放特性,最终岩体计算模型如图1c所示。地下实验室结构模型见图2。通过几何清理及重构,生成了地下实验室模型。

模型共设计了18组模拟工况(表3)。针对地下实验室模型参数,分别考虑不同震级特性、是否穿越F29-1断层、岩体动力参数特性、断层参数特性4类对计算结果有较大影响的工况因素进行组合,分别对参数进行敏感性分析。弱化带杨氏模量取两侧分区杨氏模量平均值的1/3作为破碎层力学参数。

采用随机有限断层法进行设定地震的模拟,以获取合理考虑区域地震构造环境及地震动特性的地震动。前期的断裂活动性鉴定工作表明,场址西侧距离约20 km的旧井断裂为距离场址最近的能动断裂,确定其最大潜在发震震级为6.5级,因此在旧井断裂进行震级为6.5级和7.0级的地震模拟,考虑对场址最不利的震源起始破裂情况,分别给出强地震动模拟的加速度时程和反应谱,相关区域模拟参数可见文献[20-21]。对于6.5级设定地震,新场预选地段场址的地震动持时大约在7~9 s之间,水平向基岩地震动平均值为317±130 Gal;对于7.0级设定地震,地震动持时大约在13~14 s之间,水平向基岩地震动峰值加速度平均值为565±91 Gal。反应谱中高频部分较为发育,符合硬基岩场地场址的一般地震动特性。模拟地震动见图3、图4。

为便于简化表述,本文将人员提升井定为竖井一、入风井定为竖井二、排风井定为竖井三。定义加速度放大系数为某一工况下竖井各测点加速度反应最大值与输入地震动峰值加速度的比值。图5给出了6.5级地震动作用下,不同时刻地震动在体系内的地震动传播特性。由图可知,在振动初期,地震动传播与解析特性符合较好,但由于本模型内复杂断层存在及内部开孔等几何与力学特性的特异性强,地震动在到达自由地表后,体系内的加速度出现了复杂的叠加效应,局部出现不均匀分布,且在断层地表位置加速度出现明显的集聚现象。

不同工况下竖井各测点水平方向峰值加速度放大系数如图6所示。由图可知,同一工况下,相同标高各竖井加速度反应规律较为一致:随输入地震动峰值加速度增大,人员提升竖井及入风、排风竖井峰值加速度放大系数大体呈增大的趋势,但局部位置出现了峰值加速度放大系数增大的现象。这可能是由于断层的存在使得地震动在自由地表反射后出现了体系内的复杂叠加效应,使得地下实验室结构在不同埋深位置出现了上行波与下行波的叠加,出现了局部位置的峰值加速度放大现象。在6.5级地震动作用下,C1工况的竖井各测点水平向峰值加速度放大系数最大,7. 0级地震动作用下,C10工况的竖井各测点水平方向峰值加速度放大系数最大。结果表明:各竖井同一标高处水平方向峰值加速度放大系数空间差异性小。其次,-240 m 辅助试验水平高度以深,竖井x方向放大系数先增大后减小;-240 m 辅助试验水平高度以浅,竖井x方向放大系数逐渐增大,在竖井顶部取得最大值,其加速度放大系数接近1.50。

由图6可知,同一工况下,各竖井y方向加速度反应规律较为一致;随输入地震动峰值加速度增大,竖井y方向峰值加速度放大系数呈增大的趋势。模拟结果表明:沿高度方向竖井各测点y方向加速度反应规律较为相似,除人员提升竖井个别测点外幅值均较为接近,结果表明斜坡道是否穿越F29-1断层对各竖井加速度反应影响不大。7.0级地震动作用下沿结构高度方向,各竖井y方向峰值加速度放大系数呈先减小后增大的趋势,在竖井顶部取得最大值,其加速度放大系数为1.64。

另外由图可知,地震动在目标岩体至下而上传播的过程中,峰值加速度角度出现了复杂分布特性,考虑到加速度主要反应地震动高频的传播特性,可以发现目标岩体对加速度高频的影响更为突出,分析对比得出:从地震动传播的角度,硬基岩夹软弱层的地层结构将会出现更多的界面散射、反射等特性,出现复杂的地震动在不同位置的叠加及复杂的空间分布特性;从同震级不同工况对比角度而言,相同震级条件下,地震动峰值加速度时空、幅值分布基本相同,这说明目标岩体对地震动能量的敏感性较强。

7.0 级地震动作用下地表峰值加速度分布见图7,由图可知:不同工况条件下,地表峰值加速度分布存在差异,峰值点位置也有所不同,但主要集中在断层表面所在位置,对比地震各工况地表峰值加速度分布可以发现:随着工程岩体力学参数的弱化,地表峰值加速度也随着减小,这在一定程度上体现出区域断层与岩体的相互作用,岩体弱化后地震能量在体系内重分布,相反岩体材料强化对结果影响较小(如C1对比C2 、C3、C4);当断层及岩体材料同时设置成较软材料时,地表加速度放大系数达到最大(对比C1与C5);断层材料越软,其地表的峰值加速度放大越明显,相反,若断层材料参数由周边岩体材料参数的1/3增强至1/2,地表加速度放大系数并未有显著变化(对比C1与C6、C7);是否穿越F29-1断层对地表的峰值加速度分布影响不显著(对比C1与C8)。由此可知,断层所在位置地表为加速度放大较为显著的位置。从各工况的整体分布而言,断层在地震动作用下将会对地震动产生更大的放大效应。建议地表场区建筑群集中于BS32孔地表位置周边范围内。

不同震级地震动输入条件下,地下实验室峰值见表4。由表可知:6.5级地震动作用下工况C4的地下实验室峰值加速度反应最大,工况C7、工况C1 及C8次之,工况9最小。结果表明:6.5级地震动作用下岩体弹性模量越大,地下实验室峰值加速度反应越大,加速度反应对高强度岩体更为敏感,地下实验室峰值加速度反应随岩体力学参数的强化而逐渐增大;C1的地下实验室峰值结构加速度反应与C8的峰值加速度反应基本完全相同,表明斜坡道穿越F29-1断层对地下实验室峰值加速度分布近乎没有影响。

7.0级地震动作用下,工况C16的模型结构加速度反应最大,工况C10、C11 及C17 次之,工况C18最小。结果表明:7.0级地震动作用下,断层弹性模量越大,地下实验室峰值加速度反应越大,C10与C17地下实验室峰值加速度反应结果基本相同,表明斜坡道是否穿越F29-1断层对地下实验室峰值加速度反应影响较小,可以忽略。对比6.5级地震动及7.0级地震动作用下地下实验室地震动峰值加速度反应,可以发现:不同地震动作用下,地下实验室结构的加速度响应存在差异,这是由于地震动本身的特性引起的(频率组构及时程特性),其次也体现了地下实验室地震反应对不同岩体及断层的选择性较为明显。

地震动在目标岩体内传播的谱特性见图8,具体对比观测点如图9所示。图8中各节点编号分别对应图9中的4个关键截面。由图8可知:同一工况下,人员提升竖井及入风、排风井傅里叶谱反应规律较为一致,各测点在傅里叶谱各频段振动的幅值随输入地震动峰值加速度增大呈增大的趋势;此外,不同地震动作用下加速度反应频谱组成表现出不同的特点:6.5级地震动作用下,地下实验室各测点x方向的加速度在1~10 Hz范围内较为发育,通过对比竖井y方向与x方向上各测点的傅里叶谱发现,y方向上各测点频谱的低频成分更加丰富;其次,6.5级地震动作用下,y方向上竖井顶部测点频谱振幅最大值对应的基本频率比x方向上的更高。

7.0级地震动作用下,y方向上竖井顶部测点在频段采集范围内振动的振幅出现双峰值现象,且其二阶频率较x方向上顶部测点第二幅值峰值对应的频率低。z方向上,各测点振动幅值最大值比x、y方向都大,且各频段内主要频率清晰分布较均匀。在6.5级地震动作用下,各竖井顶部测点的傅里叶谱高频成分增大。在7.0级地震动作用下,竖井顶部测点仅在10~15 Hz频段内有一个振动幅值最大值。为了研究岩体力学参数对测点傅里叶谱的影响,将C1与C4、C10与C14分别进行对比发现:随着岩体力学参数的弱化,在地震加速度较小时,各测点加速度频谱的5~12 Hz频段振动幅值显著减小;位于-240 m辅助试验水平以浅时,测点的振幅最大值基本符合先增大后减小的规律,振幅最大值对应的频率随岩体弱化越来越低;而位于-240 m辅助试验水平以深的测点振幅最大值随岩体弱化逐渐减小。

分别对比C1与C6、C3与C5、C12与C14、C10与C15,结果显示:岩体力学参数越小,各测点的振幅最大值越小;且岩体波速参数对地震加速度的大小较敏感,加速度幅值越大,振幅的变化范围越大。通过比较发现:地震动作用下斜坡道是否穿越F29-1断层对竖井高度方向上各测点的加速度反应影响不显著。

6.5级和7.0级地震动作用下斜坡道峰值加速度放大系数见表5、表6,测点编号如图8所示。由表可知,随着输入地震动峰值加速度增大,相同测点位置加速度放大系数变化并不显著,这与岩体地下实验室所在的地质环境相关。整体来看,地下实验室本身相对于地震动波长(约为300 m)对地震动的响应差异并不大,不同震级作用下地下实验室的不同响应主要源于不同地震动的频率组构差异性。同一工况下,斜坡道顶部加速度反应最大,斜坡道靠近竖井一侧及底部斜坡道与构筑物相连接处加速度反应均较大,峰值加速度整体空间分布呈现在斜坡道顶部、靠近竖井一侧及底部与构筑物连接处。

本文采用有限断层法得出位置基岩地震动,根据现有资料实现基岩面以下的地下实验室精细化建模,通过地震动输入实现地震过程的情景再现,并获取地震动传播及位移响应等特性。采用时域有限元方法开展了花岗岩硬基岩场区地下实验室的大尺度、区域性有限元模拟,工作涉及到无限域有限化、目标材料动力特性、大体系精细化建模、高精度数值算法等方面。通过开展目标区域的地震动模拟及有限元计算可以得出地下实验室地震动响应,主要结论如下:

1)因研究区发震构造(旧井断裂)距离目标场址有一定距离,传统的地震动衰减关系预测方法难以考虑近场可能存在的有限断层效应、破裂方向性效应和上盘效应等近场震源效应,本文采用随机有限断层法,通过考虑发震断层的震源破裂过程和断层几何产状,可以有效地考虑上述近场震源效应。

2)模拟得到的旧井断裂发生6.5级地震时,估计目标场址的地震动持时大约在6~9 s之间,水平向谱加速度大约在132~679 gal之间;发生7.0级地震时,估计目标场址的地震动持时大约在13~14 s之间,水平向PGA大约在444~749 gal之间。研究结果说明,目标场址场地类型为花岗岩硬基岩场地时,地下洞室群的加速度反应谱中高频部分发育明显。

3)因局部软化带的存在,花岗岩硬基岩场地的地下实验室竖井各方向加速度放大系数均存在差异,体现出明显的空间分布特征,各方向放大系数在1.5~2.0区间,并与输入地震动和软化带的材料属性相关。地下实验室巷道结构不同位置加速度分布空间效应显著,软化带对应地表峰值加速度富集明显,工程中应规避此区域,但是否穿越断层对地下硐室群的加速度反应空间分布影响较小。