LoginSignup
3
0

More than 1 year has passed since last update.

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

円柱周りの流れ解析

Posted at

はじめに

Notes on Computational Fluid Dynamics: General Principlesの第8章のサンプル問題に記載されている円柱周りの流れ(Flow Around A Cylinder)解析を行った。
OpenFOAMは2022年7月時点で最新のFoundation版バージョン10を使用した。

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

  • $FOAM_TUTORIALS/basic/potentialFoam/cylinder
  • $FOAM_TUTORIALS/incompressible/pimpleFoam/laminar/offsetCylinder

参考にしたサイト

幾何形状

円柱周りの流れ解析の幾何形状は下図の通りである。
図中の D は円柱の直径である。
fig1_geometry.png

メッシュ生成

メッシュはblockMeshで生成する。メッシュは下図のように20ブロックで構成される。

fig2_block.png

system/blockMeshDict は次の通りである。

system/blockMeshDict

/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices #codeStream
{
    codeInclude
    #{
        #include "pointField.H"
    #};

    code
    #{
        pointField points(32);
        points[0]  = point(0.001, 0, -0.5);
        points[1]  = point(0.01, 0, -0.5);
        points[2]  = point(0.1, 0, -0.5);
        points[3]  = point(0.1, 7.07107E-03, -0.5);
        points[4]  = point(7.07107E-03, 7.07107E-03, -0.5);
        points[5]  = point(7.07107E-04, 7.07107E-04, -0.5);
        points[6]  = point(0.1, 0.03, -0.5);
        points[7]  = point(7.07107E-03, 0.03, -0.5);
        points[8]  = point(0, 0.03, -0.5);
        points[9]  = point(0, 0.01, -0.5);
        points[10] = point(0, 0.001, -0.5);

        points[11] = point(-0.001, 0, -0.5);
        points[12] = point(-0.01, 0, -0.5);
        points[13] = point(-0.03, 0, -0.5);
        points[14] = point(-0.03, 7.07107E-03, -0.5);
        points[15] = point(-7.07107E-03, 7.07107E-03, -0.5);
        points[16] = point(-7.07107E-04, 7.07107E-04, -0.5);
        points[17] = point(-0.03, 0.03, -0.5);
        points[18] = point(-7.07107E-03, 0.03, -0.5);

        points[19] = point(0.1, -7.07107E-03, -0.5);
        points[20] = point(7.07107E-03, -7.07107E-03, -0.5);
        points[21] = point(7.07107E-04, -7.07107E-04, -0.5);
        points[22] = point(0.1, -0.03, -0.5);
        points[23] = point(7.07107E-03, -0.03, -0.5);
        points[24] = point(0, -0.03, -0.5);
        points[25] = point(0, -0.01, -0.5);
        points[26] = point(0, -0.001, -0.5);

        points[27] = point(-0.03, -7.07107E-03, -0.5);
        points[28] = point(-7.07107E-03, -7.07107E-03, -0.5);
        points[29] = point(-7.07107E-04, -7.07107E-04, -0.5);
        points[30] = point(-0.03, -0.03, -0.5);
        points[31] = point(-7.07107E-03, -0.03, -0.5);

        // Duplicate z points
        label sz = points.size();
        points.setSize(2*sz);
        for (label i = 0; i < sz; i++)
        {
            const point& pt = points[i];
            points[i+sz] = point(pt.x(), pt.y(), -pt.z());
        }

        os  << points;
    #};
};

n1 40;
n2 40;
n3 90;
gd1 64;
gd2 2.55;

blocks
(
    hex (5 4 9 10 37 36 41 42) ($n1 $n1 1) simpleGrading ($gd1 1 1)    // 0
    hex (0 1 4 5 32 33 36 37) ($n1 $n1 1) simpleGrading ($gd1 1 1)     // 1
    hex (1 2 3 4 33 34 35 36) ($n3 $n1 1) simpleGrading ($gd2 1 1)     // 2
    hex (4 3 6 7 36 35 38 39) ($n3 $n2 1) simpleGrading ($gd2 1 1)     // 3
    hex (9 4 7 8 41 36 39 40) ($n1 $n2 1) simpleGrading (1 1 1)        // 4

    hex (10 9 15 16 42 41 47 48) ($n1 $n1 1) simpleGrading ($gd1 1 1)  // 5
    hex (16 15 12 11 48 47 44 43) ($n1 $n1 1) simpleGrading ($gd1 1 1) // 6
    hex (13 12 15 14 45 44 47 46) ($n2 $n1 1) simpleGrading (1 1 1)    // 7
    hex (14 15 18 17 46 47 50 49) ($n2 $n2 1) simpleGrading (1 1 1)    // 8
    hex (15 9 8 18 47 41 40 50) ($n1 $n2 1) simpleGrading (1 1 1)      // 9

    hex (26 25 20 21 58 57 52 53) ($n1 $n1 1) simpleGrading ($gd1 1 1) // 10
    hex (21 20 1 0 53 52 33 32) ($n1 $n1 1) simpleGrading ($gd1 1 1)   // 11
    hex (20 19 2 1 52 51 34 33) ($n3 $n1 1) simpleGrading ($gd2 1 1)   // 12
    hex (23 22 19 20 55 54 51 52) ($n3 $n2 1) simpleGrading ($gd2 1 1) // 13
    hex (24 23 20 25 56 55 52 57) ($n1 $n2 1) simpleGrading (1 1 1)    // 14

    hex (11 12 28 29 43 44 60 61) ($n1 $n1 1) simpleGrading ($gd1 1 1) // 15
    hex (29 28 25 26 61 60 57 58) ($n1 $n1 1) simpleGrading ($gd1 1 1) // 16
    hex (27 28 12 13 59 60 44 45) ($n2 $n1 1) simpleGrading (1 1 1)    // 17
    hex (30 31 28 27 62 63 60 59) ($n2 $n2 1) simpleGrading (1 1 1)    // 18
    hex (31 24 25 28 63 56 57 60) ($n1 $n2 1) simpleGrading (1 1 1)    // 19
);

