局部差异冻结作用于既有井壁的冻胀力研究

(安徽理工大学 土木建筑学院， 安徽 淮南 232001)

0 引言

近几十年来，我国矿井建设全面开展，其中两淮矿区为我国大型煤炭资源储藏基地，是华东地区重要的能源基地，多为深厚冲积层所覆盖，煤矿均采用立井开拓[1]。目前，国内已有多例服役期矿井出现井筒破损淹水等事故。针对普通法施工矿井发生井筒淹水事故后多采用冻结法完成事故井修缮；冻结法处理事故井过程中，为确保事故井破损段井壁的安全修复，同时不使非事故段井壁由于土体冻胀产生二次破坏，既要保证冻结壁有一定的强度和厚度，又要严格控制冻结壁的发展，从而在实际工程中面临着许多技术难题。因此局部差异冻结条件下冻结壁的形成规律和研究方法至关重要。

国内外学者对于冻胀力与冻胀机理的研究颇多。EVERETT[2]根据毛细冰理论对水分迁移冻胀进行了解释，称之为第一冻胀理论。MILLER[3]为解释水分迁移过程中不连续冰透镜体的形成过程，提出了第二冻胀理论，并基于该理论提出了刚性冰模型。O’NEILL等[4]对该模型进行了完善，给出了由定量冰、土和水之间的分量总应力得出冻结锋面冰晶体产生的准则，并建立了较为完善的计算体系。国内荣传新等[5]建立了考虑冻结壁、井壁及周围土体共同作用的粘弹性计算模型，给出了外层井壁上的冻结压力解析式，以及冻结壁和井壁的应力分布规律。马茂艳等[6]基于符合深井冻土蠕变特性的改进西原体模型，利用ABAQUS软件计算分析了深部冻结井的掘砌过程，获得了作用于外层井壁的冻胀力发展变化规律，结果表明土体埋深、冻结壁温度、土体冻胀率等因素均影响冻胀力的大小。业界学者均对人工冻结温度场与应力场做出了有力探索。

在局部差异冻结的研究方面，肖朝昀等[7]针对上海地铁一区间采用四排孔局部冻结修复工程，分析了积极冻结期排内冻土壁交圈时间、发展速度，并利用公式计算出整个冻结期排外冻土壁单侧发展厚度和发展速度，得出多排局部冻结冻土壁温度场特征。根据复杂地层的冻结难易程度，提出了一种冻结管差异温度优化冻结方案，进行了差异冻结试验与数值分析[8-9]。石荣剑等[10]为解决地下有流水时盐水帷幕冻结方法无法形成封闭冻结壁的问题，采用低温液氮对局部未封闭区域进行补强冻结方式，可有效封堵地下流水，解决盐水冻结不能形成封闭冻结帷幕的难题。荣传新等[11]针对板集矿区副井井筒修复过程中的局部差异冻结条件技术进行了较为系统的分析，探究了整个井壁修复过程中冻结温度场的发展规律，为类似工程提供了参考。

综上所述，目前国内外对于局部差异冻结技术的研究成果较少，鲜有冻结修复期间既有井壁的冻胀力与安全性研究，并且关于该技术在破损井筒修复中的运用尚无成功先例，为了给类似工程提供一定的参考依据，笔者采用ABAQUS/Standard数值计算软件，对局部差异冻结条件下井壁修复过程中的冻胀力进行了模拟与研究，为指导冻结法凿井设计以及保障冻结施工安全提供帮助。

1 工程

1.1 工程背景

板集煤矿为新建的大型矿井，副井井筒主要特征参数见表1。

表1 副井井筒技术特征Tab.1 Technical characteristics of auxiliary shaft

预制井壁支护深度640 m，井筒采用变断面形式。预制钢筋混凝土井壁支护深度375.0 m，混凝土强度等级为C30～C60；预制双层钢板混凝土井壁支护深度为375.0～640.0 m，钢板厚度10～30 mm，混凝土强度等级为C60～C70。地层由老到新依次有寒武系、奥陶系、石炭系、二叠系和新生界。新生界松散层自上而下大致分为一含、一隔、二含、二隔、三含、三隔和四含共计4个含水层(组)和3个隔水层(组)。其中，四含厚度介于39.80～175.70 m，平均83.66 m左右，由上部的厚层中、细砂层和下部的砂砾层组成，其中砂层厚度大，富水性中等，全井田均有分布，且大部分与基岩直接接触，是基岩含水层的重要补给水源[12-14]。

