隧道围岩最大可移动块体求解及其应用

张敏思; 杨勇; 王述红; 查文华

doi:10.11779/CJGE20230826

隧道围岩最大可移动块体求解及其应用 English Version

张敏思^{1, 2,},
杨勇^{1, 2},
王述红^2, ,,
查文华¹

1.
东华理工大学江西省地质环境与地下空间工程研究中心, 江西南昌 330013
2.
东北大学资源与土木工程学院, 辽宁沈阳 110819

基金项目:

江西省主要学科学术和技术带头人培养计划 20212BCJ23003

江西省地质环境与地下空间工程研究中心开放基金资助项目 JXDHJJ2022-006

国家自然科学基金项目 52264003

国家自然科学基金项目 51964002

江西省“双千计划”人才资助项目 DHSQT22021002

详细信息

作者简介:
张敏思(1985—)，女，博士，副教授，主要从事非连续岩体力学方面的教学研究工作。E-mail: mscccathy@163.com

通讯作者:
王述红, E-mail: shwang@mail.neu.edu.cn

中图分类号: TU452
计量
- 文章访问数: 296
- HTML全文浏览量: 59
- PDF下载量: 71
出版历程
- 收稿日期: 2023-08-27
- 网络出版日期: 2024-03-24
- 刊出日期: 2024-10-31

Solutions and applications of maximum movable block in surrounding rock of tunnels

ZHANG Minsi^{1, 2,},
YANG Yong^{1, 2},
WANG Shuhong^2, ,,
ZHA Wenhua¹

1.
Engineering Research Center for Geological Environment and Underground Space of Jiangxi Province, East China University of Technology, Nanchang 330013, China
2.
School of Resources and Civil Engineering, Northeastern University, Shenyang 110819, China

摘要

摘要: 由于结构面的隐蔽性，目前的勘察技术无法精确获取每条结构面的定位及力学参数，仅能根据出露情况获取产状信息。而结构面的定位是求解块体的关键因素，因此有必要考虑不同定位情况下结构面与开挖面的所有组合，并找到最不利情况，为工程开挖提供预测及评价。在构建节理锥的基础上，提出了一种投影平移法，实现了任意柱状开挖面上最大可移动区域的求解。在此基础上，采用离散-切割-组合方法重构了曲面块体。首先，设置一组放射状虚拟结构面对锥形块体进行切割，完成块体的单元离散；其次，在面单元分类的基础上提出了曲面切割单元的详细算法。最后对切割后的单元进行分类组合，实现了曲面块体的重构。基于VC++和OpenGL开发了三维可视化程序，通过算例和工程实例验证了算法的精度、适用性及鲁棒性。最大可移动块体程序可用于隧道走向的选择，为工程设计提供依据。
- 块体理论 /
- 最大关键块体 /
- 隧道工程 /
- 节理岩体
Abstract: Due to the concealment of the discontinuities, the current exploration technology can not accurately obtain the location and mechanical parameters of all discontinuous, but can only obtain the occurrence information according to the exposure situation. The location of the discontinuities is a key factor for generating a block, so it is necessary to consider all the combinations of the discontinuities and excavation surface, and find the most unfavorable conditions to provide prediction and evaluation for engineering excavation. Based on the construction of joint pyramids, a projection translation method is proposed to solve the maximum movable region for cylindrical excavation surfaces. The dispersing-cutting-assembling method is used to reconstruct the curve surface block. Firstly, a set of radial virtual structures are set up to cut the pyramid block to realize the dispersion of the block. Secondly, based on the classification of face elements, a detailed surface-cutting-element algorithm is proposed. Finally, the cut elements are classified and combined to realize the reconstruction of the curve surface block. A 3D visualization program is developed based on the VC++ and OpenGL. The accuracy, applicability and robustness of the algorithm are verified by a numerical example and an engineering example. The maximum movable block program can be used to select the tunnel direction, providing reference for engineering design.
- block theory /
- maximum key block /
- tunnel engineering /
- joint rock mass

HTML全文

0. 引言

混凝土面板坝是以堆石料或砂砾石料为主要筑坝材料，以钢筋混凝土面板及其分缝止水作为防渗结构的高土石坝坝型^[1]。2000年以来，国内外面板坝进入200 m级特高坝时代，巴西于2006年建成最大坝高202 m的Campos Novos面板坝^[2]；老挝于2005年建成Nam Ngum 3水电站，其面板坝高度达到220 m^[3]；马来西亚于2009年建成Bakun面板坝，最大坝高205 m^[4]；中国则相继建成了天生桥一级、洪家渡、三板溪、水布垭、江坪河、猴子岩等一系列200 m级高混凝土面板坝^[5]。

混凝土面板是面板坝防渗体系的关键防线，其应力变形特性是影响大坝安全的关键。作为弹性薄板结构，面板上作用的法向水压力使其与垫层料之间产生相对位移，从而产生沿接触面的切向摩擦力。显然，面板应力变形主要受控于支撑其工作的堆石料在水压力作用下的应力变形特性及面板与堆石料之间的摩擦接触特性。目前，堆石料的本构模型研究取得大量成果，已能够比较全面地考虑其剪胀剪缩^[6-7]、颗粒破碎^[8]、湿化流变^[9]等复杂特性。

接触面模型研究同样长期受到关注，20世纪60年代，Goodman等^[10]就针对裂隙岩体问题提出无厚度接触面单元，成为解决复杂接触问题的基本方法；70年代初期，Wayne等^[11]提出了模拟土体与混凝土接触面剪应力与剪切位移关系的双曲线模型，并广泛应用于岩土工程^[12]。近年来，接触面模型研究主要集中在模拟单调和循环荷载作用下接触面法向胀缩、峰后软化等复杂行为^[13-17]，这类模型的本质是模拟接触面附近土体强剪切带的应力变形特性^[18]。

