OpenFoamによる熱流体解析

はじめに

OpenFOAMによる熱流体解析を紹介します。ここでは、以下のような形状・条件を例題とします。バージョンはOpenFOAM-4.1とし、buoyantBoussinesqSimpleFoamを利用します。

  • 形状
  • 物性値
    密度
    1.0e+3kg/s
    粘性係数
    1.0e−3m·s
    比熱
    4.2e+3J/(kg·K)
    熱伝導率
    0.60W/(m·K)
  • 境界条件
    inlet1
    0.5m/s, 293.0K
    inlet2
    0.8m/s, 313.0K

作業ディレクトリの作成

最初に作業用ディレクトリを作成します。

  • 実行方法
    $ mkdir -p $FOAM_RUN/emsco-jp/mixingChamber/{0,constant,system}
    $ cd $FOAM_RUN/emsco-jp/mixingChamber
    

メッシュ作成

メッシュ作成は、メッシュ作成事例を参照してください。

重力加速度の設定

強制的に発生させた流れでは、重力の影響を無視できる場合があります。ここでは重力の影響を無視し、重力加速度をゼロに設定します。

  • constant/g
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       uniformDimensionedVectorField;
        location    "constant";
        object      g;
    }
    
    dimensions      [0 1 -2 0 0 0 0];       // [kg m s K kgmol A cd]
    
    value           (0 0 0);
    

乱流モデルの設定

標準k-εモデルに設定します。

  • constant/turbulenceProperties
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        location    "constant";
        object      turbulenceProperties;
    }
    
    simulationType RAS;
    
    RAS
    {
        RASModel        kEpsilon;           // 標準k-εモデル
        turbulence      true;
        printCoeffs     true;
    }
    

物性値の設定

OpenFOAMでは非圧縮性流れ解析の場合、粘性係数のかわりに動粘性係数(粘性係数を密度で割った値)を設定します。また、重力を無視する場合は、体積膨張係数と参照温度は結果に影響しないので、適当な値とします。プラントル数(動粘性係数を温度拡散係数で割った値、または、粘性係数に比熱を乗じて熱伝導率で割った値)と乱流プラントル数(0.7~0.9)を設定します。

  • constant/transportProperties
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        object      transportProperties;
    }
    
    transportModel Newtonian;
    
    nu          [0 2 -1 0 0 0 0]    1e-6;           // 動粘性係数(nu=mu/rho)
    beta        [0 0 0 -1 0 0 0]    0;              // 体積膨張係数
    TRef        [0 0 0 1 0 0 0]     293;            // 参照温度
    Pr          [0 0 0 0 0 0 0]     7;              // プラントル数(Pr=nu/alpha=mu*Cp/k)
    Prt         [0 0 0 0 0 0 0]     0.85;           // 乱流プラントル数
    

境界条件と初期条件の設定

