热－流－固耦合岩体三轴压缩实验数值模拟

杨长德，王 鹏，2，毛金峰，李金波，张海东

（1.新疆工程学院 矿业工程与地质学院，新疆 乌鲁木齐830091；2.中国矿业大学（北京） 能源与矿业学院，北京100083；

3.山西工程技术学院 土木与建筑工程系,山西 阳泉045000）

地下煤岩体在渗流力、地应力、开采扰动及温度等各种不同组合的作用下，其内部结构与内部的物质含量会发生相应的改变，进而会导致其力学性能和渗透特性发生改变，这对煤炭的安全开采或是对矿井的稳定及维护都会产生重大的影响[1-2]。目前国内外学者对多场耦合下的煤岩体进行了大量的实验研究和数值模拟研究，并取得了较为丰富的研究成果。在多场耦合实验研究方面，谢和平等[3]采用煤岩采动卸压增透试验平台研究了不同开采方式扰动下煤岩应力场-裂隙场-渗流场的差异性行为，对指导煤与瓦斯共采提供了一定的理论依据；江宗斌[4]等进行轴压循环加卸载条件下的岩石应力-渗流耦合试验，分析了岩石加卸载过程中的渗透性变化以及力学行为变化规律；杨金保[5]等开展了单裂隙花岗岩不同围压加、卸载和不同水力梯度作用下的多场耦合试验，得出了应力历史是影响裂隙岩体渗透特性的因素之一；陈天宇[6]等对含气页岩试件进行围压循环加卸载，渗透压固定，无轴压施加方式，得到了围压敏感性和各向异性对含气页岩变形以及渗透特性有重大影响；曹亚军[7]等开展了不同围压和渗压作用下渗流-应力耦合三轴流变试验。在多场耦合数值模拟方面，赵延林[8]等建立了裂隙岩体温度场-渗流场-应力场耦合的双重介质模型，认为在一定时域下，存在起主导控制作用的双场耦合系统规律；张树光[9]等通过对裂隙岩体流-热耦合传热过程进行三维数值模拟分析，得出了岩体渗透系数的变化对其内部温度场的分布影响很大；刘泉声[10]等通过开展多场耦合作用下岩体裂隙扩展演化规律的研究，得出THM 耦合及岩体变形、失稳全过程的数值模拟算法；康永水等[11]、冯梅梅[12]等、高娟[13]等对裂隙岩体冻融损伤特性以及多场耦合作用过程中岩体的应力-应变特征、裂隙发育规律进行了研究；王军祥等[14-15]建立数值模型，利用反演的参数对隧道围岩应力场、渗流场、损伤场分布规律及衬砌结构的受力特征进行了分析。通过以上研究成果可以看出，对煤岩体在多场耦合条件下的研究逐渐在形成一种常态。

随着煤炭深部开采的日益增多，于“三高一扰动”条件下的深部煤岩体的力学作用机理和破坏失稳规律已不再适应于浅部煤岩体的发展力学规律[16-18]。为此，对温度场、应力场、渗流场三场耦合条件下的三轴压缩实验进行数值模拟研究，进而获得多场耦合条件下煤岩体的力学行为发展规律和破坏机理。

1 三场耦合下的数学模型

由于煤岩体自身具有不均质性，各向异性以及不可再生性，因此，实验室中的每1 个煤岩体试件内部成分、均质分布、结构组成以及基本力学参数都不可能完全相同，即实验室中的每1 个煤岩体试件都是不同的，并无复制性。目前体现煤岩体试件可复制性的研究手段有参数反演进行三维重构或者采用3D 打印技术，和CT 扫描进行三维重构[19]，再或者进行数值模拟研究。采用数值模拟研究建立渗流场-应力场-温度场耦合的数学模型。

建立的三场耦合模型基于如下假设：岩体属于均质连续介质体；渗流规律服从Darcy 定律；水的汽化不考虑，岩体全过程处于饱水状态；热质在固相、液相介质中传递方式以对流传热为主，各介质的比热容及热传导系数不随温度和压力而变化。

1）流体渗流控制方程。

式中：k 为渗透系数；μ 为动力黏度；▽p 为出入口速度差形成的渗透差；ρ1为流体密度；Qm为整个渗透过程流体流量。

2）温度控制方程。

式中：▽T 为岩体传热过程形成的温差；kf为岩体导热系数；Cp为比热容；u 为流体流动速度；ρ2为岩体密度；Q 为岩体本身缩传热量；Qted为流体吸收热量。

3）岩体应力控制方程。

式中：S1、S2分别为岩体轴向压力作用面面积和围压作用面面积；Fv1、Fv2分别为轴压和围压；R 为岩体半径；σz为轴向应力。

4）岩体温度-应力-渗流耦合方程。

式中：T－ext为岩体点空间温度；▽Fv为岩体耦合过程中作用力差；S 为对应力作用面积。

联立以上控制方程，得到岩体三场耦合下的数学模型，结合相应条件，进行模型的求解。

2 网格划分和边界条件

