LoginSignup
9
8

More than 3 years have passed since last update.

@ishigaki

HSMAC法によるキャビティーフロー解析のJuliaコード

Last updated at Posted at 2020-10-28

はじめに

HSMAC法（SOLA法）でキャビティーフローを解くためのコードをJuliaで作成しました。

岡本 正芳「数値計算による流体力学- ポテンシャル流,層流,そして乱流へ -」コロナ社（2016）
にあるFortranによるHSMAC法のコードを参考にしました。

補足

計算のリスタートができるようにしたコードをGithubにおきました。
https://github.com/ishigakim/2Dcavity.git

コード概要

  • 2次元のキャビティーフローの解析コード。
  • 圧力解法にはHSMAC法。
  • スタッガード格子を用いており，セル中心に圧力，uはセルの右側のフェイス，vはセルの上側のフェイスで定義。
  • スタッガード格子でu, vの位置がずれているので，境界条件には中点で速度0になるような条件が出てくる。
  • デフォルト設定では40x40の格子で，Reynolds数=100を計算。正方形領域で上側の境界(y=1)でu=1とする。それ以外の境界ではすべりなし。

物理量の配置の模式図。オレンジのところが速度0の滑りなし条件の箇所。薄い赤のところでuもしくはvが0になるようにu_{i,0}等の条件を設定する。ちなみにJuliaのコードではインデックスが1から始まる都合で，図と添字がずれます。
Picture1.png

Juliaコード

ソルバ

using Parameters
using PyPlot
using BenchmarkTools
using HDF5

@with_kw struct Param{T1,T2}
    dt::T1 = 0.01
    nu::T1 = 1e-2
    U::T1 = 1.0
    nstep::T2 = 200
    omega::T1 = 1.7
    # max iteration no. of inner iteration
    niter::T2 = 50000

    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 = "uv" 

end

function main(para)
    @unpack nx, ny, nstep, ndata = para
    # velocity
    u = zeros(Float64, nx+2, ny+2)
    v = zeros(Float64, nx+2, ny+2)
    # auxiliary velocity
    ua = zeros(Float64, nx+2, ny+2)
    va = zeros(Float64, nx+2, ny+2)
    # pressure, corrected pressure
    pa = zeros(Float64, nx+2, ny+2)
    pc = zeros(Float64, nx+2, ny+2)
    # coordinate
    x = zeros(Float64, nx+1, ny+1)
    y = zeros(Float64, nx+1, ny+1)

    uo = zeros(nx+1, ny+1)
    vo = zeros(nx+1, ny+1)
    uxy = zeros(nx+1, ny+1)
    vxy = zeros(nx+1, ny+1)
    div = zeros(Float64, nx+2, ny+2)

    init(para, x, y)
    BC(para, u, v)
    for n = 1:nstep
        solveV(para, u, v, ua, va, pa, uo, vo, uxy, vxy)
        SOLA(para, n, u, v, ua, va, pa, pc, div)
        if n%ndata == 0
            #output(para, n, u, v, pa, x, y)
            #println("u[20,20]=",u[21,21])
        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, u, v)
    @unpack nx, ny, U = para

    for i=1:nx+2
        u[i, ny+2] = 2.0U - u[i, ny+1]
        v[i, ny+1] = 0.0
        u[i, 1] = -u[i, 2]
        v[i, 1] = 0.0
    end
    for j=1:ny+2
        u[nx+1, j] = 0.0
        v[nx+2, j] = -v[nx+1, j]
        u[1, j] = 0.0
        v[1,j] = -v[2, j]
    end
    #println("BC",u[nx-1:nx+2,ny-1:ny+2])
end

function solveV(para, u, v, ua, va, pa, uo, vo, uxy, vxy)
    @unpack dt, nu, dx, dy, nx, ny, U = para

    for j=2:ny+1
        for i=2:nx+1
            uo[i,j] = (u[i,j] + u[i-1,j])*0.5
            vo[i,j] = (v[i,j] + v[i,j-1])*0.5
            uxy[i,j] = (u[i,j+1] + u[i,j])*0.5
            vxy[i,j] = (v[i+1,j] + v[i,j])*0.5
        end
    end    

    for j=2:ny+1
        for i=2:nx
            ua[i,j] = u[i,j] + dt*(
                -(uo[i+1,j]^2 - uo[i,j]^2)/dx
                -(uxy[i,j]*vxy[i,j] - uxy[i,j-1]*vxy[i,j-1])/dy
                + nu*(u[i+1,j] -2.0u[i,j] + u[i-1,j])/dx^2
                + nu*(u[i,j+1] -2.0u[i,j] + u[i,j-1])/dy^2
                - (pa[i+1,j] - pa[i,j])/dx
                )
        end
    end

    for j=2:ny
        for i=2:nx+1
            va[i,j] = v[i,j] + dt*(
                -(uxy[i,j]*vxy[i,j] - uxy[i-1,j]*vxy[i-1,j])/dx
                -(vo[i,j+1]^2 - vo[i,j]^2)/dy
                + nu*(v[i+1,j] -2.0v[i,j] + v[i-1,j])/dx^2
                + nu*(v[i,j+1] -2.0v[i,j] + v[i,j-1])/dy^2
                - (pa[i,j+1] - pa[i,j])/dy
                )
        end
    end
    BC(para, ua, va)

    #println("solveV", ua[nx-1:nx+2,ny-1:ny+2])