值得指出的是，上述接触面模型大多基于平面应变剪切试验建立，即假定剪应力和相对剪切位移是沿固定方向。实际面板坝工程大多修建于V形河谷，三维效应比较显著，面板与垫层之间接触应力和相对位移量值与方向的变化远较平面应变条件复杂。Evgin等^[19]曾开展接触面双向剪切试验，发现两个剪切方向耦合效应显著。本文简要介绍了双曲线接触面模型及其三维形式；提出了三维条件下接触面模型应满足的条件；提出了一个考虑双向剪切特性的接触面模型，并以水布垭面板坝为例，展示了接触面模型在面板坝应力变形计算中的应用。

1. 双曲线接触面剪切模型

Wayne等^[11]假定恒定法向应力σ_n条件下，接触面相对剪切位移u和平均剪应力τ之间符合双曲线关系，如图 1所示，即

$\tau = \frac{u}{{(1/{G_0}) + (u/{\tau _\infty })}} {\rm{。}}$

(1)

图 1 堆石料与混凝土接触面试验结果

Figure 1. Test results of rockfill-concrete interface

下载: 全尺寸图片幻灯片

式中： ${G_0}$ ， ${\tau _\infty }$ 分别为初始剪切模量和剪应力渐进值（图 1），

$\left. \begin{array}{l}{G}_{0}=\underset{u\to 0}{\mathrm{lim}}\left(\frac{\tau }{u}\right)\text{ }\text{, }\\ {\tau }_{\infty }=\underset{u\to \infty }{\mathrm{lim}}\tau \text{ }。\end{array} \right\}$

(2)

式（1）两侧对剪切位移u求导可得接触面切线剪切模量G_t，即

${G_{\text{t}}} = \frac{{{\text{d}}\tau }}{{{\text{d}}u}} = {\left( {1 - \frac{\tau }{{{\tau _\infty }}}} \right)^2}{G_0} = {\left( {1 - {R_{\text{f}}}\frac{\tau }{{{\tau _{\text{f}}}}}} \right)^2}{G_0} {\rm{。}}$

(3)

式中：R_f为破坏比，R_f = τ_f/τ_∞；初始剪切模量G₀与接触面法向应力σ_n有关，

${G_0} = k \cdot {\gamma _{\text{w}}} \cdot {\left( {\frac{{{\sigma _{\text{n}}}}}{{{p_{\text{a}}}}}} \right)^n} {\rm{。}}$

(4)

式中： ${\gamma _{\text{w}}}$ 为水的重度；p_a为大气压力；k，n为参数。

式（3）中， ${\tau _{\text{f}}}$ 为给定法向应力时接触面的破坏剪应力，一般采用Coulomb破坏准则，即

${\tau _{\text{f}}} = c + {\sigma _{\text{n}}} \cdot \tan \delta {\rm{。}}$

(5)

式中：c，δ分别为接触面的黏聚力和内摩擦角。

式（3）~（5）给出了单一方向剪切时接触面剪切模量表达式，对于双向剪切问题，需对其进行扩展。目前，最常用的扩展表达式为^{[12, 20]}

$\left. \begin{array}{l}{G}_{zx}={\left(1-{R}_{\text{f}}\frac{\left|{\tau }_{zx}\right|}{c+{\sigma }_{z}\mathrm{tan}\delta }\right)}^{2}\cdot k\cdot {\gamma }_{\text{w}}\cdot {\left(\frac{{\sigma }_{z}}{{p}_{\text{a}}}\right)}^{n}\text{, }\\ {G}_{zy}={\left(1-{R}_{\text{f}}\frac{\left|{\tau }_{zy}\right|}{c+{\sigma }_{z}\mathrm{tan}\delta }\right)}^{2}\cdot k\cdot {\gamma }_{\text{w}}\cdot {\left(\frac{{\sigma }_{z}}{{p}_{\text{a}}}\right)}^{n}{\rm{。}}\end{array} \right\}$

(6)

式中：τ_zx，τ_zy为接触面两个正交方向的剪应力；σ_z为接触面法向应力；x，y分别为接触面的两个切向坐标，z为接触面法向坐标。这样，接触面的应力-相对位移关系式可以表示为

$\left( \begin{array}{l} {\text{d}}{\tau _{zx}} \hfill \\ {\text{d}}{\tau _{zy}} \hfill \\ {\text{d}}{\sigma _z} \hfill \end{array} \right) = \left[ {\begin{array}{*{20}{c}} {{G_{zx}}}&0&0 \\ 0&{{G_{zy}}}&0 \\ 0&0&{{K_{zz}}} \end{array}} \right] \cdot \left( \begin{array}{l} {\text{d}}{u_x} \hfill \\ {\text{d}}{u_y} \hfill \\ {\text{d}}{u_z} \hfill \end{array} \right) {\rm{。}}$

(7)

式中：u_x，u_y，u_z分别为接触面的两个切向相对位移和法向相对位移；G_zx，G_zy为接触面的剪切模量；K_zz为法向接触模量，接触面处于压紧状态时取一高值，使接触面两侧材料近似满足不可嵌入条件；接触面处于拉开状态时取一低值，使接触面两侧材料可以相互独立变形^{[12, 20]}。

式（6）是式（3）的实用性推广，当沿x，y方向的剪应力达到破坏标准时，剪切模量趋近于零；当其中任意方向剪切位移为零时，另一方向剪应力-相对位移关系退化至式（3）。从式（7）可以看出，上述三维接触面模型的两个切线方向是相互独立的，不存在相互影响；接触面法向对切向的影响则通过剪切模量中法向应力体现。

2. 三维接触面剪切模型基本要求

接触面单元的应力位移关系一般在正交局部坐标系中建立，现设整体坐标系为XYZ，接触面单元局部坐标系为xyz，如图 2所示。下面分析说明接触面单元模型应满足的基本要求。

图 2 接触面单元分析的整体和局部坐标系

Figure 2. Global and local coordinate system for interface element analysis

下载: 全尺寸图片幻灯片

2.1 破坏准则

