CFD解析～ミニカーでドラフティングを行う場合の最適位置とその効果について～

概要

モータースポーツにおいて, 相手の背後について風避けに使うドラフティングという技術は勝負に勝つ上で重要である. スポーツカーやロードバイクにおけるドラフティングについては様々な研究がなされているが, 数cmオーダーでのドラフティングについての研究は少ない. 本研究では, 数センチメートルサイズのミニカー\(2\)台によるドラフティングを模したCFD解析を行い, ドラフティングに適した位置や, 位置による効果の違いについて調べる.

実験方法

仮定

シミュレーションを行うに当たっての仮定は以下の通りである.

仮定1

ミニカー(以下, 車)は全幅\(2cm\), 全長\(4cm\)の角柱であり, 空気の流れは2次元で考える. (車の上を通り抜ける空気の流れはないものとする.)

仮定2

相手(前の車)は高度な運転により, 一様流に対する一定の対気速度・位置を維持するものとする. (すなわち, 座標内で相手の位置は変化しない)

仮定3

自分(後ろの車)は高度な運転により, 相手との相対位置を一定に保つものとする. (すなわち, 座標内で自分の位置は変化しない)

条件設定

レイノルズ数\(Re=70\), クーラン数\(C=0.2\), 緩和法の加速パラメータ\(ω=1.0\)とし, 上流側の初期条件は\(u=1.0, v=0.0, p=0.0\)とする. イメージとしては, 全幅\(2cm\)の車が高度\(0m\)の大気中を\(5cm/s\)で進んでいることを想定している.

評価方法

Fig.1 自分(後ろの車)の9つの位置

Fig.1に示した9つの点に1つ1つ自分(後ろの車)を置いてシミュレーションを行い, 以下の2項目について評価を行う.

相手と自分(地点\(1～9\))の中心位置をそれぞれ格子番号で表すと,

相手\((100,100)\)

地点1\((150,100)\)

地点2\((200,100)\)

地点3\((250,100)\)

地点4\((150,115)\)

地点5\((200,115)\)

地点6\((250,115)\)

地点7\((150,130)\)

地点8\((200,130)\)

地点9\((250,130)\)

評価項目1

相手(前の車)の抵抗係数に対する自分(後ろの車)の抵抗係数の割合を\(R\)とおく. この\(R\)が小さいほど, 相手に対して優位な状態にいることになる.

評価項目2

カルマン渦の影響を受けて, 自分(後ろの車)には振動的な横力が働くことが予想される. この空気力の係数\(C_l\)の振幅を\(A\)とする.この\(A\)が小さいほど横風に煽られることなく運転が楽になる.

使用した言語など

プロセッサは「Intel® Core™ i5-5575R CPU @ 2.80GHz(4CPUs),~2.8GHz」, OSは「Windows10 Enterprise 2015 LTSB 64ビット(10.0, ビルド 10240)」, 言語はCythonを用い, 描画にはMatplotlibを用いた.

解法

圧力場の解析には緩和法を用いたPoisson方程式を, 速度場の解析にはKawamuraスキームを用いた. ソースコードは「7.ソースコード」に記載した.

結果

各地点において, \(n=0～6000[step]\)における圧力分布・速度ベクトル分布を描いた. そして, \(n=0～6000[step]\)までの相手, 自分それぞれの抵抗と横力の係数のグラフも描画した. なお, \(C_{d1}\)は相手の抵抗係数(\(x\)軸方向), \(C_{l1}\)は相手の横力係数(\(y\)軸方向), \(C_{d2}\)は自分の抵抗係数(\(x\)軸方向), \(C_{l2}\)は自分の横力係数(\(y\)軸方向)である. 各地点について, 圧力係数分布のアニメーション・速度ベクトル分布のアニメーション・抵抗係数の推移グラフの順に以下に示した.

地点1

地点1における圧力分布の変化

地点1における速度ベクトル分布の変化

地点1の空気力係数の推移

地点2

地点2における圧力分布の変化

地点2における速度ベクトル分布の変化

地点2の空気力係数の推移

地点3

地点3における圧力分布の変化

地点3における速度ベクトル分布の変化

地点3の空気力係数の推移

地点4

地点4における圧力分布の変化

地点4における速度ベクトル分布の変化

地点4の空気力係数の推移

地点5

地点5における圧力分布の変化

地点5における速度ベクトル分布の変化

地点5の空気力係数の推移

地点6

地点6における圧力分布の変化

地点6における速度ベクトル分布の変化

地点6の空気力係数の推移

地点7

