【OpenFOAM(円筒内の流れ)】simpleFoamで定常解析(5)

今回の記事では円筒内の流れをOpenFOAMでやってみたという内容です。
4パターンの異なるメッシュの作り方で解析を試してみました。

メッシュの作り方については↓こちらをご参考ください。

モデル作成

メッシュ作成は前回の記事で終えているので解析条件を設定します。
今回行う解析はΦ20×200の円柱内の流れで、メッシュパターンは以下の4パターンで行います。

境界条件は全て同じ条件（流入速度、面の名前）で与えています。

層流と乱流の速度分布の違いを見るために、助走区間の距離を考慮してモデルを作っています。助走区間は層流において十分に発達するまでの距離が長いので、大きくとる必要があります。

乱流において
乱流状態と同じモデルにするために$Re=4000$での乱流状態での助走区間は$d=20$mmとすると$\frac{L}{d}=10$より、$L=200$mmとなります。

層流において
助走区間の距離$L=200$mmとすると$\frac{L}{d}=0.05Re$よりレイノルズ数$Re=200$なります。念のため助走区間をもっと短くするためにレイノルズ数を小さくして$Re=100$にして解析を行います(この方がより層流側の流れになるため理論値と一致した)。

これによりモデルの円筒の長さは最低200mmにすれば良いことがわかります。

メッシュの作り方で試してみると勉強になりますね。

モデルの保存先

モデル保存先

Cylinder_Mesh_simpleFoam

流れはこんな感じです。

円筒内の流れの理論式

助走区間の距離も十分にとり、十分発達した流れの速度分布は軸方向に変化をしないから$\frac{\partial u_{z}}{\partial z}=0$,$u_{x}=u_{y}=u_{r}=0$（半径方向の流速0）となります。

\begin{align*}
u_{z}\frac{\partial u_{z}}{\partial z}+u_{r}\frac{\partial u_{z}}{\partial r}=-\frac{1}{\rho}\frac{\partial p}{\partial z}+\nu\frac{1}{r}\frac{\partial }{\partial r}\bigg(\frac{\partial u_{z}}{\partial r}\bigg)
\end{align*}

※今回は$z$方向が流路軸方向です。
$\frac{\partial u_{z}}{\partial z}=0$,$u_{x}=u_{y}=u_{r}=0$のとき、

\begin{align*}
0=-\frac{1}{\rho}\frac{\partial p}{\partial z}+\nu\frac{1}{r}\frac{\partial }{\partial r}\bigg(\frac{\partial u_{z}}{\partial r}\bigg)
\end{align*}

となるので、両辺$r$で積分して、

\begin{align*}
u_{z}&=\frac{1}{4\mu}\bigg(-\frac{\partial p}{\partial z}\bigg)\big(R^2-r^2\big)\\
&=2u_{b}\bigg(1-\big(\frac{r}{R}\big)^2\bigg)
\end{align*}

となります。
※$\nu=\frac{\mu}{\rho}$
ここで、$u_{b}=\frac{1}{\pi R^2}\int^{R}_{0}2\pi r u_{z}dr=\frac{R^2}{8\mu}\big(-\frac{\partial p}{\partial z}\big)$の平均流速を導入しています。

グラフにまとめる

Pythonのプログラムで層流と乱流での流速分布の比較を行いました。
また、層流に関しては理論式が近似的に導かれるのでその比較も行っています。

層流の場合の円筒内の流れの近似式はこちらです。

\begin{align*}
u=2u_{b}\bigg(1-\big(\frac{r}{R}\big)^2\bigg)
\end{align*}

$u_{b}$は円筒断面の平均流速です。
流入速度を一定値$1.0$[m/s]としているのでその値を与えています。

import re
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import subprocess
import os, shutil
from pathlib import Path


# 文字列を含む要素を数字で昇順にする関数
def atoi(text):
    print(text)
    return int(text) if text.isdigit() else text

def natural_keys(text):
    return [ atoi(c) for c in re.split(r'(\d+)^', text) ]

# 数字であるかどうかの判定
def is_int(s):
    try:
        int(s)
        return True
    except ValueError:
        return False

