PR

【OpenFOAMv2012】simpleFoamの解説

OpenFOAMで実装されている非圧縮粘性流れの定常ソルバーであるsimpleFoamについて解説する

スポンサーリンク

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}\)を計算する,というような手法を考える必要がある

SIMPLE法

SIMPLE(Semi-Implicit Method for Pressure-Liked Equation)法は時間発展に陰解法を用いており,反復計算を行いながら速度と圧力を少しずつ真の解に近づけていく手法である

時刻 \(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}

ただし,\(p^{(n)}\)は次の式で修正される

\begin{equation}
p^{(n)}=p^{(n-1)}+\alpha_{p}\left(p^{(n)}-p^{(n-1)}\right)
\end{equation}

\(\alpha_{u}\)と\(\alpha_{p}\)は不足緩和係数であり,計算を安定させて収束を早めるはたらきをしている

以上の式を用いて,SIMPLE法では次に示す①~⑥の手順を収束するまで繰り返すことによって真の速度場と圧力を計算する

① 時刻\(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)

⑥ 時刻を\(t_{n+1}\)に進めて,①に戻る

不足緩和係数は,通常 \(\alpha_{u}=0.3 \sim 0.8\),\(\alpha_{p}=0.3 \sim 0.8\) の値をとり,試行錯誤的に決定される

係数を小さくすると計算は安定するが収束が遅くなるので,まずは小さい値から初めて,計算が安定するぎりぎりまで大きくしていくとよい

定常流れでは時間刻みの大きさは特に考慮しなくていいので,便宜的に\(\Delta t = 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}}\) は連続の式を満たしている

↓SIMPLE法の詳しい導出はこちら

OpenFOAMにおける実装の概要

OpenFOAMでは,次の1~11のステップによってSIMPLE法を実装している

OpenFOAM - User Guide - SIMPLE algorithm

  1. 次の時間ステップ \(t=t_{n+1}\) に進む
  2. 既知の\(\mathbf{u}^{(n)}\)と\(p^{(n)}\)を用いて\(\mathbf{u}^{(n+1)}\)と\(p^{(n+1)}\)を初期化する
  3. 運動量保存則(Navier-Stokes方程式)を立てる
  4. 運動量保存則の不足緩和を行う
  5. 運動量保存則を解いて速度の予測値\(\mathbf{u}^{*}\)を求める
  6. ポアソン方程式を立てる
  7. ポアソン方程式を解いて\(p^{n+1}\)を求める
  8. 質量流束\(\phi_{n+1}\)を修正する
  9. 圧力\(p^{(n+1)}\)の不足緩和を行う
  10. 速度\(\mathbf{n}^{(n+1)}\)の修正する
  11. 収束していなければステップ2に戻る

SIMPLE法は主に次の4つのファイルに分けて記述されている

上記の11のステップは,それぞれ次のファイルに記述されている

変数の定義 → createFields.H
1,2,11 → simpleFoam.C
3,4,5 → UEqn.C
6,7,8,9,10 → pEqn.H

それでは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, simple.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)
);


// ************************************************************************* //

$FOAM_SRC/finiteVolume/cfdTools/incompressible/createPhi.H

32~35行目では,fvsolution で定義した圧力の参照点や参照値が読み込まれている

SIMPLE
{
    residualControl
    {
        p               1e-4;
        U               1e-4;
        "(k|omega|epsilon)" 1e-4;
    }
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

$FOAM_TUTORIALS/incompressible/simpleFoam/windAroundBuildings/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/simpleFoam/windAroundBuildings/constant/turbulenceProperties

simpleFoam.C

mainにあたるプログラムで,while (simple.loop()){ } でSIMPLE法の反復計算が行われている

1回のループの中では UEqn.H → pEqn.H → turbulence->correct() が呼び出され,上述した手順で速度と圧力が求められている

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "simpleControl.H"
#include "fvOptions.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Steady-state solver for incompressible, turbulent flows."
    );

    #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 (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity SIMPLE corrector
        {
            #include "UEqn.H"
            #include "pEqn.H"
        }

        laminarTransport.correct();
        turbulence->correct();

        runTime.write();

        runTime.printExecutionTime(Info);
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

laminarTransport.correct(); は何をしているのかよくわからない

runTime.write(); でユーザーが IOobject::AUTO_WRITE とした変数が書き出される

