一种用于预测砖石炭生产窑自然冷却过程的CFD（计算流体动力学）模型

《Fuel》：A CFD model to predict the natural cooling of masonry charcoal production kilns

　　蒂亚戈·帕伦特·利马（Thiago Parente Lima）| 马西奥·阿雷德斯·马丁斯（Marcio Arêdes Martins）
巴西杰基廷honha和穆库里河谷联邦大学（Universidade Federal dos Vales do Jequitinhonha e Mucuri），科学技术研究所

**摘要**
在使用砖砌窑生产木炭的过程中，冷却阶段是主要的瓶颈，因为它可能导致窑炉卸载前需要数天的停机时间，并直接影响生产效率。本研究的目的是开发并验证一个能够预测砖砌木炭窑自然冷却过程的计算流体动力学（CFD）模型。所提出的模型考虑了传导、自然对流和辐射的耦合传热效应。木炭床被建模为一个由热解气体渗透的多孔介质。通过在能量方程中加入源项来描述木炭颗粒与气体之间的热传递。窑炉通过自然对流和辐射与环境进行热交换。为了避免显式模拟碳化过程，初始温度场是根据冷却开始时的实验测量数据重建的。对于二维几何形状的模型预测结果，通过与实际矩形砖砌窑进行的实验数据进行了验证，实验数据包括了分布在木炭床、气体区域和窑炉墙壁上的十个温度测量点。模型预测的所有温度的平均绝对百分比误差均低于6.1%。其中，墙体的温度预测表现最佳，平均绝对百分比误差低于3%，R2值达到0.98。敏感性分析表明，窑炉的冷却过程主要受墙体热容量和外部热传递阻力的影响，而气体-床层间的热传递参数的影响较小。所提出的CFD框架为分析砖砌窑的自然冷却过程提供了重要工具，并为设计和优化旨在减少冷却时间和提高木炭生产效率的技术奠定了基础。

**1. 引言**
木炭是一种由生物质制成的固体生物燃料，通常来自人工种植的森林。全球木炭生产集中在非洲（65%）、亚洲（17%）和南美洲（15%）[1]。在南美洲，特别是巴西，木炭的生产和使用具有特定特征。巴西是世界上最大的木炭生产国[1][2]，其91%的木炭来自人工林，几乎全部由桉树 plantation 组成[3]。在巴西，木炭主要用于工业用途：约89%的木炭被冶金行业用于钢铁、铁合金和生铁的生产[2]。这种能源的使用为减少钢铁行业对化石燃料的依赖提供了机会，并带来相应的环境效益。
尽管具有这些潜力，但大多数木炭生产仍依赖于技术较为落后的窑炉，主要是砖砌窑。这类窑炉通常由烧制粘土砖建造而成，配备有装卸门以及侧面的小开口，用于控制窑炉的点火和氧气供应[4]。由于成本低、无需电能且操作简便，砖砌窑在巴西和非洲仍被广泛使用[5]。

木炭的生产过程包括四个阶段：装载（将原木放入窑内）、碳化（通过开口点燃原木并使其经历热解过程）、冷却（窑炉完全封闭直到内部温度降至约333 K）以及最终卸载产出的木炭[4]。使用砖砌窑生产木炭时，最大的瓶颈在于窑炉冷却所需的时间较长[5]。砖砌墙体在窑内保留热量，显著延长了冷却阶段。冷却通常依靠与环境的自然热交换来实现，未采用专门增强热传递的技术。在尺寸为32米×4米×4米的矩形窑炉中，冷却时间约为七天[6]；然而，在不利环境条件或操作故障的情况下，同一窑炉的冷却时间可能延长至17天[6]。一些作者在实验研究中提出了通过人工冷却技术（如水蒸气注入[7]和强制气体循环的热交换器[6][8]）来缩短冷却时间的方法。尽管有专利提出在冷却阶段移除窑顶的隔热层[5]，但由于电力需求或高昂的实施成本，这些快速冷却方案未能得到广泛采用。

现有的少数实验研究[6][7][8]（均在地巴西进行，鉴于该国的木炭生产和使用特点）也反映了在砖砌窑冷却方面进行实验研究的困难。这些困难源于建造原型的高成本或全尺寸实验对工业生产造成的干扰。因此，优化这些系统以降低实施和运行成本成为一项挑战。
利用计算流体动力学（CFD）等技术的计算建模和仿真为优化木炭生产系统提供了有前景的手段。一些研究[9][10]使用CFD技术模拟了木炭窑内的气体流动，但未考虑热传递过程。[11]中评估了使用双管热交换器进行冷却的效果，但由于模型简化限制了其实际应用性。
第一个专门针对木炭生产砖砌窑冷却过程的计算研究发表在[12]中。该模型采用的主要简化方法是假设内部气体处于静止状态，因此只解决了导热方程，并根据不同窑区调整了物理参数。模型采用了视角因子来考虑窑内的辐射热传递。[12]中的简化大大降低了计算成本，并且在预测窑壁温度方面与实验结果吻合良好。

