热力耦合作用下煤炭地下气化地表沉陷预测方法

李怀展，唐 超，郭广礼，陈 福，李 伟，周华安，黄建勇

（1.中国矿业大学 矿山生态修复教育部工程研究中心,江苏 徐州 221116；2.中国矿业大学 江苏省老工业基地资源利用与生态修复协同创新中心,江苏 徐州 221116）

0 引言

煤炭是我国能源安全的兜底保障，也是全球气化变暖的主要催化剂之一[1-3]。为应对全球气候变化及生态环境保护要求，各国积极进行能源转型，向着“ 清洁、低碳、安全、高效” 的现代能源体系发展[4-6]。煤炭地下气化作为一种高碳资源低碳化清洁开采技术，是煤炭行业绿色低碳高质量发展转型的重要选择[7-8]，尤其是双碳背景下，煤炭地下气化迎来了良好的发展机遇。煤炭地下气化是通过对煤炭资源进行有控制燃烧，将煤炭中的含能成分通过化学作用转换为甲烷、氢气等可燃气体，将煤基固废遗留在燃空区内，能够有效避免因采煤引起的安全和生态环境问题[9-10]。经过100 多年的发展，前苏联、中国、美国、澳大利亚、英国等许多国家的学者对煤炭地下气化工艺进行了研究[11-15]，并先后进行了40 余次的工业性试验，基本解决了煤炭地下气化工艺难题，为煤炭地下气化规模化应用奠定了坚实基础。

然而煤炭地下气化后仍会造成覆岩移动及地表沉陷，导致利用地下气化回收井工难以开采的“三下”压煤时，严重威胁地面建（构）筑物安全。如何准确地计算地表沉陷，已成为限制煤炭地下气化推广应用的瓶颈难题。目前部分学者围绕煤炭地下气化地表沉陷规律进行了研究，为煤炭地下气化地表沉陷预测模型构建奠定了坚实基础。如辛林等[16]以华亭原安口煤矿煤炭地下气化为对象，通过实测分析手段，研究了条带气化开采覆岩移动与地表沉陷规律；黄温钢等[17]利用数值模拟方法研究了条带气化开采主断面应力场、位移场以及覆岩移动规律；DERBIN等[18]基于数值模拟研究煤炭地下气化地表沉陷，明确了高温引起的岩体力学性质变化会影响地表沉陷；笔者[19]以乌兰察布煤炭地下气化场为原型，采用数值模拟方法研究了煤炭地下气化地表移动曲线形态及规律；EKNELIGODA 等[20]通过修改能量平衡方程来模拟燃煤过程，探究不同开采条件下地表沉陷变化。

同时笔者及其团队也围绕煤炭地下气化地表沉陷预测方法开展了一些工作，如2016 年，笔者采用数值模拟方法对比分析了煤炭地下气化开采和条带开采地表沉陷规律，并借鉴条带开采提出了相应的地表沉陷预测方法[21]。考虑到上述方法未考虑燃空区围岩高温效应、气化工艺特点，笔者[22]于2017 年分别提出了基于连续-随机介质理论和基于“真实采厚”的带状地下气化地表沉陷预计方法，进一步提高了煤炭地下气化地表沉陷预测精度，但提出的方法未兼顾焦化隔离煤柱的压缩变形以及简化了顶板下沉空间。总的来说，现有的煤炭地下气化地表沉陷预测方法未从煤炭地下气化引起地表沉陷的机理出发构建相应的预测方法，导致地表沉陷预测精度始终难以满足回收井工难采的“三下”压煤的技术需求。

基于此，研究分析煤炭地下气化引起地表沉陷的机理，在考虑煤炭地下气化高温热效应及其对煤岩体力学性质影响的基础上，研究建立热力耦合作用下地下气化覆岩挠曲计算方法以及焦化隔离煤柱压缩变形计算模型，进而建立热力耦合作用下煤炭地下气化地表沉陷精准预测模型。研究成果对于煤炭地下气化生产设计、难采“三下”资源回收等具有重要的理论和实践意义。

1 煤炭地下气化引起地表沉陷诱因及岩层挠曲力学解析

1.1 煤炭地下气化地表沉陷发生机理

