天然气炉窑干燥过程多孔介质传热传质耦合数值分析

同济大学机械与能源工程学院 周 宇 秦朝葵

天然气炉窑干燥过程多孔介质传热传质耦合数值分析

同济大学机械与能源工程学院 周 宇 秦朝葵

砖坯的干燥是其塑性成型的逆过程，干燥过程涉及相变及多孔介质传热传质的耦合问题，文章从非稳态传热传质和相平衡的角度出发，基于多孔介质传热传质和计算流体动力学基础理论，建立了砖坯干燥过程的非线性三维数学模型，利用流体计算软件FLUENT对模型进行求解，得到砖坯干燥过程中内部含湿量和温度的变化过程，讨论了干燥介质入口速度和温度对干燥过程的影响，结合所得的干燥数据以及临界干燥温度梯度，利用Matlab对干燥介质入口速度、温度与干燥时间的关系进行拟合，得出在小于本文所设定的临界梯度的条件下，较为经济的干燥条件是温度为550～573 K，风速为0.5～2.5 m/s。

干燥过程 多孔介质 传热传质耦合 数学模型

我国烧结砖瓦行业中，无论是新建设的砖厂还是改造扩建的砖厂，隧道窑都逐渐成为了首选的砖瓦产品干燥以及焙烧设备。作为隧道窑所附属的坯料干燥设施，干燥窑在整个生产流程中发挥着重要的作用。

干燥过程中，隧道窑干燥室内的坯料与干燥介质进行强烈的热湿交换，窑内的干燥介质沿着炉窑长度方向具有不同的温度、速度和湿度。在不同的干燥阶段，热湿交换的强度和速度各不相同，随着干燥过程的进行，坯料和干燥介质的各项参数都发生了较大的变化。干燥过程不仅消耗大量能量，而且受坯料本身物性以及外界条件的影响，干燥效果和热效率都会发生变化，由于影响干燥过程的因素众多，干燥中出现的问题也比较复杂，如干燥倒坯、坯料出现的裂纹等，然而目前对于耐材生产热力过程的研究主要集中于烧成阶段，对于干燥阶段的研究则相对较少。

砖坯的干燥过程不仅涉及相变过程，还需要考虑多孔介质传热传质的相互耦合，按照控制多孔介质传热传质过程中动量、能量和质量传递的不同模式，多孔介质的干燥过程一般可视为多孔介质中导热与质量扩散过程。多孔介质传热传质研究发展到今天己经取得了很多的进步，研究者们在过去的几十年内建立了数以百计的模型以描述多孔介质的干燥过程，并且根据湿分在在固体物质内的初始状态、驱动力和运动形式等因素相继提出液态扩散、蒸发冷凝理论、毛细理论、Luikov理论、Whitaker体积平均理论等，由于多孔介质本身结构的复杂性、多变性和随机性，使得对其传热传质过程的研究手段仍然以实验研究为主，但是实验研究相对昂贵，费时较长，相对于实验来说，数值计算的投资要小得多，具有较好的通用性、灵活性、快速性及经济性。

本文结合炉窑实测数据，建立砖坯干燥数值模型，利用Fluent软件对砖坯内部水分迁移进行数值计算，重点分析单块砖坯在不同干燥条件下的干燥情况，分析干燥介质入口速度、温度对干燥效果的影响。

1 模型及边界条件

1.1 计算模型

利用ICEM绘制单块砖坯的结构网格，砖坯模型的尺寸为150×150×300 mm，如图1所示。

图1 计算模型网格示意

干燥介质进出口为单进单出，设置左边面为进风口，右边面为出风口，墙壁为绝热面。为了保证计算精度，模型全部采用六面体结构化网格，网格数量为49 062。

1.2 边界条件

流场的数值模拟需在有限区域内进行，边界条件一般是在求解区域的边界上，求解的变量随地点和时间的变化情况。对于Fluent计算，边界条件的设置直接影响到计算结果的精度。数值模拟中，初始条件和边界条件做如下设定和处理：

