低变质煤微波热解数值模拟研究*

梁 坤 周 军， 吴 雷 周晶晶 张秋利， 宋永辉 兰新哲

(1.西安建筑科技大学化学与化工学院，710055 西安；2.陕西省冶金工程技术研究中心，710055 西安；3.西安建筑大学冶金工程学院，710055 西安)

0 引 言

煤炭资源是人类最主要的矿产资源，其最直接的利用方式是直接燃烧发电，然而，这种直接利用方式会造成严重的生态环境破坏[1]。因此，煤炭资源的高效清洁利用成为当下能源研究的重点内容之一[2-3]。低变质煤微波热解技术是微波技术与煤热解技术相结合的一种新型高效热解技术，具有加热速度快、选择性加热、控制简单、安全环保等优点[4-6]。虽然已有较多学者对煤的微波热解进行了大量实验研究，但研究工作多集中在热解产物收率及特性方面，对微波炉腔体内的电场和煤样的温度场分布的变化研究较少。微波腔体中电磁场分布是煤微波热解炉设计的关键，直接关系到微波热解炉的工业开发应用[7]，而煤样的温度场分布直接影响微波热解的能效[8]。数值模拟是一种可视化和量化的技术，借助计算机模拟，能够减少大量的中间试验，寻找最佳的操作参数，达到过程优化的目的[9-11]。CHEN et al[12]模拟了块状橡胶的微波加热行为，研究了微波频率、功率和样本大小对温度场分布的影响。结果表明，微波频率、功率和样本大小都对温度场分布存在影响，且在微波频率为2.45 GHz时，加热效率最高。块状橡胶内外温差随着功率的增加而增大；减小样品尺寸，温度分布更加均匀，但微波加热效率降低。ARSHAD et al[13]使用COMSOL软件模拟预测了空果束的微波加热行为，研究了生物质装载高度、生物质形状对温度和电场分布的影响，并对热点位置进行了预测。研究发现，终点温度很大程度上取决于样品的装载高度，模拟出的热点位置和实验观察到的位置几乎一致。

目前，有关煤微波热解数值模拟的研究较少，在仅有的研究中，研究者通常也将煤样的热容、介电系数及煤样化学反应热等关键影响因素简化为定值，使得模拟结果的准确性受到一定的影响[14]。笔者在煤样的TG-DSC分析结果的基础上，耦合电磁传热模型，研究了微波功率、热解时间和煤样位置对微波炉腔体内电场和煤样的温度场分布的影响，得出了数值模拟的最优工艺条件，并在此工艺条件下进行了低变质煤微波热解实验研究。此模型考虑了煤样化学反应热对电场和温度场分布的影响，为探索微波热解炉工业化放大难题的解决提供了新的研究思路。

1 模型建立及实验过程

1.1 几何模型和网格划分

基于COMSOL Multiphysics模拟软件对微波热解装置建立几何模型，如图1所示。波导为矩形端口，位于炉腔左侧，炉腔和波导由铜制成，内部充满空气(零介电损耗)，圆柱形煤样位于腔体中心位置。为了减小计算机运行负担，采用了对称模型。详细的几何参数如表1所示。

图1 微波加热几何模型Fig.1 Geometric model of microwave heating

表1 模型几何参数Table 1 Model geometry parameters

对建立的几何模型进行有限元网格划分，并对有限元网格划分进行质量计算，如图2所示。几何模型是由9 141个域单元、161个边界单元和22个顶点单元组成的完整网格，划分的网格元素中95%以上都满足质量大于0.6的要求。因此，模拟计算结果是可靠的，可以计算出准确的模拟结果[15]。

计算过程在Intel Core i7-9750H，3.40 GHz处理器，16 GB RAM内存和64位Windows 10 Enterprise操作系统的计算机上进行，模拟软件为COMSOL Multiphysics 5.4版本。

1.2 控制方程

1.2.1 电磁场

电磁场的计算运用麦克斯韦方程组进行精准求解，计算方程如式(1)所示。

图2 几何模型有限元网格Fig.2 Geometric model finite element mesha—Mesh map；b—Mesh quality map

(1)

式中：μ为介质磁导率，H/m；E为电场强度，V/m；ω为电磁波角频率；ε0为真空介电常数，8.85×10-12F/m；εr为相对介电系数；k0为自由空间电磁波波数；j为电流密度，A/m2；σ为电导率，S/m。

1.2.2 传热方程

