VOF法波浪模拟
理论：波浪特性与水面追踪
概述
我想用CFD计算打到海上结构物上的波浪，但造出的波走一段就变小了。
这是数值水槽最常见的烦恼。原因多半是水面附近网格相对波浪太粗，或时间步长太大。VOF法输运的是水的比例，网格粗时水面会被抹平，波浪能量随之损失。先正确求出目标波的波长，确认一个波长和波高内各有多少网格。波长由周期和水深决定，要用色散关系计算。
色散关系
$\omega = 2\pi/T$ 为角频率，$k$ 为波数，$h$ 为水深，$L$ 为波长，$c$ 为波速。深水（$kh$ 大）时 $L = gT^2/2\pi$；水变浅时，同一周期波长变短、波速变慢。周期越长的波越容易受水深影响。
VOF法
$\alpha$ 为网格内水的体积分数（1为水，0为空气）。按体积分数混合水和空气的密度，作为一个流场求解。水面位于 $\alpha = 0.5$ 处，并使用保持界面清晰的特殊对流格式（界面压缩等）。
用“体积分数”表示水面的思路
VOF法由美国洛斯阿拉莫斯实验室的赫特和尼科尔斯于1981年发表。此前的方法沿水面排列点来跟踪，或用函数表示表面形状，波浪破碎或水滴飞溅时就会失效。只给每个网格一个数——“装了多少水”——的思路，无论水面变形多复杂都能处理。此后又不断改进以保持界面清晰，如今从船舶阻力到溃坝、容器内液体晃荡都广泛使用。
算例
例1：不同周期和水深的波长
| 周期 | 深水波长 | 水深20 m波长 | 波速 | kh |
|---|---|---|---|---|
| 4 s | 25.0 m | 25.0 m | 6.2 m/s | 5.03 |
| 6 s | 56.2 m | 55.0 m | 9.2 m/s | 2.28 |
| 8 s | 99.9 m | 88.8 m | 11.1 m/s | 1.42 |
| 10 s | 156.1 m | 121.2 m | 12.1 m/s | 1.04 |
周期4秒的波在水深20 m时波长也与深水相同，但10秒的波比深水公式短两成以上。仍按深水公式搭建水槽，长周期波的网格参考和消波段长度都会出现偏差。
例2：数值水槽网格与计算量
水深20 m、波高2 m、长5个波长的二维水槽。水面±2 m带内加密，其余竖向网格粗4倍：
| 周期 | 每波长／每波高网格数 | 横向尺寸 | 竖向尺寸 | 二维网格数 | 每周期步数 |
|---|---|---|---|---|---|
| 6 s | 80 / 10 | 0.69 m | 0.20 m | 约2.1万 | 160 |
| 6 s | 100 / 20 | 0.55 m | 0.10 m | 约5.2万 | 200 |
| 10 s | 80 / 10 | 1.52 m | 0.20 m | 约2.1万 | 160 |
| 10 s | 100 / 20 | 1.21 m | 0.10 m | 约5.2万 | 200 |
二维只要几万个网格。三维会怎样？
横向也同样加密的话，网格数要再乘以横向网格数，很快就到几百万到几千万。所以三维中只在结构物附近加密，造波和消波区段放粗，或者远处用势流波浪模型求解、只在近处与CFD耦合。另外要注意10秒波时横竖网格尺寸比超过7倍。比值过大时水面计算容易不稳定，竖向加密时横向也要相应调整。
时间步长
库朗数控制在0.5左右时，每周期约100～200步。VOF法为保持界面清晰，常把库朗数取得更小（约0.25）。
数值水槽的搭建
步骤
- 由周期、波高、水深求波长，选择波浪理论（线性、斯托克斯、椭圆余弦）。
- 设置造波区段和抑制反射的消波区段（1～2个波长）。
- 先在无结构物状态下造波，确认结构物位置处的波高和周期符合目标。
- 改变网格和时间步长，确认波高衰减和结构受力不变。
“波从水槽尽头返回来”
计算打到结构物上的波浪力时，力随计算进行逐渐增大——原因是水槽出口侧消波不足，反射波返回与入射波叠加。消波区段短时，周期越长的波越容易反射。结构物反射的波也会传回造波侧，所以造波区段也需要设置吸收反射。基本的检查是在无结构物时确认水槽内各处波高不变（没有形成驻波）。
常见错误
错误与对策
| 错误 | 后果 | 对策 |
|---|---|---|
| 水面附近网格粗 | 波浪衰减 | 每波长80～100、每波高10～20 |
| 按深水公式定波长 | 长周期波长出错 | 求解色散关系 |
| 消波区段短 | 反射波叠加 | 保证1～2个波长 |
| 长宽比大的网格 | 水面不稳定 | 控制比值 |
| 一上来就带结构物计算 | 无法确认波浪质量 | 先用空水槽确认 |
我想自己试试波浪计算。