runTime.printExecutionTime(Info); は log.simplefoam に計算時間を書き出している

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}^{*}\)を求める

    // Momentum predictor

    MRF.correctBoundaryVelocity(U);

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (simple.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\)が入っていない

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

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 (simple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
    }

pEqn.H

次のステップが記述されている

6. ポアソン方程式を立てる
7. ポアソン方程式を解いて\(p^{n+1}\)を求める
8. 質量流束\(\phi_{n+1}\)を修正する
9. 圧力\(p^{(n+1)}\)の不足緩和を行う
10. 速度\(\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.makeRelative(phiHbyA);
    adjustPhi(phiHbyA, U, p);

    tmp<volScalarField> rAtU(rAU);

    if (simple.consistent())
    {
        rAtU = 1.0/(1.0/rAU - UEqn.H1());
        phiHbyA +=
            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
        HbyA -= (rAU - rAtU())*fvc::grad(p);
    }

    tUEqn.clear();

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, U, phiHbyA, rAtU(), MRF);

    // Non-orthogonal pressure corrector loop
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - pEqn.flux();
        }
    }

    #include "continuityErrs.H"

    // Explicitly relax pressure for momentum corrector
    p.relax();

    // Momentum corrector
    U = HbyA - rAtU()*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.makeRelative(phiHbyA);
    adjustPhi(phiHbyA, U, p);

    tmp<volScalarField> rAtU(rAU);

    if (simple.consistent())
    {
        rAtU = 1.0/(1.0/rAU - UEqn.H1());
        phiHbyA +=
            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
        HbyA -= (rAU - rAtU())*fvc::grad(p);
    }

    tUEqn.clear();

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, U, phiHbyA, rAtU(), 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;
}

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

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;
}


// ************************************************************************* //

src/finiteVolume/finiteVolume/fvc/fvcFlux.C

if (simple.consistent()){ } は,SIMPLE法の改良版であるSIMPLEC法のためのコードなのでこの記事では解説しない(つーか これが限界)

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(rAtU(), 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}

Programmer’s Guide(PDFダウンロード)

7. ポアソン方程式を解いて\(p^{n+1}\)を求める
    pEqn.solve();
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 (simple.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;
}

≫src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

9. 圧力\(p^{(n+1)}\)の不足緩和を行う

圧力は次の式で示すように陽的に不足緩和が行われる

\begin{equation}
p^{(n)}=p^{(n-1)}+\alpha_{p}\left(p^{(n)}-p^{(n-1)}\right)
\end{equation}

    // Explicitly relax pressure for momentum corrector
    p.relax();
10. 速度\(\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}

    // Momentum corrector
    U = HbyA - rAtU()*fvc::grad(p);
    U.correctBoundaryConditions();
    fvOptions.correct(U);

まとめ

OpenFOAMで実装されている非圧縮粘性流れの定常ソルバーであるsimpleFoamについて解説した

OpenFOAMでは,次の1~11のステップによってSIMPLE法を実装している

OpenFOAM - User Guide - SIMPLE algorithm

  1. 次の時間ステップに進む \(t=t_{n+1}\)
  2. 既知の\(\mathbf{u}^{(n)}\)と\(p^{(n)}\)を用いて\(\mathbf{u}^{(n+1)}\)と\(p^{(n+1)}\)を初期化する
  3. 運動量保存則(Navier-Stokes方程式)を立てる
  4. 運動量保存則の不足緩和を行う
  5. 運動量保存則を解いて速度の予測値\(\mathbf{u}^{*}\)を求める
  6. ポアソン方程式を立てる
  7. ポアソン方程式を解いて\(p^{n+1}\)を求める
  8. 質量流束\(\phi_{n+1}\)を修正する
  9. 圧力\(p^{(n+1)}\)の不足緩和を行う
  10. 速度\(\mathbf{n}^{(n+1)}\)の修正する
  11. 収束していなければステップ2に戻る

↓OpenFOAMにおけるRhie-Chow補間について

(準備中)

↓その他のソルバーや離散化スキームについて

参考

非圧縮性流体計算の圧力-速度連成手法 (PDF)

コメント

  1. k より:

    3. 運動量保存則(Navier-Stokes方程式)を立てる

    UEqnの中に偏差テンソルの第2項が無いのはどうしてでしょうか?
    計算過程で打ち消されるのかな…?