def graph_plot(clyndeMeshType):

    # 計算パターンのリスト
    resultDir_list = [d for d in os.listdir("./") if clyndeMeshType in d]


    # グラフのセット    
    fig, axes = plt.subplots(1, 4, figsize=(24,8), tight_layout=True)

    # colorの設定（laminar=>red, turbulence=>blue)
    color_list = ["red","blue"]

    Radius = 10/1000 # 半径 10mm = 0.01m
    ub = 1.0        # 断面平均速度(流入速度1.0m/s)
    y = np.arange(-1,1,0.01)
    uy = 2*ub*(1-y**2)

    for j, resultDir in enumerate(resultDir_list):
        resultDir = Path(resultDir)  # 計算フォルダ
        latestTime = sorted([a for a in os.listdir(resultDir) if is_int(a)], key=float)[-1] # latestTimeを抜き出す
        print(latestTime)
        sampleDictdata = resultDir/"postProcessing"/"sampleDict"/latestTime


        # ファイル名を数字部分と非数字部分に分割するための正規表現パターン
        pattern = r'(\d*\.\d+|\d+)'
        # 数字部分だけを抜き出して、数値に変換してソートする
        caseList = sorted(os.listdir(sampleDictdata), key=lambda f: [float(d) if d.isdigit() else d for d in re.split(pattern, f)])
 
        for i, zCase in enumerate(caseList):
            if j ==0 :axes[i].plot(uy, y, color="black", label="Theory")

            df_data = pd.read_table(sampleDictdata/zCase, index_col=False, names=('x[m]', 'y[m]', 'z[m]', 'Ux[m/s]', 'Uy[m/s]', 'Uz[m/s]'))
            axes[i].scatter(df_data["Uz[m/s]"], df_data["x[m]"]/Radius, color=color_list[j], label=f"{resultDir}")


            axes[i].set_ylabel("y[m]", fontsize=16)
            axes[i].set_xlabel("Uz[m/s]", fontsize=16)
            if i>=3: plt.legend(loc="upper center",bbox_to_anchor=(1.5, 1.0), fontsize=12)
            axes[i].grid("-")
            axes[i].set_title(f"z={zCase.split('_')[1]}", fontsize=20)
            
    plt.savefig(f"{clyndeMeshType}.png")
        
        
# ================ main ===================
clynderMeshTypeList = [
    "Cylinder_axialSymmetry",
    "Cylinder_blockMesh",
    "Cylinder_cfMesh",
    "Cylinder_snappyHexMesh"
    ]

for clyndeMeshType in clynderMeshTypeList:
    graph_plot(clyndeMeshType)

import re

import numpy as np

import matplotlib.pyplot as plt

import pandas as pd

import subprocess

import os, shutil

from pathlib import Path

# 文字列を含む要素を数字で昇順にする関数

def atoi(text):

print(text)

return int(text) if text.isdigit() else text

def natural_keys(text):

return [ atoi(c) for c in re.split(r'(\d+)^', text) ]

# 数字であるかどうかの判定

def is_int(s):

try:

int(s)

return True

except ValueError:

return False

def graph_plot(clyndeMeshType):

# 計算パターンのリスト

resultDir_list = [d for d in os.listdir("./") if clyndeMeshType in d]

# グラフのセット

fig, axes = plt.subplots(1, 4, figsize=(24,8), tight_layout=True)

# colorの設定（laminar=>red, turbulence=>blue)

color_list = ["red","blue"]

Radius = 10/1000 # 半径 10mm = 0.01m

ub = 1.0 # 断面平均速度(流入速度1.0m/s)

y = np.arange(-1,1,0.01)

uy = 2*ub*(1-y**2)

for j, resultDir in enumerate(resultDir_list):

resultDir = Path(resultDir) # 計算フォルダ

latestTime = sorted([a for a in os.listdir(resultDir) if is_int(a)], key=float)[-1] # latestTimeを抜き出す

print(latestTime)

sampleDictdata = resultDir/"postProcessing"/"sampleDict"/latestTime

# ファイル名を数字部分と非数字部分に分割するための正規表現パターン

pattern = r'(\d*\.\d+|\d+)'

# 数字部分だけを抜き出して、数値に変換してソートする

caseList = sorted(os.listdir(sampleDictdata), key=lambda f: [float(d) if d.isdigit() else d for d in re.split(pattern, f)])

for i, zCase in enumerate(caseList):

if j ==0 :axes[i].plot(uy, y, color="black", label="Theory")

df_data = pd.read_table(sampleDictdata/zCase, index_col=False, names=('x[m]', 'y[m]', 'z[m]', 'Ux[m/s]', 'Uy[m/s]', 'Uz[m/s]'))

