浮式液化天然气生产储卸装置(LNG floating production storage and offloading unit, 简称为FLNG)是一种新型海上天然气田开发的浮式生产装置, 具有天然气开采、装卸、预处理、液化和储存等功能。同时, FLNG打破了长距离管线与陆上生产系统的经济和成本性的限制, 适用于开发中小型边际气田。目前, 关于FLNG塔器设备的有关蒸馏或气体处理过程的设计和运行的工业经验较为缺乏, Bussy等^[1]进行了填料高度裕量的理论研究, 该研究验证了确定海上晃动条件下塔器高径比L/D临界值的必要性。Cullinane等^[2]指出晃动频率越快, 液体不均匀分布越严重, 在理论上推导了倾角与高径比对填料塔内干区面积及累计干区体积的影响, 而对晃荡条件下填料塔高径比临界值的分析较为缺乏。Iliuta等^[3]建立了填料床的三维瞬态双流体模型, 考虑了填料床内径向均匀和不均匀孔隙率分布, 比较了垂直、倾斜、对称晃荡和不对称晃荡条件下填料床内气液流动状况。Lappalaine等^[4]指出润湿面积分数可对各动量交换项进行加权。Fourati等^[5]采用相同方法对各动量交换项进行加权, 检测润湿面积分数是否影响规整填料塔内径向液体分布情况。Tsai等^[6]提出一种新的有效接触面积预测模型。Sebastia等^[7-9]研究了规整填料塔在小尺度(定位于气液相间界面)上的物理传质性能, 采用Higbie渗透理论描述物理传质, 得到了最优液体载荷, 实现液体不均匀分布的最小化及传质速率最大化。Pham等^[10]建立了胺吸收塔在宏观尺度上的CFD模型。Pham等^[11]对横摇条件下填料塔内的脱酸效率进行了预测, 对比了静态垂直和横摇条件下气相中CO₂物质的量分数沿塔高分布情况。由于Pham等仅针对一种规整填料进行研究, 对于其他类型规整填料, 需要在模型中对Ergun系数的4个修正因子进行相应调整, 为此Kim等^[12]提出了一种计算方法, 该方法须已知规整填料内持液量和压降的实验数据, 通过经验模型估算经验系数。Liu等^[13-17]提出气液相弥散力模型, 但仅有Lappalainen等进行了实验验证。这些研究者认为宏观速度与体积平均速度的偏差导致液体扩散现象发生, 由此建立了气液相弥散力表达式。鉴于目前海况对FLNG气液逆流塔内两相流动与传质特性研究不足, 笔者以应用于FLNG的气液逆流填料塔为研究对象, 采用液体不均匀分布指标对全塔内液体分布状态进行定量分析, 结合动网格数学模型, 分析倾斜角度、晃荡角度、晃荡周期对塔内各流体力学参数的定量影响, 更直观表述流场参数在横截面上分布的均匀性, 进而定量评价海况条件对全塔内的气液流场及传质性能的影响。

1 多相流逆流模型建立 1.1 几何模型

规整填料塔的三维几何结构如图 1所示, 包括顶部液体分布器、中部填料层、底部气体分配器3部分, 塔高2.295 m, 塔径0.1 m, 规整填料塔的前视图如图 2所示。其中液体分布器安装在塔顶以下约0.05 m位置处, 液体分布器上均匀分布12个滴流孔, 位于填料层上方0.006 m的位置上, 直径为0.012 m, 塔内液体通过液体分配器均匀注入, 保证液体初始均匀分布状态, 同时净化气从塔顶输出。填料结构和表征体元如图 3(图中, u为孔隙内流体速度，m/s; V为孔隙体积，m³; 下标G和L分别代表气相和液相)所示。由图 3可知, 规整填料具有几何结构规则相似性, 在填料塔内可堆砌成大孔隙的规则通道, 宏观观测时可将规则通道假设为均匀的多孔介质流道。多孔介质中某点的表征体元相当于连续介质中的质点。塔器底部空间高度为0.04 m, 气体由塔底经气体分布器均匀注入。

1.2 气液在逆流全塔板内流动数学模型

为了对气液逆流全塔内传质规律进行准确预测, 以基于液体弥散理论的填料塔内液体扩散模型为基础, 加入组分输运模型。除了闭合模型中的多孔介质阻力、气液相间拖曳力和气液相弥散力作为动量源项加入, 还需要将传质速率和化学反应速率作为源项引入, 并通过自定义函数UDF进行调用。

采用Euler-Euler模型作为多相流模型, 在孔隙尺度下采用体积平均方法对Navier-Stokes方程进行推导, 得到宏观尺度下的质量守恒方程和动量守恒方程, 模型假设包括:①各相均为牛顿流体、不可压缩相; ②气液相间无传热, 塔与环境壁面之间不考虑能量交换, 可视为绝热过程; ③孔隙率视为各相同性。针对MEA-CO₂系统, 模型假设包括:①气相中仅有CH₄和CO₂两种组分, 液相中MEA和其他组分无蒸发现象; ②在CO₂吸收过程中仅考虑一个化学反应; ③仅CO₂发生气相和液相间传质过程。

1.2.1 流动基本方程

连续性方程:

