基于FLAC3D的岩体初始高水平构造应力场反演方法

郭延辉，侯克鹏

(1.云南农业大学建筑工程学院，昆明 650201；2.昆明理工大学国土资源工程学院，昆明 650093)

地应力是地质构造和自重共同作用在地质体内形成的原始应力，初始地应力的确定是工程分析与设计的基础，对于工程的稳定与安全至关重要[1]。矿区岩体地应力测量是提供地应力场的一种最为直接的方法，然而面对广域的工程厂址，只能开展有限测点的地应力测量，而这些数据只能反映测量点处局部地应力场的情况，且有时会由于某些原因致使测量数据存在一定误差[2-3]。为了能够对工程设计和施工提供更可靠的地应力数据，根据现场的地质构造特征，结合有限的地应力实测数据，对初始应力场进行反演计算和分析，获得更为准确、适用范围较大的地应力场十分必要[4]。

近年来，研究人员提出了多种地应力场反演方法。李飞等[5]提出批数据处理(GMDH)神经网络算法，结合地应力场分布随埋深的非线性特征及局部地质构造处地应力场的非连续性特点，构建出形成复杂地质体地应力场的边界条件模式，并基于现场的稀少样本测点数据对杏山铁矿地应力场进行反演和重构；汪吉林等[6]提出以构造指标优化人工神经网络(ANN)模型的建模方案，并开展了基于构造控制的地应力人工神经网络反演方法研究；易达等[7]研究了岩体初始应力场的遗传算法与有限元联合反演法；戚蓝等[8]用系统工程灰色理论进行了地应力分析，在处理实测资料少的小样本问题取得了较好的效果；李欣等[9]采用改进的粒子群优化算法，对地应力场进行了反演计算，结果表明改进后的算法计算精度和效率更高；上述这些非线性反演方法虽然能够反演分析复杂区域地应力场，却存在收敛速度较慢、参数选择过于依赖经验等问题，并可能导致同样的研究对象而反演结果却并不唯一的现象。数值模拟方法以其三维化、可视化、成本低和能够计算复杂问题等优点广泛应用于深部地应力场的模拟计算[10-12]。然而，目前关于地应力场的反演研究大多是以自重应力场为主，针对初始高水平构造应力场的反演方法方面的研究仍然不足[13-14]。为此，以狮子山矿区为研究背景，在模拟边界讨论的基础上，采用FLAC3D反演分析各种边界条件下计算域内水平构造应力与垂直应力的分布特征，最终采用基于初始应变能理论的方法，对狮子山矿区初始高水平构造地应力场进行反演及验证，期望能对其他类似工程有一定的借鉴作用。

1 关于地应力反演方法的讨论

通常数值计算过程中赋地应力的方法主要有：①模型内部部分或某些点的地应力实测值已知，采用位移边界条件，以一定的特征和规律将已知点应力差值分布到模型中，最后对应力场反演分析并生成最终应力场；②在计算域内根据某种回归方程直接计算每个单元体的初始应力，并转化为计算单元的节点力，然后根据位移和应力边界条件，平衡计算模型系统内的应力场；③在模型位移边界条件和应力边界条件的基础上生成初始应力场，然后根据几个点实测的地应力对计算域的边界条件进行调整，并最终拟合到新的初始应力场，使之与实测点的地应力值较为吻合；④基于初始应变能理论的初始地应力场的反演方法，其基本思路为：认为将有构造应力存在的地质体作为部分具有弹性体属性的意义赋上去考虑是合理的，因此在数值模拟反演初始地应力场时，可以将构造应力作为内力施加在模型内部，即对模型内所有单元体采用三角形应力分布初值设置，而并非在边界上施加，以此来模拟构造应力场[15-16]。

2 计算模型构建及方案

