compressibleInterDyMFoamを利用した噴流シミュレーション

2,641 views
Skip to first unread message

Nakayama

unread,
Oct 18, 2012, 1:20:04 PM10/18/12
to open...@googlegroups.com
はじめて投稿させていただきます、Nakayamaと申します。
今年の8月から、海外の大学院で修士論文としてOpenFOAMを利用したliquid Jetのシミュレーションに取り組んでいます。

しかし、残念ながら、指導教官及び、同僚にOpenFOAM経験者がいないため、インターネットが教官となっています。
そのような環境の中で、日本語で情報を得られる本グループには、いつも大変お世話になっています。

さて、研究の対象は小型の噴流(医療用注射器を想定)で、以下の2つの目的があります。
1.    同僚の実験結果を模擬する事
2.    ピストン、オリフィス(針の部分)の直径や長さを変更するパラメトリックデザイン

ソルバー選定、メッシュ生成、移動格子設定、初期/境界条件等の各設定は終了して計算を実行しているのですが、途中で計算が発散してしまい、研究を進める事ができません。

以下に計算条件、質問、補足を記述しましたので、心当たりのある問題点があれば、ご指摘いただけますでしょうか。

計算条件:
使用Version: 2.1.1
ソルバー: compressibleInterDyMFoam
乱流モデル:k-epsilon
メッシュ生成: blockMesh (+checkMeshで確認)
移動格子: LayerAdditionRemoval (movingConeTopoを変更)

質問:
0.25msあたりで計算が発散してしまいます。
原因は、以下のいずれかと考えています。
ケースフォルダを添付しましたので、もし心当たりのある設定値がありましたら、ご教授いただければ幸いです。
(そもそも、ソルバーや計算ドメインの設定に問題があるといったご指摘でも構いません。)

1.    計算スキーム(fvScheme, fvSolution)
2.    初期/境界条件
3.    物性値設定 (transportproperties)

なお、液体(水)のpsi(圧縮率)を空気と同等にすると、発散せずに計算を進められます。
ただ、実験結果と比較して定性的に違う現象が起こります。(噴流が空気中で崩壊)
オリフィス内の水を押し出す力(ピストン内の圧力)が足りないのだと思います。
ちなみに移動壁の速度は、実験の平均値を用いています。
実際は、水によるダンピングの効果が加わるため、この速度は時間の関数となるよう、変更予定です。

添付ファイル
liquidJet: case folder

可視化図:time = 0- 0.2ms
計算ドメイン(単位は、mm)
左: ピストン(他端は、moving wall, constant velocity 3 m/s)
中央: オリフィス
右:大気
上図:volume fraction of water (red: water, blue: air)
下図:magnitude of velocity

liquidJet.zip
liquidJetRoccoReal9_alpha1_magU.0000.jpg
liquidJetRoccoReal9_alpha1_magU.0001.jpg
liquidJetRoccoReal9_alpha1_magU.0002.jpg

ohbuchi

unread,
Oct 21, 2012, 6:50:58 AM10/21/12
to open...@googlegroups.com
こんばんは。
dynamicMesh問題で収束性が悪い場合、最も疑わしいのはメッシュ移動ライブラリの適切性だと
思いますが、liquidJetTopoFvMeshは標準ライブラリには見当たりません。オリジナルのライブラリ
だと思いますが、これも同時に示して頂けないと問題が再現できないため、アドバイスのしようが
ないと思います。

コンタ図を見るとalphaのコンタは妥当に見えますが、Uのコンタがいかにも不自然です。
p_rghのコンタを出力してノズル出口の圧力がどうなっているか確認して下さい。
恐らく負圧になっていて、atmosphere境界から逆流が生じているためにノズル出口の流速が
直角に曲げられて上向きに早い速度になっているのではないでしょうか?
右側境界を流出境界にして、上面境界を十分に離した上で圧力規定境界にすると良いのでは
ないかと思います。
移動メッシュの詳細が不明のため、定かではありませんが、ご参考まで。



2012年10月19日金曜日 2時20分04秒 UTC+9 Nakayama:

Nakayama

unread,
Oct 21, 2012, 5:45:06 PM10/21/12
to open...@googlegroups.com
ohbuchi 様

ご回答ありがとうございました。

>dynamicMesh問題で収束性が悪い場合、最も疑わしいのはメッシュ移動ライブラリの適切性だと思いますが、liquidJetTopoFvMeshは標準ライブラリには見当たりません。
>オリジナルのライブラリだと思いますが、これも同時に示して頂けないと問題が再現できないため、アドバイスのしようが
ないと思います。

ご指摘ありがとうございました。
liquidJetTopoFvMeshを添付いたします。

ご参考までに、以下のようにして作成しました。
1.movingConeTopoをダウンロード
  http://www.openfoam.org/mantisbt/view.php?id=494

2.curLeft, left, leftExtrusionFacesに関するプログラムをコメントアウト

3.curRight, right, rightExtrusionFacesの対象座標を x = -0.01m (ピストンの可動部:movingWallの座標)に変更

4. vertexMarkup (「格子移動の有無」識別フラグ)を変更:ピストン部とそれ以外(オリフィス、大気中)の領域に識別



コンタ図を見るとalphaのコンタは妥当に見えますが、Uのコンタがいかにも不自然です。
p_rghのコンタを出力してノズル出口の圧力がどうなっているか確認して下さい。
恐らく負圧になっていて、atmosphere境界から逆流が生じているためにノズル出口の流速が
直角に曲げられて上向きに早い速度になっているのではないでしょうか?
右側境界を流出境界にして、上面境界を十分に離した上で圧力規定境界にすると良いのでは
ないかと思います。

p_rgh及び、UYを確認しました。
添付図:iquidJetRoccoReal9_UY_p_rgh.0000-0002
少ない情報にもかかわらず、詳細なアドバイスをありがとうございました。
早速試してみます。

ただ、コメントの冒頭で「最も疑わしいのはメッシュ移動ライブラリの適切性」とご指摘いただきましたように、
そもそも移動メッシュがソルバーに悪さをしていると、計算領域の変更だけではすまないかもしれないと危惧しております。

補足:
先日、試しに水の圧縮率を空気と同等にして計算を進めてみました。
すると、ピストンの速度を大きくして計算時間を進めていくと、添付図のようにピストンの可動部でalpha1の分布がおかしくなっていく現象が見受けられました。
添付図:iquidJetRoccoReal6_alpha1_magU.0004 (ピストン速度10m/s, time = 0.4 ms)

この現象は、ピストンの速度を小さくしても、計算が進むと同様に現れます。

試しにmoving meshのライブラリを movingConeTopoベースから libfvMotionSolvers.soに変更して、計算しても同様の結果となりました。
(movingConeTopoのダウンロード元に、バグらしき記述があったので、もしやと思い、別のライブラリも試してみました)
(meshは、軸対称メッシュを使用していますが、純2次元計算でも同じ結果になるか試してみようと考えています)


グループの皆様

移動メッシュとinterFoam関連ソルバーの併用で、上記のような現象を経験した方がいらっしゃれば、コメントいただければ幸いです。
また、OpenFoam経験、3ヶ月弱の初心者ですが、お役に立ちそうな情報がありましたら、投稿させていただきます。
今後もよろしくお願いいたします。
liquidJetTopoFvMesh.zip
liquidJetRoccoReal9_UY_p_rgh.0000.jpg
liquidJetRoccoReal9_UY_p_rgh.0001.jpg
liquidJetRoccoReal9_UY_p_rgh.0002.jpg
liquidJetRoccoReal6_alpha1_magU.0004.jpg

ohbuchi

unread,
Oct 22, 2012, 2:52:54 AM10/22/12
to open...@googlegroups.com
こんにちは。
liquidJetTopoFvMeshですが、私のOpenFOAM-2.1.1ではコンパイルは出来ましたが、計算実行時に
printStackエラーが出ます。

