计及气泡诱导与剪切湍流的气泡破碎、湍流相间扩散及传质模型

作者：施炜斌龙姗姗杨晓钢蔡心悦来源：《化工学报》日期：2022-08-31人气：2025

对于气液鼓泡流而言，包括相间作用力、相间传热、传质在内的相间相互作用过程一方面受到气泡动力学行为的影响，另一方面又受到液相湍流涡旋的作用。这两个方面的影响和作用是相辅相成的，不应割裂开来看待。

在数值模拟研究工作中，已有大量工作聚焦于气泡动力学行为对相间相互作用的重要影响，并有部分数学模型能够用于描述涡旋-气泡碰撞导致的聚并破碎过程[1-8]。但是，这些数学模型大多基于柯尔莫哥洛夫-5/3定律来近似获得存在于惯性子区的湍流脉动和耗散。该近似方法对剪切湍流占主导的流动较为有效，但对由气泡诱导湍流（BIT）占据大量湍流结构的气液两相流动而言，该近似或存在不甚完善之处。虽然也有部分模型如Sato模型、Rzehak and Krepper模型等试图从湍流模型的角度描述气泡对液相湍流的整体作用从而影响聚并破碎过程，但问题的根本还是在于缺乏准确的数学模型对BIT的能谱函数进行构建。因此，也难以精确刻画BIT作用下的湍流涡旋与气泡的碰撞过程。

上升气泡诱导形成的尾涡湍流对气液两相流的数值模拟工作也具有极为重要的作用。特别是如鼓泡塔这类主要由气泡上升引起的湍流流动，气泡的尾涡从产生、脱落、传递、耗散均会对气液相间相互作用产生重大影响。Risso等[9]用空间中随机分布的小球阵列诱导形成的湍流涡旋尾迹衰减类比气泡在均匀湍流中形成的尾涡，并发现即便小球体积分数上升到13%，在同样能量输入和积分尺度的情况下，这种尾迹的衰减比标准湍流流动中的湍流涡旋衰减更快。Mercado等[10]使用自主研发的相敏恒温热线（CTA）探针，并辅以激光多普勒测速仪（LDA）同步校准，在气泡群尾涡中测得了以湍流涡旋波数-3次幂斜率为特性的湍流动能能谱。Riboux等[11]使用粒子成像测速仪（PIV）测量气泡群尾涡中的湍流能谱，并肯定其波数的指数幂非常接近-3。Mendez-Diaz等[12]同时使用热膜探针和光纤探针对气泡尾迹的湍流特性进行测量，发现气含率从0.5%逐渐上升到6%的实验条件下，在脉动频率高达3000 Hz范围内湍流能谱函数的斜率始终保持在波数的-3次幂。Prakash等[13]同样使用CTA探针测量了液相的速度脉动并再次确认κ-3的特性不只适用于气泡诱导湍流而且能够普适性地描述湍动鼓泡流动。在这些实验证据的基础上，他们提出了由于气泡存在导致的能量生成和黏性耗散的能量平衡关系，与Lance等[14]在早期工作中基于热线风速仪和LDA实验得到指数幂为-8/3的结果并由此推断指数幂为-3的结论相符。还有，Roghair等[15]和Riboux等[16]使用完整解析自由上升变形气泡的直接数值模拟方法，同样清晰地指出气泡诱导湍流的湍流动能能谱具有波数的-3次指数幂的特性。基于上述研究可以认为：①气泡诱导湍流和剪切湍流在湍流动能能谱惯性子区的范围内基本处于长期并存的状态，二者共同影响气液相间相互作用过程；②尽管气泡诱导尾涡湍流能量级串过程尚未被完全理解和验证，其与单相均匀湍流二者在行为特征上的区别已经能够在湍流能谱图像上清晰识别。因此，如果能对上述实验现象加以总结归纳，并建立物理意义明确清晰的理论模型，能够为深入认识BIT作用下聚并破碎过程提供良好基础。

对于气泡聚并破碎行为发生的频率和概率都较低的缓慢流动，除了主要流动方向上的曳力作用，气泡诱导湍流引起的湍流扩散作用也不容忽视。Sommerfeld等[17]在欧拉-拉格朗日模拟中全面考察了各种相间作用力对相间作用过程的影响程度，并明确指出在数值模拟中准确描述气泡摇摆上升的行为特点也非常关键。在欧拉-拉格朗日框架下，Basset力对气泡摇摆行为的影响相对明确，但Basset力是一个瞬时的力，其作用效果会因为欧拉-欧拉模拟中的短时平均而严重弱化。在欧拉-欧拉模拟中，湍流涡旋作用在气泡上并使其左右摇摆通常可以通过湍流扩散力模型进行描述。Laviéville等[18]指出湍流扩散会导致气泡在径向上从较高相含率的区域向较低相含率的区域迁移。de Bertodano[19]提出了早期的湍流扩散力模型，模型参数取决于比例系数与气含率梯度的乘积，但是该比例系数对1 ~ 5 mm的气泡的取值由0.1 ~ 0.5变化。Drew[20]通过考虑湍流涡时间尺度和气泡弛豫时间的比率，提出了较为适用于鼓泡流动的湍流扩散力模型。Lucas等[21]指出对于在管道中的鼓泡流动，其气泡相含率径向分布的预测有必要额外考虑气泡变形诱导的湍流导致的扩散作用。但是，这些通过系综平均获得的湍流扩散力模型只是部分地反映了BIT在径向上的作用效果，还有相当一部分的作用在平均化的过程中被平滑而消失。由此可见，气泡诱导的湍流涡旋的湍流扩散作用确实是准确描述气泡在径向上的摇摆行为的关键，但如何充分反映其作用机制需要进一步的深入研究。

气液相间传质同样是气泡-涡旋相互作用的过程。研究人员对气泡尺寸、气液滑移速度、湍流耗散率、扩散速率等影响气液相间传质因素的作用进行了大量的研究，并分别发展了对流传质模型、停滞膜模型、双膜理论模型、溶质渗透模型、表面更新模型等理论模型[22-23]。随着研究的逐步深入，部分研究工作中逐渐反映出了湍流涡旋对气液相间传质过程的重要影响，并建立了Eddy Cell模型[24]、Eddy Contact模型[5, 25]等能够直接反映湍流涡旋对相间传质作用的数学模型。尤其是Eddy Contact模型通过考虑湍流能谱含能区、惯性子区、耗散区的全能谱函数，从理论上解释了从大到小的所有湍流涡旋对气液相间传质的贡献。与聚并破碎模型的问题相似，上述相间传质模型也仅考虑了剪切湍流的涡旋。但是气泡诱导湍流的湍流涡旋同样会参与相间传质，用于描述气液相间传质的模型也需要对气泡诱导湍流和剪切湍流的共同作用加以考虑。

因此，为深入认识气泡诱导湍流与剪切湍流的共同作用机理及其对相间传质过程的影响，本研究采用理论建模、大涡模拟（LES）、雷诺平均模拟、群平衡模型（PBM）等研究方法，揭示BIT在鼓泡塔气液两相流数值模拟研究工作中的重要作用。

1 数学模型

1.1 双流体模型

本研究中气液两相流采用欧拉-欧拉双流体模型，对连续相和离散相都有一组质量和动量守恒方程。LES模型可通过对动量方程的空间滤波获得。控制方程如下：

