DEM 解析

2016年12月9日

はじめに

OpenFOAM で DEM の計算を行う。

使用バージョン

OpenFOAM 4.1

ファイル

DEM ソルバー

  • icoUncoupledKinematicParcelFoam
  • DPMFoam

参考

  • $FOAM_TUTORIALS/lagrangian/icoUncoupledKinematicParcelFoam/hopper
  • $FOAM_TUTORIALS/lagrangian/DPMFoam/Goldschmidt

DEM 解析

icoUncoupledKinematicParcelFoam により DEM 計算を行うことができる。

問題

正方形領域に玉を 5 個配置して落としてみる。

玉の初期配置

constant/kinematicCloudPositions に記述する。

(
(0.025  0.075 0.005)
(0.05   0.075 0.005)
(0.075  0.075 0.005)
(0.0375 0.05  0.005)
(0.0625 0.05  0.005)
)

衝突の設定

constant/kinematicCloudProperties で衝突についての設定を行う。

玉の特性

constantProperties
{
    parcelTypeId    1;

    rhoMin          1e-15;
    minParticleMass 1e-15;

    rho0            964;
    youngsModulus   6e8;
    poissonsRatio   0.35;

    constantVolume  false;
}

粒子にかかる力

subModels
{
    particleForces
    {
        sphereDrag;
        gravity;
    }
    
    ...
    
}

粒子にかかる力には、たとえば次のようなものを指定できる。

  • gravity
  • virtualMass
  • sphereDrag
  • nonSphereDrag
  • WenYuDrag
  • ErgunWenYuDrag
  • PlessisMasliyahDrag
  • SaffmanMeiLiftForce
  • TomiyamaLift

初期配置ファイルの指定と玉の径

subModels
{

    ...

    injectionModels
    {
        model1
        {
            type manualInjection;
            massTotal       0;
            parcelBasisType fixed;
            nParticle       1;
            SOI             0;
            positionsFile   "kinematicCloudPositions";
            U0              ( 0 0 0 );
            sizeDistribution
            {
                type        fixedValue;
                fixedValueDistribution
                {
                    value   0.02;
                }
            }
        }
    }

    ...

}

sizeDistribution で玉の径を指定している。

粒子同士の衝突の設定

subModels
{

    ...
    
    collisionModel pairCollision;
    
    ...
    
    pairCollisionCoeffs
    {
        // Maximum possible particle diameter expected at any time
        maxInteractionDistance  0.02;

        writeReferredParticleCloud no;

        pairModel pairSpringSliderDashpot;

        pairSpringSliderDashpotCoeffs
        {
            useEquivalentSize   no;
            alpha               0.12;
            b                   1.5;
            mu                  0.52;
            cohesionEnergyDensity 0;
            collisionResolutionSteps 12;
        };

        wallModel wallLocalSpringSliderDashpot;

        wallLocalSpringSliderDashpotCoeffs
        {
            useEquivalentSize no;
            collisionResolutionSteps 12;
            movingWall
            {
                youngsModulus   1e10;
                poissonsRatio   0.23;
                alpha           0.12;
                b               1.5;
                mu              0.43;
                cohesionEnergyDensity 0;
            }
            fixedWalls
            {
                youngsModulus   1e10;
                poissonsRatio   0.23;
                alpha           0.12;
                b               1.5;
                mu              0.43;
                cohesionEnergyDensity 0;
            }
            frontAndBack
            {
                youngsModulus   1e10;
                poissonsRatio   0.23;
                alpha           0.12;
                b               1.5;
                mu              0.1;
                cohesionEnergyDensity 0;
            }
        };
    }
    
    ...
    
}

maxInteractionDistance は最大衝突径で、本ケースでは玉の径と同じ。wallLocalSpringSliderDashpotCoeffs ではモデルの境界分の設定が必要である (モデルが変わるたびに書き換える必要がある)。境界ごとに設定を変える必要がなければ、wallModel に wallSpringSliderDashpot を指定して設定を 1 つで済ませることもできる。

        wallModel wallSpringSliderDashpot;

        wallSpringSliderDashpotCoeffs
        {
            useEquivalentSize no;
            collisionResolutionSteps 12;
            youngsModulus   1e10;
            poissonsRatio   0.23;
            alpha           0.12;
            b               1.5;
            mu              0.43;
            cohesionEnergyDensity 0;
        };

粒子と境界の衝突の設定

subModels
{

    ...
    
    patchInteractionModel standardWallInteraction;
    
    ...
    
    standardWallInteractionCoeffs
    {
        type            rebound;
    }

}

standardWallInteraction は、すべての境界を同種の壁として設定を行う。type には以下のものを設定できる。

  • rebound (反発)
  • stick (くっつく)
  • escape (消える)

境界ごとに設定を変えたい場合は、patchInteractionModel で localInteraction を指定する。

    patchInteractionModel localInteraction;
    
    localInteractionCoeffs
    {
        patches
        (
            ".*"
            {
                type rebound;
                e    0.97;
                mu   0.09;
            }
        );
    }

実行

$ cd square
$ blockMesh
$ icoUncoupledKinematicParcelFoam

ポスト処理

VTK に変換して処理する。

$ foamToVTK

ParaView を立ち上げて、VTK/*.vtk と、VTK/lagrangian/kinematicCloud/*.vtk を読み込む。lagrangian データを Glyph で Glyph Type を Sphere で表示すると、それっぽい表示が得られる。

今回の場合、Glyph の設定を以下のようにするときれいな表示になった。

  • Glyph Type: Sphere
  • Scale Mode: scalar
  • Set Scale Factor: 1
  • Scalars: d
  • Vectors: U
  • Radius 0.52
  • Theta Resolution: 20
  • Phi Resolution: 20

計算上の注意

  • 玉は設定次第で簡単に重なるが、重なると大きな反発力が働いて玉が高速で飛んでいくことになるため、重ならないように注意する。
  • 玉をきれいに配置すると、きれいな動きしかしない。たとえば縦横きれいに並べて玉を配置して落とすと、玉は上下方向に動くだけである。より現実的な動きにするためには、配置に少々不整を入れてやる必要がある。

連続相との連成

icoUncoupledKinematicParcelFoam では連続相 (速度場) との連成が行える。連成を有効にするには constant/kinematicCloudProperties で "coupled" を "true" にする。また、ソース項の設定を追加する。

solution
{
    active          true;
    coupled         true;
    transient       yes;
    cellValueSourceCorrection off;

    interpolationSchemes
    {
        rho             cell;
        U               cellPoint;
        mu              cell;
    }

    integrationSchemes
    {
        U               Euler;
    }

    sourceTerms
    {
        schemes
	{
            U           explicit 1;
        }
    }
}

U のスキームには "explicit" と "semiImplicit" が選べる。スキームの後の数値は緩和係数らしい。

constant/transportProperties で連続相の密度と粘性係数が設定できる。

rhoInf          rhoInf [ 1 -3 0 0 0 0 0 ] 1.2;

transportModel  Newtonian;

nu              nu [ 0 2 -1 0 0 0 0 ] 1e-05;

0/U に速度場を設定することで流れの影響を考慮することができる。浮力の影響は連続相との連成を有効にしなくても考慮されるようである。

CFD-DEM 解析

DPMFoam を用いると、流体と DEM の連成計算を行える。

設定は icoUncoupledKinematicParcelFoam とほぼ同様である。実行は次のように行う。

$ cd square-cfd-dem
$ blockMesh
$ DPMFoam

古い