高校物理、大学受験、大学物理、プログラミングの記事を書きます。

流体解析フリーソフトOpenFOAMを使ったケルビンヘルムホルツ不安定性の解析計算

2019/01/19
 
この記事を書いている人 - WRITER -

 

本記事の目標

OpenFOAMを使って、ケルビンヘルムホルツ不安定性の解析計算を行う。

 

ケルビンヘルムホルツ不安定性に興味を持って、

OpenFOAMでできるかなって思ってやってみただけですので、特にメッシュの切り方やパラメータや方程式の解き方(離散化)などはこだわってやっていません。

研究に使うとか仕事に使うとかではありませんので、遊び感覚でできました。

 

※注意

本記事は、OpenFOAM を使用していますがはじめてOpenFOAMを触る人向けには書いていません。

設定ファイルは下記から入手できます。

casefile

でもOpenFOAMを触ったことがあって、使い方はちょっとは知っているよって人は上記のファイルを入手して下記のコマンドを打てば「メッシュ作成」から「解析実行」までできると思います。

 

スポンサーリンク

ケルビンヘルムホルツとは

 

概要

異なる速度(密度)をもつ2相の間で起こる不安定性現象です。

2相の界面でせん断力があるため不安定性の種(渦)が生じ、それが時間とともに増幅することで生じる2次元的な模様のことを指します。

 

 

実現象での例

ケルビンヘルムホルツ不安定性は、こんな感じ。

飛行機の窓から外を眺めているとたまに見かける渦を巻いた雲や、木星の模様もケルビンヘルムホルツ不安定性の現象が作り出した模様です。

 

 

 

 

 

 

CAE解析

 

OpenFOAMの2相流のソルバ”interFoam”を用いてケルビンヘルムホルツ不安定性を再現したいと思います。

 

ツール

OpenFOAM:ver4.x
ソルバ:inerFoam

 

解くべき方程式

 

・運動方程式(2次元)
・体積分率の保存
・表面張力
※VOFによる界面の表現

 

今回使っている「interFoam」のソルバはVolume Of Fluid(VOF)法による、二相流ソルバーです。

相の区別を各相の体積分率(0~1)で行い、密度や粘性を体積分率をかけて、ナビエストークス方程式を解いています。

ナビエストークス方程式

\begin{align*}\frac{\partial \rho_{b}\boldsymbol{u}}{\partial t}+\nabla\cdot\big(\rho_{b}\boldsymbol{u}\boldsymbol{u}\big)=-\nabla p+\nabla\cdot\mu_{b}\big(\nabla \boldsymbol{u}+\nabla \boldsymbol{u}^{T}\big)+\rho_{b}\boldsymbol{f}+\boldsymbol{F}_{s}\end{align*}

 

密度

\begin{align*}\rho_{b}=\alpha_{1}\rho_{1}+\alpha_{2}\rho_{2}\end{align*}

 

粘性

\begin{align*}\mu_{b}=\alpha_{1}\mu_{1}+\alpha_{2}\mu_{2}\end{align*}

 

表面張力

\begin{align*}\boldsymbol{F}_{s}=\gamma \kappa\big(\nabla \alpha\big)\end{align*}

\begin{align*}\kappa=\nabla\cdot\big(\frac{\nabla\alpha}{|\nabla\alpha|}\big)\end{align*}

 

非圧縮の条件

\begin{align*}\nabla\cdot \boldsymbol{u}=0\end{align*}

 

体積分率の移流方程式

\begin{align*}\frac{\partial \alpha}{\partial t}+\nabla\cdot \big(\alpha \boldsymbol{u}\big)=0\end{align*}

 

 

※\(\alpha_{2}=1-\alpha_{2}\)

※\(\kappa\):界面の曲率

※\(\gamma\):界面張力

境界条件と解析領域

 

 

その他

・2次元解析
・乱流モデル無し
・非定常計算
・時間間隔:クーラン数1

 

メッシュ作成

 

メッシュ作成はOpenFOAMのユーティリティにあるblockMeshのみで作成しています。

設定方法:blockMeshDict

 

ざっくりこんな感じで切ってみました。

界面のまわりだけちょっとだけメッシュを細かくしています。

どれくらい細かくすれば良いかは、わからないのでメッシュサイズにこだわり過ぎず、まずは荒くメッシュを切って計算がうまくいっていることを確認してからメッシュサイズを見直しても良いかもしれません。

 

初期状態の設定

初期状態の設定は、setFiledsを使います。

 

設定方法:setFiedDict

 

初期状態の速度分布として、上半分は1.0m/s、下半分は1.5m/sとします。

初期状態の体積分率として、上半分と下半分で相を分けています

 

解析結果

 

それっぽい模様ができました!(^^)!

動画速度がとてもゆっくりなので、だいぶ待ってから渦が生成され始めます。

感想

 

OpenFOAM自体はフリーで誰でも自宅PCでてきる上に、ソースコード(C++)も公開されているので、やる気になればカスタマイズもできます。

ソースコードが公開されているのは、とても勉強になります。

3次元の解析にすると、また面白いかなと思います。

おそらく奥行き方向にも渦が発生して、乱流っぽくなるんでしょうか。

そういえばOpenFOAMに「DNS(直接ナビエストークス方程式を解く)」があります。
これはスペクトル法を用いた解法で、離散化手法の中では最も精度が高いとされています。

(高波数領域のエイリアシングエラー?とかは注意が必要ですが)

いつか触ってみようかと思います。

>OpenFOAMの「dnsFoam」のドキュメント

 

おまけ

 

YouTubeでケルビンヘルムホルツ不安定の解説を実験も踏まえて易しく解説してくれています。

この記事を書いている人 - WRITER -

コメントを残す

Copyright© 宇宙に入ったカマキリ , 2018 All Rights Reserved.