movingWall を inlet境界にしてU,k,epsilonをfixedValue、p_rghをzeroGradientに変えて
compressibleInterFoamで計算したところ、それらしい結果になりました(乱流モデルは
標準kEpsilonに変えました)。
しかし、ジェット先端から巻き上がる渦芯でk、epsilonが著しく増大する現象が見られました。

もしかしたら、乱流モデルの収束に問題があって渦粘性が過大評価され、Uの結果が変になった
かも知れません。k、epsilon、nutのコンタを確認してみて下さい。

以上、ご参考まで。


2012年10月22日月曜日 6時45分06秒 UTC+9 Nakayama:

Nakayama

unread,
Oct 22, 2012, 4:13:31 PM10/22/12
to open...@googlegroups.com
ohbuchi 様

>liquidJetTopoFvMeshですが、私のOpenFOAM-2.1.1ではコンパイルは出来ましたが、計算実行時にprintStackエラーが出ます。

申し訳ありませんでした。
計算を始める前にconstant / *Zones ファイルを削除する必要があります。

movingConeTopoでは、移動faceを特定しfaceZoneとmeshModifierを作成するのですが、
constantディレクトリ以下に既にZone関連ファイルやmeshModifierがあると、それらを読み込もうとします。

今回は、blockMesh実行時に作成された*Zonesファイルを読み込もうしますが、
movingConeTopoで作成されるものと内容が違うためエラーになります。

補足事項
使用しているソルバーですが、
myCompressibleInterDyMFoamとなっていますが、中身はオリジナルと変わりありません。
(変更する可能性があったので、別名にしました)

また、KEpsilonDyM (k-epsilon乱流モデル)ですが、インターネット上でダウンロードしました。
(申し訳ありません、リンク先を探しましたが、見つけられませんでした)
オリジナル版は、(トポロジー変更を伴う)移動メッシュと併用すると問題があると書かれていましたが
私の計算結果では、大きな違いがあるようには思えませんでした。


>もしかしたら、乱流モデルの収束に問題があって渦粘性が過大評価され、Uの結果が変になったかも知れません。
>k、epsilon、nutのコンタを確認してみて下さい。

ご指摘ありがとうございます。
確かにオリフィスの出口に近い壁面のみに非常に大きな値が現れます(k=10^4, epsilon=10^10のオーダー)
(ちなみにnutは、10^-2のオーダーでした)
添付図をご参照ください。

このような場合は、epsilonやkに対する残差の設定や緩和値を変更する以外は、収束性を改善する方法はないのでしょうか。

試しにLESに変更して計算してみたのですが、同じような結果になってしまいました。
計算条件は、チュートリアルのmultiphase/interFoam/les/nozzleFlow2Dを参考にしました。

上記のチュートリアルを私の計算条件に近づけてみて、問題がない場合は、dynamic meshそのものか、ソルバーとの相性に問題があるのだと思います。(例えば、トポロジーのエラーや、ピストン内部の圧縮過程の解析など)

先日、投稿いたしました不可解なalpha1分布の件と共に、何か進展があれば、ご報告させていただきます。


2012年10月22日月曜日 15時52分54秒 UTC+9 ohbuchi:
liquidJetRoccoReal9_epsilon_k.0001.jpg
liquidJetRoccoReal9_epsilon_k.0002.jpg

ohbuchi

unread,
Oct 23, 2012, 4:56:55 PM10/23/12
to open...@googlegroups.com
本件とあまり関係ないことかも知れませんが、ご参考まで。

元々compressibleInterFoamは、等温の圧縮性自由表面流れソルバで、なぜか乱流モデルは非圧
縮のものを使っています。
状態方程式はrho=psi * pという式(psi=1/RT)で、密度が負にならない様に最小圧力Pminを指定す
る必要があります。

3ヶ月ほど前にGitサイトで大きな更新があり、一般的な状態方程式を使う形式に変更されました。transportPropertiesでのphase定義の方法が変更され、T変数の定義が必須になっています。

もしかしたら今回の様な強い圧縮性が作用する流れでは、最新ソルバの方が良いかも知れません。


2012年10月23日火曜日 5時13分31秒 UTC+9 Nakayama:

