高原多年含水冻岩爆破炮孔温度-水分场耦合分析

费鸿禄，甄 帅，李文焱，胡 刚，聂 寒

(辽宁工程技术大学 爆破技术研究院,阜新 123000)

高原高寒地区的露天矿山进行爆破作业时，爆破区域炮孔成孔后，渗流水通过岩石裂隙流入孔内，由于受当地海拔高度、气候、昼夜温差、水文地质的多种因素影响，经过一定时间导致炮孔内渗流水完全冻结，从而无法进行装药工序。研究表明，影响渗流水冻结的主要因素分别为场区温度场作用以及渗流水速度，相比而言，温度场所引起的效应更大[1,2]。

国内外许多学者对冻土温度场已经有了较为深入的研究，Noorishad等基于Biot固结理论[3]，将变分原理和伽辽金公式与有限元方法相结合，首次提出了饱和岩体固-液-热耦合方程。此后，胡向东等通过单排管、双排管、三排管对冻结壁温度影响的一系列研究[4-9]，提出“卸载状态下冻结壁-周围土体共同作用”概念并建立相关冻结壁力学模型，推导出冻结稳态温度场解析解矩阵表达式。王效宾等以南京地铁线工程为例[10]，采用ADINA有限元软件分析地铁线温度场，得出了导热系数、比热容等因素对冻土温度场的影响。张松等为研究突发渗流作用下温度场温度演变过程基于相似模型开展冻结试验并得到相关规律[11,12]。邵玉龙等建立了描述三维裂隙网络岩体渗流传热耦合模型[13]，并与二维裂隙传热模型做对比，发现大量流体流经基岩时会通过裂隙同周围岩体进行热量交换，造成裂隙周围岩体温度首先发生改变，然后向周围扩散的规律。黄峰等通过有限元软件ANSYS和TAITHERM建立了武汉某立交桥的三维温度场计算模型[14-16]，数值模拟结果与实时监测数据耦合度颇高，为类似工程温度场研究提供参考依据。邰博文等通过COMSOL有限元模拟及现场试验揭示了寒区铁路路基的水热耦合机理[17-20]，并分析其与时空演变规律。韩小妹等在高寒地区通过温控监测数据和研究分析成果[21]，总结出混凝土坝温度场和应力场的变化规律。周晓敏等在矿井工程中采用数值模拟和现场试验对渗流作用下冻结壁的温度场发展规律进行研究[22,23]。Yuchao Zeng等研究了地热田温度场、地下水初始静压力场和水密度分布特征[24]，同时分析影响温度场的主要因素，并考虑了裂缝系统的地质特征，建立了裂缝系统的概念模型。

综上所述，目前国内外学者在分析温度场-水分场时，大多是基于ANSYS、ADINA等有限元软件或现场数据及理论公式进行单一变量分析并得出相应结论，并且主要针对隧道冻结壁、深井工程、大体积混凝土等大深度工程进行温度场研究，对高寒常年冻岩地区露天矿山爆破作业过程中解决炮孔出现涌水、冰冻现象对温度-水分场多场耦合做深入系统研究的鲜见。因此，本文以拉萨市墨竹工卡县甲玛矿为背景，以光纤测温仪的试验结果与COMSOL-Multiphysics有限元模拟相结合，对比分析高原地区炮孔成孔后温度场发展规律。通过研究温度场的时空关系，从根本上解决常年冻岩条件下炮孔冰冻装药的技术难题，为高原露天矿高效开采提供有效的技术方法。

1 工程背景

甲玛露天矿位于西藏拉萨市墨竹工卡县，年度平均气温2～17℃，高寒干燥，昼夜温差大，在高寒高海拔条件下，最低温度可达-10℃。在进行爆破作业时，导致爆破场地内炮孔常年出现地表向孔内渗水和钻孔时孔壁涌(渗)水、冰冻现象。主要表现为：在炮孔成孔后3～12 h孔内渗流水经过矿区温度场的重分布在短时间内可完全冻结，导致大量完成的炮孔无法实施装药工序，当前应用最广泛的解决措施是选择传统的钻机透孔法进行破冰处理，但此方法在实际应用中存在多个不足，首先此方法是通过压缩透孔时间以满足紧后工序的装药步骤，这样就有可能出现交叉作业的风险增加，造成安全隐患。其次二次透孔的成本也是不可忽略的，增加透孔工序会严重制约爆破作业的效率。