2009年4月18日，板集煤矿副井井筒发生突水涌砂事故，专家组对事故调查后认为：突水事故的直接原因是上覆岩层受重复扰动影响，使得井筒受拉破裂，诱发井筒壁后四含水涌入井筒。通过研究分析，板集煤矿副井井筒确定采用“抛—注—冻—修—防”综合治理方案，即先向井筒内部抛填石子恢复井壁的三向受力状态，然后通过地面预注浆充填加固地层，再采用冻结法形成井筒排水修复保护帷幕，最后进行井筒套壁修复加固。根据专家对出水点位置分析的结果，出水点应位于表土段第三隔水层以下至基岩风化带，其具体位置不明确，且井壁破坏程度及其范围亦不明确。为减小冻结施工对上部完好的既有井壁的影响，并确保下部破损井壁的安全修复，工程采用局部差异冻结技术进行地层冻结[11]。

1.2 局部差异冻结井壁修复方案

板集煤矿立井修复工程井壁突水修复段高1.5 m，冻结孔布置参数见表2。

表2 冻结孔布置参数Tab.2 Freezing hole arrangement parameters

设计盐水温度为-30～-34 ℃，冻结壁厚度为5.0 m，冻结壁平均温度为-15 ℃。井壁全深采用双排管冻结方式，实际钻孔时存在一定偏斜，其中-300 m与-450 m层位冻结管布置如图1所示。为减小上部完好的强度较低的钢筋混凝土井壁遭受冻胀力影响，并确保下部破损井壁安全修复，结合对地质资料的分析结果，最终确定采取局部差异冻结，冻结深度为673 m，累深380 m以上地层采用双供液管工艺实现控制冻结。外排冻结管采用Ф159×5～7的无缝钢管内管箍连接，内下Ф75×6的聚乙烯塑料软管作供液管。内排冻结管内排冻结管立面结构如图2所示，其-380 m以浅采用Ф168×6低碳钢无缝钢管，外管箍焊接联接，内下双Ф70×5聚乙烯塑料软管；-380 m以深采用Ф159×7低碳钢无缝钢管，内管箍焊接，内下Ф70×5聚乙烯塑料软管[11]；在约-380 m位置处设置挡水隔板，以保证低温盐水仅自-380 m以深与冻结管的外管壁接触，控制有效冻结段仅覆盖突水区域所在四含段与邻近的三隔段与基岩段，以实现局部差异冻结控制效果。监测的盐水温度变化曲线如图3所示，共布置7个测温孔以实时监测冻结壁发展，内圈设置4个，内外圈之间设置2个，外圈冻结孔外设置1个。另在冻结壁内侧共布置4个水文孔，以释放早期冻胀力。

(a) -300 m层位

(b) -450 m层位

图2 内排冻结管结构Fig.2 Structure diagram of inner row freezing pipes

图3 盐水温度随时间变化曲线Fig.3 Temperature variation of brine with time

2 数值计算方法

如不考虑冻结管偏斜的影响，并假定冻结管内的盐水温度沿管壁轴向分布均匀，则立井轴向(z方向)上的土层温度梯度近似为0，计算模型可简化为平面应变问题[15]。由于土的冻胀效应影响因素众多，在进行冻胀力的数值分析时，作如下基本假设[16]：①不考虑水分迁移过程；②按平面应变问题进行计算；③对土体颗粒和冰晶，不考虑其压融效应；④以钙质黏土层和粉质黏土层为模拟对象，单一土层为均质的各向同性弹性体，但其弹性模量的取值随温度变化而变化。

2.1 温度与位移耦合方程

① 温度场控制方程

冻结温度场为带相变的瞬态热传导问题，基于Fourier热传导方程和能量守恒定律，将微元体看作刚体，其热量平衡控制微分方程[17]如下所示：

(1)

