はじめに

ブジネスク近似で浮力を考慮した2次元熱流体解析用のコードです。
底面を高温，上面を低温にしたベナール対流の解析用です。
無次元化した方程式をといてます。
「数値流体解析の基礎 - Visual C++とgnuplotによる圧縮性・非圧縮性流体解析」「数値計算による流体力学- ポテンシャル流,層流,そして乱流へ」を参考にしました。

基礎方程式

コードの概要

Paramで計算につかうパラメータを設定してます。
計算のログをlog.txtに書き出すようにしてます。
mainで変数の宣言をして，初期化init, 境界条件BCを呼び出す。solveVで予測速度を計算。SOLAで圧力の修正します。solveTで温度を計算します。
streamfunctionは流線を計算してます。
計算リスタートするときは2個目のmainを呼びます。

コード

using LinearAlgebra
using Parameters
using PyPlot
using BenchmarkTools
using HDF5

@with_kw struct Param{T1,T2}
    dt::T1 = 0.01
    # Grashoff number
    Ra::T1 = 1e4
    Pr::T1 = 1.0
    Gr::T1 = Ra/Pr
    # boundary temp
    Th::T1 = 1.0
    Tc::T1 = 0.0
    
    nstep::T2 = 200
    omega::T1 = 1.7
    # max iteration no. of inner iteration
    niter::T2 = 2000
    
    nx::T2 = 100
    ny::T2 = 20
    
    # length
    Lx::T1 = 5.0
    Ly::T1 = 1.0
    dx::T1 = Lx/nx
    dy::T1 = Ly/ny
   
    ndata::T2 = 100
     
    outfile::String = "Data/uvT" 
    logfile::String = "Data/log.txt"

end

function main(para)
    @unpack nx, ny, nstep, ndata, logfile, dt, Ra, Pr, Gr, Th, Tc, Lx, Ly = para
    log = open(logfile,"a")
    println(log, "==================")
    println(log, "parameters")
    println(log, "Ra=", Ra)
    println(log, "Pr=", Pr)
    println(log, "Gr=", Gr)
    println(log, "Th=", Th)
    println(log, "Tc=", Tc)
    println(log, "Lx=", Lx)
    println(log, "Ly=", Ly)
    println(log, "nx=", nx)
    println(log, "ny=", ny)
    println(log, "dt=", dt)
    println(log, "==================")
    close(log)
    
    # 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)
    # temperature
    T = zeros(Float64, nx+2, ny+2)
    Ta = 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)
    
    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, T)
    
    for n = 1:nstep
        solveV(para, u, v, ua, va, pa, uo, vo, uxy, vxy, T)
        SOLA(para, n, u, v, ua, va, pa, pc, T, div)
        solveT(para, u, v, T, Ta)
        if n%ndata == 0
            log = open(logfile,"a")
            println(log, n, " th step")
            close(log)
            streamfunction(para, u, v, psi)
            output(para, n, 0, u, v, pa, T, psi, x, y)
            #println("u[20,20]=",u[21,21])
        end
    end
end

function main(para, restartFile, starttime)
    @unpack nx, ny, nstep, ndata, logfile,  dt, Ra, Pr, Gr, Th, Tc, Lx, Ly = para
    
    log = open(logfile,"a")
    println(log, "==================")
    println(log, "parameters")
    println(log, "Ra=", Ra)
    println(log, "Pr=", Pr)
    println(log, "Gr=", Gr)
    println(log, "Th=", Th)
    println(log, "Tc=", Tc)
    println(log, "Lx=", Lx)
    println(log, "Ly=", Ly)
    println(log, "nx=", nx)
    println(log, "ny=", ny)
    println(log, "dt=", dt)
    println(log, "==================")
    close(log)
    
    # 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)
    # temperature
    T = zeros(Float64, nx+2, ny+2)
    Ta = 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)
    
    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)
    
    nstart = 0
    
    nstart = restart(para, restartFile, x, y, u, v, pa, T, nstart)
    println("nstart=",nstart)
    println(size(u), size(v), size(pa))
    #println(u[5,:], v[5,:], pa[5,:])
    
    BC(para, u, v, T)
    for n = 1:nstep
        solveV(para, u, v, ua, va, pa, uo, vo, uxy, vxy, T)
        SOLA(para, n, u, v, ua, va, pa, pc, T, div)
        solveT(para, u, v, T, Ta)
        if n%ndata == 0
            streamfunction(para, u, v, psi)
            output(para, n, starttime, u, v, pa, T, psi, x, y)
            #println("u[20,20]=",u[21,21])
            log = open(logfile,"a")
            println(log, n, " th step")
            close(log)
        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 restart(para, restartFile, x, y, u, v, pa, T, nstart)
    @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
    
    file = h5open(restartFile, "r") 
    u .= read(file, "u_latest")
    v .= read(file, "v_latest")
    pa .= read(file, "pa_latest")
    T .= read(file, "T_latest")
    nstart = read(file, "n")
    
    return nstart