接触面双向剪切试验结果表明^[19]，其峰值强度与应力路径无关，且在（τ_zx，τ_zy）平面上，破坏包线近似为圆形，其半径与法向应力σ_z成正比。因此，在三维应力空间中，接触面单元的破坏面是一个圆锥面，锥面顶角为2δ，如图 3所示。这一破坏准则要求，当应力状态点位于破坏面并进一步增大剪切位移时，其各向剪切模量应趋近于零，即

${G_{zx}}\left( {\frac{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }}{{c + {\sigma _{\text{n}}}\tan \delta }} = 1} \right) = {G_{zy}}\left( {\frac{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }}{{c + {\sigma _{\text{n}}}\tan \delta }} = 1} \right) \to {0_{}} {\rm{。}}$

(8)

图 3 应力空间中接触面单元破坏面

Figure 3. Failure surface of interface in stress space

下载: 全尺寸图片幻灯片

显然，式（6）并不满足式（8）的条件。在式（6）中，当且仅当其各向剪应力分别达到破坏标准时，该方向的剪切模量趋近于零。因此，在应力空间中，其破坏面是顶面为正方形的锥面，故式（6）高估了接触面的抗剪强度。

2.2 剪切耦合效应

接触面受双向剪切荷载作用时，两个剪切方向的耦合效应已由试验证实^[19]，这种剪切耦合效应可通过简单的力学分析加以说明，如图 4所示。设想一接触面受恒定法向应力作用，现对该接触面施加x方向剪应力至应力状态A，该过程不产生y向相对位移；然后保持x向剪应力不变，持续增加y向剪应力，直至应力状态达到破坏面上B点。达到B点后，y向剪应力将不再继续增加，但相对位移将持续增加。显然，应力状态B点处产生的塑性位移不可能仅有y向分量，必然还存在x向分量，接触面塑性滑移方向应与剪应力合力方向一致。y向剪应力施加过程是一个连续过程，接触面的位移响应是连续的，因此不可能在B点处突然出现x向位移增量。也就是说，即使保持x向剪应力不变，连续施加y向剪应力的过程中，x向位移增量应贯穿整个y向加载过程（应力路径AB）。如前所述，式（6），（7）给出的传统模型无法反映上述剪切耦合效应。

图 4 接触面双向剪切耦合效应示意图

Figure 4. Coupling effects in two-way shearing

下载: 全尺寸图片幻灯片

2.3 标架无差异原则

连续介质力学中，有所谓的标架无差异原则，即材料的应力应变关系是其固有属性，不应受到所选取坐标系的影响^[21]。接触面单元同样应符合标架无差异原则，即接触面的剪应力-相对位移关系与所选取的两个切向坐标的方向无关。为说明该原则对接触面剪切模量矩阵的要求，考虑两个不同的切向坐标系xy和xʹyʹ，两坐标系具有相同的法向z，xʹyʹ坐标系是由xy坐标系逆时针旋转θ角度后得到的，如图 5所示。

图 5 接触面切向坐标系旋转示意图

Figure 5. Rotation of tangential directions

下载: 全尺寸图片幻灯片

现设xy坐标系中，接触面应力位移关系可表示为

$\left. \begin{array}{l}\text{d}\mathit{\boldsymbol{\tau}} =\mathit{\boldsymbol{G}}\cdot \text{d}\mathit{\boldsymbol{u}}\text{ }\text{, }\\ \left(\begin{array}{l}\text{d}{\tau }_{zx}\\ \text{d}{\tau }_{zy}\end{array}\right)=\left[\begin{array}{cc}{G}_{xx}& {G}_{xy}\\ {G}_{yx}& {G}_{yy}\end{array}\right]\left(\begin{array}{l}\text{d}{u}_{x}\\ \text{d}{u}_{y}\end{array}\right)\text{ }{\rm{。}}\end{array} \right\}$

(9)

类似地，xʹyʹ坐标系中，接触面的应力位移关系可以表示为

$\left. \begin{array}{l}\text{d}{\mathit{\boldsymbol{\tau}} }^{\prime }={\mathit{\boldsymbol{G}}}^{\prime }\cdot \text{d}{\mathit{\boldsymbol{u}}}^{\prime }\text{ }\text{, }\\ \left(\begin{array}{l}\text{d}{{\tau }^{\prime }}_{zx}\\ \text{d}{{\tau }^{\prime }}_{zy}\end{array}\right)=\left[\begin{array}{cc}{{G}^{\prime }}_{xx}& {{G}^{\prime }}_{xy}\\ {{G}^{\prime }}_{yx}& {{G}^{\prime }}_{yy}\end{array}\right]\left(\begin{array}{l}\text{d}{{u}^{\prime }}_{x}\\ \text{d}{{u}^{\prime }}_{y}\end{array}\right)\text{ }{\rm{。}}\end{array} \right\}$

(10)

式（9），（10）中: G，Gʹ分别为两坐标系中的剪切模量矩阵。若记

$\mathit{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&{\sin \theta } \\ { - \sin \theta }&{\cos \theta } \end{array}} \right] \text{, }$

(11)

则

$\left. \begin{array}{l}\text{d}{\mathit{\boldsymbol{\tau}} }^{\prime }=\mathit{\boldsymbol{R}}\cdot \text{d}\mathit{\boldsymbol{\tau}} \text{ }\text{, }\\ \text{d}{\mathit{\boldsymbol{u}}}^{\prime }=\mathit{\boldsymbol{R}}\cdot \text{d}\mathit{\boldsymbol{u}}\text{ }{\rm{。}}\end{array} \right\}$

(12)

将式（12）代入式（10），并将其与式（9）对比可知，两局部坐标系中剪切模量矩阵应满足下面的关系：

$\left. \begin{array}{l}\mathit{\boldsymbol{G}}={\mathit{\boldsymbol{R}}}^{\text{T}}\cdot {\mathit{\boldsymbol{G}}}^{\prime }\cdot \mathit{\boldsymbol{R}}\text{ }\text{, }\\ {\mathit{\boldsymbol{G}}}^{\prime }=\mathit{\boldsymbol{R}}\cdot \mathit{\boldsymbol{G}}\cdot {\mathit{\boldsymbol{R}}}^{\text{T}}\text{ }{\rm{。}}\end{array} \right\}$

