简介：本资源面向流体动力学仿真初学者与海洋工程领域CFD实践者，聚焦二阶Stokes波理论在ANSYS Fluent平台的工程化实现，解决波浪边界条件难以精确建模的核心问题。压缩包共2个文件（136KB），含关键UDF源码文件（C语言编写）与配套二维案例设置文件（.cas），前者完整封装二阶Stokes波面位移、速度分量及非线性项的数学表达式，后者已预设网格拓扑、边界类型、求解器参数及初始条件，开箱即可加载运行。已有218人学习下载，适用于海岸防护结构受力分析、浮式平台波浪载荷预测、船舶耐波性初步评估等典型场景。用户可直接编译调用该UDF生成物理意义明确的周期性入射波边界，避免手动设置复杂时变函数，显著降低非线性波模拟门槛，并为后续耦合自由液面VOF模型或波浪-结构相互作用研究提供可靠基础。 做水动力仿真的朋友，大概率都被波浪边界条件折磨过。最近可算把 bolang.rar 这个项目里的二阶Stokes波浪模拟彻底吃透了，里面用的是 Fluent UDF 造波。网上很多教程只给一个正弦波速度入口，但那东西在 VOF 模型里根本稳不住，波面跑两步就散架。这篇我就从 UDF 公式出处讲起，把代码、Fluent 设置、调参和踩坑经验都写透，非常适合正在做无反射造波、船舶水动力评估，或者单纯想给波浪 UDF 做模板的同行。
1. 项目概述与核心需求
1.1 波浪模拟为什么需要UDF
Fluent 的内置边界条件里没有“造波”这个选项，入口的空气和水之间有一条随时间变化的相界面，你要把水面高度和速度同时按波浪理论给进去，普通的速度入口或压力入口根本做不到。很多人以为用瞬态仿真加上一个正弦速度信号就够了，但实际算出来波面很快就衰减。核心原因是入口的相体积分数没跟着波面更新，导致水相和空气相的质量通量在整个入口面上不匹配。波浪 UDF 的本质，就是在入口边界上同时控制两件事：速度分布（水平与垂直分量）和水的体积分数。VOF 模型负责捕捉界面的演化，UDF 负责把波浪解析解映射到边界上。
UDF（User-Defined Function）是 Fluent 提供的一套 C 语言扩展机制，允许用户动态设置边界条件、源项、物性参数等。在波浪模拟这个场景里，它相当于一个“会呼吸的入口”：每个时间步都根据波浪理论算出当前时刻入口面上每个网格点的速度分量和相分数，再交给 Fluent 去推进计算。如果你的目标只是看一个定常流动，根本不需要 UDF；但波浪是典型的非定常自由面流动，内置边界条件没有任何一项能覆盖这种动态过程。
1.2 二阶Stokes波浪理论基础
线性 Airy 波只保留一阶项，波面和速度都是严格正弦的，波峰和波谷完全对称。真实波浪在有限水深下会出现波峰变尖、波谷变平的现象，这时需要二阶 Stokes 理论。二阶 Stokes 波是在一阶解的基础上叠加了一个倍频项，所以波面公式和速度表达式都会多出二阶修正项。简单来说，一阶项给出主波频率的振荡，二阶项给出二倍频的振荡，这个二倍频分量正好让波面变得不对称。
设静水面为 z=0，水深为 h，海底在 z=-h，波浪沿 x 方向传播。圆频率为 ω=2π/T，波数为 k=2π/λ，色散关系为 ω²=gk·tanh(kh)。一阶波幅 A=H/2，其中 H 是波高。实际做 UDF 时，我们需要把二阶波面方程和速度方程一起写进去。经常有人只把一阶速度塞进去，结果波面形状不对、能量不平衡，所以二阶项必须加全。
这里有一套实用形式的二阶 Stokes 速度分量，取 z 轴向上，原点在静水面：
水平速度：
u(z,t) = A ω [cosh k(h+z) / sinh(kh)] cos(kx−ωt) + (3/4) A² ω k [cosh 2k(h+z) / sinh⁴(kh)] cos(2(kx−ωt))
垂直速度：
v(z,t) = A ω [sinh k(h+z) / sinh(kh)] sin(kx−ωt) + (3/4) A² ω k [sinh 2k(h+z) / sinh⁴(kh)] sin(2(kx−ωt))
二阶波面方程：
η(x,t) = A cos(kx−ωt) + (A² k / 4) [ (2 + cosh(2kh)) cosh(kh) / sinh³(kh) ] cos(2(kx−ωt))
注意，这些公式在不同文献里系数可能有差异，取决于速度势展开时的截断方式。写进 UDF 的最重要原则是：速度表达式和波面方程必须来自同一套推导，否则入口的速度场和自由面位置不闭合，算出来波面一定发散。
1.3 目标场景与软件适配
这类方案在 Fluent 2020R1、2021R1 上都实测可跑，二维和三维都能用。典型应用场景包括：波浪对固定结构的砰击载荷分析、浮式结构物响应、波浪能装置水动力评估、LNG 液舱晃荡问题等。适用前提是波浪参数不要太极端，如果波陡度 H/λ 大于 1/15，二阶 Stokes 的精度就开始下降，容易出现数值破碎。如果只是研究远场波浪传播，建议在计算域出口加一段消波区，否则出口反射波会把你的波浪场搅乱。
2. 波浪UDF的设计与编写
2.1 整体思路：从波面方程到边界条件
边界造波法是当前最常用的造波方式。在入口边界上，每个时间步每个网格节点，我们需要算出三个量：水平速度分量、垂直速度分量、水的体积分数。入口如果是垂直平面，面上的网格点就靠全局坐标 y 区分高度，所以判断当前点是在波面以下还是以上，只需要拿该点的垂直坐标和当前时刻的波面高度 η 做比较。
我把 UDF 拆成两层：第一层是辅助函数，比如波面计算、速度计算；第二层是 DEFINE_PROFILE 主函数，负责遍历入口面所有 face，把算好的值赋给 F_PROFILE。这样代码结构清晰，换波浪理论时只需要改辅助函数，不必动主循环。如果你要模拟的是不规则波，甚至可以在此基础上增加一个随机相位叠加模块。
2.2 关键公式的UDF实现
写 UDF 之前，先定义好波浪参数和全局常量。因为 Fluent 的 UDF 每次调用都在同一进程内，静态变量可以用来缓存波数和圆频率，避免每个 face 都重复计算。
#include "udf.h" #define PI 3.141592653589793 #define GRAVITY 9.81 /* 波浪参数，按实际工况修改 */ #define WAVE_HEIGHT 0.2 /* H，波高，单位m */ #define WATER_DEPTH 2.0 /* h，水深，单位m */ #define WAVE_PERIOD 3.0 /* T，周期，单位s */ #define WAVE_LENGTH 5.0 /* lambda，波长，单位m */ #define Z_SURF 0.0 /* 静水面的全局Y坐标 */ static real omega = 0.0; static real k = 0.0;接着写一个初始化函数，在第一次调用时计算 ω 和 k：
void wave_init(void) { if (omega == 0.0) { omega = 2.0 * PI / WAVE_PERIOD; k = 2.0 * PI / WAVE_LENGTH; } }这个写法可以避免在循环里频繁调用三角函数，一定程度上减少计算量。虽然对于小规模网格差别不大，但我习惯这么做，后面换到三维波浪时能明显感觉到收益。
然后写波面辅助函数。注意这里 z 是全局坐标，我要在函数内部转换成相对静水面的垂向坐标：
real wave_eta(real x, real t) { real A = WAVE_HEIGHT / 2.0; real kh = k * WATER_DEPTH; real S = sinh(kh); real C = cosh(kh); real phase = k * x - omega * t; real eta2 = A * A * k / 4.0 * ((2.0 + cosh(2.0 * kh)) * C / (S * S * S)) * cos(2.0 * phase); return A * cos(phase) + eta2; }2.3 完整UDF代码逐段解析
下面是一个可直接编译的 UDF 片段，包含水平速度、垂直速度和水相 VOF 三部分。核心思路是：先用 F_CENTROID 拿到 face 中心的全局坐标，再用坐标做判断，最后把结果赋给 F_PROFILE。
DEFINE_PROFILE(wave_x_vel, thread, index) { face_t f; real x[ND_ND]; real t = CURRENT_TIME; real z, u, phase; real A = WAVE_HEIGHT / 2.0; real kh = k * WATER_DEPTH; real sh = sinh(kh); real ch = cosh(kh); wave_init(); begin_f_loop(f, thread) { F_CENTROID(x, f, thread); z = x[1] - Z_SURF; phase = k * x[0] - omega * t; /* 一阶水平速度 */ u = A * omega / sh * cosh(k * (z + WATER_DEPTH)) * cos(phase); /* 二阶水平速度 */ u += 0.75 * A * A * omega * k / (sh * sh * sh * sh) * cosh(2.0 * k * (z + WATER_DEPTH)) * cos(2.0 * phase); F_PROFILE(f, thread, index) = u; } end_f_loop(f, thread) }垂直速度的写法类似，只是把 cosh 换成 sinh，cos 换成 sin：
DEFINE_PROFILE(wave_y_vel, thread, index) { face_t f; real x[ND_ND]; real t = CURRENT_TIME; real z, v, phase; real A = WAVE_HEIGHT / 2.0; real kh = k * WATER_DEPTH; real sh = sinh(kh); wave_init(); begin_f_loop(f, thread) { F_CENTROID(x, f, thread); z = x[1] - Z_SURF; phase = k * x[0] - omega * t; /* 一阶垂直速度 */ v = A * omega / sh * sinh(k * (z + WATER_DEPTH)) * sin(phase); /* 二阶垂直速度 */ v += 0.75 * A * A * omega * k / (sh * sh * sh * sh) * sinh(2.0 * k * (z + WATER_DEPTH)) * sin(2.0 * phase); F_PROFILE(f, thread, index) = v; } end_f_loop(f, thread) }水相 VOF 的 profile 相对简单，直接比较 face 中心的垂向位置和波面高度：
DEFINE_PROFILE(wave_vof, thread, index) { face_t f; real x[ND_ND]; real t = CURRENT_TIME; real eta; wave_init(); begin_f_loop(f, thread) { F_CENTROID(x, f, thread); eta = wave_eta(x[0], t); if (x[1] - Z_SURF <= eta) F_PROFILE(f, thread, index) = 1.0; /* 水面以下为水 */ else F_PROFILE(f, thread, index) = 0.0; /* 水面以上为空气 */ } end_f_loop(f, thread) }有人会担心，入口处的 VOF 突然从 0 跳到 1，会不会带来压力振荡？在采用 Geo-Reconstruct 的 VOF 模型里，入口边界的 alpha 值默认就是这个阶梯分布，实际运行下来问题不大。如果你实在不放心，可以在波面附近的过渡带做一个线性插值，但那样会人为引入一个模糊界面，反而影响波面精度。我更倾向于直接 0/1 跳变，然后把入口附近网格加密。
2.4 编译与加载的注意事项
UDF 的加载分为解释型和编译型两种。建议优先使用编译型（Compiled），因为代码中用到了 sinh/cosh 等数学库函数，解释型虽然也能跑，但效率偏低。在 Fluent 中，路径是 User Defined -> Functions -> Compiled，添加源文件后点击 Build，编译成功后再点击 Load。
需要特别留意几个点：
- 工程路径不能有中文，目录最好全英文，否则编译器会报找不到头文件。
- 如果改了 C 代码，必须重新 Build，不能只点 Load 旧库文件。
- 边界条件对话框里，选择对应方向的 velocity profile 时，要手动下拉选到 wave_x_vel、wave_y_vel；相分数 profile 要选 wave_vof。
- 可以在代码里加
Message("t=%f, u=%f\n", t, u);这类打印，观察控制台输出，但大规模网格下别每步都打印，否则 I/O 会很慢。
我这套代码在 Fluent 2021R1 + VS2019 环境下编译正常。如果你的 Fluent 版本较老，比如 17.0 之前，可能需要在编译前设置path环境变量，让编译器能找到 cl.exe。
3. 在Fluent中的仿真设置与实操流程
3.1 网格与求解器准备
网格质量直接决定波浪是否衰减。我的经验是，二维模型入口方向网格尺寸控制在波长的 1/50 到 1/100，垂向至少 20 个节点。以波长 5m 为例，水平方向网格取 0.1m，垂向取 0.05m 左右，入口边界附近因为速度梯度大，最好再加密一层。
以下是推荐的基础设置表，可以直接抄：
| 设置项 | 推荐值 |
|---|---|
| 求解器类型 | Pressure-Based，Transient |
| 多相流模型 | VOF，显式 Geo-Reconstruct |
| 湍流模型 | 层流 或 SST k-omega |
| 压力-速度耦合 | PISO |
| 动量离散格式 | Second Order Upwind |
| 压力离散格式 | Body Force Weighted |
| 体积分数离散 | Geo-Reconstruct |
很多波浪模拟在低雷诺数下用层流就够了。如果关注波浪破碎或者结构物绕流，再用 SST k-omega 也不迟。但要注意，湍流模型会把波面的能量耗散掉一部分，波高会略微衰减，所以验证时要留出误差余量。
3.2 边界条件与初始场设置
这个造波方案不涉及动网格，入口就是一个固定的速度入口。边界条件设置如下：
- 入口（inlet）：velocity-inlet，速度用 wave_x_vel 和 wave_y_vel，相分数用 wave_vof。
- 出口（outlet）：pressure-outlet，最好在后端加一段阻尼消波区。
- 底部：wall，无滑移。
- 顶部：如果空气区足够高，用 pressure-outlet，相对压力设为 0。
初始场是关键。很多人一开始把整个域 patch 成水，结果入口刚开始发射波浪，流域内部压力场和入口速度场不匹配，产生虚假振荡。正确做法是初始化后，用Adapt -> Region选择静水面以下的区域，patch 水相体积分数为 1；静水面以上保持空气。如果静水面以上的空气层太薄，波浪运动过程中可能撞到顶部边界，所以空气区高度建议至少 0.5 倍波长。
3.3 时间步长与松弛因子调优
时间步长的选择通常用库朗数约束：CFL = U_max * dt / dx。对 VOF 显式格式，CFL 建议小于 2。入口处最大质点速度大约为 Aω，波速 c = λ/T。举个例子，波高 0.2m，波长 5m，周期 3s，Aω ≈ 0.21 m/s，c ≈ 1.67 m/s，dx 取 0.1m，则dt < 0.1 / (1.67 + 0.21) ≈ 0.053s。实际我取 0.01s，每个周期 300 步，模拟 10 个周期就是 3000 步，计算量完全可接受。
松弛因子参考：压力设为 0.3，动量 0.7，湍流 0.5。压力速度耦合选 PISO，因为瞬态自由面问题用 PISO 比 SIMPLE 更能保证稳定性。如果计算过程中出现发散，优先把压力松弛因子继续下调到 0.2，同时把时间步长降低一个数量级，先让初始波传播出去，再逐步提高步长。
3.4 结果后处理与验证方法
算完之后，不要直接拿云图看个颜色就觉得完事了。要验证波面是否准确，我通常会在入口附近、流域中部、靠近出口处各放一个监测点，记录水体积分数随时间的变化。在 CFD-Post 或 Tecplot 里提取Water.VOF = 0.5的等值面，就是自由面。然后用这个自由面高度和时间的关系去和理论二阶 Stokes 波面公式做比较。
还可以对波高时间序列做 FFT，检查频谱中是否包含二倍频分量。二阶 Stokes 波的核心特征就是这个二倍频分量，如果 FFT 结果里只有基频，说明你的 UDF 可能只写入了一阶项，二阶项没生效。
我实测的案例是：H=0.2m，h=2m，T=3s，λ=5m，模拟 20 个周期后，入口下游第 4 个波长处波高误差在 5% 以内。如果你算出来误差超过 10%，优先检查网格尺寸和出口消波区是否足够。
4. 常见问题与排查技巧实录
4.1 UDF编译报错排查
UDF 编译报错是最常见的第一道坎。我遇到过这些典型报错：
undeclared identifier：变量未定义或大小写写错，比如CURRENT_TIME写成了current_time。- 找不到
udf.h：工程路径有中文，或者 Fluent 找不到 Visual Studio 编译器环境。 F_PROFILE未定义：可能是因为忘记包含udf.h，或者主函数里漏写了thread参数。
排查技巧很简单：把代码中分支函数逐一注释，定位挂掉的那一行；加上Message打印中间值，看是不是出现了 NaN。另外，Fluent 在 Windows 下的 UDF 编译依赖 Visual Studio 版本，比如 2021R1 通常需要 VS2019，版本不匹配会报各种莫名其妙的链接错误。
4.2 波浪衰减与数值耗散
模拟中波高沿传播方向逐渐变小，这是数值耗散的直接体现。主要原因有三个：网格太粗、时间步长太大、动量方程离散格式一阶精度。解决方法也很直接：
- 加密网格，尤其自由面附近网格；垂向至少保证波高方向有 10 个网格。
- 时间步长降到
dt < 0.1 * λ / (c + U_max)，不要为了进度硬撑。 - 动量方程用 Second Order Upwind，VOF 保持 Geo-Reconstruct。
如果波衰减还是明显，就得在出口前加阻尼消波区。最朴素的做法是把出口区域网格拉长，在 UDF 里用 DEFINE_SOURCE 给动量方程添加一个阻尼项，让波浪在进入消波区的过程中能量逐渐被吸收。这个源项的形态通常是S = -ρ * f_damp * (u - u_background)，阻尼系数从 0 逐渐增大到某个值，保证过渡平滑。
4.3 波面发散或压力不稳定的解决
波面发散经常发生在开头几个时间步。现象是入口刚给速度，流域内压力场没跟上，产生高压脉冲，波面瞬间碎掉。我的解决方法是加一个“软启动”函数：在入口速度前面乘一个平滑因子ramp = min(1.0, t / (2.0 * T))，也就是前两个周期让波幅从 0 逐渐增长到目标值。这样初始压力冲击被抹平，等波浪场建立起来后再满幅运行。
另一个常见问题是入口的 VOF 相分数和速度 UDF 不协调。如果速度沿 y 方向突变，入口界面处会产生非物理对流。实际处理时，速度 UDF 可以直接作用在整个入口面，包括空气区；但如果你在空气区给了很大的波浪速度，空气也会被带动，形成虚假剪切层。我的做法是：速度 UDF 只在水中按公式赋速度，在空气区把速度设为一个很小的值（例如 0.01 m/s），这样空气基本静止，避免干扰水面。
4.4 参数敏感性速查表
| 参数 | 影响 | 经验取值 |
|---|---|---|
| 水平网格尺寸 dx | 波高衰减随 dx 增大而上升 | λ/50 ~ λ/100 |
| 垂向网格尺寸 dy | dy 过大导致波面阶梯状变形 | ≤ H/10 |
| 时间步长 dt | dt 过大导致界面发散 | 满足 CFL ≤ 2 |
| 波陡度 H/λ | 二阶Stokes适用上限约 1/15 | 超过则换高阶或数值造波 |
| 水深比 h/λ | 有限水深范围适用二阶Stokes | 0.1 < h/λ < 0.5 |
| 消波区长度 | 太短时反射波干扰入口 | ≥ 1 个波长 |
这张表也是我后来做波浪类项目时的默认检查项。每次仿真跑飞，先对照表查一遍，多半能定位问题。
5. 一些经验与扩展建议
5.1 造波方法的选择与对比
边界速度造波适合中低波陡、计算域较短的情况，但缺点是无法处理结构物强反射。如果波浪遇到大尺度结构后反射回入口，反射波会和造波边界相互作用，导致波场失真。这种情况下，真正的主动吸收式造波需要在入口速度里叠加一个反向传播速度，计算方法基于入口监测到的实际波面，用反馈控制器实时调节。这比普通 UDF 复杂得多，但效果也更好。
另一种思路是改成源项造波，在流域内部一定高度的带状区域添加质量源项，通过源项强度控制波浪参数。源项造波的最大优势是边界上不设置入口，反射波会穿过源区，不会被边界弹回去。就是需要在 UDF 里把源项形状搞得精细一些，比如用高斯分布包络。
5.2 UDF在其他场景的复用
这套 UDF 的可复用性很强。把二阶 Stokes 公式替换成五阶 Stokes、孤立波、或者 JONSWAP 不规则波的速度谱，就能模拟不同海况。换汤不换药，核心还是“把当前时刻的相位、坐标带入波谱公式，得到边界处的速度和相分数”。
掌握这套逻辑后，你也会发现 UDF 在其他物理场景里的通用性。比如蒸发冷凝模型的 Lee 模型 UDF，本质上就是给质量传递方程加了一个源项；Presto 平台的自定义 UDF 也是类似思路，只不过调用的 API 不同。归根结底，UDF 能让用户把领域公式写进通用仿真平台，波浪模拟只是一个特别直观的例子。
最后分享一个小技巧：拿到 bolang.rar 这样的压缩包时，别急着看代码。先把二阶 Stokes 的理论公式在纸上推导一遍，重点确认坐标原点和正方向，然后对着 UDF 逐行核对。很多看似玄学的数值发散，最后都是因为某个公式的符号反了，或者二阶项系数少乘了一个 2。把这些基础工作做扎实，后面的仿真调试会顺利得多。
本文还有配套的精品资源，点击获取
