1 引言

单组元火箭发动机中的反应主要是催化反应，其催化室由随机分布的催化剂颗粒堆积而成(亦称堆积床)，推进剂在其中经过催化分解从而产生高温高压的气体。由于催化床内颗粒分布的随机性以及颗粒之间所形成的流通通道曲折而又复杂，故对其内部的流动与传热研究比较困难，很多研究都回避了这一问题^[1~3]。

催化床内的流动与传热模拟大多采用多孔介质模型，该模型是建立在试验研究基础上的，由于试验工况的限制，多孔介质模型的应用有很大局限性^[4]，所以有许多学者对堆积床内的流动开展了更为详细的实验研究^[5~8]和数值模拟研究^[9~21]。很多数值模拟直接对堆积床中的颗粒通道划分网格进行计算。这类研究中的早期计算大多是采用二维模型，为了模拟球形颗粒，堆积床内颗粒都在轴线上规则排列。例如Dalman等数值研究了二维轴对称圆管内两个球体之间的流动与传热，研究结果显示在两球之间的流体区域出现了旋涡，并且旋涡对两球之间的传热有一定的影响^[9]。Barry Lloyd等建立的二维堆积床模型中，堆积床颗粒采用等直径小球在轴线上排成一列的形式给出，研究结果表明球体之间的热导率随着两球之间距离的增大而减小^[10]。由于二维模型简化过度，难以近似三维流动，所以较近的研究大都建立三维模型。为了避免网格生成的难度，很多研究是针对颗粒规则排列的情况^{[11, 12]}，例如Min Sheng等研究了在规则、非接触排列的等直径球形堆积床通道内添加纤维的形式来改善床内的热传导^[13]。Jinliang Yuan等通过建立微型颗粒多孔介质模型和数值计算模型，对流体与固体之间的传热特性进行了分析，数值研究中颗粒也是采取了规则排列的等直径球^[14]。Alkhalaf A等通过建立规则排列堆积模型和多孔介质模型并运用商业软件FLUENT对横截面上气体组份分布进行了模拟^[15]。但是在很多实际的堆积床中，因为颗粒分布不规则，且颗粒之间为点接触，容易产生网格畸变，所以网格生成比较困难，许多研究者采用各种方法来降低难度。例如，于立章等通过对径向颗粒采用分离的方式排列，轴向的颗粒层之间采用分离与相交交替的方法排列来减小网格的划分难度^[16]；李健等通过台阶连接模型对圆球颗粒之间的接触点进行了简化^[17]；吴江权等通过将颗粒之间的点接触用短圆柱来代替以此来减小网格的划分难度^[18]；Alkhalaf A等综合了间隙排列、重叠排列和台阶连接三种方法来减小网格划分难度^[15]；郭雪岩等通过Chimera网格的方法来解决网格的划分难度^{[19, 20]}。Enrico Bianchi等建立了一种新型泡沫状催化床反应器模型，运用Openfoam软件计算分析泡沫催化床与壁面之间的距离以及流体的速度对泡沫催化床流动与传热特性的影响^[21]，但并未对网格生成过程做详细介绍。这些方法多数是对堆积床进行了简化(如颗粒有序堆积)，与实际(无序)堆积状态差别较大；个别针对无序堆积的计算网格生成工作量大，只适用于少量颗粒。

为此，本文在单组元发动机仿真研究中，针对颗粒随机分布的堆积床，提出了一种新的网格生成思路，并运用CFD方法对管道内等直径自然堆积小球的堆积床内部流动与传热特性进行了详细的分析。这种方法不仅避免多孔介质模型的过度简化，而且也避免了颗粒有序堆积与实际状况的偏差，可以直接应用于较多颗粒的堆积床仿真研究。

2 模型的建立

本文的研究对象如图 1所示，在圆柱形管道内部有一个充满颗粒的堆积床，流体从左侧进口流入，先经过一段管道，通过堆积床后再从出口流出。为了简化模型，堆积床两侧未设置挡板结构。管道直径为30mm，总长110mm，其中堆积床长度50mm，进口到堆积床的长度为40mm。

为了比较真实地模拟堆积床内的流动和传热，首先需要对其中的颗粒堆积过程进行仿真，确定颗粒自然堆积以后的空间位置；然后再划分网格进行数值计算。

2.1 堆积模型

计算空间中颗粒的堆积方法主要有：松弛算法、L-S方法和离散单元法(DEM)^[22]。本文采用离散单元法模拟颗粒在重力作用下的自然堆积过程。DEM的基本原理是将颗粒分离成离散单元的组合，通过运动学定律和接触模型，考虑颗粒间与壁面间的相互作用力，能够比较真实地反映颗粒堆积过程。堆积床和颗粒的参数如表 1所示。

堆积过程计算采用离散单元法商业软件(EDEM)：假定堆积床出口一侧向上，然后将颗粒从堆积床出口倒入，颗粒在重力作用下在堆积床中形成自然堆积。堆积完成后的颗粒分布如图 2所示，从图中可以看出，球形颗粒随机分布。

