一种近断层的隧洞应力场和地震动力响应的数值分析方法、终端及介质

本发明涉及地震分析，尤其涉及一种近断层的隧洞应力场和地震动力响应的数值分析方法、终端及介质。

背景技术：

1、断层是隧洞开挖和使用过程中较常见的不良地质构造，近断层隧洞围岩的力学行为与稳定性受构造应力场影响显著。隧洞工程的建设首先应选在较稳定的地层中，如遇断层应首选避让且避让距离要符合设计规范的要求，或使隧洞与断层的走向相互垂直。但在实际隧洞工程建设中，因地质条件、设计方案、施工方法等因素制约，部分隧洞工程不可避免地建造在距离断层较近的地方，且与断层走向斜交或平行，对围岩稳定性的控制造成较大挑战。

2、目前研究此类问题的方法有通过制备相似材料，如非连续性的断层，对近断层隧洞工程进行模拟试验。模型试验的理论依据是相似原理构建满足相似条件的模型，亦即要求模型与实体(原型)相似，模型能够反映实体的情况。这时，模型与原型，除了几何形状相似以外，同类的物理量(如应力、应变、位移)也必须按比例相似。这些比例必须满足各种力学条件。但由于材料的选择和模型的制作，使得模型试验的周期较长、工作量大，且相似条件难以完全满足，难以根据设计反应谱对相似模型施加地震动力边界条件。

3、因此，有必要采用断层接触面的模型进行模拟计算，用计算机数值平台代替真实试验室中的相似模型，缩短模型试验的准备周期，节省模型试验的工作量。

技术实现思路

1、本发明的目的在于克服现有技术中的不足，提供一种近断层的隧洞应力场和地震动力响应的数值分析方法，建立地质力学模型，模拟地质力学模型的动力边界条件，利用数值模拟平台的高效性以及对不同工程的适应性，过程方法易操作和实现，有利于提高隧洞结构抗震设计的效率。

2、为实现上述目的，本发明提供了如下技术方案：

3、第一方面，本发明提供一种近断层的隧洞应力场和地震动力响应的数值分析方法，其特征在于，包括以下步骤：

4、s1：基于flac3d6.0平台，通过预构建不含结构面的完整岩石试样的数值模型和含结构面的岩石试样的数值模型，所述含结构面的岩石试样的数值模型包括结构面，通过对比所述不含结构面的完整岩石试样的数值模型和含结构面的岩石试样的数值模型得到结构面的变形参数，根据所述结构面的变形参数，施加速度边界条件，计算结构面的法向刚度和剪切刚度参数，由flac3d软件模型库或者有限差分法得出断层接触面模型，对所述断层接触面模型标定所述结构面的法向刚度和结构面的剪切刚度；

5、s2：根据所述结构面的法向刚度和结构面的剪切刚度，构建地质力学模型，包括若干单元，每个所述单元包括若干个节点和若干个面单元；假设所述地质力学模型中的围岩材料服从线弹性强度准则，断层接触面模型服从莫尔-库伦强度准则；对所述地质力学模型赋予所述结构面的变形参数，根据所述地质力学模型的边界坐标施加法向约束条件和地应力边界条件，对所述地质力学模型进行模拟计算，并将所述地质力学模型的节点的位移和速度设置为零，保留初始地应力场的计算结果；

6、s3：以步骤s2初始地应力场的计算结果为基础，构建近断层隧洞模型并进行开挖隧洞；模拟隧洞开挖后所述近断层隧洞模型的次生应力场并计算应力集中系数；

7、s4：基于matlab平台，按照《公路隧道抗震设计规范(jtg 2232-2019)》和《地下结构抗震设计标准(gb/t 51336-2018)》中关于隧洞抗震设计的相关要求，根据地震动加速度反应谱合成设计人工地震动加速度时程，并进行滤波和基线校正，得到满足设计规范要求的人工地震波；