$ \frac{\partial }{{\partial t}}\left( {{\alpha _k}{\rho _k}} \right) + \nabla \cdot \left( {{\alpha _k}{\rho _k}{\mathit{\boldsymbol{u}}_k}} \right) = 0. $

(1)

动量方程:

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\alpha _{\rm{G}}}{\rho _{\rm{G}}}{\mathit{\boldsymbol{u}}_{\rm{G}}}} \right) + \nabla \cdot \left( {{\alpha _{\rm{G}}}{\rho _{\rm{G}}}{\mathit{\boldsymbol{u}}_{\rm{G}}}{\mathit{\boldsymbol{u}}_{\rm{G}}}} \right) = - {\alpha _{\rm{G}}}\nabla p + \\ \nabla \cdot \left( {{\alpha _{\rm{G}}}{\mathit{\boldsymbol{\tau }}_{\rm{G}}}} \right) + {\alpha _{\rm{G}}}{\rho _{\rm{G}}}\mathit{\boldsymbol{g}} - {\mathit{\boldsymbol{F}}_{{\rm{GS}}}} + {\mathit{\boldsymbol{F}}_{{\rm{IG}}}} + {\mathit{\boldsymbol{F}}_{{\rm{D,G}}}}, \end{array} $

(2)

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\alpha _{\rm{L}}}{\rho _{\rm{L}}}{\mathit{\boldsymbol{u}}_{\rm{L}}}} \right) + \nabla \cdot \left( {{\alpha _{\rm{L}}}{\rho _{\rm{L}}}{\mathit{\boldsymbol{u}}_{\rm{L}}}{\mathit{\boldsymbol{u}}_{\rm{L}}}} \right) = - {\alpha _{\rm{L}}}\nabla p + \\ \nabla \cdot \left( {{\alpha _{\rm{L}}}{\mathit{\boldsymbol{\tau }}_{\rm{L}}}} \right) + {\alpha _{\rm{L}}}{\rho _{\rm{L}}}\mathit{\boldsymbol{g}} - {\mathit{\boldsymbol{F}}_{{\rm{LS}}}} + {\mathit{\boldsymbol{F}}_{{\rm{IL}}}} + {\mathit{\boldsymbol{F}}_{{\rm{D,L}}}}. \end{array} $

(3)

组分输运方程:

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\alpha _{\rm{G}}}{\rho _{{\rm{G,}}i}}} \right) = - {\mathit{\boldsymbol{u}}_{\rm{G}}} \cdot \nabla \left( {{\alpha _{\rm{G}}}{\rho _{{\rm{G,}}i}}} \right) + {D_{{\rm{G,}}i}}{\nabla ^2}\left( {{\alpha _{\rm{G}}}{\rho _{{\rm{G,}}i}}} \right) - \\ {r_{{\rm{GL}},i}},i = {\rm{C}}{{\rm{O}}_2}. \end{array} $

(4)

$ \frac{\partial }{{\partial t}}\left( {{\alpha _{\rm{L}}}{\rho _{{\rm{L,}}k}}} \right) = - {\mathit{\boldsymbol{u}}_{\rm{L}}} \cdot \nabla \left( {{\alpha _{\rm{L}}}{\rho _{{\rm{L,}}k}}} \right) + {D_{{\rm{L,}}k}}{\nabla ^2}\left( {{\alpha _{\rm{L}}}{\rho _{{\rm{L,}}k}}} \right) + \\ {r_{{\rm{GL}},i}} + {\alpha _{\rm{L}}}{R_k},\\ k{\rm{ = MEA,C}}{{\rm{O}}_{\rm{2}}}{\rm{,HORNH}}_3^ + {\rm{,HORNHCO}}{{\rm{O}}^ - }. $

(5)

其中

$ {\mathit{\boldsymbol{F}}_{{\rm{IG}}}} + {\mathit{\boldsymbol{F}}_{{\rm{IL}}}} = 0. $

式中, α为体积分数; ρ为流体密度, kg/m³; u为孔隙内流体速度, m/s; p为压力, Pa; τ_k为平均应力矢量; F_GS和F_LS分别为气固和液固间相互作用力, N/m³; F_IG和F_IL为气液相间拖曳力, N/m³; F_{D, G}和F_{D, L}分别为气相弥散力和液相弥散力; D_i为组分i的扩散系数, m²/s; r_GL为CO₂的气液相传质速率, kg/(m³·s); R_k为均匀化学反应速率, kg/(m³·s); g为自由落体加速度, m/s²。

1.2.2 多孔介质阻力

规整填料内可以在干缝(dry slit)和湿缝(wet slit)描述气相和液相动量平衡方程, 如图 4(图中，l、s、w分别为缝隙长度、厚度和间距)所示, 图 4(b)为部分润湿规整填料的气液流动形式, 图 4(c)为缝隙几何形态, 采用润湿面积分数f_e表征规整填料的润湿状态, 即当f_e=1时, 气液共存处填料表面完全润湿; f_e＜1代表填料表面部分润湿, 进而得到液固、气固、气液界面动量交换项的权重系数分别为f_e、1-f_e、f_e。

多孔介质阻力F_GS和F_LS可用Ergun方程进行定义, 表示为

