CFDEM-PFM 境界埋め込み法ソルバー

今回はCFDEM-PFMの境界埋め込み法ソルバーであるcfdemSolverIBのチュートリアルを実施してみます。
CFDEM-PFMの紹介やインストール方法は以下の記事で紹介しています。

cfdemSolverIB

CFDEM-PFMのカップリング方法は,cfdemSolverPisoとcfdemSolverIBで２通りの方法が実装されており、DEM粒子とCFDのメッシュサイズにより使い分けます。

cfdemSolverPiso
DEM粒子のサイズがCFDのメッシュサイズよりも小さい場合に使用
cfdemSolverIB
DEM粒子のサイズがCFDのメッシュサイズよりも大きい場合に使用

また、cfdemSolverIBのカップリング手順は以下の通りです。
Step1
CFDソルバーを用いて、流体相の速度と圧力場を計算します。このとき、全領域を流体として扱い、DEM粒子の存在を無視します。
Step2
DEM計算から得られた粒子の運動情報を取り込み、流体相の速度場を修正します。この時、速度場への影響度は、CFDメッシュに対するDEM粒子の占める割合="void fraction"に応じて計算されます。
Step3
前のステップでの速度修正により、連続の式が満たされなくなるので、速度場をさらに修正してこれを考慮します。
cfdemSolverPisoの場合は、"void Fraction"に応じた運動量のソース項を追加することにより、速度場を計算するため、より単純なカップリング手順になります。ただし、Step3のような修正手順がなく、Void Fraction＝1.0となるような部分が存在すると計算が収束せず発散してしまうという特徴があります。
以下に各手法の模式図と"Void Fraction"の取扱いについての図を示します。

チュートリアル

今回はCFDEM-PFMに付属の以下のチュートリアルを実行してみます。

CFDEMcouplingPFM/tutorials/cfdemIB/falling_sphere_two_way_coupling

フォルダー構成は以下のようになっています。

falling_sphere_two_way_coupling/
        ┣  CFD/                     #Openfoam設定フォルダ
        ┃     ┝ 0/                  #初期値
        ┃     ┝　┝ p                #圧力           
        ┃     ┝　┝ phiIB            #補正計算のためのスカラーphiIB
        ┃     ┝　┝ rho              #密度
        ┃     ┝　┝ U                #流体速度
        ┃     ┝　┝ Us               #DEM粒子速度
        ┃     ┝　└ voidfraction     #ボイド率
        ┃     ┝        
        ┃     ┝ constant/
        ┃     ┝　      ┝ polymesh/               
        ┃     ┝　      ┝ couplingProperties    #カップリングパラメータの設定
        ┃     ┝　      ┝ g
        ┃     ┝　      ┝ liggghtsCommands
        ┃     ┝　      ┝ transportProperties
        ┃     ┝　      └ turblenceProperties
        ┃     ┝　
        ┃     ┝ octave/
        ┃     ┝　  └ plot_data.m
        ┃     ┝　  
        ┃     └ system/ 
        ┃     　    ┝ blockMeshDict 
        ┃     　    ┝ controlDict
        ┃     　    ┝ decomposeParDict 
        ┃     　    ┝ fvOptions 
        ┃     　    ┝ fvSchemes 
        ┃     　    └ fvSolution 
        ┃
        ┣  DEM/                     #Liggght用設定フォルダ
        ┃     ┝ post/               
        ┃     ┝ in.liggghts_init    
        ┃     └ in.liggghts_run                
        ┃ 
        ┣  Allrun.sh                #カップリング実行用スクリプトファイル
        ┣  parCFDDEMrun.sh          #CFDソルバーの実行ファイル
        ┗  parDEMrun.sh             #DEMソルバーの実行ファイル

cfdemSolverIBではphiIBというスカラー場を解きますが、これはカップリング手順Step3での補正計算で使用されるもので、連続の式の不釣り合い量に応じた分布を持ちます。
チュートリアルでは、出入口に”0”を指定し、壁面に"zeroGradient"を設定しています。

0/phiIB.c

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.6                                   |
|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      phiIB;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0;

boundaryField
{
    inlet      
    {
        //type            zeroGradient;
        type            fixedValue;
        value           uniform 0; 
    }

    wall      
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
        //type            zeroGradient;
    }

    frontAndBack    
    {
        type            empty;
    }
}

// ************************************************************************* //

また、カップリングパラメータの設定は以下のようになっています。
voidFractionModelとlocateModelでcfdemSolverIBの設定を行っています。

voidFractionModel ⇒　IB
locateModel ⇒　engineIB

constant/couplingProperties.c

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  6
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      couplingProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// sub-models & settings

modelType none;

checkPeriodicCells;

verbose;

couplingInterval 10;

depth 0;

voidFractionModel IB;    

locateModel engineIB;    

meshMotionModel noMeshMotion;

regionModel allRegion;

dataExchangeModel twoWayMPI;

IOModel basicIO;

probeModel off;

averagingModel dilute;

clockModel off;

smoothingModel off;

forceModels
(
    ShirgaonkarIB
    ArchimedesIB
);

momCoupleModels
(
);

turbulenceModelType "turbulenceProperties";

// sub-model properties

ShirgaonkarIBProps
{
    velFieldName "U";
    pressureFieldName "p";
    densityFieldName "rho";
    verbose;
    useTorque;
}

ArchimedesIBProps
{
    densityFieldName "rho";
    gravityFieldName "g";
    voidfractionFieldName "voidfractionNext";
}

twoWayMPIProps
{
    liggghtsPath "../DEM/in.liggghts_run";
}

IBProps
{
    maxCellsPerParticle 20000;
    alphaMin 0.30;
    scaleUpVol 1.0;
}

engineIBProps
{
    engineProps
    {
        treeSearch false;
    }
    zSplit 8;
    xySplit 16;
}

// ************************************************************************* //

その他の設定については以下の記事を参考にしてください。

実行~可視化

Allrun.shを実行することで、カップリング計算を開始することが可能です。

CFDEMcouplingPFM/tutorials/cfdemIB/falling_sphere_two_way_coupling/Allrun.sh

また、post.shを実行することでopenfoamの結果ファイルをvtk形式に変換できます。このファイルをparaviewに読み込んで可視化します。

CFDEMcouplingPFM/tutorials/cfdemIB/falling_sphere_two_way_coupling/post.sh

paraviewから出力した動画を示します。
DEM粒子の落下計算で、落下時の流体-DEMの相互作用を双方向の連成で解析しています。

格子線も合わせると、DEM粒子よりもメッシュサイズが小さいことが分かります。

次回以降、チュートリアル以外の解析を実施してみようと思います。