[13]中提出了一个更复杂的模型，该模型同时考虑了碳化阶段（通过加热床层实现）和冷却阶段。作者考虑了地面、墙体、气体及木炭床中的热传递，木炭床被建模为多孔介质。窑内的辐射热传递使用P1辐射模型进行模拟，模型结果与窑壁内表面的温度测量值吻合良好。

本研究旨在模拟砖砌窑在自然冷却过程中的热传递，包括传导、自然对流和辐射三种热传递机制。模型结果将与[12]中提供的实际矩形窑炉的墙体、气体和木炭床区域的温度数据进行比较。所提出的模型通过从预定义的温度场开始并避免模拟碳化阶段，从而降低了[13]中计算的复杂性。地面被近似为绝热体，从而减小了计算域，同时仍考虑了窑内的关键热传递过程（如气体流动）。预计该模型将有助于促进CFD技术在改进砖砌窑冷却过程中的应用，最终提高生产效率。

**2. 方法论**
本节介绍了用于验证数值模型结果的实验方法、提出的数学模型及用于求解该模型的数值程序。

**2.1. 实验数据**
本研究使用的实验数据来自[12]，该数据是在巴西米纳斯吉拉斯州帕拉奥佩巴市（Paraopeba）农村地区的一个全尺寸砖砌窑的自然冷却阶段获得的。实验过程在葡萄牙语中有详细描述[12]，下面简要概括以便读者理解。
该窑炉的主要外部尺寸为16.4米×3.88米×5.19米（长×高×宽），内部体积为190立方米，其构造布局如图1所示。这种带有穹顶形状屋顶的矩形窑炉通常被大型冶金公司使用，因为它允许机械化装卸木材和木炭。窑炉的墙体和屋顶由烧制粘土砖建造，并用粘土和焦油制成的砂浆粘合。墙体和屋顶的外表面也涂有这种砂浆。窑底由当地土壤构成，通常为粘土质，经过压实且未作任何内衬处理。

整个实验包括装载、碳化和冷却三个阶段，总持续时间约二十天，于1999年6月至7月期间进行。冷却阶段持续约九天，直到木炭床中心区域的温度降至333 K。仅在冷却过程的前六天内记录了窑温；这一时期结束时，热电偶显示的温度接近373 K。同期监测到的环境温度范围为284 K至313 K，但[12]未记录这些数据。实验期间未出现降雨事件，否则可能影响冷却过程。

窑内的温度监测使用了十个热电偶：三个安装在木炭床区域，五个位于床层上方的气体区域，另外两个嵌入墙体和穹顶中，如图1所示。实验中使用了K型热电偶，因为它们具有宽广的温度范围和抗氧化性能，适用于工业炉窑。热电偶T2、T3、T4、T5、T6和T8的数据由数据采集系统记录，而T9、T10、T11和T12则是在特定时间手动读取的。
窑炉从底部到顶部装载了54320公斤的Eucalyptus grandis木材，平均年龄为六年。原木的平均长度为1.8米，直径范围主要在0.10米至0.12米之间。碳化过程结束后，木炭床占据了木材初始体积的54%，据此估算出木炭块的平均直径为0.063米。冷却完成后，木炭床的高度为1.4米，其孔隙率为34.3%。

**2.2. 数学模型**
研究的问题域定义为包括窑体、墙体和内部空间（包含木炭床及热解产生的气体）。窑炉的尺寸为16.4米×3.88米×5.19米（长×高×宽）。因此假设沿纵向方向的变化可以忽略不计，仿真基于窑中央区域的横截面进行，该横截面对应于图1中热电偶的位置。这种简化降低了仿真计算成本。

预计沿纵向方向的变形在整个窑长范围内可以忽略，仅在靠近前后门的部位存在显著变化。因此，在三维仿真中，预测的冷却速率可能会高于二维模型的结果。不过，近年来矩形窑炉的长度显著增加，当前研究中的窑炉长度约为现代矩形窑炉的一半，从而减少了前后门对总体冷却时间的相对影响。

