はじめに
無次元化や座標のとり方はこのあたりと同じです。
https://qiita.com/ishigaki/items/866df91ea8f32136fa5a
https://qiita.com/ishigaki/items/60bb3c7fc9a66d278cfa
PISO法の概要についてはVersteeg "An Introduction to Computational Fluid Dynamics" http://ftp.demec.ufpr.br/disciplinas/TM702/Versteeg_Malalasekera_2ed.pdf の教科書を参考にしました。PISO法のベースとなるSIMPLE法については岡本「数値計算による流体力学- ポテンシャル流,層流,そして乱流へ -」を参考にしました。
コードのここおかしくないとかあったらコメントください。
コードの概要
- Paramで計算につかうパラメータを設定してます。
- mainで変数の宣言をして，初期化init, 境界条件BCを呼び出す。
- solveUA, solveVAで1段目の予測速度u*, v*を計算。
- solvePCで1段目の圧力修正量p'の計算をします。
- updateで2段めの予測速度u**, v**を計算します。
- solvePC2で2段めの圧力修正量p''を計算します。
- update2で3段目の予測速度u***, v***を計算します。
- solveTAで温度を計算します。
- mainのfor n = 1:niterのループ中で速度が収束したら，ループを抜けます。次の時間ステップへ。
- 計算リスタートするときは2個目のmainを呼びます。
- 対流項は2次精度の中心差分（線形補間），拡散項は2次精度の中心差分で離散化。時間微分はEuler陰解法で離散化します。
- 下面が高温，上面が低温，側面は断熱。速度については滑りなし。
コード
using Parameters
using HDF5
@with_kw struct Param{T1,T2}
dt::T1 = 0.01
# Rayleigh number, Prandtl number, Grashoff number
Ra::T1 = 1e4
Pr::T1 = 1.0
Gr::T1 = Ra/Pr
nu::T1 = 1.0/sqrt(Gr)
alpha::T1 = 1.0/Pr/sqrt(Gr)
# boundary temp
Th::T1 = 1.0
Tc::T1 = 0.0
nstep::T2 = 200
niter::T2 = 500
ipter::T2 = 100
tol::T1 = 1e-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 = 10
outfile::String = "Data-piso/uv"
end
function main(para)
@unpack nx, ny, nstep, niter, ndata, tol = para
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)
apua = zeros(Float64, nx+2, ny+2)
aeua = zeros(Float64, nx+2, ny+2)
awua = zeros(Float64, nx+2, ny+2)
anua = zeros(Float64, nx+2, ny+2)
asua = 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)
apva = zeros(Float64, nx+2, ny+2)
aeva = zeros(Float64, nx+2, ny+2)
awva = zeros(Float64, nx+2, ny+2)
anva = zeros(Float64, nx+2, ny+2)
asva = zeros(Float64, nx+2, ny+2)
bva = zeros(Float64, nx+2, ny+2)
# temperature
Ta = zeros(Float64, nx+2, ny+2)
Told = zeros(Float64, nx+2, ny+2)
apTa = zeros(Float64, nx+2, ny+2)
aeTa = zeros(Float64, nx+2, ny+2)
awTa = zeros(Float64, nx+2, ny+2)
anTa = zeros(Float64, nx+2, ny+2)
asTa = zeros(Float64, nx+2, ny+2)
bTa = zeros(Float64, nx+2, ny+2)
Ttemp = 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)
# 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)
rest = 0.0
init(para, x, y)
for l = 1:nstep
BC(para, ua, va, Ta, ua2, va2)
for n = 1:niter
solveUA(para, ua, va, pa, apua, aeua, awua, anua, asua, bua, uold, utemp)
solveVA(para, ua, va, Ta, pa, apva, aeva, awva, anva, asva, bva, vold, vtemp)
solvePC(para, ua, va, pc, ptemp, apua, aeua, awua, anua, asua, bua, apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc)
update(para, ua, va, pc, apua, apva, ua2, va2)
solvePC2(para, ua, ua2, va, va2, pc, pc2, apua, aeua, awua, anua, asua, bua, apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc, ue, uw, vn, vs, ptemp)
update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apua, aeua, awua, anua, asua, apva, aeva, awva, anva, asva)
solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
rest = converg(para, ua, va, utemp, vtemp)
#println("iter, redisual ", n,", ",rest)
if rest <=tol && n >=2
if l%10 == 0
println(l, " th step, ", n, " th iteration, residual=", rest)
end
break
end
println(rest)
end
for j=2:ny+1
for i=2:nx
uold[i,j] = ua[i,j]
end
end
for j=2:ny
for i=2:nx+1
vold[i,j] = va[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
Told[i,j] = Ta[i,j]
end
end
if l%ndata == 0
println("residual=", rest)
output(para, ua, va, pa, Ta, uold, vold, Told, x, y, l)
end
end
end
function main(para, restartFile)
@unpack nx, ny, nstep, niter, ndata, tol = 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)
apua = zeros(Float64, nx+2, ny+2)
aeua = zeros(Float64, nx+2, ny+2)
awua = zeros(Float64, nx+2, ny+2)
anua = zeros(Float64, nx+2, ny+2)
asua = 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)
apva = zeros(Float64, nx+2, ny+2)
aeva = zeros(Float64, nx+2, ny+2)
awva = zeros(Float64, nx+2, ny+2)
anva = zeros(Float64, nx+2, ny+2)
asva = zeros(Float64, nx+2, ny+2)
bva = zeros(Float64, nx+2, ny+2)
# temperature
Ta = zeros(Float64, nx+2, ny+2)
Told = zeros(Float64, nx+2, ny+2)
apTa = zeros(Float64, nx+2, ny+2)
aeTa = zeros(Float64, nx+2, ny+2)
awTa = zeros(Float64, nx+2, ny+2)
anTa = zeros(Float64, nx+2, ny+2)
asTa = zeros(Float64, nx+2, ny+2)
bTa = zeros(Float64, nx+2, ny+2)
Ttemp = 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)
# 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)
rest = 0.0
nstart = restart(para, restartFile, x, y, ua, va, pa, Ta, uold, vold, Told)
println("nstart=",nstart)
for l = nstart+1:nstart+nstep
BC(para, ua, va, Ta, ua2, va2)
for n = 1:niter
solveUA(para, ua, va, pa, apua, aeua, awua, anua, asua, bua, uold, utemp)
solveVA(para, ua, va, Ta, pa, apva, aeva, awva, anva, asva, bva, vold, vtemp)
solvePC(para, ua, va, pc, ptemp, apua, aeua, awua, anua, asua, bua, apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc)
update(para, ua, va, pc, apua, apva, ua2, va2)
solvePC2(para, ua, ua2, va, va2, pc, pc2, apua, aeua, awua, anua, asua, bua, apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc, ue, uw, vn, vs, ptemp)
update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apua, aeua, awua, anua, asua, apva, aeva, awva, anva, asva)
solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
rest = converg(para, ua, va, utemp, vtemp)
#println("iter, redisual ", n,", ",rest)
if rest <=tol && n >=2
if l%10 == 0
println(l, " th step, ", n, " th iteration, residual=", rest)
end
break
end
println(rest)
end
for j=2:ny+1
for i=2:nx
uold[i,j] = ua[i,j]
end
end
for j=2:ny
for i=2:nx+1
vold[i,j] = va[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
Told[i,j] = Ta[i,j]
end
end
if l%ndata == 0
println("residual=", rest)
output(para, ua, va, pa, Ta, uold, vold, Told, 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, Ta, ua2, va2)
@unpack nx, ny, Th, Tc = para
for i=1:nx+2
ua[i, ny+2] = 0.0
va[i, ny+1] = 0.0
ua[i, 1] = 0.0
va[i, 1] = 0.0
Ta[i, ny+2] = Tc
Ta[i, 1] = Th
# ua2[i, ny+2]をupdate2でつかうが，updateではここの値は更新されないので，BCで与える。
ua2[i, ny+2] = 0.0
va2[i, ny+1] = 0.0
ua2[i, 1] = 0.0
va2[i, 1] = 0.0
end
for j=1:ny+2
ua[nx+1, j] = 0.0
va[nx+2, j] = 0.0
ua[1, j] = 0.0
va[1, j] = 0.0
Ta[nx+2, j] = Ta[nx+1, j]
Ta[1, j] = Ta[2,j]
ua2[nx+1, j] = 0.0
va2[nx+2, j] = 0.0
ua2[1, j] = 0.0
va2[1, j] = 0.0
end
end
function solveUA(para, ua, va, pa, apua, aeua, awua, anua, asua, bua, uold, utemp)
@unpack nx, ny, dx, dy, dt, nu, = para
for j=2:ny+1
for i=1:nx+1
if i==nx+1
aeua[i,j] = nu*dy/dx
else
aeua[i,j] = nu*dy/dx - 0.25*dy*(ua[i+1, j] + ua[i,j])
end
if i==1
awua[i,j] = nu*dy/dx
else
awua[i,j] = nu*dy/dx + 0.25*dy*(ua[i, j] + ua[i-1,j])
end
if j==2
asua[i,j] = 2nu*dx/dy
else
asua[i,j] = nu*dx/dy + 0.25dx*(va[i+1,j-1] + va[i,j-1])
end
if j == ny+1
anua[i,j] = 2nu*dx/dy
else
anua[i,j] = nu*dx/dy - 0.25dx*(va[i+1,j] + va[i,j])
end
apua[i,j] = aeua[i,j] + awua[i,j] + anua[i,j] + asua[i,j] + dx*dy/dt
bua[i,j] = -dy*(pa[i+1,j] - pa[i,j]) +dx*dy/dt*uold[i,j]
if i!=1 && i!=nx+1
utemp[i,j] = (aeua[i,j]*ua[i+1,j] + awua[i,j]*ua[i-1,j] + anua[i,j]*ua[i,j+1] + asua[i,j]*ua[i,j-1] + bua[i,j])/apua[i,j]
end
end
end
for j=2:ny+1
for i=2:nx
ua[i,j] = utemp[i,j]
end
end
end
function solveVA(para, ua, va, Ta, pa, apva, aeva, awva, anva, asva, bva, vold, vtemp)
@unpack nx, ny, dx, dy, dt, nu, = para
for j=1:ny+1
for i=2:nx+1
if i==2
awva[i,j] = 2nu*dy/dx
else
awva[i,j] = nu*dy/dx + 0.25dy*(ua[i-1,j+1] + ua[i-1,j])
end
if i == nx+1
aeva[i,j] = 2nu*dy/dx
else
aeva[i,j] = nu*dy/dx - 0.25dy*(ua[i,j+1] + ua[i,j])
end
if j==ny+1
anva[i,j] = nu*dx/dy
else
anva[i,j] = nu*dx/dy - 0.25*dx*(va[i, j+1] + va[i,j])
end
if j==1
asva[i,j] = nu*dx/dy
else
asva[i,j] = nu*dx/dy + 0.25*dx*(va[i, j] + va[i,j-1])
end
apva[i,j] = aeva[i,j] + awva[i,j] + anva[i,j] + asva[i,j] + dx*dy/dt
bva[i,j] = -dx*(pa[i,j+1] - pa[i,j]) + dx*dy/dt*vold[i,j] + dx*dy*0.5*(Ta[i, j]+Ta[i, j+1])
if j!=1 && j!=ny+1
vtemp[i,j] = (aeva[i,j]*va[i+1,j] + awva[i,j]*va[i-1,j] + anva[i,j]*va[i,j+1] + asva[i,j]*va[i,j-1] + bva[i,j])/apva[i,j]
end
end
end
for j=2:ny
for i=2:nx+1
va[i,j] = vtemp[i,j]
end
end
end
function solvePC(para, ua, va, pc, ptemp, apua, aeua, awua, anua, asua, bua,
apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc)
@unpack nx, ny, dx, dy, nu, ipter = 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] = dy*dy/apua[i-1,j]
end
if i == nx+1
aepc[i,j] = 0.0
else
aepc[i,j] = dy*dy/apua[i,j]
end
if j == 2
aspc[i,j] = 0.0
else
aspc[i,j] = dx*dx/apva[i,j-1]
end
if j == ny+1
anpc[i,j] = 0.0
else
anpc[i,j] = dx*dx/apva[i,j]
end
appc[i,j] = aepc[i,j] + awpc[i,j] + anpc[i,j] + aspc[i,j]
bpc[i,j] = -dy*(ua[i,j] - ua[i-1,j]) - dx*(va[i,j] - va[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] = (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
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, apua, apva, ua2, va2)
@unpack nx, ny, dx, dy, = para
# u** = u* + u'
for j=2:ny+1
for i=2:nx
ua2[i,j] = ua[i,j] -dy*(pc[i+1,j] - pc[i,j])/apua[i,j]
end
end
for j=2:ny
for i=2:nx+1
va2[i,j] = va[i,j] -dx*(pc[i,j+1] - pc[i,j])/apva[i,j]
end
end
end
function solvePC2(para, ua, ua2, va, va2, pc, pc2, apua, aeua, awua, anua, asua, bua,
apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc, ue, uw, vn, vs, ptemp)
@unpack nx, ny, dx, dy, nu, ipter = 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] = awua[i,j]*(ua2[i-1,j]-ua[i-1,j])
else
ue[i,j] = aeua[i,j]*(ua2[i+1,j]-ua[i+1,j]) +awua[i,j]*(ua2[i-1,j]-ua[i-1,j]) +anua[i,j]*(ua2[i,j+1]-ua[i,j+1]) +asua[i,j]*(ua2[i,j-1]-ua[i,j-1])
end
if i==2
uw[i,j] = aeua[i-1,j]*(ua2[i,j]-ua[i,j])
else
uw[i,j] = aeua[i-1,j]*(ua2[i,j]-ua[i,j]) +awua[i-1,j]*(ua2[i-2,j]-ua[i-2,j]) +anua[i-1,j]*(ua2[i-1,j+1]-ua[i-1,j+1]) +asua[i-1,j]*(ua2[i-1,j-1]-ua[i-1,j-1])
end
if j==ny+1
vn[i,j] = asva[i,j]*(va2[i,j-1]-va[i,j-1])
else
vn[i,j] = aeva[i,j]*(va2[i+1,j]-va[i+1,j]) +awva[i,j]*(va2[i-1,j]-va[i-1,j]) +anva[i,j]*(va2[i,j+1]-va[i,j+1]) +asva[i,j]*(va2[i,j-1]-va[i,j-1])
end
if j==2
vs[i,j] = anva[i,j-1]*(va2[i,j]-va[i,j])
else
vs[i,j] = aeva[i,j-1]*(va2[i+1,j-1]-va[i+1,j-1]) +awva[i,j-1]*(va2[i-1,j-1]-va[i-1,j-1]) +anva[i,j-1]*(va2[i,j]-va[i,j]) +asva[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] = -dy/apua[i,j]*ue[i,j] +dy/apua[i-1,j]*uw[i,j] -dx/apva[i,j]*vn[i,j] +dx/apva[i,j-1]*vs[i,j]
end
end
# 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] = (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
rest=convergP(para, pc2, ptemp)
if rest < 1e-5
#println("rest PC2 ", rest, ", ip=", ip)
break
end
end
end
function update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apua, aeua, awua, anua, asua, apva, aeva, awva, anva, asva)
@unpack nx, ny, dx, dy, = para
# u** = u* + u'
for j=2:ny+1
for i=2:nx
ua3[i,j] = ua2[i,j] +
(aeua[i,j]*(ua2[i+1,j]-ua[i+1,j]) + awua[i,j]*(ua2[i-1,j]-ua[i-1,j]) + anua[i,j]*(ua2[i,j+1]-ua[i,j+1]) + asua[i,j]*(ua2[i,j-1]-ua[i,j-1]) -
dy*(pc2[i+1,j] - pc2[i,j]))/apua[i,j]
end
end
for j=2:ny+1
for i=2:nx
ua[i,j] = ua3[i,j]
end
end
for j=2:ny
for i=2:nx+1
va3[i,j] = va2[i,j] +
(aeva[i,j]*(va2[i+1,j]-va[i+1,j]) + awva[i,j]*(va2[i-1,j]-va[i-1,j]) + anva[i,j]*(va2[i,j+1]-va[i,j+1]) + asva[i,j]*(va2[i,j-1]-va[i,j-1]) -
dx*(pc2[i,j+1] - pc2[i,j]))/apva[i,j]
end
end
for j=2:ny
for i=2:nx+1
va[i,j] = va3[i,j]
end
end
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 solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
@unpack nx, ny, dx, dy, dt, alpha, = para
for j=2:ny+1
for i=2:nx+1
if i==2
awTa[i,j] = 2alpha*dy/dx
else
awTa[i,j] = alpha*dy/dx + 0.5dy*ua[i-1,j]
end
if i==nx+1
aeTa[i,j]=2alpha*dy/dx
else
aeTa[i,j] = alpha*dy/dx - 0.5dy*ua[i,j]
end
if j == ny+1
anTa[i,j] = 2alpha*dx/dy
else
anTa[i,j] = alpha*dx/dy - 0.5*dx*va[i,j]
end
if j == 2
asTa[i,j] = 2alpha*dx/dy
else
asTa[i,j] = alpha*dx/dy + 0.5*dx*va[i,j-1]
end
apTa[i,j] = aeTa[i,j] + awTa[i,j] + anTa[i,j] + asTa[i,j] + dx*dy/dt
bTa[i,j] = dx*dy/dt*Told[i,j]
Ttemp[i,j] = (aeTa[i,j]*Ta[i+1,j] + awTa[i,j]*Ta[i-1,j] + anTa[i,j]*Ta[i,j+1] + asTa[i,j]*Ta[i,j-1] + bTa[i,j])/apTa[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
Ta[i,j] = Ttemp[i,j]
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, Ta, uold, vold, Told, 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]+ua[1:nx, 2:ny+1])/2)
write(file, "V", (va[2:nx+1, 2:ny+1]+va[2:nx+1, 1:ny])/2)
write(file, "T", Ta[2:nx+1, 2:ny+1])
write(file, "P", 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])
write(file, "ua_latest", ua)
write(file, "va_latest", va)
write(file, "pa_latest", pa)
write(file, "Ta_latest", Ta)
write(file, "ua_old", uold)
write(file, "va_old", vold)
write(file, "Ta_old", Told)
write(file, "n", l)
end
end
function restart(para, restartFile, x, y, ua, va, pa, Ta, uold, vold, Told)
@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")
ua .= read(file, "ua_latest")
va .= read(file, "va_latest")
pa .= read(file, "pa_latest")
Ta .= read(file, "Ta_latest")
uold .= read(file, "ua_old")
vold .= read(file, "va_old")
Told .= read(file, "Ta_old")
nstart = read(file, "n")
return nstart
end
計算結果
Ra=1e8, Pr=5.3, nx=ny=512, dt=2e-3の計算結果
SIMPLE法，PISO法の理屈
有限体積法による離散化
SIMPLE法
非定常SIMPLE
PISO
もしくは
