はじめに

コロケート格子でのSIMPLE法コードが下記です。

今回はこれをPISO法に拡張します。
実行方法は上記記事とほぼ同じなので割愛します。

Re=100のキャビティーフローでテストしてて，x方向の速度はSIMPLEの計算とほぼ一致しましたが圧力が振動していました。調べてみたら，Rhie-Chow補間のように運動量に対する補間をすると不足緩和しないといけないようです。
例えば，これとかに書いてあります。

というわけで，不足緩和しながら計算するようにしたら，圧力の振動がなくなりました。それを下記のコードでも反映しました。

コード

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
    niter::T2 = 100
    ipter::T2 = 10
    alphaua::T1 = 0.5
    alphava::T1 = 0.5
    alphapa::T1 = 0.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 = 100

    outfile::String = "Data-piso-co/uv" 

end

function main(para)
    @unpack nx, ny, nstep, niter, ndata = para
    # auxiliary velocity
    ua = zeros(Float64, nx+2, ny+2)
    ua2 = zeros(Float64, nx+2, ny+2)
    ua3 = 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)
    va2 = zeros(Float64, nx+2, ny+2)
    va3 = 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)
    uW = zeros(Float64, nx+2, ny+2)
    vN = zeros(Float64, nx+2, ny+2)
    vS = zeros(Float64, nx+2, ny+2)

    uf = zeros(Float64, nx+2, ny+2)
    vf = 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)
    pc2 = 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, ua2, va2, pa, pc, pc2)
        
        for n = 1:niter
            BC(para, ua, va, ua2, va2, pa, pc, pc2)
            solveUA(para, ua, va, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, uf)
            solveVA(para, ua, va, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, vf)
            solvePC(para, ua, va, pc, ptemp, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, uf, vf)
            update(para, ua, va, pc, apuva, ua2, va2)
            solvePC2(para, ua, ua2, va, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, uE, uW, vN, vS, ptemp)
            update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva)
            
            rest = converg(para, ua, va, utemp, vtemp)
            #println("iter, redisual ", n,", ",rest)
            
            if rest <=1e-7 && n>=2
                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
    #streamfunction(para, u, v, psi)
    
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, ua2, va2, pa, pc, pc2)
    @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
        # ua2[i, ny+2]をupdate2でつかうが，updateではここの値は更新されないので，BCで与える。
        ua2[i, ny+2] = U
        va2[i, ny+2] = 0.0
        ua2[i, 1] = 0.0
        va2[i, 1] = 0.0
        pa[i, 1] = pa[i, 2]
        pa[i, ny+2] = pa[i, ny+1]
        pc[i, 1] = pc[i, 2]
        pc[i, ny+2] = pc[i, ny+1]
        pc2[i, 1] = pc2[i, 2]
        pc2[i, ny+2] = pc2[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
        ua2[nx+2, j] = 0.0
        va2[nx+2, j] = 0.0
        ua2[1, j] = 0.0
        va2[1, j] = 0.0
        pa[1, j] = pa[2, j]
        pa[nx+2, j] = pa[nx+1, j]
        pc[nx+2, j] = pc[nx+1, j]
        pc[1, j] = pc[2, j]
        pc2[nx+2, j] = pc2[nx+1, j]
        pc2[1, j] = pc2[2, j]
    end
end

function solveUA(para, ua, va, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, uf)
    @unpack nx, ny, dx, dy, dt, nu, alphaua = para
    
    for j=1:ny+2
        for i=1:nx+2
            if i==nx+2 
                aeuva[i,j] = 0.0
            elseif 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 i==1
                awuva[i,j] = 0.0
            elseif 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 j==1
                asuva[i,j] = 0.0
            elseif 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+2
                anuva[i,j] = 0.0
            elseif 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

        end
    end

    for j=2:ny+1
        for i=2:nx+1
            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

    for j=2:ny+1
        for i=2:nx
            uf[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
        uf[1,j] = 0.0
        uf[nx+1,j] = 0.0
    end
end

function solveVA(para, ua, va, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, vf)
    @unpack nx, ny, dx, dy, dt, nu, alphava = 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

    for j=2:ny
        for i=2:nx+1
            vf[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
        vf[i,1] = 0.0
        vf[i,ny+1] = 0.0
    end
end

function solvePC(para, ua, va, pc, ptemp, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, uf, vf)
    @unpack nx, ny, dx, dy, nu, ipter, alphapa = para
    
    for j=2:ny+1
        for i=2:nx+1
            pc[i,j] = 0.0
            bpc[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*(uf[i,j] - uf[i-1,j]) - dx*(vf[i,j] - vf[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
    #println("pc=",pc[15,15])
    #println("max pc=",findmax(abs.(pc)))
end    

function update(para, ua, va, pc, apuva, ua2, va2)
    @unpack nx, ny, dx, dy,  = para
    
    # u** = u* + u'
    for j=2:ny+1
        for i=2:nx+1
            ua2[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
            va2[i,j] = va[i,j] -0.5dx*(pc[i,j+1] - pc[i,j-1])/apuva[i,j]
        end
    end
end

function solvePC2(para, ua, ua2, va, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, uE, uW, vN, vS, ptemp)
    @unpack nx, ny, dx, dy, nu, ipter, alphapa = para

    for j=2:ny+1
        for i=2:nx+1
            pc2[i,j] = 0.0
            bpc[i,j] = 0.0
        end
    end

    for j=2:ny+1
        for i=2:nx+1

            if i==nx+1
                uE[i,j] = awuva[i+1,j]*(ua2[i,j]-ua[i,j])
            else
                uE[i,j] = aeuva[i+1,j]*(ua2[i+2,j]-ua[i+2,j]) +awuva[i+1,j]*(ua2[i,j]-ua[i,j]) +anuva[i+1,j]*(ua2[i+1,j+1]-ua[i+1,j+1]) +asuva[i+1,j]*(ua2[i+1,j-1]-ua[i+1,j-1])
            end

            if i==2
                uW[i,j] = aeuva[i-1,j]*(ua2[i,j]-ua[i,j])
            else
                uW[i,j] = aeuva[i-1,j]*(ua2[i,j]-ua[i,j]) +awuva[i-1,j]*(ua2[i-2,j]-ua[i-2,j]) +anuva[i-1,j]*(ua2[i-1,j+1]-ua[i-1,j+1]) +asuva[i-1,j]*(ua2[i-1,j-1]-ua[i-1,j-1])
            end

            if j==ny+1
                vN[i,j] = asuva[i,j+1]*(va2[i,j]-va[i,j])
            else
                vN[i,j] = aeuva[i,j+1]*(va2[i+1,j+1]-va[i+1,j+1]) +awuva[i,j+1]*(va2[i-1,j+1]-va[i-1,j+1]) +anuva[i,j+1]*(va2[i,j+2]-va[i,j+2]) +asuva[i,j+1]*(va2[i,j]-va[i,j])
            end

            if j==2
                vS[i,j] = anuva[i,j-1]*(va2[i,j]-va[i,j]) 
            else
                vS[i,j] = aeuva[i,j-1]*(va2[i+1,j-1]-va[i+1,j-1]) +awuva[i,j-1]*(va2[i-1,j-1]-va[i-1,j-1]) +anuva[i,j-1]*(va2[i,j]-va[i,j]) +asuva[i,j-1]*(va2[i,j-2]-va[i,j-2])
            end
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            bpc[i,j] = -0.5(dy/apuva[i+1,j]*uE[i,j] -dy/apuva[i-1,j]*uW[i,j]) -0.5(dx/apuva[i,j+1]*vN[i,j] -dx/apuva[i,j-1]*vS[i,j])           
        end
    end

    #println("ue, uw, vn, vs ", findmax(abs.(uE)), ",", findmax(abs.(uW)), ", ", findmax(abs.(vN)), ",", findmax(abs.(vS)))
    #println("apuva ", findmax(apuva))
    #println("bpc ",findmax(abs.(bpc)))
    

    # SOR法で圧力の解を求めてる
    for ip=1:ipter
        for j=2:ny+1
            for i=2:nx+1
                ptemp[i,j] = pc2[i,j]
            end
        end
        for j=2:ny+1
            for i=2:nx+1
                pc2[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
            pc2[i, 1] = pc2[i, 2]
            pc2[i, ny+2] = pc2[i, ny+1]
        end
        for j=1:ny+2
            pc2[nx+2, j] = pc2[nx+1, j]
            pc2[1, j] = pc2[2, j]
        end
        rest=convergP(para, pc2, ptemp)
        if rest < 1e-5 && ip>=2
            #println("rest PC2 ", rest, ", ip=", ip)
            break
        end
    end
    #println("max pc2=",findmax(abs.(pc2)))
end

function update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva)
    @unpack nx, ny, dx, dy,  = para
    
    #println("max ua2-ua=",findmax(abs.(ua2[2:nx,2:ny+1]-ua[2:nx,2:ny+1])))
    #println("max va2-va=",findmax(abs.(va2[2:nx+1,2:ny]-va[2:nx+1,2:ny])))
    
    #println("ua2[2,ny+2], ua[2,ny+2]", ua2[2,ny+2], ",",ua[2,ny+2])

    # u** = u* + u'
    for j=2:ny+1
        for i=2:nx+1
            ua3[i,j] = ua2[i,j] + (aeuva[i,j]*(ua2[i+1,j]-ua[i+1,j]) + awuva[i,j]*(ua2[i-1,j]-ua[i-1,j]) + anuva[i,j]*(ua2[i,j+1]-ua[i,j+1]) + asuva[i,j]*(ua2[i,j-1]-ua[i,j-1]) -
            0.5dy*(pc2[i+1,j] - pc2[i-1,j]))/apuva[i,j] 
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            ua[i,j] = ua3[i,j]
        end
    end
    
    for j=2:ny+1
        for i=2:nx+1
            va3[i,j] = va2[i,j] + (aeuva[i,j]*(va2[i+1,j]-va[i+1,j]) + awuva[i,j]*(va2[i-1,j]-va[i-1,j]) + anuva[i,j]*(va2[i,j+1]-va[i,j+1]) + asuva[i,j]*(va2[i,j-1]-va[i,j-1]) -
            0.5dx*(pc2[i,j+1] - pc2[i,j-1]))/apuva[i,j]
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            va[i,j] = va3[i,j]
        end
    end

    #println("max ua3-ua2=",findmax(abs.(ua3[2:nx,2:ny+1]-ua2[2:nx,2:ny+1])))
    #println("max va3-va2=",findmax(abs.(va3[2:nx+1,2:ny]-va2[2:nx+1,2:ny])))
    for j=2:ny+1
        for i=2:nx+1
            pa[i,j] = pa[i,j] + pc[i,j] + pc2[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+1
            restua = restua + abs(ua[i,j] - utemp[i,j])
            
        end
    end
    for j=2:ny+1
        for i=2:nx+1
            restva = restva + abs(va[i,j] - vtemp[i,j])
        end
    end
    #rest = max(restua, restva)
    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
    
    #open( outfile*time*".csv", "w")　 do f
    #    println(f, "x, y, u, v, pa")
    #    for j=2:ny+1
    #        for i=2:nx+1
    #            println(f, x[i,j], ",", y[i,j], ",",(u[i,j]+u[i-1,j])/2., ",", (v[i,j]+v[i,j-1])/2.0, ", ", pa[i,j] )
    #        end
    #    end
    #end
end

お絵かき

お絵かき用のコードはこんな感じです。

name = "uv2000"
filename = "Piso-Re100-co/"*name*".h5"
file = h5open(filename, "r") 
u = read(file, "U")
v = read(file, "V")
x = read(file, "X")
y = read(file, "Y")
pa = read(file, "P")
close(file)

skip=1

X = x[:,1]
Y = y[1,:]
fig = plt.figure(figsize = (15, 13))
ax1 = fig.add_subplot(221)
heatmap = ax1.contour(X', Y, u', cmap="jet", levels=50 )
fig.colorbar(heatmap, ax=ax1)
ax1.set_title("u")

ax2 = fig.add_subplot(222)
heatmap2 = ax2.contour(X', Y, v', cmap="jet", levels=50)
fig.colorbar(heatmap2, ax=ax2)
ax2.set_title("v")

ax3 = fig.add_subplot(223)
umag = sqrt.(u.*u .+ v.*v)
#image = ax3.quiver(x[:, 20:40], y[:, 20:40], u[:, 20:40], v[:, 20:40], umag, scale=4, cmap="jet")
image = ax3.quiver(x[1:skip:size(x)[1], 1:skip:size(x)[2]], y[1:skip:size(x)[1], 1:skip:size(x)[2]], u[1:skip:size(x)[1], 1:skip:size(x)[2]], v[1:skip:size(x)[1], 1:skip:size(x)[2]], umag[1:skip:size(x)[1], 1:skip:size(x)[2]], scale=2, cmap="jet")
#image = ax3.streamplot(x', y', u, v)
cb=fig.colorbar(image, ax=ax3)
ax3.set_title("u,v vector")

ax4 = fig.add_subplot(224)
heatmap3 = ax4.contour(X', Y, pa', cmap="magma", levels=50)
fig.colorbar(heatmap3, ax=ax4)
ax4.set_title("p")

plt.savefig("Piso-Re100-co/Fig/uv-vec"*name*"-co.pdf")

不足緩和の影響

不足緩和なし。左下が圧力のx方向の分布。振動しています。

不足緩和したもの。振動が抑制されてます。