fluid-dynamics-workflow

流体力学解析ワークフロー

数値流体力学（CFD）および解析流体力学の問題に対する、包括的で体系的なアプローチです。

概要

このワークフローでは、問題の初期設定から検証まで、流体力学解析の全工程を順を追って説明します。支配物理の厳密な取り扱い、適切なモデル選定、および結果の適切な検証を実現します。

解析ワークフロー全体

1. 問題の定義

1.1 形状と領域

• 物理領域：

  • 計算領域の範囲を定義する
  • 主要な幾何学的特徴（境界層、剥離領域、再循環）を特定する
  • 計算コストを削減するため、対称面を決定する
  • 座標系（直交座標、円筒座標、球座標）を設定する
  • 2D と 3D のどちらが必要かを評価する

• 領域サイズに関する考慮事項：

  • 流入口：代表長さの 5～10 倍上流
  • 流出口：代表長さの 20～30 倍下流
  • 横方向境界：物体から代表長さの 5～10 倍
  • 人為的な閉塞の影響がないことを確認する

1.2 流動領域の特性評価

レイノルズ数解析：

Re = ρVL/μ = VL/ν

Where:
- ρ: density
- V: characteristic velocity
- L: characteristic length
- μ: dynamic viscosity
- ν: kinematic viscosity

流動領域の分類：

• Re < 2300：層流（管内流）
• 2300 < Re < 4000：遷移流
• Re > 4000：乱流（管内流）
• 外部流：
  • Re < 5×10⁵：層流境界層
  • Re > 5×10⁵：乱流境界層

マッハ数解析：

Ma = V/c

Where:
- V: flow velocity
- c: speed of sound = √(γRT)

圧縮性の分類：

• Ma < 0.3：非圧縮性（ρ = 一定）
• 0.3 < Ma < 0.8：亜音速圧縮性
• 0.8 < Ma < 1.2：遷音速
• 1.2 < Ma < 5.0：超音速
• Ma > 5.0：極超音速

フルード数（自由表面流）：

Fr = V/√(gL)

• Fr < 1：常流
• Fr = 1：限界流
• Fr > 1：射流

1.3 その他の無次元数

ストローハル数（非定常流）：

St = fL/V

• 振動流動機構を特徴付けます
• 渦放出：円柱では St ≈ 0.2

プラントル数（熱伝達）：

Pr = ν/α = μcp/k

グラスホフ数（自然対流）：

Gr = gβΔTL³/ν²

2. 支配方程式

2.1 連続の式

非圧縮性：

∇·V = 0

圧縮性：

∂ρ/∂t + ∇·(ρV) = 0

2.2 運動量方程式（ナビエ・ストークス方程式）

非圧縮性ニュートン流体：

ρ(∂V/∂t + V·∇V) = -∇p + μ∇²V + ρg

成分形式（直交座標）：

ρ(∂u/∂t + u∂u/∂x + v∂u/∂y + w∂u/∂z) = -∂p/∂x + μ(∂²u/∂x² + ∂²u/∂y² + ∂²u/∂z²) + ρgₓ
ρ(∂v/∂t + u∂v/∂x + v∂v/∂y + w∂v/∂z) = -∂p/∂y + μ(∂²v/∂x² + ∂²v/∂y² + ∂²v/∂z²) + ρgᵧ
ρ(∂w/∂t + u∂w/∂x + v∂w/∂y + w∂w/∂z) = -∂p/∂z + μ(∂²w/∂x² + ∂²w/∂y² + ∂²w/∂z²) + ρgᵣ

圧縮性形式：

ρ(∂V/∂t + V·∇V) = -∇p + ∇·τ + ρg

ここで、応力テンソルは次のとおりです：

τᵢⱼ = μ(∂uᵢ/∂xⱼ + ∂uⱼ/∂xᵢ) + λ(∇·V)δᵢⱼ

2.3 エネルギー方程式

圧縮性流れの場合：

ρcp(∂T/∂t + V·∇T) = ∇·(k∇T) + Φ + Q

ここで：

• Φ：粘性散逸
• Q：熱源項

全エネルギー形式：

∂(ρE)/∂t + ∇·(ρEV) = -∇·(pV) + ∇·(k∇T) + ∇·(τ·V)

2.4 乱流モデル

レイノルズ平均ナビエ・ストークス（RANS）： 変数を分解します：u = ū + u'

レイノルズ応力テンソル：

τᵢⱼᴿᴬᴺˢ = -ρu'ᵢu'ⱼ

クロージャモデルの選択フロー：