(13)

显然式（6），（7）给出的剪切模量矩阵一般情况下并不符合式（13）中的转换关系，即其应力位移关系受到局部坐标系选择的影响。

为更清楚地说明上述问题，设一接触面受到恒定法向应力σ_z = 500 kPa作用，其后沿x向施加剪应力至τ_zx = 400 kPa。现采用式（6）分别在xy坐标系和xʹyʹ坐标系中研究其位移路径，相关参数如表 1所示^[20]。

表 1 堆石料-混凝土接触面计算参数

Table 1. Parameters of rockfill-concrete interface

R_f	k	n	c/ kPa	δ/(°)
0.90	6000	0.85	0.0	41.5

下载: 导出CSV

| 显示表格

图 6中给出了θ为0°（xy坐标系），30°，60°情况下接触面的应力位移关系和应力路径，在xy坐标系中考察该问题时，接触面沿x向位移达到9.37 mm；但在θ为30°，60°的局部坐标系中计算时，x向位移仅5.18 mm，可见式（6）高估了接触面剪切模量。更不合理之处在于，该案例仅沿x向施加荷载，但在θ为30°，60°的局部坐标系中计算时，会产生y向剪切位移。该案例清楚地表明，式（6），（7）预测的应力位移关系和位移路径与局部坐标系选择有关，违背了标架无差异原则。

图 6 不同局部坐标系下接触面应力位移关系与位移路径

Figure 6. Stress-displacement relations and displacement paths under different local coordinate systems

下载: 全尺寸图片幻灯片

3. 接触面的双向剪切弹塑性模型

根据前文分析，面板坝计算常用的接触面模型存在高估剪切强度、不能考虑剪切耦合效应以及违反标架无差异原则等不足，势必影响面板应力计算结果的可靠性。为解决上述问题，本文在弹塑性理论框架下，构建一个新的接触面模型。

3.1 基本假定

接触面剪切试验所揭示的剪胀剪缩特性本质上是其附近强剪切带内土体的力学行为，原则上可通过剪切带土体的本构模型模拟。因此，本文仅限于模拟所谓无厚度接触面单元在压紧状态时的剪切特性，并采用2点基本假定。

（1）不考虑接触面的剪胀与剪缩，接触面的法向和切向耦合效应仅通过法向应力对剪切模量的影响体现。

（2）接触面的剪切位移增量可分解为弹性位移增量和塑性位移增量，即

$\text{d}\mathit{\boldsymbol{u}}=\text{d}{\mathit{\boldsymbol{u}}}^{\text{e}}+\text{d}{\mathit{\boldsymbol{u}}}^{\text{p}}\text{, }\left(\begin{array}{l}\text{d}{u}_{x}\\ \text{d}{u}_{y}\end{array}\right)=\left(\begin{array}{l}\text{d}{u}_{x}^{\text{e}}\\ \text{d}{u}_{y}^{\text{e}}\end{array}\right)+\left(\begin{array}{l}\text{d}{u}_{x}^{\text{p}}\\ \text{d}{u}_{y}^{\text{p}}\end{array}\right) {\rm{。}}$

(14)

式中：上标e，p分别表示弹性和塑性。下面分别说明弹性和塑性剪切位移的计算方法。

3.2 弹性位移计算

接触面的弹性位移增量采用下式计算：

$\left( \begin{array}{l} {\text{d}}u_x^e \hfill \\ {\text{d}}u_y^{\text{e}} \hfill \end{array} \right) = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{G_{\text{e}}}}}}&0 \\ 0&{\frac{1}{{{G_{\text{e}}}}}} \end{array}} \right]\left( \begin{array}{l} {\text{d}}{\tau _{zx}} \hfill \\ {\text{d}}{\tau _{zy}} \hfill \end{array} \right) {\rm{。}}$

(15)

式中：G_e为接触面弹性剪切模量，

${G_{\text{e}}} = {k_{\text{e}}} \cdot {\gamma _{\text{w}}} \cdot {\left( {\frac{{{\sigma _{\text{n}}}}}{{{p_{\text{a}}}}}} \right)^n} {\rm{。}}$

(16)

式中：k_e为接触面弹性参数，可取k_e = (2~3)k；n为模量参数，与式（4）中的参数n一致。

3.3 塑性位移计算

接触面处于滑动状态时，其滑动方向应与剪应力合力方向一致。因此，本文假定塑性位移增量方向始终与剪应力方向一致，即

$\left( \begin{array}{l} {\text{d}}u_x^{\text{p}} \hfill \\ {\text{d}}u_y^{\text{p}} \hfill \end{array} \right) = \left( \begin{array}{l} \frac{{{\tau _{zx}}}}{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }} \hfill \\ \frac{{{\tau _{{\text{zy}}}}}}{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }} \hfill \end{array} \right){\text{d}}{u^{\text{p}}} = \left( \begin{array}{l} \frac{{{\tau _{zx}}}}{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }} \hfill \\ \frac{{{\tau _{zy}}}}{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }} \hfill \end{array} \right)\frac{1}{{{G_{\text{p}}}}}{\text{d}}\tau {\rm{。}}$

(17)

式中：du^p为塑性位移增量的合成量；G_p可称为接触面的塑性剪切模量；dτ为剪应力增量的合成量，即

$\begin{array}{l}\text{d}\tau =\text{d}\sqrt{{\tau }_{zx}{}^{2}+{\tau }_{zy}{}^{2}}\\ \ \ \ \ =\left(\frac{{\tau }_{zx}}{\sqrt{{\tau }_{zx}{}^{2}+{\tau }_{zy}{}^{2}}}\text{, }\frac{{\tau }_{zy}}{\sqrt{{\tau }_{zx}{}^{2}+{\tau }_{zy}{}^{2}}}\right)\left(\begin{array}{l}\text{d}{\tau }_{zx}\\ \text{d}{\tau }_{zy}\end{array}\right)\end{array} {\rm{。}}$