$\frac{\partial (ρ_{k} α_{k})}{\partial t} + \nabla ∙ (ρ_{k} α_{k} u_{k}) = 0$ (1) $\begin{array}{l} \frac{\partial (ρ_{k} α_{k} u_{k})}{\partial t} + \nabla ∙ (ρ_{k} α_{k} u_{k} u_{k}) = \\ - α_{k} \nabla p + \nabla ∙ τ_{k} + ρ_{k} α_{k} g + M_{F, k} \end{array}$ (2)

式中，ρ、α、 u 、 τ 、 MF 分别代表密度、相含率、速度矢量、黏性剪应力张量、相间动量交换项；k为气相g或液相l，二者相含率之和为1。

1.2 湍流模型

本研究采用Reynolds应力方程(RSM)和LES两种湍流模型对气液鼓泡流动展开研究。当采用RSM湍流模型时，BIT的各向异性特点能够通过不同方向上的源项有所反映[26-29]。本研究采用的RSM湍流模型见表1。

表1 RSM湍流模型方程

Table 1 Reynolds stress equation model

模型	方程
$ρ \bar{u_{i}^{'} u_{j}^{'}}$ 输运方程	$\frac{\partial (α_{l} ρ_{l} \bar{u_{i}^{'} u_{j}^{'}})}{\partial t} + \frac{\partial (α_{l} ρ_{l} u_{k} \bar{u_{i}^{'} u_{j}^{'}})}{\partial x_{k}} = \frac{\partial}{\partial x_{k}} [α_{l} (μ_{l} + \frac{μ_{t}}{σ_{k}}) \frac{\partial \bar{u_{i}^{'} u_{j}^{'}}}{\partial x_{k}}] + α_{l} P_{i j} + α_{l} ϕ_{i j} - \frac{2}{3} δ_{i j} α_{l} ρ_{l} ε + α_{l} S_{i j}^{B I T}$ (3)
湍动能产生项	$P_{i j} = - ρ_{l} (\bar{u_{i}^{'} u_{k}^{'}} \frac{\partial u_{j}}{\partial x_{k}} + \bar{u_{j}^{'} u_{k}^{'}} \frac{\partial u_{i}}{\partial x_{k}})$ (4)
压力-应变项	$ϕ_{i j} = ϕ_{i j, 1} + ϕ_{i j, 2} + ϕ_{i j, 1}^{W} + ϕ_{i j, 2}^{W} = - C_{1} ρ_{l} \frac{ε}{k} (\bar{u_{i}^{'} u_{j}^{'}} - \frac{2}{3} k δ_{i j}) - C_{2} ρ_{l} \frac{ε}{k} (P_{i j} - \frac{1}{3} t r (P) δ_{i j}) - C_{1}^{W} ρ_{l} \frac{ε}{k} (\bar{u_{k}^{'} u_{m}^{'}} n_{k} n_{m} δ_{i j} - \frac{3}{2} \bar{u_{k}^{'} u_{i}^{'}} n_{k} n_{j} - \frac{3}{2} \bar{u_{k}^{'} u_{j}^{'}} n_{k} n_{i}) {(\frac{k^{3 / 2}}{ε} \frac{1}{C_{l} y_{W}})}^{2} - C_{2}^{W} (ϕ_{k m, 2} n_{k} n_{m} δ_{i j} - \frac{3}{2} ϕ_{i k, 2} n_{k} n_{j} - \frac{3}{2} ϕ_{j k, 2} n_{k} n_{i}) {(\frac{k^{3 / 2}}{ε} \frac{1}{C_{l} y_{W}})}^{2}$ (5)
湍流耗散率输运方程	$\frac{\partial (α_{l} ρ_{l} ε)}{\partial t} + \frac{\partial}{\partial x_{i}} (α_{l} ρ_{l} ε u_{i}) = \frac{\partial}{\partial x_{j}} [α_{l} (μ_{l} + \frac{μ_{t}}{σ_{ε}}) \frac{\partial ε}{\partial x_{j}}] + α_{l} ρ_{l} \frac{ε}{k} (C_{1 ε} \bar{{u^{'}}_{i} {u^{'}}_{j}} \frac{\partial u_{i}}{\partial x_{k}} - C_{2 ε} ε) + α_{l} S_{ε}^{B I T}$ (6)
湍动能	$k = \frac{1}{2} (\sum_{i = 1,2, 3} u_{i}^{'} u_{j}^{'})$ (7)

新窗口打开| 下载CSV

在RSM模型中，湍流耗散率是一维标量，其源项表达式见下[30]：

$S_{ε}^{B I T} = C_{ε, B I T} \frac{S_{k}^{B I T}}{τ}$ (8)

此处的湍动能源项相应的也为标量，其表达式如下[31]：

$S_{k}^{B I T} = C_{k, B I T} F_{l}^{d r a g} (u_{g} - u_{l})$ (9)

Reynolds应力输运方程中，由于BIT的作用导致的源项 $S_{i j}^{B I T}$ 在主方向和其他方向上有所不同，其表达式为[26]：

$S_{i j}^{B I T} = [\begin{matrix} 1.0 & 0.0 & 0.0 \\ 0.0 & 0.5 & 0.0 \\ 0.0 & 0.0 & 0.5 \end{matrix}] S_{k}^{B I T}$ (10)

对于LES湍流模型，应力项可写为：

$τ_{k} = - μ_{e f f} [\nabla u_{k} + {(\nabla u_{k})}^{T} - \frac{2}{3} I (\nabla ∙ u_{k})]$ (11)

式中， $μ_{e f f}$ 代表连续相的有效黏度，由分子动力黏度 $μ_{l}$ 、湍流黏度 $μ_{t}$ 和描述BIT贡献项 $μ_{B I T}$ 共同构成。

$μ_{e f f} = μ_{l} + μ_{t} + μ_{B I T}$ (12)

其中，BIT的影响其中一部分可以通过Sato模型有所反映。

$μ_{B I T} = ρ_{l} C_{μ, B I T} α_{g} d_{b} |u_{g} - u_{l}|$ (13)

更为重要的是，考虑到气泡由周围湍流涡旋的剪切作用导致其动态响应，Long等[32]详细考察了LES模拟中气泡对液相湍流涡黏度的贡献，其修正的液相亚格子尺度湍流涡黏度模型如下：

$μ_{t, l} = ρ_{l} {(C_{s} ∆)}^{2} |S| (1 + C_{b u b b l e} \bar{α_{g}} \frac{∆}{d_{b}} {(\frac{1}{1 + S t})}^{\frac{3}{2}})$ (14)

1.3 相间作用力模型

本研究中选择曳力、升力、虚拟质量力为基本的相间作用力，表达式见表2。其中，曳力系数为Grace模型[33]、升力系数为Tomiyama模型[34]，虚拟质量力系数为0.5。

表2 相间作用力模型方程

Table 2 Interphase momentum exchange models

作用力	方程
曳力	$M_{D} = \frac{3}{4} \frac{C_{D}}{d_{b}} ρ_{l} α_{g} \|u_{g} - u_{l}\| (u_{g} - u_{l})$ (15)
升力	$M_{l i f t} = C_{l} ρ_{l} α_{g} (u_{g} - u_{l}) \times (\nabla \times u_{l})$ (16)
虚拟质量力	$M_{V M} = C_{V M} ρ_{l} α_{g} ({\frac{d u_{l}}{d t}\|}_{l} - {\frac{d u_{g}}{d t}\|}_{g})$ (17)

新窗口打开| 下载CSV

1.4 群平衡模型

