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

はじめに

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

解析モデル

図1 解析モデル

表1 解析モデル

名前	記号	設定値
タンク (長方形箱)	$L_x$ [m]	10
	$L_y$ [m]	5
移動する正方形	一辺の長さ $L$ [m]	1
	最大速度 $U_{max}$ [m/s]	1

解析条件

移動する正方形のレイノルズ数$Re$がそれぞれ 50、100 および 150 のケースを計算します。
$Re$ の設定は、$Re$ の定義式 $Re=U_{max}L/ \nu$ から $\nu$ を介して行います。

表2 解析条件1

$Re$	$\nu [m^2/s]$
50	0.02
100	0.01
150	6.67E-03

解析メッシュは、粗い、標準および詳細の3つを用意します。

表3 解析条件2

解析メッシュ	タンク (長方形箱)		移動する正方形
分割数	$n_x$	$n_y$	$n$
粗い (coarse)	267	133	53
標準 (standard)	400	200	80
詳細 (fine)	600	300	120

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

前処理 (解析準備)

解析は、次のチュートリアルを流用して行います。
$FOAM_TUTORIALS/incompressible/overPimpleDyMFoam/simpleRotor

system/blockMesh

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

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 ブロックを追加しています。

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 に変更しています。

system/decomposeParDict

numberOfSubdomains  3;

method          scotch;
//method          hierarchical;

coeffs
{
    n           (3 1 1);
}

system/fvSchemes

oversetInterpolation ブロックの「method inverseDistance;」以外はコメントにします。
コメントにしないと正常に計算ができません。

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に増やす

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

変更する必要はありません。

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 を変更します。

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

変更する必要はありません。

constant/turbulenceProperties

simulationType laminar;

constant/transportProperties

「transportModel Newtonian;」は変更しません。nu の値は表2のように変更します。
この2つ以外はコメントにします。

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 に変更します。

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

変更する必要はありません。

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

変更する必要はありません。

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

変更する必要はありません。

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

計算の実行で使用するシェルスクリプトです。
変更する必要はありません。

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 から呼び出されるシェルスクリプトです。
変更する必要はありません。

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

計算結果を削除するときに使用するシェルスクリプトです。
変更する必要はありません。

Allclean

#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
#------------------------------------------------------------------------------

cleanCase0

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

後処理 (計算結果の整理)

タンクおよび移動する正方形の簡単なアニメーションを以下に示します。

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

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$の定義は次の通りです。

C_D=\frac{\lvert F_x \rvert}{\frac{1}{2}\rho U_{max}^2 L b}

ここで、上式の変数を下表に示します。

表4 抗力係数$C_D$の算出に使用する変数

変数	説明	単位	値
$F_x$	抗力の$x$成分	$[N]$	forces(function object)の出力
$\rho$	流体の密度	$[kg/m^3]$	1
$U_{max}^2$	正方形の最大流速	$[m/s]$	1
$L$	正方形の側面の長さ	$[m]$	1
$b$	正方形の横向きの長さ	$[m]$	1

圧力場は当初、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) 抗力係数

図2 抗力係数 (左が圧力成分、右がせん断応力成分)

B) 渦度場

図3 渦度場

C) 圧力場

図4 圧力場 $p(x,y)-p_{average}$

D) 流速場

図5 流速場

まとめ

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