Skip to content
cfd-lab:~/ja/posts/2026-03-26-openfoam-back…● online
$_
Younglin's CFD Blog
毎日の記録 · 計算流体力学

この記事は著者が韓国語で書いた原文の翻訳です。 原文(한국어)を読む

NOTE #011DAY THU 유체역학DATE 2026.03.26UPDATED 2026.10.19READ 4 min read#OpenFOAM#Tutorial#simpleFoam#Backward-Facing-Step#Turbulence

【OpenFOAMチュートリアル】後方ステップ流の解析 — simpleFoamで再循環領域を捉える

下流1.5 mは格子のHで15Hです。注釈のH=0.01 mで読むと、同じ区間が150Hになります。

後方ステップの直後で、流れ方向の速度は負になります。その区間が終わる x が再付着点です。simpleFoam がその点を何 H に置くかは、格子の高さと nu の注釈の高さが同じときにだけ読めます。

Re のスライダーを 50 から 1500 へ動かすと、段の後ろの逆流区間が伸びます。粒子数のスライダー（50–500）は、その区間を埋める軌跡の本数だけを変えます。

스텝 직후 재순환 영역(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×(0.05)2=0.00375k = 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 2方程式モデル
    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);  // 2次精度
    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 での順は次のとおりです。

  1. Coloring で U を選び、Magnitude で速度の大きさを塗ります。
  2. Glyph でベクトルを撒き、逆流区間を見ます。
  3. 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: 残差（residual）が減少せずに振動する#

nNonOrthogonalCorrectors が 0 のまま非直交が大きいと、残差が揺れます。補正回数を上げます。

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

予想される結果のまとめ

項目値
レイノルズ数 (Re=UH/νRe = UH/\nu)1,000
再付着長さ (xr/Hx_r / H)6–8
最大逆流速度≈−0.1\approx -0.1 m/s
収束までの反復回数300–600回

Armaly et al. 1983 の Re=800Re = 800 は約 6H で、表の 6–8 の隣にあります。

xr≈6–8⋅H(Re∼103)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×0.01/10−5=1000Re = 1 \times 0.01 / 10^{-5} = 1000 と書くので、H を 0.01 m と置いています。表の 1000 はその注釈と一致します。格子の H を同じ式に入れると 1×0.1/10−5=1041 \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)×0.1=0.01(2/20) \times 0.1 = 0.01 m で、y が 0.005 未満のセルだけを拾うと、そのセルの半分です。

k は注釈の 1.5×(0.05)2=0.003751.5 \times (0.05)^2 = 0.00375 と一致します。epsilon の式に Cμ=0.09C_\mu = 0.09 を入れると Cμ3/4k3/2≈3.77×10−5C_\mu^{3/4} k^{3/2} \approx 3.77 \times 10^{-5} で、書かれた 0.0009 で割った長さは 0.042 m です。L≈0.1HL \approx 0.1H は 0.01 m か 0.001 m なので、そのときの epsilon は 0.00377 か 0.0377 でなければなりません。ファイルの 0.0009 はどちらでもありません。逆流 −0.1-0.1 m/s は入口 1 m/s の 10% で、強度 5% の 0.05 m/s の 2 倍です。

6Hを十倍に伸ばしてしまう段の高さ#

格子の H を残すと、表の Re 1000 は偽りです。nu を 10−410^{-4} に上げて、その 1000 が戻ります。注釈の H を残すなら scale を 0.01 に下げることになり、6H は 0.06 m まで縮んで、下流 1.5 m の手前に付きます。

どちらの高さを、このケースの H として書くでしょうか。入口の 1000 を管内遷移の Re として読むと、すでに入れた k-ε が過ぎているのか、段の剥離がその遷移より先なのかは、円管の流れが崩れる瞬間に問いが残っています。

関連記事

圧縮性多相流のための Rhie-Chow 補間法まとめ界面では水-空気、密度比 1000 で算術平均密度は 500.5、調和平均は 2.0。フィルタは駆動圧力だけにかける。
Mar 10, 2026·11 min read
歪んだ格子で勾配がずれるとき — Green–Gaussと重み付き最小二乗による再構成扁平セルでGreen–Gaussの面補間が勾配をほぼ直角にずらし、距離重みn=1が曲面誤差を0.0293から0.0222へ削る。
Jul 17, 2026·5 min read
圧力が市松模様に振動するとき — SIMPLE と Rhie-Chow の圧力-速度連成48点の周期格子では、naive の市松模様は 0.80 のまま、Rhie-Chow は 3.1e-16 まで消えます。
Jul 22, 2026·5 min read

役に立ったらシェアしてください。

GitHub でソースを見る ↗誤字・改善提案 ↗