简介：本资源面向流体力学仿真工程师、海洋工程研究人员及CFD进阶学习者，聚焦二阶Stokes波理论在ANSYS Fluent平台上的工程化实现，解决波浪非线性运动建模与边界条件动态加载的核心难点。压缩包为136KB的RAR格式，共含2个关键文件：一个C语言编写的UDF源码文件（stokes-2.c），完整实现了二阶Stokes波面位移、速度分量及周期性边界更新逻辑；一个配套的二维波浪模拟Case文件（2Dbolang.cas），已预设网格拓扑、求解器参数、边界类型及UDF挂载路径，开箱即可运行验证。已有218人学习下载，资源虽小但结构完整——既提供可直接编译调用的波浪入口函数，又包含经配置验证的仿真案例框架，显著降低从理论公式到CFD实操的转化门槛，适用于海岸防护结构受力分析、浮式平台波激响应等典型应用场景。 “bolang.rar”这个压缩包名字，在流体仿真资源站上已经快成一种符号了：解压之后里面通常是几个UDF源文件加一个现成的case，能跑，但轮到自己改水深、波高、周期时，问题就全冒出来了——编译过不了、波浪跑几步就衰减没影、入口附近碎成一团。这篇文章就把Stokes二阶波浪UDF这件事从头到尾拆开讲：理论选型、参数怎么算、代码每一行在干嘛、Fluent里怎么配、出口消波怎么写、跑完之后怎么验证，一步不落。
1. 这个压缩包里的核心：Stokes二阶造波UDF到底在干什么
1.1 一阶到二阶：UDF里多出来的那一项才是关键
先说结论：所谓“Stokes二阶波浪UDF”，本质上是在Fluent的入口边界上用C语言复现一个解析速度场。这个速度场由一个一阶项和一个二阶项叠加而成。
一阶项就是大家熟悉的线性波解，也就是Airy波理论。对速度势求导后，水质点的水平速度和垂向速度分别是：
$u^{(1)} = a\omega \frac{\cosh k(y+d)}{\sinh kd} \cos(kx-\omega t)$
$v^{(1)} = a\omega \frac{\sinh k(y+d)}{\sinh kd} \sin(kx-\omega t)$
其中 $a$ 是波幅，$\omega = 2\pi/T$ 是圆频率，$k$ 是波数，$d$ 是静水深，$y$ 是垂向坐标（静水面为零，向上为正，水底 $y=-d$）。
二阶项呢？它叠加一个频率为 $2(kx-\omega t)$ 的速度分量：
$u^{(2)} = \frac{3}{4}a^2\omega k \frac{\cosh 2k(y+d)}{\sinh^4 kd} \cos 2(kx-\omega t)$
$v^{(2)} = \frac{3}{4}a^2\omega k \frac{\sinh 2k(y+d)}{\sinh^4 kd} \sin 2(kx-\omega t)$
体现在波面上就更直观。线性波的自由水面是一条纯余弦曲线，波峰和波谷对称。二阶Stokes波的波面是：
$\eta = a\cos\theta + \frac{a^2 k}{4}\frac{\cosh kd (2+\cosh 2kd)}{\sinh^3 kd}\cos 2\theta$
其中 $\theta = kx-\omega t$。加上第二项之后，波形不再是纯余弦——波峰变尖、波谷变平，这才是真实海洋中有限振幅波的样子。UDF里多出来的那部分二阶项，做的工作就是把这种不对称性“塞”进入口边界。
1.2 速度入口和动量源造波，选哪种更适合你
实现波浪UDF入口边界，主流有三种做法，这里先给个对比，方便你根据项目阶段选择：
| 方法 | 实现复杂度 | 优点 | 缺点 |
|---|---|---|---|
| 速度入口（边界给定速度） | 低 | 代码量小，调试直观，适合固定入口规则波 | 波在入口有反射风险，需要配消波区 |
| 动量源项造波（域内加源项） | 高 | 入射波和反射波可分离，适合不规则波 | 需要额外一块造波源区，参数多 |
| 动网格/重叠网格造波 | 很高 | 适合大幅晃动、结构耦合问题 | 网格开销大，稳定性需要经验积累 |
绝大多数“bolang.rar”用的都是第一种——速度入口法。原因很简单：它不改变网格拓扑，不需要额外造波区域，把入口面的速度和一个理论解绑定就结束。代价是入口边界对反射没有抵抗力，所以后面必须做阻尼消波区，这个我在第5章专门讲。
2. 开工前的参数账本：水深、周期、波高的联动关系
2.1 选型判据：Stokes二阶波不是哪里都能用
写UDF之前先确定一个问题：你想模拟的波浪，到底适不适合用Stokes二阶理论。这个判断做错了，后面所有参数都没意义。
Stokes波理论成立的前提是小振幅、中等水深。工程上常用几个判据：
- 相对水深 $kd$：水深与波长的比值。一般要求 $kd$ 大致在 $0.1\pi$ 到 $3$ 之间，即中间水深到深水范围。
- Ursell数：$U_r = \frac{H L^2}{d^3}$，它衡量非线性和浅水效应的相对强弱。当 $U_r > 30$ 左右时，波浪的非线性太强，Stokes展开收敛很差，需要换椭圆余弦波或孤立波理论。
- 波陡 $H/L$：一般要求 $H/L$ 不超过约 0.06，超过这个值波浪接近破碎，任何基于摄动展开的理论都失效。
举个实际例子：水深10 m，周期6 s，波高0.6 m。先用色散关系算出波长 $L$ 约55 m，那么 $kd = 2\pi d / L \approx 1.14$，波陡 $H/L \approx 0.011$，Ursell数 $U_r = \frac{0.6 \times 55^2}{10^3} \approx 0.18$。这组参数用Stokes二阶完全没有问题，二阶项对波面的修正量只有厘米级，波高再大一点才会更明显。
反过来，如果水深只有2 m、周期8 s、波高1 m，Ursell数会很容易过30，Stokes二阶算出来明显失真，这种场景就不得不考虑椭圆余弦波了。选型这一步不要省，很多算不下去的case根源都在这里。
2.2 波数k的求解：色散关系与牛顿迭代
波数 $k$ 是整份UDF里最底层的参数，它由线性色散关系决定：
$\omega^2 = gk\tanh(kd)$
这是一个超越方程，没法直接写出 $k$ 的显式表达式，通常用牛顿迭代求解。在UDF里可以单独写一个求解函数，在计算初始化时调用一次。
#include "udf.h" #define PI 3.141592653589793 #define G 9.81 #define DEPTH 10.0 /* 静水深 */ #define T_WAVE 6.0 /* 波浪周期 */ static real omega; static real k_wave; real compute_wavenumber(real depth, real T) { real w = 2.0 * PI / T; real k = 0.1; /* 迭代初值 */ real f, fp; int i; for (i = 0; i < 100; i++) { f = G * k * tanh(k * depth) - w * w; fp = G * (tanh(k * depth) + k * depth / (cosh(k * depth) * cosh(k * depth))); k -= f / fp; if (fabs(f) < 1e-8) break; } return k; }初值给0.1一般就够了。如果水深特别浅或者特别深，可以把迭代初值改成 $\sqrt{\omega^2/g}$ 或 $\omega/\sqrt{gd}$，收敛会更快。但实际测试里，从0.1起步，100次内都能稳定收敛到一个合理值。
注意：迭代算出来的 $k$ 必须在计算正式开始前赋值给全局变量，否则后续profile宏里读到的 $k_wave$ 是零，速度全部变成NaN，求解器一迭代就崩。初始化用DEFINE_INIT最省事：
DEFINE_INIT(init_wave, mixture_domain) { omega = 2.0 * PI / T_WAVE; k_wave = compute_wavenumber(DEPTH, T_WAVE); Message0("Stokes2: omega=%.4f k=%.4f L=%.4f\n", omega, k_wave, 2.0*PI/k_wave); }进入计算后，控制台会输出omega和k的数值，检查一眼心里就有底。
2.3 波面高度公式与入口网格尺寸的配合
入口边界不只要给速度，还要定义水位在哪里。波面高度用二阶Stokes表达式：
$\eta(x,t) = a\cos\theta + \frac{a^2 k}{4}\frac{\cosh kd (2+\cosh 2kd)}{\sinh^3 kd}\cos 2\theta$
入口附近的网格尺寸必须能分辨这个波面。经验上：
- 沿传播方向，每个波长至少80~100个网格。
- 自由液面附近，每个波高至少10~15层网格。
- 入口处第一排网格宽高比尽量接近1。
如果入口网格太粗，波面在边界上无法被VOF的几何重构正确捕捉，出来的波马上就会产生寄生波，入口附近的水体看起来像沸腾了一样。所以参数不只是算数，还要落到网格设计上。
3. UDF代码拆解：一个可直接抄作业的完整实现
3.1 完整入口速度UDF
把前面所有讨论汇成一个完整的UDF文件。二维数值水槽，坐标原点在静水面，$x$ 向右为波传播方向，$y$ 向上。入口在计算域左端。
#include "udf.h" #define PI 3.141592653589793 #define G 9.81 #define DEPTH 10.0 #define H_WAVE 0.6 #define T_WAVE 6.0 static real omega; static real k_wave; real compute_wavenumber(real depth, real T) { real w = 2.0 * PI / T; real k = 0.1; real f, fp; int i; for (i = 0; i < 100; i++) { f = G * k * tanh(k * depth) - w * w; fp = G * (tanh(k * depth) + k * depth / (cosh(k * depth) * cosh(k * depth))); k -= f / fp; if (fabs(f) < 1e-8) break; } return k; } real wave_elevation(real x, real t) { real a = H_WAVE / 2.0; real theta = k_wave * x - omega * t; real s = sinh(k_wave * DEPTH); real c = cosh(k_wave * DEPTH); real eta = a * cos(theta) + a * a * k_wave / 4.0 * (c * (2.0 + cosh(2.0 * k_wave * DEPTH)) / (s * s * s)) * cos(2.0 * theta); return eta; } DEFINE_INIT(init_wave, mixture_domain) { omega = 2.0 * PI / T_WAVE; k_wave = compute_wavenumber(DEPTH, T_WAVE); Message0("Stokes2: omega=%.4f k=%.4f L=%.4f\n", omega, k_wave, 2.0*PI/k_wave); } DEFINE_PROFILE(wave_vel_x, thread, position) { face_t f; real xc[ND_ND]; real t = CURRENT_TIME; real a = H_WAVE / 2.0; real y, theta, u; begin_f_loop(f, thread) { F_CENTROID(xc, f, thread); y = xc[1]; theta = k_wave * xc[0] - omega * t; u = a * omega * cosh(k_wave * (y + DEPTH)) / sinh(k_wave * DEPTH) * cos(theta) + 3.0/4.0 * a * a * omega * k_wave * cosh(2.0 * k_wave * (y + DEPTH)) / pow(sinh(k_wave * DEPTH), 4.0) * cos(2.0 * theta); F_PROFILE(f, thread, position) = u; } end_f_loop(f, thread) } DEFINE_PROFILE(wave_vel_y, thread, position) { face_t f; real xc[ND_ND]; real t = CURRENT_TIME; real a = H_WAVE / 2.0; real y, theta, v; begin_f_loop(f, thread) { F_CENTROID(xc, f, thread); y = xc[1]; theta = k_wave * xc[0] - omega * t; v = a * omega * sinh(k_wave * (y + DEPTH)) / sinh(k_wave * DEPTH) * sin(theta) + 3.0/4.0 * a * a * omega * k_wave * sinh(2.0 * k_wave * (y + DEPTH)) / pow(sinh(k_wave * DEPTH), 4.0) * sin(2.0 * theta); F_PROFILE(f, thread, position) = v; } end_f_loop(f, thread) }这段代码可以直接放进记事本存成stokes2.c，在Fluent的UDF编译对话框中导入编译。代码适用于二维情况，三维的话速度入口的展向不需要变化，直接同样挂载即可。
3.2 代码里容易被忽略的细节
初看这段代码可能觉得简单，但有几个细节直接影响成败。
第一，F_CENTROID(xc, f, thread)获取的是边界面的中心坐标，Fluent默认二维时xc[0]是x，xc[1]是y。如果你的几何模型初始是把水面放在z轴方向，那对应关系就要调整，不要机械照抄。提前确定坐标轴含义，能省一个晚上的调试时间。
第二，sinh(k_wave * DEPTH)在入水口底部是有限值。当相对水深很浅时，$kd$ 小，$\sinh(kd)$ 也小，二阶项系数会急剧增大，速度剖面会出现异常大的数值。这就是之前说的“浅水不适合Stokes二阶”的数值化表现。如果编译完一运行就发散或者速度值大得离谱，先回第2章查适用条件。
第三，两个profile宏都需要依赖于初始化后的k_wave和omega。如果你不用DEFINE_INIT，而是在算例运行后才手动初始化，那k_wave就会一直是0，速度函数里的正弦余弦参数全部失真。建议每次进入计算前，在Console里确认这行输出存在：
Stokes2: omega=1.0472 k=0.1160 L=54.18看到这种输出，才说明初始化真正生效了。
3.3 VOF液面入口：光有速度不够
速度只解决一半问题。入口边界上，还得分清水相和气相的位置。典型做法是在入口边界上给体积分数：当 $y < \eta(x=0,t)$ 时，水的体积分数为1；当 $y > \eta$ 时，水的体积分数为0。
这里有一个很常见的坑：在Fluent的VOF模型里，速度入口的相分数profile和速度profile挂载的线程不一定相同。你在设置边界条件时，需要分别对水相和空气相挂载各自的速度profile，同时给水相一个体积分数profile。体积分数profile可以直接用C语言判断：
DEFINE_PROFILE(water_vof, thread, position) { face_t f; real xc[ND_ND]; real t = CURRENT_TIME; real eta; begin_f_loop(f, thread) { F_CENTROID(xc, f, thread); eta = wave_elevation(xc[0], t); if (xc[1] < eta) F_PROFILE(f, thread, position) = 1.0; else F_PROFILE(f, <p> <a href="https://download.csdn.net/download/weixin_42656416/86591163" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>