流速、圧力、温度、乱流エネルギーおよび乱流散逸率の境界条件と初期条件を設定します。buoundaryFieldに境界条件を設定し、internalFieldに各セルの初期値(以下のファイルでは初期値を一様に与えています)を設定します。また、乱流エネルギーの入口境界条件は乱流強度で設定します。乱流散逸率の入口境界条件は混合長さで指定し、水力直径の0.07倍とします。

  • 0/U
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       volVectorField;
        location    "0";
        object      U;
    }
    
    dimensions      [0 1 -1 0 0 0 0];       // [kg m s K kgmol A cd]
    
    internalField   uniform (0 0 0);
    
    boundaryField
    {
        inlet1
        {
            type            fixedValue;
            value           uniform (0.5 0 0);
        }
        inlet2
        {
            type            fixedValue;
            value           uniform (0 -0.8 0);
        }
        outlet
        {
            type            zeroGradient;
        }
        wall
        {
            type            fixedValue;
            value           uniform (0 0 0);
        }
    }
    
  • 0/p_rgh
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       volScalarField;
        object      p_rgh;
    }
    
    // 非圧縮性ソルバの場合、圧力は密度で割った単位になる
    dimensions      [0 2 -2 0 0 0 0];       // [kg m s K kgmol A cd]
    
    internalField   uniform 0;
    
    boundaryField
    {
        inlet1
        {
            type            zeroGradient;
        }
        inlet2
        {
            type            zeroGradient;
        }
        outlet
        {
            type            fixedValue;
            value           uniform 0;
        }
        wall
        {
            type            zeroGradient;
        }
    }
    
  • 0/nut
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       volScalarField;
        location    "0";
        object      nut;
    }
    
    dimensions      [0 2 -1 0 0 0 0];       // [kg m s K kgmol A cd]
    
    internalField   uniform 0;
    
    boundaryField
    {
        inlet1
        {
            type            calculated;
            value           uniform 0;
        }
        inlet2
        {
            type            calculated;
            value           uniform 0;
        }
        outlet
        {
            type            calculated;
            value           uniform 0;
        }
        wall
        {
            type            nutkWallFunction;
            value           uniform 0;
        }
    }
    
  • 0/k
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       volScalarField;
        location    "0";
        object      k;
    }
    
    dimensions      [0 2 -2 0 0 0 0];       // [kg m s K kgmol A cd]
    
    internalField   uniform 1;
    
    boundaryField
    {
        inlet1
        {
            type            turbulentIntensityKineticEnergyInlet;
            intensity       0.05;
            value           uniform 1;
        }
        inlet2
        {
            type            turbulentIntensityKineticEnergyInlet;
            intensity       0.05;
            value           uniform 1;
        }
        outlet
        {
            type            zeroGradient;
        }
        wall
        {
            type            kqRWallFunction;
            value           uniform 1;
        }
    }
    
  • 0/epsilon
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       volScalarField;
        location    "0";
        object      epsilon;
    }
    
    dimensions      [0 2 -3 0 0 0 0];       // [kg m s K kgmol A cd]
    
    internalField   uniform 1;
    
    boundaryField
    {
        inlet1
        {
            type            turbulentMixingLengthDissipationRateInlet;
            mixingLength    0.007;          // 0.07 * 0.1(m)
            value           uniform 1;
        }
        inlet2
        {
            type            turbulentMixingLengthDissipationRateInlet;
            mixingLength    0.0035;         // 0.07 * 0.1(m)
            value           uniform 1;
        }
        outlet
        {
            type            zeroGradient;
        }
        wall
        {
            type            epsilonWallFunction;
            value           uniform 1;
        }
    }
    
  • 0/T
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       volScalarField;
        location    "0";
        object      T;
    }
    
    dimensions      [0 0 0 1 0 0 0];        // [kg m s K kgmol A cd]
    
    internalField   uniform 293;
    
    boundaryField
    {
        inlet1
        {
            type            fixedValue;
            value           uniform 293;
        }
        inlet2
        {
            type            fixedValue;
            value           uniform 313;
        }
        outlet
        {
            type            zeroGradient;
        }
        wall
        {
            type            zeroGradient;
        }
    }
    
  • 0/alphat
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       volScalarField;
        location    "0";
        object      alphat;
    }
    
    dimensions      [0 2 -1 0 0 0 0];       // [kg m s K kgmol A cd]
    
    internalField   uniform 0;
    
    boundaryField
    {
        inlet1
        {
            type            calculated;
            value           uniform 0;
        }
        inlet2
        {
            type            calculated;
            value           uniform 0;
        }
        outlet
        {
            type            calculated;
            value           uniform 0;
        }
        wall
        {
            type            alphatJayatillekeWallFunction;
            Prt             0.85;
            value           uniform 0;
        }
    }
    

離散化スキームの設定

移流項の離散化スキームには、計算安定性の良い一次精度の風上差分法を選択します。

  • system/fvSchemes
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        location    "system";
        object      fvSchemes;
    }
    
    ddtSchemes
    {
        default         steadyState;
    }
    
    gradSchemes
    {
        default         Gauss linear;
    }
    
    divSchemes
    {
        default         false;
        div(phi,U)      bounded Gauss upwind;
        div(phi,T)      bounded Gauss upwind;
        div(phi,k)      bounded Gauss upwind;
        div(phi,epsilon) bounded Gauss upwind;
        div((nuEff*dev2(T(grad(U))))) Gauss linear;
    }
    
    laplacianSchemes
    {
        default         Gauss linear corrected;
    }
    
    interpolationSchemes
    {
        default         linear;
    }
    
    snGradSchemes
    {
        default         corrected;
    }
    
    wallDist
    {
        method          meshWave;
    }
    

代数方程式ソルバと圧力・速度連成方法の設定