(18)

因此式（17）可以改写为

$\left( \begin{array}{l} {\text{d}}u_x^{\text{p}} \hfill \\ {\text{d}}u_y^{\text{p}} \hfill \end{array} \right) = \frac{1}{{{G_{\text{p}}}}}\left[ {\begin{array}{*{20}{c}} {\frac{{{\tau _{zx}}^2}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}}&{\frac{{{\tau _{zx}}{\tau _{zy}}}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}} \\ {\frac{{{\tau _{zy}}{\tau _{zx}}}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}}&{\frac{{{\tau _{zy}}^2}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}} \end{array}} \right]\left( \begin{array}{l} {\text{d}}{\tau _{zx}} \hfill \\ {\text{d}}{\tau _{zy}} \hfill \end{array} \right) {\rm{。}}$

(19)

3.4 本构方程

将式（15），（19）代入式（14）可得柔度形式的接触面本构方程，即

$\left( \begin{array}{l} {\text{d}}{u_x} \hfill \\ {\text{d}}{u_y} \hfill \end{array} \right) = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{G_{\text{e}}}}} + \frac{1}{{{G_{\text{p}}}}}\frac{{{\tau _{zx}}^2}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}}&{\frac{1}{{{G_{\text{p}}}}}\frac{{{\tau _{zx}}{\tau _{zy}}}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}} \\ {\frac{1}{{{G_{\text{p}}}}}\frac{{{\tau _{{\text{zy}}}}{\tau _{{\text{zx}}}}}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}}&{\frac{1}{{{G_{\text{e}}}}} + \frac{1}{{{G_{\text{p}}}}}\frac{{{\tau _{zy}}^2}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}} \end{array}} \right]\left( \begin{array}{l} {\text{d}}{\tau _{zx}} \hfill \\ {\text{d}}{\tau _{zy}} \hfill \end{array} \right) {\rm{。}}$

(20)

由式（20）可进一步写出劲度形式的接触面本构方程，这里直接给出式（9）中剪切模量矩阵的表达式为

${\mathbf{G}} = \left[ {\begin{array}{*{20}{c}} {{G_{\text{e}}}}&0 \\ 0&{{G_{\text{e}}}} \end{array}} \right] - \frac{{{G_{\text{e}}}^2}}{{{G_{\text{e}}} + {G_{\text{p}}}}}\left[ {\begin{array}{*{20}{c}} {\frac{{\tau _{zx}^2}}{{\tau _{zx}^2 + \tau _{zy}^2}}}&{\frac{{{\tau _{zx}} \cdot {\tau _{zy}}}}{{\tau _{zx}^2 + \tau _{zy}^2}}} \\ {\frac{{{\tau _{zy}} \cdot {\tau _{{\text{zx}}}}}}{{\tau _{zx}^2 + \tau _{zy}^2}}}&{\frac{{\tau _{zy}^2}}{{\tau _{zx}^2 + \tau _{zy}^2}}} \end{array}} \right] {\rm{。}}$

(21)

式中：右侧第一项满足式（13）；第二项方括号内可表示为τ·τ^T/|τ|²，应用式（12）的转换关系可以证实其满足式（13），故上述模型满足标架无差异原则。从推导过程可以看出，本文模型满足该原则的关键是假定弹性剪切模量各向同性，且塑性位移增量方向与剪应力方向一致。

3.5 塑性模量

完成本构模型构建工作，尚需给出塑性剪切模量表达式。为此，考虑沿一固定方向施加剪切荷载的情形，即

$\frac{{{\text{d}}{u_x}}}{{{\text{d}}{u_y}}} = \frac{{{\tau _{zx}}}}{{{\tau _{zy}}}} = \frac{{{\text{d}}{\tau _{zx}}}}{{{\text{d}}{\tau _{zy}}}} = r {\rm{。}}$

(22)

式中：r为比例系数。

将式（22）代入式（20）可得

${\text{d}}u\left( \begin{array}{l} r \hfill \\ 1 \hfill \end{array} \right) = \frac{{{\text{d}}\tau }}{{{G_{\text{e}}}}}\left( \begin{array}{l} r \hfill \\ 1 \hfill \end{array} \right) + \frac{{{\text{d}}\tau }}{{{G_{\text{p}}}}}\left( \begin{array}{l} r \hfill \\ 1 \hfill \end{array} \right) {\rm{。}}$

(23)

即

$\frac{1}{{{G_{\text{t}}}}} = \frac{{{\text{d}}u}}{{{\text{d}}\tau }} = \frac{1}{{{G_{\text{e}}}}} + \frac{1}{{{G_{\text{p}}}}} {\rm{。}}$

(24)

因此

${G_{\text{p}}} = {\left( {\frac{1}{{{G_{\text{t}}}}} - \frac{1}{{{G_{\text{e}}}}}} \right)^{ - 1}} = \frac{{{G_{\text{e}}} \cdot {G_{\text{t}}}}}{{{G_{\text{e}}} - {G_{\text{t}}}}} {\rm{。}}$

(25)

注意式（25）中切线剪切模量仍采用式（3）的基本形式，但将其扩展为二维剪切情形，即

${G_{\text{t}}} = {\left( {1 - {R_{\text{f}}}\frac{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }}{{c + {\sigma _z}\tan \delta }}} \right)^2} \cdot k \cdot {\gamma _{\text{w}}} \cdot {\left( {\frac{{{\sigma _z}}}{{{p_{\text{a}}}}}} \right)^n} {\rm{。}}$

(26)

显然，无论剪应力分量如何，只要当其合成量满足破坏准则时，切线剪切模量趋近于零（R_f = 1时），因此塑性模量也趋近于零，微小的应力增量即可引起很大的塑性位移。

3.6 加卸载准则

