高温散体对流换热数值模拟Matlab编程实现

摘要：高温散体对流换热的理论求解一直是多孔介质研究的难点和热点问题，由于结构和流动的随机复杂性及非线性，许多情况需要采用数值求解。文章开发了高温散体对流换热数值模拟通用程序，实现了高温散体多孔介质对流换热离散代数方程组求解，根据高温散体对流换热过程气体、固体的流动和传热的相关参数，求解了对流换热过程的风速、风温、料温和压力场。通过可视化后处理技术绘制了不同高度的风温、料温的变化规律图、温度场图和与实验数据差异对比图。该研究为高温散体对流换热的研究和工业应用具有很好的指导意义。

关键词：Matlab编程；数值模拟；高温散体；对流换热

中图分类号：TP311.1     文献标识码：A

文章编号：1009-3044（2023）25-0004-04

开放科学（资源服务）标识码（OSID） ：

0 引言

高温散体广泛存在于化工、环保、冶金、建材、矿业等许多工业过程，而且在航空航天、核反应堆等领域也有着广泛应用，因此高温散体的高效热回收和强化传热问题近年来一直是多孔介质对流换热研究的热点和难点。人们对高温散体对流换热已经进行了一定的实验和数值模拟研究：马淑杰等[1]建立了烧结矿冷却过程的多孔介质模型，探讨了烧结矿热物性参数分别采用经验值和实验值对数值模拟的影响；夏建芳等[2]应用VB.NET语言，基于CFD平台，开发了环式鼓风冷却机烧结矿冷却过程自动仿真软件；刘立钧[3]采用有限差分法对烧结矿竖冷窑内气-固换热过程进行了数值模拟研究；倪文杰[4]将富氧和喷吹焦炉煤气与烟气循环联合使用，形成综合烧结工艺，通过数值模拟研究烟气循环条件下富氧和焦炉煤气喷吹参数对烧结工艺和产品质量的影响; 刘玮寅[5]等通过仿真模拟软件COMSOL Multiphysics对烧结矿竖罐式冷却过程进行数值模拟研究；冯军胜[6]等利用Fluent软件及其二次开发平台对烧结矿余热回收竖罐内气固传热过程进行了数值分析；果晶晶[7]等采用CFD软件对烧结矿余热罐内的气固换热过程进行了数值模拟研究; 冯军胜[8]等借助Fluent模拟计算平台，利用正交实验方法对烧结矿竖罐内气固传热过程进行了数值模拟与优化；吴礼忠[9]等基于Fluent模拟软件，通过使用UDS和UDF分别对函数方程和源项进行定义与修改，采用多孔介质模型对烧结矿层冷却过程进行模拟。本文选择了Matlab语言来实现对高温散体对流换热过程的数值模拟。

Matlab最早是专门针对矩阵简化运算开发的，所以被称为矩阵实验室（Matrix Laboratory） 。该软件具有强大高效的多维数据处理能力，可以进行数值分析、矩阵计算、科学数据可视化和非线性动态系统的建模与仿真等，具有代码简洁、编程效率高等特点，因此比其他同类软件操作简单方便[10-11]。同时软件还具有巨量的库函数，而且也可以方便定义自己库函数。

高温散体对流换热模型的数值模拟的基本步骤包括物理数学模型的建立、空间和数学方程的离散差分、代数方程的求解、计算机编程和结果分析处理。

1 高温散体对流换热模型的建立

1.1 高温散体冷却过程物理模型

目前工业高温散体对流换热常见工艺有气固逆流和气固错流两大类，一般固体为缓慢移动，本文以后者作为研究对象。假设某高温散体，温度在300～1 500℃之间，在冷空气中缓慢移动冷却直到降到某一特定温度下，冷却空气从高温散体下部垂直吹入散体料层，从上部进入烟罩或排放。整个过程可以简化为如图1所示的物理模型。

1.2 高温散体对流换热数学模型

针对图1所示的高温散体冷却过程，在料层较宽的情况下忽略边壁效应、宽度维数、长度维数，再假设高温散体由大小不同的球形颗粒随机堆积而成且各向同性，忽略对环境热损失、气体导热和辐射传热，高温散体三维数学模型就可以简化为一维非稳态模型，如式（1） 到式（6） 所示。

1） 连续性方程

[∂（ερgr）∂t+∂（ρgrwga）∂z=0]      （1）

2） 能量方程

①气体能量方程

[∂（ρgrcpgrTgr）∂tε+∂（ρgrwgacpgrTgr）∂z=Spvα（Tsr-Tgr）]  （2）

结合连续性方程可以简化为：

[ερgr∂（cpgrTgr）∂t+ρgrwga∂（cpgrTgr）∂z=Spvα（Tsr-Tgr）] （3）

②固体能量方程

[∂（ρsrcpsrTsr）∂t（1-ε）=（1-ε）∂∂z（λsr∂Tsr∂z）+Spvα（Tgr-Tsr）]     （4）

③动量方程

高温散体多采用Ergun公式计算，即：

[ΔPH=150（1-ε）2ε3μgνgd2e+1.75（1-ε）ε3ρgv2gde]  （5）

式中：[ΔP]—高温散体内部压降，Pa；

H—对应压降的高温散体高度，m；

[de]—当量直径，m。

④状态方程

[PV=NRT]  （6）

式中：[N]—气体摩尔数，[mol]；

[R]—气体常数，8.314[ J / （mol·K）]；

[M]—气体摩尔质量，[ kg / mol]。

