1. 引言

为了确保尾矿库的安全性，对尾矿库坝体要进行地震动力稳定计算，而在此之前要先对坝体进行静力分析[1]。静力状态是动力反应分析的基础，其初始静应力是影响动力状态饱和非粘性土液化的重要因素之一[2] [3]。坝体计算静力分析的目的，就是确定震前坝体与坝基中的应力状态，判断深层稳定性并为地震动力反应分析和液化分析提供基础[4] [5]。

尾砂的本构模型是表征尾砂料力学特性的基本关系，是用有限元解决尾矿坝内的应力应变及强度——变形稳定问题必不可少的关系[6]。只有选择合适的尾矿坝本构模型及其参数，才能通过尾矿坝的静动力反应分析，合理真实地给出尾矿坝的动应力、变形和孔隙水压力的发展过程和空间分布特征，对尾矿坝的抗震性能及地震稳定性做出科学的评价[7]。

尾砂料在动荷载作用下不仅有弹塑性的特点，而且还有粘性的特点，可将坝料视为弹性、塑性和粘滞性的粘弹塑性体[8]。此外，尾砂具有明显的各向异性，加上水的影响，使尾砂的动应力应变关系十分复杂。因此描述尾砂料的动应力应变关系，必须对尾砂的非线性、滞后性和变形积累性的特性均有较深入的研究[9]。

在尾矿坝、土石坝工程中，进行静力计算分析使用最多、应用最广的是邓肯$E\text{-}\nu$ 非线性弹性模型[10] [11]。该模型特点：(1) 简单清晰，物理意义明确，邓肯-E模型基于双曲线假设，能够通过常规三轴试验曲线直接获取参数。(2) 适用范围广，该模型适合描述低剪应力水平下的土体应力应变关系，适用于尾矿坝等工程中常见的应力状态。(3) 非线性特征描述清晰，能够较好地反映土体在不同应力水平下的变形特性，包括剪胀和剪缩现象。

邓肯-E模型与南水模型相比，后者更适用高土石坝复杂结构，但在尾矿坝工程中，前者更适合。邓肯-E模型与线性弹性模型相比，后者不能反映土体非线性变形，而前者通邓肯-E模型过双曲线关系能够更准确地反映应力应变关系。

邓肯-E模型与沈珠江模型相比：后者更适合计算土体的残余变形，但在初始应力状态计算中，邓肯模型适用性更高。

邓肯-E模型适用性分析：(1) 参数获取方便，实际适用性强，计算过程相对简单，计算参数可通过常规试验获得。(2) 计算结果可靠性高，与试验数据拟合性高。(3) 适应多种工况，适用性强，邓肯模型不仅适用于静力计算，还可通过改进用于动力分析。

综上所述，邓肯-E模型在尾矿坝和土石坝的静力计算中具有较高的适用性和准确性。通过与其他模型的比较和适用性分析，可以证明其在尾矿坝初始应力状态计算中的优势。

因此，在进行尾矿坝初始应力状态计算时，采用邓肯$E\text{-}\nu$ 模型计算尾矿静力状态下的应力应变，特别是剪应力状态。

2. 非线性静力模型

尾砂的基本变形特性表现为应力应变关系的非线性，邓肯$E\text{-}\nu$ 模型常用来描述这种非线性[10] [11]。尾砂料的弹性模量$E$ ，泊松比$\nu$ 是应力(或应变)的函数，根据弹塑性力学的广义虎克定律：

$\left\{\sigma \right\}=\left[D\left(\sigma \right)\right]\left\{\epsilon \right\}$ (1)

式中 $\left[D\left(\sigma \right)\right]$ ——刚度矩阵。

Kondner根据大量土的三轴试验的应力应变关系曲线，提出可以用双曲线函数拟合三轴试验的$\left({\sigma }_1-{\sigma }_3\right)~{\epsilon }_a$ 关系曲线，即：

${\sigma }_1-{\sigma }_3=\frac{{\epsilon }_a}{a+b{\epsilon }_a}$ (2)

式中 $a$ 、$b$ ——试验常数。

对于常规三轴压缩试验，${\epsilon }_a$ 等于主应变${\epsilon }_1$ ，邓肯等人根据应力应变的双曲线关系提出一种增量弹性模型。

在常规三轴压缩试验中，由于$d{\sigma }_2=d{\sigma }_3=0$ ，所以切线弹性模量：