群平衡模型是描述气泡尺寸分布的方法之一。气泡尺寸分布决定了气液相间接触面积，对相间传热、传质、化学反应等过程有重要影响。气液鼓泡流动的群平衡方程主要考虑气泡的聚并、破碎。气泡根据直径被分成不同组，di 为第i组气泡直径。因此，可将群平衡方程写为：

$\frac{\partial n_{i}}{\partial t} + \nabla ∙ (u_{i} n_{i}) = S_{i}$ (18)

式中，ni 为第i组气泡数密度; ui 为速度矢量;Si 为代表因气泡破碎聚并导致气泡生成和消亡的源项。Si 的表达式可写为：

$\begin{matrix} S_{i} = B_{c o a l e s c e n c e, i} - D_{c o a l e s c e n c e, i} + B_{b r e a k u p, i} - D_{b r e a k u p, i} = \\ \begin{array}{l} \sum_{d_{i} = d_{m i n}}^{\frac{d_{i}}{2}} Ω_{C} (d_{j} ∶ d_{i} - d_{j}) - \sum_{d_{i}}^{d_{m a x} - d_{i}} Ω_{C} (d_{j} ∶ d_{i}) + \\ \sum_{d_{j} = d_{i}}^{d_{m a x}} Ω_{B} (d_{j} ∶ d_{i}) - Ω_{B} (d_{i}) \end{array} \end{matrix}$ (19)

因此，局部气含率可通过以下表达式计算：

$α_{g} f_{i} = n_{i} V_{i}$ (20)

式中，fi 是第i组气泡占总气含率的分数；Vi 为第i组气泡的体积。Sauter平均直径d32的表达式可写为：

$\frac{1}{d_{32}} = \sum_{i = 1}^{N} \frac{f_{i}}{d_{i}}$ (21)

本文采用Luo聚并模型[35]，气泡聚并速率可表示为：

$Ω_{C} (d_{i} ∶ d_{j}) = ω_{C} (d_{i} ∶ d_{j}) p_{C} (d_{i} ∶ d_{j})$ (22)

气泡间碰撞频率的表达式为：

$ω_{C} (d_{i} ∶ d_{j}) = \frac{π}{4} {(d_{i} + d_{j})}^{2} n_{i} n_{j} {\dot{u}}_{i j}$ (23)

聚并概率的表达式为：

$p_{C} = e x p \{- c_{1} \frac{{[0.75 (1 + x_{i j}^{2}) (1 + x_{i j}^{3})]}^{1 / 2}}{{(ρ_{g} / ρ_{l} + 0.5)}^{1 / 2} {(1 + x_{i j})}^{3}} W e_{i j}^{1 / 2}\}$ (24)

其中

${\dot{u}}_{i j} = {(u_{i}^{2} + u_{j}^{2})}^{1 / 2}$ (25) $u_{i} = 1.43 {(ε d_{i})}^{1 / 3}$ (26)xij =di /dj(27) $W e_{i j} = \frac{ρ_{l} d_{i} {\dot{u}}_{i j}^{2}}{σ}$ (28)

气泡诱导湍流对气泡聚并过程有部分影响，具体表现主要有：(1) 对于湍流碰撞形式的气泡聚并，其气泡-气泡碰撞筒中两气泡所携带的湍动能都会受到周围湍流涡旋的控制，不同湍流机制作用下的湍流涡旋的平均脉动速度和耗散率很不相同，即气泡诱导湍流作用下的气泡可能具有更大的碰撞速度和更短的作用时间；(2) 对于大气泡尾涡卷吸夹带形式的气泡聚并，其尾涡作用时空范围内的小气泡在气泡诱导湍流的作用下发生聚并，但不同尺寸的气泡对气泡诱导湍流涡旋的响应又有所不同。由此可见，与气泡破碎的涡旋-气泡直接碰撞相比，虽然气泡诱导对气泡聚并过程的影响较为间接，但是同样十分复杂，也无法通过简单地调节模型参数实现。因此，为了更加深入揭示气泡诱导对破碎过程的作用机制，本研究中将问题进行了一定的简化，暂时忽略了气泡诱导对聚并过程的影响。希望在后续的研究中再建立考虑气泡诱导对气泡聚并过程复杂作用机制的理论模型。

2 BIT模型

2.1 考虑BIT涡旋的气泡破碎模型

气液鼓泡流动中气泡破碎主要是由湍流涡旋与气泡碰撞引起。因此，来流涡旋的特性尤为重要。气泡聚并破碎模型中对于来流涡旋的描述，大多基于柯尔莫哥洛夫的湍流-5/3次方律推导获得湍流平均脉动速度和耗散率等重要参数。但是，-5/3次方律的适用范围是各向同性均匀湍流，对于BIT这种在主流方向上有着更强脉动和更快耗散的湍流流动而言，其适用性值得商榷。

Shi等[28]在目前应用最为广泛的Luo and Svendsen破碎模型基础上，从BIT湍流能谱的波数κ-3特性和BIT湍动能产生和消亡的等效平衡关系出发，详细推导了BIT湍流动能能谱函数的基本形式，并建立了考虑气泡诱导湍流和剪切湍流共同作用的气泡破碎模型。关于该模型的理论求解等详细讨论，可以参考文献[29]。该模型与Luo and Svendsen模型表达式的对比见表3。

表3 气泡破碎模型方程的比较

Table 3 Comparison of bubble breakage models

项目	Luo and Svendsen模型	考虑BIT的气泡破碎模型
能谱函数	$E (κ) = C_{κ} ε^{\frac{2}{3}} κ^{- \frac{5}{3}}$ $C_{κ} \approx 1.5$	$E (κ) = δ_{l} C_{l} ε^{\frac{2}{3}} κ^{- \frac{5}{3}} + δ_{b} C_{b} \frac{α g U_{s l i p}}{ν} κ^{- 3}$ (29) $δ_{l}$ 、 $δ_{b}$ 为开关函数 $C_{l}$ 、 $C_{b}$ 取动态变化值
涡旋平均脉动速度	${\bar{u}}_{λ} = β^{\frac{1}{2}} {(ε λ)}^{\frac{1}{3}}$ $β \approx$ 2.0	${\bar{u}}_{λ} = C_{λ}^{\frac{1}{2}} \sqrt[]{C_{b} \frac{α g U_{s l i p}}{ν}} λ$ (30)
湍流涡旋数密度	$f_{λ} = \frac{C_{3} (1 - α)}{λ^{4}}$ $C_{3} = \frac{15}{2 {(2 π)}^{\frac{2}{3}} Γ (\frac{1}{3})} \approx 0.822$	$f_{λ} = \frac{C_{3} (1 - α)}{λ^{4}}$ (31) $C_{3} = \frac{12}{π C_{λ} {(2 π)}^{2}}$
涡旋-气泡碰撞概率密度函数	$ω_{B}^{T} (d_{i}, λ) = C_{4} (1 - α) n_{i} {(ε d_{i})}^{\frac{1}{3}} \frac{{(1 + ξ)}^{2}}{d_{i} ξ^{\frac{11}{3}}}$ $C_{4} = \frac{π}{4} C_{3} {(β)}^{\frac{1}{2}} \approx 0.923$	$ω_{B}^{T} (d_{i}, λ) = C_{4} (1 - α) n_{i} \frac{{(1 + ξ)}^{2}}{d_{i} ξ^{3}} \sqrt[]{C_{b} \frac{α g U_{s l i p}}{ν}}$ (32) $C_{4} = \frac{π}{4} C_{3} {(C_{λ} C_{b})}^{\frac{1}{2}}$
平均湍动能	$\bar{e} (d_{i}, λ) = \frac{π β}{12} ρ_{l} {(ε d_{i})}^{\frac{2}{3}} {d_{i}}^{3} ξ^{\frac{11}{3}}$	$\bar{e} (d_{i}, λ) = \frac{π}{12} C_{λ} ρ_{l} λ^{5} C_{b} \frac{α g U_{s l i p}}{ν}$ (33)
破碎速率	$Ω_{B} (d_{i} ∶ d_{j}) = 0.923 (1 - α)$ × $n_{i} {(\frac{ε}{d^{2}})}^{\frac{1}{3}} \int_{ξ_{m i n}}^{1} \frac{{(1 + ξ)}^{2}}{ξ^{\frac{11}{3}}} e x p (- \frac{12 σ C_{f}}{β ρ_{l} ε^{\frac{2}{3}} d_{i}^{\frac{5}{3}} ξ^{\frac{11}{3}}}) d ξ$	$Ω_{B} (d_{i} ∶ d_{j}) = 0.923 (1 - α) n_{i} {(\frac{ε}{d^{2}})}^{\frac{1}{3}}$ × $\int_{\frac{Λ}{d_{i}}}^{1} \frac{{(1 + ξ)}^{2}}{ξ^{\frac{11}{3}}} e x p (- \frac{12 σ C_{f}}{β ρ_{l} ε^{\frac{2}{3}} d_{i}^{\frac{5}{3}} ξ^{\frac{11}{3}}}) d ξ + C_{4} (1 - α)$ × (34) $n_{i} \sqrt[]{C_{b} \frac{α g U_{s l i p}}{ν}} \int_{ξ_{m i n}}^{\frac{Λ}{d_{i}}} \frac{{(1 + ξ)}^{2}}{ξ^{3}} e x p (- \frac{12 σ C_{f}}{β ρ_{l} C_{b} \frac{α g U_{s l i p}}{ν} d_{i}^{3} ξ^{5}}) d ξ$