传热只在煤样区域进行求解。利用傅立叶能量平衡方程对微波场中的热传递进行计算，计算方程如式(2)所示。

(2)

式中：ρ为介质密度，kg/m3；Cp为介质热容，J/(kg·K)；t为时间，s；k为介质热导率，W/(m2·K)；u为速度向量，m/s；T为温度，K；Q为供给煤样升温的总热量，J。

根据JERZY et al[16]对煤热解过程的热容变化研究结果，将Cp定义为式(3)的温度函数，引入计算模型。

Cp=1.01+4.19×10-3×T-2.01×10-6×
T2-5.63×10-9×T3+3.85×10-2×T4

(3)

式中：T为温度，K。

计算过程中，将煤样热解发生反应消耗的热量及热解气逸出带走的热量定义为Qloss，在计算煤样温度变化时，应从微波热源所提供的总能量中减去该部分热量。Qloss可根据热失重实验过程中获取的DSC差热曲线进行函数拟合得到，如图3所示。拟合后Qloss计算如式(4)所示，将方程式引入计算模型。

图3 DSC差热分析拟合函数曲线Fig.3 DSC differential thermal analysis fitting function curve

(4)

式中：T为温度，K。

1.3 边界条件

对于电磁场，波腔和波导的壁被定义为阻抗边界条件并在频域中求解，计算方程如式(5)所示。

(5)

式中：n为折射率，单位为1。

对称边界被描述为完美磁导体，如式(6)所示。

n×H=0

(6)

式中：H为磁场强度，A/m。

端口边界由横向电波激发，端口为矩形，频率为2 450 MHz，在TE10(m=1，n=0)模式下运行，计算方程如式(7)和式(8)所示。

(7)

(8)

式中：β为传播常数；c为传播速度，m/s；v代表微波频率；vc为截止频率；c′为电磁波在波导中的传播速度，m/s；m和n是模式编号；a和b表示矩形横截面的长度，m。

热绝缘边界条件是指在边界上没有热流，计算方程如式(9)所示。

(9)

1.4 实验过程

实验所用煤样为陕西北部某低变质煤，其工业分析及元素分析见表2。称取50 g低变质煤样，放入石英反应器中，根据实验条件连接测试仪表及石英反应器放置位置(以煤样处于微波炉腔体中心位置为0 cm位置，依次水平靠近波导以3 cm增幅点分别记为-3 cm位置、-6 cm位置、-9 cm位置，依次水平远离波导以3 cm增幅点分别记为+3 cm位置、+6 cm位置、+9 cm位置)，连接完成后检测装置的气密性，并将反应器及管路内气体排空。设置实验所需微波功率，开始热解实验，达到设定的热解时间后自动停止加热，待反应器冷却至室温。反应结束后，称量石英反应器(含有固体产物)质量和冷凝装置下集液瓶(含有液体产物)质量并记录，利用差重法计算固、液产物产量和收率，液体产物利用有机溶剂和无水硫酸钠进行油水分离，气体产物用气体分析仪进行分析检测。低变质煤微波热解产物收率计算公式如式(10)～式(14)所示。

表2 实验用煤样的工业分析与元素分析(%1))Table 2 Proximate and ultimate analysis of coal(%1))

1) Mass fraction; 2) By difference.

ms=m1-m2

(10)

(11)

mL=m3-m4

(12)

(13)

(14)

式中：m0为原煤质量，g；m1为热解结束后石英反应器(含有固体产物)质量，g；m2为空石英反应器质量，g；m3为热解结束后集液瓶(含有液体产物)总质量，g；m4为热解前空集液瓶总质量，g；m5为油水分离后焦油质量，g；ms为固体产物产量，g；ws为固体产物收率，%；mL为液体产物产量，g；wL为液体产物收率，%；wtar为焦油收率，%。

微波热解模拟部分采用单因素实验法，分别研究微波功率(280 W,420 W,560 W,700 W)、热解时间(100 s,200 s,400 s,1 800 s)和煤样位置(0 cm，±3 cm，±6 cm和±9 cm)对微波炉腔的电场分布和煤样的温度场分布的影响，确定出低变质煤微波热解模拟最优实验条件，并在模拟最优工艺条件下进行低变质煤微波热解实验，通过模拟升温曲线与热解实验升温曲线的对比，以验证数值模拟的准确性。

2 结果与讨论

2.1 微波功率对微波炉腔体内电场和煤样温度场分布的影响

