こんにちは（@t_kun_kamakiri）
本記事ではOpenFOAMを用いたスロッシングタンクの流体解析を設定の解説を行います。
- OpenFOAMでスロッシングタンクの解析設定を解説
- スロッシングタンクの実験データとの比較を行う
スロッシングタンクは、タンク内に入れた液体が周期的な振動によりうねる現象のことです。構造物には構造物特有の固有振動を持っているため、その周期と重なったときに共振を起こし小さな力でも問題になることがあります。
- OpenFOAMでダイナミックメッシュをやってみたい
- 流体解析を実験データと比較したい
ベンチマークテストの試験データとOpenFOAMとの比較も行います。
OpenFOAMで再現するのはこちらの実験です。
タンクが振動して水がバシャバシャ暴れていますね。
CAEと実機との比較検討は本記事では十分に行うことができませんが、解析設定から実験データとの比較まで一連の流れを体験することができます。
OpenFOAMv2212
Python3.8.5
解析対象
まずは解析対象を確認します。
実験の元データはこちらになります。
自身が解析するフォルダにとりあえず保存しておきましょう。
フォルダの中に試験の解説が書かれているpdfがあります。
試験の寸法はこちらにきっちり書かれています。
Sensor1~6は実験での圧力センサーの位置です。
OpenFOAMでも同じ位置にprobeを設けて圧力の時刻歴変化を出力するようにします。
試験ではz方向の厚みを31mm,62mm,124mmの3パターン行っていますが、データとして62mmのみを保存しているとのことです。
振動周期は以下のように書かれています。
$T_{1}$が振動の周期[s]です。
OpenFOAMでこれ通りに与えると、そもそも試験での振動と上記の数式による振動の周期が微妙にずれていました。
ですので、ここでは実験で得られた試験データの振動周期をOpenFOAMにも与えるようにします。試験条件が違うとなると実機とCAEのコリレーションの意味がなくなるので。
H[mm]の高さが2パターンありますが、これは液体の初期水位です。
今回示すのはH=93mmですので、93mだけ見ておけば良いです。
また、液体はwater（水）です。
物性値は以下のようになっています。
Water(水)による実験データと比較するのでwaterだけを見ておけば良いです。
以上で解析対象の説明が終わりました。
以下でOpenFOAMの解析設定の解説を行います。
チュートリアルから解析ケースをコピー
OpenFOAMにスロッシング用のチュートリアルがあるのでそれをコピーします。
|
1 |
cp -r $FOAM_TUTORIALS/multiphase/interFoam/laminar/sloshingTank3D6DoF 001_run |
※「$FOAM_TUTORIALS＝/usr/lib/openfoam/openfoam2212/tutorials」
ソルバはinterFoamを使います。
気体と液体を体積分率$\alpha$で区別するVOF法による解析ソルバです。
先ほど入手した試験データは001_runと同じ階層に保存しておきます。
- 001_run
- expdata/9449af_05f4aa0eaa1d4164a1024d07c562592a/SPHERIC_TestCase10
9449af_05f4aa0eaa1d4164a1024d07c562592aは元々のフォルダ名ですので、この辺は好みで名前やフォルダ構成を変えてください。
メッシュ作成
形状は直方体なのでblockMeshで簡単に設定することができます。
形状作成とメッシュ作成、境界面の名前はblockMeshで行います。
system/blockMeshDict
|
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 |
scale 0.001; xMin -450; xMax 450; yMin 0; yMax 508; zMin 0; zMax 62; dx 10; dy 10; dz 5; xl #calc "$xMax-($xMin)"; yl #calc "$yMax-($yMin)"; zl #calc "$zMax-($zMin)"; nx #calc "$xl/$dx"; ny #calc "$yl/$dy"; nz #calc "$zl/$dz"; vertices ( ($xMin $yMin $zMin) ($xMax $yMin $zMin) ($xMax $yMax $zMin) ($xMin $yMax $zMin) ($xMin $yMin $zMax) ($xMax $yMin $zMax) ($xMax $yMax $zMax) ($xMin $yMax $zMax) ); blocks ( hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1) ); edges ( ); boundary ( walls { type wall; faces ( (3 7 6 2) (0 4 7 3) (2 6 5 1) (1 5 4 0) (0 3 2 1) (4 5 6 7) ); } ); mergePatchPairs ( ); |
設定が終わりましたら、blockMeshを実行します。
|
1 |
blockMesh |
結果をParaViewで確認しましょう。
形状はいたってシンプルですね。
境界条件の設定
乱流モデルを使わない解析ですので、設定する物理量はU(流速)、p_rgh(圧力)、体積分率alpha.waterです。
圧力は計算の便宜上$p-\rho \boldsymbol{g}\cdot h$という値で行っている点に注意です。
いつでも初期化できるように0.origフォルダ内のファイルを編集します。
0.orig/U
|
1 2 3 4 5 6 7 8 9 10 11 12 |
dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { walls { type movingWallVelocity; value uniform (0 0 0); } } |
壁条件ではありますが振動をしているのでtypeをmovingWallVelocityとしています。
※movingWallVelocity：パッチでのフラックスが0となるよう、パッチでの値の法線方向を置き換える。
0.orig/p_rgh
|
1 2 3 4 5 6 7 8 9 10 11 |
dimensions [1 -1 -2 0 0 0 0]; internalField uniform 0; boundaryField { walls { type fixedFluxPressure; } } |
fixedFluxPressureは境界のフラックスが速度の境界条件にマッチするように圧力勾配を調整するものです。
0.orig/alpha.water
|
1 2 3 4 5 6 7 8 9 10 11 |
dimensions [0 0 0 0 0 0 0]; internalField uniform 0; boundaryField { walls { type zeroGradient; } } |
typeを法線方向の勾配0のzeroGradientとしておきます。
初期状態の設定
初期状態としては水を93mm浸した状態を作るためsetFieldsで場の設定を行います。
system/setFieldsDict
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
defaultFieldValues ( volScalarFieldValue alpha.water 0 ); regions ( boxToCell { box (-0.5 -0.1 -0.1) (0.5 0.093 0.1); fieldValues ( volScalarFieldValue alpha.water 1 ); } ); |
では初期状態を確認しましょう。
まず、先ほど境界条件で0.origで設定したので0フォルダという名前にコピーします。
|
1 |
cp -r 0.orig 0 |
setFieldsにより0フォルダ内はalpha.waterの分布を持ったファイルができてしまうので、0.origを残しておくと後で境界条件を設定しなおすときに便利です。
では、setFieldsを実行します。
|
1 |
setFields |
ParaViewで結果を見てみましょう。
ParaViewのデフォルトが0フォルダを読みこまない設定になっているので、「Skip Zero Time」のチェックを外して読み込んでください。
図のように93mmに体積分率$\alpha=1$（水）が表示されていればOKです。
移動メッシュの設定
移動メッシュはconstant /dynamicMeshDictで行います。
|
1 2 3 4 5 6 7 8 |
dynamicFvMesh dynamicMotionSolverFvMesh; motionSolver solidBody; solidBodyMotionFunction tabulated6DoFMotion; CofG (0 0 0); timeDataFileName "<constant>/1DoF_exp.dat"; |
solidBodyMotionFunctionは以下のように選択ができます。
|
1 2 3 4 5 6 7 8 |
SDA axisRotationMotion drivenLinearMotion linearMotion multiMotion oscillatingLinearMotion oscillatingRotatingMotion rotatingMotion tabulated6DoFMotion |
3次元6自由度の振動を与えるためtabulated6DoFMotionとしています。
CofGは重心位置のことです。
「timeDataFileName “<constant>/1DoF_exp.dat”;」で振動のデータを与えています。
1DoF_exp.dat”には決まったフォームでデータを書かないといけないので、エクセルで処理します。
- SPHERIC_TestCase10/data_files/lateral_oil_1x.txtをExcelにタブ区切りで貼り付ける
- 以下のように入力用のフォームにする
- コピーしてconstant/1DoF_exp.datにコピーする
|
1 |
="( " & A3 &"( (0 0 0) ( 0 0 " & F3 & ") )) "） |
constant/1DoF_exp.dat
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
167001 ( ( 0 ( (0 0 0) (0 0 0.000000001143341) ) ) ( 0.00005 ( (0 0 0) (0 0 0.001194) ) ) ( 0.0001 ( (0 0 0) (0 0 0.001194) ) ) ( 0.00015 ( (0 0 0) (0 0 0.001194) ) ) ( 0.0002 ( (0 0 0) (0 0 0.001194) ) ) ( 0.00025 ( (0 0 0) (0 0 0.001194) ) ) ( 0.0003 ( (0 0 0) (0 0 0.001194) ) ) ( 0.00035 ( (0 0 0) (0 0 0.001194) ) ) ( 0.0004 ( (0 0 0) (0 0 0.001194) ) ) (省略) |
※余談ですがチュートリアルには数式で上記の振動データを作るプログラムがありました。こちらの数式に従いプログラムを作成してタンクの振動のインプットデータとして使うことができます。試しにやってみましょう。
型定義などはOpenFOAM用のC++って感じですがわりと簡単に関数を作れます。
作ったプログラムが以下です。
|
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 |
#include "List.H" #include "vector.H" #include "Vector2D.H" #include "Tuple2.H" #include "OFstream.H" #include "fvCFD.H" using namespace Foam; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: int main(int argc, char *argv[]) { // End-time of the table const scalar endTime = 8; // Number of entries in the table const label nTimes = 100; // Amplitude of the translation [m] const vector transAmp(0, 0, 0); // Frequency of the translation [rad/s] const vector transOmega(0, 0, 0); // Amplitude of the rotation [deg] const scalar rotAmpz = 4; const vector rotAmp(0, 0, rotAmpz); // Frequency of the rotation [rad/s] const scalar L = 0.9; const scalar H = 0.093; const scalar g = 9.81; const scalar T1 = constant::mathematical::twoPi*std::pow(std::pow(constant::mathematical::pi*g/L*std::tanh(constant::mathematical::pi*H/L), 0.5),-1); const scalar rotOmegaz = constant::mathematical::twoPi /T1; const vector rotOmega(0, 0, rotOmegaz); Info << "T1 = " << T1 << endl; List<Tuple2<scalar, Vector2D<vector>>> timeValues(nTimes); forAll(timeValues, i) { scalar t = (endTime*i)/(nTimes - 1); timeValues[i].first() = t; timeValues[i].second()[0] = vector ( transAmp.x()*Foam::sin(transOmega.x()*t), transAmp.y()*Foam::sin(transOmega.y()*t), transAmp.z()*Foam::sin(transOmega.z()*t) ); timeValues[i].second()[1] = vector ( rotAmp.x()*Foam::sin(rotOmega.x()*t), rotAmp.y()*Foam::sin(rotOmega.y()*t), rotAmp.z()*Foam::sin(rotOmega.z()*t) ); } { OFstream dataFile("MyGen1DoF.dat"); dataFile << timeValues << endl; } Info<< "End\n" << endl; return 0; } |
こちらをコンパイルします。
|
1 |
wmake |
同じフォルダに「MyGen1DoF.dat」が出力されて先ほどのconstant/1DoF_exp.datとして使うことができます。
こちらのインプットデータで解析をしましたが、試験とタンクの振動周期がずれてきます。
試験が理論通りにいかないことはよくあることです。今回は振動のインプットデータは試験と同じにしたいので試験データの振動をそのままOpenFOAMのデータとして使います。
もし、振動が理論通りになっているかどうかも含めて知りたい、もしくは試験データがまだないと言った場合は理論式を使わざるを得ないのです。
物性値
物性値はconstant/transportPropertiesで行います。
こちらは資料通りに設定します。
constant/transportProperties
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
phases (water air); water { transportModel Newtonian; nu 8.96e-07; rho 998.2; } air { transportModel Newtonian; nu 1.48e-05; rho 1; } sigma 0.0728; |
離散化スキーム
離散化スキームはデフォルトの設定のままでよいです。
system/fvSchemes
|
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 |
ddtSchemes { default Euler; } gradSchemes { default Gauss linear; } divSchemes { div(rhoPhi,U) Gauss vanLeerV; div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss vanLeer; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } |
対流項は2次精度のTVDスキームvanLeerを使用しています。
代数ソルバ
代数ソルバもデフォルトのままでよいです。
system/fvSolution
|
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 |
solvers { alpha.water { nAlphaCorr 1; nAlphaSubCycles 3; cAlpha 1; } "pcorr.*" { solver PCG; preconditioner { preconditioner GAMG; tolerance 1e-05; relTol 0; smoother DICGaussSeidel; cacheAgglomeration no; } tolerance 1e-05; relTol 0; maxIter 100; } p_rgh { solver GAMG; tolerance 1e-08; relTol 0.01; smoother DIC; } p_rghFinal { solver PCG; preconditioner { preconditioner GAMG; tolerance 2e-09; relTol 0; nVcycles 2; smoother DICGaussSeidel; nPreSweeps 2; } tolerance 2e-09; relTol 0; maxIter 20; } U { solver smoothSolver; smoother GaussSeidel; tolerance 1e-06; relTol 0; nSweeps 1; } } PIMPLE { momentumPredictor no; nCorrectors 2; nNonOrthogonalCorrectors 0; correctPhi no; pRefPoint (0 0 0); pRefValue 1e5; } relaxationFactors { equations { "U.*" 1; } } |
以下で基準圧を指定しています。
- pRefPoint (0 0 0);
- pRefValue 1e5;
pRefValueの値を0Paとしていると実験と同じ基準圧になります。
計算制御
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 |
application interFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 8; deltaT 0.01; writeControl adjustable; writeInterval 0.05; purgeWrite 0; writeFormat binary;//ascii; writePrecision 6; writeCompression on; timeFormat general; timePrecision 6; runTimeModifiable yes; adjustTimeStep yes; maxCo 0.5; maxAlphaCo 0.5; maxDeltaT 1; functions { probes { type probes; libs (sampling); writeControl timeStep; writeInterval 1; probeLocations ( (-0.450 0.093 0.031) //sensor1 (-0.450 0.299 0.031) //sensor2 (-0.425 0.508 0.031) //sensor3 (0.400 0.508 0.031) //sensor4 (0.450 0.483 0.031) //sensor5 (0.450 0.222 0.031) //sensor6 ); fixedLocations false; fields ( p ); } } じ |
時間刻みは「maxCo 0.5; 」のクーラン数0.5で指定しています。
ちなみに「deltaT 0.01;」は初期値にしか使っていないので（初期値にも使ってないかも）、こちらの値を小さくしても時間刻みは変わらないので注意です。
時間刻みを変えたい場合はクーラン数を小さくしましょう。
functionsで指定した位置での圧力を出力する設定を行っています。
測定点は移動する場合と移動しない場合を考える必要があります。- メッシュ移動するけど初期の固定位置で測定したい場合はyes(or true)
- メッシュ移動すると移動した分だけ測定位置も変わる場合はno
今回は容器に圧力計測を取り付けているので、fixedLocationsはnoとする必要があります。
fixedLocations |
Do not recalculate cells if mesh moves | no | true |
以上で解析の設定が終了です。
次に計算を実行します。
計算実行
interFoamコマンドを実行して計算を実行します。
※今回は並列計算無しで計算をさせますが、8秒前の解析で数十分かかるので、待ってられない方は並列計算をさせてください。
|
1 |
interFoam > log.interFoam & |
- interFoam：ソルバの実行
- log.interFoam：log.interFoamにログを出力
- &：バックグラウンドで流す
ParaViewで確認します。
CAEは水しぶきの水がぼやけてしまっていますね。
実験の動画も載せておきます。
実験と比較するとCAEは解像度から水しぶきがぼやけているため、実機との乖離があります。今回使っているVOFによる解法はメッシュが粗いと空気と水との境界で体積分率がなまってしまうという問題があるため、ここからCAEの精度を上げるためには改善が必要です。
- メッシュを細かくして解像度を上げる
ただし計算負荷がかかる - 界面を追跡しながらダイナミックメッシュを行う
毎時刻メッシュを切り直すの計算負荷がかかる - 乱流モデルを設定する
細かい渦などを見ることを諦めて乱流状態を数理モデルで表現する
こちらのように - 界面圧縮項を大きくする
- interIsoFoamソルバを使う
今度のCAE精度アップの改善活動は必要ですが、本記事では深追いはしません。
圧力データの比較
実機とCAEの圧力データの比較を行うために簡単にPythonスクリプトを書きました。
今回はsensor1のみをグラフ化してみます。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
import pandas as pd import numpy as np import matplotlib.pyplot as plt plt.figure(figsize=(16,10)) #data_p = np.loadtxt(f'./expdata/9449af_05f4aa0eaa1d4164a1024d07c562592a/SPHERIC_TestCase10/data_files/lateral_water_1x.txt') df_p1_exp = pd.read_table(f"../expdata/9449af_05f4aa0eaa1d4164a1024d07c562592a/SPHERIC_TestCase10/data_files/lateral_water_1x.txt").iloc[:,0:2]#, columns=['time','p1','p2','p3','p4','p5','p6']) df_p1_exp['Pressure[Pa]'] = df_p1_exp['Pressure[mbar]']*100 caedata_p = np.loadtxt(f'./postProcessing/probes/0/p') df_p1_cae = pd.DataFrame(caedata_p, columns=['time','CoGp', 'p1','p2','p3','p4','p5','p6']) df_p1_cae['p1(Pa)abs'] = df_p1_cae['p1']-df_p1_cae['p1'][0] plt.title('p1 pressure',fontsize=24) plt.plot(df_p1_exp.iloc[:,0],df_p1_exp.iloc[:,2],linestyle='solid',color='black',label="Exp") plt.plot(df_p1_cae.iloc[:,0],df_p1_cae.iloc[:,8],linestyle='solid',color='red',label="OpenFOAM") plt.xlabel('time(sec)',fontsize=16) plt.ylabel('Pressure(Pa)',fontsize=16) plt.legend(fontsize=24) plt.grid() plt.xticks(list(np.arange(0,8,1)),fontsize=12) plt.yticks(list(np.arange(-1000,9000,1000)),fontsize=12) |
おおまかな圧力の挙動は再現できていますね。
ちなみにポリヘドラルメッシュに変えてみましょう。
これまた何とも言えない結果ですね。
ずいぶんノイジーな感じです。
ビタビタに実験とCAEがあることは難しいので、これをどう定量的に合っていると表現するか、どこまでの実機とCAEのコリレーションレベルで妥協するのかは納期や目的にも依るでしょう。
実験データはかなりノイジーなデータなので、n増し（複数回試験をすること）しないとノイズなのか本当にそういう振動なのかが判断付かないというのもあります。このような衝撃のある圧力データなどは1回のデータで判断せずにn増し試験を追加するのが賢明です。
追試験1
どうもOpenFOAMの結果の圧力が1.5秒3.5秒・・・あたりで負圧へ振動しているのが気になっていたのですが、原因が基準圧の座標位置の問題でした。
system/fvSolution
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
（省略） PIMPLE { momentumPredictor no; nCorrectors 2; nNonOrthogonalCorrectors 0; correctPhi no; pRefPoint (-0.45 0.508 0.031); pRefValue 0; } （省略） |
以下で基準圧を指定しています。
- pRefPoint (-0.45 0.508 0.031));
- pRefValue 0;
pRefValueの値を0Paとしていると実験と同じ基準圧になります。
これで余計な振動がなくなりました。
追試験2
実験データの「Position_original [deg]」の値が重複しているので、重複を除いたデータにして再度解析を行いました。
重複データの除去は以下のように適当にPythonで行いました。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
import pandas as pd df_data = pd.read_table("../expdata/SPHERIC_TestCase10/data_files/lateral_water_1x.txt") df_data_new = pd.DataFrame() i_list = [] data_old = "" for i, data in enumerate(df_data["Position_original [deg]"]): if data != data_old: i_list.append(i) data_old = data df_data_new = df_data.iloc[i_list] df_data_new.to_excel("lateral_water_1x_new.xlsx", index=False) |
最終的には以下のようになりました。
まとめ
今回はOpenFOAMを使ってスロッシングの解析を行いました。
また、実験データとOpenFOAMの結果を比較しおおまかな圧力の挙動は示していることを確認しました。
ただ、実験のような水しぶきを再現するにはVOF法では限界がありそうであることがわかったので、今後CAE向上のためにも改善が必要です。
おすすめ参考図書
☟こちらは、OpenFOAMの日本語書籍が無い中唯一わかりやすい参考書だと思います。
☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。
☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの？」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ（NextPublishing）)
インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。
OpenFOAMコードの辞書的な扱いとしては以下の参考書が大変役に立ちます。
本記事ではESI版のOpenFOAMを使っているため本書のFoundation版で対応していない部分がありますが、その辺を考慮しても持っていて損はないでしょう。
メッシュ作成を極めたいって方はこちらもお勧めです。
OpenFOAMでメッシュ作成
Follow @t_kun_kamakiri