砖坯的初温设置为300 K，进入干燥窑之前的含水率为6%，出窑含水率要求低于1%。湿坯的密度为2 270 kg/m3，导热系数为1.3 W/(m·K)，比热容为910 J/(kg·K)，为了简化模型以上参数在模拟干燥过程中视为常量。

砖坯设置为多孔介质，设置其多孔率为常数0.4。需要注意的是，被选为多孔介质的区域必须是Fluid类型。

在直角坐标系中定义阻力系数的办法是在二维问题中定义一个方向矢量，或在三维问题中定义两个方向矢量，然后再在每个方向上定义粘性和惯性阻力系数。本次计算中粘性阻力系数为4e+10，惯性阻力系数为4e+09。

计算过程中涉及水蒸气和液态水共存现象，因此采用两相流模型，本文对两相流的液相、气相及混合相均采用相应的udf进行定义。

除两相流外，本次计算开启能量方程，流动采用k-epsilon模型。

2 结果及分析

为了分析不同干燥条件下砖坯内温度梯度和水分梯度的变化及其对干燥效果的影响，对不同送风速度和送风温度下的干燥过程进行了数值模拟，本节选取风速为5 m/s、2 m/s，温度为573 K、600 K下的干燥过程做详细分析。

定义水分饱和度S为坯料中含湿量占总湿分的比例，初始条件水分饱和度为1，完全干燥后为0，取干燥介质入口速度垂直方向上的坯料的中心面(以下简称中心面)，观察该中心面的平均温度以及水分饱和度的变化，以得出坯料中心面的干燥曲线。

为了便于观察砖坯内部温度、水分梯度的分布，对模型做如下处理：沿着Z方向和X方向各截取8个面，如图2所示。

图2 温度梯度与水分梯度截面示意

沿Z方向截取的面与干燥介质入口风速垂直，其坐标为-0.14～-0.1 m；沿X方向截取的面与干燥介质入口风速平行，其坐标为0.1～0.14 m。温度梯度和水分饱和度梯度的算法如下：

式中：T1、T2——相邻面的面平均温度，K；

S1、S2——相邻面的面平均水分饱和度；

△L——相邻面的面的距离，m。

2.1 入口风速5 m/s，温度573 K时干燥过程

图3所示的是干燥介质入口风速为5 m/s，温度为573K时砖坯中心面温度以及水分饱和度随着干燥时间的变化曲线。该干燥曲线与文献中提及的恒定干燥条件下典型干燥速率曲线基本吻合，间接印证了本模型的正确性。

图3 入口5 m/s，573 K时坯料中心面温度以及水分饱和度曲线

从图中可以看出，在100 s之前属于升温阶段，砖坯温度升高但水分饱和度几乎不变。在100～200 s之间，温度迅速升高，水分也随之快速蒸发，该阶段属于加速干燥阶段。在干燥时间t＜200 s时，砖的表面饱和度较高，此时液态水分的毛细流起主要作用，随着干燥的进行砖坯的表面脱水使得其表面饱和度比其内部饱和度降低要快得多。200～500 s之间干燥平稳，属于等速干燥阶段，在干燥进行约500 s时，砖坯中心面的平均水分饱和度已降至S=0.4。

图4和图5所示分别是垂直风速方向和平行风速方向截面的温度以及水分饱和度在不同干燥时间的梯度。从图4可以看出，砖坯中心的温度梯度基本为0，离外表面越近，其温度梯度越大。随着干燥时间的增加，砖坯温度和水分饱和度梯度呈逐渐减小的趋势。从图5可以看出，初始阶段，由于边界处存在对流，坯料两端的温度梯度高于中间部分，随着时间的递增，内部热传导的进行以及湿度的影响，坯料的内部温差变小，温度趋于一致，并随时间逐渐上升。由于初始时刻坯料表面的湿势远高于环境，所以湿分是以水汽形式由表面向外界传递，使得边界处的湿势变化很快。同样，由于质扩散系数较小，中间部分湿度变化相对比较缓慢，需要较长时间能达到平衡。

图4 垂直风速方向温度梯度及水分梯度示意