edges
(
    arc 0 5 (9.39693E-04 3.82683E-04 -0.5)
    arc 5 10 (3.82683E-04 9.39693E-04 -0.5)
    arc 1 4 (9.23880E-03 3.82683E-03 -0.5)
    arc 4 9 (3.82683E-03 9.23880E-03 -0.5)
    arc 32 37 (9.39693E-04 3.82683E-04 0.5)
    arc 37 42 (3.82683E-04 9.39693E-04 0.5)
    arc 33 36 (9.23880E-03 3.82683E-03 0.5)
    arc 36 41 (3.82683E-03 9.23880E-03 0.5)

    arc 11 16 (-9.39693E-04 3.82683E-04 -0.5)
    arc 16 10 (-3.82683E-04 9.39693E-04 -0.5)
    arc 12 15 (-9.23880E-03 3.82683E-03 -0.5)
    arc 15 9 (-3.82683E-03 9.23880E-03 -0.5)
    arc 43 48 (-9.39693E-04 3.82683E-04 0.5)
    arc 48 42 (-3.82683E-04 9.39693E-04 0.5)
    arc 44 47 (-9.23880E-03 3.82683E-03 0.5)
    arc 47 41 (-3.82683E-03 9.23880E-03 0.5)

    arc 21 0 (9.39693E-04 -3.82683E-04 -0.5)
    arc 26 21 (3.82683E-04 -9.39693E-04 -0.5)
    arc 20 1 (9.23880E-03 -3.82683E-03 -0.5)
    arc 25 20 (3.82683E-03 -9.23880E-03 -0.5)
    arc 53 32 (9.39693E-04 -3.82683E-04 0.5)
    arc 58 53 (3.82683E-04 -9.39693E-04 0.5)
    arc 52 33 (9.23880E-03 -3.82683E-03 0.5)
    arc 57 52 (3.82683E-03 -9.23880E-03 0.5)

    arc 11 29 (-9.39693E-04 -3.82683E-04 -0.5)
    arc 29 26 (-3.82683E-04 -9.39693E-04 -0.5)
    arc 12 28 (-9.23880E-03 -3.82683E-03 -0.5)
    arc 28 25 (-3.82683E-03 -9.23880E-03 -0.5)
    arc 43 61 (-9.39693E-04 -3.82683E-04 0.5)
    arc 61 58 (-3.82683E-04 -9.39693E-04 0.5)
    arc 44 60 (-9.23880E-03 -3.82683E-03 0.5)
    arc 60 57 (-3.82683E-03 -9.23880E-03 0.5)
);

defaultPatch
{
    type empty;
}

boundary
(
    down
    {
        type wall;
        faces
        (
            (24 23 55 56)
            (23 22 54 55)
            (31 24 56 63)
            (30 31 63 62)
        );
    }
    right
    {
        type patch;
        faces
        (
            (2 3 35 34)
            (3 6 38 35)
            (19 2 34 51)
            (22 19 51 54)
        );
    }
    up
    {
        type wall;
        faces
        (
            (7 8 40 39)
            (6 7 39 38)
            (8 18 50 40)
            (18 17 49 50)
        );
    }
    left
    {
        type patch;
        faces
        (
            (14 13 45 46)
            (17 14 46 49)
            (13 27 59 45)
            (27 30 62 59)
        );
    }
    cylinder
    {
        type wall;
        faces
        (
            (10 5 37 42)
            (5 0 32 37)
            (16 10 42 48)
            (11 16 48 43)
            (0 21 53 32)
            (21 26 58 53)
            (26 29 61 58)
            (29 11 43 61)
        );
    }
);

mergePatchPairs
(
);

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

実際に生成したメッシュを下図に示す。

  • 円柱の径方向にgradingを設定することで、円柱近傍のセル高さを 0.008 D 以下とした。
  • 円柱後流の x 方向にgradingを設定することで、総メッシュを40,000とした。
    fig3_model.png

