角柱の抗力係数と揚力係数について～実験値と異なる値を示す原因～

[arima]

ユーザーの皆様へ

icoFsiFoam(DEXCS_FSI 3D) を使って流速20m/sで片持ち柱（高さH=0.6m 幅B=0.05m 奥行D=0.05m）、
動粘性係数0.001m2/s、Re=1000を　//# include "setPressure.H"、//# include
"solveSolid.H"　を読み込まずに、つまり連成させずにicoFoamと同じ条件下で、解析領域は例題の２倍の高さにした直方体(高さ
2m×奥行1m×h幅6m)で解析を行っています。片持ち柱の設置場所は解析領域幅の端から2mの位置です。

本解析が正しく表しているかどうか、他の論文と抗力係数と揚力係数を比較すると異なった値を算出しており、原因が分かりません。多くの意見を参考にして
原因を解決したいのでよろしくお願いします。

いくつかの論文を参考にするとRe=1000での実験値と解析値はともに抗力係数、揚力係数はそれぞれCd=2.0,Cl=1.5です。

本解析（icoFsiFoamの連成なし）の平均抗力係数と平均揚力係数はCd=1.24 Cl=0.06

圧力の収束により異なった結果が表れるという意見もありましたのでfvsolutionと収束状況も合わせて以下に記します。

functionは以下のように設定しました。
forces
{
type forces;
functionObjectLibs ("libforces.so"); //Lib to load
patches (consoleFluid); // change to your patch name
rhoInf 1.225; //Reference density for fluid
CofR (0.15 0 0); //Origin for moment calculations
}

forceCoeffs
{
type forceCoeffs;
functionObjectLibs ("libforces.so");
patches (consoleFluid); //change to your patch name
rhoInf 1.225;
CofR (0 0 0); //Origin of the moment calculation
liftDir (0 1 0); // direction on Cl
dragDir (1 0 0); // direction on Cd
pitchAxis (0 1 0);
magUInf 20;
lRef 0.6; //length of the object i.e. h(0.6)
Aref 0.03; //area of the object B(0.05)*h(0.6)
}


fvsolutionは以下です。

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

solvers
{
p ICCG 1e-06 0;
// p AMG 1e-06 0 500;
U BICCG 1e-05 0;
}

PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}


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

収束はこのようになっています。

Courant Number mean: 0.043256 max: 0.71997 velocity magnitude: 25.7606
Time = 0.0402

PCG: Solving for motionUx, Initial residual = 0, Final residual = 0,
No Iterations 0
PCG: Solving for motionUy, Initial residual = 0, Final residual = 0,
No Iterations 0
PCG: Solving for motionUz, Initial residual = 0, Final residual = 0,
No Iterations 0
volume continuity errors : volume = 11.9985, max error = 0, sum local
= 0, global = 0
Motion magnitude: mean = 0 max = 0
PBiCG: Solving for Ux, Initial residual = 0.00292805, Final residual
= 4.68933e-07, No Iterations 2
PBiCG: Solving for Uy, Initial residual = 0.00565219, Final residual
= 6.20721e-07, No Iterations 2
PBiCG: Solving for Uz, Initial residual = 0.00548541, Final residual
= 3.3651e-07, No Iterations 2
PCG: Solving for p, Initial residual = 0.0150603, Final residual =
9.65444e-07, No Iterations 143
PCG: Solving for p, Initial residual = 9.64287e-07, Final residual =
9.64287e-07, No Iterations 0
Moving mesh time step continuity errors : sum local = 6.30886e-11,
global = -9.04344e-13, cumulative = 8.80474e-10
PBiCG: Solving for Ux, Initial residual = 5.57864e-05, Final residual
= 1.00794e-06, No Iterations 1
PBiCG: Solving for Uy, Initial residual = 0.000178689, Final residual
= 3.81521e-06, No Iterations 1
PBiCG: Solving for Uz, Initial residual = 0.000238688, Final residual
= 8.35888e-06, No Iterations 1
PCG: Solving for p, Initial residual = 0.0101316, Final residual =
9.98591e-07, No Iterations 222
PCG: Solving for p, Initial residual = 9.99119e-07, Final residual =
9.99119e-07, No Iterations 0
Moving mesh time step continuity errors : sum local = 6.53907e-11,
global = -5.83174e-13, cumulative = 8.79891e-10
ExecutionTime = 3618.06 s