在应力场反演方法讨论分析的基础上，拟通过对比分析位移边界条件、应力边界条件、混合边界条件等情况下，计算模型内垂直应力与水平应力的分布情况，从而探讨模拟高水平构造应力场的最佳方法。在实际数值模拟过程中，为了尽可能消除模型尺寸过小所产生的边界效应，通常模型左右边界向外扩展一定距离，计算作为平面应变问题考虑。通过反复试算，当计算模型长3 000 m，高1 400 m，宽50 m时，在这一尺寸下，边界效应对模型的影响可以忽略。由于地表地形对模拟结果影响较小，因此，在模型的概化过程中，将地表概化成水平地表。计算采用Mohr-Coulomb本构模型。

为了充分研究边界条件对初始地应力场分布的影响，计算模拟采用以下边界条件[15-16]：①位移边界条件，如图1所示，即在数值模拟计算中，模型左右边界的水平方向的位移固定，模型底部的位移固定，模型顶部设置为自由表面，其中G表示自重体积力；②应力边界条件，图2所示，模型左右施加应力梯形；③混合边界条件，模型底部固定竖向位移，左右边界分别为混合的三角形分布的水平分布面力和水平向位移约束；④基于初始应变能理论的初始地应力场的反演方法，其主要为位移边界条件的基础上，将构造应力作为内力施加在模型内部，即对模型内所有单元体采用三角形应力分布初值设置。计算所采用的力学参数为青灰色白云岩的岩体力学参数，如表1所示。

图1 位移边界条件示意图

图2 应力边界条件示意图

3 应力场计算反演结果对比分析

3.1 基于位移边界条件的初始地应力场的反演

在位移边界条件下，反演得到垂直应力和水平应力分布分别如图3所示。计算结果表明：应力等值线水平展布，走向稳定，等值线间间隔几乎相等，从地表至模型底边呈线性分布。但在计算模型内，同一位置的自重应力明显大于水平应力，约为3倍左右。计算的最大主应力为竖直方向，垂直应力为最大主应力，表明此种方法仅能用于计算分析以自重应力场为主导的模型，而不能模拟高水平构造应力场。

图3 位移边界条件下的应力分布

3.2 基于应力边界条件的初始地应力场的反演

一般情况下，水平地应力场是在自重体积力存在的条件下由构造运动产生的一种应力场，尤其是在高构造应力区，岩体的原岩应力状态特征总是与显著的构造应力状态特征相联系。大量的实测结果表明，水平构造应力场一般最大主压应力σ1与最小主压应力σ3方向均接近水平方向，而中间主应力σ2一般为竖直方向，在数值上等于γh(γ为重度，h为埋深)。在数值模拟反演中，自重应力通常作为一种体积力进行考虑并赋值，然而对于构造应力条件的处理办法，特别是高水平构造应力场的反演方法存在一定的分歧。初始应力场模拟是指复原其采矿工程活动前(如开挖)的应力状态，如果将构造应力作为边界应力加载在模型边界上，正如2节讨论(图2)，计算时以三角形应力施加在模型的左右边界，并且保持左右边界三角形应力不变，通过边界单元与模型内部单元接触面的节点力向内部传递，并且使得模型外表面上的分布力与其内力平衡。

上述应力边界条件下，模型的垂直应力分布如图4(a)所示。由图4(a)可以看出，垂直应力等值线图在模型中部基本上呈水平走向，且线性增加，但模型两端自重应力场明显大于理论值，且应力等值线分布并不理想。图4(b)为应力边界条件下水平应力分布图，在模型中部，同一位置水平应力高于自重应力，但模型上部和下部均未呈水平走向分布，并且模型中间偏下部分的最大主压应力远低于理论值，而模型中间偏上部分的最大主压应力又远高于理论值。可见，左右侧边界上采用应力边界条件不能够较好地计算模拟水平高构造应力场。

图4 应力边界条件下的应力分布

3.3 基于混合边界条件的初始地应力场的反演

