李清元 赵洪涛

(上海三爱富新材料股份有限公司，上海 200241)

四氟乙烯精馏塔流体力学数值模拟

李清元 赵洪涛

(上海三爱富新材料股份有限公司，上海 200241)

通过实验对空气-水体系在填料塔内的流动行为进行了流体力学研究，运用计算流体力学方法对丝网规整填料CY建立了数值模型，研究了不同喷淋密度和气速下，压降、液膜厚度和液膜长度的关系，进而对四氟乙烯提纯填料塔进行了数值分析，使得工业生产装置数据化和可视化，为工艺技术和设备的改进提供了有力的支持。

流体力学；丝网填料CY；四氟乙烯；填料塔

0 前言

四氟乙烯(TFE)单体精馏生产装置采用规整填料塔的生产工艺基本成熟，但是单耗相对国外来说还是偏高；单体塔等设备比较大，造成操作费用和设备成本较高。由于四氟乙烯体系的特殊性以及较高的操作压力，造成难以直接通过实验获得塔内流体流动的有用数据，同时对四氟乙烯精馏先进技术也极少见公开报道。

目前,用计算流体力学方法(CFD)[1-2]研究规整填料塔内复杂的流体流动成为填料塔研究的发展趋势之一。利用流体力学软件Fluent 6.3对塔径Φ 400 mm、填料高度2 000 mm的填料塔进行流体动力学数值分析,摸清流体在填料内部的流动情况，建立正确的数学模型，然后对四氟乙烯提纯填料塔进行数值模拟，分析生产装置内流体流动情况，以便为单体精馏技术的进步提供指导。

1 规整填料塔流体流动模型

对规整填料的压降、持液量、液膜流动和传质等研究已有一定的成果,但对气液两相流动的相互作用研究较少,对精馏过程气液两相流体力学行为研究很不足,阻碍了分离效率的进一步提高。影响分离效率的因素有很多，包括填料类型、填料盘和填料段高度、安装方式、流体物性、液体喷淋密度和F因子、气液分布情况等。因此,很有必要对这些参数变化时规整填料内气液两相的流动情况进行研究。

填料塔内气液两相的流体流动对气液两相的传质和传热都会产生很大的影响。采用计算流体力学方法可以对填料塔中流体的流动行为进行非常精确的模拟。另外，规整填料的几何结构非常复杂,没有办法建立整塔的实体模型。

Raynal[3-4]等建立了Z状流动路线的二维两相流模型,模型的计算单元和边界条件见图1。

图1 填料二维模型

1.1 气液两相流数值模拟方法

在进行流体力学模拟时，选用合适的流体力学模型至关重要，包括确定流体相态和流动状态等。对于规整填料塔内气液两相逆流流动的流体力学来说，流体为多相流，流动状态含层流和湍流。

由于研究对象复杂,仅对规整填料内气液两相并流流动的动量传递性质进行了研究,没有考虑传热和传质的影响。为了模拟液膜波形变化,考虑到液膜具有自由表面,属于气液两相流动,所以选用VOF法来寻找气液相界面。由于液膜靠近规整填料壁面的外层区域为湍流流动,靠近壁面的内层区域为层流流动,拟选用RNGκ-ε两种湍流模型封闭动量方程。填料近壁处理方法采用增强函数方法。

1.1.1 流体体积函数法

VOF法是在MAC方法的基础上发展起来的。VOF方法用相函数取代了 MAC方法中虚拟无质量的彩色粒子,是该方法的改进。VOF模型是应用于固定的Euler网格上的两种或多种互不溶流体的界面追踪技术。在VOF模型中,各相流体共享一个动量方程,每一相的体积分数在整个计算域内被追踪。假定规整填料内气液两相逆流的流动过程是非稳态的气液两相分层流动的过程。而此种流动的关键是如何确定气液相界面的位置。该法适用于存在非单一表面和较陡的自由表面的情况，同时所需的计算机存储单元较少，因此，在处理含复杂自由表面流体力学的问题上具有很大的优势。

