当前位置：首页 > news >正文

OpenFOAM实战：在interFoam中植入多孔介质源项模拟复杂固壁

news 2026/5/12 9:42:41

1. 多孔介质模拟的工程需求与原理

在流体力学仿真中，我们经常遇到需要处理复杂几何边界的情况。传统方法是通过精细的网格划分来精确描述固体边界，但这会带来两个主要问题：一是计算成本急剧上升，二是对于动态变化的边界（比如相变、侵蚀等过程）难以处理。这时候，多孔介质模型就提供了一种巧妙的替代方案。

我曾在船舶螺旋桨空化模拟项目中遇到过类似问题。当需要模拟空泡溃灭对金属表面的冲击时，直接描述微观尺度的金属结构几乎不可能。后来发现通过将金属区域设置为高阻力系数的多孔介质，可以完美规避这个难题。这种方法的本质是在动量方程中添加与速度成正比的阻力源项（S_i = -D U_i），当D取值足够大时（比如1e10），该区域流体速度会被强制归零，等效于固体边界。

OpenFOAM的interFoam求解器原生支持VOF（Volume of Fluid）方法，能够精确追踪气液界面。我们通过修改其源代码，将其中一相（通常是液相）的部分区域定义为多孔介质，就能实现用流体网格模拟固体边界的效果。这种混合方法结合了VOF的界面捕捉能力和多孔介质的边界简化优势，特别适合处理以下场景：

具有复杂表面纹理的固体边界（如多孔材料、粗糙表面）
动态变化的边界形状（如融化、沉积过程）
需要快速原型验证的初期设计方案

2. interFoam求解器的定制化改造

2.1 基础环境准备

首先需要获取interFoam的源代码。在OpenFOAM 9中，求解器位于/opt/openfoam9/applications/solvers/multiphase/interFoam。我建议新建一个工作目录，将整个interFoam文件夹连同Make目录一起复制过去。同时准备一个测试案例，比如经典的溃坝案例（damBreak），可以从教程目录/opt/openfoam9/tutorials/multiphase/interFoam/laminar/获取。

目录结构应该如下：

├── damBreak # 测试案例 ├── Make # 编译规则 ├── alphaSuSp.H # 以下是interFoam源代码 ├── correctPhi.H ├── createFieldRefs.H ├── createFields.H ├── initCorrectPhi.H ├── interFoam.C ├── pEqn.H ├── rhofs.H └── UEqn.H

提示：建议使用wmake命令编译时，先执行wmRefresh更新编译环境。我在Ubuntu 20.04上测试时，曾因环境变量问题导致编译失败。

2.2 关键场变量定义

在createFields.H中添加自定义的key场，这个标量场将作为多孔介质区域的标记。以下是具体代码实现：

// 定义单位制 (kg*m^-3*s^-1) const dimensionSet mydimension(1, -3, -1, 0, 0, 0, 0); // 创建key场，初始值与alpha1相同 volScalarField key ( IOobject ( "key", runTime.timeName(), mesh, IOobject::NO_READ, // 不从文件读取 IOobject::AUTO_WRITE // 自动写入结果 ), alpha1 );

这段代码做了三件事：

定义了阻力系数的量纲（kg/m³/s）
创建了与alpha1（VOF相分数）同维度的标量场
设置该场不需要从文件读取，但会输出到结果中便于调试

在实际项目中，我发现将key场初始化为alpha1的副本是个好习惯，这样能确保多孔介质区域与液相区域初始重合。当然，后续可以通过条件判断来动态调整key场的分布。

3. 动量方程改造与源项植入

3.1 动态更新key场

在UEqn.H中，我们需要在求解动量方程前先更新key场。我的经验是添加在MRF校正之后：

MRF.correctBoundaryVelocity(U); // 根据alpha1更新key场 for (label cellI=0; cellI<mesh.C().size(); cellI++) { if(alpha1[cellI] > 0.001) { // 阈值可调 key[cellI] = alpha1[cellI]; } else { key[cellI] = 0; } }

这里使用0.001作为阈值是为了避免数值振荡带来的影响。在模拟多孔过滤器的项目中，我发现这个值能很好平衡数值稳定性和边界清晰度。

3.2 隐式源项添加

接下来修改动量方程，添加多孔介质源项。关键是要使用fvm::Sp函数实现隐式离散化：

fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + MRF.DDt(rho, U) + turbulence->divDevTau(rho, U) == phaseChange.SU(rho, rhoPhi, U) + fvModels.source(rho, U) - fvm::Sp(dimensionedScalar("tmp", mydimension, 1e10)*key, U) );

几个技术细节值得注意：

fvm::Sp会自动将源项线性化为隐式形式，显著改善收敛性
系数1e10需要根据具体问题调整，太大会导致矩阵病态，太小则不能有效抑制速度
key场作为乘数，实现了只在指定区域添加源项

在模拟船舶防污涂层时，我发现将阻力系数设为1e8~1e12范围内都能获得不错的效果，具体取决于网格尺寸和时间步长。

4. 案例配置与参数优化

4.1 输运性质设置

由于我们将"水"相区域实际当作固体使用，需要修改constant/transportProperties：

phases (water air); water { transportModel Newtonian; nu 1.48e-05; // 空气的运动粘度 rho 1; // 空气的密度 } air { transportModel Newtonian; nu 1.48e-05; rho 1; } sigma 0; // 关闭表面张力

这种设置看似违反直觉，但实际上是因为：

多孔介质区域内的流体运动已被源项抑制
将物性设为与空气相同可以避免密度差引起的虚假流动
表面张力关闭是因为气固界面不需要考虑这一效应

4.2 初始条件设置

在system/setFieldsDict中定义多孔介质区域：

defaultFieldValues ( volScalarFieldValue alpha.water 0 ); regions ( boxToCell { box (0.292 0.536 -1) (0.315999 0.584 1); fieldValues ( volScalarFieldValue alpha.water 1 ); } );

这个例子使用长方体区域，但在实际项目中，我经常用更复杂的几何定义方式：