はじめに
ESI版 (OpenCFD版) OpenFOAM に搭載されているオーバーセットメッシュの機能を使って、「SPHERIC Test06 2-D Incompressible flow around a moving square inside a rectangular box」の再現解析を行います。OpenFOAM のバージョンは、v2306 とします。
解析方法などは、オープンCAE学会と関連が深い「改訂新版 OpenFOAM の歩き方 川畑真一著 (株式会社インプレスR&D)」を参考にしています。
解析モデル
図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$ を介して行います。
| $Re$ | $\nu [m^2/s]$ |
|---|---|
| 50 | 0.02 |
| 100 | 0.01 |
| 150 | 6.67E-03 |
解析メッシュは、粗い、標準および詳細の3つを用意します。
| 解析メッシュ | タンク (長方形箱) | 移動する正方形 | |
|---|---|---|---|
| 分割数 | $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に合わせてそれぞれ変更します。
こちらも、上の行がタンクに、下の行が移動する正方形に該当します。
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 ブロックを追加しています。
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 に変更しています。
numberOfSubdomains 3;
method scotch;
//method hierarchical;
coeffs
{
n (3 1 1);
}
system/fvSchemes
oversetInterpolation ブロックの「method inverseDistance;」以外はコメントにします。
コメントにしないと正常に計算ができません。
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に増やす
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
変更する必要はありません。
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 を変更します。
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
変更する必要はありません。
simulationType laminar;
constant/transportProperties
「transportModel Newtonian;」は変更しません。nu の値は表2のように変更します。
この2つ以外はコメントにします。
//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 に変更します。
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
変更する必要はありません。
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
変更する必要はありません。
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
変更する必要はありません。
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
計算の実行で使用するシェルスクリプトです。
変更する必要はありません。
#!/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 から呼び出されるシェルスクリプトです。
変更する必要はありません。
#!/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
計算結果を削除するときに使用するシェルスクリプトです。
変更する必要はありません。
#!/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}
ここで、上式の変数を下表に示します。
| 変数 | 説明 | 単位 | 値 |
|---|---|---|---|
| $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 の結果とほぼ一致していると思われます。しかしながら、抗力係数の圧力成分が振動しており、圧力の計算の振動をより抑制できれば、より良い計算結果が得られると思われます。