Nakayama

unread,
Oct 24, 2012, 4:59:59 PM10/24/12
to open...@googlegroups.com
ohbuchi 様


> 3ヶ月ほど前にGitサイトで大きな更新があり、一般的な状態方程式を使う形式に変更されました。
>transportPropertiesでのphase定義の方法が変更され、T変数の定義が必須になっています。
> もしかしたら今回の様な強い圧縮性が作用する流れでは、最新ソルバの方が良いかも知れません。

貴重な情報、ありがとうございます。
Gitサイトから最新版をダウンロードする方法は初めてですが、説明を読んでトライしてみます。

また、移動メッシュのプログラムにバグがあるのではないか、と考え始めています。
(以下の検証計算を行いました)

添付図(sloshingTank_alpha1.0010)は、シンプルなsloshing tank(純2D)の計算結果です。
U境界条件: 上部:PressureInletOutlet、右壁:slip、下壁:non slip
P境界条件: 上部:total pressure, 上部以外:buoyantPressure
alpha1境界条件: 上部: inletOutlet, inlet value/value 0, 上部以外:zerogradient

移動メッシュとcompressibleInterDyMFoamの動作確認のために行ないました。
しかし、液相の中に空気相が混ざってしまう計算結果になってしまいました。

計算条件
ソルバー: compressibleInerDyMFoam
乱流モデル: 標準k-epsilon
メッシュ:純2次元メッシュ(z軸面:empty)
移動メッシュ: movingCone改造(左壁 moving wall、3m/s, layerAdditionRemoval使用)

ただし、以下の条件では、上記のような現状はみられませんでした。

ソルバー: interDyMFoam
乱流モデル: 標準k-epsilon
メッシュ:純2次元メッシュ(z軸面:empty)
移動メッシュ: movingCone改造(左壁 moving wall、3m/s, layerAdditionRemoval使用)
fvSolutionに以下の一文を追加。
PIMPLE
{
    correctPhi      yes;
}

上記の一文を追加すると、添付図(sloshingTank2D1_alpha1.0010)及び、計算ログ内で以下の変化が起こります。
対応するプログラムは、correctPhi.H (バグ修正用のパッチプログラム?)だと思います。

補足事項
-上記の一文なしでは、interDyMFoamのログで、Max(alpha1) が1以上の値となります。
-上記の一文は、compressibleInterDyMFoamを使用した場合、可視化結果は変わらず、おかしいままでした。

バグが本当にあるのかどうかは、わかりませんが、ソルバーとともに最新版を試してみたいと思います。
(バグレポートには、未解決のままのmoving mesh関連のバグが結構あるのが、気がかりですが・・・)

---計算ログ抜粋---
time:9.90099e-08 curMotionVel_:(3 0 0) curRight:-0.011
No topology change
Executing mesh motion
Execution time for mesh.update() = 0.17 s
time step continuity errors : sum local = 0, global = 0, cumulative = 0
GAMGPCG:  Solving for pcorr, Initial residual = 1, Final residual = 6.60522e-06, No Iterations 16
time step continuity errors : sum local = 2.70036e-05, global = 2.70034e-05, cumulative = 2.70034e-05
MULES: Solving for alpha1
Phase-1 volume fraction = 0.500003  Min(alpha1) = 0  Max(alpha1) = 1
MULES: Solving for alpha1
Phase-1 volume fraction = 0.500007  Min(alpha1) = 0  Max(alpha1) = 1
MULES: Solving for alpha1
Phase-1 volume fraction = 0.50001  Min(alpha1) = 0  Max(alpha1) = 1
MULES: Solving for alpha1
Phase-1 volume fraction = 0.500013  Min(alpha1) = 0  Max(alpha1) = 1


2012年10月24日水曜日 5時56分55秒 UTC+9 ohbuchi:
sloshingTank_alpha1.0010.jpg
sloshingTank2D1_alpha1.0010.jpg

ohbuchi

