概要・背景
物体周りの流れの解析を定常解析で実施する場合、特に形状が複雑な場合、目標残差に収束せず計算が何ステップにわたり継続し、計算がなかなか終了できない状況も起こりうる。
そのような場合、残差だけではなく、評価したい物理量の履歴をプロットして、その変動が落ち着いている（変動幅が一定範囲内にあるか、周期的であるか等）ことを確認して計算の終了を判断することも一つの手段である。
評価する物理量としては、解析する対象により様々であるが、物体周りの流れ解析の場合は、例えば以下のものが代表的であろう。
- 物体に作用する流体力（圧力、粘性力、圧力と粘性力の合成）
- 抗力
- 揚力
- モーメント
ここでは、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スクリプトでプロットした図を示す。
流体力の履歴としては以下の３種類出力される。
- トータル（圧力+粘性力）
- 圧力
- 粘性力
履歴のプロットから分かるように、流体力は圧力の寄与分が支配的である。
