こんにちは（@t_kun_kamakiri）
今回は、前回作成したメッシュを使って3次元ダムブレイクの流体解析を行い、指定した座標位置での圧力データを取得します。
そして、CAE解析結果と実験データとの比較を行います。
前回の記事では3次元ダムブレイクの障害物近辺のみメッシュの再分割を行いましたのでまだお読みでない方は下記の記事をご参考ください。
前回の記事
使用しているベースとなるチュートリアル「$FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreakWithObstacle/」よりコピーして使っています。
今回の内容
3次元ダムブレイクの実験の指定した座標位置での圧力の時刻歴をOpenFOAMの結果と比較する。
時刻歴の重ねはgnuplotを使います。
OpenFOAM初心者でチュートリアルを動かしたことがある方を対象にしています。
DEXCS2020
OpenFOAM v2006
gnuplot 5.2
3次元ダムブレイクの実験データ
こちらの実験データのより以下の時刻歴データを計測しています。
- 圧力計測点P1~P8←今回OpenFOAMと比較する対象
- 水面の高さH1~H4
本記事では実験データとの比較は圧力データのみとしますが、実験では以下のように水面の高さも計測しています(H2とH4)。
H1~H4での水面の高さは次回の記事で解説します。
解析実行
解析設定はこちらをご参考ください。
設定ファイルは再度こちらにも記載しておきます。
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で読み込みます。
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 | /*--------------------------------*- 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; // 壁面での圧力時系列モニター用のプローブ } // ************************************************************************* // |
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を使えば適当に最適化して分割してくれます。
その他の設定は特に触っていないのでそのままでいきます。
設定が終わったので以下のコマンドで解析を実行しましょう。
1 2 3 | decomposePar mpirun -np 4 interFoam -parallel reconstructPar |
計算実行中に圧力データの時刻歴を見ることができます。
別のターミナルを開き以下のコマンドを打ちます。
1 | foamMonitor postProcessing/probes/0/p |
計算が終わったらParaviewで結果を確認しましょう。
このように水が崩れて障害物で水しぶきを上げている結果が得られました。
実験データ比較(gnuplot)
こちらの実験データより圧力データを入手します。
zipファイルを解凍します。
「test_case_2_exp_data.xls」が実験データです。
gnuplotでグラフ化
gnuplotをExcelデータで読み込むことができなかったため（文字コードが違う?）ため、Excelデータを全てコピーし、新規作成のテキストファイルに貼り直してcsvファイルの名前を「test_case_2_exp_data.csv」として保存します。
treeコマンドでファイル構成を見ておきましょう。
実験データは「exp_data/3D_damBrake_experiment_data/9449af_6ce73873d6d04315853e481c337dbfca」に保存しています。
※計算結果は省略
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 | tree . ├── 0 │ ├── 0.org │ │ ├── U │ │ ├── alpha.water │ │ ├── k │ │ ├── nut │ │ ├── omega │ │ └── p_rgh │ ├── U │ ├── alpha.water │ ├── k │ ├── nut │ ├── omega │ ├── p_rgh │ └── polyMesh │ └── cellMap ├── 0.org │ ├── 0.org │ │ ├── U │ │ ├── alpha.water │ │ ├── k │ │ ├── nut │ │ ├── omega │ │ └── p_rgh │ ├── U │ ├── alpha.water │ ├── k │ ├── nut │ ├── omega │ └── p_rgh ├── Allclean ├── Allrun ├── constant │ ├── dynamicMeshDict │ ├── g │ ├── polyMesh │ │ ├── boundary │ │ ├── cellZones │ │ ├── faceZones │ │ ├── faces │ │ ├── neighbour │ │ ├── owner │ │ ├── pointZones │ │ ├── points │ │ └── sets │ │ ├── c0 │ │ ├── refineCells │ │ └── refinedCells │ ├── transportProperties │ └── turbulenceProperties ├── damBreakWithObstacle_noDynamichMesh_rasidual.OpenFOAM ├── exp_data │ └── 3D_damBrake_experiment_data │ ├── 9449af_6ce73873d6d04315853e481c337dbfca │ │ ├── test_case_2.f │ │ ├── test_case_2_exp_data.csv │ │ ├── test_case_2_exp_data.xls │ │ └── test_case_2_v1p1.pdf ├── run └── system ├── blockMeshDict ├── controlDict ├── decomposeParDict ├── fvSchemes ├── fvSolution ├── probes ├── refineMeshDict ├── setFieldsDict ├── topoSetDict1 └── topoSetDict2 16 directories, 78 files |
gnuplotを起動します。
gnuplotと打つとgnuplotが起動し対話形式になります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | $ 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> |
このようになっていればgnuplotが使えます。
以下のコマンドを打って圧力計測点P1の実験データとOpenFOAMの結果を比較してみましょう。
1 | plot "exp_data/3D_damBrake_experiment_data/9449af_6ce73873d6d04315853e481c337dbfca/test_case_2_exp_data.csv" us 1:2 with lines lt 1 lw 3 title "EXP P1", "postProcessing/probes/0/p" us 1:2 with lines lt 2 lw 3 title "OpenFOAM p1" |
さらにグラフをきれいに見せたい場合は以下のページを参考にアレンジすることができます。
あわせて読みたい
以下の記述を「gnuplot_run_P.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 | point=1 #変数 sample=sprintf("%d",point+1) #数値を文字列(整数型)に変換 # 画像の保存 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 |
スクリプトの実行は以下のコマンド行います。
1 | load "gnuplot_run_P.plt" |
以上を実行すると「graph_P1.png」という画像ファイルが作成されています。
先頭分で「point=”1″ #変数」と変数定義しています。
文字列は「.(ドット)」を使って連結できます。
例えば「set output “graph_”.sample.”.png”」にすることで「set output “graph_1.png”」となります。
では他の計測点も比較してみましょう。
全ての点で概ね良い結果ですね。
p4だけ初期のピークがあっていないですね。
gnuplotに慣れていないため今回は「point=1」の数字を１～８に変更してひとつひとつ画像作成を行いましたが、gnuplotはスクリプト言語なので繰り返し文を使えば簡単に画像作成ができるはずです。追々調べて修正します。
※追記
以下のスクリプトで各測定点の時刻歴データの画像ファイルを出力できます。
gnuplot_run_P.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が実行されてスクリプトの最初に戻る
まとめ
今回は3次元ダムブレイクとOpenFOAMの圧力データの時刻歴をgnuplotでグラフ化して比較を行いました。
概ね計算結果は実験値を再現していますね。
ただ注意点として実測データとシミュレーション結果が概ね良い結果だったとしても、映像などが結構ずれている場合もあります。今回は映像まで細かく比較できていませんが、実測データと映像の両方でシミュレーションが実験を再現できているか確認する必要がありますね。
次回は水面の高さを実験データと比較します。
こちらも実験データとOpenFOMAでの結果が概ね近しい挙動をしています。
参考書
PENGUINさんサイトを体系的に学べる書籍となっています。ネット記事でも十分勉強できるのですが、OpenFOAMの初学者でOpenFOAMをインストール済みであれば一冊持って置き、体系的に学ぶのが良いでしょう。
あとは初心者向けに丁寧に解説がされているこちらの書籍もお勧めです。最後の章にはオーバーセットメッシュ（重合メッシュ）の機能を使った解析を最後まで丁寧に解説しているので挫折することはないでしょう。
Follow @t_kun_kamakiri