一般假定气相为理想气体,液相为不可压缩流体,且都为牛顿型流体，在物理模型内,气液两相都是恒温的。对于第q相,连续性方程[5-7]为：

(1)

(2)

式(2)中的F为动量方程中的动量源项。根据填料塔中气液两相逆流流动的特点,应当考虑表面张力动量源项Fvol和气液相间作用力动量源项FLG的影响,即：

F=Fvol+FLG

(3)

1.1.1.1 表面张力动量源项

在规整填料塔内,液相的主要流动状态为液膜和液滴，而液膜在填料表面的铺展程度和液滴的大小都会受到表面张力的影响。对于气液两相的流动过程,表面张力源项可由式(4)计算：

(4)

式(4)中，κL为液相界面的曲率。同时，在Fluent 6.3中，假定流体与壁面产生的接触角常用于调整壁面附近单元表面的法向,而不是加强壁面本身的边界条件。这个接触角与一个单元正常计算的表面法向远离壁面的联合决定了表面的局部曲率,这个曲率常用于调整表面张力计算中的体积力项。在Fluent 6.3中,系统默认采用Brackbill等提出的CSF模型计算表面张力源项。表面张力动量源项可以通过软件中设定液相的表面张力直接加入。

1.1.1.2 气液相间作用力动量源项

在气液两相逆流时,气相流动会导致液相表面波动增强，从而使得气液两相的有效接触面积发生变化, Woerlee等发现气液相间的作用力将会导致部分压降的产生,造成流体能量损失。其压降为：

(5)

式(5)中，AL为液膜的横截面积，SL为液相润湿通道的周长。

由于规整填料结构复杂,假定为一系列的圆管或湿壁塔,液相的流动为环状流或膜状流,且其液膜厚度不变。通过实验获得了部分有关气液两相流动的作用力模型,如 Bird和 Hewitt。

采用Hewitt拟合出的气液相间作用力模型,其表达式如下：

(6)

(7)

式中，δ为液膜的厚度，可由Nussel降膜理论[8]计算如下：

(8)

气液相间作用力的表达式如下：

式中反映了气液相界面附近单元内的相间作用力源项。对于气液相界面以外的单元,不存在相间作用力源项。气液两相间的作用力可通过编写的用户自定义函数(UDF)从 Fluent 6.3动量源项中直接导入。

1.1.2 RNGκ-ε模型

RNGκ-ε模型是通过对瞬态的Navier-Stokes方程进行重正规化群 (Renomalization Group Theory)的数学方法进行推导求得的，其模型为：

(10)

(11)

上式中Gk是平均速率梯度产生的湍动能项,S是平均应力张量系数,Gb为由于浮力产生的湍动能项,对于恒温过程此项为0，YM为可压缩湍流中由于流体膨胀产生的湍动能,在本章计算中由于流体均假设为不可压缩流体,所以此项也可忽略。ueff为有效黏度，模型常数C1ε取1.42，C2ε取1.68。

RNGκ-ε模型使用严格的统计学方法得到。在ε方程里有一项是专门用来考虑快速变形的流动，因此对于快速变形的流动更加准确。同时考虑了湍流漩涡的影响,因此可以提高模拟涡旋流的准确性。方程中湍流 Prandtl数是一个解析的公式。模型包括了低雷诺数的影响,适用于低雷诺数流体流动的模拟。

1.1.3 处理近壁区模型

湍流模型适用于模拟处于湍流核心的流体的流动状况,但当流体在近壁区流动时,由于墙的存在会使湍流和漩涡消失,因此墙的存在会对流体的流动状况产生非常大的影响。由于κ-ε湍流模型基本上只在湍流核心区域才有效,因此，在选用这种模型时,对规整填料近壁区低雷诺数流动时要选择适当的壁面处理方法。增强壁面处理法能很好地描述近壁区的流动。本文在Fluent中设置增强壁函数法来处理这种受黏性底层影响较大的液膜流动。

