ユーザーガイド 1.5系 - 2.3 ダムの決壊
出典: OFWikiJa
[編集] 2.3 ダムの決壊
このチュートリアルでは、interFoamを用いて、単純化したダム決壊の2次元問題を解くことにします。 この問題の特徴は、くっきりとした界面や自由表面によって隔てられている2つの流体による非定常の流れ場であることです。 interFoamにおける2相流体を解くアルゴリズムは、Volume of fluid(VOF)法によるものであり、ここでは特別な輸送方程式を解いて、計算格子における2相の体積分率、もしくは相比率を決定します。 各物理量は、この相比率に(各流体の)密度をかけた平均的な値として算出されます。 個々の物質の界面は、VOF法ではその性質上明示的には解かれず、相比率場の特性として浮き上がってくるということになります。 相比率は0から1の間の任意の値をとり得るため、界面は決してくっきりと定義されませんので、本来のくっきりとした界面が存在するべき領域の周辺を、(計算上の)界面がぼんやりと占めることになります。
計算条件では、貯水池の左側に、膜で仕切られた水柱が最初存在します。 時刻t=0sに、膜が取り除かれて、水柱が崩れだします。崩壊している間、水流はは貯水槽の底にある出っ張りにぶつかり、いくつかの気泡を含む、複雑な流れ場の様相を呈します。 計算形状と初期条件は図2.20に示しました。
図2.20 ダム決壊の計算形状
[編集] 2.3.1 格子の生成
$FOAM_RUN/tutorials/interFoamにあるdamBreakのケースディレクトリに移動しましょう。 前に述べた方法でblockMeshを走らせて格子を生成してください。 このdamBreakの格子は5つのブロックで構成されます。 blockMeshDictの中身を以下に示します。
23 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 24 25 convertToMeters 0.146; 26 27 vertices 28 ( 29 (0 0 0) 30 (2 0 0) 31 (2.16438 0 0) 32 (4 0 0) 33 (0 0.32876 0) 34 (2 0.32876 0) 35 (2.16438 0.32876 0) 36 (4 0.32876 0) 37 (0 4 0) 38 (2 4 0) 39 (2.16438 4 0) 40 (4 4 0) 41 (0 0 0.1) 42 (2 0 0.1) 43 (2.16438 0 0.1) 44 (4 0 0.1) 45 (0 0.32876 0.1) 46 (2 0.32876 0.1) 47 (2.16438 0.32876 0.1) 48 (4 0.32876 0.1) 49 (0 4 0.1) 50 (2 4 0.1) 51 (2.16438 4 0.1) 52 (4 4 0.1) 53 ); 54 55 blocks 56 ( 57 hex (0 1 5 4 12 13 17 16) (23 8 1) simpleGrading (1 1 1) 58 hex (2 3 7 6 14 15 19 18) (19 8 1) simpleGrading (1 1 1) 59 hex (4 5 9 8 16 17 21 20) (23 42 1) simpleGrading (1 1 1) 60 hex (5 6 10 9 17 18 22 21) (4 42 1) simpleGrading (1 1 1) 61 hex (6 7 11 10 18 19 23 22) (19 42 1) simpleGrading (1 1 1) 62 ); 63 64 edges 65 ( 66 ); 67 68 patches 69 ( 70 wall leftWall 71 ( 72 (0 12 16 4) 73 (4 16 20 8) 74 ) 75 wall rightWall 76 ( 77 (7 19 15 3) 78 (11 23 19 7) 79 ) 80 wall lowerWall 81 ( 82 (0 1 13 12) 83 (1 5 17 13) 84 (5 6 18 17) 85 (2 14 18 6) 86 (2 3 15 14) 87 ) 88 patch atmosphere 89 ( 90 (8 20 21 9) 91 (9 21 22 10) 92 (10 22 23 11) 93 ) 94 ); 95 96 mergePatchPairs 97 ( 98 ); 99 100 // ************************************************************************* //
[編集] 2.3.2 境界条件
constant/polyMeshディレクトリのboundaryファイルを見ることでblockMeshで生成された境界の形状を確認しましょう。leftWall, rightWall, lowerWall, atmosphere, defaultFacesの5つの境界パッチがあります。パッチの種類について理解しておきましょう。 atmosphereは何の属性もなく、単に境界条件によって規定される標準のパッチです。defaultFacesは、本ケースでは2次元であるためパッチの法線方向を解析の対象としないため、emptyとします。leftWall, rightWall, lowerWallはそれぞれwallです。plainと同様にwallもメッシュについて形状や位相の情報を持ちませんが、wallとして識別することができるので、アプリケーションに特殊な壁表面のモデリングを明示するためにwallと定義しています。
interFoamのソルバーが、界面と壁面との接点における表面張力に対するモデルを含んでいる、というのがよい例です。このモデルはgamma(構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \gamma )場のgammaContactAngleの境界条件と関連付けられています。 その場合、静的な接触角度theta0 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \theta_0 や、前縁や後縁における動的な接触角度であるthetaA 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \theta_A とthetaR 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \theta_R 、そして、動的な接触角度において速度に比例する係数uThetaを指定する必要があります。
このチュートリアルでは、壁面と界面間の表面張力による効果を無視することにしたいと思います。 それは、静的な接触角度を 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \theta_0=90^{\circ} に、速度比例係数を0と設定することで可能です。 しかしながら、壁の境界条件として、通常のwallタイプの境界条件を指定する別なやり方もあります。 この場合、gammaに対してgammaContactの境界条件を設定する代わり、zeroGradientの境界条件を課すことになります。
topの境界は大気に対して開放されていることから、atmosphereタイプの境界条件を与えます。また、この2次元問題において、前後の面のようなdefaultFacesの境界条件は通常通りemptyタイプです。
[編集] 2.3.3 初期条件の設定
これまでのケースと異なり、ここでは相比率構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \gamma に対して、以下のような非一様な初期条件を与えます。
これは、setFieldsユーティリティーを実行することによって行ないます。 この実行にはsystemディレクトリ内のsetFieldsDictを必要とします。 このケースにおけるsetFieldsDictファイルの内容を以下に示します。
23 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
24
25 defaultFieldValues
26 (
27 volScalarFieldValue gamma 0
28 volVectorFieldValue U (0 0 0)
29 );
30
31 regions
32 (
33 boxToCell
34 {
35 box (0 0 -1) (0.1461 0.292 1);
36
37 fieldValues
38 (
39 volScalarFieldValue gamma 1
40 );
41 }
42 );
43
44
45 // ************************************************************************* //
ここで、defaultFieldValuesは場の規定値を設定するものであり、regionsのサブ辞書において別途指定されない場合に場に与えられる値です。 regionsのサブ辞書は、指定された領域において、規定値を上書きするfieldValuesを含んだサブ辞書のリストを含んでいます。 領域の定義は、ある位相幾何学的な制約に基づいて、点や格子、界面等の集合を生成するtopoSetSourceによって行います。 ここでは、boxToCellを使って、大きいほうと小さいほうの二つの座標点で定義されるバウンディング・ボックスを生成し、液体の領域における格子の集合を定義しています。 また、この領域における相比率構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \gamma を1と指定しています。
さて、setFieldsを他のプログラムと同様に起動してください。 そうしたら、paraFoamを用いて、初期のgamma場が図2.21のように望むような分布になっているかどうか確かめてください。
図2.21 相比率gammaの初期条件
[編集] 2.3.4 流体の物性値
constantディレクトリのtransportPropertiesファイルを確認しましょう。このファイルは、phase1とphase2に分かれて、各流体の物性値を保持しています。
各相の輸送モデルは、transportModelによって選択されます。
ここで、動粘性係数がnuというキーワードで指定され、一定値であるNewtonianモデルを選んでください。
CrossPowerLawといったその他のモデルにおける粘性係数の指定は、この例におけるCrossPowerLawCoeffsといったように、<model>Coeffsという名のサブ辞書の中で行います。
密度の指定は、rhoキーワードで行います。
二つの相の間の表面張力は、sigmaキーワードで指定します。
このチュートリアルで用いた値を表2.3に挙げます。
| phase1の物性 | |||
| 動粘性率 | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): {\rm m}^{2}{\rm s}^{-1} | nu | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): 1.0\times 10^{-6} |
| 密度 | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): {\rm kgm}^{-3} | rho | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): 1.0\times 10^{3} |
| phase2の物性 | |||
| 動粘性率 | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): {\rm m}^{2}{\rm s}^{-1} | nu | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): 1.48\times 10^{-5} |
| 密度 | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): {\rm kgm}^{-3} | rho | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): 1.0 |
| 両phaseの物性 | |||
| 表面張力 | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): N{\rm m}^{-1} | sigma | 構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): 0.07 |
表2.1: 鋼の機械的性質
environmentalPropertiesの辞書は、重力加速度ベクトルを指定していますが、このチュートリアルでは、この値は(0,9.81,0)ms-2としてください。
[編集] 2.3.5 タイムステップの制御
自由界面の捕捉においては、タイムステップの制御は重要です。というのも、界面捕捉のアルゴリズムは、通常の流体計算に比べ、クーラン数Coに対してかなり鋭敏だからです。 理想的には、界面がある領域において、およそ0.2であるCoの上限値を超えてはいけません。 伝播速度の予測が容易であるようなケースでは、Coの制限を守るような固定したタイムステップを指定することができますが、より複雑なケースの場合タイムステップの指定はずっと困難になります。 そこで、interFoamでは、controlDictにおいて、タイムステップの自動修正を指定することをお勧めします。adjustTimeStepをonにして、Coの最大値maxCoを0.2にしましょう。 タイムステップの上限maxDeltaTはこのシミュレーションでは超えようのない値、たとえば1.0等に設定すればよいでしょう。
ただし、自動タイムステップ制御を用いると、ステップ自体は決して使いやすい値に丸められることはありません。 従って、固定のタイムステップ間隔でOpenFOAMに結果を出力させた場合、その時刻はきりの良い値になりません。 ところがこの自動タイムステップ制御を用いていても、OPENFORMでは決まった時刻に結果を出力するように指定することが可能です。 この場合、OpenFOAMは、結果の出力に指定された時刻ぴったりに合うように時間刻みを補正しつつ、自動時間刻みの制御を行います。 これを行うには、controlDict辞書におけるwriteControlに対して、adjustableRunTimeオプションを選んでください。 controlDict辞書の中身は以下のようになります。
23 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 24 25 application interFoam; 26 27 startFrom startTime; 28 29 startTime 0; 30 31 stopAt endTime; 32 33 endTime 1; 34 35 deltaT 0.001; 36 37 writeControl adjustableRunTime; 38 39 writeInterval 0.05; 40 41 purgeWrite 0; 42 43 writeFormat ascii; 44 45 writePrecision 6; 46 47 writeCompression uncompressed; 48 49 timeFormat general; 50 51 timePrecision 6; 52 53 runTimeModifiable yes; 54 55 adjustTimeStep yes; 56 57 maxCo 0.5; 58 59 maxDeltaT 1; 60 61 62 // ************************************************************************* //
[編集] 2.3.6 離散化スキーム
OpenFOAMにおける自由表面の扱いは、乱流の影響を考慮しません。 これは乱流モデルに対するレイノルズ平均モデルの手法が空気と水の間のごく 薄い表面の概念と適合しない事実に基づいています。 結果として、全ての自由表面シミュレーションは流体の直接数値シミュレーション(DNS)としてみることができます。
DNSはメッシュ数に対して一定の要求があり、それは我々のテストケースにお けるメッシュの解像度を遥かに超えています。 このソルバはOpenFormで開発されたMulti dimetional universal Limiter for explicit solution(MULES)法を用いて おり、基礎を成す数値的スキームやメッシュ構造から独立な段階分数の有界性を保存するために使います。 したがって、対流スキームの選択はそれが静的であるか、また閉鎖系であるかによって制限されません。
従って、運動量の式における対流項では風上差分を使用して計算を安定させる こととします。風上差分によって導入された数値拡散は計算を安定させること でしょう。 風上の設定は、fvSchemes辞書のdivSchemesサブディレクトリ 内に作成されており、運動量方程式における構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \nabla\cdot(\rho \phi \mathbf{U}) 項に対応する div(rho*phi,U)キーワードは、Gauss upwindを読みこみます。 相構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \gamma に対する2つの方程式における対流項には幾つかの注意を必要とします。適切な精度を維持しつつ有界性を保証しなくてはならないため、対流項内での補完間のためには、有界なスカラーに対するリミッタ付きの線形なスキームであるlimitedLinear01を用います。
リミッタ付きの線形なスキームは4.4.1に記述される係数を必要とします。ここで安定性を最重視して、キーワードdiv(phi,gamma)に対応する構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \nabla\cdot(\phi \gamma) 項、キーワードdiv(phirb,gamma)に対応する構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \nabla\cdot(\phi_{rb} \gamma) 項で構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \phi=1.0 の状態を選びます。
その他の離散化項は一般に決ったスキームを使用します。 それゆえfvSchemesディクショナリーのエントリーは以下のようになるべきです。
23 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
24
25 ddtSchemes
26 {
27 default Euler;
28 }
29
30 gradSchemes
31 {
32 default Gauss linear;
33 grad(U) Gauss linear;
34 grad(gamma) Gauss linear;
35 }
36
37 divSchemes
38 {
39 div(rho*phi,U) Gauss upwind;
40 div(phi,gamma) Gauss limitedLinear01 1;
41 div(phirb,gamma) Gauss limitedLinear01 1;
42 }
43
44 laplacianSchemes
45 {
46 default Gauss linear corrected;
47 }
48
49 interpolationSchemes
50 {
51 default linear;
52 }
53
54 snGradSchemes
55 {
56 default corrected;
57 }
58
59 fluxRequired
60 {
61 default no;
62 pd;
63 pcorr;
64 gamma;
65 }
66
67
68 // ************************************************************************* //
[編集] 2.3.7 線形ソルバーの制御
fvSolutionでは、PISOサブ辞書がinterFoamに特化した要素を含んでいます。 ここには、通常と同じく運動量方程式に対する反復数だけでなく、相方程式のPISOループに対する反復数も指定します。 特に興味深いものはnGammaSubCyclesとcGammaキーワードです。 nGammaSubCyclesは構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \gamma 方程式内のsub-cycleの数を表してあり、ここで、sub-cyclesは与えられた時間ステップ内での方程式に対する付加的な解の点数です。 それは、時間ステップや計算時間の莫大な増加なしで解を安定させることができるようにするものです。 ここでは、4つのsub-cycleを指定しており、構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): \gamma 方程式は実際の各時間ステップ内で4分の1の幅の時間スッテプで4回解かれていることを意味します。
cGammaキーワードは界面の圧縮を制御する要素です。 つまり、0は無圧縮に対応し、1は保存的な圧縮に対応し、1以上は拡張された界面の圧縮を意味します。 通常はこの例題で用いられている1.0の値が推奨されます。
[編集] 2.3.8 コードの実行
コードの実行方法については、前述のチュートリアルに詳細に記述しています。 以下を試してください。teeによって計算内容をターミナルウインドウに表示しつつ、ログファイルに記録します。
cd $FOAM_RUN/tutorials/interFoam interFoam | tee log
コードは、出力のコピーをlogファイル内に記録しつつ、対話形式で実行されます。
[編集] 2.3.9 事後処理
結果の事後処理は、通常の方法で行えます。 図2.22のように参照時間の経過とともに相比率gammaの発達を見ることができます。
(a)At t=0.25s
(b)At t=0.50s
図2.22: γフェーズのスナップショット
[編集] 2.3.10 並列計算
前述の例の結果はかなり目の粗い格子を使って得られました。ここでは格子の解像度を増やして再計算します。新しいケースは、一般的に1つのプロセッサでは計算するのに数時間を要するので、複数のプロセッサにアクセスしているのであれば、OpenFOAMの並列計算という能力を試してみてもよいでしょう。
まず初めに、damBreakケースのコピーをしてください。
cd $FOAM_RUN/tutorials/interFoam
mkdir damBreakFine
cp -r damBreak/0 damBreakFine
cp -r damBreak/system damBreakFine
cp -r damBreak/constant damBreakFine
新しいケースはdamBreakFineと名づけてください。新しいケースディレクトリを開いてblockMeshDictファイル内のblocksの記述を以下の様に変更してください。
blocks
(
hex (0 1 5 4 12 13 17 16) (46 10 1) simpleGrading (1 1 1)
hex (2 3 7 6 14 15 19 18) (40 10 1) simpleGrading (1 1 1)
hex (4 5 9 8 16 17 21 20) (46 76 1) simpleGrading (1 2 1)
hex (5 6 10 9 17 18 22 21) (4 76 1) simpleGrading (1 2 1)
hex (6 7 11 10 18 19 23 22) (40 76 1) simpleGrading (1 2 1)
);
上記で、入力はblockMeshDictファイルで表示されているように、つまりは、格子の密度を変更しなければなりません。例えば46 10 1という入力や1 2 1という格子幅の勾配の入力の様にです。入力内容が正しければ、格子を生成します。
ここで格子がdamBreakの例から変更されると、時刻0のディレクトリ内のgammaというフェーズ場を再初期化しなければなりません。というのもgammaは新しい格子とは合致しないいくつかの要素を含んでいるからです。ここで、Uやpというフィールドは変更する必要が無いことに注意しましょう。それらはuniformとして明記されておりフィールド内の要素の数と独立だからです。フィールドの初期化はシャープな界面を持つように行いたいものです。つまり、その要素がγ＝１かγ＝０を持つようにです。mapFieldsでのフィールドの更新は界面に補間された値０＜γ＜１が生成されるかもしれないので、setFieldsユーティリティーを以下の様に再実行したほうがよいでしょう。その前に初期条件の一様なγのバックアップファイル0/gamma.orgを0/gammaにコピーします。
cd $FOAM_RUN/tutorials/interFoam/damBreakFine cp -r 0/gamma.org 0/gamma setFields $FOAM_RUN/tutorials/interFoam damBreakFine
OpenFOAMで用いられる並列計算の手法はいわゆる領域分割であり、幾何形状やそれに関連する場が領域毎に分解されて、解析のため個々のプロセッサに割り当てられます。 そのため、並列計算を実行するために必要な最初の段階は、decomposeParを用いて領域を分解することです。decomposeParの設定はsystem ディレクトリにある、decomposeParDictというファイルです。他のユーティリティ同様、初期状態のファイルがユーティリティのソースコードのディレクトリ（$FOAM_UTILITIES/parallelProcessing/decomposePar）にあります。
最初の入力のnumberOfSubdomainsにおいて何個のサブ領域に分割するかを指定します。通常はこのケースに利用できるプロセッサの数と対応します。
このチュートリアルでは、分解の手法はsimpleで、対応するsimpleCoeffsは以下の基準の様に編集しましょう。領域は、構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): x ,構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): y ,構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): z 方向で部分かサブ領域に分けられ、各方向へのサブ領域の数はベクトルnとして与えられます。この幾何形状は2次元なので、3次元方向の構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): z には分割され得ず、それゆえ必ず構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): n_z は1になります。構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): n_x と構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): n_y は構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): x ,構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): y 方向の領域の分割数nを構成し、構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): n_x と構文解析失敗 (texvcプログラムが見つかりません。math/READMEを読んで正しく設定してください。): n_y の積で表されるサブ領域の数がnumberOfSubdomainに指定したものと等しくなる必要があります。隣接するサブ領域間のセル面の数を最小にしたほうがよいので、正方形の幾何形状では、x,y方向を均等に分割するのが良いでしょう。deltaキーワードは0.001に設定しましょう。
例として、4つのプロセッサで計算を実行するとします。numberOfSubdomainを4に、n=(2,2,1)に設定します。decomposeParDictを閉じて、decomposeParを実行します。decomposeParのスクリーンメッセージが確認でき、分解はプロセッサ間で均等に分配されたと表示されます。
セクション3.4に並列計算の方法についての詳細があるので参照してください。このチュートリアルでは並列計算の一例を示しているに過ぎません。openMPIを用いて標準のメッセージパッシングインターフェース(MPI)を実装しています。ここでは、テストとしてローカルホストのみの単独ノードで実行します。
mpirun -np 4 interFoam -parallel > log &
3.4.2に後述しますが、ケースが実行されるマシンのホストネームを列記したファイルを作っておけばネットワーク上のより多くのノードを使って計算することも可能です。ケースはバックグラウンドで実行し、進行状況をlogファイルで監視するのがよいでしょう。
図2．33: 並列プロセスケースでのプロセッサ2のメッシュ
(a)At t=0.25s
(b)At t=0.50s
図2.34: 正確なメッシュのγフェーズのスナップショット
[編集] 2.3.11 並列計算ケースのポスト処理
一度ケースの実行が完了したら、分解されたフィールドとメッシュはreconstructParを実行してポスト処理のために再統合します。コマンドラインから容易に実行できます。細かい格子による結果はFigure 2.30に表されます。インターフェースでの結果は粗い格子のものと比較して著しく改良されたことが見てとれます。
また、単に個々のプロセッサの領域を一つのケースと扱うことで、分解された領域の部分を個々にポスト処理することもできます。
例えば、paraFoamを以下の様に起動します。
paraFoam -case processor1
そしてprocessor1はParaViesのケースモジュールとして表れます。
simpleメソッドを使った領域分割を行なったprocessor1からFigure2.31のような格子が見えます。