2.2 网格模型

由于整个堆积床中的颗粒是随机排列的，且颗粒之间的接触是点接触，如果划分贴体网格就需要人工操作，这对于大量的随机分布颗粒几乎是不能完成的，另外颗粒接触点附近的网格变形也会非常严重，使计算收敛困难。为此，本文提出了一种自动化的网格处理方法：即先不考虑颗粒的存在，首先对整个计算区域(包括堆积床区域)划分非结构网格，但是以六面体网格为主；然后逐个对自然堆积的颗粒进行扫描，通过颗粒的球心坐标和半径对网格进行计算判断，将处于颗粒内的网格标识为固体，其余网格标识为流体；这样处理后就可以采用CFD的方法对流体区域进行流动计算，而在计算传热时统一考虑流体和固体区域。

图 3显示了采用上述方法生成的固体颗粒网格。因为在这种网格处理方法中，不再对颗粒生成贴体网格，所以避免了人工操作的难度和网格变形，但是显然颗粒的近似形状与网格疏密有关。与颗粒尺度相比，网格越密则近似程度越高。本文所示算例的网格采用Gambit软件生成。

2.3 数值模拟方法及边界条件

数值模拟用FLUENT软件，采用三维N-S方程。湍流模型选用标准k-ε方程，其中各个系数的取值为默认值。计算方法为基于压力的求解方法，其中压力修正采用SIMPLE方法，其他物理量的空间差分采用二阶精度的迎风格式，松弛因子为0.5。

边界条件如下：

(1) 进口为速度边界，其速度值和温度值根据参考文献的实验条件给定。

(2) 出口边界为压力出口，其绝对压力为1个大气压。

(3) 壁面为无滑移、绝热壁面。

3 计算方法验证

网格模型中已经指出，本文提出的网格生成方法虽然自动化程度高，但是堆积床处的网格尺寸和颗粒大小相关，相对颗粒直径的网格越小，越能近似表示颗粒的形状。因此，首先对网格尺寸进行相关性验证，研究网格尺度对计算结果的影响。

计算时，堆积床内生成的网格主要为六面体网格，各方向尺寸基本一致，并且在整个堆积床内尺寸基本相同。网格尺寸分别为1/10，1/13，1/20，1/27颗粒直径(D_p)。根据参考文献[5]的实验，堆积颗粒为玻璃球，直径为5mm；流体介质为水；进口速度分别为0.1，0.2，0.3，0.4，0.5m/s。

图 4给出了不同进口流速与单位长度上的压降之间的关系，并将数值模拟值与实验值进行了比较。从图中可以看出，不同网格尺度下计算所得的规律与实验一致，当进口流速增大时，粗网格的偏差增大；随着网格尺寸的减小，计算值与实验值之间的偏差逐渐减小，当网格尺寸为1/20 D_p时，数值解已经非常接近实验值。通过对比可知，通过本文提出的计算途径可以进行颗粒随机分布条件下的堆积床流动模拟，计算结果符合自然规律；堆积床内网格为1/20颗粒直径是减小计算偏差的合适准则，再缩小网格的意义不大。

4 计算结果与分析

采用多孔介质模型计算堆积床流动时，需要设置流体在多孔介质内的摩擦系数，摩擦系数与单位压降和进口速度的关系为

$ f = \frac{{\Delta p}}{{\rho {U^2}}}\frac{{{D_{\rm{p}}}}}{L}\frac{{{\varepsilon ^3}}}{{(1-\varepsilon )}} $

(1)

式中ε=1-V_p/V为孔隙度，V_p为颗粒总体积，V为堆积管道容积；D_p为颗粒平均直径，U为来流速度。当采用Ergun方程计算压降时

$ \frac{{\Delta p}}{L} = 150\frac{{{{(1-\varepsilon )}^2}}}{{{\varepsilon ^3}}}\frac{{\mu U}}{{d_{\rm{p}}^2}} + 1.75\frac{{\rho {U^2}}}{{{d_{\rm{p}}}}}\frac{{(1-\varepsilon )}}{{{\varepsilon ^3}}} $

(2)

摩擦系数简化为

$ {f_{{\rm{Ergun}}}} = \frac{{150}}{{R{\mathit{e}_{\rm{p}}}/(1-\varepsilon )}} + 1.75 $

(3)

式中Re_p=ρuD_p/μ。

为了与多孔介质模型对比，本文采用DEM方法分别建立了粒子直径分别为5mm，6mm，7mm和8mm的自然堆积模型，并生成了网格尺度约为1/20颗粒直径的网格。经过计算，这几种不同颗粒的自然堆积模型相对应的孔隙率分别为：0.448，0.468，0.482和0.484。计算中，流体工质为100℃的液态水，颗粒初温为20℃；进口流速分别为0.1，0.2，0.3，0.4，0.5m/s。本文对堆积床内的流动进行了稳态流动模拟，对热水与堆积床的传热过程进行了非定常模拟。

