趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

SIMPLE法でもpythonならシンプルに書けるのか ⇒ かなりシンプルになるがやっぱり...(ソースコードあり)

最近科学技術計算を書いてなかったのでウォーミングアップにこんな感じの3次元キャビテーション流れをSIMPLE法で書いてみました。

f:id:SedimentHydraulics:20200610105620p:plain

f:id:SedimentHydraulics:20200610093542p:plain

大きいものはこちら


はじめに

  • SIMPLE法って複雑ですよね。そこで、pythonなら少しはわかりやすく書けるかと思いコードを書いてみました。
  • 参考になるソースコードはないかとググったけど意外とヒットしなくて仕方なく一から書きました。

計算モデル

単純な3次元のキャビテーション流れです。

計算条件

  • 計算範囲は1 x 1 x 1の直方体
  • dx=dy=dz=0.05、メッシュ数は20 x 20 x 20の8000メッシュ
  • 計算スキームは、スタッガードスキームのSIMPLE法、移流項は一次風上
  • 上面(z=max)以外はノンスリップ条件
  • 上面(z=max)はx方向流速として1を与える
  • 外力項は重力も含めて考慮しない
  • 非定常項はキャンセル(定常場を求める)
  • レイノルズ数は100

ソースコード

仮の速度の計算式
  • 3本の運動方程式をaxisで切り替えて1つの関数にしています。私の定番の書き方です。
  • 収束計算は不足緩和で係数は標準値を採用。ここは厳密じゃなくても良い。
  • 行数は150。運動方程式がこれだけで書けるのは立派だと思う。
  • スタッガードスキームつらい。コロケートならもっとスッキリするはず。
