自動車周りの解析について
suzuswift
suzuswiftと申します。いつも皆様の書き込みを参考にさせていただいております。
大学の研究でOpenFOAMを利用しております。
さて、この度どうしても自力では解決できない問題に直面したのでお知恵を貸していただければと思い書き込み致しました。
少々長くなりますがご了承下さい。
現在OpenFOAMによって自動車周りの空気の流れを解析しようとしております。
ソルバはsimpleFoam、乱流モデルはkOmegaSSTを使おうとしています。
そのため、チュートリアルのmotorbikeおよびこちらでのこれまでの書き込みを参考に境界条件を設定しました。
しかし、計算が途中で止まって(発散して)しまいうまくいきません。
どこに問題があるのか自力ではわからず、こちらの方々にチェックしていただきたく思います。
現状を整理致します。
1.メッシュは某商用メッシャにてテトラメッシュを作成。checkMeshの結果は下に載せます。最小メッシュ幅4mm、クーラン数を考慮すると時間
刻み幅は最大で0.0002s。商用メッシャとはいえ、この時点でおかしい可能性もあるかもしれません。なお、空間サイズは幅7m、高さ5m、長さ
40m。自動車サイズは幅2m、高さ1.4m、長さ4m。自動車の位置は空間内の2.5~4.5m、0~1.4m、9~13mの間に配置しておりま
す。
2.0ディレクトリ、constantディレクトリ、systemディレクトリはmotorbikeを参考にしつつ作成。こちらもそれぞれ下に載せま
す。なお、壁面は流入口、流出口、地面、空、側面、車体で考えており、地面と車体は固定値(wall)です。(本来なら地面も動かすべきなのでしょうが
今回は車体と同じ条件だとどうなるかやってみたい) 空と側面はslip壁です。
3.system内のcontrolDictに関してですが、空気抵抗係数に関するforceCoeffsのみfunctionにて追加、出力したい。
しかし、途中まで進んだ計算結果を見ると明らかにおかしい。これも下に載せます。
以上について、どこから手をつけていいのか自分ではわからないのでご助言よろしくお願い致します。
以下、各種ファイル、および出力です。
Create time
Create polyMesh for time = 0
Time = 0
Mesh stats
points: 209140
faces: 2102593
internal faces: 1946659
cells: 1012313
boundary patches: 6
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
hexahedra: 0
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 1012313
polyhedra: 0
Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).
Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface
topology
body 127066 63535 ok (closed singly
connected)
ground 9196 4758 ok (non-closed singly
connected)
inlet 858 471 ok (non-closed singly
connected)
outlet 1006 546 ok (non-closed singly
connected)
side 10308 5458 ok (non-closed singly
connected)
sky 7500 3909 ok (non-closed singly
connected)
Checking geometry...
Overall domain bounding box (0 0 -2.220446e-15) (7.0866 4.953
40.131)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (6.8826e-20 1.081551e-18 -1.628616e-19) OK.
Max cell openness = 1.868291e-16 OK.
Max aspect ratio = 606.5899 OK.
Minumum face area = 1.727455e-05. Maximum face area = 0.2641945.
Face area magnitudes OK.
Min volume = 5.200124e-08. Max volume = 0.05002422. Total volume
= 1402.506. Cell volumes OK.
Mesh non-orthogonality Max: 89.59495 average: 21.2547
*Number of severely non-orthogonal faces: 281.
Non-orthogonality check OK.
<<Writing 281 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 1.664459 OK.
Mesh OK.
End
/////////////////////////////////////////////////////////////////
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
dimensions [ 0 2 -2 0 0 0 0 ];
internalField uniform 1.483;
boundaryField
{
inlet
{
type fixedValue;
value uniform 1.483;
}
outlet
{
type inletOutlet;
inletValue uniform 1.483;
value uniform 1.483;
}
body
{
type kqRWallFunction;
}
sky
{
type slip;
}
side
{
type slip;
}
ground
{
type kqRWallFunction;
}
}
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
dimensions [0 2 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
body
{
type nutkWallFunction;
value uniform 0;
}
ground
{
type nutkWallFunction;
value uniform 0;
}
inlet
{
type calculated;
value uniform 0;
}
outlet
{
type calculated;
value uniform 0;
}
side
{
type calculated;
value uniform 0;
}
sky
{
type calculated;
value uniform 0;
}
}
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
dimensions [ 0 0 -1 0 0 0 0 ];
internalField uniform 4.492;
boundaryField
{
inlet
{
type fixedValue;
value uniform 4.492;
}
outlet
{
type inletOutlet;
inletValue uniform 4.492;
value uniform 4.492;
}
body
{
type omegaWallFunction;
}
sky
{
type slip;
}
side
{
type slip;
}
ground
{
type omegaWallFunction;
}
}
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
dimensions [ 0 2 -2 0 0 0 0 ];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
body
{
type zeroGradient;
}
sky
{
type slip;
}
side
{
type slip;
}
ground
{
type zeroGradient;
}
}
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
dimensions [ 0 1 -1 0 0 0 0 ];
internalField uniform (0 0 20);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 20);
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (0 0 20);
}
body
{
type fixedValue;
value uniform (0 0 0);
}
sky
{
type slip;
}
side
{
type slip;
}
ground
{
type fixedValue;
value uniform (0 0 0);
}
}
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object RASProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
RASModel kOmegaSST;
turbulence on;
printCoeffs on;
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
transportModel Newtonian;
nu nu [0 2 -1 0 0 0 0] 1.512e-05;
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
application simpleFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 100;
deltaT 0.00001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
adjustTimeStep yes;
maxCo 0.5;
maxDeltaT 0.1;
libs
(
"libOpenFOAM.so"
"libincompressibleRASModels.so"
"libincompressibleTurbulenceModel.so"
);
functions
(
forceCoeffs
{
type forceCoeffs;
functionObjectLibs ( "libforces.so" );
outputControl timeStep;
outputInterval 1;
patches (body);
pName p;
UName U;
rhoName rhoInf;
log true;
rhoInf 1.205;
CofR (3.5 0.66 11 );
liftDir ( 0 1 0 );
dragDir ( 0 0 1 );
pitchAxis ( 1 0 0 );
magUInf 20.0;
lRef 3.99;
Aref 2.108;
}
);
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwind grad(U);
div(phi,k) Gauss upwind;
div(phi,omega) Gauss upwind;
div((nuEff*dev(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear limited 0.5;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
fluxRequired
{
default no;
p;
}
//
************************************************************************* //
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-7;
relTol 0.1;
smoother GaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
cacheAgglomeration on;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 10;
mergeLevels 1;
};
U
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-8;
relTol 0.1;
nSweeps 1;
};
k
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-8;
relTol 0.1;
nSweeps 1;
};
omega
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-8;
relTol 0.1;
nSweeps 1;
};
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
residualControl
{
p 1e-2;
U 1e-2;
"(k|omega)" 1e-2;
}
}
potentialFlow
{
nNonOrthogonalCorrectors 10;
}
relaxationFactors
{
p 0.3;
U 0.7;
k 0.7;
omega 0.7;
}
cache
{
grad(U);
}
//
************************************************************************* //
発散直前の出力です。
Time = 0.00102
smoothSolver: Solving for Ux, Initial residual = 0.457847, Final
residual = 0.364793, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.438349, Final
residual = 0.330726, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.460352, Final
residual = 0.318314, No Iterations 2
GAMG: Solving for p, Initial residual = 1.3449e-07, Final residual =
5.79091e-8, No Iterations 1
time step continuity errors : sum local = 2.07573e+17, global =
2.28925e+13, cuulative = 6.56304e+13
smoothSolver: Solving for omega, Initial residual = 0.270269, Final
residual =0.0245962, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.309657, Final
residual = 0.072769, No Iterations 3
ExecutionTime = 74.83 s ClockTime = 122 s
forceCoeffs output:
Cd = -5.77385e+47
Cl = 9.10268e+48
Cm = 6.42402e+48
ohbuchi
運動方程式の残差の減少が少ない様です。通常、Initial ResidualとFinal Residualは相当に大きさが違う筈です。殆ど収束
していない様です。
圧力方程式の残差は小さいですが、連続の式の残差はとても大きく
このステップが正しく計算できていないことが解ります。
Uの設定ファイルを見ると、出口でInletOutletで値を指定されている
様ですが、通常はZeroGradientではないでしょうか？
出入り口で速度を規定してしまうと、連続条件を満たすことができなくなります。
恐らく、これを修正することで計算できると思います。お試しください。
k，Omegaも出口でZeroGradientにするのが無難かと思います。
尚、simpleFoamは定常計算ですのでクーラン数や時間刻みは関係ありません。
安定性を改善するにはfvSolutionのrelaxationFactorsを小さく（特にｐ）すると
良いと思います。
以上、ご参考まで。
Masashi Imano
流出口の境界条件については、ohbuchiさんが書かれている通りzeroGradientに
するほうが無難です。
また、最初は移流項の離散化スキーム(divSchemes)は、pitzDailyのfvScheme
のようにU,k,omegaに対して全てupwindにするほうが無難です。
他の条件が妥当で安定に計算できることが確認できてから、linearUpwindや
limitedLinearなどのより高次のスキームに変更したら良いと思います。
U,k,omegaの線型ソルバーについても、pitzDailyのfvSolutionのようにPBiCGに
するほうが通常速いと思います。ただし、pのソルバーはmoterBikeで用いている
GAMGのほうが並列数が大きくない限り一般に速いと思います。
suzuswift
ご助言ありがとうございます。
suzuswiftです。
まずohbuchi様の書き込みについてですが、
>圧力方程式の残差は小さいですが、連続の式の残差はとても大きく
>このステップが正しく計算できていないことが解ります。
発散直前の計算ステップを抜き出したため、やはり正しく計算できていないようですね。
そのためご指摘の通りにOutletに対して各ファイル(U,k,omega)ともzeroGradientに変更致しました。今現在、順調に計算回数
を増やせております。
>尚、simpleFoamは定常計算ですのでクーラン数や時間刻みは関係ありません。
>安定性を改善するにはfvSolutionのrelaxationFactorsを小さく（特にｐ）すると
>良いと思います。
この件は知りませんでした。てっきり時間刻みを十分小さくすれば発散しない場合もあると思っておりました。simpleFoamの場合は関係ないのです
ね。(もちろん他のソルバの場合は関係あるのでしょうが)
ちなみに、ここの「安定性を改善する」とは発散しづらくなると認識してよろしいのでしょうか?
試しに、pを0.2,それ以外を0.6にして現在計算しております。
続いて、今野様の書き込みについてですが、
> また、最初は移流項の離散化スキーム(divSchemes)は、pitzDailyのfvScheme
> のようにU,k,omegaに対して全てupwindにするほうが無難です。
チュートリアルのうち、自分がやりたいことと似たようなことをやっているmotorBikeしか見ておりませんでした。もう少し幅広い視野を持たないと
いけませんでした…
ご指摘の通り、fvScheme内のdivSchemesのU,k,omegaをupwindに切り替えました。
> U,k,omegaの線型ソルバーについても、pitzDailyのfvSolutionのようにPBiCGに
> するほうが通常速いと思います。ただし、pのソルバーはmoterBikeで用いている
> GAMGのほうが並列数が大きくない限り一般に速いと思います。
こちらもご指摘の通り、fvSolutionのU,k,omegaをPBiCGに切り替えました。pは元々GAMGだったためいじっておりません。
ちなみに、書き忘れておりましたが8コアで並列計算しております。8コアは並列数としてそれほど大きくないと思うのですがいかがでしょうか?
P.S.
現在、計算途中の結果を出力してみたところ、以下の通りになりました。
Time = 0.01079
DILUPBiCG: Solving for Ux, Initial residual = 0.00176682, Final
residual = 3.78507e-05, No Iterations 2
DILUPBiCG: Solving for Uy, Initial residual = 0.000800157, Final
residual = 1.83885e-05, No Iterations 2
DILUPBiCG: Solving for Uz, Initial residual = 0.00134856, Final
residual = 9.78541e-05, No Iterations 2
GAMG: Solving for p, Initial residual = 0.0192136, Final residual =
0.00179004, No Iterations 1
time step continuity errors : sum local = 0.00206689, global =
1.99371e-05, cumulative = 0.835609
DILUPBiCG: Solving for omega, Initial residual = 5.80788e-05, Final
residual = 9.4454e-06, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.00208878, Final
residual = 3.41734e-05, No Iterations 2
ExecutionTime = 2121.61 s ClockTime = 2741 s
forceCoeffs output:
Cd = 4.8226e+09
Cl = 2.68272e+07
Cm = -1.31267e+09
////////////////////////////////
forceCoeffsがとんでもない値を出しているように見えます…
計算回数が増えれば徐々に落ち着くものなのでしょうか?
kaji
横から失礼します。
気づいたことを簡単に述べますと、
1. 抵抗係数がおかしいとのことですが、U_inf=20 [m/s]でT=0.01079 [s]でなるとまだまだ非定常の流場であると思われます
ので、
40 [m]風洞の長さを流体が通りきる20 [s]ぐらいは様子を見た方がいいと思います。しかしこれも目安でしかないので
実際は抵抗値の推移などを見て余裕をもって定常になったと判断するべきだと思います。
つまり今はまだ非定常であり、抵抗係数が変な値でもおかしくないということです。
2. timestepが小さすぎるように感じます。これでは時間がかかると思うので0.01程度にして様子を見てはいかがでしょうか。
3. これは自信がありませんが、checkMeshの結果を見るに100万セルの計算格子であると思われますので8CPUによる並列計算ではあまり計
算速度は変わらないと思います。
4. 抵抗係数や揚力係数がforceCoeffsでうまく出力されない場合、forcesでの抵抗、揚力をもとに手計算をしてみることをお勧めしま
す。
ソースを見てもforceCoeffsにおかしい点は見つけられませんでしたが、手計算でやったらうまくいった、ということがありました。
以上です。
梶川
suzuswift
ご助言ありがとうございます。
suzuswiftです。
> 1. 抵抗係数がおかしいとのことですが、U_inf=20 [m/s]でT=0.01079 [s]でなるとまだまだ非定常の流場であると思われます
> ので、
> 40 [m]風洞の長さを流体が通りきる20 [s]ぐらいは様子を見た方がいいと思います。しかしこれも目安でしかないので
> 実際は抵抗値の推移などを見て余裕をもって定常になったと判断するべきだと思います。
> つまり今はまだ非定常であり、抵抗係数が変な値でもおかしくないということです。
やはりまだまだ計算を重ねないといけないようですね。了解致しました。
ちなみに、細かいことで申し訳ないのですが(おそらくタイプミスかと思います)、
40[m]を20[m/s]で流しているので2[s]で流れきるはずです。
> 2. timestepが小さすぎるように感じます。これでは時間がかかると思うので0.01程度にして様子を見てはいかがでしょうか。
そうですね、上述の通りsimpleFoamはクーラン数に関係ないとのことなので時間刻みをこれまでよりも大きくとって、ご指摘のように0.01でや
ってみます。
> 3. これは自信がありませんが、checkMeshの結果を見るに100万セルの計算格子であると思われますので8CPUによる並列計算ではあまり計
> 算速度は変わらないと思います。
わかりました。あまり大規模な並列計算ができない環境なので残念ながらこのまま進めていきます...
> 4. 抵抗係数や揚力係数がforceCoeffsでうまく出力されない場合、forcesでの抵抗、揚力をもとに手計算をしてみることをお勧めしま
> す。
> ソースを見てもforceCoeffsにおかしい点は見つけられませんでしたが、手計算でやったらうまくいった、ということがありました。
そのようなこともあるのですね。controlDictのfunctionに以下のようにforcesを追加致しました。
forces
{
type forces;
functionObjectLibs ("libforces.so");
patches (body);
pName p;
Uname U;
rhoName rhoInf;
rhoInf 1.205;
CofR (3.5 0.66 11);
outputControl timeStep;
outputInterval 1;
log true;
}
この状態でやってみたいと思います。
ありがとうございました。
suzuswift
追記させていただきます。
計算結果が出たのですが、やはりforceCoeffs,forcesがおかしな値を出しております。
収束判定にresidualControl(p 1e-2,U|k|omega 1e-3を設定)を使った場合の最終ステップでの出力を以下に示しま
す。
Time = 21.67
DILUPBiCG: Solving for Ux, Initial residual = 0.000518417, Final
residual = 9.12038e-06, No Iterations 2
DILUPBiCG: Solving for Uy, Initial residual = 0.000459858, Final
residual = 7.53172e-06, No Iterations 2
DILUPBiCG: Solving for Uz, Initial residual = 0.000678121, Final
residual = 8.77229e-06, No Iterations 2
GAMG: Solving for p, Initial residual = 0.0062044, Final residual =
0.000427456, No Iterations 1
time step continuity errors : sum local = 0.0354208, global =
-0.000543933, cumulative = 839.754
DILUPBiCG: Solving for omega, Initial residual = 1.53419e-05, Final
residual = 2.7426e-06, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.000999455, Final
residual = 1.54947e-05, No Iterations 2
ExecutionTime = 5358.15 s ClockTime = 5573 s
SIMPLE solution converged in 21.67 iterations
forceCoeffs output:
Cd = 7.23396e+07
Cl = -1.60102e+08
Cm = -2.45479e+06
forces output:
forces(pressure, viscous)((1.81941e+10 -8.04977e+10 3.59516e+10)
(5.71084e+08 -8.38369e+08 7.98903e+08))
moment(pressure, viscous)((-5.49548e+09 9.15458e+09 -1.15318e+10)
(5.19548e+08 2.1847e+08 -5.97667e+08))
local:
forces(pressure, viscous)((1.81941e+10 -8.04977e+10 3.59516e+10)
(5.71084e+08 -8.38369e+08 7.98903e+08))
moment(pressure, viscous)((-5.49548e+09 9.15458e+09 -1.15318e+10)
(5.19548e+08 2.1847e+08 -5.97667e+08))
End
Finalising parallel run
e+07のように指数がついていてとんでもない値になっているように見えます。
residualControlを使わずにもっと計算回数を増やすべきなのでしょうか?
それともforcesの設定にミスがあったのでしょうか?
計算エラーの類でないのでどうにも自分ではわかりません、たびたび申し訳ありませんがご助言よろしくお願い致します。
Masashi Imano
ログを見るとそもそも連続の式の誤差が大きぎます。
> time step continuity errors : sum local = 0.0354208, global =
> -0.000543933, cumulative = 839.754
これは、ohbuchiさんが以前指摘されていたように圧力の残差が十分に落ちて
いないことによると思います。
実際、pのFinal residualが0.4e-3と大きく、初期残差の1/10程度となってい
ます。
> GAMG: Solving for p, Initial residual = 0.0062044, Final residual =
> 0.000427456, No Iterations 1
恐らくfvSolutionでpのソルバーの収束判定で相対残差許容値のrelTolが
0.1程度となっていると思いますが、これは0.01程度とより小さくし、
さらに絶対残差許容値の toleranceも少なくとも1e-06程度にしたほうが
良いと思います。
なお、simpleFoamなどの定常解法ソルバーは、そもそも時間項を除いて
解いているので、時間刻みdeltaTには意味がありません。
そこで反復回数がわかりやすいようにdeltaTを1とすることが通例です。
tutorials/incompressible/simpleFoam/以下のチュートリアルケースの
controlDictを参考にしてください。
> residualControlを使わずにもっと計算回数を増やすべきなのでしょうか?
residualControlは十分残差が落ちた時に自動的に停止してくれるので便利ですが、
その許容値を大きく設定すると、自分が本当に予測したい値(Cd値等)が収束する
前に計算を止めてしまうので、危険でもあります。
そこで、residualControlの許容値は小さめにすると共に、自分が予測したい値の
収束状況をグラフを描いて確かめることが非常に重要だと思います。
ここでは、forcesのfunctionで値をファイルに出力していると思いますので、
Cd値等が妥当な値を出力するようになったら、その値をグラフに描きながら、
求める計算精度と許容できる計算負荷を考慮しながら、
妥当なresidualControlの許容値を見極めてください。
以上、ご参考まで。
2011年12月20日18:08 suzuswift <suzukis...@gmail.com>:
> --
> このメールは Google グループのグループ「OpenFOAM」の登録者に送られています。
> このグループに投稿するには、open...@googlegroups.com にメールを送信してください。
> このグループから退会するには、openfoam+u...@googlegroups.com にメールを送信してください。
> 詳細については、http://groups.google.com/group/openfoam?hl=ja からこのグループにアクセスしてください。
>
--
IMANO Masashi, Ph.D.
Assistant Professor
Department of Architecture, Graduate School of Engineering,
The University of Tokyo
7-3-1, Hongo, Bunkyo-ku, Tokyo, Japan, 113-8656
E-mail:im...@arch.t.u-tokyo.ac.jp
Phone:+81-3-5841-6164(direct), +81-3-5841-6179(Laboratory)
Facsimile:+81-3-5841-8511
suzuswift
ご助言ありがとうございます。
suzuswiftです。
> 恐らくfvSolutionでpのソルバーの収束判定で相対残差許容値のrelTolが
> 0.1程度となっていると思いますが、これは0.01程度とより小さくし、
> さらに絶対残差許容値の toleranceも少なくとも1e-06程度にしたほうが
> 良いと思います。
> なお、simpleFoamなどの定常解法ソルバーは、そもそも時間項を除いて
> 解いているので、時間刻みdeltaTには意味がありません。
> そこで反復回数がわかりやすいようにdeltaTを1とすることが通例です。
residualControlの値も小さくしました。
その結果、現在計算途中の結果が以下のようになりました。
Time = 476
DILUPBiCG: Solving for Ux, Initial residual = 0.000232139, Final
residual = 9.55796e-06, No Iterations 2
DILUPBiCG: Solving for Uy, Initial residual = 7.74334e-05, Final
residual = 2.62316e-06, No Iterations 2
DILUPBiCG: Solving for Uz, Initial residual = 6.42992e-05, Final
residual = 4.12506e-06, No Iterations 2
GAMG: Solving for p, Initial residual = 0.00242768, Final residual =
2.23265e-05, No Iterations 3
time step continuity errors : sum local = 7.27245e-06, global =
2.11368e-07, cumulative = 7.19507e-05
DILUPBiCG: Solving for omega, Initial residual = 8.1422e-06, Final
residual = 8.1422e-06, No Iterations 0
DILUPBiCG: Solving for k, Initial residual = 5.23277e-05, Final
residual = 9.99146e-06, No Iterations 1
ExecutionTime = 1293.86 s ClockTime = 1338 s
forceCoeffs output:
Cd = 0.353696
Cl = 0.25129
Cm = -0.048046
forces output:
forces(pressure, viscous)((15.1715 127.22 167.353) (-0.128114
0.442893 12.3346))
moment(pressure, viscous)((-97.5863 -6.22688 19.4512) (0.195442
-0.388035 -0.0719235))
local:
forces(pressure, viscous)((15.1715 127.22 167.353) (-0.128114
0.442893 12.3346))
moment(pressure, viscous)((-97.5863 -6.22688 19.4512) (0.195442
-0.388035 -0.0719235))
十分妥当な値を計算できているようです。
この度はありがとうございました。
>
> もっと読む ≫
suzuswift
皆様明けましておめでとうございます。
たびたびの投稿失礼致します。
現在、メッシュ数を増やして同じ解析をしようとしております。
これまでの100万セルから390万セルに増やしました。
そうすると、100万セルで十分に収束していたケースが発散してしまうようになりました。
0ディレクトリ、constantディレクトリはこれまでと変えておりません。
systemディレクトリ内のfvSolutionを変更して試しているのですが、うまくいきません。
メッシュ数を増やした場合の注意点や残差許容値および緩和係数の設定でおかしなところがあればアドバイスお願い致します。
以下にcheckMesh,fvSolution,発散直前の結果を載せます。
以下のfvSolutionの他に、pのtoleranceを1e-8,relTolを0.001として実行したものもあるのですが、これも同様に発散
してしまいました。
何度も投稿してしまい申し訳ありませんが、何卒よろしくお願い致します。
//checkMesh
Create polyMesh for time = 0
Time = 0
Mesh stats
faces: 8073314
internal faces: 7497830
cells: 3892786
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
hexahedra: 0
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).
Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface
topology
connected)
ground 403591 202974 ok (non-closed singly
connected)
inlet 1632 924 ok (non-closed singly
connected)
outlet 1705 961 ok (non-closed singly
connected)
side 18420 10311 ok (non-closed singly
connected)
sky 2638 1416 ok (non-closed singly
Checking geometry...
Overall domain bounding box (0 0 -2.220446e-15) (7.0866 4.953
40.131)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Max cell openness = 1.817578e-16 OK.
Max aspect ratio = 395.1085 OK.
Minumum face area = 1.907878e-07. Maximum face area = 0.3217578.
Face area magnitudes OK.
Min volume = 1.456293e-10. Max volume = 0.06420923. Total volume
= 1402.505. Cell volumes OK.
Mesh non-orthogonality Max: 89.03308 average: 20.4782
*Number of severely non-orthogonal faces: 11549.
Non-orthogonality check OK.
<<Writing 11549 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 3.177031 OK.
Mesh OK.
End
//fvSolution
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-7;
nPreSweeps 0;
nPostSweeps 2;
cacheAgglomeration on;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 10;
mergeLevels 1;
};
U
{
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
};
k
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
};
omega
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
};
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
residualControl
{
p 1e-4;
U 1e-4;
"(k|omega)" 1e-4;
}
}
potentialFlow
{
nNonOrthogonalCorrectors 10;
}
relaxationFactors
{
p 0.05;
U 0.15;
k 0.15;
//epsilon 0.7;
omega 0.15;
}
//
************************************************************************* //
//結果
Time = 3196
DILUPBiCG: Solving for Ux, Initial residual = 0.019562, Final
residual = 0.000183616, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.0171144, Final
residual = 0.000708975, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0187588, Final
residual = 0.000108371, No Iterations 1
GAMG: Solving for p, Initial residual = 2.95398e-07, Final residual =
3.58104e-08, No Iterations 1
time step continuity errors : sum local = 1.45497e+56, global =
2.19705e+53, cumulative = 5.00841e+54
DILUPBiCG: Solving for omega, Initial residual = 0.0171573, Final
residual = 9.18016e-05, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.0451135, Final
residual = 0.000239463, No Iterations 1
bounding k, min: -6.39487e+89 max: 9.34941e+124 average: 4.3999e+119
ExecutionTime = 34744.3 s ClockTime = 35320 s
forceCoeffs output:
Cd = -1.78204e+118
Cl = -1.75211e+119
Cm = -4.68942e+118
forces output:
forces(pressure, viscous)((1.3485e+120 -8.90178e+121 -9.03347e
+120) (2.01812e+118 5.79009e+117 -1.98057e+118))
moment(pressure, viscous)((-9.50757e+121 4.01686e+120 -5.07941e
+121) (1.96455e+118 -8.72236e+117 1.74288e+118))
local:
forces(pressure, viscous)((1.3485e+120 -8.90178e+121 -9.03347e
+120) (2.01812e+118 5.79009e+117 -1.98057e+118))
moment(pressure, viscous)((-9.50757e+121 4.01686e+120 -5.07941e
+121) (1.96455e+118 -8.72236e+117 1.74288e+118))
>
> もっと読む ≫
ohbuchi
恐らく乱流モデルの計算が原因で発散しているのではないでしょうか？
残差が大きくなりだす直前の流れをparaviewで確認して見てください。
nutがとても大きな値になっている部分がないでしょうか？
ｆｖSolutionのRelaxationは十分小さい値になっているので、あとは
・k,omegaのdivSchemeをupwindにする
・乱流モデルを標準k-εにする
・メッシュのｓkewnessやアスペクト比をもう少し改善する
くらいしか思いつきません。結局はメッシュ品質が効くと思います。
cehckMeshでは問題ない様ですが、アスペクト比が大きいのが気に
なりました。ストレッチではなく、部分細分割で段階的に物体周りを
細かくして、アスペクト比の大きなメッシュを作らない方が良い様な
気がします。
外していたらごめんなさい。
>
> もっと読む ≫
suzuswift
お世話になっております。いつも書き込みありがとうございます。
divSchemeについてはもう既にupwindにしておりましたので、乱流モデルを
標準k-epsilonで実行してみました。
それでも発散してしまいましたので、メッシュのクオリティを疑っております。
skewnessやaspect比を抑えたメッシュが生成してみようと思います。