流体力学の不安定性現象のひとつであるレイリー・テイラー不安定性の流体解析をOpenFOAMを使ってやってみました。
他に不安定性現象として有名なものとしてケルビンヘルムホルツ不安定性というのがありますが、そちらは以前記事にまとめましたのでご参考ください。
参考にした記事は以下です。
作り方として2つあります。
- 液体領域をstlファイルで作成しておきsetFiledsで液体領域を作成
- 0/alpha.waterのinternalField（内部領域）の領域をプログラムで埋め込む
今回の記事ではどちらの方法についても解説を行います。
OpenFOAMv2012(WSL)
チュートリアルのコピー
チュートリアルをコピーしてモデルを作りましょう。
まずは気液混相流ソルバのinterFoamのdamBreakのチュートリアルを自身の計算フォルダにコピーします。
|
1 |
cp -r $FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreak/damBreak 20230326_Rayleigh-Taylor_instability |
ベースのファイルを作っておきます。
今回2つの方法で液体領域を作成するので共通部分を「base」というフォルダに作ります。
|
1 |
mkdir 20230326_Rayleigh-Taylor_instability/base |
cd ..そしてフォルダ移動をします。
|
1 |
cd 20230326_Rayleigh-Taylor_instability/base |
damBreakのチュートリアルのメッシュと境界条件を変更して解析を行うことにします。
ベースメッシュの作成
ベースのメッシュはsystem/blockMeshDictで行います。
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 63 64 65 66 67 68 69 70 71 72 73 |
scale 1; vertices ( (-0.5 -2 -0.1) //0 (0.5 -2 -0.1) //1 (0.5 2 -0.1) //2 (-0.5 2 -0.1) //3 (-0.5 -2 0.1) //4 (0.5 -2 0.1) //5 (0.5 2 0.1) //6 (-0.5 2 0.1) //7 ); blocks ( hex (0 1 2 3 4 5 6 7) (200 800 1) simpleGrading (1 1 1) ); edges ( ); boundary ( topWall { type wall; // type patch; faces ( (3 7 6 2) ); } bottomWall { type wall; faces ( (1 5 4 0) ); } leftWall { type symmetryPlane; faces ( (0 4 7 3) ); } rightWall { type symmetryPlane; faces ( (2 6 5 1) ); } frontAndBack { type empty; faces ( (0 3 2 1) (4 5 6 7) ); } ); mergePatchPairs ( ); |
blockMeshを実行します。
|
1 |
blockMesh |
空の.foamファイルを作成してParaViewで読み込んでメッシュ状態を確認しておきましょう。
|
1 |
touch post.foam |
これでベースのメッシュができました。
次に境界条件を作成します。
物性値
物性値は「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 1e-06; rho 1000; } air { transportModel Newtonian; nu 1.48e-05; rho 1; } sigma 0.07; |
境界条件
damBreakのチュートリアルは初期状態の設定ファイルを元に戻せるように「0.orig」ファイルで境界状態を作成します。
あとで、setFieldsをすると内部領域の値が書き換わるので元に戻したい場合に「0.orig」で設定ファイルを作っておくと便利です。
0.orig/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 |
internalField uniform (0 0 0); boundaryField { leftWall { type symmetryPlane; // type noSlip; } rightWall { type symmetryPlane; // type noSlip; } bottomWall { type noSlip; } topWall { type noSlip; } backAndFront { type empty; } } |
0.orig/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 |
boundaryField { leftWall { type symmetryPlane; // type fixedFluxPressure; // value uniform 0; } rightWall { type symmetryPlane; // type fixedFluxPressure; // value uniform 0; } bottomWall { //type symmetryPlane; type fixedFluxPressure; value uniform 0; } topWall { //type symmetryPlane; type fixedFluxPressure; value uniform 0; } backAndFront { type empty; } } |
ここまでは共通部分ですね。
ここから液体領域の初期状態を作る必要があります。
以下液体領域を作る方法を2つ紹介します。
- 液体領域をstlファイルで作成しておきsetFiledsで液体領域を作成
- 0/alpha.waterのinternalField（内部領域）の領域をプログラムで埋め込む
液体領域を作る(alpha.water)
液体領域をstlファイルで作成しておきsetFiledsで液体領域を作成
まずは液体領域をstlファイルで作る方法を試します。
共通部分のファイルをコピーして新たに作ったフォルダで設定をしていきます。
|
1 2 3 |
cd .. cp -r base case001_stl cd base case001_stl |
modelというフォルダを作ってそちらに液体部分のモデルを保存しましょう。
|
1 |
mkdir model |
stlファイルはFreeCADで作ります。
今回は青色の領域をFreeCADで作りました（灰色部分はblockMeshで作った部分をFreeCADでも作成し、全体の位置関係の基準にしています。）
界面を拡大すると以下のように寸法を決めました。
液体領域をモデリングで決めることができるのはこの設定のメリットですね。
ではこれを「Mesh Design」>「メッシュ(M)」>「シェイプからメッシュを作成」として作成されたメッシュをエクスポートします。以下のようにするとstlファイルがbinary形式で保存されます。
ものモデルを「model.stl」という名前で「model」というフォルダに保存しましょう。
ではこのmodel.stlファイルを使ってsystem/setFiledsから読み込んで液体領域を作ります。
system/setFileds
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
defaultFieldValues ( volScalarFieldValue alpha.water 0 ); regions ( surfaceToCell { file "model/model.stl"; useSurfaceOrientation false; outsidePoints ((0 -0.1 0)); // definition of outside includeCut true; // cells cut by surface includeInside true; // cells not on outside of surf includeOutside false; // cells on outside of surf nearDistance -1; // cells with centre near surf (set to -1 if not used) curvature 0.9; // cells within nearDistance and near surf curvature (set to -100 if not used) fieldValues ( volScalarFieldValue alpha.water 1 ); } ); |
境界条件「0/alpha.water」は以下のようにします。
0/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 |
dimensions [0 0 0 0 0 0 0]; internalField uniform 0; boundaryField { leftWall { type symmetryPlane; } rightWall { type symmetryPlane; } bottomWall { type zeroGradient; } topWall { type inletOutlet; inletValue uniform 0; value uniform 0; } backAndFront { type empty; } } |
では、「0.orig」を「0」としてコピーして、
|
1 |
cp -r 0.orig 0 |
以下のコマンドで液体領域を作成します。
|
1 |
setFields |
こんな感じでできました。
今回のような境界条件では圧力の固定値がどこにもないため参照する圧力の固定値が必要なので、fvSolutionで指定します。
system/fvSolution
|
1 2 3 4 5 6 7 8 9 |
PIMPLE { momentumPredictor no; nOuterCorrectors 1; nCorrectors 3; nNonOrthogonalCorrectors 0; pRefCell 0; //Add pRefValue 0; //Add } |
では、計算させてみましょう。
|
1 |
interFoam |
なんかあまりきれいにできませんでしたね・・・
界面の形がカクカクしているからでしょうか。
もうちょっと界面を滑らかにして再度試してみました。
先ほどよりは滑らかな界面にしてみました。
これもあまり変わらない結果でした。
次はこちら・・・
こちらもあまりきれいじゃないですね。
0/alpha.waterのinternalField（内部領域）の領域をプログラムで埋め込む
次に流体領域を「0/alpha.water」に直接コードを埋め込んで設定したいと思います。
0/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 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 |
dimensions [0 0 0 0 0 0 0]; internalField #codeStream { codeInclude #{ #include "fvCFD.H" #}; codeOptions #{ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude #}; codeLibs #{ -lmeshTools \ -lfiniteVolume #}; code #{ const IOdictionary& d = static_cast<const IOdictionary&>(dict); const fvMesh& mesh = refCast<const fvMesh>(d.db()); scalarField alpha(mesh.nCells(), 0.); forAll(alpha, i) { const scalar x = mesh.C()[i][0]; const scalar y = mesh.C()[i][1]; alpha[i] = 0.5 + 0.5*tanh((y+0.1*cos(2*constant::mathematical::pi*x))/(0.007*1.414)); } alpha.writeEntry("", os); #}; }; boundaryField { leftWall { type symmetryPlane; } rightWall { type symmetryPlane; } bottomWall { type zeroGradient; } topWall { type inletOutlet; inletValue uniform 0; value uniform 0; } backAndFront { type empty; } } |
コードの埋め込みに関してはまた別の記事で解説をします。
このようにinternalFieldに直接値を代入することでalpha.waterの分布ができるので、setFieldsで設定しなくても良くなります。
同様にinterFoamで計算させてParaViewで可視化します。
界面が不安定になる現象ができましたね。
さきほどのstlファイルから界面を作るよりもきれいな結果になりました。
物性値をの変更
物性値を変えてみました。
constant/transportProperties
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
phases (water air); water { transportModel Newtonian; // nu 1e-06; // rho 1000; nu 6e-03; rho 3; } air { transportModel Newtonian; // nu 1.48e-05; // rho 1; nu 2e-03; rho 1; } // sigma 0.07; sigma 0.000333; |
物性値って結構影響するんですね。
このように界面が変化する様子の違いを作る因子をパラスタで調査したくなりますね。
それは追々やるとして本記事はこの辺で終わります。
まとめ
今回は流体力学の不安定性現象のひとつであるレイリー・テイラー不安定性をOpenFOAMでやってみました。
ケルビンヘルムホルツ不安定と並んで面白い現象ですよね。
液体と気体の界面の初期状態が曲線になるように設定するには以下の2つの方法があります。
- 液体領域をstlファイルで作成しておきsetFiledsで液体領域を作成
- 0/alpha.waterのinternalField（内部領域）の領域をプログラムで埋め込む
今回は試していませんが、interFoamは界面がぼやけてしまうという問題があり、いくつかの改善点が必要です。
- 界面圧縮項を大きくする
- メッシュを細かくする
- Adaptive mesh refinement（界面のメッシュのみ細かく）
- interIsoFoamソルバを使う
interIsoFoam
気液混相流に対してisoAdvector位相分率ベースの界面捕捉アプローチを用いたinterFoamから派生したソルバーで、適応的再メッシュを含むメッシュの動きやメッシュトポロジー変更を任意に行うことができます。
interIsoFoamを使うと界面のぼやけるのがinterFoamと比べるとかなり改善されています。
おすすめ参考図書
☟こちらは、OpenFOAMの日本語書籍が無い中唯一わかりやすい参考書だと思います。
☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。
☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの？」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ（NextPublishing）)
インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。
Follow @t_kun_kamakiri