end

function BC(para, u, v, T)
    @unpack nx, ny, Th, Tc = para
    
    for i=1:nx+2
        u[i, ny+2] = -u[i, ny+1]
        v[i, ny+1] = 0.0
        T[i, ny+2] = 2Tc - T[i, ny+1]
        #T[i, ny+2] = T[i, ny+1]
        u[i, 1] = -u[i, 2]
        v[i, 1] = 0.0
        T[i, 1] = 2Th - T[i, 2]
        #T[i, 1] = T[i, 2] 
    end
    for j=1:ny+2
        u[nx+1, j] = 0.0
        v[nx+2, j] = -v[nx+1, j]
        T[nx+2, j] = T[nx+1, j]
        #T[nx+2, j] = 2Tc - T[nx+1, j]
        u[1, j] = 0.0
        v[1, j] = -v[2, j]
        T[1, j] = T[2, j]
        #T[1, j] = 2Th - T[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, T)
    @unpack dt, Gr, dx, dy, nx, ny = para

    nu = 1/sqrt(Gr)
    
    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
                + (T[i,j] + T[i,j+1])*0.5
                )
        end
    end
    BC(para, ua, va, T)
    
    #println("solveV", ua[nx-1:nx+2,ny-1:ny+2])
end

function SOLA(para, n, u, v, ua, va, pa, pc, T, div)
    @unpack dt, dx, dy, nx, ny, niter, omega, logfile = 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-10
            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))))
                log = open(logfile, "a")
                println(log, "iter=", m, ", n=", n)
                close(log)
                
            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, T)
    
    #println("sola",u[nx-1:nx+2,ny-1:ny+2])
end

function solveT(para, u, v, T, Ta)
    @unpack dt, Gr, Pr, dx, dy, nx, ny = para
    
    alpha = 1/Pr/sqrt(Gr)
    for j=2:ny+1
        for i=2:nx+1
            Ta[i,j] = T[i,j] + dt*(
                # central differecing
                # -(u[i,j]*(T[i,j] + T[i+1,j])*0.5 - u[i-1,j]*(T[i-1,j] + T[i,j])*0.5)/dx 
                # -(v[i,j]*(T[i,j] + T[i,j+1])*0.5 - v[i,j-1]*(T[i,j-1] + T[i,j])*0.5)/dy
                # upwind differencing
                + (u[i-1,j]*(T[i-1,j] + T[i,j])*0.5 - abs(u[i-1,j])*(T[i,j] - T[i-1,j])*0.5)/dx
                - (u[i,j]*(T[i,j] + T[i+1,j])*0.5 - abs(u[i,j])*(T[i+1,j] - T[i,j])*0.5)/dx
                + (v[i,j-1]*(T[i,j-1] + T[i,j])*0.5 - abs(v[i,j-1])*(T[i,j] - T[i,j-1])*0.5)/dy
                - (v[i,j]*(T[i,j] + T[i,j+1])*0.5 - abs(v[i,j])*(T[i,j+1] - T[i,j])*0.5)/dy
                + alpha*(T[i+1,j] -2.0T[i,j] + T[i-1,j])/dx^2
                + alpha*(T[i,j+1] -2.0T[i,j] + T[i,j-1])/dy^2
                
            )
        end
    end
    for j=2:ny+1
        for i=2:nx+1
            T[i,j] = Ta[i,j]
        end
    end
    BC(para, u, v, T)
end

function output(para, n, starttime, u, v, pa, T, psi, x, y)
    @unpack outfile, dt, nx, ny = para
    time = string(dt*n + starttime) 
    #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, "P", pa[2:nx+1, 2:ny+1])
        write(file, "T", T[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, "u_latest", u)
        write(file, "v_latest", v)
        write(file, "pa_latest", pa)
        write(file, "T_latest", T)
        write(file, "n", n)
    end
    
end

function streamfunction(para, u, v, psi)
    @unpack dx, dy, nx, ny = para
    for j=2:ny+1
        for i=2:nx+1
            psi[i,j] = 0.0
        end
    end    
    
    for j=2:ny+1
        for i=2:nx+1
            for k=2:j
                psi[i,j] += u[i,k]*dy
            end
        end
    end        
end

実行

para = Param(nstep=8000, outfile="DataRa1e4/uvT", 
    logfile="DataRa1e4/log.txt")
@time main(para)

リスタートのとき

リスタートに使うデータと時刻を引数で指定します。

@time main(para, "DataRa1e4/uvT30.0.h5", 30.0)

計算

レイリー数Ra=1e7，nx=250, ny=50, dt=2.5e-3程度にすると，初期条件が少し違っただけで異なる挙動を示すのでおもしろいです。具体的には，できる渦の数が変わってきます。

計算結果の例。上から水平方向速度u, 垂直方向速度vのコンター，流速ベクトル，温度コンター。