地点7における圧力分布の変化

地点7における速度ベクトル分布の変化

地点7の空気力係数の推移

地点8

地点8における圧力分布の変化

地点8における速度ベクトル分布の変化

地点8の空気力係数の推移

地点9

地点9における圧力分布の変化

地点9における速度ベクトル分布の変化

地点9の空気力係数の推移

各地点の抵抗係数とR, A

ほぼ定常的な振動に移行しきった\(n=4000～6000[step]\)における抵抗係数\(C_{d1},C_{d2}\)の平均値を用いて \(R=\frac{C_{d2ave}}{C_{d1ave}}\) として\(R\)を計算した. また, \(n=4000～6000[step]に\)おける横力の最大値\(C_{l2max}\)と最小値\(C_{l2min}\)を用いて, \(A=\frac{C_{l2max}-C_{l2min}}{2}\)として\(A\)を計算した.

地点	\(C_{d1}\)	\(C_{d2}\)	\(R\)	\(A\)
1	1.196	-0.059	-0.049	0.003
2	1.240	0.680	0.548	2.349
3	1.249	0.739	0.592	2.116
4	1.147	1.178	1.026	0.292
5	1.214	0.866	0.714	1.421
6	1.249	0.896	0.717	1.763
7	1.220	1.313	1.076	0.183
8	1.231	1.318	1.071	1.313
9	1.254	1.254	1.000	1.396

考察

評価項目1

地点1において\(R\)が最も小さくなった. したがって, 9地点の中で相手に対して最も優位な状態にいられるのは地点1である. (しかも, 前へ押されるような空気力が働いている.) 地点4, 5を比較すると, 4の方が\(R\)が大きくなっている. これは, 4の方が相手の進行軸に対して角度があり相手があまり風避けにならないためだと考えられる. これは7, 8, 9の比較においても同様である.

評価項目2

地点1,4,7では\(A\)の値が小さいため,横風に煽られることなく運転が楽である. その他の地点は地点1,4,7と比べて相手から遠い位置にあり, 相手から発生するカルマン渦が十分に発達して襲い掛かってくるため, 大きな横風を受けることになるのだと考えられる.

ソースコード

%%cython
from __future__ import division
import numpy as np
cimport numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.animation as animation

DTYPE = np.double
ctypedef np.double_t DTYPE_t

# Change Here to change flow conditions------------------
#flow conditions
cdef double re = 70.0 #reynolds number
cdef double cfl = 0.2 #CFL number クーラン数
cdef double V=1.0 #一様流速
#SOR params for p
cdef double omegap = 1.00 #1 for gaus seidel
cdef int maxitp = 100   #ポアソン方程式における圧力解の最大更新回数
cdef double errorp = 1.0e-4
#No. of time steps
cdef int nlast = 6001   #6000回
#monitor interval
cdef int nlp = 10     #何回ごとに書き出すか


#grid
#set x-grid parameters
cdef int mx = 401   #グリッドのx範囲
cdef int i1 = 90    #角柱の前面位置
cdef int i2 = 110   #角柱の後面位置
cdef int k1 = 240   #後ろの角柱の前面位置
cdef int k2 = 260   #後ろの角柱の後面位置
#set y-grid params
cdef int my = 200   #グリッドのy範囲
cdef int j1 = 95    #角柱の左面位置
cdef int j2 = 105   #角柱の右面位置
cdef int l1 = 125    #後ろの角柱の左面位置
cdef int l2 = 135   #後ろの角柱の右面位置
cdef double dx = 1.0 / (j2-j1)  #角柱のy方向の長さを1とした時の格子幅の長さ
cdef double dy = 1.0 / (j2-j1)  #角柱のy方向の長さを1とした時の格子幅の長さ
#------------------------------------------------------

#timestep
cdef double dt = cfl*min(dx,dy)/ V  #時間間隔

def setgrd(x,y):   #角柱中心を原点とした時のx,y座標(x[i,j], y[i,j])を生成する関数
    """物体を中心とした座標の算出"""
    icent = int((i1+i2)/2)  #角柱の中心位置x
    jcent = int((j1+j2)/2)  #角柱の中心位置y
    for i in range(mx):
        for j in range(my):
            x[i,j] = dx*(i-icent)    #角柱中心を原点とした時のx座標(角柱のx軸に平行な1辺を1とする)
            y[i,j] = dy*(j-jcent)    #角柱中心を原点とした時のy座標(角柱のx軸に平行な1辺を1とする)
    return