后退式控制注气地下气化工艺的工作面布置方式类似常规条带开采，即气化1 条工作面、留设1 条煤柱，以保证气化过程的安全性及控制覆岩裂隙发育。但煤炭地下气化过程中燃空区围岩内会出现超过1 000 ℃高温，不仅会产生热应力，同时会造成围岩力学性质发生变化，形成燃空区围岩高温效应[23-24]；同时由于煤炭地下气化的工艺特点，导致气化工作面形态与常规井工开采的矩形工作面差异较大，且气化后会将粉煤灰等固体废弃物遗留在燃空区，造成煤炭地下气化与条带开采围岩变形差异较大。

总的来说，煤炭地下气化和常规条带开采的围岩裂隙特征差异较大，但均会产生一定高度的裂隙，且随着变形继续向上传递，在某一层会出现岩层实际拉应力小于其极限抗拉强度，此时其上部岩层会发生整体性与连续性的挠曲变形，引起覆岩协同沉降，进而传递到地表，造成地表移动与变形（图1a）。因此，气化工作面引起的覆岩挠曲变形会引起地表沉陷。

图1 煤炭地下气化覆岩移动与等效下沉空间Fig.1 Overburden rock movement and equivalent subsidence space of underground coal gasification

另外，由于煤炭地下气化过程中会留设隔离煤柱，当工作面气化后燃空区不能承受上覆载荷，此时覆岩的载荷将由隔离煤柱承载，进而会造成煤柱发生压缩变形。同时受气化扰动的影响，燃空区两侧煤柱会屈服软化，煤柱屈服区承载能力降低，加剧了煤柱弹性区的压缩变形。煤柱屈服区与弹性区压缩变形将诱发覆岩协同下沉，进一步增加了地表沉陷。

综上所述，煤炭地下气化引起地表沉陷的根源为覆岩挠曲与煤柱压缩变形，故可以通过计算覆岩挠曲值和隔离煤柱压缩量来预测地表沉陷，如图1b所示。

1.2 煤炭地下气化岩层挠曲力学解析

结合煤炭地下气化的工艺特点及岩层赋存特征，可以采用板模型计算覆岩的挠曲变形。根据板的厚度与板中面的最小尺寸，可以将板分为薄板与厚板，其中薄板和厚板需满足式（1）与式（2）[25]：

式中：b为 板的较短边；h为岩层厚度。

1.2.1 厚岩层板结构力学模型与解析

拉普拉斯变换法可以把对时间的偏导数从导热微分方程中消去，基于Laplace 变换，求解得燃空区顶板温度场函数的拉普拉斯变换式为

式中：s为持续温度时间；T0为初始温度；ai=λi/ρici，其中，λi为岩层的导热系数；ci为岩层的比热容；ρi为岩层密度；z为距热源的距离；Ai、Bi为待定系数。

如图2 所示，在煤炭地下气化过程中，由于岩层上下表面温差的影响，在岩层高度方向，由于温差引起的热应力为：

图2 煤炭地下气化岩板力学模型Fig.2 Mechanical model of underground coal gasification rock plate

式中：E为岩层的平均弹性模量；αi为岩层的线膨胀系数；ΔT为岩层在垂直方向上的温差。

为减小计算难度，从较为理想的角度分析岩层挠曲与破断，提出以下5 点假设：①认为岩层是连续、均质的各向同性的；②覆岩载荷以均布方式作用在岩层上；③不考虑温度场在水平方向上产生的热应力，且热应力以均布载荷方式作用在岩层上表面；④在计算岩层挠曲时，岩层弹性模量取气化后的岩层下表面值；⑤根据高温后岩层平均弹性模量计算岩层挠曲与破断。

根据符拉索夫（Vlazov）厚板理论，矩形厚板的平衡微分方程[26]为

式中：Di为第i层岩层抗弯刚度；hi第i层岩层厚度；ψx、ψy分别为x、y的转角；µi为第i层岩层泊松比；wi为第i层岩层挠度；Gi为第i层岩层剪切变形模量；q为均布载荷为岩层承载的覆岩应力。

式中：qi=γihi，γi为第i岩层的重力密度。

