COMSOL两相流水平集仿真：原理、实操与液滴生成全解析

写这篇的主要原因是我发现很多朋友卡在“两相流”这个概念上。一打开COMSOL,看到水平集、相场、动网格、移动网格一堆接口不知道哪个适合自己也不知道选完以后那一大堆默认参数到底是什么意思。尤其做微流控、液滴生成、气泡上升、喷墨打印这类问题时界面追踪的精度直接决定仿真有没有参考价值。这篇我把水平集方法从原理到实操完整捋一遍选一个最经典的T型通道液滴生成案例展开把你最容易踩的坑提前标出来。1. 为什么说水平集方法是两相流仿真里最灵活的界面捕捉器1.1 三种主流界面追踪方法我为什么推荐水平集接触过两相流仿真的人应该知道界面处理大致有三大流派动网格ALE、VOF流体体积法和水平集Level Set再加上一个COMSOL里也很常见的相场法Phase Field。它们的目标是一样的在流场中追踪两种不相溶流体之间的分界面但思路差别非常大。动网格的思路最直观直接把界面当成几何边界网格跟着界面一起变形。好处是界面非常清晰边界条件可以直接施加在界面上坏处是当界面发生拓扑变化时比如液滴从T型通道口被剪切下来、气泡发生破裂或合并网格会出现严重畸变甚至重叠仿真直接崩掉。所以动网格只适合界面变形小、不发生拓扑改变的弱两相流问题。VOF用体积分数表示每个网格单元里某相流体占的体积比例好处是质量守恒性天然很好在商业CFD软件里是绝对主流但坏处是界面需要额外的人工重构或压缩几何精度一般而且对网格方向性敏感。而COMSOL默认的VOF实现里从界面法向速度重构到压力场的耦合处理不好容易出现界面“劈叉”的假象。相场法则引入了一个连续的序参数用双阱势能描述两相自由能界面被看作一个有一定厚度的扩散区域。它的物理背景非常严谨适合研究微纳尺度下界面扩散、润湿、相分离这类本身界面就是模糊的过程。但它的控制方程里包含一个所谓的“迁移率”参数取值大小对结果影响很大而且方程阶数高非线性强算起来偏慢。水平集方法用的是隐式界面描述用一个标量场φ来标记流体界面表示为φ0.5或者φ0取决于你的符号约定的等值面。它最大的优势是天然支持界面拓扑变化液滴断掉、合并、撕裂都没有任何心理负担因为界面不是用网格追出来的而是藏在φ这个标量场里。这一点在微流控液滴生成这种“界面不断断裂”的场景里简直是量身定做。而且水平集方法在COMSOL里耦合多物理场比较自然。密度、粘度、表面张力系数这些参数可以直接写成φ的函数随φ一起输运和更新不需要手动追踪界面位置。做电润湿、热毛细、表面活性剂输运这类问题时直接在水平集方程或者流体方程里加源项就行不用操心“力该加在哪条边界上”。1.2 水平集方程逐项拆解它到底在演化什么COMSOL里层流两相流水平集接口的核心方程如下[ \frac{\partial \phi}{\partial t} \mathbf{u} \cdot abla \phi \gamma abla \cdot \left( \varepsilon abla \phi - \phi(1-\phi) \frac{ abla \phi}{| abla \phi|} \right) ]第一次看到这个方程的人会觉得全是符号。我把每一项翻译成人话。左边第一项是∂φ/∂t表示φ场随时间的变化。这个改变由两部分引起第二部分u·∇φ表示流体流动把界面“吹着走”也就是对流输运它直接来源于流场速度u这就是两相流和NS方程强耦合的桥梁。右边是重头戏γ是重初始化强度ε是控制界面厚度的参数。括号里面有两项第一项ε∇φ是一个各向同性的扩散项它让φ场从尖锐的跳跃变得平滑避免界面处数值振荡第二项是一个“压缩”或“重初始化”项它的作用是把φ的等值线重新拉回到接近阶跃函数的形状保证界面附近φ从0到1的过渡厚度大致稳定在ε量级。如果没有右边这两项φ界面会在流动中慢慢散开变得越来越模糊。如果只有扩散项没有压缩项界面肯定糊了如果只有压缩项没有扩散项数值又不稳定。所以ε和γ这两个参数必须匹配网格和流速来选。从数值角度看水平集法最大的缺点是质量守恒性不如VOF因为在界面区域φ被光滑化了计算密度和粘度时用的是光滑后的加权值导致液滴体积会在长时间仿真中出现可感知的漂移。这点我在第四节会给出具体的排查办法。不过我个人的经验是只要ε取值合理、网格能保证界面区域至少有3到5个单元体积误差在几个百分点以内是完全做得到的。相比VOF那种靠重建算法带来的几何误差水平集的整体稳定性反而更适合多物理场问题。1.3 与多物理场耦合是天然搭档标题里带“多物理场耦合”不是随口号。两相流本身就是一个流场和界面输运的耦合问题而水平集把界面写成一个场变量之后其他物理场可以非常方便地通过φ做域内耦合。举个例子某个物理属性比如粘度μ在两相区分别是μ1和μ2。你可以直接用[ \mu \mu_1 (\mu_2 - \mu_1) \cdot \phi ]把μ定义成φ的连续函数。密度同理。这样Navier-Stokes方程和水平集方程就通过ρφ和μφ双向耦合了。表面张力呢在水平集接口里表面张力是以体积力的形式施加的力的大小包含界面曲率κ、界面法向n和Dirac函数δ。其核心表达式是[ \mathbf{F}_{st} \sigma \kappa \delta(\phi) \mathbf{n} ]注意它被转写成了一个“分布在界面附近区域里的体积力”而不是一条数学意义上没有厚度的边界力。这样做的工程意义是有限元计算不需要在界面处做网格对齐力学框架不必为了一个几何边界而重构。如果你做的是电润湿或者液滴在温度梯度下运动热毛细驱动只需要在NS方程的动量方程里加上电场力或者热毛细力项这些力同样可以写成φ的函数施加。我在做某些液滴驱动器件时就因为这个耦合方式省掉了大量网格重划分的工作。2. COMSOL中的接口选择与参数标定2.1 这一堆两相流接口到底选哪个打开COMSOL的CFD模块你会看到一片“两相流”接口层流两相流水平集、层流两相流相场、层流两相流VOF还有分别对应的湍流版本。选择方式其实很简单主流道和分散相入口的雷诺数都不高用层流版本如果流动已经进入湍流考虑湍流版本如果你主要关心宏观液膜厚度不关心界面微结构和液滴的颈缩过程可以用VOF如果你关心界面处的物理化学过程比如表面活性剂吸附或者液滴在电场下的形貌演化我推荐水平集或相场。我个人在微流控液滴生成、微液滴反应器这类场景里用“层流两相流水平集”接口的次数最多。这个接口底层已经把两相流方程、水平集方程、表面张力、接触角都构建好了。你只需要设置材料属性和边界条件不需要自己手写水平集方程。这点新手容易误解以为水平集算法需要自己从零编码其实COMSOL已经做成了物理场接口你要做的是理解参数、调好数值而不是重复造轮子。2.2 界面厚度epsilon和重初始化强度gamma的玄机epsilon和gamma是水平集接口里最重要的两个参数。这两个值设不好仿真要么界面糊成一片要么振荡发散。先说epsilon它控制界面过渡区的厚度。COMSOL默认值通常是最大网格尺寸的一半但我的习惯是显式地给一个与界面附近网格尺寸匹配的值。实操中我常用的规则是如果你在界面区域打算用h0.002 mm的网格那么epsilon大约取0.002 mm到0.004 mm。界面厚度太大会让液滴看起来像个模糊的水团体积容易虚假增加太小会让界面梯度太陡导致曲率计算产生很大噪声表面张力项不稳定。再说gamma它控制重初始化过程的强度量纲是速度单位是m/s。gamma数值影响φ场被拉回阶跃形状的快慢。一个经验法则是gamma的量级取与特征速度相当或略大。假设入口速度是0.1 m/s那gamma可以先取0.01 m/s到0.1 m/s的量级。gamma太大会让界面过度刚化细小结构被抹掉太小则界面形状保持不住慢慢扩散。在实际调参时我是这么做的先跑一个不带表面张力的纯对流测试只让界面被流场运输观察界面在几十个时间步后是否还保持清晰。如果变模糊调高gamma如果界面挤出了奇怪的振铃波纹调低gamma或减小epsilon。这样把gamma定下来以后再加表面张力做完整仿真能省掉很多崩溃时间。2.3 材料属性与表面张力怎么“耦合”进方程材料定义里最关键的是把密度和粘度写对。水平集接口在材料属性设置面板上通常会自动生成一个“流体1”和“流体2”分别填两相的材料参数。如果你打算手动写变量建议这样定义密度[ ho ho_1 (ho_2 - ho_1) \cdot \phi ]动力粘度[ \mu \mu_1 (\mu_2 - \mu_1) \cdot \phi ]这里的φ在COMSOL中通常被命名为ls.phi在有些版本叫ls.phip。注意接口内部有些辅助变量会叫“水平集变量”之类的名字你需要先在“变量”中确认接口生成的变量名否则表达式求值会报错。表面张力系数的单位是N/m。COMSOL界面里会让你填σ通常直接填实测或文献值。但这里有个工程经验微流控实验中由于表面活性剂、微尺度效应有效表面张力往往比宏观值小很多因此在仿真里使用宏观表面张力常常导致液滴生成实验观测到的尺寸和仿真结果对不上。我在做微通道液滴时前期经常把σ作为扫描参数从0.005 N/m扫到0.03 N/m看哪个量级的结果能与实验吻合。扫描完以后你会发现水平集接口的灵敏度分析其实很直观很适合做这种参数化设计。另外如果是热毛细、电润湿这类问题σ就不能设成常数了需要在界面处写上σ(φ,T)之类的表达式甚至把Marangoni力加进去。这时水平集的好处就体现出来了可以非常自然地写[ \mathbf{F}_{Marangoni} abla_s \sigma(\phi,T) ]它仍然是基于φ和温度场T的域内表达式不用找界面边界去施加。3. 上手实操T型通道中的液滴生成仿真3.1 从建模到网格那些不显眼却致命的细节案例我选T型微通道液滴生成这个是微流控里最经典的液滴制备结构适合演示水平集的全流程。几何尺寸上主流道宽度100微米分散相入口通道宽度50微米两通道垂直相交。把几何单位设置为微米但计算时COMSOL会自动转换我习惯在全局定义里把特征长度设为微米量级输入物理量时仍用SI制。这样做的好处是不容易在粘度、表面张力这些量上犯量级错误。几何建模很简单画一个长方体其实是2D的矩形域或者2D轴对称模型做主流道再画一个侧臂用“并集”把它们合并再加一个圆角过渡避免尖角处压力奇异。3D模型当然更贴近真实但2D模型在早期调参阶段可以节省大量时间和内存。等你把水平集参数调稳了再升维到3D也不迟。我第一次就直接上3D结果一天大半时间都耗在网格和求解器报错上回头来还是2D起步更理智。网格划分是最能体现经验的地方。两相流仿真不是把全局网格统一下加密就行——界面激战区域需要细网格主流道远离界面的区域完全可以用粗网格。具体操作是在界面可能经过的区域也就是T型交汇口往后的-一段距离内添加一个“网格细化域”最大单元尺寸设为2微米到5微米其余区域设在20微米左右。此时还要确认一点别让最大单元尺寸超过epsilon的2倍。我在第三版模型里把细化区画长了结果网格自由度数暴涨求解时间直接翻了好几倍效果却没有更好。界面区域的网格加密原则是“恰好覆盖界面扫过通道的区域”做液滴生成时就是T型口下游一段主通道做气泡上升时就是整个液滴轨迹路径。3.2 边界条件、接触角与初始液滴的布置边界条件设置看似常规但两相流里容易出问题。主流道入口给连续相速度比如0.1 m/s分散相入口给另一速度比如0.02 m/s这样两种流的流量比大概51适合形成滴状流。出口设压力为0壁面默认无滑移。接触角设置很容易被忽略。在“层流两相流水平集”接口下壁面边界上可以选择“润湿壁”特征输入接触角。我这里要提醒一句如果壁面接触角设了但网格在壁面附近太粗这个接触角实际能实现的角度会和设定值相差很大看起来像是没设一样。要想让接触角效果出来壁面附近的网格尺寸仍要减小最好在1到2个微米以内。初始液滴怎么布置水平集接口允许你给初始φ场赋值。我会用如下表达式生成一个初始圆形液滴圆心在T型口附近[ \phi_0 0.5 - 0.5 \tanh\left( \frac{\sqrt{(x-x_0)^2 (y-y_0)^2} - r_0}{\epsilon \cdot 2} \right) ]即把半径r0内的区域设为0相比如分散相外部设为1相连续相。这个表达式的意义在于φ在界面附近光滑过渡而不是简单地从0跳到1。直接赋0或1会造成数值震荡尤其在表面张力计算时会出现尖峰力。我早期用阶梯函数初始化时第一个时间步就报了“找不到一致的初始值”换成tanh平滑表达式后问题立刻消失。3.3 瞬态求解器设置与后处理提取液滴尺寸两相流是瞬态问题求解器设置直接影响稳定性和计算时间。我用的是“瞬态”研究时间步进方式选“BDF”并使用自由步长。BDF向后差分公式对刚性问题稳定但阶数太高会出现数值扩散所以一般选BDF2即可。初始步长建议给一个很小的值比如1e-7秒让前几步把接触角和表面张力导致的瞬态尖峰消化掉再进入自适应阶段。CFL条件是这类仿真里最重要的一条经验。CFL数定义为[ CFL \frac{u \Delta t}{h} ]要让水平集计算稳定CFL通常控制在0.5以内。如果特征速度0.1 m/s、网格尺寸2微米那么Δt至少要小于1e-5秒量级。COMSOL自适应时间步的长处是它会自己调整但前提是你给初始步长和最大步长的约束方向是对的。允许的最大步长如果设得过大稳定性就失控设得太小又白白浪费时间。后处理阶段最实用的一个操作是生成液滴直径随时间的变化曲线。方法是在结果里用“表面最大值”或“体积计算”表达式统计液滴相的体积。定义一个体积计算表达式[ V_d \int_\Omega \phi , d\Omega ]并除以一个液滴的参考体积就能看到液滴从初始时刻到脱离、再到完全生成的全过程。我常把这个表达式导出成数据表和实验照片对比验证模型准确性。COMSOL里还可以用“时间-频率FFT”分析液滴生成频率这对做高通量液滴生成设计特别有用。4. 实操中躲不开的5个坑与排查思路4.1 界面糊成一片扩散了仿真正跑着界面从一条清晰的线变成了一条渐变带这通常是epsilon取太大或者gamma太小导致的。很多教程把epsilon直接设为最大网格尺寸的一半但若算法的数值扩散本身比较强这个值还不够小。我的排查步骤是先把gamma提高到特征速度的1到2倍看界面是否重新锐化如果还糊再逐步减小epsilon同时按比例加密界面网格。要记住epsilon减小以后网格尺寸也要跟着减小保持网格尺寸≤epsilon/2。如果只缩小epsilon不加密网格曲率计算在粗糙网格上会产生强烈振荡。4.2 液滴体积越算越小水平集方法不带严格守恒修正长时间仿真液滴体积缓慢缩水是很常见的。这个问题没有一劳永逸的办法但有几个工程对策。第一检查界面区域是否至少有5层网格网格越粗体积损失越严重。第二控制CFL数时间步长过大时界面输运误差会迅速累积。第三确认初始φ场确实满足重初始化条件而不是生硬跳跃。第四在COMSOL的求解器中开启守恒形式的输运方程水平集方程能以守恒形式求解。我试过在相同网格下开启这个选项后液滴体积误差从8%降到3%以内。另外强烈建议在后处理里实时跟踪总液滴体积一旦发现体积误差超过5%就应该停下调参数而不是等仿真跑完才发现数据没法用。4.3 表面张力把界面冲出了锯齿表面张力项是水平集仿真最主要的不稳定来源。当曲率计算被网格噪声污染时界面会出现明显的锯齿状伪影甚至“振荡炸开”。这时先检查epsilon与网格的比例是否合理如果网格尺寸比epsilon大曲率中来自φ二阶导数的信号会被噪声淹没。其次检查表面张力系数是否被你设得过大导致局部体积力超过粘性扩散所能平衡的范围。第三考虑给表面张力项加一点人工平滑或者用更光滑的Heaviside函数。COMSOL里界面力可以通过“弱贡献”方式加入这时可以显式控制积分位置遇到锯齿问题可以试着把弱贡献改为体积力加平滑。4.4 接触角设了没反应这个问题我被问过很多次。接触角设置以后仿真结果里液滴接触线的表观接触角与设定值不一致原因大概率在壁面网格尺度和接触线附近的曲率计算。接触线附近流体速度变化剧烈压力奇异必须用足够细的网格捕获粘性应力。经验值是接触线附近网格尺寸至少比界面epsilon小2倍以上。另外接触角模型分静态与动态两种在高速液滴生成场景里动态接触角与静态接触角差异巨大。如果实验里接触角表现出明显速度依赖可以参考COMSOL的“动态接触角”选项配合一些经验关系式而不要只填一个固定角度。4.5 非线性求解器不收敛一启动就报错多半是初始条件不平滑或者界面初始化不合理。如果跑了几步以后发散最常见的原因是时间步太大。处理办法是把最大时间步长调小并打开“非线性求解器”的阻尼参数让压力速度耦合不那么激进。如果压力基方程带表面张力出现收敛困难可以在预处理阶段关闭表面张力项先跑一个纯对流-扩散的水平集问题把流场和φ场基本稳定后再开启表面张力。这个方法看着绕路实际能省下很多排查时间。另外网格质量检查一定要做单元质量低于0.3的区域尤其是T型口和壁面接触点附近往往是发散的根源。最后再分享两个技巧做了这么多水平集仿真之后我个人最大的体会是不要把水平集接口当作傻瓜式黑箱。ε和γ这两个参数虽然看起来只是一对数字但它们背后就是数值方法的稳定性与精度之间的杠杆。每改动一个物理参数比如流速、表面张力、通道宽度都要回头审视ε和γ是否需要同步调整。另一个很实用的习惯是——在正式跑一组参数扫描之前先用粗网格、大时间步跑一遍简单的“预测试”只看液滴能不能在预期时间尺度内生成。预测试通过之后再逐步细网格、严格CFL约束这样的流程比你直接上高精度配置反复崩溃不断试错要快得多。水平集方法在COMSOL里的表达能力很强液滴生成、气泡上升、微射流破裂这些问题都能得到不错的定量结果关键是你愿意花时间把底层参数吃透而不是遇到问题就去动网格。