加高扩容尾矿库地震动力稳定性分析
Seismic Dynamic Stability Analysis of Enlarged Tailings Reservoir
摘要: 尾矿库在地震作用下,可能会出现坝体浸润线抬升、尾砂液化、坝体发生位移和变形,最后造成坝体失稳甚至酿成溃坝事故。要保证尾矿库在地震动力作用下安全稳定运行,就需要对尾矿库坝体进行地震动力稳定性分析。本文通过建立尾矿库数值模型,运用有限元法选取有代表性的剖面对尾矿库进行地震动力仿真模拟计算,根据计算结果,判断安全稳定性参数指标是否在可控范围内。得出结论:当尾矿库达到加高扩容终期标高时,尾矿库主坝在经历地震的过程中,坝体在地震动力作用下仍能满足相应等别安全滩长与安全超高的要求,均无明显成片液化区的产生,坝体永久变形的位置和位移值均处于安全范围,尾矿库抗震稳定安全系数满足工程对坝体抗滑稳定的要求。
Abstract: Under the seismic action of tailings reservoir, elevation of the seepage line of the dam body, liquefaction of the tailings, displacement and deformation of the dam may occur, and finally, the dam will become unstable or even break. In order to ensure the safe and stable operation of the tailings dam under the action of seismic dynamic, it is necessary to carry out the seismic dynamic stability analysis of the tailings dam. In this paper, based on the establishment of the numerical model of the tailing pond, the representative sections selected by the finite element method are used to carry out the simulation calculation of the seismic dynamics of the tailing pond, according to the results of the calculations, to judge whether the safety and stability parameter index is in the controllable range. It is concluded that the main dam of the tailing reservoir can still meet the requirements of equivalent safety beach length and super-high safety under the action of seismic in the process of experiencing seismic when the end of the heightening and expansion period of the tailings reservoir is reached, there is no obvious piece of liquefied area, the location and displacement of permanent deformation of the dam body are in the safe range, and the seismic stability factor of safety of the tailings reservoir can meet the requirements of the engineering for the anti-slide stability of the dam body.
文章引用:陈铁亮, 母传伟, 孙靖霄, 曹婷, 赵阳. 加高扩容尾矿库地震动力稳定性分析[J]. 矿山工程, 2025, 13(1): 87-101. https://doi.org/10.12677/me.2025.131012

1. 引言

尾矿库运行过程中,在地震作用下,可能会出现尾矿库坝体浸润线抬升、坝体尾砂液化、坝体发生位移和变形,最后造成坝体失稳甚至酿成溃坝事故[1] [2]。能否保证这种高势能危险源在地震动力作用下安全稳定运行,对尾矿库进行地震动力稳定性分析尤为重要[3]

目前国内对加高扩容尾矿库方面的研究主要侧重于在地震作用下的坝体位移、液化、渗流和安全监测,采用不同方法从不同的角度对加高扩容尾矿库进行了稳定性分析[4]-[7]。周薛淼应用Geo-studio软件对尾矿坝加高扩容前后的渗流场进行模拟分析,探讨了尾矿坝的静力稳定性和地震工况下的动力响应[8]。王峰等基于Geo-studio模拟尾矿坝加高扩容前后渗流场并优化渗透系数,分析了坝体静力稳定性和地震动力响应[9]。霍晨玮等基于饱和–非饱和渗流理论采用极限平衡法分析了持续降雨对某加高扩容尾矿坝渗流稳定性的影响[10]。吴顺川等结合理正岩土渗流分析软件对某尾矿库加高扩容工程的坝体加固设计与稳定性进行了分析[11]

借助midas软件运用有限元法建立尾矿库数值模型,选取典型剖面对尾矿库扩容工程进行地震动力稳定计算,程序自动搜索滑动面位置,并计算相应的稳定安全系数,进一步分析尾矿坝坝坡的抗震稳定性。这种方法同时考虑了在地震动力影响下的液化和渗流的作用;对数值模型实时更新,可以完成尾矿库全生命周期的动力稳定计算分析;输出彩色应力分布云图、坝体位移变形云图和液化区分布图,分析结果直观显现、一目了然。

2. 尾矿库概况和参数选取

2.1. 尾矿库基本情况