def slvflw():  #流れを解く関数
    """初期設定，流れのシミュレーションおよび出力を行う。"""
    #set initial conditions
    cdef int nbegin = 0
    cdef int nstep = 0
    cdef double time = 0.0
    cdef DTYPE_t cd, cd2, cl, cl2, cpfore, cpfore2, cpback, cpback2, cpbtm, cpbtm2, cptop, cptop2, cp1, cp2
    cdef int n, i, j, k, l, itrp
    cdef double resp
    #axis
    cdef np.ndarray[DTYPE_t, ndim=2] x = np.zeros([mx,my], dtype=DTYPE)
    cdef np.ndarray[DTYPE_t, ndim=2] y = np.zeros_like(x)
    #computation result with initial conditions
    cdef np.ndarray[DTYPE_t, ndim=2] u = np.ones([mx+4,my+1], dtype=DTYPE)
    cdef np.ndarray[DTYPE_t, ndim=2] v = np.zeros_like(u)
    cdef np.ndarray[DTYPE_t, ndim=2] p = np.zeros_like(u)
    CD1 = np.zeros(nlast,dtype=DTYPE)
    CL1 = np.zeros(nlast,dtype=DTYPE)
    CD2 = np.zeros(nlast,dtype=DTYPE)
    CL2 = np.zeros(nlast,dtype=DTYPE)

    setgrd(x,y) #グリッドを設定

    bcforp(p)   #圧力に関する境界条件を設定
    bcforv(u,v) #速度に関する境界条件を設定

    #start monitoring
    print('Step \t Res(p) \t itr \t CD1 \t CL1 \t CD2 \t CL2')  #書き出し
    # time marching
    for n in range(nlast):  #n=0～5000
        nstep = n+nbegin    #今のステップがいくつか
        time = time+dt      #今の時間を更新
        #solve poisson for p
        resp, itrp = poiseq(p, u, v)
        bcforp(p)           #圧力に関する境界条件を設定　
        #update u, v
        veloeq(p, u, v)
        bcforv(u,v)         #速度に関する境界条件を設定

        #calculate CL and CD
        cd = 0.0
        for j in range(j1, j2+1):
            cpfore = (2*p[i1,j] + 2*p[i1,j+1])/2  #角柱前面のCp(=2p)を計算
            cpback = (2*p[i2,j] + 2*p[i2,j+1])/2  #角柱後面のCp(=2p)を計算
            cd = cd+(cpfore-cpback)*dy  #合計のCdを算出
        cl = 0.0
        for i in range(i1, i2+1):
            cpbtm = (2*p[i,j1] + 2*p[i+1,j1])/2
            cptop = (2*p[i,j2] + 2*p[i+1,j2])/2
            cl = cl+(cpbtm-cptop)*dx   #合計のClを算出
        cd2 = 0.0
        for l in range(l1, l2+1):
            cpfore2 = (2*p[k1,l] + 2*p[k1,l+1])/2  #後ろの角柱前面のCp(=2p)を計算
            cpback2 = (2*p[k2,l] + 2*p[k2,l+1])/2  #後ろの角柱後面のCp(=2p)を計算
            cd2 = cd2+(cpfore2-cpback2)*dy  #合計のCdを算出
        cl2 = 0.0
        for k in range(k1, k2+1):
            cpbtm2 = (2*p[k,l1] + 2*p[k+1,l1])/2
            cptop2 = (2*p[k,l2] + 2*p[k+1,l2])/2
            cl2 = cl2+(cpbtm2-cptop2)*dx   #合計のClを算出
        CD1[n]=cd
        CL1[n]=cl
        CD2[n]=cd2
        CL2[n]=cl2

        #monitor by nlp steps
        if (n%nlp==0):
            cp1 = 2*p[i2+i2-i1, j1]
            cp2 = 2*p[i2+i2-i1, j2]
            print('%d \t %f \t %d \t %f \t %f \t %f \t %f '%(nstep,resp,itrp,cd,cl,cd2,cl2))
        if (n==0):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,0)
        if (n==150):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,150)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==300):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,300)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==450):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,450)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==600):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,600)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==750):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,750)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==900):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,900)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==1050):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,1050)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==1200):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,1200)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==1350):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,1350)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==1500):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,1500)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==1650):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,1650)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==1800):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,1800)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==1950):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,1950)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==2100):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,2100)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==2250):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,2250)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==2400):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,2400)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==2550):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,2550)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==2700):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,2700)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==2850):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,2850)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==3000):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,3000)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==3150):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,3150)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==3300):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,3300)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==3450):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,3450)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==3600):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,3600)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==3750):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,3750)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==3900):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,3900)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==4050):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,4050)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==4200):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,4200)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==4350):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,4350)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==4500):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,4500)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==4650):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,4650)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==4800):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,4800)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==4950):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,4950)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==5100):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,5100)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==5250):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,5250)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==5400):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,5400)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==5550):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,5550)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==5700):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,5700)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==5850):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,5850)
            CDCL(CD1,CL1,CD2,CL2)
        if (n==6000):
            np.savez('result_puv.npz',p=p,u=u,v=v)
            draw(p,u,v,6000)
            plt.plot(CD1,label="Cd1")
            plt.plot(CL1,label="Cl1")
            plt.plot(CD2,label="Cd2")
            plt.plot(CL2,label="Cl2")
            plt.legend()
            plt.title("Cd1,Cl1,Cd2,Cl2")
            plt.xlabel("n[step]")
            plt.ylabel("Cd1,Cl1,Cd2,Cl2[]")
            plt.show()
            print(np.average(CD1[4000:6000]),np.average(CL1[4000:6000]),np.average(CD2[4000:6000]),np.average(CL2[4000:6000]))
            print(np.average(CD2[4000:6000])/np.average(CD1[4000:6000]))
            print((np.max(CL2[4000:6000])-np.min(CL2[4000:6000]))*0.5)


