長方形箱内で移動する正方形周りの2次元非圧縮性流れ解析

# はじめに
ESI版 (OpenCFD版) OpenFOAM に搭載されているオーバーセットメッシュの機能を使って、「[SPHERIC Test06 2-D Incompressible flow around a moving square inside a rectangular box](https://www.spheric-sph.org/tests/test-06)」の再現解析を行います。OpenFOAM のバージョンは、v2306 とします。
解析方法などは、オープンCAE学会と関連が深い「改訂新版 OpenFOAM の歩き方 川畑真一著 (株式会社インプレスR&D)」を参考にしています。

# 解析モデル

![fig01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/259664de-d531-bbb6-a6d8-d97aea223702.png "図1 解析モデル")
図1  解析モデル

<table>
  <caption>表1 解析モデル</caption>
  <thead>
    <tr>
      <th>名前</th> <th>記号</th> <th>設定値</th>
    </tr>
  </thead>
  <tr>
    <td rowspan="2">タンク (長方形箱)</td> <td>$L_x$ [m]</td> <td> 10 </td>
  </tr>
  <tr>
    <td>$L_y$ [m]</td> <td> 5 </td>
  </tr>
  <tr>
    <td rowspan="2">移動する正方形</td> <td>一辺の長さ $L$ [m]</td> <td> 1 </td>
  </tr>
  <tr>
    <td>最大速度 $U_{max}$ [m/s]</td> <td> 1 </td>
  </tr>
</table>

# 解析条件
移動する正方形のレイノルズ数$Re$がそれぞれ 50、100 および 150 のケースを計算します。
$Re$ の設定は、$Re$ の定義式 $Re=U_{max}L/ \nu$ から $\nu$ を介して行います。
<table>
  <caption>表2 解析条件1</caption>
  <thead>
    <tr>
      <th>$Re$</th> <th>$\nu [m^2/s]$</th>
    </tr>
  </thead>
  <tr>
    <td> 50 </td> <td> 0.02 </td>
  </tr>
  <tr>
    <td> 100 </td> <td> 0.01 </td>
  </tr>
  <tr>
    <td> 150 </td> <td> 6.67E-03 </td>
  </tr>
</table>

解析メッシュは、粗い、標準および詳細の3つを用意します。
<table>
  <caption>表3 解析条件2</caption>
  <thead>
    <tr>
      <th> 解析メッシュ </th> <th colspan="2"> タンク (長方形箱) </th> <th> 移動する正方形 </th>
    </tr>
  </thead>
  <tr>
    <td> 分割数 </td> <td> $n_x$ </td> <td> $n_y$ </td> <td> $n$ </td>
  </tr>
  <tr>
    <td> 粗い (coarse) </td> <td> 267 </td> <td> 133 </td> <td> 53 </td>
  </tr>
  <tr>
    <td> 標準 (standard) </td> <td> 400 </td> <td> 200 </td> <td> 80 </td>
  </tr>
  <tr>
    <td> 詳細 (fine) </td> <td> 600 </td> <td> 300 </td> <td> 120 </td>
  </tr>
</table>

表2に示したレイノルズ数 $Re$ と表3に示した解析メッシュを組み合わせ、9ケースの解析を行います。

# 前処理 (解析準備)
解析は、次のチュートリアルを流用して行います。
$FOAM_TUTORIALS/incompressible/overPimpleDyMFoam/simpleRotor

## system/blockMesh
vertices ブロックを表1に合わせて変更します。
上のデータがタンクに、下のデータが移動する正方形に該当します。
blocks ブロックを表3に合わせてそれぞれ変更します。
こちらも、上の行がタンクに、下の行が移動する正方形に該当します。

```cpp:system/blockMesh
scale   1;

vertices
(
    (  0  0  0 )
    ( 10  0  0 )
    ( 10  5  0 )
    (  0  5  0 )
    (  0  0  1 )
    ( 10  0  1 )
    ( 10  5  1 )
    (  0  5  1 )

    ( 0.5  1.5  0 )
    ( 2.5  1.5  0 )
    ( 2.5  3.5  0 )
    ( 0.5  3.5  0 )
    ( 0.5  1.5  1 )
    ( 2.5  1.5  1 )
    ( 2.5  3.5  1 )
    ( 0.5  3.5  1 )
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (600 300 1) simpleGrading (1 1 1)

    hex (8 9 10 11 12 13 14 15) movingZone (120 120 1) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    overset2
    {
        type        overset;
        faces
        (
            ( 8 12 15 11)
            (10 14 13  9)
            (11 15 14 10)
            ( 9 13 12  8)
        );
    }

    walls
    {
        type wall;
        faces
        (
            (3 7 6 2)
            (1 5 4 0)
            (0 4 7 3)
            (2 6 5 1)
        );
    }
    // Populated by subsetMesh
    hole
    {
        type wall;
        faces ();
    }

    frontAndBack
    {
        type empty;
        faces
        (
            (0 3 2 1)
            (4 5 6 7)
            ( 8 11 10  9)
            (12 13 14 15)
        );
    }
);
```
## system/controlDict
endTime を 8、writeInterval を 0.1 にそれぞれ変更しています。
また、計算結果を評価するために、probe で圧力基準点 p(0,0,0) の値を、forces で抗力を、vorticity で渦度を function object で取得します。そのために、functions ブロックを追加しています。

```cpp:system/controlDict
libs            (overset fvMotionSolvers);

application     overPimpleDyMFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         8;
//endTime         0.1;

deltaT          0.00025;

//writeControl    timeStep;
//writeInterval   10;

writeControl    adjustable;
writeInterval   0.1;
//writeInterval   0.01;

purgeWrite      0;

writeFormat     ascii;

writePrecision  10;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  yes;

maxCo           1;

functions
{
    probes
    {
        type            probes;
        libs            (sampling);
        writeControl    adjustable;
        writeInterval   0.1;
        writeOption     anyWrite;
        fields          (p);
        probeLocations
        (
            (0.0001 0.0001 0.001)
        );
    }
    forces1
    {
        type            forces;
        libs            (forces);
        patches         (hole);
        rho             rhoInf;
        rhoInf          1;
        CofR            (0 0 0);
        writeControl    timeStep;
        writeInterval   1;
        porosity        no;
    }
    vorticity1
    {
        type        vorticity;
        libs        (fieldFunctionObjects);
        writeControl        adjustable;
        writeInterval       0.1;
        writeOption         anyWrite;
    }
}
```

## system/decomposeParDict
method を scotch に変更しています。

```cpp:system/decomposeParDict
numberOfSubdomains  3;

method          scotch;
//method          hierarchical;

coeffs
{
    n           (3 1 1);
}
```
## system/fvSchemes
oversetInterpolation ブロックの「method inverseDistance;」以外はコメントにします。
コメントにしないと正常に計算ができません。

```cpp:system/fvSchemes
ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(T)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss limitedLinearV 1;
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,epsilon) Gauss limitedLinear 1;
    div(phi,R)      Gauss limitedLinear 1;
    div(R)          Gauss linear;
    div(phi,nuTilda) Gauss limitedLinear 1;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
    laplacian(diffusivity,cellDisplacement)  Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

oversetInterpolation
{
    method              inverseDistance;
/*
    searchBox           (0 0 0)(0.1 0.1 0.1);
    searchBoxDivisions  (130 130 1);
    holeLayers          4;
    useLayer            2;
*/
}

fluxRequired
{
    default         no;
    pcorr           ;
    p               ;
}
```

## system/fvSolution
計算時に抗力にスパイクのような突変が見られ、そのタイミングでは圧力が収束していない (計算反復回数が上限に達している)ことがわかったため、次の対策を行いました。
・momentumPredictor を有効にする (no から yes にする)
・nCorrectors を2から3に増やす

```cpp:system/fvSolution
solvers
{
    cellDisplacement
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0;
        maxIter         100;
    }

    p
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0;
    }

    pFinal
    {
        $p;
    }
    pcorr
    {
        $p;
        solver          PCG;
        preconditioner  DIC;
    }

    pcorrFinal
    {
        $pcorr;
        relTol          0;
    }

    "(U|k|epsilon)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0;
    }

    "(U|k|epsilon)Final"
    {
        $U;
        tolerance       1e-6;
        relTol          0;
    }
}
PIMPLE
{
    momentumPredictor   yes;
    //momentumPredictor   no;
    nOuterCorrectors    1;
    nCorrectors         3;
    //nCorrectors         2;
    nNonOrthogonalCorrectors 0;

    pRefPoint           (0.0001 0.0001 0.001);
    pRefValue           0.0;
}

relaxationFactors
{
   fields
   {
   }
   equations
   {
       ".*"             1;
   }
}
```
## system/setFieldDict
変更する必要はありません。

```cpp:system/setFieldDict
defaultFieldValues
(
    volScalarFieldValue zoneID 123
);

regions
(
    // Set cell values
    // (does zerogradient on boundaries)
    cellToCell
    {
        set c0;

        fieldValues
        (
            volScalarFieldValue zoneID 0
        );
    }

    cellToCell
    {
        set c1;

        fieldValues
        (
            volScalarFieldValue zoneID 1
        );
    }
);
```
## system/topoSetDict
表1に合わせて cellSet を変更します。

```cpp:system/topoSetDict
actions
(
    {
        name    c0;
        type    cellSet;
        action  new;
        source  regionToCell;
        insidePoints ((0.001 0.001 0.001));
        //insidePoints ((0.01 0.01 0.01));
    }

    {
        name    c1;
        type    cellSet;
        action  new;
        source  cellToCell;
        set     c0;
    }

    {
        name    c1;
        type    cellSet;
        action  invert;
    }
    // Select box to remove from region 1

    {
        name    box;
        type    cellSet;
        action  new;
        source  cellToCell;
        set     c1;
    }

    {
        name    box;
        type    cellSet;
        action  subset;
        source  boxToCell;
        box     (1 2 -100)(2 3 100);
        //box     (0.025 0.045 -100)(0.075 0.055 100);
    }

    {
        name    box;
        type    cellSet;
        action  invert;
    }
);
```
## constant/turbulenceProperties
変更する必要はありません。

```cpp:constant/turbulenceProperties
simulationType laminar;
```
## constant/transportProperties
「transportModel  Newtonian;」は変更しません。nu の値は表2のように変更します。
この2つ以外はコメントにします。

```cpp:constant/transportProperties
//DT              1;

transportModel  Newtonian;

nu              0.02;  // Re=50, L=1, Umax=1
//nu              1e-05;
/*
CrossPowerLawCoeffs
{
    nu0         1e-06;
    nuInf       1e-06;
    m           1;
    n           1;
}

BirdCarreauCoeffs
{
    nu0         1e-06;
    nuInf       1e-06;
    k           0;
    n           1;
}
*/
```
## constant/dynamicMeshDict
移動する正方形の動き(変位)のデータが提供されているため、solidBodyMotionFunctionをrotatingMotion から tabulated6DoFMotion に変更します。

```cpp:constant/transportProperties
dynamicFvMesh       dynamicOversetFvMesh;

solver          multiSolidBodyMotionSolver;

multiSolidBodyMotionSolverCoeffs
{
    movingZone
    {
        solidBodyMotionFunction tabulated6DoFMotion;
        tabulated6DoFMotionCoeffs
        {
            timeDataFileName "<constant>/6DoF.dat";
            CofG (1.5 2.5 0.5);
        }
    /*
        solidBodyMotionFunction rotatingMotion;
        rotatingMotionCoeffs
        {
            origin      (0.05 0.05 0.05);
            axis        (0 0 1);
            omega       100.0;
        }
    */
    }
}
```
## constant/6DoF.dat
dynamicMeshDict で参照される、移動する正方形の動きのデータです。
提供された Motion_Body.dat を基に、dynamicMesh で読み込めるよう変換して、6DoF.dat という名前で保存したものです。
Motion_Body.dat のフォーマットは、1列目が時間[$s$]、2列目が加速度[$m^2/s$]、3列目が位置[$m$]です。
6DoF.dat のフォーマットは、1列目が時間[$s$]、2列目が変位ベクトル[$m$]、3列目が回転量ベクトル[$deg$]です。Motion_Body.dat の初期位置からの位置の差を変位として、6DoF.dat の変位ベクトルのx成分に設定します。

## 0.orig/U
変更する必要はありません。
```cpp:0.orig/U
dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    walls
    {
        type            uniformFixedValue;
        uniformValue    (0 0 0);
    }

    hole
    {
        type            movingWallVelocity;
        value           uniform (0 0 0);
    }

    overset
    {
        type            overset;
    }
}
```
## 0.orig/p
変更する必要はありません。
```cpp:0.orig/p
dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    "(walls|hole)"
    {
        type            zeroGradient;
    }

    overset
    {
        type            overset;
    }
}
```
## 0.orig/zoneID
変更する必要はありません。
```cpp:0.orig/zoneID
dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    ".*"
    {
        type            zeroGradient;
    }

    overset
    {
        type            overset;
        value           uniform 0;
    }
}
```

## 0.orig/その他
U、p、zoneID 以外のファイルは使用しません。

# 計算実行
シェルスクリプトを使用して計算を実行します。

## Allrun
計算の実行で使用するシェルスクリプトです。
変更する必要はありません。
```sh:Allrun
#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
#------------------------------------------------------------------------------

./Allrun.pre

runApplication decomposePar -cellDist

runParallel $(getApplication)

runApplication reconstructPar

#------------------------------------------------------------------------------

```
## Allrun.pre
計算の実行で Allrun から呼び出されるシェルスクリプトです。
変更する必要はありません。
```sh:Allrun.pre
#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
#------------------------------------------------------------------------------

runApplication blockMesh

# Select cellSets
runApplication topoSet

runApplication subsetMesh box -patch hole -overwrite

# Select cellSets
runApplication -s box topoSet

restore0Dir

# Use cellSets to write zoneID
runApplication setFields

#------------------------------------------------------------------------------
```
## Allclean
計算結果を削除するときに使用するシェルスクリプトです。
変更する必要はありません。
```sh:Allclean
#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
#------------------------------------------------------------------------------

cleanCase0

#------------------------------------------------------------------------------
```

# 後処理 (計算結果の整理)
タンクおよび移動する正方形の簡単なアニメーションを以下に示します。
![anim_delay_0.1s.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/f7171195-152a-66b4-431c-95e11675ac38.gif)

次の計算結果の物理量を整理します。

A) 抗力係数$C_D$ (圧力成分(pressure component)、せん断応力成分(viscous component))
B) 渦度場 ($t=5,8[s]$)
C) 圧力場 $p(x,y)-p_{average}$ ($t=5,8[s]$)
D) 流速場 ($t=5,8[s]$)