某尾矿库位于承德市,是一座生产运行中的库,因矿山扩大产能,现有库容已不能满足生产需要。尾矿库要在原设计基础上加高扩容,加高扩容后尾矿库堆积坝高度达199 m,达到了国家二等库的等别标准。

库区处于三面环山的山谷内,主坝体朝西,详见图1。尾矿库设计总库容1.45亿m3,尾矿库等别二等,初期坝为碾压碎石坝,堆积坝用上游法主坝方式。

尾矿库需要在原设计基础上加高扩容,加高扩容后总坝高增加37 m,主坝堆积标高从652 m升至689 m,尾矿库总库容量增加4335万m3,扩容后最终坝高199 m,库容量达到1.8亿m3

尾矿库加工扩容后,尾矿库坝体总高度和库容量也都大大增加,尾矿库安全风险也随之剧增,加高扩容后的尾矿库在地震工况下还是否稳定,必须要先有了可靠的结论才能实施扩容工程。需要对尾矿库加高扩容终期作坝体地震动力稳定性模拟分析,并验证坝体稳定安全系数是否满足规范要求。为完成此项任务,要建立尾矿库数值模型,选取典型剖面对尾矿库扩容工程进行地震动力稳定计算,同时还要考虑在地震动力作用下液化和渗流对坝坡的抗震稳定性影响。

Figure 1. Profile location map selected for the calculation of the total plane and stability of the tailings reservoir

1. 尾矿库总平面和稳定性计算所选取的剖面位置图

2.2. 库区地震特征参数

2.2.1. 抗震设防标准

库区抗震设防烈度6度,基本地震加速度值为0.05 g。

2.2.2. 输入地震波的选取

地震波选取的原则是使输入地震波的特性和尾矿库所在场地的条件相符合,主要包括地震烈度、场地条件和卓越周期等,同时还应使地震波满足地震动三要素——地震动峰值、地震动频谱特性和地震动持续时间。地震波目前主要分为天然波和人工拟合地震波,不同的场地类别都有实地监测到的著名的地震波,这些波在工程中应用较多,大量的工程经验表明,这种方法是合理的,能够满足工程所需的精度要求[12] [13]。尾矿库工程的动力反应,采用2条实测地震加速度记录及一条人工拟合地震加速度时程,其中两条实测地震波分别为TAR_TARZANA波和迁安波。TAR_TARZANA波、迁安波和人工波时程曲线如图2~4所示。其中,TAR_TARZANA波、迁安波和人工波总时长为30 s,峰值加速度为0.05 g,时程采样间隔均为0.01 s。

根据目标反应谱调整地震波的幅值和周期,将滤波及基线校正后的三条输入地震波反应谱与目标反应谱比较,调整后的地震波反应谱与目标反应谱在统计意义上相符,人工拟合地震波符合效果最好,可用于地震动力时程分析[14] [15]

Figure 2. TAR_TARZANA seismic wave time history curve

2. TAR_TARZANA地震波时程曲线

Figure 3. Qian’an wave time history curve

3. 迁安波时程曲线

Figure 4. Artificial fitting seismic wave time history curve

4. 人工拟合地震波时程曲线

2.3. 尾矿库材料参数

尾矿库材料参数如表1所示。

Table 1. Stability calculation parameter table

1. 稳定性计算参数表

岩土层编号

岩土层名称

天然重度

饱和重度

黏聚力

内摩擦角

γ

γsat

标准值

标准值

(kN/m3)

(kN/m3)

C (kPa)

Φ (˚)

1

碾压块石、碎石

22.00

22.50

1.0

36.0

2

碾压块石、碎石

21.50

22.00

0.5

35.0

1

尾中砂

19.30

21.60

水上:3.0

水上:30.0

水下:1.0

水下:28.0

2

尾中砂

19.80

21.20

水上:4.5

水上:33.0

水下:2.5

水下:31.0

1

尾粉砂

19.50

21.50

水上:5.0

水上:28.0

水下:2.0

水下:26.0

2

尾粉砂

20.00

21.20

水上:6.0

水上:31.0

水下:3.0

水下:29.0

尾粉质黏土

20.00

20.10

19.0

12.0

1

粉质黏土(第四系覆盖层)

19.50

19.80

26.0

18.0

2

粉土(第四系覆盖层)

19.00

19.50

12.0

20.0

含土角砾(第四系覆盖层)

19.30

19.60

1.0

33.0