式（5-24）也可以写为另一种形式：[PM=ρRT]。

3） 定解条件

初始时刻（[t]=0），料层入口处初始料温为[Ts0]，空气入口处（[z]=0），空气入口速度ug0.，空气温度为环境温度[Tg0]。

2 一维非稳态的数值求解

数值求解就是将时空离散化，把数学微分方程组转换为时空节点上的代数方程，最后解出所有节点的代数方程组，得到相应的温度时空分布。

2.1 离散代数方程组

数学方程离散方法有很多，而目前传热过程主要采取有限容积法，内部节点用有限差分，边角节点用热平衡法。针对一维非稳态模型，为了保证差分方程的稳定，采用隐式差分格式。

1） 连续方程离散

[ερn，k-ρn-1，kΔt+ρn，kwn，k-ρn，k-1wn，k-1Δz=0]   （7）

式中：[n]—时间节点编号；

[k]—空间节点编号。

2） 气体能量方程离散

[a]） 内部节点

[ερn，kcn，kTn，k-cn-1，kTn-1，kΔt+ρn，kwn，kcn，kTn，k-cn，k-1Tn，k-1Δz=Spvα（Tsn-1，k-Tn，k）]   （8）

[b]） 入口边界点：[k]=0，[Tn，k]=[Tg0]，[Tg0]为环境空气温度。

[c]） 出口边界点：[k]=[kmax]，根据能量平衡，忽略热损失，方程为：

[-Δtρn，kwn，kcn，k-1Tn，k-1+（εΔz2ρn，kcn，k+Δtρn，kwn，kcn，k+Δz2ΔtSpvα）Tn，k=εΔz2ρn，kcn-1，kTn-1，k+Δz2ΔtSpvαTsn-1，k]

（9）

3） 固体能量方程离散

①内部节点

[（1-ε）ρn，kcn，kTn，k-cn-1，kTn-1，kΔt]

[=（1-ε）λn，kTn，k+1-（λn，k+λn，k-1）Tn，k+λn，k-1Tn，k-1Δz2+Spvα（Tgn-1，k-Tn，k）]   （10）

②入口边界点：[k]=0，根据能量平衡，忽略热损失，方程为：

[ρn，kcn，kTn，k-cn-1，kTn-1，kΔt（1-ε）=（1-ε）λn，kTn，k+1-Tn，kΔz/2+Spvα（Tgn-1，k-Tn，k）]  （11）

式中，[k]=0时，[T0，k=Ts0]，[Ts0]为入料温度。

③出口边界点：[k]=[kmax]，根据能量平衡，忽略热损失，方程为：

[ρn，kcn，kTn，k-cn-1，kTn-1，kΔt（1-ε）=（1-ε）λn，kTn，k-1-Tn，kΔz/2+Spvα（Tgn-1，k-Tn，k）]  （12）

3 数值求解思路及Matlab编程实现

3.1 数值求解思路

数值求解主要包括前处理（网格划分）、代数方程求解和后处理（结果分析处理），总程序框图如图2所示。

3.2 数值求解编程实现

1） 网格和变量初始化

划分时间网格和空间网格，时间用t表示，空间高度用z表示，节点编号为1、2、3…… 设置空气、散体物性参数、运动参数和散体结构参数等，如：

...

ruos=1700;%sinter density

kxl=0.05*（15.5*a1+11.3*a2+9.1*a3+7.6*a4+6.7*a5

+6.3*a6）; %孔隙率

de=1/（a1/0.08+a2/0.07+a3/0.05+a4/0.0325+a5/0.0175+a6/0.005）;%平均直径Spv=6*（1-kxl）/de;%平均直径

ruog=ones（ml，nc）;

ruog（：，1）=MGAS*p（：，1）./（RGAS*Tg（：，1））;%air density boundary

Ts=700*ones（ml，nc）;

Ts（1，：）=811.6;%sinter inlet temperature boundary

...

2） 求解各时刻各空间节点温度和速度

①求解初始时刻各空间节点温度和速度

由于定解条件只知道[t=0]时的固体温度和[z=0]处的气体温度，初始时刻[z≠0]处的气体温度未知，所以先要假设初始温度场和压力场，进而算出速度场，采用迭代求解值反复修改初始温度场和压力场，直至温度场和压力场稳定，此时即得到初始时刻的真实温度，在此基础上就可以计算以后各个时刻的温度场了。该求解过程具体如图3所示。

②其余时刻各空间节点温度的计算

初始时刻的温度、速度和压力真值求出后，根据状态方程、初始压力、初始温度求出初始第一节点的密度，再结合初始速度利用Ergun方程求出下一个节点的压力，该压力再通过状态方程可以得到该节点的密度和速度，于是利用TDMA追赶法求出该时刻该节点的温度。就这样反复采用该方法就可以求出下一个空间节点，下下个节点直至所有空间节点的温度、压力、速度和温度。此处要注意的是物性均随温度而变化，在该节点温度未求出时物性是按上一时刻的温度来计算的，所以误差就会产生，尤其是时间间隔取得较大时更为明显。改变的方法是再加一层循环，即把现在计算出的温度再代回去重新计算，直到新旧温度间的差值小于指定精度为止。

部分代码如下：

...

%追赶法求固体温度

a（i，1）=0;

b（i，1）=（1-kxl）*ruos*cps（i，1）*deltaz/2+（1-kxl）*lamds（i，1）*deltat+Spv*h（i，1）*deltat*deltaz/2;