MRFSimpleFoamを使った計算について
Sakuma
ほぼ真ん中にファンが入っていて，そのファンを回転させて流れを計算）でトライしていますが，モデル作成の設定で分からない部分があります。
直管状のチューブは，半径45で長さが340，ファンは半径40で高さが15で，羽の枚数は37枚です。
blockMeshで基礎メッシュを切って，
vertices
(
(50 -50 180) // 0
(50 50 180) // 1
(-50 -50 180) // 2
(-50 50 180) // 3
(50 50 -200) // 4
(50 -50 -200) // 5
(-50 50 -200) // 6
(-50 -50 -200) // 7
);
blocks
(
hex (7 5 4 6 2 0 1 3) (10 10 38) simpleGrading (1 1 1)
);
edges ();
patches ();
mergePatchPairs ();
その後，snappyHexMeshでメッシュを作成しています。パッチ名の設定ですが，直管状のチューブの側面がSTATOR_Circle，流体の
流入の面をBACK_Circle，出口側の面をFRONT_ Circle，中のファンがROTOR_ROTORになっています（paraFoamで
チェックしました）。
チュートリアルでは，cellSetを行うようになっていましたので，system/cellSetDictのファイルを次のように設定しました。
name rotor;
action new;
topoSetSources ( zoneToCell { name ROTOR_ROTOR; } );
cellSetを実施すると次のようになりますが，上手くセルが設定できないようです。
Create polyMesh for time = 0
Reading cellSetDict
Backing up rotor into rotor_old
Set:rotor Size:0 Action:new
Adding all cells of cellZone ROTOR_ROTOR ...
--> FOAM Warning :
From function zoneToCell::combine(topoSet&, const bool)
in file sets/cellSources/zoneToCell/zoneToCell.C at line 87
Cannot find any cellZone named ROTOR_ROTOR
Valid names are
0
(
)
Writing rotor (size 0) to "constant/polyMesh/sets/rotor"
End
どのようにcellSetを設定すればよいのでしょうか。
また，境界条件の設定ですが，ファンにはチュートリアルのように回転を与えたいと思いますが，直管状のチューブの入口と出口の速度，圧力の境界条件です
が，入口側の速度条件が良く分かりません。
出口側：速度条件 type zeroGradient
圧力条件 type fixedValue;
value uniform 0.0;
入口：速度条件 ? どのような設定になるのでしょうか。
圧力条件 type zeroGradient;
ファンの回転により入口側の速度は変わるため，calculaedになるのでしょうか。
ご存じの方がおられましたら，アドバイスを頂きたく思います。
E.Mogura
> topoSetSources ( zoneToCell { name ROTOR_ROTOR; } );
チュートリアルの問題では、blockMeshを作成する段階で、cellの名前を区分してあるのに対して、
Sakumaさんの問題では、それが出来ていないからです。
OpenFOAM-1.6.x/applications/utileites/mesh/manipulation/cellSet/
の下に、
cellsetDict
のサンプルファイルがあるので、これを参考に、Sakumaさんの問題で、
ファンを含む回転領域を定義してやって下さい。
直管ということなので、
// Cells with centre within cylinder
cylinderToCell
{
p1 (0.0 0.0 z1); // start point on cylinder axis
p2 (0.0 0.0 ｚ2); // end point on cylinder axis
radius 42.5;
}
あたりを使えばよいと思います。
ファンの流れ（ｚ）方向位置に合わせて、ｚ１、ｚ２の値を変更してください。
> 入口側の速度条件が良く分かりません。
ファンの計算の場合は、通常、ファンのP-Q特性を参考に、
代表点を何点か選んで計算することになります。
設定した流量から、直管チューブの平均風速を算出して
その速度条件（fixedValue）で良いと思います。
> ファンの回転により入口側の速度は変わるため，calculaedになるのでしょうか。
出口圧力を固定しているので、入口圧力をzeroGradientにしておけば、
相応に変化した圧力特性を計算できることになります。
Sakuma
cellSetDictを次のように設定して
name rotor;
action new;
topoSetSources
(
// Cells with centre within cylinder
cylinderToCell
{
p1 (0.0 0.0 -13); // start point on cylinder axis
p2 (0.0 0.0 4); // end point on cylinder axis
radius 41.0;
}
);
cellSetを実行しました。
Create polyMesh for time = 0
Reading cellSetDict
Backing up rotor into rotor_old
Set:rotor Size:0 Action:new
Adding cells with centre within cylinder, with p1 = (0 0 -13), p2
= (0 0 4) and radius = 41
Writing rotor (size 403473) to "constant/polyMesh/sets/rotor"
End
OpenFOAMのユーティリティーであるMesh manipulation群の使い方を理解していないので間違っている
かもしれませんが，cellSetDictのtopoSetSourcesで指定した方法で回転部分となるファンの表面部のセル
を選択していると思います。ただ，snappyHexMeshを切った時に出来るboundaryとその中に回転部分とな
るファン部分のpatch情報が書かれていますが，その情報を利用してcellSetで回転部分となるファンの表面
部を選択して，設定する事はできないでしょうか。
ROTOR_ROTOR
{
type wall;
nFaces 70366;
startFace 1186581;
}
今回のテストケースでは直管状のチューブのほぼ真ん中にあるため，お教え頂いたtopoSetSources の設定を
cylinderToCellにして，目的とする領域を選択できましたが，実際のものではファンの近傍にシュラウドがあった
りすると思いますが，このような場合には回転部分以外を選択してしまう可能性がありますが，このような時はど
のように対処するのでしょうか。
速度や圧力の境界条件を設定する場合には，DEXCS OpenFOAM2009のlauncherSimpleFoam.pyを見せ
て頂き，pyFoamCreateBoundaryPatches.pyを使ってパッチを取り込み，初期値を設定しています。これに相
当する便利なコマンドは存在しないのでしょうか。
cfd-onlineの下記のスレッドを見ましたが良く分かりません。
http://www.cfd-online.com/Forums/openfoam-solving/58821-cellzons-mrfsimplefoam.html
ファン解析の境界条件ですが，ファンをある回転速度で回した場合，どの程度流量が見込めるかなどを計算する
場合には，陽的に流量や吸い込み側の速度を与えることは難しいと思います。
このような場合には，
入力側(inlet)の
速度条件：type zeroGradient
圧力条件：type fixedValue;
value uniform 0.0;
出口側(outlet)の
速度条件 type zeroGradient
圧力条件 type fixedValue;
value uniform 0.0;
となるのでしょうか。
まだまだ流体解析の初学者の質問の域を出なく，色々とお教え頂く事が多いと思いますが，
よろしくお願いいたします。
> > ご存じの方がおられましたら，アドバイスを頂きたく思います。- 引用テキストを表示しない -
>
> - 引用テキストを表示 -
E.Mogura
＞cellSetDictのtopoSetSourcesで指定した方法で回転部分となるファンの表面部のセル
＞を選択していると思います。ただ，snappyHexMeshを切った時に出来るboundaryとその中に回転部分とな
＞るファン部分のpatch情報が書かれていますが，その情報を利用してcellSetで回転部分となるファンの表面
＞部を選択して，設定する事はできないでしょうか。
＞
cellSetでやっていることは、表面部分のセルだけでなく、
表面から離れて、指定した円筒面までの範囲のセルを全て選択しています。これは、
・指定表面に回転速度境界を付与するだけでなく、
・指定した流体領域全体にわたって遠心力（コリオリ）力を付与
する為です。
前者だけであれば、確かにご指摘の方法が考えられなくもありませんが、
後者をどういう範囲で設定するかに関しては正解があるというものでもないようで、
解析する人の勘コツでやるしかなさそうです。
＞実際のものではファンの近傍にシュラウドがあった
＞りすると思いますが，このような場合には回転部分以外を選択してしまう可能性がありますが，このような時はど
＞のように対処するのでしょうか。
静止部分を除外する形で回転領域を定義する形状を作成し、STLファイル保存。
cellSetDictで、topoSetSourcesとして、surfaceToCellを使ってやっています。
surfaceToCell
{
file "www.avl.com-geometry.stl";
outsidePoints ((-99 -99 -59)); // definition of
outside
includeCut false; // cells cut by
surface
includeInside false; // cells not on
outside of surf
includeOutside false; // cells on outside of
surf
nearDistance -1; // cells with centre
near surf
// (set to -1 if not
used)
curvature 0.9; // cells within
nearDistance
// and near surf
curvature
// (set to -100 if not
used)
}
＞陽的に流量や吸い込み側の速度を与えることは難しいと思います。
＞
＞このような場合には，
＞入力側(inlet)の
＞速度条件：type zeroGradient
＞ 圧力条件：type fixedValue;
＞ value uniform 0.0;
＞
＞
＞出口側(outlet)の
＞速度条件 type zeroGradient
＞ 圧力条件 type fixedValue;
＞ value uniform 0.0;
＞となるのでしょうか
はい、それでもいけると思います。
収束性が思わしくないような場合には、
全圧型の境界条件（totalPressure）も試してみるのが良いかと。
＞ 速度や圧力の境界条件を設定する場合には，DEXCS OpenFOAM2009のlauncherSimpleFoam.pyを見せ
＞て頂き，pyFoamCreateBoundaryPatches.pyを使ってパッチを取り込み，初期値を設定しています。これに相
＞当する便利なコマンドは存在しないのでしょうか。
あったら、とっくにDEXCSで採用してたんですがねぇ。。。
> cfd-onlineの下記のスレッドを見ましたが良く分かりません。http://www.cfd-online.com/Forums/openfoam-solving/58821-cellzons-mrfs...
> > - 引用テキストを表示 -- 引用テキストを表示しない -
>
> - 引用テキストを表示 -
IMANO Masashi
今野です。
回転部分として、ファンの表面の第一セルのみを指定すれば良いのであれば、
system/faceSetDict
---
name faceSet_ROTOR_ROTOR;
action new;
topoSetSources
(
patchToFace
{
name "ROTOR_ROTOR";
}
);
---
として、faceSetを実行し、
system/cellSetDict
---
name cellSet_ROTOR_ROTOR;
action new;
topoSetSources
(
faceToCell
{
set faceSet_ROTOR_ROTOR;
//option owner; // ,, owner
option any; // cell with any face in faceSet
}
);
---
として、cellSetを実行すれば良いと思います。
なお、界面のfaceのownerは第一セルですから、上記のoptionはownerでも良い
のですが、後述の目的のためanyにしています。
また、上記で第一セルのcellSetを取得した後、
system/faceSetDict
---
name faceSet_ROTOR_ROTOR;
action new;
cellToFace
{
set cellSet_ROTOR_ROTOR;
option all; // All faces of cells
}
---
として、faceSet -> cellSet を実行すれば、第二セルも加えたcellSetが取得
できる気がします。
さらに、faceSet -> cellSetの実行を繰り返せば、第三、第四セル、と広げら
れると思います。
また、ファンの表面からある距離以内の格子を指定するだけで良ければ、
cellSetDict内のsurfaceToCellでnearDistanceオプションを指定すれば
良いのではないでしょうか？
E.Mogura
E.Mogura です
＞として、faceSet -> cellSet を実行すれば、第二セルも加えたcellSetが取得
＞できる気がします。
＞さらに、faceSet -> cellSetの実行を繰り返せば、第三、第四セル、と広げら
＞れると思います。
なるほど、面白いやり方もあるんですね。
cellSetの方法としては、たいへん参考になります。
ただ、今回のMRFでの使用目的を考えた場合には、
この方法ではうまくいかないと思います。
というのは、MRFでcellSet指定する（コリオリ力を付与する）領域と、
そうでない静止領域との境界面は、できるだけフラットな面になるように
構成する必要があるからです。
そうしないと、フラットでない部分で圧力不連続が生じてしまい、
ファンの圧力特性が実測よりもかなり低下してしまうことになります。
（これは実際の失敗経験例です）
以上、参考まで
Sakuma
今，お教え頂いた方法で色々と触っているところです。
E.Moguraさんの投稿の中に，
”MRFでcellSet指定する（コリオリ力を付与する）領域と、そうでない静止領域との境界面は、
できるだけフラットな面になるように”
とありますが，
前にお教え頂いた
topoSetSources
(
// Cells with centre within cylinder
cylinderToCell
{
p1 (0.0 0.0 -13); // start point on cylinder axis
p2 (0.0 0.0 4); // end point on cylinder axis
radius 41.0;
}
);
を使えば，境界面はフラットな面に近くなるのでしょうか。このcylinderToCell
で指定した部分のメッシュはparaFoamで見る事ができるのでしょうか。
OpenFOAMは色々と奥が深いですね。まだまだ知らない事が山ほどありまして
前途多難ですが，あきらめずにやっていきたいと思います。
周りにOpenFOAMを使っている方がおられない方や，聞く人がおられない方は
どのようにOpenFOAMを使っておられるのでしょうか？
> > 良いのではないでしょうか？- 引用テキストを表示しない -
>
> - 引用テキストを表示 -
E.Mogura
foamToVTK -cellSet rotor
で指定領域（rotor）のセル情報がvtk化されるので、
これで出来たvtkファイルをparafoam（またはparaview）で見て確認できます。
snappyで作成したメッシュを使う限り、
> cylinderToCell
を使っても、円柱の外周境界はどうしてもギザギザになってしまいますが、
（軸流ファンの計算の場合に主流通過面となる）円柱の端面部分は工夫次第
でフラットな面を作成できるので、そのあたりで妥協して使っています。
> > - 引用テキストを表示 -- 引用テキストを表示しない -
>
> - 引用テキストを表示 -
Sakuma
foamToVTK -cellSet rotor を行い，paraFoamで読み込んでみました。
ファンの周りの領域が選択されていましたが，結構ギサギザになっている
事を確認しました。
cellSetの次に行うsetsToZonesの引数ですが，これは何を意味しているので
しょうか。1.6のマニュアルU-90ページに簡単な説明がありますが，あまり
理解できませんし，引数の記述等はないため困っています。
グーグルで検索しても中々理解が進みません。このようなコマンドのもう少し
詳しい使い方などは，どこかのサイトにあるのでしょうか？
Sakuma
良く分からないエラーが出て，OpenFOAMが止まりました。まず，どの当たりの設
定が悪いのでしょうか。境界条件は色々とネットなどを見て設定したつもりです。
境界条件として
圧力
流れの入力側：zeroGradient
流れの出口側：fixedValueで0
ローターとステータはzeroGradient
速度
流れの入力側：fixedValueで1m/s
流れの出口側：zeroGradient
ローターとステータはfixedValueで0
乱入モデルはk-εで
kは
流れの入力側：fixedValueで0.375
流れの出口側：zeroGradient
ローターとステータはkqRWallFunction，fixedValueで0.375
εは
流れの入力側：fixedValueで14.855
流れの出口側：zeroGradient
ローターとステータはepsilonWallFunction，fixedValueで14.855
乱流粘性(nut)は
流れの入力側：calculatedでuniform 0
流れの出口側：calculatedでuniform 0
ローターとステータはnutWallFunction，uniform 0
MRFZonesの回転数は209としています。
Starting time loop
Time = 1
DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual =
1.725005e-06, No Iterations 6
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual =
1.813184e-06, No Iterations 6
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual =
9.491987e-06, No Iterations 5
DICPCG: Solving for p, Initial residual = 1, Final residual =
9.533176e-07, No Iterations 150
time step continuity errors : sum local = 0.0001180144, global =
-7.01167e-07, cumulative = -7.01167e-07
DILUPBiCG: Solving for epsilon, Initial residual = 0.229383, Final
residual = 5.014963e-06, No Iterations 6
DILUPBiCG: Solving for k, Initial residual = 1, Final residual =
2.643756e-06, No Iterations 7
ExecutionTime = 31.39 s ClockTime = 31 s
Time = 2
DILUPBiCG: Solving for Ux, Initial residual = 0.3435081, Final
residual = 2.292135e-06, No Iterations 6
DILUPBiCG: Solving for Uy, Initial residual = 0.3394074, Final
residual = 2.736647e-06, No Iterations 6
DILUPBiCG: Solving for Uz, Initial residual = 0.4738587, Final
residual = 1.599999e-06, No Iterations 8
DICPCG: Solving for p, Initial residual = 0.2902264, Final residual =
9.053265e-07, No Iterations 208
time step continuity errors : sum local = 0.0006817565, global =
-3.171762e-05, cumulative = -3.241879e-05
DILUPBiCG: Solving for epsilon, Initial residual = 0.2793774, Final
residual = 1.925245e-06, No Iterations 5
bounding epsilon, min: -100355.1 max: 1.162648e+08 average: 161309.6
DILUPBiCG: Solving for k, Initial residual = 0.9624649, Final
residual = 4.950803e-06, No Iterations 6
ExecutionTime = 59.24 s ClockTime = 59 s
Time = 3
DILUPBiCG: Solving for Ux, Initial residual = 0.4732396, Final
residual = 2.970225e-06, No Iterations 6
DILUPBiCG: Solving for Uy, Initial residual = 0.4751867, Final
residual = 2.197725e-06, No Iterations 6
DILUPBiCG: Solving for Uz, Initial residual = 0.570641, Final
residual = 3.559255e-06, No Iterations 6
DICPCG: Solving for p, Initial residual = 0.2797082, Final residual =
9.914247e-07, No Iterations 231
time step continuity errors : sum local = 0.002397115, global =
1.210766e-05, cumulative = -2.031112e-05
DILUPBiCG: Solving for epsilon, Initial residual = 0.2208734, Final
residual = 3.535682e-06, No Iterations 4
bounding epsilon, min: -5.204303e+09 max: 1.317281e+11 average:
6.959981e+07
DILUPBiCG: Solving for k, Initial residual = 0.8626755, Final
residual = 2.274892e-06, No Iterations 7
ExecutionTime = 88.2 s ClockTime = 88 s
Time = 4
DILUPBiCG: Solving for Ux, Initial residual = 0.6009585, Final
residual = 3.752826e-06, No Iterations 5
DILUPBiCG: Solving for Uy, Initial residual = 0.6253781, Final
residual = 3.372857e-06, No Iterations 5
DILUPBiCG: Solving for Uz, Initial residual = 0.7957528, Final
residual = 7.990759e-06, No Iterations 5
DICPCG: Solving for p, Initial residual = 0.09645187, Final residual
= 9.323403e-07, No Iterations 255
time step continuity errors : sum local = 0.02298909, global =
0.002187293, cumulative = 0.002166982
DILUPBiCG: Solving for epsilon, Initial residual = 0.01631645, Final
residual = 3.001472e-06, No Iterations 3
bounding epsilon, min: -1.803166e+12 max: 1.389035e+14 average:
1.942541e+10
DILUPBiCG: Solving for k, Initial residual = 0.4641816, Final
residual = 7.399794e-06, No Iterations 4
bounding k, min: -26.15993 max: 1.407792e+08 average: 144060.4
ExecutionTime = 117.62 s ClockTime = 118 s
Time = 5
DILUPBiCG: Solving for Ux, Initial residual = 0.386897, Final
residual = 3.531711e-06, No Iterations 6
DILUPBiCG: Solving for Uy, Initial residual = 0.3655837, Final
residual = 3.180235e-06, No Iterations 6
DILUPBiCG: Solving for Uz, Initial residual = 0.4326544, Final
residual = 2.530188e-06, No Iterations 6
DICPCG: Solving for p, Initial residual = 0.04111225, Final residual
= 9.296733e-07, No Iterations 351
time step continuity errors : sum local = 4.902578e+15, global =
2.797307e+14, cumulative = 2.797307e+14
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual
= 2.292685e-06, No Iterations 9
bounding epsilon, min: -4.677792e+11 max: 2.667737e+31 average:
1.052605e+27
DILUPBiCG: Solving for k, Initial residual = 1, Final residual =
8.579102e-06, No Iterations 5
ExecutionTime = 156.06 s ClockTime = 157 s
Time = 6
DILUPBiCG: Solving for Ux, Initial residual = 0.5481319, Final
residual = 7.549354e-06, No Iterations 5
DILUPBiCG: Solving for Uy, Initial residual = 0.6543282, Final
residual = 8.321401e-06, No Iterations 5
DILUPBiCG: Solving for Uz, Initial residual = 0.4962678, Final
residual = 1.122597e-06, No Iterations 6
DICPCG: Solving for p, Initial residual = 0.915788, Final residual =
0.007658455, No Iterations 1001
time step continuity errors : sum local = 3.81622e+23, global =
2.584483e+17, cumulative = 2.58728e+17
DILUPBiCG: Solving for epsilon, Initial residual = 0.1893555, Final
residual = 7.257969e-06, No Iterations 2
bounding epsilon, min: -1.88556e+39 max: 7.999178e+54 average:
2.328306e+50
DILUPBiCG: Solving for k, Initial residual = 0.8553349, Final
residual = 6.468891e-06, No Iterations 6
bounding k, min: -1.436684e+10 max: 1.012193e+43 average: 1.043081e+39
ExecutionTime = 239.88 s ClockTime = 242 s
Time = 7
DILUPBiCG: Solving for Ux, Initial residual = 2.818142e-05, Final
residual = 3.581528e-06, No Iterations 3
DILUPBiCG: Solving for Uy, Initial residual = 8.149942e-07, Final
residual = 8.149942e-07, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 1.08042e-06, Final
residual = 1.08042e-06, No Iterations 0
#0 Foam::error::printStack(Foam::Ostream&) in "/home/sakuramaru/
OpenFOAM/OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 Foam::DICPreconditioner::calcReciprocalD(Foam::Field<double>&,
Foam::lduMatrix const&) in "/home/sakuramaru/OpenFOAM/OpenFOAM-1.6.x/
lib/linuxGccDPOpt/libOpenFOAM.so"
#4 Foam::DICPreconditioner::DICPreconditioner(Foam::lduMatrix::solver
const&, Foam::dictionary const&) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#5
Foam::lduMatrix::preconditioner::addsymMatrixConstructorToTable<Foam::DICPreconditioner>::New(Foam::lduMatrix::solver
const&, Foam::dictionary const&) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#6 Foam::lduMatrix::preconditioner::New(Foam::lduMatrix::solver
const&, Foam::dictionary const&) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#7 Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&,
unsigned char) const in "/home/sakuramaru/OpenFOAM/OpenFOAM-1.6.x/lib/
linuxGccDPOpt/libOpenFOAM.so"
#8 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/
sakuramaru/OpenFOAM/OpenFOAM-1.6.x/lib/linuxGccDPOpt/
libfiniteVolume.so"
#9 main in "/home/sakuramaru/OpenFOAM/sakuramaru-1.6.x/applications/
bin/linuxGccDPOpt/MRFSimpleFoam"
#10 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#11 Foam::regIOobject::writeObject(Foam::IOstream::streamFormat,
Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const
in "/home/sakuramaru/OpenFOAM/sakuramaru-1.6.x/applications/bin/
linuxGccDPOpt/MRFSimpleFoam"
Floating point exception
ここで止まっています。詳しい方がいらしたらアドバイスをお願いします。
IMANO Masashi
今野です。
ログを見ると、最初に
> Time = 2
> DICPCG: Solving for p, Initial residual = 0.2902264, Final residual =
> 9.053265e-07, No Iterations 208
> time step continuity errors : sum local = 0.0006817565, global =
> -3.171762e-05, cumulative = -3.241879e-05
> DILUPBiCG: Solving for epsilon, Initial residual = 0.2793774, Final
> residual = 1.925245e-06, No Iterations 5
> bounding epsilon, min: -100355.1 max: 1.162648e+08 average: 161309.6
と、epsilonが発散して、最終的には
Time = 6
DICPCG: Solving for p, Initial residual = 0.915788, Final residual =
0.007658455, No Iterations 1001
time step continuity errors : sum local = 3.81622e+23, global =
2.584483e+17, cumulative = 2.58728e+17
DILUPBiCG: Solving for epsilon, Initial residual = 0.1893555, Final
residual = 7.257969e-06, No Iterations 2
bounding epsilon, min: -1.88556e+39 max: 7.999178e+54 average:
2.328306e+50
DILUPBiCG: Solving for k, Initial residual = 0.8553349, Final
residual = 6.468891e-06, No Iterations 6
bounding k, min: -1.436684e+10 max: 1.012193e+43 average: 1.043081e+39
と、kや連続の式の誤差(continuity errors)も破綻していますね。
あと、圧力に関する反復回数も多い気がします。
k,epsilonの発散は離散スキームが不安定であることも多いので、fvSchemesを
貼ってください。pの反復回数が多いのは、solverの設定のせいかもしれない
ので、fvSolutionも貼ってください。
Sakuma
計算に利用しているfvSchemesは，
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwindV Gauss linear;
div(phi,k) Gauss upwind;
div((nuEff*dev(grad(U).T()))) Gauss linear;
div(phi,epsilon) Gauss upwind;
div(phi,omega) Gauss upwind;
}
laplacianSchemes
{
default Gauss linear corrected;
// default Gauss linear limited 0.5;
// default Gauss linear limited 0.333;
}
interpolationSchemes
{
default linear;
interpolate(U) linear;
}
snGradSchemes
{
default corrected;
}
fluxRequired
{
default no;
p;
}
fvSolutionは，
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0;
}
U
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}
k
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}
epsilon
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}
omega
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
};
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
}
relaxationFactors
{
p 0.3;
U 0.7;
k 0.7;
epsilon 0.7;
omega 0.7;
}
両方とも特に考えずに，チュートリアルなどから適当に持ってきています。
（設定の基準が良く分からないため）
よろしくお願いいたします。
IMANO Masashi
今野です．
fvSchemesについてですが，現在k,epsilonの移流項スキームは一番安定な
upwindになっていますから，これが発散の原因ではないですね．
他のスキームも特に問題はないですが，発散の原因がわかるまでは，Uの
移流項スキームもupwindにしておいたほうが良いと思います．原因が取り除か
れたら，Uやk,epsilonの移流項スキームを，より高次精度なスキームに変えて
いってください．
個人的にはlinearUpwind系よりも，MRFSimpleFoam/mixerVessel2Dのチュート
リアルでも用いられている，limitedLinear系のほうがお勧めです．
次にfvSolutionを見てみますと，U,k,epsilonのrelaxationFactorsが少し大きすぎる
と思いますので，これをとりあえず0.3位にしてみてください．
あと，solversで全物理量に対してrelTolが0になっていますが，SIMPLE法で定
常解を求める場合には通常0.1程度(pは0.01～0.05)にしたほうが，スムーズに
収束すると思います．
最後に，以下の点も確認してください．
・そもそもsimpleFoam(回転数を0の状態)で解けるか?
・入口側の速度が与えられいるので風量が決まっているが，それで良いのか?
Sakuma
頂きましたアドバイスを元に計算設定を修正しました。少し長くなりますがご容赦ください。
解析設定ですが，
圧力
流れの入力側：fixedValueで0
流れの出口側：fixedValueで0
ローターとステータはzeroGradient
速度
流れの入力側：zeroGradient
流れの出口側：zeroGradient
ローターとステータはfixedValueで0
乱入モデルはk-εで
kは
流れの出口側：zeroGradient
ローターとステータはkqRWallFunction，fixedValueで0.375
εは
流れの出口側：zeroGradient
ローターとステータはepsilonWallFunction，fixedValueで14.855
乱流粘性(nut)は
流れの入力側：calculatedでuniform 0
流れの出口側：calculatedでuniform 0
ローターとステータはnutWallFunction，uniform 0
fvSolutionは，
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
U
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
}
k
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
}
epsilon
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
}
omega
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
};
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
}
relaxationFactors
{
p 0.3;
U 0.3;
k 0.3;
epsilon 0.3;
omega 0.3;
}
fvSchemesは，
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss upwind;
div(phi,k) Gauss upwind;
div((nuEff*dev(grad(U).T()))) Gauss linear;
div(phi,epsilon) Gauss upwind;
div(phi,omega) Gauss upwind;
}
laplacianSchemes
{
default Gauss linear corrected;
// default Gauss linear limited 0.5;
// default Gauss linear limited 0.333;
}
interpolationSchemes
{
default linear;
interpolate(U) linear;
}
snGradSchemes
{
default corrected;
}
fluxRequired
{
default no;
p;
}
ここで，MRFZonesの回転数は20.9とした場合は最後まで計算できます。
Time = 2000
DILUPBiCG: Solving for Ux, Initial residual = 4.793028e-06, Final
residual = 4.793028e-06, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 4.795481e-06, Final
residual = 4.795481e-06, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 1.08002e-05, Final
residual = 2.892438e-07, No Iterations 1
DICPCG: Solving for p, Initial residual = 3.520504e-05, Final
residual = 1.658821e-06, No Iterations 11
time step continuity errors : sum local = 4.828023e-06, global =
-1.466306e-06, cumulative = -0.01130605
DILUPBiCG: Solving for epsilon, Initial residual = 4.23838e-05, Final
residual = 8.182859e-07, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.0001407032, Final
residual = 1.741196e-06, No Iterations 1
ExecutionTime = 15037.91 s ClockTime = 15134 s
End
しかし，回転数を10倍の209とするとエラーが直ぐに出てしまいます。以前にご指摘を頂いたように計算が破たんしています。
Time = 15
DILUPBiCG: Solving for Ux, Initial residual = 1.648609e-14, Final
residual = 1.648609e-14, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 9.768459e-15, Final
residual = 9.768459e-15, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 6.242267e-15, Final
residual = 6.242267e-15, No Iterations 0
DICPCG: Solving for p, Initial residual = 6.319817e-25, Final
residual = 6.319817e-25, No Iterations 0
time step continuity errors : sum local = 3.397966e+121, global =
6.187454e+104, cumulative = 6.187454e+104
#0 Foam::error::printStack(Foam::Ostream&) in "/home/sakuramaru/
OpenFOAM/OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 Foam::multiply(Foam::Field<double>&, Foam::UList<double> const&,
Foam::UList<double> const&) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#4 void Foam::multiply<Foam::fvPatchField,
Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh> const&, Foam::GeometricField<double,
Foam::fvPatchField, Foam::volMesh> const&) in "/home/sakuramaru/
OpenFOAM/OpenFOAM-1.6.x/lib/linuxGccDPOpt/
libincompressibleRASModels.so"
#5 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh> > Foam::operator*<Foam::fvPatchField,
Foam::volMesh>(Foam::tmp<Foam::GeometricField<double,
Foam::fvPatchField, Foam::volMesh> > const&,
Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>
const&) in "/home/sakuramaru/OpenFOAM/OpenFOAM-1.6.x/lib/linuxGccDPOpt/
libincompressibleRASModels.so"
#6 Foam::incompressible::RASModels::kEpsilon::correct() in "/home/
sakuramaru/OpenFOAM/OpenFOAM-1.6.x/lib/linuxGccDPOpt/
libincompressibleRASModels.so"
#7 main in "/home/sakuramaru/OpenFOAM/sakuramaru-1.6.x/applications/
bin/linuxGccDPOpt/MRFSimpleFoam"
#8 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#9 Foam::regIOobject::writeObject(Foam::IOstream::streamFormat,
Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const
in "/home/sakuramaru/OpenFOAM/sakuramaru-1.6.x/applications/bin/
linuxGccDPOpt/MRFSimpleFoam"
Floating point exception
また，MRFZonesの回転数は20.9とし，2CUPで並列計算すると途中でこれも途中で破たんします。
Time = 468
DILUPBiCG: Solving for Ux, Initial residual = 0.297535, Final
residual = 0.0006827055, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.2962824, Final
residual = 0.0006729731, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.2193794, Final
residual = 0.0003174518, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.9999993, Final residual =
0.04622674, No Iterations 4
time step continuity errors : sum local = 8.368492e+28, global =
-3.611277e+28, cumulative = -3.611277e+28
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual
= 0.04408848, No Iterations 13
bounding epsilon, min: -2.578166e+57 max: 6.20081e+56 average:
-6.514203e+51
DILUPBiCG: Solving for k, Initial residual = 0.7787079, Final
residual = 0.001303681, No Iterations 4
bounding k, min: -8.035473e+47 max: 1.109702e+45 average: -2.315793e
+42
ExecutionTime = 2607.57 s ClockTime = 2653 s
Time = 469
DILUPBiCG: Solving for Ux, Initial residual = 5.289698e-06, Final
residual = 5.289698e-06, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 4.204676e-05, Final
residual = 4.630337e-06, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 1.450399e-05, Final
residual = 5.107689e-07, No Iterations 1
[0] [1] ##0 Foam::error::printStack(Foam::Ostream&)0
Foam::error::printStack(Foam::Ostream&) in "/home/sa in kuramaru/
Open"FO/hAM/OopenmFOAM-1.6.x/lib/le/sakurianuxmGccaru/OpeDPnOFpt/
libOpenFOAM.sOAM/OpenFoO"A
[1] M-#1 Foam::sigFpe::sigFpeHandler(int)1.6.x/lib/linuxGccDPOpt/
libOpenFOAM.so"
[0] #1 Foam::sigFpe::sigFpeHandler(int) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
[0] #2 Uninterpreted:
[0] #3 Foam::divide(Foam::Field<double>&, double const&,
Foam::UList<double> const&) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
[1] #2 Uninterpreted:
[1] #3 Foam::divide(Foam::Field<double>&, double const&,
Foam::UList<double> const&) in "/home/sakuramaru/OpenFOAM/
OpenFOAM-1.6.x/lib/linuxGccDPOpt/libOpenFOAM.so"
[1] #4 in "/home/sakuramaru/OpenFOAM/OpenFOAM-1.6.x/lib/
linuxGccDPOpt/libOpenFOAM.so"
[0] #4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh> > Foam::operator/<Foam::fvPatchField,
Foam::volMesh>(Foam::dimensioned<double> const&,
Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh> > const&)Foam::tmp<Foam::GeometricField<double,
Foam::fvPatchField, Foam::volMesh> > Foam::operator/
<Foam::fvPatchField, Foam::volMesh>(Foam::dimensioned<double> const&,
Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh> > const&) in "/home/sakuramaru/OpenFOAM/
sakuramaru-1.6.x/applications/bin/linuxGccDPOpt/MRFSimpleFoam"
[1] #5 in "/home/sakuramaru/OpenFOAM/sakuramaru-1.6.x/applications/
bin/linuxGccDPOpt/MRFSimpleFoam"
[0] #5 mainmain in "/home/sakuramaru/OpenFOAM/sakuramaru-1.6.x/
applications/bin/linuxGccDPOpt in "/home/sakur/aMRFSimmpleaFoarm"
u/Ope[1] #nFOAM/6 __libc_start_mainsakuramaru-1.6.x/applications/bin/
linuxGccDPOpt/MRFSimpleFoam"
[0] #6 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
[1] #7 in "/lib/tls/i686/cmov/libc.so.6"
[0] #7 Foam::regIOobject::writeObject(Foam::IOstream::streamFormat,
Foam::IOstream::versionNumber, Foam::IOstream::compressionType)
constFoam::regIOobject::writeObject(Foam::IOstream::streamFormat,
Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const
in "/home/sakuramaru/OpenFOAM/sakuramaru-1.6.x/applications/bin/
linuxGccDPOpt/MRFSimpleFoam"
in "/home/sakuramaru/O[ubuntu-vm:18555] *** Process received signal
***
penFOAM/sa[ubuntu-vm:18555] Signal: Floating point exception (8)
[ubuntu-vm:18555] Signal code: (-6)
[ubuntu-vm:18555] Failing at address: 0x487b
kuramaru-1.6.x/application[ubuntu-vm:18555] [ 0] [0xb7756440]
s/bi[ubuntu-vm:18555] [ 1] /home/sakuramaru/OpenFOAM/OpenFOAM-1.6.x/
lib/linuxGccDPOpt/libOpenFOAM.so(_ZN4Foam6sigFpe13sigFpeHandlerEi
+0x61) [0xb68bda21]
n/linuxG[ubuntu-vm:18555] [ 2] [0xb7756420]
ccDP[ubuntu-vm:18555] [ 3]
MRFSimpleFoam(_ZN4FoamdvINS_12fvPatchFieldENS_7volMeshEEENS_3tmpINS_14GeometricFieldIdT_T0_EEEERKNS_11dimensionedIdEERKS8_
+0x364) [0x808fc34]
[ubuntu-vm:18555] [ 4] MRFSimpleFoam [0x805ab24]
O[ubuntu-vm:18555] [ 5] /lib/tls/i686/cmov/libc.so.6(__libc_start_main
+0xe0) [0xb631b450]
pt/MRF[ubuntu-vm:18555] [ 6]
MRFSimpleFoam(_ZNK4Foam11regIOobject11writeObjectENS_8IOstream12streamFormatENS1_13versionNumberENS1_15compressionTypeE
+0xad) [0x80599a1]
S[ubuntu-vm:18555] *** End of error message ***
impleFoam"
[ubuntu-vm:18554] *** Process received signal ***
[ubuntu-vm:18554] Signal: Floating point exception (8)
[ubuntu-vm:18554] Signal code: (-6)
[ubuntu-vm:18554] Failing at address: 0x487a
[ubuntu-vm:18554] [ 0] [0xb7750440]
[ubuntu-vm:18554] [ 1] /home/sakuramaru/OpenFOAM/OpenFOAM-1.6.x/lib/
linuxGccDPOpt/libOpenFOAM.so(_ZN4Foam6sigFpe13sigFpeHandlerEi+0x61)
[0xb68b7a21]
[ubuntu-vm:18554] [ 2] [0xb7750420]
[ubuntu-vm:18554] [ 3]
MRFSimpleFoam(_ZN4FoamdvINS_12fvPatchFieldENS_7volMeshEEENS_3tmpINS_14GeometricFieldIdT_T0_EEEERKNS_11dimensionedIdEERKS8_
+0x364) [0x808fc34]
[ubuntu-vm:18554] [ 4] MRFSimpleFoam [0x805ab24]
[ubuntu-vm:18554] [ 5] /lib/tls/i686/cmov/libc.so.6(__libc_start_main
+0xe0) [0xb6315450]
[ubuntu-vm:18554] [ 6]
MRFSimpleFoam(_ZNK4Foam11regIOobject11writeObjectENS_8IOstream12streamFormatENS1_13versionNumberENS1_15compressionTypeE
+0xad) [0x80599a1]
[ubuntu-vm:18554] *** End of error message ***
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 18554 on node ubuntu-vm
exited on signal 8 (Floating point exception).
--------------------------------------------------------------------------
Killing PID 18547
このようなエラーの場合，どのような所かた手を付けていったらよろしいでしょうか？
MRFSimpleFoamの勉強用に適当にモデル作って計算を実施しているため，今野さんがご指摘されたようにそもそも問題が
正しいのかという問題はあるかと思いますが。
E.Mogura
--
このメールは Google グループのグループ「OpenFOAM」の登録者に送られています。
このグループに投稿するには、open...@googlegroups.com にメールを送信してください。
このグループから退会するには、openfoam+u...@googlegroups.com にメールを送信してください。
詳細については、http://groups.google.com/group/openfoam?hl=ja からこのグループにアクセスしてください。
Sakuma
k-omegaで計算をしましたが，Time = 8で発散してしまいました。
k-epsillonではTime=15で発散しましたのであまり差がないようでした。
モデルをkOmegaSSTに変更して計算を行いました。（kOmegaSSTを次に設定をしたのは，
汎用ソフトでは良く利用されると聞いた事があるためです。特に理由があってではありません）
1CUP，2CUP共に設定した最後のタイムステップ2000まで計算出来ました。
始めのうちは，残差が大きい状態でしたが。
1CUPの場合
Time = 2000
DILUPBiCG: Solving for Ux, Initial residual = 8.901176e-07, Final
residual = 8.901176e-07, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 8.341001e-07, Final
residual = 8.341001e-07, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 3.562843e-06, Final
residual = 3.562843e-06, No Iterations 0
DICPCG: Solving for p, Initial residual = 1.498378e-05, Final
residual = 9.072167e-07, No Iterations 8
time step continuity errors : sum local = 4.012237e-05, global =
-1.240391e-06, cumulative = 0.7716906
DILUPBiCG: Solving for omega, Initial residual = 9.349589e-06, Final
residual = 9.349589e-06, No Iterations 0
DILUPBiCG: Solving for k, Initial residual = 2.37624e-05, Final
residual = 2.883788e-07, No Iterations 1
bounding k, min: -2.247764e-08 max: 0.8221259 average: 0.0208806
ExecutionTime = 16372.09 s ClockTime = 16503 s
2CUPの場合
Time = 2000
DILUPBiCG: Solving for Ux, Initial residual = 2.039291e-06, Final
residual = 2.039291e-06, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 2.034913e-06, Final
residual = 2.034913e-06, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 7.32111e-06, Final
residual = 7.32111e-06, No Iterations 0
DICPCG: Solving for p, Initial residual = 3.235371e-05, Final
residual = 1.485328e-06, No Iterations 17
time step continuity errors : sum local = 6.490386e-05, global =
-8.043267e-07, cumulative = -0.09838417
DILUPBiCG: Solving for omega, Initial residual = 8.979081e-06, Final
residual = 8.979081e-06, No Iterations 0
DILUPBiCG: Solving for k, Initial residual = 3.158518e-05, Final
residual = 1.057272e-06, No Iterations 1
bounding k, min: -2.594322e-08 max: 0.8293826 average: 0.01997287
ExecutionTime = 11461.14 s ClockTime = 11763 s
設定するRANSの乱流モデルで，これほどの差が出るとは思いませんでした。
では，どのようにモデルを選択していったらよいのでしょうか？ 何か選択基準があるのでしょうか。
On 4月6日, 午前5:12, "E.Mogura" <seikun...@gmail.com> wrote:
> E.Mogura です
>
> 乱流計算で発散しているように見られます。
> k-epsillonで発散する場合、k-omegaにするというのも一方法かと。
> こちらでもよく経験するのですが、MRFSimpleFoamに限らず、
> この方法で凌いでいるケースが多いです。
>
> 2010年4月5日23:51 Sakuma <sakura.maru...@gmail.com>:
> ...
>
> もっと読む ≫- 引用テキストを表示しない -
>
> - 引用テキストを表示 -
E.Mogura
まずは計算できたようで、良かったですね。
自分の説明（k-omega）もkOmegaSSTを言及したつもりでしたが、
説明不足でした。
自分がこれを最初に選択するようになったのは、
OpenFOAMのチュートリアルの中で使っている例が多かったからです。
その後、色々調べた結果、
壁際で乱流が強くなる現象では、k-epsillonよりも
kOmegaSSTの方が良い。。。といったあたりはどこかに書いてありました。
ただ、それでも市販ソフトは、結構、k-epsillonで計算してくれてしまうという
面があって、そこんところのはっきりとした理由はわかりません。
E.Mogura
（話は脱線というか、本スレッドの少し前の投稿に対する返信になります）
回転領域を定義する形状データ（STLファイル）からcellSetを使って、
回転領域（cellZones）を作成する方法について紹介しましたが、
最近、もっと簡単に作成する方法を見つけたので紹介します。
⇒snappyHexMesh作成時に同時に作成可能であることが判りました。
たとえば、
回転領域を定義した形状データ：rotorZone.stl
ファンの形状データ：fan.stl
であるとすると、snappyHexMeshDictにて、
以下定義してメッシュ作成するだけです。
geometry
{
rotorZone.stl
{
type triSurfaceMesh;
name rotorZone;
}
fan.stl
{
type triSurfaceMesh;
name fan;
}
refinementSurfaces
{
rotorZone
{
// Surface-wise min and max refinement level
level (1 1);
faceZone ROTOR;
cellZone ROTOR;
zoneInside true;
}
fan
{
// Surface-wise min and max refinement level
level (1 1);
｝
この方法を使えば、回転領域面にスナップされ、エッジ部分を除いて
ギザギザの部分はほとんど出来ないので安心です。お試しあれ。
なお、このrefinementSurfaceで
STL内部もメッシュ作成できて（zoneInside true;）、
Zoneに名前を付ける（cellZone ROTOR;）
機能があったという点に関しては、マニュアルにはなく、
OF-1.6のtutorials中、「snappyMultiRegionHeater」の中で
使われていたものを参考にしたんですが、
実は、OF-1.5.xでも使用可能であったというオチまでついています。
IMANO Masashi
今野です。
Sakumaさん＞
kEpsilonモデルで全てのrelaxationFactorsを0.1や0.05位まで下げても発散し
ますか？
Sakuma
p 0.05;
U 0.05;
k 0.05;
epsilon 0.05;
omega 0.05;
MRFZonesの回転数を209で，k-epsilonモデルで計算するとTIME=2000（最後）まで発散せずに計算出来ましたが，
TIME=600当たりから残差は振動を繰り返してTIME=2000まで行きました。
乱流モデルのみをk-omegaにするとTIME=65で発散しました。リラクゼーションファクターが0.3の時は8回で発散しまし
たので，少しは計算が進んだと見て良いのでしょか？
乱流モデルがk-omega-sstは，MRFZonesの回転数が209でもその10倍の2090でもリラクゼーションファクターが0.3で
TIME=2000まで問題なく計算出来ました。残差は振動を繰り返すことなく下がっていきます。
ラクゼーションファクターとソルバーの設定は，どのような基準で実施したら良いでしょうか。良く理解出来ていないのが
現状です。
On 4月8日, 午前11:59, IMANO Masashi <masashi.im...@gmail.com> wrote:
> 今野です。
>
> Sakumaさん＞
> kEpsilonモデルで全てのrelaxationFactorsを0.1や0.05位まで下げても発散し
> ますか？
IMANO Masashi
今野です。
> MRFZonesの回転数を209で，k-epsilonモデルで計算するとTIME=2000（最後）まで発散せずに計算出来ましたが，
> TIME=600当たりから残差は振動を繰り返してTIME=2000まで行きました。
乱流計算の場合では、よほど流れ場が定常的かつ単純なものでないと、定常計
算では残差はきれいに下降せず、途中から振動することが多いと思います。
残差が振動する場合には、収束を判定するには、その流れ場において重要な値
や主要な位置でのサンプリング値が収束することを確認することが必要になっ
てきます。
ファン流れですと、流量が重要な値だと思いますので、
OpenFOAM-1.6.x/src/postProcessing/functionObjects/field/fieldValues/controlDict
を参考して、以下をcontrolDictに達して、毎ステップ毎に流量をモニターリング
します。
functions
(
inletFlowRate
{
type faceSource;
functionObjectLibs ("libfieldFunctionObjects.so");
enabled true;
// outputControl outputTime;
outputControl timeStep;
outputInterval 1;
log true;
valueOutput true;
source patch;
sourceName inlet;
operation sum;
fields
(
phi
);
}
);
なお、outputControlをoutputTimeにすると、計算結果を書き出す時だけ、上記の
値が出力されます。一緒の値になる筈ですが、sourceNameをoutletにして、吹
出口側の流量も一緒にモニターすることも一度試してみてください。
これで計算を走らせると、inletFlowRate/初期時間/faceSource.dat といった
ファイルが出来、その中のsum(phi)が流量なので、それをExcelやgnuplot等で
グラフ化すれば、流量の収束状況がわかると思います。
> ラクゼーションファクターとソルバーの設定は，どのような基準で実施したら良いでしょうか。良く理解出来ていないのが
> 現状です。
ソルバーというか、乱流モデルや離散化スキームの選択には、自分が解きたい
流れ場に近く、かつ実験値のような正解がわかっている流れ場をまず最初に解
いてみて、乱流モデルや離散化スキームに関するケーススタディーを行い、流
れ場に適切なモデルを選択する、いわゆるベンチーマークテストを行うことが
必要です。その後、自分が解きたい問題に進んでいきます。
緩和係数は、計算が不安定なら小さくし、安定なら大きくしますが、どの値を
用いたら良いかという基準は、ケースバイケースだと思います。
あと、この解析ではチューブの両端がinletとoutletで、その圧力条件は0なの
ですよね？こういった問題を解いたことが無いのですが、そもそも、速度分布
(動圧分布)があるチューブ両端で、静圧0の条件で解いても良いのでしょうか？
ある程度両端から離れないといけない気がするのですが。
このような境界条件の設定についても、ベンチマークテストを行なわないと、
良し悪しがはっきりしないので、まずは似たような流れ場の正解を再現する
ことから始めることを強くお勧めします。
Sakuma
りがとうございます。早速，今の計算モデルで試してみましたが，メッシュが何と言うか上手く作成
出来ませんでした。
geometry設定を次のようにしています。
FRONT_mm.stl （円管に流体の入口部）
{
type triSurfaceMesh;
name FRONT;
}
BACK_mm.stl （円管の流体の出口部）
{
type triSurfaceMesh;
name BACK;
}
ROTOR_mm.stl （ローター）
{
type triSurfaceMesh;
name ROTOR;
}
STATOR_mm.stl （円管の側面部分）
{
type triSurfaceMesh;
name STATOR;
}
ROTOR_AREA_mm.stl （ローターを囲む回転領域）
{
type triSurfaceMesh;
name ROTOR_AREA;
}
refinementSurfacesの設定を次のようにしています。
ROTOR
{
level (3 4);
}
FRONT
{
level (2 3);
}
BACK
{
level (2 3);
}
STATOR
{
level (2 3);
}
ROTOR_AREA
{
level (2 3);
faceZone ROTOR_AREA;
cellZone ROTOR_AREA;
zoneInside true;
}
refinementRegionsの設定はしませんでした。
このような設定でsnappyHexMeshを行ない，paraFoamでメッシュを確認すると，Mesh Partsには
internalMeshとROTORとROTOR_AREA のパッチしか出来ませんでした。
refinementSurfacesでROTOR_AREAを設定しないと，internalMesh，FRONT，BACK，ROTOR、
STATORの各パッチが出来ます。どこの設定がおかしいのでしょうか。
上手くメッシュ作成が出来たら，cellSetDictのtopoSetSourcesで次の設定をして
zoneToCell
{
name ROTOR_AREA
}
cellSet，setsToZones -noFlipMapを行い，計算を流し以前に教えて頂いた回転領域を
作成する方法と結果を比べたいと思っていますが。。。
E.Mogura
locationInMesh ();
の設定を確認してくれませんか？
ROTORとROTOR_AREA の中間点になっていませんか？
ROTOR_AREA の外側としてみるなど、変更してみてください。
> 上手くメッシュ作成が出来たら，cellSetDictのtopoSetSourcesで次の設定をして
また、この方法でメッシュ作成すると、
cellZonesも同時に作成してくれるので、
cellSetの処理は不要になるはずです。
E.Mogura
追記です
> > 上手くメッシュ作成が出来たら，cellSetDictのtopoSetSourcesで次の設定をして
>
> また、この方法でメッシュ作成すると、
> cellZonesも同時に作成してくれるので、
> cellSetの処理は不要になるはずです。
ROTOR_AREA
{
level (2 3);
faceZone ROTOR;
cellZone ROTOR;
zoneInside true;
}
としておけばよい、ということです。
Sakuma
ご指摘頂きましたlocationInMesh ()の設定ですが，私のモデルでは，回転領域を示すROTOR_AREAが
ファンの形状を示すROTORを取り囲む領域になっており，ROTORには軸に通す穴が空いています。このため，
locationInMesh ()設定値には，この穴の内部を指定していました。つまり，回転領域を示すROTOR_AREAの
内側に入っていました。お教え頂いたようにROTOR_AREAの外側にlocationInMesh ()の設定し，snappyHex
Meshを実施すると上手くメッシュ（internalMesh，FRONT，BACK，ROTOR、STATOR）が作成できました。
この時，polyMeshの下に作成されるboundaryですが次のようになっています。最後にある回転領域
（ROTOR_AREA_Cylinder）ですが，初期はwallとなっています。この部分は削除して6個ある境界を5としても
良いのでしょうか。またはpatchに変更するだけで良いのでしょうか。
6
(
defaultFaces
{
type empty;
nFaces 0;
startFace 1244112;
}
FRONT_Circle
{
type wall;
nFaces 1744;
startFace 1244112;
}
BACK_Circle
{
type wall;
nFaces 1744;
startFace 1245856;
}
ROTOR_ROTOR
{
type wall;
nFaces 70372;
startFace 1247600;
}
STATOR_Circle
{
type wall;
nFaces 18073;
startFace 1317972;
}
ROTOR_AREA_Cylinder
{
type wall;
nFaces 0;
startFace 1336045;
}
)
E.Mogura
ROTOR_AREA_Cylinder
{
type wall;
nFaces 0;
startFace 1336045;
}
このブロックは、（nFaces 0; ）となっているので、
実体がないということです。
ブロックごと削除して構いません。
そして、ご指摘のように、
境界数を６から５に変更する必要があります。
patchとして残しておいても実害はありませんが、
その場合は、境界条件の方にも残しておく必要があります。
以上
Sakuma
メッシュを作成した後にrenumberMeshをすると収束性が改善されると以前にアドバイスを
頂きましたのでやってみるとエラーというかワーニングが出ます。
私が行った手順は，
①blockMeshの実施
②snappyHexMeshの実施。ホルダが1,2,3と3つ出来ます。このsnappyHexMeshでは，E.Mogura さん
に教えていただいたcellZonesを同時に作成する設定をしています。
③setsToZones -noFlipMapの実施
④renumberMesh -latestTimeを実施してログを見ると
Create mesh for time = 3
Mesh size: 424064
Band before renumbering: 412685
Reading volScalarField cellLevel
--> FOAM Warning :
From function polyTopoChange::addMesh(const polyMesh&, const
labelList&,const labelList&, const labelList&,const labelList&)
in file polyTopoChange/polyTopoChange/polyTopoChange.C at line
2280
Cell:355209 centre:(-0.0129054 -0.004542283 -0.005482073) is in
two zones:ROTOR_AREA and addedCells
This is not supported. Continuing with first zone only.
--> FOAM Warning :
From function polyTopoChange::addMesh(const polyMesh&, const
labelList&,const labelList&, const labelList&,const labelList&)
in file polyTopoChange/polyTopoChange/polyTopoChange.C at line
2280
Cell:355210 centre:(0.01290542 -0.004542284 -0.005481997) is in
two zones:ROTOR_AREA and
こんな感じで256021行まで続きます。
③をぜずに④を行うと，下記のように問題なく出来ます。
Create mesh for time = 3
Mesh size: 424064
Band before renumbering: 412685
Reading volScalarField cellLevel
Band after renumbering: 15323
Writing mesh to 4
End.
何か上手い対応方法はあるのでしょうか。
E.Mogura
自分の説明不足だったかもしれませんが、、、
cellSetが不要になるという意味は、
> ③setsToZones -noFlipMapの実施
も不要になるということです。
MRFに必要なファイルはcellZonesというファイルで、
従来は、cellSet⇒setToZonesの手続きでやっていたものを、
snappyHexMeshで直接cellZonesを作成できるということです。
Sakuma
私の不理解でお手数をおかけしました。両方の方法で計算を行えるようになりました。
cellSet，setsToZones -noFlipMapを使った場合より，snappyHexMeshDictの中に
メッシュの設定を入れ込んだ方が，メッシュ作成の効率が良い事が理解できました。
OpenFOAMのチュートリアルでもこの方法を最初から紹介してくれれば，もう少しすん
なりと出来たとおもいますが。
MRFsimpleFoamの使い方もまだまだの状況ですが，なんとか動かせるまで来ました
（そのつもりです）。ただ，今野さんからアドバイスを頂いています境界条件や流量の収
束チェック(モニタリング)はまだ出来ていないため，今度はここをやってみようと思います。
まだまだ先は遠そうですが。
>
> - 引用テキストを表示 -
Sakuma
毎ステップの流入量と流出量を書き出しました。
functions
(
frontFlowRate
{
functionObjectLibs ("libfieldFunctionObjects.so");
enabled true;
// outputControl outputTime;
outputControl timeStep;
outputInterval 1;
log true;
valueOutput true;
source patch;
fields
(
phi
);
}
{
functionObjectLibs ("libfieldFunctionObjects.so");
enabled true;
// outputControl outputTime;
outputControl timeStep;
outputInterval 1;
log true;
valueOutput true;
source patch;
fields
(
phi
);
}
);
backFlowRate/0/faceSource.datのデータをexcelで図化すると，Time=800回を少し越え
る当たりで収束が確認でしました。
ご指摘を頂きましたチューブ両端の圧力条件ですが，今回のモデルは，ファンのP-Q特性
のP=0でのQの値を計算を実際のモデルで計算するための準備で，解析の流れを作れない
かと思ってやっていました。このためそこまで考えていませんでした。お教え頂きありがとう
ごさいます。
functionsですが，英語版(1.6)のユーザーガイドのU-110にチュートリアル参照のような書き
方になっていますが，operationの引数について書かれているものはないでしょうか。
また，functionsで書き出す値（今回であれば流量）を計算途中でモニタ表示（残差や連続性
のように）をさせる事は出来ないでしょうか。
On 4月10日, 午後3:27, IMANO Masashi <im...@arch.t.u-tokyo.ac.jp> wrote:
> 今野です。
>
> > MRFZonesの回転数を209で，k-epsilonモデルで計算するとTIME=2000（最後）まで発散せずに計算出来ましたが，
> > TIME=600当たりから残差は振動を繰り返してTIME=2000まで行きました。
>
> 乱流計算の場合では、よほど流れ場が定常的かつ単純なものでないと、定常計
> 算では残差はきれいに下降せず、途中から振動することが多いと思います。
>
> 残差が振動する場合には、収束を判定するには、その流れ場において重要な値
> や主要な位置でのサンプリング値が収束することを確認することが必要になっ
> てきます。
>
> ファン流れですと、流量が重要な値だと思いますので、
>
IMANO Masashi
今野です．
> 計算事体はTime=1000回まで回しました。計算後にfrontFlowRate /0/faceSource.dat，
> backFlowRate/0/faceSource.datのデータをexcelで図化すると，Time=800回を少し越え
> る当たりで収束が確認でしました。
残差が振動するkEpsilonでも流量は収束しましたか?
> functionsですが，英語版(1.6)のユーザーガイドのU-110にチュートリアル参照のような書き
> 方になっていますが，operationの引数について書かれているものはないでしょうか。
おそらく無いと思います．
こういう場合には，ソースを見るのが正攻法なのでしょうが，それも面倒な時には，
operationの引数を空にする等してわざと間違えた上でソルバーを実行をすると，
以下のように候補が出てくることが多いです．
--> FOAM FATAL IO ERROR:
is not in enumeration:
5
(
areaAverage
weightedAverage
areaIntegrate
none
sum
)
> また，functionsで書き出す値（今回であれば流量）を計算途中でモニタ表示（残差や連続性
> のように）をさせる事は出来ないでしょうか。
log true;
になっていると
faceSource outletFlowRate output:
sum(outlet) for phi = 0.000324735
のように標準出力に表示されますが，出ませんか?
Sakuma
乱流モデルをkomegaSST, kepsilon，komegaの3モデルで試してみました。
1000回まで計算をしています。
komegaSSTで
①収束 relaxationFactorsがp=0.3 U=0.3 k=0.3 omega=0.3
②発散 relaxationFactorsがp=0.3 U=0.7 k=0.7 omega=0.7
③収束 relaxationFactorsがp=0.3 U=0.5 k=0.5 omega=0.5
kepsilonで
①収束 relaxationFactorsがp=0.1 U=0.1 k=0.1 epsilon=0.1
②発散 relaxationFactorsがp=0.2 U=0.2 k=0.2 epsilon=0.2
komegaで
①収束 relaxationFactorsがp=0.1 U=0.1 k=0.1 epsilon=0.1
②発散 relaxationFactorsがp=0.2 U=0.2 k=0.2 epsilon=0.2
solversは，
0.1;
k：solver PBiCG; preconditioner DILU; tolerance 1e-05; relTol
0.1;
epsilon：solver PBiCG; preconditioner DILU; tolerance 1e-05;
relTol 0.1;
omega： solver PBiCG; preconditioner DILU; tolerance 1e-05; relTo
です。
やはりkomegaSSTが良く解いてくれるように見えますが。
また，結果の書き出しですが，計算実行のログを取るため次のようにしています。
pyFoamPlotRunner.py MRFSimpleFoam > log.MRFSimpleFoam
多分この影響で標準出力に表示されていないと思います。ログも取りながら，
標準出力に書き出す方法はあるのでしょうか。
IMANO Masashi
今野です．
なるほど，この問題ではKOmegaSSTのほうが安定のようですね．
モデル間で流量はどの位違うのでしょうか？
実験値がある問題で是非テストしてみてください．
その時には，RNGkEpsilon等他のRASモデルも検討されると良いと思います．
> また，結果の書き出しですが，計算実行のログを取るため次のようにしています。
> pyFoamPlotRunner.py MRFSimpleFoam > log.MRFSimpleFoam
> 多分この影響で標準出力に表示されていないと思います。ログも取りながら，
> 標準出力に書き出す方法はあるのでしょうか。
1. teeコマンドを使う
pyFoamPlotRunner.py MRFSimpleFoam | tee log.MRFSimpleFoam
2. 計算をバックグラウンドで走らせ，tailコマンドでログを監視する
pyFoamPlotRunner.py MRFSimpleFoam > log.MRFSimpleFoam &
tail -f log.MRFSimpleFoam
もし，流量も逐次プロットするなら，
set style data l
plot "frontFlowRate/0/faceSource.dat" using 1:3 title "Front"\
,"backFlowRate/0/faceSource.dat" using 1:3 title "Back"
reread
という内容のファイルflowRate.gpを作成して，以下のようにしてgnuplotでプ
ロットすれば良いのではないでしょうか？
pyFoamPlotRunner.py MRFSimpleFoam > log.MRFSimpleFoam &
gnuplot flowRate.gp &
tail -f log.MRFSimpleFoam
E.Mogura
Sakuma
今回のモデルでの流量ですが，次のようになっています。（ファンが円管内で単純に回転数2000rpm
している単純なものです。ベルマウスなどはありません）
最終ステップtime=1000で
komegaSST
①relaxationFactorsがp=0.3 U=0.3 k=0.3 omega=0.3 Φ= 0.004088496
③relaxationFactorsがp=0.3 U=0.5 k=0.5 omega=0.5 Φ=0.004048871
kepsilon
①relaxationFactorsがp=0.1 U=0.1 k=0.1 epsilon=0.1 Φ=0.003901237 ただグラフを書
くと収束が
まだのようですので，時間がある時にもう少し先まで計算をしようかと思っています。
komegaで
①relaxationFactorsがp=0.1 U=0.1 k=0.1 epsilon=0.1 Φ=0.004046178
というような感じになっています。
実際に実験値があるものとの比較は，次のステップとしてやって行きたいと思っています。
> > 詳細については、http://groups.google.com/group/openfoam?hl=jaからこのグループにアクセスしてください。
> --
> このメールは Google グループのグループ「OpenFOAM」の登録者に送られています。
> このグループに投稿するには、open...@googlegroups.com にメールを送信してください。
> このグループから退会するには、openfoam+u...@googlegroups.com にメールを送信してください。
>
> - 引用テキストを表示 -
Sakuma
MRFsimpleFoamで別モデルの計算をしています。
流入側，流出側の流量をモニターするために，controlDictに以前，
今野さんから教えて頂いた関数を追加しています。
functions
(
inlet_FlowRate
type faceSource;
functionObjectLibs ("libfieldFunctionObjects.so");
enabled true;
// outputControl outputTime;
outputControl timeStep;
valueOutput true;
source patch;
fields
(
phi
);
}
type faceSource;
functionObjectLibs ("libfieldFunctionObjects.so");
enabled true;
// outputControl outputTime;
outputControl timeStep;
valueOutput true;
source patch;
fields
(
phi
);
}
);
# Source : patch inlet_Mesh
# Faces : 486
# Time sum(magSf) sum(phi)
5 0 0.01336154
10 0 0.006345506
15 0 0.0002748995
20 0 -0.01378661
990 0 -0.02654588
995 0 -0.02676494
1000 0 -0.02655028
Source : patch outlet_Mesh
# Faces : 450
# Time sum(magSf) sum(phi)
5 0.02539325 -0.004955362
10 0.02539325 -0.0111422
15 0.02539325 0.001256138
20 0.02539325 0.01097203
90 0.02539325 0.02652413
995 0.02539325 0.02655748
1000 0.02539325 0.02663415
ここで良く分からないのですが，inlet_Meshのsum(magSf)がゼロになっています。
outlet_Meshのsum(magSf)は0.02539325となっていますが，流入側，流出側は
一辺が0.24mの正方形のため，面積では0.0576m2となります。このsum(magSf)
はどのような値を書き出しているのでしょうか。
また，sum(phi)は流量の単位ですが，m3/sでよいのでしょうか。
このあたりに詳しい方がおられましたら，よろしくお願いいたします。