式中,C为土的容积热容,J/(m3·K)，C=Cdρ，Cd为土的比热容,kJ/(kg·K)；λ为土的导热系数,W/(m·K)；L为水的相变潜热,kJ/kg，一般取334.56 kJ/kg；ρi为冰的密度,kg/m3；θi为土体的体积含冰量,%。

② 本构方程

在温度变化的过程中，弹性体的应变由两部分组成:①自由膨胀引起的正应变分量ε0=αΔT，其中α为膨胀系数，ΔT为温度变化值,自由膨胀引起的剪应变分量为0。②热膨胀时，土体的应力-应变关系符合广义胡克定律：

(2)

联立温度场微分方程和热力学本构关系，可以得到热力耦合控制方程。并按照弹性力学方法，建立热应力方程。其中平衡方程只与物体受力有关，而与产生力的原因无关；几何方程中应变只与位移有关，而与引起位移的原因无关，所以平衡方程和几何方程的形式不变。

2.2 冻胀力计算参数的选取及计算方法

土体假设为线弹性材料，但其弹性模量和泊松比的取值随温度变化而变化；土体的冻胀率通过原状土的冻胀试验确定；根据室内试验得出导热系数、比热以及其他指标随温度场的变化规律，在进行数值计算时，通过建立比热、导热系数随温度场变化的场变量赋予冻土的热物理参数。模型考虑热力的耦合，反之不可。进而建立了随温度场变化的热力耦合的数值计算方法。

笔者结合板集矿副井修复人工冻结工程，利用上述方法代入随温度场变化的参数指标进行计算，即首先建立温度场模型、加载，求解温度场模型并进行瞬态分析，后处理得到节点上的温度。然后采用和温度场相同的数学模型进行结构分析。结构分析中首先通过地应力平衡功能构建井壁与周围土体的初始稳定场(即消除位移场保留初始应力解)，后将求解温度场得到的节点温度作为温度荷载加载到应力场以进行耦合分析，进而得到冻胀力的分布规律。

3 数值计算

3.1 模型参数

副井-300 m与-450 m层位既有井壁为预制双层钢板混凝土支护，钢板厚度为10～30 mm，混凝土强度等级为C60，弹性模量为36 GPa，泊松比为0.167，忽略冻结过程中既有井壁的热膨胀。数值计算中各岩土层的物理力学参数取值见表3。

表3 各层位土体的物理力学参数Tab.3 Physical and mechanical parameters of soil layers at various horizons

3.2 几何模型与边界条件

① 几何模型、网格

根据井壁修复工程冻结孔造孔情况，采用ABAQUS/Standard数值分析软件，选用冻结孔终偏位置进行建模计算。数值分析模型如图4所示，考虑到冻结施工的影响范围，取冻胀影响圈径为80 m；同时为提高计算精度，对冻结管周围及井帮处进行网格加密，同时对于规则的井壁选用更为合理的结构化网格。共划分190 392个四边形单元，以及198 942个结点。

(a) 工程计算模型

(b) 井壁网格划分

② 荷载及初边值条件

温度场设置：冻结管温度按工程实测中的冻结器温度选取；初始温度场按井壁修复工程开始冻结时的实际温度场设定(约28 ℃)。同时出于简化模型计算的考虑，设置既有井壁内表面与水接触处的散热系数为0.54 W/(m2· ℃)[18]，根据工程实际抽水测温得知，水源温度为10 ℃左右。

应力场设置：根据重液公式计算永久水平地压P=0.013H，则-300 m层位的初始水平应力为3.9 MPa，-450 m层位的初始水平应力为5.85 MPa；取侧压力系数γ=μ/(1-μ)=0.538 5，故-300 m层位初始竖向应力为7.24 MPa，-450 m层位初始竖向应力为10.86 MPa。土体外边界约束法向位移。

3.3 计算结果及分析

① 瞬态温度场发展规律

为了验证数值计算的合理性，提取模型测温孔位置处的节点温度，并选取-450 m层位处测温孔的实测数据与数值计算结果进行对比(图5)，其中C1测温孔位于外圈管外侧，C2测温孔位于两圈管之间，C6测温孔位于内圈管内侧。

如图5所示，3个测温孔的温降曲线数值计算与现场实测结果基本一致，其误差绝对值均小于5%，而在冻结初期受实际工程中施工热扰动等的影响，温度偏差较大，但均不超过2.5 ℃，这说明数值计算得到的瞬态冻结温度场可较为真实地反映工程实际情况，这也将作为应力场计算的基础。