弹塑性模型的加卸载准则一般由屈服函数规定，但本文未定义屈服函数，故需补充加卸载准则。记接触面最大历史剪应力为τ_max，规定若接触面应力和位移增量同时满足下述条件

$\left. \begin{array}{l}\sqrt{{\tau }_{zx}^{2}+{\tau }_{zy}^{2}}\ge {\tau }_{\text{max}}\text{ }\text{, }\\ {\tau }_{zx}\text{d}{u}_{x}+{\tau }_{zy}\text{d}{u}_{y}\ge 0\text{ }{\rm{。}}\end{array} \right\}$

(27)

则接触面处于加载状态，剪切模量矩阵由式（21）给出；否则，接触面处于卸载状态，剪切模量矩阵仅含式（21）右端第一项。式（27）实际上规定：在剪切应力（τ_zx，τ_zy）平面上，若当前剪切应力状态点位于历史最大剪应力包线上，且剪切位移增量方向与剪应力合力方向夹锐角，则接触面处于加载状态；否则，按卸载处理。

3.7 模型验证

采用弹塑性理论，将传统平面应变接触面模型扩展至双向剪切情形，使其满足三维破坏准则和标架无差异原则，并能够考虑两个剪切方向的相互影响。除增加弹性剪切模量参数k_e外，未增加其它参数，且所有参数仍可通过恒定法向应力的单向剪切试验确定；k_e可按经验取(2~3)k。这里，首先采用Evgin等^[19]的双向剪切试验对该模型进行验证。

Evgin等针对密砂与粗糙钢板面开展了恒定法向应力的双向剪切试验，试验时首先施加法向应力至 ${\sigma _z}$ = 100 kPa；然后分别施加τ_zy至0，20，40，60 kPa；其后保持 ${\tau _{zy}}$ 恒定，以位移控制方式施加 ${\tau _{zx}}$ 直至接触面发生剪切破坏。图 7中给出了试验的应力路径和x向剪应力和位移关系，可以看出密砂与粗糙钢板接触面具有明显的峰后软化效应，但本文切线模量表达式（26）无法反映软化现象，故此处模拟接触面的峰值强度特性，表 2列出了根据试验结果确定的模型参数，图 7中曲线是应用表中参数模拟试验的结果，可以看出本文模型较好地反映了沿接触面一个切向施加剪应力后，另一切向剪切强度的降低。

图 7 接触面双向剪切试验结果^[19]

Figure 7. Two-way shear results of interface^[19]

下载: 全尺寸图片幻灯片

表 2 砂-钢板接触面计算参数

Table 2. Parameters of sand-steel interface

R_f	k/10³	k_e/10³	n	c/kPa	δ/(°)
1.0	150	300	1.0	0.0	38.7

下载: 导出CSV

| 显示表格

图 8对比了上述试验和模型计算得到的位移路径，其中不同形状点是试验结果；曲线是模拟结果。尽管模型计算的应力位移曲线与试验结果存在差异，但模型计算的位移路径与试验结果是相当相符的。由于接触面施加x向剪切位移后，很快达到峰值强度（图 7），其后进入软化阶段，故无论是x向剪切位移，还是y向剪切位移，均以塑性滑移占主导。因此，图 8中模拟结果与试验结果相符实际上表明，本文采用的塑性流动准则式（17）是可以接受的。

图 8 接触面双向剪切试验的位移路径^[19]

Figure 8. Displacement paths in two-way shear^[19]

下载: 全尺寸图片幻灯片

为进一步对比本文模型与传统模型（即式（6），（7））的差异，这里运用表 1中所列参数研究两个模型模拟等位移比剪切时的应力路径。假定接触面法向应力为1 MPa，沿7个不同的位移路径施加剪切位移至10 mm，如图 9（a）所示；运用两种模型计算的接触面剪应力路径绘制于图 9（b）中。本文模型计算结果显示，x向和y向剪切位移等比例增加时，两个方向剪应力τ_zx和τ_zy也等比例增加，应力路径与位移路径方向一致，最终应力状态点落入一圆弧上，即最终剪应力合力相等；传统模型因其两个方向剪切模量的不同步变化（式6），计算的剪应力路径明显偏离等应力比路径，最终应力状态点甚至落入强度包线之外，这显然是不合理的。

图 9 等位移比剪切应力路径

Figure 9. Stress paths in constant direction shear

下载: 全尺寸图片幻灯片

4. 接触面模型的运用

前文所述接触面模型仅涉及两个剪切方向，实际运用时仍需补充接触面法向本构关系。最简单的处理方式是按照Goodman等^[10]的方法，根据接触面压紧或张开状态，分别对法向模量赋不同量级的值。本节以水布垭面板坝为例，将本文模型运用于面板位移与应力计算。

4.1 工程概况

水布垭面板坝最大坝高233 m，是中国已建最高面板坝^[22]。大坝从上游至下游主要有5个材料分区，即垫层区（2A）、过渡区（3A）、主堆石区（3B）、次堆石区（3C）以及下游堆石区（3D），如图 10所示。实际施工时，大坝分4期填筑；面板分3期浇筑。计算模拟时对上述施工过程进行适当简化：①堆石料按图 10所示4个分期填筑；②从底至顶一次浇筑面板至405 m高程；③水库水位由176 m逐渐提升至400 m高程；④水库在400 m水位条件下运行5 a。

图 10 大坝材料分区与填筑过程

Figure 10. Dam material zones and filling processes

下载: 全尺寸图片幻灯片

4.2 模型与参数

根据坝体施工过程和面板分缝情况建立了大坝三维有限元模型，如图 11所示。面板和垫层料之间设置接触面单元；面板和趾板之间设置周边缝单元；面板和面板之间设置垂直缝单元。接触面单元参数如表 1所列，并取k_e = 2k；缝单元采用文献[12]中的模型与参数。模型的左右岸截断边界设置坝轴向（x向）位移约束；坝体与建基面接触的结点设置为三向约束。

图 11 大坝的三维有限元网格

Figure 11. Finite element mesh of dam

下载: 全尺寸图片幻灯片

