0 引 言
煤炭地下气化技术(Underground Coal Gasification,UCG)是一种采用化学开采方法将煤炭在地下原位直接进行燃烧气化产生可燃气体的高碳资源低碳化开发清洁能源新技术，也是资源与环境协调的煤炭绿色开采技术[1-2]。煤炭地下气化将建井、采煤、气化三大工艺合而为一，具有产气成本低、安全性高和环境效益好的优点[3-5]，前苏联、欧洲、美国、日本、澳大利亚、中国及亚洲多个国家都在此方面进行了大量的理论研究、试验研究和现场工业性试验[6-8]。在现阶段优先适用于不可采煤层、低品位煤层及深部煤层的原位开采与转化，可提高煤炭利用价值，带动煤炭、电力与化工等传统产业发展，符合国家节能环保及发展清洁能源政策，具有广阔的应用前景[2,9-11]。
与传统煤炭开采方式不同，煤炭地下气化采场附近顶底板岩层要经受高温热作用(最高温度可达1 000 ℃[12-13])，因此顶底板岩层移动、变形和破坏不但受采动覆岩应力影响，还受到动态变化的温度产生的热应力及岩石热损伤的多重影响[14]。在煤炭地下气化数值模拟研究方面，国内外学者做了大量富有成效的工作。PARK等[15]利用非稳态的一维模型模拟了地下煤气化过程早期燃空区发育情况。模拟结果表明早期的燃空区发育受气化剂中氧气含量及煤炭自身气化的影响。BRITTEN等[16]提出在煤炭地下气化炉中的不同位置采用不同数学模型的数值模拟方法，但该模型不能预测燃空区的高温及产气的温度，也不能反映燃空区压力场的变化。KUYPER等[17-18]在BRITTEN模型的基础上添加了层流模型，但是仍不能反映温度场和压力场的变化。文献[19-20]通过模型试验研究了带有较大倾角的煤层不同温度场的分布变化情况，利用三维非线性数值模型分析了压力场的分布情况，最后得出了不同时间的温度变化情况。席建奋等[21]对顶板应力场扩展过程及顶板稳定性进行试验研究，得出氧化区培育阶段煤层温度场扩展主方向与裂隙方向一致，煤层温度场动态扩展与顶板热弥散的双重作用使顶板应力场的扩展速率逐渐趋近于煤层温度场扩展速率的结论。MEHDI等[22]利用FLAC3D模拟软件建立模型，预测煤炭地下气化过程煤岩层应力分布，发现由于热应力的作用，在工作面前缘和煤柱边缘的应力分布增大，并在模型基础上提出了煤地下气化-煤柱稳定性分析方法。上述研究成果对煤炭地下气化覆岩高温热损伤破坏特征研究具有重要指导意义。
单独对煤炭地下气化过程中的温度场或应力场进行模拟研究，虽得到了一定的成果和应用，但是煤炭地下气化过程反应区温度可高达千度，受高温影响，煤岩层的弹性模量、导热系数、恒压热容等参数都会发生较大的变化，应力场的改变又会对煤岩层的传热产生影响，相对于单一场模拟，热固耦合模拟结果在煤炭地下气化领域更具有实际意义。基于此，笔者在简化的煤炭地下气化燃空区扩展二维模型的基础上，通过拟合燃空区边界升温曲线，进行合理的热固耦合设置。相较于一般热固耦合模拟，采用在模型边界不同分段设置不同升温函数的方式，以间接的方法模拟出燃空区边界的移动，解决了在COMSOL数值模拟软件中边界移动较难实现的问题，得到更加接近试验真实情况的温度场、第一主应力场及垂直位移场分布情况。模拟得出的温度场，应力场及位移场等可更好地了解试验炉内围岩情况，以作为后续试验的参照，模拟结果有助于从热固耦合的角度研究地下气化过程对围岩的破坏。
1 数值模拟方案
1.1 煤炭地下气化简化模型
传热设计是整个热力耦合数值模拟方案设计的基础，合理地设计出与真实情况相符的固体传热过程对整个数值模拟有着重要的意义。在数值模拟设计中，根据热弹性本构方程、热弹性运动方程以及温度场控制方程[6]，将传热分为2个阶段进行研究，对垂直于气化通道的任一截面来说，气化通道截面的传热过程可分为升温阶段传热和降温阶段传热。
1.1.1 截面升温阶段
由于注气和点火部分都在煤层底部中点位置，在升温阶段，燃空区在截面不断在径向向上及向外扩展，如图1所示。在扩展过程中，燃空区的边界上不断进行煤的燃烧放热反应，热量向煤岩层传导，直到燃空区扩展出一个最大边界，在边界上的放热反应停止。任意选取一个垂直于气化通道的截面，此阶段对选取的截面来说是升温阶段。
图1 煤炭地下气化升温阶段燃空区扩展过程示意
Fig.1 Expansion process of UCG combustion zone during warming phase
1.1.2 截面降温阶段
随着火焰工作面的定向移动，处于火焰工作面后方的某一截面，当气化燃烧达到其最大扩展半径后，其截面边界所受的温度条件将发生变化，呈现逐渐降温的趋势。这一阶段对选取的截面来说是降温阶段(图2)。
图2 煤炭地下气化降温阶段燃空区扩展过程示意
Fig.2 Expansion process of UCG combustion zone during cooling phase
1.2 燃空区边界条件设定
取煤炭地下气化相似模拟试验模型一半进行数值模拟分析，如图3所示。
图3 反应炉内沿气化方向的煤岩层截面
Fig.3 Section of coal rock strata along gasification reaction direction in reactor
模型初始温度和环境温度均设置为293.15 K，并设定燃空区最大扩展半径为0.4 m。为方便进行传热边界条件设定，将燃空区边界分为圆弧边界和下边界，如图4所示。
图4 反应炉内截面模型燃空区边界设置
Fig.4 Boundary setting of burning area in reactor
根据对移动边界条件下传热规律的研究[23]，最大扩展位置升温可近似设置为一拟合分段函数，而下边界的传热较为复杂，以下将下边界平均分为8段进行分步研究，每段长5 cm，为段1～段8(图4)。圆弧边界的升温曲线拟合函数由下边界传热的分步研究得出。圆弧边界的升温曲线拟合函数近似于下边界段8的升温曲线拟合函数。
升温阶段下边界各分段依次升温，需逐步研究，而降温阶段下边界和圆弧边界作为已扩展出的边界，边界各点降温函数几乎相同。升温阶段传热较为复杂，降温阶段传热较为简单。
1.2.1 升温阶段温度设置
下边界各段函数的拟合分步进行，首先对下边界段1的升温函数进行设定，在这里将该段函数设定为一个以试验时间为自变量的线性函数，在此次研究中，根据以往试验情况，假设到达的最高温度为1 273.15 K，在正常情况下，试验炉内发生气化反应而形成2.8 m长的气化通道，试验用时约为5 d，下边界的分段方式为5 cm一段，则假设每段从内向外燃烧用时为2.4 h。则第一段下边界温度拟合函数设置为
T1=293.15+(980/8 600)t，
(0<t≤2.4×3 600)
式中：t为试验进行时间，s；T1为燃空区下边界段1在试验进行到t时刻的温度，K。
在对下边界段1的升温函数拟合完成后，在固体传热中仅对该段温度进行设置，函数设置全部使用分段函数，瞬态研究设定时间为0～2.4×3 600 s，即在试验开始到2.4×3 600 s的时间内进行模拟，每隔60 s取值一次。在2.4 h时，截面上燃空区边界的下边界已经扩展到了段1的末尾，此时对底板的传热状态如图5所示。
图5 反应炉内截面模型在2.4 h时刻传热结果(仅设置T1)
Fig.5 Heat transfer results obtained at 2.4 hours of section model in reactor (set T1 only)
由图5可知，段2中点温度约为1 138.15 K，依次记录段2至段8各段中点温度。可以看出对下边界第二段的升温，前2.4 h升温至1 138.15 K左右，由于再过2.4 h燃空区扩展到段2末端，即到4.8 h时，该段温度上升至1 273.15 K，所以将下边界段2的升温函数拟合为一个两段的分段函数：
T2=293.15+(845/8 640)t
(0<t≤2.4×3 600)
T2=1 138.15+(135/8 640)(t-8 640)
(2.4×3 600<t≤4.8×3 600)
式中,T2为燃空区下边界段2在试验进行到t时刻的温度，K。
在固体传热中对段2设置温度为T2，同时启用段1和段2的升温函数，进行瞬态研究，设定时间为0～4.8×3 600 s，即在试验开始到4.4×3 600 s的时间内进行模拟，每隔60 s取值一次。在4.8 h时，传热状态如图6所示。
图6 反应炉内截面模型在4.8 h时刻传热结果(仅设置升温函数T1，T2)
Fig.6 Heat transfer results obtained at 4.8 hours of section model in reactor (set T1,T2 only)
由图6可知，在4 h时，段3中点温度约为1 163.15 K，依次记录段3到段8中点温度。在第一次记录中的2.4 h的时间点，段3的中点温度为974.15 K，第二次记录在4.8 h的时间点，段3的中点温度为1 163.15 K，则将下边界段3的升温函数拟合为一个三段的分段函数：
T3=293.15+(681/8 640)t
(0<t≤2.4×3 600)
T3=974.15+(189/8 640) (t-8 640)
(2.4×3 600<t≤4.8×3 600)
T3= 1 163.15+(110/8 640) (t-2×8 640)
(4.8×3 600<t≤7.2×3 600)
式中：T3为燃空区下边界第3段在试验进行到t时刻的温度，K。
以此方法依次类推，得到各段中点温度，见表1。限于篇幅，此处仅给出拟合出的最后一段的升温曲线函数图形，如图7所示。分段函数T8近似于圆弧边界升温函数，T8同样适用于固体传热的圆弧边界温度。在此基础上设定好下边界每一分段及圆弧边界的升温函数，升温阶段传热研究的设定完成。
4)各段降温函数拟合。相对于升温阶段，降温阶段用时较长，对于试验炉的试验，截面降温大约需要100.8 h，并且降温阶段对热力耦合研究影响较小，可以将其理想化为一个简单的一元一次方程。即在原先升温阶段基础上，各分段函数从升至1 273.15 K开始，保持高温到19.2 h，在各分段函数后加函数段：
Ti=1 273.15-(980/362 880)×(t-8 640×8)
(19.2×3 600<t<120×3 600)
式中：Ti为燃空区下边界第i段在试验进行到t时刻的温度，K,i=1,2,3，…,8。
表1 燃空区下边界各段中点温度
Table 1 Midpoint temperature of each segment in lower boundary of burning area
时间/s温度/K段2段3段4段5段6段7段82.4×3 6001 138.15974.15884.15821.15767.15718.15668.154.8×3 600—1 163.151 007.15921.15854.15797.15733.157.2×3 600——1 170.151 022.15935.15860.15790.159.6×3 600———1 171.151 025.15933.15846.1512×3 600————1 173.151 022.15911.1514.4×3 600—————1 172.15994.1516.8×3 600——————1 137.15
图7 反应炉内截面模型燃空区下边界段8升温拟合函数
Fig.7 Fitting function of segment 8 in lower boundary of burning area in reactor
1.2.2 降温阶段温度设置
当最大燃空区边界处于降温阶段时，其传热边界条件为第一类边界条件(即给定边界温度)，时间设为19.2×3 600～120×3 600 s，即从19.2 h开始，到120 h为止，每隔10 min取值一次。
1.3 应力边界条件设置
设置截面模型为线弹性材料，且设定为几乎不可压缩。对称设置与传热部分相同(几何图形左边界)。固定约束为试验炉底部对内部岩体的约束(即几何图形的下边界)。在右边界设定辊支撑。由于正常情况下气化所处地下环境有上覆岩层压力，在气化试验煤岩层之上覆盖带孔铁板，按试验相似比例，在上边界添加边界载荷0.008 MPa，方向向下，并在材料的参数中，对导热系数、弹性模量进行变量设置，见表2。
表2 反应炉内各层材料物理参数设置
Table 2 Physical parameters setting of the layers of materials in reactor
材料导热系数/(W·m-1·K-1)弹性模量/GPa黏聚力/MPa热膨胀系数煤层0.2510-5T2-0.016T+6.53×10-6T2-0.004T+6.0372×10-5T2-0.007T+2.131底板及顶板2材料(0.979 18-0.000 757 69T+4.667 51×10-7×T2)×0.3110-5T2-0.024T+14.49-2×10-9T3+3×10-6T2-0.003T+6.5728×10-6T2+0.394顶板1材料(1.439 16-0.000 742 995T+2.133 25×10-7T2) ×0.2410-5T2-0.024T+14.494×10-6T2-0.006T+8.6546×10-6T2+0.519
注：T为模拟过程中模型中各点温度。
2 数值模拟结果
在COMSOL中，固体力学与固体传热在热应力和热膨胀的条件下同步计算，设置时间为0～19.2×60 min，即由试验开始模拟到19.2 h，每10 min取值一次，进行热力耦合计算。
2.1 温度场分布
温度分布情况如图8所示。
1)从图8a可以看出，由于点火点位于煤层中部靠底板的位置，所以在点火之后直到燃空区边界半径扩展至最大燃空区边界半径一半左右的这段时间，对顶板的影响较小，热量主要在煤层和直接底传播，且范围逐渐扩大，这时煤层和直接底出现部分温度较高区域是受采空区移动的燃烧边界的影响，该边界温度高达1 273.15 K。
2)随着时间的增加，受固体传热的影响，试验炉内各层靠近燃烧区部分的温度都出现了小幅度的升高。随采空区高温边界向外扩展，高温区域范围逐渐扩大，各层温度逐渐上升。
3)之后直到扩展出最大燃空区，采空区高温边界扩展到最大，在试验炉内各层的任意点温度都达到了气化试验的最大值。之后一直到试验结束，采空区域温度随时间线性降低，煤岩层温度也随之降低，直到接近环境温度。
图8 反应炉内气化通道截面在11、16、19 h时温度分布情况
Fig.8 Temperature distribution at 11,16，19 h in gasification tunnel section of reactor
2.2 第一主应变分布
1)在燃空区扩展的过程中，第一主应变主要产生于燃空区附近区域，主要集中在煤层和底板部分(图9)。
图9 反应炉内气化通道截面在6、12、19 h时第一主应变分布情况
Fig.9 The first main strain distribution at 6,12，19 h in gasification tunnel section of reactor
2)从图9b、图9c可以看出，燃空区边界扩展至最大扩展半径的一半之后，随着燃空区继续向外扩展，煤岩层受主应变影响的区域范围变化不大，但燃空区圆弧形边界附近的应变量持续增加。从试验开始到最大燃空区扩展完成，燃空区圆弧形边界最大第一主应变大小约为燃空区半径的1/10，而燃空区下边界最大第一主应变约为燃空区半径的1/20。
2.3 垂直位移场分布
1)在常温未开挖情况下，垂直位移较小，最大垂直位移出现在隔热顶板中，方向向下(图10)。
图10 反应炉内气化通道截面在6、12、19 h时垂直位移场分布情况
Fig.10 Vertical displacement field distribution at 6,12，19 h in gasification tunnel section of reactor
2)随着采空区扩展，温度逐渐上升，岩体受热膨胀，垂直位移整体呈向上趋势，且在模型中由下至上数值递增，由内向外数值递减。直接底内靠近采空区部分垂直位移较大，靠近采空区下边界中心部分垂直位移为煤层厚度的1/50左右。直接顶内最大垂直位移可达直接顶厚度的1/20。整体最大垂直位移出现在隔热顶板中部位置，位移约为隔热顶板厚度的1/12。
3 结 论
1)模拟得出了煤炭地下气化过程中燃空区附近的温度场，第一主应力场及垂直位移场的部分演化规律：由于点火点位于煤层中部靠底板的位置，所以在点火之后，直到燃空区边界扩展至最大扩展半径的一半的这段时间，热量主要在煤层和直接底传播，燃空区燃烧边界的移动是煤层和直接底出现部分温度较高区域的主要原因，之后直到扩展出最大燃空区，试验炉内各层的任意点温度都达到了气化试验的最大值。到试验结束，采空区域温度随气化时间的推移线性降低，煤岩层温度也随之降低，直到接近环境温度。
2)在燃空区扩展的过程中，第一主应变主要产生于燃空区附近区域。燃空区边界扩展至最大扩展半径的一半之后，煤岩层受主应变影响的区域范围变化不大，但燃空区圆弧形边界附近的应变量持续增加。从试验开始到最大燃空区扩展完成，燃空区圆弧形边界最大第一主应变大小约为燃空区半径的1/10，而燃空区下边界最大第一主应变约为燃空区半径的1/20。
3)在常温未开挖情况下，垂直位移较小，最大垂直位移出现在隔热顶板中，方向向下。随着采空区扩展，温度逐渐上升，岩体受热膨胀，垂直位移整体呈向上趋势，且在模型中由下至上数值递增，由内向外数值递减。直接底内靠近采空区部分垂直位移较大，靠近采空区下边界中心部分垂直位移在煤层厚度的1/50左右。直接顶内最大垂直位移可达直接顶厚度的1/20。整体最大垂直位移出现在隔热顶板中部位置，位移约为隔热顶板厚度的1/12。
由于COMSOL为有限元模拟软件，难以在软件中直接设置模型的边界移动，本次模拟未能直接实现燃空区边界随时间的移动，虽然对温度场进行了间接的边界移动拟合，但应力场仍受固定边界影响，所以模拟所得结果存在各点温度较实际偏高，应力相对较大的问题，后续在移动边界的热固耦合问题的研究上，希望能通过更加合理的设置和设计，得到更加符合实际情况的模拟结果。
[1] 钱鸣高,缪协兴,许家林.资源与环境协调(绿色)开采[J].煤炭学报,2007,32(1):1-7.
QIAN Minggao,MIAO Xiexing,XU Jialin.Green mining of coal resources harmonizing with environment[J].Journal of China Coal Society,2007,32(1):1-7.
[2] 余 力.我国废弃煤炭资源的利用：推动煤炭地下气化技术发展[J].煤炭科学技术,2013,41(5):1-3.
YU Li.Utilization of abandoned coal resources to promote the development of underground coal gasification technology in China[J].Coal Science and Technology,2013,41(5):1-3.
[3] 余 力,梁 杰,余学东.煤炭资源开发与利用新方法：煤炭地下气化技术[J].科技导报,1999(4)：33-35.
YU Li,LIANG Jie,YU Xuedong.Progress in the coal underground pneumatolysis technology[J].Science and Technology Review,1999(4):33-35.
[4] 梁 杰，余 力.反向两阶段煤炭地下气化方法的研究[J].煤炭学报,1996,21(1):68-72.
LIANG Jie,YU Li.Study of two-stage underground coal gasification in counter directions[J].Journal of China Coal Society,1996,21(1):68-72.
[5] 朱 铭,徐道一,孙文鹏,等.世界煤地下气化的快速发展[J].自然杂志,2012,34(3):161-166，180.
ZHU Ming,XU Daoyi,SUN Wenpeng,et al.Rapid progress of underground coal gasification in the world[J].Chinese Journal of Nature,2012,34(3):161-166，180.
[6] 辛 林,程卫民,谢 军,等.岩石单向加热热固耦合数值模拟研究[J].煤炭科学技术,2018,46(7):145-151.
XIN Lin,CHENG Weimin,XIE Jun,et al.Study of thermal-mechnical coupling numerical of rock under uniaxial heating[J].Coal Science and Technology,2018,46(7):145-151.
[7] XIN Lin,WANG Zuotang,WANG Gang,et al.Technological aspects for underground coal gasification in steeply inclined thin coal seams at Zhongliangshan Coal Mine in China [J].Fuel,2017,191:486-494.
[8] YANG Lanhe,LIANG Jie,YU Li.Clean coal technology：study on the pilot project experiment of underground coal gasification [J].Energy,2003,28(14):1445-1460.
[9] 刘淑琴,陈 峰,庞旭林,等.煤炭地下气化反应过程分析及稳定控制工艺[J].煤炭科学技术,2015,43(1):125-128.
LIU Shuqin,CHEN Feng,PANG Xulin,et al.Analysis on reaction process of underground coal gasification and stable control technique[J].Coal Science and Technology,2015,43(1):125-128.
[10] 刘淑琴,张 军,梁 杰,等.煤炭地下气化的综合利用前景[J].煤炭科学技术,2003,31(7):50-53.
LIU Shuqin,ZHANG Jun,LIANG Jie,et al.Outlook on comprehensive utilization of underground coal gasification [J].Coal Science and Technology,2003,31(7):50-53.
[11] 刘淑琴,张尚军,牛茂斐,等.煤炭地下气化技术及其应用前景[J].地学前缘,2016,23(3):97-102.
LIU Shuqin,ZHANG Shangjun,NIU Maofei,et al.Technology process and application prospect of underground coal gasification[J].Earth Science Frontiers,2016,23(3):97-102.
[12] WANG Zuotang,DING Xueping,HUO Liwen,et al.A remining technology of underground coal gasification at Zhongliangshan Coal Mine[J].Journal of Coal Science and Engineering(China),2008,14(3):469-473.
[13] CHANDELLE V,JACQUEMIN C,LÉTOLLE R,et al.Underground coal gasification on the Thulin site:results of analysis from post-burn drillings[J].Fuel,1993,72(7):949-963.
[14] 辛 林,程卫民,王 刚,等.煤炭地下气化多层热弹性基础梁模型及其应用[J].岩石力学与工程学报,2016,35(6):1233-1244.
XIN Lin,CHENG Weimin,WANG Gang,et al.Multi-layer thermo-elastic foundation beam model of UCG and its application[J].Chinese Journal of Rock Mechanics and Engineering,2016,35(6):1233-1244.
[15] PARK K Y,EDGAR T F.Modeling of early cavity growth for underground coal gasification[J].Industrial and Engineering Chemistry Research,1987,26:237-246.
[16] BRITTEN J A,THORSNESS C B.A model for cavity growth and resource recovery during underground coal gasification[J].In Situ,1989,13:1-53.
[17] KUYPER R A,VAN Der Meer T H,HOOGENDOORN C J．Turbulent natural convection flow due to combined buoyancy forces during underground coal gasification of thin seams[J].Chemical Engineering Science,1994,49:851-861.
[18] KUYPER R A,MEER T H,BRUINING J.Simulation of underground gasification of thin coal seams[J].In Situ,1996,20:311-346.
[19] YANG Lanhe.Study on the model experiment and numerical simulation for underground coal gasification[J].Fuel,2004,83:573-584.
[20] YANG Lanhe.Three - dimensional non - linear numerical analysis on the oxygen concentration field in underground coal gasification[J].Fuel Processing Technology,2004,85:1605-1622.
[21] 席建奋,梁 杰,王张卿,等.煤炭地下气化温度场动态扩展对顶板热应力场及稳定性的影响[J].煤炭学报,2015,40(8):1949-1955.
XI Jianfen,LIANG Jie,WANG Zhangqing,et al.Effect of dynamic expansion of temperature field on thermal stress field and stability of coal underground gasification roof[J].Journal of China Coal Society,2015,40(8):1949-1955.
[22] MEHDI N,SEYED M E J.Thermal-mechanical-numerical analysis of stress distribution in the vicinity of underground coal gasification (UCG) panels[J].International Journal of Coal Geology,2014,134/135:1-16.
[23] XIN Lin,CHENG Weimin,XIE Jun,et al.Theoretical research on heat transfer law during underground coal gasification channel extension process[J].International Journal of Heat and Mass Transfer,2019,142:118409.
