fixedFluxPressure | CFD WITH A MISSION

OpenFOAM で自然対流を伴う熱流体計算や混相流 (Multiphase) 計算のように重力を考慮する計算を行う場合には，流れ場を表す変数として静圧 $p$ そのものではなく，次式\eqref{eq:prgh}で定義される変数 $p_{rgh}$ を使用します．

\begin{equation}
p_{rgh} = p\;-\;\rho{\bf g} \cdot {\bf x} \tag{1} \label{eq:prgh}
\end{equation}

ここで，$\rho$ は流体の密度を，$\boldsymbol{g}$ は重力加速度ベクトルを，$\boldsymbol{x}$ は位置ベクトルをそれぞれ表します．OpenFOAM で p_rgh と表されるこの変数の境界条件として，以前は，buoyantPressure 条件が使用されていましたが，v2.3 以降では，fixedFluxPressure 条件が使用されています．今回の投稿では，この２つの境界条件について取り上げます．

fixedFluxPressure 境界条件

セル界面の流速 $\boldsymbol{U}_f$ については，次式\eqref{eq:2}が成り立ちます．

\begin{equation}
{\bf U}_f \cdot {\bf S}_f = \left(\frac{{\bf H}}{A_p}\right)_f \cdot {\bf S}_f\;-\;\left(\frac{1}{A_p}\right)_f (\nabla p_{rgh})_f \cdot {\bf S}_f \tag{2} \label{eq:2}
\end{equation}

この式\eqref{eq:2}を変形すると，次式\eqref{eq:3}が得られます．

\begin{equation}
(\nabla p_{rgh})_f \cdot {\bf n}_f = \frac{1}{|{\bf S}_f|\left(\frac{1}{A_p}\right)_f} \left(\left(\frac{{\bf H}}{A_p}\right)_f \cdot {\bf S}_f\;-\;{\bf U}_f \cdot {\bf S}_f \right) \tag{3} \label{eq:3}
\end{equation}

fixedFluxPressure 境界条件は，その境界で課されている流速 $\boldsymbol{U}$ に対する境界条件から決まる流束を満たすように，境界の法線方向の圧力勾配を式\eqref{eq:3}の右辺により計算します．

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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 <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "fixedFluxPressureFvPatchScalarField.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedGradientFvPatchScalarField(p, iF),
    curTimeIndex_(-1)
{}


Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    fixedGradientFvPatchScalarField(p, iF),
    curTimeIndex_(-1)
{
    if (dict.found("value") && dict.found("gradient"))
    {
        fvPatchField<scalar>::operator=
        (
            scalarField("value", dict, p.size())
        );
        gradient() = scalarField("gradient", dict, p.size());
    }
    else
    {
        fvPatchField<scalar>::operator=(patchInternalField());
        gradient() = 0.0;
    }
}


Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
(
    const fixedFluxPressureFvPatchScalarField& ptf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedGradientFvPatchScalarField(p, iF),
    curTimeIndex_(-1)
{
    patchType() = ptf.patchType();

    // Map gradient. Set unmapped values and overwrite with mapped ptf
    gradient() = 0.0;
    gradient().map(ptf.gradient(), mapper);

    // Evaluate the value field from the gradient if the internal field is valid
    if (notNull(iF) && iF.size())
    {
        scalarField::operator=
        (
            //patchInternalField() + gradient()/patch().deltaCoeffs()
            // ***HGW Hack to avoid the construction of mesh.deltaCoeffs
            // which fails for AMI patches for some mapping operations
            patchInternalField() + gradient()*(patch().nf() & patch().delta())
        );
    }
}


Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
(
    const fixedFluxPressureFvPatchScalarField& wbppsf
)
:
    fixedGradientFvPatchScalarField(wbppsf),
    curTimeIndex_(-1)
{}


Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
(
    const fixedFluxPressureFvPatchScalarField& wbppsf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedGradientFvPatchScalarField(wbppsf, iF),
    curTimeIndex_(-1)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs
(
    const scalarField& snGradp
)
{
    if (updated())
    {
        return;
    }

    curTimeIndex_ = this->db().time().timeIndex();

    gradient() = snGradp;
    fixedGradientFvPatchScalarField::updateCoeffs();
}

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs

(

const scalarField& snGradp

)

{

if (updated())

{

return;

}

curTimeIndex_ = this->db().time().timeIndex();

gradient() = snGradp;

fixedGradientFvPatchScalarField::updateCoeffs();

}

この境界条件の updateCoeffs メンバ関数の 137行目の勾配値 snGradp は，最新の開発版では，constrainPressure 関数 (constrainPressure.C) で計算されます．式\eqref{eq:3}の右辺の式が，以下のソースコードの 68-73 行目に対応しています．

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2016 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 <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#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::GeometricBoundaryField& pBf = p.boundaryField();

    const volVectorField::GeometricBoundaryField& UBf = U.boundaryField();
    const surfaceScalarField::GeometricBoundaryField& phiHbyABf =
        phiHbyA.boundaryField();
    const typename RAUType::GeometricBoundaryField& rhorAUBf =
        rhorAU.boundaryField();
    const surfaceVectorField::GeometricBoundaryField& SfBf =
        mesh.Sf().boundaryField();
    const surfaceScalarField::GeometricBoundaryField& magSfBf =
        mesh.magSf().boundaryField();

    forAll(pBf, patchi)
    {
        if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))
        {
            refCast<fixedFluxPressureFvPatchScalarField>
            (
                pBf[patchi]
            ).updateCoeffs
            (
                (
                    phiHbyABf[patchi]
                  - rho.boundaryField()[patchi]
                   *MRF.relative(SfBf[patchi] & UBf[patchi], patchi)
                )
               /(magSfBf[patchi]*rhorAUBf[patchi])
            );
        }
    }
}

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::GeometricBoundaryField& pBf = p.boundaryField();