8、s5：以步骤s3次生应力场应力集中系数的计算结果为基础，对地质力学模型施加黏弹性人工边界和阻尼系数，根据震源与隧洞工程的相对位置对地质力学模型边界施加所述人工地震波作为地震动力边界条件，在反应过程中记录单元的峰值加速度和围岩峰值应力集中系数，分析地质力学模型动力响应，得到近断层隧洞的峰值应力集中系数分布应力场图，以及各监测点的水平加速度反应时程。

9、进一步地，步骤s1中岩石试样数值模型的构建包括：

10、s1.1：采用软件生成不含结构面的完整岩石试样的数值模型和含结构面的岩石试样的数值模型，按照《工程岩体试验方法标准(gb/t 50266-2013)》推荐的单轴抗压强度试验标准，将所述不含结构面的完整岩石试样的数值模型和所述含结构面的岩石试样的数值模型的高度与直径之比设置为2：0；

11、s1.2：假设所述岩石试样材料服从莫尔-库伦强度准则，对比不含结构面的完整岩石试样的数值模型和含结构面的岩石试样的数值模型得到结构面的变形参数，不含结构面的完整岩石试样的数值模型的法向变形始终小于含结构面的岩石试样的数值模型，不含结构面的完整岩石试样的数值模型的中的应力-变形曲线的斜率大于含结构面的岩石试样的数值模型，根据含结构面的岩石试样的数值模型的边界坐标向所述含结构面的岩石试样的数值模型的边界施加法向约束条件；

12、s1.3：模拟单轴压缩试验的加载速率，分别对不含结构面的完整岩石试样和含结构面的岩石试样施加速度边界条件；试样以坐标轴z为旋转轴，圆平面位于xoy坐标平面；生成上部边界和下部边界，向所述上部边界和下部边界施加速率为10-5m/s的速度边界条件；采用fish语言实现单轴压缩试验，包括以下步骤：

13、s1.3.1：生成全局映射变量surface；

14、s1.3.2：所述不含结构面的完整岩石试样的数值模型和含结构面的岩石试样的数值模型包括节点体，遍历所述节点体，如果某节点体gp位于所述上部边界，则将该节点体保存至surface；

15、s1.3.3：向全局映射变量surface中施加(1～10)×105m/s的速度边界条件，模拟单轴压缩试验的加载速率；

16、s1.3.4：遍历所述节点体，用最大不平衡力的函数gp.force.unbal.z(gp)记录全局映射变量surface上z方向不平衡力，保存至变量zforce；

17、s1.3.5：计算法向应力axiastress＝zdoece/π*r2；

18、式中axialstress表示法向应力；zforce表示z方向不平衡力；r表示不含结构面的完整岩石试样的数值模型以及含结构面的岩石试样的数值模型的半径；π表示圆周率；

19、s1.3.6：设置人工收敛条件，当当所述不含结构面的完整岩石试样的数值模型和含结构面的岩石试样的数值模型的位移达到终止应变时，模型计算收敛；

20、设完整岩石试样的变形为δvr，含结构面的岩石试样的变形为δvt，则结构面的法向闭合变形δvj为：

21、δvj＝δvt-δvr      (1)

22、s1.4：用如下方程拟合σn-δvj应力-变形关系曲线：

23、

24、式中vm表示结构面的最大法向闭合变形，为σn-δvj应力-变形关系曲线的渐近线；σn表示法向应力；表示由该方程拟合的结构面法向刚度；

25、s1.5：根据式(2)拟合的方程，求解渐近线方程在原点处的切线斜率kni，即结构面的初始法向刚度：

26、

27、s1.6：根据含结构面的岩石试样参数，标定在所述法向应力σn作用下，断层接触面模型的结构面的法向刚度和结构面的剪切刚度：

28、

29、

30、式中l为断层的倾向长度；jrc为断层结构面的粗糙度系数；jcs为结构面的壁岩强度；φr为断层结构面的残余摩擦角；kn为结构面的法向刚度；ks为结构面的剪切刚度。