def calupre(p,u,v,w,Vdir,Vbound,alphaU,dx,dy,dz,nu,axis):
# x: axis = 0, y: axis = 1, z: axis = 2
    
    zero = float(0.0)
    V0 = np.copy(Vdir)
    
    nx, ny, nz = p.shape
    
    ac  = np.zeros_like(Vdir)
    axp = np.zeros_like(Vdir)
    axm = np.zeros_like(Vdir)
    ayp = np.zeros_like(Vdir)
    aym = np.zeros_like(Vdir)
    azp = np.zeros_like(Vdir)
    azm = np.zeros_like(Vdir)
    b   = np.zeros_like(Vdir)
    
    act  = np.zeros_like(Vdir)
    
    if axis == 0:
        xs, ys, zs = 1, 0, 0
        xmb, xpb, ymb, ypb, zmb, zpb = -1, -1, 0, ny-1, 0, nz-1
    elif axis == 1:
        xs, ys, zs = 0, 1, 0
        xmb, xpb, ymb, ypb, zmb, zpb = 0, nx-1, -1, -1, 0, nz-1
    elif axis == 2:
        xs, ys, zs = 0, 0, 1
        xmb, xpb, ymb, ypb, zmb, zpb = 0, nx-1, 0, ny-1, -1, -1
        
    for i in range(xs, nx):
        for j in range(ys, ny):
            for k in range(zs, nz):
                c = (i,j,k)
                
                if axis == 0:
                    uhxp = 0.5*u[i,j,k]   + 0.5*u[i+1,j,k]
                    uhxm = 0.5*u[i,j,k]   + 0.5*u[i-1,j,k]
                    vhyp = 0.5*v[i,j+1,k] + 0.5*v[i-1,j+1,k]
                    vhym = 0.5*v[i,j  ,k] + 0.5*v[i-1,j  ,k]
                    whzp = 0.5*w[i,j,k+1] + 0.5*w[i-1,j,k+1]
                    whzm = 0.5*w[i,j  ,k] + 0.5*w[i-1,j,k]
                elif axis == 1:
                    uhxp = 0.5*u[i+1,j,k] + 0.5*u[i+1,j-1,k]
                    uhxm = 0.5*u[i,j,k]   + 0.5*u[i,j-1,k]
                    vhyp = 0.5*v[i,j,k]   + 0.5*v[i,j+1,k]
                    vhym = 0.5*v[i,j,k]   + 0.5*v[i,j-1,k]
                    whzp = 0.5*w[i,j,k+1] + 0.5*w[i,j-1,k+1]
                    whzm = 0.5*w[i,j  ,k] + 0.5*w[i,j-1,k]
                elif axis == 2:
                    uhxp = 0.5*u[i+1,j,k] + 0.5*u[i+1,j,k-1]
                    uhxm = 0.5*u[i,j,k]   + 0.5*u[i,j,k-1]
                    vhyp = 0.5*v[i,j+1,k] + 0.5*v[i,j+1,k-1]
                    vhym = 0.5*v[i,j,k] + 0.5*v[i,j,k-1]
                    whzp = 0.5*w[i,j,k] + 0.5*w[i,j,k+1]
                    whzm = 0.5*w[i,j,k] + 0.5*w[i,j,k-1]
                
                axp[c] = (max([-uhxp,0.0]) + nu/dx)*dy*dz 
                axm[c] = (max([ uhxm,0.0]) + nu/dx)*dy*dz 
                ayp[c] = (max([-vhyp,0.0]) + nu/dy)*dx*dz 
                aym[c] = (max([ vhym,0.0]) + nu/dy)*dx*dz 
                azp[c] = (max([-whzp,0.0]) + nu/dz)*dx*dy 
                azm[c] = (max([ whzm,0.0]) + nu/dz)*dx*dy 
                
                ac[c] = axp[c]+axm[c]+ayp[c]+aym[c]+azp[c]+azm[c]
                
                if axis == 0:
                    b[c] = -dy*dz*( p[i,j,k] - p[i-1,j,k] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c]
                elif axis == 1:
                    b[c] = -dx*dz*( p[i,j,k] - p[i,j-1,k] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c]
                elif axis == 2:
                    b[c] = -dx*dy*( p[i,j,k] - p[i,j,k-1] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c]
                    
                ac[c] /= alphaU
                act[c] = ac[c]
                
                if i == xmb :
                    ac[c] += axm[c]
                elif i == xpb :
                    ac[c] += axp[c]
                
                if j == ymb :
                    ac[c] += aym[c]
                elif j == ypb :
                    ac[c] += ayp[c]
                    
                if k == zmb :
                    ac[c] += azm[c]
                elif k == zpb :
                    ac[c] += azp[c]
                    b[c] += 2.0*azp[c]*Vbound
    
    def sor(u):
        u0 = np.copy(u) # deep copy
        for i in range(xs, nx):
            for j in range(ys, ny):
                for k in range(zs, nz):
                    c = (i,j,k)
                    xp,xm,yp,ym,zp,zm = (i+1,j,k),(i-1,j,k),(i,j+1,k),(i,j-1,k),(i,j,k+1),(i,j,k-1)
                    
                    if i == xmb :
                        uxp ,uxm = u0[xp], zero
                    elif i == xpb :
                        uxp ,uxm = zero, u0[xm]
                    else:
                        uxp ,uxm = u0[xp], u0[xm]
                    
                    if j == ymb :
                        uyp ,uym = u0[yp], zero
                    elif j == ypb :
                        uyp ,uym = zero , u0[ym]
                    else:
                        uyp ,uym = u0[yp], u0[ym]
                        
                    if k == zmb :
                        uzp, uzm = u0[zp], zero
                    elif k == zpb :
                        uzp, uzm = zero, u0[zm]
                    else:
                        uzp ,uzm = u0[zp], u0[zm]
                    
                    ut = (axp[c]*uxp + axm[c]*uxm \
                        + ayp[c]*uyp + aym[c]*uym \
                        + azp[c]*uzp + azm[c]*uzm + b[c])/ac[c]
                    
                    u0[c] = ut
                    
        return u0
                    
    unew = np.copy(Vdir) # deep copy
    
    for _ in range(100):
        utmp = sor(unew)
        err = np.abs(unew - utmp).max() 
        unew = np.copy(utmp)
        if err < 10**(-3) : break
        
#     print(err)
        
    return unew, act
圧力補正値の計算式
  • 収束計算は厳密にSOR
  • こちらは式も単純なのですっきりと。上の運動方程式と同じく連立一次方程式を解く形なので上手く書けば一つの関数になるかも。
