【OpenFOAM】定常流れ＋粒子追跡（カップリング無し）

本記事はOpenFOAMでの粒子追跡の備忘録です。
自分で理解したことや試したことをまとめていきます。

本記事の内容

icoUncoupledKinematicParcelFoamを使って粒子追跡の解析
流体は一定の速度、粒子は球体形状とし流体の力を受けて運動する様子をシミュレーション

別の記事で球体にはたらく重力と抗力による質点の運動とOpenFOAMのicoUncoupledKinematicParcelFoamを使って得られた粒子の運動を比較を行う予定です。OpenFOAMの結果は厳密には壁との摩擦もありますし質点ではないので、質点と考えた場合の運動とずれが生じます。

粒子追跡のソルバ

こちらを見ると公式サイトほどソルバの数はありません。

自分が使っているバージョンがOpenFOAMv2012ですが、このバージョンだと以下のソルバが用意されています。

ls /opt/OpenFOAM/OpenFOAM-v2012/applications/solvers/lagrangian/

DPMFoam            icoUncoupledKinematicParcelFoam  reactingParcelFoam    sprayFoam
coalChemistryFoam  lagrangianSolversDoc.H           simpleCoalParcelFoam  uncoupledKinematicParcelFoam

ls /opt/OpenFOAM/OpenFOAM-v2012/applications/solvers/lagrangian/

DPMFoam icoUncoupledKinematicParcelFoam reactingParcelFoam sprayFoam

coalChemistryFoam lagrangianSolversDoc.H simpleCoalParcelFoam uncoupledKinematicParcelFoam

icoUncoupledKinematicParcelFoamを使います。

モデル作成

icoUncoupledKinematicParcelFoamが使わているチュートリアルをコピーしてきます。

cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying particle_base

1	cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying particle_base

「/opt/OpenFOAM/OpenFOAM-v2012/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper」というチュートリアルには「hopperInitialState」「hopperEmptying」の2つがありますが、

hopperInitialState：初期状態
hopperEmptying：本解析

と別れています。

どちらも粒子の運動を解くチュートリアルなのでどちらを使っても良いですが、今回は「hopperEmptying」を使います。

フォルダ移動します。

cd particle_base

1	cd particle_base

blockMeshでメッシュ生成

流体領域のメッシュはblockMeshで生成します。
blockMeshで作るメッシュの頂点の番号は以下のようにして、

Xmin：inlet（流入条件）
Xmin：inlet（流出条件）
その他：滑り無し（壁条件）

とします。

system/blockMeshDict

scale   0.001;

vertices
(
    (  0   0   0)
    (500   0   0)
    (500 100   0)
    (  0 100   0)
    (  0   0 100)
    (500   0 100)
    (500 100 100)
    (  0 100 100)
);


blocks
(
    hex (0 1 2 3 4 5 6 7) (50 10 10) simpleGrading (1 1 1)
);

boundary
(

    inlet
    {
        type patch;
        faces
        (
            (0 3 7 4)
        );
    }

    outlet
    {
        type patch;
        faces
        (
            (1 5 6 2)
        );
    }

    walls
    {
        type wall;
        faces
        (
            (0 1 2 3)
            (6 7 4 5)
            (0 1 5 4)
            (2 6 7 3)
        );
    }
);

scale 0.001;

vertices

(

( 0 0 0)

(500 0 0)

(500 100 0)

( 0 100 0)

( 0 0 100)

(500 0 100)

(500 100 100)

( 0 100 100)

);

blocks

(

hex (0 1 2 3 4 5 6 7) (50 10 10) simpleGrading (1 1 1)

);

boundary

(

inlet

{

type patch;

faces

(

(0 3 7 4)

);

}

outlet

{

type patch;

faces

(

(1 5 6 2)

);

}

walls

{

type wall;

faces

(

(0 1 2 3)

(6 7 4 5)

(0 1 5 4)

(2 6 7 3)

);

}

);

以下のコマンドでblockMeshを実行します。

blockMesh

blockMesh

できたメッシュはParaViewで確認してみましょう。

問題なければ次に境界条件の設定を行います。

境界条件の設定

境界条件は後で元に戻せるように「0.orig」で設定します。

0.orig/U

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

internalField   uniform (2.0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           $internalField;
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            noSlip;
    }
}

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

internalField uniform (2.0 0 0);

boundaryField

