深圳市区域应力场有限元模拟及稳定性评价研究

张建国，马俊强，黄亚娟

（中国煤炭地质总局广东煤炭地质局勘查院，广东广州 510440）

深圳市地处广东省中部，珠江三角洲的东南部，东临大亚湾，西至珠江口，北接东莞，南隔深圳河与香港新界接壤，地理位置优越，交通发达。断裂活动性和区域（地壳）稳定性是关乎深圳规划发展的重大地质问题［1］。在珠江三角洲多个地区，均发现因强震产生的地裂缝、沙土液化、软沉积物错切和揉皱变形现象。深圳地区北东向和北西向断裂发育，至今仍有一定的活动性，根据1986—2009年对F8断层微量位移测量，表明其属现今中速活动断裂，年趋势活动速率为0.35～0.45mm［2］，1994 年F8断层附近民居曾出现地热异常，跨断层的向西村一栋六层楼房开裂［3］。深圳市自1987 年以来进行了大规模建设，由此必然改变原有地壳结构的边界平衡条件、载荷平衡条件及其平衡状态，使原本具有活动性的断裂带的弱平衡状态地质环境受到破坏［3］。目前，沿该断裂带分布的破碎岩分布区、橐状与槽状强烈风化地段、软弱地层、饱和粉细砂层等不良地基岩土分布区，已发生部分建筑物变形和产生裂缝。因此进一步加强对断裂构造活动性研究和区域稳定性评价对深圳市城市规划发展、预防和减轻地质灾害的损失、国土资源保护等方面都具有重要的现实意义。

1 断裂基本情况

以横岗-黄贝岭断裂（F1-3）、简龙村断裂（F1-3-3）、清风岭断裂（F1-3-4）、炮台山断裂（F1-3-5）、狮子石断裂（F1-3-6）、田螺坑（西坑）断裂（F1-3-7）、观澜-丹竹头断裂（F2-6-1）和大芬断裂（F2-5-1）这8条断裂分布区作为主要研究区域［4］（图1）。

图1 断裂分布位置Figure 1 Fracture distribution

2 三维有限元数值分析模型

ABAQUS 有限元数值计算软件最早由美国ABAQUS公司研制开发，经过多年发展，现已成为一款具有丰富前后处理功能与求解器模块的通用有限元软件。ABAQUS 软件体系包括前后处理模块、通用分析模块、专用分析模块与第三方软件接口模块［5］。ABAQUS 内置多种材料本构模型、单元类型、边界条件与分析手段，可以计算各类固体力学、结构力学乃至流体动力学系统，实现复杂边界条件与非线性问题处理的能力极为突出，对岩土工程领域问题具有良好适用性［6］。

2.1 模型建立

根据前人资料，深圳市最大主压应力方向以NW-NNW 为主，区内最大主应力值在30m深度以上为3.94～1.78MPa，最小主压应力值为3.33～1.18MPa，最大剪应力值为0.94MPa，各应力值基本接近平面的等压状态［7-8］。通过考虑深圳地区构造应力场特征，根据有限元数值模拟的需要，建立地质计算模型。有限元模拟区域的选取：区域范围沿经度方向取4.52×104m，沿纬度方向取3×104m，平面面积为1 356km2，竖直向下取2 500m，主要断层见表1。

表1 深圳断裂带主要断层统计Table 1 Statistics of main faults in Shenzhen Fault Zone

在模拟区域内忽略较小断裂的影响，即只考虑上述断裂及他们的切割关系，并通过简化及约定，可建立计算区域的三维地质模型（图2）。

图2 深圳主干断裂断层模型Figure 2 Shenzhen Trunk Fault model

2.2 设置边界条件和荷载

（1）模型的边界条件定义

底边界x、y、z方向全部位移约束。

南北方向边界条件为y方向位移约束。

东西方向边界条件为x方向位移约束，地面无任何约束。

边界力作用方式、方向及大小与区域构造应力场一致。

（2）荷载

区域稳定性一般受构造稳定性、岩土体以及地面稳定性等控制，其中构造是主导因素，设定模型不受外力作用，受力只考虑有限元间构造力的传导。其中重力荷载，重力加速度设置为沿z方向-9.81m/s²。

（3）断层与地层接触类型

