气液搅拌槽反应器因其气液混合性能好、传质传热效率高、相间接触面积大、操作维护方便、占地面积小等优点被广泛应用于生物制药、食品工程、水处理等工业领域^［1-2］。气液搅拌槽的性能主要与桨叶的结构、安装方式，以及搅拌容器的结构参数相关。桨叶作为搅拌槽的主要构件，对搅拌槽中气液流动状态、气体分布具有重要的影响^［3-4］。气液搅拌槽中桨叶的转动会在其周围产生一个强度巨大的剪切力场，该剪切力场对搅拌槽中气泡的聚并与破碎、气相的分布、气液两相间传质等具有重要的影响，其对搅拌槽的性能也起着决定性的作用。因此，探究气泡在剪切力下的聚并与破碎机理及运动与生长特性，对于气液两相流中桨叶的优化设计具有指导意义。

随着计算机技术的快速发展，计算流体力学(CFD)逐渐成为研究气液两相流的重要工具^［5］。在气泡模拟中研究者常常选取欧拉−欧拉多相流模型对气液两相流进行模拟。欧拉−欧拉多相流模型认为气液两相之间会相互渗透，该模型对每个气泡的运动不进行追踪，所以大大降低了计算量。但欧拉−欧拉多相流模型对气泡聚并与破碎的局部运动行为不能进行模拟，其假设气泡在运动过程中直径保持不变。气泡在高剪切力中运动时，气泡间存在剧烈的聚并与破碎行为，气泡的聚并与破碎行为对气相的局部气含率、气液相间接触面积、气液相间传质传热行为及气泡的运动有显著的影响。因此，当气泡在高剪切力中运动时不能忽略气泡间的聚并与破碎行为。

欧拉−欧拉多相流模型经常被应用于气液两相流模拟，而该模型存在一个严重的局限性，其假设每个气泡都具有相同的直径，且气泡在运动的过程中直径不发生变化。气泡直径大小是模拟时需要输入的参数，然而实际情况中气泡直径的大小很难预估且在运动过程中气泡的直径不断变化。在欧拉−欧拉多相流模型框架下耦合群体平衡模型(CFD−PBM)克服了这一局限性。近年来，CFD−PBM耦合模型可以有效地模拟气液两相间的运动行为，并能对气泡间的聚并与破碎行为进行模拟，因此该模型逐渐成为研究气泡运动行为的重要方法之一^［6］。Hulburt等^［7］和Randolph等^［8］率先提出了PBM模型，该模型通过跟踪离散相数密度函数随聚并与破碎的变化来描述离散相粒径随时间和空间的分布。Hulburt等^［7］将群平衡模型应用于旋转流场中液滴聚并与破碎的模拟中，并首次提出了聚并与破碎核函数。Lee等^［9］将群平衡模型应用于气液两相中气泡的聚并与破碎中，并用数密度函数代替概率密度分布。Luo等^［10］基于各向同性均匀湍流理论定义了气泡破碎率，认为气泡破碎是一种能量受限的现象，当湍流涡旋的湍流能量大于气泡承受的压力时，气泡发生破碎。Lehr等^［11］提出，只有当湍流涡的惯性力大于气泡表面张力引起的附加压力时，气泡才会破裂，并提出了压力约束的气泡破裂模型。

Zhu等^［12］分别应用CFD−PBM耦合模型与VOF（volume of fluid）模型对旋转桨内的流场进行模拟，研究表明，CFD−PBM耦合模型能更好地模拟旋转桨内气液两相分布。Petitti等^［13］对空气−水体系下的Rushton搅拌桨中气泡的聚并与破碎及相间传质进行了实验及数值模拟研究，发现群平衡模型的模拟结果与实验结果吻合较好。刘宝庆等^［14］采用群体平衡模型对剪切变稀体系下搅拌釜内的气液分布进行了模拟，研究发现，表观气速的增加使气泡的分布变得均匀，气泡尺寸随之增加。高颂等^［15］通过CFD−PBM耦合模型对离心泵内气泡的聚并与破碎进行了模拟研究，探究了入口体积气体率、入口气泡直径等对离心泵内气泡的聚并与破碎的影响。