因此研究炮孔成孔后温度场变化规律，并根据研究成果开发相对应的技术、装备、方法，对爆破工序的正常进行具有重要意义。为研究高原多年含水冻岩温度场的发展规律，以拉萨市墨竹工卡县甲玛露天矿5120平台为工程试验背景，如图1所示。

图 1 甲玛露天矿Fig. 1 Jiama open-pit mine

2 耦合理论分析

2.1 温度场控制方程

温度场是指某研究物体在某时间内所处空间温度的集合，它反映了温度的时空间关系，通常用空间直角坐标系以及时间变量表示,其方程如下所示

T=f(x,y,z,t)

(1)

式中：T表示温度；x、y、z为空间直角坐标系坐标；t为时间。

矿区炮孔内温度场为三维热传导问题，但由于炮孔成孔后为空心圆柱，因此可取对称结构温度变化，简化为二维热传导问题。如(2)所示

T=f(x,z,t)

(2)

假设岩体为各向同性且均匀的，通过数学物理方法中热传导方程拉普拉斯公式，基于弹性力学理论推导出其二维热传导微分方程

(3)

(4)

式中：α为常数；Q表示单位热量，J；c表示物体比热容，J/(kg·K)；ρ表示物体密度，kg/m3；2为拉普拉斯算子。

考虑相变潜热的传导方程

(5)

式中：C为热容量，J/(kg·K)；θ为体积含水量，kg/m3；λ为导热系数，W/(m·K)；对于简化二维问题为为相变潜热；ρI为冰的密度，kg/m3；θI为冰的体积含量，kg/m3。

在进行有限元软件COMSOL-Multiphysics数值建模时，炮孔内渗流水由液相转变为固相需要放热，其释放的热量会通过岩石孔壁及孔口空气流动发生热传导。

2.2 水分场控制方程

在炮孔成孔后矿区原地质构造破坏会产生裂隙，裂隙中存在的水通过裂隙渗流至炮孔，其渗流作用符合达西定律。根据Richard方程并考虑冰与水的相变可得非饱和岩石内水分场迁移方程

(6)

式中：ρW表示水的密度，kg/m3；θu为液态水体积量，kg/m3；k为渗透系数。冰冻岩土中水的扩散率如下

(7)

式中：k(θu)为土体渗透率，m/d；c(θu)为比水容量，kg；I为阻抗因子[25]，表示孔隙中冰对水产生渗流的阻碍作用。

2.3 温度场与水分场的耦合

温度场控制方程中，温度T为自变量。水分场控制方程中，液态水体积含量θu与冰体积含量θI为自变量。在有限元软件建模求解时，两个控制方程不足以将方程中三个未知变量表示，需要引入耦合项将水热方程联系求解，因此选用“固液比”的概念作为耦合项[26]。固液比即为岩土中冰的体积含量与自由水体积含量的比值，其公式可表示为

(8)

式中：Tf表示岩土体冻结温度，K；B为常数。

因此冰的体积含量为

θI=BI·θu

(9)

由岩土体的相对饱和度公式，结合VG滞水模型及Gardner渗透模型[27]，可推导出水的体积含量

(10)

θu=(θs-θr)·S+θr

(11)

式中：s为相对饱和度；θs为饱和含水率；θr为残余含水率。因此温度场控制方程中热源项可化简为

(12)

引用“固液比”的概念后，温度场控制方程为

(13)

水分场控制方程为

[D(S)S+k(S)]

(14)

此时可基于COMSOL-Multiphysics有限元软件建立温度-水分耦合模型。

3 数值模型建立

3.1 选择物理场

为分析矿区炮孔涌水冻结过程及炮孔内温度场发展规律，考虑到温度场-水分场的相互耦合作用，在COMSOL-Multiphysics的模型库中选择多孔介质传热模块，并添加流体域及固体域进行有限元分析。为以下叙述方便，称炮孔成孔后至装药前之间的时间段称之为空孔期，结合本工程背景的实际，空孔期孔内温度发展规律对于涌水冻结具有决定性影响作用，现场实际空孔期在24 h左右，因此对温度场进行装药时间最大间隔24 h内的瞬态研究。

3.2 建立几何模型