2 实验内容

在Φ 400 mm有机玻璃塔中，装有CY700丝网波纹填料或聚丙烯塑料Φ 25 mm阶梯环填料，填料高度为2 000 mm，进行实验及数据测试，流程如图2所示。气相流量通过鼓风机产生，调节排气阀开度，控制气体体积转子流量计进塔流量，塔顶气体放空；液体流量由水槽通过离心泵打入塔内，调节旁路阀，控制液体体积转子流量计进塔流量，塔底液体排入水槽循环使用，调节气液比，找出不同工况下的载点和泛点以及压降和持液量。

图2 填料塔流体力学实验流程

3 丝网规整填料CY数值处理模型

对于规整填料,Zogg提出了一个预测由层流转变为湍流的临界雷诺数公式：

(12)

式中,α为波纹两斜边的夹角,εp为填料的空隙率,θ为填料波纹与竖直方向的夹角。对于不锈钢规整丝网填料CY[9]来说, 其结构参数如表1所示。

表1 丝网规整填料CY特性参数

由表1可知，α为81°，εp为0.85，θ为45°，代入式(12)可得到Rec为140.3，因此,对于雷诺数大于140.3 的流动,可认为其流型为湍流,而雷诺数小于140.3的流动，则认为是层流。雷诺数的计算方法[10]见式(13)、式(14)：

(13)

(14)

式中,uge和ule分别为气液相的有效速率，计算方法见式(16)、式(17)。由规整丝网填料参数，取dg为3.3 mm，液相水力直径,取为规整填料比表面积的倒数为1.4 mm。

(15)

(16)

(17)

hl=δα

(18)

式中，dg为气相水力直径,对规整填料空塔水力直径由式(15)计算，其中b为规整填料波长(m)，h为波纹高度(m)，s为波纹边长(m)，hl为持液量。

3.1 丝网规整填料塔流体力学建模

对于空气-水体系填料塔流体动力学数值模拟，其物性参数如表2所示。

表2 空气-水物性参数

3.1.1 模型网格划分

规整丝网填料CY采用了与计算域相适应的四边形结构网格[11]。考虑到界面波动的复杂性,在垂直于壁面方向的液相区网格最小间距为0.1 mm,且一般要保证液膜内要不少于8个网格。从相关填料宏观波纹板结构及微观结构的研究,得到液膜流动的结构尺寸都在毫米数量级以上,其表面的液膜厚度不足1 mm,在液膜流动方向网格间距为0.1 mm,填料间网格形式如图3所示,所划分的网格总数为128 379 个。

主要研究气液两相的流动过程,模型取10个波纹进行研究,模型采用三角形波纹连接方式,波纹长度为 10.3 mm,波纹振幅为 4.3 mm，模型总长度为103 mm,包括气相出口(4.3 mm)、液相入口(1.4 mm)、液相出口(1 mm)和气相入口(3.3 mm),其余各面为壁面。

图3 CY填料二维模型

3.1.2 求解器的设置

将建立好的物理模型和网格导入Fluent 6.3之后,要对Fluent 6.3求解器进行设置。首先要对读入的网格质量进行检查并将求解器的计算单位修改为mm。Fluent 6.3求解器选用压力的分离式求解器(Pressure Based)，采用非定常计算方法；对流项采用一阶迎风格式，多相流模型选用VOF模型，湍流模型选用RNGκ-ε湍流模型,参数为默认；流场计算选用SIMPLE算法，对于基本方程的离散差分采用如下方法：压力项采用PRESTO算法；动量方程、湍动能方程和湍动能耗散方程都釆用二阶迎风格式,体积分数采用Modified HRIC格式。时间项采用隐格式，压力项的松弛因子为0.2,动量为0.5,湍动能和耗散率为0.2,其他项均取为0.5,时间步长一般可取0.000 03～0.001,在每个时间步内,当速度、湍动能和耗散率的残差都降至0.001以下时,可认为在该时间步内的流动达到稳定状态。在液相流动过程中,要考虑其重力的影响,设定重力方向为y轴负向,大小为9.81 m/s2。

