岩石统计损伤软化模型及其参数反演

，，， ， ，(.成都理工大学 地质灾害防治与地质环境保护国家重点实验室，成都 60059；.西南交通大学 地球科学与环境工程学院，成都 6003)

1 研究背景

岩石作为一种非均质的工程材料，其本构关系一直是各工程领域十分关注的内容[1-2]。自从Cook[3]首次提出岩石应力-应变全过程的概念，国内外学者对岩石本构关系进行广泛研究，提出大量岩石本构模型，但绝大多数模型难以较好地反映复杂应力状态下的岩石特性。

由于工程岩体普遍带有节理、裂隙等天然缺陷，岩体损伤理论将岩体内部的天然缺陷理想地简化为岩体的损伤，这与岩体特性十分吻合[4-5]。Lemaitre[6-7]提出应变等价性假说，Krajcinovic等[8-9]将连续损伤理论和统计强度理论有机地结合起来，提出了岩石统计损伤的概念。

目前关于岩石统计损伤本构模型的研究已取得大量成果，唐春安[10]和谢和平[11]分别从岩石微元强度以及损伤耗散能释放率分布的随机性出发进行的研究，推动了岩石统计损伤本构的研究。曹文贵等[12-13]、徐卫亚等[14]、袁小平等[15]、蒋维等[16]、Deng等[17]基于Lemaitre应变等价性假说对岩石统计损伤本构模型进行了较为全面的研究。已有相关研究中，岩石微元强度的概率分布有幂函数分布[18]、Weibull分布[12]、对数正态分布[16]和正态分布[13]等，而微元破坏准则主要利用了Mohr-Coulomb(M-C)准则和Drucker-Prager(D-P)准则。

岩石统计损伤本构模型的关键在于描述微元强度和选择其概率密度分布以及确定模型参数。已有研究确定模型参数多是采用作图法[10]或基于线性方程系数的回归求解的曲线拟合法[12,14]，但这2种方法具有一定主观性。且由于岩石材料的复杂性，所建模型未必适用于不同应力状态下的不同岩石材料，故笔者从岩石微元强度服从Weibull概率分布出发，采取参数反演分析的方法来确定模型参数，选取不同强度准则，以期建立更加符合实际的岩石统计损伤软化本构模型。

2 岩石统计损伤软化模型

2.1 损伤变量及损伤演化过程

假设岩石由n(n→∞)个微单元组成，微单元仅由k个破坏单元和(n-k)个无损单元组成。假设破坏单元无法承受荷载，无损单元在一定条件下可以转化为破坏单元，破坏单元和无损单元面积都为S，则损伤变量D可定义为

(1)

式中0≤k≤n，故D∈[0，1]。

当微元的应力水平Y达到微元强度F时，微元即刻破坏。假设微元强度F服从Weibull随机分布，则其概率密度函数P(F)为

(2)

式中m和F0为Weibull分布参数。

当加载到某应力水平Y时，已破坏的微元数目k为

(3)

联立式(1)和式(3)可得统计损伤演化方程

D=P(Y) 。

(4)

依据Lemaitre应变等价性假说[6-7]，可建立损伤本构方程为

σ*=σ/I-D=Cε/I-D。

(5)

式中：σ*为有效应力矩阵;σ为应力矩阵;C为弹性矩阵;ε为轴向应变;I为单位矩阵。

将式(3)和式(4)代入式(2)可得损伤变量D为

(6)

由式(6)可知，损伤变量D与岩石微元强度F直接相关，结合文献[12]，将微元强度F表示为

(7)

式中：α为与岩石性质相关的参数;I1为应力张量第一不变量;J2为应力偏量第二不变量。且有:

(8)

(9)

在岩石常规三轴试验中，围压σ3=σ2，假设岩石无损单元为弹性体，满足Hooke定律，无损单元向破坏单元转化瞬间完成，无时间效应，则有

(10)

式中：ε1为轴向应变；μ为泊松比；E为弹性模量。

通过式(5)、式(8)—式(10)可以得到：

(11)

(12)

2.2 统计损伤软化本构模型的建立

联立式(5)、式(6)、式(10)可得

(13)

式(13)即为本文所建立的岩石统计损伤软化本构模型。