def bcforp(np.ndarray[DTYPE_t, ndim=2] p):    #圧力に関する境界条件を設定する関数
    """pの境界条件の設定"""
    cdef int i,j

    for j in range(my+1):   #p[0,0]～p[0,my]まで圧力はゼロで, p[mx+2,0]～p[mx+2,my]まで圧力はゼロ
        #inflow condition
        p[0,j] = 0
        #downstream condition (i=mx+2)
        p[mx+2,j] = 0
    for i in range(mx+3):   #p[0,0]～p[mx+3,0]まで圧力はゼロで, p[0,my]～p[mx+2,my]まで圧力はゼロ
        #bottom condition (j=0)
        p[i,0] = 0
        #upper condition (j=my)
        p[i,my] = 0
    #wall condition    #壁面付近では圧力勾配がゼロ
    #4 corners         #4隅ではそこから斜めに離れた点と圧力は変わらない
    p[i1,j1] = p[i1-1,j1-1]
    p[i1,j2] = p[i1-1,j2+1]
    p[i2,j1] = p[i2+1,j1-1]
    p[i2,j2] = p[i2+1,j2+1]
    #sides             4辺では一つ離れた点と圧力は変わらない
    for j in range(j1+1,j2):
        p[i1,j] = p[i1-1,j]
        p[i2,j] = p[i2+1,j]
    for i in range(i1+1,i2):
        p[i,j1] = p[i,j1-1]
        p[i,j2] = p[i,j2+2]
    #後ろの角柱
    #4 corners         #4隅ではそこから斜めに離れた点と圧力は変わらない
    p[k1,l1] = p[k1-1,l1-1]
    p[k1,l2] = p[k1-1,l2+1]
    p[k2,l1] = p[k2+1,l1-1]
    p[k2,l2] = p[k2+1,l2+1]
    #sides             4辺では一つ離れた点と圧力は変わらない
    for l in range(l1+1,l2):
        p[k1,l] = p[k1-1,l]
        p[k2,l] = p[k2+1,l]
    for k in range(k1+1,k2):
        p[k,l1] = p[k,l1-1]
        p[k,l2] = p[k,l2+2]
    return

def bcforv(np.ndarray[DTYPE_t, ndim=2] u, np.ndarray[DTYPE_t, ndim=2] v):    #速度に関する境界条件を設定する関数
    """u, vの境界条件の設定"""
    cdef int i, j

    for j in range(my+1):
        #inflow condition
        u[1,j] = V           #左境界ではx方向速度がV
        v[1,j] = 0           #左境界ではy方向速度が0
        u[0,j] = V           #左境界ではx方向速度がV
        v[0,j] = 0           #左境界ではy方向速度が0
        #downstream condition
        u[mx+2,j] = 2*u[mx+1,j]-u[mx,j]   
        v[mx+2,j] = 2*v[mx+1,j]-v[mx,j]
        u[mx+3,j] = 2*u[mx+2,j]-u[mx+1,j]
        v[mx+3,j] = 2*v[mx+2,j]-v[mx+1,j]
    #bottom conditions
    for i in range(mx+4):
        u[i,1] = 2*u[i,2]-u[i,3]
        v[i,1] = 2*v[i,2]-v[i,3]
        u[i,0] = 2*u[i,1]-u[i,2]
        v[i,0] = 2*v[i,1]-v[i,2]
        u[i,my-1] = 2*u[i,my-2]-u[i,my-3]   #upper condition
        v[i,my-1] = 2*v[i,my-2]-v[i,my-3]
        u[i,my] = 2*u[i,my-1]-u[i,my-2]
        v[i,my] = 2*v[i,my-1]-v[i,my-2]
    #wall condition   #壁面(角柱内)では速度ゼロ
    for i in range(i1,i2+1):
        for j in range(j1,j2+1):
            u[i,j] = 0
            v[i,j] = 0
    #後ろの角柱について
    #wall condition   #壁面(角柱内)では速度ゼロ
    for k in range(k1,k2+1):
        for l in range(l1,l2+1):
            u[k,l] = 0
            v[k,l] = 0
    return