四周固支板结构的边界条件如下[27]：

将挠度和转角的位移函数展开成双三角级数形式，即：

在保证精度情况下，为简化计算取m=n=1，联立式（3）、（5）得

厚板的最大拉应力出现在下表面，厚板下表面最大拉应力为：

岩层破断形式主要为拉破坏，当厚板下表面拉应力极值 σmax超 过岩层的抗拉强度 σkq时，厚板将发生张拉破坏。

1.2.2 薄岩层板结构力学模型与解析

薄岩层弹性弯曲变形远小于它的厚度，符合弹性薄板的基本要求。根据薄板理论，第i层岩层在均布载荷q作用下的挠曲方程[28]为:

利用弹性薄板的应力计算公式，第i层岩层抗拉强度为

2 煤炭地下气化煤柱屈服宽度与压缩量计算方法

煤炭地下气化高温热效应不仅会使得煤柱的力学性质改变，并会产生侧向热应力。故计算煤柱屈服区宽度及压缩量时，必须考虑燃空区围岩高温效应对隔离煤柱的影响。随着气化工作面的不断推进，燃空区覆岩载荷将转移至两侧煤柱，使煤柱实际荷载增加，且煤柱屈服现象的发生使煤柱核区承载进一步增大，煤柱发生压缩变形，引起上覆岩层协同移动，进而诱发地表沉陷与变形。

2.1 基于D-P 准则的气化煤柱屈服宽度计算方法

Druck-Prager(D-P)破坏准则将材料视为理想的弹塑性体，可以考虑静水压力及中间主应力对岩石屈服特性的影响，在数值计算中得到广泛运用[29-31]。本节基于D-P 准则计算煤炭地下气化煤柱屈服区宽度。Druck-Prager 准则函数表达式为

其中，I1和J2分别为第1 应力张量不变量和第2 应力偏张量不变量。具体计算公式为

式中：λ、ζ为D-P 准则材料参数,在平面应变关联法则下,与Mohr-Coulomb (M-C) 准则中的参数围岩内摩擦角φ 和 黏聚力c之间存在如下关系

在生产实践中，3 个主应力(σ1，σ2，σ3)关系可由中间主应力系数η 表示为

η ∈[0,1]，与中间主应力的影响力成正比。

将式（18）代入式（16）得

假定：①煤体是连续、均质的各向同性弹性体；②煤柱与顶底板之间相同位置处的黏聚力与内摩擦角相同；③忽略煤柱自身重力影响；④煤柱顶部垂直应力达到煤柱支撑压力峰值极限强度时，煤柱发生塑性破坏；⑤将气化煤柱近似看作矩形。

如图3 所示，建立平面直角坐标系oxz，煤柱关于x轴对称，煤柱受水平热应力 σth，垂直热应力 σtv与覆岩垂直应力σos，其合力即为煤柱所受垂直应力σz。在高温影响下煤柱与顶底板接触面处的黏聚力和内摩擦角分别为c(T)、φ(T)，由于煤层界面处于极限平衡状态，结合煤壁温场函数，其剪应力 τzx满足

图3 塑性区微元力学模型Fig.3 Microelement mechanical model of plastic zone

屈服区内微元在x轴方向处于极限平衡状态，煤柱屈服区内宽度为dx的单元体在x轴的平衡方程为

式中：cc为待定常数。

煤炭地下气化高温热效应将改变煤柱力学性质，煤柱不同位置处微元的黏聚力与内摩擦角不同。但顾及煤柱自身力学性质变化将使公式难以计算。为提高计算结果可靠性并降低计算难度，取气化后平均值作为煤柱屈服区力学参数。在温度场的影响下，将在煤柱内部产生一个附加的水平热应力，其大小为：

在x=0处，σx=σth，即求解可得待定系数cc。

在煤柱塑性区与弹性区交界处，煤柱极限强度σzl与所受垂直应力相同，即：

根据上式求解可得煤柱屈服区宽度xp。

2.2 煤柱屈服宽度影响因素分析

