块状木材的精确燃烧模拟在建筑消防和森林消防领域具有巨大应用潜力。在建筑消防领域，通过块状木材燃烧模拟，可以评估不同建筑结构和材料在火灾事件下的响应和性能，预测火灾在建筑内部的蔓延速度和路径，为建筑物内的疏散设计提供科学依据。在森林消防领域，通过树木燃烧模拟，可以深入了解森林火灾蔓延规律和火线扩展趋势，预测火场范围和火势大小，对扑灭林火和制定森林火灾防控策略具有重要意义。
然而，FDS工具的网格划分手段与主流计算流体力学(computational fluid dynamics，CFD)软件相比存在明显差距，难以精确描述弧形表面，而且热解模型较为单一，预测结果偏差较大。
现有的热解模型大致可分为3大类：单组分单反应模型、多组分平行单反应模型和详细反应模型。单组分单反应模型是生物质热解计算中最简单的方法，提供了基于单一全局反应的对热解热释放的简单描述[7]。多组分平行单反应模型可描述多种生物质组分的脱挥发反应。与单组分单反应模型类似，该模型可以根据环境条件预测质量损失率，但未考虑生成物分布[4-5, 8-9]。大多数研究将生物质简化为3种组分——纤维素、半纤维素和木质素，并假设组分之间没有相互作用[5, 7-10]。因此，生物质热解可以表征为各组分脱挥发贡献的总和。丁彦铭[11]通过开源CFD软件包OpenFOAM对粉末状木材热解燃烧过程进行了模拟计算，使用了多组分平行热解模型详细阐释了木材的热解过程，并对热解燃烧过程进行求解计算。2008年，Ranzi等[6]提出了木材详细反应模型，该模型包含15个反应和20多种挥发性化合物。Ranzi等不断完善和增强了该模型[12-14]。Anca-Couce综述了木质纤维素生物质热解的研究进展，对生物质热解过程中的二次反应进行了讨论[15]。详细反应模型可以更全面地描述所涉及的化学反应，考虑了竞争反应、平行反应以及二次反应，但该类模型涉及更多的动力学参数，计算成本很高。
综上可以发现，现阶段对块状木材燃烧的模拟研究仍不够充分，缺少基于Fluent的木材燃烧模拟方法。本研究将基于木块引燃实验结果进行木块引燃模拟，在Fluent中复现木块引燃过程，并与实验结果进行比对，从而提出基于Fluent的木块引燃模拟方法。
1 木块引燃实验
使用体积较小的柱状木块可减小挥发物释放对整体燃烧情况的影响(挥发物产生量较低)，也可有效提高模拟计算效率，因此，本研究选用小尺寸柱状木块进行实验。
1.1 样品物性参数实验测定
选取编号3的木条，取一部分研磨成粒径小于0.18 mm的木屑，称量1.0 g的样品进行实验。首先，将样品放入热重分析仪中加热至105 ℃，维持该温度直至样品质量不再发生变化，记录此时木屑质量m1=0.911 3 g。然后，快速将气氛温度升高至500 ℃，并维持该温度不变，记录木屑最终剩余质量m2=0.232 4 g。
基于上述热重分析仪得到的数据，对木材的各组分进行初步定量分析。水分质量比例为(m-m1)/m=8.87%，与表 1中水分测试仪测得的编号3的木条含水率数据(8.68%)差距不大；在500 ℃实验环境中热解的最终转化率为(m1-m2)/m1=0.745。上述燃烧环境与开放环境差异较大。在开放环境下，受环境因素(温度、湿度等)影响，木材样品可能无法全部转化为焦炭进行固相燃烧，也就是说，木材初次分解可能会提前终止。
1.2 引燃实验设置
这里选取编号6的木条，切割成30 mm长的柱状木块，继续对木块进行表观密度和含水率测定，结果见表 2。可以看出，来自同一木条的木块同样存在表观密度和含水率差异，但差距较小。
表 2 取自同一木条的木块表观密度和含水率测定结果 |
| 样品编号 | 表观密度/(kg·m-3) | 含水率/% |
| 6.1 | 449.71 | 8.76 |
| 6.2 | 450.58 | 8.84 |
| 6.3 | 460.27 | 8.68 |
| 6.4 | 460.27 | 8.42 |
本研究采用正庚烷火对木块进行引燃，使用边长为30 cm、高度为5 cm的方形油盘。具体实验场景见图 2。木块轴线距离油盘底部45 cm，木块通过螺纹杆与电动机相连。实验过程中火焰完全包裹木块(创造低氧环境)，木块表面仅有微弱火焰出现，可有效避免氧气对木块热解过程的影响。电动机提供机械能，带动木块在燃烧过程中以稳定转速(5 r/min)旋转，以保证木块表面受热均匀。点燃油盘持续引燃木块，使用电子天平(精度0.1 g)记录木块的质量变化，直至热重(thermogravimetric，TG)数据维持稳定，实验结束。
1.3 实验结果分析
经过多次重复实验后发现，在有限的引燃时间内(约2 min)，木块的质量损失比例大致在20%~50%之间，推测该差异与木块样品所处生长环境(表现为相似含水率下不同的表观密度等)及引燃过程中的开裂程度等因素有关，具体原因较为复杂。受限于研究成本，这里并未进行细致探究，仅选取一块典型样品进行分析。
基于实验数据可大致将木块的引燃过程分为3个阶段：蒸发阶段(Ⅰ)、热解阶段(Ⅱ)和质量缓慢下降阶段(Ⅲ)。蒸发阶段持续时间较长，质量下降速率较慢，该阶段只存在少量生物质热解(木块表面区域和螺纹连接区域的浅层热解)；热解阶段持续时间较短，质量下降速率较快。从热解阶段的质量快速下降现象推测，该阶段木块热解区域的温度较高。当热解进展到一定程度后，木块质量较为稳定，整个过程进入质量缓慢下降阶段。该阶段划分可为复现热解过程中的木材内部温度分布提供重要参考。
实验结束后，木块表面覆盖了一层较为均匀的炭化层，且出现了皲裂现象，如图 5所示(木块切割过程中存在炭化层脱落现象，木块截面图中的炭化层厚度并不均匀)。与螺纹杆连接处也存在一层炭化层，但厚度可忽略不计。可以推测，木块的炭化还受到介质气密性的影响。木块表面直接接触空气，在热解过程中会产生裂隙，导致气密性较差，水蒸气、热解气体得以逃逸，从而使表面持续碳化；然而，炭化层还存在一定的气密性，当炭化层积累到一定厚度，内部产生的气体无法逃逸，反应过程无法快速进行，整个过程进入质量缓慢下降阶段。
基于上述分析可得，木块各区域的干燥、转化情况不仅依赖温度条件，还受气密性影响：当热解产物无法排出时，蒸发、热解反应也会终止；随着木块气密性降低，炭化层的转化率也随之降低。实验测量得到的木块表面各处炭化层的厚度均在5 mm左右。基于热重数据，可计算炭化层的平均转化率为19.14%。平均转化率数值明显较低，据此可以推测炭化层内部存在转化率梯度。对炭化层表面和内部区域分别采样后进行热重测试，得到的结果如图 6所示。可以看出，木块炭化层表面区域的转化率明显高于其内部区域，证明炭化层内部存在转化率梯度。
2 木块引燃模拟与修正
在本研究中，木材微元的能量释放、传质传热以及孔隙率变化等过程将使用用户自定义函数(user defined function，UDF)复现。
2.1 木材组分简化及模拟环境设置
木材主要由水分、提取物、木质素、纤维素、半纤维素和灰分组成。随着温度的升高，木材中的3种聚合组分(木质素、纤维素、半纤维素)开始热解为挥发性气体、焦油和焦炭的混合物。该热解通常被认为是单个组分分解机制的叠加：半纤维素首先分解(180~350 ℃)，然后是纤维素(275~350 ℃)和木质素(250~500 ℃)[16]。焦油在热解过程中又会分解为挥发性气体和焦炭(会有少量焦油留存，这里忽略不计)。确定上述组分比例对复现木块热解过程具有重要意义。
然而，木材每种聚合组分形成的焦炭和挥发性气体比例难以通过实验进行测定。鉴于此，本研究忽略了木材中的灰分(小于1%)，将木材的3种聚合组分简化为一种有机物组分，即将木材微元简化为水分和有机物。
火源温度场对木块引燃过程至关重要。考虑到多种CFD软件中对于油池火的模拟采用气体燃烧的方式进行，本研究采用正庚烷气体进行组分传递与燃烧反应计算，以复现实验中的引燃过程。本研究基于实验中获取的正庚烷燃烧速率，在ANSYS Fluent求解器中对气态正庚烷燃烧进行模拟，调整计算域边界条件使计算结果与实验结果之间的差值缩小到可接受范围，从而精准量化火源温度场。
此次模拟所用的几何模型和局部加密后的网格如图 7所示。木块区域的网格尺寸设置为2 mm，气相燃烧区域设置为20 mm。本研究发现，木块表面边界层区域和炭化层区域的局部网格尺寸对计算结果影响显著，经过多次尝试，定义边界层首层网格高度为1×10-5 m，增长速率为1.2，共设置12层；炭化层区域内首层网格高度同样为1×10-5 m，增长速率为1.1，共设置40层。
本研究基于编号6.3的样品引燃实验结果，使用表 2中对应的平均表观密度(460.27 kg/m3)和平均含水率(8.68%)进行木块引燃模拟。经过研究发现，不锈钢螺纹杆的传热对TG曲线的下降以及炭化层内部的温度分布影响微弱，在模拟中忽略不计。
除此之外，本研究还使用滑移网格方法实现木块旋转。木块表面区域网格较为密集，若在边界处生成滑移面易导致网格重叠，造成计算崩溃。因此，将滑移面设置在外部环境中：木块区域设置为固体，滑移面与木块表面形成的闭合区域设置为流体。设置最外层区域旋转，并设置重力方向随时间变化，间接实现木块旋转。
然而，挥发性气体的燃烧非常复杂，模拟难度较高。在模拟研究中，本文作者发现，如果在木块表面区域添加木材干燥、热解过程中释放的气体(简化组分)，易导致木块两侧(圆柱底面)受热不均匀，表现为圆柱表面温度高、底面温度低，与实验情况不符。基于正庚烷火焰强度及其营造的高温低氧环境，以及木块较小的体积占比，本研究忽略木材在引燃过程中的气体释放，仅考虑正庚烷气相燃烧，进行大涡模拟(large eddy simulation，LES)。
2.2 干燥及热解过程模拟方法
在燃烧过程中，一阶Arrhenius表达式决定了相应温度下计算单元内水分、挥发物的释放速率以及炭的生成速率。释放后的热解组分在气相中燃烧形成火焰，从而进一步加热木材，实现燃烧循环。Spearpoint等[26]通过实验研究了这一机制，并对燃烧速率进行了监测。
整体反应的模拟流程如图 8所示。木材微元受到加热后，会先进入干燥阶段释放水蒸气，在温度达到一定阈值后开始热解，生成挥发物、灰分和焦炭。由于灰分含量非常微小，在模拟过程中已将其忽略。
木材内外的能量传输的模拟至关重要。本研究考虑火焰辐射影响，将离散坐标(discrete ordinate，DO)辐射模型与能量方程相结合，对木材内外的气相辐射基于灰色气体假设进行模拟。综合考虑化学反应和热解放热引起的其他额外热源以及水的气化潜热，相关公式总结如下：
其中：Es和Ef分别是固体总能量和流体介质总能量，γ为介质孔隙率，Ji为组分i的扩散通量，keff是介质的有效导热系数，Sfh为流体焓源项，Yj为组分j的质量分数，hj组分j的显焓，kf和ks分别为流体和固体的导热系数。
本研究所使用的UDF代码最初由Andersen和Heinesen开发，用于锅炉中的生物质颗粒燃烧[27]。本文作者对初始代码进行了一定修改，采用了更为精确的求解策略，以复现块状木材干燥、热解过程。为简化计算，本研究采取的假设如下：
1) 木材各转化子过程反应速率均受一阶Arrhenius表达式的约束；
2) 木材中的灰分含量为零。
首先，在计算开始前，需要对木材的初始质量等物理量进行定量，并划分用户自定义内存空间(user defined memories，UDM)。UDM可存储不同组分质量的值、各类源项的值以及后处理中所需要的各类物理量的值。将物理量的值赋于UDM可以简化代码，提高计算效率，方便后处理。
然后，将木材粉末热重分析和热重-差示扫描量热(TG-differential scanning calorimetry，TG-DSC)测试得到的关于热解速率和吸、放热速率的拟合公式植入UDF，用以描述木材微元的热解过程。
为保证计算顺利进行并方便重复计算，使用DEFINE_ON_ DEMAND宏定义一个按需执行的命令reset_solid，用于初始化UDM；使用DEFINE_PROFILE宏定义固体介质密度等物理特性；使用DEFINE_SPECIFIC_HEAT宏定义固体介质的比热容；使用DEFINE_PROPERTY宏定义固体材料密度、不同方向的热导率，以及吸收率和发射率；使用DEFINE_EXECUTE_ AT_END宏定义一个在计算过程中可被调用的命令execute_at_end，在每个时间步长迭代结束时执行，以实现对UDM进行更新。
2.3 木块引燃模拟及其可参考性检验
本研究在模拟过程中发现，如果直接设置火焰与木块表面接触，升温速率明显较高，会导致TG曲线提前进入热解阶段。推测其原因是蒸发阶段木块产生的水蒸气在其表面形成一层隔热层，具有一定初速度的水蒸气隔绝了部分火焰，导致升温速率较慢。随着木块内部温度的上升，水蒸气的温度也会上升，水蒸气降低升温速率的能力会下降。
基于此，本研究在木块表面设置了一个虚拟换热面，并根据实验结果进行不同厚度设置，以修正木块表面升温速率。虚拟换热面的详细参数设置见表 3。
表 3 虚拟换热面参数设置 |
| 时间范围/s | 厚度/m | 材料 |
| [0, 32.2] | 0.06 | wood-pine |
| (32.2, 76] | 0.007 | (GRANTA MDS Database) |
| >76 | 0 |
木块微元释放的热量并不会被固体介质全部吸收，从微观角度看，孔隙率越高的区域，被挥发物带走的热量越多。经过多次尝试，本研究设置后续模拟所使用的热解放热阶段的热量计算公式为
其中: qeff为有效吸/放热量，qpy为热解吸/放热量。
考虑到木块炭化层内存在转化率梯度，且编号6.3的样品炭化层区域平均转化率为19.14%，设定木块微元初次分解的终止点为0.75，未炭化区域转化率为0。本研究所采用的转化率梯度建模方式如图 9所示。
根据图 9构建木块微元最终转化率αs与炭化层深度x(mm)的函数关系，得到
在应用过程中需要考虑网格尺寸带来的影响，式(5)会导致质量缓慢下降阶段的初始木块质量偏小。本研究基于现有网格尺寸，对式(5)修正后得到
3 结论
本文基于有限的计算成本，通过木材燃烧实验与模拟，探讨了基于Fluent UDF的木块引燃模拟方法，并进行了模拟方法的可参考性检验，得到的主要结论如下：
1) 实验结果表明，可将木块的引燃过程分为3个阶段，即蒸发阶段、热解阶段和质量缓慢下降阶段；燃烧结束后木块表面的炭化层分布较为均匀，厚度在5 mm左右，其形成不仅受温度影响，还受到内部固体介质气密性的影响；炭化层内部存在转化率梯度，表现为表面区域高、内部区域低。
2) 提出了基于Fluent UDF的木材引燃模拟方法，将实验数据和拟合结果以代码的形式融入UDF中，以精确描述木材引燃过程。通过设置虚拟换热面、有效热解放热量以及炭化层内部的热解终止点分布，本研究得到的模拟结果与实验数据吻合良好。
