OpenFOAM

【OpenFOAM】粒子追跡を理論式と比較してみた

こんにちは(@t_kun_kamakiri

本記事はOpenFOAMでの粒子追跡の備忘録です。
自分で理解したことや試したことをまとめていきます。
前回の記事では「icoUncoupledKinematicParcelFoamを使って粒子追跡の解析」をとりあえず試しに使ってみました。

正直、粒子についての設定ファイルの中身はよくわかっていませんが、今回は1つの粒子に対して質点の運動の理論と比較してみます。
本記事の内容
  • icoUncoupledKinematicParcelFoamを使って粒子追跡の解析
  • 流体は一定の速度、粒子は球体形状とし流体の力を受けて運動する様子をシミュレーション
  • 質点の運動理論とOpenFOAMの結果を比較する

球体にはたらく重力と抗力による質点の運動とOpenFOAMのicoUncoupledKinematicParcelFoamを使って得られた粒子の運動を比較を行うことでOpenFOAMの中身を理解しようという試みです。OpenFOAMの結果は厳密には壁との摩擦もありますし質点ではないので、質点と考えた場合の運動とずれが生じます。

OpenFOAMv2012(EDI版)
WSL2

粒子の設定

モデル作成は前回の記事で解説を行っています。

粒子の設定は「constant/kinematicCloudProperties」で行います。
constant/kinematicCloudProperties

前回の記事での設定を以下に変えました。


これで壁との反発係数が1になったのかな?

理論式

粒子には以下の力がはたらきます。

  • 抗力
  • 重力
  • 浮力
  • 地面との反発力
  • 地面との摩擦力
    ※揚力は今回はシミュレーションに入れていません。

抗力について

\begin{align*}
\boldsymbol{F_{D}}&=\frac{1}{2}C_{d}\rho_{f}A|\boldsymbol{u}_{f}-\boldsymbol{v}_{d}|^2\\
&=m_{p}\frac{3}{4}\frac{\mu_{f}C_{d}Re_{p}}{\rho_{p}d_{p}^2}|\boldsymbol{u}_{f}-\boldsymbol{v}_{d}|\tag{1}
\end{align*}

  •  $\mu_{f}=1e^{-05}$ : 流体の粘性係数
  •  $\boldsymbol{u}_{f}$:流体の速度
  •  $\boldsymbol{v}_{d}$:粒子の速度
  •  $d_{p}=10$[mm] $=0.01$[m]
  •  $\rho_{d} = 964$[kg/m3] : 粒子の密度
  •  $Re_{p}=\frac{|\boldsymbol{u}-\boldsymbol{v}|d_{p}}{\mu_{f}/\rho}$ : レイノルズ数
  •  $F_{D}=\frac{1}{2}C_{d}\rho_{f}A|\boldsymbol{u}_{f}-\boldsymbol{v}_{d}|^2=m_{p}\frac{3}{4}\frac{\mu_{f}C_{d}Re_{p}}{\rho_{p}d_{p}^2}|\boldsymbol{u}_{f}-\boldsymbol{v}_{d}|$:抗力
  •  $A=\pi\big(\frac{d_{p}}{2}\big)^2$

SphereDragForce.C

球体の抗力はレイノルズ数と関係性があります。
※以下は以前にOpenFOAMで比較したものです。

レイノルズ数は$R_{e}=1000$で場合分けしています。

抗力は以下の理論式で記述しています。

\begin{align*}
C_{d}=\left\{\begin{matrix}
\frac{24}{Re}\big(1.0 + \frac{1}{6}Re^{2/3}\big)&(Re\leq 1000)\\
0.424 &(Re>1000)
\end{matrix}\right.\tag{2}
\end{align*}

重力と浮力について

\begin{align*}
\boldsymbol{F_{g}}&=m_{d}\boldsymbol{g}-\rho_{f}\boldsymbol{g}V_{d}\\
&=m_{d}\boldsymbol{g}\bigg(1-\frac{\rho_{f}}{\rho_{d}}\bigg)\tag{3}
\end{align*}

重力と浮力はひとつにまとめることができます。

GravityForce.C

重力を定義し、


地面との反発力については以下StandardWallInteraction.Cで確認すれば良いですが難しそうだったので理解するのは保留にしました。

OpenFOAMと理論式を比較

  • $y$方向は重力と浮力
  • $x$方向は抗力

$y$方向について(重力と浮力)

$y$方向は重力(+浮力)と地面と反発した際の反発係数を適当に調整して以下のように記述できます。

\begin{align*}
y=-\frac{1}{2}g(t-t_n)^2+e^{n}v_{y}(t-t_n)+h\tag{4}
\end{align*}

地面と反発するまでは自由落下なので両者は一致していますが、地面と接触するとずれてきます。
理論式でe=0.83にするとそれっぽくなりました。
OpenFOAMで反発係数e=1としているのに、どこか理解が足りていないのか?
※あとで示しますがこの設定では粒子と地面との反発係数は「collisionModel pairCollision;」となっている場合は、ソフトスフィアモデルでpairCollisionCoeffsの設定がどうも効いているようでした。

$x$方向について(抗力)

$x$方向はレイノルズ数が$Re=2000$くらいなので抗力係数が$C_{d}=0.424$として良いでしょう。

また、流体速度が$u=2$[m/s]で粒子の速度は1秒までの計算で$v=0~0.1$の間なので、レイノルズ数はほぼ流体速度の一定値$u=2$[m/s]で決まると言ってよいでしょう。
ゆえに抗力も定数になるので、運動方程式から以下のように$x$方向の式が立ちます。

\begin{align*}
x = \frac{1}{2}\frac{F_D}{m}t^2+x_{0}\tag{5}
\end{align*}
地面と接するまでは一致していますが、やはり地面と接触して以降は徐々にずれてきます。
これは理論式が地面との摩擦力を考慮していないからだろう?と考えていますが・・・・自分の理解が進んだらアップデートします。
モデルはこちらにありますので、ご参考ください。

モデルを修正

地面との反発係数が正しいのかを再度確認するために、流体からの抗力の設定を無くし重力のみの設定に変更します。

また地面との衝突は「 patchInteractionModel standardWallInteraction;」を効かすために「collisionModel none;」としておきます。

これで地面との反発係数が1になる計算ができたようです。

ただし地面との距離が球体の直径に応じて反発していると思ったのですがどうも関係ないようです・・・詳細はソースコードStandardWallInteraction.Cを確認すれば良いかもしれません。軽く確認したところ粒子の直径に応じて地面との接触を感知しているようには見えませんでした。あくまで粒子は質点として跳ね返り位置を決めているのか?だとすると粒子の跳ね返り位置はメッシュサイズに依存する気がします。

ちなみに粒子の直径を変えても粒子の座標の軌跡は変わりませんでした。

あとは浮力の影響かもしれません。OpenFOAMでgravityを設定すると$\boldsymbol{F_{g}}=m_{d}\boldsymbol{g}-\rho_{f}\boldsymbol{g}V_{d}=m_{d}\boldsymbol{g}\bigg(1-\frac{\rho_{f}}{\rho_{d}}\bigg)$のように重力と浮力一緒にして計算しているようです。
あと気になるのは、OpenFOAMの結果は跳ね返り位置がだんだん高くなっています。出力した座標値を拾っているからなんでしょうか。
理論計算は「ymin = 0.001」で地面から跳ね返るという計算にして重ねてみました。

まとめ

本記事はicoUncoupledKinematicParcelFoamソルバを使った粒子追跡について解説しました。

球体にはたらく重力と抗力による質点の運動とOpenFOAMのicoUncoupledKinematicParcelFoamを使って得られた粒子の運動を比較してみましたが、自分の理解が足りていないので一致しませんでしたが、OpenFOAMの結果は厳密には壁との摩擦もありますし質点ではないので、質点と考えた場合の運動とずれが生じます。

シミュレーションだけしていては、内部でどのような計算が行われているかユーザはわからないのですが、本記事のような基本的な理解からはじめると少しずつ理解していきます。
理解が進んだらアップデートしたいですね。

おすすめ参考図書

☟こちらは、OpenFOAMの日本語書籍が無い中唯一わかりやすい参考書だと思います。

OpenFOAMによる熱移動と流れの数値解析(第2版)

OpenFOAMによる熱移動と流れの数値解析(第2版)

3,520円(10/30 04:54時点)
Amazonの情報を掲載しています

☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。

OpenFOAMプログラミング

OpenFOAMプログラミング

Mari´c,Tomislav, H¨opken,Jens, Mooney,Kyle
8,250円(10/30 11:26時点)
Amazonの情報を掲載しています

☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの?」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。

改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ(NextPublishing))

改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ(NextPublishing))

川畑 真一
2,640円(10/30 15:46時点)
発売日: 2022/04/15
Amazonの情報を掲載しています

インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。

OpenFOAMコードの辞書的な扱いとしては以下の参考書が大変役に立ちます。
本記事ではESI版のOpenFOAMを使っているため本書のFoundation版で対応していない部分がありますが、その辺を考慮しても持っていて損はないでしょう。

OpenFOAMライブラリリファレンス

OpenFOAMライブラリリファレンス

株式会社テラバイト, 人見 大輔
16,500円(10/30 16:56時点)
Amazonの情報を掲載しています
関連記事もどうぞ

COMMENT

メールアドレスが公開されることはありません。 が付いている欄は必須項目です