Courant Number mean: 0.0432599 max: 0.719758 velocity magnitude:
25.7521
Time = 0.0405

PCG: Solving for motionUx, Initial residual = 0, Final residual = 0,
No Iterations 0
PCG: Solving for motionUy, Initial residual = 0, Final residual = 0,
No Iterations 0
PCG: Solving for motionUz, Initial residual = 0, Final residual = 0,
No Iterations 0
volume continuity errors : volume = 11.9985, max error = 0, sum local
= 0, global = 0
Motion magnitude: mean = 0 max = 0
PBiCG: Solving for Ux, Initial residual = 0.00291522, Final residual
= 4.58361e-07, No Iterations 2
PBiCG: Solving for Uy, Initial residual = 0.00564626, Final residual
= 6.18169e-07, No Iterations 2
PBiCG: Solving for Uz, Initial residual = 0.00549071, Final residual
= 3.65492e-07, No Iterations 2
PCG: Solving for p, Initial residual = 0.0152259, Final residual =
9.66449e-07, No Iterations 143
PCG: Solving for p, Initial residual = 9.65238e-07, Final residual =
9.65238e-07, No Iterations 0
Moving mesh time step continuity errors : sum local = 6.32453e-11,
global = -1.17355e-12, cumulative = 8.78717e-10
PBiCG: Solving for Ux, Initial residual = 5.71075e-05, Final residual
= 1.01568e-06, No Iterations 1
PBiCG: Solving for Uy, Initial residual = 0.000179952, Final residual
= 3.80627e-06, No Iterations 1
PBiCG: Solving for Uz, Initial residual = 0.000238904, Final residual
= 8.26208e-06, No Iterations 1
PCG: Solving for p, Initial residual = 0.0103706, Final residual =
9.87111e-07, No Iterations 221
PCG: Solving for p, Initial residual = 9.87539e-07, Final residual =
9.87539e-07, No Iterations 0
Moving mesh time step continuity errors : sum local = 6.47464e-11,
global = -7.73899e-13, cumulative = 8.77943e-10
ExecutionTime = 3643.28 s

[ohbuchi]

DEXCS-FSIを使ったことが無いので詳しいことは解りませんが、コメントさせて頂きます。
文献が何か解らないのですが、恐らく文献は２次元角柱周りの値だと思います。
計算しているのは、角柱が領域全体に渡っておらず端面をもつ突起周りの流れなので
問題が異なると思います。突起周りの流れでは、二次元角柱よりもCdは小さくなるのでは
ないでしょうか？

また、Clはカルマン渦を放出することによる変動揚力で、流れに直交する成分の力の
振幅を無次元化した値で、froceCoeffで計算される値とは違うと思います。
forceの時間変動波形を見て、自分で計算して下さい。

以上、ご参考まで。

[arima]

Ohbuchi様

ご回答ありがとうございます。
2次元正方形(1m*1m)角柱周りも抗力係数をチェックしますとRe=1000~2000,3000...10000でCd=1.3~1.4でどのレ
イノルズ数でも変化なく、実験値Cd=1.8~2.0(Re=1000~10000)と異なって計算されてしまいます。forcesから
pressureの値と動圧＊面積を除して電卓で計算すると、たしかに抗力係数は1.3~1.4の間の値でした。
全く原因が見つからないので、ご意見よろしくお願いします。

2次元では、pimpleFoamで、　流入の圧力をzerogradient、流速を一様流, 流出の圧力をfixedvalue 0、流速を
zerogradient
乱流モデル　RNGkEpsilon、 移流スキーム　limitedLinear 1です。