除了二维域外，模型还假设环境温度恒定，以及对流和辐射热传递系数恒定。窑与地面之间的热交换被忽略了，窑内的物种传输和气相化学动力学也被忽略了。在这些假设下，研究域被划分为三个子域：(i) 窑的结构，包括墙壁和穹顶；(ii) 内部气体区域，从地面延伸到穹顶，充满热解气体；以及 (iii) 木炭层区域，包含固体木炭。如图1所示，后者部分与气体区域重叠。气体区域和木炭层之间的重叠将在后面详细描述。热解气体被建模为一种固定组成的理想气体混合物，能够吸收和发射辐射。考虑的质量分数分别为0.211%的CH4、1.05%的CO、4.33%的CO2、8.00%的H2O和86.41%的其他气体，假设这些其他气体是N2。这些组成数据来自作者在尺寸与本研究分析的类似的砖砌木炭窑中进行的测量。

简而言之，当窑被密封时，使用气体采样泵（N86KTE，KNF，德国）和流量计（G1.6，Gasometer，巴西）在颗粒间收集气体样本，这标志着冷却阶段的开始。收集的样本立即在-4°C下冷却以分离凝结物组分。气流在3升容量的洗涤器中用去离子水清洗，在硅胶柱中干燥，并通过0.1微米膜过滤以去除任何剩余颗粒。清洁后的干燥气体样本的组成通过红外气体分析仪（Gasboard 3100，武汉立方光电科技有限公司）进行量化。然后，根据收集的凝结物质量和流量计测量的气体量，将气体组成表示为质量分数（湿基）。

2.2.1 气体区域的控制方程
窑内的气体压力和速度是使用雷诺平均纳维-斯托克斯（RANS）动量守恒方程计算的，如方程(1)所示[14], [15]：
(1) ?ρu/?t + ?·ρuu = ?·τ + σt - ?p + ρg
其中ρ是流体密度，u是平均流体速度向量，t是时间，τ是平均应力张量，σt是雷诺应力张量，p是流体平均压力，g是重力加速度向量。方程(1)中的雷诺应力张量使用k-ω SST湍流模型封闭，该模型在自然对流问题中表现出良好的性能[16]。k-ω SST模型的控制方程如下[17]：
(2) ?ρk/?t + ?·ρuk = P + β?ρωk + ?·μ + σkμt
(3) ?ρω/?t + ?·ρuω = γνtP + βρω2 + ?·μ + σωμt
(4) ?ρh/?t + ?·ρuh + ?ρK/?t + ?·ρuK - ?p/?t = -?·q + ρr - ρg·u
其中k是湍流动能，P是湍流产生项，ω是特定耗散率（湍流频率），μ是动态粘度，μt是湍流涡流粘度，νt是湍流动力学粘度，F1是湍流模型函数，β*、σk、γ、β、σω、σω2是湍流模型常数。气相中的热传递使用熵形式的能量守恒方程进行建模，如方程(4)所示，其中忽略了机械能向热能的不可逆转换[18]：
(4) ?ρh/?t + ?·ρuh + ?ρK/?t + ?·ρuK - ?p/?t = -?·q + ρr

2.2.2 木炭层的控制方程
木炭颗粒的平均直径（dp）远小于层高度（Lb），即dp?Lb 或 0.063 m?1.4 m。因此，木炭层区域被建模为由两相组成的填充层：气体和固体[20]。这两相被视为相互渗透的连续介质，它们的控制方程在单独的网格上求解（图3），并通过代表气体和固体相之间对流热传递的源项（方程(7)）耦合起来。木炭层中的固体相仅通过这个源项与气体相交换热量。忽略了木炭颗粒与窑壁之间的热传递。通过气体和固体相的有效热容量将层孔隙率纳入模型。具体来说，层孔隙率?及其补数（1 - ?）被包含在计算层区域内气体和固体相的比热中。孔隙率也被考虑在计算对流热传递系数和气体与固体相之间的界面热交换面积时。

在固体相中，能量方程根据内能（方程(6)）求解[15]，包括由方程(7)给出的对流热传递源项[15]：
(6) ?ρe/?t = -?·q + ρr
其中e是内能。

添加到木炭层区域内气体（方程(4)）和固体（方程(6)）相的能量方程中的热交换源项是使用方程(7)计算的[15]：
(7) -Q = htcAVTs - Tf
其中Q是体积热传递率，htc是气体与木炭之间的对流热传递系数，AV是单位体积的热交换面积，Ts是固体相温度，Tf是气体相温度。方程(7)需要指定单位体积的热交换面积（AV）和气体与木炭之间的对流热传递系数（htc）。测量的层孔隙率为34.3%，接近于随机填充的刚性球体的36.5%的值[22]。因此，木炭层被近似为平均直径等于木炭颗粒估计平均直径的球体填充层。基于这个假设，可以使用方程(8)计算单位层体积的热交换面积[23]：
(8) AV = ApVp / (1 - ?)
其中Ap是颗粒表面积，Vp是颗粒体积，?是填充层孔隙率。

