引言

规整填料由于具有传质效率高、压降低、操作范围广等优点，广泛地应用于化工生产中。一般而言，填料传性能越好，填料层压降越高。规整填料由于具有较大的比表面积与较高的空隙率(低持液量)，在保持低压降的同时可具有较高的传质效率。自从上世纪60年代丝网规整填料投入工业应用至今，许多科研机构与商业公司均研发了新型、高效的规整填料以代替传统的X、Y型填料^[1]，以满足工业生产中的需求。

经过多年发展，人们对规整填料塔内的持液量、有效润湿面积、流动分布等进行了充分研究，归纳总结出了很多的经验、半经验半理论或理论模型，如SRP模型、Delft模型、扩散模型、单元网格模型、结点网络模型等。然而这些模型的建立一般是基于特定填料大量的实验数据及对填料结构的简化，应用范围受到很大的限制。近年来，采用计算流体力学(CFD)数值模拟研究填料的流体力学性能和传质性能已成为一个趋势^[2-3]。利用CFD方法既可得到宏观信息，也可观测到微观现象，为深入分析填料塔内的流动机理提供了依据。CFD数值模拟避免了实验设计和实验操作过程，节省了时间及操作费用。

早期的CFD模拟只涉及单相流。van Gulijk^[4]建立了“toblerones”三维简化模型，把多相流简化为单相流进行模拟，研究装有催化剂的填料内部的流体流动。其后很多研究者以“toblerones”作为几何模型来模拟填料的流体力学性能^[5-6]。Larachi等^[7]把填料层的不同位置分解成4个代表单元，分别模拟4种代表单元内的气相流动，得到每个单元的圧力损失系数，进而计算整个填料层的压降。Li等^[8]在文献[7]模型的基础上，将填料层的压降分解为5种不同结构单元产生的压降之和，研究了新型SiC规整填料内流体流动。对于更加贴近实际的两相流模拟，Szulczewska等^[9]以250Y填料几何参数为基础，建立了板波纹二维两相流模型，利用Fluent软件研究了气液两相逆流流动的过程。随后，Raynal等^[10]将填料内气液的流动轨道简化为“Z”字状，利用VOF模型得到液膜厚度，通过计算得到了填料的持液量。Shojaee等^[6]研究了气液流量对填料有效润湿面积的影响，并将实验结果与已有的关联式进行了对比。

由于填料结构复杂，计算所需网格数目庞大，受计算机计算能力限制，CFD模拟一般都局限在局部填料，不可能进行工业尺度的模拟计算。在这种情况下，填料的多尺度模拟方法提供了一种解决思路，其核心思想是在微观模拟中得到填料局部流动细节，随后在这些微观结果的基础上进行宏观模拟计算从而得到整个填料层的压降、持液量等宏观信息。Raynal等^[11]最先提出基于CFD用多尺度模型来计算大型规整填料塔的压降。Sun等^[12]与本课题组^[3]均提出了不同的多尺度模型，用来实现在计算机可接受的计算量下，模拟工业尺度填料的流体力学性能。

本文采用本课题组之前提出的多尺度方法^[3]，分别在二维、三维几何模型中对新型HSX填料的液膜厚度、有效润湿面积和塔压降等流体力学性能进行模拟研究，并将模拟计算得到的压降与实验压降进行了对比。

1 CFD控制方程

本文在模拟过程中忽略温度场及传质，假设气相为理想气体，液体不可压缩。采用流体体积(VOF)两相流模型追踪气液相界面。VOF模型在整个控制域内只需求解一组动量方程，气液两相使用同一速度场，动量方程为

$\frac{\partial }{{\partial t}}\left( {\rho \mathit{\boldsymbol{u}}} \right) + \nabla \left( {\rho \mathit{\boldsymbol{uu}}} \right) = - \nabla \mathit{\boldsymbol{p + }}\nabla \left[ {\mu \left( {\nabla \mathit{\boldsymbol{u + }}{\mathit{\boldsymbol{u}}^{\rm{T}}}} \right)} \right] + \rho \mathit{\boldsymbol{g}} + \mathit{\boldsymbol{F}}$