在四种功率下的模拟结果如图4所示。由图4a可以看出，在微波炉腔体内，四种微波功率均出现了7处强电场分布区域，其中6处以炉体中心为界呈对称分布，并且越靠近波导位置，电场强度越强。随着微波功率的增大，微波腔体中电场分布位置没有明显的变化，但强度随微波功率的增大而增大。当微波功率从280 W增加到700 W时，靠近波导的电场强度从约4.2×104V/m增加到约6.5×104V/m，增加了约1.55倍。煤样在微波场中吸收微波能转化为热能，而温度又是热能的直观表现形式。煤样的温度分布如图4b所示，煤样的高温区主要集中在煤样的右上部分，并朝煤样左下部分逐渐降低。这是因为煤样处在微波炉腔中心位置时，其右上部分处于+3 cm位置的强电场辐射范围内，而左下部分距离较强电场分布区域较远，因此导致煤样中温度分布不均匀。煤样温度分布呈羽状分布，且随着微波功率的增大，羽状面积逐渐增大。当微波功率为700 W时，煤样温度的羽状分布面积达到煤样切面的80%以上，大部分羽状区域温度超过800 ℃。微波功率增大使得微波炉腔体内电场增强，增大了煤样在微波场中的吸波放热，从而导致煤样温度分布增大，温度整体升高。在本研究中，微波功率为700 W时，煤样的高温区域面积最大，热解终温最好，因此，选择微波功率700 W为低变质煤微波热解的较优微波功率。

图4 不同功率下电场和温度场分布Fig.4 Distribution of electric field and temperature field at different powersa—Electric field；b—Temperature field

2.2 热解时间对微波炉腔体内电场分布和煤样温度场分布的影响

在微波功率700 W，煤样处于微波炉腔体中心位置(0 cm)的条件下，研究热解时间对微波炉腔体内的电场分布和煤样的温度分布的影响，结果如图5所示。由图5a可知，当微波热解时间从100 s增加到1 800 s时，微波炉腔体内的电场分布和强度几乎没有变化，都出现了6处呈对称分布的强电场，强电场的强度约为6.5×104V/m。微波功率、微波腔体、腔体内吸波物质等因素会影响微波炉腔内电场分布特性，而微波热解时间只是影响微波热解的一个因素，因此，热解时间对电场分布特性几乎没有影响，但对煤样温度分布影响十分明显。由图5b可以看出，煤样的温度分布主要集中在煤样的右上部分，但随着热解时间延长，温度分布范围越来越大，温度也越来越高。从热解时间100 s增加到400 s时，煤样的最高温度从730 ℃升高到1 100 ℃。正如前文所述，煤样的温度分布主要是热解反应热和煤样吸波转化热造成的。虽然煤样所在微波腔体内电场强度没有变化，但热解反应随着热解时间的增加而越发剧烈，因此，延长热解时间对煤样温度的分布和高低影响十分明显。当微波热解时间增加到1 800 s时，煤样的温度分布与热解时间400 s时的温度分布和最高温度相似，但整体温度略低于400 s时的温度。为了使微波热解充分反应，结合之前的研究经验，选择微波热解时间1 800 s为最优的微波热解时间。

图5 不同热解时间下电场和温度场分布Fig.5 Distribution of electric field and temperature field at different pyrolysis timea—Electric field；b—Temperature field

2.3 煤样位置对微波炉腔体内电场分布和煤样温度分布的影响

