【OpenFOAM】任意の断面での流量計算(cellZone)

こんにちは（@t_kun_kamakiri）

今回はOpenFOAMで以下のような配管内の水の流れを解析し、境界と任意の断面での流量を計算します。

境界の「inlet」「outlet1」「outlet2」での流量はOpenFOAMのfunction objectの機能を使うことで簡単に出力することができます。

前回の記事

【OpenFOAM】境界面での流量計算(volFieldValue)
※本記事のモデル作成方法はこちらをご参考ください。

しかし、c1～c5の断面での流量を出力するにはどうしたらいいのかなと思って自分なりに試してみた内容を書いています。

本記事の内容

OpenFOAMで任意の断面での流量を計算する。

【使用環境】
OpenFOAM-v2012(WSL2)

カマキリ

もっといい方法があったらコメントで教えてください_(._.)_

流量を出力する方針
②cellZoneSetを使ってfunction objectを使う
まとめ
参考書

流量を出力する方針

流量を計算する方法はパッと考えると以下の方法があるでしょう。

流量計測したい断面にあらかじめ面かセル領域を用意しておく
cellZoneSetを使ってfunction objectを使う
codedFunctionObjectという機能で直接コードを埋め込む

①の方法はあらかじめ計算する前に用意する必要があると考えて今回は試しませんでした。
今回は解析計算が終わった後でも出力できる方法として②と③を試したのでご紹介します。

③はこちら

【OpenFOAM】任意の断面での流量計算(cellZone)

②が本記事の内容

とはいえ、任意の断面での流量を直接出力するのはわからなかったので任意の領域を通過する平均流速を求めて面積を掛けて流量$Q=vA$とする方針にしました。
※任意の面を通過する流量とすると、入ってきて出ていく量の合計を計算するのことになって結局0になってしまいました。なので、いったん平均流速を求めて面積を掛けて流量を求めればいいやと思った次第です。

②cellZoneSetを使ってfunction objectを使う

まずは指定した領域を作るためにcellZoneをtopoSetで作ります。
以下のピンクの領域ですね。

topoSetの設定ファイルはデフォルトでは無いため適当なチュートリアルからコピーしてきます。

例えば以下のようなコマンドでtopoSetが使われているチュートリアルを探して適当にコピーします。

find $FOAM_TUTORIALS -name "topoSetDict"

1	find $FOAM_TUTORIALS -name "topoSetDict"

以下のように候補が出てきます。

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/system/topoSetDict
/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/system/topoSetDict
/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/driftFluxFoam/RAS/mixerVessel2D/system/topoSetDict
/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/potentialFreeSurfaceFoam/oscillatingBox/system/topoSetDict
(以下省略)

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/potentialFreeSurfaceDyMFoam/oscillatingBox/system/topoSetDict

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/system/topoSetDict

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/driftFluxFoam/RAS/mixerVessel2D/system/topoSetDict

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/potentialFreeSurfaceFoam/oscillatingBox/system/topoSetDict

(以下省略)

適当にコピーします。

cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/pimpleFoam/RAS/TJunctionFan/system/topoSetDict system/

1	cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/pimpleFoam/RAS/TJunctionFan/system/topoSetDict system/

topoSetの作り方はこちらの公式サイトを参考にします。

system/topoSetDict

actions
(
    {
        name    c1;
        type    cellZoneSet;
        action  new;  //新規追加
        source  boxToCell; //範囲の指定方法 boxToCellは六面体により指定
        box     (-0.21 -0.21 0.05) (0.21 0.21 0.06); //Min座標 Max座標
    }
    {
        name    c2;
        type    cellZoneSet;
        action  new;  //新規追加
        source  boxToCell; //範囲の指定方法 boxToCellは六面体により指定
        box     (-0.21 -0.21 0.12) (0.21 0.21 0.13); //Min座標 Max座標
    }
    {
        name    c3;
        type    cellZoneSet;
        action  new;  //新規追加
        source  boxToCell; //範囲の指定方法 boxToCellは六面体により指定
        box     (-0.21 -0.21 0.17) (0.21 0.21 0.18); //Min座標 Max座標
    }
    {
        name    c4;
        type    cellZoneSet;
        action  new;  //新規追加
        source  boxToCell; //範囲の指定方法 boxToCellは六面体により指定
        box     (-0.21 -0.21 0.25) (0.21 0.21 0.26); //Min座標 Max座標
    }
    {
        // Mandatory (inherited) entries
        name        c5;
        type        cellZoneSet;
        action      new;

        // Mandatory entries
        source      cylinderToCell; //始点終点と半径から円柱領域を指定
        p1          (0 -0.15 0.15);
        p2          (0 -0.16 0.15);
        radius      0.01;

        // Optional entries
        //innerRadius &lt;radius2&gt;;
    }
);

