环境湿度影响下粉晶白云岩的力学特征与数值模拟

2023-05-22王卫军赵延林吴秋红

煤炭学报 2023年3期

陈伟,万文,王卫军,冯涛,赵延林,3,吴秋红,3,周彧

(1. 湖南工程学院建筑工程学院,湖南湘潭 411104;2. 湖南科技大学资源环境与安全工程学院,湖南湘潭 411201;3. 湖南科技大学南方煤矿瓦斯与顶板灾害预防控制安全生产重点实验室,湖南湘潭 411201)

在非金属矿山的开采中,为了高效开采和节约成本,往往采用房柱法开采[1-2]。开挖时所留设的保安矿柱中通常含有伊利石,高岭石,蒙脱石,绿泥石,云母等亲水性矿物,这些矿物成分易水解发生软化[3-5]。一直以来,许多学者在研究水化学溶液对岩石力学性能的改造作用方面取得了巨大成就。ZHAO等[6]发现,节理岩体在水-岩作用下其强度与刚度的劣化会严重影响库岸边坡稳定性。丁梧秀等[7]利用不同的水化学溶液对灰岩进行侵蚀,发现灰岩产生损伤劣化程度受控于溶液pH与凝结核丰度。MIAO等[8]研究了花岗岩在各pH以及流速条件的化学溶液中浸泡后,其单轴、三轴的变形规律和强度损伤,并借助电子能谱和扫描电镜从微细观角度总结了其化学动力反应原理。韩铁林等[9]对砂岩进行了水化学溶液浸泡并开展了一系列力学试验,发现砂岩的劣化参数与强度呈线性递减关系。JIANG等[10]探讨了单轴压缩下砂岩经AMD溶液腐蚀后的损伤本构关系,并构建相关模型,该模型能较好的呈现砂岩在AMD蚀化下的劣化规律。ATKINSON等[11]及DUNNING等[12]发现了无机水化学溶液对石英的断裂力学性能的影响规律,并用Griffith强度理论对岩石受水化学环境的作用机制作出理论分析。

上述的诸多研究主要是针对各种水化学溶液以浸泡方式对岩石力学特性产生腐蚀作用而展开,并未涉及水化学流体以气态方式侵入岩体。近年来,深部开采进行得如火如荼,由于地温的升高,结合地下水发育,蒸发量加大,深部矿房的环境湿度平均值会超过80%RH[13]。且矿柱在长期的地质构造运动过程中其内部形成了断层、解理、裂隙等缺陷,缺陷的存在给气态水分子提供了渗透通道,为水汽侵蚀矿柱创造了有利条件[14-17]。因此,关于矿柱在高湿环境下的损伤力学特性亟待研究。

针对上述情况,笔者对湖南郴州茶山岭煤矿-850 m房柱法开采段进行了实地踏勘,其水文地质报告显示,该区常年地下水发育并呈弱酸性,相对湿度处于90%RH这一高位。为了对比不同环境湿度值对矿柱力学性能的削弱程度,笔者对现场矿柱取心制备标准试件并模拟3种高湿环境,通过MTS-815型试验系统对高湿处理60 d后的试件展开单轴压缩试验,总结白云岩微细观结构在各湿度环境下的变化特征,最后利用PFC离散元软件的二次开发虚拟实现同等条件下的湿度场水汽侵蚀数值试验,结合模拟试样的微裂纹发育特征与能量演化规律,揭示在湿度扩散效应影响下白云岩能量演化的微细观机制。

1 高湿侵蚀矿柱试样原理

1.1 试验系统组成

1.2 试样制备

对茶山岭煤矿-850 m的3号矿房5号矿柱沿沉积方向取心后加工成50 mm×100 mm的圆柱试样[23],如图2所示。试件制备结束后选取外观较为良好且无瑕疵的试样,以期获得可靠的试验结果。为了避免风化,将试件密封并以塑料薄膜封装。利用煮沸法及AiniMR-60核磁共振分析仪分别测得试件平均含有3.42%的饱水率及4.93%的孔隙率。

本次采样点矿房温度处于25～30 ℃,且湿度对试样的力学损伤是此次研究重点,故以恒定25 ℃为试验温度。表1为试验所选用白云岩试样的基本参数。

图2 岩心采取与试样制备Fig.2 Core collection and specimens preparation

表1 白云岩试样基本参数

1.3 加湿与加载程序

以湿度控制箱对试样进行加湿处理:先将酸性Na2SO4溶液按1 ℃/min速率加热到25 ℃;再利用加湿系统分别设定80%RH、90%RH、100%RH试验湿度,并持续工作60 d;最后,擦干试样表面,并遵循控制变量原则[24],在DZF-2060型真空干燥箱中以120 ℃干燥24 h。

