LoginSignup
1
2

More than 1 year has passed since last update.

@openyanagi(Mitsuhiro YANAGI)in株式会社リーディング・エッジ社

ベンチュリ管解析

Posted at

はじめに

Notes on Computational Fluid Dynamics: General Principles(以下、テキストと略す)の第8章のサンプル問題に記載されているベンチュリ管(Venturi tube)解析を行った。
OpenFOAMは2022年7月時点で最新のFoundation版バージョン10を使用した。

参考にしたチュートリアルケース

blockMesh作成では、次のケースを参考にした。

  • $FOAM_TUTORIALS/multiphase/interFoam/laminar/sloshingCylinder

参考にしたサイト

ベンチュリ管の体積流量計算

テキストの(8.2)式は、サンプル問題の$\Delta p\ $から計算した$C\ $の値が合わない。
参考資料①の体積流量計算式を使用したところ、上述の問題が解決したため、その式を使用する。

Q = \frac{C}{\sqrt {\displaystyle \left(\frac {D}{D_{t}}\right )^{2}-1}} A \sqrt{2 \Delta p} \

変数の説明は下表の通り。

変数名 説明
$Q\ $ ベンチュリ管を流れる体積流量$(m^{3}/s)$
$D\ $ ベンチュリ管入口直径$(m)$
$D_{t}\ $ ベンチュリ管スロート直径$(m)$
$A\ $ ベンチュリ管入口断面積$(m^{2})$
$\Delta p\ $ ベンチュリ管入口とスロートとの間の差圧$({m^{2}}/{s^{2}})$
$C\ $ 流量係数$(-)$

ただし、$D_{t}={D}/{2}$、$A={\pi D^{2}}/{4}$であり、圧力は密度で割ったものである。

ハーゲン・ポアズイユ流れの流速

ベンチュリ管入口での流入流速の分布は一様とはせずに、定常流（ハーゲン・ポアズイユ流れ）の速度分布とする。速度分布は参考資料②より次式で与える。

\displaystyle \frac{u}{|\mathbf{u_{max}}|} = 1 - \left(\frac{r}{r_{0}} \right)^{2} \

変数の説明は下表の通り。

■変数名■ 説明
$u\ $ ハーゲン・ポアズイユ流れの流速$(m/s)$
$|\mathbf{u_{max}}|\ $ ベンチュリ管入口の最大流速$(m/s)$
$r\ $ ベンチュリ管半径方向長さ$(m)$
$r_{0}\ $ ベンチュリ管入口半径$(m)$

ただし、$r_{0}={D}/{2}$である。

幾何形状

ベンチュリ管の幾何形状を下図に示す。
fig1.png

3D解析モデル

ベンチュリ管の3D解析モデルを下図に示す。
fig2.png
3D解析モデルはblockMeshで作成した。
system/blockMeshを以下に示す。defaultPatch がv10の新機能である。

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

r0           0.05;
r0Neg       -0.05;
rt           0.025;
rtNeg       -0.025;

box0         0.02;
box0Neg     -0.02;
box1         0.02;
box1Neg     -0.02;
box2         0.01;
box2Neg     -0.01;
box3         0.01;
box3Neg     -0.01;
box4         0.02;
box4Neg     -0.02;
box5         0.02;
box5Neg     -0.02;

x0           0;
x1           0.1;   // D
x2           0.225; // x1 + 1.25*D
x3           0.275; // x2 + 0.5*D
x4           0.525; // x3 + 2.5*D
x5           0.625; // x4 + D

nR           12;
nX1          12; // x0 <= x <= x1
nX2          18; // x1 <= x <= x2
nX3           6; // x2 <= x <= x3
nX4          32; // x3 <= x <= x4
nX5          12; // x4 <= x <= x5

gd           2.7; // delta = 0.00145, box0 = 0.02, box2 = 0.01

verbose no;

geometry
{
    cylinder0
    {
        type      searchableCylinder;
        point1    (-1 0 0);
        point2    ( 1 0 0);
        radius    $r0;
    }
    cylinder1
    {
        type      searchableCylinder;
        point1    (-1 0 0);
        point2    ( 1 0 0);
        radius    $rt;
    }
}

scale 1;