actions

(

{

name c1;

type cellZoneSet;

action new; //新規追加

source boxToCell; //範囲の指定方法 boxToCellは六面体により指定

box (-0.21 -0.21 0.05) (0.21 0.21 0.06); //Min座標 Max座標

}

{

name c2;

type cellZoneSet;

action new; //新規追加

source boxToCell; //範囲の指定方法 boxToCellは六面体により指定

box (-0.21 -0.21 0.12) (0.21 0.21 0.13); //Min座標 Max座標

}

{

name c3;

type cellZoneSet;

action new; //新規追加

source boxToCell; //範囲の指定方法 boxToCellは六面体により指定

box (-0.21 -0.21 0.17) (0.21 0.21 0.18); //Min座標 Max座標

}

{

name c4;

type cellZoneSet;

action new; //新規追加

source boxToCell; //範囲の指定方法 boxToCellは六面体により指定

box (-0.21 -0.21 0.25) (0.21 0.21 0.26); //Min座標 Max座標

}

{

// Mandatory (inherited) entries

name c5;

type cellZoneSet;

action new;

// Mandatory entries

source cylinderToCell; //始点終点と半径から円柱領域を指定

p1 (0 -0.15 0.15);

p2 (0 -0.16 0.15);

radius 0.01;

// Optional entries

//innerRadius <radius2>;

}

);

領域の作り方は1つではありません。
今回は以下の2つを使った方法を示していますが、どちらを使っても良いです。

boxToCell; //範囲の指定方法 boxToCellは六面体により指定
cylinderToCell; //始点終点と半径から円柱領域を指定

では、topoSetを実行します。

topoSet

topoSet

constant/polyMesh/setに以下のようにファイルができます。

constant/polyMesh/sets/
├── c1
├── c2
├── c3
├── c4
└── c5

constant/polyMesh/sets/

├── c1

├── c2

├── c3

├── c4

└── c5

例えば「c1」を見ると、