def caldp(p,u,v,w,acu,acv,acw,dx,dy,dz):
    
    zero = float(0.0)
    dp  = np.zeros_like(p)
    ac  = np.zeros_like(p)
    axp = np.zeros_like(p)
    axm = np.zeros_like(p)
    ayp = np.zeros_like(p)
    aym = np.zeros_like(p)
    azp = np.zeros_like(p)
    azm = np.zeros_like(p)
    b   = np.zeros_like(p)
    
    nx, ny, nz = p.shape
    
    for i in range(nx):
        for j in range(ny):
            for k in range(nz):
                c = (i, j, k)
                axp[c] =  zero if np.abs( acu[i+1,j,  k  ] ) < 10**(-8) else dy**2*dz**2/acu[i+1,j,  k  ] 
                axm[c] =  zero if np.abs( acu[i  ,j,  k  ] ) < 10**(-8) else dy**2*dz**2/acu[i  ,j,  k  ] 
                ayp[c] =  zero if np.abs( acv[i  ,j+1,k  ] ) < 10**(-8) else dx**2*dz**2/acv[i  ,j+1,k  ] 
                aym[c] =  zero if np.abs( acv[i  ,j,  k  ] ) < 10**(-8) else dx**2*dz**2/acv[i  ,j,  k  ] 
                azp[c] =  zero if np.abs( acw[i  ,j,  k+1] ) < 10**(-8) else dx**2*dy**2/acw[i  ,j,  k+1] 
                azm[c] =  zero if np.abs( acw[i  ,j,  k  ] ) < 10**(-8) else dx**2*dy**2/acw[i  ,j,  k  ] 
                
                ac[c] = axp[c]+axm[c]+ayp[c]+aym[c]+azp[c]+azm[c]
                
                b[c] = - (u[i+1,j,k] - u[i,j,k])*dy*dz \
                       - (v[i,j+1,k] - v[i,j,k])*dx*dz \
                       - (w[i,j,k+1] - w[i,j,k])*dx*dy
    
    def sor(dp0, omega):
        p0 = np.copy(dp0) # deep copy
        for i in range(nx):
            for j in range(ny):
                for k in range(nz):
                    c = (i,j,k)
                    xp,xm,yp,ym,zp,zm = (i+1,j,k),(i-1,j,k),(i,j+1,k),(i,j-1,k),(i,j,k+1),(i,j,k-1)
                    
                    if i == 0:
                        pxp ,pxm = p0[xp], zero
                    elif i == nx-1:
                        pxp ,pxm = zero, p0[xm]
                    else:
                        pxp ,pxm = p0[xp], p0[xm]
                    
                    if j == 0:
                        pyp ,pym = p0[yp], zero
                    elif j == ny-1:
                        pyp ,pym = zero , p0[ym]
                    else:
                        pyp ,pym = p0[yp], p0[ym]
                        
                    if k == 0:
                        pzp, pzm = p0[zp], zero
                    elif k == nz-1:
                        pzp, pzm = zero, p0[zm]
                    else:
                        pzp ,pzm = p0[zp], p0[zm]
                    
                    pt = (axp[c]*pxp + axm[c]*pxm \
                        + ayp[c]*pyp + aym[c]*pym \
                        + azp[c]*pzp + azm[c]*pzm + b[c])/ac[c]
                    
                    p0[c] += omega*(pt - p0[c])
    #                 p0[c] = ut
                    
        return p0
                    
    dpnew = np.copy(dp) # deep copy
    for nn in range(1000):
        dptmp = sor(dpnew, float(1.8))
        err =  np.abs(dpnew - dptmp).max() 
        dpnew = np.copy(dptmp)
        
        if err < 10**(-6) :
            print('SOR N =', nn)
            break
        
    return dpnew
全体のループ
  • 収束計算の係数は標準値を採用