end

function SOLA(para, n, u, v, ua, va, pa, pc, div)
    @unpack dt, nu, dx, dy, nx, ny, niter, omega = para

    for m=1:niter
        divmax::Float64 = 0.0
        for j=2:ny+1
            for i=2:nx+1
                div[i,j] = (ua[i,j] - ua[i-1,j])/dx + (va[i,j] - va[i,j-1])/dy
                divmax = max(abs(div[i,j]), divmax)
                pc[i,j] = - omega*div[i,j]/dt/(1.0/dx^2 + 1.0/dy^2)*0.5
                ua[i,j] = ua[i,j] + dt/dx*pc[i,j]
                ua[i-1,j] = ua[i-1,j] - dt/dx*pc[i,j]
                va[i,j] = va[i,j] + dt/dy*pc[i,j]
                va[i,j-1] = va[i,j-1] - dt/dy*pc[i,j]
                pa[i,j] = pa[i,j] + pc[i,j]
            end
        end    
        if divmax <= 1e-8
            if n%100 == 0
                println("divmax=", divmax)
                println("iter, max/min div=", m, ", ", findmax(abs.(view(div,2:nx+1, 2:ny+1))), ", ", findmin(abs.(view(div,2:nx+1, 2:ny+1))))
            end
            break
        end
    end

    for j=2:ny+1
        for i=2:nx
            u[i,j] = ua[i,j]
        end
    end
    for j=2:ny
        for i=2:nx+1
            v[i,j] = va[i,j]
        end
    end
    BC(para, u, v)

    #println("sola",u[nx-1:nx+2,ny-1:ny+2])
end

function output(para, n, u, v, pa, x, y)
    @unpack outfile, dt, nx, ny = para
    time = string(dt*n)
    println("time=", time)
    filename = outfile*time*".h5"
    h5open(filename, "w") do file
        write(file, "u", (u[2:nx+1, 2:ny+1]+u[1:nx, 2:ny+1])/2)
        write(file, "v", (v[2:nx+1, 2:ny+1]+v[2:nx+1, 1:ny])/2)
        write(file, "pa", pa[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])
    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
  • Paramで計算のパラメータを設定
  • main中に計算で使う配列を定義。
  • initが初期条件設定。BCが境界条件設定。
  • solveVが予測速度の計算。SOLAでdivV=0になるように反復計算を実行。速度を修正。
  • output，HDF5でデータ保存。CSVで保存したいときはコメントアウト部分を復活させる。

計算

para = Param(nstep=2000, ndata=100, )
@btime main(para)
10.418 s (13370 allocations: 1.13 MiB)
  • Paramでパラメータを与える。
  • ちなみに参考にしたFotranではgfotranつかって，　実行に25秒くらい。 実行に15秒くらい。Juliaのほうが早い。（Juliaのほうが早いのは早いですが，Fotranのほうの条件間違ってました。Intel FortranだとFortranのほうが少し早いくらいです）
  • 環境2.3 GHz Quad-Core Intel Core i5, メモリ16 GB 2133 MHz LPDDR3，MacBook Pro (13-inch, 2018, Four Thunderbolt 3 Ports)

実行時のメモリアロケーションが13370 allocationsとなってるが，途中にあるprintが原因。途中にprintによる結果出力をなくせば，最初の配列の宣言の分だけになる（13 allocationになる）。

ポスト処理

filename = "uv10.0.h5"
file = h5open(filename, "r") 
u = read(file, "u")
v = read(file, "v")
x = read(file, "x")
y = read(file, "y")
pa = read(file, "pa")
close(file)

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=15 )
fig.colorbar(heatmap, ax=ax1)
ax2 = fig.add_subplot(222)
heatmap2 = ax2.contour(X', Y, v', cmap="jet", levels=15)
fig.colorbar(heatmap2, ax=ax2)

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[:, :], y[:, :], u[:,:], v[:, :], umag, scale=2, cmap="jet")
#image = ax3.streamplot(x', y', u, v)
cb=fig.colorbar(image, ax=ax3)

ax4 = fig.add_subplot(224)
heatmap3 = ax4.contour(X', Y, pa', cmap="jet", levels=15)
fig.colorbar(heatmap3, ax=ax4)

fig = plt.figure(figsize = (5,5))
ax1 = fig.add_subplot(111)
image = ax1.streamplot(x'[:, :], y'[:, :], u'[:,:], v'[:, :], color=u, cmap="jet")
  • 速度コンター，速度ベクトル，圧力コンター，流線を書くコード。

uv-vec.png
stream.png

最後

見本のコードがあって，それをうつしただけなのに，まともに動かすのに1日かかりました。途中で変数を間違えていて，計算がうまくいかなくなっているのを，発見するのにだいぶてこずりました。

適宜境界条件を変えれば，他の計算もできます。

付録

ないよりはあったほうがましかな程度の式のメモ書き。

Picture2.png
Picture3.png
Picture4.png

9
8
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
9
8