新窗口打开| 下载CSV

2.2 考虑气泡涡旋响应的湍流相间扩散模型

为在LES中充分反映气泡上升过程中的摇摆行为特征[17]，并研究其作用和影响，对于亚格子尺度的湍流耗散作用模型应通过对曳力的空间滤波进行修正。

如图1所示，在亚格子尺度下，液相脉动可以形象地认为是小尺度湍流涡持续击打气泡表面，远小于气泡尺寸的湍流涡引起气泡表面变形从而引起质心振荡，同时气泡会因气泡尺寸相近的湍流涡产生动态形变从而引起质心位置移动。质心改变会在一定程度上导致气泡在非主流方向上的摇摆运动，从而影响液相的湍流运动，这种小范围的振荡行为可以在大涡模拟中有所反映。

图1

图1 亚格子尺度液相湍流脉动引起的气泡质心改变示意图

Fig.1 Schematic diagram of bubble mass centre movement caused by sub-grid scale liquid fluctuations

曳力的表达式可简化并分解写为如下形式：

$M_{D} = C_{g l} (U_{g} - U_{l}) = D_{g l} A_{g l} (U_{g} - U_{l})$ (35)

其中

$C_{g l} = \frac{3}{4} C_{d} \frac{α_{l} ρ_{g}}{d_{g}} |U_{g} - U_{l}| = \frac{1}{8} C_{d} A_{g l} ρ_{g} |U_{g} - U_{l}|$ (36) $A_{g l} = \frac{6 α_{g}}{d_{g}}$ (37)

经过空间滤波后：

