球形粒子の運動におけるマグナス揚力の効果
By M.Sato
概要
以前のブログ記事においては、 球形粒子に作用するサフマン揚力の影響を調べた。
今回は、球形粒子にマグナス揚力が作用している場合の効果を調べた。
サフマン揚力は流速の空間勾配がある場合に発生する揚力であるのに対して、 マグナス揚力は粒子自体が回転することで発生する揚力である。
解析の結果、粒子の回転運動が継続する時間は、粒子特性、流体特性値から推定される緩和時間に依存しており、 微小粒子の場合、粒子の回転は急速に減衰するためマグナス揚力の効果は、ほとんど生じないことがわかった。
また、マグナス揚力は粒子と流体の相対速度と粒子回転速度の外積で定義される方向に作用することを再現できた。
粒子の並進運動を記述する方程式
流体中を移動する球形粒子の運動方程式は次式で表される1。
$$ \begin{align} \frac{\pi d^3}{6}(\rho_p+C_M \rho) \frac{d \mathbf{v}}{d t} &= \frac{\pi d^3}{6}\rho (1+C_M)\frac{D \mathbf{u}}{D t} \\ & + \frac{1}{8}C_D\pi d^2 \rho |\mathbf{u}-\mathbf{v}|(\mathbf{u}-\mathbf{v}) \\ & + \frac{3}{2}d^2 \sqrt{\pi \rho \mu} \int_{t_0}^{t} \frac{1}{\sqrt{t-\tau}}\left(\frac{D\mathbf{u}}{D\tau}-\frac{d\mathbf{v}}{d\tau}\right)d\tau \\ & + \frac{1}{8} C_{RL}\pi d^2 \rho |\mathbf{u}-\mathbf{v}|\frac{(\mathbf{u}-\mathbf{v})\times \mathbf{\omega}_p}{|\mathbf{\omega}_p|} \\ & + 1.61 C_s d^2 \sqrt{\mu\rho} \frac{(\mathbf{u}-\mathbf{v})\times\mathbf{\omega}}{\sqrt{|\mathbf{\omega}|}} \\ & + \frac{\pi d^3}{6}(\rho_p-\rho)\mathbf{g} \end{align} $$
ここで、
- $d$：粒子径
- $\rho_p$：粒子密度
- $\rho$：流体密
- $\mu$：流体粘
- $\mathbf{u}$：流体速度
- $\mathbf{v}$：粒子速度
- $\omega$：渦度$(=\nabla\times\mathbf{u})$
- $\omega_p$：粒子回転速度
- $\mathbf{g}$：重力加速度
- $C_M$：付加質量係数(=0.5)
- $C_D$：抗力係数
- $C_{RL}$：マグナス揚力係数
- $C_s$：サフマン揚力係数
である。
粒子運動方程式の右辺項についてわかりやすいように(1)から(6)までの番号をつけている。 これらの項は、それぞれ、以下の内容を表している。
- (1) 流体の非定常運動に起因する力
- (2) 粒子の周囲を流体が流れることで発生する抵抗力
- (3) バセット項（物体周囲の境界層の発達が時間遅れを伴うことに起因する力）
- (4) マグナス項（粒子の回転に伴って発生する揚力）
- (5) サフマン項（流速勾配がある流れ場で発生する揚力）
- (6) 浮力項（流体と粒子の密度差に起因する力）
粒子運動方程式において左辺項が $d\mathbf{v}/dt$となるように変形すると以下のようになる。
$$ \begin{align} \frac{d \mathbf{v}}{d t} &= \left[ \frac{\rho}{\rho_p}(1+C_M)\frac{D \mathbf{u}}{D t} \right. \notag \\ & \qquad + \frac{\mathbf{u}-\mathbf{v}}{\tau_1} \notag \\ & \qquad + \frac{9}{d \rho_p} \sqrt{\frac{\rho \mu}{\pi}} \int_{t_0}^{t} \frac{1}{\sqrt{t-\tau}}\left(\frac{D\mathbf{u}}{D\tau}-\frac{d\mathbf{v}}{d\tau}\right)d\tau \notag \\ & \qquad + \frac{3}{4d} C_{RL} \frac{\rho}{\rho_p} |\mathbf{u}-\mathbf{v}|\frac{(\mathbf{u}-\mathbf{v})\times \mathbf{\omega}_p}{|\mathbf{\omega}_p|} \notag \\ & \qquad + 9.66 C_s \frac{ \sqrt{\mu\rho}}{\pi d \rho_p} \frac{(\mathbf{u}-\mathbf{v})\times\mathbf{\omega}}{\sqrt{|\mathbf{\omega}|}} \notag \\ & \qquad +\left. \left(1- \frac{\rho}{\rho_p} \right)\mathbf{g} \right] / \left(1+C_M \frac{\rho}{\rho_p} \right) \notag \end{align} $$
ここで、$\tau_1$は粒子緩和時間(並進運動)、$Re_p$は粒子レイノルズ数であり以下のように表される。
$$ \tau_1 = \frac{\rho_p d^2}{18\mu} \frac{24}{C_D R_{ep}} \ Re_p = \frac{\rho d |\mathbf{u}-\mathbf{v}|}{\mu} $$
以下の解析では、流体慣性力は無し($D\mathbf{u}/Dt=0$)とし、 バセット項、サフマン項を除外し、右辺には抗力項、マグナス項、浮力項のみ が含まれるようにした。
粒子の回転運動を記述する方程式
球形粒子が回転している場合、その回転運動は角運動量の保存則で記述される。
以下に球粒子の角運動量方程式を示す。右辺は球に作用しているトルクを表している。
$$ \begin{align} I_p \frac{d \vec{\omega_p}}{dt} &= T = \frac{\rho_f}{2} \left( \frac{d_p}{2}\right)^5 C_\omega |\vec{\omega_f}- \vec{\omega_p}|(\vec{\omega_f}- \vec{\omega_p}) \notag \end{align} $$
ここで、
- $\vec{\omega_p}$：球の回転角速度
- $I_p$：球の慣性モーメント
- $T$：作用トルク
- $\vec{\omega_f}(=1/2\cdot\nabla\times u_f)$：流体粒子の回転角速度
- $C_\omega$：回転に関する抵抗係数
である。 球粒子の場合、慣性モーメント$I_p$は次式で表される。
$$ \begin{align} I_p &= \frac{\pi}{60} \rho_p d_p^5 \notag \end{align} $$
粒子運動方程式において左辺項が $d \vec{\omega_p}/dt$となるように変形すると以下のようになる。
$$ \begin{align} \frac{d \vec{\omega_p}}{dt} &= \frac{ \vec{\omega_f}- \vec{\omega_p} }{\tau_2} \notag \end{align} $$
ここで、$\tau_2$は粒子緩和時間(回転運動)、$Re_\omega$は粒子回転レイノルズ数であり 以下のように表される。
$$ \begin{align} \tau_2 &= \frac{\rho_p d_p^2}{60\mu_f}\frac{16\pi}{C_\omega Re_\omega} \notag \\ Re_\omega &= \frac{\rho_f |\vec{\omega_p}| d_p^2}{4\mu_f} \notag \end{align} $$
回転抵抗係数$C_\omega$は、粒子の特性、流体の特性、流れ条件などに依存したパラメータであり最終的には 実験的に決める必要がある。
ただし流れが層流（ポテンシャル流）である場合には、$C_\omega$を理論的に導くことができる。
流れ場が層流であると仮定すると、角速度$\omega_p$で回転している球形粒子に作用しているトルク$T$は次式で表される。
$$ \begin{align} T &= 8 \pi \mu \left( \frac{d_p}{2}\right)^3 \omega_p \notag \end{align} $$
これを回転運動の方程式に代入すると、$C_\omega$は次式で表されることがわかる。
$$ \begin{align} C_\omega &= \frac{16\pi}{Re_\omega} \notag \end{align} $$
緩和時間
並進運動、回転運動している粒子は粘性抵抗のためその運動速度が徐々に低減する。
初期速度の$1/e$になるまでに要する時間が緩和時間である。
これは上記の並進運動、回転運動の運動方程式中では$\tau_1,\tau_2$と表現されている。
下表に典型的な緩和時間の例を示す。
ここでは粒子密度$\rho_p=3000[kg/m^3]$、粒子相対速度$u_f-u_p=1[m/s]$、 粒子回転角速度$\omega_p=2\pi[rad/s]$とした。
流体の物性値は後述の計算例で使用した値を使っている。
緩和時間は流体が空気の場合の方が水の場合よりも大きくなるが、 粒径$d$が1[mm]の場合でも高々$3$秒弱である。
これは、微小粒子の回転運動に関して言えば数秒で粒子の回転が減衰してしまうため、 マグナス揚力は作用しなくなることを意味している（粒子に新しく回転を与えない場合）。
並進運動の抵抗係数
並進運動に対する抵抗係数$C_D$はmorsi-alexander式を使用した。
回転運動の抵抗係数
回転運動に対する抵抗係数$C_\omega$は以下のように設定した。
抵抗係数$C_\omega$は回転レイノルズ数が20以下の場合、ポテンシャルから導かれる理論解 $16\pi/Re_\omega$を使用し、回転レイノルズ数が20以上の場合Ansys Fluentで使用している次式を 使用した。
$$ \begin{align} C_\omega &= \frac{16\pi}{Re_\omega} & (Re_\omega \le 20) \notag \\ C_\omega &= \frac{6.45}{\sqrt{Re_\omega}}+\frac{32.1}{Re_\omega} & (Re_\omega>20) \notag \end{align} $$
マグナス揚力係数
マグナス揚力は次式で表される。
$$ \begin{align} F_{RL} &= \frac{1}{8} C_{RL}\pi d^2 \rho |\mathbf{u}-\mathbf{v}|\frac{(\mathbf{u}-\mathbf{v})\times \mathbf{\omega}_p}{|\mathbf{\omega}_p|} \notag \end{align} $$
マグナス揚力係数$C_{RL}$にはいくつかの経験式が提案されている。
まず、回転強度を表すパラメータとして次に示すスピンパラメータ$S$を定義する。
$\omega_p$を粒子回転レイノルズ数$Re_{\omega}$、$\mathbf{u}-\mathbf{v}$を粒子レイノルズ数$Re_p$で表すと、 スピンパラメータは粒子回転レイノルズ数と粒子レイノルズ数の比になっていることがわかる。
スピンパラメータ
$$ \begin{align} S &= \frac{|\omega_p| d_p}{2|\mathbf{u}-\mathbf{v}|} \notag \\ &= \frac{ \frac{4\mu_f Re_\omega}{\rho_f d_p^2} d_p}{2 \frac{\mu_f Re_p}{\rho_fd_p}} \notag \\ &= 2\frac{Re_\omega}{Re_p} \notag \end{align} $$
揚力係数モデル
以下に$C_{RL}$に使われる３種類のモデルを示す。
これらはfluentで採用されているものである。
- OesterleとBui Dinhの式(OB)
$$ \begin{align} C_{RL} &= 0.45+ \left( 4\frac{Re_\omega}{Re_p} -0.45\right) \exp\left( -0.09896 Re_\omega^{0.4} Re_p^{0.3} \right) \notag \\ &= 0.45 + (2S-0.45)\exp\left( -0.075 S^{0.4} Re_p^{0.7} \right) \notag \end{align} $$
- Tsujiらの式(TJ)
$$ \begin{align} C_{RL} &= 0.4 & (for\;S\ge1) \notag \\ C_{RL} &= 0.4S & (for\;S<1) \notag \end{align} $$
- RubinowとKellerの式(RK)
$$ \begin{align} C_{RL} &= 2S \notag \end{align} $$
これらの関係式をグラフに示すと以下のようになる。
OB式については、$S$だけでなく、$Re_p$もパラメータとなるので、 $Re_p=10^{^-2},1,10^2,2\times 10^3$の場合を計算した（カッコ内に$Re_p$の値を示している）。
OB式のグラフを見ると$Re_p$が小さい場合には$C_{RL}=2S$に漸近し、 $Re_p$が大きいと$C_{RL}=0.45$に漸近することがわかる。
解析結果
解析モデル
長さ$0.2$m、高さ$0.1$mの矩形領域を解析領域とした。
左端の高さ$0.05$の位置に固体粒子を設置した。
粒子が移動することで上下の壁面に衝突した場合、 $y$方向速度の符号を逆転して計算を継続するようにした。
設定パラメータ
解析に使用した物性パラメータを示す。流体は水、空気の２種類を想定した。
| property | value | unit | description |
|---|---|---|---|
| $\rho_p$ | $3000$ | [kg/m3] | particle density |
| $\rho_w$ | $1000$ | [kg/m3] | water density |
| $\rho_a$ | $1.2$ | [kg/m3] | air density |
| $\mu_w$ | $10^{-3}$ | [Pa.s] | water viscosity |
| $\mu_a$ | $1.8\times 10^{-5}$ | [Pa.s] | air viscosity |
| $g$ | $9.81$ | [m/s2] | gravity acceleration |
計算ケース一覧
流体速度uf、粒子速度upは一方が有限の値で他方がゼロとなるように指定した。
模式図を示す。なお流体速度ufは全領域で一定とした。
case-1～8では、初期粒子回転速度を$-5\pi,-2\pi,0,2\pi,5\pi[rad/s]$の$5$種類とした。
case-9～12では粒径を$10^{-3},10^{-2},10^{-1}[m]$の３種類とした。
なお粒子追跡は最大時刻$10$秒まで実施した。
以下に各ケースの結果を示す。
最初にアニメーションによる軌跡を示す。
次に、粒子x座標、座標、回転角速度ωの経時変化、最終時刻での粒子軌跡線図を示している。
case-1
このケースでは粒子から見た相対速度($u_{f}-u_{p}$)が$+x$方向のベクトルとなることから、 $\omega_p$が負値の場合上向きの揚力$((u_f-u_p)\times\omega_p$方向）、 正値の場合下向きの揚力が作用することになる。
case-2
このケースはcase-1と同じ条件で流体が空気になったものである。粒子の軌跡線はcase-1とほぼ同じである。
これは揚力係数$C_{rl}$がほぼ同じ（スピンパラメータが同じ）であり、また、 抵抗係数$C_D$も類似した値になるためである。
しかし、水と空気では緩和時間が大きくことなるため、case-1の方が現象が速く進むことがわかる。
case-3
このケースでは粒子から見た相対速度$u_{fluid}-u_{particle}$が$-x$方向のベクトルとなることから、 $\omega_p$が負値の場合下向きの揚力$((uf-up)\times\omega_p$方向）、 正値の場合上向きの揚力が作用することになる。
case-4
このケースはcase-3の流体を空気にしたものであるが、水よりも粘度が小さいため、 粒子に作用する抵抗が小さく、粒子はそのまま下流側に移動する傾向が見られる。
case-5
これはcase-1に対して重力の効果を追加したものである。 粒子回転角速度に依存して軌跡が変化していることがわかる。
case-6
このケースでは粒子はほぼ鉛直に落下し、その後振動しながら徐々に下流側に移動している。 揚力の効果はほとんど見られない。
case-7
これはcase-3に対して重力の効果を追加したものである。 粒子回転角速度に依存して軌跡が変化していることがわかる。
case-8
これはcase-4に対して重力の効果を追加したものであるが、 case-4と同様に、揚力の効果はほとんど見られない。
case-9
case-9～12は比較のための参照ケースとして粒子回転角速度がゼロの状態で、 粒径の違いによる軌跡線への影響を調べたものである。
case-10
case-11
case-12
ソースプログラム
以下に解析に使用したjuliaプログラムを示す(case-5に対応)。
using DifferentialEquations
using Plots
struct Geometry
length::Real
height::Real
end
struct Fluid
name::String
density::Real
viscosity::Real
end
@kwdef struct Ufield{T<:Real}
umax::T
profile::Function
omegaf::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_crl(Rew,Rep)
sp = (Rep<eps() ? 0 : 2*Rew/Rep)
return 2*sp
end
#----------------------------------
function get_cwre(Rew)
if Rew <= 20
CwRe = 16*pi
else
CwRe = 6.45*sqrt(Rew) + 32.1
end
return CwRe
end
#----------------------------------
function orbit!(du,u,p,t)
dp, coefp, coeff, ratio, rhop, rhof, muf, gravity, ufield = p
xp, up, yp, vp, omgp = u
uf, vf = ufield.profile(xp, yp)
udif = sqrt((uf-up)^2 + (vf-vp)^2)
Rep = rhof * dp * udif / muf
cdre = morsi_alexander(Rep)
taup = rhop*dp^2/(18*muf) * 24/cdre
omgf = ufield.omegaf(xp, yp)
Rew = rhof*abs(omgf-omgp)*dp^2/(4*muf)
cwre = get_cwre(Rew)
tauw = rhop*dp^2/(60*muf) * (16*pi/cwre)
gravx, gravy = gravity()
crl = get_crl(Rew,Rep)
(magnx, magny) = ( (vf-vp)*sign(omgp),
-(uf-up)*sign(omgp) )
du[1] = up
du[2] = ( (uf-up)/taup
+ 0.75*crl/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*crl/dp*ratio*udif*magny
# + 9.66*cs/(pi*dp)*sqrt(muf*rhof)/rhop*saffy
+ (1-ratio)*gravy
) / coefp
du[5] = (omgf-omgp)/tauw
end
#----------------------------------
function track(params)
(rhop,fluid,dp,inipos,inivel,iniomg,gravity,ufield,geom) = params
cm = 0.5
hght = geom.height
leng = geom.length
name = fluid.name
rhof = fluid.density
muf = fluid.viscosity
ratio = rhof/rhop
coefp = 1 + ratio*cm
coeff = ratio*(1 + cm)
(xini, yini) = inipos()
(xvel, yvel) = inivel(xini, yini)
omgi = iniomg()
ntraj = length(omgi)
u0_arr = Vector{Vector{Float64}}()
for i in 1:ntraj
push!(u0_arr,[ xini, xvel, yini, yvel, omgi[i] ])
end
tspan = (0.0, 1e1)
function condition(out, u, t, integrator)
out[1] = leng-u[1] # hit outlet
out[2] = hght-u[3] # hit upper wall
out[3] = u[3] # hit lower wall
end
function affect!(integrator, idx)
if idx == 1
terminate!(integrator)
elseif idx == 2
# terminate!(integrator)
integrator.u[4] = -integrator.u[4] # negate y-velocity
elseif idx == 3
integrator.u[4] = -integrator.u[4] # negate y-velocity
end
end
cb = VectorContinuousCallback(condition, affect!, 3)
p = [ dp, coefp, coeff, ratio, rhop, rhof, muf, gravity, ufield ]
prob = ODEProblem(orbit!, u0_arr[1], tspan, p)
function prob_func(prob, i, repeat)
remake(prob, u0 = u0_arr[i] )
end
ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
sim = solve(ensemble_prob, callback=cb, trajectories=ntraj)
lasttime = 0
for i=1:ntraj
(nrow,ncol) = size(sim[i])
lasttime = max(lasttime, sim[i].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)
ops = Vector{Float64}(undef,ncol)
for j=1:ncol
time = sim[i].t[j]
xp = sim[i].u[j][1]
up = sim[i].u[j][2]
yp = sim[i].u[j][3]
vp = sim[i].u[j][4]
op = sim[i].u[j][5]
times[j] = time
xps[j] = xp
ups[j] = up
yps[j] = yp
vps[j] = vp
ops[j] = op
end
out = open( "outfile_f$(name)_d$(dp)."*string(i), "w")
println(out,"time xps ups yps vps ops")
for j=1:ncol
println(out,times[j]," ",xps[j]," ",ups[j]," ",yps[j]," ",vps[j]," ",ops[j])
end
close(out)
end
# animation
nframes = 100
anim = @animate for t in LinRange(0, lasttime, nframes)
local plt
plt = plot( plot_title = "fluid = $(name), d = $(dp) : t = $(round(t, digits = 2))")
for i in 1:ntraj
plot!(sim[i], idxs = (1,3), tspan = (0.0, t), label = "omega=$(round(omgi[i],digits=2))")
scatter!(sim[i], idxs = (1,3), tspan = (t, t), lab = nothing)
end
plot(plt, layout = (1,1),
xaxis = ("x" ,(0, leng)),yaxis = ("y" ,(0, hght)))
end
gif(anim, "trajectory_f$(name)_d$(dp).gif"; fps = 10)
# 結果をプロット
plt = plot(title = "trajectory: fluid = $(name), d = $(dp)", xlabel="x[m]", ylabel="y[m]" )
for i in 1:ntraj
plot!(sim[i],idxs=(1,3), label = "omega0="*string(round(omgi[i],digits=2)),
xlims=(0, leng), ylims=(0, hght) )
end
savefig(plt,"trajectory_f$(name)_d$(dp).png")
plt = plot(title = "history of x: fluid = $(name), d = $(dp)", xlabel="time[sec]", ylabel="x[m]" )
for i in 1:ntraj
plot!(sim[i],idxs=(0,1), label = "omega0="*string(round(omgi[i],digits=2)),
markershape=:circle, markersize=1, plotdensity=100 )
end
savefig(plt,"hist-xp_f$(name)_d$(dp).png")
plt = plot(title = "history of u: fluid = $(name), d = $(dp)", xlabel="time[sec]", ylabel="u[m/s]" )
for i in 1:ntraj
plot!(sim[i],idxs=(0,2), label = "omega0="*string(round(omgi[i],digits=2)),
markershape=:circle, markersize=1, plotdensity=100 )
end
savefig(plt,"hist-up_f$(name)_d$(dp).png")
plt = plot(title = "history of y: fluid = $(name), d = $(dp)", xlabel="time[sec]", ylabel="y[m]" )
for i in 1:ntraj
plot!(sim[i],idxs=(0,3), label = "omega0="*string(round(omgi[i],digits=2)),
markershape=:circle, markersize=1, plotdensity=100 )
end
savefig(plt,"hist-yp_f$(name)_d$(dp).png")
plt = plot(title = "history of v: fluid = $(name), d = $(dp)", xlabel="time[sec]", ylabel="v[m/s]" )
for i in 1:ntraj
plot!(sim[i],idxs=(0,4), label = "omega0="*string(round(omgi[i],digits=2)),
markershape=:circle, markersize=1, plotdensity=100 )
end
savefig(plt,"hist-vp_f$(name)_d$(dp).png")
plt = plot(title = "history of omega: fluid = $(name), d = $(dp)", xlabel="time[sec]", ylabel="omega[rad/s]" )
for i in 1:ntraj
plot!(sim[i],idxs=(0,5), label = "omega0="*string(round(omgi[i],digits=2)),
markershape=:circle, markersize=1, plotdensity=100 )
end
savefig(plt,"hist-op_f$(name)_d$(dp).png")
end
#----------------------------------
function myrun()
#--------------------------------------------------------------
#--------------------------------------------------------------
# 流路形状
leng = 0.2 # 流路長さ
hght = 0.1 # 流路半幅
# 粒子=砂
rhop = 3000
# 流体=水
rhof = 1000; muf = 1e-3; name = "water"
# 流体=空気
# rhof = 1.2; muf = 1.8e-5; name = "air"
# 粒径
dp = 1e-2
# 粒子放出位置
yini = 0.05
# 粒子初期回転角速度[rad/s]
omgi = [ -5pi, -2pi, 0, 2pi, 5pi ]
# 粒子放出速度
uini = 0
# uini = 1.0
# 流速分布タイプ
type = "uniform"
# type = "linear"
# type = "parabolla"
# 最大流速
# umax = 0
umax = 1.0
# 重力加速度
# grav = 0
grav = -9.81
#--------------------------------------------------------------
#--------------------------------------------------------------
if @isdefined yinis
xinis = fill(0.0, length(yinis))
inipos = () -> ( xinis, yinis )
else
inipos = () -> ( 0.0, yini )
end
if uini == nothing
inivel = (xp,yp) -> ( ufield.profile(xp,yp) )
else
inivel = (xp,yp) -> ( uini, 0.0 )
end
iniomg = () -> omgi
gravity = () -> ( 0.0, grav )
geom = Geometry( leng, hght )
fluid = Fluid( name, rhof, muf )
if type == "uniform"
ufield = Ufield(
umax = umax,
profile = (xp,yp) -> (umax, 0.0),
omegaf = (xp,up) -> 0,
)
elseif type == "linear"
ufield = Ufield(
umax = umax,
profile = (xp,yp) -> (umax/hght*yp, 0.0),
omegaf = (xp,up) -> -umax/(2*hght),
)
elseif type == "parabolla"
ufield = Ufield(
umax = umax,
profile = (xp,yp) -> (umax*(1.0 - ((yp-hght)/hght)^2), 0.0 ),
omegaf = (xp,up) -> umax*(yp-hght)/hght^2,
)
end
params = (rhop,fluid,dp,inipos,inivel,iniomg,gravity,ufield,geom)
track(params)
end
myrun()
まとめ
球形粒子にマグナス揚力が作用した場合の効果を計算してみた。
粒子の回転運動が継続する時間は粒子特性、流体特性値から推定される緩和時間に依存しており、 微小粒子の場合、粒子の回転は急速に減衰するためマグナス揚力の効果は、ほとんど生じないことがわかった。
また、マグナス揚力は粒子と流体の相対速度と粒子回転速度の外積で定義される方向に作用することを再現できた。
弊社の解析事例
弊社の流体解析事例については、下記のリンクからご覧ください．
-
[V. Alexiades, A. Solomon, Mathematical Modeling of Melting and Freezing Processes, Hemisphere, Washington, D.C., 1993. ] ↩︎