OpenFOAMのソルバー改良（熱拡散方程式の追加）
1.はじめに
1.1 前置き
備忘録です．OpenFOAMのverによって，書き方が異なるため，最新verで確認する．
1.2 環境
- OS: Windows10
- CPU：AMD® Ryzen 9 3900X
- RAM：64 GB
- Ubuntu 2204.2.33.0 (WSL2)
- OpenFOAM: OpenFOAM-v2306
1.3 目的
自作solverの一連の手順を学ぶ．
2. 手順
2.1 通常のsolverをそのままコンパイルしてみる
ソルバーのicoFoamを作業フォルダにコピーする．
ソルバーのディレクトリ：/OpenFOAM-v2306/applications/solvers/incompressible/icoFoam/
mkdir test
mkdir orig
cp -r /OpenFOAM-v2306/applications/solvers/incompressible/icoFoam ./orig
cp -r ./orig/icoFoam ./test
cd test/icoFoam
wmake
ここで，コンパイルが通らないとまずい．OpenFOAMを再インストールすることも視野に入れるべき．
2.2 solver名の変更
フォルダおよびファイルの名前を変更する．ここでは，フォルダ名をmyicoFoam，ファイル名をmyicoFoam.Cに変更した．
mv icoFoam myicoFoam
cd myicoFoam
mv icoFoam.C myicoFoam.C
wmake
ここで，wmakeによるコンパイルはコンパイルは通らないはず．これは，./Make/filesの書き換えができていないためで，以下のように,
./myicoFoamを変更する．
vi ./Make/files
以下のように，./Make/files変更．
myicoFoam.C
EXE = $(FOAM_APPBIN)/myicoFoam
コンパイルの確認．
wmake
フォルダ名，ファイル名が変更されただけで，コードの多くは変更されていないので，コンパイルが通るはず． これで，コンパイルの通し方が分かった．
2.3 TEqn.Hの作成
フォルダ./myicoFoam内に，ファイルTEqn.Hを作成する．
touch TEqn.H
ファイルTEqn.Hに，以下のようなFVM(Finite Volume Method)を書き込む．
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
);
TEqn.solve();
2.4 TEqn.HをcreatFields.Hに追加する
vi creatFields.H
ファイルcreatFields.Hに，以下を書き込む（書き込み場所は怒らくどこでもよい？transportPropatiesよりは下の行？）
Info<< "Reading diffusivity DT\n" << endl;
dimensionedScalar DT
(
"DT", // variable name
dimViscosity, // dimension [m^2/sec]
transportProperties // read-file where "DT" is used (tutorial/icoFoam/cavity/constant/transportProperties)
);
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
2.5 myicoFoam.C内で，TEqn.Hを読み込む
vi myicoFoam.C
ファイルmyicoFoam.Cに，#include "TEqn.H"を追加する．おそらく，PISOループの中？
// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
#include "TEqn.H" // here?
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}
2.6 myicoFoam.Cをコンパイルしてみる
wmake
以下のような，ファイル構成になるはず．
hogehoge:./test/myicoFoam$ tree
.
├── Make
│ ├── files
│ ├── linux64GccDPInt32Opt
│ │ ├── icoFoam.C.dep
│ │ ├── icoFoam.o
│ │ ├── myicoFoam.C.dep
│ │ ├── myicoFoam.o
│ │ ├── options
│ │ ├── sourceFiles
│ │ └── variables
│ └── options
├── TEqn.H
├── createFields.H
└── myicoFoam.C
3 作成したソルバーをテストする
3.1 Tutorialのコピー
tutorialのicoFoamを用いて，ソルバーが機能するか確認する．
tutorial/incompressible/icoFoam/cavity/cavity
cd test
cp -r /OpenFOAM-v2306/tutorial/incompressible/icoFoam/cavity/cavity ../orig
cp -r ../orig/cavity .
mv cavity mycavity
cd mycavity
ディレクトリ構成は以下のようになっているはず．
hogehoge:/test/cavity$ tree
.
├── 0
│ ├── U
│ └── p
├── constant
│ └── transportProperties
└── system
├── PDRblockMeshDict
├── blockMeshDict
├── controlDict
├── decomposeParDict
├── fvSchemes
└── fvSolution
3.2 0/Tの追加
0/pをコピーして，0/Tを作成する
objectやdimensionsなど注意．
cp 0/p 0/T
hogehoge:/test/cavity$ tree
.
├── 0
│ ├── U
│ └── p
│ └── T (added)
├── constant
│ └── transportProperties
└── system
├── PDRblockMeshDict
├── blockMeshDict
├── controlDict
├── decomposeParDict
├── fvSchemes
└── fvSolution
vi 0/T
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2306 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object T; // here
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// from here
dimensions [0 0 0 1 0 0 0];
internalField uniform 273.15;
boundaryField
{
movingWall
{
type fixedValue;
//value $internalField;
value uniform 373.15;
}
fixedWalls
{
// type fixedValue;
// value $internalField;
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// to here
// ************************************************************************* //
3.3 transportPropertiesに"DT"を追加
vi constant/transportProperties
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2306 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu 0.01;
DT [0 2 -1 0 0 0 0] 0.002; // here
// ************************************************************************* //
3.4 fvSchemesにTに関するスキームを追加
vi system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2306 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
div(phi,T) Gauss upwind; // here
}
laplacianSchemes
{
default Gauss linear orthogonal;
laplacian(DT,T) Gauss linear corrected; // here
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //
3.5 fvSolutionにTに関する項目を追加
vi system/fvSolution
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2306 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
// from here
T
{
//solver BICCG;
//preconditioner DILU;
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-7;
relTol 0;
}
// to here
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //
3.6 controlDictのソルバー名を変更
vi system/controlDict
applicationを自作ソルバーのmyicoFoamに変更，また解析時間を5.0，タイムステップを0.001に変更している．
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2306 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application myicoFoam; // here
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 5.0 //0.5;
deltaT 0.001 //0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
3.7 解析実行
blockMesh
myicoFoam
mkdir fig
paraFoam
cd fig
ffmpeg -r 60 -i img.%04d.jpeg -vcodec libx264 -pix_fmt yuv420p -r 24 ani.mp4
3.8 解析結果
左：温度(T)，右：速度(U)
左：温度(T)，右：速度(U)