\(\def\bm{\boldsymbol}\)
\(\def\di{\displaystyle}\)
\(\def\ve{\varepsilon_0}\)
\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)
\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)

◆ キャビティ流れ 3D

先日、「流れ解析入門」の講義を受けました。その内容を忘れないように、河村先生の著書を見ながら３次元キャビティ問題を Python の勉強も兼ねて、コードを書いてみました。配列の処理などは、スライスを使えば、かなり速くなるらしいがそのままにしています。今回はポアソン方程式の解法に、フラクショナル・ステップ法を使用しました。まず、その解説から。

非圧縮性ナビエストークス方程式は以下の形で表すことが出来る。

\(\qquad \pdr{\bm{V}}{t}+(\bm{V}\cdot\nabla)\bm{V}=-\nabla p+\dfrac{1}{\mathrm{Re}}\Delta\bm{V}\quad\cdots\,\)(1)

式(1)の左辺の時間微分項を前進差分で近似すれば時間ステップ\(\,n\,\)において以下の形になる。

\(\qquad\dfrac{\bm{V}^{n+1}-\bm{V}}{\Delta t}+(\bm{V}\cdot\nabla)\bm{V}=-\nabla p+\dfrac{1}{\mathrm{Re}}\Delta\bm{V}\quad\cdots\,\)(2)

次に式(1)で圧力項を落とした方程式を同じく前進差分で近似すれば

\(\qquad\dfrac{\bm{V}^*-\bm{V}}{\Delta t}+(\bm{V}\cdot\nabla)\bm{V}=\dfrac{1}{\mathrm{Re}}\Delta\bm{V}\quad\cdots\,\)(3)

ここで、\(\,\bm{V}^*\,\)は「仮の速度」と呼ばれ圧力項のない方程式を近似している。式(2)から式(3)を引くと、

\(\qquad\dfrac{\bm{V}^{n+1}-\bm{V}^*}{\Delta t}=-\nabla p\quad\cdots\,\)(4)

になる。この式の両辺の発散をとると、

\(\qquad\dfrac{\nabla\cdot\bm{V}^{n+1}-\nabla\cdot\bm{V}^*}{\Delta t}=-\nabla\cdot\nabla p\quad\cdots\,\)(5)

となる。連続の式\(\,\nabla\cdot\bm{V}=0\,\)を考慮すると、圧力のポアソン方程式が得られる。

\(\qquad\Delta p=\dfrac{\nabla\cdot\bm{V}^*}{\Delta t}\quad\cdots\,\)(6)

が得られる。\(\,\bm{V}^*\,\)は式(3)から計算して、圧力\(\,p\,\)を求める。さらに、式(4)から

\(\qquad\bm{V}^{n+1}=\bm{V}^*-\Delta t\nabla p\quad\cdots\,\)(7)

が得られる。

◆ 3次元直交座標での計算
(1) \(\bm{V}^*\,\)の計算

