見出し画像

【OpenFOAM用のC++勉強(5)】境界面における速度場と圧力場の取得

物理・プログラミング・Python、Web制作・自己啓発

こんにちは。
OpenFOAM用のC++の勉強のメモを自分のために残しています。

OpenFOAM v2112

前回の復習（指定した座標での速度場の取得）

前回は指定した座標での速度場や圧力場の値を取得しました。
１．ベクトルとして座標を定義します。

vector U_pos(0.05, 0.05, 0);

２．座標からセルNo.を取得します。

abel cellUNo = mesh.findCell(U_pos);

３．セルNo.から値の取得

    label cellUNo = mesh.findCell(U_pos);
    Info << "U" << U_pos << " = " << U[cellUNo] << endl;
    Info << "U" << U_pos << " = " << U.ref()[cellUNo] << endl;

結果は以下のようになります。

U(0.05 0.05 0) = (-0.196068 0.0263898 0)
U(0.05 0.05 0) = (-0.196068 0.0263898 0)

境界面の情報を取得

全体のコードを示します。

#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    Info << "================== Create mesh ==================" << endl;
    argList args(argc, argv);

    #include "createTime.H"
    
    // mesh
    fvMesh mesh
    (
        IOobject
        (
            fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );
    Info << "runTime.timeName() = " << runTime.timeName() << endl;
    
    Info << "============= volScalarField p ===============" << endl;
    // volScalarField typedef=> GeometricField => DimensionedField => Field
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1DimensionedField.html

    word t_ = "0.5";
    
    volScalarField p
    (
        IOobject
        (
            "p",
            // runTime.timeName(), // word
            t_,
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info << t_ << endl;
    vector P_pos(0.05, 0.05, 0);
    Info << "get p value at P_pos = " << P_pos << endl;

    // mesh is fvMesh class.
    // fvMesh => polyMesh
    /*
    Foam::label Foam::polyMesh::findCell
    (
        const point& p,
        const cellDecomposition decompMode
    )
    */
    label cellPNo = mesh.findCell(P_pos);
    Info << "p" << P_pos << " = " << p[cellPNo] << endl;
    Info << "p" << P_pos << " = " << p.ref()[cellPNo] << endl; //https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html#a77a3ea1ce7e2adc04d27301292b095ae
 
    Info << "============= volVectorField U ===============" << endl;

    volVectorField U
    (
        IOobject
        (
            "U",
            // runTime.timeName(), // word
            t_,
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    vector U_pos(0.05, 0.05, 0);
    Info << "get U value at U_pos = " << U_pos << endl;

    label cellUNo = mesh.findCell(U_pos);
    Info << "U" << U_pos << " = " << U[cellUNo] << endl;
    Info << "U" << U_pos << " = " << U.ref()[cellUNo] << endl;

    Info << "============= Boundary p value ===============" << endl;
    Info << "mesh.boundaryMesh()" << endl;
    Info << mesh.boundaryMesh() << endl;
    
    Info << "============================" << endl;
     forAll(mesh.boundaryMesh(), patchI)
    {
        Info << "mesh.boundaryMesh()[" << patchI  << "] = \n" << 
        mesh.boundaryMesh()[patchI] << endl;
    }
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvBoundaryMesh.html#ac89fe0aab9fe12bd4f0d3807495fe72b
    label inletPressID = mesh.boundaryMesh().findPatchID("movingWall");
    Info << "inletvelID = " << inletPressID << endl;
    const fvPatchScalarField& PinletList = p.boundaryField()[inletPressID];

    forAll(PinletList, fi)
    {
        Info << "p[" << fi  << "] =" << PinletList[fi] << endl;
    }

    Info << "============= Boundary U value ===============" << endl;
    label inletvelID = mesh.boundaryMesh().findPatchID("movingWall");
    Info << "inletvelID = " << inletvelID << endl;
    const fvPatchVectorField& UinletList = U.boundaryField()[inletvelID];

    forAll(UinletList, fi)
    {
        Info << "U[" << fi  << "] =" << UinletList[fi] << endl;
    }

    Info << "============= p of cell value near boundary ===============" << endl;
    // fvPatchScalarField = fvPatchField<scalar>
    // fvpatch have path()
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvPatch.html#acb76222e9692fae6c93d7a49125229f7
    const labelList& PfaceCell =PinletList.patch().faceCells();

    forAll(PfaceCell, ci)
    {
        Info << "cellNo = " << PfaceCell[ci] << " : pFaceCell[" << ci  << "] =" << p[PfaceCell[ci]] << endl;
    }

    Info << "============= U of cell value near boundary ===============" << endl;
    const labelList& UfaceCell =UinletList.patch().faceCells();

    forAll(UfaceCell, ci)
    {
        Info << "cellNo = " << UfaceCell[ci] << " : UFaceCell[" << ci  << "] =" << U[UfaceCell[ci]] << endl;
    }

    return 0;
}

Make/files

meshOF.C

EXE = $(FOAM_USER_APPBIN)/mesh005

Make/options

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude

EXE_LIBS = \
    -lfiniteVolume

以下のコマンドでコンパイル。

 wmake 2>&1 | tee log

※チュートリアルは以下からコピーしてくる。

cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/icoFoam/cavity/cavity mycavity

以下のコマンドでmycavityのケースファイルに移動します。

cd mycavity

そしてblockMeshを実行しておきます。

blockMesh

これはメッシュ情報には以下のファイルが必要だからです。

  • constant/polyMesh/boundary

  • constant/polyMesh/faces

  • constant/polyMesh/neighbour

  • constant/polyMesh/owner

  • constant/polyMesh/points

そして流速Uや圧力pの情報を取得するには計算を実行しておく必要があります。

icoFoam

では、mycavityの中で以下のmesh005を実行してます。
実行結果は長いのでlog.mesh005に掃き出します。

mesh005 > log.meesh005

log.mesh007

================== Create mesh ==================
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _7bdb509494-20201222 OPENFOAM=2012
Arch   : "LSB;label=32;scalar=64"
Exec   : mesh007
Date   : Dec 11 2022
Time   : 21:47:31
Host   : DESKTOP-CT961AV
PID    : 4762
I/O    : uncollated
Case   : /mnt/d/openfoam/20221011_wolf_benchmark/OpenFOAM_Training/OF9/101programming/myWork/mesh007/mycavity
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

runTime.timeName() = 0
============= volScalarField p ===============
0.5
get p value at P_pos = (0.05 0.05 0)
p(0.05 0.05 0) = -0.00133453
p(0.05 0.05 0) = -0.00133453
============= volVectorField U ===============
get U value at U_pos = (0.05 0.05 0)
U(0.05 0.05 0) = (-0.196068 0.0263898 0)
U(0.05 0.05 0) = (-0.196068 0.0263898 0)
============= Boundary p value ===============
mesh.boundaryMesh()
3
(
    movingWall
    {
        type            wall;
        inGroups        1(wall);
        nFaces          20;
        startFace       760;
    }
    fixedWalls
    {
        type            wall;
        inGroups        1(wall);
        nFaces          60;
        startFace       780;
    }
    frontAndBack
    {
        type            empty;
        inGroups        1(empty);
        nFaces          800;
        startFace       840;
    }
)
============================
mesh.boundaryMesh()[0] = 
type            wall;
inGroups        1(wall);
nFaces          20;
startFace       760;

mesh.boundaryMesh()[1] = 
type            wall;
inGroups        1(wall);
nFaces          60;
startFace       780;

mesh.boundaryMesh()[2] = 
type            empty;
inGroups        1(empty);
nFaces          800;
startFace       840;

inletPressID = 0
p.boundaryField()[0][5] = -0.240143
p[0] = -4.36666
p[1] = -2.04463
p[2] = -0.973538
p[3] = -0.563177
p[4] = -0.357842
p[5] = -0.240143
p[6] = -0.170115
p[7] = -0.123779
p[8] = -0.087477
p[9] = -0.0535059
p[10] = -0.0166499
p[11] = 0.0276653
p[12] = 0.0846122
p[13] = 0.161734
p[14] = 0.272416
p[15] = 0.442103
p[16] = 0.718223
p[17] = 1.22571
p[18] = 2.41104
p[19] = 4.84854
============= Boundary U value ===============
inletvelID = 0
U[0] = (1 0 0)
U[1] = (1 0 0)
U[2] = (1 0 0)
U[3] = (1 0 0)
U[4] = (1 0 0)
U[5] = (1 0 0)
U[6] = (1 0 0)
U[7] = (1 0 0)
U[8] = (1 0 0)
U[9] = (1 0 0)
U[10] = (1 0 0)
U[11] = (1 0 0)
U[12] = (1 0 0)
U[13] = (1 0 0)
U[14] = (1 0 0)
U[15] = (1 0 0)
U[16] = (1 0 0)
U[17] = (1 0 0)
U[18] = (1 0 0)
U[19] = (1 0 0)
============= p of cell value near boundary ===============
PinletPatch = movingWall
cellNo = 380 : pFaceCell[0] = -4.36666
cellNo = 381 : pFaceCell[1] = -2.04463
cellNo = 382 : pFaceCell[2] = -0.973538
cellNo = 383 : pFaceCell[3] = -0.563177
cellNo = 384 : pFaceCell[4] = -0.357842
cellNo = 385 : pFaceCell[5] = -0.240143
cellNo = 386 : pFaceCell[6] = -0.170115
cellNo = 387 : pFaceCell[7] = -0.123779
cellNo = 388 : pFaceCell[8] = -0.087477
cellNo = 389 : pFaceCell[9] = -0.0535059
cellNo = 390 : pFaceCell[10] = -0.0166499
cellNo = 391 : pFaceCell[11] = 0.0276653
cellNo = 392 : pFaceCell[12] = 0.0846122
cellNo = 393 : pFaceCell[13] = 0.161734
cellNo = 394 : pFaceCell[14] = 0.272416
cellNo = 395 : pFaceCell[15] = 0.442103
cellNo = 396 : pFaceCell[16] = 0.718223
cellNo = 397 : pFaceCell[17] = 1.22571
cellNo = 398 : pFaceCell[18] = 2.41104
cellNo = 399 : pFaceCell[19] = 4.84854
============= U of cell value near boundary ===============
cellNo = 380 : UFaceCell[0] = (0.299604 0.146091 0)
cellNo = 381 : UFaceCell[1] = (0.396272 0.127503 0)
cellNo = 382 : UFaceCell[2] = (0.55742 0.0760895 0)
cellNo = 383 : UFaceCell[3] = (0.674487 0.0437233 0)
cellNo = 384 : UFaceCell[4] = (0.743889 0.0255721 0)
cellNo = 385 : UFaceCell[5] = (0.787154 0.0152422 0)
cellNo = 386 : UFaceCell[6] = (0.815486 0.00921237 0)
cellNo = 387 : UFaceCell[7] = (0.833956 0.00550944 0)
cellNo = 388 : UFaceCell[8] = (0.845386 0.00302719 0)
cellNo = 389 : UFaceCell[9] = (0.851372 0.00114935 0)
cellNo = 390 : UFaceCell[10] = (0.852667 -0.000499394 0)
cellNo = 391 : UFaceCell[11] = (0.849383 -0.00221509 0)
cellNo = 392 : UFaceCell[12] = (0.841 -0.00433859 0)
cellNo = 393 : UFaceCell[13] = (0.826141 -0.00742411 0)
cellNo = 394 : UFaceCell[14] = (0.802078 -0.0125323 0)
cellNo = 395 : UFaceCell[15] = (0.763621 -0.0217729 0)
cellNo = 396 : UFaceCell[16] = (0.699108 -0.0392829 0)
cellNo = 397 : UFaceCell[17] = (0.583952 -0.0726601 0)
cellNo = 398 : UFaceCell[18] = (0.415827 -0.127868 0)
cellNo = 399 : UFaceCell[19] = (0.308819 -0.149468 0)

以下コードの解説です。

    Info << "============= Boundary p value ===============" << endl;

からコードを追加しています。

mesh.boundaryMesh()はpolyMeshクラスの中のメンバ関数boundaryMeshで以下のように、polyBoundaryMeshを返します。
const polyBoundaryMesh & boundaryMesh()const

https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1polyBoundaryMesh.html
    Info << "============= Boundary p value ===============" << endl;
    Info << "mesh.boundaryMesh()" << endl;
    Info << mesh.boundaryMesh() << endl;

これにより結果は境界面の情報をリストとして取得することができます。

============= Boundary p value ===============
mesh.boundaryMesh()
3
(
    movingWall
    {
        type            wall;
        inGroups        1(wall);
        nFaces          20;
        startFace       760;
    }
    fixedWalls
    {
        type            wall;
        inGroups        1(wall);
        nFaces          60;
        startFace       780;
    }
    frontAndBack
    {
        type            empty;
        inGroups        1(empty);
        nFaces          800;
        startFace       840;
    }
)

for文でリストからひとつひとつ要素を取得することもできます。

    Info << "============= Boundary p value ===============" << endl;
    forAll(mesh.boundaryMesh(), patchI)
    {
        Info << "mesh.boundaryMesh()[" << patchI  << "] = \n" << 
        mesh.boundaryMesh()[patchI] << endl;
    }

境界面の名前から境界面の速度と圧力の値を取得

１．境界面の名前からパッチのIDを取得

fvMeshはpolyMeshを継承しており、polyMeshのメンバ関数であるboundaryMesh()はpolyBoundaryMeshを返しています。polyBoundaryMeshはfindPatchIDというメンバ関数があり、これが境界面パッチのIDを取得する関数です。

    label inletPressID = mesh.boundaryMesh().findPatchID("movingWall");
    Info << "inletPressID = " << inletPressID << endl;

inletPressID = 0とパッチ面のIDを取得できました。
２．pはvolScalarField型（GeometricFieldのスカラー版がvolScalarField型）なので、GeometricFieldのメンバ関数boundaryFIeldを使うことで速度や圧力の値を取得できます。

     const fvPatchScalarField& PinletList = p.boundaryField()[inletPressID];

    forAll(PinletList, fi)
    {
        Info << "p[" << fi  << "] =" << PinletList[fi] << endl;
    }

結果は以下のようになります。

inletPressID = 0
p[0] =-4.36666
p[1] =-2.04463
p[2] =-0.973538
p[3] =-0.563177
p[4] =-0.357842

ちなみにp.boundaryField()だけだと以下のように境界面の定義のリストが出力されます。

p.boundaryField()

3
(
    type            zeroGradient;

    type            zeroGradient;

    type            empty;

)

以下のようにp.boundaryField()[境界面ID]から続けて[5]のようにindexを入れると値の取得ができます。

    Info << "p.boundaryField()[0][5] = " << p.boundaryField()[0][5] << endl;

結果は以下のようになります。

p.boundaryField()[0][5] = -0.240143

同様に速度についても以下のように境界面での値を取得できます。

    Info << "============= Boundary U value ===============" << endl;
    label inletvelID = mesh.boundaryMesh().findPatchID("movingWall");
    Info << "inletvelID = " << inletvelID << endl;
    const fvPatchVectorField& UinletList = U.boundaryField()[inletvelID];

    forAll(UinletList, fi)
    {
        Info << "U[" << fi  << "] =" << UinletList[fi] << endl;
    }

境界面近傍のセルの値での速度と圧力の取得

fvPatchVectorFieldやfvPatchVectorFieldは以下のように、

  • fvPatchScalarField = fvPatchField<scalar>

  • fvPatchVectorField = fvPatchField<vector>

と定義されているので実態はfvPatch型です。

https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvPatchField.html
     Info << "============= p of cell value near boundary ===============" << endl;
    // fvPatchScalarField = fvPatchField<scalar>
    // fvpatch have path()
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvPatch.html#acb76222e9692fae6c93d7a49125229f7
    const fvPatch& PinletPatch = PinletList.patch();
    Info << "PinletPatch = " << PinletPatch.name() << endl;

    const labelList& PfaceCell =PinletList.patch().faceCells();

    forAll(PfaceCell, ci)
    {
        Info << "cellNo = " << PfaceCell[ci] << " : pFaceCell[" << ci  << "] = " << p[PfaceCell[ci]] << endl;
    }

    Info << "============= U of cell value near boundary ===============" << endl;
    const labelList& UfaceCell =UinletList.patch().faceCells();

    forAll(UfaceCell, ci)
    {
        Info << "cellNo = " << UfaceCell[ci] << " : UFaceCell[" << ci  << "] = " << U[UfaceCell[ci]] << endl;
    }
const fvPatch& PinletPatch = PinletList.patch()

とするとfvPatch型を返すのでメンバ関数findCellsにより値のリストを返します。

https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvPatch.html#a05ff7aca2b2d901acaec7624079f415e

返ってきたのがリストなのでfor文でひとつひとつ要素を取り出します。

    forAll(PfaceCell, ci)
    {
        Info << "cellNo = " << PfaceCell[ci] << " : pFaceCell[" << ci  << "] = " << p[PfaceCell[ci]] << endl;
    }

    Info << "============= U of cell value near boundary ===============" << endl;
    const labelList& UfaceCell =UinletList.patch().faceCells();

    forAll(UfaceCell, ci)
    {
        Info << "cellNo = " << UfaceCell[ci] << " : UFaceCell[" << ci  << "] = " << U[UfaceCell[ci]] << endl;
    }

結果は以下のようになります。

============= p of cell value near boundary ===============
PinletPatch = movingWall
cellNo = 380 : pFaceCell[0] = -4.36666
cellNo = 381 : pFaceCell[1] = -2.04463
cellNo = 382 : pFaceCell[2] = -0.973538
cellNo = 383 : pFaceCell[3] = -0.563177
cellNo = 384 : pFaceCell[4] = -0.357842
cellNo = 385 : pFaceCell[5] = -0.240143
cellNo = 386 : pFaceCell[6] = -0.170115
cellNo = 387 : pFaceCell[7] = -0.123779
cellNo = 388 : pFaceCell[8] = -0.087477
cellNo = 389 : pFaceCell[9] = -0.0535059
cellNo = 390 : pFaceCell[10] = -0.0166499
cellNo = 391 : pFaceCell[11] = 0.0276653
cellNo = 392 : pFaceCell[12] = 0.0846122
cellNo = 393 : pFaceCell[13] = 0.161734
cellNo = 394 : pFaceCell[14] = 0.272416
cellNo = 395 : pFaceCell[15] = 0.442103
cellNo = 396 : pFaceCell[16] = 0.718223
cellNo = 397 : pFaceCell[17] = 1.22571
cellNo = 398 : pFaceCell[18] = 2.41104
cellNo = 399 : pFaceCell[19] = 4.84854
============= U of cell value near boundary ===============
cellNo = 380 : UFaceCell[0] = (0.299604 0.146091 0)
cellNo = 381 : UFaceCell[1] = (0.396272 0.127503 0)
cellNo = 382 : UFaceCell[2] = (0.55742 0.0760895 0)
cellNo = 383 : UFaceCell[3] = (0.674487 0.0437233 0)
cellNo = 384 : UFaceCell[4] = (0.743889 0.0255721 0)
cellNo = 385 : UFaceCell[5] = (0.787154 0.0152422 0)
cellNo = 386 : UFaceCell[6] = (0.815486 0.00921237 0)
cellNo = 387 : UFaceCell[7] = (0.833956 0.00550944 0)
cellNo = 388 : UFaceCell[8] = (0.845386 0.00302719 0)
cellNo = 389 : UFaceCell[9] = (0.851372 0.00114935 0)
cellNo = 390 : UFaceCell[10] = (0.852667 -0.000499394 0)
cellNo = 391 : UFaceCell[11] = (0.849383 -0.00221509 0)
cellNo = 392 : UFaceCell[12] = (0.841 -0.00433859 0)
cellNo = 393 : UFaceCell[13] = (0.826141 -0.00742411 0)
cellNo = 394 : UFaceCell[14] = (0.802078 -0.0125323 0)
cellNo = 395 : UFaceCell[15] = (0.763621 -0.0217729 0)
cellNo = 396 : UFaceCell[16] = (0.699108 -0.0392829 0)
cellNo = 397 : UFaceCell[17] = (0.583952 -0.0726601 0)
cellNo = 398 : UFaceCell[18] = (0.415827 -0.127868 0)
cellNo = 399 : UFaceCell[19] = (0.308819 -0.149468 0)

ちなみにOpenFOAMのコードの中で出てくるforAllはC++の文法というわけではなく、OpenFOAMが独自で定義したforループです。

#define forAll(list, i) \
    for (Foam::label i=0; i<(list).size(); ++i)

と定義しています。

まだまだわかっていない部分が多いですが、引き続き勉強します。

参考にしているOpenFOAMの辞書です。
Foundation版をベースにしていますがわからないことはこちらで調べています。

引き続きC++の勉強を続けます。

Twitter➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ（物理ブログ）

いいなと思ったら応援しよう！