球形粒子の運動におけるサフマン揚力の効果

#(parameter.jl)

# 流路形状
    leng = 0.1  # 流路長さ
    hb = 0.05   # 流路半幅

# 粒子=空気、流体=水
    rhop = 1.2; rhof = 1000; muf = 1e-3
# 粒子=水、流体=空気
#   rhop = 1000; rhof = 1.2; muf = 1.8e-5

# 粒子回転角速度(マグナス効果)
    omgp = 0

# 粒径
    dp = 1e-5

# 粒子放出位置
    xini = -0.04

# 粒子放出速度
# (1) directly specify velocity value
    vini = 0
# (2) set same value as velocity-profile
#   vini = nothing

# 流速分布タイプ
    type = "uniform"
#   type = "linear"
#   type = "parabolla"

# 最大流速
#   vmax = 1.0
    vmax = 0.0
#   vmax = -1.0

#(eqn.jl)

using DifferentialEquations
using Plots

struct Geometry{T<:AbstractFloat}
    length::T
    hwidth::T
end

@kwdef struct Ufield{T<:AbstractFloat}
    vmax::T
    profile::Function
    saffman::Function
end

#----------------------------------
function morsi_alexander(Re)
    if     Re < 0.1
        CdRe = 24;
    elseif Re < 1.0
        CdRe = (3.690 + (22.73/Re) + (0.0903/(Re*Re)))*Re;
    elseif Re < 10.0
        CdRe = (1.222 + (29.1667/Re) + (-3.8889/(Re*Re)))*Re;
    elseif Re < 100.0
        CdRe = (0.6167 + (46.50/Re) + (-116.67/(Re*Re)))*Re;
    elseif Re < 1000.0
        CdRe = (0.3644 + (98.33/Re) + (-2778/(Re*Re)))*Re;
    elseif Re < 5000.0
        CdRe = (0.357 + (148.62/Re) + (-47500/(Re*Re)))*Re;
    elseif Re < 10000.0
        CdRe = (0.46 + (-490.546/Re) + (578700/(Re*Re)))*Re;
    else
        CdRe = (0.5191 + (-1662.5/Re) + (5416700/(Re*Re)))*Re;
    end
    return CdRe
end

#----------------------------------
function get_cl(dp,udif,omgp)
    sp = (udif<eps() ? 0 : 0.5*dp*abs(omgp)/udif)
    return 2*sp
end

#----------------------------------
function get_cs(dp,udif,rhof,muf)
    re = 0.5*dp*udif*rhof/muf
    cs = exp(-re/20.0)
    return cs
end


#----------------------------------
function particle!(du,u,p,t)
    dp, coefp, coeff, ratio, rhop, rhof, muf, omgp, gravity, ufluid = p

#  u[1] --> x-coord
#  u[2] --> x-velocity
#  u[3] --> y-coord
#  u[4] --> y-velocity

    xp, up, yp, vp = u
    uf, vf = ufluid.profile(xp, yp)
    udif = sqrt((uf-up)^2 + (vf-vp)^2)
    Re = rhof * dp * udif / muf
    cdre = morsi_alexander(Re)
    taup = rhop*dp^2/(18*muf) * 24/cdre

    gravx, gravy = gravity()

    cl = get_cl(dp,udif,omgp)
    cs = get_cs(dp,udif,rhof,muf)

    (saffx, saffy) = ufluid.saffman(xp,up, yp,vp, uf,vf)
    (magnx, magny) = ( (vf-vp)*sign(omgp),
                      -(uf-up)*sign(omgp) )

    du[1] = up
    du[2] = (
               (uf-up)/taup
             + 0.75*cl/dp*ratio*udif*magnx
             + 9.66*cs/(pi*dp)*sqrt(muf*rhof)/rhop*saffx
             + (1-ratio)*gravx
            ) / coefp
    du[3] = vp
    du[4] = (
               (vf-vp)/taup
             + 0.75*cl/dp*ratio*udif*magny
             + 9.66*cs/(pi*dp)*sqrt(muf*rhof)/rhop*saffy
             + (1-ratio)*gravy
            ) / coefp

end



#----------------------------------
function track(params)
(cm,rhop,rhof,muf,dp,omgp,inipos,inivel,gravity,ufluid,geom) = params

hb = geom.hwidth
leng = geom.length

ratio = rhof/rhop

coefp = 1 + ratio*cm
coeff = ratio*(1 + cm)

(xini, yini) = inipos()
(xvel, yvel) = inivel(xini, yini)

u0 = [ xini, xvel, yini, yvel ]
tspan = (0.0, 1e4)

# when particle hit outlet face, terminate calculation
# also, particle will reflect when hitting side wall
function condition(out, u, t, integrator)
    out[1] = leng-abs(u[3])
    out[2] = (hb-u[1])*(hb+u[1])
end
function affect!(integrator, idx)
    if idx == 1
        terminate!(integrator)
    elseif idx == 2
        integrator.u[2] = -integrator.u[2]
    end
