OpenFOAMを用いたサイクロン流れ解析のメモ

はじめに

本記事では、OpenFOAMを用いてサイクロン内の流れ解析を行った結果を紹介します。サイクロンは粒子を気流から分離するための装置であり、工業的にも広く利用されています。数値流体力学（CFD）による詳細な流れの解析は、設計改善や性能評価に大きな価値を持ちます。

モデル

stl作成ソフト：FREECAD
使用stlファイル：wall.slt, buffle.stl, inlet.stl, outlet_top.stl, outlet_bottom.stl
形状: サイクロン（集塵機）※形状は化学工学便覧にも掲載している標準サイズとしている

形状

設定概要

使用ソルバーと設定概要

ソルバー: simpleFoam
乱流モデル: k-εモデル（ファイル k, epsilon を使用）
格子生成: snappyHexMesh
流体物性: 圧縮性流体（ファイル rho, T を使用）

初期条件の確認

シミュレーションデータをもとに、代表的なフィールドである速度 U と圧力 p を調査しました。

フィールド	内容
速度 `U`	全セルで `(0 0 0)`（静止状態） inletの面から20m/sの速度で空気が流入し続ける。
圧力 `p`	全セルで `101325` Pa（一様な大気圧）

初期状態として、流れのない状態からサイクロンが駆動される様子を再現することを意図しています。

使用された乱流モデル

シミュレーション種別: RAS（Reynolds-Averaged Simulation、時間平均型乱流モデル）
使用モデル: kOmegaSST（SST k-ω モデル）

SST（Shear Stress Transport）k-ω モデルは、壁面近傍の精度と自由流中の乱流モデルの安定性を両立するために設計されており、サイクロンのような旋回流や分離流が発生する流れに対して非常に有効です。

物理モデルと数値スキームの設定

1. 乱流モデルの選択：SST k-ω モデル

この解析では、RAS型乱流モデル（時間平均）が選択されており、具体的には kOmegaSST モデルが使用されています。これは、次のような特徴を持ちます：

壁面近傍では k-ω モデルとしての精度
主流域では k-ε モデルとしての安定性
渦や分離流が強いサイクロン流に適した性能

2. 熱物性の設定：空気・理想気体モデル

熱力学モデルには perfectGas（理想気体）が採用され、流体は空気（分子量28.966 g/mol）を想定しています。その他の代表的な設定は以下の通りです：

項目	値	単位
動粘性係数 `ν`	1.5 × 10⁻⁵	m²/s
粘性係数 `μ`	1.831 × 10⁻⁵	Pa·s
比熱 `Cp`	1005	J/(kg·K)
プラントル数 `Pr`	0.705	—

このような設定により、流体の熱的性質や粘性の効果を考慮した高精度な流れ解析が可能になります。

各種設定のファイル抜粋

０/U

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    wall
    {
        type            noSlip;
    }

    buffle
    {
        type            noSlip;
    }

    outlet_bottom
    {
        type pressureInletOutletVelocity;
        value uniform (0 0 0);
    }

    outlet_top
    {
        type pressureInletOutletVelocity;
        value uniform (0 0 0);

    }

    inlet
    {
        type            fixedValue;
        value           uniform (20 0 0); 
    }
}

0/ｐ

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      ;

internalField   uniform 101325;

boundaryField
{
    wall
    {
        type            zeroGradient;
    }

    buffle
    {
        type            zeroGradient;
    }

    outlet_bottom
    {
        type            zeroGradient;
    }
    outlet_top
    {
        type            zeroGradient;
    }
    inlet
    {
        type            fixedValue;
        value           uniform 121325;
    }
}

system/controlDict

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     rhoPimpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         1;

deltaT          0.001;

writeControl    timeStep;

writeInterval   100;

purgeWrite      200;  //古い出力データの削除

writeFormat     ascii;

writePrecision  10; //数値の精度（小数点以下の桁数）

writeCompression off;  //出力ファイルの圧縮設定

timeFormat      general;

timePrecision   6;  //時間データの精度

runTimeModifiable true;

adjustTimeStep  yes;

maxCo          0.5;

maxDeltaT       0.1;

snappyHexMeshDictの留意事項（設定概要と注意）

featuresのレベル設定

    features
    (
        {
            file "wall.eMesh";
            level 1;
        }
        {
            file "inlet.eMesh";
            level 1;
        }
        {
            file "outlet_bottom.eMesh";
            level 2;
        }
        {
            file "outlet_top.eMesh";
            level 2;
        }
        {
            file "buffle.eMesh";
            level 1;
        }

    );

refinementSurfacesの設定

wallとbuffleのメッシュレベルの最大を境界（inlet、outlet_top、outlet_bottom）とそろえて5としている。

refinementSurfaces
    {
        wall
        {
            level (2 5);
            patchInfo
            {
                type wall;
            }
        }
        inlet
        {
            level (4 5);
            patchInfo
            {
                type patch;
            }
        }

        buffle
        {
            level (2 5);
            patchInfo
            {
                type wall;
            }

        }

        outlet_top
        {
            level (4 5);
            patchInfo
            {
                type patch;
            }
        }

        outlet_bottom
        {
            level (4 5);
            patchInfo
            {
                type patch;
            }
        }
    }

snapControlsの設定