unread,
Oct 24, 2012, 8:38:48 PM10/24/12
to open...@googlegroups.com
こんにちは。
確かに変な結果ですね。Mantisを見るとtopoChangerFvMeshに関するバグレポート0000494が4月に
報告されていながら半年以上も何のリアクションもなしに放置されていますね。

私はsliding meshは良く使うのですが、大きな問題は感じていません。しかし、移動メッシュ全般となると
確かに、正規版は弱いという印象を持っています。
移動メッシュをメインに使うのであれば、Extend ProjectのOpenFOAM-1.6-extの方が、例題も多く
様々な移動メッシュクラスが利用できるのでおすすめです。
http://www.extend-project.de/

ご参考まで。


2012年10月25日木曜日 5時59分59秒 UTC+9 Nakayama:

ohbuchi

unread,
Oct 24, 2012, 9:37:24 PM10/24/12
to open...@googlegroups.com
少し調べてみました。
correctPhi.Hは移動メッシュの場合に、境界パッチのU,phiを補正する処理と、
圧力補正方程式を解いてphiを修正する処理を含んでいます。
correctNonOrthogonalが有効の時だけ圧力補正を行う様です。
メッシュ移動が生じるとセル体積が変化するのでphiの補正を行わないと
alphaが0~1に収まらなくなるのではないかと思います。
問題は、圧縮性流れの場合にphiが正しく補正されるかですね。

correctPhi.Hのphi補正処理は以下の通りです。
ぱっと見では問題なさそうなのですが。。。

***** compressibleInterDyMFoamのcorrectPhi.H *********

dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);

adjustPhi(phi, U, pcorr); ←圧力境界がないときに圧力基準を決める処理?

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU  ←divUが圧縮性補正?
);

pcorrEqn.solve();

if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();   ←修正圧力からフラックスを求めPhiを補正
}
}

***** pimpleDyMFoam ********

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rAU, pcorr) == fvc::div(phi)  ←divU項(=0)なし
);

pcorrEqn.setReference(pRefCell, pRefValue);  ←圧力基準設定(adjustPhi代わり?)
pcorrEqn.solve();

if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();  ←修正圧力からフラックスを求めPhiを補正
}
}





2012年10月25日木曜日 9時38分48秒 UTC+9 ohbuchi:

Nakayama

unread,
Oct 26, 2012, 4:19:55 PM10/26/12
to open...@googlegroups.com
ohbuchi 様

correctPhi.Hについて、ご説明ありがとうございました。
解説いただきましたとおり、

fvSolution:
PIMPLE
{
correctPhi yes;
correctNonOrthogonal yes;
}

として、以下の条件で再計算してみました。
しかし、数time Stepsのうちに、alpha1が発散してしまい計算が続けられませんでした。

ソルバー:compressibleInterDyMFoam
乱流モデル: LES (RANS: k-epsilonは、オリフィス付近でepsilonが増大し、Uが発散しました)
メッシュ: 軸対称+ 移動メッシュ(movingConeTopo改造)

残念ですが、version 2.1.1での計算は諦める事にします。
Gitサイトから最新版をダウンロードして、状態方程式を使うcompressibleInterDyMFoamにトライします。
また、ご紹介いただきましたext版もダウンロードしたので、movingConeTopoよりも最適なlibを探してみます。

貴重なアドバイスを多数いただき、ありがとうございました。


2012年10月25日木曜日 10時37分25秒 UTC+9 ohbuchi:

Nakayama

unread,
Oct 27, 2012, 6:35:30 PM10/27/12
to open...@googlegroups.com

早速、OpenFOAM-2.1.xをダウンロードして、"最"最新版のcompressibleInterDyMFoamを実行してみました。
しかし、計算を始めると以下のエラーがでて、最初のtime stepで計算がストップしてしまいます。

おそらく、具体的なケースの話ではなく、ソルバーそのもの(ライブラリ?)のエラーだと思うのですが、
どなたか、以下のエラーの内容及び、修正方法をご存知の方がいましたら、ご連絡いただけないでしょうか?

公式サイトのbugレポートで、以下のキーワードを入れて探してみましたが、残念ながら該当するバグは見つかりませんでした。
"compressibleInterDyMFoam", "libOpenFOAM.so"

