GW前半 京都旅行（１）
このGW前半（4/27〜4/29）に2泊３日で京都に旅行に行ってきた。
まず、初日（4/27）から。
京都駅に13:30頃に到着して、河原町のホテルに移動し14:00 頃にホテルにチェックイン。
チェックイン後は、四条大橋や祇園周辺を散策しながら、三十三間堂まで歩いて、そのまま三十三間堂を見学した。
千体もの仏像が並んでいる様相はやはり圧巻であった。
思い出してみると、三十三間堂は高校の修学旅行以来なので、20年以上ぶりになる。その時も鮮明に覚えているが、それだけ歴史的な重みや手入れがちゃんと行き届いているからであろう。 残念ながら、三十三間堂は内部は撮影できず。
夜は、祇園の裏道に入ったところの和食居酒屋で夕食を取った。京都らしく味付けは上品であるが、出汁がちゃんとしていて、美味しかった。
gmshインストールメモ
gmshインストールメモ
これは、gmshをインストールした際のメモである。
インストール環境
インストーラーの取得&インストール
gsmhの公式ホームページ中のmacOS(ARM)をダウンロード。
gmsh公式ホームページ
https://gmsh.info
2024年5月3日現在、安定版最新バージョンは、4.12.2である。
ダウンロードしたインストーラー「gmsh-4.12.2-MacOSARM.dmg」をダブルクリックしてインストール。
下図のようなウィンドウが現れる。「Gmsh.app」をApplicationsフォルダに移動する。
Gmshの起動
アプリケーション中の「Gmsh.app」アイコンをダブルクリックすると、以下のようにgmshが起動する。
Gmshのインストール（Homebrew経由）
Homebrew経由でもgmshをインストール可能である。 ターミナルを開いて、以下のコマンドを入力する。
brew install gmsh
which gmshを入力すると、以下のフォルダにインストールされているようである。
/opt/homebrew/bin/gmsh
こちらのインストールでは、CUIも使えるので、コマンドを使用しての高度な利用や自動化においては、こちらの方が良いかもしれない。
simpleFoam流体力履歴のプロット
概要・背景
物体周りの流れの解析を定常解析で実施する場合、特に形状が複雑な場合、目標残差に収束せず計算が何ステップにわたり継続し、計算がなかなか終了できない状況も起こりうる。
そのような場合、残差だけではなく、評価したい物理量の履歴をプロットして、その変動が落ち着いている（変動幅が一定範囲内にあるか、周期的であるか等）ことを確認して計算の終了を判断することも一つの手段である。
評価する物理量としては、解析する対象により様々であるが、物体周りの流れ解析の場合は、例えば以下のものが代表的であろう。
- 物体に作用する流体力（圧力、粘性力、圧力と粘性力の合成）
- 抗力
- 揚力
- モーメント
ここでは、OpenFOAMのsimpleFomaソルバー（定常乱流場解析）で出力される物体に作用する流体力（圧力、粘性力、圧力と粘性力の合成）の履歴をプロットするPythonスクリプトを作成したメモである。
使用環境
使用環境は以下になる。
- OpenFOAM ver2312
- Python 3.10
OpenFOAMにおける設定
simpleFoamの計算で物体に作用する流体力を出力するには、functionObjectの機能を用いる。具体的には、/system/controlDictファイルに以下の出力指示を追加する。設定項目の詳細説明は追々調査して記事にしてみようと思う。
なお、ここではmotorbikeのチュートリアルで実施した。
functions
{
forces1
{
type forces;
libs (forces);
writeControl timeStep;
timeInterval 1;
patches (motorBikeGroup);
rho rhoInf;
log true;
rhoInf 1;
CofR (0.72 0 0);
pitchAxis (0 1 0);
}
}
流体力の出力
流体力は以下のとおり、/postProcessing/ディレクトリに出力される。
- 流体力の履歴 /postProcessing/forces1/0/force.dat
流体力のプロット
上記で出力された流体力の履歴データファイルを読み込み、プロットする。 以下、プロットするのに作成したPythonスクリプトを示す。
import pandas as pd import numpy as np import matplotlib.pyplot as plt import os # Forceダータの読み込み def read_force_dat(): df_force = np.loadtxt("./postProcessing/forces1/0/force.dat") df_force = pd.DataFrame(df_force, columns=['Time','total_x','total_y','total_z','pressure_x','pressure_y','pressure_z','viscous_x','viscous_y','viscous_z']) return df_force # 読み込んだForceデータのプロット def force_plot(df): x = ['Time'] y1 = ['total_x','total_y','total_z'] y2 = ['pressure_x','pressure_y','pressure_z'] y3 = ['viscous_x','viscous_y','viscous_z'] for i in range(len(y1)): plt.plot(df[x[0]], df[y1[i]],label=y1[i]) plt.grid() plt.legend(loc='best', bbox_to_anchor=(1, 1)) plt.xlabel(x[0]) plt.ylabel('force_total') plt.savefig("force_total.png") plt.show() for i in range(len(y2)): plt.plot(df[x[0]], df[y2[i]],label=y2[i]) plt.grid() plt.legend(loc='best', bbox_to_anchor=(1, 1)) plt.xlabel(x[0]) plt.ylabel('force_pressure') plt.savefig("force_pressure.png") plt.show() for i in range(len(y3)): plt.plot(df[x[0]], df[y3[i]],label=y3[i]) plt.grid() plt.legend(loc='best', bbox_to_anchor=(1, 1)) plt.xlabel(x[0]) plt.ylabel('force_viscous') plt.savefig("force_viscous.png") plt.show() # メイン処理 df = read_force_dat() force_plot(df)
以下に上記Pythonスクリプトでプロットした図を示す。
流体力の履歴としては以下の３種類出力される。
- トータル（圧力+粘性力）
- 圧力
- 粘性力
履歴のプロットから分かるように、流体力は圧力の寄与分が支配的である。
simpleFoamの残差プロット
概要
本記事は、OpenFOAMのsimpleFomaソルバー（定常乱流場解析）の残差をプロットするPythonスクリプトを作成したメモである。
使用環境
使用環境は以下になる。
- OpenFOAM ver2312
- Python 3.10
OpenFOAMにおける設定
simpleFOAMの計算において、logファイルに残差が出力されるが、logファイルだけでは残差を抽出するには、少々フォーマットが煩雑で扱いにくい。幸い、OpemFOAMでは残差だけを抽出するような設定が準備されている。以下のように、/system/controlDictファイルに残差を出力するように命令を追記する。
continuityError1は連続式の残差（保存則）、residualsは速度、乱流諸量の残差の出力設定である。
functions
{
continuityError1
{
type continuityError;
libs
(
fieldFunctionObjects
); // Mandatory enties (unmodifiable) // Optional entries (runtime modifiable)
phi phi; //Optional (inrerited) entries
writePrecision 8;
writeToFile yes;
useUserTime yes;
region region0;
enabled yes;
log yes;
timeStart 0;
timeEnd 500;
executeControl timeStep;
executeInterval 1;
writeControl timeStep;
writeInterval 50;
}
residuals
{
type solverInfo;
libs ("libutilityFunctionObjects.so");
fields (U p k omega);
}
}
残差の出力
残差は以下のとおり、/postProcessing/ディレクトリに出力される。
- 連続式の残差 /postProcessing/continuityError1/0/continuityError.dat
- 流速、乱流諸量の残差 /postProcessing/residuals/0/solverInfo.dat
残差のプロット
上記で出力された残差データファイルを読み込み、プロットする。 以下、プロットするのに作成したPythonスクリプトを示す。
流速および乱流諸量の残差をプロットするPythonスクリプト例。
import pandas as pd import matplotlib.pyplot as plt import os # 残差ファイルの読み込み def read_residual(): df_residual = pd.read_table('./postProcessing/residuals/0/solverInfo.dat', skiprows=1) df_residual = pd.DataFrame(df_residual) return df_residual #print(df_residual.columns) # 残差プロット def residual_plot(df): init_res = [data for data in df.columns if "initial" in data] final_res = [data for data in df.columns if "final" in data] # initial residualプロット df.plot(x="# Time ", y=init_res) plt.grid() plt.legend(loc='best', bbox_to_anchor=(1, 1)) plt.yscale('log') plt.xlabel('time') plt.ylabel('initial residual') plt.show() # final residualプロット df.plot(x="# Time ", y=final_res) plt.grid() plt.legend(loc='best', bbox_to_anchor=(1, 1)) plt.yscale('log') plt.xlabel('time') plt.ylabel('final residual') plt.show() df = read_residual() residual_plot(df)
連続の式の残差（保存則）をプロットするPythonスクリプト例。
連続の式の残差データファイルは、区切りスペースが統一されていないため、まずnumpyのloadtxt関数で読み込み、numpy配列とした後にpandasでデータフレーム化している。
import pandas as pd import numpy as np import matplotlib.pyplot as plt import os # 連続の式の残差ファイルを読み込む def read_continuity_error(): df_continuity = np.loadtxt("./postProcessing/continuityError1/0/continuityError.dat") df_continuity = pd.DataFrame(df_continuity, columns=['Time','Local','Global','Cumulative']) return df_continuity # 連続の式の残差プロット def continuity_error_plot(df): x = ['Time'] y = ['Local', 'Global', 'Cumulative'] for i in range(len(y)): plt.plot(df[x[0]], df[y[i]]) plt.grid() plt.legend(loc='best', bbox_to_anchor=(1, 1)) plt.xlabel(x[0]) plt.ylabel(y[i] + ' ' + 'continuity') plt.show() df = read_continuity_error() continuity_error_plot(df)
以下に上記Pythonスクリプトでプロットした図を示す。
流速、乱流諸量の残差は、initialとfinalの2種類、連続式の残差はLocal、Global、Cumulativeの3種類ある。これらの違いは追々調査してまとめる。
reactingFoam機能確認その1
概要
OpenFOAMでは、化学反応や燃焼を計算するソルバーとして、reactingFoamが用意されている。本記事は、reactingFoamの基本的な機能を確認した際のメモである。
実施した例題
reactingFoam内にあるチュートリアルのうち、以下の例題を実施した。
- counterFlowFlame2D
燃焼計算を行う場合
チュートリアルファイルを変更せず、そのまま実行する。
blockMesh reactingFoam
以下結果を示す。
上記の結果より、左側より燃料（メタン）、右側より空気（酸素）が流入し、中央部で燃料と空気が混合、燃焼反応が生じて温度の上昇が見られる。
コールドフロー計算を行う場合
次は、反応を伴わない、つまりコールドフローの計算を行うことを考える。例えば、化学種の移流拡散計算を行うことに相当する。
反応を伴わない計算を実行する場合、constant/chemistryPropertiesファイルとconstant/combustionPropertiesファイルを編集する必要がある。
具体的には以下を変更する。
- constant/chemistryPropertiesの場合
chemistry on → offに変更.
- constant/combustionPropertiesの場合
active true → fasleに変更.
下記に結果を示す。
左側からメタン、右側から酸素が流入しているが、燃焼反応がせず温度上昇は見られないことより、コールドフローの計算がされていることが確認できる。
