\(\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}}\)
非圧縮性流体の解析で使用する方程式は以下の2式である。(\(\,\bm{f_B}\,\)は重力や電磁力などの体積力である)
\(\left\{\begin{array}{l}\nabla\cdot\bm{V}=0\hspace{20mm}\cdots\,(1)\\
\pdr{\bm{V}}{t}+(\bm{V}\cdot\nabla)\bm{V}=-\nabla p+\dfrac{1}{\mathrm{Re}}\nabla^2\bm{V}+\bm{f_B}\qquad\cdots\,(2)\end{array}\right.\)
ここで(2)式を時刻について差分化して、仮の速度場\(\,\bm{V}^{\,*}\,\)を求める。
\(\quad\dfrac{\bm{V}^{\,*}-\bm{V}}{\Delta t}=-(\bm{V}\cdot\nabla)\bm{V}-\nabla p+\dfrac{1}{\mathrm{Re}}\nabla^2\bm{V}+\bm{f_B}\qquad\cdots\,(3)\)
この速度場\(\,\bm{V}^{\,*}\,\)は連続の式((1)式)を満足していないので、速度と圧力を補正して時間発展を求める。
\(\quad \bm{V}^{\,n+1}=\bm{V}^{\,*}+\bm{V}^{\,\prime}\quad,\qquad p^{\,n+1}=p^{\,n}+p^{\,\prime}\qquad\cdots\,(4)\)
この(4)式を(1),(2)式に代入して、式(3)を考慮すると以下の式が得られる。
\(\left\{\begin{array}{l}\nabla\cdot\bm{V}^{\,\prime}=-\nabla\cdot\bm{V}^{\,*}\hspace{18mm}\cdots\,(5)\\
\dfrac{\bm{V}^{\prime}}{\Delta t}=-\dfrac{1}{\mathrm{Re}}\nabla p^{\,\prime}\hspace{24mm}\cdots\,(6)\end{array}\right.\)
(6)式を(5)式に代入すると圧力補正\(\,p^{\,\prime}\,\)に対するポアソン方程式が得られる。
\(\quad\nabla^2 p^{\,\prime}=\dfrac{\mathrm{Re}}{\Delta t}\,\nabla\cdot\bm{V}^{\,*}\qquad\cdots\,(7)\)
◆ 手順
- １．式(3)を使って\(\,\bm{V}^{\,*}\,\)を求める
- ２．求めた\(\,\bm{V}^{\,*}\,\)から式(7)をSOR法で解いて、\(p^{\,\prime}\)を求める
- ３．式(6)を使って、\(\bm{V}^{\,\prime}\,\)を求める
- ４．式(4)を用いて、\(p^{\,n+1}\,,\,\bm{V}^{\,n+1}\,\)を求める
◆ 問題
図のような流路を考える。
計算領域は高さ\(\,L_y\)、横幅\(\,L_x\,\)の横長の長方形であり、左境界の上半分の部分(\(\,x=0\,,\,L_y\,/\,2\ge y\ge L_y\,\))から流体が流入し、右境界から流出する。
メッシュはスタガード格子を用いて解く。
◆ 境界条件
・流入境界
入口の外側(\(\,i=0\,\)のセル)では、\(u_{0,j}=U_{in}\,\,,\,v_{0,j}=2V_{in}-v_{1,j}\,\,,\,p^{\,\prime}_{0,j}=p^{\,\prime}_{1,j}\,\)となる。
圧力補正値\(\,p^{\,\prime}\,\)の境界条件は法線方向の勾配\(\,0\,\)のノイマン条件となる。
・流出境界
出口の外側(\(\,i=N_i+1\,\)のセル)では、各速度成分および物理量が移流によってスムーズに流出すると考え、境界法線方向の勾配を\(\,0\,\)とすると、\(\,u_{N_i+1,j}=u_{N_i,j}\,\,,\,v_{N_i+1,j}=v_{N_i,j}\,\)となる。
また、流出境界で断面内の圧力が一様で一定値と見なせるとし、その出口圧力\(\,P_{out}\,\)を与える。条件は
\(\quad p^{\,\prime}_{N_i+1,j}=p^{\,\prime}_{N_i,j}\,\,,\,p_{N_i+1,j}=2P_{out}-p_{N_i,j}\,\)となる。
◆ 計算結果
\(\quad L_y=1\,\,,\,L_x=5\,\,,\,U_{in}=2\,\,,\,V_{in}=0\,\,,\,P_{out}=0\,\)で、格子を\(\,65\,\times\,13\,\)で分割した結果を示す。
一番目は流れのベクトルプロットである。入口部分に剥離領域が生じているのが解る。次が圧力と流線である。３番目に横方向の\(\,4/5\,\)の位置の断面の流速分布である。流入流速で正規化してあるので、平均流速は\(\,1/2\,\)となっている。
◆ 参考文献
- 数値流体解析の基礎：肖 鋒、長崎 孝夫(2020)、コロナ社
- 流体力学第２版：杉山 弘、遠藤 剛、新井 隆景(2014)、森北出版
import numpy as np
from matplotlib import pyplot as plt
#from mpl_toolkits.mplot3d import axes3d
import japanize_matplotlib
# 事前に > pip install japanize_matplotlib が必要
import warnings
# 警告非表示
#warnings.simplefilter('ignore')
# 条件設定
Lx = 5.
Ly = 1.
nx = 65 # 分割数 X
ny = 13 # 分割数 Y
dx = Lx / float(nx) # 格子サイズ X
dy = Ly / float(ny) # 格子サイズ Y
re = 150 # レイノルズ数
vis = 1./re # 動粘性係数
idx = 1./dx
idy = 1./dy
dx2 = dx * dx
dy2 = dy * dy
idx2 = 1./ dx2
idy2 = 1./ dy2
uin = 2. # 流入速度
cn = 0.4 # 時間刻み係数
eps= 1.e-6 # SORの収束判定基準
alp = 1.8 # SORの加速係数
nx1 = nx+1
ny1 = ny+1
nx2 = nx+2
ny2 = ny+2
ja = ny // 2 # 入口のサイズ
dt = cn/(uin/dx+2.*vis*(idx2+idy2)) # 時間刻み
# 配列
u = np.zeros((ny2,nx2)) # 流れ x方向成分 u
v = np.zeros((ny2,nx2)) # 流れ y方向成分 v
p = np.zeros((ny2,nx2)) # 圧力
pc = np.zeros((ny2,nx2)) # 圧力補正
un = np.zeros((ny2,nx2)) # 流れ更新分 u
vn = np.zeros((ny2,nx2)) # 流れ更新分 v
fe = np.zeros((ny2,nx2)) # 検査体積の右界面の流速
fn = np.zeros((ny2,nx2)) # 検査体積の上界面の流速
# グラフ条件設定
nxh = nx // 2
nxx = nx // 5 * 2
uc = np.zeros((ny,nxh)) # ベクトル表示用 u
vc = np.zeros((ny,nxh)) # ベクトル表示用 v
x, y = np.meshgrid(np.linspace(0, Lx, nxh), np.linspace(0, Ly, ny))
pd = np.zeros((ny,nx)) # 圧力表示用
xp, yp = np.meshgrid(np.linspace(0, Lx, nx), np.linspace(0, Ly, ny))
uu = np.zeros(ny)
# 圧力ポアソン方程式の係数
a_e = np.zeros((ny2,nx2)) # 方程式の係数 east
a_w = np.zeros((ny2,nx2)) # 方程式の係数 west
a_n = np.zeros((ny2,nx2)) # 方程式の係数 north
a_s = np.zeros((ny2,nx2)) # 方程式の係数 south
a_p = np.zeros((ny2,nx2)) # 方程式の係数 圧力
b = np.zeros((ny2,nx2)) # 方程式の係数 div V
for j in range(1,ny1):
for i in range(1,nx1):
a_w[j][i]=idx2
a_e[j][i]=idx2
a_s[j][i]=idy2
a_n[j][i]=idy2
a_p[j][i]=2.*(idx2+idy2)
if j==1:
a_p[j][i]=a_p[j][i]-a_s[j][i]
a_s[j][i]=0
if j==ny:
a_p[j][i]-a_p[j][i]-a_n[j][i]
a_n[j][i]=0
if i==1:
a_p[j][i]=a_p[j][i]-a_w[j][i]
a_w[j][i]=0
if i==nx:
a_p[j][i]-a_p[j][i]-a_e[j][i]
a_e[j][i]=0
# 流入口の流れ設定
for j in range(ja+1,ny1):
u[j][0]=uin # 流入
un[j][0]=uin # 流れ更新分も設定
#
# 計算ループ
#
nt = 120 # 計算ループ回数
time = 0.
for n in range(1,nt+1):
# 境界条件設定
for i in range(nx1): # 境界条件 上下
u[0][i] = -u[1][i] # 下 u=0
u[ny1][i] = -u[ny][i] # 上 u=0
for j in range(ny1): # 境界条件 左右
v[j][0] = -v[j][1] # 左 v=0
u[j][nx1] = u[j][nx] # 右 du/dx=0
v[j][nx1] = v[j][nx] # dv/dx=0
p[j][nx1] = -p[j][nx] # p=0
# 仮の 流速 V* の計算
for j in range(ny1): # 流束 u
for i in range(nx1):
u1 = (u[j][i]+u[j][i+1])/2.
v1 = (v[j][i]+v[j][i+1])/2.
fe[j][i] = u1*(u[j][i+1]+u[j][i])/2.-vis*idx*(u[j][i+1]-u[j][i])
fn[j][i] = v1*(u[j+1][i]+u[j][i])/2.-vis*idy*(u[j+1][i]-u[j][i])
for j in range(1,ny1): # u の時間発展
for i in range(1,nx1):
dp = (p[j][i]-p[j][i+1])*idx
dudt=(fe[j][i-1]-fe[j][i])*idx+(fn[j-1][i]-fn[j][i])*idy+dp
un[j][i]=u[j][i]+dudt*dt
for j in range(ny): # 流束 v
for i in range(nx1):
u1 = (u[j][i]+u[j+1][i])/2.
v1 = (v[j][i]+v[j+1][i])/2.
fe[j][i] = u1*(v[j][i+1]+v[j][i])/2.-vis*idx*(v[j][i+1]-v[j][i])
fn[j][i] = v1*(v[j+1][i]+v[j][i])/2.-vis*idy*(v[j+1][i]-v[j][i])
for j in range(1,ny): # v の時間発展
for i in range(1,nx1):
dp = (p[j][i]-p[j+1][i])*idy
dvdt=(fe[j][i-1]-fe[j][i])*idx+(fn[j-1][i]-fn[j][i])*idy+dp
vn[j][i]=v[j][i]+dvdt*dt
for j in range(1,ny1): # div V の計算
for i in range(1,nx1):
b[j][i]=-((un[j][i]-un[j][i-1])*idx+(vn[j][i]-vn[j-1][i])*idy)/dt
pc[j][i]=0
# 圧力のポアソン方程式の計算 ( SOR法 )
for itx in range(1000):
err=0
for i in range(1,nx1):
for j in range(1,ny1):
a1=a_w[j][i]*pc[j][i-1]+a_e[j][i]*pc[j][i+1]
a2=a_s[j][i]*pc[j-1][i]+a_n[j][i]*pc[j+1][i]
pcnew=(a1+a2+b[j][i])/a_p[j][i]
err1=np.abs(pcnew-pc[j][i])
pc[j][i]=pc[j][i]+alp*(pcnew-pc[j][i])
if err1 > err:
err=err1
continue
if err < eps:
break
continue
# V と 圧力の更新
for j in range(1,ny1):
for i in range(1,nx1):
u[j][i]=un[j][i]+(pc[j][i]-pc[j][i+1])*idx*dt
if j != ny:
v[j][i]=vn[j][i]+(pc[j][i]-pc[j+1][i])*idy*dt
p[j][i]=p[j][i]+pc[j][i]
time=time + dt
# 結果の表示
print("n=",n," time=",'{:.2f}'.format(time)," itx=",itx," err=",'{:.2e}'.format(err))
for j in range(ny): # V 速度ベクトル
for i in range(nxh): # x 方向の表示は半分に
m = i*2
uc[j][i]=(u[j+1][m+1]+u[j+1][m])/2. # セルの中心
vc[j][i]=(v[j+1][m+1]+v[j][m+1])/2.
for j in range(ny): # 圧力
for i in range(nx):
pd[j][i]=p[j+1][i+1]
for j in range(ny):
uu[j]=uc[j][nxx]/uin
# グラフ表示
fig=plt.figure(figsize=(12,4))
ax1 = fig.add_subplot(1, 1, 1)
ax1.set(xlabel='x',ylabel='y')
ax1.set_title('流れ',fontsize=20)
ax1.quiver(x,y,uc,vc,color='blue',scale=80,headwidth=4,headaxislength=4,width=0.0015)
plt.show()
fig=plt.figure(figsize=(12,4))
ax2 = fig.add_subplot(1, 1, 1)
ax2.set(xlabel='x',ylabel='y')
ax2.set_title('圧力',fontsize=20)
ax2.contourf(xp, yp, pd, cmap='cool')
ax2.streamplot(x,y,uc,vc,color='k',density=0.5,linewidth=1)
plt.show()
plt.figure(figsize=(5,5))
plt.title('流れの速度 at Lx = 4',fontsize=20)
plt.ylabel('速度(入口比)',fontsize=16)
plt.xlabel('縦方向',fontsize=16)
plt.plot(uu)
plt.show()
n= 120 time= 1.57 itx= 149 err= 9.44e-07
