ANSYS FLUENT UDF实战：自定义波浪边界与衰减模型开发指南

1. 项目概述：当计算流体力学遇上自定义边界

在船舶与海洋工程、海岸防护以及近海结构物设计领域，波浪与结构的相互作用是核心的物理问题。我们通常使用ANSYS FLUENT这类计算流体动力学软件进行数值模拟，但软件自带的波浪模型，比如VOF（Volume of Fluid）方法结合波浪边界，有时难以满足一些特殊的、非标准的波浪工况需求。比如，你想模拟一个特定谱型的随机波，或者需要在计算域入口处精确生成一个非线性波（如斯托克斯五阶波），又或者你的实验数据表明波浪在传播过程中存在某种独特的衰减特性，而标准模型无法描述。这时候，FLUENT的用户自定义函数（User-Defined Function, UDF）就成了打通任督二脉的关键。

这个项目，就是一次深入FLUENT内核的实战：利用C/C++编写UDF，实现边界上的自定义造波与域内的自定义波浪衰减。这不仅仅是调用几个API那么简单，它要求你同时具备流体力学原理、数值方法以及C语言编程的交叉能力。你需要告诉FLUENT：在我的计算域入口，每一时刻每个网格面上的速度和压力（或体积分数）是多少；在域内，如何根据我的理论或经验公式，对波浪的动能或波高进行人为的、可控的衰减，以模拟海绵层、物理阻尼或数值耗散效应。

为什么是C/C++？因为FLUENT的求解器核心是用C写的，UDF作为其扩展，自然需要用C语言接口（尽管FLUENT也支持某些方式的C++编译，但核心仍是C风格）。这个过程，涉及到对FLUENT数据结构（如线程Thread、面face、单元cell）的理解，对求解过程（初始化、迭代步、结束）的把握，以及对并行计算（如果你的案例是并行的）数据交换的认知。这就像给一台精密的发动机编写新的燃油喷射控制程序，你需要知道发动机的工况（迭代步、时间）、每个气缸的位置（网格面坐标），然后精确地注入你的“燃料”（边界条件值）。

2. 核心需求解析：为什么标准模型不够用？

在开始敲代码之前，我们必须厘清需求：我们到底要解决什么标准模型解决不了的问题？这决定了UDF的复杂度和编写方向。

2.1 边界造波的定制化需求

标准FLUENT的波浪模型，如使用“速度入口”或“压力入口”结合波浪理论公式，对于线性波（小振幅波）是方便的。但面临以下情况时，就显得力不从心：

  1. 复杂波浪谱：模拟JONSWAP谱、PM谱等描述的随机波浪场。标准界面难以直接输入一个连续的谱函数并实时生成对应的时域信号。
  2. 高阶非线性波：如斯托克斯二阶、三阶乃至五阶波。这些波的波面方程、水质点速度势表达式复杂，标准输入框无法容纳冗长的公式。
  3. 聚焦波或畸形波：需要在特定时间、特定位置产生一个巨大的波峰，这需要精确控制波幅随时间的变化历程。
  4. 与外部数据耦合：入口边界条件来自另一款软件的计算结果（如SWAN波浪模型）或物理实验测量数据，需要动态读取数据文件并赋值。

UDF通过DEFINE_PROFILE宏，可以让我们在每一个迭代步（或时间步），根据当前仿真时间和面的空间坐标，动态计算并返回一个标量值（如速度分量、压力、体积分数），完美解决上述所有定制化需求。

2.2 波浪衰减的物理与数值需求

在数值波浪水槽中，为了避免波浪在边界反射回计算域干扰流场，我们通常在出口或侧边界设置“消波区”（或称海绵层）。标准方法可能是简单地增加粘性，但这通常不够有效或物理意义不明确。自定义衰减UDF可以更优雅地实现：

  1. 物理阻尼模拟：模拟多孔介质、植被层对波浪能量的耗散。衰减系数可能与水深、波高、频率相关。
  2. 主动吸收式造波：这是更高级的技术。在造波板/边界处，不仅生成入射波，还通过测量临近区域的波面信号，实时计算并叠加一个出射波以抵消反射波，实现“无反射”边界。这需要UDF能够读取计算域内的实时解（如某监测点的波高）。
  3. 数值耗散控制：在特定区域（如自由表面附近）增加人工粘性，以抑制高波陡时可能出现的数值振荡，保持计算稳定。这可能需要通过DEFINE_SOURCE宏为动量方程添加源项。