> > ExecutionTime = 3643.28 s- 引用テキストを表示しない -
>
> - 引用テキストを表示 -

[ohbuchi]

試しに、大きさ１ｍの角柱周りの二次元流れをRe=1000で計算してみました。
標準ｋ－ε、simpleFoamで計算したところ、非常に粗いメッシュ（70ｘ３０）で
Cd=2.1となりました。

メッシュとforceCoeffsの設定が合っていないのが原因ではないでしょうか？

ちなみに私のテストしたケースは以下の通りです。


******* transportProperties **************************
transportModel Newtonian;
nu nu [0 2 -1 0 0 0 0] 1e-03;

******* RASProperties ********************************
RASModel kEpsilon;
turbulence on;

******** blockMeshDict ********************************
convertToMeters 1.0;

vertices
(
(-3 -3 -0.5) // 0
(-0.5 -3 -0.5) // 1
( 0.5 -3 -0.5) // 2
(10 -3 -0.5) // 3
(-3 -0.5 -0.5) // 4
(-0.5 -0.5 -0.5) // 5
( 0.5 -0.5 -0.5) // 6
(10 -0.5 -0.5) // 7
(-3 0.5 -0.5) // 8
(-0.5 0.5 -0.5) // 9
( 0.5 0.5 -0.5) // 10
(10 0.5 -0.5) // 11
(-3 3 -0.5) // 12
(-0.5 3 -0.5) // 13
( 0.5 3 -0.5) // 14
(10 3 -0.5) // 15

(-3 -3 0.5) // 16
(-0.5 -3 0.5) // 17
( 0.5 -3 0.5) // 18
(10 -3 0.5) // 19
(-3 -0.5 0.5) // 20
(-0.5 -0.5 0.5) // 21
( 0.5 -0.5 0.5) // 22
(10 -0.5 0.5) // 23
(-3 0.5 0.5) // 24
(-0.5 0.5 0.5) // 25
( 0.5 0.5 0.5) // 26
(10 0.5 0.5) // 27
(-3 3 0.5) // 28
(-0.5 3 0.5) // 29
( 0.5 3 0.5) // 30
(10 3 0.5) // 31
);

blocks
(
hex ( 0 1 5 4 16 17 21 20) (15 10 1) simpleGrading (0.2 0.2 1)
hex ( 1 2 6 5 17 18 22 21) (10 10 1) simpleGrading (1.0 0.2 1)
hex ( 2 3 7 6 18 19 23 22) (45 10 1) simpleGrading (5.0 0.2 1)
hex ( 4 5 9 8 20 21 25 24) (15 10 1) simpleGrading (0.2 1.0 1)
hex ( 6 7 11 10 22 23 27 26) (45 10 1) simpleGrading (5.0 1.0 1)
hex ( 8 9 13 12 24 25 29 28) (15 10 1) simpleGrading (0.2 5.0 1)
hex ( 9 10 14 13 25 26 30 29) (10 10 1) simpleGrading (1.0 5.0 1)
hex (10 11 15 14 26 27 31 30) (45 10 1) simpleGrading (5.0 5.0 1)
);

edges
(
);

patches
(
patch inlet
(
(0 16 20 4)
(4 20 24 8)
(8 24 28 12)
)
patch outlet
(
( 3 7 23 19)
( 7 11 27 23)
(11 15 31 27)
)
wall blockWall
(
(5 21 22 6)
(9 10 26 25)
(5 9 25 21)
(6 22 26 10)
)
wall walls
(
(0 1 17 16)
(1 2 18 17)
(2 3 19 18)
(12 28 29 13)
(13 29 30 14)
(14 30 31 15)
)
empty back
(
(0 4 5 1)
(1 5 6 2)
(2 6 7 3)
(4 8 9 5)
(6 10 11 7)
(8 12 13 9)
(9 13 14 10)
(10 14 15 11)
)
empty front
(
(16 17 21 20)
(17 18 22 21)
(18 19 23 22)
(20 21 25 24)
(22 23 27 26)
(24 25 29 28)
(25 26 30 29)
(26 27 31 30)
)
);