def poiseq(np.ndarray[DTYPE_t, ndim=2] p, np.ndarray[DTYPE_t, ndim=2] u, np.ndarray[DTYPE_t, ndim=2] v):
    """緩和法によるポアソン方程式の解法"""
    cdef np.ndarray[DTYPE_t, ndim=2] rhs = np.zeros_like(p)
    cdef double ux, uy, vx,vy, res, dp
    cdef int i,j, itr, itrp

    #右辺の計算
    for i in range(2,mx+2):
        for j in range(2,my-1):
            if (i<i1 or i>i2 or j<j1 or j>j2)and(i<k1 or i>k2 or j<l1 or j>l2):  #角柱外部の場合
                ux = (u[i+1,j]-u[i-1,j]) / (2*dx) #uのx偏微分
                uy = (u[i,j+1]-u[i,j-1]) / (2*dy) #uのy偏微分
                vx = (v[i+1,j]-v[i-1,j]) / (2*dx) #vのx偏微分
                vy = (v[i,j+1]-v[i,j-1]) / (2*dy) #vのy偏微分
                rhs[i,j] = (ux+vy)/dt - (ux**2+2*uy*vx+vy**2)  #Right Hand Sideの計算
    #iterations
    for itr in range(maxitp):
        res = 0.0 #残差
        #緩和法
        for i in range(2,mx+2):
            for j in range(2,my-1):
                if (i<i1 or i>i2 or j<j1 or j>j2)and(i<k1 or i>k2 or j<l1 or j>l2):
                    dp = (p[i+1,j]+p[i-1,j])/dx**2 + (p[i,j+1]+p[i,j-1])/dy**2 - rhs[i,j]
                    dp = dp/(2/dx**2+2/dy**2) - p[i,j]   #この点における残差を求めた
                    res = res + dp**2 #残差    #各点の残差の2乗和を一旦resとする.
                    p[i,j] = p[i,j]+omegap*dp #緩和して修正 p[i,j]をアップデート
        #境界条件
        bcforp(p) #新たにもう一度圧力に関する境界条件を設定
        res = np.sqrt(res/mx/my) #正しい残差にする
        itrp = itr   #現在のイタレーション回数を入れておく
        if (res<errorp): #残差が許容範囲内になったら計算終了
            break

    return res, itrp  #残差とイタレーション回数を出力