现有的研究主要集中在对搅拌釜内气泡的聚并与破碎及运动状态进行研究，而关于气泡在高剪切力下的聚并与破碎特性的研究较少。因此，本文采用CFD−PBM耦合模型对高剪切力下气泡的聚并与破碎特性及运动状态进行模拟研究，探究剪切力强度、高剪切力下进气速度及气泡塔高度对气泡的聚并与破碎特性的影响规律，并结合聚并破碎模型对气泡在高剪切力下的聚并与破碎机理进行探究。

1 数学模型 1.1 双流体模型

在欧拉−欧拉双流体模型中离散相被看作连续的拟流体，并认为离散相与连续相之间可以相互渗透。该模型能较好地描述气液两相流，且计算量较小，能极大提高计算速度。双流体模型的基本方程主要包括质量守恒方程、动量守恒方程和能量守恒方程。在温度恒定且不考虑换热的情况下方程描述为

质量守恒方程(连续性方程)

$\qquad \frac{\partial \left({\alpha }_{i}{\rho }_{i}\right)}{\partial t} + \nabla \cdot \left({\alpha }_{i}{\rho }_{i}{\boldsymbol{u}}_{i}\right)=0 $

(1)

式中：$ {\rho }_{i} $为密度；$ {\alpha }_{i} $为相含率；${\boldsymbol{u}}_{i}$为速度；$ i $为g、$ {\rm{l}} $时分别代表气体、液体；t为时间。

动量守恒方程

$\qquad \begin{split} &\frac{\partial \left({\alpha }_{i}{\rho }_{i}{\boldsymbol{u}}_{i}\right)}{\partial t} + \nabla \cdot \left({\alpha }_{i}{\rho }_{i}{\boldsymbol{u}}_{i}{\boldsymbol{u}}_{i}\right)=\\&-{\alpha }_{i}\nabla p + \nabla \cdot \left({\alpha }_{i}{\tau }_{i}\right) + {\alpha }_{i}{\rho }_{i}\boldsymbol{g} + {\boldsymbol{M}}_{} \end{split} $

(2)

式中：$ p $为作用于两相间的压力；$ {\tau }_{i} $为应力张量；$\boldsymbol{g}$为重力加速度；${\boldsymbol{M}}_{}$为动量交换项，相间动量交换主要是由于相间力的作用，其包括曳力、升力、虚拟质量力、湍流分散力、壁面润滑力。

1.2 湍流模型

气泡在剪切力场中的流动状态为湍流。本文综合考虑计算精度及计算量后选择雷诺时均法中的标准$k-\varepsilon$湍流模型。标准$k-\varepsilon$湍流方程描述为

$ k $方程

$\quad\begin{array}{l}\dfrac{\partial \left({\alpha }_{{\rm{l}}}{\rho }_{{\rm{l}}}k\right)}{\partial t} + \nabla \cdot \left({\alpha }_{{\rm{l}}}{\rho }_{{\rm{l}}}{\boldsymbol{u}}_{{\rm{l}}}k\right)=\nabla \cdot \left[{\alpha }_{{\rm{l}}}\left({\mu }_{{\rm{L}},{\rm{l}}} + \dfrac{{\mu }_{{\rm{T}},{\rm{l}}}}{{\sigma }_{k}}\right)\nabla k\right] +\\ {\alpha }_{{\rm{l}}}\left({G}_{k} + {\rho }_{{\rm{l}}}\varepsilon \right)\\[-15pt] \end{array}$

(3)

$\varepsilon$方程

$\quad \begin{array}{l} \dfrac{\partial \left({\alpha }_{{\rm{l}}}{\rho }_{{\rm{l}}}\varepsilon \right)}{\partial t} + \nabla \cdot \left({\alpha }_{{\rm{l}}}{\rho }_{{\rm{l}}}{\boldsymbol{u}}_{{\rm{l}}}\varepsilon \right)=\nabla \cdot \left[{\alpha }_{{\rm{l}}}\left({\mu }_{{\rm{L}},{\rm{l}}} + \dfrac{{\mu }_{{\rm{T}},{\rm{l}}}}{{\sigma }_{\varepsilon }}\right)\nabla \varepsilon \right] +\\ {\alpha }_{{\rm{l}}}\left({C}_{1\varepsilon }\dfrac{\varepsilon }{k}{G}_{k} + {C}_{2\varepsilon }{\rho }_{{\rm{l}}}\dfrac{{\varepsilon }^{2}}{k}\right)\\[-15pt] \end{array} $