def SIMPLE(p,u,v,w,Uslip,Vslip,alphaU,alphaP,dx,dy,dz,nu):
    
    zero = float(0.0)
    nx, ny, nz = p.shape
    
    for nn in range(1000):
        
        u, acu = calupre(p,u,v,w,u,Uslip,alphaU,dx,dy,dz,nu,axis=0)
        v, acv = calupre(p,u,v,w,v,Vslip,alphaU,dx,dy,dz,nu,axis=1)
        w, acw = calupre(p,u,v,w,w,zero,alphaU,dx,dy,dz ,nu,axis=2)
        dp = caldp(p,u,v,w,acu,acv,acw,dx,dy,dz)
        
        p += alphaP*dp
        
        unew = np.copy(u) # deep copy
        vnew = np.copy(v) # deep copy
        wnew = np.copy(w) # deep copy
        
        for i in range(nx):
            for j in range(ny):
                for k in range(nz):
                    if i !=0 : unew[i,j,k] = u[i,j,k] + (dp[i-1,j,k] - dp[i,j,k])*dy*dz/acu[i,j,k]
                    if j !=0 : vnew[i,j,k] = v[i,j,k] + (dp[i,j-1,k] - dp[i,j,k])*dx*dz/acv[i,j,k]
                    if k !=0 : wnew[i,j,k] = w[i,j,k] + (dp[i,j,k-1] - dp[i,j,k])*dx*dy/acw[i,j,k]
        
        err =  max([np.abs(unew - u).max(), np.abs(vnew - v).max(), np.abs(wnew - w).max()])
        print('N', nn, 'error =', err)
        
        u = np.copy(unew) # deep copy
        v = np.copy(vnew) # deep copy
        w = np.copy(wnew) # deep copy
        
        if err < 10**(-6) : break
        
    return p,u,v,w
コード全体

こちら

結論(タイトルの回収とか)

  • SIMPLE法の解法を理解していることを前提とするとpythonを使うことにより随分スッキリとコーディングできると思います。少なくともFortranやCのグローバル変数まみれのものよりはましだと思います。でも、やっぱりSIMPLE法は煩雑です。
  • 連立一次方程式の解法の部分を関数化すればもっとスッキリする。
  • ちなみに計算時間は私のノートPC(スペックはこちら)で10分程度。フル3Dの計算が、ノートPCで、スクリプト言語で、10分でできるっていい時代になりましたね。

SIMPLE法について思うこと

  • SIMPLE法は、高速に定常解を求めることが可能であり、OpenFOAMに採用される等実務的な解法だと思う。
  • 更に本領を発揮するのはコロケート格子や非構造格子の場合であり、Rhie and Chowの補間を用いることにより複雑な形状でも解析可能となる。 その場合、連立一次方程式の解法に非対称行列に適用可能なものを用いる必要がある。近年は共役勾配法ベースのもの(BiCGStab法とか)が一般的だと思う。 また、プログラミングとしてもオブジェクト指向に適している。
  • 欠点は高次精度のスキームが使えないこと。乱流問題とかには向いてないのかな。もちろん矩形格子なら可能です。
  • 非定常問題はPISO等もあるが基本的には向いてないのかな。でも、非定常問題で厳密な計算を行う機会は少ないのであまり考えなくてもいいかなと。

参考書籍

プログラミングは以下を参考にしました。

レンダリングについて

Holoviews

二次元の図はHoloviewsで書いてます。 ソースはこちら

関連記事

computational-sediment-hyd.hatenablog.jp

ParaView

三次元の図はParaViewで書きました。 XMLベースのvtk形式に関する覚書 - Qiita を参考にxmlベースのvtkを作成しました。 ソースはこちら

ParaView Glance

話題のParaView Glance(web版ParaView)でレンダリングしてみました。

以下のサイトを参考にParaViewで作った図をvtkjsでエクスポートしました。以下ではhtmlも作成していますが上手くいきませんでした。

描画はParaView Glanceにアクセスしてvtkjsファイル(サンプル)をロードするだけです。 手順をgif化しました。

f:id:SedimentHydraulics:20200610144316g:plain

GitHub

github.com

おわりに

  • 次はコロケート版を書こうかなと思ってます。
  • その前に高速化の話を書くかも。
  • 詳細な解説などリクエストがあればコメントをお願いします。