Comments
Description
Transcript
天井駆動のキャビティー流れ
OpenFOAM ユーザーガイド 1.5 系
天井駆動のキャビティー流れ
出典：OFwikiJa
■天井駆動のキャビティー流れ
二次元正方形領域の等温非圧縮性流れに関して、
プリプロセス(前処理)
計算
ポストプロセス(後処理)する方法
を解説していきます。
図 1 に示すのは正方形のすべての境界が壁面境界であるジオメトリです。
上部の壁面境界は x 軸方向に 1m/s の速度で動き、他の 3 つの壁面境界は静止です。
まずは層流を仮定し、層流等温非圧縮性流れのためのicoFoam ソルバー(ソルバー：計算ソフト)を使用し
均一メッシュ上で解きます。また、メッシュの解像度の増加や壁方向への勾配の影響を調べます。これ
によりレイノルズ数を増加させ、turbFoam ソルバーを乱流、等温、非圧縮性流れに使用します。
図 1. 天井駆動キャビティーのジオメトリ
■１-1. 前処理
ケースは、OpenFOAM でケースファイルを編集することで設定します。ケースファイルは emacs や vi、
gedit、kate、nedit などのテキストエディタで作成・編集します。それは、入出力が初心者でもわかりやす
い辞書(Dict)形式が使われているからです。
※以前のヴァージョン(1.4 系)では FoamX という GUI(GUI：グラフィカルユーザーインターフェース。コンピュー
タグラフィックスとポインティングデバイスを用いて直感的な操作が可能なユーザーインターフェース）ケースエデ
ィタがありましたが、利用者がエディタによって編集できるファイルを好み、メンテナンスを怠る可能性
があったため、1.5 系ではなくなりました。
解析ケースはメッシュ、物理量、物性、コントロールパラメータなどの要素を含んでいます。しかし多くの
CFD ソフトが 1 つのファイルにこれらのデータを格納するのに対し、OpenFOAM は一連のファイルセット
と して解析ケースディレ ク トリに 格納する仕組みとなっ てい ます。わかりやすい よう に、以下に
OpenFOAM のケースのファイル構造を示します。
◆OpenFOAM のケースのファイル構造
◇system ディレクトリ
解析の手順そのもののためのパラメータの設定に関するディレクトリです。最低でも以下の 3 つのフ
ァイルを含んでいます。
・controlDict
パラメータが開始/終了する時間や時間ステップ及びデータのアウトプットのためのパラメータを含ん
でいるようにコントロールを実行します。
・fvSchemes
実行時に選択される解析に使われるスキームを記述しています。
・fvSolution
方程式のソルバー、許容誤差及びその他のアルゴリズム制御を実行のために設定します。
◇constant ディレクトリ
サブディレクトリの polyMesh 内のケースメッシュと物理的性質を定めるファイルの完全な説明と関連
するアプリケーション、例えば tranceportProperties を含みます。
◇time ディレクトリ
特定領域用にデータの個別のファイルを持っているディレクトリです。
データは、問題を定義するためにユーザーが指定する初期値と境界条件、または書き込まれた
OpenFOAM のファイルの結果として存在します。OpenFOAM のフィールドは、定常状態の問題のよう
に厳密に解く必要がない場合でも、常に初期化する必要があります。
各タイムディレクトリの名称は、データが書き込まれた時点のシミュレーションが行われた時間に基
づいています。私達は通常、時間$t=0$でシミュレーションを開始し、初期の条件は指定された名前
のフォーマットに依存して 0、または 0.000000e+00 と名付けられたディレクトリに収納されます。
例えば今回の天井駆動のキャビティー流れで用いる cavity の場合、速度場の U と圧力場の p、それ
ぞれファイル 0/U と 0/p から初期化されます。
1-1.1. メッシュの生成
これから解析の編集・実行していきます。まずは解析対象のディレクトリに移動するために、ターミナル
に以下のコマンドを入力します。
cd $FOAM_RUN/tutorials/icoFoam/cavity
OpenFOAM は、常に 3 次元デカルト座標系(デカルト座標系は直交座標。この場合は x、y に z を加え
た 3 次元の座標系のことを指す)で動くので、全てのジオメトリを 3 次元で生成することになります。
OpenFOAM はデフォルトの設定において問題を 3 次元で解きますが、2 次元を解く場合は、解決が必
要でない(第 3)次元方向に垂直な境界に、特別な『empty』という境界条件を指定します。
x-y 平面上の一辺の長さ d=0.1m の正方形から成る cavity の領域に、まずは 20×20 セルの均一なメッ
シュを設定します(図 3 参照)。
図 3. キャビティーのメッシュのブロック構造
OpenFOAM で提供されるメッシュ・ジェネレーターblockMesh は、constant/polyMesh ディレクトリに格納
されている blockMeshDict で指定された記述よりメッシュを生成します。このケースの blockMeshDict は
以下の通りです(要確認)。
ファイルの最初はバナー(1-7 行)形式のヘッダー情報で、FoamFile サブディレクトリはファイル情報を含
み、({…})で囲まれています。
※以後の説明においては、簡便化とスペースの都合上、バナーと FoamFile サブディレクトリを含むファイルヘッダ
ーはケースファイルの引用の際に省きます。
まずファイルは初めにブロックの頂点の座標 vertices を指定します。それに続き、頂点名とセル番号か
ら blocks(ここでは１つのみ)を定義します。そして最後に境界パッチ(patches)を設定します。
ケースディレクトリ内からターミナル上に以下のコマンドを入力すると、blockMeshDict 上で blockMesh
が実行され、メッシュが生成されます。
blockMesh
blockMesh の実行状況はターミナルウインドウに表示されます。
blockMesh ファイルに誤りがあった場合にはエラーメッセージが表示され、ファイルのどの行に問題が
あるかを教えてくれますが、現段階においてエラーメッセージが出ることはないでしょう。
1-1.2. 境界条件と初期条件
メッシュの生成が完了すると、物理的条件の初期状態を確認することが出来ます。
このケースは開始時刻が t=0 に設定されているので、解析領域の初期状態のデータは cavity ディレク
トリの 0 というサブディレクトリに格納されています。
0 ディレクトリ内には圧力(p)と速度(U)の 2 つのファイルがあり、それらの初期値と境界条件を設定する
必要があります。
P のファイルを参照に、説明します。
物理的条件のデータファイルには、3 つの主要な項目があります。
◇dimensions
物理量の次元を定義します。
◇internalField
内部領域の値。内部の物理量は単一の値で記述すれば一様となり、一様でない場合は全ての値を
指定する必要があります。
◇boundaryField
境界領域。境界面の物理量は境界条件と境界パッチに与えるデータを記述します。
この cavity 流れの解析において境界は壁面のみですが、2 つのパッチが用意されています。
(1) キャビティの固定された側面と底面用の fixedWalls
(2) キャビティの駆動天井用の movingWall
どちらも p の値が zeroGradient ですが、これは圧力の勾配が 0 だということです。frontAndBack は 2
次元の問題の場合の表裏の面を示しています。本ケースでは empty となっています。
このケースで最もよく目にするものですが、物理量の初期条件が uniform(一様)となっています。ここで
は動圧のみの非圧縮ケースであるため、絶対値は解析と関係ないため便宜上 uniform0 としています。
0/U の速度ファイルにおいても同様に、dimensions は速度であり、内部の初期条件はベクトル量で 3 成
分とも 0 であり、uniform(0 0 0)となっています。
速度の境界条件は frontAndBack パッチと同じです。fixedWalls に関しては、滑りが無いため uniform(0
0 0)となります。上面は 1m/s で移動するので uniform(1 0 0)で固定値を設定します。
1-1.3. 物理量
ケースの物理量の名前には「…Properties」という語尾が与えられ、ディレクトリに保存されます。
icoForm ケースでは、transportProperties 辞書に保存される動粘性係数を指定します。
transportProperties 辞書を開き、動粘性係数が正確にセットされていることを確認します。動粘性係数
は nu(ν)で表されます。
このケースではレイノルズ数を 10 として計算します。レイノルズ数とは流体に作用する慣性力と粘性力
の比を示します。
d と U はそれぞれ特性長さと速度を表します。ここで、d=0.1m、U=1m/s、Re=10 として計算すると、レイノ
ルズ数はν=0.01m2/s となります。動粘性係数の適切な設定を以下に示します。
1-1.4 制御
計算時間の制御を行います(system ディレクトリ内 controlDict 参照)。
◇スタート時間の設定
OpenFoam では柔軟性の高い時間制御を提供しています。
まずは、時間 t=0 から実行を開始します。つまり、OpenFoam は 0 というディレクトリからフィールドデー
タを読み込む必要があります。
startFrom に「startTime」と設定したら、startTime は 0 に指定します。
◇終了時間の設定
終了時間(endTime)には、流れがキャビティー周りを循環している定常解に達することを目標とします。
概して、流体は層流で定常状態に到達するために領域を 10 回通り抜けなければなりませんが、この
cavity ケースでは出入口がないので解析領域を通り抜けしません。
代わりに、上部の壁面境界がキャビティーを 10 回移動する時間(＝1秒)を終了時間としてセットします。
しかし、実際は後の知見により移動時間は 0.5 秒で十分だとわかるので、こちらの値を採用しましょう。
この終了時間を指定するためには stopAt キーワードを「endTime」にして、endTime に 0.5 を設定しま
す。
◇時間ステップの設定
deltaT によって設定されます。icoFoam を動かす時、時間の精度と安定性を達成するため、1 未満のク
ーラン数(流体解析の数値解析に必要となる無次元数)が必要です。クーラン数を定義します。
Δt は時間ステップ、¦U¦はセルを通る流速の大きさ、そしてδx は流速方向のセルサイズです。
流速が領域内で変化しても、Co<1 が成立しなければなりません。だから、大きな流速と小さなセルサイ
ズの組み合わせによる最大の Co を元にδt を決定します。
今回セルサイズは解析領域内全域で固定されているので最大 Co は上部の流速 1m/s に近似します。
従って、領域内で 1 以下のクーラン数を達成するために、時間ステップ deltaT を次のように設定しま
す。
シミュレーションが進行する時、後でポスト処理パッケージで見ることができるように、ある一定の時
間間隔での結果の書き出しをもとめる際、 writeControl キーワードは結果が書かれる時間を決める
ためのいくつかのオプションを提示します。
timeStep オプションは、結果が
回時間ステップ毎に結果を書き出すということを意味し、その時の
値は writeInterval キーワードで指定されます。
0.1, 0.2, ･･･, 0.5 秒で結果を書きたいとしましょう。従って、0.005 秒の時間ステップなので、時間ステップ
20 回毎に結果を出力する必要があります。よって writeInterval に 20 を設定します。
icoFoam ソルバーでは、U や p の各項目ごとに結果を時間ディレクトリに書き込みます。
1-1.5. 離散化と線形ソルバーの設定
ユーザーは fvSchemes 辞書（system ディレクトリ）内で有限体積離散化法を選択するかどうか指定しま
す。線形方程式ソルバーとトレランスおよび他のアルゴリズムコントロールの指定は fvSolution 辞書内
に作られています。ユーザーは自由にこれらの辞書を見ることができます。
キャビティーのような閉じた非圧縮系では、圧力は相対的であり、重要なのは絶対値ではなく圧力範囲
です。 このような場合では、ソルバーはセル pRefCell に pRefValue による参照レベルをセットします。
この例では、両方が 0 に設定されます。しかし、これらの値のどちらかを変えると、速度と相対圧力では
なく、絶対圧力が変化します。
※有限体積法…FVM とも言う。多くの商用流体解析コードに標準的な離散化解析手法として用いられ
る。解析領域をセルに分割し、セルの格子点を中心とする領域(コントロールボリューム)あるいは検査
領域 De を定義するもの。
※離散化…音などの連続した信号を時間や振幅で区切ること。
■1-2. メッシュの確認
解析を実行する前に、正しくメッシュが出来ているか確認しましょう。確認には OpenFoam が提供する
paraFoam というポスト処理ソフトをします。
解析ケースのディレクトリ上で、ターミナルに次のコマンドを打ち込み、paraFoam を起動します。
paraFoam
あるいは、オプションに ‒case をつけることで他のディレクトリからでも起動することが出来ます。
paraFoam ‒case $FOAM_RUN/tutorials/icoFoam/cavity
図 4 に ParaView のウインドウが開いた状態を示します。
図 4 paraFoam でのメッシュ表示
PipelineBrowser にて ParaView が cavity.OpenFOAM、キャビティケースのモデュールを開いていること
が確認できます。
Apply ボタンをクリックする前に Region Status とパネルから表示する要素を選択する必要があります。
解析ケースが単純なので Region Status パネルのチェックボックスで全てのデータを選択することが簡
単です。パネル内の全要素を自動的にチェックすることが出来ます。
ParaView でジオメトリを読込むためには Apply ボタンをクリックします。Display タブを開き選択したモジ
ュールの表示形式を調整します。
ここで色の設定をします。まずは Color というブラウザにおいて、Color by Solid Color を設定し、次にそ
の下の Set Solid Color でメッシュの色を選択します(例えば背景が白なら黒など)。背景色は、トップメニ
ューの edit から View Settings...を選択して設定することができます。最後に、Color ブラウザの下、Style
ブラウザにおいて Representation メニューから Wireframe を選択します。
paraFoam における視点操作を試してみましょう。Edit メニューの View Settings の General パネルで Use
Parallel Projection を選択してください。マウスを左クリックしながら、あるいはローラーを押さえながら動
かすと、軸の操作が出来ます。
一通りの操作を終えたら、ParaView を閉じます。
■1-3. アプリケーションの実行
あらゆる UNIX/Linux の実行ファイル同様、OpenFOAM アプリケーションは 2 つの方法で実行すること
が出来ます。
1 つ目はフォアグラウンドでのプロセスで、コマンドプロンプトを与えるのにシェルが命令終了まで待つ
ものです。2 つ目はバックグラウンドプロセスで、シェルがさらなる命令を受け入れるのに命令完了の必
要がないものです。
ここではフォアグラウンドにて icoFoam を起動します。cavity ケースディレクトリ内で、ターミナルに
icoFoam
と入力すると、icoFoam が実行されます。あるいはオプションに ‒case をつけることで、他のディレクト
リからでも起動可能です。
icoFoam ‒case $FOAM_RUN/tutorials/icoFoam/cavity
ジョブの進行はターミナルウインドウに表示されます。現在の時刻、最大クーラン数、すべてのフィール
ドの初期値と最終的結果を表示します。
■1-4. ポスト・プロセッシング
結果が時間ディレクトリに書かれるので、paraFoam を起動します。
cavity.foam ケースが黄色でハイライトされているか/それと並んだ eye ボタンにグラフィックス可能であ
ることを示すスイッチは入っているかを確認してください。
今回は最初に必須の実行時間として 0.5 秒分のデータを読み込まなければなりません。
ケースの実行中に ParaView を開いていても、時間ディレクトリの出力データは ParaView に自動的にロ
ードされません。ロードするには、Properties 画面で Update GUI を選択し、Apply ボタンを押します。
1-4.1. コンタープロット
図 5 キャビティーケースでの圧力等圧線のディスプレイ
圧力を見るには Display パネルを開き、選択したモジュールの表示形式を調整します。
圧力分布を見るには図 5 に示すように StylePanel の Representation メニューを surface にして、
ColorPanel の Set Color by で「・p」を選択します。
そして Rescale to Data Range ボタンをクリックし、メニューバーの下のツールバーにある Current
Time(または VCR Control)を t=0.5 にして、t=0.5 における解析結果を表示します。
Active Variable Control ツールバーの Toggle Color Legend Visibility ボタンをクリックするか、View メニ
ューから Show Color Legend を選択することでカラーバーの表示が可能となります。
Active Variable Control toolbar か Display ウインドウの Color panel にある Edit Color Map button をク
リックすると、フォントの大きさや種類、スケールの番号付けの形式などの設定変更が出来ます。
カラーバーはドラッグ&ドロップでイメージウインドウに置くことも出来ます。
イメージを回転させることによって全ての表面に圧力分布で着色されていることが確認できます。
図 6 cavity ケースでの圧力分布とカラーバー
正しい断面のコンターを得るために、slice フィルターを用いてジオメトリを slice します。
断面の中心座標は(0.05, 0.05, 0.005)、基準点は(0, 0, 1)とします。断面を作成後、上部のメニューバー
の Filter メニューから Contour を選択することによりコンターを描画することができます。
図 7 slice 後の状態
1-4.2. ベクトルプロット
流速ベクトルを描画する前に、作成した断面やコンターなどのモジュールは不要なので取り除きましょ
う。Pipeline Browser でそれらのモジュールを選択して Properties Panel の Delete をクリックして削除す
るか、Pipeline Browser で eye ボタンをクリックしてモジュールを非表示にすればよいでしょう。
各格子の中心におけるベクトルグラフを作成します。
ベクトルは、デフォルトによりセルの頂点上に作成されますが、今回はセルの中心上にプロットデータを
作成します。Pipeline Browser で cavity.OpenFOAM のモジュールが選択されているか確認し、Filter メ
ニューから Cell Centers を選択します。
すると Centers が強調表示されるので、Filter メニューから Glyph を選択します。図 8 のように Properties
パネルが表示されます。パネルを見ると、vectors メニューにおいて、ベクトル場は速度のみなので速
度場 U が選択されています。Scale Mode は速度の Vector Magnitude が初期値として設定されていま
すが off にします。そして Set Scale Factor の Edit に 0.005 と入力して全体の速度を見ます。
図 8 Glyph フィルターのパラメータパネル
Apply を押すと、単色(例えば白)のベクトルが表示されます。
通常は Display パネルで Color by U を選択して大きさに合わせた色付けをします。Edit Color Map を
Show Color Legend に設定し、速度の凡例(カラーバー)を表示させます(図 9 参照)。
背景色は View Settings の General Panel で白に設定されています。
カラーバーのフォントの変更等は前述にならいます。
図 9 cavity ケースの速度ベクトル
1-4.3. 流線プロット
ParaView でポストプロセスを続ける前に、前述のベクトルプロットは不要なので除去します。
Pipeline Browser で cavity.foam のモジュールを選択した状態で、Filter メニューから Stream Tracer を
選択し、Apply をクリックします。すると図 10 のように Properties ウインドウが表示されます。
Seed Type はジオメトリの中心を垂直に通って Line Source に従うように指定します(例えば、(0.05, 0.0,
0,005)から(0.05, 1.0, 0.005)まで)。
今回は Point Resolution を 21、Max.Propagation を Length で 0.5 に、Init.Step Len を Cell Length で 0.01
に、Integration Direction を BOTH という設定とします。
また、Runge-Kutta 2 Integrator Type はデフォルトパラメータを使います(変更しない)。
Apply をクリックすると、Tracer(トレーサー)が生成されます。
図 10 Stream Tracer のパラメータパネル
図 11 生成したトレーサー
Filter メニューから Tube を選択すると、高品質の流線図を作ることが出来ますが、画面が非常に重くな
るので焦らず操作を行いましょう。おそらくモジュールを動かすだけで数秒かかります。
このときの設定は、Num.sides を 20、Radius を 0.003、Radius factor を 10 とします。Apply をクリックする
と図 12 のような流線図が出来上がります。
図 12 cavity ケースの流線
■1-5. メッシュの解像度の増加
さて、流線図まで描画し終えたなら、次はメッシュの解像度を増やします。
問題の初期条件として使うため、粗いメッシュでの結果を、細かいメッシュ上に写像します。そして、細
かいメッシュの解を粗いメッシュの解と比較します。
手順を説明していきます。
1-5.1. 既存ケースを用いた新しいケースの作成
cavity ケースを複製・修正することで新しい解析ケース cavityFine を作成します。
icoFoam に入り、cavityFine というケースを作ります。
cd $FOAM_RUN/tutorials/icoFoam
mkdir cavityFine
しかしここで「 cavityFine を作成出来ません：File exists」というメッセージが出ます。つまり既にこのデ
ィレクトリは存在しているという事になります。
ls
とコマンドを打つと、icoFoam ディレクトリの内容がリスト表示され、その中に cavityFine ケースがあるこ
とを確認できますので、そのまま作業に戻ります。
基本となる解析ケース cavity の内容を解析ケース cavityFine にコピーし、cavityFine ケースに移動しま
す。
cp ‒r cavity/constant cavityFine
cp ‒r cavity/system cavityFine
cd cavityFine
※コマンド解説※
cd…カレントディレクトリ(ユーザーが現在作業を行っているディレクトリ)の移動を行う。
mkdir…指定したディレクトリ名でディレクトリの作成を行う。
cp -r…ファイルのみ(例えばテキストファイルなど)をコピーする。cp だけだとディレクトリごとコピーする
ことになる。
1-5.2. 細かいメッシュの作成
cavityFine 内において、blockMesh を使って計算格子数を増やします。
そのためには、blockMeshDict ファイルをエディターで開いて blocks に関する記述の修正を行います。
ここで、メッシュのブロックの頂点とブロックについて説明します。
図 13 に示すのは、cavity ケースの blockMeshDict における vertices と blocks です。
図 13 cavity ケースの blockMeshDict
メッシュのブロックの頂点は vertices と名付けられた標準のリストとして与えられます。
そしてブロックの定義は blocks と名付けられたリストに含まれています。各ブロックの定義は図 13 で示
された順序を持つ頂点ラベルのリストから成る複合入力です。ベクトルが各方向に必要なセルの数、タ
イプ、及び各方向のセル拡大比のリストを与えます。
ブロックの定義を示します。
◇hex (0 1 2 3 4 5 6 7) → vertex numbers
◇(20 20 1) → numbers of cells in each direction
◇simpleGrading (1 1 1) → cell expansion ratios
次に、それぞれのブロックの定義を示します。
◆Vertex numbering…ブロックの形状識別子。いつもブロックが六面体なので hex となっています。
◆Number of cells…そのブロックの x1 と x2 と x3 それぞれの方向のセルの数を与えます。
◆Cell expansion ratios…ブロックにおける各方向へのセルの拡大比を与えます。
以上のように定義がなされているので、今回は hex が最初の頂点リストで、各方向の計算格子の番号
リストがあるという事を知っておけば良いでしょう。
そして、先の cavity ケースでは(20 20 1)となっていますが、これを(41 41 1)に変えて上書き保存します。
ここで blockMesh を run させる(走らせる)ことで、新しく、より細かいメッシュを生成することが出来ます。
1-5.3. 粗いメッシュの結果を細かなメッシュにマッピングする
mapFields ユーティリティは、他のジオメトリに対応するフィールドの上へジオメトリに関した一つ以上の
フィールドをマッピング(例：A と B を対応付けること)します。
今回の例では、入力フィールドと求める結果のフィールド両方のジオメトリ・境界の種類・境界条件が同
一なので、フィールドは首尾一貫していると考えられます。
mapFields maps のフィールドデータは、目的のケースの cotrolDict 内の startFrom/stsrtTime で指定さ
れる時間ディレクトリから読まれます。
今回は cavityFine ケースの細かいメッシュ上に、cavity ケースの粗いメッシュの最終結果をマッピングし
ます(この最終結果は cavity の 0.5 というディレクトリに格納されています)。
cavityFine ケースの system ディレクトリ内にある cotrolDict を編集します。startTime を 0.5 に、
startFrom を startTime にセットして保存します。そして blockMesh を走らせます。
これで mapFields を実行する準備が出来ました。mapFields の実行には入力ケースのディレクトリ指定
が必要です。「フィールドが首尾一貫している」ので ‒consistent オプションを使います。
mapFields $FOAM_RUN/tutorials/icoFoam/cavity ‒consistent
mapFields が実行され、次のように出力されます。
1-5.4. 設定の調整
全てのセルの寸法が半分になったので、1 より小さいクーラン数を維持するために、時間ステップを半
分にします。
まずは controlDict を開き、deltaT を 0.0025 にセットします。今まではフィールドデータを固定のステップ
回数の下での時間間隔で出力する方法を示していましたが、今回は固定の計算時間でデータ出力を
指定する方法を示します。controlDict の writeControl を runTime に、writeInterval を 0.1 にセットします。
このようにすることでケースは粗いメッシュでの解を入力条件として計算を始めるので、定常状態に集
束するには適切な短い時間だけ動かします。
従って endTime は 0.7 秒でよいことになります。セッティング確認を行い、ケースを保存します。
1-5.5. バックグラウンドプロセスとしてコードを動かす
icoFoam をバックグラウンドプロセスとして動かしてみます。そしてその最終結果を後で見ることが出来
るように log ファイルに出力します。
そのためには、cavityFine ディレクトリにおいて、ターミナルにコマンドを入力します。
icoFoam > log &
cat log
1-5.6. 精密なメッシュによるベクトルプロット
各々の新しいケースは本質的には単なる Pipline Browser に現れる他のモジュールですので、
ParaView で同時に複数のケースを開くことができます。ただし、ParaView で新しいケースを開くときに
は選んだデータが拡張子を含むファイル名である必要があります。
しかし OpenFOAM において、各々のケースは特定のディレクトリ構造の中に拡張子無しの状態で、複
数のファイルに保存されます。
解決法としては、ParaFoam スクリプトが自動的に拡張子「.OpenFoam」と付けたダミーファイルを作成す
ることになっています。故に、cavity ケースモジュールは cavity.OpenFOAM と呼ばれています。
ParaView 内から他のケースディレクトリをあきたい場合、ダミーファイルを作成しなければいけません。
例えば cavityFine ケースの読み込みには、ターミナルに次のコマンドを入力し、ファイルの作成をしま
す。
cd $FOAM_RUN/tutorials/icoFoam
touch cavityFine/cavityFine.foam
※コマンド解説※
touch…指定したファイル名が存在しない場合、空のファイルを作成する。
File メニューから Open Data を選んでディレクトリツリーを辿り、cavityFine.foam を選ぶことで cavityFine
ケースを ParaView に読み込むことが可能となります。
そして ParaView で精密なメッシュの結果のベクトルプロットを作成することが出来ます。同時に両方の
ケースの Glyph を見られるようにすることによって、cavityFine ケースのプロットを cavity ケースと比較
することが出来ます。
図 14 グラフ作成のためのフィールド選択
1-5.7. グラフの描画
OpenFOAM は速度のスカラー値を抽出して 2 次元のグラフに描画したい場合のデータの取り扱いに長
けています。
データ操作のための特別なユーティリティが多数あり、単純な計算を foamCalc によって組み合わせる
ことができます。次のようにユーティリティを指定して実行します。
foamCalc <calcType> <fieldName1 … fieldNameN>
処理を規定する<calcType>には、div、components、mag、magGrad、magSpr を指定可能です。
calcType のリストを見るには、意図的に無効な操作をします。すると、エラーメッセージとともにリストを
見ることが出来ます。
例えばディレクトリの指定で、「icoFoam」というワードを入れずにコマンドを打ってみます。
seminar@ubuntu-vm:~/OpenFOAM/seminar-1.5.x/run/tutorials/icoFoam/cavityFine$ foamCa
lc coponents U -case $FOAM_RUN/tutorials/(本当はココにicoFoamと入る)cavityFine
するとエラーメッセージとCalcTypeのリストが表示される。
Selecting calcType coponents
unknown calcType type coponents, constructor not in hash table
Valid calcType selections are:
5
(
div
components
mag
magGrad
magSqr
)
#0
Foam::error::printStack(Foam::Ostream&) in
"/home/seminar/OpenFOAM/OpenFOAM-1.5.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#1
Foam::error::abort() in
"/home/seminar/OpenFOAM/OpenFOAM-1.5.x/lib/linuxGccDPOpt/libOpenFOAM.so"
#2
Foam::calcType::New(Foam::word const&) in
"/home/seminar/OpenFOAM/OpenFOAM-1.5.x/lib/linuxGccDPOpt/libfoamCalcFunctions.so"
#3
__gxx_personality_v0 in
"/home/seminar/OpenFOAM/OpenFOAM-1.5.x/applications/bin/linuxGccDPOpt/foamCalc"
#4
__libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#5
__gxx_personality_v0 in
"/home/seminar/OpenFOAM/OpenFOAM-1.5.x/applications/bin/linuxGccDPOpt/foamCalc"
From function calcType::New()
in file calcType/newCalcType.C at line 43.
FOAM aborting
Aborted (core dumped)
ケースにて「foamCalc components U」を動かすと、各時刻のディレクトリから速度のベクトル場を読み
込み、各ディレクトリに各軸方向成分のスカラー場 Ux、Uy、Uz を書き出します。foamCalc は cavity、
cavityFine のどちらに対しても実行可能です。例として cavity に対して実行します。
foamCalc components U ‒case $FOAM_RUN/tutorials/icoFoam/cavity
それぞれのコンポーネントが ParaView 内でグラフとして描画されます。容易に形式やラベルの調整が
可能です。描画をする前に新しく生成された Ux、Uy、Uz のデータを ParaView に読み込ませる必要があ
ります。
作業をしている基本のモジュール、この場合は cavityFine.foam の Properties パネル下部にある Update
GUI ボタンにチェックをします。Apply をクリックすると Vol Field Status ウィンドウに新しいデータが読み
込まれます。境界領域が Region Status パネルで選択されていると境界部分のデータ補完が不適切に
行われることになるので、Region Status パネルの「movingwall・fixedwall・frontAndBack」という、3 つの
パッチの選択を解除しもう一度 Apply を押します。
次に ParaView の図を描画します。
まずは描画したいモジュールの選択をし、Filter から Data Analisys を選択し、Plot Over Line フィルター
を選択します。すると 3D View ウィンドウの側に新しい XY Plot ウインドウが開きます。Properties ウィン
ドウで線の終点を指定すると Pｒobe Line モジュールが作成されるので、Apply を押すと XY Plot ウイン
ドウに図が描画されます。Display パネルの Plot Type を Scatter Plot に設定し、Attribute Mode Point
Data にします。Use Data Array オプションの arc_length で X 軸データが cavity の基部からの距離にな
ります。
Displays ウインドウの Line Series パネルから表示するデータを選択することができます。表示されてい
るスカラー場のリストから、ベクトルの大きさや成分を初期値とすることもできます。つまり Ux を
foamCalc から計算する必要はありませんが、今回は Ux 以外の系列の選択はすべて解除します。
選択した系列の上の四角形の色が線の色です。この上でダブルクリックをすれば容易に変更可能です。
グラフの初期化には XY Plot 自体を動かします。カーソルがグラフ上にある状態で右クリックをして、メ
ニューから Properties を選択します。各軸のタイトルや形式の設定をする Chart Options ウインドウが
表示されます。各軸のメニューはダブルクリックをして Layout and Tytle にすることで拡大することがで
きます。フォントや色、軸名の位置、軸の値の範囲や線形か対数かといった設定が行うことができま
す。
図 15 作成した cavity.OpenFOAM ケースのグラフ
■1-6. 勾配メッシュ
天井駆動キャビティー問題のために壁に向かって勾配をつけた 20×20 セルのメッシュを作成し、1-5.2.
の細かいメッシュの結果を初期条件として勾配付けされたメッシュに適用します。そして、勾配付けされ
たメッシュの結果を前のメッシュの結果と比較します。blockMeshDict の書き換えはとても重要です。今
回この部分を使ったケースは cavityGrade として icoFoam ディレクトリの中に入っています。
1-6.1. 勾配メッシュの作成
ここで、4 つの異なるメッシュ間隔の計算メッシュが計算領域となる上下左右のブロックに必要となりま
す。
図 16 cavity ケースの勾配メッシュのブロック構造(ブロック番号)
cavityGrade の constant/polyMesh サブディレクトリで blockMeshDict ファイルを見ることができます。
blockMeshDict の中身を示しますので確認してください。それぞれのブロックは x 方向、y 方向に 10 セ
ルを有し、もっとも大きなセルともっとも小さなセルとの大きさの比は 2 です。
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
convertTOMeters 0.1;
vertics
(
(0 0 0)
(0.5 0 0)
(1 0 0)
(0 0.5 0)
(0.5 0.5 0)
(1 0.5 0)
(0 1 0)
(0.5 1 0)
(1 1 0)
(0 0 0.1)
(0.5 0 0.1)
(1 0 0.1)
(0 0.5 0.1)
(0.5 0.5 0.1)
(1 0.5 0.1)
(0 1 0.1)
(0.5 1 0.1)
(1 1 0.1)
);
blocks
(
hex (0 1 4 3 9 10 13 12) (10 10 1) simpleGrading (2 2 1)
hex (1 2 5 4 10 11 14 13) (10 10 1) simpleGrading (0.5 2 1)
hex (3 4 7 6 12 13 16 15) (10 10 1) simpeGrading (2 0.5 1)
hex (4 5 8 7 13 14 17 16) (10 10 1) simpleGrading (0.5 0.5 1)
);
edges
(
);
patches
(
wall movingWall
(
(6 15 16 7)
(7 16 17 8)
)
wall fixedWalls
(
(3 12 15 6)
(0 9 12 3)
(0 1 10 9)
(1 2 11 10)
(2 5 14 11)
(5 8 17 14)
)
empty frontAndBack
(
(0 3 4 1)
(1 4 5 2)
(3 6 7 4)
(4 7 8 5)
(9 10 13 12)
(10 11 14 13)
(12 13 16 15)
(13 14 17 16)
)
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
blockMeshDict を確認したら、blockMesh を実行しましょう。そして paraFoam を実行すれば、勾配づけされたメッシ
ュを見ることが出来ます。
1-6.2. 計算時間と時間ステップの変更
もっとも速い速度と小さいセルが上蓋に面することになり、したがって、セクション 1-1.4.で示したように、
もっとも高いクーラン数が上蓋に面するセルに生じます。このようなことから上蓋に面するセルの大きさ
を見積もることは、本ケースにて適当な時間ステップを計算する上で有効です。
不一様なメッシュ勾配を使用している場合、blockMesh は形状に関する数列をもちいてセルの大きさを
算出します。長さ l に沿って、最初と最後のセルとの間に、比 R の n 個の計算セルが必要であるならば、
もっとも小さいセルの大きさ
は、次のように与えられます。
ここで、r はあるセルの大きさとその隣のセルの大きさとの比であり、次式で表されます。
そして、
cavityGrade ケースにおいては、各方向のセルの数は 10 であり、もっとも大きなセルと小さなセルとの
比は 2、ブロックの縦横は 0.05m です。
したがって、最小のセルサイズは 3.45mm となります。よって時間ステップは、クーラン数を１以下に抑
えるために 3.45ms 以下にしなければなりません。
有意な解析結果を得るために、時間ステップ deltaT を 0.0025 まで短くし、writeInterval を 40 とします。
これより解析結果は 0.1s 毎に書き出されることとなります。
startTime はその cavityFine ケースの最終的な条件、すなわち 0.7 に設定される必要があります。
cavity と cavityFine が規定されたランタイムの中でよく収束させるためには、cavityGrade ケースのラン
タイムを 0.1 に設定、すなわち endTime を 0.8 とします。
1-6.3. 解析場のマッピング
1-5.3.と同様 mapFields を使用して、cavityFine ケースの最終的な結果を cavityGrade ケースのメッシュ
にマッピングします。
cd $FOAM_RUN/tutorials/icoFoam/cavityGrade
mapFields $FOAM_RUN/tutorials/icoFoam/cavityFine ‒consistent
そして次に cavityGrade ケースディレクトリから icoFoam を実行します。ランタイム情報をモニタリングし
て、このケースの完全に収束した結果を見たら、1-5.6.と 1-5.7.で説明した後処理ツールを用いて他の
結果と比較します。以下にランタイム情報を示します。
Create time
Create mesh for time = 0.7
Reading transportProperties
Reading field p
Reading field U
Reading/calculating face flux field phi
Starting time loop
Time = 0.7025
Courant Number mean: 0.0501319 max: 0.434876
DILUPBiCG:
Solving for Ux, Initial residual = 0.040115, Final residual = 1.98324e-06,
No Iterations 4
DILUPBiCG:
Solving for Uy, Initial residual = 0.0596489, Final residual = 4.15659e-06,
No Iterations 4
DICPCG: Solving for p, Initial residual = 0.376256, Final residual = 8.94975e-07, No
Iterations 30
time step continuity errors : sum local = 3.71185e-09, global = -2.83723e-20, cumulative
= -2.83723e-20
DICPCG: Solving for p, Initial residual = 0.0807109, Final residual = 6.47566e-07, No
Iterations 30
time step continuity errors : sum local = 2.28593e-09, global = 2.38373e-19, cumulative
= 2.1e-19
ExecutionTime = 0.04 s
ClockTime = 0 s
.
.
.
Time = 0.8
Courant Number mean: 0.0500944 max: 0.42312
DILUPBiCG:
Solving for Ux, Initial residual = 6.81511e-07, Final residual = 6.81511e-07,
No Iterations 0
DILUPBiCG:
Solving for Uy, Initial residual = 1.15523e-06, Final residual = 1.15523e-06,
No Iterations 0
DICPCG: Solving for p, Initial residual = 1.25334e-06, Final residual = 7.8711e-07, No
Iterations 4
time step continuity errors : sum local = 3.41885e-09, global = 4.35779e-19, cumulative
= 2.42458e-18
DICPCG: Solving for p, Initial residual = 1.22684e-06, Final residual = 7.7321e-07, No
Iterations 2
time step continuity errors : sum local = 3.39185e-09, global = 2.70157e-19, cumulative
= 2.69474e-18
ExecutionTime = 0.18 s
End
ClockTime = 0 s
そして、グラフ化したものを示します。
図 17 cavityGrade ケースのグラフ結果
図 18 cavityFine のグラフ
■1-7. レイノルズ数の増大
これまで解いた cavity、cavityFine、cavityGrade というケースらはレイノルズ数が一様に 10 でした。これ
は大変に低い条件であり、cavity の底部中央に小さな二次渦を伴うのみだったので迅速に安定解を解
くことが可能でした。しかしここで、例えばレイノルズ数を 50 に上げるとしたら、収束解を得るのにより長
い時間を要します。
そこで cavity ケースのメッシュを初期条件として使用します。Cavity ケースディレクトリを、cavityHighRe
という名前でコピーします。icoFoam ディレクトリに移動してから行います。
cd $FOAM_RUN/tutorials/icoFoam
cp ‒r cavity cavityHighRe
1-7.1. プリプロセッシング
コピーを終えたら cavityHighRe ケースに入ります。そして transportProperties 辞書を編集します。
レイノルズ数を 10 倍増加させるには動粘性係数を 10 分の 1 にしなくてはなりません(今回なら 0.01→
0,001)。
これで cavity ケースの実行結果からリスタートしてこのケースを実行することが出来ます。
startFrom キーワードを latestTime に設定し直すことによって icoFoam は最新の時間ディレクトリを所期
データとして使用します(例えば 0.5)。endTime は 2 秒に設定し、本ケースを保存します。
1-7.2. コード実行
まずはケースディレクトリから icoFoam を実行し、ランタイム情報を見ます。
このとき、バックグラウンドでジョブを実行するときには以下の UNIX コマンドを用います。
◇ nohup…ユーザーがログアウトしても稼動し続けるコマンド。
◇ nice…カーネル・スケジューラのジョブの優先順位を変えるコマンド。-20 が最優先で 19 は最も低
い優先度。
これらのコマンドは、例えばユーザーがリモートマシンでケースを実行できるよう設定し、頻繁にモニタ
リングしなくてもいいという場合、リモートマシンではケース実行をあまり優先させたくは無いでしょうが、
そのような場合に便利です。
その場合、ユーザーは nohup コマンドで稼動しているリモートマシンをログアウトしてジョブを実行し続
けることが出来ます。一方 nice は優先度を 19 と設定します。
cd $FOAM_RUN/tutorials/icoFoam
nohup nice ‒n 19 icoFoam > log &
cat log
前述の解析方法では icoFoam は速度 U の計算が止まっても、それよりもずっと長い間、もしくは解析が
終わるまで圧力 p の計算をし続けていました。実際には icoFoam が一旦 U の計算を止め、p の初期の
残差、寛容範囲以下は fvSolution 辞書で決められた公差(一般的には 10^-6)を下回ると結果が効率的
に収束するので、フィールドデータを一旦時間ディレクトリに書き出して計算を止めることが出来ます。
例として、cavityHighRe ケースの収束の log ファイルを以下に示します、1.62 秒後に速度はすでに収束
し、初期圧力の残差は小さくなります。
log において No Iterations 0 は、U の計算の停止を示します。
1
2
Time = 1.63
3
4
Courant Number mean: 0.108642 max: 0.818175
5
DILUPBiCG:
6
No Iterations 0
7
DILUPBiCG:
8
No Iterations 0
9
DICPCG: Solving for p, Initial residual = 3.54721e-06, Final residual = 7.13506e-07,
Solving for Ux, Initial residual = 7.86044e-06, Final residual = 7.86044e-06,
Solving for Uy, Initial residual = 9.4171e-06, Final residual = 9.4171e-06,
10
No Iterations 4
11
time step continuity errors : sum local = 6.46788e-09, global = -9.44516e-19,
12
cumulative = 1.04595e-17
13
DICPCG:
14
No Iterations 3
15
time step continuity errors : sum local = 8.67501e-09, global = 7.54182e-19,
16
cumulative = 1.12136e-17
17
ExecutionTime = 1.02 s ClockTime = 1 s
Solving for p, Initial residual = 2.15824e-06, Final residual = 9.95068e-07,
18
19
Time = 1.635
20
21
Courant Number mean: 0.108643 max: 0.818176
22
DILUPBiCG:
23
No Iterations 0
24
DILUPBiCG:
25
No Iterations 0
26
DICPCG:
27
No Iterations 4
28
time step continuity errors : sum local = 8.15435e-09, global = -5.84817e-20,
29
cumulative = 1.11552e-17
30
DICPCG:
Solving for Ux, Initial residual = 7.6728e-06, Final residual = 7.6728e-06,
Solving for Uy, Initial residual = 9.19442e-06, Final residual = 9.19442e-06,
Solving for p, Initial residual = 3.13107e-06, Final residual = 8.60504e-07,
Solving for p, Initial residual = 2.16689e-06, Final residual = 5.27197e-07,
31
No Iterations 14
32
time step continuity errors : sum local = 3.45666e-09, global = -5.62297e-19,
33
cumulative = 1.05929e-17
34
ExecutionTime = 1.02 s ClockTime = 1 s
■1-8. 高レイノルズ数流れ
上記が終わったら praFoam を起動させて結果を確認した後、速度ベクトルを示します。
図 19 HighRe の速度ベクトル
計算領域の角における二次渦が幾分増大していることがわかります。 このような時、ユーザーは粘性
係数を下げることによりレイノルズ数を増大させた計算ケースを再度実行できます。 渦の数が増加す
るにともない、より複雑な流れを解くために当該領域でのメッシュ解像度を上げる必要が出てきます。
更に、レイノルズ数は収束に要する時間を増加させます。 このような場合、残差をモニターし、解を収
束させるために endTime を延長したほうがよいでしょう。
空間および時間解像度の増加を要することは、流れが乱流域に移行するという非現実的な状態となり、
解法の安定性の問題が生じることとなります。 もちろん、多くのエンジニアリングな問題は極めて高い
レイノルズ数条件となっており、乱流挙動を直接解くのに多くのコストを負担することとなり実行不可能
であります。
Reynolds-averaged stress(RAS)にかわり乱流モデルが平均流れの挙動を解くのに用いられ、ゆらぎの
統計値が計算されています。 壁関数を伴う標準
ャビティケース（レイノルズ数
◇乱流エネルギー
◇乱流消散速度
モデルが本チュートリアルの上面が移動するキ
）を解くのに用いられています。 二つの追加変数が解かれています。
乱流のための追加の方程式およびモデルは turbFoam と呼ばれる OpenFOAM ソルバーにおいて実行
されます。
1-8.1. プリプロセッシング
$FOAM_RUN/tutorials/turbFoam ディレクトリの cavity ケースに移動します。
今まで通り、blockMesh を走らせ、メッシュを生成します。壁関数付き標準
モデルを用いる場合は、
壁近傍のセルにおける流れがモデル化されることにより、壁方向へのメッシュ勾配は必ずしも必要では
ありません。
と のファイル(k と Epsilon )を開き、境界条件を確かめます。壁タイプの境界条件の選択には、 に
ついては zeroGradient 境界条件を、
については fixedValue 0 を指定します。
と の初期値は解
法アルゴリズムにてゼロ割を避けるために、正の値を与えます。
、 の適当な初期条件は、速度変動
ここで
では
は
と乱流長さスケール を用いて設定することができます。
モデルの定数であり、その値は 0.09 です。カーデシアン座標系(＝デカルト座標系)
は、
で表されます。 各項は x、y、z 方向速度ゆらぎ成分です。
ここで、初期乱流が等方的であると仮定します。例えば、
となり、これら速度
は上面速度の5%に等しく、また、乱流長さスケール はボックス幅 0.1m の 20%に等しいとすると、次のよ
うに表されます。
上記のとおり
と 0 です。
、 を設定してください。
と
に対する初期条件は前と同じように、それぞれ(0,0,0)
次いで、transportProperties ディクショナリの層流動粘度を設定します。
レイノルズ数
を実現するために動粘度を
m にする必要があります。
RASProperties ディクショナリを開き、RAS 乱流モデルを選択します。乱流モデルは RASModel エントリ
ーで選択されます。 多くの使用可能なモデルが与えられていますので、ユーザーは標準
モデル
を表す kEpsilon をここでは選択します。 そして、turbulence のスイッチを on にします。 乱流モデルに
関 す る 係 数 は 標 準 デ ィ ク シ ョ ナ リ の kEpsilonCoeffs 以 下 に 、 ま た 、 同 デ ィ ク シ ョ ナ リ に
wallFunctionCoeffs の設定もあります。
クーラン数の制限を満たすために deltaT を 0.005s に設定し、endTime は 10s とします。
しっかりと確認しましょう。
1-8.2. コード進行
ケースディレクトリに入り、"turbFoam"とタイプすることで turbFoam を実行します。10 秒までのランタイ
ム情報を示します。
.
Time = 10
Courant Number mean: 0.0643285 max: 0.251781
DILUPBiCG:
Solving for Ux, Initial residual = 3.87417e-05, Final residual = 4.99184e-06, No
Iterations 1
DILUPBiCG:
Solving for Uy, Initial residual = 3.76683e-05, Final residual = 2.72331e-07, No
Iterations 1
DICPCG: Solving for p, Initial residual = 5.08098e-05, Final residual = 4.49886e-07, No Iterations
19
time step continuity errors : sum local = 6.58961e-10, global = -2.88024e-19, cumulative =
-1.15313e-17
DICPCG: Solving for p, Initial residual = 3.85791e-06, Final residual = 9.94715e-07, No Iterations
4
time step continuity errors : sum local = 1.61121e-09, global = -5.87629e-19, cumulative =
-1.21189e-17
DILUPBiCG:
Solving for epsilon, Initial residual = 5.52647e-05, Final residual = 7.15286e-08, No
Iterations 1
DILUPBiCG:
Solving for k, Initial residual = 0.0001052, Final residual = 1.27527e-07, No Iterations
1
ExecutionTime = 11.9 s ClockTime = 27 s
End
粘性が小さいこの計算ケースでは、移動している上面近傍の境界層は極めてうすく、そして上面に面
するセルは比較的大きいことから、上面速度よりもそれらセル中心の流体速度は極めて小さいです。
事実、100 時間ステップ後、上面に隣接したセルにおける速度は、上限である 0.2m/s 程度ですので最
大クーラン数は 0.2 以上にはなりません。
クーラン数がより１に近接するように時間ステップを大きくし、解析時間を増やすことは理にかなってい
ます。従って deltaT を 0.02s にセットし直し、startFrom を latestTime にセットします。
この操作は、turbFoam が最新のディレクトリ、例えば 10.0 からスタートデータを読み込むように指示す
るものです。
endTime は層流条件よりも収束に時間を要するため、20s にセットします。
従来通り計算をリスタートし、解析の収束をモニターするため turbFoam を実行します。
解析が進行したら、連続した時間における結果を見てください。そして解析が安定状態に収束するか、
もしくは周期的に振動しているか確認してください。後者の場合には、収束は決して起こりませんが、結
果が不正確であるという意味ではありません。
DILUPBiCG:
Solving for k, Initial residual = 0.000213662, Final residual = 6.14921e-06, No
Iterations 1
ExecutionTime = 3.16 s ClockTime = 7 s
Time = 19.26
Courant Number mean: 0.259588 max: 1.00692
DILUPBiCG:
Solving for Ux, Initial residual = 0.000345851, Final residual = 7.62684e-06, No
Iterations 1
DILUPBiCG:
Solving for Uy, Initial residual = 0.000143404, Final residual = 5.96648e-06, No
Iterations 1
DICPCG: Solving for p, Initial residual = 0.00769462, Final residual = 9.46176e-07, No Iterations
25
time step continuity errors : sum local = 1.84618e-08, global = -3.83812e-19, cumulative =
-2.47942e-17
DICPCG: Solving for p, Initial residual = 0.000961491, Final residual = 8.3159e-07, No Iterations
13
time step continuity errors : sum local = 1.79586e-08, global = -3.91753e-19, cumulative =
-2.5186e-17
DILUPBiCG:
Solving for epsilon, Initial residual = 0.000174142, Final residual = 3.72344e-06, No
Iterations 1
DILUPBiCG:
Solving for k, Initial residual = 0.000216359, Final residual = 1.03873e-07, No
Iterations 2
ExecutionTime = 3.16 s ClockTime = 7 s
.
.
.
Time = 20
Courant Number mean: 0.259597 max: 1.00659
DILUPBiCG:
Solving for Ux, Initial residual = 0.00036408, Final residual = 3.70552e-07, No
Iterations 2
DILUPBiCG:
Solving for Uy, Initial residual = 0.000152561, Final residual = 9.39551e-06, No
Iterations 1
DICPCG: Solving for p, Initial residual = 0.00731517, Final residual = 9.0445e-07, No Iterations
26
time step continuity errors : sum local = 1.76682e-08, global = -2.062e-18, cumulative = -2.70389e-17
DICPCG: Solving for p, Initial residual = 0.000901978, Final residual = 9.61862e-07, No Iterations
12
time step continuity errors : sum local = 2.19134e-08, global = 1.44525e-18, cumulative =
-2.55936e-17
DILUPBiCG:
Solving for epsilon, Initial residual = 0.000168746, Final residual = 3.05665e-06, No
Iterations 1
DILUPBiCG:
Solving for k, Initial residual = 0.000208709, Final residual = 6.02202e-06, No
Iterations 1
ExecutionTime = 3.46 s ClockTime = 8 s
End
■1-9. ケース形状の変更
計算ケースの形状を変更して新たな解析を行いたい場合、そのスタート条件としてオリジナルの解析
の全てないし一部を保持しておくことは有効でしょう。しかしこれは少し複雑になります。なぜなら、オリ
ジナルの解析の物理量が、新しい解析ケースの物理量と一致しないからです。
しかし mapFields ユーティリティは、形状や境界のタイプもしくはその両者が不一致な場を位置づけるこ
とが出来ます。
例であるように、icoFoam ディクショナリ内にある cavityClipped ケースを開きます。このケースは、標準
的 な cavity ケ ー ス か ら 成 り ま す が 、 底 部 右 側 、 長 さ 0.04m の 正 方 形 を 除 い た も の で あ り 、
blockMeshDict は以下のようになっています。
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.1;
vertices
(
(0 0 0)
(0.6 0 0)
(0 0.4 0)
(0.6 0.4 0)
(1 0.4 0)
(0 1 0)
(0.6 1 0)
(1 1 0)
(0 0 0.1)
(0.6 0 0.1)
(0 0.4 0.1)
(0.6 0.4 0.1)
(1 0.4 0.1)
(0 1 0.1)
(0.6 1 0.1)
(1 1 0.1)
);
blocks
(
hex (0 1 3 2 8 9 11 10) (12 8 1) simpleGrading (1 1 1)
hex (2 3 6 5 10 11 14 13) (12 12 1) simpleGrading (1 1 1)
hex (3 4 7 6 11 12 15 14) (8 12 1) simpleGrading (1 1 1)
);
edges
(
);
patches
(
wall lid
(
(5 13 14 6)
(6 14 15 7)
)
wall fixedWalls
(
(0 8 10 2)
(2 10 13 5)
(7 15 12 4)
(4 12 11 3)
(3 11 9 1)
(1 9 8 0)
)
empty frontAndBack
(
(0 2 3 1)
(2 5 6 3)
(3 6 7 4)
(8 9 11 10)
(10 11 14 13)
(11 12 15 14)
)
);
mergePatchPairs
(
);
// ************************************************************************* //
blockMesh を実行してメッシュを生成します。パッチ(境界)は cavity と同様に設定されています。物理量
の適用の過程を明確にするために、元となるケース cavity で movingWall だった upper wall(天井)は lid
という名前に変更されています。
パッチが一致しない場合、全ての物理量のデータが元のケースからマップされるという保証はありませ
ん。残っているデータは元のケースと同一であるべきです。
よってマッピングする前に時間のディレクトリに物理量のデータが存在している必要があります。
controlDict の startTime が 0.5 秒に設定されているので cavityClipped ケースにおけるマッピングは時
刻 0.5 秒に予定されています。よって、初期状態の物理量のデータ、例えば時刻 0 秒からコピーする必
要があります。
cavityClipped ケースに入り、時刻のコピーを行いますのでターミナルに以下を入力します。
cd $FOAM_RUN/tutorials/icoFoam/cavityClipped
cp ‒r 0 0.5
データのマッピング前に 0.5 秒における形状と物理量の状況を見ておきます。
速度場と圧力場を cavity から cavityClipped にマップしようとしています。パッチが一致しないため
system ディレクトレリの mapFieldsDict を編集します。patchMap と cuttingPatches という 2 つの項目が
あります。patchMap リストは元となる物理量のパッチとマッピング対象となる 物理量のパッチを含みま
す。
対象物理量のパッチに元となる物理量のパッチの値を引き継ぎたいときに利用します。
cavityClipped において lid の境界値を cavity の movingWall から引き継ぎたいので次のように patchMap
に入力します。
patchMap
(
lid movingWall
)：
cuttingPatches リストは対象パッチを削除した、元の場の内部の値を写像した対象のパッチを含みます。
今回は fixedWalls を実例に用います。
cuttingPatches
(
fixedWalls
)：
そして mapFields を実行します。
cd mapFields $FOAM_RUN/tutorials/icoFoam/cavity
paraFoam で速度ベクトルを示します。境界パッチは元のケースからの値が引き継がれているはずで
す。
図 20 cavity ケースで解いた速度場を cavityClipped 上にマッピングした図
比較のため、次に 0.6 秒の時の速度場を示します。修正されたパッチに従って、流速も変化します。
図 21 速度フィールドの cavityClipped の解法
右下の切り取られた壁面付近で渦が発生しているのがわかります。
これらをしっかりと確認しましょう。