方程(7)中使用的对流热传递系数是使用Whitaker关联式计算的[23]，该关联式适用于随机填充的球体和圆柱体。平均努塞尔特数（Nu）的定义在方程(9)和(10)中给出[23]：
(9) Nu = 0.5Re^(1/2) + 0.2Re^(2/3)Pr^(1/3)
(10) Nu = htc / (dkpκ?^(1 - ?))
其中Re = (ρUdp) / (μ(1 - ?))是雷诺数，U是表面（达西）速度，Pr = 0.7是普朗特数，κ是热导率。

在层区域内，当气体流经多孔介质时，它们会受到与木炭颗粒相互作用的粘性和惯性力的影响。这些效应通过在整个层区域内添加动量源项到方程(1)中来考虑。这个源项是使用Darcy-Forchheimer模型计算的，如方程(11)所示[15]：
(11) ?p = -μd + ρu^2f
其中d是达西系数，f是Forchheimer系数。方程(11)中的达西和Forchheimer系数是使用Dos Santos等人提出的实验关联式计算的[21]，用于通过木炭层的流动压力梯度[12]：
(12) -?p = 21.58 + 343.65uu

2.2.3 窑墙的控制方程
窑墙被建模为连续固体，其温度场使用方程(6)计算。窑墙通过传导和对流在内部表面与窑内的气体交换热量，通过传导和对流在与外部的表面与环境交换热量，具体边界条件如下所述。

2.2.4 边界条件
作为边界条件，假设外部窑表面通过自然对流和辐射与环境交换热量。采用了两个平均对流热传递系数来考虑窑与环境空气之间的对流。对于垂直墙，使用方程(13)定义的平均努塞尔特数[24]，并通过方程(14)给出的关联式计算[24]，其中墙高度被视为特征长度：
(13) Nu = hcLκ
(14) Nu^(1/2) = 0.825 + 0.387Ra^(1/6) + 0.492/(Pr^(9/16)8/27)
其中hc是对流热传递系数，L是特征长度，Ra = [gβf(Tw-Ta)L3]/(να)是瑞利数，βf是体积热膨胀系数，Tw是墙体表面温度，Ta是环境温度，ν是层流运动粘度，α是热扩散率。

对于窑穹顶，使用方程(15)估计努塞尔特数[25]，特征长度定义为L = As/Pd，其中As是窑穹顶表面积，Pd是窑穹顶周长：
(15) Nu = 0.15Ra^(1/3)
假设墙和穹顶的对流热传递系数在整个模拟过程中保持不变。

外部窑表面与环境之间的辐射热交换使用方程(16)中的线性化辐射系数计算[15]：
(16) hr = εsσTa^2 + Tw^2/Ta
其中hr是辐射热传递系数，εs是表面发射率。对流和辐射系数结合成一个有效热传递系数，如方程(17)所示，从而在整个窑外部表面再现了第三类（Robin型）边界条件[15]：
(17) heff = 1/hc + 1/hr^(-1)
其中heff是有效对流热传递系数。在内部流体-固体界面，应用了耦合边界条件以确保能量方程（方程(4)）中的温度和能量通量的连续性。在动量方程（方程(1)）中，对速度施加无滑移条件，对压力采用零梯度条件。对于湍流模型方程，使用了对y+不敏感的墙边界条件，适用于低和高质量的雷诺数流动[15]。

2.2.5 热物理性质
表1总结了模型中使用的一组热物理性质和输入参数。在方程(13)、(14)、(15)中使用的物理性质是在垂直墙和穹顶的膜温度407 K和448 K下评估的。假设的环境温度为298 K，对应于实验期间观察到的最小和最大值的平均值[12]。