DEFINE_SOURCE宏允许我们为指定的输运方程（如X动量、Y动量、湍流方程）添加一个源项。对于衰减，我们通常是在动量方程中添加一个与速度方向相反的力源项，其大小与当地速度成正比，即S = -ρ * C * U，其中C是衰减系数，可以是常数，也可以是空间（甚至时间）的函数。

3. 开发环境搭建与UDF编译基础

工欲善其事，必先利其器。FLUENT UDF开发的环境配置是第一个小挑战，尤其是对于习惯了现代IDE（如Visual Studio Code）的开发者来说，回到FLUENT内置的文本编辑器和命令行编译环境会有些复古。

3.1 编译器配置：Visual Studio的核心地位

FLUENT在Windows上依赖于Microsoft Visual C++编译器。你必须安装与你的ANSYS/FLUENT版本匹配的Visual Studio版本。例如，ANSYS 2022 R2通常需要Visual Studio 2019。这不是可选的，是必须的。安装时，务必勾选“使用C++的桌面开发”工作负载，确保MSVC编译器、Windows SDK等组件齐全。

注意：切勿尝试使用MinGW或Cygwin的GCC编译器来编译FLUENT UDF，它们与FLUENT内部的数据结构不兼容，会导致链接错误。FLUENT的UDF编译系统是紧密绑定MSVC的。

安装后，你不需要打开庞大的Visual Studio IDE来写UDF。你可以用任何轻量级文本编辑器（如VS Code、Notepad++、Sublime Text）来编写.c源文件。编译环节由FLUENT在内部调用MSVC完成。

3.2 UDF代码结构与编译流程详解

一个最基本的UDF源文件结构如下：

#include "udf.h" // 必须包含的头文件，定义了所有FLUENT宏和数据结构 DEFINE_PROFILE(my_wave_velocity, thread, position) { real t = CURRENT_TIME; // 获取当前物理仿真时间 real x[ND_ND]; // 用于存储面中心坐标的数组 face_t f; // 面标识符 begin_f_loop(f, thread) // 循环遍历该边界线程上的所有面 { F_CENTROID(x, f, thread); // 获取面f的中心坐标，存入数组x real x_coord = x[0]; // 假设波浪沿x方向传播，x[0]即x坐标 // 根据时间t和坐标x_coord，计算波浪速度u // 例如，线性波水质点水平速度：u = A * g * k / omega * cosh(k*(z+h))/cosh(k*h) * cos(k*x - omega*t) // 这里需要你实现具体的波浪理论公式 real wave_u = ...; // 你的计算逻辑 F_PROFILE(f, thread, position) = wave_u; // 将计算值赋给该面的指定位置（如速度分量） } end_f_loop(f, thread) }

编译流程：

  1. 在FLUENT中，通过Define -> User-Defined -> Functions -> Compiled打开编译对话框。
  2. 添加你的.c源文件。
  3. 点击“Build”。FLUENT会在后台调用nmake（MSVC的命令行构建工具）来编译UDF，生成一个共享库（.dll文件）。
  4. 如果代码有语法错误，FLUENT的TUI（文本用户界面）或一个弹出的控制台窗口会显示编译错误信息。这是调试的第一步，也是最常见的一步。
  5. 编译成功后，点击“Load”将共享库载入当前FLUENT会话。

一个关键的心得：在编写UDF时，务必保持FLUENT案例文件（.cas）的路径和名称不要包含中文或特殊字符，且路径不宜过深。编译生成的临时文件和.dll文件会放在案例文件所在目录，路径复杂有时会引发难以排查的权限或文件访问问题。

4. 边界造波UDF的深度实现

让我们深入一个具体的例子：实现一个二阶斯托克斯波的入口速度边界。

4.1 波浪理论到代码的映射

首先，我们需要二阶斯托克斯波的理论公式。以水深h，波高H，波数k（k=2π/L，L为波长），圆频率ω（ω=2π/T，T为周期）为例。其波面η和水平速度u的近似解为：

η = (H/2)*cos(θ) + (H²k/16)*cosh(kh)*(2+cosh(2kh))/sinh³(kh) * cos(2θ) u = (Hω/2)*cosh(k(z+h))/sinh(kh)*cos(θ) + (3/64)*(H²ωk)*cosh(2k(z+h))/sinh⁴(kh)*cos(2θ) 其中 θ = kx - ωt

注意，z是垂直坐标（向上为正），原点在静水面。在FLUENT中，我们通常将静水面设为z=0。