代数方程式ソルバと圧力・速度連成方法を設定します。buoyantBoussinesqSimpleFoamでは、圧力と速度の連成方法としてSIMPLE法が使われています。

  • system/fvSolution
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        location    "system";
        object      fvSolution;
    }
    
    solvers
    {
        p_rgh
        {
            solver          PCG;
            preconditioner  DIC;
            tolerance       1e-8;
            relTol          0;
        }
        U
        {
            solver          PBiCG;
            preconditioner  DILU;
            tolerance       1e-6;
            relTol          0;
        }
        "(k|epsilon)"
        {
            solver          PBiCG;
            preconditioner  DILU;
            tolerance       1e-6;
            relTol          0;
        }
        T
        {
            solver          PBiCG;
            preconditioner  DILU;
            tolerance       1e-6;
            relTol          0;
        }
    }
    
    SIMPLE
    {
        nNonOrthogonalCorrectors 0;
    
        //pRefCell        0;
        //pRefValue       0;
    
        //consistent      false;
    
        residualControl
        {
            //p_rgh           1e-4;
            //U               1e-4;
            //"(k|epsilon)"   1e-4
            //T               1e-4;
        }
    }
    
    relaxationFactors
    {
        fields
        {
            p_rgh           0.3;
        }
        equations
        {
            U               0.7;
            "(k|epsilon)"   0.8;
            T               1;
        }
    }
    

制御パラメータとモニタの設定

計算制御パラメータと残差モニタを設定します。buoyantBoussinesqSimpleFoamなどの定常ソルバでは、通例としてdeltaTを1にして、イテレーション数と時間が一致するようにします。計算終了時間のendTimeは500とします。

  • system/controlDict
    FoamFile
    {
        version 2.0;
        format      ascii;
        class       dictionary;
        location    "system";
        object      controlDict;
    }
    
    application     buoyantBoussinesqSimpleFoam;
    
    startFrom       latestTime;
    
    startTime       0;
    
    stopAt          endTime;
    
    endTime         500;
    
    deltaT          1;
    
    writeControl    timeStep;
    
    writeInterval   100;
    
    purgeWrite      0;
    
    writeFormat     ascii;
    
    writePrecision  6;
    
    writeCompression false;
    
    timeFormat      general;
    
    timePrecision   6;
    
    runTimeModifiable true;
    
    functions
    {
        residuals
        {
            type        residuals;
    
            libs
            (
                "libutilityFunctionObjects.so"
            );
    
            writeControl timeStep;
    
            writeInterval 1;
    
            fields
            (
                U
                k
                epsilon
                T
                p_rgh
            );
        }
        outletNut
        {
            type        surfaceRegion;
    
            libs
            (
                "libfieldFunctionObjects.so"
            );
    
            writeControl timeStep;
    
            writeInterval 1;
    
            log         true;
    
            writeFields false;
    
            regionType  patch;
    
            name        outlet;
    
            operation   areaAverage;
    
            fields
            (
                nut
            );
        }
        temperatureProbes
        {
            type        probes;
    
            libs
            (
                "libsampling.so"
            );
    
            writeControl timeStep;
    
            writeInterval 1;
    
            fields
            (
                T
            );
    
            probeLocations
            (
                (      .2           0.           0.)
                (      .5           0.           0.)
                (      .7          -.2           0.)
                (      .7          -.5           0.)
            );
        }
    }
    

領域分割数の設定

並列計算のために領域分割数を設定します。ここでは並列数は4とします。

  • system/decomposeParDict
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        location    "system";
        object      decomposeParDict;
    }
    
    numberOfSubdomains 4;
    
    method      simple;
    
    simpleCoeffs
    {
        n           (2 2 1);
        delta       1e-6;
    }
    

計算の実行

領域分割してから、計算を実行します。また計算実行中はfoamMonitorを利用して残差などを確認します。

  • 実行方法
    $ decomposePar
    $ mpirun -np 4 buoyantBoussinesqSimpleFoam -parallel > log.buoyantBoussinesqSimpleFoam &
    $ foamMonitor -logscale -refresh 1 postProcessing/residuals/0/residuals.dat &
    $ foamMonitor -refresh 1 postProcessing/outletNut/0/surfaceReion.dat &
    $ foamMonitor -refresh 1 postProcessing/temperatureProbes/0/T &
    $ tailf log.buoyantBoussinesqSimpleFoam
    

解析結果

計算が完了したらParaViewで結果を確認します。

  • 温度コンター