\(\left\{\begin{array}{@{\\,}l}\,\,u^*=u+\Delta t\left\{-u\pdr{u}{x}-v\pdr{u}{y}-w\pdr{u}{z}+\dfrac{1}{\mathrm{Re}}\left(\ppdr{u}{x^2}+\ppdr{u}{y^2}+\ppdr{u}{z^2}\right)\right\}\quad\cdots\,(10)\\
\,\,v^*=v+\Delta t\left\{-u\pdr{v}{x}-v\pdr{v}{y}-w\pdr{v}{z}+\dfrac{1}{\mathrm{Re}}\left(\ppdr{v}{x^2}+\ppdr{v}{y^2}+\ppdr{v}{z^2}\right)\right\}\quad\cdots\,(11)\\
\,\,w^*=w+\Delta t\left\{-u\pdr{w}{x}-v\pdr{w}{y}-w\pdr{w}{z}+\dfrac{1}{\mathrm{Re}}\left(\ppdr{w}{x^2}+\ppdr{w}{y^2}+\ppdr{w}{z^2}\right)\right\}\quad\cdots\,(12)\end{array}\right.\)

(2) 圧力のポアソン方程式の計算

\(\qquad\ppdr{p}{x^2}+\ppdr{p}{y^2}+\ppdr{p}{z^2}=\dfrac{1}{\Delta t}\left(\pdr{u^*}{x}+\pdr{v^*}{y}+\pdr{w^*}{z}\right)\quad\cdots\,\)(13)

(3) \(\bm{V}^*\,\)の計算

\(\left\{\begin{array}{@{\\,}l}\,\,u^{n+1}=u^*-\Delta t\pdr{p}{x}\quad\cdots\,(16)\\
\,\,v^{n+1}=v^*-\Delta t\pdr{p}{y}\quad\cdots\,(17)\\
\,\,w^{n+1}=w^*-\Delta t\pdr{p}{z}\quad\cdots\,(18)\end{array}\right.\)

◆ 参考文献

流体解析の基礎：河村哲也(2014)、朝倉書店
数値流体解析の基礎：肖鋒、長崎孝夫(2020)、コロナ社
流体力学第２版：杉山弘、遠藤剛、新井隆景(2014)、森北出版

In [8]:

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import japanize_matplotlib
import warnings
# 警告非表示
warnings.simplefilter('ignore')

# 条件設定
nx = 21
ny = 21
nz = 21
re = 40      # レイノルズ数
dt = 0.01    # 時間刻み
eps= 0.00001 # ポアソン方程式計算誤差限度
lm = 20      # ポアソン方程式の計算回数
nx1 = nx-1
ny1 = ny-1
nz1 = nz-1
nx2 = nx-2
ny2 = ny-2
nz2 = nz-2
dx = 1. / float(nx1)
dy = 1. / float(ny1)
dz = 1. / float(nz1)
td= 1. / dt
xdh=.5 * float(nx1)
ydh=.5 * float(ny1)
zdh=.5 * float(nz1)
dx2 = dx * dx
dy2 = dy * dy
dz2 = dz * dz
idx2 = 1./ dx2
idy2 = 1./ dy2
idz2 = 1./ dz2
idxyz = .5*dx2*dy2*dz2/(dy2*dz2+dx2*dz2+dx2*dz2)
ire = 1. / float(re)
nyh = ny // 2+1

# 配列
u = np.zeros((nx,ny,nz))   # 流れ x方向成分 u
ut= np.zeros((nx,ny,nz))   # 仮の流れ x方向成分 u
v = np.zeros((nx,ny,nz))   # 流れ y方向成分 v
vt= np.zeros((nx,ny,nz))   # 仮の流れ y方向成分 v
w = np.zeros((nx,ny,nz))   # 流れ z方向成分 w
wt= np.zeros((nx,ny,nz))   # 仮の流れ z方向成分 w
uh = np.zeros((nz,nx))     # 断面表示用 u
vh = np.zeros((nz,nx))     # 断面表示用 v

p = np.zeros((nx,ny,nz))   # 圧力
q = np.zeros((nx,ny,nz))   # ポアソン方程式の右辺の値

# グラフ条件設定
x, y, z = np.meshgrid(np.linspace(0, 1, nx),
                      np.linspace(0, 1, ny),
                      np.linspace(0, 1, nz))  # ３次元
X, Y = np.meshgrid(np.linspace(0, 1, nx),
                      np.linspace(0, 1, nz))  # ２次元

#
# 計算ループ
#
nt = 300     # 計算ループ回数

for n in range(1,nt+1):

#　境界条件設定

    for i in range(nx):  # 境界条件 上下
        for j in range(ny):
            u[i][j][0] = ut[i][j][0] = 0.       # 下
            v[i][j][0] = vt[i][j][0] = 0.
            w[i][j][0] = wt[i][j][0] = 0.
            u[i][j][nz1] = ut[i][j][nz1] = 0.95  # 上 ux=0.95  
            v[i][j][nz1] = vt[i][j][nz1] = 0.15  #    uy=0.15
            w[i][j][nz1] = wt[i][j][nz1] = 0.
                  
    for j in range(ny):  # 境界条件 左右
        for k in range(nz):
            u[0][j][k] = ut[0][j][k] =0.       # 左
            v[0][j][k] = vt[0][j][k] =0.
            w[0][j][k] = wt[0][j][k] =0.
            u[nx1][j][k] = ut[nx1][j][k] =0.   # 右
            v[nx1][j][k] = vt[nx1][j][k] =0.
            w[nx1][j][k] = wt[nx1][j][k] =0.

    for k in range(nz):  # 境界条件 前後
        for i in range(nx):
            u[i][0][k] = ut[i][0][k] =0.       # 前
            v[i][0][k] = vt[i][0][k] =0.
            w[i][0][k] = wt[i][0][k] =0.
            u[i][ny1][k] = ut[i][ny1][k] =0.   # 奥
            v[i][ny1][k] = vt[i][ny1][k] =0.
            w[i][ny1][k] = wt[i][ny1][k] =0.

#　仮の　V*　の計算

    for k in range(1, nz1):
        for j in range(1, ny1):
            for i in range(1, nx1):
                rx = u[i][j][k]*(u[i+1][j][k]-u[i-1][j][k])*xdh \
                    +v[i][j][k]*(u[i][j+1][k]-u[i][j-1][k])*ydh \
                    +w[i][j][k]*(u[i][j][k+1]-u[i][j][k-1])*zdh
                vx = (u[i+1][j][k]-2*u[i][j][k]+u[i-1][j][k])*idx2 \
                    +(u[i][j+1][k]-2*u[i][j][k]+u[i][j-1][k])*idy2 \
                    +(u[i][j][k+1]-2*u[i][j][k]+u[i][j][k-1])*idz2
                ut[i][j][k]=u[i][j][k]+dt*(-rx+vx*ire)  # 式(10)

                ry = u[i][j][k]*(v[i+1][j][k]-v[i-1][j][k])*xdh \
                    +v[i][j][k]*(v[i][j+1][k]-v[i][j-1][k])*ydh \
                    +w[i][j][k]*(v[i][j][k+1]-v[i][j][k-1])*zdh
                vy = (v[i+1][j][k]-2*v[i][j][k]+v[i-1][j][k])*idx2 \
                    +(v[i][j+1][k]-2*v[i][j][k]+v[i][j-1][k])*idy2 \
                    +(v[i][j][k+1]-2*v[i][j][k]+v[i][j][k-1])*idz2
                vt[i][j][k]=v[i][j][k]+dt*(-ry+vy*ire)  # 式(11)
                
                rz = u[i][j][k]*(w[i+1][j][k]-w[i-1][j][k])*xdh \
                    +v[i][j][k]*(w[i][j+1][k]-w[i][j-1][k])*ydh \
                    +w[i][j][k]*(w[i][j][k+1]-w[i][j][k-1])*zdh
                vz = (w[i+1][j][k]-2*w[i][j][k]+w[i-1][j][k])*idx2 \
                    +(w[i][j+1][k]-2*w[i][j][k]+w[i][j-1][k])*idy2 \
                    +(w[i][j][k+1]-2*w[i][j][k]+w[i][j][k-1])*idz2
                wt[i][j][k]=w[i][j][k]+dt*(-rz+vz*ire)  # 式(12)
                
#　圧力のポアソン方程式の計算
 
    for k in range(1, nz1):
        for j in range(1, ny1):
            for i in range(1, nx1):
                q[i][j][k]=( (ut[i+1][j][k]-ut[i-1][j][k])*xdh \
                            +(vt[i][j+1][k]-vt[i][j-1][k])*ydh \
                            +(wt[i][j][k+1]-wt[i][j][k-1])*zdh )*td # 式(13) の右辺

    for l in range(1,lm):  # ポアソン方程式の解法

        for k in range(1, nz): # 圧力設定 左右 内側からコピー
            for j in range(1,ny):
                p[0][j][k]=p[1][j][k]
                p[nx1][j][k]=p[nx2][j][k]

        for k in range(1, nz): # 圧力設定 前後 内側からコピー
            for i in range(1, nx):
                p[i][0][k]=p[i][1][k]
                p[i][ny1][k]=p[i][ny2][k]

        for j in range(1, ny): # 圧力設定 上下 内側からコピー
            for i in range(1, nx):
                p[i][j][0]=p[i][j][1]
                p[nx1][j][nz1]=p[i][j][nz2]

        g2=0.
        
        for k in range(1, nz1):  # 解法ループ
            for j in range(1, ny1):
                for i in range(1, nx1):
                    rhs=(p[i+1][j][k]+p[i-1][j][k])*idx2+(p[i][j+1][k]+p[i][j-1][k])*idy2 \
                        +(p[i][j][k+1]+p[i][j][k-1])*idz2
                    uli=(rhs-q[i][j][k])*idxyz-p[i][j][k]  # 式(13)
                    g2+=uli*uli
                    p[i][j][k]+=uli

        if(g2<=eps*nz):  # 誤差以下ならループを抜ける
            break;

# 流れ$\,\bm{V}^{n+1}\,$の計算
    for k in range(1, nz1):
        for j in range(1, ny1):
            for i in range(1, nx1):
                u[i][j][k]=ut[i][j][k]-dt*(p[i+1][j][k]-p[i-1][j][k])*xdh
                v[i][j][k]=vt[i][j][k]-dt*(p[i][j+1][k]-p[i][j-1][k])*ydh
                w[i][j][k]=wt[i][j][k]-dt*(p[i][j][k+1]-p[i][j][k-1])*zdh

                
    # グラフ表示
    if n%100 == 0 or n==1 or n==10:
        print("loop= ", n)
        for k in range(nz):
            for i in range(nx):
                uh[k][i]=u[i][nyh][k]
                vh[k][i]=v[i][nyh][k]
        uh[ny1][::]=vh[ny1][::]=uh[ny2][::]=vh[ny2][::]=0 # 飛び出すベクトル消す
        fig=plt.figure(figsize=(10,5))
        ax1 = fig.add_subplot(1, 2, 1,projection='3d')
        ax1.set(xlabel='x',ylabel='y',zlabel='z')
        ax1.set_title('流れ3D')
        ax1.quiver(x,y,z,u,v,w,length=0.04,linewidth=1,color='blue')
        ax2 = fig.add_subplot(1, 2, 2)
        ax2.set(xlabel='x',ylabel='z')
        ax2.set_title('y=0.5 断面')
        ax2.quiver(X,Y,uh,vh,linewidth=1,color='blue')
        plt.subplots_adjust(left=0.125, bottom=0.2, 
                    right=0.9, top=0.9, wspace=0.3, hspace=0.3)
        plt.show() 
    if n==nt:
        for k in range(nz):
            for i in range(nx):
                uh[k][i]=p[i][nyh][k]

        fig=plt.figure(figsize=(4,4))
        plt.xlabel("X")
        plt.ylabel(("Z"))
        plt.title('圧力(y=0.5 断面)')
        plt.contourf(X, Y, uh, cmap='cool')
        plt.show

loop=  1

loop=  10

loop=  100

月別アーカイブ: 2022年8月

流れ解析の練習-1 キャビティ流れ 3D