概要・背景

物体周りの流れの解析を定常解析で実施する場合、特に形状が複雑な場合、目標残差に収束せず計算が何ステップにわたり継続し、計算がなかなか終了できない状況も起こりうる。

そのような場合、残差だけではなく、評価したい物理量の履歴をプロットして、その変動が落ち着いている（変動幅が一定範囲内にあるか、周期的であるか等）ことを確認して計算の終了を判断することも一つの手段である。

評価する物理量としては、解析する対象により様々であるが、物体周りの流れ解析の場合は、例えば以下のものが代表的であろう。

物体に作用する流体力（圧力、粘性力、圧力と粘性力の合成）
抗力
揚力
モーメント

ここでは、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 スクリプトを示す。

流体力の履歴をプロットする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)