START
  │
  ├─ Simple geometry, equilibrium turbulence?
  │  └─ YES → k-ε Standard
  │
  ├─ Wall-bounded flow, adverse pressure gradient?
  │  └─ YES → k-ω SST (Menter)
  │
  ├─ Separated flow, vortex shedding?
  │  └─ YES → k-ω SST or Transition SST
  │
  ├─ Swirling flow, rotation?
  │  └─ YES → Reynolds Stress Model (RSM)
  │
  ├─ High accuracy needed, large-scale unsteadiness?
  │  └─ YES → Large Eddy Simulation (LES)
  │
  └─ Research, turbulence structure critical?
     └─ YES → Direct Numerical Simulation (DNS)

モデルの比較：

モデル	利点	制約	最適な用途
k-ε Standard	頑健で、十分に検証されている	壁面近傍や剥離流で精度が低い	産業流れ、単純な形状
k-ε Realizable	剥離流に対する精度が高い	壁関数が必要	複雑な形状、剥離
k-ω Standard	壁面近傍で良好	自由流に敏感	壁面拘束流
k-ω SST	全般的に優れている	計算コストが高い	空気力学、逆圧力勾配
Spalart-Allmaras	1 方程式で効率的	航空宇宙分野向けに調整されている	翼型、外部空気力学
RSM	異方性を捉える	5～7 方程式で複雑	旋回流、回転流
LES	大規模渦を解像する	計算コストが非常に高い	非定常流、音響
DNS	モデル化誤差がない	計算コストが極めて高い	基礎研究、Re < 10⁴

3. 境界条件

3.1 流入境界条件

速度流入：

• 指定：u、v、w（速度成分）
• 乱流の場合は、次の値も指定します：
  • 乱流強度：I = u'/U（通常は1～10％）
  • 乱流長さスケール：l = 0.07L
  • または直接指定：k、ε（またはω）

乱流量の計算：

k = 3/2(UI)²
ε = C_μ^(3/4) k^(3/2) / l
ω = k^(1/2) / (C_μ^(1/4) l)

ここで、C_μ = 0.09です。

圧力流入：

• 指定：全圧p₀
• 流れの方向
• 温度（圧縮性流れの場合）

質量流量流入：

• 指定：ṁ = ρVA
• 流れの方向

3.2 流出境界条件

圧力流出：

• 静圧を指定します
• その他すべての量にゼロ勾配を設定します：∂φ/∂n = 0
• 最も一般的で安定した方法です

流出：

• 十分に発達した流れを仮定します
• ∂u/∂x = 0（すべての変数）
• 流出口が十分下流にある場合に使用します

流出口までの距離要件：

• 下流方向に最低10～15油圧直径
• はく離流れの場合：20～30直径

3.3 壁面境界条件

すべりなし壁面（標準）：

• 壁面でV = 0
• 温度（等温）または熱流束を指定します

すべり壁面：

• 接線速度 ≠ 0
• 法線速度 = 0
• V·n = 0、∂V_tangential/∂n = 0

壁関数アプローチ： 乱流の場合、次の条件を満たすときに壁関数を使用します：

• 30 < y⁺ < 300

y⁺ = ρu_τy/μ
u_τ = √(τ_w/ρ)  [friction velocity]

壁面近傍の解像度要件：

アプローチ	第1セルのy⁺	プリズム層の合計数	成長率
壁関数	30-300	3-5	1.2-1.3
低レイノルズ数モデル	y⁺ < 1	10-20	1.1-1.2
拡張壁面処理	y⁺ < 5	8-15	1.15-1.25

移動壁面：

• 壁面速度を指定します
• 用途：回転機械、コンベヤー

3.4 対称境界と周期境界

対称境界：

• 法線速度 = 0
• その他の量にゼロ勾配を設定します
• ∂φ/∂n = 0

周期境界：

• 並進：φ(x₁) = φ(x₂)
• 回転：回転ジオメトリの場合
• 計算領域の縮小に使用します

3.5 自由表面（混相流）

VOF（流体体積法）：

• 体積分率で界面を追跡します
• α = 1（相1）、α = 0（相2）

レベルセット法：

• 符号付き距離関数を使用します
• 鮮明な界面を追跡します

4. 離散化

4.1 メッシュ要件

メッシュの種類：

メッシュの種類	利点	欠点	最適な用途
構造格子（六面体）	効率的、高精度	複雑なジオメトリには不向き	単純な領域、流路
非構造格子（四面体）	ジオメトリに柔軟に対応	より多くのセルが必要	複雑なCADジオメトリ
ハイブリッド（プリズム＋四面体）	両者の長所を兼備	慎重な設定が必要	ほとんどのCFDアプリケーション
多面体	セル数が少なく、安定	専用のメッシャーが必要	産業用CFD