图5(a)、图5(b)分别为分布应力施加在模型左侧时混合边界条件下垂直和水平应力分布，图6(a)、图6(b)分别为分布应力施加在模型右侧时混合边界条件下垂直和水平应力分布。从图5(a)、图6(a)可以看出，远离模型施加三角形分布面力的端部，自重应力等值线呈线性增加，一旦靠近该端部，自重应力等值线便发生变异，尤其是在施加分布面力一侧，竖直应力等值线起伏较大，自重应力远大于理论值。图5(b)、图6(b)中模型施加分布面力一侧的较小范围内，模拟结果与预期结果较好，一旦离开该端部，水平应力等值线迅速变异。具体表现为在水平位移约束侧浅部等值线过于密集，但在深部等值线又过于稀疏。矿体附近的水平应力低于预期结果，水平应力等值线越往深部越偏斜，等值线间隔在不同位置变化大，表明模拟得到的水平应力梯度变化不均匀。综上分析可知，混合边界条件既不适于模拟自重应力场，也不适于模拟水平构造应力场。

图5 混合边界条件下应力场分布(左侧应力边界、右侧位移边界)

图6 混合边界条件下应力场分布左位移边界、右应力边界

3.4 基于初始应变能理论的初始地应力场的反演分析

在地质力学中，构造应力场通常是指地壳内一定范围内某一瞬时的应力状态。构造应力场中的应力分布和变化是连续并且有规律的。根据传统的弹性理论和弹塑性理论，在自重体积力和构造运动的共同作用下，产生了位移并且做了相应功，在这一过程中，除塑性变形功不能转化为应变能外，地质体主要以弹性变形的方式储存应变能。近年来，诸多学者从能量的角度和观点出发，对岩石变形破坏过程中发生的能量积聚，耗散与释放进行了分析[15]，使得能量观点在岩体变形失稳和破坏分析中的应用愈加广泛。

假若将岩体视为弹性体，那么岩体中的某一点的微小单元的单位体积具有的应变能可表示为[16-17]

(1)

式(1)可改写为

Gεijεij=WV+WF

(2)

式中:E为弹性模量；ν为泊松比；G为剪切模量；W为应变能；WV为体积变形应变能；WF为形状改变应变能；θ为应变，θ=ε11+ε22+ε33。

如式(1)、式(2)所示，工程开挖后，初始应力场发生变化，二次应力场重新分布形成过程中，受到开挖影响的地质体中某单元体将会同时出现体积改变和形状改变。即微小单元体体积改变弹性比能WV和形状改变弹性比能WF一般都会随之变化。也就是说，岩体发生的变形破坏，其本质是能量转化的过程[13-15]，在工程开挖过程中，受影响的地质体中所储存的弹性应变能必然会随之变化。

通过前面的讨论，将有构造应力存在的地质体作为部分具有弹性体属性的意义上去考虑，在理论上是比较合理的。因此，在数值法反演初始地应力场时，可以将构造应力作为内力施加在模型内部，即对模型内所有单元体采用梯形应力分布初值设置，以此来模拟高水平构造应力场，并采用位移边界条件的情况下进行初始高水平构造应力场的反演计算。

图7(a)为初始地应力场中垂直应力等值线分布。从图7(a)可以看出，垂直应力从地表向下呈线性增加，应力迹线水平展布，走向稳定，自重应力计算值与理论值较为吻合。图7(b)为初始地应力场中水平压应力等值线分布，水平应力从地表向下呈线性增加，应力迹线较平直，水平展布，等值线间隔也很均匀，走向稳定，水平应力模拟结果与预期吻合较好。

图7 初始应地应力场中的应力分布

图8为初始地应力场主压应力迹线分布。从局部放大图上看到最大主应力均为水平方向，铅直方向为中间主应力，满足预期要实现的水平构造应力为最大主压应力的初始应力场特征量的条件。

图8 模型中部的主压应力边迹线分布

4 基于FLAC3D的狮子山矿区高水平构造应力场的反演与验证

4.1 狮子山矿区地应力场测试结果

分别在矿区4个不同中段对应的4个测点进行了实测，通过对实测数据综合分析，得出如下结果。

(1)4个测点的最大主应力均位于近水平方向，各测点的最大主应力方向均接近为NNW-SSE向(近垂直于矿体走向方向)，说明矿区以近水平构造应力场为主。