もう少し調べてみて、私のインプットファイルの設定ミスやコンパイルミスでなければ、バグレポートに上げようと思います。

---以下、エラー内容---
Starting time loop

Courant Number mean: 0 max: 0
deltaT = 9.90099e-08
Time = 9.9009901e-08


time:9.90099e-08 curMotionVel_:(3 0 0) curRight:-0.01
No topology change
Executing mesh motion
Execution time for mesh.update() = 0.16 s
time step continuity errors : sum local = 2.97039e-05, global =
-2.97039e-05, cumulative = -2.97039e-05

GAMGPCG:  Solving for pcorr, Initial residual = 1, Final residual =
8.08114e-06, No Iterations 16
GAMGPCG:  Solving for pcorr, Initial residual = 1.00758e-08, Final
residual = 1.00758e-08, No Iterations 0
time step continuity errors : sum local = 2.40041e-10, global =
1.76098e-11, cumulative = -2.97038e-05
--> FOAM Warning :
    From function cubeRootVolDelta::calcDelta()
    in file cubeRootVolDelta/
cubeRootVolDelta.C at line 52
    Case is 2D, LES is not strictly applicable

DILUPBiCG:  Solving for k, Initial residual = 1, Final residual =
1.01666e-15, No Iterations 2
MULES: Solving for alphawater
Liquid phase volume fraction = 0.500005  Min(alpha1) = 0  Min(alpha2)
= -2.22045e-16
MULES: Solving for alphawater
Liquid phase volume fraction = 0.50001  Min(alpha1) = 0  Min(alpha2) =
-2.22045e-16
MULES: Solving for alphawater
Liquid phase volume fraction = 0.500015  Min(alpha1) = 0  Min(alpha2)
= -2.22045e-16
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0,
No Iterations 0
DILUPBiCG:  Solving for T, Initial residual = 0.512853, Final residual
= 6.89042e-09, No Iterations 68
#0  Foam::error::printStack(Foam::Ostream&) in
"/home/haruka/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::sigFpe::sigHandler(int) in
"/home/haruka/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::divide(Foam::Field<double>&, double const&,
Foam::UList<double> const&) in
"/home/haruka/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#4
 in "/home/haruka/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/compressibleInterDyMFoam"
#5  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#6
 in "/home/haruka/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/compressibleInterDyMFoam"
Floating point exception (core dumped)


2012年10月27日土曜日 5時19分55秒 UTC+9 Nakayama:

ohbuchi

unread,
Oct 28, 2012, 6:18:55 AM10/28/12
to open...@googlegroups.com
こんばんは。
発散時のメッセージを見ると、TEqnの計算のあと、pressure corrector loopで
ゼロ割による浮動小数点例外が発生している様です。
pEqn.Hでゼロ割の発生しそうな箇所は
     volScalarField rAU("rAU", 1.0/UEqn.A());
        solve
        (
            (
                (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
              + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
            )
          + p_rghEqnIncomp,
            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
        );

        if (pimple.finalNonOrthogonalIter())
        {
            dgdt =
            (
                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
            );

            phi = phiHbyA + p_rghEqnIncomp.flux();

            U = HbyA
              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
            U.correctBoundaryConditions();
        }
    }
のあたりだと思われます。どこかまではわかりませんが。

私の環境で、compressibleInterFoam/les/depthCharge2Dに何もしないdynamicMeshDict
を作成してcompressibleInterDyMFoamで計算したところ、compressibleInterFoamと
同じ結果になりましたので、ソルバのバグということでもなさそうです。
設定ファイルのどこかに問題があってゼロ割が生じているのではないでしょうか。



2012年10月28日日曜日 7時35分30秒 UTC+9 Nakayama:

Nakayama

unread,
Oct 28, 2012, 2:03:26 PM10/28/12
to open...@googlegroups.com
詳細に調べていただき、ありがとうございました。
「ゼロ割」をキーワードにp_rghの内の初期値を0から1気圧(101,325Pa)変更したら、計算が流れました。