$E_t=\frac{d\left({\sigma }_1-{\sigma }_3\right)}{d{\epsilon }_1}=\frac{a}{{\left(a+b{\epsilon }_1\right)}^2}$ (3)

试验的起始点时，${\epsilon }_1=0$ ，$E_t=E_i$ ，则：

$E_i=\frac{1}{a}$ (4)

式(2)中，令${\epsilon }_1\to \infty$ ，则：

${\left({\sigma }_1-{\sigma }_3\right)}_{ult}=\frac{1}{b}$ (5)

即参数$a$ 等于试验中初始变形模量$E_i$ 的倒数，参数$b$ 等于双曲线渐近线对应的极限偏差应力${\left({\sigma }_1-{\sigma }_3\right)}_{ult}$ 的倒数。

一般在土的试样中，若应力应变曲线接近于双曲线关系，则根据一定的应变值(如${\epsilon }_1=15\%$ )来确定土的强度${\left({\sigma }_1-{\sigma }_3\right)}_f$ ，因此定义破坏比$R_f$ 为：

$R_f=\frac{{\left({\sigma }_1-{\sigma }_3\right)}_f}{{\left({\sigma }_1-{\sigma }_3\right)}_{ult}}$ (6)

经推导可以得出：

$E_t=E_i{\left[1-R_f\frac{{\sigma }_1-{\sigma }_3}{{\left({\sigma }_1-{\sigma }_3\right)}_f}\right]}^2$ (7)

根据摩尔——库伦(Mohr-Coulomb)强度准则：

${\left({\sigma }_1-{\sigma }_3\right)}_f=\frac{2c\cos \phi +2{\sigma }_3\sin \phi }{1-\sin \phi }$ (8)

式中 $c$ 、$\phi$ ——材料的内聚力和内摩擦角。

Janbu发现三轴试验的初始变形模量$E_i$ 与围压${\sigma }_3$ 有关，$\lg \left(E_i/Pa\right)$ 与$\lg \left({\sigma }_3/Pa\right)$ 近似成线性关系。则：

$E_i=KPa{\left(\frac{{\sigma }_3}{Pa}\right)}^n$ (9)

其中$Pa$ 为大气压，量纲同${\sigma }_3$ ；$K$ 、$n$ 为试验常数，分别为$\lg \left(E_i/Pa\right)$ 与$\lg \left({\sigma }_3/Pa\right)$ 函数直线的截距和斜率。将式(8)和式(9)代入式(7)中得到切线变形模量：

$E_t=KPa{\left(\frac{{\sigma }_3}{Pa}\right)}^n{\left[1-\frac{R_f\left({\sigma }_1-{\sigma }_3\right)\left(1-\sin \phi \right)}{2c\cos \phi +2{\sigma }_3\sin \phi }\right]}^2$ (10)

同时邓肯等人根据实验资料，假定常规三轴压缩试验中的轴向应变${\epsilon }_1$ 与侧向应变$-{\epsilon }_3$ 之间也成双曲线关系，经推导可得切线泊松比：

${\nu }_t=\frac{G-F\lg \left({\sigma }_3/Pa\right)}{{\left\{1-\frac{D\left({\sigma }_1-{\sigma }_3\right)}{KPa{\left(\frac{{\sigma }_3}{Pa}\right)}^n\left[1-\frac{R_f\left({\sigma }_1-{\sigma }_3\right)\left(1-\sin \phi \right)}{2c\cos \phi +2{\sigma }_3\sin \phi }\right]}\right\}}^2}$ (11)

式中$D$ 为${\epsilon }_1~-{\epsilon }_3$ 关系渐近线的倒数，$G$ 、$F$ 为试验常数。根据弹性理论，$0<{\nu }_t<0.5$ 。

所以邓肯$E\text{-}\nu$ 模型中共包括$K$ 、$n$ 、$\phi$ 、$c$ 、$R_f$ 、$D$ 、$G$ 、$F$ 八个材料常数。

邓肯$E\text{-}\nu$ 模型通过常规三轴压缩试验的卸载——再加载曲线确定其卸载模量。由于此过程中应力应变表现为滞回圈曲线，所以用平均斜率表示$E_{ur}$ ：

$E_{ur}=K_{ur}Pa{\left(\frac{{\sigma }_3}{Pa}\right)}^n$ (12)