煤柱温度场分布因素较多，限于篇幅，仅考虑热源温度与侧压力对煤柱屈服区宽度的影响。温度场改变煤柱力学性质的同时给予煤柱水平方向的侧向应力，但气化作业结束后煤柱所受侧向热应力急剧下降，而煤柱高温变质无法还原，因此需分别探究煤岩高温变质与热应力作用对煤柱屈服宽度的影响。

如图4 所示，煤柱屈服区宽度与煤柱高度、热源温度及煤柱抗压强度总体呈线性正相关，与侧压力呈负相关。结合图4b 和图4d，煤炭地下气化温度场对煤柱屈服区宽度有较大影响。其中，侧向应力对煤柱屈服区影响相对较大，因此必须考虑气化结束后无侧向热应力支撑条件下煤柱屈服区宽度的变化。

图4 不同因素对煤柱屈服区宽度的影响Fig.4 Influence of different factors on width of coal pillar yield zone

2.3 气化煤柱最大压缩量计算

煤炭地下气化煤柱载荷主要来源于覆岩应力与热应力，在两者共同作用下气化煤柱形成屈服区与塑性区。煤柱塑性区实际承载为

煤柱塑性区本构关系呈非线性关系，具有应变软化特性，煤柱塑性区内载荷与煤柱压缩量的关系可以表示[32]为

式中，wzy为煤柱压缩量；u0为煤柱最大压缩量；Ecp为煤柱塑性区平均弹性模量。

煤炭地下气化高温-热效应不影响上覆岩层容重，因此煤炭地下气化高温效应不直接影响煤柱实际载荷。煤柱实际载荷为：

式中：H为开采高度；cz为 燃空区宽度；γ为上覆岩层平均重力密度。

煤柱弹性区承载为

弹性区内，煤体载荷与煤柱压缩量wcp的关系为

3 煤炭地下气化场地地表沉陷预测模型

结合煤炭地下气化地表沉陷发生机理，煤炭地下气化的地表沉陷量为煤柱压缩和覆岩挠曲变形引起的地表下沉量之和。因此，可以分别计算煤柱压缩引起的地表下沉和覆岩挠曲变形造成的下沉沉降，进而基于叠加原理可以计算得到煤炭地下气化的地表沉陷[33-34]。

3.1 基于煤柱压缩的地表下沉预测

根据概率积分法预测模型，计算区域为2s0×2t0×wcp的煤柱弹性区压缩区引起的地表任意点的下沉计算公式为

区域为 2s1×2t1×wzy的煤柱塑性区压缩区引起的地表任意点的下沉计算公式为

3.2 基于岩层挠曲下沉的地表沉陷计算

岩层挠曲下沉的等效开采空间为煤柱压缩引起的岩层下沉wcp与 岩层本身的挠曲变形wf，岩层挠曲下沉的等效开采空间为

从影响函数的理论出发，三维条件下开采空间内某微元引起的地表沉陷为：

计算区域为 2u0×2v0×wos的岩层挠曲下沉引起的地表任意点的下沉计算公式为

煤炭地下气化等效开采空间引起地表任一点沉陷计算公式为

4 实例应用

4.1 研究区域概况

乌兰察布煤炭地下气化试验区位于内蒙古自治区中部，隶属于察哈尔右翼前旗玫瑰营子镇。气化试验区地表较为平坦，所在区域大多被第四系覆盖，少量第三系玄武岩覆盖，岩石裸露于地表。新奥气化采煤技术有限公司与中国矿业大学（北京）联合开展“煤炭地下气化产业化关键技术”研究与工业性实验，采用的气化工艺为目前国际先进的“条采-面采”气化炉后退式控制注气地下气化工艺，布置4 个气化工作面与3 条隔离煤柱。煤炭地下气化工作面平均埋深为275.75 m；煤层厚度为5 m，该层煤主要为褐煤。单条工作面气化时间为90 d，形成长170 m、宽16 m 的气化燃空区，设计气化工作面四条，工作面之间留设宽度为24 m 的煤柱，如图5 所示。

图5 “条采-面采”气化炉后退式控制注气地下气化工艺示意Fig.5 Schematic diagram of “strip mining-regional mining”gasifier backward controlled gas injection underground gasification technology

4.2 岩层破断判别与挠曲变形计算