{

inlet

{

type fixedValue;

value $internalField;

}

outlet

{

type zeroGradient;

}

walls

{

type noSlip;

}

次に流体の密度を設定します。

0.orig/rho

dimensions      [1 -3 0 0 0 0 0];

internalField   uniform 1.2;

boundaryField
{
    ".*"
    {
        type            calculated;
        value           uniform 1.2;
    }
}

dimensions [1 -3 0 0 0 0 0];

internalField uniform 1.2;

boundaryField

{

".*"

{

type calculated;

value uniform 1.2;

}

流体の粘性係数を設定します。

0.orig/mu

cp -r 0.orig 0

1	cp -r 0.orig 0

乱流モデル（laminar）

constant/turbulenceProperties

simulationType  laminar;

1	simulationType laminar;

乱流モデルは使わない設定です。

物性値の設定

constant/transportProperties

rhoInf          1.2;

transportModel  Newtonian;

nu              1e-05;

rhoInf 1.2;

transportModel Newtonian;

nu 1e-05;

「0/rho」や「0/mu」でも密度$rho=1.2$、粘性係数$mu=\rho \nu = 1.2\times 10^{-5}$としているのに「constant/transportProperties」でも設定していて注意しないと二重に設定してしまっていないか、両者で食い違った設定になっていないかが気になりますね。

ソースコードを少し確認してどうなっているのかを理解しておきましょう。

ソースコードを確認

こちらより確認します。
もしくはこちら( icoUncoupledKinematicParcelFoam.C)。

icoUncoupledKinematicParcelFoam.C

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "basicKinematicCollidingCloud.H"

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

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Transient solver for the passive transport"
        " of a single kinematic particle cloud"
    );
    argList::addOption
    (
        "cloud",
        "name",
        "specify alternative cloud name. default is 'kinematicCloud'"
    );

    #include "postProcess.H"

    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createFields.H"

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

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        Info<< "Evolving " << kinematicCloud.name() << endl;

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

        kinematicCloud.evolve();

        runTime.write();

        runTime.printExecutionTime(Info);
    }

    Info<< "End\n" << endl;

    return 0;
}

#include "fvCFD.H"

#include "singlePhaseTransportModel.H"

#include "turbulentTransportModel.H"

#include "basicKinematicCollidingCloud.H"

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

int main(int argc, char *argv[])

{

argList::addNote

(

"Transient solver for the passive transport"

" of a single kinematic particle cloud"

);

argList::addOption

(

"cloud",

"name",

"specify alternative cloud name. default is 'kinematicCloud'"

);

#include "postProcess.H"

#include "addCheckCaseOptions.H"

#include "setRootCaseLists.H"

#include "createTime.H"

#include "createMesh.H"

#include "createControl.H"

#include "createFields.H"

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

Info<< "\nStarting time loop\n" << endl;

while (runTime.loop())

{

Info<< "Time = " << runTime.timeName() << nl << endl;

Info<< "Evolving " << kinematicCloud.name() << endl;

laminarTransport.correct();

mu = laminarTransport.nu()*rhoInfValue;

kinematicCloud.evolve();

runTime.write();

runTime.printExecutionTime(Info);

}

Info<< "End\n" << endl;

return 0;

}

メインのファイルはたったこれだけ？って感じです。

物理量の設定は「createFields.H」というファイル名で行っているので確認してみましょう。

createFields.H

#include "readGravitationalAcceleration.H"

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

#include "createPhi.H"

singlePhaseTransportModel laminarTransport(U, phi);

dimensionedScalar rhoInfValue
(
    "rhoInf",
    dimDensity,
    laminarTransport
);

volScalarField rhoInf
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    rhoInfValue
);

autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

