含层面岩体初始地应力场的应力函数解析法
吴锦亮，陈胜宏，何吉
(武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉，430072)
摘要：从应力函数解析法出发，考虑初始地应力场在岩层分界面处的应力不连续性，采用分片四次多项式应力函数对含层面岩体的初始地应力场进行分区拟合。针对多项式中的待定系数，引入地应力测点，并考虑地表面和岩层分界面处的应力边界条件，建立矛盾方程组，利用最小二乘法进行求解。为验证其合理性，假设某简化岩质边坡的边界条件和力学参数均为已知，利用有限单元法计算出“真实初始地应力场”；然后，取每个片区内若干节点处的应力作为实测值；最后，应用应力函数解析法求得“反演初始地应力场”。研究结果表明：在弹性范围内，使用该方法可快速有效地反演含层面岩体的三维初始地应力场。
关键词：初始地应力场；应力不连续性；分区拟合
中图分类号：P553 文献标志码：A 文章编号：1672-7207(2014)11-3922-08
Stress function analytic method for initial stress fields of rock masses with strata interfaces
WU Jinliang, CHEN Shenghong, HE Ji
(State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University,
Wuhan 430072, China)
Abstract: Based on stress function analytic method, piecewise quartic polynomial functions were established to analyze three dimensional initial stress fields of rock masses with strata interfaces in elastic range. Inconsistent equations were set up with measuring points and stress boundary conditions, and solved with the least square algorithm to determine the unknown parameters in the equations. A case study of a simplified slope was conducted to validate the method’s rationality. The slope was divided into four regions and each of them had some measuring points and some stress boundary conditions given. The slope “initial stress field” was analyzed with the proposed method. The method’s effectiveness and speediness is proved by comparing the analyzed results and the “real” initial stress field.
Key words: initial stress field; discontinuity strain; partition fitting
初始地应力对岩体位移以及其破坏形式有重要影响，是重要的地质参数[1-2]。岩体初始地应力最直接的获取途径是现场实测地应力，但由于受到测量费用和地形等种种因素的限制，在实际工程中只能对少数几个测点进行地应力测量。实测应力仅能反映测点附近岩体的局部应力特征[3]，因此，为了利用有限的测点来获取工程区域较大范围内的地应力场分布，进行岩体初始地应力场的拟合研究十分必要。此外，不同测量方法和测量误差使得实测应力结果具有一定的离散 性[3]，利用拟合方法可以予以修正。目前，利用实测应力资料进行初始地应力场拟合的常用方法主要有：1) 回归分析法[4-5]；2) 边界荷载法[6-7]和边界位移法[8-9]；3) 位移函数解析法[10]；4) 应力函数解析法[11]。在这些方法中，前3种方法都建立在有限元分析的基础上，地应力场拟合结果精确，但同时存在需要建立复杂的有限元模型、计算量大、耗时长的缺点。与前3种方法相比，应力函数解析法无需建立有限元模型，计算量小，且初始应力为构造的解析公式，使用方便，它特别适用于工程前期的地应力场定性分析，能根据地质勘查资料的调整变化而快速对初始地应力场进行重新拟合计算。应力函数解析法由张有天等[11]提出，即在计算域内假定一个带未知参数的应力函数(四次多项式)，使通过应力函数计算所得应力与已知实测值及边界条件相吻合，以确定未知参数。该法所求应力自然满足静力平衡方程和变形协调方程，并且在测点和表面点处与初始应力保持最小差值平方和的误差。张有天等[11]在地应力场拟合时仅使用一个应力函数，要求所研究区域岩体岩性必须单一，且必须满足应力场和位移场连续的条件。对于含多种介质的岩体，岩体初始地应力场在软弱岩层分界面处存在应力不连续(1个法向应力及2个沿分界面方向的切向应力连续，其余3个应力分量不连续)。为解决这个问题，张有天等[11]提出利用多个应力函数进行分区拟合的思想，但还未实现。据张有天等[11]的思路，本文作者通过引入分片四次多项式应力函数对含层面岩体的三维初始地应力场进行分区拟合，同时保证分界面上3个应力分量连续，从而得到一个统一却不连续的三维地应力场。
1 理论基础
根据弹性力学[12]理论，对于一个岩性单一且均匀的岩体区域，在弹性范围内，其内部各点都处在静力平衡状态，并且在变形后保持连续。也即该区域的应力函数必须满足平衡微分方程及变形协调方程。
对于空间问题，采用Maxwell应力函数，区域内的应力分量为
(1)
式中：
，
和
为应力函数；
为岩体容重，MN/m3。
上式自动满足空间问题的力平衡方程：
(2)
式中：
。
在体积力
为常数的情况下，上述应力分量应满足Beltrami-Michell方程，即用应力表示的变形协调方程：
(3)
式中：
为泊松比；
；
。
2 分片应力函数的构造
利用分片四次多项式应力函数来解决岩体初始地应力场在岩层分界面处的应力不连续问题。假定研究区域存在nd种介质，由ns个岩层分界面切割而成，包围区域的有ne个边界面(其中包括nb个地形面)。图1所示为岩体区域，nd=4，ns=3，nb=4。
实际工程区域中，地形面和结构面往往不规则，复杂多变，可根据它们的大致变化趋势，将地形面或结构面进行曲面的拟合，也可以是多个折面的组合。
图1 边坡岩体分区
Fig. 1 Partition of rock mass
对任一子域i(i=1, 2, 3, …, nd)，选取3个三元四次多项式应力函数
，
，
(式(4))，考虑到6个变形协调方程(式(3))，可消去6个待定系数，再考虑应力函数与应力分量的特定关系(式(1))，又可消去若干待定系数，最后待定的未知待定系数为36个(式(5))。以下式子为分片应力函数的一般表达式：
(4)
(5)
(6)
根据检验，式(4)~(6)能自动满足式(3)。将式(4)~(6)代入式(1)，可得子域i任意一点的6个应力分量：
(7)
(8)
(9)
(10)
(11)
(12)
3 分片应力函数待定系数的求解
由式(7) ~(12)可知，任意一个子域i(i=1, 2, 3, …, nd)的地应力场解析表达式存在36个待定系数，需要建立至少含有36个方程的线性方程组才能进行求解。以某子域二次初始地应力场x方向的正应力和剪应力为例，正应力表达式(7)分别由x，y和z 3个方向的二次多项式组成，其待定系数的确定需在3个方向上分别进行二次曲线拟合，同理，确定剪应力表达式(10)的待定系数则应在x方向进行二次曲线拟合和y，z等2个方向分别进行线性拟合。因此，选取的地应力测点和边界点的空间位置布置及数量能否满足以上条件是求解分片应力函数待定系数的关键，从数学角度容易发现，随着选点数目的增多将有益于曲线拟合精度的提高，同时又能使方程数目达到要求，但是，当子域内部无测点时，单独利用边界条件却很难满足曲线拟合的条件。
根据以上条件建立线性方程组，当方程个数大于待定系数个数时，可以利用最小二乘法进行求解，所拟合的地应力场能反映三维应力变化的大致趋势，且不要求反演所得地应力值与测值完全吻合，同时，由于考虑了岩层分界面处的应力协调条件，将使不同子域的地应力场成为相互影响的统一整体。
3.1 地应力测点
若子域i的地应力测点为m个，而每个地应力监测点j(j=1, 2, 3, …, m)可得到pj个应力分量。不同地应力测试方法所确定的应力分量不同，例如传统水压致裂法[13]只能测得水平主应力σhmax和σhmin的大小和方向(
, 
, 
这3个应力分量可由2个水平主应力转换得到)，而铅直向主应力的值则是由岩石的密度按静岩压力估算得到，不参与反演；孔径变形法[14]在多孔交汇条件下可测得6个应力分量，，
，，
和
。由该子域测点应力资料共得到以下
个方程：
(13)
3.2 地形面
在子域i的地形面选取n个点，已知这些点的法向应力和切应力为零，设选定的第j个地面点的法向方向余弦为
，
和
，与其共同构成笛卡尔坐标系的另外2个坐标轴方向的方向余弦分别为：
，
；
，
和
。则在该子域地形面处，有
(14)
3.3 岩层分界面
在岩层分界面j上选取sj个点，根据分界面处应力不连续特性，即1个法向应力及2个沿分界面方向的切向应力连续，其余3个应力分量不连续，存在以下关系式：
；
；
(15)
式中：“+”代表从坐标正向趋近分界面；“-”代表从坐标负向趋近分界面。
从ns个分界面可以总共得到
个方程式。
3.4 联立方程与求解
把式(7)~(12)代入式(13)~(15)后，可联立方程组：
(16)
式中：[A]为
阶矩阵；[X]为系数矩阵；[B]为右端项。
当R=L时，式(16)有唯一解；当R＞L时，可用最小二乘法求出系数矩阵[X]。当系数矩阵求解后，每个岩体子域内的地应力即可由式(7)~(12)解析表示。
4 算例考核与分析
假定某研究区域为一简化边坡(如图1)，边坡由片麻岩构成，包含4个子域和3个岩层分界面，模型以东西向为x轴，指向东为正；以南北向为y轴，指向北为正；以铅直向为z轴，指向上为正。边坡有限元网格如图2所示，3个层面ABC统一采用0.4 m厚的薄层单元，并顺沿薄层单元分别内置一组层理，边坡其余岩体单元统一内置两组节理(非贯通微小裂隙)。各子域节理岩体单元及层面处的薄层单元皆采用陈胜宏等[15-16]提出节理岩体等效流变模型进行等效模拟。边界条件及参数大小如图3所示，其中O点为三维空间坐标原点，所有边界约束点受到活动铰链约束，且边坡受到Z向重力荷载和X向构造作用，其中X向构造应力随深度的变化取线性表达式P=2+0.01h (h为埋深)。边坡岩体及结构面力学参数与产状如表1和表2所示。其中表2中的结构面倾向计算结果与地质记法结果一致。
由已知参数及边界条件，利用有限元弹性计算得到边坡“真实初始地应力场”。假定采用孔径变形法在多孔交汇条件下测量地应力(在每个测点处可测得6个应力分量)，每个子域分别选取3个测点，测点位置如表3所示。考虑到各子域的地形面和分界面均为规则的矩形面，采用大致均匀的取点原则，在这些面上统一选取9个控制点(见图4)。此外，在每个子域内分别预留3个点作为反演精度后验的校核点，不参与反演。
图2 边坡有限元网格
Fig. 2 Finite element mesh of slope
表1 各子域岩体力学参数
Table 1 Mechanical parameters of rock masses
表2 结构面力学参数与产状
Table 2 Mechanical parameters and attitude of joint
图3 几何参数及边界条件
Fig. 3 Geometry dimension and boundary conditions
图4 地形面及分界面取点位置
Fig. 4 Position of control points
利用以上资料对边坡的初始地应力场进行反演。表4和表5所示为各子域测点及校核点处3个坐标轴方向的正应力对比结果。此外，为直观反映边坡地应力场的变化趋势，给出边坡有限元初值与反演值在y=10 m断面处的图形化比较结果，如图5所示。
表3 各子域测点及校核点位置坐标
Table 3 Coordinates of measuring points and checking points
表4 各子域测点正应力及反演误差
Table 4 Normal stress values and normal stress errors of measuring points
表5 各子域校核点正应力及反演误差
Table 5 Normal stress values and normal stress errors of checking points
图5 y=10 m断面处正应力等值线对比图
Fig. 5 Comparison of normal stress isolines for y=10 m
根据表4和表5的对比结果可知：正应力测值和反演值很接近，最大反演误差在15.78%。此外，从图5可见：反演结果反映了边坡初始地应力场的大致变化趋势，而且在岩层分界面处体现了应力的不连续性。因此，对于本算例的边坡，采用分片多项式应力函数反演初始地应力场的效果令人满意。
5 结论
1) 引入分片四次多项式应力函数来反演含层面岩体的三维初始地应力场。该法原理简单，计算速度快，无需建立有限元模型，在弹性范围内，可用于反演含层面岩体的三维初始地应力场。反演结果可反映岩体地应力场在岩层分界面处的应力不连续特性，且子域内任意一点的应力都可由解得的分片应力函数解析表示。
2) 四次多项式应力函数只能描述二次应力场，若要解析高次应力区域，则需采用更高阶的应力函数。
3) 该法特别适用于工程前期的地应力场定性分析，能根据地质勘查资料的调整变化而快速对初始地应力场进行重新拟合计算。但由于实际工程中的岩体地应力场受到地下水位、塑性以及地形地貌等众多因素的影响，复杂多变，若要精确拟合，还需要借助传统的有限元方法。
参考文献：
[1] 杨林德. 岩土工程问题的反演理论与工程实践[M]. 北京: 科学出版社, 1996: 64-147.
YANG Linde. Back analysis theory and application to geotechnical engineering[M]. Beijing: Science Press, 1996: 64-147.
[2] 陈胜宏. 计算岩体力学与工程[M]. 北京: 中国水利水电出版社, 2006: 71-79.
CHEN Shenghong. Computational rock mechanics and engineering[M]. Beijing: China Water Power Press, 2006: 71-79.
[3] 葛修润, 侯明勋. 三维地应力BWSRM测量新方法及其测井机器人在重大工程中的应用[J]. 岩石力学与工程学报, 2011, 30(11): 2161-2180.
GE Xiurun, HOU Mingxun. A new 3D in situ rock stress measurement method: borehole wall stress relief method (BWSRM) and development of geostress measuring instrument based on BWSRM and its primary applications to engineering[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(11): 2161-2180.
[4] 郭怀志, 马启超, 薛玺成, 等. 岩体初始应力场的分析方法[J]. 岩土工程学报, 1983, 5(3): 64-75.
GUO Huaizhi, MA Qichao, XUE Xicheng, et al. The analytical method of the initial stress field for rock masses[J]. Chinese Journal of Geotechnical Engineering, 1983, 5(3): 64-75.
[5] 朱伯芳. 岩体初始地应力反分析[J]. 水利学报, 1994(10): 30-35.
ZHU Bofang. Back analysis of initial stresses in rock masses[J]. Journal of Hydraulic Engineering, 1994(10): 30-35.
[6] 白世伟, 李光煜. 二滩水电站坝区岩体应力场研究[J]. 岩石力学与工程学报, 1982, 1(1): 45-56.
BAI Shiwei, LI Guangyu. Study on stress field of Ertan hydroelectric station[J]. Chinese Journal of Rock Mechanics and Engineering, 1982, 1(1): 45-56.
[7] 易达, 陈胜宏, 葛修润. 岩体初始应力场的遗传算法与有限元联合反演法[J]. 岩土力学, 2004, 25(7): 1077-1080.
YI Da, CHEN Shenghong, GE Xiurun. Method of back analysis combining genetic algorithm with finite element method to initial stress field of rock mass[J]. Rock and Soil Mechanics, 2004, 25(7): 1077-1080.
[8] 金伟良, 喻军华, 邹道勤. 模拟高边坡初始地应力场的边界位移法[J]. 土木工程学报, 2003, 36(10): 72-75.
JIN Weiliang, YU Junhua, ZOU Daoqin. A new method of simulating initial ground stress field in high slope[J]. China Civil Engineering Journal, 2003, 36(10): 72-75.
[9] 郭明伟, 李春光, 王水林, 等. 优化位移边界反演三维初始地应力场的研究[J]. 岩土力学, 2008, 29(5): 1269-1274.
GUO Mingwei, LI Chunguang, WANG Shuilin, et al. Study on inverse analysis of 3D initial geostress field with optimized displacement boundaries[J]. Rock and Soil Mechanics, 2008, 29(5): 1269-1274.
[10] 喻军华, 金伟良, 邹道勤. 分析初始地应力场的位移函数法[J]. 岩土力学, 2003, 24(3): 417-419.
YU Junhua, JIN Weiliang, ZOU Daoqin. Displacement function method for analyzing initial earth stress[J]. Rock and Soil Mechanics, 2003, 24(3): 417-419.
[11] 张有天, 胡惠昌. 地应力场的趋势分析[J]. 水利学报, 1984(4): 31-38.
ZHANG Youtian, HU Huichang. Trend analysis of residual stress distribution in rock mass[J]. Journal of Hydraulic Engineering, 1984(4): 31-38.
[12] 徐芝纶. 弹性力学简明教程[M]. 北京: 高等教育出版社, 2002: 157-169.
XU Zhilun. Concise course of elasticity mechanics[M]. Beijing: Higher Education Press, 2002: 157-169.
[13] 蔡美峰. 地应力测量原理和方法的评述[J]. 岩石力学与工程学报, 1993, 12(3): 275-283.
CAI Meifeng. Commentary of principle and techniques of rock stress measurement[J]. Chinese Journal of Rock Mechanics and Engineering, 1993, 12(3): 275-283.
[14] 李方全, 廖椿庭. 地应力测量实例[J]. 地质科技, 1973(3): 67-73.
LI Fangquan, LIAO Chunting, Cases of geostress measurement[J]. Geological Science, 1973(3): 67-73.
[15] 陈胜宏, Peter Egger, 熊文林. 加锚节理岩体流变模型及三维弹粘塑性有限元分析[J]. 水利学报, 1998(9): 41-47.
CHEN Shenghong, Peter Egger, XIONG Wenlin. Rheological model and three dimensional elasto-viscoplastic finite element analysis of jointed rock masses reinforced by bolts[J]. Journal of Hydraulic Engineering, 1998(9): 41-47.
[16] CHEN Shenghong, Peter Egger. Three dimensional elasto- viscoplastic finite element analysis of reinforced rock masses and its application[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1999, 23(1): 61-78.
(编辑 邓履翔)
收稿日期：2013-12-25；修回日期：2014-03-02
基金项目(Foundation item)：国家自然科学基金资助项目(51079109) (Project(51079109) supported by the National Natural Science Foundation of China)
通信作者：陈胜宏(1957-)，男，浙江宁波人，博士，教授，从事水工结构及岩土工程的教学与研究；电话：027-68772086； E-mail: chensh@whu.edu.cn