4.2 DEFINE_PROFILE宏的实战编码

我们的目标是实现上述u的速度剖面。假设波浪沿x正方向传播，入口边界是x=0的平面。

#include "udf.h" #define H 0.1 // 波高 (m) #define T 2.0 // 周期 (s) #define h 1.0 // 水深 (m) #define g 9.81 // 重力加速度 (m/s^2) DEFINE_PROFILE(stokes2nd_u_velocity, thread, position) { real t = CURRENT_TIME; real x[ND_ND]; face_t f; // 根据线性色散关系计算波数和频率 (可预先算好，这里演示动态计算) real omega = 2.0 * M_PI / T; // 线性色散关系: omega^2 = g * k * tanh(k*h) // 需要迭代求解k，这里为了简化，假设已知波长L，则 k = 2*M_PI / L; // 更严谨的做法是预先用脚本算好k，或实现一个简单的迭代求解器。 real L = 5.0; // 示例波长，实际应根据T和h通过色散关系求得 real k = 2.0 * M_PI / L; real A = H / 2.0; // 波幅 begin_f_loop(f, thread) { F_CENTROID(x, f, thread); real z_coord = x[2]; // 假设z是第三个坐标（索引2），根据你的模型设置确认 // 注意：FLUENT中坐标索引：0-x, 1-y, 2-z。确保你的模型坐标系一致。 real theta = k * 0.0 - omega * t; // 入口边界x=0，所以x坐标为0？不，F_CENTROID获取的是面的实际中心坐标。 // 更通用的写法：theta = k * x[0] - omega * t; 但入口边界上x[0]应该是常数（如0）。这里用x[0]更稳妥。 theta = k * x[0] - omega * t; // 双曲函数 real sinh_kh = sinh(k*h); real cosh_kh = cosh(k*h); real sinh_2kh = sinh(2*k*h); real cosh_2kh = cosh(2*k*h); // 一阶速度项系数 real u1_coef = A * omega * cosh(k*(z_coord + h)) / sinh_kh; real u1 = u1_coef * cos(theta); // 二阶速度项系数 (简化公式，不同文献系数略有差异) real u2_coef = (3.0/64.0) * H * H * omega * k * cosh(2*k*(z_coord + h)) / (sinh_kh * sinh_kh * sinh_kh * sinh_kh); real u2 = u2_coef * cos(2*theta); real total_u = u1 + u2; F_PROFILE(f, thread, position) = total_u; } end_f_loop(f, thread) }

关键点解析：

  1. CURRENT_TIME：获取的是当前迭代步对应的物理时间，这对于非定常模拟至关重要。
  2. F_CENTROID：获取每个面的中心坐标。对于入口边界，x[0]通常是固定的（如0），但使用它可以使代码更通用（例如用于斜向波）。
  3. 坐标轴确认：这是最大的坑之一。你必须清楚你的FLUENT模型坐标系：哪个轴是波浪传播方向？哪个轴是垂直方向？x[0],x[1],x[2]分别对应什么？在2D模型中，ND_ND=2，只有x[0]和x[1]。通常2D水槽，x是传播方向，y是垂直方向。代码中的x[2]需要改为x[1]。
  4. 静水面基准：公式中的z坐标原点在静水面，向上为正。你的FLUENT几何中，静水面对应的y坐标值是多少？如果静水面在y=0，那么z_coord就直接是x[1]。如果静水面在y=1.5，那么公式中的(z + h)应替换为(x[1] - 1.5 + h)。坐标转换必须极其小心，否则生成的波会完全不对。

4.3 压力入口与VOF多相流耦合

如果使用VOF方法模拟自由表面，入口边界通常设置为“速度入口”并指定水的体积分数（alpha）分布。这时，你需要另一个DEFINE_PROFILE来定义alpha的分布。