def veloeq(np.ndarray[DTYPE_t, ndim=2] p, np.ndarray[DTYPE_t, ndim=2] u, np.ndarray[DTYPE_t, ndim=2] v):    #圧力を踏まえたうえで, 速度u[i,j],v[i,j]を求める関数
    """Kawamura スキームを用いた移流方程式の解法"""
    cdef np.ndarray[DTYPE_t, ndim=2] urhs = np.zeros_like(p)
    cdef np.ndarray[DTYPE_t, ndim=2] vrhs = np.zeros_like(p)
    cdef int i, j

    #pressure gradient   #圧力勾配の項について計算
    for i in range(2,mx+2):
        for j in range(2,my-1):
            if (i<i1 or i>i2 or j<j1 or j>j2)and(i<k1 or i>k2 or j<l1 or j>l2):
                urhs[i,j] = -(p[i+1,j]-p[i-1,j]) / (2*dx)
                vrhs[i,j] = -(p[i,j+1]-p[i,j-1]) / (2*dy)
    #viscous term        #粘性項を加算
    for i in range(2,mx+2):
        for j in range(2,my-1):
            if (i<i1 or i>i2 or j<j1 or j>j2)and(i<k1 or i>k2 or j<l1 or j>l2):
                urhs[i,j] += (u[i+1,j]-2*u[i,j]+u[i-1,j])/(re*dx**2) + (u[i,j+1]-2*u[i,j]+u[i,j-1])/(re*dy**2)
                vrhs[i,j] += (v[i+1,j]-2*v[i,j]+v[i-1,j])/(re*dx**2) + (v[i,j+1]-2*v[i,j]+v[i,j-1])/(re*dy**2)
    #Advection term in x direction with Kawamura scheme
    #物体内部を外挿
    for j in range(j1+1,j2):
        u[i1+1,j] = 2*u[i1,j]-u[i1-1,j]
        u[i2-1,j] = 2*u[i2,j]-u[i2+1,j]
        v[i1+1,j] = 2*v[i1,j]-v[i1-1,j]
        v[i2-1,j] = 2*v[i2,j]-v[i2+1,j]
    #後ろの角柱について
    #物体内部を外挿
    for l in range(l1+1,l2):
        u[k1+1,l] = 2*u[k1,l]-u[k1-1,l]
        u[k2-1,l] = 2*u[k2,l]-u[k2+1,l]
        v[k1+1,l] = 2*v[k1,l]-v[k1-1,l]
        v[k2-1,l] = 2*v[k2,l]-v[k2+1,l]
    for i in range(2,mx+2):
        for j in range(2,my-1):
            if (i<i1 or i>i2 or j<j1 or j>j2)and(i<k1 or i>k2 or j<l1 or j>l2):
                urhs[i,j] += -u[i,j]*(-u[i+2,j]+8*(u[i+1,j]-u[i-1,j])+u[i-2,j])/(12*dx) - abs(u[i,j])*(u[i+2,j]-4*u[i+1,j]+6*u[i,j]-4*u[i-1,j]+u[i-2,j])/(4*dx)
                vrhs[i,j] += -u[i,j]*(-v[i+2,j]+8*(v[i+1,j]-v[i-1,j])+v[i-2,j])/(12*dx) - abs(u[i,j])*(v[i+2,j]-4*v[i+1,j]+6*v[i,j]-4*v[i-1,j]+v[i-2,j])/(4*dx)
    #Advection term in y direction with Kawamura scheme
    #物体内部を外挿  #物体内部の速度を外の速度から決める
    for i in range(i1+1,i2):
        u[i,j1+1] = 2*u[i,j1]-u[i,j1-1]
        u[i,j2-1] = 2*u[i,j2]-u[i,j2+1]
        v[i,j1+1] = 2*v[i,j1]-v[i,j1-1]
        v[i,j2-1] = 2*v[i,j2]-v[i,j2+1]
    #物体内部を外挿
    for k in range(k1+1,k2):
        u[k,l1+1] = 2*u[k,l1]-u[k,l1-1]
        u[k,l2-1] = 2*u[k,l2]-u[k,l2+1]
        v[k,l1+1] = 2*v[k,l1]-v[k,l1-1]
        v[k,l2-1] = 2*v[k,l2]-v[k,l2+1]
    for i in range(2,mx+2):
        for j in range(2,my-1):
            if (i<i1 or i>i2 or j<j1 or j>j2)and(i<k1 or i>k2 or j<l1 or j>l2):
                urhs[i,j] += -v[i,j]*(-u[i,j+2]+8*(u[i,j+1]-u[i,j-1])+u[i,j-2])/(12*dx) - abs(v[i,j])*(u[i,j+2]-4*u[i,j+1]+6*u[i,j]-4*u[i,j-1]+u[i,j-2])/(4*dx)
                vrhs[i,j] += -v[i,j]*(-v[i,j+2]+8*(v[i,j+1]-v[i,j-1])+v[i,j-2])/(12*dy) - abs(v[i,j])*(v[i,j+2]-4*v[i,j+1]+6*v[i,j]-4*v[i,j-1]+v[i,j-2])/(4*dy)
    #update
    for i in range(2,mx+2):
        for j in range(2,my-1):
            if (i<i1 or i>i2 or j<j1 or j>j2)and(i<k1 or i>k2 or j<l1 or j>l2):
                u[i,j] = u[i,j]+dt*urhs[i,j]
                v[i,j] = v[i,j]+dt*vrhs[i,j]

    return

def minicar():  #ミニカーの長方形を描画する
    ax=plt.axes()
    r1=patches.Rectangle(xy=(i1,j1),width=(i2-i1),height=(j2-j1),ec='#000000',fill=False)
    r2=patches.Rectangle(xy=(k1,l1),width=(k2-k1),height=(l2-l1),ec='#000000',fill=False)
    ax.add_patch(r1)
    ax.add_patch(r2)
    plt.axis('scaled')
    ax.set_aspect('equal')


