[OpenFOAM 教程] 后步台阶流动分析 — 利用 simpleFoam 捕捉再循环区域

下游1.5 m按网格H是15H。若按注释里的H=0.01 m来读，同一段就变成150H。

后向台阶紧后面，流向速度会变成负的。这段区间结束的 x 就是再附着点。simpleFoam 把这个点记成几个 H，只有在网格高度和 nu 注释用的是同一个 H 时才读得动。

Re 滑块从 50 拉到 1500，台阶后面的逆流区会变长。粒子数滑块（50–500）只改变填进这块区域的轨迹条数。

Re 400입자 220

스텝 직후 재순환 영역(Lr)이 형성된다. Re가 커지면 Lr가 길어지다가 Re ≳ 1000부근부터 비정상이 된다 — Armaly et al. 1983 실험과 일치.

画面下方的 Lr 会写成几倍台阶高。Re 停在低端时，泡贴在台阶上；滑到顶端，它就往下游拖。

前提条件

OpenFOAM v2312 或更高。.com 和 .org 上这套算例的目录是一样的。
终端，以及算例的三格（0、constant、system）。
图在 ParaView 5.x 里打开。

第 1 步：创建案例目录#

在空目录里放上 0、constant、system。

# 加载 OpenFOAM 环境 (取决于安装路径)
source /opt/openfoam2312/etc/bashrc
 
# 复制教程目录或直接创建
mkdir -p backwardFacingStep
cd backwardFacingStep

目录树应该分成这样。

tree backwardFacingStep/

backwardFacingStep/
├── 0/
│   ├── U          # 速度初始/边界条件
│   ├── p          # 压力初始/边界条件
│   ├── k          # 湍动能
│   └── epsilon    # 湍流耗散率
├── constant/
│   ├── transportProperties
│   └── turbulenceProperties
└── system/
    ├── blockMeshDict
    ├── fvSchemes
    └── fvSolution

第 2 步：网格生成 — blockMeshDict#

system/blockMeshDict 把台阶切成两块。顶点的 y 间隔是 1，scale 是 0.1，所以注释里的 H = 0.1 m 和这张网格是一致的。

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2312                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
scale   0.1;   // 整体比例 (m 单位 → 乘以 0.1 得到最终尺寸)
 
vertices
(
    // 入口通道块 (台阶上方)
    (0  1  0)  // 0
    (5  1  0)  // 1
    (5  2  0)  // 2
    (0  2  0)  // 3
    (0  1  1)  // 4
    (5  1  1)  // 5
    (5  2  1)  // 6
    (0  2  1)  // 7
 
    // 台阶下游块 (台阶全高)
    (5  0  0)  // 8
    (20 0  0)  // 9
    (20 2  0)  // 10
    (5  2  0)  // 11 (= vertex 2)
    (5  0  1)  // 12
    (20 0  1)  // 13
    (20 2  1)  // 14
    (5  2  1)  // 15 (= vertex 6)
);
 
blocks
(
    // 入口通道
    hex (0 1 2 3 4 5 6 7)   (50 10 1) simpleGrading (1 1 1)
    // 台阶下游 (更密的网格)
    hex (8 9 10 11 12 13 14 15) (150 20 1) simpleGrading (1 1 1)
);
 
boundary
(
    inlet
    {
        type patch;
        faces ((0 3 7 4));
    }
    outlet
    {
        type patch;
        faces ((9 10 14 13));
    }
    wall
    {
        type wall;
        faces
        (
            (0 1 5 4)   // 上端入口壁
            (1 8 12 5)  // 台阶垂直面
            (8 9 13 12) // 下端底部
            (3 2 6 7)   // 上端壁 (入口)
            (11 10 14 15) // 上端壁 (下游)
        );
    }
    frontAndBack
    {
        type empty;  // 2D 分析
        faces
        (
            (0 1 2 3)
            (4 5 6 7)
            (8 9 10 11)
            (12 13 14 15)
        );
    }
);

# 执行网格生成
blockMesh

预期输出:

Create time
 
Creating block mesh from "system/blockMeshDict"
...
Writing polyMesh with 3600 cells
End

第 3 步：设置边界条件#

0/U — 速度边界条件#

入口是 x 方向 1 m/s，壁面是 noSlip，出口是 zeroGradient。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
dimensions      [0 1 -1 0 0 0 0];  // m/s
 
internalField   uniform (0 0 0);   // 初始值：静止状态
 
boundaryField
{
    inlet
    {
        type        fixedValue;
        value       uniform (1 0 0);  // x 方向 1 m/s 均匀流入
    }
    outlet
    {
        type        zeroGradient;  // 出口：自由流出
    }
    wall
    {
        type        noSlip;        // 壁面：无滑移条件
    }
    frontAndBack
    {
        type        empty;
    }
}