const volVectorField::GeometricBoundaryField& UBf = U.boundaryField();

const surfaceScalarField::GeometricBoundaryField& phiHbyABf =

phiHbyA.boundaryField();

const typename RAUType::GeometricBoundaryField& rhorAUBf =

rhorAU.boundaryField();

const surfaceVectorField::GeometricBoundaryField& SfBf =

mesh.Sf().boundaryField();

const surfaceScalarField::GeometricBoundaryField& magSfBf =

mesh.magSf().boundaryField();

forAll(pBf, patchi)

{

if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))

{

refCast<fixedFluxPressureFvPatchScalarField>

(

pBf[patchi]

).updateCoeffs

(

phiHbyABf[patchi]

- rho.boundaryField()[patchi]

*MRF.relative(SfBf[patchi] & UBf[patchi], patchi)

)

/(magSfBf[patchi]*rhorAUBf[patchi])

);

}

それ以前のバージョンでは，各ソルバーの pEqn.H において，境界での勾配計算が行われています．例えば，v3.0.x の interFoam の場合には，以下の部分でこの計算を行っています．

{
    volScalarField rAU("rAU", 1.0/UEqn.A());
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    MRF.makeRelative(phiHbyA);
    adjustPhi(phiHbyA, U, p_rgh);

    surfaceScalarField phig
    (
        (
            mixture.surfaceTensionForce()
          - ghf*fvc::snGrad(rho)
        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

    // Update the fixedFluxPressure BCs to ensure flux consistency
    setSnGrad<fixedFluxPressureFvPatchScalarField>
    (
        p_rgh.boundaryField(),
        (
            phiHbyA.boundaryField()
          - MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())
        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
    );

// Update the fixedFluxPressure BCs to ensure flux consistency

setSnGrad<fixedFluxPressureFvPatchScalarField>

(

p_rgh.boundaryField(),

(

phiHbyA.boundaryField()

- MRF.relative(mesh.Sf().boundaryField() & U.boundaryField())

)/(mesh.magSf().boundaryField()*rAUf.boundaryField())

);

buoyantPressure 境界条件

この境界条件は，2.3.x 以降のバージョンには実装されていませんが，ここで簡単に触れたいと思います．この境界条件は，変数 p_rgh の境界法線方向 ${\bf n}$ の勾配を次式\eqref{eq:4}により計算します．

$$ \frac{\partial (p_{rgh})}{\partial {\bf n}} = -({\bf g} \cdot {\bf x}) \frac{\partial \rho}{\partial {\bf n}} \tag{4} \label{eq:4} $$

したがって，静圧 $p$ の境界法線方向の勾配は次式\eqref{eq:5}で与えられます．

$$ \frac{\partial p}{\partial {\bf n}} = \frac{\partial (p_{rgh})}{\partial {\bf n}}\;+\;\frac{\partial (rgh)}{\partial {\bf n}} = \rho ({\bf g} \cdot {\bf n}) \tag{5} \label{eq:5} $$

よって，法線方向 ${\bf n}$ が重力ベクトル ${\bf g}$ に直交するような境界においては，静圧 $p$ に関して勾配0条件を意味します．

$$ \frac{\partial p}{\partial {\bf n}} = 0 \tag{6} \label{eq:6} $$

対して法線方向が重力ベクトルに平行な境界においては，静圧の勾配を次式により計算します．

$$ \frac{\partial p}{\partial {\bf n}} = \pm \rho |{\bf g}| \tag{7} \label{eq:7} $$

式\eqref{eq:4}は，buoyantPressureFvPatchScalarField クラスのメンバ関数 updateCoeffs に実装されています．

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 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 <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "buoyantPressureFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "uniformDimensionedFields.H"


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::buoyantPressureFvPatchScalarField::
buoyantPressureFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedGradientFvPatchScalarField(p, iF),
    rhoName_("rho")
{}


Foam::buoyantPressureFvPatchScalarField::
buoyantPressureFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    fixedGradientFvPatchScalarField(p, iF),
    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
{
    if (dict.found("value") && dict.found("gradient"))
    {
        fvPatchField<scalar>::operator=
        (
            scalarField("value", dict, p.size())
        );
        gradient() = scalarField("gradient", dict, p.size());
    }
    else
    {
        fvPatchField<scalar>::operator=(patchInternalField());
        gradient() = 0.0;
    }
}


Foam::buoyantPressureFvPatchScalarField::
buoyantPressureFvPatchScalarField
(
    const buoyantPressureFvPatchScalarField& ptf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
    rhoName_(ptf.rhoName_)
{}


Foam::buoyantPressureFvPatchScalarField::
buoyantPressureFvPatchScalarField
(
    const buoyantPressureFvPatchScalarField& ptf
)
:
    fixedGradientFvPatchScalarField(ptf),
    rhoName_(ptf.rhoName_)
{}


Foam::buoyantPressureFvPatchScalarField::
buoyantPressureFvPatchScalarField
(
    const buoyantPressureFvPatchScalarField& ptf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedGradientFvPatchScalarField(ptf, iF),
    rhoName_(ptf.rhoName_)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::buoyantPressureFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const uniformDimensionedVectorField& g =
        db().lookupObject<uniformDimensionedVectorField>("g");

    const fvPatchField<scalar>& rho =
        patch().lookupPatchField<volScalarField, scalar>(rhoName_);

    // If the variable name is "p_rgh", "ph_rgh" or "pd"
    // assume it is p? - rho*g.h and set the gradient appropriately.
    // Otherwise assume the variable is the static pressure.
    if
    (
        dimensionedInternalField().name() == "p_rgh"
     || dimensionedInternalField().name() == "ph_rgh"
     || dimensionedInternalField().name() == "pd"
    )
    {
        gradient() = -rho.snGrad()*(g.value() & patch().Cf());
    }
    else
    {
        gradient() = rho*(g.value() & patch().nf());
    }

    fixedGradientFvPatchScalarField::updateCoeffs();
}

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

void Foam::buoyantPressureFvPatchScalarField::updateCoeffs()

{

if (updated())

{

return;

}

const uniformDimensionedVectorField& g =

db().lookupObject<uniformDimensionedVectorField>("g");

const fvPatchField<scalar>& rho =

patch().lookupPatchField<volScalarField, scalar>(rhoName_);

// If the variable name is "p_rgh", "ph_rgh" or "pd"

// assume it is p? - rho*g.h and set the gradient appropriately.

// Otherwise assume the variable is the static pressure.

(

dimensionedInternalField().name() == "p_rgh"

|| dimensionedInternalField().name() == "ph_rgh"

|| dimensionedInternalField().name() == "pd"

)

{

gradient() = -rho.snGrad()*(g.value() & patch().Cf());

}

else

{

gradient() = rho*(g.value() & patch().nf());

}

fixedGradientFvPatchScalarField::updateCoeffs();

}

終わりに

最新の開発版では fixedFluxExtrapolatedPressure という名前の境界条件が新たに実装されています．fixedFluxPressure 境界条件が，重力という体積力を伴う計算において，圧力勾配を適切に設定するための条件であるのに対し，この新しい境界条件は体積力の有無に関わらずより汎用的な圧力境界条件として設計されているようです．まだテストの段階であり，実用レベルに至るには更なる検証計算が必要とのことです．この境界条件についても調査して別の投稿で取り上げたいと思います．

Tag: fixedFluxPressure

fixedFluxPressure 境界条件