DEFINE_PROFILE(wave_volume_fraction, thread, position) { real t = CURRENT_TIME; real x[ND_ND]; face_t f; real k, omega, eta; // 波数，频率，波面高度 // ... 计算k, omega (同上) ... begin_f_loop(f, thread) { F_CENTROID(x, f, thread); real y_coord = x[1]; // 假设y是垂直轴 real theta = k * x[0] - omega * t; // 计算波面eta (以静水面y=0为基准) eta = A * cos(theta) + ... ; // 加上二阶项 // 根据面的y坐标与波面eta的关系，判断是水还是空气 if (y_coord <= eta) // 如果面中心低于波面，则为水 { F_PROFILE(f, thread, position) = 1.0; // 水的体积分数为1 } else { F_PROFILE(f, thread, position) = 0.0; // 空气的体积分数为0 (水的体积分数为0) } // 注意：这是一种简化的“阶梯”赋值，在波面穿过网格面时不够光滑。 // 更精细的做法可以设置一个过渡层，但会复杂很多。 } end_f_loop(f, thread) }

在FLUENT界面中，你需要将入口边界的“速度”和“体积分数”都设置为“udf”，并分别选择对应的UDF函数名（如stokes2nd_u_velocity和wave_volume_fraction）。

5. 波浪衰减UDF的实现策略

衰减UDF通常通过源项（DEFINE_SOURCE）实现，作用于动量方程，在指定的消波区内施加一个与速度反向的力。

5.1 定义衰减区域与系数

首先，我们需要在FLUENT中通过“Adapt -> Region...”或直接在网格划分时，标记出消波区（例如，计算域最后2米长的区域）。假设我们标记了这个区域为一个“Cell Zone”，并命名为sponge_zone。

我们的UDF需要判断一个单元是否位于该区域，并计算其衰减系数。衰减系数C可以是常数，也可以随进入消波区的深度d（从消波区起点开始算的距离）增加而增大，常用线性或二次函数。

#include "udf.h" DEFINE_SOURCE(momentum_source_x, c, t, dS, eqn) { real source = 0.0; real C = 0.0; real x[ND_ND]; Thread *sponge_thread = NULL; // 1. 获取名为“sponge_zone”的细胞线程指针 // 注意：这需要在FLUENT中提前定义好该区域并命名。 sponge_thread = Lookup_Thread(Get_Domain(1), "sponge_zone"); // Get_Domain(1)获取第一个域 // 2. 判断当前单元c是否属于消波区线程 if (t == sponge_thread) // 如果当前单元线程就是消波区线程 { C_CENTROID(x, c, t); real x_coord = x[0]; // 假设消波区沿x方向布置 // 定义消波区起点和终点坐标 real x_sponge_start = 8.0; // 消波区从x=8m开始 real x_sponge_end = 10.0; // 消波区在x=10m结束（计算域出口） if (x_coord >= x_sponge_start && x_coord <= x_sponge_end) { // 计算归一化的衰减强度，从0到1线性增加 real d = (x_coord - x_sponge_start) / (x_sponge_end - x_sponge_start); real C_max = 5.0; // 最大衰减系数 (1/s)，需要根据网格尺寸和时间步长调试 C = C_max * d * d; // 使用二次函数，末端衰减更强 // 获取当前单元的x方向速度 real u_vel = C_U(c, t); // 计算源项：S = -ρ * C * u source = - C_R(c, t) * C * u_vel; // 为雅可比矩阵dS提供导数项 (可选，但能提高收敛性) // dS[eqn] = ∂S/∂u = -ρ * C dS[eqn] = - C_R(c, t) * C; } } return source; }

关键点解析：

  1. Lookup_Thread：这是一个非常重要的函数，用于通过区域名称获取线程指针。确保在FLUENT中设置的区域名称与代码中的字符串完全一致（包括大小写）。
  2. 判断逻辑：if (t == sponge_thread)是判断当前单元c所在的线程t是否就是消波区线程。这是最直接的判断方法。也可以使用THREAD_ID(t) == THREAD_ID(sponge_thread)。
  3. 源项公式：source = -ρ * C * u。负号表示力与速度方向相反，起阻尼作用。C的量纲是[1/时间]，C越大，衰减越快。
  4. 雅可比项dS：这是可选的，但强烈建议提供。它告诉求解器源项相对于求解变量（这里是速度u）的导数，有助于牛顿迭代法的收敛。对于线性源项S = -K * u，导数dS/du = -K。
  5. 系数C的调试：C_max的值需要调试。太小则衰减效果不足，反射波依然明显；太大则可能使方程刚性过大，导致计算不稳定或发散。通常从较小的值（如0.1~1.0）开始试算，观察消波区末端的速度场是否平稳接近零。

5.2 在FLUENT中设置源项

编写好UDF并编译加载后，在FLUENT中：

  1. 进入Define -> Boundary Conditions。
  2. 选择Cell Zone Conditions，选中你的流体区域（通常是整个计算域）。
  3. 点击Edit...，在弹出的对话框中，找到Momentum选项卡。
  4. 在X Momentum Source Terms或相应的方向动量源项中，选择udf并从下拉列表中选择你定义的momentum_source_x。
  5. 如果你的衰减是各向同性的（也衰减y方向速度），需要为Y Momentum也添加一个类似的源项UDF（公式中的u_vel需改为C_V(c,t)）。

6. 调试技巧与常见问题实录

UDF开发调试过程如同侦探破案，需要耐心和系统的方法。以下是我踩过无数坑后总结的实战经验。

6.1 编译与加载阶段的“拦路虎”

  1. “找不到 udf.h” 或编译错误：

    • 原因：udf.h路径未包含。FLUENT编译时自动设置，但如果你在外部用IDE编译可能会遇到。
    • 解决：永远使用FLUENT内置的编译对话框进行编译。这是最可靠的方式。确保你的.c文件路径无中文、无空格。

  2. “error LNK2001: 无法解析的外部符号”：

    • 原因：这是最常见的链接错误。意味着你的UDF中声明了一个函数（比如DEFINE_PROFILE），但FLUENT的编译环境找不到它的实现。99%的情况是你的UDF代码有语法错误，导致编译器没有生成该函数的对象文件。
    • 解决：仔细查看FLUENT TUI窗口或弹出的控制台中的编译输出信息，而不是最后的“链接错误”。往上翻，通常会有更早的C语法错误提示，比如“missing ';' before 'type'”。修正这些语法错误。

  3. UDF加载成功，但勾选后无效果：

    • 原因A：UDF函数名与FLUENT界面中选择的名称不匹配。区分大小写。
    • 检查：在FLUENT控制台输入define/user-defined/function-hooks可以列出所有已加载的UDF。核对名字。
    • 原因B：边界条件类型设置错误。例如，你的UDF是DEFINE_PROFILE用来定义速度，但你将边界条件类型设为了“压力入口”。
    • 解决：确保边界条件类型（速度入口、压力入口等）与UDF的预期用途一致，并在该边界条件的相应字段（如速度分量、压力、体积分数）中选择UDF。

6.2 运行时逻辑错误排查

当UDF能加载并能被调用，但模拟结果明显不对（如波浪没生成、波形畸变、计算发散），就需要进行运行时调试。

  1. 使用Message宏输出调试信息：

    #if !RP_NODE // 确保只在主机进程上打印，避免并行时每个进程都打印造成刷屏 Message("Time = %f, My calculated value = %f\n", CURRENT_TIME, some_variable); #endif

    将关键变量（如计算出的速度、坐标、衰减系数）打印到FLUENT控制台。这是最直接的调试手段。注意在并行计算时，用#if !RP_NODE包裹，否则每个计算节点都会打印，信息会极多。

  2. 利用外部文件记录数据：

    FILE *fp; fp = fopen("udf_debug.log", "a"); fprintf(fp, "Time: %f, Cell ID: %d, Coord: (%f, %f), Velocity: %f\n", CURRENT_TIME, c->id, x[0], x[1], calculated_velocity); fclose(fp);

    将数据写入日志文件，可以更详细地分析UDF在每个单元、每个时间步的行为。注意：在并行计算中，每个进程都会试图创建/写入同一个文件，会导致冲突。需要为每个进程创建不同的文件名，例如使用PRINCIPAL_HOST_P和MY_PROCESSOR_ID宏。

  3. 检查坐标与物理量单位：

    • “幽灵波”或“反重力波”：大概率是坐标转换错误。反复检查你的波浪理论公式中的z（垂直坐标）与F_CENTROID或C_CENTROID获取的x[1]或x[2]之间的换算关系。画个简单的草图，标出FLUENT全局坐标系原点、静水面位置、水深方向。
    • 量纲不一致：FLUENT内部使用SI单位制（米、秒、千克）。确保你公式中的所有常数（如重力加速度g=9.81 m/s²）和输入参数（波高H、周期T）单位一致。

  4. 并行计算特有问题：

    • UDF只在部分区域生效：在并行计算中，计算域被分割。Lookup_Thread查找的线程可能只在主机（Host）进程上有效。在节点（Node）进程上，该指针可能为NULL。安全的做法是，在初始化宏（如DEFINE_ON_DEMAND或DEFINE_INIT）中查找线程指针，并将其存储在全局变量中，供其他UDF使用。
    Thread *sponge_thread_global = NULL; DEFINE_INIT(my_init, domain) { sponge_thread_global = Lookup_Thread(domain, "sponge_zone"); Message("Sponge zone thread ID found: %d\n", THREAD_ID(sponge_thread_global)); } // 然后在DEFINE_SOURCE中使用 sponge_thread_global
    • 数据不同步：如果你在UDF中修改了某个全局变量（非FLUENT求解变量），这个修改不会自动在其他进程间同步。需要用到FLUENT提供的并行通信宏，如PRF_CSEND_INT等，这属于高级话题，初期尽量规避。

6.3 稳定性与收敛性问题

  1. 源项导致发散：

    • 现象：添加衰减源项后，计算在几个迭代步内就发散。
    • 原因：衰减系数C设置过大，导致源项-ρ*C*U的值巨大，使得动量方程失衡。
    • 解决：大幅减小C_max（比如从5.0降到0.5）。同时，务必提供源项的雅可比项dS[eqn]，这能显著改善含有源项的方程的收敛行为。也可以尝试使用隐式松弛因子。

  2. 造波边界引发初始瞬态冲击：

    • 现象：模拟开始时，入口突然从静止变为一个有限振幅的波浪，产生一个非物理的冲击波在域内传播。
    • 解决：在造波UDF中实现一个“缓启动”函数。让波浪的振幅在最初几个周期内从0逐渐增加到目标值。
    real ramp_time = 2.0 * T; // 缓启动时间，例如2个波周期 real ramp_factor; if (t < ramp_time) { ramp_factor = 0.5 - 0.5 * cos(M_PI * t / ramp_time); // 使用余弦函数平滑过渡 } else { ramp_factor = 1.0; } real actual_A = A * ramp_factor; // 将缓启动因子应用到波幅上

  3. 时间步长与波长的匹配：

    • 经验法则：每个波周期内至少要有50-100个时间步，才能较好地解析波浪运动。即Δt <= T / 50。同时，库朗数（CFL）条件也必须满足，对于VOF模拟，通常要求更小的时间步。

7. 从理论到实践：一个完整案例的搭建思路

假设我们要模拟一个长20米，深1米的水槽，水深0.6米，生成一个波高0.1米，周期1.2秒的波浪，并在最后3米设置消波区。

  1. 前处理（SpaceClaim/DesignModeler + Meshing）：

    • 创建2D矩形，长20m，高1m。
    • 划分结构化网格。在自由水面附近（y=0附近）进行局部加密，垂直方向至少布置20层网格以分辨波面。水平方向网格尺寸应小于波长的1/20。
    • 命名边界：velocity_inlet（左边界），pressure_outlet（右边界），bottom_wall（下边界），top（上边界，设为压力出口以模拟大气）。
    • 创建一个名为sponge的Face Zone（用于后续标记消波区单元），覆盖x从17m到20m的区域。

  2. FLUENT设置：

    • 启用瞬态求解器。
    • 启用多相流模型（VOF），主相为水，次相为空气。
    • 操作密度设为水的密度（减轻浮力计算负担）。
    • 重力加速度y方向设为-9.81。
    • 边界条件：
      • velocity_inlet：速度指定方法选“udf”，选择stokes2nd_u_velocity；体积分数选“udf”，选择wave_volume_fraction。
      • pressure_outlet：回流体积分数设为“air”（即次相），防止水从出口回流。
      • top：设为压力出口，回流体积分数也为“air”。
    • 初始化：使用标准初始化，从velocity_inlet补丁初始化。
    • 动网格？本例不需要，是固定边界造波。

  3. UDF准备：

    • 将前面章节的造波UDF和衰减源项UDF写在一个或两个.c文件中。
    • 在DEFINE_INIT宏中，查找并存储sponge区域的线程指针到全局变量。
    • 修正所有坐标索引和单位。
    • 编译并加载。

  4. 计算与监控：

    • 设置时间步长，例如0.01秒（满足T/120）。
    • 设置总物理时间，例如20秒（约16个波周期）。
    • 创建监测点：在消波区前（如x=15m， y=0m）和消波区后（x=19.5m， y=0m）设置波高监测点。
    • 计算并观察：入射波是否稳定生成？消波区后的波高是否显著减小（理想情况接近零）？计算域中部（x=10m）的波高时程曲线是否稳定、周期性良好？

这个过程充满了试错。第一次运行很可能不成功，需要你结合第6章的调试技巧，反复检查UDF逻辑、模型设置和网格质量。当看到稳定的波浪从入口生成，平滑地传播，并在消波区逐渐消失，几乎没有反射时，那种成就感是对所有调试工作的最好回报。这不仅仅是完成了一个CFD模拟，更是真正意义上将物理理论、数值方法和编程实践融会贯通的一次深度实战。