axes[i].scatter(df_data["Uz[m/s]"], df_data["x[m]"]/Radius, color=color_list[j], label=f"{resultDir}")

axes[i].set_ylabel("y[m]", fontsize=16)

axes[i].set_xlabel("Uz[m/s]", fontsize=16)

if i>=3: plt.legend(loc="upper center",bbox_to_anchor=(1.5, 1.0), fontsize=12)

axes[i].grid("-")

axes[i].set_title(f"z={zCase.split('_')[1]}", fontsize=20)

plt.savefig(f"{clyndeMeshType}.png")

# ================ main ===================

clynderMeshTypeList = [

"Cylinder_axialSymmetry",

"Cylinder_blockMesh",

"Cylinder_cfMesh",

"Cylinder_snappyHexMesh"

]

for clyndeMeshType in clynderMeshTypeList:

graph_plot(clyndeMeshType)

こちらを実行すると以下のようにグラフ化されます。

軸対称モデルでの結果
【OpenFOAM(円筒内の流れ)】軸対称モデルでメッシュ作成(4)

blockMeshでの結果
【OpenFOAM(円筒内の流れ)】blockMeshでメッシュ作成(1)

cfMeshでの結果
【OpenFOAM(円筒内の流れ)】cfMeshでメッシュ作成(3)

snappyHexMeshでの結果
【OpenFOAM(円筒内の流れ)】snappyHexMeshでメッシュ作成(2)

cfMeshとsnappyHexMeshで層流での点がガタついている感じですね。
もう少しメッシュ数増やすなりすれば改善するかもしれません。

また、乱流と層流で速度分布が全然違うというのも円筒内流れの特徴ですね。

まとめ

今回はメッシュの作り方違いで円筒内流れをOpenFOAMで解析してみました。
また、層流と乱流での速度分布の違いをPythonでグラフ化して確認しました。

基本的な解析内容ですが、教科書的な内容で良い勉強になりました。

【OpenFOAM(円筒内の流れ)】simpleFoamで定常解析(5)

モデル作成

モデルの保存先

円筒内の流れの理論式

グラフにまとめる

まとめ

おすすめ参考図書

COMMENT コメントをキャンセル

モデル作成

モデルの保存先

円筒内の流れの理論式

グラフにまとめる

まとめ

おすすめ参考図書

【OpenFOAM(混相流)】特定領域から時間変化で条件を与える

【OpenFOAM球体周りの抗力係数(4)】simpleFoamで定常解析でシミュレーション

【オープンソースCAE解析】OpenFOAMで流体解析！球体周りの流れ：渦の生成（第1部）

【回転するバスケットボールまわりの流れ(3)】外部領域をblockMeshでメッシュ作成

「実験とCAE 強制流」の勉強会参加！自宅でOpenFOAMを使って流体解析

【OpenFOAM(cyclicAMI)】異なるメッシュの境界で円筒内の流れ

【OpenFOAM球体周りの抗力係数(5)】Pythonでパラメータスタディ

【オープンソース流体解析】OpenFOAMカスタマイズしたいので準備をしています。

【gnuplot実験の圧力データ比較(OpenFOAM)】3次元ダムブレイクの解析(7)

OpenFOAMの熱流体固体連成のチュートリアル（snappyMultiRegionHeater）

【OpenFOAM(混相流)】円の液相率で初期状態を作る

【OpenFOAM】移流方程式で色々な離散化スキームをPythonで自動実行して試してみた

【OpenFOAM】定常流れ＋粒子追跡（カップリング無し）

【OpenFOAM】snappyHexMeshを使ったメッシュ境界層がうまくいかない

【OpenFOAM】移流方程式で色々な離散化スキーム(解の振る舞いと計算時間)

【OpenFOAM(円筒内の流れ)】snappyHexMeshでメッシュ作成(2)

【OpenFOAMの流入条件】 質量流の時刻歴変化をcsvファイルから読み込む方法

【OpenFOAM(円筒内の流れ)】軸対称モデルでメッシュ作成(4)

【OpenFOAMカスタマイズ】空間依存する温度拡散係数に変更

【OpenFOAM】境界面での流量計算(volFieldValue)

COMMENT コメントをキャンセル

【OpenFOAMの流入条件】質量流の時刻歴変化をcsvファイルから読み込む方法