POF：浮力气泡捕获疏水性颗粒的数值模拟研究

1年前浏览3063

摘要

当一个浮力气泡在疏水性颗粒悬浮液中上升时，与气泡碰撞的颗粒通常会附着在气泡上，从而随气泡一同向上运动。在此，我们运用完全解析的数值模拟方法，研究在含有颗粒悬浮液的区域内，单个气泡的捕获速率。对气泡进行足够长时间的跟踪，以收集多个颗粒，我们探究了在适中的气泡雷诺数下，控制参数对颗粒积累速率的影响。结果表明，捕获速率对体积分数的依赖程度最强，随着气泡上积累的颗粒阻碍其他颗粒的附着，颗粒的捕获速率会迅速下降。随着气泡的可变形性增强，捕获效率略有降低，这是因为气泡底部的颗粒会使气泡拉长，减小其横截面积；而捕获效率会随着颗粒尺寸的增大而提高。对于此处模拟的参数范围而言，其他参数的影响较小。

1、引言

泡沫浮选是矿物加工中广泛应用的一种技术，通过利用矿物与脉石在疏水性和亲水性上的差异来实现二者的分离。例如，在铜矿开采中，矿石在球磨机中被研磨成细粉，理想情况下，粉末足够细，使得每个颗粒主要由单一纯净物质构成。这些颗粒与水混合形成矿浆，通常还会向混合物中添加一种称为 “捕收剂” 的化学试剂，使富含铜的颗粒具有疏水性。注入浮选槽底部的气泡会捕获疏水性颗粒，并将它们带到浮选槽表面，在那里可以将其撇去，而亲水性颗粒则沉淀在底部。通常，气泡尺寸范围在0.2至2毫米之间，而颗粒尺寸范围在10至200微米之间。通过泡沫浮选生产的铜精矿按重量计含铜10% - 40%，而尾矿（即提取矿物后剩余的物质）含铜量在0.05%至0.5%之间（Neethling和Cillier，2012）。在泡沫中收集到的高品位精矿，通常会通过熔炼或电解进一步提纯以获得纯金属。除了矿物加工领域，泡沫浮选还应用于废物处理行业，用于从水中分离悬浮的油类、脂肪和其他污染物；在回收行业也有应用，比如在再生纸脱墨以及不同类型塑料的分离中。例如，King（2001年）、Fuerstenau和Han（2003年）以及Fuerstenau（2007年）等已详细阐述了泡沫浮选的基本原理。

浮选过程中涉及的微观过程已通过大量实验展开研究。Huang（2011）、Verrelli（2011）以及Xing（2017）利用高速摄像机研究了颗粒在单个气泡上的附着情况，而Ding（2023）则详细考察了颗粒从气泡上脱离的现象。颗粒与气泡的附着在很大程度上取决于其疏水性，疏水性由三相线上固体表面与气 - 液界面之间的接触角决定，并且只有当接触角超过某一临界值时，气泡与颗粒才会发生附着（Kruszelnicki，2024）。通过测量接触角，能够确定矿物颗粒的表面性质，并选择合适的化学试剂来提高浮选的选择性和回收率。尽管接触角是影响选择性的关键因素，但由于系统的复杂性以及其他多种可能影响结果的因素，目前还无法直接依据接触角来预测回收率。一个重要的变量是颗粒与气泡的相对尺寸。如果颗粒过小，它们会顺着气泡周围的流线流动而不发生碰撞，从而进入尾矿。反之，如果颗粒过大，那么每颗颗粒都需要多个气泡才能托起，这使得从脉石材料中分离变得困难。关于泡沫浮选中颗粒尺寸影响的研究包括Dobby和Finch（1987）以及Feng和Aldrich（1999）。近期综述方面，关于细颗粒的研究可参考Farrokhpay（2021），关于粗颗粒的研究可参考Anzoom（2024）。气泡的可变形性也很重要，因为可变形性更强的气泡具有更大的扫掠面积，因而有可能收集更多颗粒，并且不同的气泡形状会导致不同的上升速度。Yan（2021）预测了负载颗粒气泡的终端上升速度，其中考虑了颗粒覆盖率和气泡可变形性的影响，Wang（2019）研究了负载颗粒气泡的行为，并得出了其阻力系数的修正因子。此外，搅拌产生的湍流程度对于促使气泡与颗粒接触至关重要。关于湍流对浮选影响的实验研究有Brady（2006）、Evans（2008）以及Darabi（2020），而Ngo - Cong（2018）从理论上证明了各向同性湍流有利于气泡 - 颗粒碰撞。然而，如果湍流程度过高，颗粒可能会从气泡上脱落，从而影响分离效果。Nguyen（2016）综述了湍流在浮选中的影响，并指出目前关于气泡 - 颗粒湍流脱附的模型仍不完善。

