OpenFOAM v2006で低Re数型非線形RANSモデルと壁関数を実装する
計算に使ったファイルはGithubにアップしてある
≫mtkbirdman.com/OpenFOAM/
はじめに
参考にしたサイトは以下の通り
≫OpenFOAMカスタマイズ入門：レイノルズ平均型乱流モデルのカスタマイズに挑戦 - Qiita
≫Turbulence Models in OpenFOAM
今回実装する低Re数型非線形RANSモデルはこれ
≫Ken-ichi ABE, A Hybrid LES/RANS Approach Using an Anisotropy-Resolving Algebraic Turbulence Model, International Journal of Heat and Fluid Flow, Vol. 26, pp. 204-222, 2005.
≫Ken-ichi ABE, Yong-Jun JANG and Michael A. LESCHZINER, An Investigation of Wall-Anisotropy Expressions and Length-Scale Equations forNon-Linear Eddy-Viscosity Models, International Journal of Heat and Fluid Flow, Vol. 24, pp. 181-198, 2003.
2003年の論文のモデルを改良した2005年の論文のモデル(AJL2005)を実装する．記述は2003年の論文の方が詳しい
乱流モデルAJL2005についての説明は面倒なので紙面の都合上割愛させていただく
C++については次のサイトを見て完全に理解した
≫C++入門
OpenFOAMのテンソル計算については次のサイトが詳しい
≫Tensor Operations in OpenFOAM _ CFD WITH A MISSION
OpenFOAMにおけるRANSについてのドキュメントはこちら
≫OpenFOAM User Guide Reynolds Averaged Simulation (RAS)
乱流モデル実装までの手順
乱流モデルの実装は以下の手順で行う
- 実装したい乱流モデルを見つける
- OpenFOAMに実装されている乱流モデルの中から実装したいモデルに近いものを見つける (このモデルをベースに実装する)
- ベースとなる乱流モデルのファイルをコピーして乱流モデルの名前を書き換える
- 名前を書き換えた乱流モデルをコンパイルして計算してみる
- 乱流モデルのファイルの中身を書き換える
- 書き換えた乱流モデルをコンパイルして計算する
手順1はすでに完了しているので，手順2から説明していく
ベースとなる乱流モデルの決定
AJL2005は低Re数型非線形RANSモデルである
OpenFOAMで実装されているNon-linear eddy viscosity modelsには
の2つがあるが，低Re数型であるのはLienCubicKEだけなのでこれをベースとする
ちなみになんとLienCubicKEの参考論文はインターネット上では手に入らない (筆者調べ2020年9月時点)
どうしても読みたい人は近くの図書館で探してみよう
≫CiNii 図書 - Engineering turbulence modelling and experiments 3 _ proceedings of the Third International Symposium on Engineering Turbulence Modelling and Measurements, Heraklion-Crete, Greece, 27-29.html
ベースとなる乱流モデルのコピー・名前の書き換え
ベースとなる乱流モデルをコピーしてくる
LienCubicKEに関するファイルは$WM_PROJECT_DIR/src/TurbulenceModels/incompressible/turbulentTransportModels/RASというディレクトリにあるので，これを$WM_PROJECT_USER_DIR/src/myTurbulenceModelsにコピーする
myTurbulenceModelsにあたるディレクトリ名はわかりやすければ何でもいい
>> mkdir -p $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> cp -r $WM_PROJECT_DIR/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKE $WM_PROJECT_USER_DIR/src/myTurbulenceModelsちなみに$WM_PROJECT_DIRや$WM_PROJECT_USER_DIRは$WM_PROJECT_DIR/etc/bashrcで定義されている変数で，echoを使って中身を確認できる
>> echo $WM_PROJECT_DIR
/opt/OpenFOAM/OpenFOAM-v2006
>> echo $WM_PROJECT_USER_DIR
/home/mtkbirdman/OpenFOAM/mtkbirdman-v2006コピーしてきたすべてのファイルでLienCubicKEをAJL2005に書き換える
>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> mv LienCubicKE AJL2005
>> cd AJL2005
>> mv LienCubicKE.C AJL2005.C
>> mv LienCubicKE.H AJL2005.H
>> sed -i -e "s/LienCubicKE/AJL2005/g" AJL2005.*乱流モデルのコンパイル(1回目)
名前を書き換えたAJL2005（実際はLienCubicKE）をコンパイルするために$WM_PROJECT_USER_DIR/src/myTurbulenceModelsに2つのファイルMake/files，Make/optionsを作成する
>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> mkdir Make
>> cd Make
>> touch files
>> touch optionsそれぞれのファイルに以下を書きこむ
Make/files
AJL2005/AJL2005.C
LIB = $(FOAM_USER_LIBBIN)/libmyTurbulenceModelsMake/options
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
LIB_LIBS = \
-lincompressibleTurbulenceModels \
-lturbulenceModels \
-lfvOptions \
-lfiniteVolume \
-lmeshToolsここまででディレクトリ構成は以下のようになる
myTurbulenceModels
├── AJL2005
│ ├── AJL2005.C
│ └── AJL2005.H
└── Make
├── files
└── optionsAJL2005をコンパイルする
>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> wmake libso
wmake libso (myTurbulenceModels)
Make/options:9:12: warning: backslash-newline at end of file
LIB_LIBS = \
ln: ./lnInclude
Making dependency list for source file AJL2005.C
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas -O3 -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude -fPIC -c AJL2005/AJL2005.C -o Make/linux64Gcc63DPInt32Opt/AJL2005/AJL2005.o
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas -O3 -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude -fPIC -shared -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64Gcc63DPInt32Opt/AJL2005/AJL2005.o -L/opt/OpenFOAM/OpenFOAM-v2006/platforms/linux64Gcc63DPInt32Opt/lib \
-lincompressibleTurbulenceModels -lturbulenceModels -lfvOptions -lfiniteVolume -lmeshTools -o /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/platforms/linux64Gcc63DPInt32Opt/lib/libmyTurbulenceModels.soプログラムに問題があればここで大量のエラーを吐き出してくれる
※ちなみに，Make/filesとMake/optionsの編集をVS codeでやりたいときはWSLのhomeディレクトリをエクスプローラーで開く必要がある．そのときのhomeディレクトリのパスを以下に示す
C:\Users\<ユーザー名>\AppData\Local\Packages\CanonicalGroupLimited.Ubuntu20.04onWindows_79rhkp1fndgsc\LocalState\rootfs\home\<WSLのユーザー名>テストケースの準備
テストケースとして/tutorials/incompressible/pimpleFoam/LES/channel395を$WM_PROJECT_USER_DIRにコピーしてくる
>> cp -r $WM_PROJECT_DIR/tutorials/incompressible/pimpleFoam/LES/channel395 $WM_PROJECT_USER_DIR
もともとあった0を削除して0/nuTildaを0/epsilonに変更し，0.origをコピーする
>> cd $WM_PROJECT_USER_DIR/channel395
>> rm -r 0
>> mv 0.orig/nuTilda 0.orig/epsilon
>> cp -r 0.orig 0いくつかのファイルに修正を加える
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -3 0 0 0 0];
internalField uniform 1;
boundaryField
{
bottomWall
{
type epsilonWallFunction;
value $internalField;
}
topWall
{
type epsilonWallFunction;
value $internalField;
}
sides1_half0
{
type cyclic;
}
sides2_half0
{
type cyclic;
}
inout1_half0
{
type cyclic;
}
inout2_half0
{
type cyclic;
}
sides2_half1
{
type cyclic;
}
sides1_half1
{
type cyclic;
}
inout1_half1
{
type cyclic;
}
inout2_half1
{
type cyclic;
}
}
// ************************************************************************* //
内部領域の値は1を指定し，壁関数はepsilonWallFunctionを用いる
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "1";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 1;
boundaryField
{
bottomWall
{
type fixedValue;
value uniform 1e-32;
}
topWall
{
type fixedValue;
value uniform 1e-32;
}
sides1_half0
{
type cyclic;
}
sides2_half0
{
type cyclic;
}
inout1_half0
{
type cyclic;
}
inout2_half0
{
type cyclic;
}
sides2_half1
{
type cyclic;
}
sides1_half1
{
type cyclic;
}
inout1_half1
{
type cyclic;
}
inout2_half1
{
type cyclic;
}
}
// ************************************************************************* //
内部領域の値は1を指定し，壁面上ではできるだけ小さい値(1e-32)を指定する
本来は滑りなし条件のため壁面上でk=0であるが，0を指定するとプログラム上のどこかで0割りが発生してしまう
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "1";
object nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
bottomWall
{
type nutLowReWallFunction;
value uniform 0;
}
topWall
{
type nutLowReWallFunction;
value uniform 0;
}
sides1_half0
{
type cyclic;
}
sides2_half0
{
type cyclic;
}
inout1_half0
{
type cyclic;
}
inout2_half0
{
type cyclic;
}
sides2_half1
{
type cyclic;
}
sides1_half1
{
type cyclic;
}
inout1_half1
{
type cyclic;
}
inout2_half1
{
type cyclic;
}
}
// ************************************************************************* //
壁関数はnutLowReWallFunctionを使う
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType RAS;
RAS
{
RASModel AJL2005;
turbulence on;
printCoeffs on;
}
// ************************************************************************* //
RASモデルのAJL2005を指定する
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application simpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 5000;
deltaT 1;
writeControl timeStep;
writeInterval 5000;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
functions
{
#includeFunc residuals
#includeFunc "sample"
#includeFunc "turbulenceFields"
#includeFunc "wallShearStress"
}
libs
(
"libmyTurbulenceModels.so"
);
// ************************************************************************* //ALJ2005のライブラリを読み込むために，libs("libmyTurbulenceModels.so");を加えている
functionsについてはGithubのコードを参照
≫mtkbirdman.com/OpenFOAM/channel395/
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) bounded Gauss linearUpwind grad(U);
div(phi,k) bounded Gauss limitedLinear 1;
div(phi,epsilon) bounded Gauss limitedLinear 1;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
div(nonlinearStress) Gauss linear;
}
laplacianSchemes
{
default Gauss linear uncorrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default uncorrected;
}
wallDist
{
method meshWave;
}
// ************************************************************************* //div(nonlinearStress) Gauss linear;を加えるなど
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-10;
relTol 0.1;
smoother DICGaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
cacheAgglomeration on;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 10;
mergeLevels 1;
}
"(U|k|epsilon|R|nuTilda)"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-10;
relTol 0.1;
}
omega
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-10;
relTol 0.1;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
consistent yes;
pRefCell 1001;
pRefValue 0;
}
relaxationFactors
{
equations
{
U 0.9;
k 0.8;
epsilon 0.7;
nuTilda 0.8;
omega 0.8;
}
}
cache
{
grad(U);
}
// ************************************************************************* //PIMPLE→SIMPLE（あんまりわかっていない）
Allrun
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
runApplication blockMesh
runApplication decomposePar -cellDist
runParallel $(getApplication)
runApplication reconstructPar
runApplication postChannel
#------------------------------------------------------------------------------
計算の実行(1回目)
計算を実行して，コンパイルがうまくできていることを確認する
ちなみに計算はこんな感じ
結果が正確かどうかは置いといて，とりあえず正常に計算は回っているので良しとする
プログラムの修正
プログラムの全体はGithubにあるコードを確認してもらうとして，ここでは重要な部分のみ説明する
≫mtkbirdman.com/OpenFOAM/src/myTurbulenceModels/AJL2005/
AJL2005.C
メインの部分がここ
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<volScalarField> AJL2005::Rt() const
{
return sqr(k_)/(nu()*epsilon_);
}
tmp<volScalarField> AJL2005::fw(const dimensioned<scalar> xi) const
{
return exp(-sqr((pow(nu()*epsilon_,0.25)*y_/nu())/xi));
}
tmp<volScalarField> AJL2005::fMu(const volScalarField& Rt,const volScalarField& fw) const
{
return (scalar(1) + (35.0/pow(Rt,0.75))*exp(-pow(Rt/30.0,0.75)))
*(scalar(1) - fw);
}
tmp<volScalarField> AJL2005::f2(const volScalarField& Rt,const volScalarField& fw) const //feps
{
return (scalar(1) - 0.3*exp(-sqr(Rt/6.5)))
*(scalar(1) - fw);
}
void AJL2005::correctNut()
{
correctNonlinearStress(fvc::grad(U_));
}
void AJL2005::correctNonlinearStress(const volTensorField& gradU)
{
nut_ = Cmu_*fMu(Rt(),fw(26.0))*sqr(k_)/epsilon_;
nut_.correctBoundaryConditions();
volSymmTensorField S(symm(gradU));
volTensorField W(-skew(gradU));
volScalarField tau(nut_/k_);
volScalarField fB(scalar(1) + Ceta_*(CD_*tau)*(mag(W) - mag(S)));
volScalarField CB(1.0/(scalar(1) + (22.0/3.0)*sqr(CD_*tau)*magSqr(W) + (2.0/3.0)*sqr(CD_*tau)*(magSqr(W) - magSqr(S))*fB));
volScalarField fr1((magSqr(W) - magSqr(S))/max(magSqr(W) + magSqr(S),dimensionedScalar("minWS",dimensionSet(0,0,-2,0,0,0,0),SMALL)));
volScalarField fr2(magSqr(S)/max(magSqr(W) + magSqr(S),dimensionedScalar("minWS",dimensionSet(0,0,-2,0,0,0,0),SMALL)));
volScalarField fs1(fr1*fr2*(0.15*Ceta_)*sqr(CD_*tau)*(magSqr(W) - magSqr(S)));
volScalarField fs2(-fr1*fr2*(scalar(1)+(0.07*Ceta_)*(CD_*tau)*(mag(W) - mag(S))));
volScalarField nStar((nu()*epsilon_,0.25)*y_/nu());
volVectorField N(fvc::grad(nStar));
volVectorField d(N/mag(N));
volSymmTensorField dd(symm(d*d));
volScalarField taud1((scalar(1) - fw(15.0))*(k_/epsilon_) + fw(15.0)*(1.0)*sqrt(nu()/epsilon_));
volSymmTensorField bw1(
fw(26.0)*(
- (1.0)*(1.0/2.0)*dev(dd)
+ (scalar(1) - sqr(fr1))*sqr(taud1)*(
((0.25*0.5)/(scalar(1) + 0.5*sqr(taud1)*mag(S)*mag(W)))*twoSymm(S&W)
+ ((1.5*0.5)/(scalar(1) + 0.5*sqr(taud1)*magSqr(S)))*dev(innerSqr(S))
)
)
);
volScalarField taud2((scalar(1)-fw(15.0))*(k_/epsilon_) + fw(15.0)*(3.0)*sqrt(nu()/epsilon_));
volSymmTensorField bw2(
fw(26.0)*(
- (0.0)*(1.0/2.0)*dev(dd)
+ (scalar(1) - sqr(fr1))*sqr(taud2)*(
(((13.0/30.0)*1.0)/(scalar(1) + 1.0*sqr(taud2)*mag(S)*mag(W)))*twoSymm(S&W)
+ ((0.6*1.0)/(scalar(1) + 1.0*sqr(taud2)*magSqr(S)))*dev(innerSqr(S))
)
)
);
volScalarField ftau(exp(-3.0*(scalar(1) - fw(26.0))*(k_/epsilon_)*mag(S)));
volSymmTensorField bw(ftau*bw1 + (scalar(1) - ftau)*bw2);
nonlinearStress_ =
2.0*nut_*S // Cancel linear term
+2.0*(k_/CD_)*(
- CB*(tau*CD_)*S //b1
+ (scalar(1) - fw(26.0))*(
2.0*CB*sqr(tau*CD_)*(twoSymm(S&W) + dev(innerSqr(S))) //b2
- CB*fs1*(tau*CD_)*S + 2.0*fs2*sqr(tau*CD_)*dev(innerSqr(S)) //bs
)
+ CD_*bw //bw
);
}
注意すべき点は以下の通り
- \({\bf \Omega}={\bf -}{\sf skew(gradU)}\)
nonlinearStress_はvolSymmTensorFieldnonlinearStress_はレイノルズ応力の非線形項のみ
\({\bf \Omega}={\bf -}{\sf skew(gradU)}\)についてはTensor Operations in OpenFOAMに書いてある
レイノルズ応力テンソルRはeddyViscosity.Cで次のように定義されている
template<class BasicTurbulenceModel>
Foam::tmp<Foam::volSymmTensorField>
Foam::eddyViscosity<BasicTurbulenceModel>::R() const
{
tmp<volScalarField> tk(k());
// Get list of patchField type names from k
wordList patchFieldTypes(tk().boundaryField().types());
(中略)
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
IOobject::groupName("R", this->alphaRhoPhi_.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
((2.0/3.0)*I)*tk() - (nut_)*dev(twoSymm(fvc::grad(this->U_))),
patchFieldTypes
)
);
}この時点でRは次の形になっている
\begin{eqnarray}
{\bf R_{ij}}={\overline {u'_{i} u'_{j}}}=\frac{2}{3}k\delta_{ij}-\nu_{t}(\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial U_{j}}{\partial x_{i}})
\end{eqnarray}
そして，nonlinearStress_はnonlinearEddyViscosity.Cで次のように定義され，レイノルズ応力テンソルRに加えられている
nonlinearStress_
(
IOobject
(
IOobject::groupName("nonlinearStress", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_
),
this->mesh_,
dimensionedSymmTensor(sqr(dimVelocity), Zero)
)template<class BasicTurbulenceModel>
Foam::tmp<Foam::volSymmTensorField>
Foam::nonlinearEddyViscosity<BasicTurbulenceModel>::R() const
{
tmp<volSymmTensorField> tR
(
eddyViscosity<BasicTurbulenceModel>::R()
);
tR.ref() += nonlinearStress_;
return tR;
}このようにnonlinearStress_はvolSymmTensorFieldでレイノルズ応力の非線形項のみなので，nonlinearStress_ = ...の右辺のテンソルはすべてvolSymmTensorFieldである必要があるし，AJL2005では線形項をキャンセルする必要がある
輸送方程式の部分はほとんど変更していない
void AJL2005::correct()
{
if (!turbulence_)
{
return;
}
nonlinearEddyViscosity<incompressible::RASModel>::correct();
tmp<volTensorField> tgradU = fvc::grad(U_);
const volTensorField& gradU = tgradU();
volScalarField G
(
GName(),
(nut_*twoSymm(gradU) - nonlinearStress_) && gradU
);
// Update epsilon and G at the wall
epsilon_.boundaryFieldRef().updateCoeffs();
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
Ceps1_*G*epsilon_/k_
- fvm::Sp(Ceps2_*f2(Rt(),fw(3.3))*epsilon_/k_, epsilon_)
);
epsEqn.ref().relax();
epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
solve(epsEqn);
bound(epsilon_, epsilonMin_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::laplacian(DkEff(), k_)
==
G
- fvm::Sp(epsilon_/k_, k_)
);
kEqn.ref().relax();
solve(kEqn);
bound(k_, kMin_);
// Re-calculate viscosity and non-linear stress
correctNonlinearStress(gradU);
}AJL2005.H
大事なのはここらへん
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff() const
{
return tmp<volScalarField>
(
new volScalarField(
"DkEff",
nut_/(sigmak_/(scalar(1)+5.0*fw(5.0))) + nu()
)
);
}
//- Return the effective diffusivity for epsilon
tmp<volScalarField> DepsilonEff() const
{
return tmp<volScalarField>
(
new volScalarField(
"DepsilonEff",
nut_/(sigmaEps_/(scalar(1)+5.0*fw(5.0))) + nu()
)
);
}ヘッダーファイルの中でもfwは使えた
低Re数型RANSモデル用の壁関数
OpenFOAMに実装されているepsilonの壁関数はepsilonWallFunctionしかないが，なんか無駄な機能がごちゃごちゃしているのでシンプルな低Re数型RANSモデル用の壁関数も実装しておく
>> mkdir -p $WM_PROJECT_USER_DIR/src/myWallFunctions
>> cp -r $WM_PROJECT_DIR/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction $WM_PROJECT_USER_DIR/src/myWallFunctions>> cd $WM_PROJECT_USER_DIR/src/myWallFunctions
>> mv epsilonWallFunction epsilonLowReWallFunction>> cd $WM_PROJECT_USER_DIR/src/myWallFunctions/epsilonLowReWallFunction
>> mv epsilonWallFunctionFvPatchScalarField.C epsilonLowReWallFunctionFvPatchScalarField.C
>> mv epsilonWallFunctionFvPatchScalarField.H epsilonLowReWallFunctionFvPatchScalarField.H
>> sed -i s/epsilonWallFunctionFvPatchScalarField/epsilonLowReWallFunctionFvPatchScalarField/g epsilonLowReWallFunctionFvPatchScalarField.*
>> mkdir $WM_PROJECT_USER_DIR/src/myWallFunctions/Make
>> cd $WM_PROJECT_USER_DIR/src/myWallFunctions/Make
>> touch files
>> touch optionsepsilonLowReWallFunctionFvPatchScalarField.C
大事な部分だけ
void Foam::epsilonLowReWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbModel,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& epsilon0
)
{
const label patchi = patch.index();
const scalarField& y = turbModel.y()[patchi];
const tmp<volScalarField> tk = turbModel.k();
const volScalarField& k = tk();
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const tmp<scalarField> tnutw = turbModel.nut(patchi);
const scalarField& nutw = tnutw();
const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
const scalarField magGradUw(mag(Uw.snGrad()));
forAll(nutw, facei)
{
const label celli = patch.faceCells()[facei];
const scalar w = cornerWeights[facei];
epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
}
}epsilonLowReWallFunctionFvPatchScalarField.H
Githubのコードを参照
≫mtkbirdman.com/OpenFOAM/src/myWallFunctions/epsilonLowReWallFunction/
Make/files
epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
LIB = $(FOAM_USER_LIBBIN)/libmyWallFunctionsMake/options
もしかしたら必要ないものが混じっているかも
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
LIB_LIBS = \
-lincompressibleTurbulenceModels \
-lturbulenceModels \
-lfvOptions \
-lfiniteVolume \
-lmeshToolsコンパイル・計算の実行(2回目)
コンパイル
AJL2005とepsilonLowReWallFunctionをそれぞれコンパイルする
>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> wclean libso
>> wmake libso
wmake libso (myTurbulenceModels)
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas -O3 -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude -fPIC -c AJL2005/AJL2005.C -o Make/linux64Gcc63DPInt32Opt/AJL2005/AJL2005.o
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas -O3 -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude -fPIC -shared -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64Gcc63DPInt32Opt/AJL2005/AJL2005.o -L/opt/OpenFOAM/OpenFOAM-v2006/platforms/linux64Gcc63DPInt32Opt/lib \
-lincompressibleTurbulenceModels -lturbulenceModels -lfvOptions -lfiniteVolume -lmeshTools -o /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/platforms/linux64Gcc63DPInt32Opt/lib/libmyTurbulenceModels.so>> cd $WM_PROJECT_USER_DIR/src/myWallFunctions
>> wmake libso
wmake libso (myWallFunctions)
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas -O3 -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude -fPIC -c epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C -o Make/linux64Gcc63DPInt32Opt/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.o
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas -O3 -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude -fPIC -shared -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64Gcc63DPInt32Opt/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.o -L/opt/OpenFOAM/OpenFOAM-v2006/platforms/linux64Gcc63DPInt32Opt/lib \
-lincompressibleTurbulenceModels -lturbulenceModels -lfvOptions -lfiniteVolume -lmeshTools -o /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/platforms/linux64Gcc63DPInt32Opt/lib/libmyWallFunctions.soいくつかのファイルを変更する
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -3 0 0 0 0];
internalField uniform 1;
boundaryField
{
bottomWall
{
type epsilonLowReWallFunction;
value $internalField;
}
topWall
{
type epsilonLowReWallFunction;
value $internalField;
}
sides1_half0
{
type cyclic;
}
sides2_half0
{
type cyclic;
}
inout1_half0
{
type cyclic;
}
inout2_half0
{
type cyclic;
}
sides2_half1
{
type cyclic;
}
sides1_half1
{
type cyclic;
}
inout1_half1
{
type cyclic;
}
inout2_half1
{
type cyclic;
}
}
// ************************************************************************* //
壁関数をepsilonLowReWallFunctionに変更する
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application simpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 5000;
deltaT 1;
writeControl timeStep;
writeInterval 5000;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
functions
{
#includeFunc residuals
#includeFunc "sample"
#includeFunc "turbulenceFields"
#includeFunc "wallShearStress"
}
libs
(
"libmyTurbulenceModels.so"
"libmyWallFunctions.so"
);
// ************************************************************************* //
libs("libmyWallFunctions.so");を追加する
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scale 1;
vertices
(
(0 0 0)
(4 0 0)
(0 1 0)
(4 1 0)
(0 2 0)
(4 2 0)
(0 0 2)
(4 0 2)
(0 1 2)
(4 1 2)
(0 2 2)
(4 2 2)
);
blocks
(
hex (0 1 3 2 6 7 9 8) (40 25 30) simpleGrading (1 25 1)
hex (2 3 5 4 8 9 11 10) (40 25 30) simpleGrading (1 -25 1)
);
edges
(
);
boundary
(
bottomWall
{
type wall;
faces ((0 1 7 6));
}
topWall
{
type wall;
faces ((4 10 11 5));
}
sides1_half0
{
type cyclic;
neighbourPatch sides1_half1;
faces ((0 2 3 1));
}
sides1_half1
{
type cyclic;
neighbourPatch sides1_half0;
faces ((6 7 9 8));
}
sides2_half0
{
type cyclic;
neighbourPatch sides2_half1;
faces ((2 4 5 3));
}
sides2_half1
{
type cyclic;
neighbourPatch sides2_half0;
faces ((8 9 11 10));
}
inout1_half0
{
type cyclic;
neighbourPatch inout1_half1;
faces ((1 3 9 7));
}
inout1_half1
{
type cyclic;
neighbourPatch inout1_half0;
faces ((0 6 8 2));
}
inout2_half0
{
type cyclic;
neighbourPatch inout2_half1;
faces ((3 5 11 9));
}
inout2_half1
{
type cyclic;
neighbourPatch inout2_half0;
faces ((2 8 10 4));
}
);
mergePatchPairs
(
);
// ************************************************************************* //
メッシュのGradingを変えてy+の値を小さくする
計算の実行
計算を実行する
結果はこんな感じ
きちんと論文に書いてあるのと同じような分布を得ることができた
まとめ
今回は乱流モデルの実装に挑戦してみた
低Re数型RANSモデルの実装についての資料はあまりなかったのでどこかの誰かの参考になればうれしい
この記事でグラフを書くのに使ったExcelファイルはこちら
そのほかの乱流モデルについてはこちら
コメント