非圧縮性ブシネスク近似乱流モデル用のライブラリにLanderSharmaKEを追加
shin
OpenFOAM初心者、かつCFD初心者の山田と申します。
今井雅さんが作成した、
非圧縮性ブシネスク近似乱流モデルのライブラリに
LanderSharmaKEを追加しようと試みております。
Giroさんのサイトを参考にさせて頂きながら、
どうにかこうにかしてやってみて、
計算が走りだすまで至ったのですが、
自分のやり方であっているのかどうか、
まったく自信がありません。
そこで今回
私のLanderSharmaKE追加までの一連の流れについて
皆さんのご意見を聞かせて頂ければと思い、
投稿させて頂きました。
一連の流れを以下に示します。
１ incompressble内のLaunderSharmaKEをboussinesqにコピー
cp -r $HOME/OpenFOAM/OpenFOAM-1.5/src/turbulenceModels/RAS/
incompressible/LaunderSharmaKE $HOME/OpenFOAM/OpenFOAM-1.5/src/
turbulenceModels/RAS/boussinesq
２ $HOME/OpenFOAM/OpenFOAM-1.5/src/turbulenceModels/RAS/boussinesq/
LaunderSharmaKE/LaunderSharmaKE.C を以下の箇所のみを変更する。
\*---------------------------------------------------------------------------
*/
#include "LaunderSharmaKE.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
namespace Foam
{
namespace boussinesq ←変更箇所
{
namespace RASModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * *
* * * //
defineTypeNameAndDebug(LaunderSharmaKE, 0);
addToRunTimeSelectionTable(RASModel, LaunderSharmaKE, dictionary);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * *
* * * //
tmp<volScalarField> LaunderSharmaKE::fMu() const
{
return exp(-3.4/sqr(scalar(1) + sqr(k_)/(nu()*epsilonTilda_)/
50.0));
}
tmp<volScalarField> LaunderSharmaKE::f2() const
{
return
scalar(1)
- 0.3*exp(-min(sqr(sqr(k_)/(nu()*epsilonTilda_)), scalar
(50.0)));
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * *
* * * //
LaunderSharmaKE::LaunderSharmaKE
(
const volVectorField& U,
const surfaceScalarField& phi,
const volScalarField& T, ←変更箇所
transportModel& lamTransportModel
)
:
RASModel(typeName, U, phi, T, lamTransportModel), ←変更箇
所
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
coeffDict_,
0.09
)
),
C1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C1",
coeffDict_,
1.44
)
),
C2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C2",
coeffDict_,
1.92
)
),
Cb_ ←変更箇所
( ←変更箇所
dimensioned<scalar>::lookupOrAddToDict ←変更箇所
( ←
変更箇所
"Cb",
←変更箇所
coeffDict_,
←変更箇所
1.44
←変更箇所
)
←変更箇所
),
←変更箇所
alphaEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaEps",
coeffDict_,
0.76923
)
),
transportProperties
←変更箇所
( ←
変更箇所
IOobject
←変更箇所
( ←
変更箇所
"transportProperties",
←変更箇所
"constant",
←変更箇所
mesh_,
←変更箇所
IOobject::MUST_READ,
←変更箇所
IOobject::NO_WRITE
←変更箇所
)
←変更箇所
),
←変更箇所
beta(transportProperties.lookup
("beta")),
←変更箇所
Prt(transportProperties.lookup
("Prt")),
←変更箇所
environmentalProperties
←変更箇所
( ←
変更箇所
IOobject
←変更箇所
( ←
変更箇所
"environmentalProperties",
←変更箇所
"constant",
←変更箇所
mesh_,
←変更箇所
IOobject::MUST_READ,
←変更箇所
IOobject::NO_WRITE
←変更箇所
)
←変更箇所
),
←変更箇所
g(environmentalProperties.lookup
("g")),
←変更箇所
k_
(
IOobject
(
"k",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
epsilonTilda_
(
IOobject
(
"epsilon",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
nut_(Cmu_*fMu()*sqr(k_)/(epsilonTilda_ + epsilonSmall_))
{
printCoeffs();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * *
* * * //
tmp<volSymmTensorField> LaunderSharmaKE::R() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"R",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
k_.boundaryField().types()
)
);
}
tmp<volSymmTensorField> LaunderSharmaKE::devReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"devRhoReff",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
-nuEff()*dev(twoSymm(fvc::grad(U_)))
)
);
}
tmp<fvVectorMatrix> LaunderSharmaKE::divDevReff(volVectorField& U)
const
{
return
(
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))
);
}
bool LaunderSharmaKE::read()
{
if (RASModel::read())
{
Cmu_.readIfPresent(coeffDict_);
C1_.readIfPresent(coeffDict_);
C2_.readIfPresent(coeffDict_);
Cb_.readIfPresent
(coeffDict_);
←変更箇所
alphaEps_.readIfPresent(coeffDict_);
return true;
}
else
{
return false;
}
}
void LaunderSharmaKE::correct()
{
transportModel_.correct();
if (!turbulence_)
{
return;
}
RASModel::correct();
volScalarField S2 = 2*magSqr(symm(fvc::grad(U_)));
volScalarField G = nut_*S2;
volScalarField E = 2.0*nu()*nut_*fvc::magSqrGradGrad(U_);
volScalarField D = 2.0*nu()*magSqr(fvc::grad(sqrt(k_)));
volScalarField B = nut_ * beta * (g & (fvc::grad(T_))) /
Prt ; ←変更箇所
volScalarField Bepsilon = 0.5*(B + mag
(B));
←変更箇所
// Dissipation rate equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(epsilonTilda_)
+ fvm::div(phi_, epsilonTilda_)
- fvm::laplacian(DepsilonEff(), epsilonTilda_)
==
C1_*G*epsilonTilda_/k_
- fvm::Sp(C2_*f2()*epsilonTilda_/k_, epsilonTilda_)
+ E
+ Cb_ * Bepsilon *epsilonTilda_/
k_
←変更箇所
);
epsEqn().relax();
solve(epsEqn);
bound(epsilonTilda_, epsilon0_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::laplacian(DkEff(), k_)
==
G - fvm::Sp((epsilonTilda_ + D)/k_, k_) +
B ←変更箇所
);
kEqn().relax();
solve(kEqn);
bound(k_, k0_);
// Re-calculate viscosity
nut_ = Cmu_*fMu()*sqr(k_)/epsilonTilda_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
//
************************************************************************* //
３ ＄HOME/OpenFOAM/OpenFOAM-1.5/src/turbulenceModels/RAS/boussinesq/
wmake/file の変更。
/* RAS turbulence models */
RASModel/RASModel.C
RASModel/newRASModel.C
kEpsilon/kEpsilon.C
LaunderSharmaKE/LaunderSharmaKE.C
←変更箇所
/* Wall functions */
wallFunctions/nutWallFunctions/nutWallFunction/
nutWallFunctionFvPatchScalarField.C
wallFunctions/nutWallFunctions/nutStandardWallFunction/
nutStandardWallFunctionFvPatchScalarField.C
wallFunctions/nutWallFunctions/nutStandardRoughWallFunction/
nutStandardRoughWallFunctionFvPatchScalarField.C
/* Patch fields */
derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/
turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/
turbulentMixingLengthFrequencyInletFvPatchScalarField.C
LIB = $(FOAM_LIBBIN)/libboussinesqRASModels
４ cd ＄HOME/OpenFOAM/OpenFOAM-1.5/src/turbulenceModels/RAS
wmake libso boussinesq
とコマンド入力。
以上のような手順を行ったのですが、
どこか、不備や訂正箇所などありましたら、
ご指摘して頂けませんでしょうか？
shin
申し訳ありません。
先日投稿した内容に不備がありましたので、
再度、投稿致します。
１ incompressble内のLaunderSharmaKEをboussinesqにコピー
cp -r $HOME/OpenFOAM/OpenFOAM-1.5/src/turbulenceModels/RAS/
incompressible/LaunderSharmaKE $HOME/OpenFOAM/OpenFOAM-1.5/src/
turbulenceModels/RAS/boussinesq
２ $HOME/OpenFOAM/OpenFOAM-1.5/src/turbulenceModels/RAS/boussinesq/
LaunderSharmaKE/LaunderSharmaKE.C
それと同じ階層にあるLaunderSharmaKE.Hを
以下の箇所のみを変更する。
#######LaunderSharmaKE.C#######
\*---------------------------------------------------------------------------
*/
#include "LaunderSharmaKE.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
namespace Foam
{
namespace boussinesq //変更箇所
{
namespace RASModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * *
* * * //
defineTypeNameAndDebug(LaunderSharmaKE, 0);
addToRunTimeSelectionTable(RASModel, LaunderSharmaKE, dictionary);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * *
* * * //
tmp<volScalarField> LaunderSharmaKE::fMu() const
{
return exp(-3.4/sqr(scalar(1) + sqr(k_)/(nu()*epsilonTilda_)/
50.0));
}
tmp<volScalarField> LaunderSharmaKE::f2() const
{
return
scalar(1)
- 0.3*exp(-min(sqr(sqr(k_)/(nu()*epsilonTilda_)), scalar
(50.0)));
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * *
* * * //
LaunderSharmaKE::LaunderSharmaKE
(
const volVectorField& U,
const surfaceScalarField& phi,
const volScalarField& T, //変更箇所
transportModel& lamTransportModel
)
:
RASModel(typeName, U, phi, T, lamTransportModel), //変更箇所
Cb_ //変更箇所
(
dimensioned<scalar>::lookupOrAddToDict //変更箇所
( //変更箇所
"Cb", //変更箇所
coeffDict_, //変更箇所
1.44 //変更箇所
) //変更箇所
), //変更箇所
alphaEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaEps",
coeffDict_,
0.76923
)
),
transportProperties //変更箇所
( //変更箇所
IOobject //変更箇所
( //変更箇所
"transportProperties", //変更箇所
"constant", //変更箇所
mesh_, //変更箇所
IOobject::MUST_READ, //変更箇所
IOobject::NO_WRITE //変更箇所
) //変更箇所
), //変更箇所
beta(transportProperties.lookup("beta")), //変更箇所
Prt(transportProperties.lookup("Prt")), //変更箇所
environmentalProperties //変更箇所
( //変更箇所
IOobject //変更箇所
( //変更箇所
"environmentalProperties", //変更箇所
"constant", //変更箇所
mesh_, //変更箇所
IOobject::MUST_READ, //変更箇所
IOobject::NO_WRITE //変更箇所
) //変更箇所
), //変更箇所
g(environmentalProperties.lookup("g")), //変更箇所
nut_(Cmu_*fMu()*sqr(k_)/(epsilonTilda_ + epsilonSmall_))
{
printCoeffs();
}
Cb_.readIfPresent(coeffDict_); //変更箇所
alphaEps_.readIfPresent(coeffDict_);
return true;
}
else
{
return false;
}
}
void LaunderSharmaKE::correct()
{
transportModel_.correct();
if (!turbulence_)
{
return;
}
RASModel::correct();
volScalarField S2 = 2*magSqr(symm(fvc::grad(U_)));
volScalarField G = nut_*S2;
volScalarField E = 2.0*nu()*nut_*fvc::magSqrGradGrad(U_);
volScalarField D = 2.0*nu()*magSqr(fvc::grad(sqrt(k_)));
volScalarField B = nut_ * beta * (g & (fvc::grad(T_))) /Prt ; //変更箇
所
volScalarField Bepsilon = 0.5*(B + mag(B)); //変更箇所
// Dissipation rate equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(epsilonTilda_)
+ fvm::div(phi_, epsilonTilda_)
- fvm::laplacian(DepsilonEff(), epsilonTilda_)
==
C1_*G*epsilonTilda_/k_
- fvm::Sp(C2_*f2()*epsilonTilda_/k_, epsilonTilda_)
+ E
+ Cb_ * Bepsilon *epsilonTilda_/k_ //変更箇所
);
epsEqn().relax();
solve(epsEqn);
bound(epsilonTilda_, epsilon0_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::laplacian(DkEff(), k_)
==
G - fvm::Sp((epsilonTilda_ + D)/k_, k_) + B //変更箇所
);
kEqn().relax();
solve(kEqn);
bound(k_, k0_);
// Re-calculate viscosity
nut_ = Cmu_*fMu()*sqr(k_)/epsilonTilda_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
//
************************************************************************* //
#######LaunderSharmaKE.H#######
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
namespace Foam
{
namespace boussinesq //変更箇所
{
namespace RASModels
{
/
*---------------------------------------------------------------------------
*\
Class LaunderSharmaKE Declaration
\*---------------------------------------------------------------------------
*/
class LaunderSharmaKE
:
public RASModel
{
dimensionedScalar Cmu_;
dimensionedScalar C1_;
dimensionedScalar C2_;
dimensionedScalar Cb_; //変更箇所
dimensionedScalar alphaEps_;
IOdictionary transportProperties; //変更箇所
dimensionedScalar beta; //変更箇所
dimensionedScalar Prt; //変更箇所
IOdictionary environmentalProperties; //変更箇所
dimensionedVector g; //変更箇所
volScalarField k_;
volScalarField epsilonTilda_;
volScalarField nut_;
tmp<volScalarField> fMu() const;
tmp<volScalarField> f2() const;
public:
TypeName("LaunderSharmaKE");
LaunderSharmaKE
(
const volVectorField& U,
const surfaceScalarField& phi,
const volScalarField& T, //変更箇所
transportModel& transport
);
~LaunderSharmaKE()
{}
tmp<volScalarField> nut() const
{
return nut_;
}
tmp<volScalarField> DkEff() const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", nut_ + nu())
);
}
//- Return the effective diffusivity for epsilon
tmp<volScalarField> DepsilonEff() const
{
return tmp<volScalarField>
(
new volScalarField("DepsilonEff", alphaEps_*nut_ + nu
())
);
}
//- Return the turbulence kinetic energy
tmp<volScalarField> k() const
{
return k_;
}
//- Note that epsilonTilda is returned as epsilon.
// This is the appropriate variable for most purposes.
tmp<volScalarField> epsilon() const
{
return epsilonTilda_;
}
//- Return the Reynolds stress tensor
tmp<volSymmTensorField> R() const;
//- Return the effective stress tensor including the laminar
stress
tmp<volSymmTensorField> devReff() const;
//- Return the source term for the momentum equation
tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
//- Solve the turbulence equations and correct the turbulence
viscosity
void correct();
//- Read turbulenceProperties dictionary
bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
} // End namespace RASModels
} // End namespace boussinesq //変更箇所
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
#endif
//
************************************************************************* //
３ ＄HOME/OpenFOAM/OpenFOAM-1.5/src/turbulenceModels/RAS/boussinesq/
wmake/file の変更。
###### file #######
/* RAS turbulence models */
RASModel/RASModel.C
RASModel/newRASModel.C
kEpsilon/kEpsilon.C
LaunderSharmaKE/LaunderSharmaKE.C //変更箇所
IMANO Masashi
山田さん
今野です。
この実装が正しいかどうかは、kやepsilonの輸送方程式に浮力生産項をどのように
付加した基礎式を用いるかによるので、実装したいkやepsilonの浮力生産項を書く
か、掲載してあるWebページのURLを書いてもらえますか？
例えば、以下の論文
http://ci.nii.ac.jp/naid/110006973283
では、k方程式へ付加するの浮力生産項は、
Gk=nut_*beta*(g & (fvc::grad(T_)))/Prt
とこの実装のままですが、epsilon方程式へ付加するの浮力生産項は、安定・
不安定により係数を切り変える、いわゆるViollet型
Ref. http://ci.nii.ac.jp/naid/110004180970/
を用いていないので、それを行なっている以下の文が不要になり、
> volScalarField Bepsilon = 0.5*(B + mag(B)); //変更箇所
epsilon方程式(epsEqn)での変更は
> + Cb_ * Bepsilon *epsilonTilda_/k_ //変更箇所
でなく、
> + Cb_ * B *epsilonTilda_/k_ //変更箇所
になると思います。
ただ、最初の文献のような既往のベンチマークテストを追試してみて、コード
のチェックをすることも必要です。がんばってください。
shin
返信遅くなりまして申し訳ありません。
山田です。
ご回答して頂きありがとうございます。
kやepsilonの輸送方程式に浮力生産項を付加した基礎式は
今井様が載せていたURL（http://ci.nii.ac.jp/naid/110006973283 ）にあるものと同様です。
お手数御掛け致しまして、申し訳ありませんでした。
また、Viollet型を用いていないものでしたので、
今井様のご指摘どおり、
プログラムを変更させていただきました。
誠にありがとうございました。
引き続きベンチマークを行い、
コードのチェックを行って参ります。
では失礼致します。
On 12月21日, 午後8:40, IMANO Masashi <im...@arch.t.u-tokyo.ac.jp> wrote:
> 山田さん
>
> 今野です。
>
> この実装が正しいかどうかは、kやepsilonの輸送方程式に浮力生産項をどのように
> 付加した基礎式を用いるかによるので、実装したいkやepsilonの浮力生産項を書く
> か、掲載してあるWebページのURLを書いてもらえますか？
>
> 例えば、以下の論文
>
> http://ci.nii.ac.jp/naid/110006973283
>
> では、k方程式へ付加するの浮力生産項は、
>
> Gk=nut_*beta*(g & (fvc::grad(T_)))/Prt
>
> とこの実装のままですが、epsilon方程式へ付加するの浮力生産項は、安定・
> 不安定により係数を切り変える、いわゆるViollet型
>
> Ref.http://ci.nii.ac.jp/naid/110004180970/