2 试验结果分析

2.1 粉晶白云岩物理特性变化规律

2.1.1 不同湿度作用下白云岩孔隙率的变化规律

气态水分子侵入试样后,水解、溶蚀部分亲水性矿物,并通过外渗将反应物运移,岩石的孔隙率因其内部产生的溶洞连通为渗透通道而增加,其孔隙压力和渗透特性因此受到影响[21],试样结构疏松度、比表面积和次生孔隙等增大。利用NMR核磁共振技术获取白云岩在不同湿度环境下60 d后的孔径分布,如图3所示。T2为饱和样品孔隙中含水原子核的横向弛豫时间,利用外部静磁场所测得。鉴于试样横向弛豫时间T2与其孔隙尺寸成正比[25],其表达式为

图3 4种湿度环境下放置60 d后的白云岩T2谱规律Fig.3 T2 spectra of dolomites under four humidity conditions for 60 d

(1)

式中,ρ为饱和样品T2表面弛豫强度;S为岩石孔隙表面积;V为流体体积。

基于孔隙的半径与喉道大小的正比关系,式(1)可表达为

(2)

式中,r为孔隙半径;Fs为样品孔隙喉道的几何形状因子。

令1/(ρFs)=N(转换系数),式(2)可化简为

T2=Nr

(3)