使用COMSOL-Multiphysics组件中几何模型建立板块。在空孔期中孔内温度主要受到自由水在炮孔壁裂隙的渗流作用影响，同时，炮孔内的温度场通过空气对流以及孔内与地表面的热交换效应，也会对炮孔内的温度变化造成一定影响。简化炮孔成孔后孔内y方向温度变化，建立相应二维几何模型，模型总分为2部分，中间部分为直径0.12 m、深度18 m炮孔，炮孔周边为待爆矿区，左右区域各设5 m影响范围，如图2所示。

图 2 矿区几何模型Fig. 2 The geometric model of the mining area

3.3 模型边界条件的确定

(1)力学边界条件的确定

模型两边为辊支撑，底面为固定约束，顶面为自由约束。

(2)温度边界条件的确定

炮孔成孔后温度主要受炮孔周边岩体的热传导作用影响，使用SG-DTS-84U光纤测温仪进行矿区稳定状态下孔内温度测量作为炮孔成孔后的初始温度。取不同地表温度分析初始地表温度对炮孔温度发展的影响，并以外部自然对流的方式向矿区底部传热，矿区底面平均温度为炮孔底部温度。

(3)水分场边界条件的确定

炮孔在成孔后裂隙水会随着孔壁流入孔内，由于渗流速度在不同平台下各不相同，因此分析不同渗流速度下裂隙水对炮孔内温度场的影响，并给予流体添加相变材料，由液态转变为固态[28]，相变温度273.15[K],相变转变间隔2[K]，相变潜热为333[kJ/kg]。见图3。

图 3 相变潜热Fig. 3 Latent heat of phase change

3.4 模型材料定义

模型材料分为流体、固体、冰3种。具体属性如表1所示。

表 1 模型材料

3.5 划分网格

使用COMSOL-Multiphysics内置Mesh对模型网格划分，全局划分定义超细化，矿区边界部分定义为极细化，炮孔内部定制单元大小参数，最大面积单元0.06 m2，最小面积单元大小0.0004 m2。总计算单元16392个，计算网格面积202.4 m2。如图4所示。

图 4 几何模型网格划分Fig. 4 Meshing of geometric model

4 模拟结果与分析

拉萨市墨竹工卡县冬季平均气温为-3～13℃，夏季平均气温6～20℃，由于甲玛矿区处在高寒高海拔地区，因此温度相对于墨竹工卡县更低。考虑矿区地表区域温度受外部环境影响大，炮孔成孔后渗流水的流速对温度的扩散起决定性作用，因此本模型分别对外界温度为0℃(273.15 K)、5℃(278.15 K)、10℃(283.15 K)渗流水在炮孔内的涌出速度1 m/h、2 m/h、3 m/h、5 m/h、10 m/h的情况进行模拟分析，当渗流水溢满炮孔时水位达到相对稳定状态不再上升。使用COMSOL-Multiphysics有限元分别分析成孔后矿区温度场变化情况，并导出各个时间段的温度数据。

4.1 炮孔成孔后24 h内温度分析

结合现场工程情况，对地表温度为5℃(278.15 K)炮孔内渗流速度为3 m/h时成孔后24 h内温度发展规律进行模拟分析。t=0时矿区温度分布如图5所示，渗流水速度3 m/h时炮孔成孔后1 h、2 h、4 h、8 h、12 h、16 h、20 h、24 h孔内温度场随时间分布图如图6所示。通过有限元软件模拟结果，对比分析炮孔成孔后24 h内温度随深度变化规律，如图7所示。

图 5 成孔后矿区初始温度分布图Fig. 5 Initial temperature distribution in the mining area after blast hole formation

图 6 渗流水速度3 m/h温度随时间变化图Fig. 6 The temperature changes with time with the seepage water velocity of 3 m/h

由图5可见，炮孔成孔后t=0时矿区温度分布均匀，且温度随深度递减。

由图6模拟结果显示：在炮孔成孔后，孔内温度与矿区内部温度随时间变化明显，炮孔内0～8 h温度变化快，8 h后逐步趋于稳定，其相同深度的温度差异不大。炮孔深度0～5 m区域温度受对流热通量影响大，在深度3～6 m时呈现出温度分层。深度5～18 m区域温度变化由炮孔周边岩体温度决定，温度随深度整体呈下降趋势。

