OpenFOAMで実装されている非圧縮粘性流れの非定常ソルバーであるpisoFoamについて解説する
CFDの基礎方程式
密度\(\rho\)，粘性係数\(\mu\)が一定な非圧縮性流体の基礎方程式は，連続の式（質量保存式）およびNavier-Stokes方程式（運動量保存式）である
（連続の式）
\begin{equation}
\frac{\partial u_{i}}{\partial x_{i}}=0
\end{equation}
（Navier-Stokes方程式）
\begin{equation}
\frac{\partial u_{i}}{\partial t}+\frac{\partial \left(u_{i}u_{j}\right)}{\partial x_{j}}=-\frac{1}{\rho}\frac{\partial P}{\partial x_{i}}+\nu_{e}\frac{\partial}{\partial x_{j}}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)
+f_{i}
\end{equation}
\begin{equation}
\nu_{e}=\nu+\nu_{t}, \qquad {P}={p}+\frac{2}{3}\rho k
\end{equation}
※モデル化のためのフィルター操作を表すオーバーライン（\(\overline{u}_{i}\)，\(\overline{p}\)）は省略している
密度\(\rho\)と動粘性係数\(\nu\)は水や空気などの物性値であり，渦粘性係数\(\nu_{t}\)は乱流の影響による拡散を表している
未知数は速度\(\mathbf{u}\)の3つの成分\(u_{i}\)，圧力\(p\)，渦粘性係数\(\nu_{t}\)の5つであり，渦粘性係数を乱流モデルによって与えれば，連続の式とNavier-Stokes方程式を解くことですべての物理量を求めることができる
また，これ以降は簡単のため \(P \rightarrow p\) と表記することとする
圧力-速度連成
連続の式とNavier-Stokes方程式には次のような問題がある
- 連続の式には時間微分項が含まれておらず，速度場に対して強い制約（発散なし：divergence free）を課しているだけである
- 変数である圧力\(p\)に対する時間微分項\(\frac{\partial p}{\partial t}\)が明示的に示されていない
流れ場の速度と圧力を計算するためには，基礎方程式で明示的に示されていない圧力\(p\)を計算しつつ，連続の式による強い制約を満たすような速度\(u_{i}\)を計算する，というような手法を考える必要がある
PISO法
PISO（Pressure Implicit with Splitting Operators）法は，非定常計算が行えるようにpiso法を拡張した手法である
piso法の導出において，速度と圧力の修正量を求める際に\(\mathbf{u'}_{nb}\)の項を無視したが，それが原因でpiso法では収束までの時々刻々における速度と圧力の値は不正確なものになっている
上記の問題を解決するために，PISO法では速度と圧力の推定値と修正量の計算を複数回（通常は2回）行うことで，時々刻々における速度と圧力の精度を確保して非定常流れの計算を可能にしている
時刻 \(t_{n}\) における速度の推定値は次の式で求められる
\begin{equation}
\frac{1}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{*}
+\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}
=-\nabla p^{(n-1)}+\mathbf{s}_{_{P}}+\frac{1-\alpha_{u}}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{(n-1)}
\tag{1}
\end{equation}
ここで，\(\mathbf{\tilde{u}}^{*}_{_{P}}\)を次のように定義する
\begin{equation}
\mathbf{\tilde{u}}^{*}_{_{P}}
=\frac{\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}}{\mathbf{A}_{_{P}}}
\tag{2}
\end{equation}
時刻\(t_{n}\)における圧力\(p^{(n)}\)と速度\(\mathbf{u}^{(n)}\)は次のように求められる
\begin{align}
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n)}\right)
=\nabla\cdot\mathbf{\tilde{u}}^{*}_{_{P}}
\tag{3} \\ \\
&\mathbf{u}^{(n)}_{_{P}}=\mathbf{\tilde{u}}^{*}_{_{P}}-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n)}
\tag{4}
\end{align}
\(\alpha_{u}\)は不足緩和係数であり，計算を安定させて収束を早めるはたらきをしている
PISO法では圧力の不足緩和は必要ない
以上の式を用いて，PISO法では次に示す①～⑦の手順を収束するまで繰り返すことによって真の速度場と圧力を計算する
① 時刻\(t_{n}\)において，\(\mathbf{u}^{(n-1)}\) と \(p^{(n-1)}\) から \(\nu_{t}\) を計算し，偏微分方程式を離散化して係数\(\mathbf{A}_{_{P}}\)と\(\mathbf{A}_{nb}\)を求める
② \(\mathbf{u}^{(n-1)}\) と \(p^{(n-1)}\) から時刻 \(t_{n}\) における速度の推定値 \(\mathbf{u}^{*}\) を求める (1)
③ 速度の推定値\(\mathbf{u}_{_{P}}^{*}\)から\(\mathbf{\tilde{u}}^{*}_{_{P}}\)を求める (2)
④ \(\mathbf{\tilde{u}}^{*}_{_{P}}\)から時刻\(t_{n}\)における圧力\(p^{(n)}\)を求める (3)
⑤ \(\mathbf{\tilde{u}}^{*}_{_{P}}\)と圧力\(p^{(n)}\)から時刻\(t_{n}\)における速度\(\mathbf{u}^{(n)}\)を求める (4)
⑥ ⑤で得られた\(\mathbf{u}^{(n)}_{_{P}}\)を\(\mathbf{u}^{*}_{_{P}}\)として，③に戻る
⑦ 時刻を\(t_{n+1}\)に進めて，①に戻る
不足緩和係数は，通常 \(\alpha_{u}=0.3 \sim 0.8\) の値をとり，試行錯誤的に決定される
係数を小さくすると計算は安定するが収束が遅くなるので，まずは小さい値から初めて，計算が安定するぎりぎりまで大きくしていくとよい
ちなみに，(3)と(4)より
\begin{equation}
\nabla\cdot\left(\mathbf{\tilde{u}}^{*}_{_{P}}-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n)}\right)
=\nabla\cdot\mathbf{u}^{(n)}_{_{P}}=0
\tag{A}
\end{equation}
となり，\(\mathbf{u}^{(n)}_{_{P}}\) は連続の式を満たしている
↓PISO法の詳しい導出はこちら
OpenFOAMにおける実装の概要
OpenFOAMでは，次の1~11のステップによってPISO法を実装している
- 次の時間ステップ \(t=t_{n+1}\) に進む
- 既知の\(\mathbf{u}^{(n)}\)と\(p^{(n)}\)を用いて\(\mathbf{u}^{(n+1)}\)と\(p^{(n+1)}\)を初期化する
- 運動量保存則（Navier-Stokes方程式）を立てる
- 運動量保存則の不足緩和を行う
- 運動量保存則を解いて速度の予測値\(\mathbf{u}^{*}\)を求める
- ポアソン方程式を立てる
- ポアソン方程式を解いて\(p^{n+1}\)を求める
- 質量流束\(\phi_{n+1}\)を修正する
- 速度\(\mathbf{n}^{(n+1)}\)の修正する
- 指定した回数だけステップ6に戻る
- 設定した計算終了時間に達していなければステップ2に戻る
PISO法は主に次の4つのファイルに分けて記述されている
- pisoFoam.C
- UEqn.H
- pEqn.H
- createFields.H
上記の11のステップは，それぞれ次のファイルに記述されている
変数の定義 → createFields.H
1，2，10, 11 → pisoFoam.C
3，4，5 → UEqn.C
6，7，8，9 → pEqn.H
ソースコードはほぼpisoFoamと同じだが，1つずつ説明していく
createFields.H
圧力\(p\)と速度\(\mathbf{u}\)の宣言などが行われている
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
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"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, piso.dict(), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
#include "createMRF.H"
#include "createFvOptions.H"
createPhi.H では質量流束が定義されている
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::flux(U)
);
// ************************************************************************* //32~35行目では，fvsolution で定義した圧力の参照点や参照値が読み込まれている
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}tutorials/incompressible/pisoFoam/RAS/cavity/system/fvSolution
40~43行目で，turbulenceProperties で定義した乱流モデルが読み込まれている
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2106 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType RAS;
RAS
{
RASModel kEpsilon;
turbulence on;
printCoeffs on;
}
// ************************************************************************* //$FOAM_TUTORIAL/incompressible/pisoFoam/windAroundBuildings/constant/turbulenceProperties
pisoFoam.C
mainにあたるプログラムで，while ( runtime.loop() ){ } でpiso法の時間発展の計算が行われている
1回の時間ステップの中では UEqn.H → PISO loop → turbulence->correct() が呼び出され，while ( piso.correct() ){ #include "pEqn.H" } で，ユーザーが nCorrectors で指定した回数だけ速度と圧力の修正を行う（通常は2回）
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
-------------------------------------------------------------------------------
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
#include "fvOptions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"Transient solver for incompressible, turbulent flow,"
" using the PISO algorithm."
);
#include "postProcess.H"
#include "addCheckCaseOptions.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "initContinuityErrs.H"
turbulence->validate();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "CourantNo.H"
// Pressure-velocity PISO corrector
{
#include "UEqn.H"
// --- PISO loop
while (piso.correct())
{
#include "pEqn.H"
}
}
laminarTransport.correct();
turbulence->correct();
runTime.write();
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
PISO loop の回数は fvSolution の中で nCorrectors として指定する
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
}tutorials/incompressible/pisoFoam/LES/pitzDaily/system/fvSolution
laminarTransport.correct() は何をしているのかよくわからない
runTime.write() でユーザーが IOobject::AUTO_WRITE とした変数が書き出される
runTime.printExecutionTime(Info) は log.pisofoam に計算時間を書き出している
Time = 1
smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 0.0204448, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 1, Final residual = 0.0369173, No Iterations 1
Pressure gradient source: uncorrected Ubar = 0.999835, pressure gradient = 0.0149157
GAMG: Solving for p, Initial residual = 1, Final residual = 0.0866042, No Iterations 3
time step continuity errors : sum local = 0.000351997, global = -5.26749e-18, cumulative = -5.26749e-18
Pressure gradient source: uncorrected Ubar = 0.999881, pressure gradient = 0.0107589
ExecutionTime = 0.4 s ClockTime = 1 s // <- this!turbulence->correct() ではユーザーが設定した乱流モデルが呼び出され，渦粘性係数\(\nu_{t}\)や乱流応力の非線形項\(N_{ij}\)が計算される
↓OpenFOAMで実装されている乱流モデルについてはこちら
UEqn.H
次のステップが記述されている
3. 運動量保存則（Navier-Stokes方程式）を立てる
4. 運動量保存則の不足緩和を行う
5. 運動量保存則を解いて速度の予測値\(\mathbf{u}^{*}\)を求める
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
fvVectorMatrix UEqn
(
fvm::ddt(U) + fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
fvOptions.correct(U);
}3. 運動量保存則（Navier-Stokes方程式）を立てる
Ueqn は次のように定義され，離散化されている
\begin{align}
&\nabla\cdot\left(\mathbf{uu}\right)
+\frac{\partial \mathbf{u}}{\partial t}
-\nu_{e}\nabla\cdot\left\{\nabla\mathbf{u}+\left(\nabla\mathbf{u}\right)^{T}\right\}
= 0 \\
\\ \rightarrow \quad
&\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{*}
+\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}
=\mathbf{s}_{_{P}}
\end{align}
ここではまだ圧力勾配 \(\nabla p\)が入っていない
fvVectorMatrix UEqn
(
fvm::ddt(U) + fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
);turbulence->divDevReff(U) は IncompressibleTurbulenceModel.C で定義されている
template<class TransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::IncompressibleTurbulenceModel<TransportModel>::divDevReff
(
volVectorField& U
) const
{
return divDevRhoReff(U);
}src/TurbulenceModels/incompressible/IncompressibleTurbulenceModel/IncompressibleTurbulenceModel.C
divDevRhoReff(U) は linearViscousStress.C で定義されている
template<class BasicTurbulenceModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
volVectorField& U
) const
{
return
(
- fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))
- fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
);
}src/TurbulenceModels/turbulenceModels/linearViscousStress/linearViscousStress.C
コードを見ると，粘性項のうち\(-\nu_{e}\nabla\cdot\left(\nabla\mathbf{u}\right)^{T}\) の項は陽的に処理され，\(-\nu_{e}\nabla\cdot\left(\nabla\mathbf{u}\right)\) の項は陰的に処理されていることがわかる
粘性係数 \(\nu_{e}=\nu+\nu_{t}\) は RASModel.H と LESModel.H でそれぞれ同じように定義されている
//- Return the effective viscosity
virtual tmp<volScalarField> nuEff() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject::groupName("nuEff", this->alphaRhoPhi_.group()),
this->nut() + this->nu()
)
);
}乱流応力の非線形項\(N_{ij}\)は nonlinearEddyViscosity.C で divDevReff に付け加えられる形で陽的に定義されている
template<class BasicTurbulenceModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::nonlinearEddyViscosity<BasicTurbulenceModel>::divDevRhoReff
(
volVectorField& U
) const
{
return
(
fvc::div(this->rho_*nonlinearStress_)
+ eddyViscosity<BasicTurbulenceModel>::divDevRhoReff(U)
);
}src/TurbulenceModels/turbulenceModels/nonlinearEddyViscosity/nonlinearEddyViscosity.C
4. 運動量保存則の不足緩和を行う
UEqn.relax() で運動量保存則の式に不足緩和が施される
\begin{align}
&\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{*}
+\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}
=\mathbf{s}_{_{P}}\\
\\ \rightarrow \quad
&\frac{1}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{*}
+\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}
=\mathbf{s}_{_{P}}
+\frac{1-\alpha_{u}}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{(n-1)}
\end{align}
UEqn.relax();5. 運動量保存則を解いて速度の予測値\(\mathbf{u}^{*}\)を求める
UEqn に圧力項が加えられ，離散化方程式を解くことで\(\mathbf{u}^{*}\)が求められる
\begin{align}
&\frac{1}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{*}
+\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}
=\mathbf{s}_{_{P}}
+\frac{1-\alpha_{u}}{\alpha_{u}}\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{(n-1)}
-\nabla p
\end{align}
if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
fvOptions.correct(U);
}pEqn.H
次のステップが記述されている
6. ポアソン方程式を立てる
7. ポアソン方程式を解いて\(p^{n+1}\)を求める
8. 質量流束\(\phi_{n+1}\)を修正する
9. 速度\(\mathbf{n}^{(n+1)}\)の修正する
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi))
);
MRF.makeRelative(phiHbyA);
adjustPhi(phiHbyA, U, p);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU, MRF);
// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}
#include "continuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);変数の宣言
UEqn.H で宣言した UEqn を次のように書き換える（\(\mathbf{A}_{_{P}}\)と\(s_{_{P}}\)に不足緩和の影響を含めている）
\begin{align}
&\mathbf{A}\mathbf{u}
= \mathbf{H} \\
\\ \rightarrow \quad
&\mathbf{A}=\mathbf{A}_{_{P}}, \quad
\mathbf{u}=\mathbf{u}_{_{P}}^{*}, \quad
\mathbf{H}=\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}
\end{align}
\(\mathbf{A}\)は UEqn.A() で取得でき，\(\mathbf{H}\) は UEqn.H() で取得できる
上記の式を用いて，いくつかの変数が定義される
\begin{align}
\mathrm{rAU}&=\frac{1}{\mathbf{A}}=\frac{1}{\mathbf{A}_{_{P}}} & \\
\mathrm{HbyA}&=\frac{\mathbf{H}}{\mathbf{A}}=\frac{\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}}{\mathbf{A}_{_{P}}}
=\mathbf{\tilde{u}}^{*}_{_{P}} & \\
\mathrm{phiHbyA}&=\left(\frac{\mathbf{H}}{\mathbf{A}}\right)_{f}
=\left(\frac{\mathbf{s}_{_{P}}-\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}}{\mathbf{A}_{_{P}}}\right)_{f}
=\left(\mathbf{\tilde{u}}^{*}_{_{P}}\right)_{f} &
\end{align}
ここで，添え字\(_{f}\)はセルの境界面上を意味し， surfaceScalarField として定義される
\(\mathrm{HbyA}\)は\(\mathbf{\tilde{u}}^{*}_{_{P}}\)であり，\(\mathrm{phiHbyA}\) は \(\mathbf{\tilde{u}}^{*}_{_{P}}\) による質量流束である
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi))
);
MRF.makeRelative(phiHbyA);
adjustPhi(phiHbyA, U, p);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU, MRF);いちおう A() と H() は fvMatrix.C で定義されているが，何が書いてあるかはさっぱりわからない
template<class Type>
Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
{
tmp<volScalarField> tAphi
(
new volScalarField
(
IOobject
(
"A("+psi_.name()+')',
psi_.instance(),
psi_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
psi_.mesh(),
dimensions_/psi_.dimensions()/dimVol,
extrapolatedCalculatedFvPatchScalarField::typeName
)
);
tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
tAphi.ref().correctBoundaryConditions();
return tAphi;
}
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
Foam::fvMatrix<Type>::H() const
{
tmp<GeometricField<Type, fvPatchField, volMesh>> tHphi
(
new GeometricField<Type, fvPatchField, volMesh>
(
IOobject
(
"H("+psi_.name()+')',
psi_.instance(),
psi_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
psi_.mesh(),
dimensions_/dimVol,
extrapolatedCalculatedFvPatchScalarField::typeName
)
);
GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi.ref();
// Loop over field components
for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
{
scalarField psiCmpt(psi_.primitiveField().component(cmpt));
scalarField boundaryDiagCmpt(psi_.size(), Zero);
addBoundaryDiag(boundaryDiagCmpt, cmpt);
boundaryDiagCmpt.negate();
addCmptAvBoundaryDiag(boundaryDiagCmpt);
Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
}
Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
addBoundarySource(Hphi.primitiveFieldRef());
Hphi.primitiveFieldRef() /= psi_.mesh().V();
Hphi.correctBoundaryConditions();
typename Type::labelType validComponents
(
psi_.mesh().template validComponents<Type>()
);
for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
{
if (validComponents[cmpt] == -1)
{
Hphi.replace
(
cmpt,
dimensionedScalar(Hphi.dimensions(), Zero)
);
}
}
return tHphi;
}constrainHbyA() は constrainHbyA.C で定義されており，境界面上での HbyA の値に何かしらの制限を加えていると考えられる
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
\*---------------------------------------------------------------------------*/
#include "constrainHbyA.H"
#include "volFields.H"
#include "fixedFluxExtrapolatedPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::tmp<Foam::volVectorField> Foam::constrainHbyA
(
const tmp<volVectorField>& tHbyA,
const volVectorField& U,
const volScalarField& p
)
{
tmp<volVectorField> tHbyANew;
if (tHbyA.isTmp())
{
tHbyANew = tHbyA;
tHbyANew.ref().rename("HbyA");
}
else
{
tHbyANew = new volVectorField("HbyA", tHbyA);
}
volVectorField& HbyA = tHbyANew.ref();
volVectorField::Boundary& HbyAbf = HbyA.boundaryFieldRef();
forAll(U.boundaryField(), patchi)
{
if
(
!U.boundaryField()[patchi].assignable()
&& !isA<fixedFluxExtrapolatedPressureFvPatchScalarField>
(
p.boundaryField()[patchi]
)
)
{
HbyAbf[patchi] = U.boundaryField()[patchi];
}
}
return tHbyANew;
}
// ************************************************************************* //src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C
fvc::flux() は fvcFlux.C で定義されている．名前の通り流束を計算していると思われる
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
\*---------------------------------------------------------------------------*/
#include "fvcFlux.H"
#include "surfaceInterpolate.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::tmp<Foam::surfaceScalarField> Foam::fvc::flux
(
const volVectorField& vvf
)
{
return scheme<vector>
(
vvf.mesh(),
"flux(" + vvf.name() + ')'
)().dotInterpolate(vvf.mesh().Sf(), vvf);
}
Foam::tmp<Foam::surfaceScalarField> Foam::fvc::flux
(
const tmp<volVectorField>& tvvf
)
{
tmp<surfaceScalarField> Flux(fvc::flux(tvvf()));
tvvf.clear();
return Flux;
}
// ************************************************************************* //phiHbyA の中にsimpleFoamにはなかった MRF.zeroFilter(...) という項が入っているが，何をしているのかはわからない
constrainPressure() は constrainPressure.C で定義されており，境界面上での p の値に何かしらの制限を加えていると考えられる
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
\*---------------------------------------------------------------------------*/
#include "constrainPressure.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "geometricOneField.H"
#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class RhoType, class RAUType, class MRFType>
void Foam::constrainPressure
(
volScalarField& p,
const RhoType& rho,
const volVectorField& U,
const surfaceScalarField& phiHbyA,
const RAUType& rhorAU,
const MRFType& MRF
)
{
const fvMesh& mesh = p.mesh();
volScalarField::Boundary& pBf = p.boundaryFieldRef();
const volVectorField::Boundary& UBf = U.boundaryField();
const surfaceScalarField::Boundary& phiHbyABf =
phiHbyA.boundaryField();
const typename RAUType::Boundary& rhorAUBf =
rhorAU.boundaryField();
const surfaceVectorField::Boundary& SfBf =
mesh.Sf().boundaryField();
const surfaceScalarField::Boundary& magSfBf =
mesh.magSf().boundaryField();
forAll(pBf, patchi)
{
if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))
{
refCast<fixedFluxPressureFvPatchScalarField>
(
pBf[patchi]
).updateSnGrad
(
(
phiHbyABf[patchi]
- rho.boundaryField()[patchi]
*MRF.relative(SfBf[patchi] & UBf[patchi], patchi)
)
/(magSfBf[patchi]*rhorAUBf[patchi])
);
}
}
}
template<class RAUType>
void Foam::constrainPressure
(
volScalarField& p,
const volScalarField& rho,
const volVectorField& U,
const surfaceScalarField& phiHbyA,
const RAUType& rAU
)
{
constrainPressure(p, rho, U, phiHbyA, rAU, NullMRF());
}
template<class RAUType, class MRFType>
void Foam::constrainPressure
(
volScalarField& p,
const volVectorField& U,
const surfaceScalarField& phiHbyA,
const RAUType& rAU,
const MRFType& MRF
)
{
constrainPressure(p, geometricOneField(), U, phiHbyA, rAU, MRF);
}
template<class RAUType>
void Foam::constrainPressure
(
volScalarField& p,
const volVectorField& U,
const surfaceScalarField& phiHbyA,
const RAUType& rAU
)
{
constrainPressure(p, U, phiHbyA, rAU, NullMRF());
}
// ************************************************************************* //src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C
6. ポアソン方程式を立てる
宣言した変数を使ってポアソン方程式 pEqn を立てる
\begin{align}
&\nabla\cdot\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n)}\right)
=\nabla\cdot\mathbf{\tilde{u}}^{*}_{_{P}} \\
\\ \rightarrow \quad
&\nabla\cdot\left(\frac{1}{\mathbf{A}}\nabla p^{(n)}\right)
=\nabla\cdot\left(\frac{\mathbf{H}}{\mathbf{A}}\right)_{f} \\
\end{align}
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);OpenFOAMのProgrammer’s Guideより，laplacian() と div() は次のようにセル境界面上での積分として定義されるので，pEqn の両辺はともに surfaceScalarField である
\begin{align}
&\int_{V}{\nabla\cdot\left(\Gamma\nabla\phi\right)}dV
=\int_{S}{d\mathbf{S}\cdot\left(\Gamma\nabla\phi\right)}
=\sum_{f}{\Gamma_{f}\mathbf{S}_{f}\cdot\left(\nabla\phi\right)_{f}} \\ \\
&\int_{V}{\nabla\cdot\phi}dV
=\int_{S}{d\mathbf{S}\cdot\phi}
=\sum_{f}{\mathbf{S}_{f}\cdot\phi_{f}}
\end{align}
7. ポアソン方程式を解いて\(p^{n+1}\)を求める
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));8. 質量流束\(\phi_{n+1}\)を修正する
有限体積法において，連続の式はセルの境界面上で定義される
(A)より，\(\left(\mathbf{\tilde{u}}^{*}\right)_{f}\)だけでは境界面上で連続の式を満たさないので，境界上で連続の式を満たすために質量流束\(\phi\)を次のように修正する
\begin{align}
\phi = \left(\mathbf{\tilde{u}}^{*}_{_{P}}\right)_{f}
-\left(\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n)}\right)_{f}
\end{align}
if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}OpenFOAMのProgrammer’s Guideより，fvm::laplacian(rAU(),p) は次のように計算される
\begin{align}
&\int_{V}{\nabla\cdot\left(\frac{1}{\mathbf{A}}\nabla p\right)}dV
=\int_{S}{d\mathbf{S}\cdot\left(\frac{1}{\mathbf{A}}\nabla p\right)}
\end{align}
上式より，\(\frac{1}{\mathbf{A}}\nabla p\) は pEqn において流束として取り扱われていることがわかる
flux() は fvMatrix.C で定義されているが，何が書いてあるかはさっぱりわからない
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
Foam::fvMatrix<Type>::
flux() const
{
if (!psi_.mesh().fluxRequired(psi_.name()))
{
FatalErrorInFunction
<< "flux requested but " << psi_.name()
<< " not specified in the fluxRequired sub-dictionary"
" of fvSchemes."
<< abort(FatalError);
}
// construct GeometricField<Type, fvsPatchField, surfaceMesh>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tfieldFlux
(
new GeometricField<Type, fvsPatchField, surfaceMesh>
(
IOobject
(
"flux("+psi_.name()+')',
psi_.instance(),
psi_.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
psi_.mesh(),
dimensions()
)
);
GeometricField<Type, fvsPatchField, surfaceMesh>& fieldFlux =
tfieldFlux.ref();
fieldFlux.setOriented();
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
{
fieldFlux.primitiveFieldRef().replace
(
cmpt,
lduMatrix::faceH(psi_.primitiveField().component(cmpt))
);
}
FieldField<Field, Type> InternalContrib = internalCoeffs_;
forAll(InternalContrib, patchi)
{
InternalContrib[patchi] =
cmptMultiply
(
InternalContrib[patchi],
psi_.boundaryField()[patchi].patchInternalField()
);
}
FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_;
forAll(NeighbourContrib, patchi)
{
if (psi_.boundaryField()[patchi].coupled())
{
NeighbourContrib[patchi] =
cmptMultiply
(
NeighbourContrib[patchi],
psi_.boundaryField()[patchi].patchNeighbourField()
);
}
}
typename GeometricField<Type, fvsPatchField, surfaceMesh>::
Boundary& ffbf = fieldFlux.boundaryFieldRef();
forAll(ffbf, patchi)
{
ffbf[patchi] = InternalContrib[patchi] - NeighbourContrib[patchi];
}
if (faceFluxCorrectionPtr_)
{
fieldFlux += *faceFluxCorrectionPtr_;
}
return tfieldFlux;
}9. 速度\(\mathbf{n}^{(n+1)}\)の修正する
時刻\(t_{n}\)における速度\(\mathbf{u}^{(n)}\)は次のように求められる
\begin{align}
\mathbf{u}^{(n)}_{_{P}}=\mathbf{\tilde{u}}^{*}_{_{P}}-\frac{1}{\mathbf{A}_{_{P}}}\nabla p^{(n)}
\end{align}
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);まとめ
OpenFOAMで実装されている非圧縮粘性流れの定常ソルバーであるpisoFoamについて解説した
OpenFOAMでは，次の1~11のステップによってPISO法を実装している
- 次の時間ステップに進む \(t=t_{n+1}\)
- 既知の\(\mathbf{u}^{(n)}\)と\(p^{(n)}\)を用いて\(\mathbf{u}^{(n+1)}\)と\(p^{(n+1)}\)を初期化する
- 運動量保存則（Navier-Stokes方程式）を立てる
- 運動量保存則の不足緩和を行う
- 運動量保存則を解いて速度の予測値\(\mathbf{u}^{*}\)を求める
- ポアソン方程式を立てる
- ポアソン方程式を解いて\(p^{n+1}\)を求める
- 質量流束\(\phi_{n+1}\)を修正する
- 速度\(\mathbf{n}^{(n+1)}\)の修正する
- 指定した回数だけステップ6に戻る
- 設定した計算終了時間に達していなければステップ2に戻る
↓OpenFOAMにおけるRhie-Chow補間について
（準備中）
↓その他のソルバーや離散化スキームについて
参考
≫非圧縮性流体計算の圧力-速度連成手法 (PDF)
コメント