选取模型尺寸为25 mm(直径）×100 mm（高），采用3 节点三棱柱体单元进行网格剖分，共剖分为94 个节点，156 个三角形单元。实体单元网格剖分如图1。

图1 实体单元网格剖分图Fig.1 Mesh generation diagram of entity cells

上边界由温度梯度传热计算出热-流边界条件，下边界初始温度为75 ℃，外部温度为25 ℃，四周为对流传热边界。上边界为渗流入口，流速0.009 m/s，下边界为出水口，流速为0.005 m/s。左右两边界为封闭状态，无水流速，整个过程选值分别为1×10-10、1×10-9m2、1×10-8m2的渗透系数进行模拟。右边界为围压4、5、6 MPa 应力载荷边界，上边界为轴压载荷和指定位移边界，轴压载荷边界通过指定位移产生的应力积分函数定义，下边界为固定约束边界。计算中所用到的计算参数和数据见表1。

表1 计算模拟参数表Table 1 Calculation of simulation parameters table

3 数值模拟结果与分析

采用有限元软件对数学模型进行求解[20]，渗流-应力-热耦合下的岩体试件温度等温线分布、表面温度分布、应力分布以及渗流速度如图2～图5。

图2 岩体温度等值线图Fig.2 Temperature contour map of rock mass

图3 岩体温度表面图Fig.3 Temperature surface of rock mass

图4 应力分布图Fig.4 Stress distribution diagram

图5 渗流速度表面图Fig.5 Surface diagram of seepage velocity

从图2 和图3 可以看出，在热-应力-渗流三场耦合作用下，岩体热传递并非是均匀传递的，温度变化由岩体底部向上部非线性传递，因而岩体各局部位置温度不一，但整体温降趋势比较显著。同时由温度最大值与最小值可以看出，由于渗流场的介入，整个多场控制下的耦合系统存在热量散失，这主要是因为渗透液吸收了岩体部分热量所致。

从图4 可以看出，在热-应力-渗流三场耦合控制系统下，岩体应力分布呈现端部应力集中，但其应力分布是呈对称分布的。这也就表明了在多场耦合条件下，岩体试件实验存在场耦合作用的显著影响以及明显的端部效应现象。

从图5 可以看出，在多场耦合控制机制下，由于受到应力作用，岩体内部孔隙结构发生改变，即渗流通道发生改变，从而形成了连续均质岩体的非均质渗流，并且在岩体内部产生了流速阀值，且已远远超过了初始渗流速度。

3.1 渗流场分析

由不同围压梯度下岩体各指定位移节点下渗透速率表面图如图6。可知，在应力场和温度场耦合作用下，原本均质岩体均匀流动转变成非均匀流动，流体在岩体内流速阀值受围压影响而不同。围压越大，流体在岩体中的流速阀值越低。流体渗流速度增值按围压值升高梯度增值成正向比例增加。围压越高，岩体内流体流速阀值点数越多，而且在各围压梯度下的流体速度阀值均远超过流体的初始速度，造成这现象的原因可能是三轴应力下的渗透岩体内部均匀结构发生了改变，岩体内部产生了新生裂隙，进而贯通了岩体内部的新增渗流通道，从而改变了流经岩体的渗透速率，形成了以下岩体渗流场的速度分布状态。由此可以看出作用于岩体的应力场对流经岩体的渗流场产生的影响是比较明显的。

图6 指定位移节点下渗透速率表面图Fig.6 Surface diagram of penetration rate under specified displacement node

3.2 温度场分析

通过渗流场、温度场、应力场三场耦合下的三轴压缩实验进行模拟，在模型中心位置每10 mm 截取1个点，共计11 个空间点，获得了在不同围压，不同渗透系数下的温度数据列，不同渗透系数下温度变化图如图7，不同围压梯度下温度变化图如图8。

图7 不同渗透系数下温度变化图Fig.7 Temperature changes under different permeability coefficients

图8 不同围压梯度下温度变化图Fig.8 Temperature changes under different confining pressure gradients

由图7 和图8 可以看出在同一围压不同渗透系数下，由于水的渗流作用，各工况下的岩体温度均发生非线性降低，而且岩体的渗透系数与岩体温降梯度成正向比例关系。另外，作用于岩体的围压越大，岩体温降梯度越大。而在同一渗透系数不同围压梯度下岩体的降温梯度是相差不大的，但其温度降低趋势均是一致的。尽管岩体在应力场和渗流场综合作用下温度产生了非线性降低，但岩体的最低温度还是高于外界温度，这表明岩体内部温度场传热产生了一定的耦合作用。综合上述分析表明，渗流场和应力场的耦合作用会影响岩体温度场的分布，同时温度场对渗流场和应力场也会起到一定的附和耦合作用。从以上分析可以得出渗流场对温度场的影响是起主导作用的，而应力场对温度场的影响作用并非十分显著。

3.3 应力场分析

指定点应力表面图如图9。

由图9 可知，岩体在渗流场和温度场以及应力场三场耦合作用下，岩体应力分布呈现出一定规律：

1）岩体在三轴压缩过程中，端部效应显著，边缘应力分布集中，这意味着在实验室做岩体三轴实验时，应注意采取一定措施排除端部效应影响。

图9 指定点应力（ pc=5 MPa, K=1×10-3 m2）Fig.9 Stress surface at specified point（ pc=5 MPa,K=1×10-3 m2）

2）岩体在多场耦合下进行三轴压缩实验，岩体应力分布呈现上下对称分布状态。

3）应力场中轴压和围压以及渗流场和温度场综合作用下，岩体内部应力相互抵消，进而形成一定的有效应力场。形成的有效应力场才是造成岩体产生形变的主要原因。

4 结 论

1）围压是影响岩体渗透性的因素之一。在其他各场条件一致的工况下，围压越大，煤岩体的渗透性越低。

2）作用于岩体的应力场对流经岩体的渗流场产生的影响比较明显。而流经岩体的渗流场对岩体内部分布的温度场产生的影响占主导作用，但作用于岩体的应力场对岩体内部分布的温度场产生的影响并不十分显著。

3）岩体在多场耦合条件下进行三轴压缩实验模拟，端部效应显著，边缘应力分布集中，岩体应力分布呈现上下对称分布状态。