积鼎流体仿真VirtualFlow：基于粗粒连续接触角模型的微流控装置充注模拟

为了设计和开发新的微流控装置，必须提供可靠的建模和仿真工具。需要对传统计算流体动力学进行扩展，特别是多相流体动力学的模拟能力。VirtualFlow提供一种新的动态接触角模型，该模型不需要确定接触角或接触线速度。自由表面的捕捉采用level-set方法。对VirtualFlow中该模型进行了单元问题的测试，如接触线的松弛到平衡。将建立的模型应用于典型微器件的流体填充问题，以验证VirtualFlow作为设计CFD工具的实用性。

1 绪论

近年来，微流控技术发展迅速。然而，在建模和仿真工具方面很难有成熟的设计支持。其核心的难点在于需要对传统计算流体动力学(CFD)工具进行许多扩展，例如表面力、三相接触线动力学、非均相化学反应、表面活性剂等小规模的物理效应等。

具有移动接触线的不可压缩两相流在各种应用中很常见，例如涂层和生物过程。模拟这种流动的困难在于两种流体的N-S方程结合无滑移边界条件会导致接触线上存在剪切应力奇点。此外，还需要指定自由表面与固体表面形成的接触角。

在过去几十年里，人们曾多次尝试在连续流体动力学框架内对接触线动力学进行建模。然而，这些模型几乎很少考虑到作用于三线区域的作用力。VirtualFlow中提供的模型是基于接触线的分子动力学建立的，从而不需要给定接触角。在基于杨氏定律的动量方程中加入不平衡杨氏力，而非指定接触角。采用滑移长度对接触线上的应力奇点进行正则化，并对直至滑移长度的对数奇点进行积分，从而解决这一问题。

将模型与Spelt提出的单元模拟进行了比较，并使用其分析具有不同表面特性的液体对原型微器件的填充过程。气泡夹带是微流控器件设计中的一个常见问题。通过精细设计器件的几何形状、控制填充过程和选择材料特性(疏水或亲水性)，可以避免这些气泡的形成。在这一过程中，可以使用可靠的仿真工具对几何形状进行设计修改，以确保设计是有效可行的。

2 仿真框架

2.1 VirtualFlow微流体计算模块

VirtualFlow是一个多物理场有限体积法多相流CFD软件。VirtualFlow采用单流体方法，即双流体流动被视为一种流体，其材料特性根据区分两种流体之间的边界或界面函数而变化。具体来说，Level-Set和VOF接口跟踪方法(ITM)都可以用于跟踪不断变化的自由界面。

2.2接口跟踪

当分离两种流体的界面的确切形状不确定时，有时可以采用平均双流体方法，其中每个相需要单独的守恒方程，并具有恰当的界面作用力。当需要精确地识别界面时，例如在气泡、液滴或液体射流破碎过程中，可能会调用ITM。关键是使用一组具有可变材料性质和表面力的单流体守恒方程。这种方法提供了一种比双流体方法更准确的模型，同时最小化了建模假设。

2.3输运方程

用单流体形式表示的不可压缩流体动力学方程采用如下形式:

其中动量方程(方程(2))中的RHS项分别为体力、表面张力（方程(4)）及其接触线壁面贡献。不考虑相变时，为牛顿应力张量。在level-set方法中，两种流体之间的界面用连续函数φ表示，φ表示到界面的距离，该距离在一侧为正，在另一侧为负。通过这种方式，可以识别两种流体，并且将物理自由界面的位置与0等值面相关联。Level-set控制方程由方程3给出：

材料性质，如密度和粘度是基于φ的局部物理量，随时间变化，并使用改进的Heaviside函数在界面上做平滑处理。此外，φ是跨越自由界面的连续函数，有助于确定界面的法向量，从而确定定义表面张力所需的表面曲率，如方程(4)所示

其中γ为表面张力，kappa为界面曲率，为以界面为中心的光滑狄拉克函数。方程在一次迭代之后，level-set函数不再是与界面的有符号距离。为了恢复其在界面附近的正确分布，需要重新赋值距离，将下式积分为稳态:

其中sgn(x)是sgum函数。在VirtualFlow中方程(6)在方程(3)的每个迭代步后，使用非振荡(WENO)三阶格式求解。

2.4动态接触角模型

VirtualFlow采用与三线相关物理力的润湿动力学数值处理方法，可以提供对润湿动力学的物理充分描述，从而消除了指定接触角的任何特定边界条件的需要。在本研究中使用的三线力是基于考虑界面自由能的。因此，它只包含两个参数:流体之间的界面张力γ和平衡接触角。

其中，是瞬时动态接触角，是狄拉克函数。

3 结果

3.1松弛到平衡

在这里，我们考虑Spelt（Spelt P. D. M., J. Comp. Phys., 207, 389-404 (2005)）提出的两种情况，即液滴浸入不同的流体中，并附着在矩形域(2 × 1，由2N × N离散)的长边界上。在这些情况下使用的流体性质列于表1中。在Case1中，密度相等(= 1)，而粘度不同(液滴为4.95×10−2，周围的相为4.95 × 10−3)。然而，在case 2中液滴与周围流体的密度比为1:20，粘度恒定为4.95 × 10−2。两个Case的表面张力分别为γ = 0.11和2.21，网格尺寸为128x64。

(a) t=0s

(b) t=3.7s

(c) t=19.8s

图1:非润湿液体恢复平衡; = 1200

图2:接触角随时间的变化。图例:圆形-情况I;三角形-情况II

首先，液滴呈现圆帽形(半径为1.66)，图1和图3中的蓝色相的接触角分别为300度和150度，这与平衡接触角=120度和600度不同。当前模型与Spelt模型的主要区别在于，后一种情况需要额外规定接触线速度。

(a) t=0s

(b) t=0.6s

(c) t=1.0s

(d) t=3.6s

图3:润湿液恢复平衡; = 60°

图1和图3分别显示了case 1和2情况下平衡时的状态。接触角在流动趋于平衡时的变化如图2所示。对于case 2，由于较高的表面张力以及由于较高的外相密度而导致其可以更快地恢复平衡。

3.2微流控装置的灌装

为了说明VirtualFlow在实际应用中的实用性，对=70°、90°和110°三种不同液体的微型泵原型装置的充液过程进行了仿真。图4显示了最终结果，该结果表明，对于疏水性情况，必须特别小心，才能设计出没有气泡截留的设计。

(a) 几何与网格

(b) 接触角=70°

(c) 接触角=90°

(d) 接触角=110°

图4:微器件原型液充

图中蓝色代表液体，灰色代表空气，实线表示相间界面，白色表示固体区域

4 总结

VirtualFlow中搭载全新的动态接触角模型，该模型通过在动量方程中增加力项来消除指定接触角或接触线速度的要求。本文对液滴恢复平衡状态的过程进行了模拟，并与前人所做的可靠数据进行了对比。对具有不同润湿特性的液体填充微器件的仿真结果表明，VirtualFlow中搭载的基于粗粒连续接触角模型是一种有效的设计工具。