メッシュ品質指標：

• アスペクト比 < 100（< 20が望ましい）
• スキューネス < 0.85（< 0.5が望ましい）
• 直交性 > 0.15（> 0.7が望ましい）

格子独立性の検討：

Refinement ratio: r = 1.5 to 2.0
Test at least 3 mesh levels:
- Coarse: N cells
- Medium: r²N cells
- Fine: r⁴N cells

Check convergence of key outputs:
- Forces (drag, lift)
- Pressures at critical points
- Flow rates

リチャードソン外挿：

f_exact ≈ f_fine + (f_fine - f_coarse)/(r^p - 1)

ここで、pは精度の次数です。

4.2 乱流のY⁺値

重要なy⁺範囲：

Viscous sublayer:     y⁺ < 5
Buffer layer:         5 < y⁺ < 30
Log-layer:            30 < y⁺ < 300
Outer layer:          y⁺ > 300

第1セル高さの計算：

y = y⁺μ/(ρu_τ)

For flat plate:
u_τ ≈ U∞√(C_f/2)
C_f = 0.058Re_L^(-0.2)  [turbulent]

目標のy⁺に対するyの推定：

Example: Airfoil at Re = 1×10⁶, V = 50 m/s
1. Estimate C_f ≈ 0.0037
2. u_τ = 50√(0.0037/2) = 2.15 m/s
3. For y⁺ = 1: y = 1×1.5e-5/(1.225×2.15) ≈ 5.7e-6 m

First cell: 5.7 μm

4.3 時間刻みの選定（非定常）

CFL条件（クーラン・フリードリヒス・ルーウィ条件）：

CFL = VΔt/Δx ≤ CFL_max

For explicit schemes: CFL < 1
For implicit schemes: CFL < 10-20 (stability)
                      CFL < 1-5 (accuracy)

時間刻みの計算：

Δt = CFL × Δx_min / V_max

非定常現象の時間解像度：

• 渦放出を捉える：Δt < T/100（T = 放出周期）
• 過渡解析：Δt < τ/50（τ = 特性時間）

デュアルタイムステップ法：

• 物理時間刻み：精度のため
• 内部反復：収束のため

5. 解法戦略

5.1 ソルバーの選択

圧力・速度カップリング：

アルゴリズム	種類	利点	使用例
SIMPLE	分離型	メモリ効率が高い	定常状態
SIMPLEC	分離型	収束が速い	定常／非定常
PISO	分離型	非定常計算の精度が高い	非定常性が非常に強い流れ
Coupled	完全連成型	堅牢で高速	圧縮性流れ、複雑な流れ

運動量の離散化：

• 1 次精度風上差分：安定していますが、数値拡散が大きくなります
• 2 次精度風上差分：精度が良好です
• QUICK：精度が高く、値が有界になります
• 中心差分：数値拡散が小さいものの、安定性が必要です

圧力補間：

• Standard：ほとんどのケース
• PRESTO：旋回流、高速回転流
• Body-force-weighted：浮力駆動流

5.2 収束基準

残差のモニタリング：

Continuity:   < 1e-4
Velocity:     < 1e-4
Turbulence:   < 1e-4
Energy:       < 1e-6

追加の収束チェック：

• 力係数（C_L、C_D）をモニタリングします
• 質量流量の収支（流入量＝流出量）を追跡します
• 積分量を確認します
• 流れ場を可視化し、妥当性を確認します

収束プロットには次の特徴が必要です：

• 残差が滑らかに減少している
• 定常値に漸近している
• すべての残差が同時に減少している

5.3 アンダーリラクゼーション係数

目的：解を安定させ、発散を防ぎます

一般的な値（SIMPLE）：

Pressure:         0.3
Momentum:         0.7
Turbulence (k):   0.8
Turbulence (ε,ω): 0.8
Energy:           0.9

調整方法：

• 不安定な問題では値を下げます
• 解が進展するにつれて、徐々に値を上げます
• 連成ソルバーでは、より高い係数を使用できます

6. ポスト処理と検証

6.1 流れの可視化

必須の可視化項目：

1. 速度コンター／ベクトル：

   • 再循環領域を特定します
   • 境界層の発達を確認します
   • 流れの方向を検証します

2. 圧力分布：

   • 表面圧力係数：C_p = (p - p_∞)/(½ρV²)
   • 圧力勾配を特定します
   • 非物理的な振動がないか確認します