根据地质钻孔资料，煤炭地下气化试验区的覆岩层结构及其岩性力学参数见表1[35]。

表1 岩层分布及其力学参数Table 1 Rock strata distribution and its mechanical parameters

温度场的拉普拉斯变换式计算复杂，且温度场在上覆岩层传递高度十分有限，在计算温度场传递时，可将上覆岩层简化为12 m 泥岩与15 m 砂岩。假定热源温度为1 000 ℃；原岩初始温度为25 ℃，岩层热物理参数见表2，岩层弹性模量随温度变化规律的拟合曲线函数表达式见表3[35-36]。

表2 岩层热物理参数Table 2 Thermal physical parameters of rock strata

表3 岩层弹性模量随温度的变化规律Table 3 Variation law of elastic modulus of rock strata with temperature

计算得5 m 厚泥岩未发生破断，但其不能承载上覆岩层所有载荷，岩层破断层数为5 层。15 m 厚粉砂岩可承载上覆岩层载荷，该岩层将控制上覆岩层移动与变形。经计算该层粉砂岩最大挠曲变形为14 cm，岩层三维挠曲如图6 所示。

图6 岩体三维挠曲示意Fig.6 Three-dimensional diagram of rock flexure

4.3 气化煤柱最大压缩量计算

煤柱热物理参数：热源温度1 000 ℃，原岩初始温度25 ℃，导热系数0.5 W/(m·K)，比热容1 670 J/(kg·K)，时间t为5 184 000 s，密度1 400 kg/m-3，热扩散系数9.76×10-8m-2/s；常温状态下煤柱力学参数：煤层与顶板、底板间的黏聚力700 kPa，煤层与顶板、底板间的内摩擦角28°，煤柱黏聚力1 300 kPa，煤柱内摩擦角30°，煤柱极限强度16 MPa；煤柱力学性质随温度变化见式（42）-式（44）[37-38]：

式中，c0为 常温状态下煤柱与顶底板的黏聚力；φ0为常温状态下煤柱与顶底板的内摩擦角；Ten0为常温状态下煤柱的极限强度；E0为常温状态下煤柱的弹性模量。

计算得煤柱屈服区宽度为4.4 m；煤柱屈服区压缩高度为5 cm；煤柱弹性区压缩高度为8 cm。

4.4 试验区地表沉陷预测

参考《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》中依据覆岩岩性综合评价系数计算公式，计算得下沉系数为0.6。同时，根据4.1 节与4.2 节对煤炭地下气化等效开采空间计算，结合式（35）计算得乌兰察布煤炭地下气化试验区地表最大下沉为45 mm，地表下沉等值线如图7 所示。试验区CG05 监测点地表下沉为36 mm，距离工作面中部为10 m。由此可知，新提出的方法符合实际情况，相对于目前基于“真实采厚”地表沉陷预测方法得到的地表最大下沉为95 mm 而言，预测结果精度更高，可以更好地满足回收井工难采的“三下”压煤的技术需求。

图7 煤炭地下气化地表下沉等值线Fig.7 Surface subsidence contour of underground coal gasification

5 结论

1）煤炭地下气化引起地表沉陷的根源为覆岩挠曲与煤柱压缩变形，可以通过计算覆岩挠曲值和隔离煤柱压缩量来预测地下气化地表沉陷。同时基于板理论，建立了热力耦合条件下煤炭地下气化覆岩挠曲计算方法及岩层破断判别方法。

2）基于D-P 破坏准则，构建了热力耦合作用下煤炭地下气化隔离煤柱屈服区宽度计算模型，并根据煤柱屈服区与弹性区承载特征，建立了煤炭地下气化隔离煤柱最大压缩计算模型。另外，煤柱高度、温度场与煤柱极限强度均对煤柱屈服宽度有较大影响，在计算气化煤柱屈服区宽度时需顾及热应力时序变化的影响。

3）基于叠加原理，构建了热力耦合作用下煤炭地下气化地表沉陷预测新方法。该方法在乌兰察布煤炭地下气化试验场进行了应用，得出试验场地表最大下沉为45 mm，相对于现有方法而言，论文提出的方法与实测结果更切合，证明了方法的准确性。