Pythonによる流体シミュレーション(円柱周りの流れ)
はじめに
使用環境
・MacOS Mojave バージョン 10.14.3
・anaconda3
・Python3.6
・Jupyter notebook
今回は友人からのリクエスト「流体力学」について書いていきます。「流体力学についてなんか書いて!」と言われましても自分自身専門外で、特に面白いネタも持ってないので、円柱周りの流れのシミュレーションをしてみることにしました。
NS式を数値的に解く
流体力学において最も有名な式の1つ「ナビエ・ストークス方程式(NS式)」を数値的に解く方法として、MAC法というのがあります。大雑把にいうと、MAC法は連続の式を満たすような圧力Pを予め決定し(ポアソン方程式)、その値を用いて次のステップの速度を計算する(NS式)といった方法です。詳しい説明は調べたら出てくると思うので省略。実際に使ったコードは最後にあります。
シミュレーション結果
色の薄い部分が速度の早いところ、濃い部分が速度の遅いところです。中心の濃い円が円柱(速度0)を表しています。
何となく後方に渦が2つ伸びていくのがわかる?かな?(微妙)あまり綺麗になりませんでした...もっと細かく区切った方がいいんでしょうがオーバーフローするし、ただでさえ激遅な計算がさらに遅くなってしまいます。
でもまぁ色々勉強になったので満足(アニメーションの描画とか)
コード全文
ライブラリをインポート
%matplotlib nbagg import numpy as np import matplotlib.pyplot as plt import sys,random import matplotlib.animation as anm
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推奨。