简介Fluent二维波浪水槽的UDF造波是流体仿真中常见且易踩坑的需求这套材料围绕“1005_udf波”主题面向海洋工程、船舶设计及海岸防护方向的CFD学习者适合已掌握Fluent基础操作、希望深入UDF二次开发的读者。压缩包共3个文件包含C语言UDF源码、msh网格文件和TXT说明文档源码对应波浪生成算法网格用于划分计算域并分配物理属性说明文字系统梳理UDF的编译、链接与调用逻辑整体仅1.4MB轻量精悍便于对照学习。已有787人学习具备一定参考热度。内容覆盖线性/非线性波浪理论在UDF中的表达、Fluent源项与边界条件设置、求解器时间步长与迭代配置以及波高、速度场、压力分布等后处理思路通过实例代码与网格文件可直接作为实际造波模拟项目的改造起点也便于后续修改波浪参数、边界条件做二次开发。1. 用 UDF 造波等于给 Fluent 入口边界装了一个物理上正确的速度解析器抛一个反直觉的结论Fluent 里做波浪模拟最难的往往不是你写不出那个波浪公式而是波浪进了计算域之后怎么不让它“出去又弹回来”。很多刚接触 UDF 造波的人第一天就能把线性波的速度剖面写进DEFINE_PROFILE但跑到第 20 秒就发现出口边界把波反射回来整个自由液面变成一团乱麻。所以数值波浪水槽里真正决定成败的是入口怎么给出波、出口怎么吃掉波、中间怎么让 VOF 界面不发散这三件事。本文只讨论一条主线用 UDF 在 Fluent 里实现波浪模拟覆盖“操作”这个词背后大多数人真正需要的做法——速度入口造波附带推板式动边界造波作为替代方案。你会看到线性波和 Stokes 波在入口边界上的代码长什么样、阻尼消波区的 UDF 怎么写、libudf平台不匹配的报错怎么处理以及为什么网格分辨率、时间步和初始化的选择能毁掉一个完全正确的波。读者定位为已经在用 Fluent、知道边界条件面板长什么样、但第一次碰“用 UDF 造波”这条路的工程师。新手照着做能跑通老手可以跳过基础代码直接看第 4 章和第 6 章的坑。2. 造波 UDF 的前提选对波浪理论再谈速度表达式2.1 波浪理论选型的三个判据不管 UDF 怎么写入口边界给的速度场必须来自某个波浪理论。物理上没错数值上才有意义。工程上最常用的三个选择是线性波Airy 波、Stokes 五阶波、椭圆余弦波。绝大多数 Fluent 造波场景落在前两个。判断用哪个理论三个无量纲数就够了相对水深h/Lh/L 1/20是浅水1/20 ~ 1/2是中等水深 1/2是深水波陡H/L小于 0.02 时线性波误差尚可接受超过 0.05 建议直接上 Stokes 二阶以上相对波高H/h浅水中超过 0.4 就要小心非线性已不是微扰。有人会在深水条件h/L 0.5下依然用线性波这是安全区但在近岸防波堤、滩涂波浪这种相对水深明显小于 0.1 的场景线性波的底部水质点速度趋势完全不对这时候 UDF 写得再漂亮波高沿传播方向的变化也是失真的。判定量线性波AiryStokes 五阶椭圆余弦波Cnoidal适用相对水深h/L 0.1较安全h/L 0.1h/L 0.1适用波陡H/L 0.02H/L 0.1H/L 0.1UDF 表达式复杂度低闭式解中需要多个系数高需要雅可比椭圆函数Fluent 里的常见程度最常见教学和工程起步首选进阶海洋工程经常用极少见通常用 OpenFOAM我一般对刚上手的人只推荐线性波代码短、变量少、出问题好排查。Stokes 五阶波在 Fluent 里要处理 5 个待定系数其中涉及求解超越方程组放在 UDF 里不是不能写但调试成本翻倍。如果模型确实需要 Stokes 五阶建议先拿 MATLAB 把系数解出来再以静态数组方式嵌进 UDF而不是让 Fluent 每次迭代都解一遍。2.2 线性波速度场在入口边界上的数学形式线性波理论中取静水面为z 0向上为正水深为h即水底在z -h波面方程为η(x,t) A cos(kx - ωt)水平与垂向速度分量为u(x,z,t) Aω · cosh(k(zh)) / sinh(kh) · cos(kx - ωt) v(x,z,t) Aω · sinh(k(zh)) / sinh(kh) · sin(kx - ωt)其中波数k由色散关系ω² gk·tanh(kh)决定ω 2π/T。注意入口边界不是一个点而是一条线二维或一个面三维。UDF 要做的就是对入口边界上的每一个面根据它自己的z坐标算出对应的(u, v)。这也解释了为什么F_CENTROID这一步不可省略边界面上不同位置的 z 值不同速度也不同。底部附近速度趋近于零水面附近速度最大这是线性波速度剖面的基本形态后处理时可以用这个特点判断 UDF 有没有写反。2.3 UDF 里求波数 k不要手算代常数k和ω之间不是线性关系tanh让色散方程没有简单闭式解。UDF 里最稳的做法是用牛顿迭代在每次计算开始时求一次而不是在代码里硬编码一个近似值。硬编码的问题是当你把波高 A 从 0.05 改到 0.15水深从 0.5 改到 0.3原来那个 k 值可能偏离真实值几个百分点波形相位就会逐渐漂移。迭代初值建议取深水近似k₀ ω²/g对大多数中等水深情况都能在两三次迭代内收敛。临界点是避免在迭代里出现除零。df g·tanh(kh) g·kh·(1 - tanh²(kh))在kh很小的时候趋近于g·kh·(1) g·kh 2gkh不会为 0所以牛顿迭代在这类问题上是安全的。3. 用 DEFINE_PROFILE 在 Fluent 入口实现 UDF 造波从代码到边界参数3.1 搭建数值波浪水槽的边界布局二维为例左侧为造波入口右侧为出流边界顶部为压力出口底部为壁面。入口边界高度要覆盖整个波浪可能起伏的范围不要只画到静水面。很多人一开始把入口画到z 0齐平波峰高于静水面后就从入口顶部溢出去了速度场被截断波高直接就错了。几何尺寸的经验值水槽长度取 5~8 个波长水深方向网格在自由液面附近加密到波高的 1/10~1/20。入口前 1~2 个波长的范围内不要放任何结构物让波先充分发展。VOF 模型设置上选择“Volume of Fluid”的显式格式勾选“Sharp Interface Modeling”下的“Sharp/Dispersed”中的 Sharp 类型。隐式 VOF 格式在瞬态造波中时间精度较差而且对库朗数敏感不建议用于波浪模拟。3.2 完整可编译的线性波速度入口 UDF下面是一个二维线性波速度入口的完整 UDF把波数求解、水平速度、垂向速度、水体积分数四个部分放在一个文件中。竖直方向坐标按x[1]处理如果你的模型竖直方向是Z轴把下面所有x[1]改成x[2]。#include udf.h #include math.h #define PI 3.141592653589793 #define G 9.81 #define H_DEPTH 0.5 /* 水深单位 m */ #define WAVE_A 0.05 /* 波幅 波高/2单位 m */ #define WAVE_T 1.5 /* 波周期单位 s */ #define X_INLET 0.0 /* 入口处的 x 坐标 */ /* 色散关系用牛顿迭代求波数 k */ static real wave_number(real omega) { real k omega * omega / G; int i; for (i 0; i 50; i) { real kh k * H_DEPTH; real f G * k * tanh(kh) - omega * omega; real df G * (tanh(kh) kh * (1.0 - tanh(kh) * tanh(kh))); k k - f / df; } return k; } /* 入口水平速度 u */ DEFINE_PROFILE(inlet_u, thread, index) { real x[ND_ND]; real z, u, phase; real omega 2.0 * PI / WAVE_T; real k wave_number(omega); face_t f; begin_f_loop(f, thread) { F_CENTROID(x, f, thread); z x[1]; phase k * X_INLET - omega * CURRENT_TIME; u WAVE_A * omega * cosh(k * (z H_DEPTH)) / sinh(k * H_DEPTH) * cos(phase); F_PROFILE(f, thread, index) u; } end_f_loop(f, thread) } /* 入口垂向速度 v */ DEFINE_PROFILE(inlet_v, thread, index) { real x[ND_ND]; real z, v, phase; real omega 2.0 * PI / WAVE_T; real k wave_number(omega); face_t f; begin_f_loop(f, thread) { F_CENTROID(x, f, thread); z x[1]; phase k * X_INLET - omega * CURRENT_TIME; v WAVE_A * omega * sinh(k * (z H_DEPTH)) / sinh(k * H_DEPTH) * sin(phase); F_PROFILE(f, thread, index) v; } end_f_loop(f, thread) } /* 入口水体积分数z 低于波面给 1高于波面给 0 */ DEFINE_PROFILE(inlet_vof, thread, index) { real x[ND_ND]; real z, eta; real omega 2.0 * PI / WAVE_T; real k wave_number(omega); face_t f; eta WAVE_A * cos(k * X_INLET - omega * CURRENT_TIME); begin_f_loop(f, thread) { F_CENTROID(x, f, thread); z x[1]; if (z eta) F_PROFILE(f, thread, index) 1.0; else F_PROFILE(f, thread, index) 0.0; } end_f_loop(f, thread) }3.3 UDF 的编译与边界参数设置对照三个DEFINE_PROFILE分别对应入口边界在 Fluent 里需要挂接的三个量X Velocity、Y Velocity、Water Volume Fraction。操作路径User-Defined → Functions → Compiled点击Add选择源文件后按Build成功后再点Load。运行目录和 UDF 源文件路径都不能出现中文和空格这是老生常谈但高频踩坑的点。边界条件面板的对应关系如下Fluent 边界面板位置挂接的 UDF 名称说明Velocity Specification Method → Componentsinlet_u入口水平速度分量Velocity Specification Method → Componentsinlet_v入口垂向速度分量Multiphase → Volume Fractioninlet_vof水相体积分数必须与速度同步给注意当你手动指定了速度分量后Fluent 不再做法向/切向分解UDF 给出的u、v就是绝对坐标系下的速度。若入口边界是倾斜的斜坡堤前造波可能会遇到需要在 UDF 里做旋转坐标变换不要直接沿用这里的cos(phase)表达式。3.4 初始化混合初始化更适合入口边界设置完成后初始化方式直接决定了第一个波是否成立。混合初始化Hybrid Initialization会尝试预测一个流场但问题在于它会把 VOF 界面拍平成一个大的过渡带——如果计算域的水面不在你预设的z0位置初始化的气液界面会挤出莫名的速度脉冲。标准初始化的操作路径是先初始化全场Standard Initialization 填好水的速度、压力再用Adapt → Region框选出静水面以下的区域把Water Volume Fractionpatch 为 1。这个环节省不掉因为压力场和密度场不匹配会在前几步产生虚假的重力波叠加到你造的波上波高监测曲线前 1~2 秒就不能用。时间步方面线性波造波的初步估算Δt 0.2·Δx_min / (A·ω·coth(kh))分母是水面附近最大水平速度的近似。自由面附近的网格尺寸按波高的 1/15 取一般能得到可接受的波形。4. UDF 造波的另一半阻尼消波区代码与反射控制4.1 出口反射为什么是 UDF 造波的头号杀手速度入口把波送进计算域波传到出口边界时如果遇到固定壁面或简单的压力出口水体的垂向运动被约束动能就转化为势能反弹回来——这就是反射波。反射波会与入射波叠加形成驻波形态的伪波浪场波高监测数据直接报废。压力出口并非不能吸波但你几乎不可能调出一个在宽频范围内都有效的透射边界阻尼消波区是行业内的标准做法。常见的替代方案是 Fluent 自带的多孔介质区域或数值海滩Numerical Beach但前者要额外定义粘性阻力系数后者相关的 UDF 写法并不直观。自己写动量源项阻尼是最可控的做法。4.2 海绵层阻尼源项 UDF对垂向速度施加衰减阻尼消波区的原理很简单在出口前一段区域内对动量方程加一个与速度成正比的负源项让波浪在到达出口前被逐步耗散。阻尼系数从消波区起点到出口逐渐增大避免物理不连续产生的新反射。这里给出对垂向速度分量生效的 UDF 版本更适合抑制重力波#include udf.h #define DAMP_START 24.0 /* 消波区起点 x 坐标单位 m */ #define DAMP_END 30.0 /* 消波区终点水槽出口单位 m */ #define DAMP_COEF 12.0 /* 阻尼强度系数需按波周期调试 */ DEFINE_SOURCE(damping_v_source, c, t, dS, eqn) { real x[ND_ND]; real damp 0.0; real C_V_value; C_CENTROID(x, c, t); if (x[0] DAMP_START) { real xi (x[0] - DAMP_START) / (DAMP_END - DAMP_START); damp DAMP_COEF * xi * xi; /* 二次曲线起点处光滑过渡 */ } C_V_value C_V(c, t); dS[eqn] -C_R(c, t) * damp; return -C_R(c, t) * damp * C_V_value; }这个源项要挂到Cell Zone Conditions → 水区域 → Source Terms下选择Y Momentum二维或Z Momentum三维取决于轴向。DAMP_COEF的量纲是 1/s物理含义是速度衰减的速率。取值太小波穿过去还在出口弹取值太大会在消波区起点形成新的反射。DAMP_COEF 的调试顺序先定为4π/T观察波高衰减再逐步加密。4.3 消波区长度和阻尼系数的联动经验消波区长度取 1~2 倍波长阻尼系数从4~6/T下限试起。判定是否有效看出口前 0.5 波长处的波高实测值与理论值偏差是否在 5% 以内。第 4.2 节的二次权重是平滑的关键线性权重会在消波区起点处产生导数突变形成微弱但持续的二次波。反射率的粗略检测方法在水槽中段设一个波高探针取前 3 个周期的平均波高再用出口附近探针的值做对比。如果出口反射波传回中段后波高增加超过 8%先加消波区长度而不是盲目加DAMP_COEF。Fluent 后处理里看出入口的流量正负也有参考意义正常情况下入口Mass Flow Rate为负流入计算域如果看到周期性的正负交替说明有反射波到达入口区域边界条件已经受到污染。5. 推板造波用 DEFINE_CG_MOTION 驱动边界运动5.1 推板造波的运动学基础速度入口造波适合规则波的研究但在做浮体与波浪相互作用时有时候入口处的均匀速度场与物理造波机并不等价。推板造波Piston Wave Maker用一块竖直平板在底部做水平简谐运动通过板面推动水体产生波浪更贴近物理波浪水槽的实验数据。推板冲程 S 与目标波高 H 的关系用 Biesel 公式估算H / S 4 · sinh²(kh) / (sinh(2kh) 2kh)知道目标波高和周期先用色散关系求出 k再用公式反解冲程 S。注意这是在推板运动频率固定、波高较小的前提下成立的一阶近似波高大时非线性会导致实际波高低于公式估算值需要做两次试算修正。5.2 DEFINE_CG_MOTION 驱动刚性运动推板在 Fluent 里建模为一个薄矩形壁面与水体接触的一侧设为壁面边界整个推板区域设置为刚体运动由 UDF 控制#include udf.h #include math.h #define PI 3.141592653589793 #define STROKE 0.12 /* 推板冲程单位 m来自 Biesel 公式估算 */ #define WAVE_FREQ 0.667 /* 圆频率 ω 2π/T 对应的频率系数单位 Hz */ DEFINE_CG_MOTION(piston_motion, dt, vel, omega, time, dtime) { real w 2.0 * PI * WAVE_FREQ; real amplitude STROKE / 2.0; /* 推板沿 X 方向做正弦运动 */ vel[0] amplitude * w * cos(w * time); vel[1] 0.0; vel[2] 0.0; omega[0] 0.0; omega[1] 0.0; omega[2] 0.0; }vel是刚体中心的线速度omega是角速度。这里推板只有水平平动所以只设置了vel[0]。5.3 入口造波与推板造波的使用边界对比维度速度入口造波推板DEFINE_CG_MOTION造波自由液面捕捉由 VOF 计算入口只给速度剖面由 VOF 计算动边界推动水体非线性波支持需要换用高阶波理论UDF 复杂度大推板运动本身会产生非线性波动网格开销无需要 Dynamic Mesh计算量明显增加反射控制仍需出口消波推板处也会产生二次反射适合场景规则波、教学验证、参数扫描浮体耦合实验对标、非规则波做浮体运动响应时推板造波更接近实验但动网格与浮体六自由度运动耦合后网格质量退化速度快需要设置合理的网格重构阈值。一个折中方案是先用入口造波完成参数敏感性分析最后用推板造波做一两组验证性计算避免全程承担动网格的计算成本。6. UDF 造波结果验证探针、FFT 与三个高频坑6.1 用探针和 FFT 验证波高、波周期造波计算跑了 10 个周期不能只看云图上“有波的样子”。规范的验证流程是在水槽中段距离入口 2~3 个波长处布置监测点监测静水面位置的水体积分数变化时间间隔设为T/100。后处理里把该点的Volume Fraction曲线转成波面高程η (C_VOF - 0.5)·Δz_local再对时间序列做 FFT主频应该精确落在你设置的1/T上谐波能量占比不超过主峰的 10%。波高验证不能只看一个周期取第 3 到第 8 个周期的平均波峰谷差略去前两个周期波形尚未稳定和最后两个周期可能受消波区反射影响。如果 FFT 主频正确但波高偏小 5% 以上优先检查入口体积分数 UDF 的波面公式是否与速度场公式一致两者相位不同步是常见原因。6.2 三个一碰就翻车的 UDF 细节第一个坑是libudf库平台报错。错误信息形如error: the udf library you are trying to load (libudf) is not compiled for p...实际含义是当前 Fluent 进程与现存的 libudf 库针对不同处理器架构编译。解决办法不要直接Load外部拷来的 .so 或 .dll在 Fluent 内用Compiled UDFs重新Build。如果 Build 失败检查环境变量AWP_ROOT和VS路径是否指向当前版本的 Visual Studio 编译器Fluent 2020 之后的版本对编译器的版本匹配要求更严格。第二个坑是 VOF 初始化把界面抹平。混合初始化对气液两相流场的预测会把界面拉伸成几个网格宽的过渡带前几个周期的波高偏小几乎一定与此有关。标准做法是Standard Initialization后手动 patch 水的区域。第三个坑是自由度方向写错。二维模型竖直方向是 yx[1]没问题但从外部导入的模型如果一开始就在其他平面内建模可能在 XY 平面内水平方向和竖直方向搞混速度入口的 u、v 全给颠倒。排查方法很朴素把inlet_u临时改为常数 0.1inlet_v改为 0跑 10 步看入口附近是水平流动还是垂直流动。6.3 网格分辨率与时间步的匹配技巧救活一个“看起来完全正确但就是跑不动”的造波案例最后一步几乎都是统一网格与时间步的尺度。自由液面处相邻两层网格的高度比别超过 1.2网格最大尺寸与最小尺寸比控制在 5 以内避免数值波在粗细网格交界处发生伪反射。时间步的 CFL 数按 0.3 起步每个波周期至少 500 步这样等值线云图和时间序列曲线都会平滑很多。把这些排查顺序固定成流程下次造波项目就能直接照着走一遍。本文还有配套的精品资源点击获取 Fluent UDF造波实战：从速度入口到阻尼消波
简介Fluent二维波浪水槽的UDF造波是流体仿真中常见且易踩坑的需求这套材料围绕“1005_udf波”主题面向海洋工程、船舶设计及海岸防护方向的CFD学习者适合已掌握Fluent基础操作、希望深入UDF二次开发的读者。压缩包共3个文件包含C语言UDF源码、msh网格文件和TXT说明文档源码对应波浪生成算法网格用于划分计算域并分配物理属性说明文字系统梳理UDF的编译、链接与调用逻辑整体仅1.4MB轻量精悍便于对照学习。已有787人学习具备一定参考热度。内容覆盖线性/非线性波浪理论在UDF中的表达、Fluent源项与边界条件设置、求解器时间步长与迭代配置以及波高、速度场、压力分布等后处理思路通过实例代码与网格文件可直接作为实际造波模拟项目的改造起点也便于后续修改波浪参数、边界条件做二次开发。1. 用 UDF 造波等于给 Fluent 入口边界装了一个物理上正确的速度解析器抛一个反直觉的结论Fluent 里做波浪模拟最难的往往不是你写不出那个波浪公式而是波浪进了计算域之后怎么不让它“出去又弹回来”。很多刚接触 UDF 造波的人第一天就能把线性波的速度剖面写进DEFINE_PROFILE但跑到第 20 秒就发现出口边界把波反射回来整个自由液面变成一团乱麻。所以数值波浪水槽里真正决定成败的是入口怎么给出波、出口怎么吃掉波、中间怎么让 VOF 界面不发散这三件事。本文只讨论一条主线用 UDF 在 Fluent 里实现波浪模拟覆盖“操作”这个词背后大多数人真正需要的做法——速度入口造波附带推板式动边界造波作为替代方案。你会看到线性波和 Stokes 波在入口边界上的代码长什么样、阻尼消波区的 UDF 怎么写、libudf平台不匹配的报错怎么处理以及为什么网格分辨率、时间步和初始化的选择能毁掉一个完全正确的波。读者定位为已经在用 Fluent、知道边界条件面板长什么样、但第一次碰“用 UDF 造波”这条路的工程师。新手照着做能跑通老手可以跳过基础代码直接看第 4 章和第 6 章的坑。2. 造波 UDF 的前提选对波浪理论再谈速度表达式2.1 波浪理论选型的三个判据不管 UDF 怎么写入口边界给的速度场必须来自某个波浪理论。物理上没错数值上才有意义。工程上最常用的三个选择是线性波Airy 波、Stokes 五阶波、椭圆余弦波。绝大多数 Fluent 造波场景落在前两个。判断用哪个理论三个无量纲数就够了相对水深h/Lh/L 1/20是浅水1/20 ~ 1/2是中等水深 1/2是深水波陡H/L小于 0.02 时线性波误差尚可接受超过 0.05 建议直接上 Stokes 二阶以上相对波高H/h浅水中超过 0.4 就要小心非线性已不是微扰。有人会在深水条件h/L 0.5下依然用线性波这是安全区但在近岸防波堤、滩涂波浪这种相对水深明显小于 0.1 的场景线性波的底部水质点速度趋势完全不对这时候 UDF 写得再漂亮波高沿传播方向的变化也是失真的。判定量线性波AiryStokes 五阶椭圆余弦波Cnoidal适用相对水深h/L 0.1较安全h/L 0.1h/L 0.1适用波陡H/L 0.02H/L 0.1H/L 0.1UDF 表达式复杂度低闭式解中需要多个系数高需要雅可比椭圆函数Fluent 里的常见程度最常见教学和工程起步首选进阶海洋工程经常用极少见通常用 OpenFOAM我一般对刚上手的人只推荐线性波代码短、变量少、出问题好排查。Stokes 五阶波在 Fluent 里要处理 5 个待定系数其中涉及求解超越方程组放在 UDF 里不是不能写但调试成本翻倍。如果模型确实需要 Stokes 五阶建议先拿 MATLAB 把系数解出来再以静态数组方式嵌进 UDF而不是让 Fluent 每次迭代都解一遍。2.2 线性波速度场在入口边界上的数学形式线性波理论中取静水面为z 0向上为正水深为h即水底在z -h波面方程为η(x,t) A cos(kx - ωt)水平与垂向速度分量为u(x,z,t) Aω · cosh(k(zh)) / sinh(kh) · cos(kx - ωt) v(x,z,t) Aω · sinh(k(zh)) / sinh(kh) · sin(kx - ωt)其中波数k由色散关系ω² gk·tanh(kh)决定ω 2π/T。注意入口边界不是一个点而是一条线二维或一个面三维。UDF 要做的就是对入口边界上的每一个面根据它自己的z坐标算出对应的(u, v)。这也解释了为什么F_CENTROID这一步不可省略边界面上不同位置的 z 值不同速度也不同。底部附近速度趋近于零水面附近速度最大这是线性波速度剖面的基本形态后处理时可以用这个特点判断 UDF 有没有写反。2.3 UDF 里求波数 k不要手算代常数k和ω之间不是线性关系tanh让色散方程没有简单闭式解。UDF 里最稳的做法是用牛顿迭代在每次计算开始时求一次而不是在代码里硬编码一个近似值。硬编码的问题是当你把波高 A 从 0.05 改到 0.15水深从 0.5 改到 0.3原来那个 k 值可能偏离真实值几个百分点波形相位就会逐渐漂移。迭代初值建议取深水近似k₀ ω²/g对大多数中等水深情况都能在两三次迭代内收敛。临界点是避免在迭代里出现除零。df g·tanh(kh) g·kh·(1 - tanh²(kh))在kh很小的时候趋近于g·kh·(1) g·kh 2gkh不会为 0所以牛顿迭代在这类问题上是安全的。3. 用 DEFINE_PROFILE 在 Fluent 入口实现 UDF 造波从代码到边界参数3.1 搭建数值波浪水槽的边界布局二维为例左侧为造波入口右侧为出流边界顶部为压力出口底部为壁面。入口边界高度要覆盖整个波浪可能起伏的范围不要只画到静水面。很多人一开始把入口画到z 0齐平波峰高于静水面后就从入口顶部溢出去了速度场被截断波高直接就错了。几何尺寸的经验值水槽长度取 5~8 个波长水深方向网格在自由液面附近加密到波高的 1/10~1/20。入口前 1~2 个波长的范围内不要放任何结构物让波先充分发展。VOF 模型设置上选择“Volume of Fluid”的显式格式勾选“Sharp Interface Modeling”下的“Sharp/Dispersed”中的 Sharp 类型。隐式 VOF 格式在瞬态造波中时间精度较差而且对库朗数敏感不建议用于波浪模拟。3.2 完整可编译的线性波速度入口 UDF下面是一个二维线性波速度入口的完整 UDF把波数求解、水平速度、垂向速度、水体积分数四个部分放在一个文件中。竖直方向坐标按x[1]处理如果你的模型竖直方向是Z轴把下面所有x[1]改成x[2]。#include udf.h #include math.h #define PI 3.141592653589793 #define G 9.81 #define H_DEPTH 0.5 /* 水深单位 m */ #define WAVE_A 0.05 /* 波幅 波高/2单位 m */ #define WAVE_T 1.5 /* 波周期单位 s */ #define X_INLET 0.0 /* 入口处的 x 坐标 */ /* 色散关系用牛顿迭代求波数 k */ static real wave_number(real omega) { real k omega * omega / G; int i; for (i 0; i 50; i) { real kh k * H_DEPTH; real f G * k * tanh(kh) - omega * omega; real df G * (tanh(kh) kh * (1.0 - tanh(kh) * tanh(kh))); k k - f / df; } return k; } /* 入口水平速度 u */ DEFINE_PROFILE(inlet_u, thread, index) { real x[ND_ND]; real z, u, phase; real omega 2.0 * PI / WAVE_T; real k wave_number(omega); face_t f; begin_f_loop(f, thread) { F_CENTROID(x, f, thread); z x[1]; phase k * X_INLET - omega * CURRENT_TIME; u WAVE_A * omega * cosh(k * (z H_DEPTH)) / sinh(k * H_DEPTH) * cos(phase); F_PROFILE(f, thread, index) u; } end_f_loop(f, thread) } /* 入口垂向速度 v */ DEFINE_PROFILE(inlet_v, thread, index) { real x[ND_ND]; real z, v, phase; real omega 2.0 * PI / WAVE_T; real k wave_number(omega); face_t f; begin_f_loop(f, thread) { F_CENTROID(x, f, thread); z x[1]; phase k * X_INLET - omega * CURRENT_TIME; v WAVE_A * omega * sinh(k * (z H_DEPTH)) / sinh(k * H_DEPTH) * sin(phase); F_PROFILE(f, thread, index) v; } end_f_loop(f, thread) } /* 入口水体积分数z 低于波面给 1高于波面给 0 */ DEFINE_PROFILE(inlet_vof, thread, index) { real x[ND_ND]; real z, eta; real omega 2.0 * PI / WAVE_T; real k wave_number(omega); face_t f; eta WAVE_A * cos(k * X_INLET - omega * CURRENT_TIME); begin_f_loop(f, thread) { F_CENTROID(x, f, thread); z x[1]; if (z eta) F_PROFILE(f, thread, index) 1.0; else F_PROFILE(f, thread, index) 0.0; } end_f_loop(f, thread) }3.3 UDF 的编译与边界参数设置对照三个DEFINE_PROFILE分别对应入口边界在 Fluent 里需要挂接的三个量X Velocity、Y Velocity、Water Volume Fraction。操作路径User-Defined → Functions → Compiled点击Add选择源文件后按Build成功后再点Load。运行目录和 UDF 源文件路径都不能出现中文和空格这是老生常谈但高频踩坑的点。边界条件面板的对应关系如下Fluent 边界面板位置挂接的 UDF 名称说明Velocity Specification Method → Componentsinlet_u入口水平速度分量Velocity Specification Method → Componentsinlet_v入口垂向速度分量Multiphase → Volume Fractioninlet_vof水相体积分数必须与速度同步给注意当你手动指定了速度分量后Fluent 不再做法向/切向分解UDF 给出的u、v就是绝对坐标系下的速度。若入口边界是倾斜的斜坡堤前造波可能会遇到需要在 UDF 里做旋转坐标变换不要直接沿用这里的cos(phase)表达式。3.4 初始化混合初始化更适合入口边界设置完成后初始化方式直接决定了第一个波是否成立。混合初始化Hybrid Initialization会尝试预测一个流场但问题在于它会把 VOF 界面拍平成一个大的过渡带——如果计算域的水面不在你预设的z0位置初始化的气液界面会挤出莫名的速度脉冲。标准初始化的操作路径是先初始化全场Standard Initialization 填好水的速度、压力再用Adapt → Region框选出静水面以下的区域把Water Volume Fractionpatch 为 1。这个环节省不掉因为压力场和密度场不匹配会在前几步产生虚假的重力波叠加到你造的波上波高监测曲线前 1~2 秒就不能用。时间步方面线性波造波的初步估算Δt 0.2·Δx_min / (A·ω·coth(kh))分母是水面附近最大水平速度的近似。自由面附近的网格尺寸按波高的 1/15 取一般能得到可接受的波形。4. UDF 造波的另一半阻尼消波区代码与反射控制4.1 出口反射为什么是 UDF 造波的头号杀手速度入口把波送进计算域波传到出口边界时如果遇到固定壁面或简单的压力出口水体的垂向运动被约束动能就转化为势能反弹回来——这就是反射波。反射波会与入射波叠加形成驻波形态的伪波浪场波高监测数据直接报废。压力出口并非不能吸波但你几乎不可能调出一个在宽频范围内都有效的透射边界阻尼消波区是行业内的标准做法。常见的替代方案是 Fluent 自带的多孔介质区域或数值海滩Numerical Beach但前者要额外定义粘性阻力系数后者相关的 UDF 写法并不直观。自己写动量源项阻尼是最可控的做法。4.2 海绵层阻尼源项 UDF对垂向速度施加衰减阻尼消波区的原理很简单在出口前一段区域内对动量方程加一个与速度成正比的负源项让波浪在到达出口前被逐步耗散。阻尼系数从消波区起点到出口逐渐增大避免物理不连续产生的新反射。这里给出对垂向速度分量生效的 UDF 版本更适合抑制重力波#include udf.h #define DAMP_START 24.0 /* 消波区起点 x 坐标单位 m */ #define DAMP_END 30.0 /* 消波区终点水槽出口单位 m */ #define DAMP_COEF 12.0 /* 阻尼强度系数需按波周期调试 */ DEFINE_SOURCE(damping_v_source, c, t, dS, eqn) { real x[ND_ND]; real damp 0.0; real C_V_value; C_CENTROID(x, c, t); if (x[0] DAMP_START) { real xi (x[0] - DAMP_START) / (DAMP_END - DAMP_START); damp DAMP_COEF * xi * xi; /* 二次曲线起点处光滑过渡 */ } C_V_value C_V(c, t); dS[eqn] -C_R(c, t) * damp; return -C_R(c, t) * damp * C_V_value; }这个源项要挂到Cell Zone Conditions → 水区域 → Source Terms下选择Y Momentum二维或Z Momentum三维取决于轴向。DAMP_COEF的量纲是 1/s物理含义是速度衰减的速率。取值太小波穿过去还在出口弹取值太大会在消波区起点形成新的反射。DAMP_COEF 的调试顺序先定为4π/T观察波高衰减再逐步加密。4.3 消波区长度和阻尼系数的联动经验消波区长度取 1~2 倍波长阻尼系数从4~6/T下限试起。判定是否有效看出口前 0.5 波长处的波高实测值与理论值偏差是否在 5% 以内。第 4.2 节的二次权重是平滑的关键线性权重会在消波区起点处产生导数突变形成微弱但持续的二次波。反射率的粗略检测方法在水槽中段设一个波高探针取前 3 个周期的平均波高再用出口附近探针的值做对比。如果出口反射波传回中段后波高增加超过 8%先加消波区长度而不是盲目加DAMP_COEF。Fluent 后处理里看出入口的流量正负也有参考意义正常情况下入口Mass Flow Rate为负流入计算域如果看到周期性的正负交替说明有反射波到达入口区域边界条件已经受到污染。5. 推板造波用 DEFINE_CG_MOTION 驱动边界运动5.1 推板造波的运动学基础速度入口造波适合规则波的研究但在做浮体与波浪相互作用时有时候入口处的均匀速度场与物理造波机并不等价。推板造波Piston Wave Maker用一块竖直平板在底部做水平简谐运动通过板面推动水体产生波浪更贴近物理波浪水槽的实验数据。推板冲程 S 与目标波高 H 的关系用 Biesel 公式估算H / S 4 · sinh²(kh) / (sinh(2kh) 2kh)知道目标波高和周期先用色散关系求出 k再用公式反解冲程 S。注意这是在推板运动频率固定、波高较小的前提下成立的一阶近似波高大时非线性会导致实际波高低于公式估算值需要做两次试算修正。5.2 DEFINE_CG_MOTION 驱动刚性运动推板在 Fluent 里建模为一个薄矩形壁面与水体接触的一侧设为壁面边界整个推板区域设置为刚体运动由 UDF 控制#include udf.h #include math.h #define PI 3.141592653589793 #define STROKE 0.12 /* 推板冲程单位 m来自 Biesel 公式估算 */ #define WAVE_FREQ 0.667 /* 圆频率 ω 2π/T 对应的频率系数单位 Hz */ DEFINE_CG_MOTION(piston_motion, dt, vel, omega, time, dtime) { real w 2.0 * PI * WAVE_FREQ; real amplitude STROKE / 2.0; /* 推板沿 X 方向做正弦运动 */ vel[0] amplitude * w * cos(w * time); vel[1] 0.0; vel[2] 0.0; omega[0] 0.0; omega[1] 0.0; omega[2] 0.0; }vel是刚体中心的线速度omega是角速度。这里推板只有水平平动所以只设置了vel[0]。5.3 入口造波与推板造波的使用边界对比维度速度入口造波推板DEFINE_CG_MOTION造波自由液面捕捉由 VOF 计算入口只给速度剖面由 VOF 计算动边界推动水体非线性波支持需要换用高阶波理论UDF 复杂度大推板运动本身会产生非线性波动网格开销无需要 Dynamic Mesh计算量明显增加反射控制仍需出口消波推板处也会产生二次反射适合场景规则波、教学验证、参数扫描浮体耦合实验对标、非规则波做浮体运动响应时推板造波更接近实验但动网格与浮体六自由度运动耦合后网格质量退化速度快需要设置合理的网格重构阈值。一个折中方案是先用入口造波完成参数敏感性分析最后用推板造波做一两组验证性计算避免全程承担动网格的计算成本。6. UDF 造波结果验证探针、FFT 与三个高频坑6.1 用探针和 FFT 验证波高、波周期造波计算跑了 10 个周期不能只看云图上“有波的样子”。规范的验证流程是在水槽中段距离入口 2~3 个波长处布置监测点监测静水面位置的水体积分数变化时间间隔设为T/100。后处理里把该点的Volume Fraction曲线转成波面高程η (C_VOF - 0.5)·Δz_local再对时间序列做 FFT主频应该精确落在你设置的1/T上谐波能量占比不超过主峰的 10%。波高验证不能只看一个周期取第 3 到第 8 个周期的平均波峰谷差略去前两个周期波形尚未稳定和最后两个周期可能受消波区反射影响。如果 FFT 主频正确但波高偏小 5% 以上优先检查入口体积分数 UDF 的波面公式是否与速度场公式一致两者相位不同步是常见原因。6.2 三个一碰就翻车的 UDF 细节第一个坑是libudf库平台报错。错误信息形如error: the udf library you are trying to load (libudf) is not compiled for p...实际含义是当前 Fluent 进程与现存的 libudf 库针对不同处理器架构编译。解决办法不要直接Load外部拷来的 .so 或 .dll在 Fluent 内用Compiled UDFs重新Build。如果 Build 失败检查环境变量AWP_ROOT和VS路径是否指向当前版本的 Visual Studio 编译器Fluent 2020 之后的版本对编译器的版本匹配要求更严格。第二个坑是 VOF 初始化把界面抹平。混合初始化对气液两相流场的预测会把界面拉伸成几个网格宽的过渡带前几个周期的波高偏小几乎一定与此有关。标准做法是Standard Initialization后手动 patch 水的区域。第三个坑是自由度方向写错。二维模型竖直方向是 yx[1]没问题但从外部导入的模型如果一开始就在其他平面内建模可能在 XY 平面内水平方向和竖直方向搞混速度入口的 u、v 全给颠倒。排查方法很朴素把inlet_u临时改为常数 0.1inlet_v改为 0跑 10 步看入口附近是水平流动还是垂直流动。6.3 网格分辨率与时间步的匹配技巧救活一个“看起来完全正确但就是跑不动”的造波案例最后一步几乎都是统一网格与时间步的尺度。自由液面处相邻两层网格的高度比别超过 1.2网格最大尺寸与最小尺寸比控制在 5 以内避免数值波在粗细网格交界处发生伪反射。时间步的 CFL 数按 0.3 起步每个波周期至少 500 步这样等值线云图和时间序列曲线都会平滑很多。把这些排查顺序固定成流程下次造波项目就能直接照着走一遍。本文还有配套的精品资源点击获取