def draw(p,u,v,n):
    minicar()
    result = np.load('result_puv.npz')
    plt.imshow(result['p'].T)  #0ステップ目の圧力を描画
    plt.colorbar(shrink=.92)
    plt.xlim(0,400)
    plt.ylim(0,200)
    plt.gca().set_aspect('equal', adjustable='box')
    plt.savefig("result_p_"+str(n)+".png",format='png',dpi=300)
    plt.close()
 
    minicar()
    plt.imshow(result['u'].T) #0ステップ目のx方向速度を描画
    plt.colorbar(shrink=.92)
    plt.xlim(0,400)
    plt.ylim(0,200)  
    plt.gca().set_aspect('equal', adjustable='box')
    plt.savefig("result_u_"+str(n)+".png",format='png',dpi=300)
    plt.close()

    minicar()
    plt.imshow(result['v'].T) #0ステップ目のy方向速度を描画
    plt.colorbar(shrink=.92)
    plt.xlim(0,400)
    plt.ylim(0,200)
    plt.gca().set_aspect('equal', adjustable='box')
    plt.savefig("result_v_"+str(n)+".png",format='png',dpi=300)
    plt.close()

    x=range(0,401)
    y=range(0,201)
    X,Y=np.meshgrid(x,y)
    Z=p[X,Y]
    plt.contour(X,Y,Z)
    minicar()
    plt.colorbar(shrink=.92)
    plt.gca().set_aspect('equal',adjustable='box')
    plt.savefig("result_isobars_"+str(n)+".png",format='png',dpi=300)
    plt.close()
  
    X=np.zeros((80,41))
    Y=np.zeros((80,41))
    U=np.zeros((80,41))
    V=np.zeros((80,41))
    plt.figure()
    for i in range(0,80):
        for j in range(0,41):
            X[i,j]=5*i
            Y[i,j]=5*j
            U[i,j]=u[5*i,5*j]
            V[i,j]=v[5*i,5*j]
    plt.quiver(X,Y,U,V,angles='xy',scale_units='xy',scale=0.1)
    minicar()  
    plt.xlim(0,400)
    plt.ylim(0,200)
    plt.grid()
    plt.draw()
    plt.savefig("result_vectors_"+str(n)+".png",format='png',dpi=300)
    plt.close()
    
    
def CDCL(CD1,CL1,CD2,CL2):
    plt.plot(CD1,label="Cd1")
    plt.plot(CL1,label="Cl1")
    plt.plot(CD2,label="Cd2")
    plt.plot(CL2,label="Cl2")
    plt.legend()
    plt.title("Cd1,Cl1,Cd2,Cl2")
    plt.xlabel("n[step]")
    plt.ylabel("Cd1,Cl1,Cd2,Cl2[]")
    plt.close()




slvflw()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

%%cython

from __future__ import division

import numpy as np

cimport numpy as np

import matplotlib.pyplot as plt

import matplotlib.patches as patches

import matplotlib.animation as animation

DTYPE = np.double

ctypedef np.double_t DTYPE_t

# Change Here to change flow conditions------------------

#flow conditions

cdef double re = 70.0 #reynolds number

cdef double cfl = 0.2 #CFL number クーラン数

cdef double V=1.0 #一様流速

#SOR params for p

cdef double omegap = 1.00 #1 for gaus seidel

cdef int maxitp = 100 #ポアソン方程式における圧力解の最大更新回数

cdef double errorp = 1.0e-4

#No. of time steps

cdef int nlast = 6001 #6000回

#monitor interval

cdef int nlp = 10 #何回ごとに書き出すか

#grid

#set x-grid parameters

cdef int mx = 401 #グリッドのx範囲

cdef int i1 = 90 #角柱の前面位置

cdef int i2 = 110 #角柱の後面位置

cdef int k1 = 240 #後ろの角柱の前面位置

cdef int k2 = 260 #後ろの角柱の後面位置

#set y-grid params

cdef int my = 200 #グリッドのy範囲

cdef int j1 = 95 #角柱の左面位置

cdef int j2 = 105 #角柱の右面位置

cdef int l1 = 125 #後ろの角柱の左面位置

cdef int l2 = 135 #後ろの角柱の右面位置

cdef double dx = 1.0 / (j2-j1) #角柱のy方向の長さを1とした時の格子幅の長さ

cdef double dy = 1.0 / (j2-j1) #角柱のy方向の長さを1とした時の格子幅の長さ

