こんにちは。
本記事では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) ); |
本記事は以下の記事を参考にしています。
以下の記事の内容は、数年前に書かれているため年月が経つと今のOpenFOAMのバージョンでのコードと変わってきています。
今のバージョンOpenFOAM-v2112で解説をし直したのがこの記事です。
本記事がOpenFOAMのソースコードの理解の助けになれば幸いです_(._.)_
ラプラシアンソルバの編集
ラプラシアンソルバを自身のディレクトリにコピーします。
1 | cp -R $FOAM_SOLVERS/basic/laplacianFoam/ $WM_PROJECT_USER_DIR/solvers/ |
フォルダ名を変えておきます。
1 | mv laplacianFoam laplacianMultiRegionFoam |
以下のコマンドでフォルダ移動をします。
1 | cd laplacianMultiRegionFoam/ |
ファイル構成は以下のようになっています。
1 2 3 4 5 6 7 8 9 10 11 12 13 | ├── Make │ ├── files │ └── options ├── createFields.H ├── laplacianFoam.C ├── overLaplacianDyMFoam │ ├── Make │ │ ├── files │ │ └── options │ ├── createFields.H │ ├── overLaplacianDyMFoam.C │ └── write.H └── write.H |
「laplacianMultiRegionFoam/Make/files」が実行コマンドを作るファイルですので、以下のように変えておきます。
laplacianMultiRegionFoam/Make/files
1 2 3 | laplacianMultiRegionFoam.C EXE = $(FOAM_USER_APPBIN)/laplacianMultiRegionFoam |
今回新しく作るファイルは「laplacianMultiRegionFoam.C」とします。
環境変数は以下のように設定されています。
1 | $FOAM_USER_APPBIN = /home/ユーザ名/OpenFOAM/kamakiri-v2012/platforms/linux64Gcc63DPInt32Opt/bin |
作成されたlaplacianMultiRegionFoamという実行コマンドは「$FOAM_USER_APPBIN」以下に保存されます。
環境変数は以下のようにすると確かめることができます。
1 2 3 | echo $FOAM_USER_APPBIN #結果 #/home/kamakiri/OpenFOAM/kamakiri-v2012/platforms/linux64Gcc63DPInt32Opt/bin |
では、メインとなるファイル名を変更します。
ここでは元のファイルをコピーして、念のため元のファイルはoldという名前のフォルダに保存しておきます。
※本記事ではlaplacianMultiRegionFoam.Cは変更しないのですが、ソルバの一部を変更した場合は混乱しないようにファイル名も変えておきましょう。
1 2 3 | cp laplacianFoam.C laplacianMultiRegionFoam.C mkdir old mv laplacianFoam.C old |
ここまででファイル構成は以下のようになりました。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | . ├── Make │ ├── files │ └── options ├── createFields.H ├── laplacianMultiRegionFoam.C ├── old │ └── laplacianFoam.C ├── overLaplacianDyMFoam │ ├── Make │ │ ├── files │ │ └── options │ ├── createFields.H │ ├── overLaplacianDyMFoam.C │ └── write.H └── write.H |
ここまでで作業に問題がないかwmakeでコンパイルして確かめてみましょう。
1 | wmake 2>&1 | tee log |
コンパイルをlogファイルに書き出すようにします。
logファイルに出力結果を書き出しておく方が、もしエラーがあった場合にじっくり眺めることができるので便利です。
以下のコマンドで「$FOAM_USER_APPBIN」にlaplacianMultiRegionFoamがあることを確認します。
1 | ls $FOAM_USER_APPBIN |
これで実行コマンドが作成できました。
温度拡散係数の設定
温度拡散係数の設定は「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 | Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading diffusivity DT\n" << endl; 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" |
まず以下のようにvolScalarField型でDTという名前を付けてインスタンス化しています。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | Info<< "Reading diffusivity DT\n" << endl; volScalarField DT ( IOobject ( "DT", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), mesh, dimensionedScalar(dimViscosity, Zero) ); |
volScalarFieldのGeometricFieldを以下のように定義しなおした型なので、volScalarFieldを知りたい場合はGeometricFieldを確認すれば良いです。
1 | typedef GeometricField<scalar, fvPatchField, volMesh> volScalarField; |
何はともあれvolScalarFieldはコントロールボリュームにおけるスカラー値（物理量など）を持った場という認識で良いかと思います。
上のvolScalarFieldをDTという名前にすることで、温度拡散係数DTをコントロールボリューム内の値として保持することになります。
1 | DT = dimensionedScalar("DT", dimViscosity, transportProperties); |
dimensionedクラスを見るとわかりますが、１つ目の引数に変数名、2つ目の引数に次元、3つ目の引数にdictのnameキーワードで保存された値で初期化しています。
1 | dimensioned (const word &name, const dimensionSet &dims, const dictionary &dict) |
つまり、transportPropertiesというdict型から”DT”を探して値で初期化しているということです。
ESI版でもバージョン違いやFoundation版などを使っている場合は以下のようにtransportPropertiesからlookup関数で値と次元を取得している方法も見られますが、いずれの記述も dimensionedScalar DTのようにDTは次元を持った変数だということです。
1 2 3 4 5 | Info<< "Reading diffusivity DT\n" << endl; dimensionedScalar DT ( transportProperties.lookup("DT") ); |
このように次元を持つ変数にすることで異なる次元同士の演算（足し算引き算など）を行うとでエラー表示してくれます。
「IOobject::MUST_READ_IF_MODIFIED」としているのでDTというファイルがあれば読み込みますが、ない場合は読み込まないという設定になっています。
はじめは初期値としてDTは全てのセルで0値ですが、
1 | DT = dimensionedScalar("DT", dimViscosity, transportProperties); |
の記述の後はDTは全てのセルで値が0.1になっています。
しかし、このままだとDTは毎時間固定値を使うことになります。
後で（この記事では扱いませんが）、DTが時間変化に対応する変数にするために下のように変更します。DTが時間変化する場合というのは、例えばDTが温度依存する場合などがあります。温度が時間によって変化するため$DT(T(t))$という具合にDTが変化する場合などに対応します。
1 2 3 4 5 6 7 8 9 10 11 12 13 | volScalarField DT ( IOobject ( "DT", runTime.timeName(), mesh, IOobject::MUST_READ,//READ_IF_PRESENT, IOobject::AUTO_WRITE ), mesh // dimensionedScalar(dimViscosity, Zero) ); |
変更したのは「IOobject::MUST_READ」の部分です。
「runTime.timeName()」が毎時間という意味なので、時間が進むにつれて作成されるフォルダ(0, 10, 20など）の中のDTというファイルを読み込んで、「IOobject::AUTO_WRITE」データを出力していきます。
※ここでは、DTは時間変化する変数ではないので「「IOobject::NO_WRITE」でも良いです。
「IOobject::MUST_READ」としている場合は初期値「0」フォルダにDTという設定ファイルが必要なので用意する必要がありますので、後で用意します。
再度コンパイルします。
1 | wmake 2>&1 | tee log |
エラーがなければOKです。
チュートリアルの編集
ラプラシアンソルバを使った適当なチュートリアルをコピーしてきます。
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とDTを用意します。
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; } } // ************************************************************************* // |
元々のチュートリアルにはありませんでしたので用意しました。
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 | /*--------------------------------*- 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; location "10"; object DT; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 -1 0 0 0 0]; internalField uniform 0.1; boundaryField { inlet { type zeroGradient; } outlet { type zeroGradient; } bottom { type zeroGradient; } top { type zeroGradient; } frontAndBack { type empty; } } // ************************************************************************* // |
以下設定をします。
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; } // ************************************************************************* // |
DTを場所に応じて値を変えたいのでsetFieldを使います。
その前にsetFieldをすると値が更新されてしまうので、DTをDT.orgというオリジナルファイルとして別で保存しておきます。
1 | cp -r 0/DT 0/DT.org |
setFieldの設定ファイルは適当なところからコピーしてきます。
1 | cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/laminar/damBreak/damBreak/system/setFieldsDict system/ |
setFieldは以下のように設定します。
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 | /*--------------------------------*- 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 setFieldsDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // defaultFieldValues ( volScalarFieldValue DT 0.1 ); regions ( boxToCell { box (0.25 0 0) (0.75 1 1); fieldValues ( volScalarFieldValue DT 0.2 ); } ); // ************************************************************************* // |
x方向[0.25,0.75]を温度拡散係数DT=0.2としています。
節点での値を表示すると補完されているのか以下のようになっています。
セルでの値を表示すると以下のようになります。
先ほど作った実行コマンドを実行します。
1 | laplacianMultiRegionFoam > log.laplacianMultiRegionFoam |
温度拡散係数DTは初期状態で値が異なります。
それが時間が経過しても値を変えるような方程式になっていないので、初期状態のままずっと同じ値をとります。
一方温度に関しては以下の方程式を解くため、場所に依存した温度拡散係数に応じて温度分布が時刻変化します。
\frac{\partial T}{\partial t} = \nabla \cdot (D\nabla T)\tag{1}
\end{align*}
初期状態はセル内の温度の値は273Kでしたが、t>0以降で左側の境界条件は273Kのままで、右側は0Kとします。初期状態で右側境界条件をT=0としているため、補完して分布が表示されています。
t=10となるとほぼ定常になっています。
まとめ
温度拡散係数を空間分布させて熱拡散方程式を解くようにカスタマイズしました。
特に意味のないカスタマイズかもしれませんが、OpenFOAMのカスタマイズの練習としては良い例題かなと思います。
参考書
Foundation版で書かれていますが、OpenFOAMの型定義や関数などを調べるのによくまとまった辞書です。
Follow @t_kun_kamakiri