基于实验研究数据，人们开发了各种模型来辅助设计和优化浮选过程。Gharai和Venugopal（2016）将这些模型分为基于化学类比的动力学模型、概率模型、群体平衡模型和确定性模型。已提出一些动力学模型，这些模型由带有经验拟合参数的常微分方程组成，用于描述颗粒浓度的变化。Mao和Yoon（1997）提出了一阶动力学模型，但浮选动力学可能受其他因素影响，Yalcin和Kelebek（2011）发现二阶动力学模型拟合效果更佳。Lee（1987）和Bascur（2000）提出了基于群体平衡模型的浮选框架。Schuhmann（1942）基于一系列事件的成功概率提出了回收率的概率模型。Duan（2003）根据碰撞效率、附着效率和稳定效率计算了浮选速率常数，而Yoon和Mao（1996）考虑了颗粒与气泡之间的表面力。Williams和Meloy（1983）、Persechini（2000）以及Popli（2015）研究了其他描述浮选过程随时间动态变化的模型。在Fu（2020）的研究中，机器学习和动力学研究也被用于探究颗粒大小和浮选时间的影响。尽管已有多种浮选模型，但Gharai和Venugopal（2016）指出，仍需更深入理解整个浮选过程，以助力设计与优化工作。

在过去二十年里，计算机模拟已成为研究浮选过程的有力工具。其大致可分为三类：双流体模型，该模型将所有相都视为相互贯穿的连续介质；一类是对气相进行完全解析，而将固相视为点粒子的模型；还有一类是对气相和固相都进行完全解析的模型。Wang（2018）已完成了对浮选数值模拟的全面综述。

双流体模型的应用实例包括：Koh（2000年）和Kostoglou（2007年）采用欧拉 - 欧拉（E - E）湍流模型，将气体和颗粒都视为连续介质，以此估算气泡 - 颗粒的碰撞速率；Koh和Schwarz（2008年）运用E - E多相流方法预测浮选槽中的附着速率；Karimi（2014年）同样采用带有k−ε湍流模型的E - E框架来预测浮选速率，并通过实验数据进行验证。

有几位学者对气相进行了完全解析，并对固相采用点粒子模型。Liu和Luo（2018）结合流体体积法（VOF）与离散粒子法（DPM），模拟了一个和两个气泡与固体颗粒的碰撞，不过这是针对二维流动的模拟。van Sint Annaland（2005）运用前沿追踪法和DPM模拟了气泡在固体颗粒悬浮液中的上升过程。Liu和Schwarz（2009a）以及Liu和Schwarz（2009b）采用点粒子法和k−ε湍流模型，研究了单个气泡与多个颗粒之间的相互作用，而Bogner（2017）则采用格子玻尔兹曼方法和点粒子进行研究。在诸如Maxwell（2012）和Gao（2014）的研究中，运用离散元法（DEM）对疏水性颗粒与单个气泡进行建模。Xia（2023）结合前沿追踪法和DEM来模拟颗粒与气泡之间的附着与脱离。Washino（2023）最近的一项研究使用DEM方法预测了作用于具有广泛接触角范围的非球形颗粒上的毛细力。Fayed和Ragab（2013）将气泡和颗粒都视为点质量，并对均匀各向同性湍流中流体的速度波动进行了完全解析。

对气泡和颗粒都进行完全解析的模拟相对较少，并且大多数此类模拟假定颗粒是亲水性的，因此不会发生附着现象。Baltussen（2017）结合了前沿追踪法和浸入边界法，发现了有效阻力与相体积分数之间的相关性。Sasic（2014）对气泡采用流体体积法（VOF），对颗粒采用浸入边界法，而Qin和Suckale（2017）对气泡使用水平集方法。对于颗粒附着到气泡的情况，Dinariev和Evseev（2018）使用他们所谓的密度泛函流体动力学来解析三相混合物，不过其结果仅局限于一些图片，并未包含定量分析。

在本研究中，我们对泡沫浮选进行了数值研究，重点关注单个大气泡在含有多个重质疏水性颗粒的流体中移动的情况。颗粒会附着在气泡上，并且我们采用Qian（2003年）提出的广义纳维边界条件（GNBC）连续介质模型来考虑移动接触线，该模型允许颗粒表面存在滑移。尽管主要出于计算原因，我们大多使用相对较大的颗粒（因为它们更容易解析），但由于在研磨过程中节能以及其他经济和环境方面的原因，粗颗粒浮选在工业上受到了相当大的关注（Anzoom，2024年）。本研究的主要目标是探究收集速率如何依赖于各种控制参数，包括气泡的可变形性、固体颗粒的体积分数以及气泡与颗粒的尺寸比。我们还研究了一系列的接触角和接触线滑移速度。为了研究控制参数变化的影响，我们引入了一个简单模型来预测气泡的雷诺数和颗粒的捕获速率，该模型包含我们通过拟合数据确定的系数。通过比较各种情况下的拟合系数，我们得到了一种简洁的比较方式。尽管只研究单个气泡是一种显著的简化，但我们认为这不仅是一个重要的极限情况，能够提供重要的见解，而且对于更大、更复杂系统的研究来说，它也是一个必要的基准。

2、数值方法

三相不可压缩流动由N-S方程和连续性方程控制。