当微波功率为700 W，热解时间为1 800 s，煤样处于0 cm，±3 cm，±6 cm和±9 cm七个位置时微波炉腔的电场分布和煤样的温度分布如图6所示。由图6a可以清晰看出，当改变煤样位置时，微波炉腔中的电场分布和强度也会随之发生变化。当煤样处于-9 cm位置时，除波导位置外，微波炉腔内电场没有出现较强的分布区域，这可能是由于煤样太靠近波导，除小部分微波被煤样吸收，大部分微波被煤样阻挡反射，导致微波炉腔内其他地方电场分布强度较弱。当煤样处于-6 cm位置时，微波炉腔体内在-9 cm位置出现了一个强电场分布区域，强电场分布区域中心位置强度约为1.5×104V/m，其他区域电场强度较煤样处于-9 cm时电场分布均匀，电场强度也有所增强。当煤样处于-3 cm位置时，微波炉腔内明显出现了6个强电场分布区域，且呈对称分布形态，与煤样处于-6 cm位置相似，同样在电场分布区域-9 cm靠左位置出现了最强电场，电场强度约为1.1×105V/m，较煤样处于-6 cm位置时相同位置的电场强度提高了一个数量级。其他位置的电场强度增强，强电场分布区域面积增大、形状规则。 当煤样处于0 cm位置时，与煤样位置处于-3 cm最强电场相同，微波炉腔体内出现6处强电场分布区域，最强电场分布区域在-9 cm位置；与之不同的是，强电场分布区域面积增大，但电场强度有所降低，最强电场区域的最大电场强度约为6.2×104V/m。当煤样处于+3 cm位置时，微波炉腔体内强电场分布区域数量上与之前一样，存在6个强电场分布区域，但不同的是，强电场分布区域位置有所差异，从之前的竖向分布变为横向分布，分布区域面积有所减小，但强电场分布区域的电场强度较为均匀。当煤样处于-6 cm和-9 cm时，强电场分布区域较为相似，分布区域位置大致在一条线上，形状无规则，也均在-9 cm靠左位置出现了最强电场分布区域，最强电场分别约为3.5×104V/m和5.2×104V/m。从电场分布和强度来看，当煤样位置处于-3 cm到+3 cm时，微波炉腔内强电场分布呈规则对称分布，强电场分布区域面积较大，电场强度也较高，尤其是煤样处于-3 cm位置。

由图6b可以看出，煤样位置显著影响了煤样的温度分布，这是因为煤样位置影响电场分布和强度，而电场变化又直接影响煤样的吸波转化热。当煤样处于-9 cm位置时，煤样温度分布主要集中在煤样的左侧中心位置，呈核壳型分布，中心最高温度为590 ℃，由电场分布图可知，煤样左侧处在一个弱电场分布区域，因此，煤样的温度集中在左侧位置。当煤样位置处于-6 cm时，煤样的温度分布区域面积增大，呈圆环型分布，整体温度较高，最高温度达到860 ℃，但煤样中心位置温度较低，仅约为470 ℃，与煤样左侧的最高温相差近400 ℃。由电场分布图同样可以看出，煤样在-6 cm位置时，煤样左侧处在一个强电场分布区域，而右侧处在一个弱电场分布区域，因此，煤样温度分布呈圆环型分布，且煤样左侧温度高于右侧温度。当煤样处于-3 cm位置时，煤样的温度分布最为均匀，温度也最高，煤样整体温度约在1 000 ℃。尽管煤样在-3 cm时，煤样整体处于强电场的边缘区域和弱电场区域，但强电场边缘区域和弱电场区域的电场强度分别约为6×104V/m和2.5×104V/m，均高于煤样处于-6 cm位置时的电场强度，这也导致了煤样温度高于之前所处位置煤样温度，温度分布也较为均匀。前文已经介绍了煤样处在中心(0 cm)位置的微波炉腔体的电场分布和煤样温度分布变化规律和原因，此处不再重复介绍。当煤样处在+3 cm位置时，煤样温度分布主要在煤样的底部，此处煤样的最高温度只有约530 ℃。对照此处电场分布可知，煤样处在+3 cm位置时，没有处于任何强电场区域，此处的煤样吸波转化热较低。当煤样处于+6 cm位置时，煤样温度分布主要集中在煤样上半部，分布较为均匀，最高温度到达940 ℃。当煤样处于+9 cm位置时，煤样分布在左上部区域，但由于所处位置远离强电场分布区域，煤样整体温度较低。由以上分析可知，煤样作为弱吸波介质，会使微波炉腔体内的电场分布发生改变，当煤样处于强电场分布区域时，煤样的温度分布均匀，最高温度较高。因此，选择煤样处于-3 cm位置为低变质煤的最优热解反应位置。

图6 不同煤样位置下电场和温度场分布Fig.6 Distribution of electric field and temperature field at different coal sample locationsa—Electric field；b—Temperature field

2.4 实验验证及产物分析

为了保证模拟的准确性，在微波功率700 W，加热时间30 min，煤样位置处于-3 cm位置的条件下，对低变质煤微波热解的模拟热解温度进行了验证实验，实验采用与模拟模型相同尺寸的实验装置，结果见图7。如图7所示，模拟热解温度曲线较实际测得热解温度曲线平滑，没有出现温度波动，在热解升温阶段，模拟热解温度曲线与实际热解温度曲线吻合度较高，但热解终温存在差异，模拟终温均高于实际监测终温。在煤样热解温度的升高过程中，煤样放出的热量主要由煤样热解的化学热和煤样吸波的转化热构成，在热解升温过程中，煤样热解的化学热起主要作用，但随着化学反应的完成，煤样吸波的转化热逐渐发挥作用。然而，随着热解温度的增高，煤样的介电损耗正切值呈现先增大后减小的变化，而本模型主要考虑了化学反应热，假定煤样的介电损耗是一个常数，因此，出现了升温阶段实验和模拟吻合度较高，而热解终温存在差异。在接下来的研究中将进一步修正模型，考虑煤样介电损耗与热解温度的关系，以进一步优化模型。