1

强风化花岗岩(基岩)

22.30

22.30

60.0

35.0

2

中风化花岗岩(基岩)

24.20

24.20

400.0

37.0

2.4. 尾砂动参数

根据项目勘察报告并参考同类尾矿库的尾砂动力试验资料,尾矿库工程动力计算采用的尾砂动参数如表2表3所示。

Table 2. Dynamic strength parameter table

2. 动强度参数表

土样名称

干密度

固结比

10周

20周

30周

C d

φ d

C d

φ d

C d

φ d

g/cm3

Kc

kPa

kPa

kPa

尾中砂

1.68

1.0

0

10.0

0

9.8

0

9.6

1.5

0

20.3

0

19.7

0

19.4

2.0

0

26.6

0

26.3

0

26.0

尾粉砂

1.65

1.0

2

9.0

2

8.9

2

8.8

1.5

2

17.8

2

17.6

2

17.4

2.0

2

25.5

2

25.2

2

25.0

尾粉质黏土

1.54

1.0

4

9.4

4

9.2

4

9.0

1.5

4

13.4

4

13.2

4

13.0

2.0

4

19.4

4

19.2

4

19.0

注: E dmax :最大动弹性模量, G dmax :最大动剪切模量。

Table 3. Maximum modulus data table

3. 最大模量数据表

土类

名称

干密度

固结比

围压

a

E dmax

G dmax

g/cm3

kPa

MPa

MPa

尾中砂

1.68

1.0

100

0.0049

204

68

150

0.0040

251

84

200

0.0027

373

124

1.5

100

0.0036

276

92

150

0.0028

355

118

200

0.0026

384

128

2.0

100

0.0031

324

108

150

0.0027

376

125

200

0.0024

416

139

尾粉砂

1.65

1.0

100

0.0055

183

61

150

0.0042

236

79

200

0.0037

271

90

1.5

100

0.0050

200

67

150

0.0036

276

92

200

0.0034

298

99

2.0

100

0.0044

227

76

150

0.0038

265

88

200

0.0028

360

120

尾粉质黏土

1.54

1.0

100

0.0052

192

64

150

0.0044

228

76

200

0.0038

263

88

1.5

100

0.0053

190

63

150

0.0041

242

81

200

0.0035

288

96

2.0

100

0.0046

217

72

150

0.0040

251

84

200

0.0032

313

104

注: E dmax :最大动弹性模量, G dmax :最大动剪切模量。

3. 扩容终期主坝剖面动力计算

3.1. 地震结束时有效大小主应力分布

用midas软件建立尾矿库数值模型,运用有限元法计算出尾矿坝单元地震作用下的动应力,选取有代表性的剖面对尾矿库扩容工程在地震作用下进行仿真模拟计算,根据计算结果,判断安全稳定性参数指标是否在可控范围内。在动力计算中,假定滑动面形状,给定搜索范围,由程序自动寻找滑动面的位置,并计算相应的稳定安全系数,进一步分析尾矿坝坝坡的抗震稳定性。

图5~7为尾矿库扩容终期689 m标高主坝剖面模型在三种地震波作用结束后的有效大小主应力分布云图。TAR_TARZANA波作用下的有效大主应力最大值为3506 kPa,最小值为62.84 kPa;迁安波作用下的有效大主应力最大值为3491 kPa,最小值为63.9 kPa;人工波作用下的有效大主应力最大值为3573 kPa,最小值为61.9 kPa。在三种地震波作用下,尾矿库最大主应力大小均在合理范畴内,应力随尾矿埋深的增加逐渐增大,应力变化分布合理,尾矿库坝体没有应力集中现象,坝内主应力是压应力,不存在受拉区。

Figure 5. Effective major principal stress distribution of the main dam section at the end of the TAR wave expansion period

5. TAR波结束时扩容终期主坝剖面有效大主应力分布

Figure 6. Effective major principal stress distribution of the main dam section at the end of Qian’an wave expansion period

6. 迁安波结束时扩容终期主坝剖面有效大主应力分布

Figure 7. Effective major principal stress distribution of main dam section at the end of artificial wave expansion period

7. 人工波结束时扩容终期主坝剖面有效大主应力分布

3.2. 地震结束时峰值动剪应力分布