其中$K_{ur}$ 、$n$ 分别为$\lg \left(E_{ur}/Pa\right)~\lg \left({\sigma }_3/Pa\right)$ 直线的截距和斜率。

1984年邓肯等人提出加载函数(13)来作为加载卸载的判断准则：

$S_S=S\sqrt[4]{{\sigma }_3/Pa}$ (13)

其中$S$ 为应力水平：

$S=\frac{{\sigma }_1-{\sigma }_3}{{\left({\sigma }_1-{\sigma }_3\right)}_f}$ (14)

如果加载过程中加载函数的最大值为$S_{sm}$ ，那么临界应力水平为：

$S_m=\frac{S_{sm}}{\sqrt[4]{{\sigma }_3/Pa}}$ (15)

当$S>S_m$ 时，则为加载，切线变形模量按(3)式计算；

当$S<0.75S_m$ 时，则为卸载或再加载，切线弹性模量按式(12)计算；

当$0.75S_m<S<S_m$ 时，则用$E_t$ 与$E_{ur}$ 线性插值计算。

3. 静力计算参数及模型

坝体的静力计算参数根据相关勘察资料，结合类似工程实践共同选取，静力分析参数如表1和表2所示。

Table 1. Material parameter table for dam body static analysis

表1. 坝体静力分析材料参数表

Table 2. Initial mechanical index of dam materials

表2. 初期坝材料力学指标

静动力稳定性计算所选取的剖面位置如图1所示。

Figure 1. Selected section position for static and dynamic stability calculation of tailing dam

图1. 尾矿库坝动力稳定性计算所选取的剖面位置

图2~4为尾矿库原设计终期652 m标高、扩容终期689 m标高主坝和扩容终期689 m标高副坝模型图。

Figure 2. Model of 652 m elevation in the final design stage of tailing reservoir

图2. 尾矿库原设计终期652 m标高模型图

Figure 3. Model of main dam with 689 m elevation at the end of expansion of tailing reservoir

图3. 尾矿库扩容终期689 m标高主坝模型图

Figure 4. Model of auxiliary dam with 689 m elevation at the end of tailing reservoir expansion

图4. 尾矿库扩容终期689 m标高副坝模型图

4. 静力计算结果与分析

选取尾矿库原设计终期652 m标高、扩容终期689 m标高主坝模型和扩容终期689 m标高副坝模型进行有限元分析。

4.1. 原设计终期652 m标高静力结果

(1) 有效主应力

图5、图6分别代表尾矿库原设计终期652 m标高下尾矿坝内的有效大主应力${{\sigma }^{\prime }}_1$ 及有效小主应力${{\sigma }^{\prime }}_3$ 等值线云图，计算结果显示，坝内有效大主应力${{\sigma }^{\prime }}_1$ 最大值为3230 kPa，有效小主应力${{\sigma }^{\prime }}_3$ 最大值为1178 kPa。

Figure 5. Effective principal stress isoline at 652 m elevation at the end of original design

图5. 原设计终期652 m标高有效大主应力等值线图

Figure 6. Effective small principal stress isoline at 652 m elevation at the end of original design

图6. 原设计终期652 m标高有效小主应力等值线图

(2) 位移

图7和图8分别代表尾矿库原设计终期652 m标高下尾矿坝内水平位移、垂直位移分布状况，考虑尾砂在库内堆积的长期固结沉降过程，其中水平方向最大固结变形位于堆积坝顶，总体向初期坝方向位移为0.27 m；最大固结沉降亦位于堆积坝顶，最大位移0.86 m。

Figure 7. Horizontal displacement isoline of 652 m elevation at the end of original design

图7. 原设计终期652 m标高水平位移等值线图

Figure 8. 652 m elevation vertical displacement isoline in the end of original design

图8. 原设计终期652 m标高垂直位移等值线图

(3) 剪应力

图9、图10分别代表尾矿库原设计终期652 m标高下尾矿坝内的最大剪应力${\tau }_{\max }$ 及XY剪应力${\tau }_{xy}$ 分布状况，坝内最大剪应力${\tau }_{\max }$ 为1036 kPa，XY剪应力${\tau }_{xy}$ 最大值为461.5 kPa。剪应力较大值均分布于底部基岩处，堆积尾砂区域剪应力最小，无剪应力突变现象，尾矿坝是静力稳定的。

Figure 9. Maximum shear stress isoline at 652 m elevation at the end of original design of tailing reservoir

图9. 尾矿库原设计终期652 m标高最大剪应力等值线图

