回転ドラム内の混相流解析

2017年12月24日

はじめに

回転ドラム内の混相流解析の手順を示す。

使用バージョン

OpenFOAM v1706

ファイル

解析条件

モデルは直径 0.3 m、奥行き 0.05 m の円筒とする。中に流体を入れて側面をある回転速度で回転させて、流体の挙動を見る。流体は 2 つのパターンを考える。

  • 水-空気
  • 粉体-空気

円筒側面を "side"、表と裏を "frontAndBack" という名前にする。

メッシュ

メッシュは snappyHexMesh で作成する (参考)。

ジオメトリ

ジオメトリ (STL) は FreeCAD などで作成する (参考)。ただの円筒なので、ParaView でも作れる。

バックグラウンドセル (blockMesh)

blockMesh の範囲をジオメトリの範囲より少し大きめに作る。

system/blockMeshDict

minx -0.151;
maxx  0.151;
miny -0.151;
maxy  0.151;
minz -0.001;
maxz  0.051;

分割数は適当に設定する。

blockMesh を実行。

$ blockMesh

snappyHexMesh

simpleFoam のチュートリアルケースの motorBike をベースにする。

system/snappyHexMeshDict

geometry
{
    drum.stl
    {
        type triSurfaceMesh;
        name drum;
    }
};

...

castellatedMeshControls
{
    ...

    maxGlobalCells 1000;

    features
    (
        {
            file "drum.eMesh";
            level 6;
        }
    );

    refinementSurfaces
    {
        drum
        {
            level (0 0);
        }
    }

    refinementRegions
    {
    }

    locationInMesh (0 0 0.001);

    ...
}

...

addLayersControls
{
    ...

    layers
    {
        drum_side
        {
            nSurfaceLayers 3;
        }
    }

    expansionRatio 1.1;

		...
}

system/surfaceFeatureExtractDict

drum.stl
{
    ...
}

surfaceFeatureExtract、snappyHexMesh を実行。

surfaceFeatureExtract
snappyHexMesh -overwrite

境界名の設定 (createPatch)

createPatch で境界の名前とタイプを設定する (手で修正してもよいが、やり直しが面倒なので)。

system/createPatchDict

pointSync false;

patches
(
    {
        name side;
        patchInfo
        {
            type wall;
        }
        constructFrom patches;
        patches (drum_side);
    }

    {
        name frontAndBack;
        patchInfo
        {
            type wall;
        }
        constructFrom patches;
        patches (drum_front drum_back);
    }
);

STL から作られた境界名 "drum_side"、"drum_front"、"drum_back" を "side" と "frontAndBack" に変更している。境界タイプは "wall" としている。

水-空気系

水-空気系の計算は interFoam で行う。damBreak ケースをベースにする。メッシュの分割数は 50x50x10 とした。

setFields で水を溜める。

system/setFieldsDict

defaultFieldValues
(
    volScalarFieldValue alpha.water 0
);

regions
(
    boxToCell
    {
        box (-0.15 -0.15 0) (0.15 -0.05 0.05);
        fieldValues
        (
            volScalarFieldValue alpha.water 1
        );
    }
);

円筒側面 (side) の回転速度を rotatingWallVelocity で与える。ここでは回転速度 10 [rad/s] としている。表と裏 (frontAndBack) はすべり条件とする。

0/U

internalField   uniform (0 0 1e-6);

boundaryField
{
    side
    {
        type            rotatingWallVelocity;
        origin          (0 0 0);
        axis            (0 0 1);
        omega           constant 10;
    }

    frontAndBack
    {
        type            slip;
    }
}

setFields を実行。

$ cp 0/alpha.water.org 0/alpha.water
$ setFields

計算実行。

$ interFoam > log &

計算結果。

反時計回りに回転する側面によって水の膜が左端に移動し水面に突入するので、空気を巻き込んでいる。

粉体-空気系

粉体-空気系の計算は twoPhaseEulerFoam で行う。fluidizedBed ケースをベースにする。メッシュの分割数は 50x50x3 とした。

相の設定を constant/phaseProperties で行う。

constant/phaseProperties

particles
{
    diameterModel constant;
    constantCoeffs
    {
        d               10e-6;
    }

    alphaMax        0.62;
    residualAlpha   1e-6;
}

...

drag
(
    (particles in air)
    {
        type            SchillerNaumann;
        residualRe      1e-3;
        swarmCorrection
        {
            type        none;
        }
    }
);

粉体の径を 10 [μm] とし、抗力モデルを Schiller-Naumann モデルに設定している。

粉体の挙動を解くためのモデルを constant/turbulenceProperties.particle で設定する。

RAS
{
    RASModel kineticTheory;

    turbulence      on;
    printCoeffs     on;

    kineticTheoryCoeffs
    {
        equilibrium             on;

        e                       0.8;
        alphaMax                0.62;
        alphaMinFriction        0.5;
        residualAlpha           1e-6;

        viscosityModel          Gidaspow;
        conductivityModel       Gidaspow;
        granularPressureModel   Lun;
        frictionalStressModel   JohnsonJackson;
        radialModel             SinclairJackson;

        JohnsonJacksonCoeffs
        {
            Fr                      0.05;
            eta                     2;
            p                       5;
            phi                     28.5;
            alphaDeltaMin           0.05;
        }
    }

    phasePressureCoeffs
    {
        preAlphaExp     500;
        expMax          1000;
        alphaMax        0.62;
        g0              1000;
    }
}

RASModel を kineticTheory に設定している。これは粉体同士の接触を扱うためのものである。いろいろ細かいモデルがあるが、チュートリアルケースで指定してあるままとした。

setFields で粉体を溜める。

system/setFieldsDict

defaultFieldValues
(
    volScalarFieldValue alpha.particles 0
);

regions
(
    boxToCell
    {
        box (-0.15 -0.15 0) (0.15 -0.045 0.05);
        fieldValues
        (
            volScalarFieldValue alpha.particles 0.6
        );
    }
);

円筒側面 (side) の回転速度を rotatingWallVelocity で与える。ここでは回転速度 0.28 [rad/s] としている。表と裏 (frontAndBack) はすべり条件とする。

0/U.particles

internalField   uniform (0 0 0);

boundaryField
{
    side
    {
        type            rotatingWallVelocity;
        origin          (0 0 0);
        axis            (0 0 1);
        omega           constant 0.28;
    }

    frontAndBack
    {
        type            slip;
    }
}

0/U.air

internalField   uniform (0 0 0);

boundaryField
{
    side
    {
        type            rotatingWallVelocity;
        origin          (0 0 0);
        axis            (0 0 1);
        omega           constant 0.28;
    }

    frontAndBack
    {
        type            slip;
    }
}

setFields を実行。

$ cp 0/alpha.particles.org 0/alpha.particles
$ setFields

計算実行。

$ twoPhaseEulerFoam > log &

計算結果。

反時計回りに回転する側面によって粉体が右に持ち上がろうとするが、すべり落ちるので、傾いた状態になる。

参考文献を参考にして、回転速度や粉体の量をいろいろ変えると挙動が変わって面白い。

参考文献