$\bar{M_{k}} = D_{g l} [\bar{A_{g l}} (\bar{U_{g}} - \bar{U_{l}}) + \bar{A_{g l}^{'} (U'_{g} - U'_{l})}]$ (38)

对后一项使用涡流扩散假定：

$D_{g l} \bar{A_{g l}^{'} (U'_{g} - U'_{l})} = - \bar{C_{g l}} (\frac{ν_{t g}}{σ_{A g}} - \frac{ν_{t l}}{σ_{A l}}) \frac{\nabla \bar{A_{g l}}}{\bar{A_{g l}}}$ (39)

Favre平均后的速度可写为：

$\tilde{u} = \bar{u} + \frac{\bar{α_{k}^{'} u_{k}^{'}}}{\bar{α_{k}}}$ (40)

其中 $\bar{α_{k}^{'} u_{k}^{'}}$ 为与湍流扩散相关项，将式(39)代入式(38)并与湍流扩散力Favre平均后的形式进行比对，可写为：

$\bar{M_{k}} = C_{g l} (\tilde{U_{g}} - \tilde{U_{l}}) + M_{k}^{T u r b u l e n t D i f f u s i o n}$ (41)

值得注意的是，与RANS模拟中湍流扩散力的系综平均不同，此处的湍流耗散力相当于空间平均和相平均后的量，其扰动信息在每个单位时间内都需要实时计算。因此，再经过一定的数学变换后，式(41)可写为：

$\bar{M_{k}} = C_{g l} (\bar{U_{g}} - \bar{U_{l}}) - \bar{C_{g l}} \frac{ν_{t}}{σ_{a}} (\frac{1}{\bar{α_{l}}} + \frac{1}{\bar{α_{g}}}) \nabla \bar{α_{l}}$ (42)

再考虑式(14)气泡对涡旋相应的液相亚格子尺度湍流涡黏度模型，可将亚格子尺度湍流相间扩散模型写为：

$\begin{array}{l} \bar{M_{T D}} = C_{T D} \frac{3 {\bar{α}}_{g}}{4} \frac{ρ_{l}}{d_{b}} ({\bar{u}}_{l} - {\bar{u}}_{g}) \times \\ \frac{{(C_{s} Δ)}^{2} |S| (1 + C_{b u b b l e} {\bar{α}}_{G} \frac{λ}{d} {(\frac{1}{1 + S t_{S G S}})}^{\frac{3}{2}})}{σ_{A}} (\frac{1}{\bar{α_{l}}} + \frac{1}{\bar{α_{g}}}) \nabla \bar{α_{l}} \end{array}$ (43)

3 模型求解方法

本研究采用三维非稳态压力求解器求解模型方程。入口边界条件为速度入口，出口边界条件为压力出口，壁面条件为无滑移。模拟的时间步长根据模拟的体系和采用的模型在0.001~0.005 s范围内动态取值，并认为在100~120 s后达到了准稳态。当鼓泡塔内的模拟达到稳定后，获取了充分发展段处截面的时均气含率径向分布和气泡直径分布等数据。PBM模型采用离散方法将气泡大小按照等体积比方法离散求解。其中，BIT气泡破碎模型通过用户自定义函数 (UDF)植入时，涉及破碎速率表达式的一重积分和子气泡尺寸分布概率密度函数表达式的二重积分。本研究中采用的重积分数值积分算法经过严格的检验，与计算平台内置的通过Gamma函数对Luo and Svendsen破碎模型求解所得的气泡尺寸分布结果非常相近。

3.1 模拟的气液体系

表4为模拟的气液体系的各项参数。实验选取的气相为空气，液相为水，环境温度为室温。为了能够清晰地比较模拟参数的改变，表4中还简要列出了各模拟中主要模型的选择。

表4 模拟的气液体系实验参数及主要模型选择

Table 4 Experimental details of gas-liquid systems and model selections

Case	实验	塔径/m	高度/m	表观气速/(m/s)	静液位高度/m	测量高度/m	湍流模型	湍流扩散	破碎模型
Case 1	Gemello等[36-37]	0.4	3.6	0.16	1.6	1.5	RSM	—	Luo and Sevndsen
Case 2	Gemello等[36-37]	0.4	3.6	0.16	1.6	1.5	RSM	—	Luo and Sevndsen ΩB’(di ∶dj )=10ΩB(di ∶dj )
Case 3	Gemello等[36-37]	0.4	3.6	0.16	1.6	1.5	RSM	—	式(34)
Case 4	Guan等[38]	0.15	1.6	0.08	1.2	0.8	RSM	—	Luo and Sevndsen
Case 5	Guan等[38]	0.15	1.6	0.08	1.2	0.8	RSM	—	式(34)
Case 6	Sommerfeld等[17]	0.14	1.4	0.0029	0.65	0.325	LES	Burns等[39]	—
Case 7	Sommerfeld等[17]	0.14	1.4	0.0029	0.65	0.325	LES	式(43)	—

新窗口打开| 下载CSV

3.2 网格合理性验证

图2(a)为CFD-PBM模拟采用的网格，该网格在半径、圆周、高度方向上分别采用28(r)×64(θ)×100(z)个均匀分布节点。图2(b)为LES所采用的网格，与RANS模拟采用的划分方式整体相似，但考虑到LES对于网格密度的要求，该鼓泡塔网格划分在高度方向上保持均匀的Δz+ = 100，在横截面中心处网格大小采用dmean/ Δ = 0.6375，在径向上由Δr+ = 5开始以1. 2的增长率逐渐变大。该网格的分辨率合理地接近Milleli限制标准，并被认为足够精细到可以解析湍流的大尺度结构。采用该网格设置，该鼓泡塔的总网格数量约为95400个。

图2

图2 网格设置示意图

Fig.2 Schematic diagram of mesh set-up

如图3所示，分别增加了两种不同的网格密度对如图2(a)中所示的Grid 2网格进行网格合理性验证：即Grid 1采用了20(r)×40(θ)×80(z)个节点，总网格数量减半；Grid 3采用了36(r)×72(θ)×126(z)个节点，总网格数量加倍。基于同样范围的时间平均和同样径向位置的空间平均，三组网格的敏感性基本相似。图3(b)中的纵坐标为归一化后的轴向液速，此处用的归一化方法为“线性函数归一化(Min-Max scaling)”，即Ynorm=(Y-Ymin)/(Ymax-Ymin)。其中，Grid 2和Grid 3无论在气含率还是液相轴向速度的径向分布结果从趋势上看都较为接近实验值，仅有Grid 1的结果与实验值呈现略大一些的偏离。可以认为，在后续模拟中采用的网格设置不仅有更高的计算效率，还能够有效捕捉较为可靠的模拟结果。

图3

图3 网格合理性验证结果

Fig.3 Grid sensitivity test results

4 结果与讨论

4.1 BIT气泡破碎模型的影响

气泡破碎速率对于气泡尺寸分布的准确预测具有十分重要的意义。对于气液鼓泡流动，采用Luo and Svendsen破碎模型往往低估了破碎速率[40]，其原因很可能是该模型低估了与气泡碰撞的湍流涡旋所携带的湍动能。这也导致了Case 1该模型模拟的气含率和液相轴向速度远低于实验值，如图4所示。为此，在数值计算中采用了Chen等[41]提出的将Luo and Svendsen模型的破碎速率人为提高10倍的办法之后，Case 2的模拟结果与实验值非常接近。这说明在某些特定工况下，Luo and Svendsen原模型确实无法准确地描述破碎速率，需要进行一定的参数调节。

图4

图4 不同破碎模型模拟结果与实验结果对比

Fig.4 Comparison of experimental results with simulation results by different breakage models

在鼓泡反应器的气液流动中，气体鼓入塔内前，体系处于没有任何流动的静止状态。当气泡通过气体分布器鼓入塔内，液相湍流由气泡或气泡群诱导而成。尤其是当表观气速较大，塔内气泡十分密集时，这种气泡诱导形成的液相湍流会占据大量的湍流结构。因此，如果根据该情况修正原模型，将单相湍流能谱函数以具有波数κ-3特性的BIT能谱函数式(29)进行替代，并且同时修正湍流脉动速度、涡旋数密度等相关参数，就能够得到充分考虑BIT的气泡破碎模型。该模型无须人为调整破碎速率参数，但Case 3中其模拟结果却与Chen等 [40]的模拟结果非常相近。尤其在鼓泡塔中心区域，BIT气泡破碎模型的结果更为准确。这说明了在这种气泡诱导湍流占据主导的工况下，在模拟中采用BIT湍流能谱，更能准确描述湍流涡旋的动态行为特征。

图5比较了不同破碎模型获得的气泡尺寸分布模拟结果。原模型的模拟结果中，大气泡组占据了大量的体积分率。采用BIT气泡破碎模型后，气泡尺寸分布预测结果与人为将破碎速率提高10倍的气泡尺寸分布结果较为相似。但是，不能简单认为二者起作用的机制相似。从子气泡尺寸分布概率密度函数可以看出，人为提高破碎速率并不改变子气泡尺寸分布概率密度；而采用BIT湍流能谱对相关参数进行修正，子气泡尺寸分布概率密度会随诱导气泡的特征长度尺度在每个计算网格中动态改变，从而改变整体的计算结果。该BIT气泡破碎模型的子气泡尺寸分布概率密度特性的详细论述可参见文献[28-29]。

图5

图5 不同破碎模型模拟的气泡尺寸分布对比

Fig.5 Comparison of bubble size distributions by different breakage models

当鼓泡反应器的塔径和表观气速都减小时，如Case 4和Case 5的工况下，该BIT气泡破碎模型的准确性也同样得到验证。与Ug = 0.16 m/s的工况下鼓泡塔几乎处于完全湍动状态不同，在Ug = 0.08 m/s的工况下，鼓泡塔处于由均匀鼓泡向完全湍动转化的中间状态，湍流强度低于完全湍动状态。此时，BIT并不完全占据主导，剪切湍流在气泡较为稀疏的局部时空中也有部分贡献。但总体而言，气泡诱导湍流的作用依然不容忽视。图6比较了该工况下原模型和BIT气泡破碎模型对气含率、Sauter平均气泡尺寸的模拟结果。从模拟结果来看，气泡尺寸分布对于气含率预测结果的准确性有着较强的相关性。由于原模型估算的剪切湍流涡旋与该工况下湍流涡旋的湍动能差距较大，导致聚并速率远大于破碎速率，并因此大大高估了平均气泡尺寸。其气含率预测结果在鼓泡塔中心区域与实验值的偏差极为突出，暗示中心区存在大量大于20 mm的大气泡，但事实上这一预测结果并不准确。而BIT气泡破碎模型在气含率和平均气泡尺寸两个参数的模拟结果都较好地与实验结果吻合。从表面上看，这是因为BIT破碎模型预测的破碎速率和聚并速率较为匹配而取得的结果。但是，从更深的层次讲，BIT破碎模型的准确性得益于其所采用的能谱函数既考虑了剪切湍流又考虑了气泡诱导湍流在惯性子区的发展规律，较好地捕捉了该工况下参与气泡碰撞和破碎过程的湍流涡旋的脉动与耗散的特性。在BIT破碎模型中，分为剪切湍流和气泡诱导湍流两部分计算破碎速率，湍流耗散率ε仅影响剪切湍流涡旋导致的破碎，对气泡诱导湍流造成的破碎并不直接以湍流耗散率计算湍流涡旋的湍动能，不受湍流模型计算所得的湍流耗散率的影响。因此，模拟结果中气泡尺寸分布的改进基本上得益于BIT气泡破碎模型对湍流涡旋特性的准确描述。其影响机制可简要概括为：（1）BIT改变了与气泡碰撞的湍流涡旋的平均脉动速度uλ 的大小和数密度fλ 的分布情况，从而动态地修正了气泡破碎速率（对不同尺寸气泡破碎速率的修正比例不同，不是都×10），这对于气泡诱导湍流仅在部分时空范围占据主导的工况而言具有更加重要且明确的物理意义；（2）BIT破碎模型还隐性地改变了子气泡尺寸分布，这对最终通过数值模拟获得的气泡尺寸分布结果也有一定影响，但影响程度还需要进一步研究。

图6

图6 不同破碎模型模拟的气含率和气泡尺寸分布对比

Fig.6 Comparison of gas holdup and bubble size distributions by different breakage models

4.2 湍流相间扩散模型的影响

当鼓泡反应器的表观气速较小时，体系中气泡数量相对较少。同时，气泡与气泡、气泡与涡旋碰撞发生聚并、破碎的频率很低，此时采用PBM模拟的必要性较低。虽然湍流能量的整体来源还是鼓入液相的气泡，但是由于气泡诱导湍流快速耗散的特性，气泡诱导湍流主要存在于气泡尾涡中。气泡尾涡伴随着气泡上升快速移动变化，其影响的局部时空范围受限于气泡上升路径附近。虽然气泡诱导湍流的各项异性特征并没有发生变化，但这种局部的各项异性特征对整体流动而言可能并不占主导。考虑到这些因素，对于Case 6和Case 7工况的模拟可以采用LES湍流模型，而BIT的作用可以在SGS湍流涡黏度修正及SGS湍流相间扩散模型中进行反映。

图7为LES模拟的鼓泡塔纵截面液相速度矢量图，颜色深浅由气含率分布标示。无论是从时间平均值、瞬态分布还是不同的时间序列来看，LES都能较好地解析大部分的大尺度湍流涡。尤其是从图7(b)中能够在气含率较高的区域清晰分辨出气泡诱导的局部湍流涡，并且这些局部涡结构随着气泡摇摆上升动态变化。因此可以认为，除主方向外，其他方向上也存在着较为强烈的气液相间相互作用力。Sommerfeld等[17]在欧拉-拉格朗日模拟中全面考察了曳力、壁面润滑力、升力、浮力、虚拟质量力、Basset力的重要性，并认为Basset力是造成气泡及其尾涡摇摆上升的关键。但是在欧拉-欧拉模拟中，由于其短时平均的特性，Basset力的作用被严重弱化。Zhang等[42]在CFD-PBM模拟中全面考察了除曳力外的其他相间作用力对模拟结果准确性的影响，并认为湍流扩散力对气含率预测结果具有重要影响。

图7

图7 液相速度矢量图

Fig.7 Axial liquid velocity vectors

在大涡模拟中以不同湍流扩散模型预测的时间平均气泡上升速度径向分布对比如图8所示。在表观气速Ug = 0.0029 m/s、平均气泡尺寸为2.55 mm的工况下，鼓泡塔整体气含率较低，主方向上曳力的作用相对明显，气含率径向分布的结果一般较为准确，但是气泡上升速度这类受到局部作用影响较大的参数预测显得尤为困难。由图8可见，采用Burns湍流扩散力模型虽然能一定程度上提高预测的准确性，但依然较为偏离实验值；而采用式(43)的湍流相间扩散模型进行模拟，其结果与实验值更为接近。可以认为，该模型通过对曳力的空间滤波，较为精准地反映了气泡和湍流涡旋在亚格子尺度的动态响应，如图1所示。因其捕捉到气泡及其尾涡摇摆上升这一特性，所以能够更准确地反映在气泡上升速度这一重要参数。值得注意的是，因为气泡发生聚并、破碎的频率较低，体系中的气泡尺寸分布几乎完全取决于由于气体分布器形成的气泡尺寸分布。而且，气泡诱导的尾涡在局部的影响范围取决于跟气泡尺寸高度相关的特征长度尺度Λ = d / CD。因此，如果再把体系中的气泡尺寸分布加以考虑，则模拟结果更为接近实验值。

图8

图8 气泡上升速度径向分布

Fig.8 Time-averaged radial distribution of bubble rise velocity

从图8可以看出，在Euler/Euler-LES的工作中，无论是采用标准的Burns湍流扩散力模型还是修正后的亚格子尺度湍流扩散模型（SGS-TD），在预测气泡速度曲线时，在入口处使用多气泡尺寸模型(considering BSD)比对单一气泡尺寸(db=2.55 mm)表现更好。研究结果充分表明了不同尺寸的气泡的瞬态行为实际上与指定一个等效的平均气泡尺寸所描述的行为不同。针对Case 6可以看出，采用Burns 湍流扩散力所预测的气泡速度曲线显示出与实验结果的较大差异。相比之下，Case 7模拟预测的气相轴向速度分布的一致性虽然在反应器的中心区域仍与实验数据稍有差距，但在近壁区域更为贴合，这表明在LES模拟中加入修正的SGS-TD模型对气泡的径向迁移有显著影响。这可能是由于气泡在中心区域更为集中，导致该区域等效气泡直径db变化较大。在考虑多气泡尺寸模型的算例中，在中心区域分配了大约4 mm的气泡尺寸，相当于db/Δ∈(0.830,1.025)，该范围与最佳的LES网格标准db/Δ∈(0.670, 0.830)略有偏差，所以BSD的使用可能高估了主流动方向上的气泡波动。相对较小的气泡更倾向于聚集在壁面区域附近，对周围的湍流涡流响应更为敏感，从而更能准确模拟实验中壁面附近的真实情况。这进一步表明，使用修正的SGS-TD模型具有调节气泡分散行为的功能，从而更准确地预测了气含率梯度以及气泡横向分散行为。基于与Sommerfeld等[17]的欧拉/拉格朗日LES结果（Case B）的比较，表明欧拉/欧拉LES加上修正的SGS-TD模型在与实验数据比较时仍能提供一致的气泡动力学结果，并减少了计算量以及复杂的相间作用力封闭模型。

鼓泡塔中心区域某一点在主方向上的湍动能统计规律可以通过LES中获取该点轴向脉动速度，并进行时间关联和快速傅里叶变换，从而获得其湍流动能能谱。如图9所示，若在LES中采用Smargorinsky亚格子湍流涡黏度模型，其所得到的能谱曲线更接近由Pope[45]提出的标注能谱函数形式；若在LES中考虑BIT的作用，即采用考虑气泡涡旋响应的亚格子湍流涡黏度模型，修正的能谱曲线在惯性子区内有一部分的斜率更加接近-3。虽然-5/3定律是针对高Reynolds数的各向同性均匀湍流，但是Pope[45]的model spectrum指出，对于积分尺度Reynolds数Rel 在100~100000范围内，其惯性子区的斜率基本都符合-5/3。而对0.0029 m/s工况下的模拟结果以量纲分析大致估算，Rel 在100~1000范围内。因此，在不考虑BIT作用的前提下对LES结果进行处理得到了斜率接近-5/3的湍流能谱曲线是较为合理的。

图9

图9 通过LES获得的湍流动能能谱

Fig.9 Turbulent energy spectrum at middle point of z=325 mm plane by LES

若将气泡信号的特征频率定义为 $f_{b} = |u_{g} - u_{l}| / 2 π d_{b}$ ，则当气泡尺寸为2.55mm时，该特征频率为12.48 Hz。通过观察图9可知，修正的能谱曲线在惯性子区内斜率由-5/3到-3的转变点f1在14.80 Hz左右，与气泡信号的特征频率极为接近。可以认为，气泡所具有的湍动能在fb、f1频率附近对液相湍动能进行了部分能量传递。但是，这部分的能量如何沿着能量级串以及湍流拟序结构进行传递，尚未有准确结论。

4.3 BIT对气液相间传质的作用

基于上述分析与讨论，在气液两相鼓泡流动的几种具有代表性的工况下，BIT对于气泡-涡旋相互作用的数学描述起着非常关键的作用。而气液相间传质同样也是气泡-涡旋相互作用的过程，因此在描述气液相间传质系数kla的模型中考虑BIT的作用也有其必要性。本研究将kl 和a分别进行考虑，其中气液相间接触面积可以通过PBM模拟获得气泡尺寸分布，甚至可以在模拟中进一步对不同尺寸的气泡形状及其破碎过程中的形状及表面积变化的影响进行考虑，Shi等[43-44]对该部分工作进行了详尽的讨论。而对于kl，有较多因素会对其产生影响，如气泡尺寸、气液滑移速度、湍流耗散率、扩散速率等。研究人员对这些影响因素的作用进行了大量的研究，并分别发展了对流传质模型、停滞膜模型、双膜理论模型、溶质渗透模型、表面更新模型等理论模型[22-23]。随着研究的逐步深入，部分研究工作中逐渐反映出了湍流涡旋对气液相间传质过程的重要影响。Lamont等[24]在表面更新模型的基础上，通过假定湍流中从黏性涡到惯性涡都会对相间传质过程起作用，并以通过柯尔莫哥洛夫时间尺度 ${(ε_{l} / ν)}^{- 1 / 2}$ 计算了气泡表面更新速率，并以此建立了Eddy Cell模型，其表达式为：

$k_{l} = K D_{l}^{1 / 2} {(\frac{ε_{l}}{ν})}^{1 / 4}$ (44)

从这一角度看来，表面更新速率的计算，即湍流涡旋特性的合理描述是kl 的决定性因素。Han等[5, 25]在此基础上，建立了Eddy Contact模型。该模型通过考虑含能区、惯性子区、耗散区全范围的能谱函数，把湍流涡旋的尺寸、数量、脉动速度等因素与相间传质过程联系起来，其表达式可写为：

$k_{l} = 0.682 D_{l}^{1 / 2} \int_{κ_{m i n}}^{κ_{m a x}} [E {(κ)}^{5 / 4} κ^{3 / 4} d κ] / \int_{κ_{m i n}}^{κ_{m a x}} [E (κ) d κ]$ (45)

Eddy Contact模型中使用了全范围的能谱函数，意味着从相当于柯尔莫哥洛夫长度尺度η的最小湍流涡旋到与鼓泡塔尺寸相当的最大湍流涡旋都对相间传质过程产生了贡献。该机制的明确对于建立描述气液相间传质过程的理论模型具有十分重要的价值。但在该模型中，湍流涡旋对气泡表面的作用机制并不完全清晰。而且，从Pope[45]提出的能谱函数表达式 $E {(κ)}_{w i d e} = E {(κ)}_{i n e r t i a l} f_{L} (κ L) f_{η} (κ η)$ 中可以发现其关键核心在于惯性子区的能谱函数。若依然沿用柯尔莫哥洛夫-5/3次方律作为惯性子区的能谱，对于气液鼓泡流动而言，则无法正确体现出BIT的动态行为特点。因此，基于上述考虑，本研究中建立了考虑BIT的多尺度气液相间传质模型。

受到BIT在不同工况下占据不同时空范围内主导地位的思路启发，与气泡进行碰撞并发生物质传递的湍流涡旋大致可以分为三类：远大于气泡尺寸、与气泡尺寸相当、远小于气泡尺寸。对尺寸远大于气泡的涡旋而言，其与气泡的相互作用更多的是起到输运作用，这类接触导致的表面更新更多在于由剪切导致的气泡表面的扭曲变形或小气泡脱落，如图10(a)所示。而对于尺寸与气泡相当的湍流涡旋而言，其表面更新来源于气泡破碎生成的子气泡与母泡的面积差异，如图10(b)所示。因此，当表观气速较大、气泡聚并破碎活动较为剧烈时，可以假定这类型的表面更新对气液相间传质作主要贡献。而对于远小于气泡尺寸的湍流涡旋而言，其在惯性力的作用下会击打在气泡表面并与气泡表面张力或内部压力对抗。虽然其携带的湍动能往往不足以将气泡击碎，但很容易在气泡表面造成一个个“凹坑”，这类型的涡旋气泡接触引起的表面更新来源于气泡表面的局部拉伸变形和质心振荡导致的晃动变形，如图10(c)所示。

图10

图10 不同尺度的气泡与湍流涡旋接触示意图

Fig.10 Schematic diagram for mass transfer due to bubble-eddy contact

若从能谱函数开始推导，式(44)中的液相湍流耗散率可写为：

$ε = 2 ν \int_{κ_{m i n}}^{κ_{m a x}} κ^{2} E (κ) d κ$ (46)

若直接采用式(29)中的BIT湍流能谱函数，则BIT气液相间传质系数表达式可写为：

$k_{l} = K D_{l}^{1 / 2} {[2 \int_{κ_{m i n}}^{κ_{m a x}} (γ_{S T} C_{l} ε^{\frac{2}{3}} κ^{\frac{1}{3}} + γ_{B I T} C_{b} \frac{α_{g} U_{s l i p}}{ν} κ^{- 1}) f_{L} (κ L) f_{η} (κ η) d κ]}^{1 / 4}$ (47)

式中， $γ_{S T}$ 和 $γ_{B I T}$ 分别代表剪切湍流涡旋和气泡诱导湍流涡旋对整体相间系数的贡献，其取值可以通过对应的湍动能和总湍动能之比近似估算：

$γ_{S T} = \frac{\int_{κ_{m i n}}^{2 π / Λ} (C_{l} ε^{\frac{2}{3}} κ^{- \frac{5}{3}}) f_{L} (κ L) f_{η} (κ η) d κ}{\int_{κ_{m i n}}^{2 π / Λ} (C_{l} ε^{\frac{2}{3}} κ^{- \frac{5}{3}}) f_{L} (κ L) f_{η} (κ η) d κ + \int_{2 π / Λ}^{κ_{m a x}} (C_{b} \frac{α_{g} U_{s l i p}}{ν} κ^{- 3}) f_{L} (κ L) f_{η} (κ η) d κ}$ (48) $γ_{B I T} = \frac{\int_{2 π / Λ}^{κ_{m a x}} (C_{b} \frac{α_{g} U_{s l i p}}{ν} κ^{- 3}) f_{L} (κ L) f_{η} (κ η) d κ}{\int_{κ_{m i n}}^{2 π / Λ} (C_{l} ε^{\frac{2}{3}} κ^{- \frac{5}{3}}) f_{L} (κ L) f_{η} (κ η) d κ + \int_{2 π / Λ}^{κ_{m a x}} (C_{b} \frac{α_{g} U_{s l i p}}{ν} κ^{- 3}) f_{L} (κ L) f_{η} (κ η) d κ}$ (49)

图11中溶质对渗透模型、Eddy Cell模型、Eddy Contact模型以及BIT气液相间传质模型的模拟结果与Terasaka等[46]在不同表观气速下测量的单位体积传质系数进行了比较。其中，气液相间接触面积a的预测结果通过CFD-PBM模拟获得。为体现在相间传质模型中考虑湍流涡旋作用的重要性，同时也将溶质渗透模型(slip penetration model)的计算结果也列入比较。渗透模型的表达式可写为kl = 2F(Dl /πτe)0.5，其中τe为渗透的特征时间，可以用db/Uslip进行估算。溶质渗透模型和Eddy Cell模型都需要通过调整模型参数才能够获得与实验值接近的模拟结果。因此，这两个模型中模型参数K与Zhang等[47]建议的取值保持一致。对比结果显示，随着表观气速增大，溶质渗透模型和Eddy Cell模型预测获得的相间传质系数都低于实验值。其原因有可能是这两个模型都只考虑了与气泡发生接触的湍流涡旋的部分作用。与此相反，Eddy Contact模型和BIT气液相间传质模型通过能谱函数考虑了所有湍流涡旋对相间传质过程的贡献。当表观气速低于0.1 m/s时，Eddy Contact模型和BIT传质模型的模拟结果在整体上都与实验值较为接近；但更加值得注意的是，随着表观气速的进一步增大，BIT在气液相间传质过程中的重要性逐步显现，特别是当表观气速高于0.1 m/s后，BIT传质模型的模拟结果依然保持与实验值的良好吻合。

图11

图11 气液相间传质系数模拟结果对比

Fig.11 Comparison of the predicted and measured volumetric mass transfer coefficients

BIT传质模型的物理意义可以认为是参与相间传质的湍流涡旋的湍动能来源于剪切湍流和气泡诱导湍流两部分。当采用该模型进行数值计算时，困难之处在于式(47) 计算的相间传质系数是一个体系平均值，而由于气泡诱导湍流的强烈时空不均匀性，模型参数的取值从理论上来说应该随着BIT在不同局部时空作用范围变化。因此，此处BIT传质模型中模型参数的取值只好采用较为折中的方式，引入两种湍流机制对相间传质的贡献因子，但是这一估算方式可能带来一定的误差。当表观气速较低时，剪切湍流对相间传质的贡献较大，气泡诱导湍流的影响范围和对相间传质的作用效果逐渐加强直到与剪切湍流不相上下，式(49) 在较低气速的工况下可能略微过高估计了BIT的贡献因子。而Eddy Contact模型主要考虑剪切湍流的作用，其模拟结果在较低气速工况下较为准确；当表观气速高于0.1 m/s时，气泡诱导湍流的作用开始逐渐占据主导，BIT传质模型的作用逐渐凸显出来。但是，从BIT传质模型的理论框架而言，如果能够在数值计算中全面考虑到BIT从局部时空对整个体系的作用，在较低气速工况下也能够提高传质模型的准确性。

5 结论

通过对不同直径的鼓泡反应器不同工况的数值模拟，对气泡诱导湍流在气液鼓泡流动中的作用机理和气液相间传质特性展开了深入的探讨。在气液鼓泡流动中，气泡诱导湍流的影响不容忽视，但其作用方式随着BIT在不同时空中的主导地位有所变化，具体结论如下。

（1）当气泡诱导湍流在整个体系中占据主导时，其作用可以在CFD-PBM模拟中通过描述气泡-涡旋碰撞的BIT气泡破碎模型体现。

（2）当气泡诱导湍流仅在局部时空占据主导时，瞬态的气泡-涡旋对湍流的调制作用可以在LES中通过亚格子尺度的BIT湍流涡黏度模型体现。此外，通过LES空间滤波捕捉的瞬态扰动信息，并考虑气泡-涡旋动态响应作用可以隐式地在修正后的亚格子尺度湍流相间扩散模型体现气泡变形所引起的质心振荡行为，进一步揭示BIT在湍流能谱中所作贡献的重要性。

（3）虽然不同尺度的湍流涡旋对气液相间传质的作用方式不尽相同，但随着表观气速逐渐上升，BIT对气液相间传质过程的贡献越发突出。

（4）通过正确描述气泡诱导湍流的动态行为特征，在不同工况下对气含率、液相轴向速度、气泡尺寸分布、平均气泡尺寸、气泡上升速度、湍流动能、气液相间传质系数等关键流体力学参数的数值模拟均取得了与实验测量非常接近的结果，充分体现了BIT在气液鼓泡流动中的重要性。

符号说明

a	界面面积浓度，m-1
CD	等效曳力系数
D	鼓泡塔直径，m
Dl	扩散系数，m2/s
d	气泡直径，m
d32	Sauter平均直径，m
es	表面能增量，kg·m2/s2
ē	平均湍动能，kg·m2/s2
MD	曳力，N/m3
Mlift	升力，N/m3
MVM	虚拟质量力，N/m3
fi	第i组气泡占总气含率的分数
g	重力加速度，m/s2
H	高度位置，m
k	湍动能，m2/s2
kl	液相传质系数，m/s
l	长度尺度，m
n	单位体积数密度，m-3
t	时间，s
Re	Reynolds数
U	表观气速，m/s
Uslip	滑移速度, m/s
u	速度矢量，m/s
u′	u方向脉动速度，m/s
ūλ	湍流涡旋平均脉动速度，m/s
V	体积，m3
v′	v方向脉动速度，m/s
α	相含率
$γ$	部分湍动能与总湍动能之比
ε	湍流耗散率，m2/s3
η	柯尔莫哥洛夫长度尺度，m
κ	波数，m-1
λ	湍流涡旋特征长度尺度，m
Λ	气泡诱导湍流特征尺度，m
μeff	有效动力黏度，Pa·s
μl	分子动力黏度，Pa·s
μt	湍流涡动力黏度，Pa·s
υ	运动黏度，m2/s
ξ	涡旋长度尺度与气泡直径之比
ρ	密度，kg/m3
σ	表面张力，N/m
τ	黏性剪切应力，Pa
下角标
B	破碎
BIT	气泡诱导湍流
b	气泡
C	聚并
g	气体
i,j,k	气泡组
l	液体
ST	剪切湍流

关键字：优秀论文

上一篇：并行微通道内液液两相流及介尺度效应
下一篇：基于条件生成对抗网络的书法字笔画分割

栏目分类

热门排行

推荐信息

期刊知识

计及气泡诱导与剪切湍流的气泡破碎、湍流相间扩散及传质模型

1 数学模型

1.1 双流体模型

1.2 湍流模型

1.3 相间作用力模型

1.4 群平衡模型

2 BIT模型

2.1 考虑BIT涡旋的气泡破碎模型

2.2 考虑气泡涡旋响应的湍流相间扩散模型

图1

3 模型求解方法

3.1 模拟的气液体系

3.2 网格合理性验证

图2

图3

4 结果与讨论

4.1 BIT气泡破碎模型的影响

图4

图5

图6

4.2 湍流相间扩散模型的影响

图7

图8

图9

4.3 BIT对气液相间传质的作用

图10

图11

5 结 论

符 号 说 明

栏目分类

5 结论

符号说明