$ {\mathit{\boldsymbol{F}}_{{\rm{LS}}}} = - {f_{\rm{e}}}\left( {\frac{{\alpha {E_1}}}{{36}} \cdot \frac{{a_{\rm{g}}^2}}{{\varepsilon \theta _{\rm{L}}^2}}{\mu _{\rm{L}}} + \frac{{b{E_2}}}{6} \cdot \frac{{{a_{\rm{g}}}}}{{{\theta _{\rm{L}}}}}{\rho _{\rm{L}}}\left\| {{\mathit{\boldsymbol{u}}_{\rm{L}}}} \right\|} \right){\mathit{\boldsymbol{u}}_{\rm{L}}} = - {K_{{\rm{LS}}}}{\mathit{\boldsymbol{u}}_{\rm{L}}}. $

(6)

$ \begin{array}{l} {\mathit{\boldsymbol{F}}_{{\rm{GS}}}} = - \left( {1 - {f_{\rm{e}}}} \right)\left( {\frac{{a{E_1}}}{{36}} \cdot \frac{{a_{\rm{g}}^2}}{\varepsilon }{\mu _{\rm{G}}} + \frac{{b{E_2}}}{6}{a_{\rm{g}}}{\rho _{\rm{G}}}{\theta _{\rm{G}}}\left\| {{\mathit{\boldsymbol{u}}_{\rm{G}}}} \right\|} \right){\mathit{\boldsymbol{u}}_{\rm{G}}} = \\ {K_{{\rm{GS}}}}{\mathit{\boldsymbol{u}}_{\rm{G}}}. \end{array} $

(7)

其中

$ {f_{\rm{e}}} = \frac{{{a_{\rm{e}}}}}{{{a_{\rm{s}}}}} = 1.34{\left[ {\frac{{{\rho _{\rm{L}}}}}{\sigma }{g^{1/3}}{{\left( {\frac{Q}{{{L_p}}}} \right)}^{4/3}}} \right]^{0.116}},Q = 5.1159Ah_{\rm{L}}^{2.854}. $

式中, K_LS和K_GS分别为液固及气固相间拖曳系数, kg/(m³·s); E₁和E₂为Ergun系数; a和b为Ergun系数的修正因子, 不同规整填料类型对应的a和b值也不同; ε为孔隙率量; θ为波纹角, (°); μ为相动力黏度, Pa/s; f_e为润湿面积分数, 采用Tsai提出的预测模型^[6]求得。对于规整填料Mellapak 250.X, E₁和E₂分别为160和0.16;对于规整填料Mellapak 500.X, E₁和E₂需要进行修正; a_e为有效接触面积, m²; a_s为填料的表面积, m²; ρ_L为输送液体密度, kg/m³; σ为溶液的表面张力, N/m; Q为输送液体流速, m³/s; L_p为湿周, m; A为塔器横截面积, m²; h_L为持液率。

1.2.3 气液相间拖曳力

气液相间拖曳力封闭方程^[18-20]为

$ \begin{array}{l} {\mathit{\boldsymbol{F}}_{{\rm{IG}}}} = - {f_{\rm{e}}}\left( {\frac{{c{E_1}}}{{36}} \cdot \frac{{a_{\rm{g}}^2}}{{{\varepsilon ^2}{\theta _{\rm{G}}}}}{\mu _{\rm{G}}} + \frac{{d{E_2}}}{6}\frac{{{a_{\rm{g}}}}}{\varepsilon }{\rho _{\rm{G}}}\left\| {{\mathit{\boldsymbol{u}}_{\rm{G}}} - {\mathit{\boldsymbol{u}}_{\rm{L}}}} \right\|} \right) \times \\ \left( {{\mathit{\boldsymbol{u}}_{\rm{G}}} - {\mathit{\boldsymbol{u}}_{\rm{L}}}} \right) = - {K_{{\rm{IG}}}}\left( {{\mathit{\boldsymbol{u}}_{\rm{G}}} - {\mathit{\boldsymbol{u}}_{\rm{L}}}} \right). \end{array} $

(8)

式中, K_IG为气液相间拖曳系数, kg/(m³·s); c和d为Ergun系数的修正因子, 不同规整填料类型对应的c和d值也不同。

1.2.4 气液相弥散力

为了实现对气液逆流填料塔内液体分布的准确预测, 引入液体弥散机制, 需要对液体弥散机制进行合理建模。经过分析, 描述规整填料内液体扩散时可忽略毛细管弥散作用, 以机械弥散作用为主。采用Lappalainen模型^[21-22]描述气液相弥散项, 表示为

$ {\mathit{\boldsymbol{F}}_{{\rm{D,L}}}} = {\theta _{\rm{L}}}{\mathit{\boldsymbol{K}}_{{\rm{LS}}}}{\mathit{\boldsymbol{u}}_{{\rm{D,L}}}} - \varepsilon {\mathit{\boldsymbol{K}}_{{\rm{IG}}}}\left( {{\mathit{\boldsymbol{u}}_{{\rm{D,G}}}} - {\mathit{\boldsymbol{u}}_{{\rm{D,L}}}}} \right), $

(9)

$ {\mathit{\boldsymbol{F}}_{{\rm{D,G}}}} = {\theta _{\rm{G}}}{\mathit{\boldsymbol{K}}_{{\rm{LS}}}}{\mathit{\boldsymbol{u}}_{{\rm{D,G}}}} - \varepsilon {\mathit{\boldsymbol{K}}_{{\rm{IG}}}}\left( {{\mathit{\boldsymbol{u}}_{{\rm{D,G}}}} - {\mathit{\boldsymbol{u}}_{{\rm{D,L}}}}} \right), $