4.1 颗粒直径对流阻的影响

图 5给出了四种不同颗粒直径的堆积床在相同进口流速情况下的压力云图，该图清楚地表明，流阻随着颗粒直径的减小而增大。图 6给出了不同颗粒直径的堆积床阻力特性与进口流速之间的关系。从中可以观察到：随着进口流速的提高，堆积床内的阻力压降逐渐增大；同一进口流速下，流阻随着颗粒直径的增大而减小；采用多孔介质的Ergun方程计算得到的流阻偏大，进口速度越高，Ergun方程的偏差越大。颗粒直径对流阻的影响可以从孔隙率来分析，大颗粒的堆积床孔隙率较大，流通面积随之增大，故流阻降低。Ergun方程是关于进口速度的二次多项式，由实验确定的系数缺乏通用性，随着进口流速增大，堆积床中湍流增强，偏差增大。

4.2 摩擦系数

为了比较本文方法与多孔介质方法的差异，图 7给出了根据本文结果计算的摩擦系数(应用公式(1))与Ergun方程所计算的摩擦系数(公式(3))的对比。从图 7所示的摩擦系数与雷诺数之间的关系曲线可以看到：摩擦系数随着雷诺数的增大呈递减趋势，当雷诺数大于3×10³以后，摩擦系数都趋于稳定值；颗粒的直径越小，堆积床的孔隙度就越小，所对应的摩擦系数就越大；应用Ergun方程计算出的摩擦系数只有在极低雷诺数时才会体现出颗粒直径的影响，而在雷诺数增大后不同颗粒直径的摩擦系数相差不大，这说明Ergun方程在流速较高时具有较大偏差。

4.3 孔隙率分布的影响

图 8是对颗粒直径为5mm，进口速度为0.5m/s的堆积床上某横截面的速度幅值云图和涡量强度云图，从图中可以看出，堆积床内部通道内出现了漩涡，并且涡量大的地方流速也高，分析其原因是孔喉作用使流速增大。

图 9给出的是不同颗粒直径堆积床某横截面上的平均速度分布。其取值是沿着圆周方向选取该截面上若干不同的半径，对相同半径位置的有效流体速度进行算术平均。图中横坐标(R-r)/D_p为无量纲壁面距离(R为管道半径，r为到管壁的距离，D_p为颗粒直径)。从图中可以看出，速度幅值的分布趋势是左边比右边高，即靠近壁面附近的平均流速比其他地方的流速稍大，这是由于堆积床内孔隙率的分布不均匀所造成的。因为等直径颗粒的堆积床在壁面附近出现了较大孔隙率，所以流阻减小，速度增大，这种现象为“壁面效应”，局部孔隙率较大所引起的较大流通量称之为沟流^[20]。从图 9中还可以观察到，颗粒直径越小，速度极差(最大值减去最小值)越小，分布越均匀，因此直径小的堆积床可以改善孔隙度分布的不均性，同样也可减小壁面效应和沟流。

4.4 管道内的传热过程

本文计算了催化剂颗粒直径为5mm，管道的进口速度为0.1m/s，进口水温为100℃，初温为20℃条件下，水流过堆积床并与之换热的过程。图 10显示了不同时刻的管道纵剖面温度分布。从图中可以观察到，颗粒的温度随着时间推移逐渐升高，且堆积床下游中心处的温度要稍低于四周的温度。这是因为堆积床内孔隙度分布不均匀，壁面附近的“沟流”使热流体先到达，造成传热未平衡时壁面附近的温度比中心温度高。

图 11给出了t=0.5s时，堆积床某截面上的速度幅值云图和温度云图。图 12给出了t=0.5s时，不同颗粒直径堆积床某截面上温度沿径向的分布。这两个图中可以清晰地看到由于横截面上速度分布的不均匀，壁面效应使得壁面附近的速度明显较高，温度也较高。但是随着时间推移，当流体和颗粒之间逐渐达到热平衡时，这种不均匀性会逐渐消失。

5 结论

本文通过对自然堆积的催化床建模，运用CFD方法直接求解了堆积床内部的流动和传热过程，并与多孔介质的Ergun方程进行了对比分析，得到以下结论：

(1) 本文提出的针对颗粒随机分布的催化床内流动与传热计算思路，可以有效地减小网格生成难度，适合进行大规模颗粒堆积的堆积床仿真，在网格尺度小于颗粒直径的1/20时，与实验值符合较好。

(2) 管道进口流速增大、颗粒直径减小都会使流阻增大；颗粒直径减小还会使摩擦系数增大；而Ergun方程在高雷诺数下计算的流阻偏大。

(3) 本文计算方法可以有效模拟出孔隙率不均匀所引起的堆积床壁面效应和沟流，减小颗粒的直径可以削弱孔隙度分布的不均匀性；由于壁面效应的影响，当高温流体流过堆积床非定常传热时，壁面附近的温度会高于中心温度。