Python

【Pythonプログラム】衝撃波菅の理論とOpenFOAMの結果比較

こんにちは(@t_kun_kamakiri

OpenFOAMの衝撃波菅内の解析と厳密解との比較を行います。
本記事では衝撃波管内の厳密解についてのプログラムを載せておきます。

衝撃波管に関して関係する内容を以下に記載しておきます。
OpenFOAMでは、rhoCentralFoamとsonicFoamソルバにshock tubeという名前のチュートリアルがあります。
rhoCentralFoamの離散化スキームにKurganovとTadmorが実装されています。

コードは自分のメモ用としてPythonでべた書きしています。
後日OpenFOAMのチュートリアルと比較するためのものです。

コードの実行にはjupyter labで順番に実行していきます。

衝撃波官

  • 1次元波動のショックチューブ(衝撃波管)
  • 20世紀後半に研究したゲーリー A. ソッド(Gary A. Sod)にちなんでsod shock tubeと呼ばれる
  • 1次元の非粘性圧縮方程式(オイラー方程式)と理想気体を仮定したの厳密解があり、数値流体の精度や安定性の評価に使われる
  • OpenFOAMにはrhoCentralFoamとsonicFoamにチュートリアルがある
  • 非粘性の圧縮流れの方程式(1次元オイラー方程式)で記述
  • 衝撃波前後でランキン・ユゴニオの関係式と理想気体の断熱過程の式より厳密解が求まる

必要なライブラリと初期状態

まずは適当に初期状態を与えます。

反復計算用の関数

以下は$\frac{p_{2}}{p_{1}}=p_{21}$を求めるための関数です。
式で書くと

\begin{align*}
\frac{p_{2}}{p_{1}}=\frac{p_{R}}{p_{L}}\bigg(\,\,\,1+ \frac{\gamma-1}{2c_{L}}\bigg(u_{L}-u_{R}-\frac{c_{R}(\frac{p_{2}}{p_{1}}-1)}{\gamma\sqrt{\frac{\gamma+1}{2\gamma}}(\frac{p_{2}}{p_{1}}-1)+1}\bigg)\,\,\,\bigg)^{\frac{2\gamma}{\gamma-1}}
\end{align*}
と少々複雑です。

厳密解を解くための関数

こちらが厳密解を解くためのプログラム部分です。

  1. 領域1 $ x>x_{0}+V_{s}t$
    $\rho_{1}=\rho_{R}$, $p_{1}=p_{R}$, $u_{1}=u_{R}$
    衝撃波速度$ V_{s}=u_{1}+c_{1}\sqrt{\frac{\gamma + 1}{2\gamma}\big(\frac{p_{2}}{p_{1}-1}+1\big)}$
  2. 領域2 $ x_{0}+V_{c}t < x \leq x_{0}+V_{s}t$
    衝撃波の前後でランキン・ユゴニオの関係式より、$\frac{\rho_{2}}{\rho_{R}}=\frac{\frac{p_{2}}{p_{R}}+\frac{\gamma-1}{\gamma+1}}{\frac{\gamma-1}{\gamma+1}\frac{p_{2}}{p_{R}}+1}$
    $u_{2}=u_{R}+\frac{c_{R}\big(\frac{p_{2}}{p_{R}}-1\big)}{\sqrt{\frac{\gamma}{2}\big(\gamma-1+\frac{p_{2}}{p_{R}}(\gamma+1)\big)}}$
  3. 領域3 $x_{0}+V_{rt}t < x \leq x_{0}+V_{c}t$ $(V_{rt}=u_{2}-\frac{\gamma-1}{2}(u_{L}-u_{2})+c_{L})=u_{3}-c_{3}$
    等エントロピー過程よりポアソンの関係式を使って、$p_{3}\rho_{3}^{\gamma}=p_{L}\rho_{L}^{\gamma}$,
    流速と圧力は変化がないので$u_{3}=u_{2}=V_{c}$, $p_{3}=p_{2}$
    $c_{3}=\sqrt{\frac{\gamma p_{3}}{\rho_{3}}}=\frac{\gamma-1}{2}(u_{L}-u_{3})+c_{L}$
  4. 領域4 $x_{0}+V_{rh}t < x \leq x_{0}+V_{rt}t$
    流速$u_{4}=\frac{2}{\gamma+1}\big(\frac{x-x_{0}}{t}+c_{L}+\frac{\gamma-1}{2}u_{L}\big)$
    音速$c_{4}=c_{L}-\frac{\gamma-1}{2}(u_{4}-u_{L})$
    等エントロピー過程よりポアソンの関係式を使って、$p_{4}\rho_{4}^{\gamma}=p_{L}\rho_{L}^{\gamma}$ or $p_{4}c_{4}^{\frac{2\gamma}{\gamma-1}}=p_{L}c_{L}^{\frac{2\gamma}{\gamma-1}}$
  5. 領域5 $x\leq x_{0}+V_{rh}t$
    $p_{5}=p_{L}$, $u_{5}=u_{L}$, $\rho_{5}=\rho_{L}$

まずは$\frac{p_{2}}{p_{1}}=p_{21}$を求めるための反復計算を行っています。

\begin{align*}
\frac{p_{2}}{p_{1}}=\frac{p_{R}}{p_{L}}\bigg(\,\,\,1+ \frac{\gamma-1}{2c_{L}}\bigg(u_{L}-u_{R}-\frac{c_{R}(\frac{p_{2}}{p_{1}}-1)}{\gamma\sqrt{\frac{\gamma+1}{2\gamma}}(\frac{p_{2}}{p_{1}}-1)+1}\bigg)\,\,\,\bigg)^{\frac{2\gamma}{\gamma-1}}
\end{align*}
は「F(p21, p1, p5, rho1, rho5, u1, u5)」の関数から求めています。

※理想気体の状態方程式から、温度と密度と圧力は$p=\rho R_{MW} T$
ただし、$R_{MW}$は気体定数$R=8.31$[J/k mol]とモル質量$MW=28.9$[mol](空気のモル質量)から、$R_{MW}=\frac{R}{MW}$で計算される値であることに注意

時間発展の計算

座標に対する物理量をリストに格納しています。
さらに、各時刻の物理場(リスト)をリストに入れています。

アニメーションの作成

まずは密度のアニメーション

続いて流速のアニメーション

続いて圧力のアニメーション

続いて温度のアニメーション

各時刻のcsvデータを出力

ついでに各時刻のcsvデータも出力しておきます。

以下のようにcsvファイルが出力されます。

  • time0.001.csv
  • time0.001.csv
  • time0.002.csv
  • time0.003.csv
  • time0.004.csv
  • time0.005.csv
  • time0.006.csv
  • time0.007.csv

time0.007.csv

OpenFOAMの結果と比較

OpenFOAMの衝撃波菅はshock tubeというチュートリアルとして用意されています。
ソルバはsonicFoamとrhoCentralFoamがあり、どちらも超音速近辺での圧縮性のソルバです。

今回はrhoCentralFoamのshock tubeのチュートリアルの結果と上で示した理論解との比較を行います。

モデルを少し編集しています。
必要とする流速はx方向のみであるためmag(U)を出力しておきます。

system/controlDict

物理量の取得のためにsampleDictを用意します。

system/sampleDict

ここから適当にPythonで厳密解とOpwnFOAMを比較するためのアニメーションを作成します。密度を出力するプログラムのみを示しますが、流速と圧力も同様にアニメーション作成ができます。

以下は自身のフォルダ構成に合わせて変更

  • datafilename:OpenFOAMの結果ファイル名
  • filenamepath :相対パス

以下は密度を出力するプログラムです。
同様に流速と圧力も出力してみます。

OpenFOAMのrhoCentralFoamのチュートリアルはデフォルトでかなり振動していますね。
厳密解と比較するとよくわかります。
この振動を抑えるためには、

  • クーラン数を小さくする(時間刻みを小さくする)
  • メッシュを細かくする
  • 安定する離散化スキームを使う

などが挙げられます。

例えばクーラン数を0.09にすると以下のように改善されます。

こちらに解析ケースとParaViewでのグラフ作成、gnuplotでのグラフ作成のデータをアップしました。

参考書

圧縮性流体についてはこちらの参考書がおすすめです。

圧縮性流体力学(第2版): 内部流れの理論と解析

圧縮性流体力学(第2版): 内部流れの理論と解析

松尾 一泰
3,960円(07/21 22:00時点)
Amazonの情報を掲載しています

本記事の理論解はこちらの参考書のプログラムを参考にしています。
C言語(タイトルにC++とあるが・・)で書かれたプログラムをPythonに書き換えたのが本記事で内容です。

数値流体解析の基礎 - Visual C++とgnuplotによる圧縮性・非圧縮性流体解析 -

数値流体解析の基礎 – Visual C++とgnuplotによる圧縮性・非圧縮性流体解析 –

鋒, 肖, 孝夫, 長崎
3,960円(07/22 18:00時点)
Amazonの情報を掲載しています

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

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

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

3,520円(07/22 17:01時点)
Amazonの情報を掲載しています

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

OpenFOAMプログラミング

OpenFOAMプログラミング

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

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

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

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

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

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

関連記事もどうぞ

COMMENT

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