mergePatchPairs
(
);

********* boundary **********************
6
(
inlet
{
type patch;
nFaces 30;
startFace 3880;
}
outlet
{
type patch;
nFaces 30;
startFace 3910;
}
blockWall
{
type wall;
nFaces 40;
startFace 3940;
}
walls
{
type wall;
nFaces 140;
startFace 3980;
}
back
{
type empty;
nFaces 2000;
startFace 4120;
}
front
{
type empty;
nFaces 2000;
startFace 6120;
}
)

********* 0/U ************************************
dimensions [0 1 -1 0 0 0 0];

internalField uniform (1 0 0);

boundaryField
{
inlet
{
type fixedValue;
phi phi;
value uniform (1.0 0.0 0.0);
}

outlet
{
type zeroGradient;
}

blockWall
{
type fixedValue;
value uniform (0.0 0.0 0.0);
}

walls
{
type fixedValue;
value uniform (0.0 0.0 0.0);
}

front
{
type empty;
}

back
{
type empty;
}
}

******** 0/p *************************************
dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{
inlet
{
type zeroGradient;
}

outlet
{
type fixedValue;
value uniform 0;
}

walls
{
type zeroGradient;
}

blockWall
{
type zeroGradient;
}

front
{
type empty;
}

back
{
type empty;
}
}

******** controlDict ****************************
application simpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 2000;
deltaT 1;
writeControl timeStep;
writeInterval 500;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;

functions
{

forces
{
type forces;
functionObjectLibs ( "libforces.so" );

outputControl timeStep;
outputInterval 500;
patches (blockWall);
pName p;
UName U;
rhoName rhoInf;
log true;
rhoInf 1;
CofR (0.0 0.0 0.0);

}
forceCoeffs
{
type forceCoeffs;
functionObjectLibs ("libforces.so");

patches (blockWall);
rhoName rhoInf;
rhoInf 1.0;
CofR (0 0 0);
liftDir (0 1 0);
dragDir (1 0 0);
pitchAxis (0 0 1);
magUInf 1;
lRef 1;
Aref 1;
outputControl timeStep;
outputInterval 500;
}
}


***** forceCoeffs/0/forceCoeffs.dat ********************************
# Time Cd Cl Cm
500 2.11932 -0.0820795 0.00985084
1000 2.12105 0.0816687 -0.0104678
1500 2.11935 -0.0569716 0.00799756
2000 2.11851 0.0123197 -0.00272281

[arima]

ohbuchi様

ケースの公開本当にありがとうございます。
貴氏のメッシュで私の設定条件(pimpleFoam)で動かしてみると、Re=1000で実験値と同等の抗力係数Cd≒2.2がでました。つまり、単
純に私のメッシュが適切でなかったようです。さらに細かくメッシュを切ると他のRe数でも実験値と同様の結果を得ることができました。私のメッシュは、
物体を中心に放射状にメッシュを切っていましたので、物体境界面のメッシュはどうしても台形もしくは平行四辺形で切れており、メッシュが歪んだため？正
確な値がでなかったと考えられます。
私の素人知見ですが、適切な抗力（揚力）係数を出すためには
①対象物体の周りはなるべくメッシュを細かく切ること
②正方形もしくは長方形のように直角に切ってあげること
だと思います。
しかしながら、流速分布や圧力分布は私のメッシュの切り方でも期待通りの結果が出ており、少々歪んでも大丈夫なようです。ただ、抗力係数のように物体周
りの圧力がシビアに要求される場合は、適切にメッシュを作成する必要があるようです。抗力係数がこれほどメッシュに依存するとは、驚きとともに非常に勉
強になりました。本当にありがとうございます。
これでなんとか卒業研究が進められます。＾＾；

> ...
>
> もっと読む ≫