円柱周りの流れ解析

# はじめに
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

## 参考にしたサイト
- [Notes on Computational Fluid Dynamics: General Principles](https://doc.cfd.direct/notes/cfd-general-principles/flow-around-a-cylinder)
- [Tutorial Guide 2.2 flow-around-a-cylinder](https://www.openfoam.com/documentation/tutorial-guide/2-incompressible-flow/2.2-flow-around-a-cylinder)

# 幾何形状
円柱周りの流れ解析の幾何形状は下図の通りである。
図中の _D_ は円柱の直径である。
![fig1_geometry.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/7883e4a3-f9e5-09b2-3110-57e3e2ca0a95.png)

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

![fig2_block.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/eb2c4401-61fd-6555-58e1-1586b60ca14a.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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/fd66bc9a-b6da-0c48-7eb1-cd51a861141b.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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/f1081909-2fc4-26ae-b466-bbf41544fd08.png)

## _U_
速度分布のアニメーションを以下に示す。
![animation.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/cd4c145a-d8c7-07a0-c217-055acd92f084.gif)

## _CL_
揚力係数 _CL_ の時刻歴を下図に示す。
![fig6_CL.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/02aff9bf-667f-1b0e-4d0a-fb82d48dfbdc.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されているものとする」としか記載がなく、メッシュを再現するのに苦労した。