在这里，是速度，是密度，是压力，是重力加速度。气泡和液体中的应力为，其中是粘度，是表面张力项。是是一个施加的罚力项，用于保持固相的刚性，并且取值为：

其中

满足固体颗粒运动方程

。这里，

是颗粒质心速度，

是角速度，r是距颗粒质心的距离。为了识别不同的相并指定它们的属性，我们使用两组指示函数，

用于识别气相，

用于识别固体颗粒相，

我们使用与从固定流体网格插值得到的流体速度一同移动的连接标记点，来追踪流体 - 流体界面。表面张力也在标记点处确定，并且它与指示函数的梯度一起，利用Peskin（1977）引入的平滑函数，从标记点平滑到固定网格上。然后通过求解泊松方程（Unverdi和Tryggvason，1992），由其梯度构建指示函数，并将表面张力添加到离散的纳维 - 斯托克斯方程中。这种方法通常被称为前沿追踪法（Tryggvason等人，2011）。球形颗粒通过平流其质心来移动，并且指示函数根据质心的位置通过解析方式构建（Lu等人，2023）。由于前沿追踪法不能精确守恒体积，我们每隔十个时间步长，通过沿垂直于界面的方向移动表面点，来校正气泡的体积，以维持气泡的原始体积。这种校正量非常小，而对于颗粒，无需进行校正。

纳维 - 斯托克斯方程采用有限体积法在交错网格上求解，并使用二阶投影法进行时间积分。对流项采用QUICK格式，扩散项采用二阶中心差分格式，压力方程通过HYPRE方法求解（Falgout和Yang，2002）。在每个时间步结束时，颗粒内部的变形率并不完全为零，为了在下一个时间步强制实现刚体运动，我们将式（2）定义的体力加到纳维 - 斯托克斯方程中。为了得到和，我们对每个颗粒进行积分：

其中，是固体颗粒的指示函数。使用指示函数的平方能实现更快的收敛（Lu等人，2023）。此处采用的策略与常用于固体颗粒模拟的浸入边界法类似（Uhlmann，2005），但标准浸入边界法通常处理单相纳维 - 斯托克斯方程，而我们在此使用双流体版本。因此，对固体颗粒的处理与对气泡的处理非常相似，这遵循了Yabe等人（2001）提出的思路。这里使用的代码基于一个内部代码，该内部代码用于大量气泡流模拟（例如，近期应用可参考Lu和Tryggvason，2019；Liu等人，2023），以及一个为模拟固体颗粒而开发的代码，如Lu等人（2023）所述。

当气泡与疏水性颗粒碰撞时，我们假设气泡与固体颗粒之间的液膜瞬间排干，使得气泡能够立即捕获颗粒。目前，如果代表气泡界面的一个前沿单元与固体表面相互靠近且它们的相对切向速度较小时，我们通过将该前沿单元移动到固体表面，来使固体附着到气泡上。在实际情况中，附着过程更为复杂，因为气泡与固体之间的液膜排干需要时间，这可能取决于气泡和固体的表面条件。尽管这将涉及额外的代码开发，但可以添加更复杂的模型来描述这些复杂情况。接触线基于广义纳维边界条件（GNBC）模型移动。接触线的滑移速度由粘性应力和不平衡的杨氏应力决定，由以下方程给出：

在此，，其中是液相和气相的平均黏度，k是滑移长度。不平衡的杨氏应力是作用在三相点处的表面力的结果，而黏性应力可以忽略不计，因为它远小于杨氏应力（Yamamoto 等人，2013 年）。因此，滑移速度可写为

在此，σ是液-气表面张力， θd是动态接触角，θo是静态接触角，θc是接触线狄拉克函数，t是沿固体表面的切向向量。方程（6）表明，滑移速度与不平衡力成正比，因此我们写成

。如果我们用

来近似狄拉克函数（其中h是网格间距），那么

，由此可以计算出k。一般来说，滑移长度是未知的，在这里我们根据经验选择一个

的值，使得颗粒能平稳地部分嵌入气泡。如后文所示，结果对

的具体取值并不敏感。在附录中，我们展示了当固体颗粒附着在气泡上时，经过较长时间后静态接触角会恢复，关于该实现方法的其他测试，可参考 Zeng 等人（2024）的研究。此前已有几篇论文使用前沿追踪法，将广义纳维边界条件（GNBC）用于静止表面的模拟，包括 Yamamoto 等人（2013、2014、2016），Shin 等人（2018），以及 Shang 等人（2018、2019）。关于将 GNBC 用于尖锐界面问题的全面讨论，可参考 Fullana 等人（2024）的研究。

为了帮助维持接触线处的动量守恒，我们采用“虚拟界面”的方法，即将分隔气泡和液体的界面延伸至固体内部。这种方法是 Fujita 等人（2013 年、2015 年）在使用水平集方法模拟漂浮在自由表面上的颗粒时引入的。虚拟界面的表面张力与真实界面相同，但其在固体内的斜率由静态接触角确定。原则上，界面的其余部分可以是任意形状，但为了使施加在固体上的力均匀分布，让界面尽可能平滑似乎是合理的。例如，我们可以使曲率的平方最小化。对于二维流动，这相当于最小化以下离散形式的