堆石料采用“南水”双屈服面弹塑性本构模型^[7]，并考虑流变影响^[23]，其中垫层料、过渡料、主堆石料、下游堆石料采用相同参数；次堆石料采用另一组参数，如表 3，4所列。

表 3 筑坝堆石料“南水”模型参数^[22]

Table 3. 'Nanshui' constitutive model parameters of damming rockfill materials^[22]

材料	ρ/(g·cm^-3)	φ₀/(°)	Δφ/(°)	R_f	k	n	c_d/ %	n_d	R_d
2A, 3A, 3B, 3D	2.20	54.7	10.4	0.81	994	0.33	0.29	0.84	0.72
3C	2.20	51.3	10.4	0.83	602	0.25	0.28	0.98	0.75

下载: 导出CSV

| 显示表格

表 4 筑坝堆石料流变模型参数^[23]

Table 4. Creep model parameters of damming rockfill materials^[23]

材料	t₀/月	c₁/ %	n₁	c₂/ %	n₂	c₃/ %	n₃	n₄
2A, 3A, 3B, 3D	5.5	0.160	0.45	0.006	0.88	0.032	0.66	0.50
3C	5.5	0.180	0.43	0.007	0.90	0.035	0.55	0.50

下载: 导出CSV

| 显示表格

对于混凝土趾板和面板，采用各向同性线弹性模型模拟，弹性模量为E = 28 GPa；弹性泊松比为 $\nu$ = 0.167。

4.3 计算结果分析

图 12，13中分别绘制了用本文接触面模型计算的蓄水期和运行期面板轴向位移、法向位移（挠度）、轴向应力和顺坡向应力等值线分布。水库蓄水至400 m高程时，两岸面板均向河谷中央位移，因此在河床中部面板中产生轴向挤压应力，最大值约12.71 MPa，两岸坡部位面板则出现轴向拉应力，最大值约2.20 MPa；水压力作用下面板挠度从坝顶至坝底先增大后减小，最大挠度发生于1/3坝高处；面板顺坡向几乎处于全断面受压状态，最大压应力为14.84 MPa。

图 12 本文模型计算的蓄水后面板位移和应力分布等值线

Figure 12. Displacements and stresses of concrete slabs after impounding predicted with proposed model

下载: 全尺寸图片幻灯片

图 13 本文模型计算的运行期面板位移和应力分布等值线

Figure 13. Displacements and stresses of concrete slabs during operation predicted by proposed model

下载: 全尺寸图片幻灯片

运行期因坝体堆石料产生流变，两岸坡部位面板进一步向河床中央挤压，中部面板轴向压应力显著增大至35.15 MPa，峰值区上移至320 m高程以上；因堆石体沉降，面板挠度峰值区上移至坝顶部位，最大值达到84.5 cm，与实测结果基本接近^[24]，与此同时，面板顺坡向压应力也大幅增长至28.00 MPa。在该案例中，堆石流变产生的面板坝轴向压应力增幅超过顺坡向压应力增幅。该坝于2008年11月首次蓄水至正常蓄水位附近，运行过程中河床中部面板出现垂直缝两侧面板表层混凝土挤压破损，部位与本文所计算的压应力极值区基本吻合^[24]，表明本文接触面模型较好地反映了堆石坝体与面板之间的荷载传递机制。

5. 结论

分析了现有三维接触面模型的不足，提出了一个模拟混凝土面板与堆石料双向摩擦接触特性的无厚度接触面模型，该模型与传统模型及同类弹塑性模型相比，具有3个特点。

（1）满足接触面剪切强度准则和力学响应的标架无差异原则，且可以反映两个剪切方向的相互影响。满足标架无差异原则，使该模型实际应用时，计算结果与局部坐标系选择无关。

（2）由于未考虑剪切软化等复杂特性，模型参数仍比较简单，与传统接触面模型相比，仅增加一个弹性参数，只需通过常规的恒定法向应力的接触面剪切试验即可确定全部参数。

（3）模型的数值实现方便，只需简单修改传统接触面模型的剪切模量矩阵即可，且模型模拟结果与实际案例出现的情况定性上相符。

图 1 节理锥有效性判断

Figure 1. Validity judgment of joint pyramids

下载: 全尺寸图片幻灯片

图 2 切点的确定

Figure 2. Determination of tangential points

下载: 全尺寸图片幻灯片

图 3 锥体的构建

Figure 3. Construction of joint pyramid

下载: 全尺寸图片幻灯片

图 4 曲面块体的重构

Figure 4. Reconstruction of surface blocks

下载: 全尺寸图片幻灯片

图 5 面的类型示意图

Figure 5. Diagram of face type

下载: 全尺寸图片幻灯片

图 6 曲面切割单元块体

Figure 6. Curve surface cutting element block

下载: 全尺寸图片幻灯片

图 7 锥形块体的构建

Figure 7. Construction of pyramid blocks

下载: 全尺寸图片幻灯片

图 8 曲面块体重构

Figure 8. Reconstruction of curve surface blocks

下载: 全尺寸图片幻灯片

图 9 重构后和重构前参数的比值

Figure 9. Ratios of parameters before and after reconstruction

下载: 全尺寸图片幻灯片

图 10 最大块体形态

Figure 10. Shape of maximum block

下载: 全尺寸图片幻灯片

图 11 夹角与计算精度间的关系

Figure 11. Relationship between intersecting angle and calculation

下载: 全尺寸图片幻灯片

图 12 最大块体体积与隧道走向间的关系

Figure 12. Relationship between volume of maximum block and tunnel direction

下载: 全尺寸图片幻灯片

图 13 曲线极值处四面体形态

Figure 13. Tetrahedral morphologies at extreme of curve

下载: 全尺寸图片幻灯片

图 14 曲线极值处五面体形态

Figure 14. Pentahedral morphologies at extreme of curve

下载: 全尺寸图片幻灯片

图 15 存在交集的两锥体

Figure 15. Two pyramids with intersection