31、进一步地，步骤s2中隧洞模型和断层接触面模型的构建包括：

32、s2.1：隧洞模型的构建，对于形状较为简单的隧洞断面，如圆形、半圆拱形隧洞，采用结构化网格建立模型；对于形状相对复杂的隧洞断面，如直墙三心拱形、马蹄形隧洞，采用非结构化网格建立模型；

33、按照真实尺寸建立隧洞断面的xoz平面模型，使模型的几何中心位于坐标原点；设置模型的厚度为1.0；

34、s2.2：围岩模型的构建，在已建立的隧洞模型周边设置范围为4～6倍隧洞半径的围岩，分别对围岩的上盘和下盘进行分组命名；

35、s2.3：断层接触面模型的构建，基于有限差分数值平台，根据断层的产状和性质建立接触面模型；对于结构化网格，在断层的上下盘之间创建面单元生成接触面模型；对于非结构化网格，在上下盘分组之间直接生成接触面模型；对接触面模型赋予kn、ks等断层的力学参数。

36、s2.4:根据所述地质力学模型的边界坐标施加法向约束条件和地应力边界条件，将模型边界面单元的法向速度设置为0，对模型施加的水平应力和竖直应力。

37、进一步地，步骤s3中采用fish语言实现隧洞开挖后围岩应力场的模拟：

38、s3.1：遍历模型单元，提取每个单元体zp的应力状态，使用fish语言实现：

39、

40、式中σx为数值计算模型中单元体的水平应力；σz为数值计算模型中单元体的垂直应力；τxz为数值计算模型中单元体的剪应力；zp为每个单元体的指针变量；

41、s3.2：计算每个单元体中心到坐标原点/隧洞几何中心点的距离r：

42、

43、s3.3：不失一般性，以直墙半圆拱为例，根据隧洞断面的形状，计算每个单元体中心到坐标原点连线的正弦和余弦，使用fish语言实现：

44、s3.3.1：单元体位于圆弧拱部位，即zone.pos.z(zp)>0.0，则

45、

46、s3.3.2：单元体位于直墙部位，即-1.2<zone.pos.z(zp)<0.0，则

47、

48、s3.3.3：单元体位于墙脚部位，即zone.pos.z(zp)<-r，且math.abs(zone.pos.x(zp))>r，则

49、

50、s3.3.4：单元体位于底板部位，即zone.pos.z(zp)<-r，且math.abs(zone.pos.x(zp))<r，则

51、

52、s3.4：根据转轴公式计算隧洞周边围岩的径向应力σr和环向应力σθ：

53、

54、s3.5：定义应力集中系数k为环向应力σθ和垂直地应力pv之比，求解围岩单元体zp的k：

55、k＝σθ/pv     (13)

56、将计算结果保存至变量zone.extra(zp)中。

57、进一步地，步骤s4采用三角级数法合成人工地震动时程，使用matlab平台实现：

58、s4.1：创建合成人工地震动反应时程加速度的主程序，设置变量用以保存反应持续时间t、反应间隔dt，频谱上限ωu、频谱下限ωl、频域采样间隔δω；生成空数组变量acc用以保存反应加速度；

59、s4.2：根据随机振动理论，输入平稳随机过程的单质点体系的最大反应分布：

60、

61、式中ag(t)为平稳随机过程，初相位为(0,2π)内均匀分布的随机相角，ωk为第k个频率的圆频率，ak(ωk)为第k个频率的幅值，由下式计算：

62、

63、式中s(ωk)为功率谱密度函数，δω为频率采样间隔，由下式计算：

64、δω＝2π/t   (16)

65、式中t为反应持续时间；

66、s4.3：创建强度包络函数子程序，在合成人工地震动反应时程加速度的主程序中进行调用，向子程序传入时间变量t，返回强度包络函数值f(t)；根据等效震级m和等效震中距d，计算强度包络函数中的加速度平稳段开始时间t1、加速度平稳段结束时间t2和衰减系数c：

67、

68、构建强度包络函数如下：