天然岩石具备多孔介质属性,因此目前还未能形成统一的标准对其进行孔径划分,文献[26]以10 μm为界,将岩石中孔径区分为大孔(r>10 μm)、小孔(r≤10 μm);针对低渗透岩石,文献[27]划分为大孔(1～5 μm)、中孔(0.1～1.0 μm)、小孔(r<0.1 μm)。笔者参考毛细管压力测量孔隙半径分级方法,结合试验用样的孔径集中分布的范围,将孔径划分为小孔(r≤ 0.1 μm)、中孔(0.1 μm

T2=10r

(4)

由图3结合式(4)可将粉晶白云岩的T2谱中的横向驰豫时间分布在0～1 ms划分为小孔隙、1～10 ms为中孔隙、10 ms以上为大孔隙。试件内不同孔径的谱面积分布及占比统计见表2。

表2 试件内不同孔径、不同环境湿度的谱面积分布及占比

干燥试样主要分布在0.1～10.0 ms弛豫时间内,且振幅远低于经过高湿处理后的样本,其T2谱总面积为3 255.26。对比来看,80%RH环境下试样的T2谱总面积为160 631.61,增幅达49倍之多,且在湿度逐渐接近100%RH过程中,波峰不断向右偏移,表明在高湿环境下与试件内部亲水化合物发生剧烈的物理、化学反应,孔隙结构持续劣化。值得说明的是,100%RH环境下试样谱线出现3峰,这是由于试样内部黏土矿物吸水膨胀内部挤压,导致一部分中、大孔隙体积减小转化成小孔隙,小孔谱面积增大,即形成左侧第一波峰。具体反应如下所述。

粉晶白云岩中主要碎屑矿物为正长石和云母,在酸性气态水分子侵蚀下产生溶蚀孔洞,形成黏土矿物:

(5)

(6)

生成的黏土矿物高岭石(Al2Si2O5(OH)4),其粒径不大于2 μm,具备较大表面能,包括强亲水性与膨胀性。高岭石吸水膨胀化学式为

(7)

测定干燥状态及3种高湿环境作用下60 d并烘干后白云岩试样的总孔隙率(包括开口与闭合两部分孔隙),计算公式为

n=(ρp-ρd)/ρp×100%

(8)

式中,n、ρd、ρp分别为试件的总孔隙率、干视密度、真密度。

从统计结果,以及白云岩表面的高清照片(50倍率,图4)可知,岩样矿物颗粒间的孔隙结构随环境湿度的升高而愈发松散。

具体来看,白云岩孔隙率的增幅呈现明显的阶段性变化趋势,当试样所处环境从干燥状态转为80%RH时,其平均孔隙率增幅较大,由4.96%增至7.41%,增幅高达49.39%。100%RH高湿作用后的白云岩平均孔隙率较80%RH湿度作用后累计增幅达8.09%,符合负指数函数增长。

2.1.2 不同环境湿度作用下白云岩质量的变化规律

高湿作用带来白云岩的孔隙率增加,因此有必要对各湿度环境作用前后的质量差异进行分析[28]。称得3种经高湿环境处理60 d并冲洗外表松散脱落物后烘干的白云岩试件的质量,计算出质量差。根据前后质量差,可求出试样的相对质量变化率,即质量损失因子,定义为D,即

图4 4种湿度环境下放置60 d后的白云岩试样孔隙率Fig.4 Porosity of dolomite specimens under 4 humidity conditions for 60 d

(9)

其中,M0、M(h)、ΔM(h)分别为岩样的初始质量、经各湿度环境作用60 d并烘干后的质量、经各湿度环境作用60 d前后质量差。不同高湿作用前后各岩样的质量变化与质量损失因子分布如图5所示。

图5 3种高湿环境下放置60 d的白云岩试样的质量变化对比Fig.5 Comparison of mass changes of dolomite specimens under 3 humidity conditions for 60 d

从图5(a)可以看出,白云岩试件进行各湿度作用后其质量皆有不同的减少,湿度越大质量损失越严重。试样的平均质量损失在100%RH环境下为5.23 g,是80%RH环境下的5.03倍及90%RH环境下的2.41倍。同时,观察图5(b)、(c)发现,环境湿度越大,试样的平均质量损失因子越大,且增幅不断上升,符合幂函数增长趋势。试样在100%RH湿度下其平均质量损失因子是0.95%,远高于80%、90%RH湿度下的0.19%和0.39%。

2.2 粉晶白云岩力学特性劣化规律

本次试验中20组试件的单轴压缩力学参数见表3。

表3 白云岩单轴压缩力学参数

图6展示了各湿度环境处理60 d后所有试样的单轴压缩应力-应变曲线。

不同湿度条件下,试件的单轴压缩应力-应变曲线可分为(Ⅰ)～(Ⅳ)阶段,分别为裂纹闭合压密阶段(Ⅰ),弹性变形至裂纹稳定扩展阶段(Ⅱ),裂纹非稳定扩展阶段(Ⅲ)与峰后扩容阶段(IV)[4]。分析各阶段的曲线特征以便进一步解读白云岩矿柱(岩心)的渐进破坏过程。

(Ⅰ)为裂纹闭合压密阶段,在此阶段,白云岩内部的微裂隙结构出现了闭合。此时应力-应变曲线呈下凹型,试样的变形模量逐渐增大但侧向膨胀较小。90%RH和100%RH湿度处理过的试样甚至出现了微弱的侧向收缩。主要是由于这2种环境下,空气中的水分子密度较高,水岩接触界面将会产生湿度梯度,迫使气态水分子由高质量浓度区跃迁至低质量浓度区。湿度梯度越大,水分子对内部矿物组分中的亲水性矿物溶蚀随之加大,形成空隙。加载初期,孔隙被压实。随着载荷的增加试样体积逐渐减少,此时轴向压缩占主体地位。随着环境湿度的增加,其裂纹的闭合应力降低,由35.97 MPa(干燥)降为13.56 MPa(100%RH),降幅达62.31%(图7(a))。

试样内部裂纹闭合完毕后进入弹性变形至裂纹稳定扩展阶段(Ⅱ)。白云岩在此阶段的应力-应变曲线近似直线段。参数泊松比μ和弹性模量E则用来反映试样这一阶段的变形特性。结合式(10)、(11)求得平均弹性模量和泊松比[29],结果如表3和图7(b)所示。

(10)

(11)

式中,参考试样应力-应变曲线近似直线段(图6),起点A和终点B的应力分别为σA和σB;轴向应变分别为εhA和εhB;径向应变分别为εdA和εdB。

应力-应变曲线上直线段的E越小说明水分子对白云岩的侵蚀程度越深[22]。由图7(b)可知,白云岩在3种湿度环境处理过后的弹性模量均低于干燥状态,随湿度值呈负指数函数形式降低并趋于稳定。相对湿度为80%、90%、100%环境下60 d后,试样平均弹性模量分别为10.39、9.34和8.41 GPa,与干燥试样相比分别降低了54.08%、58.76%和62.88%。由此可见,化学作用引起的腐蚀软化效应会随着环境湿度的上升而增强。此外,白云岩的变形参数泊松比与环境湿度值呈正相关。与干燥状态相比,在100%RH条件下的白云岩泊松比由0.314增至0.826,累积增加量为163.41%。

当轴向应力达到了损伤应力时,表明试样进入了裂纹非稳定扩展阶段(Ⅲ)。弹性变形至裂纹稳定扩展阶段(Ⅱ)和裂纹非稳定扩展阶段(Ⅲ)以扩容起始应力点为界(损伤应力),即体积压缩与体积膨胀的拐点。各湿度环境下试样的损伤应力在31.25～63.22 MPa,为峰值应力的51.86%～58.49%。该阶段的轴向应力-应变曲线改变直线状转而呈现上凸形态,即进入塑性变形阶段。干燥、80%RH与90%RH环境下,轴向应力-应变曲线在峰值处产生应力平台并波动(图6(b)、(c))。此阶段,新生新裂纹持续不断的扩展致使试件的侧向应变呈非线性增加。试样内部微裂纹相互交汇连通导致试件出现宏观破坏。在每种湿度下放置60 d后,干燥状态和100%RH环境下试样的平均峰值应力分别为121.91和53.44 MPa,明显降低了41.05%。相反,随着环境湿度值的增加,其试样峰值点的轴向和侧向应变不断增加:干燥状态和100%RH环境下试样的平均轴向应变分别为0.72×10-2和0.96×10-2,增幅达32.78%;此外,这2种状态下的平均侧向应变分别为0.43×10-2和1.03×10-2,增加了137.97%。

图6 不同相对湿度作用60 d后白云岩的单轴压缩应力-应变曲线Fig.6 Uniaxial compressive stress-strain curves of dolomite specimens under different humidity actions for 60 d

这里,分别定义试样峰值点应力的劣化度Q(h)和峰值点轴、侧向应变的增长率P(h)、R(h),用来定量分析环境湿度对白云岩强度变形的影响规律,分别为

(12)

(13)

(14)

式中,σc(h)为岩样经不同高湿环境处理60 d后峰值强度;ε1c(h)、ε2c(h)分别为其峰值点轴向和侧向应变;σc0、ε1c0和ε2c0为干燥环境岩样峰值点强度、轴向应变和侧向应变。

图7(c)展示了60 d的处理时间过后,白云岩单轴压缩下峰值点的应力、轴向和侧向应变随环境湿度值而发生的典型阶段变化率。白云岩试样在峰值点的强度的劣化度和轴向、侧向应变的增长率均随着环境湿度的升高而增加。与干燥状态相比,同样的湿度作用时间60 d后,经80%RH、90%RH和100%RH处理后的白云岩峰值强度平均劣化度随湿度呈负指数函数上升,分别为42.37%、54.14%和60.28%,分布曲线的增幅具有显著非均匀性,即由干燥环境进入80%RH环境后,两者的劣化度增幅较大,由80%RH环境进入100%RH环境导致的两者劣化度的增幅相对较小。80%RH、90%RH和100%RH处理后的白云岩较干燥试样峰值点轴向、侧向应变对应的平均增长率随湿度值线性增加,轴向应变增长率分别为18.76%、21.52%和33.01%,侧向应变增长率分别为45.58%、92.79%和140.32%。

图7 不同相对环境湿度下60 d后白云岩的单轴压缩力学参数Fig.7 Uniaxial compression mechanical parameters of dolomite under different humidity conditions for 60 d

应力达到峰值强度后,试样进入峰后扩容阶段(IV)。此时,轴向应力-应变曲线分别呈断崖式跌落(图6(a))和峰值后继续发生部分残余变形这2种状态(图6(b)～(d))。第1种下跌状态是由于干燥试样的结构具有致密性,呈现出较强的脆性,应力出现骤降现象的同时往往伴随脆裂声,且产生短时间的峰值软化。另一种下跌中试样峰后残余变形会随着湿度的增加而增加。这一阶段,裂纹发展速度加快,侧向应变增长加剧,且发生多次应力降,每次应力降都会导致侧向应变急剧增加而呈延性破坏。因此,白云岩随环境湿度的升高会由脆性破坏向延性破坏过渡,且加剧试样峰后的累积残余变形。

3 基于微细观损伤的数值模型

3.1 SEM电镜扫描分析

对各湿度环境处理后的白云岩侧面进行SEM电镜扫描拍摄(500倍率),选取部分典型图片展示(图8)。观察可知,干燥的白云岩试样在切割打磨后表面结构相对稳定致密,基本无微裂隙发育,不同晶体之间存在层理面,矿物颗粒整齐排列,完全被胶结物覆盖,形态呈现明晰片状及块状结晶(图8(a))。这说明在进入高湿环境前,白云岩试件的宏观力学性能良好。经80%RH环境60 d后,原规则紧密的块状、片状结晶结构部分胶结物在水化作用下溶解,展露出海绵状的团絮形态,层间隆起,形成陡峭不一的“山脊”(图8(b))。90%RH环境下,外表面部分矿物颗粒由于水化侵蚀溶解了胶结物,削弱了颗粒间的粘结力,部分颗粒呈凸出状。与此同时,试样表面发育有蜂窝状的次生孔隙(图8(c))。100%RH环境下,表层矿物颗粒间的胶结物继续水化分解,片状晶体在溶蚀作用下产生孔洞。相较80%RH的试件,矿物颗粒间胶结物缺失严重且出现颗粒大面积掉落现象,孔隙通道交汇连通产生微裂隙,整体微细结构趋于疏松,规则排列的晶体颗粒消失殆尽,白云岩微细观结构损伤愈发严重。

图8 不同湿度环境下60 d白云岩表观电镜扫描Fig.8 Scanning electron microscopy of dolomite surfaces under different humidity conditions for 60 d

由于矿柱上、下端面与顶底板连接,并未与水蒸气直接接触。因此,为了得到湿度扩散引起的白云岩矿柱腐蚀深度变化规律,将试件上、下端面用薄膜覆盖,仅侧面用来接受气态水分子的侵蚀作用。为了定量分析白云岩在高湿环境下的微细观损伤规律,将微细观损伤定义为电镜观测面中的微裂隙和微孔隙。对干燥和3种湿度环境中处理60 d后的白云岩试件的顶底面利用扫描电镜进行拍摄。为了拍摄精度,已对仪器内置的载物台进行了坐标网格划分(即以载物台中心为圆心,半径为2.5 cm圆形)共计78张电镜扫描图,作为损伤原始面(图9)。

图9 载物台坐标网格划分Fig.9 Coordinate gridding of loading platform

载物台的尺寸经后期厂家改装为半径5 cm的圆盘,可容纳高度13 cm的试样进入拍摄。对各试件顶、底面按图9中的编号进行依次逐点扫描。根据前期尝试的效果,以及采集清晰度的需求,最终选取1 000倍作为本次定点定倍的拍摄方式。在SEM图中,相对于完整部位来说,损伤部位的灰度为0。表4中利用Coreldraw软件对每张扫描图进行二值化处理。根据软件内置的柱状图工具得到最佳灰度分割阈值为97,从而得出原始面每个点的SEM的“损伤因子”O[30]为

(15)

式中,以二值化图片为参考,LD为数值为0的面积;L为总面积。

表4展示了各湿度环境下白云岩试样上、下端面的扫描图以及二值化图片。

平均损伤因子选取试件上、下端面对应两点的损伤因子平均值,图10展示了试件上、下端面平均损伤因子的频率分布。考虑到干燥试件存在天然缺陷以及加工制备期间端面产生的初始损伤,因此其电镜扫描图的损伤因子不为0。观察可知,每个试件上、下端面的平均损伤因子分布具有明显差别,以图10中的数据经三参数峰值函数拟合后来直观对比试件置于不同湿度下的损伤程度,结果表明白云岩上、下端面的平均损伤因子频率直方图符合Gauss函数分布。式(16)为Gauss函数拟合公式:

表4 白云岩试样上、下端面的扫描图以及二值化图片

(16)

式中,J为比例参数;w为形状参数;xl为位置参数。

干燥试件损伤因子在0.02～0.04内占总拍摄点数40.08%,而100%RH环境下试件的损伤因子在0.04～0.06内占总拍摄点数29.21%,反映出高斯曲线的波峰会因环境湿度升高而右移(图10)。可以看出,损伤因子分布逐渐由低损区往高损区过渡,环境湿度增加直接加剧了微观损伤的形成。其内部会在长期的高湿环境下产生一定的细观损伤累积,这种岩体的强度渐进劣化作用会影响深部岩石工程的稳定性。

基于图9中的原始面二值化图片分布坐标网格,可得到各试件上、下原始面平均损伤因子分布(图11)。由图11可知,对比干燥试样,经高湿侵蚀作用后,试件边缘区多见高损伤区,中心区偶见高损伤区。随着湿度的增加,细观损伤由边缘侧逐渐向中心发育,不断形成的高损区以“损伤-联结-侵入”往试样中心延展,且湿度越高,损伤侵入的深度越深。

图10 白云岩试样上、下端面的平均损伤因子频率直方图与拟合曲线Fig.10 Average damage factor frequency histogram and fitting curves of dolomite specimen top and bottom end faces

3.2 高湿侵蚀深度

上述微观表征显示高湿环境下白云岩试样产生了由外及内的侵蚀作用。对各湿度环境下处理60 d后的试样冲洗表层脱落颗粒后,进行尺寸测定,图12展示了各环境下5个试样的平均直径变化曲线。

图12显示,白云岩试样经高湿作用60 d后由于表层受剥蚀脱落其直径减小,然而随着湿度值越高,下降速度随之平缓。进入80%RH环境后,试件直径呈最大降幅,为1.06 mm,缩减率达2.31%。100%RH对比80%RH环境下的试样直径累计减小量为0.15 mm,缩减率为0.31%。由图12可知,80%RH环境中的试样较干燥试样而言,气态水分子会对白云岩产生强烈侵蚀作用,致使其外层发生较为严重的剥蚀脱落,80%RH～90%RH侵蚀速率相对较缓。

3.3 湿度场数值模型的建立

根据湿度传质理论可知,存在湿度差的2种物体一旦产生接触面就会形成湿度梯度,迫使气态水分子由高浓度组分移至低浓度组分,通过调节最终趋于平衡[31]。这个过程称为湿度扩散,本次试验中的白云岩与湿度环境通过扩散效应进行能量交换[32]。能量交换的主要方式是分子扩散和对流扩散:分子运动引起的能量交换被称为分子扩散,主要存在于流体和固体之间;而对流运动引起的物质转移被称为对流扩散,主要存在于2种流体之间。2种传递方式相比,分子扩散的物质传递比较平稳[33]。在深部矿房中,主要存在空气和矿柱之间水分的迁移,以分子扩散进行湿度的传递。基于传质理论,空气和白云岩矿柱之间的水分能量交换如式(17)所示,表示单位时间内交换的对流水量。交换对流水量的影响因素众多,比如对流交换系数,能量交换时间和接触面积等。其中,不稳定因素是对流交换系数。

Hf=mC(vw-vf)

(17)

式中,Hf为矿房中的湿交换量,kg/(m2·s);mC为湿度交换系数,m/s;vw为矿柱水岩接触面的饱和水汽质量浓度,kg/m3;vf为矿房中对流层的水汽质量浓度,kg/m3。

水汽流动本质上属于多孔介质渗流,气体在多孔介质附近发生渗流的阻力公式[34]为

图11 白云岩试样上、下端面的平均损伤因子分布Fig.11 Average damage factor distribution images of dolomite specimen top and bottom end faces

图12 试样的直径和直径变化量随环境湿度变化曲线Fig.12 Diameter and diameter variation curves of specimen with environment humidity

(18)

式中,Δp为压降,N/m2;l为矿房中渗流带长度,m;я为矿房中流体黏度,N/(m2·s);ο为矿柱离散元颗粒平均直径,m;g为矿柱孔隙率,%;ρ0为矿房水汽质量浓度,kg/m3;v为矿房风流流速,m/s。

利用PFC2D模拟水汽传输,动量方程中的源项分别表示为黏性、惯性损失,如

(19)

式中,等号右侧首项、次项分别为黏性、惯性损失系数项;Ki为x、y、z方向上矿房中多孔介质在动量方程中的源项;Aij为黏性损失的预定义矩阵;Bij为惯性损失的预定义矩阵;vj为x、y、z方向上的风流流速的分速度。

矿房水汽流动方向的黏性阻力系数q和惯性阻力系数I2为

(20)

(21)

式中,Ap为颗粒粒径。

在上述由外及内的细观损伤基础上,引入应力腐蚀理论建立了高湿作用下矿柱的颗粒离散元接触模型。在恒载地应力下,水分子持续腐蚀矿柱内部结构,矿柱的强度变形随着时间的推移不断劣化。编写了FISH函数将PBM的模型进行改造,结合细观侵蚀深度对矿柱试样模型进行二次开发,从而实现了应力腐蚀模拟:基于PBM模型和断裂力学中KI型张拉裂纹和KⅡ型剪切裂纹(图13)的亚临界扩展速率的概念定义损伤速率v(计算见式(22)),不断缩减颗粒的黏结半径,得到试样时效损伤变形的腐蚀数值模型。

图13 湿度扩散引发颗粒间接触力损伤的应力腐蚀模型Fig.13 Stress corrosion model for moisture diffusion-induced interparticle contact force damage

(22)

其中,rb为PBM线性平行黏结本构模型的黏结半径;N1、N2为岩石材料的相关性系数;σb0、σb1、σb2分别为PBM线性平行黏结本构模型的法向应力、法向抗拉强度与起裂应力。采用PFC建模前,首先需要调试试样的微观力学参数,即依据室内试验的应力-应变曲线,标定模型中的宏微观参数。根据应力-应变曲线特征,确定干燥试样的起裂强度σci为26.51 MPa左右。表5为最终选用的细观力学参数。

值得说明的是,本次建模未能考虑非饱和渗流过程对岩样力学行为的影响。我们通过简化模型,直接以试样在非饱和吸水状态下的水汽入侵深度引入非饱和渗透系数。将粉晶白云岩试样颗粒的润湿性直接定义成微细观损伤的损伤程度,认为颗粒之间的微张拉力和微剪切力的削弱是引起粉晶白云岩试样微细观时效损伤演化的内在驱动力[35-36]。

表5 数值模型的微观和宏观性质

结合图11的湿度扩散导致的白云岩损伤速率,并通过自编程FISH语言来实现两侧加载板的湿度传递程序,计算得出3种环境湿度下试样60 d内湿度扩散损伤区分布云图(图14),此处的大气压力是以标准大气压为参考。

由图14可知,在矿柱(岩心)上、下隅角处湿度扩散速度最小,损伤最小,越靠近矿柱中段位置,腐蚀深度越深,且损伤分布面呈“扇型”对称分布。由于水汽流的不断涌入,矿柱中段容易接触的水分子面积较大,且渗透压力随着垂直方向开始不断减小。相同的作用时间下,更高湿度下水汽侵蚀的深度较深。

图14 不同湿度环境下60 d白云岩湿度扩散分布云图Fig.14 Moisture dispersion distribution clouds of dolomite under different humidity conditions for 60 d

3.4 数值模拟结果分析

图15对比了室内试验和PSC模型所得到的单轴压缩应力-应变曲线。由于模型所赋参数是基于试验数据的平均值,而不是单个试样,因此2者的曲线之间存在一些误差。但仍可以注意到,模拟试样的应力-应变曲线和破裂形态与试验结果基本一致。

单轴压缩下,试样只受轴向应力,因此岩样的侧向变形往往没有约束。在本试验中,试样的最终破坏主要包括2种类型:拉伸破坏和混合剪切/拉伸破坏(图15)。

图15 白云岩受压数值模型与室内试验结果对比:应力-应变曲线与破裂形态Fig.15 Comparison of numerical modeling with laboratory results of dolomite under compression:Stress-strain curves and failure modes

干燥状态下的白云岩主要呈现拉伸破坏,并且破坏形态会随着环境湿度的增加转变为混合剪切/拉伸破坏。天然岩样破裂面较顺直,起伏不明显,用手触摸摩擦感不强,且破裂面上有少量岩屑分布;100%RH处理后的岩样破裂面凹凸不平、起伏明显,用手触摸摩擦感较强且有明显湿润感,破裂面上岩屑分布多、颗粒大。由外及内的结合水膜的润滑作用以及孔隙水压力使得颗粒间胶结强度减弱,降低了岩样外层原生裂隙尖端的应力强度因子。在单轴压缩下,试样外层最先发生起裂并呈斜对角线连接和贯通,最终导致宏观剪切裂隙数量增加[37],致使白云岩宏观破坏由拉伸破坏发展为剪切/拉伸混合破坏。

此外,观察图15与未示出的试样发现,岩样破坏时其表面发育新裂纹的数量会随湿度值的增大而增多。这是由于干燥试样受载前在长期的地应力作用下,内部缺陷已经形成稳定的力学状态,因此荷载作用下表面只萌生导致试样形成最终破坏的宏观裂纹。另一方面,水分子侵蚀改变了试样的内部结构的稳定性。因此,在相同的加载速率下,新生裂纹由于松散结构的重组充分扩展,试样表面的宏观裂纹急剧增多。试样数值模拟结果也显示两侧的裂纹数量随湿度的升高而不断增加,与室内试验结果基本吻合,证明所构建的腐蚀离散元颗粒模型(PSC)是可以反映现场高湿环境对白云岩矿柱力学特性由外及内的侵蚀损伤。

4 讨论

4.1 微裂纹发育特征

对PSC模型进行加载时,一旦其颗粒黏结力低于内部接触强度,黏结力链即发生断裂,导致平行黏结接触转化为线性接触。这标志着PSC模型内部开始萌生微裂纹。试样微裂纹的数量与类型极大地影响了颗粒的黏结半径与颗粒平均半径的比值。当黏结力小于接触面颗粒的抗拉强度时,会产生张拉微裂纹。反观黏结剪切强度大于黏结力时,剪切微裂纹随之发育。图16比较了4种不同湿度状态下模拟试样的微裂纹发育数量特征。

图16 模拟试样微裂纹发育特征Fig.16 Microcrack developmental characteristics of samples

可见,当试样发生破坏时,内部微裂纹总数量随着环境湿度的升高而增加,但张拉微裂纹数的增幅远低于剪切微裂纹的增幅。如干燥状态下的拉伸微裂纹和剪切微裂纹数分别是124和98,以拉伸微裂纹为主,而100%RH环境下的两者数量基本接近,分别为689和736,涨了4.5和6.5倍。试验过程中白云岩微裂纹的发育特征与数值模拟结果基本类似,说明在湿度扩散效应下由外及内的渐进侵蚀作用劣化了试样的强度。同时,湿度越高,白云岩试样压缩破坏后更为破碎,呈现明显的剪切特性和较弱的张拉特性。根据模拟结果可以得出,不同湿度状态下试样的微裂纹集中萌生于塑性变形段并在荷载达到峰值强度后数量陡增。这也进一步说明了峰后宏观裂缝断面的贯通和摩擦滑移促使矿物颗粒间的黏结断裂,导致裂纹数量剧增。通过试验所得出的耗散能演化特征与PSC模拟所得出的微裂纹发展规律相一致,揭示了高湿环境对白云岩侵蚀后受载所产生裂纹的能量机制。

4.2 能量演化的微细观机制

单轴压缩下PSC离散元模型的数值模拟,所有颗粒之间的能量关系为

Eboundary=Estrain+Eslip+Epbstrain+Ekinetic

(23)

式中,Eboundary为颗粒边界能;Ekinetic为颗粒动能;Eslip为颗粒间摩擦滑移能;Estrain为颗粒间应变能;Epbstrain为颗粒间黏结能。

PSC模型颗粒间的能量形式主要有2种,分别是动能和接触能。颗粒动能主要是模型受外力作用所吸收的边界能,类似于总吸收能U。接触能主要是储存线性弹簧中的应变能以及平行黏结弹簧中的黏结能,类似于弹性应变能Ue;其次为耗散于微裂纹形成与剪切滑移中的摩擦能,与耗散能Ud类同。所以,受载时PSC模型的边界能约等于应力达到峰值前的颗粒应变能、摩擦能和黏结能的总和。

为了方便对比,图17给出了干燥状态和100%RH环境下白云岩PSC模型在单轴压缩下的能量演化监测曲线。可知:第1,加载过程中所做的功大部分转化为颗粒的应变能。当试样中的微裂纹开始出现时,应变能开始释放并转化为颗粒耗散能和局部阻尼能。第2,轴向荷载达到峰值以后,应变能释放并转化为摩擦断裂能,表明断裂后颗粒的摩擦作用直接导致峰后残余强度的形成。图17还显示,在应力达到峰值之前,应变能和黏结能的总和在边界能中的占比超过了90%。当总值超过能量累积的极限,储存在模型内部的能量会随试样破坏而释放,这个储能极值受控于岩石材料的强度劣化度。此外,这与试验中高湿作用降低白云岩的储能极限相吻合。随着微裂隙界面摩擦作用的加剧,摩擦能也在不断地攀升。与此同时,在应力达到峰值时,宏观裂隙发生贯通,导致摩擦能出现陡增的趋势。这与室内试验中摩擦耗散能随岩样塑性变形的发展而增加,峰值后其能量值急剧增加相一致。

图17 PSC模型试样的能量演化曲线Fig.17 Energy evolution curves of PSC model specimens

图17中能量演化特征的计算结果表明:高湿环境下白云岩的边界能、黏结能和应变能远小于干燥状态下试样所对应的能量值。由此可知,环境湿度的升高降低了白云岩的总吸收能U与弹性应变能Ue,即高湿作用使得试样的吸能性能减弱。对比干燥状态和高湿环境下试样的摩擦能,白云岩在高湿侵蚀后的能量耗散略有增加,说明湿度扩散效应导致内部裂隙发育,其微裂纹数量的增加提升了颗粒间摩擦滑移耗散能。

在单轴压缩条件下,根据白云岩模型在干燥和3种湿度状态下的模拟试验数据,分别从能量演化角度进一步揭示白云岩受环境湿度影响的微细观机制:① 白云岩在轴向应力的持续做功下不断吸收能量,并以颗粒黏结能与应变能结合的形式储存和积累,形成边界能;高湿侵蚀迫使颗粒间的黏结发生断裂,因此进入塑性阶段时首先需要克服的是颗粒间的黏结能,然后微裂纹才能随着应变能的增加而继续扩展;摩擦能对裂纹滑移有着决定性的作用,其占边界能的比例随着裂纹扩展和滑移程度的加剧而增加。② 应力达到峰值后,储存的黏结能和应变能迅速释放;高湿处理后试件的微裂纹数量以及发育、贯通程度都在增加,使摩擦能在总能量中占据主导地位;然后,剪切微裂纹迅速扩展并贯通,成为宏观裂纹的主要表现形式。