垂直高度函数方法在液滴模拟中的创新与应用
1. 垂直高度函数方法的核心原理与创新
在计算流体力学（CFD）领域，液滴与固体表面的相互作用模拟一直是极具挑战性的课题。传统水平高度函数（HHF）方法在处理与坐标轴平行的固体表面时表现良好，但当固体表面与坐标平面接近平行（如接近水平的xy平面）且接触角θ较小时，HHF方法会完全失效。这种情况对应于所谓的"几何退化"场景——在任何笛卡尔方向都无法重建有效界面。
1.1 水平高度函数方法的局限性
HHF方法的核心假设是界面重建必须与某个特定坐标方向对齐。这种假设在以下两种典型情况下会导致严重问题：
- 当接触角θ < 45°或θ > 135°时
- 当固体表面法向量与坐标轴接近平行时（如接近水平的表面）
在这些情况下，HHF方法会产生显著的界面重建误差，进而影响接触角计算的准确性。例如，在模拟水滴在超疏水表面（接触角>150°）的行为时，传统方法往往无法得到物理合理的结果。
1.2 垂直高度函数方法的创新设计
针对HHF的局限性，垂直高度函数（VHF）方法通过三个关键创新解决了这些问题：
PLIC界面表示法：采用分段线性界面构造（Piecewise Linear Interface Construction, PLIC）技术，不再要求界面与任何特定坐标方向对齐。这使得液/气界面Γₗ可以在任意方向与固体边界Γₛ相交，精确满足预设的接触角θ条件。
接触角依赖的法向量：在单元(i,j,k)中重建界面时，使用基于接触角θ的法向量nₗ,θ来确定界面方向（图10中的红色虚线面）。这个法向量同时考虑了表面张力和固液相互作用的影响。
双高度值定义：
- 流体区域高度值h_f由重建界面质心（图10绿色点）定义
- 接触角约束点由交线段MM'的中点（图10蓝色点）确定
关键提示：VHF方法保留了HHF方法中抛物面拟合的核心步骤（算法2的步骤3除外），包括基于流体区域HF点的预拟合（式38）、接触线相对位置的四个辅助HF点计算，以及结合接触角条件的非线性系统求解（式39）。这使得VHF能够充分利用HHF的成熟框架，同时扩展其适用范围。
2. VHF方法的实现细节与技术要点
2.1 界面重建与高度函数计算
VHF方法的界面重建过程可分为以下步骤：
接触线单元识别：在混合单元（同时包含流体和固体的单元）中，根据体积分数和相邻单元状态确定潜在的接触线位置。
法向量计算：根据预设接触角θ和固体表面法向量nₛ，计算液/气界面的法向量nₗ,θ：
# 伪代码：接触角依赖的法向量计算 def compute_interface_normal(n_s, contact_angle): # n_s: 固体表面单位法向量 # contact_angle: 预设接触角（弧度） n_l = rotate_vector(n_s, contact_angle) # 绕接触线旋转接触角度 return normalize(n_l)PLIC界面定位：使用计算得到的nₗ,θ确定界面位置，使得：
- 界面与单元边界的交线满足体积分数约束
- 界面与固体表面的交线满足接触角条件
高度函数值确定：
- h_f = 界面质心到参考平面的垂直距离
- h_c = 接触点到参考平面的垂直距离
2.2 抛物面拟合与曲率计算
为了准确计算界面曲率和表面张力，VHF方法采用二次抛物面拟合技术：
预拟合阶段：仅使用流体区域内的HF点拟合初始抛物面（式38）：
z = ax^2 + bxy + cy^2 + dx + ey + f辅助点计算：在接触线附近确定四个辅助HF点（图11b中的红点），这些点的位置由接触线法向量nₗ,θ在xy平面的投影nₗ,θ-xy指导。
最终拟合：求解包含接触角条件的非线性系统（式39），得到最终抛物面参数。这个系统通过牛顿迭代法求解，通常3-5次迭代即可收敛。
表1对比了HHF与VHF在抛物面拟合关键步骤中的差异：
| 步骤 | HHF方法 | VHF方法 |
|---|---|---|
| 界面方向确定 | 对齐坐标轴 | 由nₗ,θ决定，任意方向 |
| 辅助点位置确定 | 使用nₛ-yz投影 | 使用nₗ,θ-xy投影 |
| 高度函数扩展方向 | 从混合单元向相邻固体单元扩展 | 从接触线液体侧向气体侧扩展 |
2.3 接触角滞后效应的建模
在实际应用中，接触角往往不是单一值，而是在前进角θₐ和后退角θᵣ之间变化的。VHF方法通过以下策略模拟这种滞后效应：
滞后窗口定义：为每个接触线单元预设[θᵣ, θₐ]范围，反映表面粗糙度或化学异质性。
二分搜索算法：在体积分数更新后，通过迭代寻找使接触线位置与上一时间步相同的θ值：
- 固定接触线切线向量τₛ保持几何一致性
- 体积分数在迭代过程中保持不变
- 通常需要5-7次迭代达到收敛
界面重构：一旦确定θ值，应用标准抛物面拟合程序强制执行接触角条件。
图12展示了这种迭代过程的示意图，其中虚线表示上一时间步的界面Γₗⁿ⁻¹/²，实线表示更新后的界面Γₗⁿ⁺¹/²，两者在固体表面相交于同一接触线位置。
3. 数值验证与性能分析
3.1 体积守恒测试
通过同心球体模型（图13）验证VHF方法的体积守恒性能。测试配置包括：
- 静止固体球（半径rₛ=0.4m）
- 旋转液体壳层（外半径r_d=0.6m）
- 刚性旋转速度场（式43）
表2对比了传统几何VOF平流方案与修正方案在19.2 cpr分辨率下的表现：
| 指标 | 传统方案 | 修正方案 |
|---|---|---|
| L∞误差 | 6.1×10⁻² | 7.9×10⁻¹⁵ |
| L₂误差 | 1.4×10⁻² | 7.0×10⁻¹⁶ |
| 全局体积损失 | ~10⁻¹⁴ | ~10⁻¹⁴ |
| CPU时间 | 950s (CFL约束) | 1238s (+30%) |
虽然修正方案增加了30%的计算成本，但将局部体积误差降低到机器精度级别，这对于长期模拟的准确性至关重要。
3.2 表面张力驱动液滴铺展
3.2.1 平板上液滴铺展
验证不同接触角（15°-165°）和表面取向（nₛ=[0,0,1]和[1,1,1]）下的铺展行为。关键发现：
嵌入边界不敏感性：如图18所示，抛物面拟合方法在不同嵌入深度d（0.2-0.8）下保持一致的铺展半径，相对误差<1%，而线性拟合方法误差随d增大而显著增加。
接触线对称性：图19显示抛物面拟合方法在倾斜表面（nₛ=[1,1,1]）上仍能保持完美的圆形接触线，而线性拟合产生明显三角形畸变。
体积守恒：全局体积变化保持在3×10⁻⁸以内（图20b），与无嵌入参考案例相当。
3.2.2 球面上液滴铺展
测试更复杂的曲面几何（图21），比较3D VHF结果与2D轴对称参考解。表3总结了关键指标：
| 接触角θ | 半径相对误差 | 体积损失E_V |
|---|---|---|
| 30° | 1.0×10⁻³ | 9.4×10⁻⁷ |
| 60° | 2.6×10⁻³ | 1.0×10⁻⁶ |
| 90° | 3.5×10⁻³ | 1.4×10⁻⁷ |
| 120° | 5.4×10⁻³ | 5.8×10⁻⁸ |
| 150° | 5.5×10⁻³ | 4.0×10⁻⁷ |
图22展示了不同接触角下的接触线轮廓，3D结果与理论预测（式48）高度一致，验证了方法在曲面上的准确性。
3.3 剪切流驱动液滴运动
3.3.1 平板上滞后效应
模拟剪切流驱动下接触角滞后现象（图23），两个测试案例参数见表4。关键观察：
- 案例A（θᵣ=40°, θₐ=90°）：后接触线立即滑动，前接触线在θ达到90°前保持钉扎，导致液滴前部隆起（图24a）
- 案例B（θₐ=90°, θᵣ=50°）：前接触线立即移动，后接触线保持钉扎直到θ降至50°，形成流线型液滴（图24b）
表5对比了稳态液滴尺寸与文献结果，显示良好一致性：
| 参数 | 案例A (本文) | 案例A [35] | 案例B (本文) | 案例B [36] |
|---|---|---|---|---|
| L | 2.33 | 2.36 | 2.51 | 2.49 |
| W | 2.96 | 2.72 | 2.00 | 2.00 |
| H | 0.77 | 0.73 | 0.88 | 0.86 |
3.3.2 正弦波状基底上的液滴运动
更复杂的几何测试（图25），正弦表面由式50定义。图26展示了液滴在六个时间点的形态演变：
- t₀-t₁：毛细主导的对称铺展，接触角从初始~90°趋向设定值60°
- t₁后：剪切流效应显现，液滴开始变形和迁移
- t₂,t₄：液滴遇到波峰时速度降低，流体在波谷积聚形成局部隆起
- t₃,t₅：液滴越过波峰后尾部拉长变尖
这些结果与Zhu等[37]的实验观察定性一致，证实了方法处理复杂几何的能力。
4. 工程应用与实施建议
4.1 实际应用场景选择指南
VHF方法特别适用于以下场景：
- 极端接触角情况（θ < 45°或θ > 135°）
- 具有倾斜或复杂曲率的固体表面
- 需要精确模拟接触线动力学的应用（如液滴冲击、微流体器件）
- 涉及接触角滞后的多相流问题
对于简单几何和中等接触角（45°-135°），传统HHF方法可能仍具计算效率优势。
4.2 计算参数选择经验
基于大量测试，推荐以下实践参数：
网格分辨率：
- 一般应用：25-35 cpr（cells per droplet radius）
- 高精度需求：50+ cpr
- 注意：在接触线区域需要局部加密
时间步长：
# 伪代码：自适应时间步计算 def compute_time_step(CFL, dx, u_max, sigma, rho): dt_advection = CFL * dx / u_max dt_capillary = sqrt(rho * dx**3 / (2*pi*sigma)) return min(dt_advection, dt_capillary)- 典型CFL数：0.25-0.5
- 高表面张力情况需考虑毛细时间步限制
抛物面拟合参数：
- 牛顿迭代容差：1×10⁻⁶
- 最大迭代次数：10
- 正则化参数：1×10⁻⁸（防止矩阵奇异）
4.3 常见问题排查
表6总结了典型问题及解决方案：
| 问题现象 | 可能原因 | 解决方案 |
|---|---|---|
| 接触线位置振荡 | 时间步过大 | 减小CFL数至0.25以下 |
| 体积损失逐渐累积 | 平流方案不守恒 | 启用体积重分配方案 |
| 小接触角下界面破裂 | 曲率计算不准确 | 增加抛物面拟合迭代次数 |
| 滞后模拟中接触线不钉扎 | 二分搜索容差过大 | 将角度容差减小到0.1°以下 |
| 高剪切流下液滴畸变 | 表面张力系数过低 | 增加σ或减小剪切率 |
4.4 性能优化技巧
自适应网格细化：在接触线区域动态加密网格，其他区域保持粗网格，可节省30-50%计算资源。
混合并行化：
- 使用MPI进行节点间并行
- 在节点内使用OpenMP加速关键循环
- 对抛物面拟合等计算密集型部分考虑GPU加速
接触线检测优化：
# 伪代码：快速接触线检测 def detect_contact_cells(phi, solid_frac): contact_mask = (phi > 0) & (phi < 1) & (solid_frac > 0) # 仅处理接触线附近3×3×3邻域 return dilate(contact_mask, 3x3x3)内存访问优化：对高度函数数组使用结构体封装，提高缓存利用率，实测可加速15-20%。
5. 扩展应用与未来方向
5.1 工业应用案例
喷墨打印：VHF方法可精确模拟微小液滴（pL级）在多种基材上的沉积过程，包括：
- 多孔介质渗透
- 非均匀表面润湿
- 高精度图案形成
微流体器件设计：用于优化：
- 液滴生成芯片的几何参数
- 表面处理方案选择
- 流动控制策略评估
涂层工艺优化：模拟复杂表面上的液膜扩展，预测：
- 涂层均匀性
- 缺陷形成风险
- 干燥过程应力分布
5.2 方法扩展可能性
动态接触角模型：将当前静态接触角条件扩展为速度依赖的动态模型：
θ_d = f(θ_s, Ca, ...)其中Ca为毛细数，θ_s为静态接触角。
多物理场耦合：整合：
- 相变（蒸发/凝结）
- 溶解/沉淀过程
- 电润湿效应
机器学习加速：使用神经网络：
- 替代昂贵的抛物面拟合迭代
- 预测接触线动态行为
- 实时参数调优
5.3 开源实现建议
对于希望实现VHF方法的研究者，推荐以下开发路径：
基础框架选择：
- 已有VOF代码：集成VHF模块
- 从头开发：考虑Basilisk或OpenFOAM架构
关键数据结构：
struct VHFData { double h_fluid; // 流体区域高度值 double h_contact; // 接触点高度值 vec3d n_l; // 界面法向量 vec3d contact_line; // 接触线位置 double curvature; // 拟合曲率 };验证流程：
- 单元测试：验证各组件（法向量计算、抛物面拟合等）
- 基准测试：比较文献案例（如本节5.1-5.3）
- 网格收敛研究：确认理论精度阶数
在实际应用中，我们发现保持严格的代码模块化至关重要——将几何计算、界面重建、高度函数更新等分离为独立模块，大幅提高了代码可维护性和扩展性。对于大规模3D模拟，建议采用层次化网格数据结构来高效管理不同尺度的界面特征。