ここで、本題からははずれてしまうのですが、圧力の初期/境界条件について、お聞きしたい事があります。

1. p_rghの意味
以下のCFD Online
http://www.cfd-online.com/Forums/openfoam/80454-p_rgh-1-7-a.html
によれば、p_rghは、静圧からhydrostatic pressure (静水圧?)を引いた値とあります。

---抜粋---
p_rgh = p - rho*gh;
So it is the pressure without the hydrostatic pressure and is intialized from the pressure field in the p-file.
In the pEqn.H file, the pressure equation is written and solved for the p_rgh, so the boundary conditions important for the pressure solution are the p_rgh conditions. After the pressure solution, the p is calculated with:
p = p_rgh + rho*gh;
---

また、以下の本グループの投稿によれば、p_rghは、ゲージ圧力との表現もあります。
https://groups.google.com/forum/?fromgroups=#!msg/openfoam/_Bi06fDJBhg/fr3oTlnrUtsJ

この2つの情報を参考に、今回の計算では(静圧)1気圧下での計算であったため、
初期値を0に設定したのですが、今回は計算がながれませんでした。
(version2.1.1のcompressibleInterDyMFoamでは、計算できました)

何か私の理解(設定)が間違っていたのでしょうか?
それとも最新版のcompressibleInterDyMFoamでは、状態方程式を解いている事が何か関係しているのでしょうか?


2. buoyantPressure
「fixedGradientベースの境界条件で、密度勾配値から圧力勾配を計算する」と理解しています。

以下を参照しました。
http://www.openfoam.org/docs/user/damBreak.php
http://www.cfd-online.com/Forums/openfoam/67155-interfoam-mass-conservation.html

これらをもとに、今回は以下に設定しました。

boundaryField
{
    movingWall
    {
        type            buoyantPressure;
    }
}

この条件と、下記に違いはあるのでしょうか?
boundaryField
{
    movingWall
    {
        type            buoyantPressure;
    uniform     0; (or 101325;)
    }
}

また、私の計算は、VOF(水と空気の2相)を利用しており、密度はこの2相の線形結合(平均)から計算されるので、
圧力p_rghの壁境界には、すべてbuoyantPressureを使用しました。
しかし、チュートリアル:multiphase/interFoam/les/nozzleFlow2Dでは、壁面にzeroGradientを使用しています。
これは、何か理由があるのでしょうか?


2012年10月28日日曜日 19時18分55秒 UTC+9 ohbuchi:

Nakayama

unread,
Oct 28, 2012, 2:45:14 PM10/28/12
to open...@googlegroups.com
compressibleInterDyMFoam(v2.1.x)+movingConeTopoFvMeshk(改造)による
sloshing tankとliquid jetの計算結果をご報告します。

残念ながら、v2.1.1のcompressibleInterDyMFoamと同じく、物理的な結果を得る事ができませんでした。
内容
・水相の中に空気相が混じる
・alpha1の値が、0-1にならない

よって、やはりmovingConeTopoFvMeshにバグがあると判断いたしました。
そこで、OpenFOAM-1.6-extのmovingConeTopoチュートリアルを調べたところ、
dynamicFvMeshとして、movingBodyTopoFvMeshが使用されていました。

ちなみにソースコードには、作者のコメントがありました。
(おそらくバグの修正の事だと思います)

  // In order to do a correct update on a mask on processor boundaries,
    // Detection of moving cells should use patchNeighbourField for
    // processor (not coupled!) boundaries.  This is done by expanding
    // a moving cell set into a field and making sure that processor patch
    // points move in sync.  Not done at the moment, probably best to do
    // using parallel update of pointFields.  HJ, 19/Feb/2011

そこで、compressibleInterDyMFoam(v2.1.x)+movingBodyTopoFvMeshを使用して、
今回のケースを実行したいと考えています。

この場合、最も簡便なコンパイル、及び実行方法は、どのような方法でしょうか。
といいますのも、上記のソルバーとライブラリは、バージョン違い(v2.1.xと1.6-ext)のため、
どちらかのプロジェクトフォルダに、もう一方のソースコードをコピーして、コンパイルする必要があります。