计算采用Dell电脑进行运算,主要配置如下:core 2核,4 GB内存,500 GB硬盘,CPU 2.66 GHz。求解过程所需的时间较长,计算过程中每个模拟条件下的时间都在72 h以上,有时超过120 h。

3.1.3 初始条件

初始时,假定计算域内不存在液相,只存在静止的气相，αL=0，αG=1。操作压力为常压和加压。

3.1.4 边界条件

规整填料内气液两相逆流流动过程,对于气液两相,均有2个进口和2个出口。设置液相入口为速度入口(Velocity Inlet),对于填料塔，湍动能κ和湍动能耗散率ε采用经验关联式[12]计算如下：

κ=(0.003～0.005)u2

(19)

(20)

式中，u为进口速率(m/s),l为填料特征尺寸(m)，一般为入口直径或长度的0.005～0.05。液相出口为压力出口(Pressure Outlet)；气相入口也为速度入口(Velocity Inlet),气相出口为压力出口(Pressure Outlet)；壁面为无滑移壁面(Wall)。

4 结果分析

4.1 单液相数值模拟

采用流体动力学软件Fluent 6.3分别对喷淋密度为20、28.2、36、44和55 m3/(m2·h)的液相进行了数值模拟，其液相在填料表面流动情况如图4所示。

图4 填料表面挂膜情况

从图4可以看出，喷淋密度在20～55 m3/(m2·h)范围内，液体都能在填料上成膜状流动。液相的流动由于受表面张力的作用在起初以液膜状态流动，但随着流动过程的进行,当液相所受重力高于表面张力时,液膜发生断裂,液相不再是连续的液膜。

模拟的液膜厚度和喷淋密度的关系见表3。从表3可以看出，液膜厚度随着喷淋密度的增加而增加。对气液两相逆流流动的情况，学术界认为液相在规整填料内以液膜的形式流动，液相是以溪流的形态流动的,而溪流的形状只和液体本身所受的作用力有关,与填料自身无关,而与填料的材质有关。通过分析可知，模拟值和计算值存在的误差为4.5%，模型基本反映了实验填料塔的运行工况。

表3 喷淋密度与液膜厚度关系

另外，不同的喷淋密度，液膜竖直高度也不一样，喷淋密度和液膜高度关系见表4。由液膜高度可以推断，其填料盘高也可以按照同样的液膜高度分段，然后上下层填料盘旋转90°交叉堆放。

表4 喷淋密度与液膜高度及填料盘高关系

4.2 气液两相流数值模拟

对喷淋密度为28.2 m3/(m2·h)、气速为0.16～0.49 m/s范围内的气液两相进行了数值模拟，其液体成膜和两相速度以及压降分布见图5。喷淋密度为28.2 m3/(m2·h)时，液膜高度为3.5个填料波段，故流体力学建模取3个填料波段进行模拟。

图5 液相挂膜情况和气液两相流速

从图5可以看出，在喷淋密度为28.2 m3/(m2·h)时，不同空气的流速产生的模拟压降和计算值的偏差约在10%，其中液泛压降模拟值与实验值差别较大，主要是由于实验中液泛只是局部液泛，气体还是可以顺利通过塔填料从塔顶排出，造成液泛压降较小；而模拟液泛过程是液体不能顺利流下，气体不能顺利从顶部排出，也就是全塔液泛，造成气液在塔内滞留，压降较大。