表1. 模型中使用的热物理性质和输入参数
| 参数/性质 | 值/注释 | 参考文献 |
|--------|---------|---------|
| 气体 | μ = AS ± TS/TH | [26] |
| | AS = 2.978 × 10^-7 kg m^-1 s^-1 K^-1/2, TS = 1064 K |
| CH4 | AS = 1.0194 × 10^-6 kg m^-1 s^-1 K^-1/2, TS = 171.06 K |
| CO2 | AS = 1.503 × 10^-6 kg m^-1 s^-1 K^-1/2, TS = 222 K |
| CO | AS = 1.512 × 10^-6 kg m^-1 s^-1 K^-1/2, TS = 120 K |
| N2 | AS = 1.401 × 10^-6 kg m^-1 s^-1 K^-1/2, TS = 107 K |
| cp/R | = a4T + a3T + a2T + a1T + a0 | [27] |
| 木炭层 | κ = μcv^1.32 + 1.77R/cv | [15] |
| | ρ = 1/RT | [15] |
| | a | [19] |
| | cp = 1017 J kg^-1K^-1 | [21] |
| | κ = 0.03 W m^-1K^-1 | [21] |
| | ρ = 344.6 kg m^-3 | [21] |
| | φ = 0.34 | [12] |
| | AV = 62.6 m^2 m^-3 | [23] |
| | Pr = 0.70 | [12] |
| | dp = 0.063 m | [12] |
| | ?p = 方程(12) | [21] |
| 墙壁 | cv = 920 J kg^-1K^-1 | [28] |
| | κ = 2.0321 × 10^-1 + 1.9473 × 10^-3T^-5.8656 × 10^-7T^2 | [29] |
| | ρ = 2000 kg m^-3 | [28] |
| | ε = 0.95 | [28] |
| (垂直墙) | hc = 6.56 | [24] |
| (穹顶墙) | hc = 9.10 | [25] |

