简介sdfibm 是一套基于 C 与 OpenFOAM v6 实现的浸入边界方法求解器面向从事流固耦合、多颗粒流数值模拟的研究生、科研人员与工程开发者用于在复杂几何与运动边界条件下开展 CFD 仿真。资源包共 153 个文件约 28.24MB以 44 个 .h 头文件与 12 个 .cpp 源文件构成核心求解器代码并配套 blockMeshDict、controlDict、fvSchemes、fvSolution、transportProperties 等 OpenFOAM 算例字典以及 Makefile 编译脚本、sh 运行脚本和 gif、png 演示图覆盖从网格生成、求解设置到结果展示的完整流程。目前已有 787 人学习下载。读者可据此获得可直接编译运行的求解器源码与示例算例理解符号距离场与浸入边界耦合的实现思路掌握在笔记本或 HPC 集群上通过 git clone、make 三步完成部署的方法并借助字典文件快速搭建自己的多颗粒流仿真案例。1. 从一张“卡在网格里”的圆柱绕流图说起如果你做过浸入边界法Immersed Boundary Method, IBM大概率经历过这种场景几何模型明明是个光滑圆柱可一映射到笛卡尔网格上边界就变成了锯齿状的台阶。力项在网格节点上跳来跳去算出来的阻力系数振荡得让人怀疑人生。sdfibm 这个项目就是冲着这个痛点来的——它把符号距离场Signed Distance Field, SDF和 OpenFOAM 的有限体积求解器缝在一起用 SDF 提供亚网格精度的边界距离和法向让 IBM 的力项不再依赖“网格是否恰好切到边界”。项目主体用 C 写成核心是一个 OpenFOAM 的 solver 加上一套 SDF 采样与插值工具。适合谁适合已经跑过 icoFoam 或 pimpleFoam、想往流固耦合或动边界方向走又不想从零手写 IBM 力项的人。下面我按“它怎么算 → 怎么编译跑起来 → 怎么改”的顺序把这份源码拆一遍。2. 架构拆解SDF 与 IBM 力项是怎么接进 OpenFOAM 的2.1 顶层目录与模块职责拿到源码包先别急着wmake。顶层目录大致是这么分的src/下放 SDF 的构造与查询类solvers/下是 OpenFOAM 的求解器入口cases/里是几个能直接跑的算例scripts/里是网格生成和场初始化的辅助脚本。这个分法很“OpenFOAM 原生”——求解器只负责组装方程SDF 作为独立库被链接进去算例目录里用system/controlDict挂载。我一般会先看solvers/里的主文件因为那里决定了整个时间步的推进逻辑。sdfibm 的求解器继承自fvCFD那套基础类在动量方程组装之前插入了一个sdfForce的源项。这个源项不是简单的fvm::SuSp而是显式地把 IBM 力加到UEqn的源项里因为力的大小依赖上一时刻的边界位置和 SDF 梯度。2.2 SDF 类的接口设计SDF 部分的核心是一个sdfField类它对外暴露三个关键方法distance(point)返回点到最近边界的符号距离normal(point)返回该点的单位外法向gradient(point)返回距离场的梯度。这三个方法在 IBM 力项里被反复调用所以实现上做了空间哈希加速——把边界三角形按包围盒分桶查询时只遍历邻近桶。// src/sdf/sdfField.C scalar sdfField::distance(const point p) const { // 先查哈希桶拿到候选三角形索引 labelList candidates hash_.query(p, searchRadius_); scalar minDist GREAT; forAll(candidates, i) { const triangle tri tris_[candidates[i]]; // 点到三角形的最短距离带符号 scalar d tri.signedDistance(p); if (mag(d) mag(minDist)) minDist d; } return minDist; }这段代码里searchRadius_是个可调参数默认取网格最小边长的两倍。如果几何特别细长这个半径要手动调大否则会漏掉邻近三角形导致距离场出现“空洞”。signedDistance的符号由三角形法向和点相对位置决定法向朝外为正朝内为负——这个约定在 IBM 力项里必须一致否则力会反向。2.3 IBM 力项的组装逻辑力项的计算分两步先把流体速度插值到边界点再根据边界点与流体的相对速度算反馈力。sdfibm 用的是直接力法direct forcing力的大小正比于速度差除以时间步长再乘上一个松弛因子。// solvers/sdfIBM/UEqn.H 片段 forAll(boundaryPoints, i) { point bp boundaryPoints[i]; vector Uf interpolateU(U, bp); // 流体速度插值到边界点 vector Ub boundaryVelocity(bp); // 边界自身运动速度 vector dU Ub - Uf; // 反馈力rho为密度dt为时间步relax为松弛因子 vector f rho * dU / dt * relax; // 把力扩散到周围网格节点 spreadForce(f, bp, forceField); }interpolateU用的是三线性插值边界点落在哪个网格单元就取周围八个节点加权。spreadForce是反向操作把边界点上的力按核函数分配到邻近节点。核函数的选择直接影响力的光滑程度——项目里默认用余弦核宽度取两倍网格间距。如果算出来力场有高频振荡可以把核宽度调到三倍网格间距代价是边界附近的有效分辨率会下降。2.4 时间推进与耦合方式sdfibm 用的是显式耦合每个时间步先根据当前边界位置算 SDF再算 IBM 力然后组装动量方程最后用 PISO 或 PIMPLE 算法修正压力。边界位置更新在压力修正之后用的是边界点上的力积分得到的速度。这种显式方式对时间步长有要求CFL 数一般要压到 0.5 以下否则边界附近的力会发散。注意如果算例里固体密度远大于流体密度显式耦合会非常僵硬时间步得压到很小。常见做法是改用隐式或半隐式耦合但 sdfibm 的源码里没有内置需要自己改。3. 编译与跑通从 wmake 到第一个圆柱绕流算例3.1 环境依赖与编译顺序sdfibm 依赖 OpenFOAM 的 foundation 版本建议用 OpenFOAM 8 或 9因为源码里用到了fvCFD的新接口。编译前先确认$WM_PROJECT_DIR指向正确然后按顺序编先编src/下的 SDF 库再编solvers/下的求解器。# 先编 SDF 库 cd src/sdf wmake libso # 再编求解器 cd ../../solvers/sdfIBM wmake如果wmake报找不到sdfField.H说明库没编成功或者LD_LIBRARY_PATH没包含$FOAM_USER_LIBBIN。我一般会在~/.bashrc里加一行export LD_LIBRARY_PATH$FOAM_USER_LIBBIN:$LD_LIBRARY_PATH省得每次手动设。3.2 算例目录结构与关键字典cases/下有个cylinderRe100的算例可以直接拿来验证。目录结构是标准的 OpenFOAM 布局0/放初始场constant/放网格和物性system/放控制字典。多出来的东西是constant/triSurface/下的圆柱 STL 文件以及system/sdfDict这个 SDF 专用字典。sdfDict里几个关键参数参数含义典型值surfaceFile边界几何文件路径cylinder.stlsearchRadius哈希查询半径2 * minCellSizekernelWidth力扩散核宽度2 * minCellSizerelax反馈力松弛因子0.5relax取 1.0 时力最“硬”边界能严格贴合但容易振荡取 0.3 到 0.5 之间比较稳代价是边界会有轻微“渗透”。我一般先用 0.5 跑通再根据残差曲线微调。3.3 运行与监控跑算例就是标准的sdfIBM命令cd cases/cylinderRe100 sdfIBM log.run 21 监控残差用foamMonitor看postProcessing/residuals/0/residuals.dat。如果Ux残差在 1e-3 附近震荡不降多半是relax太大或者时间步太大。先看Co数如果超过 0.8把deltaT砍半。如果Co没问题但残差还是震荡把relax降到 0.3 再试。提示第一次跑建议先用potentialFoam初始化流场再跑sdfIBM能省不少时间。sdfibm 的算例里没有内置这一步需要手动在0/之前插一个potentialFoam的 run。3.4 后处理与验证算完之后用foamToVTK导出在 ParaView 里看涡量场。圆柱尾迹的卡门涡街应该清晰可见如果边界附近出现明显的力“条纹”说明核宽度太小或者searchRadius不够。阻力系数可以用forces函数对象监控在system/controlDict里加一段functions { forces { type forces; libs (libforces.so); patches (cylinder); rho rhoInf; rhoInf 1.0; CofR (0 0 0); } }Re100 的圆柱阻力系数应该在 1.3 到 1.4 之间升力系数幅值在 0.3 左右。如果差得远先检查 STL 的单位——很多 STL 是毫米而 OpenFOAM 默认是米差一千倍。4. 避坑与排查跑 sdfibm 时最容易翻车的五个地方4.1 现象编译通过但运行时报“segmentation fault”原因SDF 库的哈希表在构造时没有正确初始化或者 STL 文件读入后三角形数量为零。常见于 STL 文件路径写错但 OpenFOAM 不会报“文件不存在”而是静默返回空列表。解决在sdfDict里把surfaceFile写成绝对路径或者在sdfField构造后加一行Info triangles: tris_.size() endl;确认读入数量。如果为零检查 STL 是不是二进制格式——sdfibm 的读入器只支持 ASCII STL。4.2 现象边界附近的力场出现棋盘状振荡原因力扩散核宽度小于网格间距导致力只落在单个节点上形成奇偶失稳。解决把kernelWidth调到至少2 * minCellSize如果还振荡就调到3 * minCellSize。代价是边界会变“模糊”但稳定性大幅提升。4.3 现象残差在 1e-2 附近卡住不降原因IBM 力项是显式加的如果relax太大力项会主导方程压力修正跟不上。解决先把relax降到 0.3如果还不行把deltaT降到原来的四分之一。另外检查fvSchemes里div(phi,U)是不是用了Gauss linearUpwind如果是Gauss linear改成linearUpwind能改善稳定性。4.4 现象边界点“穿透”固体流体进入圆柱内部原因searchRadius太小SDF 在圆柱内部某些点算出了正值导致力方向反了。解决把searchRadius调到3 * minCellSize并在sdfField::distance里加一个断言如果点在 STL 包围盒内部但距离为正打印警告。这个断言在调试阶段非常有用。4.5 现象并行跑的时候结果和单核不一致原因SDF 的哈希表在并行时没有做域分解每个核都持有完整的三角形列表但边界点插值时跨核通信没做对。解决sdfibm 的并行支持不完整建议先用单核跑通再考虑并行。如果一定要并行把decomposeParDict里的method改成hierarchical并且把searchRadius调大减少跨核查询。5. 二次开发把 sdfibm 改成你自己的动边界求解器5.1 替换边界运动逻辑sdfibm 默认的边界速度是零也就是静止固体。要改成动边界最直接的办法是继承boundaryMotion类重写velocity(point)方法。比如让圆柱做横向振荡// myMotion.H class oscillatingCylinder : public boundaryMotion { scalar A_; // 振幅 scalar omega_; // 角频率 public: virtual vector velocity(const point p) const { // 只给x方向速度y和z为零 return vector(A_ * omega_ * cos(omega_ * runTime.value()), 0, 0); } };然后在sdfDict里加一行motionType oscillatingCylinder;并在求解器初始化时根据这个字段构造对应的运动对象。注意runTime.value()返回的是当前物理时间如果算例里deltaT是自适应的这个速度会跳变最好用runTime.timeName()对应的整步时间。5.2 加入流固耦合的隐式迭代显式耦合的稳定性瓶颈在固体密度大时特别明显。一个可行的改法是每个时间步内做两到三次子迭代先算流体得到边界上的力用这个力更新固体速度再用新速度重算流体。sdfibm 的源码里UEqn.H和pEqn.H是分开的可以在两者之间插一个solidUpdate.H。// 在pEqn.H之后下一时间步之前 for (int subIter 0; subIter 3; subIter) { vector F computeForceOnBoundary(); vector Unew Uold F / mass * dt; updateBoundaryVelocity(Unew); // 重新组装动量方程并求解 solve(UEqn -fvc::grad(p)); }子迭代次数不是越多越好三次以上收敛就饱和了反而增加计算量。判断收敛可以用边界速度的变化量如果两次子迭代之间速度差小于 1e-4就提前退出。5.3 验证方法用静止圆柱的解析解对表改完之后怎么知道对不对最稳的办法是拿静止圆柱的阻力系数对表。Re100 时文献值大约在 1.32 到 1.38 之间。如果你的结果在这个区间外先别怀疑算法检查网格分辨率——圆柱直径上至少要有 20 个网格否则 IBM 的亚网格精度优势发挥不出来。另一个验证是看边界附近的压力场。在圆柱前缘点压力系数应该接近 1.0在尾缘点应该接近 -0.5。如果前缘压力系数只有 0.5说明力扩散核太宽把压力抹平了。5.4 一个我踩过的坑有一次我把relax设成 1.0想看看“硬”边界的效果。结果跑了 200 步之后圆柱尾迹里出现了非物理的涡阻力系数直接跳到 2.0。查了半天才发现relax1.0时力项的时间尺度小于网格的声学时间尺度压力修正根本来不及响应。从那以后我每次改relax都强制走一遍“先 0.3 跑 100 步看残差曲线再逐步加到 0.5”的流程再也不敢直接上 1.0 了。希望这份拆解能帮你把 sdfibm 跑起来并且敢动手改它的力项和边界运动逻辑。源码包里那几个算例是最好的起点先让圆柱转起来再换成你自己的几何。本文还有配套的精品资源点击获取
