こんにちは（@t_kun_kamakiri）
今回も前回に引き続きCAE解析結果と実験データとの比較を行います。
前回の記事
今回の内容
3次元ダムブレイクの実験の指定した座標位置での水面の高さの時刻歴をOpenFOAMの結果と比較する。
時刻歴の重ねはgnuplotを使います。
OpenFOAM初心者でチュートリアルを動かしたことがある方を対象にしています。
DEXCS2020
OpenFOAM v2006
gnuplot 5.2
3次元ダムブレイクの実験データ
こちらの実験データのより以下の時刻歴データを計測しています。
- 圧力計測点P1~P8（前回の記事)
- 水面の高さH1~H4←今回OpenFOAMと比較する対象
前回の記事では圧力の時刻歴データと実験とOpenFOAMで比較を行いました。
本記事ではH1~H4の水面の高さの時刻歴をOpenFOAMで計測し実験データとの比較を行います。
解析実行
解析設定はこちらをご参考ください。
設定ファイルは再度こちらにも記載しておきます。
constant/turbulenceProperties
今回乱流モデルは使いません。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2006 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object turbulenceProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // simulationType laminar; // ************************************************************************* // |
今回こちら「$FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreakWithObstacle/」のチュートリアルを利用していますが、乱流モデルなしのケースファイルでも「k,nut,omega」ファイルがありますが、こちらは乱流モデルの際に使うので今回は特に設定をしなくも良いです（使っていないので）。
0.org/U
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2006 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volVectorField; location "0"; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { region_maxZ { type pressureInletOutletVelocity; value uniform (0 0 0); } ".*" { type noSlip; } } // ************************************************************************* // |
0.org/p_rgh
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2006 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object p_rgh; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [1 -1 -2 0 0 0 0]; internalField uniform 0; boundaryField { ".*" { type fixedFluxPressure; phi phiAbs; value uniform 0; } region_maxZ { type totalPressure; p0 uniform 0; } } // ************************************************************************* // |
0.org/alpha.water
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2006 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; location "0"; object alpha.water; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 0 0 0 0 0 0]; internalField nonuniform List<scalar> boundaryField { region_maxZ { type inletOutlet; inletValue uniform 0; value uniform 0; } ".*" { type zeroGradient; } } |
constant/dynamicFvMeshDict
チュートリアルは空気と水の界面でのメッシュを時刻歴で再分割するようにdynamicMeshが採用されています。
今回このような設定では解析時間がとてもかかるのと、結果処理がものすごく重いためノートPCで耐えられる計算ではないため、「dynamicFvMesh staticFvMesh;」を使って無効にしておきます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2006 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object dynamicMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dynamicFvMesh staticFvMesh; // |
constant/g
重力加速度はz方向の下向きに設定します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2006 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class uniformDimensionedVectorField; location "constant"; object g; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 -2 0 0 0 0]; value (0 0 -9.81); // ************************************************************************* // |
system/probes
今回はこちらの実験データの圧力計測点と同じ位置の圧力データを取得する設定を行います。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | /*--------------------------------*- C++ -*----------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Version: v2006 \\ / A nd | Website: www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- Description Writes out values of fields from cells nearest to specified locations. \*---------------------------------------------------------------------------*/ #includeEtc "caseDicts/postProcessing/probes/probes.cfg" //probes 設定の挿入 fields (p p_rgh); // プローブする場 probeLocations ( (0.8245001 0 0.0205) // プローブ座標 P1 (0.8245001 0 0.0605) // プローブ座標 P2 (0.8245001 0 0.1005) // プローブ座標 P3 (0.8245001 0 0.1405) // プローブ座標 P4 (0.8040 0 0.161 ) // プローブ座標 P5 (0.7640 0 0.161 ) // プローブ座標 P6 (0.7240 0 0.161 ) // プローブ座標 P7 (0.6840 0 0.161 ) // プローブ座標 P8 ); interpolationScheme cellPoint; // 補間スキーム // ************************************************************************* // |
- 「fields (p p_rgh);」で出力する物理量を指定します。
- probeLocationsで出力したい座標位置を指定します。
- 「interpolationScheme cellPoint」はセル値による線形重み付け補間スキーム
probesをsystem/controlDictで読み込みます。
水面高さを計測するためinterfaceHeight1を設定します。
system/controlDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2006 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application interFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 8; deltaT 0.001; writeControl adjustable; writeInterval 0.1; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable yes; adjustTimeStep yes; maxCo 0.5; maxAlphaCo 0.5; maxDeltaT 1; functions // 関数定義 { #includeFunc solverInfo(U,p,alpha.water); residuals { type solverInfo; libs ("libutilityFunctionObjects.so"); fields (".*"); } #includeFunc probes; // 壁面での圧力時系列モニター用のプローブ interfaceHeight1 { // Mandatory entries (unmodifiable) type interfaceHeight; //libs (“libfieldFunctionObjects.so”); libs (fieldFunctionObjects); // Mandatory entries (runtime modifiable) locations ( (0.496 0 0) (0.992 0 0) (1.488 0 0) (2.638 0 0) ); // Optional entries (runtime modifiable) alpha alpha.water; liquid true; direction (0 0 -1); interpolationScheme cellPoint; // Optional (inherited) entries writePrecision 8; writeToFile true; useUserTime true; region region0; enabled true; log true; timeStart 0; timeEnd 1000; executeControl timeStep; executeInterval 1; writeControl timeStep; writeInterval 1; } } // ************************************************************************* // |
「direction (0 0 -1);」は重力の方向を正しく指定する必要があります。
今回はz方向下向きに重力がはたらいているためこのような設定としました。
また、「postProcessing/interfaceHeight1/0/height.dat」には一つの座標に対して以下の出力結果が出されます。
- # hB : Interface height above the boundary：水面の高さ
- # hL : Interface height above the location：波の高さ
いずれも重力の方向を正しく指定することが重要です。
system/decomposeParDict
計算時間もかかるため4並列にして計算を実行します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2006 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object decomposeParDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // numberOfSubdomains 4; method scotch; // ************************************************************************* // |
scotchを使えば適当に最適化して分割してくれます。
その他の設定は特に触っていないのでそのままでいきます。
※並列計算で実行すると水面高さのデータが出力されませんでした。
原因は調査中ですので今回は並列計算は諦めて並列なしで計算実行します。
（追記：OpenFOAM v2106で改善されていましたので、v2106以降をお使いの方は並列計算でも水面高さは出力されます。）
設定が終わったので以下のコマンドで解析を実行しましょう。
1 | interFoam & |
計算実行中に圧力データの時刻歴を見ることができます。
別のターミナルを開き以下のコマンドを打ちます。
1 | foamMonitor postProcessing/interfaceHeight1/0/height.dat |
計算が終わったらParaviewで結果を確認しましょう。
このように水が崩れて障害物で水しぶきを上げている結果が得られました。
解析結果
指定した位置での水面の他高さはの時刻歴の記事同様gnuplotスクリプトで可視化します。
以下の記述を「gnuplot_run_H.plt」というファイル名で保存しgnuplotを起動してスクリプトを読み込みみます。
gnuplot_run_H.plt
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | if (exist("n")==0 || n<0) n=0 #変数の初期化 point=n+1 sample=sprintf("%d",point+1) #数値を文字列(整数型)に変換 # exp_data="exp_data/3D_damBrake_experiment_data/9449af_6ce73873d6d04315853e481c337dbfca/test_case_2_exp_data.csv" CAE_data="postProcessing/probes/0/p" # 画像の保存 set terminal png set output "graph_P".point.".png" # 画枠のマージン set lmargin 15 set tmargin 5 set rmargin 10 set bmargin 5 set mxtics #set mxtics 10 # グリッドを表示 set grid # 座標軸レンジとラベル set yrange[0:12000] set xlabel font "Arial,16" set ylabel font "Arial,16" set tics font "Arial,20" set ylabel 'Pressure(Pa)' offset -2,0 set xlabel 'Time(sec)' offset 0,-1 # グラフタイトル set title "Pressure P".point offset -5,0 font "Arial,24" plot "exp_data/3D_damBrake_experiment_data/9449af_6ce73873d6d04315853e481c337dbfca/test_case_2_exp_data.csv" us 1:@sample with lines lt 1 lw 3 title "EXP P".point, "postProcessing/probes/0/p" us 1:@sample with lines lt 2 lw 3 title "OpenFOAM P".point if (n<3) pause 0.1; n=n+1; \ reread # スクリプトの再読み込み n=-1 #end of script |
- 1行目で変数nを初期化
rereadが実行されてスクリプトの最初に戻る
ターミナル上で以下のコマンドを打ちます。
1 2 3 4 5 6 7 8 9 10 11 12 13 | $ gnuplot G N U P L O T Version 5.2 patchlevel 8 last modified 2019-12-01 Copyright (C) 1986-1993, 1998, 2004, 2007-2019 Thomas Williams, Colin Kelley and many others gnuplot home: http://www.gnuplot.info faq, bugs, etc: type "help FAQ" immediate help: type "help" (plot window: hit 'h') Terminal type is now 'qt' |
gnuplotが起動したら以下のように打ってスクリプトを実行します。
1 | gnuplot> loat "gnuplot_run_H.plt" |
各位置での水面の高さについて実験データとOpenFOAMを比較すると以下のようになります。
概ね良い結果ですね。
まとめ
今回は実験データの指定した位置での水面の高さとOpenFOAMの結果を比較を行い、概ね良い結果を得ることができました。
しかし、実は解析上は問題がありなぜか並列計算をした際にinterfaceHeight1を認識してくれませんでした。
interfaceHeightについてなぜmpirunにすると値出てこなくなるのかわかってない_(._.)_ https://t.co/fIjsKW5Cyb
— カマキリ🐲Python頑張る昆虫 (@t_kun_kamakiri) May 27, 2022
なので、仕方なく今回はシングルで計算を行ったため物凄い計算時間がかかってしまいました。
どなたか原因がわかる方いましたらご教授いただければ幸いです_(._.)_
参考書
PENGUINさんサイトを体系的に学べる書籍となっています。ネット記事でも十分勉強できるのですが、OpenFOAMの初学者でOpenFOAMをインストール済みであれば一冊持って置き、体系的に学ぶのが良いでしょう。
あとは初心者向けに丁寧に解説がされているこちらの書籍もお勧めです。最後の章にはオーバーセットメッシュ（重合メッシュ）の機能を使った解析を最後まで丁寧に解説しているので挫折することはないでしょう。
Follow @t_kun_kamakiri