2.2.6 初始温度场
在冷却阶段开始时，窑具有在碳化阶段发展出的温度场。因此，为了模拟自然冷却过程而无需显式建模碳化阶段，必须估计一个初始温度场。这个初始温度场是根据冷却阶段开始时的实验数据得出的。数学模型被拟合到这些数据上，以重现已计算域内的初始温度分布。最初假设木炭层上方区域的气体温度场是均匀的，温度为964 K，对应于位于该区域的温度计T3、T6、T8、T9和T10的读数的平均值（图1）。在模拟过程中观察到，在该区域采用不同的初始温度剖面不会显著影响结果，因为由于辐射模型考虑的吸收和发射效应，温度场会迅速均匀化。在木炭层区域，气体和固体相采用了相同的初始温度场。该场是通过将多项式模型拟合到位于层内的温度计T2、T4和T5测量的温度数据并受到位于其边界的温度计T6、T3和T11测量的温度限制而构建的。因此，床层区域中气体和固体相的初始温度场由方程（18）定义。(18)T=849.06-0.057395x3+209.70y3+34.048x2y-0.19853xy2-52.046x2-181.66y2-152.00xy+233.90x-178.00y，其中x和y分别是x方向和y方向的坐标。窑墙的初始温度场是通过初步模拟确定的。由于墙壁的高热容量，假设在冷却阶段开始时，墙区域呈现出近似线性的准稳态温度分布，这是由于窑内表面温度几乎恒定以及外部表面的对流和辐射热交换造成的。这一假设与[13]中报告的窑冷却早期阶段的结果一致。因此，为了定义初始墙温度场，在稳态条件下模拟了墙区域的热传递，考虑了热电偶T2、T3和T10在这些传感器附近的内表面测量的恒定温度边界条件，以及使用方程（13）、（14）、（15）、（16）、（17）计算出的与环境的相同对流和辐射热交换条件。基于这些假设，墙的初始温度场如图2所示。下载：下载高分辨率图像（70KB）下载：下载全尺寸图像图2. 初始温度场。初始速度场和入射辐射场被设置为零。压力场是根据相应的静水压力分布初始化的。湍流模型所需的其余场被初始化为均匀的，其值分别为k = 3 × 10-4 m2 s?2，ω = 0.1 s?1，νt = 0 m2 s?1，以及湍流热扩散率αt = 0 m2s?1。下载：下载高分辨率图像（194KB）下载：下载全尺寸图像图3. 计算域的离散化：气体和墙子域（左）；以及竹炭填充床的多孔子域（右）。2.3. 数值程序数学模型使用OpenFOAM 13软件包[15]实现和解决。OpenFOAM是一个开源的计算流体动力学（CFD）框架，由大量的用户和开发者社区开发和维护了二十多年[30]。它应用于耦合热传递问题和多孔介质流动在文献中得到了很好的验证[31]、[32]、[33]。计算域使用hexahedral和polyhedral单元的组合进行离散化，这些单元是用blockMesh和snappyHexMesh工具生成的。最终离散化的域（即计算域）由两个不同的网格组成（图3）。一个网格用于离散化气体和墙区域，分别求解方程（1）、（2）、（3）、（4）、（5）及其相应的源项，以及方程（6）。第二个网格用于离散化竹炭床的固体相，其中求解包括方程（7）中定义的源项的方程（6）。使用三种不同的网格分辨率N1、N2和N3进行了网格独立性研究。对于墙和气体区域，网格分别包含15386、8250和4579个单元，而对于竹炭床的固体相，相应的网格分别包含7003、3069和1386个单元。在所有情况下，如[34]中定义的网格细化比率r均保持在1.3以上。图1中显示的十个热电偶位置测量的温度被用作网格独立性研究中的感兴趣量。在最终模拟时间130小时时，网格N1和N2之间的最大近似相对误差[34]出现在热电偶T2处，为0.6160%。基于这一结果，认为网格N1对于本工作中的模拟来说是足够细化的。最终的计算网格（图3）总共包含21,973个单元，分布如下：墙区域7419个单元，气体区域7967个单元，竹炭床的固体相7003个单元。控制方程使用二阶精确方案进行空间离散化，包括对于散度项的二阶迎风方案和对于梯度、拉普拉斯和其他空间算子的线性方案。瞬态项使用Euler方案[15]离散化。结果方程组使用PISO算法求解压力-速度耦合。模拟时间步长被调整以保持最大Courant数为4[35]。3. 结果和讨论图4展示了冷却阶段热电偶T2、T4和T5的温度随时间的变化，这些热电偶位于竹炭床区域内。在数值模型中，这些热电偶与竹炭床内的气体温度相关联，因为在实验中，热电偶测量的是通过床层渗透的气体温度，而不是固体竹炭颗粒的温度。下载：下载高分辨率图像（248KB）下载：下载全尺寸图像图4. 热电偶（a）T2、（b）T4、（c）T5的温度，以及（d）它们在窑中的位置。模型预测的温度与热电偶T2和T4的测量结果吻合良好，这两者相对于窑底的高度大致相同。相反，模型低估了热电偶T5测量的温度，T5的位置略高于T2和T4，位于床层内部。热电偶T5的实验数据表明，在冷却阶段开始时温度略有上升，这是由于热量从床层的较热核心向周边区域传递造成的。模型能够定性捕捉到这种行为；然而，它低估了该区域的热惯性，导致温度上升和下降的速度比实验观测到的要快。确定竹炭床的有效热性能本质上具有挑战性。空间变化可能来源于构成床层的材料本身的差异，如竹炭颗粒的热性能和几何特性，以及它们的随机堆积方式，这可能会局部增强或阻碍气体流动。这种不确定性是建模此类系统的主要挑战之一，因为这些变化通常是事先不知道的。使用图5中的身份图比较了位于床层内的热电偶的数值和实验结果。还显示了一条虚线带，对应于与身份线相差±5%的偏差。在床层区域，热电偶T2的模型性能最好，决定系数R2 = 0.9920，平均绝对误差（MAE）为9 K。由于需要估计初始温度场，因此在冷却过程开始时模型缺乏碳化阶段的考虑导致较大的误差，如图5所示。随着冷却的进行，这些误差趋于减小。下载：下载高分辨率图像（276KB）下载：下载全尺寸图像图5. （a）T2、（b）T4、（c）T5和（d）它们在窑中的位置的模拟和实验温度比较。图6展示了冷却阶段热电偶T3、T6、T8、T9和T10的温度随时间的变化，这些热电偶位于竹炭床上方的气体区域。模型预测与热电偶T9和T10的测量结果吻合良好。相比之下，模型低估了热电偶T3、T6和T8测量的温度。当分析气体区域的实验数据时，热电偶T9和T10一致报告的温度低于同一区域的其他热电偶。应该注意的是，热电偶T9和T10的测量是手动获得的，而热电偶T3、T6和T8是使用自动数据采集记录的。除了这种测量方法的差异外，[12]中没有报告热电偶之间的其他区别特征。下载：下载高分辨率图像（340KB）下载：下载全尺寸图像图6. （a）T3、（b）T6、（c）T8、（d）T9、（e）T10的温度，以及（d）它们在窑中的位置。图6显示，在冷却阶段的最初几个小时内气体温度急剧下降。模型捕捉到了这种行为，尽管它略微高估了初始气体冷却速率。应强调在辐射模型中考虑气体作为参与介质的重要性。在初步模拟中，假设气体是透明的，气体温度在过程开始几小时后出现了一个人为峰值。在这种条件下，这个峰值是由离开床层并到达窑上部的热气体上升引起的。当在辐射模型中包含参与介质的假设时，这种行为消失了，冷却曲线显示了实验中观察到的趋势，如图6所示。这一结果表明，当考虑气体流动时，基于视图因子的辐射模型不适合模拟竹炭窑。图7展示了模型预测的床层上方体积平均气体温度随时间的变化，以及同一区域热电偶的实验测量结果。热电偶T9和T10一致报告的温度低于热电偶T3、T6和T8，尽管这种差异在冷却过程中逐渐减小。热电偶T3和T6位于床层更接近的位置，因此更多地暴露在热竹炭床的辐射下，这可能解释了它们记录的温度较高。尽管热电偶T8与热电偶T9位于相似的高度，但它位于窑的中心附近，模型预测那里有一个上升的热气体羽流（图8）。这种流动结构可能解释了热电偶T8记录的温度比T9高的原因，后者与床层的距离相似。实验数据与模拟的体积平均温度之间的比较得出决定系数R2 = 0.9732，表明这个指标合理地代表了该区域气体的平均热行为。下载：下载高分辨率图像（92KB）下载：下载全尺寸图像图7. 竹炭床上方气体的温度。下载：下载高分辨率图像（230KB）下载：下载全尺寸图像图8. 10小时冷却后的速度矢量场和背景温度场。图9(a)显示了安装在窑的垂直壁和圆顶上的热电偶测量的模拟和实验温度变化。模型预测的温度与两个热电偶的实验数据吻合良好。只有两个计算的温度点落在±5%相对误差带之外，如图9(b)中的身份图所示。初步模拟中采用的边界条件足以提供热电偶T11的初始温度的良好估计。尽管在初始时间没有为热电偶T12达到相同的准确度，模型仍然与实验数据整体上吻合良好。下载：下载高分辨率图像（385KB）下载：下载全尺寸图像图9. （a, c）T11、（b, d）T12和（e）它们在窑中的位置的模拟和实验温度。窑的圆顶显示的温度显著高于垂直壁。这些升高的温度是由于热气体在窑中央区域的持续上升形成的热羽流造成的，如图8所示。这种流动模式将床层侧面的气体带入中心，促进它们上升到圆顶。到达窑的上部后，气体被重新导向侧面并沿墙壁下降，冷却直到重新进入床层，这导致垂直壁的温度较低。壁温预测的良好一致性也在[12]和[13]中报告，这一点尤为重要，因为这些热电偶由于受到的辐射效应较少，因此通常提供更可靠的测量结果。表2总结了所有分析的热电偶的R2、MAE和MAPE值。总体而言，模型准确预测了竹炭床内的气体温度，特别是对于热电偶T2和T4，而热电偶T5的预测也可以认为是令人满意的。床层内的气体温度是一个重要的操作参数，因为它指示了开启窑的适当时间，通常与333 K的温度相关。对于床层上方的气体区域，模型对于距离床层较远的熱电偶T9和T10的表现更好。对于窑壁内的温度预测，T11和T12获得了最佳的整体结果，结合了高R2值和低MAE和MAPE。根据完整的实验结果，可以得出结论：所提出的模型能够 satisfactoryly reproduce the thermal behavior of natural cooling in masonry kilns used for charcoal production（即能够令人满意地再现用于木炭生产的砖砌窑的自然冷却过程），并可以作为通过计算模拟开发 and 评估 new kiln cooling technologies（新的窑冷却技术）的基础。表2显示了所有热电偶的R2值、MAE值和MAPE值。