3. 流線／流跡線：

   • 流れのパターンを可視化します
   • はく離／再付着を特定します
   • 出口で逆流がないか確認します

4. 渦度：

   • ω = ∇ × V
   • 渦構造を特定します
   • 乱流を可視化します

5. 壁面せん断応力：

   • τ_w = μ(∂u/∂y)_wall
   • はく離：τ_w = 0
   • 摩擦係数：C_f = τ_w/(½ρV²)

6. 乱流量：

   • 乱流運動エネルギー（k）
   • 乱流粘度比（μ_t/μ）
   • y⁺分布

6.2 力の計算

圧力力と粘性力：

F_pressure = ∫∫ p·n dA
F_viscous = ∫∫ τ·n dA
F_total = F_pressure + F_viscous

無次元係数：

C_D = D/(½ρV²A)    [Drag coefficient]
C_L = L/(½ρV²A)    [Lift coefficient]
C_M = M/(½ρV²AL)   [Moment coefficient]

圧力係数：

C_p = (p - p_∞)/(½ρV²)

摩擦係数：

C_f = τ_w/(½ρV²)

6.3 検証と妥当性確認

検証（方程式を正しく解いているか）：

1. 格子収束指数（GCI）：

GCI = F_s |ε|/(r^p - 1)

Where:
- F_s = safety factor (1.25 for 3+ grids)
- ε = (f_fine - f_coarse)/f_fine
- r = refinement ratio
- p = order of accuracy

2. 時間収束：

• Δt/2 で再計算します
• 時間刻み幅への依存性がないことを確認します

3. 反復収束：

• 残差が十分に減少していることを確認します
• 解の変数をモニタリングします

妥当性確認（適切な方程式を解いているか）：

1. 解析解：

   • ポアズイユ流れ（チャネル）
   • クエット流れ
   • ブラジウス境界層
   • 円柱周りのポテンシャル流れ

2. 実験データ：

   • 風洞計測
   • PIV（粒子画像流速測定法）
   • 熱線流速計測

3. ベンチマークケース：

   • 蓋駆動キャビティ流れ
   • 円柱周りの流れ
   • 後向きステップ流れ
   • Ahmed ボディ

妥当性確認の指標：

Mean Absolute Error:
MAE = (1/N)Σ|φ_CFD - φ_exp|

Root Mean Square Error:
RMSE = √[(1/N)Σ(φ_CFD - φ_exp)²]

Coefficient of Determination:
R² = 1 - Σ(φ_exp - φ_CFD)²/Σ(φ_exp - φ̄_exp)²

デシジョンツリー

流動様式のデシジョンツリー

START: Define Problem
  │
  ├─ Calculate Re = ρVL/μ
  │
  ├─ Re < 1 (Stokes flow)
  │  └─ Use: Creeping flow equations
  │
  ├─ 1 < Re < 100 (Low Re)
  │  └─ Use: Laminar NS, full domain
  │
  ├─ Re < 2300 (Internal)
  │  └─ Use: Laminar NS
  │
  ├─ 2300 < Re < 4000 (Transitional)
  │  └─ Use: Transition models (γ-Re_θ)
  │
  ├─ Re > 4000 (Turbulent)
  │  ├─ Calculate Ma = V/c
  │  │
  │  ├─ Ma < 0.3
  │  │  └─ Incompressible RANS
  │  │
  │  ├─ 0.3 < Ma < 1
  │  │  └─ Compressible RANS
  │  │
  │  └─ Ma > 1
  │     └─ Compressible + shock capturing
  │
  └─ Unsteady features?
     ├─ YES → Calculate St = fL/V
     │  ├─ St > 0.1 → Use URANS or LES
     │  └─ St < 0.1 → RANS may suffice
     └─ NO → RANS sufficient

乱流モデル選択ツリー

START: Turbulent Flow Confirmed
  │
  ├─ Computational Budget?
  │
  ├─ VERY HIGH (Research)
  │  └─ Use DNS (if Re allows)
  │
  ├─ HIGH
  │  ├─ Massively separated?
  │  │  └─ YES → LES (WALE, Smagorinsky)
  │  │
  │  └─ NO → URANS with fine mesh
  │
  └─ MODERATE (Engineering)
     │
     ├─ Flow Type?
     │
     ├─ EXTERNAL (Aero/Hydro)
     │  ├─ Pressure gradient?
     │  │  ├─ ADVERSE → k-ω SST
     │  │  └─ FAVORABLE → Spalart-Allmaras
     │  │
     │  └─ Separation likely? → k-ω SST
     │
     ├─ INTERNAL (Pipes, Channels)
     │  ├─ Simple geometry → k-ε Standard
     │  ├─ Complex, bends → k-ε Realizable
     │  └─ Heat transfer → k-ω SST
     │
     ├─ ROTATING/SWIRLING
     │  └─ RSM (Reynolds Stress Model)
     │
     └─ TRANSITIONAL
        └─ γ-Re_θ Transition SST