(4)

$\qquad {\mu }_{{\rm{T}},{\rm{l}}}={\rho }_{{\rm{l}}}{C}_{\mu }\dfrac{{k}^{2}}{\varepsilon } $

(5)

式中：k为湍动能；$ \varepsilon$为湍动能耗散率；u_l为液体速度；μ_L,l、μ_T,l分别为流体黏度和涡旋湍流黏度；G_k为平均速度梯度产生的湍动能 ；${C}_{1\varepsilon }=1.44$，${C}_{2\varepsilon }=1.92$，${C}_{\mu}=0.09$，$ {\sigma }_{k}=1.0 $，${\sigma }_{\varepsilon }=1.3$。

1.3 群体平衡模型 1.3.1 群体平衡方程

群平衡模型比较通用的描述方程为

$\qquad \frac{\partial n\left(V,x,t\right)}{\partial t} + \nabla \cdot \left[\boldsymbol{u}n\left(V,x,t\right)\right]=S $

(6)

式中：$ n\left(V,x,t\right) $为气泡分布的数密度函数，其值与气泡体积$ V $、空间变量$ x $及时间t有关；$ S $为气泡聚并与破碎的源项；u为气泡速度。

$ \begin{split} &S=\frac{1}{2}{\int }_{0}^{V}c\left(V-{V}^{'},{V}^{'},t\right)n\left(V-{V}^{'},x,t\right)n\left({V}^{'},x,t\right){\rm{d}}{V}^{'}-\\&{\int }_{0}^{\infty }c\left(V,{V}^{'},x,t\right)n\left(V,x,t\right)n\left(\left({V}^{'},x,t\right)\right){\rm{d}}{V}^{'} +\\& {\int }_{{\varOmega }\mathrm{v}}^{}b\left({V}^{'}\right)\beta \left(V|{V}^{'},x,t\right)n\left({V}^{'},x,t\right){\rm{d}}{V}^{'}-b\left(V,x,t\right)n\left(V,x,t\right) \end{split} $

(7)

式中：等式右侧4项依次分别表示两个体积较小的气泡聚并成体积为V气泡的源项，体积为V的气泡与其他体积的气泡聚并消灭的源项，体积较大的气泡破碎生成体积为V气泡的源项，体积为V的气泡破碎成其他体积气泡的源项，V'为 子气泡体积；c为气泡聚并核函数；b为破碎核函数 ；Ω_V为破碎气泡的体积 ；β为体积为V的气泡破裂成体积为V'的子气泡分布函数。

1.3.2 聚并核函数

气泡在湍流流动中的聚并可能是由多种机制引起的，如湍流波动、黏性剪切应力、湍流捕获、浮力和尾迹相互作用^［16］。本文选取Luo^［17］的聚并模型，直径为$ {d}_{i} $的气泡与直径为$ {d}_{j} $的气泡聚并表示为

$\qquad c\left({d}_{i},{d}_{j}\right)={\varpi }_{{\rm{c}}}{P}_{{\rm{c}}} $

(8)

$\qquad {\varpi }_{{\rm{c}}}=\frac{{\text{π}}}{4}\sqrt{2}{\varepsilon }^{1/3}{\left({d}_{i} + {d}_{j}\right)}^{2}{\left({d}_{i}^{2/3} + {d}_{j}^{2/3}\right)}^{1/2} $

(9)

$ {P}_{{\rm{c}}}\left({d}_{i},{d}_{j}\right)={\rm{exp}}\left(-0.4\frac{{\left[0.75\left(1 + {\xi }_{ij}^{2}\right)\left(1 + {\xi }_{ij}^{3}\right)\right]}^{1/2}}{{\left(\dfrac{{\rho }_{{\rm{g}}}}{{\rho }_{{\rm{l}}}} + 0.5\right)}^{1/2}{\left(1 + {\xi }_{ij}\right)}^{3}}{We}_{ij}^{1/2}\right) $

(10)

式中：$ {\varpi }_{{\rm{c}}} $为直径为$ {d}_{i}{\text{、}}{d}_{j} $的两气泡碰撞的频率；$ {P}_{{\rm{c}}} $为碰撞后聚并的效率；韦伯数${We}_{ij}=\dfrac{{\rho }_{{\rm{l}}}{d}_{i}{u}_{ij}}{\sigma }$；${\xi }_{ij}=\dfrac{{d}_{i}}{{d}_{j}}$；σ为液体表面张力函数。

1.3.3 破碎核函数

气泡的破碎过程非常复杂，且与液体的流动及气液之间的相互作用密切相关。现广泛使用的破碎模型为Luo & Svendsen模型与Leher模型。本文采用Luo & Svendsen模型，该模型认为当气泡与湍流涡发生碰撞，且湍流涡的能量大于气泡破碎所需能量时气泡发生破碎，一个体积为${V}_{i}$的气泡破碎为两个体积分别为${V}_{j}$和$ \left({V}_{i}-{V}_{j}\right) $的子气泡，其破碎模型描述为

$\qquad b\left({V}_{i}\right)=\frac{1}{2}{\int }_{0}^{{V}_{i}}{\varOmega }\left(V,{V}_{i}\right){\rm{d}}V $

(11)

$\qquad \begin{array}{l} {\varOmega }\left({V}_{j},{V}_{i}\right)=0.923\left(1-{\alpha }_{{\rm{g}}}\right){\left(\dfrac{{\varepsilon }_{{\rm{l}}}}{{d}_{i}^{2}}\right)}^{1/3}{\int }_{{\xi }_{{\rm{min}}}}^{1}\dfrac{{\left(1 + \xi \right)}^{2}}{{\xi }^{11/3}}\\{\rm{exp}}\left(-\dfrac{12{c}_{{\rm{f}}}\sigma }{2{\rho }_{{\rm{l}}}{\varepsilon }_{{\rm{l}}}^{2/3}{d}_{i}^{5/3}{\xi }^{11/3}}\right){\rm{d}}\xi\\[-15pt] \end{array} $

(12)

式中：${\varOmega }\left({V}_{j},{V}_{i}\right)$为破碎频率，$ {\xi }_{{\rm{min}}}={\lambda }_{{\rm{min}}}/{d}_{i} $，$ {\lambda }_{{\rm{min}}}=\left(11.4-31.4\right)\eta $，$\eta ={\left({\mu }_{{\rm{l}}}/{\rho }_{{\rm{l}}}\right)}^{0.75}/{\varepsilon }_{{\rm{l}}}^{0.25}$，${c}_{{\rm{f}}}= {\left(\dfrac{{V}_{j}}{{V}_{i}}\right)}^{2/3} + {\left(1-\dfrac{{V}_{j}}{{V}_{i}}\right)}^{2/3}-1$；ξ为涡流大小与气泡大小之比 ；ε_l为气泡的当地分数。

2 数值方法 2.1 高剪切力环境的构建

在气泡塔底部两侧分别通入不同速度的液体，由于两侧液体间存在速度差，因此在气泡塔中形成了一个剪切力场，通过改变两侧液速的差值进而改变剪切力的大小。为了避免因侧壁面的存在所产生涡流影响剪切力场，将所有的侧壁面均设置为速度进口，沿流向通入一定速度的液体，从而实现气泡在高剪切力中的运动。

2.2 几何模型

本文采用的几何模型及边界条件如图1所示。长方体气泡塔高为H，底面为0.2 m × 0.2 m，底面中央开有0.2 m × 0.01 m的进气口。初始时气泡塔内充满水，在底面Y<0侧通入的水的速度较大，在底面Y>0侧水以0.2 m·s⁻¹的速度通入塔内。在塔Y<0侧壁面沿Z轴方向通入速度较大的液体，在Y>0侧壁面沿Z轴方向以0.2 m·s⁻¹的速度通入液体，液体速度如表1所示。所有气体、液体进口均为速度进口，塔上部出口为压力出口。

2.3 数值模拟方案

本文利用ANSYS Fluent软件对气泡在高剪切力中的运动状态进行模拟，群平衡方程采用离散方法(class method，CM)进行求解，将气泡按直径大小分为11组子气泡，各子气泡组的直径如表2所示。气液流动状态为湍流，在综合考虑计算精度及计算量的情况下本文选用雷诺时均法中的标准k−ε湍流模型。为了获得一个稳定的剪切力场，模拟的前5 s不通入气体，塔底部及侧面仅通入液体，模拟时间步长为0.001 s，5 s时塔内液体速度分布如图2所示。5 s时，底面入口通入体积分数为1，气泡直径为4 mm，速度为0.02 m·s⁻¹的气体，时间步长改为0.000 5 s。通过监测发现，通入气体5 s后，塔内气液流动达到准稳态，在后续的研究中无特别说明时所取数据均为通气后15 s、Z=0.5 m、X=0 m处的数据。

3 计算结果与分析 3.1 剪切力对气泡聚并与破碎的影响

通过改变液体速度来构建不同强度的剪切力场。图3为不同液体速度下液体速度梯度的径向分布。由图中可知，当两侧通入液体速度差值增大时，液体速度梯度随之增加，液体间的剪切力增大。在侧面处，液体的速度梯度较大，这是因为侧面沿Z轴方向通入一定速度的液体，而底部液体速度方向与重力方向相反，在流动过程中内部液体速度减小；内部液体速度梯度主要分布在−0.32<r/R<0.70范围内，其中：r 为该位置至Y=0的距离；R为气泡塔半径。在后续研究中主要探究该处剪切力对气泡聚并与破碎的影响。

图4为不同液体速度下气含率径向分布。由图可知，当液体速度较大时，气体主要分布在中央，且气含率较大；当液体速度减小时，气体的分布向Y轴右侧移动，气含率的最大值随之减小。因为液体速度较大，气体将随液体以较大速度流过塔体，导致气体没有足够的时间向两侧运动，所以气体主要分布在r/R=0附近；液体速度减小，气体流经塔体的时间变长，气体在相间作用力升力的作用下，向液体速度较小处流动，使气体分布范围变大，气含率的最大值随之减小。

图5为气泡数密度随剪切力的变化。由图中可知；当剪切强度较小时，气泡尺寸分布较宽，气相主要由直径为3～9 mm的气泡组成；随着剪切力增大，气泡尺寸分布逐渐变窄，小直径气泡的数量显著增多。因此，气泡在较小的剪切力场中运动时，气泡间主要以聚并为主，随着剪切强度的增加，气泡聚并逐渐减弱，气泡破碎逐渐增强，气泡直径逐渐减小。

3.1.1 剪切力对气泡聚并的影响

图6为不同剪切力下气泡组径向分布。图6(a)、(c)、(e)为不同工况下聚并产生的子气泡径向分布。由图6(e)可知，工况C下聚并产生的子气泡主要分布在−0.34<r/R<0.85范围内，聚并产生的气泡体积约占气相的87 %，且以大直径子气泡组(Bin−0、Bin−1、Bin−2)为主，可见在该剪切力强度下气泡以聚并为主。由图6(a)、(c)可知，工况A、工况B下随着剪切力强度的增大，聚并产生的子气泡体积占气相的比例逐渐减小，且主要以小直径子气泡组(Bin−4)为主。随着剪切力强度的增大，气泡聚并行为逐渐减弱，聚并产生的气泡直径逐渐减小。当气泡在剪切流场中运动时，气泡在剪切力作用下沿剪切力方向拉伸，气泡直径变大，当剪切强度较小时，剪切力对气泡的作用小于气泡具有的表面能，气泡保持稳定的临界直径增大；当剪切强度较大时，剪切力对气泡的作用增强，气泡保持稳定的临界直径减小。

3.1.2 剪切力对气泡破碎的影响

图6(b)、(d)、(f)分别为工况 A、工况 B、工况 C下破碎产生的子气泡组径向分布。由图6(f)可知，工况 C下破碎产生的小气泡约占气相的5 %，且以大直径的子气泡组(Bin−6)居多；工况 B下破碎产生的小气泡体积约占气相的60 %，通入的子气泡组Bin−5体积约占气相的20 %，破碎产生的子气泡组中Bin−6～Bin−10的体积分数逐渐减小；在工况 A下破碎产生的小气泡体积占气相的95 %，其中Bin−6～Bin−10占比逐渐增大。气泡在剪切力场中运动时，气泡的破碎主要由黏性剪切作用引起。气泡在剪切力作用下发生形变，当气泡直径大于临界直径时，气泡将发生破碎形成小气泡；当气泡直径较大时，气泡尾部在剪切力作用下会发生脱离并产生小气泡，所以剪切力增大时气泡临界直径将会减小。随着剪切力强度的增大，由黏性剪切引起的气泡破碎效应逐渐增强，气相中小气泡的体积分数逐渐增大。

3.2 进气速度对气泡聚并与破碎影响

在3.1节的基础上，本文研究在工况 B的剪切强度下进气速度v_g对气泡聚并与破碎特性的影响，进气速度分别为0.02、0.06、0.12 m·s⁻¹。图7为不同进气速度下液体速度梯度的径向分布。从图中可以看出，进气速度的变化对液体速度梯度分布几乎无影响，在该剪切流中液体速度与进气速度之间为单项耦合，气体在流场中跟随液体运动，流场中液体速度对气体速度起着决定性作用。图8为气泡数密度随不同进气速度的变化。由图中可知，在相同剪切强度下，随着气体速度增大，气泡的尺寸分布变宽，子气泡组Bin−2～Bin−8的数量逐渐增多。因为当气体速度增大时，塔内气体含量逐渐增多，各直径子气泡组的数量相应增多，但小气泡组间的直径跨度较小，所以小气泡的数量变化不明显。

图9为不同剪切力下气泡组径向分布。由图9(a)、(c)、(f)中可以看出，在该剪切强度下，随着进气速度的增大，气相中大直径子气泡组的含量逐渐增加，例如，Bin−2子气泡组的体积分数最大值由2%增加至20%，所以气泡间的聚并效应随气速的增大而增强。从图9(b)、(d)、(e)中可得，随着进气速度的增大，气泡破碎产生的小气泡逐渐减小，如子气泡组Bin−8 的体积分数最大值由22%减小至13%，所以气泡间的破碎随着进气速度的增大而减弱。因为剪切力场的强度主要由液体之间的速度差决定，进气速度对剪切强度无明显的影响，进气速度的增大使液体中的气含率增加，随着气含率的增加，气泡间碰撞的频率增大，使气泡间聚并效应增加，大直径子气泡组的含量随之增大。

3.3 气泡塔高径比对气泡聚并破碎的影响

图10为不同气泡塔高径比（H/W）下液体速度梯度的径向分布，其中：H、W分别为气泡塔的高度、直径。从图中可以看出，不同高径比下液体速度梯度的径向分布基本相同。因为不同高径比并未改变液体之间的速度差，剪切强度主要由液体之间的速度差决定，所以不同高径比下液体速度梯度不变。图11为气泡数密度随气泡塔高径比的变化。由图中可得，不同高径比下气泡的尺寸分布基本相同，各子气泡组的数量基本相同，气泡聚并与破碎不随气泡塔高径比的变化而变化。因为高径比的改变仅造成液体压力的轴向分布变化，对气含率大小、分布，以及液体速度的大小、分布等均没有影响，所以高径比的变化对气泡的聚并与破碎未有显著影响。

4 结　论

为了研究气泡在搅拌槽及离心泵等旋转桨叶作用下的聚并与破碎特性，根据旋转桨叶的工作特点，将其简化为气泡在高剪切力场中的运动。应用CFD−PBM耦合模型对气泡在不同剪切强度下的运动进行数值模拟，探究了剪切力强度、高剪切力下气速及气泡塔高径比对气泡聚并与破碎的影响，主要结论为：

（1）气泡在不同剪切强度中运动时，当剪切强度较小时，气泡间主要存在聚并行为，气泡尺寸分布较宽，且气相主要由大直径的气泡组成；随着剪切强度的增强，气泡间的破碎行为逐渐增强，气泡尺寸分布变窄，气相主要由破碎产生的小气泡组成。

（2）在高剪切强度下对0.02、0.06和0.12 m·s⁻¹三种进气速度下的气泡运动进行数值模拟。结果表明，随着进气速度的增大，气泡尺寸分布逐渐变宽，大直径子气泡体积分数增加，气泡间的聚并效应增强。

（3）在高剪切力下气泡塔高径比增大时，气泡塔内液体压力的轴向分布发生改变，而气泡塔内的剪切强度几乎不变，气相的尺寸分布及组成气相的各直径气泡数量均未发生明显变化，其对气泡聚并与破碎行为无明显影响。