ためしに、ライブラリをv2.1.xのプロジェクトディレクトリにコピーして、コンパイルしました。
しかし、リンクしているその他のプログラムをコピーしていなかったので、コンパイルエラーとなりました。
この場合、やはり地道にすべてのリンク先を調べて、そのリンク先にあるプログラムもすべてコピーするしか方法はないのでしょうか?

ちなみに今のディレクトリ構造は、以下のようになっています。
HOME/OpenFOAM/
OpenFOAM-1.6-ext, OpenFOAM-2.1.x, (USER NAME)-2.1.1, (USER NAME)-2.1.x


2012年10月29日月曜日 3時03分26秒 UTC+9 Nakayama:

ohbuchi

unread,
Oct 28, 2012, 9:55:23 PM10/28/12
to open...@googlegroups.com
こんにちは。
p_rghはおっしゃり通りの定義です。絶対値から静水圧を除いた値になります。
非圧縮の場合には圧力の絶対値そのものは意味を持たず、圧力勾配のみが流れに影響しますから
基準となる値を決め、そこからの相対値で設定できますが、圧縮性の場合には圧力絶対値そのものが
重要なので0はダメだと思います。

buoyantPressureは静水圧分布と整合する壁面圧力条件です。重力がなければzeroGradientと同じ
になります。nozzleFlow2Dチュートリアルでは、gが0なので違いはないと思います。


2012年10月29日月曜日 3時03分26秒 UTC+9 Nakayama:

Nakayama

unread,
Nov 5, 2012, 5:12:30 PM11/5/12
to open...@googlegroups.com
ohbuchi 様

返信が遅れて申し訳ありませんでした。

p_rghの定義について、ご回答いただき、ありがとうございました。

重力加速度が0の場合は、確かに密度差は圧力計算に影響しないですね。
チュートリアルのケースを詳細に見るべきでした。
ご指摘ありがとうございました。


さて、一週間ほど以下のケースについて取り組んだのですが、物理的な結果が全く得られませんでした。
何度もご教授いただきたのに、よい結果をご報告できず、申し訳ありません。

スロッシングタンク2D
ケース名:sloshingTank2D2
ソルバー:compressibleInterDyMFoam (2.1.x)
移動格子:motion solver
格子:純2次元格子
乱流モデル:RANS(k-epsilon), LES
備考:PIMPLEの移動格子補正あり(correctPhi, correctNonOrthogonal)

extバージョンもトライしましたが、うまくいきませんでした。
といいますのも、1.6-extのcompressibleInterDyMFoamは、圧力計算がp_rghではなく、pのままだったためか、計算が発散してしまいました。
また、1.6-extのdynamic mesh関連のlibを2.1.xに移植しようとしましたが、コンパイルエラーが盛り沢山で、使用できませんでした。(リンク先のファイルもすべて移植するか、make/optionのリンク先を変更する必要があるのだと思いました)

最後に、以下の2つの方法を試してみて、dynamic meshとVOFのコンビネーションにバグがあるのか、
それとも、私の設定ミスなのか判断しようと思います。

1.dynamic meshの代替:a piston wave maker BCという方法(標準lib? swak4Foam?)
   https://www.youtube.com/watch?v=MLXCw_9kSyY

2.liquid jetをcompressibleInterFoam(2.1.x)で計算できるか
  dynamic meshを使わずに、(ohbuchi 様にも計算いただいた)moving wallをinletに変えてみます
     最新版のソルバーで、RANS, LESの両ケースを試そうと思います。
    (RANSは、2.1.1版では、オリフィス領域で速度が発散してしまいます)


なにかわかりましたら、再度ご報告させていただきます。
全く望みがなさそうであれば、残念ですが、OpenFOAMからfluentに変更する予定です。


2012年10月29日月曜日 10時55分23秒 UTC+9 ohbuchi:
sloshingTank2D2.0000.jpg
sloshingTank2D2.0014.jpg
sloshingTank2D2.zip
Reply all
Reply to author
Forward
0 new messages