法向行为设置为硬接触，切向行为设置为penalty 接触（点面接触），依据断层角砾岩内摩擦角13°，切向接触摩擦系数取tan13°=0.230 868 181 1。因断层角砾岩黏聚力较小，模型z 方向厚度较大，除地面下浅层土外，黏聚力对抗滑力影响较小，可忽略。

2.3 有限元网格划分

对地质模型进行有限元计算网格剖分，岩石和断裂均采用4 节点线性四面体实体单元，整个计算模型共剖分44 140 个节点，195 484 个单元，有限元网格剖分结果见图3。

图3 三维有限元剖分网格Figure 3 Three-dimensional finite element mesh

2.4 计算参数设置

本次有限元模拟的岩土介质类型可简化为5种，根据试验成果，参考深圳地区经验参数，经反算对比等综合确定其物理力学参数，见表2。

表2 深圳断裂带三维地应力场有限元计算参数Table 2 Three-dimensional finite element calculation parameters of in-situ stress field in Shenzhen fault zone

2.5 计算结果分析

本文通过建立的断层-基岩接触单元模型来模拟地质体中断层的存在和其一定摩擦机制下的滑动特征，其实质是将断层与基岩之间的力学作用以接触摩擦代替，断层的活动则通过不同地质体间的接触滑动来模拟［9］。基于强度折减有限元法进行断裂稳定性计算，分析断裂带发生错动前的应力场分布及其稳定性系数。三维弹塑性有限元数值计算中以拉应力、拉应变为正，压应力、压应变为负，与岩土工程定义相反［10］。但为了和岩土工程领域一致，下文讨论中，所提及最大、最小主应力按岩土工程领域的确定方式，但数值的正负没有改变，即按有限元计算结果进行分析。

2.6 断裂稳定性结果分析

断裂构造对地应力场的影响是十分复杂的［11］。由于测量技术和资料的不足，加之岩体中的应力分布又受多种因素的制约，目前尚难总结出定量的规律性特征。因此，通过系统的数值模拟，可以研究断裂构造对地应力状态的影响特征，分析断裂构造附近应力状态的变化规律，进而较为准确地评估工程场地的地应力状态，具有重要的实用价值［12］。

（1）断层附近的应力方位

从模拟结果可以看出，断层附近的应力方位均发生不同程度的变化，在断层附近应力方位与边界应力方位呈一定夹角，而远离断层，应力方位逐渐趋于与边界应力方位一致。断层附近的主应力方位与断层可以斜交，也可以平行或垂直，而且断层两盘的应力方位可以不同，即断层附近的应力方位的变化是十分复杂的。分析表明，上述断层附近应力方位变化的幅度和发生变化的范围因断层的力学性质、围岩的物理力学性质及边界条件等的不同而不同。断层带两侧应力大小不连续，说明断层对应力的传递有较明显的阻隔作用［13］。

（2）断层附近的应力量值

模拟结果表明，断层附近的应力量值有明显的变化，而且这种变化也是十分复杂的。最大主应力和最大剪应力在断裂附近表现为应力集中，断层两侧的主应力量值有差别。应力集中的幅度和发生应力集中的范围均受断层的力学性质、围岩的物理力学性质及边界条件等影响［14］。

（3）断层附近的最大剪应力

最大剪应力也受到断层的明显影响，断层带两侧最大剪应力大小往往不连续，多在断层交汇处形成剪应力集中。

为了便于观察横岗-黄贝岭断裂北东段与南西段、简龙村断裂、清风岭断裂、炮台山断裂5 组断裂在进行强度折减后的位移变化情况，在结果模型中截取横截面A-A′，B-B′（图4），分别展示出它们带在不同折减系数下的位移云图，如图5、图6所示，随着强度折减程度的增大，五组断裂带均出现了一定程度的位移，其中最大位移出现在简龙村断裂带区域，当折减系数增至6.71 时，简龙村断裂带与邻近基岩之间出现了明显的错动，其余断裂带并未出现明显位移，表明简龙村断裂带的稳定性应弱于其余四组断裂带。选取简龙村断裂带中位移最大的节点作为监测点，得到该节点处的位移与折减系数关系曲线如图5f所示。

图4 模型界面选取位置Figure 4 Model interface selection location

图6 B-B′截面处位移云图Figure 6 Displacement cloud on section B-B′