vertices
(
    // x = x0
    // Inner
    ($x0 $box0Neg $box0Neg) // 0
    ($x0 $box0    $box0Neg) // 1
    ($x0 $box0Neg $box0)    // 2
    ($x0 $box0    $box0)    // 3
    // Outer block points
    project ($x0 $r0Neg $r0Neg) (cylinder0) // 4
    project ($x0 $r0    $r0Neg) (cylinder0) // 5
    project ($x0 $r0Neg $r0) (cylinder0)    // 6
    project ($x0 $r0    $r0) (cylinder0)    // 7
    // x = x1
    // Inner
    ($x1 $box1Neg $box1Neg) // 8
    ($x1 $box1    $box1Neg) // 9
    ($x1 $box1Neg $box1)    // 10
    ($x1 $box1    $box1)    // 11
    // Outer block points
    project ($x1 $r0Neg $r0Neg) (cylinder0) // 12
    project ($x1 $r0    $r0Neg) (cylinder0) // 13
    project ($x1 $r0Neg $r0) (cylinder0)    // 14
    project ($x1 $r0    $r0) (cylinder0)    // 15
    // x = x2
    // Inner
    ($x2 $box2Neg $box2Neg) // 16
    ($x2 $box2    $box2Neg) // 17
    ($x2 $box2Neg $box2)    // 18
    ($x2 $box2    $box2)    // 19
    // Outer block points
    project ($x2 $rtNeg $rtNeg) (cylinder1) // 20
    project ($x2 $rt    $rtNeg) (cylinder1) // 21
    project ($x2 $rtNeg $rt) (cylinder1)    // 22
    project ($x2 $rt    $rt) (cylinder1)    // 23
    // x = x3
    // Inner
    ($x3 $box3Neg $box3Neg) // 24
    ($x3 $box3    $box3Neg) // 25
    ($x3 $box3Neg $box3)    // 26
    ($x3 $box3    $box3)    // 27
    // Outer block points
    project ($x3 $rtNeg $rtNeg) (cylinder1) // 28
    project ($x3 $rt    $rtNeg) (cylinder1) // 29
    project ($x3 $rtNeg $rt) (cylinder1)    // 30
    project ($x3 $rt    $rt) (cylinder1)    // 31
    // x = x4
    // Inner
    ($x4 $box4Neg $box4Neg) // 32
    ($x4 $box4    $box4Neg) // 33
    ($x4 $box4Neg $box4)    // 34
    ($x4 $box4    $box4)    // 35
    // Outer block points
    project ($x4 $r0Neg $r0Neg) (cylinder0) // 36
    project ($x4 $r0    $r0Neg) (cylinder0) // 37
    project ($x4 $r0Neg $r0) (cylinder0)    // 38
    project ($x4 $r0    $r0) (cylinder0)    // 39
    // x = x5
    // Inner
    ($x5 $box5Neg $box5Neg) // 40
    ($x5 $box5    $box5Neg) // 41
    ($x5 $box5Neg $box5)    // 42
    ($x5 $box5    $box5)    // 43
    // Outer block points
    project ($x5 $r0Neg $r0Neg) (cylinder0) // 44
    project ($x5 $r0    $r0Neg) (cylinder0) // 45
    project ($x5 $r0Neg $r0) (cylinder0)    // 46
    project ($x5 $r0    $r0) (cylinder0)    // 47
);