69、

70、s4.4：创建功率密度函数子程序，在合成人工地震动反应时程加速度的主程序中进行调用；向子程序传入频率变量ωk，返回功率谱密度函数值s(ωk)；设置变量用以保存结构自振周期t、场地特征周期tg、阻尼比ζ、加速度反应谱最大值smax、下降段衰减指数γ、反应超越概率p；用于合成设计加速度时程的反应谱sa为：

71、

72、功率谱密度函数由下式计算：

73、

74、s4.5：运行合成人工地震动反应时程加速度的主程序，调用上述子程序，合成人工地震动加速度时程，保存至数组变量acc：

75、

76、s4.6：创建低通滤波主程序，生成数组变量input_acc保存步骤s5.5中的人工地震动加速度时程acc；生成变量保存滤波器阶数n和截止频率ωn；

77、使用matlab平台的butterworth低通滤波器功能，生成滤波器的传递函数系数a和b，返回长度为(n+1)的行矢量：

78、

79、通过传递函数系数进行滤波，返回滤波后的加速度反应时程filered_acc；

80、s4.7：创建基线校正主程序，生成数组变量iput_acc保存步骤s5.6中滤波后的加速度反应时程filtered_acc；使用采样滑动窗口算法计算运行平均加速度：设置滑动窗口大小k，返回本地k点均值的数组，其中每个均值都是在iput_acc的相邻元素上通过长度为k的滑动窗口计算得出的；当k为偶数时，窗口以当前和前一个元素为中心；当没有足够的元素填充窗口时，窗口大小会在端点处自动截断；当窗口被截断时，只对填满窗口的元素取平均值；返回运行平均加速度smooth ed_a，则基线校正后的加速度反应时程为：

81、modified_a＝iput_acc-smooth ed_a     (23)

82、式中modified_a为加速度反应时程；input_acc为数组变量；smooth ed_a为返回运行平均加速度；

83、s4.8：创建时程转换子程序，可任意主程序中进行调用；将加速度反应时程转换为速度反应时程和位移反应时程；向子程序传入加速度反应时程保存为数组变量iput_acc，初始化速度变量vel和位移变量disp：

84、

85、采用数值积分计算速度反应时程和位移反应时程作为返回值：

86、

87、s4.9：创建傅里叶变换子程序，可在任意主程序中进行调用；采用矩阵计算数值平台的快速傅里叶变换功能，将时域地震动信号转换为频域信号进行分析；向子程序传入加速度反应时程保存为数组变量iput_acc，返回加速度幅值频谱。

88、第二方面，本发明提供了一种电子终端，包括处理器及存储介质；所述存储介质用于存储指令；所述处理器用于根据所述指令进行操作以执行根据权利要求1～6任一项所述方法的步骤。

89、第三方面，本发明提供了一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现权利要求1～6任一项所述方法的步骤。

90、与现有技术相比，本发明所达到的有益效果：

91、(1)本发明基于有限差分数值模拟平台，计算结构面的变形参数，对断层接触面模型参数进行标定并建立地质力学模型，以地质力学模型的计算结果为基础，构建近断层隧洞模型模拟隧洞开挖后近断层隧洞模型的次生应力场并计算应力集中系数，可根据实际工程中断层的产状和性质，充分考虑由断层形成的构造应力边界条件，有利于对此类隧洞围岩的稳定性控制。

92、(2)本发明可根据规范反应谱的参数，合成设计人工地震动加速度时程，并通过滤波、基线调整和傅里叶变换等信号处理方式，模拟地质力学模型的动力边界条件，生成满足设计要求的人工地震动加速度反应时程，有利于提高隧洞结构抗震设计的效率。

93、(3)本发明充分发挥了数值模拟平台的高效性以及对不同工程的适应性，根据设计反应谱生成的人工地震波作为动力边界条件，计算地震作用下近断层隧洞围岩应力场和动力响应，过程方法易操作和实现，计算结果具体可视，易于理解和判断。