コロケート格子上でのSIMPLE法によるキャビティーフロー解析コード（Julia）

はじめに

「数値流体解析の基礎 - Visual C++とgnuplotによる圧縮性・非圧縮性流体解析」を参考にしました。

HSMAC法で，スタッガード格子のコードは下記にあります。

ここでは，コロケート格子使って，SIMPLE法で解析するコードを示します。コロケート格子というのはx方向の速度u，y方向の速度v，圧力pがともにセルの中心で定義された格子です。スタッガード格子はu, v, pの定義位置がずれた格子です。HSMAC法のコードの説明に図があるので参考にしてください。

コロケート格子を使って，なんの対処もしないと，チェッカーボード問題（圧力が格子点ごとに振動する問題）が発生します。それを避けるために，Rhie-Chow補間という手法があります。セル界面の速度に圧力勾配の補正をいれることで，圧力が振動するのを抑制する方法です。今回それを行ったコードになってます。コードではue, vnを計算しているところです。

コロケート格子のため境界での圧力の条件が必要になります。ここでは圧力勾配0としてます。dp/dn=0 (nは壁垂直方向）。壁隣接セルでの圧力をその内側と同じ値にすることで境界条件としています。

コロケート格子なので，格子が変形してもわかりやすいというのがあります。このコードの改良として，一般座標系に対応させる等があります。一般座標系に対応できれば，円筒周りの流れの計算なんかもできるようになります（多分）。

コード

using Parameters

using HDF5

@with_kw struct Param{T1,T2}
    dt::T1 = 0.01
    # viscosity
    nu::T1 = 1e-2
    # boundary velocity
    U::T1 = 1.0
    nstep::T2 = 2000
    alphaua::T1 = 0.5
    alphava::T1 = 0.5
    alphapa::T1 = 0.5
    niter::T2 = 5000
    ipter::T2 = 5
    
    # mesh size
    nx::T2 = 40
    ny::T2 = 40
    
    # length
    Lx::T1 = 1.0
    Ly::T1 = 1.0
    dx::T1 = Lx/nx
    dy::T1 = Ly/ny

    ndata::T2 = 50

    outfile::String = "Data-simple/uv" 

end

function main(para)
    @unpack nx, ny, nstep, niter, ndata = para
    # auxiliary velocity
    ua = zeros(Float64, nx+2, ny+2)
    utemp = zeros(Float64, nx+2, ny+2)
    uold = zeros(Float64, nx+2, ny+2)
    apuva = zeros(Float64, nx+2, ny+2)
    aeuva = zeros(Float64, nx+2, ny+2)
    awuva  = zeros(Float64, nx+2, ny+2)
    anuva = zeros(Float64, nx+2, ny+2)
    asuva = zeros(Float64, nx+2, ny+2)
    bua  = zeros(Float64, nx+2, ny+2)
    va = zeros(Float64, nx+2, ny+2)
    vtemp = zeros(Float64, nx+2, ny+2)
    vold = zeros(Float64, nx+2, ny+2)
    bva = zeros(Float64, nx+2, ny+2)
    ue = zeros(Float64, nx+2, ny+2)
    vn = zeros(Float64, nx+2, ny+2)
    
    # pressure, corrected pressure
    pa = zeros(Float64, nx+2, ny+2)
    pc = zeros(Float64, nx+2, ny+2)
    ptemp = zeros(Float64, nx+2, ny+2)
    appc = zeros(Float64, nx+2, ny+2)
    aepc = zeros(Float64, nx+2, ny+2)
    awpc = zeros(Float64, nx+2, ny+2)
    anpc = zeros(Float64, nx+2, ny+2)
    aspc = zeros(Float64, nx+2, ny+2)
    bpc = zeros(Float64, nx+2, ny+2)
    # coordinate
    x = zeros(Float64, nx+1, ny+1)
    y = zeros(Float64, nx+1, ny+1)
    # stream function
    psi = zeros(Float64, nx+2, ny+2)
    
    rest = 0.0
    init(para, x, y)

    for l = 1:nstep
        BC(para, ua, va, pc, pa)
        
        for n = 1:niter
            BC(para, ua, va, pc, pa)
            solveUA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, ue)
            solveVA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, vn)
            #println("ua, va=", ua[2,2],", ", va[2,2])
            solvePC(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, ptemp, ue, vn)
            update(para, ua, va, pc, pa, apuva, utemp, vtemp)
            rest = converg(para, ua, va, utemp, vtemp)
            #println("iter, redisual ", n,", ",rest)
            if rest <=1e-7
                if l%10 == 0
                    println(l, " th step, ", n, " th iteration, residual=", rest)
                end
                break
            end
        end
        for j=2:ny+1
            for i=2:nx+1
                uold[i,j] = ua[i,j]
            end
        end
    
        for j=2:ny+1
            for i=2:nx+1
                vold[i,j] = va[i,j]
            end
        end
        if l%ndata == 0
            println("residual=", rest)
            output(para, ua, va, pa, x, y, l)
        end
    end
    
end