0/p — 压力边界条件#

这里的压力是运动压力。只有出口固定为 0。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
 
dimensions      [0 2 -2 0 0 0 0];  // 运动压力 (p/rho), m²/s²
 
internalField   uniform 0;
 
boundaryField
{
    inlet
    {
        type        zeroGradient;
    }
    outlet
    {
        type        fixedValue;
        value       uniform 0;  // 出口参考压力 0
    }
    wall
    {
        type        zeroGradient;
    }
    frontAndBack
    {
        type        empty;
    }
}

0/k — 湍动能#

入口湍流强度 5% 时， $k = 1.5 \times (0.05)^2 = 0.00375$ ，和文件里写的值一样。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      k;
}
 
dimensions      [0 2 -2 0 0 0 0];  // m²/s²
 
// 入口湍流强度 5%, 流速 1 m/s → k = 1.5*(I*U)² = 1.5*(0.05*1)² ≈ 0.00375
internalField   uniform 0.00375;
 
boundaryField
{
    inlet
    {
        type        fixedValue;
        value       uniform 0.00375;
    }
    outlet
    {
        type        zeroGradient;
    }
    wall
    {
        type        kqRWallFunction;
        value       uniform 0.00375;
    }
    frontAndBack
    {
        type        empty;
    }
}

第 4 步：设置物性参数和湍流模型#

constant/transportProperties#

ν 写在这个文件里。凑出 Re 1000 的那个 H 是否等于网格高度，放到结果表后面再对。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      transportProperties;
}
 
transportModel  Newtonian;
nu              1e-5;   // 运动粘度系数 (m²/s), Re = U*H/nu = 1*0.01/1e-5 = 1000

constant/turbulenceProperties#

RAS 打开，模型用 kEpsilon。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      turbulenceProperties;
}
 
simulationType  RAS;
 
RAS
{
    RASModel        kEpsilon;   // k-epsilon 二方程模型
    turbulence      on;
    printCoeffs     on;
}

第 5 步：设置数值方案和求解器#

system/fvSchemes#

时间项是 steadyState，速度对流是 linearUpwind，k 和 epsilon 用 upwind。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
 
ddtSchemes      { default steadyState; }  // 稳态分析
 
gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;  // 梯度限制器
}
 
divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);  // 二阶精度
    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}
 
laplacianSchemes { default Gauss linear corrected; }
 
interpolationSchemes { default linear; }
 
snGradSchemes { default corrected; }

system/fvSolution#

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
 
solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-7;
        relTol          0.01;
    }
    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
    }
    "(k|epsilon)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
    }
}
 
SIMPLE
{
    nNonOrthogonalCorrectors    0;
    pRefCell                    0;
    pRefValue                   0;
 
    residualControl
    {
        p           1e-4;
        U           1e-4;
        "(k|epsilon)" 1e-4;
    }
}
 
relaxationFactors
{
    fields      { p 0.3; }           // 压力松弛因子 (防止发散)
    equations
    {
        U       0.7;
        k       0.7;
        epsilon 0.7;
    }
}

第 6 步：执行分析#

压力松弛 0.3，加上 U、k、epsilon 的 0.7，就是把日志在 500 步带到残差 1e-4 以下的那组耦合。这种压力插值怎么把格子振荡藏起来，写在 SIMPLE 与 Rhie–Chow 里。

# 检查网格质量 (可选，建议执行)
checkMesh
 
# 运行 simpleFoam
simpleFoam 2>&1 | tee log.simpleFoam

收敛日志:

Time = 1
 
smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0234, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0.412, Final residual = 0.00891, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.987, Final residual = 0.00821, No Iterations 15
...
 
Time = 500
 
smoothSolver:  Solving for Ux, Initial residual = 1.23e-05, ...
GAMG:  Solving for p, Initial residual = 8.91e-05, ...
 
SIMPLE solution converged in 500 iterations

停机标准是每个变量的残差都低于 1e-4。这份日志在第 500 步跨过那条线。

第 7 步：结果可视化 (ParaView)#

用 paraFoam 打开。

# 创建 ParaView 用的空文件并运行
touch backwardFacingStep.foam
paraFoam

ParaView 里的顺序是这样。

在 Coloring 里选 U，用 Magnitude 涂速度大小。
用 Glyph 撒出矢量，看逆流区。
用 Stream Tracer 画流线。

转成 VTK 之后只取数字，用这条命令。

foamToVTK -latestTime

把近壁 ux 变号的 x 当成再附着。脚本把 H 钉在 0.01 m。

import numpy as np
 
# 从 foamToVTK 输出文件中读取底部壁面 (y=0) 数据
# (实际路径请根据环境调整)
data = np.loadtxt("VTK/backwardFacingStep_500/internal.vtk",
                  skiprows=10)  # 排除 VTK 文件头
 