图7 实验升温曲线和模拟升温曲线Fig.7 Experimental and simulated heating curves

对低变质煤微波热解最优条件(模拟所得)产物进行了较为详细的表征与分析。利用差量法计算微波热解四种产物的收率，结果如图8所示。由图8可知，低变质煤在微波场中吸收微波能，转化为热能，发生了热分解反应，产生了兰炭、热解水、焦油和煤气等产物，其收率分别为70.25%，5.31%，3.99%和20.45%。在微波热解过程中，热解水是由结合水、吸附水和反应生成水三部分组成，其中反应生成水占主要部分。表3所示为兰炭的工业分析和元素分析，与表2中数据对比可知，兰炭中灰分增大，而挥发分显著降低，这是由于煤中灰分在高温热解过程后，依然存在于固体产物兰炭中，而挥发分由于高温作用从煤中逸出，造成了兰炭中灰分占比增大，而挥发分占比减少。兰炭的SEM照片(见图9)也能证明这一说法。由图9可以看出，兰炭表面存在大小不一的孔洞和裂纹，这些孔洞和裂纹就是挥发分、煤气和煤焦油气在逸出过程造成的。兰炭中氢元素和氧元素显著降低，是由于在微波热解过程中，煤中的氢元素由于高温和微波作用，产生了大量的自由氢，自由氢在热解过程中发生加氢热解反应，使得煤中的羟基、羰基等含氧官能团断裂，生成了焦油、热解水和煤气。

图8 低变质煤微波热解产物分布Fig.8 Distribution of microwave pyrolysis products of low-rank coal

表3 兰炭的工业分析与元素分析(%1))Table 3 Proximate and ultimate analysis of bluecoke(%1))

1) Mass fraction; 2) By difference.

图9 固体产物兰炭的SEM照片Fig.9 SEM image of solid product bluecoke

按照焦油中单个成分的碳原子数，将焦油成分划分为三类：碳原子数小于10的成分称为轻质油，碳原子数在10～20之间的成分称为中质油，碳原子数大于20的成分称为重油。微波热解所得焦油的组成见图10。由图10a可知，低变质煤微波热解的焦油成分组成十分复杂，焦油中大部分成分在15 min～38 min被检测出来，也就是说焦油中的主要组分是轻质油和中质油。对GC-MS分析的焦油组分进行统计分类，如图10b所示。焦油中轻质油、中质油和重油成分分别占37.43%，59.58%和2.99%，轻质油和中质油的质量分数占焦油总量的97.01%，统计结果与GC-MS光谱一致。低变质煤微波热解煤气组分如表4所示，煤气中氢气体积分数达到35.65%，甲烷体积分数达到17.37%，这两种“富氢”气体体积分数占煤气总体积的53.02%，高价值气体总量超过了70%，煤气的热值约为14.70 MJ/m3。

图10 微波热解所得焦油组成Fig.10 Composition of tar from microwave pyrolysisa—GC-MS spectra；b—Classification of tar components

表4 微波热解所得煤气的主要成分Table 4 Main components in gas from microwave pyrolysis

3 结 论

1) 在构建模型时考虑了煤样热容和煤样化学反应热的影响，模拟结果与实验吻合度较高。

2) 微波功率和煤样位置对腔体内电场和煤样温度场分布影响十分明显，而热解时间对电场分布影响较小，对煤样温度分布影响较大。微波腔体中电场强度、煤样的升温速率和热解终温都随微波功率的增大而增大。当煤样处于-3 cm位置时，微波炉腔体内电场分布较为均匀，电场强度较高，煤样的温度分布较为均匀，处于高温区的面积最大。

3) 模拟的最优低变质煤微波热解条件为：微波功率700 W、热解时间1 800 s和煤样处于-3 cm位置。在此条件下进行煤的热解实验，其产物兰炭、焦油、热解水和煤气的收率分别为70.25%，3.99%，5.31%和20.45%，其中，焦油中轻质油含量为37.43%，煤气中“富氢”气体组分含量为53.02%。