解析条件

主な解析条件を下表に示す。

項目 設定内容
使用ソフトウェアとバージョン OpenFOAM ver.10
ソルバ pimpleFoam
総メッシュ数 40,000
円柱直径 D 2E-03 (m)
円柱近傍 dx 0.007～0.462 (D)
円柱近傍 dy 0.010 (D)
その他 dx 0.088～0.794 (D)
その他 dy 0.088～0.287 (D)
時間刻み dt 2.5E-05 (s)
解析終了時間 0.4 (s)
自由流れ速度 u∞ 1 (m/s)
動粘性係数 ν 1E-05 (m2/s)

system/controlDict と揚力係数CLを出力するための system/forceCoeffs 、y+を出力するための system/yPlus を以下に示す。

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     pimpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.4;

deltaT          2.5e-5;

writeControl    runTime;

writeInterval   0.001;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

functions
{
    #include "forceCoeffs"
    #include "yPlus"
}

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

system/forceCoeffs

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/

forceCoeffs1
{
    type            forceCoeffs;

    libs            ("libforces.so");

    writeControl    timeStep;
    timeInterval    1;

    log             no;

    patches         (cylinder);
    rho             rhoInf;      // Indicates incompressible
    rhoInf          1;           // Redundant for incompressible
    liftDir         (0 1 0);
    dragDir         (1 0 0);
    CofR            (0 0 0);     // Axle midpoint on ground
    pitchAxis       (0 0 1);
    magUInf         1;
    lRef            2e-3;        // Diameter of cylinder
    Aref            2e-3;        // Estimated
    /*
    binData
    {
        nBin        20;          // output data into 20 bins
        direction   (1 0 0);     // bin direction
        cumulative  yes;
    }
    */
}


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

system/yPlus

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
-------------------------------------------------------------------------------
Description
    Calculates the turbulence y+, outputting the data as a yPlus field.

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

yPlus1
{
    type            yPlus;
    libs            ("libfieldFunctionObjects.so");

    executeControl  writeTime;
    writeControl    writeTime;
}

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

system/fvSchmes および system/fvSolution は次の通りである。

system/fvSchmes

/*--------------------------------*- 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         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

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

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}


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

system/fvSolution

/*--------------------------------*- 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-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    "U.*"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-5;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor   no;
    nOuterCorrectors    5;
    nCorrectors         1;
    nNonOrthogonalCorrectors 0;
}

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

動粘性係数 ν の設定は、control/transportProperties で行う。

control/transportProperties

/*--------------------------------*- 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              1E-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;

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

境界条件

  • 流入(left)および流出(right)の流速は x 方向に 1 (m/s)、圧力は 0 (Pa) の自由流れ境界条件とする。
  • 上部(up)、下部(down)および円柱(cylinder)の流速はスリップなし条件、圧力は勾配なし条件とする。

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 (1 0 0);

boundaryField
{
    left
    {
        type            freestreamVelocity;
        freestreamValue $internalField;
        value           $internalField;
    }

    right
    {
        type            freestreamVelocity;
        freestreamValue $internalField;
        value           $internalField;
    }

    down
    {
        type            noSlip;
    }

    up
    {
        type            noSlip;
    }

    cylinder
    {
        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
{
    left
    {
        type            freestreamPressure;
        freestreamValue $internalField;
        value           $internalField;
    }

    right
    {
        type            freestreamPressure;
        freestreamValue $internalField;
        value           $internalField;
    }

    down
    {
        type            zeroGradient;
    }

    up
    {
        type            zeroGradient;
    }

    cylinder
    {
        type            zeroGradient;
    }

    defaultFaces
    {
        type            empty;
    }
}

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

解析結果

y+

計算全体を通して、円柱の y+ は y+ < 1 を維持できていた。
これは本解析の必要条件である（らしい）。
fig4_yPlusCylinder.png

U

速度分布のアニメーションを以下に示す。
animation.gif

CL

揚力係数 CL の時刻歴を下図に示す。
fig6_CL.png

振動は t = 0.08 (s) 頃に始まる。t = 0.15 (s) 頃に渦のパターンが安定するため、CL の振幅および周期(振動数)も安定した。
CL の振幅および周期(周波数)が安定する t = 0.15(s)からt = 0.4(s)までのデータを用いて、LibreOffice calc でフーリエ変換を行った。フーリエ変換により CL の周波数が95.99(Hz)と求まり、周期 T は0.0104(s) であった。
ストローハル数は、Sr = D / (T u∞) = 0.192 でほぼ0.2であった。

まとめ

  • Notes on Computational Fluid Dynamics: General Principles に記載されているサンプル問題の結果をよく再現できた。
  • 同書には「円柱周りは適切にgradingされているものとする」としか記載がなく、メッシュを再現するのに苦労した。
3
0
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
3
0