end
cb = VectorContinuousCallback(condition, affect!, 2)

p = [ dp, coefp, coeff, ratio, rhop, rhof, muf, omgp, gravity, ufluid ]
prob = ODEProblem(particle!, u0, tspan, p)

sol = solve(prob, callback=cb )
(nrow,ncol) = size(sol)
lasttime = sol.t[ncol]

times = Vector{Float64}(undef,ncol)
xps   = Vector{Float64}(undef,ncol)
ups   = Vector{Float64}(undef,ncol)
yps   = Vector{Float64}(undef,ncol)
vps   = Vector{Float64}(undef,ncol)

for j=1:ncol
    time = sol.t[j]
    xp   = sol.u[j][1]
    up   = sol.u[j][2]
    yp   = sol.u[j][3]
    vp   = sol.u[j][4]

    times[j] = time
    xps[j]   = xp
    ups[j]   = up
    yps[j]   = yp
    vps[j]   = vp
end

out = open( "outfile", "w")
println(out,"time  xps  ups  yps  vps")

for j=1:ncol
    println(out,times[j],"  ",xps[j],"  ",ups[j],"  ",yps[j],"  ",vps[j])
end
close(out)

# animation (trajectory)
nframes = 100
anim = @animate for t in LinRange(0, lasttime, nframes)
    plt = plot(sol, idxs = (1,3), tspan = (0.0, t), label = false)
    scatter!(sol, idxs = (1,3), tspan = (t, t), lab = nothing)

    plot(plt, layout = (1,1),
            xaxis = ("x" ,(-0.05, 0.01)),yaxis = ("y" ,(0, 0.1)),
            plot_title = "t = $(round(t, digits = 2))")
end
gif(anim, "trajectory.gif"; fps = 10)


# trajectory plot
plot(sol,idxs=(1,3),   title="trajectory",
        markershape=:circle, markersize=1, plotdensity=100,
        xlabel="x[m]", ylabel="y[m]", label=nothing )
savefig("trajectory.png")


# history plot
plot(sol,idxs=(0,1),   title="history of x-position",
    markershape=:circle, markersize=1, plotdensity=100,
    xlabel="time[sec]", ylabel="x[m]", label="xp")
savefig("hist-xp.png")

plot(sol,idxs=(0,2),   title="history of x-velocity",
    markershape=:circle, markersize=1, plotdensity=100,
    xlabel="time[sec]", ylabel="u[m/s]", label="up")
savefig("hist-up.png")

plot(sol,idxs=(0,3),   title="history of y-position",
    markershape=:circle, markersize=1, plotdensity=100,
    xlabel="time[sec]", ylabel="y[m]", label="yp")
savefig("hist-yp.png")

plot(sol,idxs=(0,4),   title="history of y-velocity",
    markershape=:circle, markersize=1, plotdensity=100,
    xlabel="time[sec]", ylabel="v[m/s]", label="vp")
savefig("hist-vp.png")

end

#----------------------------------
function run()

# all simulation parameters are imported from parameter.jl
include("parameter.jl")

if @isdefined xinis
    yinis = fill(0.0, length(xinis))
    inipos = () -> ( xinis, yinis )
else
    inipos = () -> ( xini, 0.0 )
end

if vini == nothing
    inivel = (xp,yp) -> ( ufluid.profile(xp,yp) )
else
    inivel = (xp,yp) -> ( 0.0, vini )
end

cm = 0.5  # virtual mass coef. for sphere
geom = Geometry( leng, hb )
gravity = () -> ( 0.0, -9.81 )

if type == "uniform"
    ufluid = Ufield(
                vmax = vmax,
                profile = (xp,yp) -> (0.0, vmax),
                saffman = (xp,up, yp,vp, uf,vf) -> (0,0),
              )
elseif type == "linear"
    ufluid = Ufield(
                vmax = vmax,
                profile = (xp,yp) -> (0.0, vmax/(2*hb)*(xp+hb)),
                saffman = (xp,up, yp,vp, uf,vf) -> begin
                    omega = 0.5*vmax/hb
                    bunbo = sqrt(abs(omega))
                    return ( (vf-vp)*omega,
                            -(uf-up)*omega )
                end
              )
elseif type == "parabolla"
    ufluid = Ufield(
                vmax = vmax,
                profile = (xp,yp) -> (0.0, vmax*(1.0 - (xp/hb)^2) ),
                saffman = (xp,up, yp,vp, uf,vf) -> begin
                    bunbo = sqrt(2*abs(vmax*xp))/hb
                    if bunbo < eps()
                        return ( 0.0, 0.0 )
                    else
                        return ( -2*vmax*(vf-vp)*xp/hb^2 /bunbo,
                                  2*vmax*(uf-up)*xp/hb^2 /bunbo ) 
                    end
                end
              )
end

params = (cm,rhop,rhof,muf,dp,omgp,inipos,inivel,gravity,ufluid,geom)
track(params)

end

run()