1216
(
86065
86066
81972
(以下省略)

1216

(

86065

86066

81972

(以下省略)

1216個のセルがあるということでしょうか。
この数字は次のcodedFunctionObjectでも出てきますのでちょっと記憶の片隅に残しておきましょう。

topoSetで作ったcellZoneはParaViewで確認することができます。

次にcontrolDictのfunction objectに以下を追加します。

    //Add
    volFieldValue1
    {
        // Mandatory entries (unmodifiable)
        type            volFieldValue;
        libs            (fieldFunctionObjects);

        // Mandatory entries (runtime modifiable)
        fields          (U p);
        operation       volAverage;
        regionType      cellZone;
        name            c1;

        // Optional entries (runtime modifiable)
        postOperation   none;
        weightField     alpha1;

        // Optional (inherited) entries
        writeFields     false;
        scalingFactor   1.0;
        writePrecision  8;
        writeToFile     true;
        useUserTime     true;
        region          region0;
        enabled         true;
        log             true;
        timeStart       0;
        timeEnd         10000;
        executeControl  timeStep;
        executeInterval 1;
        writeControl    timeStep;
        writeInterval   100;
    }
    volFieldValue2
    {
        $volFieldValue1
        name            c2;
    }
    volFieldValue3
    {
        $volFieldValue1
        name            c3;
    }
    volFieldValue4
    {
        $volFieldValue1
        name            c4;
    }
    volFieldValue5
    {
        $volFieldValue1
        name            c5;
    }

//Add

volFieldValue1

{

// Mandatory entries (unmodifiable)

type volFieldValue;

libs (fieldFunctionObjects);

// Mandatory entries (runtime modifiable)

fields (U p);

operation volAverage;

regionType cellZone;

name c1;

// Optional entries (runtime modifiable)

postOperation none;

weightField alpha1;

// Optional (inherited) entries

writeFields false;

scalingFactor 1.0;

writePrecision 8;

writeToFile true;

useUserTime true;

region region0;

enabled true;

log true;

timeStart 0;

timeEnd 10000;

executeControl timeStep;

executeInterval 1;

writeControl timeStep;

writeInterval 100;

}

volFieldValue2

{

$volFieldValue1

name c2;

}

volFieldValue3

{

$volFieldValue1

name c3;

}

volFieldValue4

{

$volFieldValue1

name c4;

}

volFieldValue5

{

$volFieldValue1

name c5;

}

cellZoneによる領域の物理量は平均流速以外にも以下の量を出力することができます。

none : no operation
sum : summation
sumMag : sum of magnitude
average: average
weightedAverage: weighted average
volAverage: volume average
weightedVolAverage: weighted volume average
volIntegrate: volume integral
min: minimum
max: maximum
CoV: coefficient of variation

none : no operation

sum : summation

sumMag : sum of magnitude

average: average

weightedAverage: weighted average

volAverage: volume average

weightedVolAverage: weighted volume average

volIntegrate: volume integral

min: minimum

max: maximum

CoV: coefficient of variation

参考記事

実行すればcellZoneでの平均流速と平均圧力が出力されます。

simpleFoam > log.simpleFoam &

1	simpleFoam > log.simpleFoam &

計算が終わった後にcontrolDictのfunction objectだけを実行したい場合は、

simpleFoam -postProcess

1	simpleFoam -postProcess

とします。

例えばc1に対しては、

postProcessing/volFieldValue1/0/volFieldValue.dat

# Region        : cellZone c1
# Cells         : 1216
# Volume        : 1.25374567e-05
# Time          	volAverage(U)	volAverage(p)
100             	(-6.91239995e-06 -5.72589460e-06 1.00545029e-01)	1.80967825e+01
200             	(-6.90135638e-06 -5.60361338e-06 1.00541477e-01)	1.80968310e+01
（省略)
2000            	(-6.71016485e-06 -4.74256383e-06 1.00544414e-01)	1.81051521e+01
2100            	(-6.84017444e-06 -5.02236067e-06 1.00543984e-01)	1.81051721e+01

# Region : cellZone c1

# Cells : 1216

# Volume : 1.25374567e-05

# Time volAverage(U) volAverage(p)

100 (-6.91239995e-06 -5.72589460e-06 1.00545029e-01) 1.80967825e+01

200 (-6.90135638e-06 -5.60361338e-06 1.00541477e-01) 1.80968310e+01

（省略)

2000 (-6.71016485e-06 -4.74256383e-06 1.00544414e-01) 1.81051521e+01

2100 (-6.84017444e-06 -5.02236067e-06 1.00543984e-01) 1.81051721e+01

流速が「-6.84017444e-06 -5.02236067e-06 1.00543984e-01」なので断面積を掛けると流量$Q=0.000126284$になりますね。

理論値との誤差は0.54%です。

その他にもc2~c5も同様に平均流速がpostProcessingに出力されています。

まとめ

今回は流量を出力する方法をまとめました。

任意の断面の流量計算(topoSet)

lob.simpleFoamのログファイルに、以下同じc1領域を計算したものを出力しています。

topoSetDictで作った領域からvolFieldValueで計算したのが上
codedFunctionObjectで計算したのが下

volFieldValue volFieldValue1 write:
    volAverage(c1) of U = (-6.71016e-06 -4.74256e-06 0.100544)
    volAverage(c1) of p = 18.1052
(省略)
============= codevolFieldValue ===============
volFieldValue codevolFieldValue write:
volAverage(codevolFieldValue) of U =(-6.71016e-06 -4.74256e-06 0.100544)
volAverage(codevolFieldValue) of p =18.1052
volAverage(codevolFieldValue) of phi = U*A = 0.000126284
Vi 1216
End