こんにちは。
本記事ではOpenFOAMのカスタマイズのためのメモを残しておきます。
今回は熱拡散方程式（熱伝導方程式とも呼ぶ）の以下の式を扱います。
\frac{\partial T}{\partial t} = \nabla \cdot (D\nabla T)\tag{1}
\end{align*}
OpenFOAMのソースコードでも以下のようになっており上記のイメージと結びつきやすいでしょう。
1 2 3 4 5 6 | fvScalarMatrix TEqn ( fvm::ddt(T) - fvm::laplacian(DT, T) == fvOptions(T) ); |
今回は熱拡散方程式の温度拡散係数に温度依存性をさせたものにします。
温度拡散係数の温度依存性を追加
DT=A_{1}T^2+A_{2} T +A_{3}\tag{2}
\end{align*}
本記事は以下の記事を参考にしています。
以下の記事の内容は、数年前に書かれているため年月が経つと今のOpenFOAMのバージョンでのコードと変わってきています。
今のバージョンOpenFOAM-v2112で解説をし直したのがこの記事です。
本記事がOpenFOAMのソースコードの理解の助けになれば幸いです_(._.)_
ラプラシアンソルバの編集
まずはlaplacianFoamソルバのベースとなるコードを「laplacianFoamDepT」という名前のフォルダにしてコピーします。
1 | cp -R $FOAM_SOLVERS/basic/laplacianFoam/ $WM_PROJECT_USER_DIR/solvers/laplacianFoamDepT |
「$FOAM_APP」は環境変数で「/opt/OpenFOAM/OpenFOAM-v2012/applications」とパスの設定がされているので、各自の環境に合わせて読みかえてください。
※カスタマイズするものによっては「$WM_PROJECT_USER_DIR = /home/kamakiri/OpenFOAM/kamakiri-v2012」などの環境変数に応じたフォルダでカスタマイズソースを作成する方がOpenFOAMのファイル構成で扱いやすいかもしれません。今回はその辺は気にしなくても良いです。
以下のコマンドで「laplacianFoamDepT」フォルダに移動します。
1 | cd solvers/laplacianFoamDepT |
フォルダとファイルの構成は以下のようになっています。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | . ├── Make │ ├── files │ └── options ├── createFields.H ├── laplacianFoam.C ├── overLaplacianDyMFoam │ ├── Make │ │ ├── files │ │ └── options │ ├── createFields.H │ ├── overLaplacianDyMFoam.C │ └── write.H └── write.H |
以下のコマンドでファイル名を変えておきます。
1 | mv laplacianFoam.C laplacianFoamDep.C |
温度拡散係数の温度依存の追加
温度拡散係数はデフォルトでは温度依存の定数ですが、温度依存する物性値というのはおおくあります。
今回は温度拡散係数の温度依存を以下のように2次の項までの展開した形に変更を行います。
DT=A_{1}T^2+A_{2} T +A_{3}\tag{2}
\end{align*}
「createFields.H」の「DT」を以下のように変更します。
createFields.H
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 | Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); IOdictionary transportProperties ( IOobject ( "transportProperties", runTime.constant(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE ) ); Info<< "Reading diffusivity DT\n" << endl; dimensionedScalar A1(transportProperties.lookup("A1") ); dimensionedScalar A2(transportProperties.lookup("A2") ); dimensionedScalar A3(transportProperties.lookup("A3") ); volScalarField DT ( IOobject ( "DT", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), A1*T*T+A2*T+A3 ); // volScalarField DT // ( // IOobject // ( // "DT", // runTime.timeName(), // mesh, // IOobject::READ_IF_PRESENT, // IOobject::AUTO_WRITE // ), // mesh, // dimensionedScalar(dimViscosity, Zero) // ); // if (!DT.headerOk()) // { // IOdictionary transportProperties // ( // IOobject // ( // "transportProperties", // runTime.constant(), // mesh, // IOobject::MUST_READ_IF_MODIFIED, // IOobject::NO_WRITE // ) // ); // DT = dimensionedScalar("DT", dimViscosity, transportProperties); // } #include "createFvOptions.H" |
まずはtransportPropertiesという名前でIOdictionary型を用意します。
1 2 3 4 5 6 7 8 9 10 11 | IOdictionary transportProperties ( IOobject ( "transportProperties", runTime.constant(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE ) ); |
「runTime.constant()」は文字列のconstantを返します。つまり、OpenFOAMのケースのconstantフォルダに当たります。constantの中の「transportProperties」ファイルを読み込むようにしており、「IOobject::MUST_READ」でファイルの読み込みを、「IOobject::NO_WRITE」ファイルの書き込みを行わないように設定しています。
つぎに、温度拡散係数の各項の係数を定義します。
1 2 3 | dimensionedScalar A1(transportProperties.lookup("A1")); dimensionedScalar A2(transportProperties.lookup("A2")); dimensionedScalar A3(transportProperties.lookup("A3")); |
A1, A2, A3はスカラーを持った量として定義する必要があります
- $A_{1}T^2$
- $A_{2} T$
- $A_{3}$
の各項は温度拡散係数と同じ次元$m^2/s$を持つようにするため、
- $A_{1}$：$m^2/K^2 s$
- $A_{2}$：$m^2/Ks$
- $A_{3}$：$m^2/s$
の単位を持つように設定する必要があります。
ここで単純に値を設定して「1.0T*T + 2.0*T +3.0」などとしても上手くいきそうな感じはしますが、実際は異なる物理量の足し算や引き算は間違いのモノですのでOpenFOAMでは単位がそろったもの同士の足し算引き算を行うようになっています。
この単位の設定は後ほどケースフォルダの「constant/transportProperties」で行います。
「volScalarField」という型で変数名DTとして定義します。
1 2 3 4 5 6 7 8 9 10 11 12 | volScalarField DT ( IOobject ( "DT", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), A1*T*T+A2*T+A3 ); |
ここでは初期値を設定しているだけです。
createFields.HはlaplacianFoamDep.Cで時間発展のループの前にインクルードで読み込んでいるため、 A1*T*T+A2*T+A3と設定しても初期値の温度でDTの値を計算しているにすぎません。
ちなみに、GeometricFieldをtypedefでvolScalarFieldに定義しなおしたものなので詳細を知りたい場合はGeometricFieldを見ることになります。
引数的に以下のコンストラクタが呼び出されているのかと思います。
1 | GeometricField(const IOobject& io, const GeometricField< Type, PatchField, GeoMesh >& gf) |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 | /*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation Copyright (C) 2019 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. Application laplacianFoam Group grpBasicSolvers Description Laplace equation solver for a scalar quantity. \heading Solver details The solver is applicable to, e.g. for thermal diffusion in a solid. The equation is given by: \f[ \ddt{T} = \div \left( D_T \grad T \right) \f] Where: \vartable T | Scalar field which is solved for, e.g. temperature D_T | Diffusion coefficient \endvartable \heading Required fields \plaintable T | Scalar field which is solved for, e.g. temperature \endplaintable \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "fvOptions.H" #include "simpleControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { argList::addNote ( "Laplace equation solver for a scalar quantity." ); #include "postProcess.H" #include "addCheckCaseOptions.H" #include "setRootCaseLists.H" #include "createTime.H" #include "createMesh.H" simpleControl simple(mesh); #include "createFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nCalculating temperature distribution\n" << endl; while (simple.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; while (simple.correctNonOrthogonal()) { DT = A1*T*T + A2*T + A3; //add fvScalarMatrix TEqn ( fvm::ddt(T) - fvm::laplacian(DT, T) == fvOptions(T) ); fvOptions.constrain(TEqn); TEqn.solve(); fvOptions.correct(T); } #include "write.H" runTime.printExecutionTime(Info); } Info<< "End\n" << endl; return 0; } // ************************************************************************* // |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | while (simple.correctNonOrthogonal()) { DT = A1*T*T + A2*T + A3; //add fvScalarMatrix TEqn ( fvm::ddt(T) - fvm::laplacian(DT, T) == fvOptions(T) ); fvOptions.constrain(TEqn); TEqn.solve(); fvOptions.correct(T); } |
では、できたら実行コマンドを作ります。
Make/files
1 2 3 | laplacianFoamDep.C EXE = $(FOAM_USER_APPBIN)/laplacianFoamDep |
Make/options
1 2 3 4 5 6 7 8 | EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ -lfiniteVolume \ -lfvOptions \ -lmeshTools |
以下のコマンドでコンパイルを行います。
1 2 | wclean wmake 2>&1 | tee log |
「$FOAM_USER_APPBIN」にlaplacianFoamDepができていたら成功です。
チュートリアルの編集
ケースファイルは前回の記事同様です。
以下に設定ファイルを載せておきます。
ラプラシアンソルバを使った適当なチュートリアルをコピーしてきます。
1 | cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/basic/laplacianFoam/flange myCase |
このチュートリアルはblockMeshが無いので適当に持ってきます。
1 | cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/basic/scalarTransportFoam/pitzDaily/system/blockMeshDict ../myCase/system/ |
blockMeshDictは適当に設定します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object blockMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // scale 1.0; vertices ( ( 0.0 0.0 0.0) ( 1.0 0.0 0.0) ( 1.0 0.0 1.0) ( 0.0 0.0 1.0) ( 0.0 1.0 0.0) ( 1.0 1.0 0.0) ( 1.0 1.0 1.0) ( 0.0 1.0 1.0) ); blocks ( hex (0 1 5 4 3 2 6 7) (5 1 1) simpleGrading (1 1 1) ); edges ( ); patches ( patch inlet ( (0 4 7 3) ) patch outlet ( (1 5 6 2) ) patch bottom ( (0 1 5 4) ) patch top ( (3 2 6 7) ) empty frontAndBack ( (0 1 2 3) (4 5 6 7) ) ); mergePatchPairs ( ); // ************************************************************************* // |
では、blockMeshを実行します。
1 | blockMesh |
ParaViewでモデルを読み込んで確認しましょう。
今回はわかりやすく1×1×1(m3)でx方向に5分割したモデルにしました。
初期条件として「0」フォルダにはTを用意します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object T; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 0 0 1 0 0 0]; internalField uniform 273; boundaryField { inlet { type fixedValue; value uniform 1; } outlet { type fixedValue; value uniform 0; } top { type slip; } bottom { type slip; } frontAndBack { type empty; } } // ************************************************************************* // |
※前回の記事でDT設定ファイルが必要でしたがcreateField.HでDTの定義は以下のようになっているため、0/DTは必要ありません。
1 2 3 4 5 6 7 8 9 10 11 12 13 | Info<< "Reading diffusivity DT\n" << endl; volScalarField DT ( IOobject ( "DT", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), A1*T*T+A2*T+A3 ); |
注目は「IOobject::NO_READ」の部分です。
「runTime.timeName()」は毎時間出力したフォルダを指しますが、「IOobject::NO_READ」となっていることで、そこにDTが無くてもエラーとなりません。「A1*T*T+A2*T+A3」で次元を持った初期値を設定しているイメージです。
温度拡散係数の係数をconstant/transportPropertiesで行います。
constant/transportProperties
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object transportProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // DT 0.1; A1 A1 [ 0 2 -1 -2 0 0 0 ] 1e-8; A2 A2 [ 0 2 -1 -1 0 0 0 ] 1e-5; A3 A3 [ 0 2 -1 0 0 0 0 ] 0.01; // ************************************************************************* // |
DTは値だけですが、A1,A2,A3は次元と値を設定しています。
少し気持ちの悪い感じがしますので、A1,A2,A3も値だけを設定すれば計算ができるように後ほど変更します。実際は、A1,A2,A3の次元を頭に入れて計算するなんてユーザー目線からすると面倒なのでプログラムの中で自動的に次元も考慮するようにしたいものですよね。
後ほど変更するとして、今はこのまま進めていきます。
またプログラムが正常に動作しているかを確認するために、A1,A2,A3の値はわかりやすい数値にしました。
以下設定をします。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application laplacianFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 10; deltaT 0.1; writeControl runTime; writeInterval 1; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; // ************************************************************************* // |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default Euler; } gradSchemes { default none; grad(T) Gauss linear; } divSchemes { default none; } laplacianSchemes { default none; laplacian(DT,T) Gauss linear corrected; //laplacian(DT,T) Gauss harmonic corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } // ************************************************************************* // |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { T { solver PCG; preconditioner DIC; tolerance 1e-06; relTol 0; } } SIMPLE { nNonOrthogonalCorrectors 2; } // ************************************************************************* // |
以下のコマンドで計算を実行します。
1 | laplacianFoamDep > log.laplacianFoamDep & |
ParaViewで結果を確認しましょう。
特に温度分布から温度拡散係数を手計算で求めた値と一致しているかを確認します。
問題はなさそうです。
コードを編集する
createFields.Hが以下のように書いていたせいでconstant/transportPropertiesにA1, A2, A3の次元を書く必要がありました。
createFields.H(1部抜粋）
1 2 3 | dimensionedScalar A1(transportProperties.lookup("A1") ); dimensionedScalar A2(transportProperties.lookup("A2") ); dimensionedScalar A3(transportProperties.lookup("A3") ); |
constant/transportProperties
1 2 3 4 | DT 0.1; A1 A1 [ 0 2 -1 -2 0 0 0 ] 1e-8; A2 A2 [ 0 2 -1 -1 0 0 0 ] 1e-5; A3 A3 [ 0 2 -1 0 0 0 0 ] 0.01; |
これではDTは値だけ書けばよく、次元は意識しなくても良いのにA1,A2,A3は次元と値も書く必要があり少々面倒に感じます。
できれば次元を意識することなくA1,A2,A3の値のみを設定したいです。
そのためには、以下のようにcreateFields.HでA1,A2,A3の次元も一緒に定義しておけば良いです。
createFields.H(1部抜粋）
1 2 3 | dimensionedScalar A1("A1", dimensionSet(0, 2, -1, -2, 0, 0, 0), transportProperties.lookup("A1")); dimensionedScalar A2("A2", dimensionSet(0, 2, -1, -1, 0, 0, 0), transportProperties.lookup("A2")); dimensionedScalar A3("A3", dimensionSet(0, 2, -1, 0, 0, 0, 0), transportProperties.lookup("A3")); |
※dimensionSet(kg, m, s, K, mol, A, Cd)のように関数の引数に次数を記述します。
※ちなみにdimensionedScalarは以下のようにdimensioned型にテンプレートでscalar(実数）を定義しなおしています。
1 | typedef dimensioned<scalar> dimensionedScalar; |
以下のようにコンストラクタで引数を設定します。
1 | dimensioned (const word &name, const dimensionSet &dims, const Type &val) |
変更したら再度wmakeでコンパイルしなおします。
そうするとconstant/transportProperties以下のようにシンプルな記述で済みます。
constant/transportProperties
1 2 3 4 | DT 0.1; A1 1e-8; A2 1e-5; A3 0.01; |
以上のように設定できたら再度計算をして結果を確認してみましょう。
先ほどと同じ結果だったらOKです。
補足（コンパイル時のwarning)
※コンパイルすると以下のようなwarningが出ます。
1 | ‘Foam::dimensioned<Type>::dimensioned(const Foam::word&, const Foam::dimensionSet&, Foam::Istream&) [with Type = double]’ is deprecated: Since 2018-11 |
どうやら
1 | dimensionedScalar A1("A1", dimensionSet(0, 2, -1, -2, 0, 0, 0), transportProperties.lookup("A1")); |
のような書き方は非推奨となっているようです。
以下のような書き方で良いようです。
1 2 3 | dimensionedScalar A1("A1", dimensionSet(0, 2, -1, -2, 0, 0, 0), transportProperties); dimensionedScalar A2("A2", dimensionSet(0, 2, -1, -1, 0, 0, 0), transportProperties); dimensionedScalar A3("A3", dimensionSet(0, 2, -1, 0, 0, 0, 0), transportProperties); |
以上のように書き換えて再度コンパイルするとwarningは消えます。
解析結果も先ほどの結果と同じです。
まとめ
温度拡散係数に温度依存性を持たせた熱拡散方程式を解くようにカスタマイズしました。
特に温度依存の係数は適当に設定したので意味のないカスタマイズかもしれませんが、OpenFOAMのカスタマイズの練習としては良い例題かなと思います。
参考書
Foundation版で書かれていますが、OpenFOAMの型定義や関数などを調べるのによくまとまった辞書です。
Follow @t_kun_kamakiri