function init(para, x, y)
    @unpack nx, ny, dx, dy = para
    
    for j in 2:ny+1
        for i in 2:nx+1
            x[i,j] = (i-2+0.5)*dx
        end
    end
    for j in 2:ny+1
        for i in 2:nx+1
            y[i,j] = (j-2+0.5)*dy
        end
    end
end

function BC(para, ua, va, pc, pa)
    @unpack nx, ny, U = para
    for i=1:nx+2
        ua[i, ny+2] = U
        va[i, ny+2] = 0.0
        ua[i, 1] = 0.0
        va[i, 1] = 0.0
        pc[i, 1] = pc[i, 2]
        pc[i, ny+2] = pc[i, ny+1]
        pa[i, 1] = pa[i, 2]
        pa[i, ny+2] = pa[i, ny+1]
    end
    
    for j=1:ny+2
        ua[nx+2, j] = 0.0
        va[nx+2, j] = 0.0
        ua[1, j] = 0.0
        va[1, j] = 0.0
        pc[nx+2, j] = pc[nx+1, j]
        pc[1, j] = pc[2, j]
        pa[nx+2, j] = pa[nx+1, j]
        pa[1, j] = pa[2, j]
    end
end

function solveUA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, ue)
    @unpack nx, ny, dx, dy, dt, nu, alphaua, alphava, alphapa, niter, ipter = para
    
    for j=2:ny+1
        for i=2:nx+1
            if i==2
                awuva[i,j] = 2nu*dy/dx
            else
                awuva[i,j] = nu*dy/dx + 0.25*dy*(ua[i, j] + ua[i-1,j])

            end
            if i==nx+1
                aeuva[i,j] = 2nu*dy/dx
            else
                aeuva[i,j] = nu*dy/dx - 0.25*dy*(ua[i+1, j] + ua[i,j])
            end
            
            if j==2
                asuva[i,j] = 2nu*dx/dy
            else
                asuva[i,j] = nu*dx/dy + 0.25dx*(va[i,j] + va[i,j-1])
            end
            if j == ny+1
                anuva[i,j] = 2nu*dx/dy
            else
                anuva[i,j] = nu*dx/dy - 0.25dx*(va[i,j+1] + va[i,j])
            end
            
            apuva[i,j] = aeuva[i,j] + awuva[i,j] + anuva[i,j] + asuva[i,j] + dx*dy/dt
            bua[i,j] = -0.5dy*(pa[i+1,j] - pa[i-1,j]) +dx*dy/dt*uold[i,j]
            
            utemp[i,j] = ua[i,j]*(1.0-alphaua) + alphaua*(
                aeuva[i,j]*ua[i+1,j] + awuva[i,j]*ua[i-1,j] + anuva[i,j]*ua[i,j+1] + asuva[i,j]*ua[i,j-1] + bua[i,j])/apuva[i,j]
        end
    end
    for j=2:ny+1
        for i=2:nx+1
            ua[i,j] = utemp[i,j]
        end
    end
    # ue: interface velocity
    for j=2:ny+1
        for i=2:nx
            ue[i,j] = 0.5(ua[i+1,j]+ua[i,j]) + 0.5(dy/apuva[i,j] + dy/apuva[i+1,j])*(pa[i,j]-pa[i+1,j]) -0.25dy/apuva[i,j]*(pa[i-1,j]-pa[i+1,j]) -0.25dy/apuva[i+1,j]*(pa[i,j]-pa[i+2,j])   
        end
    end
    for j=2:ny+1
        ue[1,j] = 0.0
        ue[nx+1,j] = 0.0
    end
end

function solveVA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, vn)
    @unpack nx, ny, dx, dy, dt, nu, alphaua, alphava, alphapa, niter, ipter = para
    
    for j=2:ny+1
        for i=2:nx+1
            bva[i,j] = -0.5dx*(pa[i,j+1] - pa[i,j-1]) + dx*dy/dt*vold[i,j]
            
            vtemp[i,j] = va[i,j] * (1.0 - alphava) + alphava * (
                aeuva[i,j]*va[i+1,j] + awuva[i,j]*va[i-1,j] + anuva[i,j]*va[i,j+1] + asuva[i,j]*va[i,j-1] + bva[i,j])/apuva[i,j]
        end
    end
    for j=2:ny+1
        for i=2:nx+1
            va[i,j] = vtemp[i,j]
        end
    end
    # vn: interface velocity
    for j=2:ny
        for i=2:nx+1
            vn[i,j] = 0.5(va[i,j+1]+va[i,j]) + 0.5(dx/apuva[i,j]+dx/apuva[i,j+1])*(pa[i,j]-pa[i,j+1]) -0.25dx/apuva[i,j]*(pa[i,j-1]-pa[i,j+1]) -0.25dx/apuva[i,j+1]*(pa[i,j]-pa[i,j+2])   
        end
    end
    for i=2:nx+1
        vn[i,1] = 0.0
        vn[i,ny+1] = 0.0
    end
end