由图7可见：成孔1 h内与炮孔深度5 m左右位置存在明显温度拐点，且温度拐点的深度随时间的增加逐渐降低。当成孔12 h时孔内温度已达到0℃(273.15 K)以下，此时可认为炮孔内部完全冰冻。在成孔2 h内深度在6m以内的炮孔温度下降速率大，成孔2～24 h深度在2 m以内炮孔温度变化剧烈。这是因为随着时间的增加，外界地表温度对炮孔内的温度场影响逐渐减小，孔内渗流水因温度的降低逐渐冻结为冰相，冰相冻结体的导热系数大于空气介质导热系数，因此炮孔底部的低温会迅速向孔外传导，但由于炮孔顶部的地表温度相对固定，所以造成在孔深6 m以内温度显著变化。

图 7 炮孔成孔24 h内温度随深度变化图Fig. 7 Temperature change with depth within 24 hours of blasthole formation

4.2 不同渗流水速度作用下炮孔内温度分析

为进一步探讨不同渗流水速度作用下炮孔内温度的发展规律，在此模型上改变渗流水速度，分析流速为1 m/h、2 m/h、5 m/h、10 m/h时炮孔成孔3 h、8 h、12 h、24 h后孔内温度变化，如图8所示。

图 8 不同渗流速度作用下温度分布Fig. 8 Temperature distribution under different seepage velocities

由图8可见：①孔内渗流水速度对炮孔内24 h后的整体温度分布有较大影响。在一定流速范围内，渗流水速度越快，炮孔内整体温度越低。由于矿区底部温度低且较为恒定，因此渗流水速度主要影响深度较浅的孔区温度，对于深度较深的炮孔温度影响不大。②随着炮孔内温度的降低，炮孔周边岩体会随渗流水的热传导作用温度也相应降低。③从二维角度分析，炮孔成孔后不同速度的渗流水对矿区表面初始温度的影响范围≤2 m，对炮孔周边岩体的影响范围≤3 m，渗流速度越快，影响范围越大。

分析不同渗流速度作用下孔内温度随时间变化规律及24 h时孔内温度，不同渗流水速度作用下温度对比如图9所示，不同渗流水速度作用下24 h时孔内温度如表2所示。

图 9 不同渗流水速度作用下温度对比Fig. 9 Temperature comparison under different seepage water velocities

表 2 不同渗流水速度作用下24 h时孔内温度

由图9可见：①在炮孔成孔后的任意时刻，渗流水速度越快，孔内温度下降速度越快。②在不同渗流速度作用下，渗流速度越慢炮孔顶部温度下降幅度越大，反之越小。当渗流速度1 m/h时，炮孔顶端温度在24h内由277.92 K降为271.22 K，温度下降幅度6.70K。渗流速度10 m/h时，炮孔顶端温度在24 h内由267.16 K降为263.10 K，温度下降幅度4.06 K。

由表2可见：①当渗流水为1 m/h时，孔内整体温度都已达到0℃(273.15 K)以下，此温度极有可能导致孔内裂隙水产生冰冻。②炮孔内最低温度位于炮孔最深处且温度恒定，当孔内裂隙渗流水速度高于5 m/h时24 h后炮孔内最大温差将小于3 K，平均温度为262.66 K，此时渗流水流速的增加对炮孔温度变化影响不明显。

4.3 不同初始地表温度状态下炮孔内温度分析

由于不同地表温度对炮孔深度较浅区域温度影响较大，因此研究初始地表温度对炮孔温度的影响范围及效果。依次改变地表初始温度为0℃(273.15 K)及10℃(283.15 K)，对裂隙水渗流速度为3 m/h的场区温度进行模拟，如图10所示。对相同初始地表温度状态下不同时间炮孔内温度变化进行对比，如图11、图12所示。

图 10 不同初始地表温度对炮孔温度影响Fig. 10 The influence of different initial surface temperature on blasthole temperature

图 11 地表初始温度0℃(273.15 K)孔内温度变化Fig. 11 The initial temperature of the ground surface is 0℃(273.15 K) and the temperature change in the hole

图 12 地表初始温度10℃(283.15 K)孔内温度变化Fig. 12 The initial temperature of the ground surface is 10℃(283.15 K) and the temperature change in the hole

由图10可见：①随着初始地表温度的上升，矿区深度在6 m内区域温度明显升高，深度6 m以上区域温度变化不明显，即初始地表温度对矿区温度影响限一定范围。②从地表初始温度为0℃(273.15 K)及10℃(283.15 K)温度影响范围对比可发现，地表温度的升高会减弱渗流水对矿区温度的影响范围。

