April 29, 2016 | CFD WITH A MISSION

OpenFOAM のソースコードを眺めていると，fvc::reconstruct という演算子を目にします．例えば，buoyantSimpleFoam の UEqn.H ではこの演算子が下記のように使用されています．

    // Solve the Momentum equation

    MRF.correctBoundaryVelocity(U);

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

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (simple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                )*mesh.magSf()
            )
        );

        fvOptions.correct(U);
    }

if (simple.momentumPredictor())

{

solve

(

UEqn

fvc::reconstruct

(

- ghf*fvc::snGrad(rho)

- fvc::snGrad(p_rgh)

)*mesh.magSf()

)

);

fvOptions.correct(U);

}

これは，具体的にどのような演算を行っているのでしょうか？こんな時は，ソースコードを確認しましょう．ソースコードを参照して実装方法をすぐに確認できるのが，オープンソースの１つのメリットですね．この演算子のソースファイルは，fvcReconstruct.C です．このファイルの 76 行目に演算の定義が記述されています．

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-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 "fvcReconstruct.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "fvcSurfaceIntegrate.H"
#include "extrapolatedCalculatedFvPatchFields.H"

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

namespace Foam
{

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

namespace fvc
{

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

template<class Type>
tmp
<
    GeometricField
    <
        typename outerProduct<vector,Type>::type, fvPatchField, volMesh
    >
>
reconstruct
(
    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
    typedef typename outerProduct<vector, Type>::type GradType;

    const fvMesh& mesh = ssf.mesh();

    surfaceVectorField SfHat(mesh.Sf()/mesh.magSf());

    tmp<GeometricField<GradType, fvPatchField, volMesh>> treconField
    (
        new GeometricField<GradType, fvPatchField, volMesh>
        (
            IOobject
            (
                "volIntegrate("+ssf.name()+')',
                ssf.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            inv(surfaceSum(SfHat*mesh.Sf()))&surfaceSum(SfHat*ssf),
            extrapolatedCalculatedFvPatchField<GradType>::typeName
        )
    );

    treconField.ref().correctBoundaryConditions();

    return treconField;
}

template<class Type>

tmp

GeometricField

typename outerProduct<vector,Type>::type, fvPatchField, volMesh

reconstruct

(

const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf

)

{

typedef typename outerProduct<vector, Type>::type GradType;

const fvMesh& mesh = ssf.mesh();

surfaceVectorField SfHat(mesh.Sf()/mesh.magSf());

tmp<GeometricField<GradType, fvPatchField, volMesh>> treconField

(

new GeometricField<GradType, fvPatchField, volMesh>

(

IOobject

(

"volIntegrate("+ssf.name()+')',

ssf.instance(),

mesh,

IOobject::NO_READ,

IOobject::NO_WRITE

inv(surfaceSum(SfHat*mesh.Sf()))&surfaceSum(SfHat*ssf),

extrapolatedCalculatedFvPatchField<GradType>::typeName

)

);

treconField.ref().correctBoundaryConditions();

return treconField;

}

私は一見しただけでは何を行っているのか分かりませんでした．ヘッダーファイル fvcReconstruct.H も確認してみます．

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

InNamespace
    Foam::fvc

Description
    Reconstruct volField from a face flux field.

SourceFiles
    fvcReconstruct.C

InNamespace

Foam::fvc

Description

Reconstruct volField from a face flux field.

SourceFiles

fvcReconstruct.C

冒頭の Description の項目に，Reconstruct volField from a face flux field という記述があります．最初に挙げた buoyantSimpleFoam での使用例を確認すると，fvc::reconstruct(引数) の引数は，フェイス（セル界面）で値をもっています（face flux field）．具体的には，

snGrad(rho) は，フェイスの法線方向の密度の勾配を，
mesh.magSf() は，フェイスの面積を

それぞれ計算するので，両方とも値はフェイスで保持しています．この例では，OpenFOAM の用語を使用すると，surfaceScalarField を引数をとして受け取り，セル中心で値をもつ volVectorField を返しています．

次回予告

次回の投稿では，簡単な例でこの演算子を書き下して，数式を使用した説明にチャレンジしたいと思います．最終的には，重力のような体積力を考慮した計算において，圧力の変数として \(p\) ではなく \(p_{rgh}\) を使用することの数値計算上の利点に関する説明につなげられたらと考えています．では！

M	T	W	T	F	S	S
				1	2	3
4	5	6	7	8	9	10
11	12	13	14	15	16	17
18	19	20	21	22	23	24
25	26	27	28	29	30

Day: April 29, 2016

reconstruct 演算子について #1