在砖坯的边角处，由于同时存在两个方向上的剧烈传热传质，其温度上升和湿度下降的速度都要较其他地方快，因此边界处也是最易造成干燥不均匀并产生裂缝的地方。砖坯的不同侧面之间温度有一定的差距，但差别并不大。干燥的开始阶段侧面干燥得较快，因而温度首先上升。而正面温度则上升较晚，但正面在表面干燥后温度上升得较快，并在干燥的最后阶段超过侧面。传质也有类似情况，干燥开始阶段侧面响应快，但是当表面干了之后，速率就有明显的下降。与之相反，正面虽然开始干燥速率慢，但是干燥速率变化平稳。砖坯在干燥过程中，由于水分的剧烈蒸发以及蒸发面的移动会造成局部温度梯度和湿度梯度过于集中，从而使得局部的水分饱和度呈上升状态。温度和湿度梯度过于集中均会导致应力集中，而且这种集中又特别容易发生在边角处，如果这一应力超过砖坯的承受极限就会产生缺陷。

图5 平行风速方向温度梯度及水分梯度示意

图6 入口2 m/s，573 K时坯料中心面温度以及水分饱和度曲线

2.2 入口风速2 m/s，温度573 K时干燥过程

图6所示的是干燥介质入口风速为2 m/s，温度为573 K时砖坯中心面的温度以及水分饱和度随着干燥时间的变化曲线。

与风速为5 m/s时相比可以看出，温升显然小很多，在干燥完全时的平均温度为330 K，而5 m/s风速下的为440 K。然而，水分的蒸发速率并没有因为入口风速的减小而减小，相反，与较大的入口风速相比，2m/s风速下干燥至S=0.1时所需要的时间为450 s，而风速为5 m/s时需要600 s，很显然风速并不是决定干燥时间的唯一因素。

图7和图8所示分别是垂直风速方向和平行风速方向截面温度及水分饱和度在不同干燥时间的梯度。

图7 垂直风速方向温度梯度及水分梯度示意

图8 平行风速方向温度梯度及水分梯度示意

从图中可以看出温度梯度的变化不是很大，并且水分梯度从一开始就很小。

对比5 m/s风速下迎风处和背风处的水分梯度相差不大，风速为2 m/s时，Z方向背风处的水分梯度明显要大于迎风处，原因是风速较小时背风面的回流以及扰流使得总的传热传质系数下降。与入口风速为5 m/s时相比，2 m/s风速下的平均温度梯度平均减小了40%，平均水分梯度平均减少了50%，较小的温度以及水分梯度可以有效防止热湿应力过大而产生的裂缝。

2.3 风速5 m/s，温度600 K时干燥过程

图9所示是入口风速为5 m/s，温度为600 K时中心面平均温度以及水分饱和度随着干燥时间的变化。

图9 入口5 m/s，600 K时坯料中心面温度以及水分饱和度曲线

从图中可以看出来，风速为5 m/s，温度为600 K时的温升很快，水分饱和度的降低也很快，在200 s之内水分饱和度S就降至了0.2，从水分饱和度变化曲线可以看出来，200 s之前的干燥属于快速干燥阶段。这种情况下虽然干燥速度很快，但是由于温度梯度和水分梯度产生的热湿应力也很大，砖坯很容易开裂。

图10和图11所示分别是垂直风速方向和平行风速方向截面温度以及水分饱和度在不同干燥时间的梯度。

与风速为5 m/s，温度为573 K相比，该条件下垂直风向和平行风向的平均梯度并不是一直都大。在200 s左右是一个临界点，在此之前，600 K时的梯度大于573 K时；200 s之后则慢慢变小，原因是600 K的条件下200 s后已经基本干燥完毕，这种情况下很容易出现大的热湿应力，从而引起裂缝。

图10 垂直风速方向温度梯度及水分梯度示意

图11 垂直风速方向温度梯度及水分梯度示意

3 干燥过程拟合