下载: 全尺寸图片幻灯片

图 16 锥体的合并

Figure 16. Conflation of pyramids

下载: 全尺寸图片幻灯片

表 1 面的类型

Table 1 Types of face

类别	判断方法	图示
类别	判断方法	A型面	B型面
Ⅰ类面	除落在曲面上的点外，其余顶点均在开挖面外侧
Ⅱ类面	除落在曲面上的点外，其余顶点均在开挖面内侧
Ⅲ类面	顶点分布在曲面内外两侧
Ⅳ类面	位于开挖面外，以边的方式与开挖面接触
Ⅴ类面	位于开挖面内，以边的方式与开挖面接触

下载: 导出CSV

表 2 结构面参数

Table 2 Parameters of the discontinuities

节理组	产状代表值	走向
F1	NE15°∠52°	NW105
F2	SW242°∠43°	NW152°
F3	SW203°∠39°	NW113°
F4	NW347°∠45°	NE77°

下载: 导出CSV

表 3 F1-F2-F3-F4构建的有效节理锥棱线向量

Table 3 Vectors of effective edges of joint pyramid constructed by F1-F2-F3-F4

交线	棱线向量（x, y, z）
交线	JP-A1	JP-A6	JP-A11	JP-A16
F1F2	—	-N₁	+N₁	—
F1F3	+N₂	-N₂	+N₂	-N₂
F1F4	-N₃	—	—	+N₃
F2F3	-N₄	—	—	+N₄
F2F4	+N₅	+N₅	-N₅	-N₅
F3F4	—	+N₆	-N₆	—
注	体积小			体积小

下载: 导出CSV

表 4 φ为5°，0°的体积比

Table 4 Volumetric ratios under φ of 5°and 0°

隧道走向θ	0°	45°	90°	135°
φ=5°	0.4438	0.3395	0.4134	0.3784
φ=0°	0.4404,	0.3357	0.4094	0.3748
差值百分比/%	0.77	1.13	0.978	0.96

下载: 导出CSV

参考文献(14)

[1]	GOODMAN R E, SHI G H. Block theory and its application to rock engineering[M]. Englewood Cliff: Prentice Hall, 1985.
[2]	PRIEST S D. Discontinuity Analysis for Rock Engineering[M]. London: Chapman & Hall, 1993.
[3]	DELPORT J L, MARTIN D H. A multiplier method for identifying keyblocks in excavations through jointed rock[J]. SIAM Journal on Algebraic and Discrete Methods, 1986, 7(2): 321-330. doi: 10.1137/0607035
[4]	LIN D, FAIRHURST C, STARFIELD A M. Geometrical identification of three-dimensional rock block systems using topological techniques[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1987, 24(6): 331-338.
[5]	YU Q, OHNISHI Y, XUE G, et al. A generalized procedure to identify three-dimensional rock blocks around complex excavations[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2009, 33(3): 355-375. doi: 10.1002/nag.720
[6]	ZHANG Z X, WANG S F, HUANG X, et al. Application of block theory for evaluating face stability under disc cutters loading of TBM, case study of a water-conveyance tunnel project[J]. Tunnelling and Underground Space Technology, 2019, 90: 249-263. doi: 10.1016/j.tust.2019.05.002
[7]	ZHANG Y T, XIAO M, CHEN J T. A new methodology for block identification and its application in a large scale underground cavern complex[J]. Tunnelling and Underground Space Technology, 2010, 25(2): 168-180. doi: 10.1016/j.tust.2009.10.005
[8]	ZHENG Y H, XIA L, YU Q C. A method for identifying three-dimensional rock blocks formed by curved fractures[J]. Computers and Geotechnics, 2015, 65: 1-11. doi: 10.1016/j.compgeo.2014.11.005
[9]	张子新, 孙钧. 块体理论赤平解析法及其在硐室稳定分析中的应用[J]. 岩石力学与工程学报, 2002, 21(12): 1756-1760. doi: 10.3321/j.issn:1000-6915.2002.12.002 ZHANG Zixin, SUN Jun. Stereoanalytic method for block theory and its application in stability analysis of a cave[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(12): 1756-1760. (in Chinese) doi: 10.3321/j.issn:1000-6915.2002.12.002
[10]	张奇华, 邬爱清, 石根华. 关键块体理论在百色水利枢纽地下厂房岩体稳定性分析中的应用[J]. 岩石力学与工程学报, 2004, 23(15): 2609-2614. doi: 10.3321/j.issn:1000-6915.2004.15.024 ZHANG Qihua, WU Aiqing, SHI Genhua. Application of key block theory to analysis of rock stability for underground plant in Baise hydraulic project[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(15): 2609-2614. (in Chinese) doi: 10.3321/j.issn:1000-6915.2004.15.024
[11]	杨勇, 张敏思, 王述红, 等. 钓鱼台隧道复杂围岩设计施工过程关键块确定技术[J]. 现代隧道技术, 2017, 54(3): 105-111. YANG Yong, ZHANG Minsi, WANG Shuhong, et al. Key block determination in the progress of design and construction of the Diaoyutai tunnel in complex surrounding rock[J]. Modern Tunnelling Technology, 2017, 54(3): 105-111. (in Chinese)
[12]	ELMOUTTIE M, POROPAT G, KRÄHENBÜHL G. Polyhedral modelling of underground excavations[J]. Computers and Geotechnics, 2010, 37(4): 529-535. doi: 10.1016/j.compgeo.2010.02.009
[13]	LI M C, ZHANG Y, ZHOU S B, et al. Refined modeling and identification of complex rock blocks and block-groups based on an enhanced DFN model[J]. Tunnelling and Underground Space Technology, 2017, 62: 23-34. doi: 10.1016/j.tust.2016.11.002
[14]	LIU J, LI Z K, ZHANG Z Y. Stability analysis of block in the surrounding rock mass of a large underground excavation[J]. Tunnelling and Underground Space Technology, 2004, 19(1): 35-44. doi: 10.1016/S0886-7798(03)00084-1