我们通过一个三步迭代过程来实现这一点。首先，我们对相连点的位置进行平均；然后，对这些平均值再次进行平均；最后，用两倍的第一次平均值减去第二次平均值，从而计算出最终位置。对于三维流动，我们对与我们所考虑的前沿点直接相连的所有点进行平均。尽管这并不完全等同于最小化总曲率，但它能得到一个足够平滑的表面。有关该方法的详细描述，包括收敛性测试，请参阅 Zeng 等人（2024）的研究。

使用虚拟界面将气泡界面延伸至颗粒内部意味着，当颗粒附着在气泡上时，会存在一个和重叠的区域。因此，我们通过以下公式来设定局部属性，如密度和黏度：

固体指示函数会覆盖气泡指示函数，因此，即便对于处于气泡内部的那部分固体，颗粒的属性也总是能被正确设定。

3、问题设置

我们研究了单个大气泡和多个较小的疏水性颗粒在液体槽中的运动情况。模拟在一个矩形区域内进行，该区域在水平方向采用周期性边界条件，顶部和底部采用壁面边界条件。垂直尺寸是水平尺寸的四倍。初始时，使用随机数生成器为颗粒质心坐标赋值，将颗粒插入到区域中包含颗粒的部分，但会舍弃与先前插入的颗粒重叠的颗粒。气泡被放置在静止液体中靠近区域底部的位置，模拟会一直运行到气泡接近区域顶部的时刻。重力方向向下，使得气泡上升，颗粒相对于液体下落。每个固体颗粒的网格分辨率约为 10 个网格点/直径，这是通过网格细化研究确定的，详见附录。具体的物理和数值参数列于表 1。颗粒相对于未变形气泡的相对尺寸为

，颗粒总体积分数为

其中Lx、Ly和Lz是区域的尺寸。这些参数使得上升气泡的雷诺数达到适中的约 80。为了确定每个控制参数对颗粒捕获率的影响，我们每次改变一个参数，同时保持其他参数固定，以此进行测试。

表1

在黏性液体中上升的气泡运动可用阿基米德数和Eötvös数来表征，其定义为

其中。莫顿数，有时会用它替代阿基米德数。在所有工况下，阿基米德数Ar=3110.4，但Eötvös数和Morton数会随表面张力而变化，如表II所示。除了对附着在气泡上的疏水性颗粒进行模拟外，我们还针对不附着在气泡上的亲水性颗粒重复进行了部分模拟。

表2

A. 一个简单模型

对于此处所研究的简单情况，我们可以通过相对基础的建模来深入了解其演变过程。许多作者都提出过此类模型，比如 Bloom 和 Heindel（1997 年）分析了附着有累积颗粒的气泡所受的力；Yoon 和 Mao（1996 年）给出了气泡上颗粒累积的一阶速率方程：，其中 P 是收集概率，Vg是空气的表观速度。更近一些，Gharai 和 Venugopal（2016 年）回顾了浮选动力学模型，其中颗粒的累积由给出，这里α是速率常数，m是过程的阶数。我们的简单模型借鉴了这些论文中提出的思路。

图 1. 建模示意图。在时间间隔Δt 内，只有气泡上方圆柱体内的固体颗粒能够被捕获。图中绘制了几条流线以供参考。

图 1 的示意图展示了一个气泡穿过颗粒床上升的过程，在这个过程中，气泡所经过的圆柱形体积内的部分颗粒最终会附着在气泡上。气泡捕获颗粒的速率直接与固体颗粒数密度（n）的某一比例、气泡的横截面积（Ab）以及相对速度（Dub）的乘积成正比。气泡的上升速度可以通过多种方式来测量，比如相对于固定区域、相对于流体，或者如我们此处所采用的，相对于颗粒来测量。颗粒与气泡相互靠近的速率决定了颗粒与气泡的碰撞情况，因此我们在此使用相对速度（Dub）。所以，颗粒在气泡上的累积速率为

这里，我们运用了

以及

。Ca是一个“捕获系数”，它描述了气泡正上方圆柱体内被捕获的颗粒比例。对于那些不会移开的小而重的颗粒，我们预计Ca≈1；对于随流体一起移动的非常轻的颗粒，我们预计Ca<1；而如果颗粒又大又重，我们可能会有Ca>1，这是因为一个球形气泡穿过球形颗粒床上升时的有效扫掠面积是气泡半径与颗粒半径之和。这对于大颗粒而言很重要，但鉴于我们分析中存在的其他简化情况，以及我们所研究的气泡是可变形的，我们选择仅使用等效球形气泡的横截面积。不过，在所有情况下，我们都预期该系数接近 1。通过以下方式对时间进行无量纲化

我们可以将方程 (11) 改写为

其中雷诺数定义为

。

我们假设气泡的上升速度处于准稳态，并且颗粒的重量会降低气泡的浮力。通过对浮力和阻力进行平衡，我们得到：