blocks
(
    // x0 <= x <= x1
    hex (4 12 8 0 6 14 10 2) ($nX1 $nR $nR) simpleGrading (1 $gd 1)      // 0
    hex (6 14 10 2 7 15 11 3) ($nX1 $nR $nR) simpleGrading (1 $gd 1)     // 1
    hex (7 15 11 3 5 13 9 1) ($nX1 $nR $nR) simpleGrading (1 $gd 1)      // 2
    hex (5 13 9 1 4 12 8 0) ($nX1 $nR $nR) simpleGrading (1 $gd 1)       // 3
    hex (0 8 9 1 2 10 11 3) ($nX1 $nR $nR) simpleGrading (1 1 1)         // 4
    // x1 <= x <= x2
    hex (12 20 16 8 14 22 18 10) ($nX2 $nR $nR) simpleGrading (1 $gd 1)  // 5
    hex (14 22 18 10 15 23 19 11) ($nX2 $nR $nR) simpleGrading (1 $gd 1) // 6
    hex (15 23 19 11 13 21 17 9) ($nX2 $nR $nR) simpleGrading (1 $gd 1)  // 7
    hex (13 21 17 9 12 20 16 8) ($nX2 $nR $nR) simpleGrading (1 $gd 1)   // 8
    hex (8 16 17 9 10 18 19 11) ($nX2 $nR $nR) simpleGrading (1 1 1)     // 9
    // x2 <= x <= x3
    hex (20 28 24 16 22 30 26 18) ($nX3 $nR $nR) simpleGrading (1 $gd 1) // 10
    hex (22 30 26 18 23 31 27 19) ($nX3 $nR $nR) simpleGrading (1 $gd 1) // 11
    hex (23 31 27 19 21 29 25 17) ($nX3 $nR $nR) simpleGrading (1 $gd 1) // 12
    hex (21 29 25 17 20 28 24 16) ($nX3 $nR $nR) simpleGrading (1 $gd 1) // 13
    hex (16 24 25 17 18 26 27 19) ($nX3 $nR $nR) simpleGrading (1 1 1)   // 14
    // x3 <= x <= x4
    hex (28 36 32 24 30 38 34 26) ($nX4 $nR $nR) simpleGrading (1 $gd 1) // 15
    hex (30 38 34 26 31 39 35 27) ($nX4 $nR $nR) simpleGrading (1 $gd 1) // 16
    hex (31 39 35 27 29 37 33 25) ($nX4 $nR $nR) simpleGrading (1 $gd 1) // 17
    hex (29 37 33 25 28 36 32 24) ($nX4 $nR $nR) simpleGrading (1 $gd 1) // 18
    hex (24 32 33 25 26 34 35 27) ($nX4 $nR $nR) simpleGrading (1 1 1)   // 19
    // x4 <= x <= x5
    hex (36 44 40 32 38 46 42 34) ($nX5 $nR $nR) simpleGrading (1 $gd 1) // 20
    hex (38 46 42 34 39 47 43 35) ($nX5 $nR $nR) simpleGrading (1 $gd 1) // 21
    hex (39 47 43 35 37 45 41 33) ($nX5 $nR $nR) simpleGrading (1 $gd 1) // 22
    hex (37 45 41 33 36 44 40 32) ($nX5 $nR $nR) simpleGrading (1 $gd 1) // 23
    hex (32 40 41 33 34 42 43 35) ($nX5 $nR $nR) simpleGrading (1 1 1)   // 24
);

edges
(
    // x = x0
    project  4  5 (cylinder0)
    project  7  5 (cylinder0)
    project  6  7 (cylinder0)
    project  4  6 (cylinder0)
    // x = x1
    project 12 13 (cylinder0)
    project 13 15 (cylinder0)
    project 12 14 (cylinder0)
    project 14 15 (cylinder0)
    // x = x2
    project 20 21 (cylinder1)
    project 21 23 (cylinder1)
    project 20 22 (cylinder1)
    project 22 23 (cylinder1)
    // x = x3
    project 28 29 (cylinder1)
    project 29 31 (cylinder1)
    project 28 30 (cylinder1)
    project 30 31 (cylinder1)
    // x = x4
    project 36 37 (cylinder0)
    project 39 37 (cylinder0)
    project 38 39 (cylinder0)
    project 36 38 (cylinder0)
    // x = x5
    project 44 45 (cylinder0)
    project 47 45 (cylinder0)
    project 46 47 (cylinder0)
    project 44 46 (cylinder0)
);

defaultPatch
{
    name tube;
    type wall;
}

boundary
(
    inlet
    {
        type    patch;
        faces
        (
            (4 5 1 0)
            (4 0 2 6)
            (1 5 7 3)
            (2 3 7 6)
            (0 1 3 2)
        );
    }
    outlet
    {
        type    patch;
        faces
        (
            (44 45 41 40)
            (44 40 42 46)
            (41 45 47 43)
            (42 43 47 46)
            (40 41 43 42)
        );
    }
);

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

解析条件

本解析は simpleFoam を使用した定常計算である。
主な解析条件を下表に示す。
テキストによれば、速度分布の解決のため、ベンチュリ管壁に最も近いセルの高さは$1.5(mm)$以下とする必要がある。
径方向のセルサイズの最小値を示すセルは管壁に最も近いセルであり、この条件を満たしている。

