はじめに

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で結果を確認します。

• 温度コンター