x = data[:, 0]
ux = data[:, 3]  # x 方向速度
 
# 过滤 y ≈ 0 的单元 (底部壁面附近)
wall_mask = data[:, 1] < 0.005
x_wall = x[wall_mask]
ux_wall = ux[wall_mask]
 
# 符号变化位置 = 再附着点
sign_changes = np.where(np.diff(np.sign(ux_wall)))[0]
if len(sign_changes) > 0:
    x_reattach = x_wall[sign_changes[0]]
    H = 0.01  # 台阶高度 (m)
    print(f"再附着长度: {x_reattach/H:.1f}H")

常见错误及解决办法

错误 1: `FOAM FATAL ERROR: Cannot find file "0/epsilon"`#

--> FOAM FATAL ERROR:
    cannot find file
    "backwardFacingStep/0/epsilon"

缺了 0/epsilon。按 k 的骨架补一个。

// 0/epsilon
dimensions      [0 2 -3 0 0 0 0];  // m²/s³
 
// epsilon = C_mu^(3/4) * k^(3/2) / L, L ≈ 0.1*H
internalField   uniform 0.0009;
 
boundaryField
{
    inlet       { type fixedValue; value uniform 0.0009; }
    outlet      { type zeroGradient; }
    wall        { type epsilonWallFunction; value uniform 0.0009; }
    frontAndBack { type empty; }
}

错误 2: 分析发散 (Floating point exception)#

Floating point exception (core dumped)

松弛太猛，或者网格太稀，计算就会发散。把 fvSolution 里的松弛系数调低。

relaxationFactors
{
    fields      { p 0.2; }   // 0.3 → 0.2
    equations   { U 0.5; }   // 0.7 → 0.5
}

错误 3: 残差不下降而是振荡#

nNonOrthogonalCorrectors 还是 0，网格又很不正交时，残差会来回晃。把修正次数加上去。

SIMPLE
{
    nNonOrthogonalCorrectors    2;  // 0 → 2
    ...
}

预期结果汇总

项目	数值
雷诺数 ( $Re = UH/\nu$ )	1,000
再附着长度 ( $x_r / H$ )	6–8
最大反向速度	$\approx -0.1$ m/s
收敛迭代次数	300–600 次

Armaly et al. 1983 在 $Re = 800$ 时约为 6H，就挨在表里的 6–8 旁边。

x_r \approx 6\text{–}8 \cdot H \quad (Re \sim 10^3)

顶点里台阶的 y 从 0 到 1。乘上 scale 0.1，H 就是 0.1 m。nu 注释写成 $Re = 1 \times 0.01 / 10^{-5} = 1000$ ，等于把 H 记成 0.01 m。表里的 1000 和这句注释一致。把网格的 H 代进同一个式子，得到 $1 \times 0.1 / 10^{-5} = 10^{4}$ 。表上的 6–8H，用 0.1 m 是 0.6–0.8 m，用 0.01 m 是 0.06–0.08 m。

下游块的 x 从 5 到 20。台阶面在 0.5 m，出口在 2.0 m，中间 1.5 m。按网格 H 是 15H，0.6–0.8 m 落在这段里面。用注释的 H 去除，同一段变成 150H，6H 只在台阶面下游 0.06 m，坐标是 0.56 m。提取脚本打印时用的 H 也是 0.01 m。底部格子高度是 $(2/20) \times 0.1 = 0.01$ m，只留 y 低于 0.005 的格子，拿到的是这一格的一半。

k 与注释里的 $1.5 \times (0.05)^2 = 0.00375$ 一致。把 $C_\mu = 0.09$ 代进 epsilon 的式子， $C_\mu^{3/4} k^{3/2} \approx 3.77 \times 10^{-5}$ ，再除以写下的 0.0009，混合长度是 0.042 m。 $L \approx 0.1H$ 要么是 0.01 m，要么是 0.001 m，对应的 epsilon 应是 0.00377 或 0.0377。文件里的 0.0009 两个都不是。逆流 $-0.1$ m/s 是入口 1 m/s 的 10%，是湍流强度 5% 所对应的 0.05 m/s 的两倍。

把6H拉长十倍的那个台阶高度#

若留下网格的 H，表里的 Re 1000 就是假的。要把 nu 提高到 $10^{-4}$ ，这个 1000 才会回来。若留下注释的 H，scale 就得降到 0.01，6H 缩成 0.06 m，贴在下游 1.5 m 的最前头。

该把哪一个高度写成这个算例的 H？把入口 1000 读成圆管转捩的雷诺数，已经打开的 k-ε 是不是过重，台阶分离是不是先于那次转捩，问题留在圆管流动开始散掉的那个雷诺数。