りぽログ

twitterアカウント→@nakayoshi0505 とあるChemical Engineerによる徒然日記。好物はリポビタンD。たまに趣味(ゲーム・ドライブ)のことも?

Pythonによる流体シミュレーション(円柱周りの流れ)

はじめに

使用環境
MacOS Mojave バージョン 10.14.3
・anaconda3
・Python3.6
・Jupyter notebook

今回は友人からのリクエスト「流体力学」について書いていきます。「流体力学についてなんか書いて!」と言われましても自分自身専門外で、特に面白いネタも持ってないので、円柱周りの流れのシミュレーションをしてみることにしました。

NS式を数値的に解く

流体力学において最も有名な式の1つ「ナビエ・ストークス方程式(NS式)」を数値的に解く方法として、MAC法というのがあります。大雑把にいうと、MAC法は連続の式を満たすような圧力Pを予め決定し(ポアソン方程式)、その値を用いて次のステップの速度を計算する(NS式)といった方法です。詳しい説明は調べたら出てくると思うので省略。実際に使ったコードは最後にあります。

シミュレーション結果

f:id:cheme_nd:20190216031815g:plain
色の薄い部分が速度の早いところ、濃い部分が速度の遅いところです。中心の濃い円が円柱(速度0)を表しています。

何となく後方に渦が2つ伸びていくのがわかる?かな?(微妙)あまり綺麗になりませんでした...もっと細かく区切った方がいいんでしょうがオーバーフローするし、ただでさえ激遅な計算がさらに遅くなってしまいます。

でもまぁ色々勉強になったので満足(アニメーションの描画とか)

コード全文

ライブラリをインポート

%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
import sys,random
import matplotlib.animation as anm

MAC法でナビエストークス方程式を数値的に解く

Re = 100
####
#grid making
igr, ityp = 0, 1
Jmax = 70
Kmax = 35
AA  = 0.75
BB = 1.0
HH = 0.025
RA  = 1.15
RR = np.zeros(Kmax+1)
RR[1] = 1.0
for k in range(2, Kmax+1):
 RR[k]  = RR[k-1] + HH*RA**(k-1)

x, y = np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1])
for k in range(1, Kmax+1):
 for j in range(1, Jmax+1):
  TT  =2.0*np.pi * float(j-2)/float(Jmax-2)
  x[j, k] = AA*RR[k]*np.cos(TT)
  y[j, k] = AA*RR[k]*np.sin(TT)
        
#input data
Jm = Jmax - 1
Km = Kmax - 1
nsteps = 4000
istep0 = 20
#dt = 0.0004
dt = 0.0003*np.sqrt(Re)
print("Re =",Re)
print("dt =",dt)
eps = 0.001
const = 1.0
Rei = 1.0/Re
dti = 1.0/dt

#caluculating metrics
xx, xy, yx, yy = np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1])
c1, c2, c3 = np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1])
c4, c5, aj= np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1])

for k in range(1, Kmax+1):
 for j in range(1, Jmax+1):
        
  if(k == 1):
   xe = 0.5*(- x[j, 3] + 4.0*x[j, 2] - 3.0*x[j, 1])
   ye = 0.5*(- y[j, 3] + 4.0*y[j, 2] - 3.0*y[j, 1])
  elif(k == Kmax):
   xe = 0.5*(x[j, Kmax-2] - 4.0*x[j, Kmax-1] + 3.0*x[j, Kmax])
   ye = 0.5*(y[j, Kmax-2] - 4.0*y[j, Kmax-1] + 3.0*y[j, Kmax])
  else:
   xe = 0.5*(x[j, k+1] - x[j, k-1])
   ye = 0.5*(y[j, k+1] - y[j, k-1])
            
  if(j == 1):
   xxi = 0.5*(- x[3, k] + 4.0*x[2, k] - 3.0*x[1, k])
   yxi = 0.5*(- y[3, k] + 4.0*y[2, k] - 3.0*y[1, k])
   if(ityp == 1):
    xxi = 0.5*(x[2, k] - x[Jmax -2, k])
    yxi = 0.5*(y[2, k] - y[Jmax -2, k])
  elif(j == Jmax):
   xxi = 0.5*(x[Jmax-2, k] - 4.0*x[Jmax-1, k] + 3.0*x[Jmax, k])
   yxi = 0.5*(y[Jmax-2, k] - 4.0*y[Jmax-1, k] + 3.0*y[Jmax, k])
   if(ityp == 1):
    xxi = 0.5*(x[3, k] - x[Jmax -1, k])
    yxi = 0.5*(y[3, k] - y[Jmax -1, k])
  else:
   xxi = 0.5*(x[j+1, k] - x[j-1, k])
   yxi = 0.5*(y[j+1, k] - y[j-1, k])   
            
  if(ityp == 1 and j ==1):
   xxi = 0.5*(x[j+1, k] - x[Jm -1, k])
   yxi = 0.5*(y[j+1, k] - y[Jm -1, k])
            
  if(ityp == 1 and j == Jmax):
   xxi = 0.5*(x[3, k] - x[j-1, k])
   yxi = 0.5*(y[3, k] - y[j-1, k])
            
  ajj = xxi*ye - xe*yxi
  xx[j, k] = ye/ajj
  yx[j, k] = - yxi/ajj
  xy[j, k] = - xe/ajj
  yy[j, k] = xxi/ajj
  aj[j, k] = ajj
        