#------------------------------------------------------

#timestep

cdef double dt = cfl*min(dx,dy)/ V #時間間隔

def setgrd(x,y): #角柱中心を原点とした時のx,y座標(x[i,j], y[i,j])を生成する関数

"""物体を中心とした座標の算出"""

icent = int((i1+i2)/2) #角柱の中心位置x

jcent = int((j1+j2)/2) #角柱の中心位置y

for i in range(mx):

for j in range(my):

x[i,j] = dx*(i-icent) #角柱中心を原点とした時のx座標(角柱のx軸に平行な1辺を1とする)

y[i,j] = dy*(j-jcent) #角柱中心を原点とした時のy座標(角柱のx軸に平行な1辺を1とする)

return

def slvflw(): #流れを解く関数

"""初期設定，流れのシミュレーションおよび出力を行う。"""

#set initial conditions

cdef int nbegin = 0

cdef int nstep = 0

cdef double time = 0.0

cdef DTYPE_t cd, cd2, cl, cl2, cpfore, cpfore2, cpback, cpback2, cpbtm, cpbtm2, cptop, cptop2, cp1, cp2

cdef int n, i, j, k, l, itrp

cdef double resp

#axis

cdef np.ndarray[DTYPE_t, ndim=2] x = np.zeros([mx,my], dtype=DTYPE)

cdef np.ndarray[DTYPE_t, ndim=2] y = np.zeros_like(x)

#computation result with initial conditions

cdef np.ndarray[DTYPE_t, ndim=2] u = np.ones([mx+4,my+1], dtype=DTYPE)

cdef np.ndarray[DTYPE_t, ndim=2] v = np.zeros_like(u)

cdef np.ndarray[DTYPE_t, ndim=2] p = np.zeros_like(u)

CD1 = np.zeros(nlast,dtype=DTYPE)

CL1 = np.zeros(nlast,dtype=DTYPE)

CD2 = np.zeros(nlast,dtype=DTYPE)

CL2 = np.zeros(nlast,dtype=DTYPE)

setgrd(x,y) #グリッドを設定

bcforp(p) #圧力に関する境界条件を設定

bcforv(u,v) #速度に関する境界条件を設定

#start monitoring

print('Step \t Res(p) \t itr \t CD1 \t CL1 \t CD2 \t CL2') #書き出し

# time marching

for n in range(nlast): #n=0～5000

nstep = n+nbegin #今のステップがいくつか

time = time+dt #今の時間を更新

#solve poisson for p

resp, itrp = poiseq(p, u, v)

bcforp(p) #圧力に関する境界条件を設定　

#update u, v

veloeq(p, u, v)

bcforv(u,v) #速度に関する境界条件を設定

#calculate CL and CD

cd = 0.0

for j in range(j1, j2+1):

cpfore = (2*p[i1,j] + 2*p[i1,j+1])/2 #角柱前面のCp(=2p)を計算

cpback = (2*p[i2,j] + 2*p[i2,j+1])/2 #角柱後面のCp(=2p)を計算

cd = cd+(cpfore-cpback)*dy #合計のCdを算出

cl = 0.0

for i in range(i1, i2+1):

cpbtm = (2*p[i,j1] + 2*p[i+1,j1])/2

cptop = (2*p[i,j2] + 2*p[i+1,j2])/2

cl = cl+(cpbtm-cptop)*dx #合計のClを算出

cd2 = 0.0

for l in range(l1, l2+1):

cpfore2 = (2*p[k1,l] + 2*p[k1,l+1])/2 #後ろの角柱前面のCp(=2p)を計算

cpback2 = (2*p[k2,l] + 2*p[k2,l+1])/2 #後ろの角柱後面のCp(=2p)を計算

cd2 = cd2+(cpfore2-cpback2)*dy #合計のCdを算出

cl2 = 0.0

for k in range(k1, k2+1):

cpbtm2 = (2*p[k,l1] + 2*p[k+1,l1])/2

cptop2 = (2*p[k,l2] + 2*p[k+1,l2])/2

cl2 = cl2+(cpbtm2-cptop2)*dx #合計のClを算出

CD1[n]=cd

CL1[n]=cl

CD2[n]=cd2

CL2[n]=cl2

#monitor by nlp steps

if (n%nlp==0):

cp1 = 2*p[i2+i2-i1, j1]

cp2 = 2*p[i2+i2-i1, j2]

print('%d \t %f \t %d \t %f \t %f \t %f \t %f '%(nstep,resp,itrp,cd,cl,cd2,cl2))

if (n==0):