OpenFOAMで実装されている非圧縮層流流れの非定常ソルバーであるicoFoamについて解説する
CFDの基礎方程式
密度\(\rho\)，粘性係数\(\mu\)が一定な非圧縮性流体の基礎方程式は，連続の式（質量保存式）およびNavier-Stokes方程式（運動量保存式）である
層流流れ \(\nu_{t}=0\) であることを仮定すると
（連続の式）
\begin{equation}
\nabla\cdot\mathbf{u}=0
\end{equation}
（Navier-Stokes方程式）
\begin{equation}
\frac{\partial \mathbf{u}}{\partial t}+\nabla\cdot\left(\mathbf{uu}\right)
=-\frac{1}{\rho}\nabla P
+\nu\nabla\cdot\left(\nabla \mathbf{u}\right)
+\mathbf{f}
\end{equation}
\begin{equation}
{P}={p}+\frac{2}{3}\rho k
\end{equation}
密度\(\rho\)と動粘性係数\(\nu\)は水や空気などの物性値である
未知数は速度\(\mathbf{u}\)の3つの成分\(u_{i}\)，圧力\(p\)の4つであり，連続の式とNavier-Stokes方程式を解くことですべての物理量を求めることができる
また，これ以降は簡単のため \(P \rightarrow p\) と表記することとする
圧力-速度連成
連続の式とNavier-Stokes方程式には次のような問題がある
- 連続の式には時間微分項が含まれておらず，速度場に対して強い制約（発散なし：divergence free）を課しているだけである
- 変数である圧力\(p\)に対する時間微分項\(\frac{\partial p}{\partial t}\)が明示的に示されていない
流れ場の速度と圧力を計算するためには，基礎方程式で明示的に示されていない圧力\(p\)を計算しつつ，連続の式による強い制約を満たすような速度\(u_{i}\)を計算する，というような手法を考える必要がある
icoFoam
icoFoamは，PISO法を用いて非定常層流流れを計算するソルバーである
時刻 \(t_{n}\) における速度の推定値は次の式で求められる
\begin{equation}
\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{*}
+\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}
=-\nabla p^{(n-1)}+\mathbf{s}_{_{P}}
\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}
以上の式を用いて，層流流れを仮定したPISO法では次に示す①～⑦の手順を収束するまで繰り返すことによって真の速度場と圧力を計算する
① 時刻\(t_{n}\)において偏微分方程式を離散化して係数\(\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}\)に進めて，①に戻る
ちなみに，(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~10のステップによってicoFoamを実装している
- 次の時間ステップ \(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に戻る
icoFoamは主に次の2つのファイルに分けて記述されている
変数の定義は createFields.H で行い，上記の10のステップはすべて icoFoam.C に記述されている
createFields.H
圧力\(p\)と速度\(\mathbf{u}\)の宣言などが行われている
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
dimensionedScalar nu
(
"nu",
dimViscosity,
transportProperties
);
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, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());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)
);
// ************************************************************************* //55~58行目では，fvsolution で定義した圧力の参照点や参照値が読み込まれている
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}tutorials/incompressible/icoFoam/cavity/cavity/system/fvSolution
icoFoam.C
mainにあたるプログラムである
while ( runTime.loop() ){...} でPISO法の時間発展が計算され，1回の時間ステップの中では UEqn → PISO loop が呼び出され，while ( piso.correct() ){...} で，ユーザーが nCorrectors で指定した回数だけ速度と圧力の修正を行う（通常は2回）
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "pisoControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"Transient solver for incompressible, laminar flow"
" of Newtonian fluids."
);
#include "postProcess.H"
#include "addCheckCaseOptions.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
pisoControl piso(mesh);
#include "createFields.H"
#include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "CourantNo.H"
// Momentum predictor
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);
if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}
// --- PISO loop
while (piso.correct())
{
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*icoFoam.C(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);
// 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();
}
runTime.write();
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
ico loop の回数は fvSolution の中で nCorrectors として指定する
ico
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
}tutorials/incompressible/icoFoam/LES/pitzDaily/system/fvSolution
laminarTransport.correct() は何をしているのかよくわからない
runTime.write() でユーザーが IOobject::AUTO_WRITE とした変数が書き出される
runTime.printExecutionTime(Info) は log.icofoam に計算時間を書き出している
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!3. 運動量保存則（Navier-Stokes方程式）を立てる
Ueqn は次のように定義され，離散化されている
\begin{align}
&\frac{\partial \mathbf{u}}{\partial t}
+\nabla\cdot\left(\mathbf{uu}\right)
-\nu\nabla\cdot\left(\nabla\mathbf{u}\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)
- fvm::laplacian(nu, U)
);4. 運動量保存則を解いて速度の予測値\(\mathbf{u}^{*}\)を求める
UEqn に圧力項が加えられ，離散化方程式を解くことで\(\mathbf{u}^{*}\)が求められる
\begin{align}
\mathbf{A}_{_{P}}\mathbf{u}_{_{P}}^{*}
+\sum{\mathbf{A}_{nb}\mathbf{u}_{nb}^{*}}
=\mathbf{s}_{_{P}}
-\nabla p
\end{align}
if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}変数の宣言
ここからは while ( piso.correct() ){...} の中のプログラムである
103～108行目で宣言した UEqn を次のように書き換える
\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)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);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 の中にpisoFoamにはなかった fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) という項が入っているが，何をしているのかはわからない
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
5. ポアソン方程式を立てる
宣言した変数を使ってポアソン方程式 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}
6. ポアソン方程式を解いて\(p^{n+1}\)を求める
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));7. 質量流束\(\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;
}8. 速度\(\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();まとめ
OpenFOAMで実装されている非圧縮粘性流れの定常ソルバーであるicoFoamについて解説した
OpenFOAMでは，次の1~10のステップによって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)
コメント