其中，

是阻力系数，

。对单个气泡在不同重力条件的连续流体中上升进行的模拟显示，阻力系数可以近似为与雷诺数成反比。因此，对于带有颗粒的气泡，我们采用一个简单的模型，即

，其中 k是一个常数（称为“阻力定律常数”）。将此方程代入式 (14)，我们得到：

这里，β是单个颗粒的“浮力”与气泡的“浮力”之比。将式 (15) 代入式 (13)，可得：

这个方程很容易进行积分，从而得到：

其中，C是积分常数。最后，将式 (17) 代入式 (15)，结果为：

式 (17) 和式 (18) 中的未知常数分别是阻力定律常数 (k)、捕获系数 (Ca) 和积分常数 (C)。对于每个具体情况，这些常数通过将数据与方程进行拟合来确定。

在气泡上累积的固体颗粒会使其上升速度减慢。当附着颗粒的重量与气泡的浮力达到平衡时，也就是当

时，气泡将停止移动。我们还可以通过以下方法估算能够覆盖一个气泡的最大颗粒数量：用气泡的表面积

除以单个颗粒的最大横截面积

，再乘以堆积效率Sp，由此可得

。通常情况下，球体上圆形的堆积效率低于平面上的堆积效率，但随着圆形尺寸变小，堆积效率会提高。对于无限多个小圆形而言，其堆积效率上限约为 0.9（Clare 和Kepert，1991 年）。

四、结果

我们首先详细研究一个案例的非定常演化过程，然后改变控制参数。对于此处研究的参数范围，结果具有相似性，因此第一部分展示的随时间变化的结果也代表了其他案例的情况。

A. 时间演变

图2展示了在埃特沃什数Eo=1、阿基米德数Ar=3110.4、c=1/4、静态接触角以及滑移系数

的条件下，一个气泡（黄色）和五百个颗粒

在三个不同时刻的情况。图2（a）呈现的是初始状态；图2（b）中，气泡已上升至区域略超一半的位置；图2（c）里，气泡接近区域顶部。随着气泡上升，颗粒下落，聚集在底部，在顶部形成一个无颗粒区。气泡正前方的许多颗粒与气泡碰撞并附着其上。附着的颗粒会较快地被卷入气泡，直至接触角接近静态接触角，随后沿气泡表面滑动，并在气泡后部附近聚集。为了表明颗粒对流场造成的扰动，我们在后平面绘制涡量

。尽管在图中并非总是清晰可见，但高涡量区域与颗粒所处位置紧密相关。下落的颗粒挡住了我们观察气泡的视线，因此在图3中，我们展示了气泡上升至区域约一半高度时，从两个不同角度看到的气泡及附着在上面的颗粒。图3（a）为侧面视角，图3（b）是从底部观察气泡的视角。侧面视角表明，由于颗粒向气泡后部滑动，气泡顶部没有颗粒，但气泡下半部分被颗粒覆盖。底部视角显示，颗粒分布相对均匀。我们让颗粒保持略微分离的状态，因为对于相互接触的颗粒，很难维持规则的表面网格。随着网格细化，颗粒会靠得更近。

图2. 对于埃特沃什数Eo=1、阿基米德数Ar=3110.4、α= 2.8%、c=1/4、静态接触角θo=100，以及滑移系数Cslip=0.2的情况，展示1个气泡（黄色）和500个疏水性颗粒（灰色）在不同时刻的状态：(a) τ= 0，(b) τ= 7，(c) τ= 14 。涡量显示在区域后部的一个平面上。

图 3. 图 2 在 τ= 7时的细节：(a) 侧面视图，展示了穿过气泡中部平面的速度矢量和涡量；(b) 底部视图。

颗粒最初随机分布在区域内。为确保结果对初始颗粒分布不敏感，我们使用四种不同（但仍为随机）的初始分布进行了模拟。在图4中，我们绘制了（a）气泡的上升速度和（b）附着在气泡上的颗粒数量随时间的变化，展示了四种初始条件下的情况。雷诺数最初像单个气泡那样增加，但并未达到稳定状态，而是随着颗粒附着在气泡上立即开始下降。由于颗粒要么附着在气泡上，要么没有附着，图（b）中显示附着颗粒数量的曲线是不连续的，但很明显，随着更多颗粒附着在气泡上且气泡速度减慢，颗粒的累积速率也会减慢。不同初始颗粒分布的结果表明，对于相同的控制参数，结果几乎相同，曲线的整体形状保持不变。因此，在本文其余的模拟中，我们仅研究一种初始条件。对于这个特定案例,

和

的计算结果均为57.6（如果我们使用），所以在模拟接近结束时，气泡聚集的颗粒数量接近理论最大值的一半。一般来说，我们预期

和

是不同的。

图4. 对于对于埃特沃什数Eo=1、阿基米德数Ar=3110.4、α= 2.8%、c=1/4、静态接触角θo=100，以及滑移系数Cslip=0.2的情况，展示四种不同初始颗粒分布下：(a) 雷诺数随时间的变化；(b) 附着在气泡上的颗粒数量随时间的变化。