抗力係数$C_D$の定義は次の通りです。

```math
C_D=\frac{\lvert F_x \rvert}{\frac{1}{2}\rho U_{max}^2 L b}
```
ここで、上式の変数を下表に示します。
<table>
  <caption>表4 抗力係数$C_D$の算出に使用する変数</caption>
  <thead>
    <tr>
      <th> 変数 </th> <th> 説明 </th> <th> 単位 </th> <th> 値 </th>
    </tr>
  </thead>
  <tr>
    <td> $F_x$ </td> <td> 抗力の$x$成分 </td> <td> $[N]$ </td> <td> forces(function object)の出力 </td>
  </tr>
  <tr>
    <td> $\rho$ </td> <td> 流体の密度 </td> <td> $[kg/m^3]$ </td> <td> 1 </td>
  </tr>
  <tr>
    <td> $U_{max}^2$ </td> <td> 正方形の最大流速 </td> <td> $[m/s]$ </td> <td> 1 </td>
  </tr>
  <tr>
    <td> $L$ </td> <td> 正方形の側面の長さ </td> <td> $[m]$ </td> <td> 1 </td>
  </tr>
  <tr>
    <td> $b$ </td> <td> 正方形の横向きの長さ </td> <td> $[m]$ </td> <td> 1 </td>
  </tr>