メッシュ戦略のデシジョンツリー

START: Domain Defined
  │
  ├─ Geometry Complexity?
  │
  ├─ SIMPLE (Rectangular, cylindrical)
  │  └─ Structured hexahedral
  │     └─ O-grid for cylinders
  │
  ├─ MODERATE
  │  └─ Hybrid mesh
  │     ├─ Prism layers at walls
  │     └─ Hex/tet in core
  │
  └─ COMPLEX (CAD geometry)
     └─ Unstructured
        ├─ Prism layers at walls (always)
        └─ Tet or poly in core

  │
  ├─ Boundary Layer Resolution?
  │
  ├─ WALL FUNCTIONS (y⁺ = 30-300)
  │  ├─ First cell: y⁺ ≈ 50
  │  ├─ Prism layers: 3-5
  │  └─ Growth rate: 1.2-1.3
  │
  └─ NEAR-WALL RESOLUTION (y⁺ < 1)
     ├─ First cell: y⁺ < 1
     ├─ Prism layers: 10-20
     ├─ Growth rate: 1.1-1.2
     └─ At least 10 cells in boundary layer

  │
  └─ Refinement Zones?
     ├─ Wake regions: 2-3× density
     ├─ Shear layers: graded refinement
     ├─ Recirculation: fine mesh
     └─ Free stream: coarse OK

よくある落とし穴と解決策

問題：解が収束しない

診断：

• 残差が振動している
• 力係数が変動している
• 解が発散している

解決策：

1. アンダーリラクゼーション係数を下げます
2. メッシュ品質を改善します（アスペクト比、歪み）
3. 最初は 1 次精度スキームを使用します
4. 境界条件の整合性を確認します
5. 時間刻みを小さくします（非定常の場合）
6. 連成ソルバーに切り替えます

問題：非物理的な結果

診断：

• 負の圧力（絶対圧）
• 入口での逆流
• 質量不均衡が 1% 超
• 非現実的な剥離

解決策：

1. 境界条件を確認します
2. 物性値を確認します
3. 適切な基準値が設定されていることを確認します
4. 座標系を確認します
5. メッシュにエラーがないか確認します
6. 乱流モデルの選択を見直します

問題：高い y⁺ 値

診断：

• 低レイノルズ数モデルを使用しているときに y⁺ > 5
• 壁関数を使用しているときに y⁺ < 30

解決策：

1. 最初のセルの高さを再計算します
2. インフレーションパラメーターを調整します
3. 拡張壁面処理を使用します
4. 壁関数と壁面解像アプローチを切り替えます

ベストプラクティスのまとめ

1. 常に問題の物理から始めます：

   • 無次元数を計算します
   • 支配的な現象を特定します
   • 適切な簡略化を選択します

2. 体系的にメッシュを構築します：

   • 粗いメッシュから始め、反復的に細分化します
   • 境界層の解像を優先します
   • 格子独立性試験を実施します

3. シンプルに始め、複雑さを加えていきます：

   • 定常 RANS から始めます
   • 安定性のために一次風上スキームを使用します
   • より高次のスキームへと進みます
   • 必要に応じて非定常性を加えます

4. あらゆる項目を監視します：

   • 残差
   • 力係数
   • 質量流量の収支
   • y⁺ の分布

5. 厳密に検証します：

   • 実験や理論と比較します
   • 物理的に実現可能か確認します
   • 不確かさの定量化を実施します
   • 仮定と制約を文書化します

6. 不確かさを報告します：

   • 格子収束指数
   • 検証データとの比較
   • 境界条件に対する感度
   • モデル形式の不確かさ

参考文献とリソース

基礎的な文献：

• Anderson 著「Computational Fluid Dynamics」
• Versteeg と Malalasekera 著「An Introduction to Computational Fluid Dynamics」
• Pope 著「Turbulent Flows」

検証データベース：

• NASA Turbulence Modeling Resource
• ERCOFTAC Database
• NPARC Alliance Validation Archive

ベストプラクティスのガイドライン：

• CFD における検証と妥当性確認のための AIAA ガイド
• ASME V&V 規格
• ERCOFTAC ベストプラクティスガイドライン