blockMeshとsnappyHexMeshで二次元の円柱後方の流れ場モデルを作成する方法~その8~controlDict

ヒョウ

今までにメッシュを作成し初期条件と境界条件を設定しました。計算実行までの実施項目としてあとは”system”ディレクトリ内のいくつかのディクショナリファイルを編集します。ここではcontrolDictについて説明します。icoFoamの選択,⊿tの設定,出力頻度,エンドタイムなどの設定を行います。計算が発散しないように数値を入力すれば最低限よいのかもしれませんが,作成したモデルに対して最適な設定方法に対する考え方を解説します。

controlDictの中身

以下の内容です。コピーして利用する場合は,赤字は説明文なのでコピーしないか,コピー後に削除してください。
*****************************************
application icoFoam;     //icoFoamでの計算を想定しています。

startFrom startTime;    //後述しますがstatTimeをlatestTimeと書き換えてもokです。

startTime 0;    //startFormをstartTimeとした場合に参照されます。

stopAt endTime;    //計算終了のタイムステップをendTimeとします。

endTime 400;    //計算終了のタイムステップを入力します。

deltaT 0.05;    //1ステップ当たりの⊿tを入力します。

writeControl timeStep;    //計算結果の自動保存のタイミング設定です。

writeInterval 200;    //タイムステップ200回ごとに自動保存として
                             //そのタイムステップでの計算結果を保存します。

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable yes;
*****************************************

endTime

計算をタイムステップで何秒まで実行すればよいか?答えは「前提や目的による。」わけですが,今回対象としている円柱後方の流れ場モデルでは,初期条件(速度分布ゼロ)から流れが十分発達した際の流れ場をシミュレーションにより可視化したいので,流れ場が十分に発達すると思われるタイムステップまで計算します。

滞在時間θから設定する

どれだけのタイムステップまで計算したらモデル内の流れ場が十分発達するか?を検討初期の段階で考えるならば,入口面流速とモデルの長さから求まる滞在時間 θ を使って2θとしてみます。
θ = L / U  [s]

θ; おおよその見込みとしてのモデル内滞在時間 [s]
L; 流路方向のモデル長さ [m]
U; 入口面流速 [m/s]

ここでは流路方向のモデル長さ L = 1.05 [m], 入口面流速 U = 0.005 [m/s]ですので
θ = 1.05 / 0.005 = 210 [s]
であり,θ ≈ 200 [s]とします

したがってendTimeは 2θ = 400 [s]とします。

2θとすることに深い根拠はありません。また,ここで言っているθは厳密な滞在時間ではなく,おおよその見込みとしての滞在時間ということになります。実際の流れはプラグフローではないですし,流体要素の滞在時間には分布があります。

deltaT

⊿tの設定についてです。与える数値が大きすぎると発散しますし,小さすぎるといつまでたっても計算が終わりません。クーラン数が1より小さくなるように⊿tを与えます。

クーラン数 Cr

クーラン数は「1タイムステップでの流体の移動距離」と「1メッシュの長さの比」で,ここでは次式で与えます。
Cr = ⊿t * U / x [-]

Cr; クーラン数 [-]
⊿t; タイムステップ [s]
U; 入口面流速 [m/s]
x; 1メッシュの長さ [m]

モデル内では局所的に流速が大きくなる領域がありますので,Cr = 0.1 とすることを考えます。U = 0.005 [m/s](0/Uにて設定), x = 0.003 [m](blockMeshにて設定)であり,
⊿t = Cr * x / U = 0.1 * 0.003 / 0.005 = 0.06 [s]
であり,⊿t ≈ 0.05 [s]とします。

以上より,deltaTは0.05 [s]です。

Cr > 1 となるよう⊿tにはするべきではない

隣接するメッシュ間での速度の移送を解くことでモデル全体としての速度場・圧力場を求めます。差分計算における⊿tが大きいために,1タイムステップでの速度の移送が1メッシュの距離を越えてしまう場所があると計算の精度が低下するか発散します。計算の初期での設定ではクーラン数は0.1程度になるようにすればよいと考えます。
VOF法で混相系の計算を行う場合などではもっと小さいクーラン数が求められることがあります。

writeInterval

startTime 0秒からendTime 400秒まで,⊿t 0.05秒で計算をするよう設定しました。endTimeまで計算を実施する中で,途中経過を保存する設定がwriteIntervalです。endTimeまでの計算には,実際の時間で半日から1日かかる場合もあり,場合によっては3日間程度かかる計算もあります。時間(という資源)がかかるので計算の途中経過を残すことで,計算の記録を振り返れるようにします。

途中経過を保存する目的を考える

途中経過を残す目的はいくつかあると思います。
①途中経過を確認することで,計算途中で計算に設定ミスがないか確認する。
②アニメーション(動画)作成のために,途中経過を出力する。
③長期間の計算のため,途中での予期せぬ計算停止に備えて定期的に保存する。
などあると考えます。

ここでのwriteInterval設定の考え方

今回の計算は私が使っているようなロースペックのノートPCでも半日程度で完了します。半日というのも間違った設定で計算し続け,浪費してよい時間ではないでしょう。また,計算が途中で停止することも心配なので,上記の①と③の目的でwriteIntervalを考えます。あまりたくさん途中経過を保存してもディスク容量を浪費してしまいますのでendTimeのデータを含めて10~20個データが保存できていればよいでしょう。
ここではwriteIntervalは400とします。entTime 400秒まで⊿t 0.05秒で計算する場合,タイムステップは8000ステップです。データを20個保存するには400ステップごとに保存すればよいことになります。

startFoamについて

計算をスタートする際の開始時刻の設定です。初期条件を0/Uと0/pに保存していますが,初期のタイムステップは0です。一度停止した計算を途中のタイムステップから開始したり,計算完了したデータをもう少し進めたい場合にstartFoamやstartTimeを設定します。

startFoam latestTime; とした場合

controlDictのstartFoamをstartTime;からlatestTime;に変更することで,wtireIntervalで保存された各タイムステップのディレクトリ中の最新のタイムステップから計算を再スタートできます。

startFoam startTime; で途中から計算する場合

何らかの理由で途中のタイムステップを指定して計算を再スタートしたい場合は,startFoamをstartTimeのままとし,次の行に書かれているstartTimeに指定したいタイムステップを入力します。ただし,指定したタイムステップのデータが保存されている必要があります。

計算の実行まで残りの説明としてmeshQualityDict, fvShemes, fvSolutionがありますがmeshQualityDictは例題からコピーして持ってくればよく,fvShemes, fvSolutionについては数値流体力学の本題であり,説明すると1年くらいかかると思いますので,こちらも例題からコピーして利用するという説明にとどめます。次回これらの説明をして,その後,openFoamのディレクトリ構造を説明したのちにicoFoamの実行となります。

スポンサーリンク
プロフィール
この記事を書いた人

openfoamを使って流体力学のシミュレーションをやっています。openfoamはまだまだ始めたばかりの人にとってわかりにくさがあると思います。そういったわかりにくさを説明するtipsを書いていきながら,他のこともいろいろやっていきたいと思ます。

ヒョウをフォローする
ヒョウ 流体力学
ヒョウをフォローする
mister Big5  サファリ系サラリーマンblog

コメント

タイトルとURLをコピーしました