在干燥的几个主要阶段中，等速干燥过程会由于热质传递的不匹配而产生大的热湿应力，从而产生裂纹，而降速干燥阶段不会产生裂纹。因此存在一个干燥的临界点和干燥的最大温度梯度。根据文献中提及的干燥窑内截面水平温度梯度要求。本节将结合前文对不同干燥条件下水分及温度梯度的讨论，对特定温度梯度及水分梯度要求下的干燥进行最优化分析。本文规定温度梯度小于5000 K/m，即50 ℃/cm，时不会产生裂纹。因为干燥过程中有多个变量，现在控制一个量不变，即水分饱和度，观察干燥至水分饱和度为0.4时的干燥时间与干燥介质入口温度以及风速的关系。表1为不同风速以及温度下的最大梯度，灰色部分为临界干燥条件，即在此干燥条件下其温度梯度等于50 ℃/cm。

表1 不同干燥条件下温度梯度

根据表1可以绘制相应的二维图，横坐标为干燥时间，纵坐标为送风速度，如图12所示，图中灰点代表其温度梯度大于临界值50 ℃/cm，小五星点表示其温度梯度小于50℃/cm，黑点为临界点。阴影部分的温度梯度小于50 ℃/cm，该部分所对应的送风速度、送风温度为合适的干燥条件。

图12 不同干燥条件下的温度梯度

为得出干燥时间较小，并且温度梯度不大于临界点的送风风速以及温度，利用Matlab绘制三维曲面图如图13所示，用来描述温度，水分和干燥时间之间的关系，图中所用的干燥条件都满足其温度梯度小于临界点。

图13 不同干燥条件下干燥至S=0.4所需的干燥时间

从图13中可以看出，干燥时间最低的区域是温度为550～573 K，风速为0.5 ～2.5 m/s。

4 结语

本文针对多孔介质耦合传热传质，建立相关的数学模型，利用数值模拟软件FLUENT模拟了单块砖坯在不同边界条件下的干燥效果。

(1)详细考察在送风温度为573K、600K，送风风速分别为5m/s，2m/s时砖坯干的燥曲线以及不同方向的温度梯度和水分梯度。结果显示，相同温度下，随着干燥介质入口风速的增大，干燥时间呈减少的趋势，但是砖坯内部截面温度和水分梯度变大。同样，风速相同时，随着温度的升高，干燥时间减少，砖坯内部截面温度和水分梯度增大。

(2)在干燥过程中由于正面和侧面所处的边界条件不同，温度上升和水分下降的过程也不相同。在砖坯的边角处，由于同时存在两个方向上与外界的剧烈传热传质，温度上升和湿度下降的速度都要比其他地方快，因此边界处也是最易造成干燥不均匀并产生裂缝的地方。

(3)根据所得的干燥数据以及所确定的临界干燥温度梯度，利用Matlab绘制相关的三维图用来描述送风风速、温度与干燥时间的关系，并得出在小于本文所设定的临界梯度的条件下，即温度梯度小于50 ℃/cm时，较为经济的干燥区域是温度为550～573 K，风速为0.5～2.5 m/s。

Numerical Analysis of Coupled Heat and Mass Transfer in Porous Media within Adobe Drying Process in Industrial Furnace

Tongji University College of Mechanical and Energy Engineering Zhou Yu Qin Chaokui

Adobe drying is the inverse process of plastic molding and involves phase change and coupling of heat and mass transfer. From the perspective of unsteady heat and mass transfer and phase equilibrium, a three dimensional mathematical and computational model based on porous heat and mass transfer and computational fluid dynamics theories was established in this paper. By using the fluid calculation software FLUENT, the changing process of moisture content and temperature in the adobe was obtained. The influence of drying medium inlet velocity and temperature on drying process was discussed. Combining the drying datum and the critical drying temperature gradient, the relationship of drying medium inlet velocity, temperature and drying time was fitted by Matlab. Results show that the more economical drying condition is that the temperature is 550～573 K,and the wind speed is 0.5～2.5 m/s under the condition of less than the critical gradient set in this paper.

rdrying process, porous media, coupled heat and mass transfer, mathematical model