【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 x}\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】snappyHexMeshを使ったメッシュ境界層がうまくいかない

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

【OpenFOAM】キャビティ流れをマッピングして粒子追跡

【OpenFOAM】流体解析用のベンチマークのサイト

【OpenFOAM球体周りの抗力係数(3)】snappyHexMeshで球体周りのメッシュ作成

【OpenFOAM】レイリー・テーラー不安定性の流体解析

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

【OpenFOAM(円筒内の流れ)】境界条件をcodeStremでコードを埋め込む

【OpenFOAM】円筒形状の異なるメッシュの境界を結合する

【水面高さ(OpenFOAM)】ParaViewで表示してPythonスクリプトにする(11)

【圧力データ(OpenFOAM)】ParaViewで表示してExcelデータと比較(9)

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

【3次元ダムブレイク流体解析(OpenFOAM)】圧力データと水面高さの実機比較

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

【refineMesh(OpenFOAM)】メッシュの再分割(6)

【gnuplot実験の水面高さデータ比較(OpenFOAM)】3次元ダムブレイクの解析(8)

【圧力時刻歴データ(OpenFOAM)】ParaViewとPythonを使って自動化する方法(10)

【3次元ダムブレイク(OpenFOAM)】クーラン数違いによる結果比較

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

【OpenFOAM入門】熱交換器の軸対称モデルでは失敗した

COMMENT コメントをキャンセル