3 模型参数的确定

将式(7)、式(11)、式(12)代入式(13)并求偏微分可得

(14)

设岩石全应力-应变曲线的峰值点为b，并引入峰值条件，即

(15)

将式(15)代入式(14)整理可得

式中：σ1b为峰值应力；ε1b为峰值应力对应的应变；Fb为峰值点对应的F值。

再将式(13)变形为

(17)

联立式(16)和式(17)可得：

(18)

(19)

由式(18)可看出，m的取值通过E，μ，σ3，σ1b，ε1b即可确定，与微元强度F无关。由式(7)和式(19)可看出，F0的取值取决于α，而式(7)中，不同的α对应不同的破坏准则。图1为不同α值所对应的岩石应力-应变曲线。

图1 不同α值对应的应力-应变曲线Fig.1 Stress-strain curves in the presence of different α values

由图1可看出，不同α值的峰后应力-应变曲线形态有较为明显的差异，随着α值的增大，峰后曲线由平坦逐渐变得陡直，峰后的应力衰减速度逐渐加剧，岩石更趋于表现出脆性特征。

不同强度准则的α值都不相同，现列出几种常用的强度准则[18](如表1)。

表1 不同强度准则下的α值Table 1 Values of α under different strength criteria

注：ψ为内摩擦角

图3 试验曲线与理论曲线的对比Fig.3 Comparison between test curve and theoretical curves

常用的模型参数确定方法多是采用作图法[10]或基于线性方程系数的回归求解的曲线拟合法[12,14]。由于岩石材料的复杂性，事先选定好α值所建模型未必适用于不同应力状态下的不同岩石材料，笔者认为这2种方法具有一定的主观性，故采用优化反演分析[19]方法对试验数据进行参数反演,确定参数α值。

在对弹塑性问题建立的优化反演分析法中，可在各参数可能的变化范围内找到一组使误差为最小的最佳参数，目标函数建立为

(20)

式中：N为某种工况下应变测试点总数;σi为应力理论计算值;σi*为应力测试值。

对参数α进行约束，文献[4]对此进行了具体论述，约束条件为

(21)

反演过程中α在约束范围内不断变换直至J(α)取得最小值。当α确定后，通过已解出的m再由式(7)、式(11)、式(12)、式(19)即可求得F0，最终代入式(13)即可得到本文统计损伤软化模型。

4 模型验证

4.1 常规三轴压缩试验

为了验证本文模型，采用YSJ-01-00岩石三轴流变试验机(如图2)对红层泥岩进行常规三轴压缩试验。岩样取自川西盆地某红层边坡[20-21]，密封后运回实验室加工成Φ50 mm×100 mm的圆柱样;试验围压设置为1，2，3 MPa，试样过程以0.5 MPa/s的速率施加围压，围压稳定后以0.1 mm/s的轴向位移速率施加轴向荷载。

图2 试验仪器Fig.2 Test apparatus

4.2 模型验证及参数求解

由本文试验得，所取红层泥岩的E=10.24 MPa，μ=0.22，内摩擦角ψ=28.72°。通过式(20)和式(21)确定α的反演值，将该反演值和利用表1中不同强度准则得到的α值进行理论曲线与试验曲线数据的对比。由于利用D-P2与D-P3准则确定的α值十分接近，故在图3中略去D-P2准则。基于反演α值通过式(18)和式(19)可得到统计损伤软化模型参数,如表2所示。

表2 统计损伤软化模型参数Table 2 Parameters of statistical damage model

表1归纳了不同强度准则下的α值，图3是基于表1的不同α值和反演α值拟合得到的应力-应变曲线。由图3可知，红层泥岩在σ3=1 MPa和2 MPa的情况下利用D-P3准则确定的α值更接近于反演结果，红层泥岩在σ3=3 MPa的情况下利用D-P1准则确定的α值更接近于反演结果。由此可知，α在不同的围压下取值有所差异，将α提前取为一个定值有待商榷,对试验曲线进行理论辨识时，不建议提前选取强度准则，通过反演分析的方法可有效提高辨识精度。

5 结 论