(1)

式中，ρ为流体密度，kg/m³；u为流体速度；μ为流体黏度，Pa·s；p为压力。

对于气液两相，出现在动量方程中的物理参数由每个控制体内的分相决定

$\rho = {\alpha _{\text{L}}}{\rho _{\text{L}}} + \left( {1 - {\alpha _{\text{L}}}} \right)\rho $

(2)

$\mu = {\alpha _{\text{L}}}{\mu _{\text{L}}} + \left( {1 - {\alpha _{\text{L}}}} \right){\mu _{\text{L}}}$

(3)

式中，α为体积分数；下标L、G表示液相、气相。在每一控制体积内，体积分数总和为1

$\sum\limits_{i = 1}^n {{\alpha _i}} = 1$

(4)

在VOF模型中，气、液相界面的追踪是通过求解连续性方程来完成的，连续性方程为

$\frac{\partial }{{\partial t}}\left( {{\alpha _i}{\rho _i}} \right) + \nabla \left( {{\alpha _i}{\rho _i}{\mathit{\boldsymbol{v}}_i}} \right) = 0$

(5)

在模拟过程中，方程(1) 中源项F为表面张力, Fluent基于Brackbill等^[13]提出的连续表面力(CSF)模型来计算表面张力。

多孔介质模型中的阻力通过在动量方程中添加源项来实现。源项由黏性阻力相和惯性损失相构成，即

${S_i} = - \left( {\sum\limits_{j = 1}^3 {{\mathit{\boldsymbol{D}}_{ij}}\mu {v_j}} + \sum\limits_{j = 1}^3 {{\mathit{\boldsymbol{C}}_{ij}}\frac{1}{2}\rho \left| v \right|{v_j}} } \right)\quad \left( {i = x,y,z} \right)$

(6)

其中，D、C分别为黏性阻力及惯性阻力系数矩阵，即

${S_i} = - \frac{\mu }{r}{u_i} + {C_2}\frac{1}{2}\rho \left| u \right|{u_i}\quad \left( {i = x,y,z} \right)$

(7)

其中，r为渗透率；C₂为惯性阻力系数。在模型中也可以用幂函数来表示阻力，即

${S_i} = - {C_0}{\left| v \right|^{{C_1}}} = - {C_0}{\left| v \right|^{{C_1} - 1}}{v_i}\quad \left( {i = x,y,z} \right)$

(8)

式中，C₀和C₁为经验系数。

选用RNG-k-ε模型为湍流计算模型

