はじめに（着手日：210916, HP更新日：221011）
OpenFOAMのver.が古いので，最新のに対応しているかわからない．自作solverの制作 (interPhaseChangeFoamの利用)
公開solverが古いverのOpenFOAMで制作されており、(自分の)現在の環境では利用できなさそうであった。 そこで自作solverを制作する。 solverのオリジナルは"interPhaseChangeFoam"を利用する。 (解析対象はチュートリアルのcavitation bullet )。 解析環境OpenFOAM 8 Fundation版 OS：Ubuntu 20.04.3 CPU：AMD® Ryzen 9 3900x 12-core processor × 24 メモリ：64 GB参考サイト
▶ A interphaseChangeFoam tutorial http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2011/MartinAndersen/Tutorial_interPhaseChangeFoam.pdf
0 一連の流れ (solver非編集)
1 ディレクトリの準備（interPhaseChangeFoamをコピー）
$ cd run
$ cp -r $FOAM_SOLVERS/multiphase/interPhaseChangeFoam myiPCF_2
$ cd myiPCF
$ mv interPhaseChangeFoam.C myiPCF.C
自分のワークフォルダに既存のsolverをコピーしてくる (チュートリアルの"interPhaseChangeFoam")。
今回の自作solverの名前を"myiPCF"とした(my inter Phase Change Foam の略)。
2 "myiPCF.C" (自作solver.C)内のヘッダーの修正
.
├── Allwclean
├── Allwmake
├── Make
│ ├── files
│ └── options
├── UEqn.H
├── alphaControls.H
├── alphaEqn.H
├── alphaEqnSubCycle.H
├── correctPhi.H
├── createFields.H
├── initCorrectPhi.H
├── myiPCF.C (もともとはinterPhaseChangeFoam.C)
├── pEqn.H
└── phaseChangeTwoPhaseMixtures
├── Kunz (略)
├── Make (略)
├── Merkle (略)
├── SchnerrSauer (略)
├── lnInclude (略)
└── phaseChangeTwoPhaseMixture (略)
■ OpenFOAM8の場合 ("myiPCF.C"内のヘッダーを確認)
（上略）
#include "CorrectPhi.H"
（下略）
絶対参照がないのでOpenFoam8では修正必要なさそう
■ 古いOpenFOAMの場合
"interFoam"内の"correctPhi.H"をコピーする。
$ cp $FOAM_SOLVERS/multiphase/interFoam/correctPhi.H .
"myiPCF.C"内のヘッダーを修正する
修正前：(参照ディレクトリが記載している場合)
（上略）
#include "../interFoam/correctPhi.H"
（下略）
修正後：
（上略）
#include "correctPhi.H"
（下略）
3 "Make/files"の修正
修正前：
interPhaseChangeFoam.C
EXE = $(FOAM_APPBIN)/interPhaseChangeFoam
修正後：
myiPCF.C
EXE = $(FOAM_APPBIN)/myiPCF
古いverのOpenFOAMだと"Make/files"の記載内容が異なるが、基本は下記の通り。
①どのファイルを読むか (〇〇.C)
②どんな名前のファイルを出力するか (EXE = $(FOAM_APPBIN)/〇〇)
4 コンパイル
$ wclean
$ wmake 2>&1 | tee log.wmake
■ エラー（outputを書き込み権限がない）
usr/bin/ld.bfd: cannot open output file openfoam
書き込む場所の権限を与える。
$ sudo chown -R tokudy /opt/openfoam8/platforms/
■ 他のエラー "$ less log.wmake"で各自解決してほしい。
5 solverの出力確認
$ cd $FOAM_APPBIN
$ ls
ここに制作した"myiPCF"が出力されているはず。
6 solverの動作確認
$ run
$ cp -r $FOAM_TUTORIALS/multiphase/interPhaseChangeFoam/cavitatingBullet cavitatingBullet_test
$ cd cavitatingBullet_test
$ vi system/controlDict
$ ./Allrun
チュートリアルの"cavitatingBullet"で確認を行う。
"system/controlDict"を制作したsolverで解くように修正する。
"application"を"myiPCF"(自作solver)とし、"endTime"と"writeInterval"を適切な範囲に設定する。
(下記の条件の計算時間：ExecutionTime = 4020.92 s ClockTime = 4025 s)
system/controlDict
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application myiPCF;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 0.02;
deltaT 1e-8;
writeControl adjustableRunTime;
writeInterval 0.001;
purgeWrite 0;
writeFormat binary;
writePrecision 6;
writeCompression off;
timeFormat general;
runTimeModifiable yes;
adjustTimeStep on;
maxCo 5;
maxAlphaCo 2;
// ************************************************************************* //
7 結果
自作ソルバのコンパイルは通り解析結果まで出力できたのでOK。
これが物理的にどうのこうの、定性的、定量的にあーだこーだの議論は保留する。
一連の流れが確認できた、ヨシ！。
deltaT 0.001; endTime 0.02;
paraViewで確認後、jpegで出力し動画制作。コマンドは以下
$ ffmpeg -r 30 -i myiPCF_alpha_water.%04d.jpeg -vcodec libx264 -pix_fmt yuv420p -r 60 210916_myiPCF_alpha_water.mp4
8 おまけ 解析log (最初と最後)
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
Build : 8-1c9b5879390b
Exec : myiPCF
Date : Sep 16 2021
Time : 17:36:24
Host : "tokudy-desktop"
PID : 10259
I/O : uncollated
Case : /home/tokudy/OpenFOAM/tokudy-8/run/cavitatingBullet_test
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModification
Skew 10)
allowSystemOperations : Allowing user-supplied system call operations
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
Create mesh for time = 0
PIMPLE: No convergence criteria found
PIMPLE: Operating solver in transient mode with 1 outer corrector
PIMPLE: Operating solver in PISO mode
Reading field p_rgh
Reading field U
Reading/calculating face flux field phi
Creating phaseChangeTwoPhaseMixture
Selecting phaseChange model SchnerrSauer
Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
Selecting turbulence model type laminar
Selecting laminar stress model Stokes
Reading g
Reading hRef
Calculating field g.h
No finite volume options present
GAMG: Solving for pcorr, Initial residual = 1, Final residual = 9.71108e-09, No Iterations 19
time step continuity errors : sum local = 2.36612e-20, global = -7.53186e-21, cumulative = -7.53186e-21
Courant Number mean: 0.00013526 max: 0.00764258
Starting time loop
Courant Number mean: 0.00013526 max: 0.00764258
Interface Courant Number mean: 0 max: 0
deltaT = 1.19999e-08
Time = 1.19999e-08
smoothSolver: Solving for alpha.water, Initial residual = 1, Final residual = 1.3889, No Iterations 10
Phase-1 volume fraction = 1 Min(alpha.water) = 1 Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 1 Min(alpha.water) = 1 Max(alpha.water) = 1
GAMG: Solving for p_rgh, Initial residual = 1, Final residual = 0.0343803, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 0.572606, Final residual = 2.87102e-05, No Iterations 1
GAMGPCG: Solving for p_rgh, Initial residual = 3.26599e-05, Final residual = 9.21097e-08, No Iterations 2
ExecutionTime = 3.1 s ClockTime = 3 s
-----------------中略-----------------
smoothSolver: Solving for alpha.water, Initial residual = 6.44222e-06, Final residual = 9.08914e-10, No Iterations 2
Phase-1 volume fraction = 0.970299 Min(alpha.water) = 3.13574e-05 Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 0.970221 Min(alpha.water) = 3.13634e-05 Max(alpha.water) = 1.00039
GAMG: Solving for p_rgh, Initial residual = 6.11211e-05, Final residual = 2.00492e-07, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 1.63396e-05, Final residual = 7.53242e-08, No Iterations 1
GAMGPCG: Solving for p_rgh, Initial residual = 5.73702e-06, Final residual = 1.72077e-08, No Iterations 1
ExecutionTime = 4019.57 s ClockTime = 4023 s
Courant Number mean: 0.0908173 max: 2.12591
Interface Courant Number mean: 0.00579412 max: 1.98598
deltaT = 7.0893e-06
Time = 0.02
smoothSolver: Solving for alpha.water, Initial residual = 6.45889e-06, Final residual = 9.11336e-10, No Iterations 2
Phase-1 volume fraction = 0.970295 Min(alpha.water) = 3.1736e-05 Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 0.970218 Min(alpha.water) = 3.17365e-05 Max(alpha.water) = 1.00035
GAMG: Solving for p_rgh, Initial residual = 6.1257e-05, Final residual = 1.99878e-07, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 1.64261e-05, Final residual = 7.32739e-08, No Iterations 1
GAMGPCG: Solving for p_rgh, Initial residual = 5.80019e-06, Final residual = 1.77596e-08, No Iterations 1
ExecutionTime = 4020.92 s ClockTime = 4025 s
End
1 solverの中身の理解
状況の確認
前章までに一連の流れ (ディレクトリの作製〜solverのコンパイル〜テストケースの検証まで)を確認した。本章では、メイントピック(新たなヘッダーファイルを作成し、solverのコードを修正し、オリジナルのsolverをコンパイルする)の大前提となる"ソースコードの内容を理解"を試みる。まず、ソースコード "myiPCF.C" (自作コード)を見てみる。しかし、何が書いてあるかはよくわからない (一応解釈を試みるが各自調べてほしい)。
Useful Reference
▶ Description and utilization of interFoam multiphase solver
http://infofich.unl.edu.ar/upload/3be0e16065026527477b4b948c4caa7523c8ea52.pdf
▶ interFoam
http://penguinitis.g1.xrea.com/study/OpenFOAM/tankentai/22-interFoam.html
▶ 【OpenFOAM】手法を変えて処理を追加② : ソルバーに処理を追加する
https://zenn.dev/inabower/articles/bd0406a50919e6
▶ OpenFOAMライブラリリファレンス (株式会社テラバイト, 2020)
https://www.amazon.co.jp/dp/4627691610/
▶ OpenFOAMプログラミング (森北出版, 2017)
https://www.amazon.co.jp/dp/4627670915/
▶ 液体食品の充填製造プロセスにおける流動操作条件の最適化
file:///tmp/LID201704201006.pdf
https://repo.lib.tokushima-u.ac.jp/ja/110074
myiPCF.C
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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
" while (pimple.run(runTime))"が数値解析の本体。中を見てみる。
まず、計算を開始する前に、メッシュ、時間、クーラン数などの基本的な値を確認する。
#include "readDyMControls.H"
"correctPhi", "checkMeshCourantNo", "moveMeshOuterCorrectors"の取得。
flux variableである\(\,\phi\)、クーラン数、移動メッシュ？を取得するっぽいがよくわからない。
型はboolだと思う。
▶ readDyMControls.H
https://cpp.openfoam.org/v8/readDyMControls_8H_source.html
readDyMControls.H
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault
(
"correctPhi",
correctPhi
);
checkMeshCourantNo = pimple.dict().lookupOrDefault
(
"checkMeshCourantNo",
checkMeshCourantNo
);
moveMeshOuterCorrectors = pimple.dict().lookupOrDefault
(
"moveMeshOuterCorrectors",
moveMeshOuterCorrectors
);
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
"volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));"で値をコピーする。補正された\(\,\phi\)が同じ発散を有するように、前のメッシュから\(\,{\nabla}U\)を保存し、それをマッピングしてcorrectPhiで使用できるようにする。
"volScalarField"に関しては、下記に簡潔に書いてある。
▶ スカラーフィールド
http://penguinitis.g1.xrea.com/study/OpenFOAM/tankentai/07-sfield.html
▶ OpenFOAM プログラミングメモ
http://penguinitis.g1.xrea.com/study/OpenFOAM/programming_memo.html
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
ここもよくわからないが、クーラン数とタイムステップを調整している？
▶ 液体食品の充填製造プロセスにおける流動操作条件の最適化
file:///tmp/LID201704201006.pdf
https://repo.lib.tokushima-u.ac.jp/ja/110074
次に、時間の更新。
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
次に、速度と圧力の連成解析を行う。ここではPIMPLE法を使用している。PIMPLE法については後述。
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
{
mesh.update();
if (mesh.changing())
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
if (correctPhi)
{
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & Uf();
#include "correctPhi.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
}
mixture.correct();
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
}
次に、上記の連成解析で得られた圧力を用いて、仮の速度を出す。
divU = fvc::div(fvc::absolute(phi, U));
次に、密度、flux？を更新する。
surfaceScalarField rhoPhi
(
- 略 -
);
"alpha"VOFパラメータ？に関してあれやこれやする。
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
たぶん混合密度、混合粘性の修正？
mixture.correct();
速度について解く。
#include "UEqn.H"
圧力について解く？修正？。
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
よくわからん。乱流に関して修正？
if (pimple.turbCorr())
{
turbulence->correct();
}
以上を繰り返す。
PIMPLE法は"PIS0"と"SIMPLE"を組み合わせた方法で、アルゴリズムは、
① 現在の速度 (n-step)から仮の速度 (*-step)を算出する。
\begin{align}
u_* = u_n + c_1(A_{AD,*} + A_{NU,*} + A_{p,n})
\end{align}
② 圧力の補正量を算出する。
\begin{align}
{\nabla}{\cdot}u = 0
\end{align}
\begin{align}
u_{n+1} = u_* - c_2(P_{correct})
\end{align}
\begin{align}
{\nabla}{\cdot}(A_{p,*}) = \cfrac{1}{c_2}{\nabla}{\cdot}u_*
\end{align}
③ 圧力と速度を補正補正する。
\begin{align}
u_{n+1} = u_* - c_2(P_{correct})
\end{align}
④ 値が収束するまで繰り返す。
なお、実際の(数学的な)アルゴリズムは下記のようになる。OpenFOAMライブラリリファレンス (第15章 自由表面流れの解析手法, p648~650)を参考にした。
支配方程式
非圧縮性流体の運動保存式は、気液界面での表面張力と重量を考慮すると、
\begin{align}
\cfrac{{\partial}({\rho}u_i)}{{\partial}t}
+ \cfrac{{\partial}({\rho}{u_i}{u_j})}{{\partial}x_j}
=
- \cfrac{{\partial}p_r}{{\partial}x_i}
- g_j(x_j-r_j)\cfrac{{\partial}{\rho}}{{\partial}x_j}
+ \cfrac{\partial}{{\partial}x_j}
\left[
{\mu} (\cfrac{{\partial}u_i}{{\partial}x_j}
+ \cfrac{{\partial}u_j}{{\partial}x_i})
\right]
+ F_{CSF,i}
\end{align}
時間変化 + 移流項 = 圧力項 + 重力項 + 粘性項 + 表面張力項
\(\, r_i\)は参照座標。\(\, g_j\)は重力。
圧力\(\, p_r\)は静水圧の寄与を取り除いた値。
\begin{align}
p_r = p - {\rho}g_j(x_j-r_j)
\end{align}
右辺の最終項は表面張力項(単相流には無い)。CSF (Continuum SUrface Force)モデルとして扱われる。
\begin{align}
F_{CSF,i} = {\sigma}{\kappa}\cfrac{{\partial}F}{{\partial}x_i}
\end{align}
\(\, \sigma\)は表面張力係数。界面の曲率\(\,\kappa\)は、
\begin{align}
\kappa = - {\nabla}{\cdot}{\bf{n}} = {\nabla}{\cdot}
\left(
{\cfrac{{\nabla}F}{|{\nabla}F|}}
\right)
\end{align}
混合密度\(\,\rho\)と混合粘度\(\,\mu\)は、VOF関数体積率\(\,F\)を用いて、
\begin{align}
\rho = \rho_2 + F(\rho_1 - \rho_2)
\end{align}
\begin{align}
\mu = \mu_2 + F(\mu_1 - \mu_2)
\end{align}
interFoamで実装されているPIMPLE方での計算手順
(1) 流速の初期値、体積率\(\,F\)を設定し、セル境界面のフラックス\(\,\phi(u_i^0)\)を計算。
(2) VOF関数の移流方程式を解いて、体積率\(\,F\)を算出 (MULES or CMULES)。
\begin{align}
\cfrac{{\partial}F}{{\partial}t}
+ \cfrac{{\partial}(F{u_j})}{{\partial}x_j}
= 0
\end{align}
(3) 運動量保存式を解く。
\begin{align}
A_P u_{iP}^{m*}
+ {\sum\limits_{l}^{NC}} A_l u_{il}^{m*}
= Q_{iP}^{m-1}
- \left(
\cfrac{{\delta}p^{m-1}} {{\delta} x_i}
\right)_P \Delta V
\end{align}
(注：上式を解かない場合、fvSolutionファイルのPIMPLE辞書のmomentunPredictorキーワードにfalseを指定する。)
(4) 運動量保存式から得られた流速\(\,u_i^{m*}\)を用いて暫定流速を獲得 ((3)を解かない場合は前タイムステップの\(\,u_i^n\))。
\begin{align}
\tilde{u}_{iP}^{**} = \cfrac{1}{A_P}
\left(
- \sum\limits_{l}^{NC} A_l u_{il}^{m*} + Q_{iP}^{m-1}
\right)
\end{align}
(5) 暫定流速を用いて、Rhie-Chow補間により暫定境界面流束を計算。
\begin{align}
\phi_{\alpha}(\tilde{u}_i^{**})
= {\langle \tilde{u}_i^{**} \rangle}_{\alpha} S_{i \alpha}
+ \left\langle \cfrac{\rho \Delta V}{A_P} \right\rangle _{\alpha}
\left(
\cfrac{\phi _\alpha (u_i^{n-1}) - {\langle u_i^{n-1} \rangle}_\alpha S_{i\alpha}}{\Delta t}
\right)
L \left( \phi(u_i^{n-1}) \right)
- \left\langle \cfrac{\Delta V}{A_P} \right\rangle _{\alpha}
\left(
F_{CSF,i\alpha} - g_j (x_j - r_j) \cfrac{\delta \rho}{\delta x_i}
\right)_\alpha
S_{i \alpha}
\end{align}
ここで、制限関数\(\, L \)は、
\begin{align}
L \left( \phi(u_i^{n-1}) \right)
= 1 - \min
\left(
\cfrac{ \left\lvert \phi (u_i^{n-1}) - {\langle u_i^{n-1} \rangle} S_i \right\rvert \Delta t}{d_{CF} \left\lvert \bf{S} \right\rvert}
, 1
\right)
\end{align}
(6) 圧力方程式をlaplacianスキーム (反復計算)で解く。
\begin{align}
\sum\limits_{\alpha}^{NS}
\left\langle \cfrac{\Delta V}{A_P} \right\rangle _{\alpha}
\left( \cfrac{ \delta p^m }{ \delta x_i } \right)_{\alpha}
S_{i \alpha}
=
\sum\limits_{\alpha}^{NS}
\phi_{\alpha}(\tilde{u}_i^{**})
\end{align}
(7) 次サイクルで必要となるセル境界面流束を計算する。
\begin{align}
\phi _\alpha (u_i^m) = \left\langle \tilde{u}_i^{**} \right\rangle _\alpha S_{i\alpha}
- \left\langle \cfrac{\Delta V}{A_P} \right\rangle _{\alpha}
\left( \cfrac{ \delta p^m }{ \delta x_i } \right)_{\alpha}
S_{i \alpha}
\end{align}
(8) 流速を更新する。
\begin{align}
u_i^m = \tilde{u}_i^{**}
- \cfrac{\Delta V}{A_P}
\left(
F_{CSF,i} - g_j (x_j - r_j) \cfrac{\delta \rho}{\delta x_i}
\right)_\alpha
- \cfrac{\Delta V}{A_P}
\cfrac{\delta p^m}{\delta x_i}
\end{align}
(9) 圧力補正と流速補正を繰り返す(PISOループ)ため、
\(\,u_i^{m*} = u_i^m\)
として、STEP(4)に戻り、STEP(8)まで繰り返す。
(10)
運動量保存式の係数行列を更新する外部反復ループを実行するために、STEP(2)に戻り、
STEP(9)まで反復計算する。収束したと判断された場合、
\(\,u_i^{n+1} = u_i^m\)
として次の計算ステップに進む。
実は移流方程式を解く方法
http://t2r2.star.titech.ac.jp/rrws/file/CTT100759558/ATD100000413/outline_15D31012_%E6%B1%A0%E7%AB%AF%E6%98%AD%E5%A4%AB.pdf
http://i-rep.emu.edu.tr:8080/xmlui/bitstream/handle/11129/1610/EmadVahid.pdf?sequence=1
https://sites.google.com/site/trujillolabgroup/understanding-mules
http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2017/ElinOlsson/Report_ElinOlsson_Updated050118.pdf
2 熱拡散方程式の追加
一連の実行コマンド (solver側)
$ run
$ cd myiPCF
$ wclean
$ vi createFields.H
$ touch TEqn.H
$ vi TEqn.H
$ vi myiPCF.C
$ wmake
追加内容: 熱拡散項・支配方程式の実装
熱拡散 (thermal diffusivity)を実装する。通常は\(\,\alpha\)で表されるがコーディングの都合上\(\, DT\)とする (\(\,\alpha\)はVOF値に取られてる)。
\(\, DT = \cfrac{k}{\rho c_p}\)
\(\, DT\): thermal diffusivity [m^2/S]
\(\, k\): thermal conductivity [W/(m*K)]
\(\, \rho\): density [kg/m^3]
\(\, c_p\): specific heat capacity [J/kg*K]
上式の内容を"createFields.H"の一番上 (おそらくどこでもよい)に書き足す。
createFields.H
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportPropertiesDict
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar DT
(
transportPropertiesDict.lookup("DT")
);
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr mixturePtr
(
phaseChangeTwoPhaseMixture::New(U, phi)
);
phaseChangeTwoPhaseMixture& mixture = mixturePtr();
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Construct incompressible turbulence model
autoPtr turbulence
(
incompressible::momentumTransportModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "createFvOptions.H"
2~17行はthermal diffusion (DTとして)の新しいtransport prppetyの記述。
18~30行は温度分布Tの生成。
支配方程式 (熱拡散方程式)
\begin{align}
\cfrac{{\partial}T}{{\partial}t} + {\nabla}{\cdot}({\phi}T) = {\nabla}{\cdot}(DT{\nabla}T)
\end{align}
上式をコードで表すと下記のようになる。新しく"TEqn.H"というファイルを作り、このコードを書き込む。
TEqn.H
{
fvScalarMatrix TEqn
(s
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
);
TEqn.solve();
}
上記ファイルを作っただけではだめなので、自作solverの"myiPCF.C"にて #include をする。おそらく、#include のタイミングが重要であり、下記の箇所に書き込む ("UEqn.H"の後なのは速度の項を流用するからか？)。
myiPCF.C
- 略 -
while (pimple.run(runTime))
{
- 略 -
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
- 略 -
}
divU = fvc::div(fvc::absolute(phi, U));
surfaceScalarField rhoPhi
(
- 略 -
);
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
mixture.correct();
#include "UEqn.H"
// ここに追加する↓
#include "TEqn.H
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
- 略 -
古いverのOpenFOAMは下記のように編集するっぽい。
flux variable \(\phi\)に関しては、 圧力のループ計算 (ポアソン方程式)のにおいて使用されている (それを流用するものだと思われる)。熱輸送は速度分布に依存するため、"Pressure-velocity PIMPLE corrector loop"のループ解析後に、TEqn.Hのファイルはincludeされる。phase changeはこのエネルギー方程式に依存し、"twoPhaseProperties->correct();:"(昔のOpenFOAMのコード)の前にincludeされなければならない。
一連の実行コマンド (case 側)
$ cp -r $FOAM_TUTORIALS/multiphase/interPhaseChangeFoam/cavitatingBullet $FOAM_RUN/temperatureTest
$ cd $FOAM_RUN/temperatureTest
$ vi Allrun
$ rm -r constant/triSurface system/snappyHexMeshDict
$ sed -i "/bullet/,/}/d" 0/*
$ vi constant/transportProperties
$ cp 0/p_rgh 0/T
$ vi 0/T
$ vi system/fvSchemes
$ vi system/fvSolution
$ vi sysetm/controlDict
$ ./Allrun
$ paraFoam
$ ffmpeg -r 30 -i 210919_cavitationBullet_myiPCF_T.%04d.jpeg -vcodec libx264 -pix_fmt yuv420p -r 60 210919_cavitationBullet_myiPCF_T.mp4
"cavitationBullet"のケースディレクトリをコピーしてくる。しかし、簡単に確認するために、bulletの部分は取り除いて自作solverのコンパイルの確認を行う。
まず、"Allrun"から、"potentialFoam"と"snappyHexMesh"の記述を取り除き下記のようにする。
Allrun
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Generate the base block mesh
runApplication blockMesh
# Run the solver
runApplication $(getApplication)
#------------------------------------------------------------------------------
次に、bulletのジオメトリやコードの記述を取り除く。
$ rm -r constant/triSurface system/snappyHexMeshDict
$ sed -i "/bullet/,/}/d" 0/*
メモ：
「sed」は「Stream EDitor」の略。「sed スクリプトコマンド ファイル名」で、指定したファイルをスクリプトコマンドに従って行単位に処理し、標準出力へ出力する。ファイル名を省略した場合は、標準入力からのデータを処理する。
「-i」 ファイルを直接編集する。
「s/置換前/置換後/」
［置換前］で指定した文字列にマッチした部分を［置換後］に置き換える。複数マッチした場合は先頭のみ置換、全てを置換したい場合は、「s/置換前/置換後/g」のように「g」オプションを指定する。
"transportProperties"に"DT"に関して定義する。
場所はどこでもよく、最後に追加した。
constant/transportProperties
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
phases (water vapour);
phaseChangeTwoPhaseMixture SchnerrSauer;
pSat 2300; // Saturation pressure
sigma 0.07;
water
{
transportModel Newtonian;
nu 9e-07;
rho 1000;
}
vapour
{
transportModel Newtonian;
nu 4.273e-04;
rho 0.02308;
}
KunzCoeffs
{
UInf U20.0;
tInf t0.005; // L = 0.1 m
Cc C1000;
Cv C1000;
}
MerkleCoeffs
{
UInf 20.0;
tInf 0.005; // L = 0.1 m
Cc 80;
Cv 1e-03;
}
SchnerrSauerCoeffs
{
n 1.6e+13;
dNuc 2.0e-06;
Cc 1;
Cv 1;
}
DT DT [0 2 -1 0 0 0 0] 0.002; // [m^2/s]
// ************************************************************************* //
"T"の設定。"0/p_rgh"をコピーして下記のように記述する。
0/T
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0]; // [K]
internalField uniform 363.15; // 90 deg. C
boundaryField
{
inlet
{
type fixedValue;
value uniform 383.15; // 110 deg. C
}
outlet
{
type fixedValue;
value $internalField;
}
walls
{
type symmetry;
}
}
// ************************************************************************* //
"system/fvSchemes"の設定。支配方程式を見ても明らかなように、"divSchemes"と"laplacianSchemes"の部分を修正する必要がある。
system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
interpolationSchemes
{
default linear;
}
divSchemes
{
default none;
div(rhoPhi,U) Gauss linearUpwind grad(U);
div(phi,omega) Gauss linearUpwind grad(omega);
div(phi,k) Gauss linearUpwind grad(k);
div(phi,alpha) Gauss interfaceCompression vanLeer 1;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
div(phi,T) Gauss upwind;
}
gradSchemes
{
default Gauss linear;
}
laplacianSchemes
{
default Gauss linear limited corrected 0.5;
laplacian(DT,T) Gauss linear corrected;
}
snGradSchemes
{
default limited corrected 0.5;
}
// ************************************************************************* //
"system/fvSolution"の設定。最初は"T"で設定していたがコンパイル通らず、"TFinal"にしたらうまくいった。solverはとりあえず"PBiCG"にしたがどれが良いかは不明。
system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
TFinal
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-7;
relTol 0;
};
"alpha.water.*"
{
nAlphaCorr 2;
nAlphaSubCycles 1;
MULESCorr yes;
nLimiterIter 5;
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-8;
relTol 0;
maxIter 10;
};
"U.*"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-6;
relTol 0;
};
p_rgh
{
solver GAMG;
tolerance 1e-8;
relTol 0.1;
smoother DICGaussSeidel;
maxIter 50;
};
p_rghFinal
{
solver PCG;
preconditioner
{
preconditioner GAMG;
tolerance 1e-6;
relTol 0;
nVcycles 2;
smoother DICGaussSeidel;
};
tolerance 1e-7;
relTol 0;
maxIter 50;
};
"pcorr.*"
{
$p_rgh;
relTol 0;
};
Phi
{
$p_rgh;
relTol 0;
};
}
potentialFlow
{
nNonOrthogonalCorrectors 3;
}
PIMPLE
{
momentumPredictor no;
nOuterCorrectors 1;
nCorrectors 3;
nNonOrthogonalCorrectors 0;
}
relaxationFactors
{
equations
{
".*" 1;
}
}
// ************************************************************************* //
"controlDict"で自作solverを設定することを忘れないように注意する。また、今回はbulletなしの簡単な解析のため、"endTime"は小さく、非定常解析のため"writeInterval"は細かく設定する。
system/controlDict
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application myiPCF; // 追加したところ
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 1e-3; // 追加したところ
deltaT 1e-8;
writeControl adjustableRunTime;
writeInterval 1e-5; // 追加したところ
purgeWrite 0;
writeFormat binary;
writePrecision 6;
writeCompression off;
timeFormat general;
runTimeModifiable yes;
adjustTimeStep on;
maxCo 5;
maxAlphaCo 2;
// ************************************************************************* //
Results
とりあえず熱拡散はうまくいってそう。
deltaT 1e-5; endTime 1e-3;
3 飽和蒸気圧曲線の追加 Psat(T)
August-Roche-Magnus formula
solver内にて飽和蒸気圧は２つのループで計算されている (\alpha -loopと pressure-loop)。
計算の不必要なとき (高速化や非線形方程式を再研鑽するとき)、温度場が更新されたとき飽和蒸気圧だけで計算され、圧力場が再構築される。ただし、必ずしもすべてのタイムステップにおいてそれが実行されるとは限らない。
$ touch calcPSatField.H
$ vi calcPSatField.H
calcPSatField.H
{
const dimensionedScalar t30_11("30.11", dimensionSet(0,0,0,1,0,0,0), 30.11);
const dimensionedScalar t273_15("273.15", dimensionSet(0,0,0,1,0,0,0), 273.15);
const dimensionedScalar t1("1", dimensionSet(0,0,0,1,0,0,0), 1);
const dimensionedScalar p610_94("610.94", dimensionSet(1,-1,-2,0,0,0,0), 610.94);
// dimensionSet( [kg], [m], [s], [K], [kg*mol], [A], [cd]), [kg/(m*S^2)]=[Pa]
// August-Roche-Magnus formula
pSat = p610_94 * exp( 17.625*(T-t273_15) / max(t1, T-t30_11) );
//max(1,...) is included to avoid problems with devision by 0
}
単位が必ず必要。ここではまだPsatを決めてない。
$ vi createFields.H
以下の一文を取り除く。
createFields.H
const dimensionedScalar& pSat = twoPhaseProperties->pSat();
そして、"volScalarField T"の後にincludeする文を付け加える。
"volScalarField p_rgh"の後。(飽和蒸気圧場の初期化のために"volScalarField p_rgh"が"p_rgh"として用いられるので。)
"phaseChangeTwoPhaseMixture"の前。(飽和蒸気圧が"phaseChangeTwoPhaseMixture"に用いられるので。)
createFields.H
volScalarField pSat
(
IOobject
(
"pSat",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh // initial value will be overwritten in calcPSatField.H
);
#include "calcPSatField.H"
温度"T"に依存する飽和蒸気圧を更新するには、"myiPCF.C"内に｢#include "TEqn.H"｣の後に、｢#include "calcPSatField.H"｣を記載する。
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportPropertiesDict
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar DT
(
transportPropertiesDict.lookup("DT")
);
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField pSat
(
IOobject
(
"pSat",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh // initial value will be overwritten in calcPSatField.H
);
#include "calcPSatField.H"
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr mixturePtr
(
phaseChangeTwoPhaseMixture::New(U, phi)
);
phaseChangeTwoPhaseMixture& mixture = mixturePtr();
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Construct incompressible turbulence model
autoPtr turbulence
(
incompressible::momentumTransportModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "createFvOptions.H"
/myiPCF/phaseChangeTwoPhaseMixtures
.
├── Kunz
│ ├── Kunz.C
│ └── Kunz.H
├── Make
│ ├── files
│ └── options
├── Merkle
│ ├── Merkle.C
│ └── Merkle.H
├── SchnerrSauer
│ ├── SchnerrSauer.C
│ └── SchnerrSauer.H
├── lnInclude
│ ├── Kunz.C -> ../Kunz/Kunz.C
│ ├── Kunz.H -> ../Kunz/Kunz.H
│ ├── Merkle.C -> ../Merkle/Merkle.C
│ ├── Merkle.H -> ../Merkle/Merkle.H
│ ├── SchnerrSauer.C -> ../SchnerrSauer/SchnerrSauer.C
│ ├── SchnerrSauer.H -> ../SchnerrSauer/SchnerrSauer.H
│ ├── phaseChangeTwoPhaseMixture.C -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
│ ├── phaseChangeTwoPhaseMixture.H -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
│ └── phaseChangeTwoPhaseMixtureNew.C -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
└── phaseChangeTwoPhaseMixture
├── phaseChangeTwoPhaseMixture.C
├── phaseChangeTwoPhaseMixture.H
└── phaseChangeTwoPhaseMixtureNew.C
6 directories, 20 files
stacic saturation pressure を取り除く。
$ vi phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
次の記述
//- Saturation vapour pressure
dimensionedScalar pSat_;
■ OpenFOAM8
修正前：
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
修正後：
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
return pSat_;
}
■ 古いOpenFOAM
修正前：
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
修正後：
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
return pSat;
}
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 New
(
const volVectorField& U,
const surfaceScalarField& phi
);
//- Destructor
virtual ~phaseChangeTwoPhaseMixture()
{}
// Member Functions
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
//- Return the mass condensation and vaporisation rates as a
// coefficient to multiply (1 - alphal) for the condensation rate
// and a coefficient to multiply alphal for the vaporisation rate
virtual Pair> mDotAlphal() const = 0;
//- Return the mass condensation and vaporisation rates as coefficients
// to multiply (p - pSat)
virtual Pair> mDotP() const = 0;
//- Return the volumetric condensation and vaporisation rates as a
// coefficient to multiply (1 - alphal) for the condensation rate
// and a coefficient to multiply alphal for the vaporisation rate
Pair> vDotAlphal() const;
//- Return the volumetric condensation and vaporisation rates as
// coefficients to multiply (p - pSat)
Pair> vDotP() const;
//- Correct the phaseChange model
virtual void correct() = 0;
//- Read the transportProperties dictionary and update
virtual bool read() = 0;
// Member Operators
//- Disallow default bitwise assignment
void operator=(const phaseChangeTwoPhaseMixture&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
pSatに関係するすべての記述を消去。
以下の3行目と2行目最後のコンマを消す。
$ vi phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
■ OpenFOAM8
修正前：
immiscibleIncompressibleTwoPhaseMixture(U, phi),
phaseChangeTwoPhaseMixtureCoeffs_(optionalSubDict(type + "Coeffs")),
pSat_("pSat", dimPressure, lookup("pSat"))
{}
修正後：
immiscibleIncompressibleTwoPhaseMixture(U, phi),
phaseChangeTwoPhaseMixtureCoeffs_(optionalSubDict(type + "Coeffs"))
{}
■ 古いOpenFOAM
修正前：
twoPhaseMixture(U, phi, alpha1Name),
phaseChangeTwoPhaseMixtureCoeffs_(subDict(type + "Coeffs")),
pSat_(lookup("pSat"))
{}
修正後：
twoPhaseMixture(U, phi, alpha1Name),
phaseChangeTwoPhaseMixtureCoeffs_(subDict(type + "Coeffs"))
{}
次の記述を消す。
lookup("pSat") >> pSat_;
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 >
Foam::phaseChangeTwoPhaseMixture::vDotAlphal() const
{
volScalarField alphalCoeff(1.0/rho1() - alpha1()*(1.0/rho1() - 1.0/rho2()));
Pair> mDotAlphal = this->mDotAlphal();
return Pair>
(
alphalCoeff*mDotAlphal[0],
alphalCoeff*mDotAlphal[1]
);
}
Foam::Pair>
Foam::phaseChangeTwoPhaseMixture::vDotP() const
{
dimensionedScalar pCoeff(1.0/rho1() - 1.0/rho2());
Pair> mDotP = this->mDotP();
return Pair>(pCoeff*mDotP[0], pCoeff*mDotP[1]);
}
void Foam::phaseChangeTwoPhaseMixture::correct()
{
immiscibleIncompressibleTwoPhaseMixture::correct();
}
bool Foam::phaseChangeTwoPhaseMixture::read()
{
if (immiscibleIncompressibleTwoPhaseMixture::read())
{
phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
lookup("pSat") >> pSat_;
return true;
}
else
{
return false;
}
}
// ************************************************************************* //
以上で修正終了。
$ wclean
$ wmake
躓いたエラー
エラーログ
Making dependency list for source file myiPCF.C
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -c myiPCF.C -o Make/linux64GccDPInt32Opt/myiPCF.o
In file included from myiPCF.C:49:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H: In member function ‘const volScalarField& Foam::phaseChangeTwoPhaseMixture::pSat() const’:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:124:42: error: reference to ‘alpha1_’ is ambiguous
124 | const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
| ^~~~~~~
In file included from myiPCF.C:48:
/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude/interfaceProperties.H:69:31: note: candidates are: ‘const volScalarField& Foam::interfaceProperties::alpha1_’
69 | const volScalarField& alpha1_;
| ^~~~~~~
In file included from /opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude/incompressibleTwoPhaseMixture.H:40,
from /opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude/immiscibleIncompressibleTwoPhaseMixture.H:38,
from phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:38,
from myiPCF.C:49:
/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude/twoPhaseMixture.H:58:24: note: ‘Foam::volScalarField Foam::twoPhaseMixture::alpha1_’
58 | volScalarField alpha1_;
| ^~~~~~~
In file included from myiPCF.C:49:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:124:82: error: expected primary-expression before ‘>’ token
124 | const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
| ^
make: *** [/opt/openfoam8/wmake/rules/General/transform:26: Make/linux64GccDPInt32Opt/myiPCF.o] エラー 1
たぶん、‘alpha1_’ の定義がはっきりしてない。
"Foam::interfaceProperties::alpha1_"や"Foam::twoPhaseMixture::alpha1_"が候補として出ているが、ここではひとまず、"Foam::twoPhaseMixture::alpha1_"とした。
修正内容
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = Foam::twoPhaseMixture::alpha1_.db().lookupObject("pSat");
return pSat;
}
以上で修正終了。
$ wclean
$ wmake
コンパイルログ
tokudy@tokudy-desktop:~/OpenFOAM/tokudy-8/run/myiPCFPsat$ wmake
Making dependency list for source file myiPCF.C
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -c myiPCF.C -o Make/linux64GccDPInt32Opt/myiPCF.o
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -fuse-ld=bfd -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPInt32Opt/myiPCF.o -L/opt/openfoam8/platforms/linux64GccDPInt32Opt/lib \
-lphaseChangeTwoPhaseMixtures -ltwoPhaseMixture -linterfaceProperties -ltwoPhaseProperties -lincompressibleTransportModels -lmomentumTransportModels -lincompressibleMomentumTransportModels -lfiniteVolume -ldynamicFvMesh -ltopoChangerFvMesh -lfvOptions -lmeshTools -lOpenFOAM -ldl \
-lm -o /opt/openfoam8/platforms/linux64GccDPInt32Opt/bin/myiPCF
コンパイル通っているぽい。おｋ．
状況の確認
前章までに一連の流れ (ディレクトリの作製〜solverのコンパイル〜テストケースの検証まで)を確認した。本章では、メイントピック(新たなヘッダーファイルを作成し、solverのコードを修正し、オリジナルのsolverをコンパイルする)の大前提となる"ソースコードの内容を理解"を試みる。まず、ソースコード "myiPCF.C" (自作コード)を見てみる。しかし、何が書いてあるかはよくわからない (一応解釈を試みるが各自調べてほしい)。Useful Reference
▶ Description and utilization of interFoam multiphase solver
http://infofich.unl.edu.ar/upload/3be0e16065026527477b4b948c4caa7523c8ea52.pdf
▶ interFoam
http://penguinitis.g1.xrea.com/study/OpenFOAM/tankentai/22-interFoam.html
▶ 【OpenFOAM】手法を変えて処理を追加② : ソルバーに処理を追加する
https://zenn.dev/inabower/articles/bd0406a50919e6
▶ OpenFOAMライブラリリファレンス (株式会社テラバイト, 2020)
https://www.amazon.co.jp/dp/4627691610/
▶ OpenFOAMプログラミング (森北出版, 2017)
https://www.amazon.co.jp/dp/4627670915/
▶ 液体食品の充填製造プロセスにおける流動操作条件の最適化
file:///tmp/LID201704201006.pdf
https://repo.lib.tokushima-u.ac.jp/ja/110074
myiPCF.C
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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
" while (pimple.run(runTime))"が数値解析の本体。中を見てみる。
まず、計算を開始する前に、メッシュ、時間、クーラン数などの基本的な値を確認する。
#include "readDyMControls.H""correctPhi", "checkMeshCourantNo", "moveMeshOuterCorrectors"の取得。 flux variableである\(\,\phi\)、クーラン数、移動メッシュ？を取得するっぽいがよくわからない。 型はboolだと思う。
▶ readDyMControls.H
https://cpp.openfoam.org/v8/readDyMControls_8H_source.html
readDyMControls.H
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault
(
"correctPhi",
correctPhi
);
checkMeshCourantNo = pimple.dict().lookupOrDefault
(
"checkMeshCourantNo",
checkMeshCourantNo
);
moveMeshOuterCorrectors = pimple.dict().lookupOrDefault
(
"moveMeshOuterCorrectors",
moveMeshOuterCorrectors
);
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
"volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));"で値をコピーする。補正された\(\,\phi\)が同じ発散を有するように、前のメッシュから\(\,{\nabla}U\)を保存し、それをマッピングしてcorrectPhiで使用できるようにする。
"volScalarField"に関しては、下記に簡潔に書いてある。
▶ スカラーフィールド
http://penguinitis.g1.xrea.com/study/OpenFOAM/tankentai/07-sfield.html
▶ OpenFOAM プログラミングメモ
http://penguinitis.g1.xrea.com/study/OpenFOAM/programming_memo.html
#include "CourantNo.H" #include "alphaCourantNo.H" #include "setDeltaT.H"ここもよくわからないが、クーラン数とタイムステップを調整している？
▶ 液体食品の充填製造プロセスにおける流動操作条件の最適化
file:///tmp/LID201704201006.pdf
https://repo.lib.tokushima-u.ac.jp/ja/110074
次に、時間の更新。
runTime++; Info<< "Time = " << runTime.timeName() << nl << endl;
次に、速度と圧力の連成解析を行う。ここではPIMPLE法を使用している。PIMPLE法については後述。
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
{
mesh.update();
if (mesh.changing())
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
if (correctPhi)
{
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & Uf();
#include "correctPhi.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
}
mixture.correct();
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
}
次に、上記の連成解析で得られた圧力を用いて、仮の速度を出す。
divU = fvc::div(fvc::absolute(phi, U));
次に、密度、flux？を更新する。
surfaceScalarField rhoPhi
(
- 略 -
);
"alpha"VOFパラメータ？に関してあれやこれやする。
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
たぶん混合密度、混合粘性の修正？
mixture.correct();
速度について解く。
#include "UEqn.H"
圧力について解く？修正？。
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
よくわからん。乱流に関して修正？
if (pimple.turbCorr())
{
turbulence->correct();
}
以上を繰り返す。 PIMPLE法は"PIS0"と"SIMPLE"を組み合わせた方法で、アルゴリズムは、
① 現在の速度 (n-step)から仮の速度 (*-step)を算出する。 \begin{align} u_* = u_n + c_1(A_{AD,*} + A_{NU,*} + A_{p,n}) \end{align}
② 圧力の補正量を算出する。 \begin{align} {\nabla}{\cdot}u = 0 \end{align}
\begin{align} u_{n+1} = u_* - c_2(P_{correct}) \end{align}
\begin{align} {\nabla}{\cdot}(A_{p,*}) = \cfrac{1}{c_2}{\nabla}{\cdot}u_* \end{align}
③ 圧力と速度を補正補正する。 \begin{align} u_{n+1} = u_* - c_2(P_{correct}) \end{align} ④ 値が収束するまで繰り返す。
なお、実際の(数学的な)アルゴリズムは下記のようになる。OpenFOAMライブラリリファレンス (第15章 自由表面流れの解析手法, p648~650)を参考にした。
実は移流方程式を解く方法 http://t2r2.star.titech.ac.jp/rrws/file/CTT100759558/ATD100000413/outline_15D31012_%E6%B1%A0%E7%AB%AF%E6%98%AD%E5%A4%AB.pdf http://i-rep.emu.edu.tr:8080/xmlui/bitstream/handle/11129/1610/EmadVahid.pdf?sequence=1 https://sites.google.com/site/trujillolabgroup/understanding-mules http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2017/ElinOlsson/Report_ElinOlsson_Updated050118.pdf支配方程式
非圧縮性流体の運動保存式は、気液界面での表面張力と重量を考慮すると、 \begin{align} \cfrac{{\partial}({\rho}u_i)}{{\partial}t} + \cfrac{{\partial}({\rho}{u_i}{u_j})}{{\partial}x_j} = - \cfrac{{\partial}p_r}{{\partial}x_i} - g_j(x_j-r_j)\cfrac{{\partial}{\rho}}{{\partial}x_j} + \cfrac{\partial}{{\partial}x_j} \left[ {\mu} (\cfrac{{\partial}u_i}{{\partial}x_j} + \cfrac{{\partial}u_j}{{\partial}x_i}) \right] + F_{CSF,i} \end{align}
時間変化 + 移流項 = 圧力項 + 重力項 + 粘性項 + 表面張力項
\(\, r_i\)は参照座標。\(\, g_j\)は重力。 圧力\(\, p_r\)は静水圧の寄与を取り除いた値。 \begin{align} p_r = p - {\rho}g_j(x_j-r_j) \end{align} 右辺の最終項は表面張力項(単相流には無い)。CSF (Continuum SUrface Force)モデルとして扱われる。 \begin{align} F_{CSF,i} = {\sigma}{\kappa}\cfrac{{\partial}F}{{\partial}x_i} \end{align} \(\, \sigma\)は表面張力係数。界面の曲率\(\,\kappa\)は、 \begin{align} \kappa = - {\nabla}{\cdot}{\bf{n}} = {\nabla}{\cdot} \left( {\cfrac{{\nabla}F}{|{\nabla}F|}} \right) \end{align} 混合密度\(\,\rho\)と混合粘度\(\,\mu\)は、VOF関数体積率\(\,F\)を用いて、 \begin{align} \rho = \rho_2 + F(\rho_1 - \rho_2) \end{align} \begin{align} \mu = \mu_2 + F(\mu_1 - \mu_2) \end{align}
interFoamで実装されているPIMPLE方での計算手順
(1) 流速の初期値、体積率\(\,F\)を設定し、セル境界面のフラックス\(\,\phi(u_i^0)\)を計算。
(2) VOF関数の移流方程式を解いて、体積率\(\,F\)を算出 (MULES or CMULES)。 \begin{align} \cfrac{{\partial}F}{{\partial}t} + \cfrac{{\partial}(F{u_j})}{{\partial}x_j} = 0 \end{align}
(3) 運動量保存式を解く。 \begin{align} A_P u_{iP}^{m*} + {\sum\limits_{l}^{NC}} A_l u_{il}^{m*} = Q_{iP}^{m-1} - \left( \cfrac{{\delta}p^{m-1}} {{\delta} x_i} \right)_P \Delta V \end{align} (注：上式を解かない場合、fvSolutionファイルのPIMPLE辞書のmomentunPredictorキーワードにfalseを指定する。)
(4) 運動量保存式から得られた流速\(\,u_i^{m*}\)を用いて暫定流速を獲得 ((3)を解かない場合は前タイムステップの\(\,u_i^n\))。 \begin{align} \tilde{u}_{iP}^{**} = \cfrac{1}{A_P} \left( - \sum\limits_{l}^{NC} A_l u_{il}^{m*} + Q_{iP}^{m-1} \right) \end{align}
(5) 暫定流速を用いて、Rhie-Chow補間により暫定境界面流束を計算。 \begin{align} \phi_{\alpha}(\tilde{u}_i^{**}) = {\langle \tilde{u}_i^{**} \rangle}_{\alpha} S_{i \alpha} + \left\langle \cfrac{\rho \Delta V}{A_P} \right\rangle _{\alpha} \left( \cfrac{\phi _\alpha (u_i^{n-1}) - {\langle u_i^{n-1} \rangle}_\alpha S_{i\alpha}}{\Delta t} \right) L \left( \phi(u_i^{n-1}) \right) - \left\langle \cfrac{\Delta V}{A_P} \right\rangle _{\alpha} \left( F_{CSF,i\alpha} - g_j (x_j - r_j) \cfrac{\delta \rho}{\delta x_i} \right)_\alpha S_{i \alpha} \end{align} ここで、制限関数\(\, L \)は、
\begin{align} L \left( \phi(u_i^{n-1}) \right) = 1 - \min \left( \cfrac{ \left\lvert \phi (u_i^{n-1}) - {\langle u_i^{n-1} \rangle} S_i \right\rvert \Delta t}{d_{CF} \left\lvert \bf{S} \right\rvert} , 1 \right) \end{align}
(6) 圧力方程式をlaplacianスキーム (反復計算)で解く。 \begin{align} \sum\limits_{\alpha}^{NS} \left\langle \cfrac{\Delta V}{A_P} \right\rangle _{\alpha} \left( \cfrac{ \delta p^m }{ \delta x_i } \right)_{\alpha} S_{i \alpha} = \sum\limits_{\alpha}^{NS} \phi_{\alpha}(\tilde{u}_i^{**}) \end{align}
(7) 次サイクルで必要となるセル境界面流束を計算する。 \begin{align} \phi _\alpha (u_i^m) = \left\langle \tilde{u}_i^{**} \right\rangle _\alpha S_{i\alpha} - \left\langle \cfrac{\Delta V}{A_P} \right\rangle _{\alpha} \left( \cfrac{ \delta p^m }{ \delta x_i } \right)_{\alpha} S_{i \alpha} \end{align}
(8) 流速を更新する。 \begin{align} u_i^m = \tilde{u}_i^{**} - \cfrac{\Delta V}{A_P} \left( F_{CSF,i} - g_j (x_j - r_j) \cfrac{\delta \rho}{\delta x_i} \right)_\alpha - \cfrac{\Delta V}{A_P} \cfrac{\delta p^m}{\delta x_i} \end{align}
(9) 圧力補正と流速補正を繰り返す(PISOループ)ため、 \(\,u_i^{m*} = u_i^m\) として、STEP(4)に戻り、STEP(8)まで繰り返す。
(10) 運動量保存式の係数行列を更新する外部反復ループを実行するために、STEP(2)に戻り、 STEP(9)まで反復計算する。収束したと判断された場合、 \(\,u_i^{n+1} = u_i^m\) として次の計算ステップに進む。
2 熱拡散方程式の追加
一連の実行コマンド (solver側)
$ run
$ cd myiPCF
$ wclean
$ vi createFields.H
$ touch TEqn.H
$ vi TEqn.H
$ vi myiPCF.C
$ wmake
追加内容: 熱拡散項・支配方程式の実装
熱拡散 (thermal diffusivity)を実装する。通常は\(\,\alpha\)で表されるがコーディングの都合上\(\, DT\)とする (\(\,\alpha\)はVOF値に取られてる)。
\(\, DT = \cfrac{k}{\rho c_p}\)
\(\, DT\): thermal diffusivity [m^2/S]
\(\, k\): thermal conductivity [W/(m*K)]
\(\, \rho\): density [kg/m^3]
\(\, c_p\): specific heat capacity [J/kg*K]
上式の内容を"createFields.H"の一番上 (おそらくどこでもよい)に書き足す。
createFields.H
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportPropertiesDict
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar DT
(
transportPropertiesDict.lookup("DT")
);
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr mixturePtr
(
phaseChangeTwoPhaseMixture::New(U, phi)
);
phaseChangeTwoPhaseMixture& mixture = mixturePtr();
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Construct incompressible turbulence model
autoPtr turbulence
(
incompressible::momentumTransportModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "createFvOptions.H"
2~17行はthermal diffusion (DTとして)の新しいtransport prppetyの記述。
18~30行は温度分布Tの生成。
支配方程式 (熱拡散方程式)
\begin{align}
\cfrac{{\partial}T}{{\partial}t} + {\nabla}{\cdot}({\phi}T) = {\nabla}{\cdot}(DT{\nabla}T)
\end{align}
上式をコードで表すと下記のようになる。新しく"TEqn.H"というファイルを作り、このコードを書き込む。
TEqn.H
{
fvScalarMatrix TEqn
(s
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
);
TEqn.solve();
}
上記ファイルを作っただけではだめなので、自作solverの"myiPCF.C"にて #include をする。おそらく、#include のタイミングが重要であり、下記の箇所に書き込む ("UEqn.H"の後なのは速度の項を流用するからか？)。
myiPCF.C
- 略 -
while (pimple.run(runTime))
{
- 略 -
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
- 略 -
}
divU = fvc::div(fvc::absolute(phi, U));
surfaceScalarField rhoPhi
(
- 略 -
);
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
mixture.correct();
#include "UEqn.H"
// ここに追加する↓
#include "TEqn.H
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
- 略 -
古いverのOpenFOAMは下記のように編集するっぽい。
flux variable \(\phi\)に関しては、 圧力のループ計算 (ポアソン方程式)のにおいて使用されている (それを流用するものだと思われる)。熱輸送は速度分布に依存するため、"Pressure-velocity PIMPLE corrector loop"のループ解析後に、TEqn.Hのファイルはincludeされる。phase changeはこのエネルギー方程式に依存し、"twoPhaseProperties->correct();:"(昔のOpenFOAMのコード)の前にincludeされなければならない。
一連の実行コマンド (case 側)
$ cp -r $FOAM_TUTORIALS/multiphase/interPhaseChangeFoam/cavitatingBullet $FOAM_RUN/temperatureTest
$ cd $FOAM_RUN/temperatureTest
$ vi Allrun
$ rm -r constant/triSurface system/snappyHexMeshDict
$ sed -i "/bullet/,/}/d" 0/*
$ vi constant/transportProperties
$ cp 0/p_rgh 0/T
$ vi 0/T
$ vi system/fvSchemes
$ vi system/fvSolution
$ vi sysetm/controlDict
$ ./Allrun
$ paraFoam
$ ffmpeg -r 30 -i 210919_cavitationBullet_myiPCF_T.%04d.jpeg -vcodec libx264 -pix_fmt yuv420p -r 60 210919_cavitationBullet_myiPCF_T.mp4
"cavitationBullet"のケースディレクトリをコピーしてくる。しかし、簡単に確認するために、bulletの部分は取り除いて自作solverのコンパイルの確認を行う。
まず、"Allrun"から、"potentialFoam"と"snappyHexMesh"の記述を取り除き下記のようにする。
Allrun
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Generate the base block mesh
runApplication blockMesh
# Run the solver
runApplication $(getApplication)
#------------------------------------------------------------------------------
次に、bulletのジオメトリやコードの記述を取り除く。
$ rm -r constant/triSurface system/snappyHexMeshDict
$ sed -i "/bullet/,/}/d" 0/*
メモ：
「sed」は「Stream EDitor」の略。「sed スクリプトコマンド ファイル名」で、指定したファイルをスクリプトコマンドに従って行単位に処理し、標準出力へ出力する。ファイル名を省略した場合は、標準入力からのデータを処理する。
「-i」 ファイルを直接編集する。
「s/置換前/置換後/」
［置換前］で指定した文字列にマッチした部分を［置換後］に置き換える。複数マッチした場合は先頭のみ置換、全てを置換したい場合は、「s/置換前/置換後/g」のように「g」オプションを指定する。
"transportProperties"に"DT"に関して定義する。
場所はどこでもよく、最後に追加した。
constant/transportProperties
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
phases (water vapour);
phaseChangeTwoPhaseMixture SchnerrSauer;
pSat 2300; // Saturation pressure
sigma 0.07;
water
{
transportModel Newtonian;
nu 9e-07;
rho 1000;
}
vapour
{
transportModel Newtonian;
nu 4.273e-04;
rho 0.02308;
}
KunzCoeffs
{
UInf U20.0;
tInf t0.005; // L = 0.1 m
Cc C1000;
Cv C1000;
}
MerkleCoeffs
{
UInf 20.0;
tInf 0.005; // L = 0.1 m
Cc 80;
Cv 1e-03;
}
SchnerrSauerCoeffs
{
n 1.6e+13;
dNuc 2.0e-06;
Cc 1;
Cv 1;
}
DT DT [0 2 -1 0 0 0 0] 0.002; // [m^2/s]
// ************************************************************************* //
"T"の設定。"0/p_rgh"をコピーして下記のように記述する。
0/T
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0]; // [K]
internalField uniform 363.15; // 90 deg. C
boundaryField
{
inlet
{
type fixedValue;
value uniform 383.15; // 110 deg. C
}
outlet
{
type fixedValue;
value $internalField;
}
walls
{
type symmetry;
}
}
// ************************************************************************* //
"system/fvSchemes"の設定。支配方程式を見ても明らかなように、"divSchemes"と"laplacianSchemes"の部分を修正する必要がある。
system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
interpolationSchemes
{
default linear;
}
divSchemes
{
default none;
div(rhoPhi,U) Gauss linearUpwind grad(U);
div(phi,omega) Gauss linearUpwind grad(omega);
div(phi,k) Gauss linearUpwind grad(k);
div(phi,alpha) Gauss interfaceCompression vanLeer 1;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
div(phi,T) Gauss upwind;
}
gradSchemes
{
default Gauss linear;
}
laplacianSchemes
{
default Gauss linear limited corrected 0.5;
laplacian(DT,T) Gauss linear corrected;
}
snGradSchemes
{
default limited corrected 0.5;
}
// ************************************************************************* //
"system/fvSolution"の設定。最初は"T"で設定していたがコンパイル通らず、"TFinal"にしたらうまくいった。solverはとりあえず"PBiCG"にしたがどれが良いかは不明。
system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
TFinal
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-7;
relTol 0;
};
"alpha.water.*"
{
nAlphaCorr 2;
nAlphaSubCycles 1;
MULESCorr yes;
nLimiterIter 5;
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-8;
relTol 0;
maxIter 10;
};
"U.*"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-6;
relTol 0;
};
p_rgh
{
solver GAMG;
tolerance 1e-8;
relTol 0.1;
smoother DICGaussSeidel;
maxIter 50;
};
p_rghFinal
{
solver PCG;
preconditioner
{
preconditioner GAMG;
tolerance 1e-6;
relTol 0;
nVcycles 2;
smoother DICGaussSeidel;
};
tolerance 1e-7;
relTol 0;
maxIter 50;
};
"pcorr.*"
{
$p_rgh;
relTol 0;
};
Phi
{
$p_rgh;
relTol 0;
};
}
potentialFlow
{
nNonOrthogonalCorrectors 3;
}
PIMPLE
{
momentumPredictor no;
nOuterCorrectors 1;
nCorrectors 3;
nNonOrthogonalCorrectors 0;
}
relaxationFactors
{
equations
{
".*" 1;
}
}
// ************************************************************************* //
"controlDict"で自作solverを設定することを忘れないように注意する。また、今回はbulletなしの簡単な解析のため、"endTime"は小さく、非定常解析のため"writeInterval"は細かく設定する。
system/controlDict
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application myiPCF; // 追加したところ
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 1e-3; // 追加したところ
deltaT 1e-8;
writeControl adjustableRunTime;
writeInterval 1e-5; // 追加したところ
purgeWrite 0;
writeFormat binary;
writePrecision 6;
writeCompression off;
timeFormat general;
runTimeModifiable yes;
adjustTimeStep on;
maxCo 5;
maxAlphaCo 2;
// ************************************************************************* //
Results
とりあえず熱拡散はうまくいってそう。
deltaT 1e-5; endTime 1e-3;
3 飽和蒸気圧曲線の追加 Psat(T)
August-Roche-Magnus formula
solver内にて飽和蒸気圧は２つのループで計算されている (\alpha -loopと pressure-loop)。
計算の不必要なとき (高速化や非線形方程式を再研鑽するとき)、温度場が更新されたとき飽和蒸気圧だけで計算され、圧力場が再構築される。ただし、必ずしもすべてのタイムステップにおいてそれが実行されるとは限らない。
$ touch calcPSatField.H
$ vi calcPSatField.H
calcPSatField.H
{
const dimensionedScalar t30_11("30.11", dimensionSet(0,0,0,1,0,0,0), 30.11);
const dimensionedScalar t273_15("273.15", dimensionSet(0,0,0,1,0,0,0), 273.15);
const dimensionedScalar t1("1", dimensionSet(0,0,0,1,0,0,0), 1);
const dimensionedScalar p610_94("610.94", dimensionSet(1,-1,-2,0,0,0,0), 610.94);
// dimensionSet( [kg], [m], [s], [K], [kg*mol], [A], [cd]), [kg/(m*S^2)]=[Pa]
// August-Roche-Magnus formula
pSat = p610_94 * exp( 17.625*(T-t273_15) / max(t1, T-t30_11) );
//max(1,...) is included to avoid problems with devision by 0
}
単位が必ず必要。ここではまだPsatを決めてない。
$ vi createFields.H
以下の一文を取り除く。
createFields.H
const dimensionedScalar& pSat = twoPhaseProperties->pSat();
そして、"volScalarField T"の後にincludeする文を付け加える。
"volScalarField p_rgh"の後。(飽和蒸気圧場の初期化のために"volScalarField p_rgh"が"p_rgh"として用いられるので。)
"phaseChangeTwoPhaseMixture"の前。(飽和蒸気圧が"phaseChangeTwoPhaseMixture"に用いられるので。)
createFields.H
volScalarField pSat
(
IOobject
(
"pSat",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh // initial value will be overwritten in calcPSatField.H
);
#include "calcPSatField.H"
温度"T"に依存する飽和蒸気圧を更新するには、"myiPCF.C"内に｢#include "TEqn.H"｣の後に、｢#include "calcPSatField.H"｣を記載する。
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportPropertiesDict
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar DT
(
transportPropertiesDict.lookup("DT")
);
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField pSat
(
IOobject
(
"pSat",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh // initial value will be overwritten in calcPSatField.H
);
#include "calcPSatField.H"
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr mixturePtr
(
phaseChangeTwoPhaseMixture::New(U, phi)
);
phaseChangeTwoPhaseMixture& mixture = mixturePtr();
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Construct incompressible turbulence model
autoPtr turbulence
(
incompressible::momentumTransportModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "createFvOptions.H"
/myiPCF/phaseChangeTwoPhaseMixtures
.
├── Kunz
│ ├── Kunz.C
│ └── Kunz.H
├── Make
│ ├── files
│ └── options
├── Merkle
│ ├── Merkle.C
│ └── Merkle.H
├── SchnerrSauer
│ ├── SchnerrSauer.C
│ └── SchnerrSauer.H
├── lnInclude
│ ├── Kunz.C -> ../Kunz/Kunz.C
│ ├── Kunz.H -> ../Kunz/Kunz.H
│ ├── Merkle.C -> ../Merkle/Merkle.C
│ ├── Merkle.H -> ../Merkle/Merkle.H
│ ├── SchnerrSauer.C -> ../SchnerrSauer/SchnerrSauer.C
│ ├── SchnerrSauer.H -> ../SchnerrSauer/SchnerrSauer.H
│ ├── phaseChangeTwoPhaseMixture.C -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
│ ├── phaseChangeTwoPhaseMixture.H -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
│ └── phaseChangeTwoPhaseMixtureNew.C -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
└── phaseChangeTwoPhaseMixture
├── phaseChangeTwoPhaseMixture.C
├── phaseChangeTwoPhaseMixture.H
└── phaseChangeTwoPhaseMixtureNew.C
6 directories, 20 files
stacic saturation pressure を取り除く。
$ vi phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
次の記述
//- Saturation vapour pressure
dimensionedScalar pSat_;
■ OpenFOAM8
修正前：
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
修正後：
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
return pSat_;
}
■ 古いOpenFOAM
修正前：
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
修正後：
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
return pSat;
}
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 New
(
const volVectorField& U,
const surfaceScalarField& phi
);
//- Destructor
virtual ~phaseChangeTwoPhaseMixture()
{}
// Member Functions
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
//- Return the mass condensation and vaporisation rates as a
// coefficient to multiply (1 - alphal) for the condensation rate
// and a coefficient to multiply alphal for the vaporisation rate
virtual Pair> mDotAlphal() const = 0;
//- Return the mass condensation and vaporisation rates as coefficients
// to multiply (p - pSat)
virtual Pair> mDotP() const = 0;
//- Return the volumetric condensation and vaporisation rates as a
// coefficient to multiply (1 - alphal) for the condensation rate
// and a coefficient to multiply alphal for the vaporisation rate
Pair> vDotAlphal() const;
//- Return the volumetric condensation and vaporisation rates as
// coefficients to multiply (p - pSat)
Pair> vDotP() const;
//- Correct the phaseChange model
virtual void correct() = 0;
//- Read the transportProperties dictionary and update
virtual bool read() = 0;
// Member Operators
//- Disallow default bitwise assignment
void operator=(const phaseChangeTwoPhaseMixture&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
pSatに関係するすべての記述を消去。
以下の3行目と2行目最後のコンマを消す。
$ vi phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
■ OpenFOAM8
修正前：
immiscibleIncompressibleTwoPhaseMixture(U, phi),
phaseChangeTwoPhaseMixtureCoeffs_(optionalSubDict(type + "Coeffs")),
pSat_("pSat", dimPressure, lookup("pSat"))
{}
修正後：
immiscibleIncompressibleTwoPhaseMixture(U, phi),
phaseChangeTwoPhaseMixtureCoeffs_(optionalSubDict(type + "Coeffs"))
{}
■ 古いOpenFOAM
修正前：
twoPhaseMixture(U, phi, alpha1Name),
phaseChangeTwoPhaseMixtureCoeffs_(subDict(type + "Coeffs")),
pSat_(lookup("pSat"))
{}
修正後：
twoPhaseMixture(U, phi, alpha1Name),
phaseChangeTwoPhaseMixtureCoeffs_(subDict(type + "Coeffs"))
{}
次の記述を消す。
lookup("pSat") >> pSat_;
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 >
Foam::phaseChangeTwoPhaseMixture::vDotAlphal() const
{
volScalarField alphalCoeff(1.0/rho1() - alpha1()*(1.0/rho1() - 1.0/rho2()));
Pair> mDotAlphal = this->mDotAlphal();
return Pair>
(
alphalCoeff*mDotAlphal[0],
alphalCoeff*mDotAlphal[1]
);
}
Foam::Pair>
Foam::phaseChangeTwoPhaseMixture::vDotP() const
{
dimensionedScalar pCoeff(1.0/rho1() - 1.0/rho2());
Pair> mDotP = this->mDotP();
return Pair>(pCoeff*mDotP[0], pCoeff*mDotP[1]);
}
void Foam::phaseChangeTwoPhaseMixture::correct()
{
immiscibleIncompressibleTwoPhaseMixture::correct();
}
bool Foam::phaseChangeTwoPhaseMixture::read()
{
if (immiscibleIncompressibleTwoPhaseMixture::read())
{
phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
lookup("pSat") >> pSat_;
return true;
}
else
{
return false;
}
}
// ************************************************************************* //
以上で修正終了。
$ wclean
$ wmake
躓いたエラー
エラーログ
Making dependency list for source file myiPCF.C
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -c myiPCF.C -o Make/linux64GccDPInt32Opt/myiPCF.o
In file included from myiPCF.C:49:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H: In member function ‘const volScalarField& Foam::phaseChangeTwoPhaseMixture::pSat() const’:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:124:42: error: reference to ‘alpha1_’ is ambiguous
124 | const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
| ^~~~~~~
In file included from myiPCF.C:48:
/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude/interfaceProperties.H:69:31: note: candidates are: ‘const volScalarField& Foam::interfaceProperties::alpha1_’
69 | const volScalarField& alpha1_;
| ^~~~~~~
In file included from /opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude/incompressibleTwoPhaseMixture.H:40,
from /opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude/immiscibleIncompressibleTwoPhaseMixture.H:38,
from phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:38,
from myiPCF.C:49:
/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude/twoPhaseMixture.H:58:24: note: ‘Foam::volScalarField Foam::twoPhaseMixture::alpha1_’
58 | volScalarField alpha1_;
| ^~~~~~~
In file included from myiPCF.C:49:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:124:82: error: expected primary-expression before ‘>’ token
124 | const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
| ^
make: *** [/opt/openfoam8/wmake/rules/General/transform:26: Make/linux64GccDPInt32Opt/myiPCF.o] エラー 1
たぶん、‘alpha1_’ の定義がはっきりしてない。
"Foam::interfaceProperties::alpha1_"や"Foam::twoPhaseMixture::alpha1_"が候補として出ているが、ここではひとまず、"Foam::twoPhaseMixture::alpha1_"とした。
修正内容
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = Foam::twoPhaseMixture::alpha1_.db().lookupObject("pSat");
return pSat;
}
以上で修正終了。
$ wclean
$ wmake
コンパイルログ
tokudy@tokudy-desktop:~/OpenFOAM/tokudy-8/run/myiPCFPsat$ wmake
Making dependency list for source file myiPCF.C
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -c myiPCF.C -o Make/linux64GccDPInt32Opt/myiPCF.o
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -fuse-ld=bfd -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPInt32Opt/myiPCF.o -L/opt/openfoam8/platforms/linux64GccDPInt32Opt/lib \
-lphaseChangeTwoPhaseMixtures -ltwoPhaseMixture -linterfaceProperties -ltwoPhaseProperties -lincompressibleTransportModels -lmomentumTransportModels -lincompressibleMomentumTransportModels -lfiniteVolume -ldynamicFvMesh -ltopoChangerFvMesh -lfvOptions -lmeshTools -lOpenFOAM -ldl \
-lm -o /opt/openfoam8/platforms/linux64GccDPInt32Opt/bin/myiPCF
コンパイル通っているぽい。おｋ．
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportPropertiesDict
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar DT
(
transportPropertiesDict.lookup("DT")
);
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr mixturePtr
(
phaseChangeTwoPhaseMixture::New(U, phi)
);
phaseChangeTwoPhaseMixture& mixture = mixturePtr();
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Construct incompressible turbulence model
autoPtr turbulence
(
incompressible::momentumTransportModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "createFvOptions.H"
2~17行はthermal diffusion (DTとして)の新しいtransport prppetyの記述。
18~30行は温度分布Tの生成。
TEqn.H
上記ファイルを作っただけではだめなので、自作solverの"myiPCF.C"にて #include をする。おそらく、#include のタイミングが重要であり、下記の箇所に書き込む ("UEqn.H"の後なのは速度の項を流用するからか？)。
myiPCF.C
古いverのOpenFOAMは下記のように編集するっぽい。
"cavitationBullet"のケースディレクトリをコピーしてくる。しかし、簡単に確認するために、bulletの部分は取り除いて自作solverのコンパイルの確認を行う。 まず、"Allrun"から、"potentialFoam"と"snappyHexMesh"の記述を取り除き下記のようにする。
Allrun
次に、bulletのジオメトリやコードの記述を取り除く。
メモ：
"transportProperties"に"DT"に関して定義する。 場所はどこでもよく、最後に追加した。
constant/transportProperties
"T"の設定。"0/p_rgh"をコピーして下記のように記述する。
0/T
"system/fvSchemes"の設定。支配方程式を見ても明らかなように、"divSchemes"と"laplacianSchemes"の部分を修正する必要がある。
system/fvSchemes
"system/fvSolution"の設定。最初は"T"で設定していたがコンパイル通らず、"TFinal"にしたらうまくいった。solverはとりあえず"PBiCG"にしたがどれが良いかは不明。
system/fvSchemes
支配方程式 (熱拡散方程式)
\begin{align} \cfrac{{\partial}T}{{\partial}t} + {\nabla}{\cdot}({\phi}T) = {\nabla}{\cdot}(DT{\nabla}T) \end{align} 上式をコードで表すと下記のようになる。新しく"TEqn.H"というファイルを作り、このコードを書き込む。TEqn.H
{
fvScalarMatrix TEqn
(s
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
);
TEqn.solve();
}
上記ファイルを作っただけではだめなので、自作solverの"myiPCF.C"にて #include をする。おそらく、#include のタイミングが重要であり、下記の箇所に書き込む ("UEqn.H"の後なのは速度の項を流用するからか？)。
myiPCF.C
- 略 -
while (pimple.run(runTime))
{
- 略 -
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
- 略 -
}
divU = fvc::div(fvc::absolute(phi, U));
surfaceScalarField rhoPhi
(
- 略 -
);
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
mixture.correct();
#include "UEqn.H"
// ここに追加する↓
#include "TEqn.H
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
- 略 -
古いverのOpenFOAMは下記のように編集するっぽい。
flux variable \(\phi\)に関しては、 圧力のループ計算 (ポアソン方程式)のにおいて使用されている (それを流用するものだと思われる)。熱輸送は速度分布に依存するため、"Pressure-velocity PIMPLE corrector loop"のループ解析後に、TEqn.Hのファイルはincludeされる。phase changeはこのエネルギー方程式に依存し、"twoPhaseProperties->correct();:"(昔のOpenFOAMのコード)の前にincludeされなければならない。
一連の実行コマンド (case 側)
$ cp -r $FOAM_TUTORIALS/multiphase/interPhaseChangeFoam/cavitatingBullet $FOAM_RUN/temperatureTest $ cd $FOAM_RUN/temperatureTest $ vi Allrun $ rm -r constant/triSurface system/snappyHexMeshDict $ sed -i "/bullet/,/}/d" 0/* $ vi constant/transportProperties $ cp 0/p_rgh 0/T $ vi 0/T $ vi system/fvSchemes $ vi system/fvSolution $ vi sysetm/controlDict $ ./Allrun $ paraFoam $ ffmpeg -r 30 -i 210919_cavitationBullet_myiPCF_T.%04d.jpeg -vcodec libx264 -pix_fmt yuv420p -r 60 210919_cavitationBullet_myiPCF_T.mp4
"cavitationBullet"のケースディレクトリをコピーしてくる。しかし、簡単に確認するために、bulletの部分は取り除いて自作solverのコンパイルの確認を行う。 まず、"Allrun"から、"potentialFoam"と"snappyHexMesh"の記述を取り除き下記のようにする。
Allrun
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Generate the base block mesh
runApplication blockMesh
# Run the solver
runApplication $(getApplication)
#------------------------------------------------------------------------------
次に、bulletのジオメトリやコードの記述を取り除く。
$ rm -r constant/triSurface system/snappyHexMeshDict $ sed -i "/bullet/,/}/d" 0/*
メモ：
「sed」は「Stream EDitor」の略。「sed スクリプトコマンド ファイル名」で、指定したファイルをスクリプトコマンドに従って行単位に処理し、標準出力へ出力する。ファイル名を省略した場合は、標準入力からのデータを処理する。 「-i」 ファイルを直接編集する。 「s/置換前/置換後/」 ［置換前］で指定した文字列にマッチした部分を［置換後］に置き換える。複数マッチした場合は先頭のみ置換、全てを置換したい場合は、「s/置換前/置換後/g」のように「g」オプションを指定する。
"transportProperties"に"DT"に関して定義する。 場所はどこでもよく、最後に追加した。
constant/transportProperties
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
phases (water vapour);
phaseChangeTwoPhaseMixture SchnerrSauer;
pSat 2300; // Saturation pressure
sigma 0.07;
water
{
transportModel Newtonian;
nu 9e-07;
rho 1000;
}
vapour
{
transportModel Newtonian;
nu 4.273e-04;
rho 0.02308;
}
KunzCoeffs
{
UInf U20.0;
tInf t0.005; // L = 0.1 m
Cc C1000;
Cv C1000;
}
MerkleCoeffs
{
UInf 20.0;
tInf 0.005; // L = 0.1 m
Cc 80;
Cv 1e-03;
}
SchnerrSauerCoeffs
{
n 1.6e+13;
dNuc 2.0e-06;
Cc 1;
Cv 1;
}
DT DT [0 2 -1 0 0 0 0] 0.002; // [m^2/s]
// ************************************************************************* //
"T"の設定。"0/p_rgh"をコピーして下記のように記述する。
0/T
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0]; // [K]
internalField uniform 363.15; // 90 deg. C
boundaryField
{
inlet
{
type fixedValue;
value uniform 383.15; // 110 deg. C
}
outlet
{
type fixedValue;
value $internalField;
}
walls
{
type symmetry;
}
}
// ************************************************************************* //
"system/fvSchemes"の設定。支配方程式を見ても明らかなように、"divSchemes"と"laplacianSchemes"の部分を修正する必要がある。
system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
interpolationSchemes
{
default linear;
}
divSchemes
{
default none;
div(rhoPhi,U) Gauss linearUpwind grad(U);
div(phi,omega) Gauss linearUpwind grad(omega);
div(phi,k) Gauss linearUpwind grad(k);
div(phi,alpha) Gauss interfaceCompression vanLeer 1;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
div(phi,T) Gauss upwind;
}
gradSchemes
{
default Gauss linear;
}
laplacianSchemes
{
default Gauss linear limited corrected 0.5;
laplacian(DT,T) Gauss linear corrected;
}
snGradSchemes
{
default limited corrected 0.5;
}
// ************************************************************************* //
"system/fvSolution"の設定。最初は"T"で設定していたがコンパイル通らず、"TFinal"にしたらうまくいった。solverはとりあえず"PBiCG"にしたがどれが良いかは不明。
system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
TFinal
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-7;
relTol 0;
};
"alpha.water.*"
{
nAlphaCorr 2;
nAlphaSubCycles 1;
MULESCorr yes;
nLimiterIter 5;
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-8;
relTol 0;
maxIter 10;
};
"U.*"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-6;
relTol 0;
};
p_rgh
{
solver GAMG;
tolerance 1e-8;
relTol 0.1;
smoother DICGaussSeidel;
maxIter 50;
};
p_rghFinal
{
solver PCG;
preconditioner
{
preconditioner GAMG;
tolerance 1e-6;
relTol 0;
nVcycles 2;
smoother DICGaussSeidel;
};
tolerance 1e-7;
relTol 0;
maxIter 50;
};
"pcorr.*"
{
$p_rgh;
relTol 0;
};
Phi
{
$p_rgh;
relTol 0;
};
}
potentialFlow
{
nNonOrthogonalCorrectors 3;
}
PIMPLE
{
momentumPredictor no;
nOuterCorrectors 1;
nCorrectors 3;
nNonOrthogonalCorrectors 0;
}
relaxationFactors
{
equations
{
".*" 1;
}
}
// ************************************************************************* //
"controlDict"で自作solverを設定することを忘れないように注意する。また、今回はbulletなしの簡単な解析のため、"endTime"は小さく、非定常解析のため"writeInterval"は細かく設定する。
system/controlDict
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application myiPCF; // 追加したところ
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 1e-3; // 追加したところ
deltaT 1e-8;
writeControl adjustableRunTime;
writeInterval 1e-5; // 追加したところ
purgeWrite 0;
writeFormat binary;
writePrecision 6;
writeCompression off;
timeFormat general;
runTimeModifiable yes;
adjustTimeStep on;
maxCo 5;
maxAlphaCo 2;
// ************************************************************************* //
Results
とりあえず熱拡散はうまくいってそう。
deltaT 1e-5; endTime 1e-3;
3 飽和蒸気圧曲線の追加 Psat(T)
August-Roche-Magnus formula
solver内にて飽和蒸気圧は２つのループで計算されている (\alpha -loopと pressure-loop)。
計算の不必要なとき (高速化や非線形方程式を再研鑽するとき)、温度場が更新されたとき飽和蒸気圧だけで計算され、圧力場が再構築される。ただし、必ずしもすべてのタイムステップにおいてそれが実行されるとは限らない。
$ touch calcPSatField.H
$ vi calcPSatField.H
calcPSatField.H
{
const dimensionedScalar t30_11("30.11", dimensionSet(0,0,0,1,0,0,0), 30.11);
const dimensionedScalar t273_15("273.15", dimensionSet(0,0,0,1,0,0,0), 273.15);
const dimensionedScalar t1("1", dimensionSet(0,0,0,1,0,0,0), 1);
const dimensionedScalar p610_94("610.94", dimensionSet(1,-1,-2,0,0,0,0), 610.94);
// dimensionSet( [kg], [m], [s], [K], [kg*mol], [A], [cd]), [kg/(m*S^2)]=[Pa]
// August-Roche-Magnus formula
pSat = p610_94 * exp( 17.625*(T-t273_15) / max(t1, T-t30_11) );
//max(1,...) is included to avoid problems with devision by 0
}
単位が必ず必要。ここではまだPsatを決めてない。
$ vi createFields.H
以下の一文を取り除く。
createFields.H
const dimensionedScalar& pSat = twoPhaseProperties->pSat();
そして、"volScalarField T"の後にincludeする文を付け加える。
"volScalarField p_rgh"の後。(飽和蒸気圧場の初期化のために"volScalarField p_rgh"が"p_rgh"として用いられるので。)
"phaseChangeTwoPhaseMixture"の前。(飽和蒸気圧が"phaseChangeTwoPhaseMixture"に用いられるので。)
createFields.H
volScalarField pSat
(
IOobject
(
"pSat",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh // initial value will be overwritten in calcPSatField.H
);
#include "calcPSatField.H"
温度"T"に依存する飽和蒸気圧を更新するには、"myiPCF.C"内に｢#include "TEqn.H"｣の後に、｢#include "calcPSatField.H"｣を記載する。
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportPropertiesDict
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar DT
(
transportPropertiesDict.lookup("DT")
);
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField pSat
(
IOobject
(
"pSat",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh // initial value will be overwritten in calcPSatField.H
);
#include "calcPSatField.H"
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr mixturePtr
(
phaseChangeTwoPhaseMixture::New(U, phi)
);
phaseChangeTwoPhaseMixture& mixture = mixturePtr();
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Construct incompressible turbulence model
autoPtr turbulence
(
incompressible::momentumTransportModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "createFvOptions.H"
/myiPCF/phaseChangeTwoPhaseMixtures
.
├── Kunz
│ ├── Kunz.C
│ └── Kunz.H
├── Make
│ ├── files
│ └── options
├── Merkle
│ ├── Merkle.C
│ └── Merkle.H
├── SchnerrSauer
│ ├── SchnerrSauer.C
│ └── SchnerrSauer.H
├── lnInclude
│ ├── Kunz.C -> ../Kunz/Kunz.C
│ ├── Kunz.H -> ../Kunz/Kunz.H
│ ├── Merkle.C -> ../Merkle/Merkle.C
│ ├── Merkle.H -> ../Merkle/Merkle.H
│ ├── SchnerrSauer.C -> ../SchnerrSauer/SchnerrSauer.C
│ ├── SchnerrSauer.H -> ../SchnerrSauer/SchnerrSauer.H
│ ├── phaseChangeTwoPhaseMixture.C -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
│ ├── phaseChangeTwoPhaseMixture.H -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
│ └── phaseChangeTwoPhaseMixtureNew.C -> ../phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
└── phaseChangeTwoPhaseMixture
├── phaseChangeTwoPhaseMixture.C
├── phaseChangeTwoPhaseMixture.H
└── phaseChangeTwoPhaseMixtureNew.C
6 directories, 20 files
stacic saturation pressure を取り除く。
$ vi phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
次の記述
//- Saturation vapour pressure
dimensionedScalar pSat_;
■ OpenFOAM8
修正前：
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
修正後：
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
return pSat_;
}
■ 古いOpenFOAM
修正前：
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
修正後：
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
return pSat;
}
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 New
(
const volVectorField& U,
const surfaceScalarField& phi
);
//- Destructor
virtual ~phaseChangeTwoPhaseMixture()
{}
// Member Functions
//- Return const-access to the saturation vapour pressure
const dimensionedScalar& pSat() const
{
return pSat_;
}
//- Return the mass condensation and vaporisation rates as a
// coefficient to multiply (1 - alphal) for the condensation rate
// and a coefficient to multiply alphal for the vaporisation rate
virtual Pair> mDotAlphal() const = 0;
//- Return the mass condensation and vaporisation rates as coefficients
// to multiply (p - pSat)
virtual Pair> mDotP() const = 0;
//- Return the volumetric condensation and vaporisation rates as a
// coefficient to multiply (1 - alphal) for the condensation rate
// and a coefficient to multiply alphal for the vaporisation rate
Pair> vDotAlphal() const;
//- Return the volumetric condensation and vaporisation rates as
// coefficients to multiply (p - pSat)
Pair> vDotP() const;
//- Correct the phaseChange model
virtual void correct() = 0;
//- Read the transportProperties dictionary and update
virtual bool read() = 0;
// Member Operators
//- Disallow default bitwise assignment
void operator=(const phaseChangeTwoPhaseMixture&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
pSatに関係するすべての記述を消去。
以下の3行目と2行目最後のコンマを消す。
$ vi phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
■ OpenFOAM8
修正前：
immiscibleIncompressibleTwoPhaseMixture(U, phi),
phaseChangeTwoPhaseMixtureCoeffs_(optionalSubDict(type + "Coeffs")),
pSat_("pSat", dimPressure, lookup("pSat"))
{}
修正後：
immiscibleIncompressibleTwoPhaseMixture(U, phi),
phaseChangeTwoPhaseMixtureCoeffs_(optionalSubDict(type + "Coeffs"))
{}
■ 古いOpenFOAM
修正前：
twoPhaseMixture(U, phi, alpha1Name),
phaseChangeTwoPhaseMixtureCoeffs_(subDict(type + "Coeffs")),
pSat_(lookup("pSat"))
{}
修正後：
twoPhaseMixture(U, phi, alpha1Name),
phaseChangeTwoPhaseMixtureCoeffs_(subDict(type + "Coeffs"))
{}
次の記述を消す。
lookup("pSat") >> pSat_;
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 >
Foam::phaseChangeTwoPhaseMixture::vDotAlphal() const
{
volScalarField alphalCoeff(1.0/rho1() - alpha1()*(1.0/rho1() - 1.0/rho2()));
Pair> mDotAlphal = this->mDotAlphal();
return Pair>
(
alphalCoeff*mDotAlphal[0],
alphalCoeff*mDotAlphal[1]
);
}
Foam::Pair>
Foam::phaseChangeTwoPhaseMixture::vDotP() const
{
dimensionedScalar pCoeff(1.0/rho1() - 1.0/rho2());
Pair> mDotP = this->mDotP();
return Pair>(pCoeff*mDotP[0], pCoeff*mDotP[1]);
}
void Foam::phaseChangeTwoPhaseMixture::correct()
{
immiscibleIncompressibleTwoPhaseMixture::correct();
}
bool Foam::phaseChangeTwoPhaseMixture::read()
{
if (immiscibleIncompressibleTwoPhaseMixture::read())
{
phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
lookup("pSat") >> pSat_;
return true;
}
else
{
return false;
}
}
// ************************************************************************* //
以上で修正終了。
$ wclean
$ wmake
躓いたエラー
エラーログ
Making dependency list for source file myiPCF.C
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -c myiPCF.C -o Make/linux64GccDPInt32Opt/myiPCF.o
In file included from myiPCF.C:49:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H: In member function ‘const volScalarField& Foam::phaseChangeTwoPhaseMixture::pSat() const’:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:124:42: error: reference to ‘alpha1_’ is ambiguous
124 | const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
| ^~~~~~~
In file included from myiPCF.C:48:
/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude/interfaceProperties.H:69:31: note: candidates are: ‘const volScalarField& Foam::interfaceProperties::alpha1_’
69 | const volScalarField& alpha1_;
| ^~~~~~~
In file included from /opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude/incompressibleTwoPhaseMixture.H:40,
from /opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude/immiscibleIncompressibleTwoPhaseMixture.H:38,
from phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:38,
from myiPCF.C:49:
/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude/twoPhaseMixture.H:58:24: note: ‘Foam::volScalarField Foam::twoPhaseMixture::alpha1_’
58 | volScalarField alpha1_;
| ^~~~~~~
In file included from myiPCF.C:49:
phaseChangeTwoPhaseMixtures/lnInclude/phaseChangeTwoPhaseMixture.H:124:82: error: expected primary-expression before ‘>’ token
124 | const volScalarField& pSat = alpha1_.db().lookupObject("pSat");
| ^
make: *** [/opt/openfoam8/wmake/rules/General/transform:26: Make/linux64GccDPInt32Opt/myiPCF.o] エラー 1
たぶん、‘alpha1_’ の定義がはっきりしてない。
"Foam::interfaceProperties::alpha1_"や"Foam::twoPhaseMixture::alpha1_"が候補として出ているが、ここではひとまず、"Foam::twoPhaseMixture::alpha1_"とした。
修正内容
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
//- Return const-access to the saturation vapour pressure
const volScalarField& pSat() const
{
const volScalarField& pSat = Foam::twoPhaseMixture::alpha1_.db().lookupObject("pSat");
return pSat;
}
以上で修正終了。
$ wclean
$ wmake
コンパイルログ
tokudy@tokudy-desktop:~/OpenFOAM/tokudy-8/run/myiPCFPsat$ wmake
Making dependency list for source file myiPCF.C
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -c myiPCF.C -o Make/linux64GccDPInt32Opt/myiPCF.o
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3 -DNoRepository -ftemplate-depth-100 -I. -I/opt/openfoam8/src/transportModels/lnInclude -I/opt/openfoam8/src/twoPhaseModels/twoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/interfaceProperties/lnInclude -I/opt/openfoam8/src/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude -I/opt/openfoam8/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude -IphaseChangeTwoPhaseMixtures/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/momentumTransportModels/lnInclude -I/opt/openfoam8/src/MomentumTransportModels/incompressible/lnInclude -I/opt/openfoam8/src/finiteVolume/lnInclude -I/opt/openfoam8/src/dynamicFvMesh/lnInclude -I/opt/openfoam8/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam8/src/OpenFOAM/lnInclude -I/opt/openfoam8/src/OSspecific/POSIX/lnInclude -fPIC -fuse-ld=bfd -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPInt32Opt/myiPCF.o -L/opt/openfoam8/platforms/linux64GccDPInt32Opt/lib \
-lphaseChangeTwoPhaseMixtures -ltwoPhaseMixture -linterfaceProperties -ltwoPhaseProperties -lincompressibleTransportModels -lmomentumTransportModels -lincompressibleMomentumTransportModels -lfiniteVolume -ldynamicFvMesh -ltopoChangerFvMesh -lfvOptions -lmeshTools -lOpenFOAM -ldl \
-lm -o /opt/openfoam8/platforms/linux64GccDPInt32Opt/bin/myiPCF
コンパイル通っているぽい。おｋ．