COMSOL 软体建模教程：如何模拟自由液面 (一)
COMSOL Multiphysics? 软体提供了四种可用于模拟自由液面的方法：水平集、相场、动网格和稳态自由表面。作为系列的第一节，我们将讨论水平集和相场法，这两种基于场的方法几乎可以描述任何类型的自由液面。在第二节中，我们计划将本文的求解结果与通过动网格介面获得的结果进行比较。
何为水平集和相场法？
水平集和相场法都是基于场的方法，这类方法将自由液面表征为水平集或相场函数的等值面。自由液面对应固定网格框架下的液体和气体之间的相分界面。
下图为管道内两颗液滴的表面，摘自附加产品「微流体模块」的「案例库」所提供的液滴破碎模型。从这张图中可以看出，尽管液滴的表面非常明显，但液滴周围的单元并没有贴合到液滴表面上。
水平集和相场函数都是由纳维-斯托克斯方程计算的速度矢量进行对流传输的。在水平集和相场法中，对应公式为：
需要注意的是，水平集和相场函数都使用了 Φ。二者的不同在于方程右侧的 F。在初始水平集方法中 F = 0，因此得到纯对流传输方程。然而当 F = 0 时，数值解不仅不稳定，而且大部分情况下实用性很小。所以为了保持相界面清晰，我们在水平集方法的 F 中添加了高阶导数项 Φ 。
在相场法中，F 代表设法将系统的自由能最小化的项。此项也引入了高阶导数 Φ。实际上，相场方程中的源项中包含了四阶项。这意味著，出于实用性考虑，方程经常被分解为两个方程，与此同时，辅助因变数被定义为 Φ 的二阶导数函数形式。COMSOL Multiphysics 中也采取了这种做法。
两种方法均将自由液面的表面张力引入到纳维-斯托克斯方程的源项中。水平集方法利用表征自由边界的水平集等值面的曲率来描述表面张力。相场法根据化学势计算出表面张力对纳维-斯托克斯方程的源项贡献，因为正是化学势导致了界面附近产生的表面张力和相场函数梯度的产生。
在给定的水平集或相场函数值范围内，即自由表面的数值范围内，流体特性从液体平滑地过渡到气体。
穿过自由表面，水平集函数 Φ 在 0 和 1 之间变化，在两种流体内部分别为常数 0 或 1。比如在液相中为 0，气相中为 1。至于自由液体面处，也就是液体和气体之间的分界面上，对应的水平集函数值为 Φ = 0.5。根据下方公式，密度是水平集函数的函数：
动力粘度也是如此：
从这两个方程中我们可以看出，当 Φ = 0 时，可得到流体 1（比如液体）的属性，当 Φ = 1 时，可得到流体 2（比如气体）的属性。由此，流体属性在整个自由表面上实现了平滑过渡，像水平集函数一样平滑。
相场函数 Φ 在 -1 和 1 之间变化，当 Φ = 0 时的等值面为自由液面。粘度和密度的计算是通过与水平集法相似的方式完成的，只不过因为相场函数在 -1 和 1 之间变化，所以表达式不同。
水平集与相场介面的设置
水平集和相场介面中包含若干个参数，为了使计算取得最优效果，我们需要对这些参数进行调整。参数的值取决于被模拟的系统和方程的数值离散。几何结构的尺寸（长宽高）、流体的特性、壁的润湿条件、初始条件、工作条件以及网格的单元尺寸都会对模型设置选项产生影响。
水平集介面
我们从水平集方法的设置开始。用户界面中的重新初始化参数 γ 可以确保水平集函数中的梯度随时间推移逐渐集中到自由表面上，也就是下方示例中水和空气之间的交界面。它没有设定表面厚度，但能够确保水平集函数的变化在厚度范围之内。如果该参数值设置过低，水平集函数的变化可能被截留在一种流体域内，然后凭空生成液-气界面。过高的值往往造成时间步过小、计算时间过长。
ε 表示的界面厚度参数，它理解起来更加直观，唯一的作用是控制水平集函数发生变化的区域的厚度。过厚会模糊自由表面的形状。虽然这让问题变得更容易求解，但是损失了自由表面的形状精度。较小的 ε 值有利于生成清晰的表面，但这也要求网格具有相应的解析度。设定小于网格单元尺寸的值是没有意义的，因为这会使得界面变得不规则，而无法解析的表面张力可能会导致速度失真。ε 的设定值应与表面附近的最大单元大约相同。
相场介面
接下来是相场法，界面厚度参数 ε 与水平集方法基本一致，前文中的论证同样适用于此设定值。
相场法的迁移调整参数 χ 在某种程度上决定了相场的扩散系数。它的值必须足够大，才能保证相场方程稳定，但是为了保证表面尖锐，同时也要足够小。合理值与表面的速度成正比，与表面张力系数成反比。
这意味著一旦确定了某组工作条件的 χ 值，就可以利用上述关系为一组新的条件设定 χ 的值。
自由液面建模
通过下图中的示例问题，我们将演示如何在 COMSOL Multiphysics 中利用水平集和相场法进行自由液面建模。
实心条一半浸没在小型管道内的水中。沿与水面相切的方向来回移动实心条，使水面产生柔和的波浪。为了保持层流状态，管道和实心条的尺寸必须足够小。
在此例中，我们利用动网格功能指定小方形条在水面上反复运动。不过，在水平集和相场法中，网格不会随液体表面移动。
下面是关于在 COMSOL Multiphysics 中创建模型的几点备注：
- 添加重力作为动量方程的源
- 添加参考压力和压力约束点
- 为了方便比较结果，对壁应用 Navier 滑移条件，并使滑移长度等於单元长度
水平集与相场法的模拟结果
下方的结果图显示了 0.07 秒、0.57 秒和 1.0秒的流场和水面形状。两种方法在各个时间点上的流线和水面形状均大致相同。最大速度和水面高度稍有差异，原因可能在于两种方法处理表面张力的不同方式。
再经过一段时间，流线也出现了差别，例如 t = 1.0 s 秒时的回流区。通过相场法得到的液面相对更加平静。水平集方法利用基于水平集函数的梯度获得的表面曲率来计算表面张力。所以，与相场法中更加平滑的力相比，水平集方法得到的表面力更「尖锐」。
我们可以通过动画完整地展示 0~2 秒之间的结果。可以看到，尽管长方形条的运动引起了剧烈搅动，但是表面从未断裂。
水平集和相场法还能计算自由液面上的空气域流场。从图片中可以看出，长方形条的运动在整个相界面上方引起了持续显著变化的流场。
如果增强搅动能使液面运动更剧烈，那么液面可能先破裂再合并，如下方动画所示。这也是水平集与相场法的一个优势：我们能够更简单地处理自由表面的拓扑变化。
尽管水平集和相场法类似，但表面张力的处理对二者的稳定性产生很大的影响，至少在 COMSOL Multiphysics 中存在区别。在处理对表面张力非常敏感的问题时，相场法在求解时间方面比水平集方法表现得更好，原因在于利用水平集方法计算表面曲率时，瞬态求解器需要采用比相场法小得多的时间步长。
在此例中，水平集方法的平均时间步长是相场法的六倍，所以水平集方法需要六倍的计算时间。因此，对于较小尺度自由液面问题和对表面张力极其敏感的层流（例如微流体）问题，相场法通常是更好的选择。
自适应网格细化与自由表面
针对本文探讨的二维示例，在采用了水平集和相场法的整个自由表面区域中，网格足够密集。不过，对于三维模型而言，这种水平的解析度的计算成本太高。一种替代方案是使用自适应网格细化功能，它能够自动根据我们选择的任何函数创建更密集的网格。
举例来说，下方动画展示了基于最大速度梯度的位置（上图）和相场函数的最大梯度（下图）生成的网格自适应结果。绘图展示了气相和水相。需要注意的是，与剪切率相关的自适应网格在气相中生成了细化网格，而水相的细化程度相对较低。在此例中，我们对水相更感兴趣，因此需要通过缩放相场函数来修改自适应网格的误差指示器。
00:06
利用速度梯度作为误差指示器
00:06
利用相场梯度作为误差指示器
动网格与自由表面
本文讨论了两种基于场的自由表面建模方法：相场法和水平集方法。在此系列的第二篇中，我们将使用动网格为自由表面建模，并比较两篇文章介绍的方法。敬请关注！
后续操作
单击下方按钮，了解 COMSOL Multiphysics 的附加「CFD模块」中与流体流动模拟相关的专业功能。
使用 CFD 模块模拟流体流动应用
延伸阅读
- 动手操作：在「案例下载」页面中下载 MPH 文件
- 阅读博客文章，了解水平集与相场法的其他实际应用：
- 模拟掉落在水面上的球体
- 设计可实现材料精确沉积的喷墨列印头
经授权转载自cn.comsol.com/blogs，原作者Ed Fontes
推荐阅读：