本文基于岩石微元强度服从Weibull随机分布，建立统计损伤软化模型，其中参数m可通过应力-应变曲线的峰值点确定，参数F0取决于α，而α通过反演分析得到，并根据反演结果选取适用于不同围压条件的强度准则。本文通过开展常规三轴压缩试验对此进行验证，证明了本文所建模型和通过反演分析确定参数的方法合理可行。

本文所建模型对于岩石应力-应变曲线有较强的辨识能力，模型参数较少，应用方便。模型中的微元强度参数α并非一个定值，α不仅可以反映岩石的不同应力状态，而且对于辨识峰后曲线形态起到了关键性作用。

参考文献：

[1] 黄海峰, 巨能攀, 李 路, 等. 基于改进Harris函数的岩石统计损伤软化模型[J].工程地质学报,2018,26(2):510-517.

[2] WANG Z L, LI Y C, WANG J G. A Damage-softening Statistical Constitutive Model Considering Rock Residual Strength[J]. Computers & Geosciences, 2007, 33(1):1-9.

[3] COOK N G W. The Failure of Rock[J]. International Journal of Rock Mechanics & Mining Sciences & Geomechanics Abstracts, 1965, 2(4):389-403.

[4] 郑颖人, 沈珠江, 龚晓南. 岩土塑性力学原理[M]. 北京: 中国建筑工业出版社, 2002.

[5] 温 韬, 唐辉明, 刘佑荣,等. 考虑裂纹闭合效应的岩石损伤力学模型及耗散能量分析[J]. 长江科学院院报, 2016, 33(5):69-75.

[6] LEMAITRE J. A Continuous Damage Mechanics Model for Ductile Fracture[J]. Journal of Engineering Materials & Technology, 1985, 107(1):83-89.

[7] LEMAITRE J. How to Use Damage Mechanics[J]. Nuclear Engineering & Design, 1984, 80(2):233-245.

[8] KRAJCINOVIC D, SILVA M A G. Statistical Aspects of the Continuous Damage Theory[J]. International Journal of Solids & Structures, 1982, 18(7):551-562.

[9] KRAJCINOVIC D. Continuum Damage Mechanics: When and How?[J]. International Journal of Damage Mechanics, 1995, 4(3):217-229.

[10] 唐春安. 岩石破裂过程中的灾变[M]. 北京:煤炭工业出版社, 1993.

[11] 谢和平. 岩石、混凝土损伤力学[M]. 徐州: 中国矿业大学出版社, 1990.

[12] 曹文贵, 赵明华, 刘成学. 基于 Weibull分布的岩石损伤软化模型及其修正方法研究[J]. 岩石力学与工程学报, 2004, 23(19):3226-3231.

[13] 曹文贵, 莫 瑞, 李 翔,等. 基于正态分布的岩石软硬化损伤统计本构模型及其参数确定方法探讨[J]. 岩土工程学报, 2007, 29(5):671-675.

[14] 徐卫亚, 韦立德. 岩石损伤统计本构模型的研究[J]. 岩石力学与工程学报, 2002, 21(6):787-791.

[15] 袁小平, 刘红岩, 王志乔. 基于Drucker-Prager准则的岩石弹塑性损伤本构模型研究[J]. 岩土力学, 2012, 33(4):148-153.

[16] 蒋 维, 邓 建, 李 隐,等. 基于对数正态分布的岩石损伤本构模型研究[J]. 地下空间与工程学报, 2010, 6(6):1190-1194.

[17] DENG J, GU D. On a Statistical Damage Constitutive Model for Rock Materials[J]. Computers & Geosciences, 2011, 37(2):122-128.

[18] 刘成学, 杨林德, 曹文贵. 岩石统计损伤软化本构模型及其参数反演[J]. 地下空间与工程学报, 2007, 3(3):453-457.

[19] 杨林德. 岩土工程问题的反演理论与工程实践[M]. 北京: 科学出版社, 1995.

[20] 巨能攀, 黄海峰, 郑 达,等. 考虑含水率的红层泥岩蠕变特性及改进伯格斯模型[J]. 岩土力学,2016, 37(增2):67-74.

[21] 黄海峰, 巨能攀, 周 新,等. 红层滑坡滑带土经验型蠕变模型研究[J]. 人民长江, 2017, 48(5):91-95.