图8~10为尾矿库扩容终期主坝剖面在三种地震波作用结束后的峰值动剪应力分布云图,可以看出峰值动剪应力值从滩顶及坝体外坡面以下向库底基岩方向呈带状分布,应力集中的部位都在远离坝体的库中和库尾的下部,对尾矿库的安全性没有影响。TAR_TARZANA波作用下的峰值动剪应力最大值为79.1 kPa;迁安波作用下的峰值动剪应力最大值为81.8 kPa;人工波作用下的峰值动剪应力最大值为95.2 kPa。

Figure 8. Peak dynamic shear stress distribution of the main dam section at the end of the TAR wave expansion period

8. TAR波结束时扩容终期主坝剖面峰值动剪应力分布

Figure 9. Peak dynamic shear stress distribution of the main dam section at the end of Qian’an wave expansion period

9. 迁安波结束时扩容终期主坝剖面峰值动剪应力分布

Figure 10. Peak dynamic shear stress distribution of the main dam section at the end of artificial wave expansion period

10. 人工波结束时扩容终期主坝剖面峰值动剪应力分布

3.3. 地震结束时水平向位移分布

图11~13为尾矿库扩容终期主坝剖面在三种地震波作用结束时的水平向位移分布云图,可以看出三种地震波结束时位移最大值均分布在滩顶位置。TAR_TARZANA波作用下的水平向最大值为0.13 m;迁安波作用下的水平向位移最大值为0.14 m;人工波作用下的水平向位移最大值为0.16 m。

Figure 11. Horizontal displacement distribution of the main dam section at the end of the TAR wave expansion period

11. TAR波结束时扩容终期主坝剖面水平向位移分布

Figure 12. Horizontal displacement distribution of the main dam section at the end of Qian’an wave expansion period

12. 迁安波结束时扩容终期主坝剖面水平向位移分布

Figure 13. Horizontal displacement distribution of the main dam section at the end of the artificial wave expansion period

13. 人工波结束时扩容终期主坝剖面水平向位移分布

3.4. 堆积坝顶监测点加速度反应曲线

图14~16为尾矿库扩容终期主坝剖面在三种地震波作用下堆积坝顶监测点的加速度反应曲线。计算结果表明,加速度反应曲线和输入的加速度时程曲线吻合程度较好,仅在幅值上有所不同,表明了输入地震动引起的地震动过程不存在过大的放大效应或者不合理的阻尼效应。

Figure 14. Acceleration response curve of TAR wave accumulation dam crest at the end of expansion period

14. 扩容终期TAR波堆积坝顶加速度反应曲线

Figure 15. Acceleration response curve of Qian’an wave accumulation dam crest at the end of expansion period

15. 扩容终期迁安波堆积坝顶加速度反应曲线

Figure 16. Acceleration response curve of artificial wave accumulation dam crest at the end of expansion period

16. 扩容终期人工波堆积坝顶加速度反应曲线

3.5. 液化区分布

图17~19为三种地震波作用下尾矿库扩容终期689 m标高主坝剖面的液化结果,三图显示均无明显成片液化区。尾矿库浸润线埋深及干滩长度处于安全范围,坝体安全程度较高,在6度地震作用下各种计算工况下均无明显成片液化区的产生。

Figure 17. Liquefaction zone of main dam section under the action of TAR wave at the end of tailings reservoir expansion period

17. 尾矿库扩容终期主坝剖面TAR波作用下液化区

Figure 18. Liquefaction zone of the main dam section under the action of Qian’an wave at the end of tailings reservoir expansion period

18. 尾矿库扩容终期主坝剖面迁安波作用下液化区

Figure 19. Liquefaction zone of main dam section under the action of artificial wave at the end of tailings reservoir expansion period

19. 尾矿库扩容终期主坝剖面人工波作用下液化区

3.6. 永久变形分布

图20~22为三种地震波作用下尾矿库扩容终期689 m标高主坝剖面的水平方向永久变形云图。可以看出,尾矿库在经历地震的过程中,扩容终期689 m标高主坝剖面在TAR_TARZANA波、迁安波和人工波作用下,水平向永久变形值分别为0.22 m、0.23 m和0.20 m。

Figure 20. Horizontal permanent deformation cloud diagram under the action of TAR wave at the end of tailings reservoir expansion period

20. 尾矿库扩容终期TAR波作用下水平方向永久变形云图