(a) C1测点实测与模拟对比

(b) C2测点实测与模拟对比

(c) C6测点实测与模拟对比

井壁冻结修复期间-300 m和-450 m层位冻结壁发展状况如图6所示。①-450 m层位于94 d交圈，交圈时井帮(井壁与土体交界面处)平均温度为-1.28 ℃；-300 m层位滞后22 d于116 d交圈，交圈时井帮平均温度为3.15 ℃。冻结60、172 d时，-450 m层位井帮平均温度为4.43 ℃、-17.75 ℃，同时期-300 m层位井帮平均温度为11.62 ℃、-2.74 ℃。②冻结94 d时，-450 m层位冻结壁平均温度为-13.97 ℃，冻结116 d时，-300 m层位冻结壁平均温度为-11.64 ℃；冻结172 d时，-450 m层位冻结壁平均温度为-21.53 ℃，同时期-300 m层位冻结壁平均温度为-16.33 ℃，较-450 m层位高5.2 ℃。③冻结94 d时，-450 m层位冻结壁平均厚度为5.11 m；冻结116 d时，-300 m层位冻结壁平均厚度4.13 m。冻结172 d时，-450 m层位冻结壁平均厚度为6.02 m，同时期-300 m层位冻结壁平均厚度为5.07 m，较-450 m层位薄0.95 m。同时期，-300 m层位的井帮温度、冻结壁平均厚度及平均温度均明显低于-450 m层位，这表明差异冻结技术可以显著减小冻结壁向内发展的速度，达到了预期的冻结控制目的。

图6 冻结温度场云图Fig.6 Freezing temperature field nephogram

② 应力场发展规律

将节点应力计算结果减去初始地应力即为冻胀力，后文所述冻胀力均为既有井壁与土体交界面处节点的径向应力。分别提取-300 m钙质黏土层以及-450 m粉质黏土层距离该立井井筒中心r=4.65 m(井帮)界面处节点的冻胀力，绘制不同方位上的冻胀力分布状况如图7所示，及平均冻胀力时间历程曲线如图8所示；同时将两层位在东、北偏东30°、北偏西6°、北偏西45°、南偏西45°、南偏东15°六个特征方位上的冻胀力数据整理见表5。

图7及表5结果表明：①粉质、钙质黏土层冻胀力呈现出类同心圆的发展趋势，且由于冻结孔在实际成孔的过程中存在偏斜，两层冻胀力存在显著的不均匀性，最大冻胀力主要出现在东北侧、西北侧、南侧、西南侧。②在整个积极冻结期内，-300 m层位井帮处的冻胀力较小，最大不超过0.70 MPa，占初始地应力的17.94%，而-450 m层位的冻胀力较大，最大值为2.43 MPa，占初始地应力的41.54%，可见局部差异冻结技术可以有效避免冻土冻胀力对既有井壁的不利影响。③冻结94 d时，差异冻结层位与全深冻结层位在北偏东30°、北偏西6°、北偏西45°、南偏西45°、南偏东15°方位上的井帮冻胀力差值分别为 0.167、 0.138、 0.133、 0.215、 0.322 MPa。冻结172 d时，30°、北偏西6°、北偏西45°、南偏西45°、南偏东15°方位上的井帮冻胀力差值分别为1.757、1.771、1.672、1.803、1.779 MPa。可见在整个积极冻结期内，-300 m层位的冻胀力始终小于-450 m层位，冻结控制效果显著，采取局部差异冻结技术能有效减弱土体冻胀对既有井壁的影响，避免完好的井壁由于冻结施工产生二次破坏。

差异冻结层位与全深冻结层位在北偏东

(a) -300 m层位

(b) -450 m层位

图8 既有井壁与土体交界面处冻胀力的时间历程曲线Fig.8 Time-history curve of frost heaving force at the interface between the existing shaft wall and soil

表5 特征方向上的冻胀力Tab.5 Frost heaving force in characteristic positions