function solvePC(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, 
        bva, appc, aepc, awpc, anpc, aspc, bpc, ptemp, ue, vn)
    @unpack nx, ny, dx, dy, nu, alphaua, alphava, alphapa, niter, ipter = para

    for j=2:ny+1
        for i=2:nx+1
            pc[i,j] = 0.0
        end
    end
    
    for j=2:ny+1
        for i=2:nx+1
            if i == 2
                awpc[i,j] = 0.0
            else
                awpc[i,j] = 0.5dy*(dy/apuva[i-1,j]+dy/apuva[i,j])
            end
            
            if i == nx+1
                aepc[i,j] = 0.0
            else          
                aepc[i,j] = 0.5dy*(dy/apuva[i,j]+dy/apuva[i+1,j])
            end
            
            if j == 2
                aspc[i,j] = 0.0
            else
                aspc[i,j] = 0.5dx*(dx/apuva[i,j-1]+dx/apuva[i,j])
            end
                    
            if j == ny+1
                anpc[i,j] = 0.0
            else
                anpc[i,j] = 0.5dx*(dx/apuva[i,j]+dx/apuva[i,j+1])
            end
                        
            appc[i,j] = aepc[i,j] + awpc[i,j] + anpc[i,j] + aspc[i,j]
            bpc[i,j] = -dy*(ue[i,j] - ue[i-1,j]) - dx*(vn[i,j] - vn[i,j-1])
        end
    end
    
    # SOR法で圧力の解を求めてる
    for ip=1:ipter
        for j=2:ny+1
            for i=2:nx+1
                ptemp[i,j] = pc[i,j]
            end
        end
        for j=2:ny+1
            for i=2:nx+1
                pc[i,j] = ptemp[i,j]*(1.0 - alphapa) + alphapa*(
                    aepc[i,j]*ptemp[i+1,j] + awpc[i,j]*ptemp[i-1,j] + anpc[i,j]*ptemp[i,j+1] + aspc[i,j]*ptemp[i,j-1] + bpc[i,j])/appc[i,j]
            end
        end
        for i=1:nx+2
            pc[i, 1] = pc[i, 2]
            pc[i, ny+2] = pc[i, ny+1]
        end
        for j=1:ny+2
            pc[nx+2, j] = pc[nx+1, j]
            pc[1, j] = pc[2, j]
        end
        rest=convergP(para, pc, ptemp)
        if rest < 1e-5 && ip >=2
            println("rest PC ", rest, ", ip=", ip)
            break
        end
    end
end    

function update(para, ua, va, pc, pa, apuva, utemp, vtemp)
    @unpack nx, ny, dx, dy, nu, alphaua, alphava, alphapa, niter, ipter = para
    
    # u = u* + u'
    for j=2:ny+1
        for i=2:nx+1
            utemp[i,j] = ua[i,j]
            ua[i,j] = ua[i,j] -0.5dy*(pc[i+1,j] - pc[i-1,j])/apuva[i,j]
            
        end
    end
    
    for j=2:ny+1
        for i=2:nx+1
            vtemp[i,j] = va[i,j]
            va[i,j] = va[i,j] -0.5dx*(pc[i,j+1] - pc[i,j-1])/apuva[i,j]
            
        end
    end
 
    
    for j=2:ny+1
        for i=2:nx+1
            pa[i,j] = pa[i,j] + alphapa*pc[i,j]
        end
    end
    
    pref = pa[2,2]
    for j=2:ny+1
        for i=2:nx+1
            pa[i,j] -= pref
        end
    end
end

function converg(para, ua, va, utemp, vtemp)
    @unpack nx, ny = para
    
    restua = 0.0
    restva = 0.0
    
    for j=2:ny+1
        for i=2:nx
            restua = restua + abs(ua[i,j] - utemp[i,j])
            
        end
    end
    for j=2:ny
        for i=2:nx+1
            restva = restva + abs(va[i,j] - vtemp[i,j])
        end
    end
    rest = (restua + restva)/2.0
    return rest/(nx*ny)
end

function convergP(para, pc, ptemp)
    @unpack nx, ny = para
    
    rest = 0.0
    
    for j=2:ny+1
        for i=2:nx+1
            rest = rest + abs(pc[i,j] - ptemp[i,j])
        end
    end

    return rest/(nx*ny)
end

function output(para, ua, va, pa, x, y, l)
    @unpack outfile, dt, nx, ny = para
    #println("time=", time)
    filename = outfile*string(l)*".h5"
    h5open(filename, "w") do file
        write(file, "U", ua[2:nx+1, 2:ny+1])
        write(file, "V", va[2:nx+1, 2:ny+1])
        write(file, "P", pa[2:nx+1, 2:ny+1])
        #write(file, "Psi", psi[2:nx+1, 2:ny+1])
        write(file, "X", x[2:nx+1, 2:ny+1])
        write(file, "Y", y[2:nx+1, 2:ny+1])
        write(file, "ua_latest", ua)
        write(file, "va_latest", va)
        write(file, "pa_latest", pa)
    end
    
end

実行方法

さきほどのコードをSIMPLE-colocate.jlとして保存したとします。Simple-Re100ディレクトリにデータを書くようにするには以下のようにします。データがHDF5ファイルで保存されます。

julia> include("SIMPLE-colocate.jl")
julia> para = Param(outfile="Simple-Re100/uv")
julia> @time main(para)