for k in range(1, Kmax+1):
 for j in range(1, Jmax+1):
  c1[j, k] = xx[j, k]**2 + xy[j, k]**2
  c3[j, k] = yx[j, k]**2 + yy[j, k]**2
  c2[j, k] = 2.0*(xx[j, k]*yx[j, k] + xy[j, k]*yy[j, k])
        
for k in range(2, Km+1):
 for j in range(2, Jm+1):
  c77 = xx[j, k]*(xx[j+1, k] - xx[j-1, k]) + yx[j, k]*(xx[j, k+1] - xx[j, k-1]) + xy[j, k]*(xy[j+1, k] - xy[j-1, k]) + yy[j, k]*(xy[j, k+1] - xy[j, k-1])
  c88 = xx[j, k]*(yx[j+1, k] - yx[j-1, k]) + yx[j, k]*(yx[j, k+1] - yx[j, k]) + xy[j, k]*(yy[j+1, k] - yy[j-1, k]) + yy[j, k]*(yy[j, k+1] - yy[j, k-1])
  c4[j, k] = c77*0.5
  c5[j, k] = c88*0.5
            
#init condtion
alp = 0.
TTT = alp*np.pi/180.0
u, v, p = np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1])
us, vs, ps = np.zeros([Jmax+1, Kmax+1, nsteps]), np.zeros([Jmax+1, Kmax+1, nsteps]), np.zeros([Jmax+1, Kmax+1, nsteps])
q, d  = np.zeros([Jmax+1, Kmax+1]), np.zeros([Jmax+1, Kmax+1])
random.seed(1)
for k in range(1, Kmax+1):
 for j in range(1, Jmax+1):
  u[j, k] = np.cos(TTT)+0.05*(random.random()-0.5)
  v[j, k] = np.sin(TTT)+0.05*(random.random()-0.5)
  p[j, k] = 0.0 - 0.5*(u[j, k]**2 + v[j, k]**2)

#main
for n in range(1, nsteps+1):
#calculating RHS of poisson eq
 for k in range(2, Km+1):
  for j in range(2, Jm+1):
   uxd = xx[j, k]*(u[j+1, k] - u[j-1, k]) + yx[j, k]*(u[j, k+1] - u[j, k-1])
   uyd = xy[j, k]*(u[j+1, k] - u[j-1, k]) + yy[j, k]*(u[j, k+1] - u[j, k-1])
   vxd = xx[j, k]*(v[j+1, k] - v[j-1, k]) + yx[j, k]*(v[j, k+1] - v[j, k-1])
   vyd = xy[j, k]*(v[j+1, k] - v[j-1, k]) + yy[j, k]*(v[j, k+1] - v[j, k-1])
   q[j, k] = - 0.25*(uxd*uxd + 2.0*uyd*vxd + vyd*vyd) + 0.5*(uxd + vyd)*dti
            
 err = 0.0
 for i in range(1, istep0+1):
    #calculating Pressure 
  for k in range(2, Km+1):
   for j in range(2, Jm+1):
    cc = 0.5/(c1[j, k] + c3[j, k])
    pa1 = c1[j, k]*(p[j+1, k] + p[j-1, k]) + c3[j, k]*(p[j, k+1] + p[j, k-1]) + 0.25*c2[j, k]*(p[j+1, k+1] - p[j-1, k+1] - p[j+1, k-1] + p[j-1, k-1])
    pa2 = 0.5*c4[j, k]*(p[j+1, k] - p[j-1, k]) + 0.5*c5[j, k]*(p[j, k+1] - p[j, k-1])
    pa = pa1 + pa2
    pp = (pa - q[j, k])*cc
    err = err + (pp - p[j, k])**2
    p[j, k] = p[j, k]*(1.0 - const) + pp*const
                
  #Pressure boundary condition
  for j in range(2, Jm+1):
   ul = 0.5*c2[j, 1]*(u[j+1, 2] - u[j-1, 2]) + c5[j, 1]*u[j, 1]
   vl = 2.0*c3[j, 1]*v[j, 2]
   p[j, 1] = p[j, 2] - Rei*aj[j, 1]*(xx[j, 1]*vl - xy[j, 2]*ul)
   p[j, Kmax] = - 0.5*(u[j, Kmax]**2 + v[j, Kmax]**2)
            
  if(ityp == 1):
   for k in range(1, Kmax+1):
    p[1, k] = p[Jm, k]
    p[Jmax, k] = p[2, k]
  else:
   for k in range(1, Kmax+1):
    p[1, k] = p[2, k]
    p[Jmax, k] = p[Jm, k]
        
  #print(err)
  if(err < eps ):
   count = i
   break
 sys.stdout.write("\r%d" %n +" / "+"%d" % nsteps +"  "+"%d" % i +"  "+"%e" % err +"  "+"%e" % eps+"    ")
 sys.stdout.flush()
    
 #convergence check
