高承压水地层冻结温度场数值模拟研究＊

张 瀚 李栋伟 于 奇

（1.矿山地下工程教育部工程研究中心，安徽省淮南市，232001；2.安徽理工大学土木建筑学院，安徽省淮南市，232001）

高承压水地层冻结温度场数值模拟研究＊

张 瀚1，2李栋伟1，2于 奇2

（1.矿山地下工程教育部工程研究中心，安徽省淮南市，232001；2.安徽理工大学土木建筑学院，安徽省淮南市，232001）

通过高承压水裂隙岩体室内土工试验原始热物理参数、冻结钻孔的实际偏斜资料和测温孔的测温数据，利用ANSYS软件对裂隙岩体的物理力学参数进行反演，得到冻结岩体导热系数、比热容的最优参数。并以最优参数为各岩层的物理参数，利用ANSYS软件模拟各岩层的冻结壁厚度。

冻结法凿井 高承压水 冻结温度场 裂隙岩体 数值模拟 ANSYS

高承压水地层在深井多圈冻结管条件下冻结壁场形成规律十分复杂，裂隙岩体深井冻结壁的冻结温度场在形成过程中受到渗流场和应力场的相互影响。渗流场的变化将引起地层含水率的变化，而地层含水率的变化又将引起裂隙岩体导热系数、比热容等热物理力学参数的变化，进而影响整个冻结温度场的传热及分布。同时应力场的变化又将引起裂隙岩体密度、孔隙率等变化，从而影响到整个温度场的分布。本文以平煤一矿回风井3个不同的冻结岩层段为原型，结合现场实测的测温孔温度，通过数值计算的方法模拟冻结温度场的发展规律，为高承压水层的冻结工程的设计和施工提供一定的参考。

1 工程概况

平顶山平煤一矿回风井位于郏县李口乡李口村南，属于山前洪坡积平原地带。井筒设计采用立井开拓方式，冻结法施工，井筒全深1075 m。平煤一矿回风井井筒净直径6.5 m，井壁最大厚度0.8 m，井筒最大掘进直径8.1 m，采用主排孔＋防片帮孔冻结方案进行施工。主排孔采用全深冻结方式，其深度为660 m。主排冻结管布置圈径14 m，管数28根；防片帮孔冻结深度为120 m，防片帮冻结管布置圈径11.1 m，管数14根。深度300 m以内主排孔采用ø159 mm×6 mm低碳钢无缝钢管内管箍对焊联接，深度超过300 m选用ø159 mm×7 mm低碳钢无缝钢管内管箍对焊联接。防片帮孔选用ø159 mm×5 mm低碳钢无缝钢管内管箍对焊。设计布置3个测温孔，分别布置在冻结壁中间薄弱位置、冻结壁内侧界面位置、冻结壁外侧最大孔间距位置且冻结面外侧界面位置，并且为准确测量裂隙岩体冻结温度场温度，测温元件放置在测温管外侧，将不同的测温线用塑料管绑扎使得下放过程中对电缆线无损坏。冻结管及测温管施工布置如图1所示。冻结壁厚度表层由拉麦公式计算得2.3 m。基岩段冻结壁厚度按维亚洛夫-扎列茨基公式计算得3.9 m，故最终确定平煤一矿回风井的冻结壁厚度为3.9 m，并且要求深度超过300 m后井帮温度不超过-4℃。

图1 冻结孔与测温孔布置图

2 冻结壁温度场热物理参数的数值反演

在进行冻结壁温度场热物理参数反演的过程中，由于影响冻结壁温度场的因素较多，所以在这里只考虑土体的密度 （ρ）、比热 （C）、导热系数（λ）。并以92 m卵石层热物理参数的反演为例进行说明。

根据平煤一矿冻土试验报告选取92 m卵石层处，土体的密度包括2400 kg／m3、2500 kg／m3、2600 kg／m3；土体的比热由于受温度的影响这里取等效比热，92 m卵石层冻结状态时等效比热分别取0.678 kJ／ （kg·K），0.786 kJ／ （kg·K）、0.834 kJ／（kg·K），而在未冻结状态下，等效比热几乎相同，故在未冻结状态等效比热取0.925 kJ／（kg·K）；导热系数同比热一样受温度的影响，且在冻结状态时导热系数几乎均为2.582 W／（m·K），而在未冻结状态下导热系数分别为2.146 W／（m·K），2.237 W／（m·K），2.358 W／（m·K）。故最后运用ANSYS进行数值反演时取密度、冻结状态下的等效比热和未冻结状态下的导热系数共9个变量进行反演。

因为土体的密度、冻结状态下的比热、未冻结下的导热系数均有3个影响因素，若每个排列组合均进行ANSYS反演，则有3×3×3次试验，工程量太大，这里采用正交试验，按L9（3＾3）排列进行反演试验，具体安排如表1。

表1 92 m处卵石层试验方案

利用ANSYS软件对各试验进行模拟，再将模拟温度与实测温度进行对比，最后利用最小二乘拟合得出最优解。其中最小二乘拟合是数学上的一种近似和优化，利用数值模拟得出的曲线，使之在坐标系上与已知数据之间的距离的平方和最小，是离散情形下的最佳平方逼近。其式为：

式中：Δδ2——最小二乘拟合误差；

S（xi）——测温孔模拟值 （i＝5，10，15，…150）；

yi——实测测温孔温度值。

如果取每天的实测温度与模拟温度进行比较，计算起来很繁琐，为了简化计算过程取有效冻结期0～150 d每隔5 d的冻结温度，并从第5 d开始最小二乘拟合计算。在这里取最优解第五组试验进行说明，由式 （1）计算得：Δδ2＝236.023。第五组试验实测与模拟值对比见图2。