由图11及图12可见：①当矿区地表初始温度为0℃(273.15 K)时，炮孔成孔3 h时整体温度可达到零下，孔内平均温度265.39 K，距离炮孔深度4 m以外时温度变化相对于4 m以内较为明显。炮孔成孔8 h后孔内整体温度变化减缓，最高温度267.39 K，孔内平均温度264.44 K。12 h、24 h时炮孔平均温度分别为264.14 K、263.74 K，温度达到相对稳定值。当矿区地表初始温度为10℃时，矿区地表温度与矿区底部温度梯度大，因此地表传热现象较为显著。②从炮孔成孔后的4个时间节点可发现，随着时间的增加，深度0～6 m处孔内温度受外界环境影响大，孔深6 m以下区域温度基本保持稳定。

5 现场试验与模拟对比分析

试验地点位于拉萨市墨竹工卡县甲玛露天矿海拔高度5120 m，天气阴，温度8～20℃，矿区顶平均气温6℃，最低温度-2℃。

5.1 试验仪器及使用

分布式光纤测温系统(DTS)也称为光纤测温，主要原理以及喇曼(Raman)散射效应和光时域反射(OTDR)对光纤温度的感度实现温度监测。本试验采用测温仪器为SG-DTS-84U光纤测温系统，测温系统需使用熔纤机将分布式光纤与SG-DTS-84U光纤测温系统串联，待测前利用红光笔检测光纤通路，确认联通后方可进行现场试验。

5.2 试验结果

通过炮孔成孔后24 h内现场试验结果与模拟对比分析。SG-DTS-84U光纤测温系统如图13所示，炮孔内部图如图14所示，孔内渗流水溢出图如图15所示。

图 13 SG-DTS-84U光纤测温系统Fig. 13 SG-DTS-84U Fiber Optic Temperature Measurement System

图 14 炮孔内部图Fig. 14 Internal view of the blasthole

图 15 孔内渗流水溢出Fig. 15 Seepage water overflow in the hole

通过现场试验及COMSOL-Multiphysics有限元数值模拟，得到24 h内不同时间段炮孔内温度变化情况，对比图如图16所示，孔内温度误差如表3所示。

图 16 24 h内不同时间段炮孔内温度对比图Fig. 16 Comparison of the temperature in the blasthole at different time periods in 24 h

表 3 孔内温度误差

由图16可见：①炮孔成孔后，最低温度一般为炮孔底部，孔底低温通过孔内渗流水逐渐向炮孔顶部扩散，在成孔2 h时0℃深度约为4 m，成孔8 h左右0℃深度约为2～3 m，成孔16 h左右0℃深度约为2 m，成孔24 h后孔内温度可达到0℃(273.15 K)以下。②炮孔成孔后孔内温度在8 h后下降速率明减慢并逐渐趋于稳定。

由表3可见：①实测温度与模拟温度在深度为0～4 m时温差较大，最大温差可达5 K，深度为4～18 m时误差较小，一般不会高于2 K。这是因为在此深度的孔温受外界环境影响较大，且实验地点为高海拔地区，气温并不稳定，而炮孔内部温度相对恒定，受外界因素干扰较小，因此接近地平面处孔温会发生较大变化。

6 结论

以拉萨市墨竹工卡县甲玛矿区为研究对象，使用有限元模拟及现场试验分析了高原高寒地区炮孔成孔后孔内温度场变化规律，所得结论如下：

1)炮孔成孔后温度变化主要受地层温度影响，导致孔内积水冻结，且炮孔温度随深度整体呈下降趋势，初始地表温度对炮孔上部的影响趋于一定范围(≤6 m)，与炮孔上部(0～5 m)区域在受空气对流热通量耦合作用影响致使孔内积水未达到冻结临界状态。

2)孔内24 h后温度变化主要由渗流水速度决定，但当流速达到温度影响流速阈值(5 m/h)时，渗流水速度对炮孔内温度变化影不再显著。

3)炮孔成孔8 h后温度场会达到相对稳定值，装药更为合理，若因特殊原因无法在此时间内装药，笔者会在今后的研究方案中提出应对措施。

(感谢成远矿业开发股份有限公司及罗乃鑫、姚毅、唐玲彪、白宇在科研现场试验期间给予的支持！)