(10)

$ {\mathit{\boldsymbol{u}}_{{\rm{D,G}}}} = \left( { - S\left\| {{\mathit{\boldsymbol{u}}_{\rm{G}}}/{\alpha _{\rm{G}}}} \right\|/{\alpha _{\rm{G}}}} \right)\nabla {\alpha _{\rm{G}}}, $

(11)

$ {\mathit{\boldsymbol{u}}_{{\rm{D,L}}}} = \left( { - S\left\| {{\mathit{\boldsymbol{u}}_{\rm{L}}}} \right\|/{\alpha _{\rm{L}}}} \right)\nabla {\alpha _{\rm{L}}}. $

(12)

式中, u_{D, G}和u_{D, L}为气液相漂移速度, m/s; S为扩散系数, mm。

1.2.5 传质速率

一般情况下, 整体传质系数采用双膜理论进行描述, 表示为

$ \frac{1}{K} = \frac{1}{{E{k_1}}} + \frac{1}{{H{k_{\rm{g}}}}}. $

(13)

式中, k_l和k_g分别为液膜和气膜内传质系数, m/s; K为整体液体传质系数, m/s; E为化学反应的增强因子, cal/gmol; H为亨利常数, Pa·m³/mol。

在CO₂-MEA系统中, 整体传质阻力(1/K)主要来自液体侧阻力(1/k_l)。因此, 整体传质系数可以简化为

$ K \approx E{k_1}. $

(14)

增强因子E是MEA浓度函数, 如果MEA浓度在反应区域几乎不变时, 增强因子E可视为常数。实际上, 在该胺吸收塔内MEA浓度为2.6~3.0 kmol/m³。

传质速率r_GL可表示为液膜传质系数和液体浓度驱动力函数, 表示为