从图2可以看出，温度在冻结初期下降较快，随着冻结天数的增加，到一定的时间后曲线趋于平缓。1＃测温孔位于外圈孔之外，2＃测温孔位于内圈孔以内，3＃测温孔位于内圈孔和外圈孔之间。各测温孔的降温速度为：3＃测温孔＞2＃测温孔＞1＃测温孔。说明两圈冻结管之间的冻结速度最快，内侧冻结管内侧部分处冻结速度稍慢，而在外圈冻结管外侧部分处的冻结速度最慢。同理也可获得207 m处砂质泥岩层和373 m处中粒砂岩处各热物理参数，如表2所示。

图2 第五组试验实测与模拟值对比

表2 热物理计算参数取值

3 冻结壁温度场的超前预测及分析

以实际掌握的冻结孔偏斜资料并结合表2中所得到的各个层位的最优热物理参数，利用ANSYS软件对各个层位进行模拟预测。以92 m卵石层为例，得到其冻结壁温度场的分布规律，冻结100 d的温度场冻结云图见图3，冻结温度场100 d时，-2℃等温线分布见图4。

图3 冻结100 d卵石层温度场云图分布

从图3和图4中可以看出在92 m卵石层处冻结100 d时冻结壁已基本交圈。同时可利用ANSYS软件进一步模拟得出3个层位冻结壁厚度和冻结壁平均温度随时间的变化关系，见图5和图6。

图4 冻结100 d卵石层-2℃冻结温度等温线分布

从图5中可以看出，冻结壁厚度初期增长较快，之后随着冻结时间的延长，当冻结时间达到150 d后增速就趋于缓慢，而且对于不同的3个岩层变化规律基本相同，但各自的冻结速度不同，其中373 m的中粒砂岩层的冻结壁厚度增长速度最快，最快发展速度达到0.056 m／d，相对于一般岩层的冻结壁厚度发展速度，高承压水层的发展速度明显偏快。当冻结时间达到250 d时，卵石层、砂质泥岩层和中粒砂岩层的冻结壁厚度分别为：4.4 m，5 m和4.7 m。从图6可以看出冻结壁平均温度变化情况大致与冻结壁厚度变化相同，初期增长速度快，当达到一定的冻结温度后，冻结壁的平均温度趋于稳定，此后再增加冻结时间对于降低冻结壁的平均温度作用不甚明显。

最后利用ANSYS软件对各个层位的井帮温度进行预测，以373 m中粒砂岩层为例，沿井帮周长分为120段等分，井帮温度变化曲线图如图7所示。图7显示井帮温度均低于-4℃，满足工程需要。

图7 373 m中粒砂岩层处井帮温度沿周长方向变化曲线

4 结论

（1）以实际冻结孔的偏斜资料建立的有限元数值模型，进而进行各个层位的热物理力学参数的反演拟合并，最终预测各层位的冻结壁厚度，不仅可行，对于工程的设计和施工也是极其必要的。

（2）通过对冻结壁厚度的数值模拟预测可知，冻结壁厚度初期发展快，当冻结时间达到一定程度时，增速趋于缓慢，但最大冻结壁发展速度为0.056 m／d，相对于一般岩层的冻结壁发展速度明显偏大。

（3）通过本文的数值模拟方法，可以及时判断冻结壁的薄弱环节，进而强化薄弱部位的冻结，为深井裂隙岩体冻结施工提供一定的指导作用。

［1］陈军浩，汪仁和.多圈管冻结下冻结壁温度场融化特性的实测研究 ［J］.中国煤炭，2009（12）

［2］吴王青，李栋伟，马乾坤.白垩系地层冻结壁温度和压力实测分析 ［J］.中国煤炭，2012（3）

［3］刘建国.高承压水地层冻结温度场形成特性现场实测研究 ［J］.工程技术，2011 （18）

［4］王衍森，杨维好，任彦龙.冻结法凿井冻结温度场的数值反演与模拟 ［J］.中国矿业大学学报，2005（5）

［5］张朝晖.ANSYS热分析教程与实例解析 ［M］.北京：中国铁道出版社，2007

［6］沈科伟，荣传新，余文文.查干卓尔矿副井冻结温度场有限元分析 ［J］.广西大学学报 （自然科学版），2011 （4）

［7］李栋伟，旺仁和，胡璞.多圈管冻结瞬态温度场有限元数值分析 ［J］.煤田地质月勘探，2007（2）

［8］刘建国.深井裂隙岩体冻结温度场数值模拟研究［J］.工程技术，2011 （19）

［9］马芹永.人工冻结法的理论与施工技术 ［M］.北京：人民交通出版社，2007

Numerical simulation of freezing temperature field in strata with high confined water

Zhang Han1，2，Li Dongwei1，2，Yu Qi2
（1.Research Center of Mine Underground Engineering of the Ministry of Education，Huainan，Anhui 232001，China；2.School of Civil Engineering and Architecture，Anhui University of Science and Technology，Huainan，Anhui 232001，China）

Based on the original thermophysical parameters of fractured rock with high confined water obtained from lab soil tests，the practical deflection data of the freezing borehole and the temperature data from thermometry hole，the paper uses the ANSYS software to conduct an inversion analysis on the physics-mechanical parameters of the fractured rock，and obtains the optimal parameters of specific heat capacity and thermal conductivity coefficient of the freezing rock.And taking the optimal parameters as the physical parameters of every stratum，the paper simulates the frozen wall thicknesses of the strata with ANSYS software.

freezing sinking，high confined water，freezing temperature field，fractured rock，numerical simulation，ANSYS

TD353

A

国家自然科学基金 （41271071）；教育部新世纪优秀人才支持计划 （NCET-11-0887）

张瀚 （1985-），男，硕士研究生，主要从事冻土力学与工程研究工作。

（责任编辑 张毅玲）