</table>

圧力場は当初、SPHERIC Test06 と同じく $p(x,y)-p(0,0)$ で表示しようとしましたが、それぞれの解析メッシュの $p(0,0)$ のバラツキが大きく、図化すると差異が大きくなりました。
そこで、$p(0,0)$ ではなく圧力場の平均値 $p_{average}$ を用いて、$p(x,y)-p_{average}$ で図化したところ、差異が小さくなったため、こちらを採用しています。
なお、圧力場の平均値 $p_{average}$ は paraview の plotDataOverTime で取得しています。

## A) 抗力係数
![fig02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/cdf06080-23aa-cff0-5b95-81d9d6422049.png "図2 抗力係数")
図2 抗力係数 (左が圧力成分、右がせん断応力成分)

## B) 渦度場
![fig03_coarse.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/350cc640-a8e2-db3a-99d3-eb7b74d304af.png "図3 渦度場(粗いメッシュ)")
![fig03_standard.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/57f0cdd2-59cd-6dd8-6bc7-c3347692c163.png "図3 渦度場(標準メッシュ)")
![fig03_fine.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/b926a02b-8e28-9053-0b50-05d6be5969bd.png "図3 渦度場(詳細メッシュ)")
図3 渦度場

## C) 圧力場
![fig04_coarse.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/235ea109-b7d9-0208-5bd1-7432c6bcae9b.png "図4 圧力場(粗いメッシュ)")
![fig04_standard.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/21507e60-12e9-cc37-8fe7-cd61181f0c50.png "図4 圧力場(標準メッシュ)")
![fig04_fine.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/050c4a55-8e40-2d1a-8528-3fc1c30ab16f.png "図4 圧力場(詳細メッシュ)")
図4 圧力場 $p(x,y)-p_{average}$

## D) 流速場
![fig05_coarse.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/c19e6a88-a810-1b41-ba37-2f137461c679.png "図5 流速場(粗いメッシュ)")
![fig05_standard.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/039d590e-ef7d-2584-0250-fa1a8753186b.png "図5 流速場(標準メッシュ)")
![fig05_fine.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1784513/0931a57b-b799-f812-6a89-49043a72f68d.png "図5 流速場(詳細メッシュ)")
図5 流速場

# まとめ
今回の解析結果は、SPHERIC Test06 の結果とほぼ一致していると思われます。しかしながら、抗力係数の圧力成分が振動しており、圧力の計算の振動をより抑制できれば、より良い計算結果が得られると思われます。