Figure 10. XY shear stress isoline at 652 m elevation at the end of original design of tailing reservoir

图10. 尾矿库原设计终期652 m标高XY剪应力等值线图

4.2. 扩容终期689 m标高主坝静力结果

(1) 有效主应力

图11、图12分别代表尾矿库扩容终期689 m标高主坝剖面尾矿坝内的有效大主应力${{\sigma }^{\prime }}_1$ 及有效小主应力${{\sigma }^{\prime }}_3$ 等值线云图，计算结果显示，坝内有效大主应力${{\sigma }^{\prime }}_1$ 最大值为3505 kPa，有效小主应力${{\sigma }^{\prime }}_3$ 最大值为1253 kPa。

Figure 11. Effective principal stress isoline of 689 m elevation main dam profile at the end of capacity expansion

图11. 扩容终期689 m标高主坝剖面有效大主应力等值线图

Figure 12. Effective small principal stress isoline of main dam section at 689 m elevation at the end of capacity expansion

图12. 扩容终期689 m标高主坝剖面有效小主应力等值线图

(2) 位移

图13和图14分别代表尾矿库扩容终期689 m标高主坝剖面水平位移、垂直位移分布状况，考虑尾砂在库内堆积的长期固结沉降过程，其中水平方向最大固结变形位于堆积坝顶，总体向初期坝方向位移为0.25 m；最大固结沉降亦位于堆积坝顶，最大位移1.04 m。

Figure 13. Horizontal displacement isoline of main dam section at 689 m elevation at the end of capacity expansion

图13. 扩容终期689 m标高主坝剖面水平位移等值线图

Figure 14. Vertical displacement isoline of main dam section at 689 m elevation at the end of capacity expansion

图14. 扩容终期689 m标高主坝剖面垂直位移等值线图

(3) 剪应力

图15、图16分别代表尾矿库扩容终期689 m标高主坝剖面坝内的最大剪应力${\tau }_{\max }$ 及XY剪应力${\tau }_{xy}$ 分布状况，坝内最大剪应力${\tau }_{\max }$ 为1139 kPa，XY剪应力${\tau }_{xy}$ 最大值为482.4 kPa。剪应力较大值均分布于底部基岩处，堆积尾砂区域剪应力最小，无剪应力突变现象，尾矿坝是静力稳定的。

Figure 15. Maximum shear stress isoline of main dam section at 689 m elevation at the end of expansion

图15. 扩容终期689 m标高主坝剖面最大剪应力等值线图

Figure 16. XY shear stress isoline of main dam section at 689 m elevation at the end of capacity expansion

图16. 扩容终期689 m标高主坝剖面XY剪应力等值线图

4.3. 扩容终期689 m标高副坝静力结果

(1) 有效主应力

图17、图18分别代表尾矿库扩容终期689 m标高副坝剖面尾矿坝内的有效大主应力${{\sigma }^{\prime }}_1$ 及有效小主应力${{\sigma }^{\prime }}_3$ 等值线云图，计算结果显示，坝内有效大主应力${{\sigma }^{\prime }}_1$ 最大值为1520.7 kPa，有效小主应力${{\sigma }^{\prime }}_3$ 最大值为548.6 kPa。

Figure 17. Effective major principal stress isoline of auxiliary dam profile at 689m elevation at the end of capacity expansion

图17. 扩容终期689 m标高副坝剖面有效大主应力等值线图

Figure 18. Effective small principal stress isoline of auxiliary dam profile at 689 m elevation at the end of capacity expansion

图18. 扩容终期689 m标高副坝剖面有效小主应力等值线图

(2) 位移

图19和图20分别代表尾矿库扩容终期689 m标高副坝剖面水平位移、垂直位移分布状况，考虑尾砂在库内堆积的长期固结沉降过程，其中水平方向最大固结变形位于堆积坝顶，总体向初期坝方向位移为0.07 m；最大固结沉降亦位于堆积坝顶，最大位移0.51 m。

Figure 19. Horizontal displacement isoline of auxiliary dam section at 689 m elevation at the end of capacity expansion

图19. 扩容终期689 m标高副坝剖面水平位移等值线图

Figure 20. Isogram of vertical displacement of auxiliary dam at 689 m elevation at the end of capacity expansion

图20. 扩容终期689 m标高副坝剖面垂直位移等值线图

(3) 剪应力