另外，图6为填料塔内气液两相流速情况。由图6可以看出，气液流速在填料拐角处变大，相应压降也变大。对气速和压降关系作图，如图7所示。从图7中以看出，压降随气速增加而增加；模拟压降和实验压降值在液泛以下，其值差别不大，变化趋势基本一致。在液泛时，模拟压降与实验值差别较大，实际液泛时压降应该参考模拟值。

图6 填料塔内气液两相流速情况

图7 压降与气速关系

5 四氟乙烯精馏塔体系数值模拟

由空气和水在常温常压情况下的单相和两相流体流动的实验和流体力学模拟分析可知，对不锈钢丝网规整填料流体力学建模可以较好地模拟实验过程。

另外，对四氟乙烯1#和2#精馏塔在中压低温生产过程中，不能看到其塔内的流体流动情况，也不能做实验，也未见公开材料报道，这样就造成对塔内流体的流动情况不清楚，不能对这两个塔的生产运行情况作出准确的评价，进而不能有效作出技术上进一步的升级改造。因此，通过流体力学建模和模拟，可以比较真实地再现塔内流体的流动情况，找出技术瓶颈，以便为生产提质增效、节能降耗提出可行性方案。

四氟乙烯1#和2#精馏塔气液两相物性参数见表5，其中气液两相在填料表面的接触角参考文献报道取值得到。

表5 四氟乙烯精馏塔体系物性参数

5.1 液相单相流动情况

1#和2#塔在中压低温情况下，喷淋密度分别为26 m3/(m2·h)和20 m3/(m2·h)，其液体在填料表面成膜情况见图8。从图8可以看出，在不同喷淋密度情况下，液相成膜厚度不一样，其膜厚度见表6。

图8 液相挂膜情况

1#塔在喷淋密度为26 m3/(m2·h)时，液相膜厚度均值为0.48 mm，并且在波纹竖直方向上，上表面成膜厚度薄，下表面成膜厚度大。另外，从图8可以看出，1#塔喷淋密度有点大，液膜厚度有点过大，其操作工况还不是最佳运行方案，可以通过增大塔径或提高填料效率或增加塔高甚至用更加高效的填料等来适当减小喷淋密度，使得该塔操作处于最佳操作范围。

2#塔喷淋密度为20 m3/(m2·h)时，液相膜厚度均值为0.36 mm，其成膜机理同1#塔。其次,2#塔喷淋密度适中，液膜厚度还有增加的空间，其操作工况应该是最佳运行方案。此外,1#和2#塔在喷淋密度分别为26 m3/(m2·h)和20 m3/(m2·h)时，液膜在填料表面连续成膜高度不一样，具体见表6。

表6 喷淋密度和液膜厚度关系

从表6可以看出，1#和2#塔在喷淋密度分别为26 m3/(m2·h)和20 m3/(m2·h)时，液体在丝网规整填料表面挂膜高度分别为51.5 mm和72.1 mm，由此建议填料盘高分别为52 mm和72 mm。

由前述可知，1#和2#塔液相挂膜高度接近空气-水体系喷淋密度在28.2～36 m3/(m2·h)时的值。从不同喷淋密度时液体分布情况的实验结果可知，塔径在Φ 400 mm左右时，1#和2#塔填料分段高度建议在 1～2 m。

5.2 2#塔气液两相流体动力学模拟分析

通过2#塔气液两相流体动力学模拟分析，对喷淋密度分别为20 m3/(m2·h)和33 m3/(m2·h),气速分别为0.135 m/s和0.223 m/s这两个操作工况进行比较，模拟结果见图9。

图9 液相挂膜和气液两相流速关系

从图9可以看出，在不同喷淋密度和气速情况下，塔内流体流动不同，在喷淋密度为20 m3/(m2·h)时，塔内流体流动和压降在正常范围。当喷淋密度为33 m3/(m2·h)时，塔的操作基本在液泛状态。这时应该减小喷淋密度至接近20 m3/(m2·h)。塔内不同喷淋密度和气速与压降的关系见表7，其中液泛时的模拟压降远远大于实际操作值，这主要是由于实际操作是局部液泛状态，而模拟液泛是全塔液泛状态，从而压降过大。