$ {r_{{\rm{GL}}}} = {{k'}_1}{a_{\rm{e}}}\left( {\rho _{{\rm{L,C}}{{\rm{O}}_2}}^*{\rho _{{\rm{L,C}}{{\rm{O}}_2}}}} \right). $

(15)

式中, ρ^*_{L, CO2}为气液界面处CO₂的平衡液体质量浓度, kg/m³; k′_l为将液膜阻力和化学反应增强因子综合考虑的参数, m/s, 其值取决于系统本身, 可以通过理论或经验关联式得到。

ρ^*_{L, CO2}通过亨利定律得到:

$ \rho _{{\rm{L,C}}{{\rm{O}}_2}}^* = {M_{{\rm{w,C}}{{\rm{O}}_2}}}\frac{{{p_{{\rm{C}}{{\rm{O}}_2}}}}}{{{H_{{\rm{C}}{{\rm{O}}_2} - {\rm{MEA}}}}}}. $

(16)

式中, M_{w, CO2}为CO₂的分子质量, kg/kmol; p_CO2为气相中CO₂分压, Pa; H_CO2-MEA为MEA溶液中CO₂的亨利常数, 可以通过N₂O-CO₂的气液平衡类比得到。

$ {H_{{\rm{C}}{{\rm{O}}_2} - {\rm{MEA}}}} = {H_{{{\rm{N}}_{\rm{2}}}{\rm{O - MEA}}}}\frac{{{H_{{\rm{C}}{{\rm{O}}_{\rm{2}}}{\rm{ - water}}}}}}{{{H_{{{\rm{N}}_{\rm{2}}}{\rm{O - water}}}}}}. $

(17)

CO₂和N₂O在纯水中的亨利常数是温度函数, 为得到物质的量分数为6%MEA(或质量分数为31.7%的MEA)条件下的H_CO2-MEA, 需要借助PENTTILÄ^[8]的实验数据, 如图 5所示。操作温度为304 K, 则H_CO2-MEA=3 690 Pa·m³/mol。

1.2.6 化学反应速率

对于CO₂-MEA系统, 化学反应速率可以用二阶动力学方程^[10]表示为

$ {\rm{C}}{{\rm{O}}_{\rm{2}}}\left( {\rm{L}} \right){\rm{ + 2OHORN}}{{\rm{H}}_{\rm{2}}} \leftrightarrow {\rm{HORNH}}_{\rm{3}}^ + {\rm{ + HORNHCO}}{{\rm{O}}^ - }{\rm{.}}\\ {R_{{\rm{C}}{{\rm{O}}_2}}} = \frac{{\rm{d}}}{{{\rm{d}}t}}{\rho _{{\rm{C}}{{\rm{O}}_2}}} = - {k_{\rm{c}}}{\rho _{{\rm{L,C}}{{\rm{O}}_2}}}{\rho _{{\rm{L,MEA}}}}{k_{\rm{c}}} = {k_0}{{\rm{e}}^{ - {E_{\rm{a}}}/\left( {RT} \right)}}. $

(18)

式中, k_c为化学反应速率常数, m³/(kmol·s), 取决于温度T, 可表示为Arrhenius方程形式。

液相中共有5种组分:MEA、H₂O、HORNH₃⁺、HORNHCOO^-和CO₂(L)。由于摩尔浓度(C_i=ρ_{L, i}/M_w)优于质量浓度(ρ_{L, i}), 因此将式(18)转化为

$ {R_{{\rm{C}}{{\rm{O}}_2}}} = \frac{{\rm{d}}}{{{\rm{d}}t}}{C_{{\rm{C}}{{\rm{O}}_2}}} = - {k_{\rm{c}}}{C_{{\rm{C}}{{\rm{O}}_2}}}{C_{{\rm{MEA}}}}. $

(19)

其中, 传质速率r_GL和化学反应速率R均通过UDF加入组分输运方程中。当MEA-CO₂系统进行反应时, 液体的物理性质(扩散系数、黏度、表面张力)均会改变, 这将对有效接触面积和传质系数带来一定影响, 但在简单反应模型中可忽略不计。

1.3 网格划分及边界条件

胺吸收塔的操作压力为1.013×10⁵ Pa, 操作温度为304 K。液体分布器为液体入口, 采用质量流量入口边界; 塔顶为气体出口, 采用质量流量边界; 塔底为气体入口和液体出口, 采用压力出口边界; 壁面为无滑移边界。

塔顶处的液体分布器包含12个滴流孔, 结构较为复杂, 网格质量要求较高, 因此采用六面体网格进行精细划分, 如图 6所示。中间填料层为多孔介质区域, 结构较为简单, 网格质量要求较低, 因此采用六面体网格进行粗糙划分, 横截面上的网格结构如图 7所示。将非共形网格技术应用于顶部液体分布器和中间填料层, 以及底部气体分布器和中间填料层间的界面, 实现不同网格区域边界上的网格节点的位置不一致的情况。非共形界面通过传递通量使得相邻网格区域更为容易地连接, 从而减少了界面的网格数量, 液体分布器和气体分布器与多孔区域界面的网格划分如图 8所示。同时, 在塔壁处的边界层内进行网格加密。多孔介质区域选取网格数量30 500、57 820、103 200进行网格独立性检验, 塔器顶部和底部持液率均略低于填料区域内持液率, 网格数量为57 820和103 200时的持液率变化趋势较为接近, 但网格数量为30 500时的持液率变化趋势差异较大, 塔器顶部持液率低至0.092, 底部持液率低至0.093, 而填料区域内持液率接近0.095。经网格独立性检验满足计算的精度和收敛要求。

1.4 模型参数及求解

模拟介质为MEA-CO₂系统, 塔顶注入质量分数为31.7%的MEA溶液, 塔底注入物质的量分数为85%的CH₄和物质的量分数为15%的CO₂混合气体; 规整填料采用M500.X, 比表面积为500 m²/m³, 孔隙率为0.91, Ergun系数E₁=160, E₂=0.16, Ergun系数修正因子分别为a=1.9, b=1.37, c=0.17, d=1.16, MEA溶液的表面张力0.062 N/m, 填料塔湿周为4.05 m, 扩散系数取7.4 mm。模拟考虑传质和化学反应, 液体传质系数为0.00 104 m/s, 指前因子为5.928×10⁴ m³/(kmol·s), 化学反应的活化能为1.788×10⁴ J/mol, 气相和液相扩散系数分别取2.8×10^-7和2.88×10^-9 m²/s。

采用基于压力的非稳态、隐式方法求解欧拉多相流模型, 对非线性偏微分方程进行离散时, 采用一阶迎风格式离散格式对动量方程和体积分数进行空间离散, 时间项采用隐格式; 压力-速度耦合方程的求解采用phase coupled SIMPLEC算法。闭合模型中的多孔介质阻力、气液相间拖曳力和气液相弥散力均作为动量源项加入, 传质速率和化学反应速率均需加入组分输运方程中, 以上闭合模型均通过自定义函数UDF进行调用。规整填料采用多孔介质模型进行简化, 研究全塔内气液流场和传质性能, 但由于塔内存在逆向气流流动的两相流, 流动过程中两相间相互作用是一个非常复杂的动态过程, 对计算机配置要求高, 计算时间长。非稳态时间步长取0.001 s, 收敛残差为10^-5。

2 模型计算结果验证

模型主要通过塔内持液率、单位长度湿压降、气相中CO₂物质的量分数进行实验验证。在液体载荷q_L为10 m³/(m²·h), 气体动能因子为0.21 Pa^0.5条件下, 规整填料选用Mellapak 500.X。

图 9为液体载荷q_L为24.4 m³/(m²·h)时的单位长度湿压降随气体动能因子的变化情况。由图 9可得, 随着气体动能因子增加, 塔内单位长度湿压降也随之增加。规整填料单位长度湿压降的模拟结果与实验数据吻合较好。

图 10为气体动能因子为0.21 Pa^0.5(即气体载荷q_G为38.5 kmol/(m²·h))时的平均持液率随液体载荷的变化情况。其中持液率h_L为整个填料区域内的平均持液率, 填料为规整填料Mellapak 500.X。由图 10可知, 模拟结果与实验结果偏差在合理范围内。

图 11为气相中CO₂物质的量分数的径向平均值与实验结果对比。由图 11可得, 所建模型在塔内气相中CO₂物质的量分数和实验结果吻合较好, 可以进行全塔内气液流场及传质规律分析。

3 海况边界条件的动网格数学模型

通过对模拟对象使用六自由度定义方法, 利用研究对象的力和时间计算物体重心平移和角度的运动。在惯性坐标系中, 研究对象重心平移运动的控制方程可以表示为

$ {{\mathit{\boldsymbol{\dot v}}}_{\rm{G}}} = \frac{1}{m}\sum {{\mathit{\boldsymbol{f}}_{\rm{G}}}} . $

(20)

式中, ${\mathit{\boldsymbol{\dot v}}_{\rm{G}}}$为中心位移运动加速度; f_G为重力。

模拟对象的角度运动, ${\mathit{\boldsymbol{\dot \omega }}_B}$可表示为

$ {{\mathit{\boldsymbol{\dot \omega }}}_B} = {I^{ - 1}}\left( {\sum {{\mathit{\boldsymbol{M}}_B} - {\mathit{\boldsymbol{\omega }}_B} \times L{\mathit{\boldsymbol{\omega }}_B}} } \right). $

(21)

式中, I为惯性矢量; M_B为时间矢量; ω_B为刚体角速度向量; L为晃荡位移, m。

采用弦等效波方法简化复杂的海况条件, Hoerner等验证了该简化方法的适用性:

$ {P_0} = {a_{\rm{w}}}\sin \left( {{\omega _{\rm{e}}}t + \varphi } \right). $

(22)

式中, P₀为研究对象在时间t的位置; a_w为波幅; ω_e为遭遇频率, Hz; φ为相位角, (°)。

对式(22)关于时间求导得到角速度^{[1, 22]}为

$ {\omega _I} = \frac{{{\theta _{\max }}\pi }}{{180}} \cdot \frac{{2\pi }}{T}\cos \left( {\frac{{2\pi t}}{T}} \right). $

(23)

式中, θ_max为晃荡幅度, (°); T为晃动周期, s; t为时刻值, s; 下标I代表x、y和z轴任一方向。

线速度可以表示为

$ {v_I} = L \times \frac{{2\pi }}{T}\cos \left( {\frac{{2\pi t}}{T}} \right). $

(24)

在一些实际问题中, 流体流动区域通常会因控制域边界随时间变化而造成形状改变和位置改变, 因此需要采用动网格技术对流体域进行建模。

采用动网格技术后, 共有两种速度, 一种是静止坐标系中的流体速度u_L和u_G, 另一种是新坐标系中的网格速度v_mesh, 表达形式如下:

$ {\mathit{\boldsymbol{v}}_{{\rm{mesh}}}} = \mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{r}}\mathit{\boldsymbol{.}} $

(25)

式中, ω为角速度, s^-1; r为位置向量。

倾斜条件下的单位方向向量(u_d)_tilt表示为

$ {\left( {{\mathit{\boldsymbol{u}}_{\rm{d}}}} \right)_{{\rm{tilt}}}} = \frac{\mathit{\boldsymbol{r}}}{{\left| \mathit{\boldsymbol{r}} \right|}} = \sin \theta \cdot \mathit{\boldsymbol{i}} + \cos \theta \cdot \mathit{\boldsymbol{j}}\mathit{\boldsymbol{.}} $

(26)

对称晃荡是沿着对称轴做旋转周期摆动运动, 可将三维晃动模型简化为二维晃动模型, 则横摇和纵摇条件下的单位方向矢量(u_d)_roll表示为

$ {\left( {{\mathit{\boldsymbol{u}}_{\rm{d}}}} \right)_{{\rm{roll}}}} = \frac{\mathit{\boldsymbol{r}}}{{\left| \mathit{\boldsymbol{r}} \right|}} = \sin \theta \cdot \mathit{\boldsymbol{i}} + \cos \theta \cdot \mathit{\boldsymbol{j,}} $

(27)

其中

$ \theta \mathit{\boldsymbol{ = }}{\theta _{\max }}\cos \left( {\frac{{2\pi t}}{T} + \frac{\pi }{2}} \right). $

式中, i为x轴方向向量; j为y轴方向向量。将两个单位方向向量应用于填料塔的晃动中心处, 由船体倾斜或晃动引起的动量源项表示为

$ {({\mathit{\boldsymbol{S}}_{\rm{G}}})_{{\rm{ship}}}} = {{S'}_{{\rm{Gx}}}} \cdot \mathit{\boldsymbol{i'}} + {{S'}_{{\rm{Gy}}}} \cdot \mathit{\boldsymbol{j'}} + {{S'}_{{\rm{Gz}}}} \cdot \mathit{\boldsymbol{k',}} $

(28)

$ {({\mathit{\boldsymbol{S}}_{\rm{L}}})_{{\rm{ship}}}} = {{S'}_{{\rm{Lx}}}} \cdot \mathit{\boldsymbol{i'}} + {{S'}_{{\rm{Ly}}}} \cdot \mathit{\boldsymbol{j'}} + {{S'}_{{\rm{Lz}}}} \cdot \mathit{\boldsymbol{k'}}\mathit{\boldsymbol{.}} $

(29)

式中, i′、j′和k′为新坐标系的方向向量。

以气相为例, 将(S_G)_ship用原来x-y坐标系表示如图 12所示。

$ {{S'}_{{\rm{Gx}}}} \cdot \mathit{\boldsymbol{i'}} = \left| {{{S'}_{{\rm{Gx}}}}} \right| \cdot \cos \theta \cdot \mathit{\boldsymbol{i}} - \left| {{{S'}_{{\rm{Gx}}}}} \right| \cdot \sin \theta \cdot \mathit{\boldsymbol{j}}, $

(30)

$ {{S'}_{{\rm{Gy}}}} \cdot \mathit{\boldsymbol{j'}} = \left| {{{S'}_{{\rm{Gy}}}}} \right| \cdot \sin \theta \cdot \mathit{\boldsymbol{i}} + \left| {{{S'}_{{\rm{Gy}}}}} \right| \cdot \cos \theta \cdot \mathit{\boldsymbol{j}}, $

(31)

$ {{S'}_{{\rm{Gz}}}} \cdot \mathit{\boldsymbol{k'}} = {S_{{\rm{Gz}}}} \cdot \mathit{\boldsymbol{k}}\mathit{\boldsymbol{.}} $

(32)

4 液体在塔内不均匀分布评价指标

填料塔内液体不均匀分布状态采用均匀系数R_UI作为评价指标。R_UI表示流域内某个特定参数在截面上的变化情况, 当R_UI越接近于1代表该参数在截面上分布越均匀, R_UI越接近于0代表该参数在截面上分布越不均匀。

$ {R_{{\rm{UI}}}} = 1 - \frac{{\sum\limits_{i = 1}^n {\left[ {\left( {{\varphi _i} - \bar \varphi } \right){A_i}} \right]} }}{{2\left| {\bar \varphi } \right|\sum\limits_{i = 1}^n {{A_i}} }}, $

(33)

其中

$ \bar \varphi = \sum\limits_{i = 1}^n {{\varphi _i}{A_i}} /\sum\limits_{i = 1}^n {{A_i}} . $

式中, φ_i为第i个面索引的参数值; φ为在该截面上参数的平均值; A_i为第i个面索引的截面面积, m²。

均匀系数R_UIα为填料塔内液体体积分数在截面上分布的评价指标, 计算公式为

$ {R_{{\rm{UI}}\alpha }} = 1 - \frac{{\sum\limits_{i = 1}^n {\left[ {\left( {\left| {{\alpha _{{\rm{L,}}i}} - {{\bar \alpha }_{\rm{L}}}} \right|} \right){A_i}} \right]} }}{{2\left| {{{\bar \alpha }_{\rm{L}}}} \right|\sum\limits_{i = 1}^n {{A_i}} }}, $

(34)

其中

$ {{\bar \alpha }_{\rm{L}}} = \sum\limits_{i = 1}^n {{\alpha _{{\rm{L,}}i}}{A_i}} /\sum\limits_{i = 1}^n {{A_i}} . $

均匀系数R_UIu为填料塔内液体速度在截面上分布的评价指标, 计算公式为

$ {R_{{\rm{UIu}}}} = 1 - \frac{{\sum\limits_{i = 1}^n {\left[ {\left( {\left| {{u_{{\rm{L,}}i}} - {{\bar u}_{\rm{L}}}} \right|} \right){A_i}} \right]} }}{{2\left| {{{\bar u}_{\rm{L}}}} \right|\sum\limits_{i = 1}^n {{A_i}} }}, $

(35)

其中

$ {{\bar u}_{\rm{L}}} = \sum\limits_{i = 1}^n {{u_{{\rm{L}},i}}{A_i}} /\sum\limits_{i = 1}^n {{A_i}} . $

5 海况条件对塔内气液流场影响 5.1 倾斜角度影响

当填料塔处于倾斜状态时, 填料内部的流体流动可假设为垂直, 导致部分填料内部无足够的液体通过, 而另一侧液体流量过大, 此时气液相无法充分均匀接触导致塔器效率降低。

图 13为不同倾斜角度下塔内液体体积分数的均匀系数R_UIα。由图 13可得, 倾斜角度为1°时均匀系数R_UIα的变化范围为0.8~ 0.9, 均匀系数随着填料高度降低而缓慢降低; 在倾斜角度为3°时均匀系数R_UIα的变化范围为0.5~0.8, 与静止状态时相比降低幅度约达30%, 均匀系数亦随着填料高度降低而缓慢降低, 但接近最底层填料时均匀系数快速降低; 在倾斜角度为6°时R_UIα集中于0.5左右, 与静止状态时相比降低幅度约达50%, 均匀系数随着填料高度降低而快速降低, 在塔器顶部前两层填料内均匀系数降低约0.3, 在塔器底部填料内均匀系数降低约0.12。因此, 倾斜角度3°以上均对塔内液体体积分数分布影响较大。

图 14为不同倾斜角度下液体速度大小的均匀系数。

由图 14可知, 倾斜角度为1°时均匀系数几乎未变, 倾斜角度为3°和6°时均匀系数变化范围为0.65~0.99, 与静止状态相比降低幅度达5%~10%;越接近塔底均匀系数降低趋势越明显, 在塔底最底层填料内均匀系数变化范围为0.65~0.75。倾斜角度3°以上时对塔内液体速度影响较大。

图 15为不同倾斜角度下气液相间接触面积的纵截面分布。由图 15可知:倾斜角度为1°时, 倾斜侧沿-x方向气液相间接触面积a_e较高, 塔壁处a_e约为400~ 410 m²·m^-3, 另一侧气液相间a_e较低, 塔壁处a_e约为290~300 m²·m^-3, 且填料高度越接近塔底, 填料内气液相间接触面积增加趋势越明显; 随着倾斜角度增加, 倾斜侧气液相间接触面积呈增加趋势; 在倾斜角度为3°时, 倾斜侧塔壁处a_e约为420~450 m²·m^-3; 在倾斜角度为6°时, 倾斜侧塔壁处a_e约为440~460 m²·m^-3。

5.2 晃荡角度影响

填料塔的晃动中心为塔底中心(x=0.0 m, y=0.0 m, z=0.0 m), 以晃荡周期6 s, 晃荡角度1°、3°、6°为例, 分析晃荡角度对全塔内气液流场参数的影响。图 16为晃荡角度1°, 晃荡周期6 s时均匀系数R_UIα随时间的变化。由图 16可知, 各纵截面的均匀系数变化很小, 填料高度z≥1.96 m时的均匀系数为0.9~1.0, 而z=2.205 m处的均匀系数约为0.776。图 17为不同晃荡角度下均匀系数R_UIα和R_UIu, 以填料高度z=1.225 m为例进行分析。由图 17(a)可知, 当填料塔进入晃荡的第一个周期内, 各晃荡角度条件下的均匀系数R_UIα略低于垂直状态时的均匀系数R_UIα, 但从第二个周期开始, 各晃荡角度条件下的均匀系数R_UIα均高于垂直状态的, 且晃荡角度越高, 均匀系数R_UIα增加越明显。由此可得在晃荡条件下, 填料塔内液体分布状态呈现前期不均匀化加剧, 一定时间后趋于均匀化的趋势, 这将有利于塔内气液相间接触面积增加, 在一定程度上将弥补由晃荡条件造成的气液相偏流现象。

由图 17(b)可知, 晃荡角度为1°时均匀系数R_UIu变化幅度非常小, 晃荡角度为3°时均匀系数R_UIu在前10 s均低于垂直状态的, 在10 s之后略高于垂直状态的。晃荡角度为6°时均匀系数R_UIu在前12 s均低于垂直状态的, 在12 s之后略高于垂直状态的。晃荡角度越高, 均匀系数R_UIu的变化趋势越明显。

5.3 晃荡周期影响

填料塔的晃动中心为塔底中心(x=0.0 m, y=0.0 m, z=0.0 m), 以晃荡角度3°, 晃荡周期为6和10 s为例, 分析晃荡周期对全塔内气液流场参数的影响。图 18为晃荡角度3°、晃荡周期10 s时均匀系数随时间的变化。均匀系数为0.9~1.0, 各纵截面的均匀系数变化很小。

图 19为不同晃荡周期下均匀系数随时间的变化。以填料高度z=1.225 m进行分析。晃荡周期为6 s时, 在第一个周期内均匀系数均低于垂直状态时的均匀系数, 之后其均高于垂直状态的; 而晃荡周期为10 s时, 第一个周期内均匀系数均低于垂直状态的, 在12之后均匀系数与垂直状态时的均匀系数相差不大。

图 20为不同晃荡条件下净化气中CO₂物质的量分数的变化情况。由图 20可得, 净化气中CO₂物质的量分数y_CO2以塔器垂直静止状态时的数值为轴上下波动。在同一横摇周期下, 净化气中CO₂物质的量分数的波动周期与填料塔的横摇周期近似相同; 横摇角度为3°和6°时, y_CO2的波动幅度分别为0.032和0.071, 即填料塔的横摇角度越大, 净化气中CO₂物质的量分数的波动幅度越大。在同一横摇角度下, 净化气中CO₂物质的量分数的波动周期与填料塔的横摇周期近似相同, 横摇周期分别为6 s和10 s时, y_CO2的波动幅度分别为0.038和0.071, 即填料塔的横摇周期越长, 净化气中CO₂物质的量分数的波动幅度越大。相比于6 s时, 波动幅度增加了97%。

6 结论

(1) 塔器在倾斜状态下, 随着倾斜角度增加, 塔内各纵截面各参数不均匀分布趋势逐渐增强; 倾斜角度3°以上对塔内液体速度和液体体积分数分布影响较大。在倾斜角度为3°时, 倾斜侧塔壁处a_e约为420~450 m²·m^-3; 倾斜角度大于6°时均匀系数集中于0.5左右, 其范围为0.9~0.99。

(2) 在晃荡条件下填料塔内液体分布状态呈现前期不均匀化加剧, 一定时间后趋于均匀化的趋势。

(3) 在晃荡条件下净化气中CO₂物质的量分数以塔器垂直静止状态时的数值为轴上下波动。净化气中CO₂物质的量分数的波动周期与填料塔的横摇周期近似相同。在同一横摇周期下, 横摇角度为3°时, CO₂物质的量分数的波动幅度为0.032;横摇角度为6°时, CO₂物质的量分数的波动幅度为0.071, 即填料塔的横摇角度越大, CO₂物质的量分数的波动幅度越大; 在同一横摇角度下, 填料塔的横摇周期越长, CO₂物质的量分数的波动幅度越大。