模型的冷却时间受到表1中列出的多个输入参数的影响。因此，有必要评估模型对这些参数不确定性的敏感度。为此，我们进行了敏感性分析，考虑了五个参数：气体与木炭层之间的对流热传递系数（公式（7））、Darcy–Forchheimer系数（公式（11）、窑壁和穹顶处的对流热传递系数（公式（14）、（15））、外部窑表面的发射率（公式（16）以及窑壁的热容量。前两个参数直接影响气体与木炭层之间的热传递阻力，而后三个参数则影响窑与周围环境之间的热传递阻力。敏感性分析是通过在预定义的范围内单独改变这些参数并评估它们对模拟窑温的影响来进行的。

所考虑的五个参数具有不同的不确定性水平。例如，应用于窑穹顶和木炭层的努塞尔特数的不确定性分别在大约±7% [25] 和 ±25% [23] 之间，即使是在受控的实验室条件下也是如此。其他参数，如窑壁的热性质，也表现出显著的变化性。巴西关于建筑热性能的标准推荐固体烧制黏土砖的发射率值在0.85到0.95之间 [28]，相当于平均值约±5%的变化范围。用于窑建造的砖的密度在文献中显示出更广泛的变化范围：[36] 对来自29个不同制造商的砖进行的测量报告了1610到2120 kg/m3之间的密度，相当于平均值约±13%的变化；而同一作者汇编的文献数据则指出密度在1200到2200 kg/m3之间。还有研究报道了更高的密度，达到2400到2550 kg/m3 [37]。文献中报告的砖的比热值变化相对较小，大约在±3% [28]、[29] 的范围内。因此，通过改变砖的密度来评估模型对墙热容量的敏感性。