Figure 21. Vertical permanent deformation cloud diagram under the action of TAR wave at the end of tailings reservoir expansion period

21. 尾矿库扩容终期TAR波作用下竖直方向永久变形云图

Figure 22. Horizontal permanent deformation cloud diagram under the action of Qian’an wave at the end of tailings reservoir expansion period

22. 尾矿库扩容终期迁安波作用下水平方向永久变形云图

4. 计算结果分析与讨论

4.1. 计算结果分析

对尾矿库扩容终期689 m标高主坝剖面模型进行震后动力稳定分析,计算条件为震后孔隙水压力及震后应力重分布的结果。

Figure 23. Safety factor curve of TAR wave dynamic time history method of main dam section at the end of expansion period

23. 扩容终期主坝剖面TAR波动力时程线法安全系数曲线

Figure 24. Safety factor curve of Qian’an wave dynamic time history method of main dam section at the end of expansion period

24. 扩容终期主坝剖面迁安波动力时程线法安全系数曲线

Figure 25. Safety factor curve of artificial wave dynamic time history method of main dam section at the end of expansion period

25. 扩容终期主坝剖面人工波动力时程线法安全系数曲线

Figure 26. Safety factor curve of Quasi-static simplified Bishop method of main dam section at the end of expansion period

26. 扩容终期主坝剖面拟静力简化毕肖普法安全系数

图23~25分别示三种剖面坝体动力时程线法安全系数曲线及拟静力法计算结果。图23显示扩容终期主坝剖面TAR波动力时程线法安全系数曲线,最小值为1.72,最大值为1.97;图24显示扩容终期主坝剖面迁安波动力时程线法安全系数曲线,最小值为1.79,最大值为1.95;图25显示扩容终期主坝剖面人工波动力时程线法安全系数曲线,最小值为1.76,最大值为1.97;图26显示扩容终期主坝剖面拟静力简化毕肖普法计算结果,结果显示安全系数为1.576。

计算结果显示的安全系数均满足工程对坝体抗滑稳定的要求。

4.2. 计算结果讨论

1) 运用有限元法计算尾矿坝单元地震作用下的动应力,采用三条特征地震波,其中两条为实测的TAR_TARZANA波和迁安波,另一条为根据目标反应谱人工拟合的地震波,作为坝体动力时程分析的输入地震波。选取代表性剖面对尾矿库扩容工程进行地震动力稳定计算,得出结果如下:

设计地震作用下计算剖面的坝内浸润线抬升幅度不明显,未出现超孔隙水压力集中区,超孔隙水压力最大值远小于有效应力计算值,计算标高情况下无明显成片液化区的产生。原因是随着坝体的升高,前期排放的尾砂固结程度较好,地震时坝体内部超孔隙水压力上升幅度有限,尾矿库整体动力安全度较大。由此认为地震作用对坝体的稳定影响较小,坝体在地震动力作用下仍能满足相应等别安全滩长与安全超高的要求,坝体稳定状态处于安全范围。

2) 采用Duncan的E-ν模型研究尾矿库工程运行过程中各阶段的应力应变状态、变形大小、应力水平以及坝体永久变形和坝坡的抗震稳定性,得出坝体地震残余变形和动力稳定性的有关规律和结论:

① 尾矿库工程在抗震设防烈度地震作用下,分别在初期坝顶和堆积坝顶选取监测点监测水平加速度反应,加速度反应曲线和输入的加速度时程曲线吻合程度较好,仅在幅值上有所不同,表明输入地震动引起的地震动过程不存在过大的放大效应或者不合理的阻尼效应。

② 尾矿坝的典型剖面在加高扩容终期标高内,其最大主应力及最小主应力大小均在合理范畴内,坝内主应力是压应力,不存在受拉区,具备静力稳定特征。

③ 尾矿坝的位移变形主要表现在垂直方向,服务期内尾矿坝的变形量并不大,总体上属于一个漫长平稳的沉降过程。

④ 尾矿库内剪应力近似呈层状分布,最大剪切应力主要分布在库底基岩层。库内尾砂绝大部分区域剪应力值较小,不会产生相对错动滑移。

⑤ 坝体在地震波作用下,水平向地震反应大于竖向地震反应,震后坝体水平位移大于竖向位移。终期坝体的最大水平向位移为0.16 m,坝体沉降位移量较小。尾矿坝整体剪应力水平较小,无明显剪切变形区域,不会发生剪切滑动现象。