$\begin{gathered} \frac{\partial }{{\partial t}}\left( {\rho \varepsilon } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho {\mu _i}\varepsilon } \right) = \frac{\partial }{{\partial {x_j}}}\left( {{\alpha _\varepsilon }{\mu _{{\text{eff}}}}\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right) - \hfill \\ {C_{1\varepsilon }}\frac{\varepsilon }{k}\left( {{G_{\text{k}}} + {C_{3\varepsilon }}{G_{\text{b}}}} \right) - {C_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k} - {R_\varepsilon } + {S_\varepsilon } \hfill \\ \end{gathered} $

(9)

其中，ε为湍流耗散率，G_k为由平均速度引起的湍动能产生项，G_b为由浮力引起的湍动能产生项， ${R_\varepsilon } = \frac{{{C_\mu }\rho {\eta ^3}\left( {1 - \left( {\eta /{\eta _0}} \right)} \right)}}{{1 + \beta {\eta ^3}}}\frac{{{\varepsilon ^2}}}{k},\eta = \sqrt {2{S_{ij}}{S_{ij}}} \frac{k}{\varepsilon }$ 。模型中的常数C_1ε, C_2ε, C_μ, η₀和β分别为1.42, 1.68, 0.085, 4.38及0.012。

本文使用Ansys15.0软件进行模拟计算。模拟过程中，水与空气作为气液两相。模拟为非稳态，对流项采用二阶迎风格式，压力项采用PRESTO!算法，压力-速度耦合方程的求解采用PISO方法。VOF模型中采用Geo-Reconstruct界面重构技术对气液相界面进行追踪。模拟中时间步长均设置为1×10^-5 s。计算过程中，设置监视器检测气液相出口流量，当气液两相出口流量围绕某一值变化不大时，可以认为计算达到拟稳态状态。

2 几何模型及边界条件

将HSX丝网填料的流体力学性能与传质性能的实验结果与Mellapak125X型填料进行了对比^[14]，实验表明HSX填料相对于传统填料有更高的传质性能和更低的压降。但是，对于填料内液相的分布情况、流体的流动轨迹等详细信息仍然未知。本文对HSX125板波纹填料建立几何模型，利用多尺度方法对其进行CFD模拟，并与实验数值进行对比。HSX125填料的几何结构见图 1，几何参数见表 1。

建立填料多尺度模拟所需几何结构(图 2)。建模过程中，结构参数应与实际填料尽量保持一致。在2D模型(图 2(a))中利用VOF方法追踪气液接触面，得到液膜厚度δ；3D模型(图 2(b))取相邻的3个流道作为计算域，模拟得到填料的有效润湿面积；用膜厚与有效润湿面积来计算填料持液量，利用多孔介质模型(图 2(c))在Fluent中得到填料塔的湿塔压降。2D模型中，液相入口和气相出口宽度分别为2 mm和3 mm。3D模型中，液相入口是垂直于填料片方向宽2 mm沿填料片横向延伸的狭缝，以确保液相是以均匀液膜形式进入模型计算域的。

采用结构化网格，对近壁层网格进行加密。Raynal等^[10]建议液膜内最少有8~10层网格，本文中第一层网格的厚度仅为0.005 mm。随后进行网格独立性验证，为了节省计算时间，最终确定2D与3D网格数量为36587与503284。

设置气、液出口的边界条件为压力出口(0 Pa)，气、液进口均为速度进口边界条件：u_L=u_L0、u_G=u_G0，无滑移壁面。

3 模拟结果与讨论 3.1 干塔压降

干塔压降是湿塔压降模拟的基础，采用多孔介质模型对全塔进行模拟计算，得到干塔压降。模拟值与实验值如图 3所示。模拟值与实验值误差平均在5%之内，吻合度较好。利用CFD方法可以比较准确地得到填料塔干塔压降。

本文多尺度模拟使用虚拟单相流计算湿塔压降，虚拟流的有效速度为

${u_{\text{e}}} = {u_{{\text{e}},{\text{G}}}} + {u_{{\text{e}},{\text{L}}}}$

(10)

其中 ${u_{{\text{e}},{\text{G}}}} = \frac{{{u_{\text{G}}}}}{{\varepsilon \left( {1 - {h_{\text{L}}}} \right){\text{sin}}\theta }},{u_{{\text{e}},{\text{L}}}} = \frac{{{u_{\text{L}}}}}{{\varepsilon {h_{\text{L}}}{\text{sin}}\theta }}$ 分别代表气相与液相的有效速度，θ是填料变化角度的算术平均值，h_L是填料的持液量，其值由液膜厚度δ及有效润湿面积a_e求得

${h_{\text{L}}} = {a_{\text{e}}}\delta $

(11)

3.2 平均液膜厚度的影响因素 3.2.1 液相速度

首先研究液相速度对平均液膜厚度的影响，设定气速为0。图 4显示了不同时间液膜在填料表面的流动状况。

图 4的模拟结果显示，液膜的厚度在整个模拟域内具有不均一性。模拟过程中发现入口处液膜较厚，在液相向下流动的过程中，液膜厚度逐渐减小；在填料折线处，液体有部分累积。进一步利用Tecplot360软件对流场相分率图进行处理，得到气液两相的相界面。图 5给出了不同液相雷诺数下的气液相界面的位置，从而可得到液膜厚度。

由图 5可以看出，随着液速的增大，液膜厚度逐渐增大。利用得到的液膜厚度与液相雷诺数作图(图 6)，可以看出在不考虑气相作用时，平均液膜厚度呈线性增长。

3.2.2 气相速度

固定液相速度为0.01 m/s，模拟不同气相速度下液膜厚度的变化如图 7所示。

从图 7中可以看出，液膜厚度随着气速的增大逐渐增大，在气速较小时增加平缓，接近液泛区域(F>2.5(m/s)·(kg/m³)^0.5)时，液膜厚度急剧增加。这是由于气速较小时，气体对液体影响很小，液膜的流速较大，膜厚较小；随着气速增大，气体对液体向下流动产生了较大的阻碍作用，从而导致液体累积，液膜增厚。

3.3 有效润湿面积的影响因素

利用CFD模拟可以较为精确地得到填料的有效润湿面积，从而降低了在以往研究过程中把气液两相有效润湿面积与实际填料表面积之比(a_e/a_p)当作1而产生的较大误差。

3.3.1 液相速度

设定气速为0，研究液相速度对有效润湿面积的影响。

由图 8可以看出，液体以膜状流动进入填料层，在重力与表面张力的作用下沿填料表面向下流动。当喷淋密度较小时，液体以溪流状态沿填料的流道底部流动, a_e/a_p远小于1。随着喷淋密度的增加，填料表面被液体覆盖区域逐渐增大，虽然还是以溪流为主，但已经有膜状流动区域存在，并且液体有向相邻流道流动的趋势；喷淋密度继续增加，整个填料表面逐渐被覆盖，a_e/a_p接近1。图 9给出了a_e/a_p的值随液相雷诺数的变化趋势。

3.3.2 气相速度

固定液相速度为0.01 m/s，模拟不同气相速度下有效润湿面积的变化，如图 10所示。

有效润湿面积与平均液膜厚度受气速影响变化趋势一致，在低气速区域增长缓慢，在液泛区域(F>2.5(m/s)·(kg/m³)^0.5)增加明显。其原因也是由于随着气速的增加，气液交互作用逐渐增强，从而增加了有效润湿面积。

3.4 湿塔压降模拟结果

将3.2节与3.3节模拟计算得到的平均液膜厚度、有效润湿面积代入式(10)、(11) 中，计算求取虚拟单相流的有效速度。在全塔几何模型中，采用多孔介质模型对全塔湿塔压降进行模拟计算。图 11描述了喷淋密度L=39.34 m³/(m²·h)时模拟结果与实验结果对比情况。

由图 11可以看出，在整个测量范围内，模拟值与实验值趋势相同，模拟结果与实验结果平均误差在15%左右，表明此多尺度模型可用来对填料的压降进行预测。分析造成误差的原因，可能包括如湍流模型不适用于低速流动状态、VOF模型对于湍动较大的流场有较大误差等。

4 结论

(1) 当气相速度为0时，液膜厚度随液相速度呈线性递增变化。低气速对液膜厚度影响较小，当气速增加时，接近液泛区域，液膜厚度随着气速的增加明显增大。

(2) 有效润湿面积随着气速与液速的增大均增大。当不考虑气相影响时，低喷淋密度下气液两相有效润湿面积与实际填料表面积之比(a_e/a_p)仅有0.3左右，远小于1。固定喷淋密度研究气相对有效润湿面积影响，结果表明，气相对其影响与对液膜厚度影响一致，高气速(泛点附近)下由于气液交互作用增强，有效润湿面积增加显著。

(3) 湿塔压降模拟结果表明，模拟计算值与实验值吻合度较高，整体平均误差基本在15%以内。