項■目 ■■■■■■設■定■内■容■■■■■■
使用ソフトウェアとバージョン OpenFOAM v10
ソルバ simpleFoam
総セル数 $57,600$
径方向セルサイズ $1.45×10^{-3}〜7.69×10^{-3} (m)$
軸方向セルサイズ $6.94×10^{-3}〜8.34×10^{-3} (m)$
ベンチュリ管入口直径$D$ $0.1 (m)$
ベンチュリ管入口最大流速$|\mathbf{u_{max}}|$ $0.4 (m)$
動粘性係数$\nu$ $4×10^{-5} (m^{2}/s)$
計算収束条件$r_{abs}$ $1×10^{-4} (-)$

解析条件に関係するファイルを以下に示す。

system/controlDict

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     simpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         500;

deltaT          1;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

functions
{
    probes
    {
        type            probes;
        libs            ("libsampling.so");
        writeControl    writeTime;
        probeLocations
        (
            (0.05 0 0)
            (0.25 0 0)
        );
        fixedLocations  false;
        fields
        (
            p
        );
    }
}

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

system/fvSchemes

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}


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

system/fvSolution
SIMPLE/residualControlで指定した値で計算の収束を判定する。

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0.01;
    }

    U
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-8;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent no; // using SIMPLE method
    residualControl
    {
        p               1e-4;
        U               1e-4;
    }
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
    }
}

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

control/transportProperties
動粘性係数$\nu$を設定する。

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

transportModel  Newtonian;

nu              4E-05;

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

control/momentumTransport
シミュレーションタイプは層流とする。

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      momentumTransport;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType laminar;

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

境界条件

境界条件を下表に示す。

  • 流入境界$inlet$の流速$U$は、前述のハーゲン・ポアズイユ流れを模擬するために、$codedFixedValue$を使用している。
  • 流出境界$outlet$では、スロートの下流の壁近傍で再循環する流れにより、逆流が発生するため、流速$U$および圧力$p$の両方で、逆流の計算ができる境界条件としている。
パッチ $U\ $ $p\ $
流入境界$inlet$ $codedFixedValue$ $zeroGradient$
流出境界$outlet$ $pressureInletOutletVelocity$ $totalPressure$
ベンチュリ管壁$tube$ $noSlip\ $ $zeroGradient$

以下に実際の境界条件のファイルを示す。
0/U

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            codedFixedValue;
        value           uniform (0 0 0);
        name            quadratic;
        code
        #{
            const fvPatch& boundaryPatch = patch();
            const vectorField& Cf = boundaryPatch.Cf();
            vectorField& field = *this;
            
            const scalar r0 = 0.05;
            const scalar uMax = 0.4;
            
            forAll(Cf, faceI)
            {
                scalar r = pow((pow(Cf[faceI].y(), 2.0) + pow(Cf[faceI].z(), 2.0)), 0.5);
                scalar u = uMax * (1.0 - pow(r / r0, 2.0));
                field[faceI] = vector(u, 0, 0);
            }
        #};
    }

    outlet
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }

    tube
    {
        type            noSlip;
    }

    defaultFaces
    {
        type            empty;
    }
}

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

0/p

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            totalPressure;
        p0              uniform 0;
    }

    tube
    {
        type            zeroGradient;
    }

    defaultFaces
    {
        type            empty;
    }
}

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

解析結果

速度分布

計算収束時の速度分布を下図に示す。
境界条件で想定した通り、ベンチュリ管入口の速度分布は、最大流速を頂点とした曲面が設定され、スロートより下流では逆流が起こっていることが分かる。なお、図中の流速は最大流速である。
fig3_simple.png

テキストのサンプル問題との比較

テキストのサンプル問題と本解析の計算結果の比較を下表に示す。
テキストのサンプル問題と本解析の計算結果には大きな差はなく、本解析はサンプル問題をよく再現していると言える。

■■■■ サンプル問題 本解析
反復回数 $292$ $225$
$\Delta p\ $ $0.434$ $0.420$
$C\ $ $0.831$ $0.845$

まとめ

  • Notes on Computational Fluid Dynamics: General Principles に記載されているサンプル問題の結果をよく再現できた。
  • 「円柱周りの流れ」と同様に、同書にはメッシュの作り方の記載が少なく、メッシュを作成するのに苦労した。
  • fvSolutions/SIMPLE/consistentをyesとすると、SIMPLE法からSIMPLEC法に切り替えることができる。
  • SIMPLEC法では収束性が改善し、反復回数が132回と少なくなるが、計算結果はSIMPLE法での計算結果とほぼ変わらない。
1
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2