为了研究系统行为如何随控制参数变化，我们使用由方程（17）和（18）给出的简单模型。该模型包含三个未知系数：捕获系数（Ca），表示气泡前方最终附着在气泡上的颗粒比例；阻力定律常数（k），用于近似阻力系数与雷诺数的关系；以及积分常数（C）。我们定义：

并将方程（17）和（18）改写为：

然后通过拟合来确定这些系数。该模型并未描述雷诺数最初的增加情况，因此我们从通过直观观察估计的雷诺数对时间的二阶导数变号的大致时刻开始拟合数据，一直到气泡到达含颗粒区域顶部的时刻。系数C 、k和Ca是使用MATLAB的“lsqnonlin”函数拟合方程（20）得到的，该函数通过求解非线性最小二乘方程来最小化一个组合损失函数。利用这些拟合值，再依次通过方程（19）计算C、k和Ca。图5表明，尽管在开始和结束时存在一些细微偏差，但拟合结果与数据的吻合度相当不错。我们对其他案例的结果也进行了拟合（稍后讨论），通常发现这个简单模型能够很好地描述其演变过程。由于一般来说每个案例的拟合系数都不同，观察这些系数如何随控制参数变化，是一种简洁的方式来了解系统演变如何随控制参数改变。

图5. 方程（17）和（18）对数值结果的拟合情况：(a) 雷诺数；(b) 附着在气泡上的颗粒数量（Na）。

我们收集了其他数据来量化这一演变过程。在图6中，我们展示了在运行结束时，针对四种不同初始条件，计算包含悬浮颗粒的区域部分液体中的平均湍动能和平均耗散率。这意味着要排除底部上方和顶部下方高度约为四个固体颗粒直径的两个区域，因为在运行结束时，底部区域包含静止颗粒，而顶部区域没有颗粒。最初，随着气泡开始上升、颗粒开始下落，动能和耗散率会增加，但随着气泡速度减慢，动能增加的速度放缓，耗散率也会略有下降。随着气泡上升并聚集颗粒，气泡前方和后方的流动条件有所不同。为了研究这些变化，我们在图7中绘制了气泡前方和后方两个区域内的局部平均液体湍动能（a）和局部平均湍流耗散率（b）。这两个区域高度为两个气泡直径，分别位于气泡前方和后方半个气泡直径处。通过将每个区域的局部值除以整个区域（不包括顶部和底部附近的小区域）的平均值进行缩放，我们发现这两个值在时间上几乎保持恒定。该图表明，当气泡经过时，局部湍动能和局部湍流耗散率显著增加。显然，我们的数据仅限于气泡至少上升了两个气泡直径且尚未过于接近顶部的时间段。在此，我们仅展示τ= 6 - 9.6这一时间段的数据。

图6. （a）平均液体湍动能（Ek）和（b）平均耗散率（ϵ）随时间的变化，数据来自整个区域，但不包括顶部和底部附近的小区域。

图7. （a）气泡前方（Ekt）和后方（Ekb）的局部平均湍动能，以及（b）气泡前方（ϵt）和后方（ϵb）的局部平均湍流耗散率，均通过整个区域（不包括顶部和底部区域）的平均动能进行归一化。

为研究附着在气泡上的颗粒如何改变气泡后方的颗粒分布，图8（a）展示了在τ= 7-14时间间隔内，以与气泡固定的参考系下颗粒质心的运动轨迹。左边的颗粒是疏水性的，其中一些会附着在气泡上；右边的颗粒是亲水性的，不会附着。附着的疏水性颗粒部分被卷入气泡，其轨迹因此终止在气泡表面，而亲水性颗粒则从气泡上反弹并滑过。很明显，气泡正前方的大多数疏水性颗粒最终附着在了气泡上，因此与右边的尾流相比，左边的尾流中颗粒相对较少。注意，由于气泡是完全三维且对称的，这些轨迹会以不同的方位角经过气泡，所以离气泡轴越远，我们看到的轨迹就越多。图8（b）展示了在高度为两个气泡直径、分别位于气泡前方和后方半个气泡直径处的区域之间，颗粒到气泡质心的相应平均水平距离的归一化变化。由于平均体积较小，即我们对相对较少的颗粒进行平均，两条曲线波动都很明显。然而，很清楚的是，亲水性颗粒到气泡对称轴的平均距离没有变化，而疏水性颗粒的平均距离增加了。

图8. （a）对于Eo=1，固体颗粒相对于气泡的轨迹。左图：疏水性颗粒；右图：亲水性颗粒。（b）颗粒到气泡对称轴的平均水平距离的变化，通过初始值进行归一化。

B. 控制参数的影响

在本节中，我们研究结果如何依赖于各种控制参数。其演变过程与第四节A部分所研究的情况类似，尽管细节有所不同。在图9中，我们展示了在τ= 7时的两种情况，一种是具有更大可变形性的气泡（Eo=4），另一种是更小的颗粒(c=1/5)，但其他所有参数与第四节A部分保持一致。