由图8可以看出,① 0～30 d时，-300 m层位与-450 m层位的冻胀力几乎为0，这主要是因为冻结管内的盐水温度还未降至指定温度，冻结壁正在缓慢形成。② 30 d至各层位冻结壁交圈前(-300 m层位为30～116 d, -450 m层位为30～94 d)，两层位的冻胀力呈现出线性增长趋势。在各层位交圈后的14 d内有一段近似阶跃性的增长，其后冻胀力增长速度放缓，最终基本趋于稳定。③冻结172 d时，-450 m层位冻胀力达2.35 MPa，-300 m层位冻胀力为0.58 MPa，约为-450 m层位的1/4。

3.4 井壁安全性评价

冻结法施工会直接影响混凝土的强度发展、变形与裂缝，因此既有井壁的安全性综合冻胀力和混凝土应变测值进行评价。其中混凝土应变是井壁安全性评价的最直接指标。

对数值计算得到的结果进行后处理，井壁等效应力结果如图9所示。-300 m层位和-450 m层位的最值结果统计见表6。将既有井壁混凝土等效应力与其混凝土强度进行比较，对既有井壁的安全性进行评价。

图9 井壁等效应力结果Fig.9 Von Mises stress of the shaft wall

表6 两层位井壁安全性评价指标结果Tab.6 Results of safety evaluation indicators at two layers

由图9及表6可知：①井壁等效应力随冻结时间的发展规律与冻胀力发展规律完全吻合。在冻结初期，井壁的应力未受到冻结施工的影响，两层位井壁应力均由初始地应力引起；随着冻结施工的进行，-300 m层位井壁混凝土等效应力增长缓慢，而-450 m层位增长较快。②对比两个层位井壁混凝土的等效应力可以发现，至冻结172 d，-300 m层位的最大等效应力为9.31 MPa，约为同时期-450 m层位的49.97%。且-300 m层位的井壁应力增长明显小于-450 m层位，约为3.28 MPa。然两层位混凝土等效应力均未超过C60混凝土的抗压强度标准值38.5 MPa，既有井壁是安全的。这表明本工程局部差异冻结达到了预期目的，局部差异冻结可以有效的减弱土体冻胀对完好井壁的影响。③在板集副井内层钢板高强钢筋混凝土井壁结构中，混凝土环向应变较大。-300 m层位混凝土的环向应变为-125.54 με，-450 m层位为-355.83 με，均远小于C60的极限应变值，说明在冻结修复过程中，既有井壁始终处于弹性受力状态，符合设计要求，冻结修复期间井壁结构并未开裂，安全可靠。

4 结论

① 板集煤矿副井井壁破损和突水修复工程采用了沿着深度方向进行局部差异冻结的修复方案。根据现场测温孔的实测数据及土体的物理力学参数，采用冻结孔的终偏位置进行数值计算，对冻结过程中的温度场和应力场进行分析。温度场计算结果表明，-300 m层位交圈时间较-450 m层位晚22 d；冻结172 d时，两层位冻结壁平均厚度均达到了设计标准，且-300 m层位冻结壁平均温度较-450 m层位高5.2 ℃，井帮温度较-450 m层位高15.01 ℃，这表明差异冻结技术可以显著减小冻结壁向内发展速度。

② 在整个冻结期内，冻胀力随着冻结施工的进行而增大，土体的冻胀对既有井壁会产生一定的影响。两层位冻胀力呈现出类同心圆发展形式，且冻胀力的发展有明显的不均匀性。其中，-300 m层位的冻胀力增长缓慢，且增长较小，172 d时冻胀力占初始地应力的17.94%。-450 m层位井帮处的冻胀力增长迅猛，172 d时冻胀力占初始地应力的41.54%。井帮处冻胀力自冻结30 d时呈现出线性增长趋势，在交圈后的14 d内有阶跃性增长，其后冻胀力增长速度放缓，最终基本趋于稳定。

③ 数值计算结果表明，在板集煤矿副井井壁结构冻结修复过程中，既有井壁混凝土的环向应变较大。-300 m层位混凝土的环向应变为-125.54 με，-450 m层位为-355.83 με，远小于C60混凝土极限应变值；井壁混凝土最大等效应力为18.63 MPa，低于C60混凝土的抗压强度标准值，冻结修复过程中井壁结构安全可靠。