#    print(n,err, i)
    
 #calculating velocity
 for k in range(2, Km+1):
  for j in range(2, Jm+1):
   unl1 = xx[j, k]*(u[j+1, k]**2 - u[j-1, k]**2)*0.5 + yx[j, k]*(u[j, k+1]**2 - u[j, k-1]**2)*0.5
   unl2 = xy[j, k]*(u[j+1, k]*v[j+1, k] - u[j-1, k]*v[j-1, k])*0.5 + yy[j, k]*(u[j, k+1]*v[j, k+1] - u[j, k-1]*v[j, k-1])*0.5
   unl = unl1 + unl2
   vnl1 = xx[j, k]*(u[j+1, k]*v[j+1, k] - u[j-1, k]*v[j-1, k])*0.5 + yx[j, k]*(u[j, k+1]*v[j, k+1] - u[j, k-1]*v[j, k-1])*0.5
   vnl2 = xy[j, k]*(v[j+1, k]**2 - v[j-1, k]**2)*0.5 + yy[j, k]*(v[j, k+1]**2 - v[j, k-1]**2)*0.5
   vnl = vnl1 + vnl2
   uvs1 = c1[j, k]*(u[j+1, k] - 2.0*u[j, k] + u[j-1, k]) + c2[j, k]*(u[j+1, k+1] - u[j+1, k-1] - u[j-1, k+1] + u[j-1, k-1])*0.25
   uvs2 = c3[j,k]*(u[j, k+1] - 2.0*u[j, k] + u[j, k-1]) + 0.5*c4[j, k]*(u[j+1, k] - u[j-1, k]) + 0.5*c5[j, k]*(u[j, k+1] - u[j, k-1])
   uvs = uvs1 + uvs2
   vvs1 = c1[j, k]*(v[j+1, k] - 2.0*v[j, k] + v[j-1, k]) + c2[j, k]*(v[j+1, k+1] - v[j+1, k-1] - v[j-1, k+1] + v[j-1, k-1])*0.25
   vvs2 = c3[j,k]*(v[j, k+1] - 2.0*v[j, k] + v[j, k-1]) + 0.5*c4[j, k]*(v[j+1, k] - v[j-1, k]) + 0.5*c5[j, k]*(v[j, k+1] - v[j, k-1])
   vvs = vvs1 + vvs2
   d[j, k] = u[j, k] + dt*(- unl - 0.5*xx[j, k]*(p[j+1, k] - p[j-1, k]) - 0.5*yx[j, k]*(p[j, k+1] - p[j, k-1]) + Rei*uvs)
   q[j, k] = v[j, k] + dt*(- vnl - 0.5*xy[j, k]*(p[j+1, k] - p[j-1, k]) - 0.5*yy[j, k]*(p[j, k+1] - p[j, k-1]) + Rei*vvs)
            
 for k in range(2, Km+1):
  for j in range(2, Jm+1):
   u[j, k] = d[j, k]
   v[j, k] = q[j, k]

 #boundary condition
 TTT  = alp*np.pi/180.0
 for j in range(1, Jmax+1):
  u[j, 1] = 0.0
  v[j, 1] = 0.0
  u[j, Kmax] = np.cos(TTT)
  v[j, Kmax] = np.sin(TTT)
        
 if(ityp == 1):
  for k in range(1, Kmax+1):
   u[1, k] = u[Jm, k]
   v[1, k] = v[Jm, k]
   u[Jmax, k] = u[2, k]
   v[Jmax, k] = v[2, k]
 else:
  for k in range(1, Kmax+1):
   u[1, k] = u[2, k]
   v[1, k] = 0.0
   u[Jmax, k] = u[Jm, k]
   v[Jmax, k] = 0.0

 for k in range(Kmax):
  for j in range(Jmax):
   us[j,k,n-1]=u[j,k]
   vs[j,k,n-1]=v[j,k]
   ps[j,k,n-1]=p[j,k]

print("finished")

アニメーションの作成

#make plot(animation)

fig=plt.figure(figsize = (3.5, 3.5))

def update(ns, fig_title, A):

 if i != 0:
  plt.cla()                      # 現在描写されているグラフを消去
    
 tplot=A*dt*(ns)
 u=us[:,:,A*(ns)-1].copy()
 v=vs[:,:,A*(ns)-1].copy()
 z=np.sqrt(u**2+v**2)
 plt.title("Re="+str(Re)+":   "+"t="+str(tplot))
 start, finish = 1, 63
 plt.pcolor(x, y, z)
 plt.xlim(-3, 3)
 plt.ylim(-3, 3)
 plt.tick_params(labelbottom='off')
 plt.tick_params(labelleft='off')
 plt.axes().set_aspect('equal', 'datalim')
    
ani = anm.FuncAnimation(fig, update, fargs = ('Animation ', 40),interval = 1, frames =100 )

ani.save("flow.gif", writer = 'imagemagick')

コードが折りたたまれて見にくい...
色々試したけどスマホ画面では表示がうまく変わってくれなかった。コード見るときはPC推奨。