图9. 侧视图：(a) 埃特沃什数Eo=4（Eötvös）变为时；(b) 尺寸比变为c=1/5时。两种情况均对应τ= 7。

在图10中，针对三种改变一个参数而其他参数保持不变的情况，绘制了气泡雷诺数（左列）和附着在气泡上的颗粒数量（右列）随时间的变化。图10（a）展示了通过改变表面张力进而改变埃特沃什数（Eötvös number）所带来的可变形性影响。作为参考，我们还绘制了气泡穿过不粘附气泡的亲水性颗粒悬浮液时的雷诺数。尽管由于颗粒的存在会出现一些波动，但平均而言，气泡以近乎稳定的雷诺数向上移动。虽然气泡穿过亲水性颗粒悬浮液时的雷诺数取决于埃特沃什数，但在后期，当颗粒为疏水性时，基本上没有差异。当亲水性颗粒与可变形性更强的气泡碰撞时，气泡更大程度的变形会导致稳态雷诺数显著降低。气泡捕获疏水性颗粒后，随着接触线移动，它会平稳地吞没这些颗粒。对比图3（a）和图9（a），不同埃特沃什数下气泡的形状略有变化，这减弱了可变形性对雷诺数的影响。

图10. 左列：气泡雷诺数随时间的变化。右列：附着在气泡上的颗粒数量随时间的变化。在（a）中，我们研究埃特沃什数Eo变化带来的影响；在（b）中，改变的是体积分数；在（c）中，改变的是颗粒的相对大小。其他参数与第四节A部分的参考案例相同。在（a）中，虚线表示气泡在亲水性颗粒悬浮液中上升的情况。

对于此处研究的参数范围，如图10（b）所示，改变固体体积分数对演化过程有显著影响。我们通过增加颗粒数量但保持其尺寸不变来提高固体体积分数。方程（17）和（18）中指数函数的参数与体积分数（α）直接成正比，这表明较高的体积分数会导致颗粒在气泡上更快地累积，且雷诺数下降得更快。对于最低体积分数（α=1.4%），气泡捕获的颗粒较少，上升速度更快，所以当气泡到达颗粒床顶部时，我们提前终止了模拟。在图10（b）中，我们还给出了最高和最低体积分数的拟合曲线。尽管简单模型[方程（17）和（18）]所给出的拟合能够准确描述此处进行的大多数模拟中的演化情况，但对于较大的体积分数，误差会更高。特别是，如果我们仅使用早期数据来确定拟合系数，后期的累积速度比预测的要慢。

在泡沫浮选过程中，颗粒大小是一个重要因素。在图10（c）中，我们通过未变形气泡直径与颗粒直径之比c=Ds/Db来改变颗粒大小。由于需要确保有足够的网格点来准确解析颗粒运动，颗粒大小的减小程度受到限制。不过，附录中的网格细化研究表明，最小颗粒的整体演化过程能够得到较为准确的捕捉（尽管可能不如较大颗粒那么精确）。由于体积分数保持不变，所以小颗粒数量更多，大颗粒数量更少，因此我们可以看到，随着颗粒尺寸减小，气泡上附着的颗粒数量会增加。然而，颗粒越小，单个颗粒的质量就越小，每个颗粒对气泡上升速度的降低作用也就越小。在图11中，我们重新绘制了图10（c）右侧的图像，但将气泡上的颗粒数量乘以它们的质量，从而绘制出气泡收集到的质量（Ma）随时间的变化。最初，颗粒大小对累积的总质量影响不大，但在后期，我们可以看到显著差异。尽管随着颗粒尺寸减小，气泡上的颗粒数量增加，但累积质量却显著下降。对于给定质量，较小的颗粒会覆盖气泡表面更大的部分，从而阻碍其他颗粒与气泡碰撞并附着。

图11展示了不同颗粒尺寸下，附着在气泡上的固体总质量（颗粒数量乘以单个颗粒重量）随时间的变化情况。

在图12（a）中，我们对比了在τ= 7-14时间间隔内，体积分数α=1.4%（左图）和α=5.5%（右图）时颗粒的运动轨迹。结果表明，在较高体积分数下，已附着在气泡上的颗粒会阻止其他颗粒附着[图12（b）]，因为颗粒累积速度足够快，当另一个颗粒与气泡碰撞时，新附着的颗粒仍在气泡表面。另一方面，对于较低体积分数，颗粒有时间立即“让开”，从而为后续接近的颗粒呈现出无颗粒附着的气泡表面。

图12. （a）体积分数α=1.4%（左图）和α=5.5%（右图）时颗粒的运动轨迹。（b）α=5.5%体积分数，在τ= 11时刻的侧视图。