表7 喷淋密度和气速、压降关系

6 结论

1)通过丝网规整填料塔流体力学实验和Fluent 6.3数值模拟，建立了规整填料CY的二维模型，对气液两相流体在塔内的压降和液泛也进行了模拟分析，误差小于5%。

2)运用流体力学软件清楚地对空气-水体系填料塔内液体挂膜进行了数据化和透明化，包括液膜厚度和液膜长度以及液膜发展情况，进而对不同喷淋密度的填料盘高度给出了最佳值。

3)利用丝网规整填料CY数值模型分别对四氟乙烯体系1#和2#精馏塔进行了流体力学分析，对 1#和2#精馏塔液相成膜厚度和高度以及填料分段高度都给出了参考值。

4)对四氟乙烯体系2#精馏塔进行了两种工况的流体力学模拟分析，给出了正常操作和泛点操作的喷淋密度和压降值，为精馏塔的改进提供了可靠的数据。

[1] 陈江波, 刘春江. 计算流体力学方法模拟规整填料塔内流体流动行为的研究进展[J].石油化工, 2007, 36(9): 961- 965.

[2] 尹晔东, 王运东, 费维扬. 计算流体力学(CFD)在化学工程中的应用[J]. 石化技术, 2000, 7(3): 166-169.

[3] Raynal L, Royon L A. A multi-scale approach for CFD calculations of gas-liquid flow within large size columns equipped with structured packing[J]. Chem. Eng. Sci., 2007, 62: 7196-7204.

[4] Raynal L, Boyer C, Ballaguet J P. Liquid holdup and pressure drop determinationin structured packing with CFD simulations[J]. Can. J. Chem., 2004, 82: 871-879.

[5] 王福军. 计算流体动力学分析-CFD软件原理与应用[M]. 北京: 清华大学出版社, 2010.

[6] 丁欣硕, 焦楠. FLUENT14.5流体仿真计算从入门到精通[M]. 北京: 清华大学出版社, 2014.

[7] 李鹏飞, 徐敏义, 王飞飞. 精通CFD工程仿真与案例实战[M]. 北京: 人民邮电出版社, 2014.

[8] Zhao L, Cerro R L. Experimental characterization of viscous film flows over complex surfaces[J]. Int. J. Multiphase Flow, 1992, 18: 495-516.

[9] 王树楹. 现代填料塔技术指南[M]. 北京: 中国石化出版社, 1998.

[10] Fernandes J, Simos P C, Mota J P B, et al. Application of CFD in the study of supercritical fluid extraction with structured packing: dry pressure drop calculations[J]. J. Supercrit. Fluids, 2008, 47(1): 17-24.

[11] 刘春江, 成洁, 袁希钢. 规整填料内鼓泡流动的二维CFD模拟[J]. 中国科技论文在线, 2008, 3(12): 933-940.

[12] 余国琮, 袁希钢. 化工计算传质学导论[M]. 天津: 天津大学出版社, 2011.

Numerical Simulation of Fluid Dynamics of Tetrafluoroethylene Distillation Tower

Li Qingyuan, Zhao Hongtao

(Shanghai 3F New Materials Co., Ltd., Shanghai 200241, China)

By means of the experiment on the flow behavior of air-water system in a packed tower taking the fluid dynamics research, and the method of computational fluid dynamics, a numerical model of gauze structured packing CY was established. It was studied for the relationship between liquid film thickness and film length in different spray density and gas velocity, pressure drop. Thus, the similar numerical analysis for those packing towers in purification of tetrafluoroethylene was made,which makes it possible to provide strong support for improving process technology and equipment to the industrial production scale device and setting data visualization.

fluid dynamics; CY wire packing; TFE; packing tower

李清元(1978—)，男，硕士，工程师，从事化工过程开发、工程技术等方面研究。