图21、图22分别代表尾矿库扩容终期689 m标高副坝剖面坝内的最大剪应力${\tau }_{\max }$ 及XY剪应力${\tau }_{xy}$ 分布状况，坝内最大剪应力${\tau }_{\max }$ 为507.5 kPa，XY剪应力${\tau }_{xy}$ 最大值为995.4 kPa。剪应力较大值均分布于底部基岩处，堆积尾砂区域剪应力最小，无剪应力突变现象，尾矿坝是静力稳定的。

Figure 21. Maximum shear stress isoline of auxiliary dam section at 689 m elevation at the end of expansion

图21. 扩容终期689 m标高副坝剖面最大剪应力等值线图

Figure 22. XY shear stress isoline of auxiliary dam section at 689 m elevation at the end of capacity expansion

图22. 扩容终期689 m标高副坝剖面XY剪应力等值线图

5. 分析和结论

5.1. 结果分析

(1) 尾矿坝在各个计算标高内，其有效大主应力、有效小主应力的分布及大小均在正常范围内，全部为压应力，坝体不存在拉应力区。这种结果的形成与以下几个方面有关。① 尾矿坝的筑坝材料具有良好的抗压性能。尾矿颗粒之间通过物理咬合方式相互支撑，使得坝体在正常工况下主要承受压应力。② 合理的颗粒级配和较高的密实度。尾矿颗粒大小分布均匀，且经过良好的压实处理，坝体内部的孔隙率较低。良好的颗粒级配和密实度能够使坝体内部的应力分布较为均匀，不能产生拉应力的条件，主要以压应力存在。③ 良好的排水系统能够有效降低坝体内部的孔隙水压力，坝体内部的孔隙水能够顺利排出，降低了坝体内部孔隙水压力。坝体在承受荷载时，主要是尾矿颗粒之间的有效应力在起作用，而且这种有效应力主要是压应力。因为排水系统降低了孔隙水压力对坝体内部应力分布的干扰，使得坝体内部应力分布更加符合材料的抗压特性。④ 在正常运行条件下，尾矿坝所受到的荷载主要是尾矿自身重力、水压力等。这些荷载在正常范围内时，坝体的应力响应也在正常范围内。

(2) 尾矿坝的位移变形主要表现在垂直方向，即沉积尾砂的沉降，服务期内尾矿坝的变形量并不大，总体上属一个漫长变形沉降过程。尾砂在库内堆积的长期固结沉降过程中，水平方向最大固结变形和最大固结沉降均位于堆积坝顶，尾矿堆积坝顶部位的固结沉降和水平变形最大，是由于自重应力集中、尾矿颗粒特性、固结过程的不均匀性、渗流压力以及放矿方式等多方面因素共同作用的结果。

(3) 坝内最大剪切应力主要分布在初期坝底及库底基岩层，是由于坝体自重和尾矿压力的传递、基岩的力学特性、渗流作用、坝体结构与施工过程以及时间效应等多方面因素共同作用的结果。这些因素导致坝底和基岩接触面处的应力集中，从而产生了较大的剪切应力。尾矿库内尾砂区域绝大部分剪应力值较小且不会产生相对错动滑移，主要是由于尾砂的物理力学特性(如颗粒特性、低压缩性、高内摩擦角)、堆积方式(均匀堆积、应力传递路径)、固结过程(固结沉降、孔隙压力均匀分布)、尾矿库的运行与管理(放矿管理、渗流控制)以及时间效应(长期稳定性、次固结变形)等多方面因素共同作用的结果。这些因素共同维持了尾砂区域的稳定性，使得剪应力值保持在较低水平，且不易发生相对错动滑移。

5.2. 结论

通过静力有限元计算分析，得出结论如下：

(1) 尾矿坝在各个计算标高内，其有效大主应力、有效小主应力的分布及大小均在正常范围内，全部为压应力，坝体不存在拉应力区。

(2) 尾矿坝的位移变形主要表现在垂直方向，即沉积尾砂的沉降，服务期内尾矿坝的变形量并不大，总体上属于一个漫长变形沉降过程。

(3) 坝内最大剪切应力主要分布在初期坝底及库底基岩层。库内尾砂区域绝大部分剪应力值较小，不会产生相对错动滑移。

综上分析，静力状态下不存在深部滑动及较大变形区域，坝体是静力稳定的。

NOTES

^*第一作者。

^#通讯作者。

参考文献