お詫び
Jacobianの計算がおかしかったのと，速度の収束判定がおかしかったのを修正しました。
ついでに格子点の生成も修正しました。
普通有限体積法のコードだと，一般直交座標とかは使わずに，四面体とかで複雑形状に適用すると思います。それはややこしいので，ここでは一般直交座標を適用します。コロケート格子のPISO法のコードのベースはここにあります。
これに一般直交座標に適用します。カルマン渦の計算しました。下記の結果を出します。
基礎式
説明ややこしいので，下の画像を見てください。基本的には一般直交座標でのSIMPLE法の説明です。RhieとChowの論文が参考になります。
Numerical study of the turbulent flow past an airfoil with trailing edge separation
C. M. Rhie and W. L. Chow
AIAA Journal 1983 21:11, 1525-1532
グリッドの生成
まずはグリッドを作ります。グリッド生成のコード。周期境界条件を適用するため，j方向のグリッドが一部重なります。
高見，河村「偏微分方程式の差分解法」を参考にしました。
using Parameters
using HDF5
using PyPlot
@with_kw struct Param{T1,T2}
R1::T1 = 0.5
R2::T1 = 10.0
xr::T1 = 3.
ytop::T1 = 3.
ε::T1 = 1e-6
itermax::T2 = 100
nx::T2 = 150
ny::T2 = 146
outputname::String = "Grid/grid-largecirc2"
end
function gridgen(para)
@unpack nx, ny, = para
x = zeros(nx, ny)
y = zeros(nx, ny)
xT = zeros(ny, nx)
yT = zeros(ny, nx)
# control function
P = zeros(nx, ny)
Q = zeros(nx, ny)
initialize(para, x, y)
xT .= x'
yT .= y'
output(para, x, y, xT, yT)
# solve elliptic equation
ellip(para, x, y, P, Q)
xT .= x'
yT .= y'
output(para, x, y, xT, yT)
end
function initialize(para, x, y)
@unpack nx, ny, R1, R2, xr, ytop = para
# initial grid
for j = 1:ny
r = R1
θ = 2pi*(j-1)/(ny-2)
x[1, j] = r*cos(θ)
y[1, j] = r*sin(θ)
end
for j = 1:ny
for i = 2:nx
r = R1 + (R2-R1)*(i-1.5)/(nx-1)
θ = 2pi*(j-1)/(ny-2)
x[i, j] = r*cos(θ)
y[i, j] = r*sin(θ)
end
end
end
function output(para, x, y, xT, yT)
@unpack nx, ny, outputname = para
filename = outputname*".h5"
h5open(filename, "w") do file
write(file, "nx", nx)
write(file, "ny", ny)
write(file, "x", x)
write(file, "y", y)
write(file, "xT", xT)
write(file, "yT", yT)
end
end
function ellip(para, x, y, P, Q)
@unpack nx, ny, itermax, ε = para
itermax3 = div(itermax, 3)
itermax6 = itermax3*2
println(x[1:3,1:3])
cons = 0.05
for i = 1:itermax
error = 0.0
if i >= itermax3
cons = 0.25
elseif i>=itermax6
cons = 1.0
end
for l = 2:ny-1
for k = 2:nx-1
xξ = 0.5(x[k+1, l] - x[k-1, l])
yξ = 0.5(y[k+1, l] - y[k-1, l])
xη = 0.5(x[k, l+1] - x[k, l-1])
yη = 0.5(y[k, l+1] - y[k, l-1])
α = xη^2 + yη^2
β = xξ*xη + yξ*yη
γ = xξ^2 + yξ^2
J = xξ*yη - xη*yξ
PP = P[k, l]
QQ = Q[k, l]
ag = 0.5/(α + γ)
sx = ( α*(x[k+1,l] + x[k-1,l]) -
0.5β*(x[k+1,l+1] - x[k-1,l+1] - x[k+1,l-1] + x[k-1,l-1]) +
γ*(x[k,l+1] + x[k,l-1]) + J^2*(PP*xξ + QQ*xη) )*ag - x[k,l]
sy = (α*(y[k+1,l] + y[k-1,l]) -
0.5β*(y[k+1,l+1] - y[k-1,l+1] - y[k+1,l-1] + y[k-1,l-1]) +
γ*(y[k,l+1] + y[k,l-1]) + J^2*(PP*yξ + QQ*yη) )*ag - y[k,l]
x[k,l] = x[k,l] + cons*sx
y[k,l] = y[k,l] + cons*sy
error += sx^2 + sy^2
end
end
for i = 1:nx
x[i, 1] = x[i, ny-1]
y[i, 1] = y[i, ny-1]
x[i, ny] = x[i, 2]
y[i, ny] = y[i, 2]
end
if error < ε
println("iteration, error=", i, ", ", error)
break
elseif i==itermax
println("error= ",error)
end
end
end
実行と可視化
para = Param(nx=150, ny=146)
gridgen(para)
@unpack outputname = para
filename = outputname*".h5"
file = h5open(filename, "r")
x = read(file, "x")
y = read(file, "y")
xT = read(file, "xT")
yT = read(file, "yT")
nx = read(file, "nx")
ny = read(file, "ny")
close(file)
fig = plt.figure(figsize = (13, 13))
ax1 = fig.add_subplot(111)
ax1.plot(xT[:,:], yT[:,:], "*" )
#ax1.plot(x[:,:], y[:,:], color="black")
#ax1.plot(xT[:,:], yT[:,:], color="black")
#ax1.set_xlim(-1,1)
#ax1.set_ylim(-1,1)
こんなグリッドです。
流体コード
流体の解析は以下です。グリッドの点数を表す，nx, nyと数字がずれてます。これは流体は1ーnx+2, 1-ny+2まで点数があるためです。
できたばかりなので，間違っているところが含まれる可能性があるので，注意してください。
using Parameters
using HDF5
@with_kw struct Param{T1,T2}
dt::T1 = 0.01
# viscosity
nu::T1 = 0.2e-2
Uext::T1 = 1.0
nstep::T2 = 2000
alphaua::T1 = 0.9
alphava::T1 = 0.9
alphapa::T1 = 0.5
niter::T2 = 500
ipter::T2 = 25
ipter2::T2 = 100
start::T2 =1
# mesh size
nx::T2 = 147
ny::T2 = 144
ndata::T2 = 50
gridfile::String = "Grid/grid-rectcirc4.h5"
outfile::String = "Data-piso-Kar/uv"
end
function main(para)
@unpack nx, ny, nstep, niter, ndata, Uext = para
# auxiliary velocity
ua = ones(Float64, nx+2, ny+2)
ua2 = ones(Float64, nx+2, ny+2)
ua3 = ones(Float64, nx+2, ny+2)
utemp = ones(Float64, nx+2, ny+2)
uold = ones(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)
Vn = 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)
omega = zeros(Float64, nx+2, ny+2)
beta2 = zeros(Float64, nx+2, ny+2)
gam = zeros(Float64, nx+2, ny+2)
U = zeros(Float64, nx+2, ny+2)
V = zeros(Float64, nx+2, ny+2)
#dξ = zeros(Float64, nx+2, ny+2)
#dη = zeros(Float64, nx+2, ny+2)
# pressure, corrected pressure
pa = zeros(Float64, nx+2, ny+2)
pc = zeros(Float64, nx+2, ny+2)
pc2 = 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+2, ny+2)
y = zeros(Float64, nx+2, ny+2)
ξx = zeros(Float64, nx+2, ny+2)
ξy = zeros(Float64, nx+2, ny+2)
ηx = zeros(Float64, nx+2, ny+2)
ηy = zeros(Float64, nx+2, ny+2)
J = zeros(Float64, nx+2, ny+2)
rest = 0.0
readgrid(para, x, y)
jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
for l = 1:nstep
BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
for n = 1:niter
BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
makesource(bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J)
solveUA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
solveVA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, Vn, ξx, ξy, ηx, ηy, J, U, V)
RhieChowInterpolation(para, utemp, vtemp, pa, apuva, U, V, Ue, Vn, ξx, ξy, ηx, ηy, J)
#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, ξx, ξy, ηx, ηy, J,)
update(para, utemp, vtemp, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
solvePC2(para, utemp, ua2, vtemp, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J)
update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)
rest = converg(para, ua, va, ua2, va2)
if rest <=1e-7 && n >=2
if l%1 == 0
println(l, " th step, ", n, " th iteration, residual=", rest)
end
break
end
end
for j=1:ny+2
for i=1:nx+2
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)
vor(para, omega, ua, va, ξx, ξy, ηx, ηy)
output(para, ua, va, pa, uold, vold, x, y, l, omega)
end
end
#streamfunction(para, u, v, psi)
end
function main(para, restartFile)
@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)
Vn = zeros(Float64, nx+2, ny+2)
omega = 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)
beta2 = zeros(Float64, nx+2, ny+2)
gam = zeros(Float64, nx+2, ny+2)
U = zeros(Float64, nx+2, ny+2)
V = zeros(Float64, nx+2, ny+2)
dξ = zeros(Float64, nx+2, ny+2)
dη = zeros(Float64, nx+2, ny+2)
# pressure, corrected pressure
pa = zeros(Float64, nx+2, ny+2)
pc = zeros(Float64, nx+2, ny+2)
pc2 = 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+2, ny+2)
y = zeros(Float64, nx+2, ny+2)
ξx = zeros(Float64, nx+2, ny+2)
ξy = zeros(Float64, nx+2, ny+2)
ηx = zeros(Float64, nx+2, ny+2)
ηy = zeros(Float64, nx+2, ny+2)
J = zeros(Float64, nx+2, ny+2)
rest = 0.0
readgrid(para, x, y)
jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
nstart = restart(para, restartFile, ua, va, pa, uold, vold, )
println("nstart=",nstart)
for l = nstart+1:nstart+nstep
BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
for n = 1:niter
BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
makesource(bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J)
solveUA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
solveVA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, Vn, ξx, ξy, ηx, ηy, J, U, V)
RhieChowInterpolation(para, utemp, vtemp, pa, apuva, U, V, Ue, Vn, ξx, ξy, ηx, ηy, J)
#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, ξx, ξy, ηx, ηy, J,)
update(para, utemp, vtemp, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
solvePC2(para, utemp, ua2, vtemp, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J)
update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)
rest = converg(para, ua, va, ua2, va2)
if rest <=1e-7 && n >=2
if l%1 == 0
println(l, " th step, ", n, " th iteration, residual=", rest)
end
break
end
end
for j=1:ny+2
for i=1:nx+2
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)
vor(para, omega, ua, va, ξx, ξy, ηx, ηy)
output(para, ua, va, pa, uold, vold, x, y, l, omega)
end
end
#streamfunction(para, u, v, psi)
end
function BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
@unpack nx, ny, Uext, dt, start = para
for j=1:ny+2
ua[1, j] = 0.0
va[1, j] = 0.0
ua2[1,j] = 0.0
va2[1,j] = 0.0
pc[1, j] = pc[2,j]
pc2[1, j] = pc2[2,j]
pa[1, j] = pa[2,j]
end
# outlet
for j=1:2div(ny,8) +1
ua[nx+2,j] = ua[nx+1,j]
va[nx+2,j] = va[nx+1,j]
ua2[nx+2,j] = ua2[nx+1,j]
va2[nx+2,j] = va2[nx+1,j]
pa[nx+2, j] = 0 #pa[nx+1, j]
pc[nx+2, j] = 0 #pc[nx+1, j]
pc2[nx+2, j] = 0
end
# inlet
Ubound = min(Uext*l/start, Uext)
for j=2div(ny,8)+2:ny -2div(ny,8)
ua[nx+2,j] = Ubound
va[nx+2,j] = 0.0
ua2[nx+2,j] = Ubound
va2[nx+2,j] = 0.0
pa[nx+2,j] = pa[nx+1, j]
pc[nx+2,j] = pc[nx+1, j]
pc2[nx+2,j] = pc2[nx+1, j]
end
# outlet
for j=ny+1-2div(ny,8):ny+2
ua[nx+2,j] = ua[nx+1,j]
va[nx+2,j] = va[nx+1,j]
ua2[nx+2,j] = ua2[nx+1,j]
va2[nx+2,j] = va2[nx+1,j]
pa[nx+2, j] = 0 #pa[nx+1, j]
pc[nx+2, j] = 0 #pc[nx+1, j]
pc2[nx+2, j] = 0
end
# periodic condition
for i=1:nx+2
ua[i, ny+2] = ua[i, 2]
va[i, ny+2] = va[i, 2]
pc[i, ny+2] = pc[i, 2]
pa[i, ny+2] = pa[i, 2]
ua2[i, ny+2] = ua2[i, 2]
va2[i, ny+2] = va2[i, 2]
pc2[i, ny+2] = pc2[i, 2]
ua[i, 1] = ua[i, ny+1]
va[i, 1] = va[i, ny+1]
pc[i, 1] = pc[i, ny+1]
pa[i, 1] = pa[i, ny+1]
ua2[i, 1] = ua2[i, ny+1]
va2[i, 1] = va2[i, ny+1]
pc2[i, 1] = pc2[i, ny+1]
end
end
function readgrid(para, x, y)
@unpack gridfile = para
file = h5open(gridfile, "r")
x .= read(file, "x")
y .= read(file, "y")
close(file)
end
function jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
@unpack nx, ny = para
xξ = zeros(Float64, nx+2, ny+2)
yξ = zeros(Float64, nx+2, ny+2)
xη = zeros(Float64, nx+2, ny+2)
yη = zeros(Float64, nx+2, ny+2)
for j = 1:ny+2
for i = 2:nx+1
xξ[i,j] = 0.5(x[i+1,j] - x[i-1,j])
yξ[i,j] = 0.5(y[i+1,j] - y[i-1,j])
end
end
for j=1:ny+2
xξ[1,j] = -(3*x[1,j] - 4*x[2,j] + x[3,j])*0.5
yξ[1,j] = -(3*y[1,j] - 4*y[2,j] + y[3,j])*0.5
xξ[nx+2,j] = (3*x[nx+2,j] - 4*x[nx+1,j] + x[nx,j])*0.5
yξ[nx+2,j] = (3*y[nx+2,j] - 4*y[nx+1,j] + y[nx,j])*0.5
end
for j = 2:ny+1
for i = 1:nx+2
xη[i,j] = 0.5(x[i,j+1] - x[i,j-1])
yη[i,j] = 0.5(y[i,j+1] - y[i,j-1])
end
end
for i=1:nx+2
xη[i,1] = 0.5(x[i,2] - x[i,ny])
yη[i,1] = 0.5(y[i,2] - y[i,ny])
xη[i, ny+2] = 0.5(x[i,3] - x[i,ny+1])
yη[i, ny+2] = 0.5(y[i,3] - y[i,ny+1])
end
for j = 1:ny+2
for i = 1:nx+2
J[i,j] = 1.0/(xξ[i,j]*yη[i,j] - xη[i,j]*yξ[i,j])
ξx[i,j] = J[i,j]*yη[i,j]
ηx[i,j] = -J[i,j]*yξ[i,j]
ξy[i,j] = -J[i,j]*xη[i,j]
ηy[i,j] = J[i,j]*xξ[i,j]
end
end
println("metrics ", ξx[1:3,1:2], ξy[1:3,1:2], ηx[1:3,1:2], ηy[1:3,1:2], J[1:3,1:2])
end
function makesource(bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J)
@unpack nx, ny, dt, nu, = para
for j=2:ny+1
for i=2:nx+1
beta2[i,j] = nu/J[i,j]*(ξx[i,j]*ηx[i,j] + ξy[i,j]*ηy[i,j])
end
end
# bua
for j=2:ny+1
for i=2:nx+1
bua[i,j] = 0.125( (uold[i,j+1] - uold[i,j-1])*(beta2[i+1,j] - beta2[i-1,j]) + (uold[i+1,j] - uold[i-1,j])*(beta2[i,j+1] - beta2[i,j-1]) + uold[i+1,j+1]*(2beta2[i,j] + beta2[i+1,j] +beta2[i,j+1]) - uold[i-1,j+1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j+1]) - uold[i+1,j-1]*(2beta2[i,j] + beta2[i+1,j] + beta2[i,j-1]) + uold[i-1,j-1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j-1]))
bua[i,j] += -0.5(ξx[i+1,j]*pa[i+1,j]/J[i+1,j] - ξx[i-1,j]*pa[i-1,j]/J[i-1,j] + ηx[i,j+1]*pa[i,j+1]/J[i,j+1] - ηx[i,j-1]*pa[i,j-1]/J[i,j-1]) + uold[i,j]/J[i,j]/dt
end
end
#bva
for j=2:ny+1
for i=2:nx+1
bva[i,j] = 0.125( (vold[i,j+1] - vold[i,j-1])*(beta2[i+1,j] - beta2[i-1,j]) + (vold[i+1,j] - vold[i-1,j])*(beta2[i,j+1] - beta2[i,j-1]) + vold[i+1,j+1]*(2beta2[i,j] + beta2[i+1,j] +beta2[i,j+1]) - vold[i-1,j+1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j+1]) - vold[i+1,j-1]*(2beta2[i,j] + beta2[i+1,j] + beta2[i,j-1]) + vold[i-1,j-1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j-1]))
bva[i,j] += -0.5(ξy[i+1,j]*pa[i+1,j]/J[i+1,j] - ξy[i-1,j]*pa[i-1,j]/J[i-1,j] + ηy[i,j+1]*pa[i,j+1]/J[i,j+1] - ηy[i,j-1]*pa[i,j-1]/J[i,j-1]) + vold[i,j]/J[i,j]/dt
end
end
end
function makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)
@unpack nx, ny, = para
for j=1:ny+2
for i=1:nx+2
U[i,j] = ξx[i,j]*ua[i,j] + ξy[i,j]*va[i,j]
V[i,j] = ηx[i,j]*ua[i,j] + ηy[i,j]*va[i,j]
end
end
end
function solveUA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
@unpack nx, ny, dt, nu, alphaua, alphava, alphapa, niter, ipter, Uext = para
makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)
for j=1:ny+2
for i=1:nx+2
if i==1
awuva[i,j] = 0
elseif i==2
awuva[i,j] = nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i-1,j]^2 + ξy[i-1,j]^2)/J[i-1,j])
else
awuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i-1,j]^2 + ξy[i-1,j]^2)/J[i-1,j]) + 0.25*(U[i,j]/J[i,j] + U[i-1,j]/J[i-1,j])
end
if i==nx+2
aeuva[i,j] = 0
elseif i==nx+1
aeuva[i,j] = nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i+1,j]^2 + ξy[i+1,j]^2)/J[i+1,j]) - 0.25*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j])
else
aeuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i+1,j]^2 + ξy[i+1,j]^2)/J[i+1,j]) - 0.25*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j])
end
if j==1
asuva[i,j] = 0
elseif j==2
asuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,ny+1]^2 + ηy[i,ny+1]^2)/J[i,ny+1]) + 0.25*(V[i,j]/J[i,j] + V[i,ny+1]/J[i,ny+1])
else
asuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j-1]^2 + ηy[i,j-1]^2)/J[i,j-1]) + 0.25*(V[i,j]/J[i,j] + V[i,j-1]/J[i,j-1])
end
if j==ny+2
anuva[i,j] = 0
elseif j == ny+1
anuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,2]^2 + ηy[i,2]^2)/J[i,2]) - 0.25*(V[i,2]/J[i,2] + V[i,j]/J[i,j])
else
anuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j+1]^2 + ηy[i,j+1]^2)/J[i,j+1]) - 0.25*(V[i,j+1]/J[i,j+1] + V[i,j]/J[i,j])
end
apuva[i,j] = aeuva[i,j] + awuva[i,j] + anuva[i,j] + asuva[i,j] + 1/dt/J[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
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=1:ny+2
utemp[1,j] = ua[1,j]
utemp[nx+2,j] = ua[nx+2,j]
end
for i=1:nx+2
utemp[i,1] = ua[i,1]
utemp[i,ny+2] = ua[i,ny+2]
end
#for j=2:ny+1
# for i=2:nx+1
# dξ[i,j] = (ξx[i,j]^2+ξy[i,j]^2)/apuva[i,j]/J[i,j]
# dη[i,j] = (ηx[i,j]^2+ηy[i,j]^2)/apuva[i,j]/J[i,j]
# end
#end
end
function solveVA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, Vn, ξx, ξy, ηx, ηy, J, U, V)
@unpack nx, ny, dt, nu, alphaua, alphava, alphapa, niter, ipter, = para
for j=2:ny+1
for i=2:nx+1
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=1:ny+2
vtemp[1,j] = va[1,j]
vtemp[nx+2,j] = va[nx+2,j]
end
for i=1:nx+2
vtemp[i,1] = va[i,1]
vtemp[i,ny+2] = va[i,ny+2]
end
end
function RhieChowInterpolation(para, utemp, vtemp, pa, apuva, U, V, Ue, Vn, ξx, ξy, ηx, ηy, J)
@unpack nx, ny, = para
makeCotravariant(para, utemp, vtemp, U, V, ξx, ξy, ηx, ηy)
# ue: interface velocity
for j=2:ny+1
for i=2:nx
#Ue[i,j] = 0.5(U[i+1,j]+U[i,j]) + 0.5(dξ[i,j]+dξ[i+1,j])*(pa[i,j]-pa[i+1,j]) -0.25dξ[i,j]*(pa[i-1,j]-pa[i+1,j]) -0.25*dξ[i+1,j]*(pa[i,j]-pa[i+2,j])
Ue[i,j] = 0.5(U[i+1,j]+U[i,j]) -0.5((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])/apuva[i,j] + (ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])/apuva[i+1,j])*(pa[i+1,j]-pa[i,j]) +
0.25((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])*pa[i+1,j] -(ξx[i,j]*ξx[i-1,j]/J[i-1,j]+ξy[i,j]*ξy[i-1,j]/J[i-1,j])*pa[i-1,j] )/apuva[i,j] +
0.25((ξx[i+1,j]*ξx[i+2,j]/J[i+2,j]+ξy[i+1,j]*ξy[i+2,j]/J[i+2,j])*pa[i+2,j] -(ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])*pa[i,j] )/apuva[i+1,j]
end
end
for j=2:ny+1
Ue[1,j] = U[1,j]
end
for j=2:ny+1
Ue[nx+1,j] = U[nx+2,j]
end
for j=2:ny
for i=2:nx+1
#Vn[i,j] = 0.5(V[i,j+1]+V[i,j]) + 0.5(dη[i,j]+dη[i,j+1])*(pa[i,j]-pa[i,j+1]) -0.25dη[i,j]*(pa[i,j-1]-pa[i,j+1]) -0.25dη[i,j+1]*(pa[i,j]-pa[i,j+2])
Vn[i,j] = 0.5(V[i,j+1]+V[i,j]) -0.5((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])/apuva[i,j] + (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])/apuva[i,j+1] )*(pa[i,j+1]-pa[i,j]) +
0.25((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])*pa[i,j+1] - (ηx[i,j]*ηx[i,j-1]/J[i,j-1]+ηy[i,j]*ηy[i,j-1]/J[i,j-1])*pa[i,j-1])/apuva[i,j] +
0.25((ηx[i,j+1]*ηx[i,j+2]/J[i,j+2]+ηy[i,j+1]*ηy[i,j+2]/J[i,j+2])*pa[i,j+2] - (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])*pa[i,j])/apuva[i,j+1]
end
end
for i=2:nx+1
j=ny+1
#Vn[i,j] = 0.5(V[i,2]+V[i,j]) + 0.5(dη[i,j]+dη[i,2])*(pa[i,j]-pa[i,2]) -0.25dη[i,j]*(pa[i,j-1]-pa[i,2]) -0.25dη[i,2]*(pa[i,j]-pa[i,3])
Vn[i,j] = 0.5(V[i,2]+V[i,j]) -0.5((ηx[i,j]*ηx[i,2]/J[i,2]+ηy[i,j]*ηy[i,2]/J[i,2])/apuva[i,j] + (ηx[i,2]*ηx[i,j]/J[i,j]+ηy[i,2]*ηy[i,j]/J[i,j])/apuva[i,2] )*(pa[i,2]-pa[i,j]) +
0.25((ηx[i,j]*ηx[i,2]/J[i,2]+ηy[i,j]*ηy[i,2]/J[i,2])*pa[i,2] - (ηx[i,j]*ηx[i,j-1]/J[i,j-1]+ηy[i,j]*ηy[i,j-1]/J[i,j-1])*pa[i,j-1])/apuva[i,j] +
0.25((ηx[i,2]*ηx[i,3]/J[i,2]+ηy[i,2]*ηy[i,3]/J[i,3])*pa[i,3] - (ηx[i,2]*ηx[i,j]/J[i,j]+ηy[i,2]*ηy[i,j]/J[i,j])*pa[i,j])/apuva[i,2]
end
for i=2:nx+1
Vn[i,1] = Vn[i,ny+1]
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, ξx, ξy, ηx, ηy, J)
@unpack nx, ny, nu, alphaua, alphava, alphapa, niter, ipter, nu = 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
else
awpc[i,j] = 0.5((ξx[i-1,j]*ξx[i,j]/J[i,j]+ξy[i-1,j]*ξy[i,j]/J[i,j])/apuva[i-1,j] + (ξx[i,j]*ξx[i-1,j]/J[i-1,j]+ξy[i,j]*ξy[i-1,j]/J[i-1,j])/apuva[i,j])/(0.5(J[i-1,j]+J[i,j]))
end
if i == nx+1
aepc[i,j] = 0
else
aepc[i,j] = 0.5((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])/apuva[i,j] + (ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])/apuva[i+1,j])/(0.5(J[i,j]+J[i+1,j]))
end
if j == 2
aspc[i,j] = 0.5((ηx[i,ny+1]*ηx[i,j]/J[i,j]+ηy[i,ny+1]*ηy[i,j]/J[i,j])/apuva[i,ny+1] + (ηx[i,j]*ηx[i,ny+1]/J[i,ny+1]+ηy[i,j]*ηy[i,ny+1]/J[i,ny+1])/apuva[i,j] )/(0.5(J[i,ny+1]+J[i,j]))
else
aspc[i,j] = 0.5((ηx[i,j-1]*ηx[i,j]/J[i,j]+ηy[i,j-1]*ηy[i,j]/J[i,j])/apuva[i,j-1] + (ηx[i,j]*ηx[i,j-1]/J[i,j-1]+ηy[i,j]*ηy[i,j-1]/J[i,j-1])/apuva[i,j] )/(0.5(J[i,j-1]+J[i,j]))
end
if j == ny+1
anpc[i,j] = 0.5((ηx[i,j]*ηx[i,2]/J[i,2]+ηy[i,j]*ηy[i,2]/J[i,2])/apuva[i,j] + (ηx[i,2]*ηx[i,j]/J[i,j]+ηy[i,2]*ηy[i,j]/J[i,j])/apuva[i,2] )/(0.5(J[i,j]+J[i,2]))
else
anpc[i,j] = 0.5((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])/apuva[i,j] + (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])/apuva[i,j+1] )/(0.5(J[i,j]+J[i,j+1]))
end
appc[i,j] = aepc[i,j] + awpc[i,j] + anpc[i,j] + aspc[i,j]
bpc[i,j] = -(Ue[i,j]/(0.5(J[i,j]+J[i+1,j])) - Ue[i-1,j]/(0.5(J[i,j]+J[i-1,j]))) - (Vn[i,j]/(0.5(J[i,j]+J[i,j+1])) - Vn[i,j-1]/(0.5(J[i,j]+J[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] = pc[i,j]*(1.0 - alphapa) + alphapa*(
aepc[i,j]*pc[i+1,j] + awpc[i,j]*pc[i-1,j] + anpc[i,j]*pc[i,j+1] + aspc[i,j]*pc[i,j-1] + bpc[i,j])/appc[i,j]
end
end
for j=1:ny+2
pc[1,j] = pc[2,j]
end
for j=1:2div(ny,8) +1
pc[nx+2,j] = 0 #pc[nx+1, j]
end
for j=2div(ny,8)+2:ny -2div(ny,8)
#for i=1:nx+2
pc[nx+2,j] = pc[nx+1, j]
end
for j=ny+1-2div(ny,8):ny+2
pc[nx+2,j] = 0 #pc[nx+1, j]
end
for i=1:nx+2
pc[i, ny+2] = pc[i, 2]
pc[i, 1] = pc[i, ny+1]
end
rest=convergP(para, pc, ptemp)
if rest < 1e-5 && ip >=2
#println("converged rest PC ", rest, ", ip=", ip)
break
end
if ip ==ipter
println("rest PC ", rest, ", ip=", ip)
end
end
end
function update(para, utemp, vtemp, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
@unpack nx, ny, alphapa, = para
# u = u* + u'
for j=2:ny+1
for i=2:nx+1
ua2[i,j] = utemp[i,j] -0.5(ξx[i+1,j]*pc[i+1,j]/J[i+1,j] - ξx[i-1,j]*pc[i-1,j]/J[i-1,j] + ηx[i,j+1]*pc[i,j+1]/J[i,j+1] - ηx[i,j-1]*pc[i,j-1]/J[i,j-1])/apuva[i,j]
va2[i,j] = vtemp[i,j] -0.5(ξy[i+1,j]*pc[i+1,j]/J[i+1,j] - ξy[i-1,j]*pc[i-1,j]/J[i-1,j] + ηy[i,j+1]*pc[i,j+1]/J[i,j+1] - ηy[i,j-1]*pc[i,j-1]/J[i,j-1])/apuva[i,j]
end
end
#pref = pa[2,2]
#for j=1:ny+2
# for i=1:nx+2
# pa[i,j] -= pref
# end
#end
end
function solvePC2(para, utemp, ua2, vtemp, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J)
@unpack nx, ny, ipter2, alphapa, nu = para
makeCotravariant(para, ua2 .- utemp, va2 .- vtemp, U, V, ξx, ξy, ηx, ηy)
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]*(U[i,j])
else
UE[i,j] = aeuva[i+1,j]*(U[i+2,j]) +awuva[i+1,j]*(U[i,j]) +anuva[i+1,j]*(U[i+1,j+1]) +asuva[i+1,j]*(U[i+1,j-1])
end
if i==2
UW[i,j] = aeuva[i-1,j]*(U[i,j])
else
UW[i,j] = aeuva[i-1,j]*(U[i,j]) +awuva[i-1,j]*(U[i-2,j]) +anuva[i-1,j]*(U[i-1,j+1]) +asuva[i-1,j]*(U[i-1,j-1])
end
if j==ny+1
VN[i,j] = asuva[i,j+1]*(V[i,j])
else
VN[i,j] = aeuva[i,j+1]*(V[i+1,j+1]) +awuva[i,j+1]*(V[i-1,j+1]) +anuva[i,j+1]*(V[i,j+2]) +asuva[i,j+1]*(V[i,j])
end
if j==2
VS[i,j] = anuva[i,j-1]*(V[i,j])
else
VS[i,j] = aeuva[i,j-1]*(V[i+1,j-1]) +awuva[i,j-1]*(V[i-1,j-1]) +anuva[i,j-1]*(V[i,j]) +asuva[i,j-1]*(V[i,j-2])
end
end
end
for j=2:ny+1
for i=2:nx+1
bpc[i,j] = -0.5(UE[i,j]/apuva[i+1,j]/(0.5(J[i+1,j]+J[i,j])) - UW[i,j]/apuva[i-1,j]/(0.5(J[i,j]+J[i-1,j]))) -0.5(VN[i,j]/apuva[i,j+1]/(0.5(J[i,j]+J[i,j+1])) -VS[i,j]/apuva[i,j-1]/(0.5(J[i,j]+J[i,j-1])))
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:ipter2
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] = pc2[i,j]*(1.0 - alphapa) + alphapa*(aepc[i,j]*pc2[i+1,j] + awpc[i,j]*pc2[i-1,j] + anpc[i,j]*pc2[i,j+1] + aspc[i,j]*pc2[i,j-1] + bpc[i,j])/appc[i,j]
end
end
for i=1:nx+2
pc2[i, 1] = pc2[i, ny+1]
pc2[i, ny+2] = pc2[i, 2]
end
for j=1:ny+2
pc2[1, j] = pc2[2, j]
end
for j=1:2div(ny,8) +1
pc2[nx+2, j] = 0
end
for j=2div(ny,8)+2:ny -2div(ny,8)
pc2[nx+2,j] = pc2[nx+1, j]
end
for j=ny+1-2div(ny,8):ny+2
pc2[nx+2, j] = 0
end
rest=convergP(para, pc2, ptemp)
if rest < 1e-7 && ip>=2
#println("converged rest PC2 ", rest, ", ip=", ip)
break
end
if ip == ipter2
println("rest PC2 ", rest, ", ip=", ip)
end
end
end
function update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)
@unpack nx, ny, = para
# 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]-utemp[i+1,j]) + awuva[i,j]*(ua2[i-1,j]-utemp[i-1,j]) + anuva[i,j]*(ua2[i,j+1]-utemp[i,j+1]) + asuva[i,j]*(ua2[i,j-1]-utemp[i,j-1]) -
0.5(ξx[i+1,j]*pc2[i+1,j]/J[i+1,j] - ξx[i-1,j]*pc2[i-1,j]/J[i-1,j] + ηx[i,j+1]*pc2[i,j+1]/J[i,j+1] - ηx[i,j-1]*pc2[i,j-1]/J[i,j-1]))/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]-vtemp[i+1,j]) + awuva[i,j]*(va2[i-1,j]-vtemp[i-1,j]) + anuva[i,j]*(va2[i,j+1]-vtemp[i,j+1]) + asuva[i,j]*(va2[i,j-1]-vtemp[i,j-1]) -
0.5(ξy[i+1,j]*pc2[i+1,j]/J[i+1,j] - ξy[i-1,j]*pc2[i-1,j]/J[i-1,j] + ηy[i,j+1]*pc2[i,j+1]/J[i,j+1] - ηy[i,j-1]*pc2[i,j-1]/J[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
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 = (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, uold, vold, x, y, l, omega)
@unpack outfile, dt, nx, ny = para
#println("time=", time)
filename = outfile*string(l)*".h5"
h5open(filename, "w") do file
write(file, "U", ua)
write(file, "V", va)
write(file, "P", pa)
write(file, "omega", omega)
#write(file, "Psi", psi[2:nx+1, 2:ny+1])
write(file, "X", x[1:nx+2, 1:ny+2])
write(file, "Y", y[1:nx+2, 1:ny+2])
write(file, "ua_old", uold)
write(file, "va_old", vold)
write(file, "n", l)
end
end
function restart(para, restartFile, ua, va, pa, uold, vold, )
@unpack nx, ny, = para
file = h5open(restartFile, "r")
ua .= read(file, "U")
va .= read(file, "V")
pa .= read(file, "P")
uold .= read(file, "ua_old")
vold .= read(file, "va_old")
nstart = read(file, "n")
close(file)
return nstart
end
function vor(para, omega, ua, va, ξx, ξy, ηx, ηy)
@unpack nx, ny, = para
for j=2:ny+1
for i=2:nx+1
omega[i,j] = (ξx[i,j]*0.5*(va[i+1,j] -va[i-1,j]) + ηx[i,j]*0.5*(va[i,j+1] -va[i,j-1])) -
(ξy[i,j]*0.5*(ua[i+1,j] -ua[i-1,j]) + ηy[i,j]*0.5*(ua[i,j+1] -ua[i,j-1]))
end
end
for i=2:nx+1
j=1
omega[i,j] = (ξx[i,j]*0.5*(va[i+1,j] -va[i-1,j]) + ηx[i,j]*0.5*(va[i,j+1] -va[i,ny])) -
(ξy[i,j]*0.5*(ua[i+1,j] -ua[i-1,j]) + ηy[i,j]*0.5*(ua[i,j+1] -ua[i,ny]))
j=ny+2
omega[i,j] = (ξx[i,j]*0.5*(va[i+1,j] -va[i-1,j]) + ηx[i,j]*0.5*(va[i,3] -va[i,j-1])) -
(ξy[i,j]*0.5*(ua[i+1,j] -ua[i-1,j]) + ηy[i,j]*0.5*(ua[i,3] -ua[i,j-1]))
end
end
#para = Param(nu=0.02, nstep=500, outfile="Data-simple/uv", nx=30, ny=30, ndata=50, alphapa=0.5)
#para = Param(nu=0.02, nstep=1, outfile="Data-simple/uv", nx=30, ny=30, ndata=1, niter=5000)
#@time main(para)
実行と可視化
include("Unsteady-PISO_colocation-gene-KarmanBC.jl")
para = Param(nu=0.002, nx=148, ny=144, gridfile="Grid/grid-largecirc2.h5", outfile="Data-piso-Kar-largeCircle/uv", ndata=50, niter=500)
@time main(para)
using HDF5
using PyPlot
dir="Data-piso-Kar-largeCircle/"
name = "uv7000"
filename = dir*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")
omega = read(file, "omega")
close(file)
skipx = 4
skipy = 1
X = x[:,1]
Y = y[1,:]
fig = plt.figure(figsize = (15, 20))
ax1 = fig.add_subplot(321)
heatmap = ax1.pcolormesh(x, y, u, cmap="jet", )
fig.colorbar(heatmap, ax=ax1)
ax1.set_title("u")
#heatmap.set_clim(vmin=0, vmax=1)
ax1.set_xlim(-3,9)
ax1.set_ylim(-6,6)
ax2 = fig.add_subplot(322)
heatmap2 = ax2.pcolormesh(x, y, v, cmap="jet", )
fig.colorbar(heatmap2, ax=ax2)
ax2.set_title("v")
#heatmap2.set_clim(vmin=-1, vmax=1)
ax2.set_xlim(-3,9)
ax2.set_ylim(-6,6)
ax3 = fig.add_subplot(323)
umag = sqrt.(u.*u .+ v.*v)
#image = ax3.quiver(x[:,40:40], y[:,40:40], u[:,40:40], v[:,40:40], scale=4, cmap="jet")
image = ax3.quiver(x[1:skipx:size(x)[1], 1:skipy:size(x)[2]], y[1:skipx:size(x)[1], 1:skipy:size(x)[2]], u[1:skipx:size(x)[1], 1:skipy:size(x)[2]], v[1:skipx:size(x)[1], 1:skipy:size(x)[2]],
umag[1:skipx:size(x)[1], 1:skipy:size(x)[2]], scale=20, cmap="jet")
#image = ax3.streamplot(x', y', u, v)
cb=fig.colorbar(image, ax=ax3)
ax3.set_title("u,v vector")
ax3.set_xlim(-3,9)
ax3.set_ylim(-6,6)
ax4 = fig.add_subplot(324)
heatmap3 = ax4.pcolormesh(x, y, pa, cmap="viridis", )
fig.colorbar(heatmap3, ax=ax4)
ax4.set_title("p")
#heatmap3.set_clim(vmin=0.35, vmax=0.5)
ax4.set_xlim(-3,9)
ax4.set_ylim(-6,6)
ax5 = fig.add_subplot(325)
heatmap4 = ax5.pcolormesh(x, y, omega, cmap="bwr", )
fig.colorbar(heatmap4, ax=ax5)
ax5.set_title("vorticity")
heatmap4.set_clim(vmin=-1, vmax=1)
ax5.set_xlim(-3,9)
ax5.set_ylim(-6,6)
plt.savefig(dir*"Fig/"*name*".pdf")