5. 结论

通过建立尾矿库数值模型,运用有限元法选取代表性剖面对尾矿库扩容工程进行地震动力稳定计算,计算出相应的稳定安全系数,进一步分析尾矿坝坝坡的抗震稳定性,并得出如下结论:

1) 当尾矿库达到加高扩容终期标高时,坝体在地震动力作用下仍能满足相应等别安全滩长与安全超高的要求,库区范围内均无明显成片液化区的产生,库内尾砂固结程度较好,地震时坝体内部超孔隙水压力上升幅度有限,尾矿库整体动力安全度较高,坝体稳定状态处于安全范围。

2) 尾矿库主坝在经历地震的过程中,坝体永久变形的位置和位移值均处于安全范围,尾矿库抗震稳定安全系数满足工程对坝体抗滑稳定的要求。

借助midas软件运用有限元法建立数值模型实现了对尾矿库在地震工况下稳定分析计算,不同于其他手段的坝体位移分析、液化分析、静力稳定分析和安全监测,也不同于理正岩土软件的渗流分析。这种方法同时考虑了在地震动力影响下的液化和渗流的作用;尾矿库的运行是动态过程,库容和坝高每天都在变化,借助软件可以对数值模型实时更新,实现尾矿库全生命周期的稳定计算模拟;输出彩色应力分布云图、坝体位移变形云图和液化区分布图,分析结果直观显现。

NOTES

*第一作者。

#通讯作者。

参考文献

[1] 刘佳浩, 刘红岩, 邹宗山. 竖向及水平地震作用下尾矿库动力响应及液化稳定性分析[J]. 金属矿山, 2023(11): 98-100.
[2] 张树茂. 西北地区某尾矿库地震动力时程响应分析[J]. 中国矿业, 2022, 31(S1): 114-118.
[3] 李全明, 杨曌, 张红. 基于QUAKE/W模型的尾矿库地震动力响应研究[J]. 中国矿业, 2019, 28(4): 163-167.
[4] 郑彬彬. 高浓度尾矿上游式堆坝基础性问题研究及坝体稳定性分析[D]: [博士学位论文]. 重庆: 重庆大学, 2017.
[5] 肖文杰. 西马架子尾矿坝稳定性研究[D]: [硕士学位论文]. 长春: 吉林大学, 2023.
[6] 郄永波, 周汉民, 崔旋. 强震条件下陡底坡尾矿坝动力稳定性分析[J]. 中国矿业, 2018, 27(z1): 390-393, 399.
[7] 张水兵. 高堆磷石膏尾矿坝扩容稳定性分析及加固措施研究[D]: [硕士学位论文]. 昆明: 云南大学, 2018.
[8] 周薛淼. 三元洞尾矿库加高扩容后坝体的稳定性研究[J]. 山西建筑, 2024, 50(7): 175-181.
[9] 王峰, 杨勇, 王宇驰, 等. 基于Geo-studio的尾矿库加高扩容后坝体稳定性分析[J]. 有色金属(矿山部分), 2022, 74(5): 90-98.
[10] 霍晨玮, 沈振中. 降雨持续时间对某加高扩容尾矿坝渗流稳定性影响分析[J]. 有色金属(矿山部分), 2021, 73(2): 16-20.
[11] 吴顺川, 赵志强, 张小强, 等. 某尾矿库加高扩容工程的坝体加固设计与稳定性分析[J]. 昆明理工大学学报(自然科学版), 2021, 46(1): 36-44.
[12] 高德志, 毕鹏, 刘艳. 各国规范关于时程分析中地震波选取方法的对比[J]. 建筑结构, 2022, 52(16): 67-73.
[13] 冯伟栋, 曹均锋, 彭刘亚, 等. 地震安评中天然地震波选取方法研究[J]. 华南地震, 2021, 41(4): 97-101.
[14] 王文松, 尹光志, 魏作安, 等. 基于时程分析法的尾矿坝动力稳定性研究[J]. 中国矿业大学学报, 2018, 47(2): 271-101.
[15] 尹光志, 王文松, 魏作安, 等. 尾矿库加高扩容坝体动力反应与抗震性能分析[J]. 岩石力学与工程学报, 2018(S1): 3132-3142.