图1 渗流求解流程
Fig.1 Seepage solving procedure
郭 影1,2，姜忻良1，曹东波3，白铁钧2，朱广轶2，冯 春4
(1. 天津大学建筑工程学院，天津 300072；2. 沈阳大学建筑工程学院，沈阳 110044；3. 沈阳工程学院公外部，沈阳 110136；4. 中国科学院力学研究所，北京 100190)
摘 要：传统的渗流数值计算方法难以较真实地描述岩石材料性能劣化和渗透性演化机制。该文提出了一种渗流吸水诱发岩体强度弱化的有限体积数值计算方法。利用高斯散度定理采用有限体积法求解水在岩体中的渗流过程及岩体的变形破坏过程，建立了基质吸水引起岩体模量及强度弱化的理论模型及相应的数学表述。基于考虑吸水弱化算法的各向同性孔隙渗流模型，模拟了低、高两种边界流体压力下某粉砂岩试样的吸水软化过程和不同吸水时间下试样的单轴压缩过程。数值算例表明：边界流体压力越高，试样达到整体饱和状态的时间越快；渗流初期以自由水渗流填充孔隙为主，渗流后期以孔隙内自由水向基质吸水的转化为主；边界流体压力对渗流速率具有明显控制作用，但对基质吸水速度无影响；随着吸水时间的增加，试样的强度(黏聚力、内摩擦角)逐渐减小至残余值，得到的基质吸水含量随时间变化的数值解与理论解基本一致，表明了数值算法的计算精度，可以用于隧道突水、围岩稳定性等实际岩体工程问题的渗流-应力耦合效应分析。
关键词：数值分析；岩体；渗流特性；强度弱化；有限体积法
岩石中流体不仅影响有效应力，而且可造成岩石部分物质软化或溶解，改变应力、应变的条件，同时又会改变渗流场特征[1―3]。加强岩石介质渗流规律及吸水弱化效应数值技术研究对于水库、堤坝、隧道建设、水下采矿意义重大。
大量文献调研表明[4―9]数值计算方法是进行岩石介质渗流规律及吸水弱化效应研究的有效手段。然而传统方法仅根据吸水弱化实验曲线对岩体的强度及模量进行唯像地折减，未考虑岩体渗流过程、吸水过程及岩土弱化过程的相互耦合作用。目前罕有针对渗流吸水弱化的耦合算法，鉴于此，该文基于连续-非连续数值模拟方法(CDEM)，将最常用的基于连续介质力学的有限元法与基于非连续介质力学的离散元法进行耦合，利用高斯散度定理采用有限体积法，将场量梯度的求解转化为单元外表面场量的积分求解，从而简化了计算过程，且渗流计算时体积守恒将自动满足；并求解水在岩体中的渗流过程及岩体的变形破坏过程，并在数值模型中引入基于单元饱水持续时间的岩体模量及强度弱化本构，实现了岩体吸水弱化动态过程的模拟。
针对粉砂岩石饱水易软化、水对岩石力学性质具有较大影响、饱水状态下岩石强度及弹性模量均大幅度降低的力学性质变化规律，国内外学者做了大量研究工作。Burshtein[10]、Dyke和Dobereiner[11]研究了水对砂岩强度的影响及其水敏感性；赵延林[12]深入研究裂隙水对岩体弹性模量和强度的影响。黄宏伟和车平[13]、张永安等[14]研究了泥岩遇水作用后发生软化的微观机制和水对软化效应的影响。陈栋、文献[15―19]中开展了水对不同类型粉砂岩物理力学性能影响的试验研究，测定了试样的单轴抗压强度、弹性模量等参数，结果都表明，随着吸水时间的增加，试样单轴抗压强度和弹性模量均减小，呈现负指数函数变化规律。综上，该文在数值模拟中选用粉砂岩为研究对象。
岩体中的渗流模型采用各向同性的孔隙渗流模型，渗流过程的求解方式为显式求解。每一时步的求解过程包括单元渗流速度及流量求解、节点总压力求解两个方面[20―23]。具体的求解流程如图1所示。
图1 渗流求解流程
Fig.1 Seepage solving procedure
1.1.1 单元渗流速度及节点流量求解
设岩体中水的输运满足达西定律为：
式中：vi为单元第 i个方向的流体流速；k/(m2·s·Pa−1)为流体的渗透系数；ks为相对渗透系数；P为节点流体的总压力；xj为节点的空间坐标。
ks的表达式为：
为单元的平均饱和度，表达式为：
式中：N为组成单元的节点数量；si为第i个节点的饱和度。
采用高斯散度定理对式(1)进行化简，可表示为：
式中：V为单元体积；M为单元内总面数；
为节点总压力P在第j个面内的平均值；ni为第j个面的单位外法向在 i方向的分量；∆Sj为第 j个面的面积。
基于单元流速，即可计算单元施加给节点i的流量为：
式中：qi为第i个节点获得的流量；Mi为与第i个节点相关的单元面的个数；v为单元流速矢量；nj为与节点i相关的第j个面的单位外法向；∆Sj为第 j个面的面积；Nj为第 j个面的总节点数。
当有多个单元时，需对公共节点处的流速、流量进行叠加，设某节点是Ne个单元的公共节点，则该节点的平均流速及总流量可表述为：
式中：Vi为单元叠加后的节点平均流速；
为第j个单位在该节点处的流速；Q为单元叠加后的节点总流量；qj第j个单元在该节点处的流量。
1.1.2 节点压力的显式求解
如果某节点为压力边界条件施加节点，则节点孔隙水压力即为外界给定的压力。如果该节点不是压力边界条件节点，基于式(7)计算获得的节点流量计算当前时步的节点饱和度s为：
式中：n为孔隙率；Vn为节点总体积；Qapp为流量边界；∆t为计算时步。如果饱和度小于 1，则令该节点的孔隙水压力为0，如果饱和度等于1，则根据式(9)计算节点孔隙水压力pp为：
式中，Kf为流体体积模量。
基于单元的平均饱和度
，单元任意节点的总压力可表述为：
式中：x、y、z为单元某节点全局坐标的3个分量；gx、gy、gz为全局重力加速度的3个分量。
在不考虑水对岩体软化效应的前提下，一般采用考虑应变软化效应的Mohr-Coulomb模型及最大拉应力模型[20―25]描述岩体在外载荷作用下的塑性变形过程及损伤破裂过程。
首先利用增量形式的有限元法计算本时步单元的应力增量为：
式中：∆σij为当前时步的应力增量；∆εij为当前时步的应变增量；∆θ为当前时步的体应变增量；K为体积模量；G为剪切模量；δij为Kronecker记号。
而后计算本时步单元的试探应力为：
式中：σij为本时步的试探应力；σij-o为上一时步的应力。
根据试探应力张量 σij计算当前时步的主应力σ1、σ2及 σ3，根据式(13)判断该应力状态是否已经达到或超过 Mohr-Coulomb准则及最大拉应力准则：
式中：c(t)、φ、σt(t)为块体当前时步的黏聚力、内摩擦角及抗拉强度；Nφ、αP、σP为常数。
如果fs≥0且h≤0，单元发生剪切破坏；如果ft≥0且h＞0，则发生拉伸破坏。
当单元发生剪切破坏时，采用式(15)进行主应力的修正为：
式中，λs、Nψ、α1、α2为常数，其表达式为：
式中，ψ表示剪胀角。
当单元发生拉伸破坏时，采用式(17)进行主应力的修正为：
式中，
将经过Mohr-Coulomb准则及最大拉应力准则修正后的主应力转换至整体坐标系，根据有限元法计算由单元应力贡献出的节点力。
同时，根据当前时步的等效塑性剪应变及等效塑性体应变，对黏聚力及抗拉强度值进行折减，为
式中：c(t +∆t)及 σt(t +∆t )为下一时步的黏聚力及抗拉强度值；∆t为计算时步；c及σt为初始时刻的黏聚力及抗拉强度值；γp及εp为当前时刻等效塑性剪应变及等效塑性体应变；γlim及εlim为剪切断裂应变及拉伸断裂应变。
基于拉剪复合应变软化模型[20]，可以定义3类损伤因子，分别为拉伸损伤因子α，剪切损伤因子β及联合损伤因子χ，为：
岩体中的部分矿物质具有一定的吸水特性，当岩体遇水后，这些矿物质将捕获岩体中的自由水，以物理或化学的方式将这些水存储起来，从而导致岩体中渗流特性的改变，并导致岩体内部微观结构、胶结程度及孔隙特征发生变化，最终导致岩体的宏观弹性模量及强度出现不同程度的降低。
大量的岩石试样吸水实验表明[15―19]，岩石试样的强度、弹模与吸水时间呈现指数下降型关系，具体如图2所示。
图2 粉砂岩试样吸水试验结果[19]
Fig.2 Water absorption experimental results of siltstone samples[19]
当然，通过物理实验只能给出岩石强度、弹模与吸水时间的关系，并不能给出岩石吸水量与吸水时间，以及岩石强度与吸水量之间的关系。因此，该文设岩石中矿物质的吸水满足如下假定：
1) 吸水速率具有先快后慢的指数衰减型特征；
2) 对于特定的岩体，存在一个最大吸水量；
3) 岩体中的渗流压力对吸水速率及最大吸水量无影响。
基于以上假设，可以给出饱水状态下代表性体积单元的吸水量随时间的变化规律：
式中：V为t时刻代表性体积单元的吸水体积；Ve-max为该代表性体积单元的最大吸水体积；t为代表性体积单元的吸水累积时间；∆T为达到63.2%吸水体积所对应的时间。
设代表性体积单元的体积为V0，则式(20)可表述为：
式中：nt为某时刻的吸水率；ne为最大吸水率为：
设ne=10%，∆T=104 s，则可绘制出饱和状态下不同时刻的吸水率随时间的变化，如图3所示。
图3 表征元吸水率随时间的变化
Fig.3 Characterize change of element water absorption with time
岩体吸水后将对表征元的渗流过程发生影响，当前时步的节点饱和度公式s及节点孔隙水压力公式 pp需要调整为[23―25]：
式中，Qs为某时刻单元的吸水速率，可通过对式(20)求导获得：
设单元的体积模量、剪切模量随着单元吸水率的增加线性减小，表达式为：
式中：K′及G′为吸水率是nt时的体积模量及剪切模量；K、G为初始体积模量及剪切模量；KR及GR为达到最大吸水率 ne后的残余体积模量及残余剪切模量。
同理，假设黏聚力、内摩擦角、剪胀角、抗拉强度随着单元吸水率的增加也呈线性减小趋势为：
式中：cR、Rφ、ψR及σtR分别为表征元达到最大吸水率后黏聚力、内摩擦角、剪胀角及抗拉强度的残余值。
将式(21)代入式(26)及式(27)，可获得表征元在任意吸水时刻的模量变化及强度变化为：
设某粉砂岩的密度为2100 kg/m3，弹性模量为500 MPa，泊松比为 0.27，黏聚力为 80 kPa，内摩擦角为22°，剪胀角为10°，抗拉强度为20 kPa，剪切断裂应变为0.3%，拉伸断裂应变为0.1%。
该粉砂岩的孔隙率为 0.5%，渗透系数为5 × 10-13m2·Pa-1·s-1(即 5 × 10-9m/s)，最大吸水率为 10%，特征吸水时间为 104 s，达到最大吸水率后的残余模量是初始模量的0.5倍，残余强度是初始强度的0.1倍。
取直径为5 cm、高为10 cm的试样进行吸水实验，分别观察该试样在低、高两种水压下的渗流过程，试样内含水量、孔隙压力及饱和度分布随时间的变化，试样平均模量及强度的弱化情况。其中，低压值取10 kPa(1 m水深)、高压值取250 kPa (25 m水深)。
考虑模型的对称性，建立1/2圆柱模型，并采用64378个四面体进行剖分(如图4)。渗流计算时，在模型的四周(对称面除外)施加渗流压力边界条件，计算时步为5 ms。
图4 半圆柱体网格部分
Fig.4 Semi-cylindrical mesh generation
2.2.1 低压不同工况下岩石试样吸水演化规律
低压情况下，考虑及不考虑基质吸水特性时的饱和度随时间的演化规律如图5所示，图中浅色表示饱和度为1，深色表示饱和度为0。由图5及图6可得，不考虑吸水时的渗透速度明显快于考虑吸水时的渗透速度；从图5(e)及图6(d)的对比中可以看出，不考虑吸水时试样达到完全饱和的时间比考虑吸水时试样达到完全饱和的时间快30s左右。
图5 低压考虑基质吸水特性的饱和度随时间的变化
Fig.5 Change in the saturation over time at low pressure taking into account matrix’s water absorption characteristic
图6 低压不考虑基质吸水特性的饱和度随时间的变化
Fig.6 Change in saturation over time at low pressure without taking into account matrix’s water absorption characteristics
定义试样中水的总体积含量为：
式中：Vs为自由水体积；Va为基质吸水体积；V0为试样总体积。
则考虑低压吸水及不考虑吸水两种工况下，w随时间的变化规律如图7所示。由图可得，随着渗流时间的增大，试样中水的体积含量逐渐增大。刚开始，两种工况下w随时间的变化规律基本一致；50 s以后，两种工况下的曲线逐渐分开。不考虑基质吸水的工况，由于试样中的水全部来自参与渗流计算的自由水，因此一旦渗流场达到稳定，水的体积含量基本保持不变，约为0.5%。考虑基质吸水的工况，随着渗流时间的增大，一方面外界的水逐渐渗透到试件内部，另一方面，试件内部孔隙中的自由水将被基质捕获，并以某种方式存储在基质中，因此即使渗流场达到稳定，试样中水的体积含量仍然在逐渐增加；当然，当试样达到最大吸水量后，试样中水的总体积含量将保持不变。从图7可以看出，当试样浸入水中的初期，渗流特性占主导，试样中的水大都为自由水；随着时间的推移，渗流效应逐渐减弱，而基质吸水效应将逐渐增强。
图7 低压下试样中水的总体积含量随时间的变化
Fig.7 Total volume of water in sample at low pressure versus time
令饱和单元的体积占比=某时刻达到饱和的单元体积总和/试样体积，则该占比随时间的变化如图8所示。由图可得，不考虑基质吸水时，所有单元均可达到饱和状态；考虑基质吸水时，仅有约90%的体积长期处于饱和状态，且因基质吸水速度与单元内含水量间的动态平衡关系，饱和单元体积占比将在总体平衡的形势下出现局部范围的波动。
图8 低压下饱和单元体积占比随时间的变化
Fig.8 Change of volume fraction of saturated unit with time under low pressure
试样中心点的孔隙压力随时间的变化如图9所示。由图可得，不考虑基质吸水时，渗流稳定后的孔隙压力维持不变，为10 kPa，到达稳态压力的时间约为170 s。考虑基质吸水时，在短期内，无法达到渗流稳定状态，监测点孔隙压力随时间呈周期性变化，变化幅度约为60 kPa。
图9 低压下试样中心点孔隙压力随时间的变化规律
Fig.9 Variation of pore pressure with time at center of sample under low pressure
2.2.2 高水压条件下数值解与理论解对比分析
高压(250 kPa)且考虑粉砂岩吸水的情况下，试样中饱和度的变化规律如图 10所示。由图可得，随着渗流时间的增大，试样迅速达到饱和状态，与图4低压下考虑基质吸水的工况相比，高压吸水工况下试样达到完全饱和的时间仅为8 s，而低压吸水工况下试样达到完全饱和的时间约为200 s。
试样中水的总体积含量随渗流时间的变化如图 11所示。由图可得，渗流刚开始，试样迅速达到饱和状态，试样中的水含量迅速上升至 0.5%(与试样的孔隙率一致)，此部分水为自由水；随着渗流时间的持续增加，试样中的水含量继续增加，但增加趋势逐渐变缓，当渗流时间为 3×104s时，试样中的水的总体积含量已达 9.99%(该试样设置的最大吸水量为10%)，此部分水体积主要由基质吸水贡献。
图10 高水压下考虑基质吸水特性的饱和度随时间的变化
Fig.10 Change of saturation of matrix with time under high water pressure
图11 高压吸水情况下试样中水的总体积含量随渗流时间的变化
Fig.11 Total volume of water in sample under high pressure with water uptake varying with seepage time
扣除图 11中的自由水体积，可得数值计算获得的基质吸水含量随时间的变化曲线如图12。由图可得，数值计算的结果与理论解式(21)趋势基本一致，表明了数值计算的正确性。
图12 高压吸水情况下试样吸水体积含量随渗流时间的变化
Fig.12 Variation of water absorbent volume in case of high pressure water absorption
低压及高压状态下，在渗流初期，试样内水的总体积含量随时间的变化曲线如图 13所示。由图可得，渗流开始阶段，高压下试样中水的体积增加速率明显快于低压下的速率；但到了渗流后期，两种压力下的体积增加速率基本一致(基质吸水占主导)；高压下的渗流过程与基质吸水过程分界线较为明显(曲线上存在转折点)，而低压下两者的分界线不明显。
图13 渗流初期试样中水的总体积含量变化规律
Fig.13 Variation law of total volume of water in initial sample of seepage
对高压(250 kPa)渗流下 500 s、2500 s、7500 s、15000 s及30000 s这5个典型时刻试样的平均黏聚力、内摩擦角进行统计，获得试样平均强度与吸水时间的关系，如图 14所示。由图可得，随着吸水时间的增加，试样的强度(黏聚力、内摩擦角)逐渐减小至残余值；数值解与理论解(式(29))基本一致，表明了数值算法的计算精度。
以某粉砂岩试样的数值模型(图 4)及材料参数为基础，对其进行单轴压缩分析。对数值模型的一端施加法向位移约束，另一端施加垂直于试样的准静态速度载荷，模型的对称面法向约束。基于前面高压吸水的数值结果，分别取吸水时间为0 s、500 s、2500 s、7500 s、15000 s及30000 s时的试样进行数值实验，研究不同吸水时间下试样的峰值强度及破坏模式的变化规律。
图14 黏聚力及内摩擦角随吸水时间的变化
Fig.14 Cohesion and internal friction angle with water absorption time
不同吸水时间下试样单轴压缩的应力-应变曲线如图15所示。由图可得，随着吸水时间的增加，试样的峰值强度、弹性模量均逐渐减小，但减小趋势逐渐变缓；随着吸水时间的增加，试样达到峰值后的理想弹塑性段逐渐增长，但峰值下降段的长度基本不变。
图15 不同吸水时刻试样的单轴压缩应力-应变曲线
Fig.15 Stress-strain curves of uniaxial compressive samples at different water absorption time
对图 15的曲线进行分析，可得试样的单轴抗压强度随吸水时间的变化规律(如图16所示)。由图可得，随着吸水时间的增大，试样的单轴抗压强度逐渐减小，但减小趋势逐渐变缓。干燥试样的单轴抗压强度约为237 kPa，吸水3×104s后试样的单轴抗压强度为24.6 kPa。
对图 16的曲线进行分析，可得试样的弹性模量随吸水时间的变化规律(如图17所示)。由图可得，随着吸水时间的增加，弹性模量逐渐减小，但减小趋势逐渐变缓。干燥试样的弹性模量为500 MPa，吸水3×104s后试样的弹性模量为263 MPa。
图16 试样单轴抗压强度随吸水时间的变化
Fig.16 Uniaxial compressive strength of specimen varying with time of water absorption
图17 试样的弹性模量随吸水时间的变化
Fig.17 Elastic modulus of specimen varying with time of water absorption
不同吸水时间下试样的单轴压缩破坏模式如图 18所示。由图可得，干燥情况下，单轴压缩实验的压剪破坏发生于试件的下部；不同吸水情况下，试件的剪切破坏均发生于试样的上部，仅裂缝的发展形态有所不同。
该文提出了一种渗流吸水诱发岩体强度弱化的有限体积数值计算方法。该方法建立了基质吸水量与吸水时间之间的指数衰减型函数关系，给出了基质吸水引起岩体模量及强度弱化的理论模型及相应的数学表述；通过在各向同性孔隙渗流模型中引入基质吸水弱化算法，实现了渗流作用下孔隙内自由水向基质吸附水转化过程的模拟，并实现了岩体模量与强度随渗流时间逐渐弱化过程的模拟。
(1) 基于考虑吸水弱化算法的各向同性孔隙渗流模型，模拟了低(10 kPa)、高(250 kPa)两种边界流体压力下某粉砂岩试样的渗流过程及吸水弱化过程。数值结果表明，边界流体压力越高，试样达到整体饱和状态的时间越快；渗流初期，以自由水渗流填充孔隙为主，渗流后期，以孔隙内自由水向基质吸水的转化为主。边界流体压力对渗流速率具有明显控制作用，但对基质吸水速度无影响。
图18 不同吸水时刻试样的单轴压缩破坏模式(塑性剪应变云图)
Fig.18 Uniaxial compression failure mode of different samples at different water absorption time(plastic shear strain image)
(2) 基于考虑应变软化的 Mohr-Coulomb模型及最大拉应力模型，模拟了不同吸水时间下某粉砂岩试样的单轴压缩过程。数值结果表明，随着吸水时间的增加，试样的峰值强度、弹性模量均逐渐减小，但减小趋势逐渐变缓；随着吸水时间的增加，试样达到峰值后的理想弹塑性段逐渐增长，但峰值下降段的长度基本不变。不同吸水情况下，试样的剪切破坏均发生于试样的上部，仅裂缝的发展形态有所不同。
该文所述的岩体吸水弱化理论表达式是基于某一特定岩石(粉砂岩)提出的。在后续的研究中，需要深入探讨该吸水方程在其他岩石中的适用性以及用该方法研究工程尺度的问题。
参考文献：
[1] 佘诗刚, 林鹏. 中国岩石工程若干进展与挑战[J]. 岩石力学与工程学报, 2014, 33(3): 433―457.She Shigang, Lin Peng. Some developments challenging issues in rock engineering field in China [J]. Chinese Journal of Rock Mechanics and Engineering, 2014,33(3): 433―457. (in Chinese)
[2] She S G, Dong L J. Statistics and analysis of academic publications for development of rock mechanics in China[J]. Chinese Journal of Rock Mechanics and Engineering,2013, 32(3): 442―464.
[3] Lin P, Zhou Y, Liu H Y, et al. Reinforcement design and stability analysis for large-span tailgated tunnels with irregular geometry [J]. Tunneling and Underground Space Technology, 2013, 38(9): 189―204.
[4] Dafalias Y F, Manzari M T. Simple plasticity sand model accounting for fabric change effects [J]. Journal of Engineering Mechanics, 2004, 130(6): 622―634.
[5] Li X S, Dafalias Y F. Dilatancy for cohesionless soils [J].Géotechnique, 2000, 50(4): 449―460.
[6] Bathe K J, Khoshgoftaar M R. Finite element free surface seepage analysis without mesh iteration [J].International Journal for Numerical and Analytic Methods in Geo-mechanics, 1979, 3(1): 13―22.
[7] Zheng H, Liu D F, Lee C F, et al. A new formulation of Signorinip's type for seepage problems with free surfaces[J]. International Journal for Numerical Methods in Engineering, 2005, 64(1): 1―16.
[8] 李腊梅, 冯春. 一种非连续介质中热传导过程的数值模拟方法[J]. 工程力学, 2016, 33(1): 25―31.Li Lamei, Feng Chun. A numerical silulation method for heat conduction in discontinuous media [J]. Engineering Mechanics, 2016, 33(1): 25―31. (in Chinese)
[9] Hatzor Y H. Fundamentals of discrete element methods for rock engineering: theory and applications [J].International Journal of Rock Mechanics & Mining Sciences, 2008, 45(8): 1536―1537.
[10] Burshtein L S. Effect of moisture on the strength and deformability of sandstone [J]. Journal of Mining Science, 1969, 5(5): 573―576.
[11] Dyke C G, Dobereiner L. Evaluating the strength and deform-ability of sandstone [J]. Quarterly Journal of Engineering Geology and Hydro-geology, 1991, 24(1):123―134.
[12] 赵延林. 裂隙岩体渗流-损伤-断裂锅合的理论研究和工程应用研究[D]. 长沙: 中南大学, 2009: 43―57.Zhao Yanlin. Coupling the of seepage-damage-fracture in fractured rock masses and its application [D]. Changsha:Central South University, 2009: 43―57. (in Chinese)
[13] 黄宏伟, 车平. 泥岩遇水软化微观机理研究[J]. 同济大学学报(自然科学版), 2007, 35(7): 866―870.Huang Hongwei, Che Ping. Research on micromechanism of softening and argillitization of mudstone[J]. Journal of Tongji University (Natural Science), 2007,35(7): 866―870. (in Chinese)
[14] 张永安, 李峰, 陈军. 红层泥岩水岩作用特征研究[J].工程地质学报, 2008, 16(1): 22―26.Zhang Yongan, Li Feng, Chen Jun. Analysis of the interaction between mudstone and water [J]. Journal of Engineering Geology, 2008, 16(1): 22―26. (in Chinese)
[15] 陈栋. 温湿环境下深井软岩强度软化试验研究[D]. 哈尔滨: 黑龙江科技大学, 2014: 80―95.Chen Dong. Experimental study on strength softening of soft rock in deep mine under temperature humidity environment [D]. Harbin: Heilongjiang University of Science and Technology, 2014: 80―95. (in Chinese)
[16] 周莉, 韩朝龙, 梦祥民, 等. 砂岩单面吸水强度软化实验[J]. 黑龙江科技学院学报, 2012, 22(3): 320―324+329+208.Zhou Li, Han Zhaolong, Meng Xiangmin, et al.Softening experiment on single-side absorption strength for sandstone [J]. Journal of Heilongjiang Institute of Science & Technology, 2012, 22(3): 320 324+329+208.(in Chinese)
[17] 赵春雷, 赵成刚, 张卫华. 饱和砂土基于相变状态的不排水本构模型[J]. 工程力学, 2015, 32(12): 68―76.Zhao Chunlei, Zhao Chenggang, Zhang Weihua. An undrained constitutive model of saturated sands based on the phase transformation [J]. Engineering Mechanics,2015, 32(12): 68―76. (in Chinese)
[18] 李铮, 何川, 丁建军, 等. 矿山法城市隧道运营期排水量与水压力关系的一种预判方法[J]. 工程力学, 2017,34(1): 14―21.Li Zheng, He Chuan, Ding Jianjun, et al. A method to predict the relationship between water discharge and pressure during operational period of city tunnels constructed using the mining method [J]. Engineering Mechanics, 2017, 34(1): 14―21. (in Chinese)
[19] 郝耐, 张秀莲, 王淑鹏, 等. 敦煌石窟砂岩吸水特性及力学效应试验研究[J]. 科学技术与工程, 2017, 17(12):21―26.Hao Nai, Zhang Xiulian, Wang Shupeng, et al.Experimental study on water absorption characteristic and mechanical effect of sandstone in Mogao Grottoes,Dunhuang, China [J]. Science Technology and Engineering, 2017, 17(12) : 21―26. (in Chinese)
[20] Itasca Consulting Group Inc. Fast Lagrangian Analysis of Continua in 3 Dimensions, Version 3.0, Users Manual[R]. Minneapolis, Minnesota: Itasca Consulting Group Inc, 2005.
[21] 陆晶晶, 李承亮. 基于 CDEM 的高桩码头在地震作用下破坏模式数值模拟研究[C]// 第十五届中国海洋工程学术讨论会论文集, 2011.Lu Jingjing, Li Chengliang. Numerical Simulation Research on CDEM high pile wharf damage in the earthquake mode [C]// The 15th China Ocean Engineering Symposium Proceedings, 2011. (in Chinese)
[22] 冯春, 李世海, 刘晓宇. 基于颗粒离散元法的连接键应变软化模型及其应用[J]. 力学学报, 2016, 48(1): 76―85.Feng Chun, Li Shihai, Liu Xiaoyu. Particle-Dem based linked bar strain softening model and its application [J].Chinese Journal of Theoretical and Applied Mechanics,2016, 48(1): 76―85. (in Chinese)
[23] 冯春, 李世海, 王理想. 一种基于单元局部坐标系求解二维孔隙渗流问题的数值方法[J]. 岩土力学, 2014,35(2): 584―590.Feng Chun, Li Shihai, Wang Lixiang. A numerical method to solve pore seepage problems based on element local coordinate system [J]. Rock and Soil Mechanics,2014, 35(2): 584―590. (in Chinese)
[24] 李勇. 岩质边坡动力放大系数及近似计算方法的研究[D]. 成都: 西南交通大学, 2013: 36―47.Ling Yong. The study of the dynamic magnification factor of the rock slope and approximate calculation method [D]. Chengdu: Southwest Jiaotong University,2013: 36―47. (in Chinese)
[25] Feng C, Li S H, Liu X Y, et al. A semi-spring and semi-edge combined contact model in CDEM and its application to analysis of Jiweishan landslide [J]. Journal of Rock Mechanics and Geotechnical Engineering, 2014,6(1): 26―35.
A FINITE VOLUME NUMERICAL SIMULATION METHOD FOR ROCK MASS STRENGTH WEAKENING BY SEEPAGE WATER ABSORBING
GUO Ying1,2, JIANG Xin-liang1, CAO Dong-bo3, BAI Tie-jun2, ZHU Guang-yi2, FENG Chun4
(1. School of Civil Engineering, Tianjin University, Tianjin 300072, China; 2. School of Civil Engineering, Shenyang University, Shenyang 110044, China;3. Public English Department, Shenyang Institute of Engineering, Shenyang 110136, China;4. Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China)
Abstract:It is difficult to describe the deterioration mechanism and permeability evolution mechanism of rock materials by the traditional numerical method of seepage flow. A finite volume numerical method is proposed for the weakening strength of rock mass induced by seepage water absorption. The seepage process of water and the deformation and destruction process of rock mass are solved with the finite volume method by using the Gauss divergence theorem, and a theoretical model for modulus and strength of the weakening rock mass induced by the matrix suction is established and the corresponding mathematical expressions are given. Based on the isotropic pore seepage model considering a water weakening algorithm, the softening process of a powder sandstone sample and the uniaxial compression process of specimens with different water absorption time are simulated under both low and high boundary fluid pressures. The numerical examples show that the higher the pressure of boundary fluid, the faster the sample reaches the overall saturation state. In the early stage of seepage, the free water flow fills the pores, while it is dominated by the transformation from the free water in the pore to the matrix suction in the later stage. The boundary fluid pressure has obvious control effect on the flow rate, but has no effect on the water absorption speed of the matrix suction. With the increase of the absorption time, the strength(cohesion and internal friction angle) of the sample is gradually reduced to the residual value, and the matrix water absorption content obtained by using this method is basically consistent with the theoretical solution. This demonstrates the calculation precision of the numerical algorithm and that it can be used to analyze the seepage and stress coupling effect on actual rock mass engineering problems such as tunnel water breakthrough and surrounding rock stability.
Key words:numerical analysis; rock mass; seepage characteristic; strength weakening; finite volume method
中图分类号：TU452
文献标志码：A
doi:10.6052/j.issn.1000-4750.2017.03.0209
文章编号：1000-4750(2018)07-0139-11
收稿日期：2017-03-17；修改日期：2017-10-24
基金项目：国家青年科学基金项目(51208356)；辽宁省自然科学基金指导计划项目(20170540649，20170540651)；国家留学基金项目(201708210323)
通讯作者：郭 影(1979―)，女，辽宁人，副教授，博士，主要从事地下结构和岩土工程方面的研究(E-mail: guoying0829@163.com).
作者简介：
姜忻良(1951―)，男，浙江人，教授，博士，博导，主要从事地下工程及结构与土相互作用的研究(E-mail: jiangxinliang@126.com);
曹东波(1980―)，女，辽宁人，讲师，硕士，主要从事外语教学和数值软件翻译研究(E-mail: 2352052703@qq.com);
白铁钧(1962―)，男，辽宁人，教授，博士，主要从事工程力学方面研究(E-mail: tiejun1009@126.com);
朱广轶(1962―)，男，辽宁人，教授，硕士，主要从事开采沉陷与治理方面的研究(E-mail: yxj2038@163.com);
冯 春(1982―)，男，浙江人，助理研究员，硕士，主要从事岩土力学领域数值方法的研究(E-mail: fengchun@imech.ac.cn).