鉴于这些不确定性，在敏感性研究中，参数的变化范围为?25%到+25%，除了外部窑表面的发射率，因为这种变化会导致物理上不现实的值；对于Darcy–Forchheimer系数，即使变化±25%，模拟也无法达到数值收敛。

图10显示了使用公式（19）计算得出的热电偶T11在冷却阶段开始48小时后的温度变化。选择热电偶T11是因为它对分析的参数最为敏感。

图10中观察到的趋势对所有热电偶都是相似的。模型对与气体和木炭层之间的热传递相关的参数敏感性较低。窑的冷却过程主要由窑与环境之间的热传递阻力控制，这包括自然对流热传递系数、外部窑表面的发射率以及窑材料的热容量。如图10所示，发射率的影响与自然对流系数的影响相当。然而，应该注意的是，与估计砖的发射率相关的不确定性通常大于与估计对流热传递系数相关的不确定性。

决定砖砌窑冷却的主要因素及其建模的是窑壁的热容量。正如预期的那样，更高的墙热容量会导致更大的热惯性，从而延长冷却时间。在本研究中，采用了2000 kg/m3的砖密度，这个值在[36]报告的范围内，并且符合巴西建筑标准[28]的推荐值。在研究过程中发现，1800到2000 kg/m3的砖密度足以使表2中报告的所有决定系数R2保持高于0.95。

图11展示了冷却阶段开始12小时后热电偶T2的模拟温度变化。该图显示，气体-床层对流热传递系数和Darcy–Forchheimer系数对木炭层内气体温度的影响更为显著。这种增加的敏感性是由于冷却初期存在更大的温度梯度造成的。随着冷却的进行和温度梯度的减弱，这两个参数的影响也减小，达到与图10中观察到的相似程度。

图10和图11表明，降低窑壁的热容量可以加速冷却过程。在[5]中提交的专利调查中报告了几种这方面的策略，例如在金属容器中生产木炭，这些容器在碳化阶段后可以从砖砌窑中移除；或者构建在穹顶具有可拆卸隔热层的窑。另一种方法是通过在冷却阶段将热解气体引导到外部热交换器来绕过窑壁的热阻。这种策略在[16]、[19]中被证明是有前景的。然而，在包含此类修改的窑中进行实验的成本和技术挑战阻碍了其优化。

提出了一个CFD模型来模拟用于木炭生产的砖砌窑的自然冷却阶段。模型预测的木炭层、气体区域和窑壁的温度与全尺寸实验数据进行了验证，显示出良好的一致性。对于位于窑壁内的热电偶，最佳预测结果的R2值高于0.98，MAE值低于16 K。

初始温度场是根据实验测量结果估计的，这消除了明确模拟碳化阶段的需要。同样，忽略了窑与地面之间的热传递，从而可以在不显著影响模型性能的情况下减少网格尺寸和计算成本。

模型对窑壁的物理性质以及窑与周围环境之间的对流热传递系数更为敏感，其中墙热容量是最具影响力的参数。因此，在估计这些性质时应特别小心，尽可能使用实验测量的值。

敏感性分析表明，有效缩短窑冷却时间需要降低墙的热容量或减少外部窑表面与环境之间的热传递阻力。然而，降低这些参数也会降低窑在碳化阶段储存热能的能力，这是不希望看到的。因此，简单的改变，如更换建筑材料或减少墙厚，并不是吸引人的解决方案。更有前景的策略涉及在碳化阶段结束和冷却阶段开始之间对窑的配置进行修改，例如去除墙的隔热层或使用热交换器。

本文提出的CFD模型可以用来研究这类配置更改，或者与热交换器模型结合使用，以优化窑的设计和操作。通过这种方式，它的应用有助于提高用于木炭生产的砖砌窑的生产效率。

Thiago Parente Lima 的作者贡献声明：
- 编写 – 审稿与编辑
- 编写 – 原始草稿
- 验证
- 方法论
- 调查
- 形式分析
- 构想

Marcio Arêdes Martins 的作者贡献声明：
- 编写 – 审稿与编辑
- 编写 – 原始草稿
- 监督
- 方法论
- 调查
- 构想