图5f表明，当强度折减系数增至5.52时，节点位移随折减系数的进一步增大而呈现出急剧提升的趋势，此后断裂带出现了显著的失稳滑动，表明在该折减系数下断裂带达到临界稳定状态，因此包含横岗-黄贝岭断裂北东段与南西段、简龙村断裂、清风岭断裂、炮台山断裂的断裂组最小安全系数为5.52。

同样得出狮子石断裂、田螺坑（西坑）断裂最小安全系数为10.55（图6f）；观澜-丹竹头断裂、大芬断裂最小安全系数为4.04。

3 区域稳定性评价

根据地震等烈度线、震源机制解、前期地应力测量测量及地形变特征等均说明深圳市现今以水平地应力作用为主，最大主压应力的优势方向为NW-SE 向（图7）。构造现今应力场有限单元数值计算结果与实际资料基本一致。

广东省区域地壳运动发展趋势总体上是由强到弱，到貌近时期以来以间歇性上升为主。区域性主要断裂带以NE向为主，等间距展布。断裂带间夹持有相对稳定地块。深圳位于现今活动量相对较小的五华-深圳断裂带西南段及其旁侧的相对稳定地块上。从现今地形变资料得知，深圳位于由东南相对上升西北相对下降的斜坡带。从区域地震活动规律看，其主要发震构造为NE 向断裂带，深圳所处的莲花山断裂带在我国东南沿海仅为中等水平。因此，从广东省区域地壳活动看，深圳市处于相对稳定—相对较稳定的部位。从区内几个中、强震发展地区计算的最大剪切应力与深圳地区的实际情况获知，深圳地区未来发生中、强地震的可能性较小。

稳定性评价因子主要选取地震活动性、断层活动特征、深部构造特征、构造应力场特征、地壳形变特征和岩体性质6 种，根据有限元模拟的结果和区域地质条件，本文将地壳稳定性划分为稳定和较稳定等两个等级（图8）。

图8 深圳市区域稳定性分布Figure 8 Regional stability distribution of Shenzhen city

3.1 稳定区地质特征

1）断裂构造特征：断裂构造较发育—不发育，断裂规模较小，构造岩胶结较好。

2）现今构造活动性与地质灾害特征：现今构造活动性微弱，岩体完整性好，各类地质灾害均不发育。

3）岩体质量特征：区内岩体质量较好—好，岩石强度高。

4）水文地质环境及特殊岩土体特征：低山丘陵地貌，沟谷不发育，水文地质环境中等及以上，基岩裂隙水丰富，风化层厚度不大。

3.2 较稳定区地质特征

1）断裂构造特征：断裂构造发育—较发育，断裂规模较大，断层破碎带宽，构造岩胶结差。

2）现今构造活动性与地质灾害特征：现今构造活动性微弱或有一定活动性，崩塌及滑坡等坡地类地质灾害较发育，局部岩溶塌陷较发育。

3）岩体质量特征：区内岩体质量较差，灰岩发育，沿断裂带裂隙发育，岩体质量较差。

4）水文地质环境及特殊岩土体特征：水文地质环境复杂，基岩裂隙水较丰富，风化层厚度不均匀，软土发育。

4 结论

在深圳市城市地质调查工作中开展深圳市断裂活动性和区域稳定性评价研究，既是贯彻落实国务院和广东省的决策部署，也是夯实深圳市基础地质、水文地质、工程地质、地质灾害现状工作的基础，编制国土空间规划、健全完善城市地质安全监测预警体系的重要依据。构造稳定性是确定区域稳定性的主导因素，只有在构造稳定性相同的情况下，进一步评价地面稳定性，外动力地质作用才会有条件地上升为重要因素。

1）有限元模拟情况与区域构造运动演化发展特征、构造现今活动特征、现今地应力状态特征等基本一致，可信程度较高。

2）断裂构造对地应力场的影响是十分复杂的。由于测量技术和资料的不足，加之岩体中的应力分布又受多种因素的制约，目前尚难总结出定量的规律性特征。通过系统的数值模拟，可以研究断裂构造对地应力状态的影响特征，分析断裂构造附近应力状态的变化规律，进而较为准确地评估工程场地的地应力状态，具有重要的实用价值。

3）基于深圳市区域应力场有限元模拟，将深圳市陆域的区域稳定性分区分为稳定区和较稳定区，其中稳定区面积1 107km2，占55.43%，较稳定区面积890km2，占44.57%。