volScalarField mu
(
    IOobject
    (
        "mu",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    laminarTransport.nu()*rhoInfValue
);

const word kinematicCloudName
(
    args.getOrDefault<word>("cloud", "kinematicCloud")
);

Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
basicKinematicCollidingCloud kinematicCloud
(
    kinematicCloudName,
    rhoInf,
    U,
    mu,
    g
);

IOobject Hheader
(
    "H",
    runTime.timeName(),
    mesh,
    IOobject::MUST_READ,
    IOobject::AUTO_WRITE
);

autoPtr<volVectorField> HPtr;

if (Hheader.typeHeaderOk<volVectorField>(true))
{
    Info<< "\nReading field H\n" << endl;

    HPtr.reset(new volVectorField (Hheader, mesh));
}

IOobject HdotGradHheader
(
    "HdotGradH",
    runTime.timeName(),
    mesh,
    IOobject::MUST_READ,
    IOobject::AUTO_WRITE
);

autoPtr<volVectorField> HdotGradHPtr;

if (HdotGradHheader.typeHeaderOk<volVectorField>(true))
{
    Info<< "Reading field HdotGradH" << endl;

    HdotGradHPtr.reset(new volVectorField(HdotGradHheader, mesh));
}

#include "createNonInertialFrameFields.H"す

100

101

102

103

104

105

106

107

108

109

110

111

#include "readGravitationalAcceleration.H"

Info<< "Reading field U\n" << endl;

volVectorField U

(

IOobject

(

"U",

runTime.timeName(),

mesh,

IOobject::MUST_READ,

IOobject::AUTO_WRITE

mesh

);

#include "createPhi.H"

singlePhaseTransportModel laminarTransport(U, phi);

dimensionedScalar rhoInfValue

(

"rhoInf",

dimDensity,

laminarTransport

);

volScalarField rhoInf

(

IOobject

(

"rho",

runTime.timeName(),

mesh,

IOobject::NO_READ,

IOobject::AUTO_WRITE

mesh,

rhoInfValue

);

autoPtr<incompressible::turbulenceModel> turbulence

(

incompressible::turbulenceModel::New(U, phi, laminarTransport)

);

volScalarField mu

(

IOobject

(

"mu",

runTime.timeName(),

mesh,

IOobject::NO_READ,

IOobject::AUTO_WRITE

laminarTransport.nu()*rhoInfValue

);

const word kinematicCloudName

(

args.getOrDefault<word>("cloud", "kinematicCloud")

);

Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;

basicKinematicCollidingCloud kinematicCloud

(

kinematicCloudName,

rhoInf,

mu,

);

IOobject Hheader

(

"H",

runTime.timeName(),

mesh,

IOobject::MUST_READ,

IOobject::AUTO_WRITE

);

autoPtr<volVectorField> HPtr;

if (Hheader.typeHeaderOk<volVectorField>(true))

{

Info<< "\nReading field H\n" << endl;

HPtr.reset(new volVectorField (Hheader, mesh));

}

IOobject HdotGradHheader

(

"HdotGradH",

runTime.timeName(),

mesh,

IOobject::MUST_READ,

IOobject::AUTO_WRITE

);

autoPtr<volVectorField> HdotGradHPtr;

if (HdotGradHheader.typeHeaderOk<volVectorField>(true))

{

Info<< "Reading field HdotGradH" << endl;

HdotGradHPtr.reset(new volVectorField(HdotGradHheader, mesh));

}

#include "createNonInertialFrameFields.H"す

少し詳しく見てみます。
以下で流速の「U」の定義を行っていますね。

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field U\n" << endl;

volVectorField U

(

IOobject

(

"U",

runTime.timeName(),

mesh,

IOobject::MUST_READ,

IOobject::AUTO_WRITE

mesh

);

「 runTime.timeName()」が時間に応じて読み込むファイルを指定している部分です。
「IOobject::MUST_READ」で毎回時間に応じて読み込むという意味で、「IOobject::AUTO_WRITE」が書き込みの指定です。
つまり、「U」は値をファイルに書き込んで出力するということです。

続いて以下で何やら定義しています。

singlePhaseTransportModel laminarTransport(U, phi);

1	singlePhaseTransportModel laminarTransport(U, phi);

これが良くわかりませんが、調べると「singlePhaseTransportModel.H」のコンストラクタのようです。
コンストラクタはクラスのインスタンス化した際に呼び出される関数で初期設定のようなものです。

    // Constructors
 
        //- Construct from components
        singlePhaseTransportModel
        (
            const volVectorField& U,
            const surfaceScalarField& phi
        );

// Constructors

//- Construct from components

singlePhaseTransportModel

(

const volVectorField& U,

const surfaceScalarField& phi

);

.Hのようなヘッダーファイルは関数の定義しか書いていないので、「singlePhaseTransportModel.C」で内容を確認してみます。

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::singlePhaseTransportModel::singlePhaseTransportModel
(
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    IOdictionary
    (
        IOobject
        (
            "transportProperties",
            U.time().constant(),
            U.db(),
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    ),
    viscosityModelPtr_(viscosityModel::New("nu", *this, U, phi))
{}

// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

Foam::singlePhaseTransportModel::singlePhaseTransportModel

(

const volVectorField& U,

const surfaceScalarField& phi

)

IOdictionary

(

IOobject

(

"transportProperties",

U.time().constant(),

U.db(),

IOobject::MUST_READ_IF_MODIFIED,

IOobject::NO_WRITE

)

viscosityModelPtr_(viscosityModel::New("nu", *this, U, phi))

{}

これを見ると「U.time().constant()」は「constant」フォルダの意味で「”transportProperties”」が読み込むファイル名なので、結局「constant/transportProperties」の「nu」の値を使っているように思えます。

続いて「rhoInfValue」を次元付きで定義しています。

dimensionedScalar rhoInfValue
(
    "rhoInf",
    dimDensity,
    laminarTransport
);

dimensionedScalar rhoInfValue

(

"rhoInf",

dimDensity,

laminarTransport

);

先ほど「laminarTransport」は「constant/transportProperties」から定義されたsinglePhaseTransportModelクラスなので、おそらくここでの「rhoInf」は「constant/transportProperties」で定義された「rhoInf」のことでしょう。
注意すべきは値は「constant/transportProperties」で定義された「rhoInf」ですが、変数名はrhoInfValueと定義している点です。

ここまででわかったことは、

nu
rhoInfValue

「constant/transportProperties」で定義された値を使っているということです。

続いて、

volScalarField rhoInf
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    rhoInfValue
);

autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

volScalarField mu
(
    IOobject
    (
        "mu",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    laminarTransport.nu()*rhoInfValue
);

volScalarField rhoInf

(

IOobject

(

"rho",

runTime.timeName(),

mesh,

IOobject::NO_READ,

IOobject::AUTO_WRITE

mesh,

rhoInfValue

);

autoPtr<incompressible::turbulenceModel> turbulence

(

incompressible::turbulenceModel::New(U, phi, laminarTransport)

);

volScalarField mu

(

IOobject

(

"mu",

runTime.timeName(),

mesh,

IOobject::NO_READ,

IOobject::AUTO_WRITE

laminarTransport.nu()*rhoInfValue

);

「rhoInf」も「mu」も「IOobject::NO_READ」となっていることから、「0」フォルダに「rho」「mu」があるのですがそれは読み込まないという意味になっています。
つまり「0」フォルダに「rho」と「mu」もあって、「constant/transportProperties」にもrhoIhfとnuを読み込む設定だと2重に定義していることになっていないかという疑問に対しては、結局「constant/transportProperties」でのrhoIhfとnuの値しか使っていないということです。

ですので「0」フォルダのrhoとmuの設定ファイルは消しても計算はできます。
なんでrhoとmuの設定ファイルを用意しているのかわかりませんが（解析結果には影響ないと思いますがどうでしょう…）

また、muについては「laminarTransport.nu()*rhoInfValue」で定義しています。
どちらも「constant/transportProperties」でのrhoIhfとnuの値を使って「mu」を定義しているということですね。

粒子の設定

粒子の設定は「constant/kinematicCloudProperties」で行います。

constant/kinematicCloudProperties

solution
{
    active          true;
    coupled         false;
    transient       yes;
    cellValueSourceCorrection off;
    maxCo           0.3;

    sourceTerms
    {
        schemes
        {
        }
    }

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

    integrationSchemes
    {
        U               Euler;
    }
}


constantProperties
{
    rho0            964;
    youngsModulus   6e8;
    poissonsRatio   0.35;
}


subModels
{
    particleForces
    {
        sphereDrag;
        gravity;
    }

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

    dispersionModel none;

    patchInteractionModel standardWallInteraction;//none;

    surfaceFilmModel none;

    stochasticCollisionModel none;

    collisionModel pairCollision;

    pairCollisionCoeffs
    {
        // Maximum possible particle diameter expected at any time
        maxInteractionDistance  0.006;

        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;
            walls
            {
                youngsModulus   1e10;
                poissonsRatio   0.23;
                alpha           0.12;
                b               1.5;
                mu              0.43;
                cohesionEnergyDensity 0;
            }
        };
    }
    standardWallInteractionCoeffs
    {
        type            rebound;
    }
}


cloudFunctions
{}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

solution

{

active true;

coupled false;

transient yes;

cellValueSourceCorrection off;

maxCo 0.3;

sourceTerms

{

schemes

{

}

interpolationSchemes

{

rho cell;

U cellPoint;

mu cell;

}

integrationSchemes

{

U Euler;

}

constantProperties

{

rho0 964;

youngsModulus 6e8;

poissonsRatio 0.35;

}

subModels

{

particleForces

{

sphereDrag;

gravity;

}

injectionModels

{

model1

{

type manualInjection;

massTotal 0;

parcelBasisType fixed;

nParticle 1;

SOI 0;

positionsFile "kinematicCloudPositions";

U0 ( 0 0 0 );

sizeDistribution

{

type fixedValue;

fixedValueDistribution

{

value 0.01;

}

dispersionModel none;

patchInteractionModel standardWallInteraction;//none;

surfaceFilmModel none;

stochasticCollisionModel none;

collisionModel pairCollision;

pairCollisionCoeffs

{

// Maximum possible particle diameter expected at any time

maxInteractionDistance 0.006;

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;

walls

{

youngsModulus 1e10;

poissonsRatio 0.23;

alpha 0.12;

b 1.5;

mu 0.43;

cohesionEnergyDensity 0;

}

};

}

standardWallInteractionCoeffs

{

type rebound;

}

cloudFunctions

{}

今回は粒子に働く力として以下を設定しています。

y方向に重力
粒子を球体とみなした抗力

粒子の初期位置は「constant/kinematicCloudPositions」で行い「constant/kinematicCloudProperties」で読み込んでいます。

constant/kinematicCloudPositions

(
(0.25  0.05 0.05)
)

(

(0.25 0.05 0.05)

)

複数の粒子を設定したい場合は、リストとして座標を複数設定します。

球体に働く抗力

粒子には球体に働く抗力が設定されています。
その理論式は以下で確認することができます。

SphereDragForce< CloudType >

SphereDragForce.C

template<class CloudType>
Foam::scalar Foam::SphereDragForce<CloudType>::CdRe(const scalar Re) const
{
    // (AOB:Eq. 35)
    if (Re > 1000.0)
    {
        return 0.424*Re;
    }
 
    return 24.0*(1.0 + (1.0/6.0)*pow(Re, 2.0/3.0));
}

template<class CloudType>

Foam::scalar Foam::SphereDragForce<CloudType>::CdRe(const scalar Re) const

{

// (AOB:Eq. 35)

if (Re > 1000.0)

{

return 0.424*Re;

}

return 24.0*(1.0 + (1.0/6.0)*pow(Re, 2.0/3.0));

}

template<class CloudType>
Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled
(
    const typename CloudType::parcelType& p,
    const typename CloudType::parcelType::trackingData& td,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const
{
    // (AOB:Eq. 34)
    return forceSuSp(Zero, mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d())));
}

template<class CloudType>

Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled

(

const typename CloudType::parcelType& p,

const typename CloudType::parcelType::trackingData& td,

const scalar dt,

const scalar mass,

const scalar Re,

const scalar muc

) const

{

// (AOB:Eq. 34)

return forceSuSp(Zero, mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d())));

}

KinematicCloud< CloudType >

解析実行

解析実行とVTKファイルの出力を行います。

icoUncoupledKinematicParcelFoam;foamToVTK

1	icoUncoupledKinematicParcelFoam;foamToVTK

「;(セミコロン)」でつなぐとコマンドを続けて実行してくれます。

解析が終わったら以下の手順でParaViewで可視化してみましょう。

「.foam」の空ファイルの読み込み
「VTK/lagrangian/kinematicCloud/kinematicCloud_*」の読み込み
Glyph>Glyph TypeをSphere

あとはアニメーションなどで挙動を確認。

まとめ

本記事はicoUncoupledKinematicParcelFoamソルバを使った粒子追跡について解説しました。

別の記事で球体にはたらく重力と抗力による質点の運動とOpenFOAMのicoUncoupledKinematicParcelFoamを使って得られた粒子の運動を比較を行う予定です。OpenFOAMの結果は厳密には壁との摩擦もありますし質点ではないので、質点と考えた場合の運動とずれが生じます。シミュレーションも内部でどのような計算が行われているかユーザはわからないのですが、本記事のような基本的な理解からはじめるのが良いでしょう。

粒子追跡で使えるソルバ（おまけ）

ESI版の粒子追跡によるソルバの種類を確認します。

coalChemistryFoam	Transient solver for compressible, turbulent flow, with coal and limestone particle clouds, an energy source, and combustion	エネルギー源である石炭や石灰岩の粒子雲、燃焼を伴う圧縮性乱流のトランジェントソルバ
DPMFoam	Transient solver for the coupled transport of a single kinematic particle cloud including the effect of the volume fraction of particles on the continuous phase	連続相に対する粒子の体積分率の影響を含む単一の運動学的粒子雲の連成輸送のための過渡ソルバー
DPMDyMFoam	Transient solver for the coupled transport of a single kinematic particle cloud including the effect of the volume fraction of particles on the continuous phase, with optional mesh motion and mesh topology changes	連続相に対する粒子の体積分率の影響を含む単一の運動学的粒子雲の連成輸送のための過渡ソルバーで、オプションでメッシュの動きとメッシュのトポロジーを変更することができます。
MPPICDyMFoam	Transient solver for the coupled transport of a single kinematic particle cloud including the effect of the volume fraction of particles on the continuous phase. Multi-Phase Particle In Cell (MPPIC) modeling is used to represent collisions without resolving particle-particle interactions, with optional mesh motion and mesh topology changes	連続相に対する粒子の体積分率の影響を含む、単一の運動学的粒子雲の連成輸送のための過渡ソルバーです。粒子間相互作用を解決せずに衝突を表現するMPPIC（Multi-Phase Particle In Cell）モデリングは、オプションでメッシュの動きとメッシュのトポロジーを変更することができます。
MPPICFoam	Transient solver for the coupled transport of a single kinematic particle cloud including the effect of the volume fraction of particles on the continuous phase. Multi-Phase Particle In Cell (MPPIC) modeling is used to represent collisions without resolving particle-particle interactions	連続相に対する粒子の体積分率の影響を含む、単一の運動学的粒子雲の連成輸送のための過渡ソルバーです。粒子間相互作用を解決せずに衝突を表現するためにMPPIC（Multi-Phase Particle In Cell）モデリングが使用されています。
icoUncoupledKinematicParcelFoam	Transient solver for the passive transport of a single kinematic particle cloud	単一運動論的粒子雲の受動輸送のための過渡ソルバー
icoUncoupledKinematicParcelDyMFoam	Transient solver for the passive transport of a single kinematic particle cloud, with optional mesh motion and mesh topology changes	単一の運動学的粒子雲の受動的輸送のためのトランジェントソルバ、オプションでメッシュの動きとメッシュのトポロジーを変更可能
reactingParcelFoam	Transient solver for compressible, turbulent flow with a reacting, multiphase particle cloud, and surface film modelling	反応する多相粒子雲を伴う圧縮性乱流の過渡ソルバー、および表面膜モデリング
reactingHeterogenousParcelFoam	Transient solver for the coupled transport of a single kinematic particle cloud including the effect of the volume fraction of particles on the continuous phase. Multi-Phase Particle In Cell (MPPIC) modeling is used to represent collisions without resolving particle-particle interactions	連続相に対する粒子の体積分率の影響を含む、単一の運動学的粒子雲の連成輸送のための過渡ソルバーです。粒子間相互作用を解決せずに衝突を表現するためにMPPIC（Multi-Phase Particle In Cell）モデリングが使用されています。
simpleReactingParcelFoam	Steady-state solver for compressible, turbulent flow with reacting, multiphase particle clouds and optional sources/constraints	反応する多相の粒子雲とオプションのソース/制約を持つ圧縮性乱流の定常状態ソルバー
simpleCoalParcelFoam	Steady-state solver for compressible, turbulent flow with coal particle clouds and optional sources/constraints	石炭粒子雲とオプションのソース/制約を含む圧縮性乱流の定常状態ソルバー
sprayFoam	Transient solver for compressible, turbulent flow with a spray particle cloud	噴霧粒子雲を有する圧縮性乱流の過渡ソルバー
engineFoam	Transient solver for compressible, turbulent engine flow with a spray particle cloud	噴霧粒子雲を有する圧縮性乱流エンジン流れの過渡ソルバー
simpleSprayFoam	Steady state solver for compressible, turbulent flow with a spray particle cloud and optional sources/constraints	スプレー粒子雲とオプションのソース/制約を持つ圧縮性乱流の定常ソルバー
sprayDyMFoam	Transient solver for compressible, turbulent flow with a spray particle cloud, with optional mesh motion and mesh topology changes	スプレー粒子雲を伴う圧縮性乱流のトランジェントソルバ、オプションでメッシュの動きやメッシュトポロジーの変更も可能
uncoupledKinematicParcelFoam	Transient solver for the passive transport of a particle cloud	粒子雲の受動輸送のためのトランジェントソルバ
uncoupledKinematicParcelDyMFoam	Transient solver for the passive transport of a particle cloud	粒子雲の受動輸送のためのトランジェントソルバ