图13以简洁的形式展示了演化过程对控制参数的依赖关系。在此图中，我们绘制了方程（17）和（18）的归一化拟合系数与控制参数的关系。归一化处理是将第四节A部分所呈现案例的系数（Ca=0.8；k=55.7，以及C=0.95）作为标准案例，用其他案例的系数除以这些标准值。这样我们就能在一张图中展示所有案例的拟合参数。图13（a）展示了不同埃特沃什数下的归一化拟合系数。随着气泡可变形性增强，阻力定律常数（k）和积分常数（C）略有增加，但捕获系数（Ca)在Eo=1时，随着气泡可变形性稍有增强先略有上升，随后下降。对于可变形性更强的气泡，捕获系数下降是因为附着在气泡底部的颗粒导致气泡变形，这些颗粒将气泡表面向下拉，使气泡呈扁长形，从而减小了扫掠面积。图13（b）展示了体积分数的影响。很明显，尽管附着在气泡上的颗粒数量随体积分数增加，但捕获系数会下降，因为更快的累积速度导致气泡上已有的颗粒会阻止新到达颗粒的附着。阻力定律常数基本保持不变，但积分常数略（C）有下降，因为最初捕获的颗粒更多，导致该值变小。图13（c）研究了颗粒相对于气泡的大小的影响。相对于气泡的颗粒大小对积分常数基本没有影响，但阻力定律常数随颗粒大小略有下降。另一方面，捕获系数确实依赖于颗粒大小，对于可能沿流线运动的小颗粒，捕获系数会下降；而对于即使在标称扫掠区域之外也可能与气泡碰撞的大颗粒，捕获系数会上升。

图13. 捕获系数Ca、阻力定律常数k和积分常数C的归一化值与控制参数的关系。(a) 埃特沃什数Eo；(b) 体积分数α；(c) 尺寸比c。此处Ca=0.8，k=55.7，C=0.95 。

我们还研究了静态接触角

和滑移系数

的影响，并在图14中展示了改变这些参数如何影响拟合系数。归一化方法与图13相同，采用标准案例的值。气 - 液界面与固体表面的接触角对泡沫浮选的选择性至关重要。然而，一旦颗粒附着，接触角的主要影响在于确定颗粒嵌入气泡的深度（随着接触角增大，颗粒在气泡内的深度越深）。对于此处研究的取值范围，我们发现基本上不受静态接触角

或滑移速度

的影响。

图14. 捕获系数Ca、阻力定律常数kC和积分常数的归一化值与控制参数的关系。(a) 静态接触角θo。(b) 滑移系数Cslip。此处，Ca=0.8，k=55.7且C=0.95。

5、结论

浮力驱动的气泡对疏水性颗粒的捕获是一个复杂问题，众多物理过程会影响最终的选择性。这些过程包括流体湍流、气泡 - 气泡相互作用、气泡破裂与聚并，以及表面活性化学物质的存在。尽管这些过程之间的相互作用错综复杂，但大多数基本过程已得到相对较好的理解。因此，能够精确捕捉所有连续空间和时间尺度的全解析模拟，应有助于极大地增进我们对这一完整问题的理解。在此，我们聚焦于一个相对简单的设定，仅考虑单个气泡在固体颗粒悬浮液中移动的情况，并研究了各种控制参数的影响。为简化计算，我们将雷诺数维持在相对适中的水平，采用非常简单的附着和脱离标准，且仅研究相对较大的颗粒。为了理解动力学如何随控制参数变化，我们使用了一个带有自由参数的简单模型，并针对每组数据进行拟合。该模型能够较好地描述颗粒在气泡上的初始累积情况，但在后期，已附着的颗粒会干扰新颗粒的附着，实际的附着过程比模型所描述的更为复杂。即便如此，研究拟合系数如何随控制参数变化，为我们审视这些参数的影响提供了一种简洁的方式。

结果表明，在本文所考察的参数范围内，气泡的可变形性对结果的影响相对较小，不过在较高的埃特沃什数下，捕获系数会有所降低。虽然较低的表面张力会使纯净气泡变得更扁平，但颗粒在气泡底部的聚集会使其在流动方向上拉长，从而抵消了因表面张力降低而产生的 “扁平化” 效应。尺寸比有一定影响，对小颗粒而言会降低捕获系数，而对较大颗粒则会提高捕获系数。

影响最为显著的参数是体积分数。随着气泡负载的颗粒增多，附着速率会下降，且在较高体积分数下，这种下降会更快。由于我们的简单模型仅考虑了捕获颗粒的重量，而未考虑它们对降低附着速率的作用，因此发现在本文模拟的时间段内，不仅捕获系数随体积分数降低，拟合的质量也会下降。最高体积分数的结果表明，即使对于较低体积分数，如果跟踪气泡的时间足够长，该模型的适用性最终也会变差。在所考察的范围内，静态接触角和滑移速度的影响都较小。尽管完整的问题要复杂得多，但我们认为本文的研究是重要的第一步，其结果将为进一步的研究提供基准。

翻译转自：《Physics of Fluids》：“A numerical study of the capture of hydrophobic particles by a buoyant bubble”

来源：多相流在线

Maxwell 碰撞非线性多相流化学群体平衡 MATLAB UG UM 理论材料

著作权归作者所有，欢迎分享，未经许可，不得转载

首次发布时间：2025-03-13

水务污水处理仿真技术创新与应用

免费

还没有评论