multiRegionFeatureSnapの設定は「true」にしておく必要があります。初期のfalseのままにしておくと境界パッチが分離されて、発散してエラーの原因になってしまいました。（これに気が付くのにすごく時間がかかってしまった、、、）

snapControls
{
    nSmoothPatch 5;

    tolerance 2.0;

    nSolveIter 40;

    nRelaxIter 5;

        nFeatureSnapIter 10;

        implicitFeatureSnap false;

        explicitFeatureSnap true;

        multiRegionFeatureSnap true;
}

ちなみに各設定の詳細も以下に説明を載せておきます。

各パラメータの詳細解説

`nSmoothPatch 5`

吸着前にパッチの法線ベクトルを平均化するスムージング処理の回数
値を大きくすると吸着がなめらかになりやすいが、形状が丸まりやすくなる
通常 5〜10 くらいが妥当

`tolerance 2.0`

吸着の際、点が特徴点・エッジに引き寄せられる距離の係数
実際の距離は「この値 × ローカルなセルエッジ長さ」
値を大きくすると、遠くのポイントも吸着されやすくなるが、メッシュが歪む可能性もある
推奨範囲：1.0〜2.5

`nSolveIter 40`

メッシュ変形を解く反復回数
スナッピングでポイントが動く際のメッシュの滑らかさを保つための内部ソルバのステップ数
精度が欲しいときは大きめ（30〜100）、時間短縮したいときは小さめ

`nRelaxIter 5`

スナッピング処理の中でのリラクゼーション反復の最大回数
通常は 5〜10程度でOK

特徴エッジのスナッピング設定

この部分は、STL表面に加えて、エッジ（角）にも吸着するかどうかを制御します。

`nFeatureSnapIter 10`

特徴エッジにスナップする反復回数
これを省略・コメントアウトすると、エッジへの吸着が無効になります
10〜20あれば十分な場合が多いです

`implicitFeatureSnap false`

エッジを自動で検出（STLの角度に基づく）するかどうか
false にすると、自動検出せず、.eMesh で指定されたものだけを使う

`explicitFeatureSnap true`

castellatedMeshControls::features セクションで指定した .eMesh ファイルを使って、 明示的に定義した特徴エッジにスナップする。
通常は true にします（STLに基づいたエッジ指定が正確なため）

`multiRegionFeatureSnap false`

複数のSTLサーフェスが交わる点での吸着を有効にするか
一般的な単一領域のケースでは false でOK
多領域（multi-region）での接触に対して吸着したい場合は true　←とても重要

system/fvSchemes

チュートリアルから持ってきてほぼそのまま使用している。

div(phi,T)のみ Gauss upwind;に変更している。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}
d2dt2Schemes
{
    default         Euler;
}
gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(phi,U)      bounded Gauss upwind;
    div(U)          Gauss linearUpwind grad(U);

    div(phi,p)     Gauss linearUpwind grad(p);
    div(phi,rho)    Gauss linearUpwind grad(rho);
    div(phi,K)      Gauss linearUpwind grad(K);
    div(phi,T)      Gauss upwind;
    div(phi,h)      Gauss linearUpwind grad(h);
    div(phi,k)      Gauss linearUpwind grad(k);
    div(phi,epsilon) Gauss linearUpwind grad(epsilon);
    div(phi,R)      Gauss linearUpwind grad(R);
    div(phi,omega)  Gauss linearUpwind grad(omega);
    div((rho*R))    Gauss linear;
    div(R)          Gauss linear;

    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}

wallDist
{
    method          meshWave;
}

system/fvSolution

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        smoother        DIC;
        tolerance       1e-8;
        relTol         0.05;
        nPreSweeps     0;
        nPostSweeps    2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    rho
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.01;
    }

    rhoFinal
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-10;  // 最終的な収束条件を厳しくする
        relTol          0;
    }

    "(U|h|R|k|epsilon|omega)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }
    "(U|h|R|k|epsilon|omega)Final"
    {
        $U;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor yes;
    nOuterCorrectors 2;
    nCorrectors     3;
    finalOnLastPimpleIterOnly   true;
    nNonOrthogonalCorrectors 3;
    pRefCell 42001;
    pRefValue 101325;
    //pRefValue 0;
    residualControl
    {
        p               1e-6;
        U               1e-6;
        h               1e-6;
        "(k|omega|epsilon)" 1e-6;
    }

    //fvConstraints に pressure limits を追加（使用していないので無視してください）
    //pMax            1.2e5;
    //pMin            0.8e5;
}

constant/thermophysicalProperties

このファイルはとても重要ですので全部載せておきます。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermophysicalProperties;
}

viscosityModel Newtonian;
nu              1.5e-05;  // [m?/s] 運動粘性係数（動粘性）

// 熱物性モデル
thermoType
{
    type             hePsiThermo;
    mixture         pureMixture;
    transport       const;
    thermo         hConst;
    //thermo         eConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
    //energy          sensibleInternalEnergy;
}

// 分子量
mixture
{
    specie
    {
        molWeight   28.966;
    }
    // 熱力学パラメータ
    thermodynamics
    {
       Cp          1005;
       Hf          0;
    }

  transport
  {
      mu  1.831e-05;
      Pr   0.705;
  }
}