(2)各测点最大主应力、垂直主应力和最小主应力的大小均随深度基本成线性增长，并且垂直主应力和最小主应力的随埋深的增加速率要远小于最大主应力的增加速率，进一步证实了矿区最大主应力的影响具有明显的主导性。通过实测结果分析，得到最大主应力、中间主应力和最小主应力三个方向地应力(单位：MPa)随深度H(单位：m)的线性回归方程：

σh,max=-0.016 3+0.051 1H

(3)

σh,min=-0.388 3+0.021H

(4)

σv=-0.601 8+0.030 7H

(5)

式中：σh,max为最大主应力；σh,min为最小主应力；σv为垂直主应力。

(3)结合室内岩石力学试验结果和三维地应力实测结果，按照高地应力定量标准[18]，狮子山铜矿属于高地应力区。

4.2 狮子山铜矿地应力场反演数值模型的构建与计算参数

计算模型以矿体为中心，计算域上边界取至地表。计算模型X方向(垂直矿体走向方向)长度1 400 m；模型Y方向(沿矿体走向方向)，长度1 500 m；模型Z方向为竖直方向，模型底部标高687 m，模型最大高度1 430 m。模型共包括394 713个单元体和411 825个节点，计算模型和矿体形态分别如图9、图10所示。采用Mohr-Coulomb弹塑性本构模型[15]。计算模型X方向两端约束X方向位移，模型Y方向两端约束Y方向位移，模型底部固定位移，模型顶部为自由边界。地应力按3.4节讨论的高水平构造应力场的施加方法进行反演，即在位移边界条件下，对模型内所有单元体采用三角形应力分布初值设置，最大主应力为X方向，最小主应力为Y方向，竖直方向为中间主应力。矿区各岩性宏观岩体力学参数如表1所示。

图9 矿区三维有限差分数值模型

图10 矿体形态

图11分别为初始应力平衡后模型最大主应力和最小主应力分布。初始应力场反演计算表明，X方向(垂直于矿体走向方向)为最大主应力方向，Y方向(沿矿体走向方向)为最小主应力，Z方向(竖直方向)为中间主应力方向，与现场三维地应力实测主应力方向一致。根据三维地应力实测的4个测点位置，在模型内部相应位置布设4个监测点，以分析初始应力平衡后各监测点的主应力数值与地应力实测值之间的关系。

图11 初始应力平衡后模型主应力分布

为检验狮子山矿区地应力反演结果的可靠性，绘制实测地应力与反演地应力对比曲线，如图12所示。各测点最大主应力反演值与实测地应力值之间的最大相对误差为-8.10%，为测点D3；4个测点中，D1测点最小主应力反演值与实测地应力值之间的相对误差最大，为-18.29%。各测点中，中间主应力反演值与实测地应力之间最大相对误差为9.25%，为测点D3。综上可见，该地应力场的反演方法具有一定的可靠性。

图12 狮子山矿区4个测点实测地应力与反演地应力对比

5 结论

(1)在地应力边界条件讨论的基础上，采用FLAC3D对4种地应力场反演方法和边界条件进行数值模拟对比分析，结果表明位移边界条件在模拟自重应力场中效果较好，但不能较好地模拟高水平构造应力场。混合边界条件既不适于模拟自重应力场，也不适于模拟水平构造应力场。

(2)采用基于初始应变能理论的初始地应力场的反演方法，即在位移边界条件下，对计算域内所有单元体采用梯形应力分布初值设置，即设定初始弹性应变能状态，使满足或拟合应力随深度变化的实测值，模拟效果与实际吻合最好，能够较好地模拟初始高水平构造应力场。

(3)建立了狮子山矿区复杂三维力学模型，以初始应变能理论的初始地应力场的反演方法进行设置，数值反演表明，各测点主应力的方向与实测地应力方向一致，各测点主应力的反演值与实测地应力回归值之间的误差在允许范围内，表明矿区水平构造应力场反演具有一定的可靠性。