侧边栏壁纸
博主头像
CyberWEI博主等级

We are the history

  • Wrote 28 Articles
  • Made 9 Tags
  • Got 24 Comments

OutlineCONTENT

Outline

自由落下水流における振動現象のモデリング計算

CyberWEI
2022-08-09 / 1 Comments / 1 Likes / 1,378 Views / 5,138 Words

中文/CN:https://www.cyberwei.com/archives/220810cn

1. 研究背景

図1

蛇口から出る水流を小さく調整して観察すると,その水流は重力の影響を受けて自由落下していることがわかる.水流の断面は円になっており,落下すると同時にその円半径が縮まっている.十分高い場所から落下させると最終的に水流は千切れ,粒状の水滴になる.また,水流が千切れる前に障害物などに当たると,水流は粒状の水滴にならず,直接障害物の表面に沿って流れる.観察してみると,障害物と衝突する直前の水流の部分が小さな波状の凸凹が生じていることがわかる.その凸凹が生成している同時に,水流量を小さく調整すると,凸凹の細いところから水流が千切れ.つまり,この現象は水流が水滴になる時の原因と同じと推定できる.今回の研究では,この凸凹現象の生成についてのモデリング・計算を行う.

2. 現象のモデリング・計算

図1の波動は定常波のように見られるので,液体の落下速度と液体柱の波動速度の和がゼロと考えられる。そのため,液体柱の波動を計算する必要がある.この論文では流体力学の連続方程式及びナビエストークス方程式を使い計算を行う.

2.1. 自由落下水流の半径変化

液体柱振動の計算を行う前にまずは自由落下水流の落下高度と断面半径の関係を計算する.

図2

vdS=const(2.1.1)vdS=const\tag{2.1.1}

p+0.5ρv2+ρgh=const(2.1.2)p+0.5ρv^2+ρgh=const\tag{2.1.2}

水は非圧縮性の連続流体で,ベルヌーイの定理により以下の方程式が得られる.v0v_0vvは液体の初期落下速度と高度hhの時の落下速度で,r0r_0SSは水流の半径と断面積.

水流は自由落下なのでp=0p=0,初期条件を代入して

vdS=vS=const=v0S0(2.1.3)∫vdS=vS=const=v_0 S_0\tag{2.1.3}

0.5ρv2+ρgh=const=0.5ρv02(2.1.4)0.5ρv^2+ρgh=const=0.5ρv_0^2\tag{2.1.4}

連立すると

S=S0v02v02+2gh(2.1.5)S=S_0\sqrt{\frac{v_0^2}{v_0^2+2gh}}\tag{2.1.5}

r=r0v02v02+2gh4(2.1.6)r=r_0\sqrt[4]{\frac{v_0^2}{v_0^2+2gh}}\tag{2.1.6}

落下高度と半径の関係がわかる.

2.2. ナビエストークス方程式を基づいた計算

2.2.1. 連続性方程式

ここでは半径R0R_0,密度ρρ,表面張力係数σσの長い円柱で流れている流体(非粘性)を考える.まずは液体の連続性を使って計算する.

ρt+ρv=0{\frac{∂ρ}{∂t}}+∇ρv=0

ρt=0{\frac{∂ρ}{∂t}}=0

図3

振動を分析するために,この円柱の液体の表面が微小な摂動をしていると考える.ここではαR0α≪R_0である.ωは角速度で,複素数と考えられる.

R~=R0+αei(ωt+kx)(2.2.1)\widetilde{R}=R_0+αe^{i(ωt+kx)}\tag{2.2.1}

urS2+uxS3ux+dxS3=0u_r S_2+u_x S_3-u_{x+dx} S_3=0

urrdθdxur+dr(r+dr)dθdx+(uxux+dx)(12(r+dr)212r2)dθ=0u_rrdθdx-u_{r+dr} (r+dr)dθdx+(u_x-u_{x+dx})(\frac{1}{2}(r+dr)^2-\frac{1}{2}r^2)dθ=0

(12(r+dr)212r2)(\frac{1}{2}(r+dr)^2-\frac{1}{2}r^2)項をテーラー展開で線形化する.

12(r+dr)212r2rdr  (drr0)\frac{1}{2}(r+dr)^2-\frac{1}{2}r^2≈rdr\space\space(\frac{dr}{r}→0)

u(r+dr)urdrr+ur+dr+u(x+dx)uxdxr=0(2.2.2)\frac{u_(r+dr)-u_r}{dr} r+u_{r+dr}+\frac{u_(x+dx)-u_x}{dx} r=0\tag{2.2.2}

ur+druru_{r+dr}-u_rux+dxuxu_{x+dx}-u_xを偏微分として考える.

2.2.2. ナビエストークス方程式

次は非粘性流体のナビエストークス方程式を使って計算する.

ρvt+(ρvv)=p+ρg\frac{∂ρ\textbf{\textit{v}}}{∂t}+∇(ρ\textbf{\textit{v}}\cdot\textbf{\textit{v}})=∇\textbf{\textit{p}}+ρg

今回分析する流体は非圧縮性流体なので,以下の式が得られる.

vt+(vv)=pρ(2.2.3)\frac{∂\textbf{\textit{v}}}{∂t}+∇(\textbf{\textit{v}}\cdot\textbf{\textit{v}})=\frac{∇\textbf{\textit{p}}}{ρ}\tag{2.2.3}

円柱座標でのナビエストークス方程式は以下の形式になる.θθ方向の速度は0.

uθt=0\frac{∂u_θ}{∂t}=0

{urt+ururr+uxurx=1ρpruxt+uruxr+uxuxx=1ρpx(2.2.4)\begin{cases} \frac{∂u_r}{∂t}+u_r\frac{∂u_r}{∂r}+u_x\frac{∂u_r}{∂x} =-\frac{1}{ρ}\frac{∂p}{∂r}\\ \frac{∂u_x}{∂t}+u_r\frac{∂u_x}{∂r}+u_x\frac{∂u_x}{∂x}=-\frac{1}{ρ}\frac{∂p}{∂x}\\ \end{cases}\tag{2.2.4}

上記の式は微小振幅ααの二乗の式が含まれているため,これを0と近似して捨てる.

{urt=1ρpruxt=1ρpx(2.2.5)\begin{cases} \frac{∂u_r}{∂t}=-\frac{1}{ρ}\frac{∂p}{∂r}\\ \frac{∂u_x}{∂t}=-\frac{1}{ρ}\frac{∂p}{∂x}\\ \end{cases}\tag{2.2.5}

ここでは速度と圧力の摂動をRRと同じような形式で表す.

{R~=R0+αei(ωt+kx)ur=R(r)αei(ωt+kx)ux=X(r)αei(ωt+kx)p=P(r)ei(ωt+kx)(2.2.6)\begin{cases} \widetilde{R}=R_0+\alpha e^{i(\omega t+kx)}\\ u_r=R(r)\alpha e^{i(\omega t+kx)}\\ u_x=X(r)\alpha e^{i(\omega t+kx)}\\ p=P(r)e^{i(\omega t+kx)}\\ \end{cases}\tag{2.2.6}

2.2.3. 振動方程式の導出及びその解

式(2.2.6)を(2.2.5)と(2.2.2)に代入する.

{1ρpr=iωR1ρipk=iωXdRdr+Rr+ikX=0(2.2.7)\begin{cases} -\frac{1}{ρ}\frac{∂p}{∂r}=i\omega R\\ -\frac{1}{ρ}ipk=i\omega X\\ \frac{dR}{dr}+\frac{R}{r}+ikX=0\\ \end{cases}\tag{2.2.7}

式(2.2.7)の式2と式3をそれぞれXとRに対して微分する.

{1ρpr=iωR1ρikpr=iωdXdrd2Rdr2Rr2+1rdRdr+ikdXdr=0(2.2.8)\begin{cases} -\frac{1}{ρ}\frac{∂p}{∂r}=i\omega R\\ -\frac{1}{ρ}ik\frac{∂p}{∂r}=i\omega \frac{dX}{dr}\\ \frac{d^2R}{dr^2}-\frac{R}{r^2}+\frac{1}{r}\frac{dR}{dr}+ik\frac{dX}{dr}=0\\ \end{cases}\tag{2.2.8}

上記の式を連立して以下の微分方程式が得られる

r2d2Rdr2+rdRdrR(1+k2r2)=0(2.2.9)r^2\frac{d^2R}{dr^2}+r\frac{dR}{dr}-R(1+k^2 r^2)=0\tag{2.2.9}

この微分方程式は1次の修正ベッセル関数である.方程式の解は第一種修正ベッセル関数で表すことができる.

R(r)=C1I1(kr)(2.2.10)R(r)=C_1 I_1 (kr)\tag{2.2.10}

次はこの方程式の境界条件を考える.

境界条件1:円柱自身の半径の変化率

この時,r=R0r=R_0ur=dRdtu_r=\frac{dR}{dt}が適用される.

ur=R(R0)αei(ωt+kx)=iαωei(ωt+kx)u_r=R(R_0)αe^{i(ωt+kx)}= iαωe^{i(ωt+kx)}

R(R0)=C1I1(kR0)=iαωR(R_0 )=C_1 I_1 (kR_0 )=iαω

C1=iαωI1(kR0)C_1=\frac{iαω}{I_1 (kR_0)}

R(r)=iαωI1(kR0)I1(kr)(2.2.11)R(r)=\frac{iαω}{I_1 (kR_0)}I_1(kr)\tag{2.2.11}

境界条件2:円柱表面の表面張力による応力

図4

P0P_0は振動が発生していないときの円柱の表面の応力で,ppは振動する時の摂動.(図4)

P0+p=σ(1R1+1R2)P_0+p=σ(\frac{1}{R_1}+\frac{1}{R_2})

{P0=σR01R1=cosθR0+αei(ωt+kx)1R0(1αR0ei(ωt+kx))1R2=2Rx2=αk2ei(ωt+kx)\begin{cases} P_0=\frac{σ}{R_0} \\ \frac{1}{R_1} =\frac{cosθ}{R_0+αe^{i(ωt+kx)}}≈\frac{1}{R_0} (1-\frac{α}{R_0} e^{i(ωt+kx)} )\\ \frac{1}{R_2} =-\frac{∂^2R}{∂x^2}=αk^2 e^{i(ωt+kx)} \\ \end{cases}

1R1\frac{1}{R_1} の境界条件は,液体表面を傾きも計算しなくてはならないが,微小振動なのでテーラー展開した結果,その影響は微小量ααの2乗なので1として近似する.上記の式と連立して解くと以下の式が得られる.

p=σ(1R1+1R2)P0=ασR02(k2R021)ei(ωt+kx)p=σ(\frac{1}{R_1} +\frac{1}{R_2})-P_0=\frac{ασ}{R_0^2} (k^2 R_0^2-1) e^{i(ωt+kx)}

次は解いた微分方程式(2.2.11)を式(2.2.7)1ρpr=iωR-\frac{1}{ρ}\frac{∂p}{∂r}=iωRに代入してP(r)P(r)に対して積分する.

P(r)=iiαρω2I1(kR0)I1(kr)drP(r)=-\frac{i\cdot iαρω^2}{I_1 (kR_0)}∫{I_1 (kr)dr}

=αρω2I1(kR0)0rI1(kr)dr=\frac{αρω^2}{I_1 (kR_0)} \int_0^{r}{I_1 (kr)dr}

=αρω2kI0(kr)I1(kR0)(2.2.12)=\frac{αρω^2}{k}\frac{I_0 (kr)}{I_1(kR_0)}\tag{2.2.12}

よって,

p=αρω2kI0(kr)I1(kR0)ei(ωt+kx)(2.2.13)p=\frac{αρω^2}{k}\frac{I_0 (kr)}{I_1(kR_0)}e^{i(ωt+kx)}\tag{2.2.13}

境界条件で得たppの式と微分方程式の解を積分して得たppを連立すると,以下の式が得られる.

ασR02(k2R21)=αρω2kI0(kR0)I1(kR0)\frac{ασ}{R_0^2}(k^2R^2-1)=\frac{αρω^2}{k}\frac{I_0 (kR_0)}{I_1(kR_0)}

ωωR0kR_0 kの関係が得られる.

ω2=σρR03kR0(k2R21)I1(kR0)I0(kR0)(2.2.14)ω^2=\frac{σ}{ρR_0^3}kR_0(k^2R^2-1)\frac{I_1(kR_0)}{I_0(kR_0)}\tag{2.2.14}

液体の下落の速度と円柱液体の振動が上に遡る速度が一致するとき,液体柱の表面の凸凹が定常波のように見える,この時v=ωkv=\frac{ω}{k}の関係が存在する.

v2=ω2k2=σρkR02kR0(k2R21)I1(kR0)I0(kR0)(2.2.15)v^2=\frac{ω^2}{k^2}=\frac{σ}{ρkR_0^2}kR_0(k^2R^2-1)\frac{I_1(kR_0)}{I_0(kR_0)}\tag{2.2.15}

また,微小振動の波長はλ=2πk=2πR0kR0λ=\frac{2π}{k}=\frac{2πR_0}{kR_0}であるので,波長と流れの速度と液体柱の関係は以下の方程式の解となる.

{v2=ω2k2=σρkR02kR0(k2R21)I1(kR0)I0(kR0)λ=2πk=2πR0kR0(2.2.16)\begin{cases} v^2=\frac{ω^2}{k^2}=\frac{σ}{ρkR_0^2}kR_0(k^2R^2-1)\frac{I_1(kR_0)}{I_0(kR_0)}\\ λ=\frac{2π}{k}=\frac{2πR_0}{kR_0}\\ \end{cases}\tag{2.2.16}

2.2.4. 数値計算

次はPythonを使い式(2.2.15)の描画を行う.修正ベッセル関数の値はscipy.specialで計算する.式(2.2.15)で計算したい未知数はkなので,ここではkR0kR_0を変数として関数を描画し,極値を計算する.

f(kR0)=kR0(k2R21)I1(kR0)I0(kR0)f(kR_0)=kR_0(k^2R^2-1)\frac{I_1(kR_0)}{I_0(kR_0)}

Pythonのコードは以下である.

from scipy import special
import numpy as np
import matplotlib.pyplot as plt

x = []
y = []
for i in range(100000):
	xi=i / 100000
    IV1=special.iv(1,xi)
    IV0=special.iv(0,xi)
    y.append(IV1/IV0 * (1 - xi**2) * xi)
    x.append(xi)

plt.plot(x,y)

上記のコードで計算を行い、f(kR0)f(kR_0)のグラフを得ることができる(図5).またこの関数f(kR0)f(kR_0)kR0=0.697kR_0=0.697で最小値を取ることがわかる.

図5

ω2ω^2が負の値になると,R~\widetilde{R}の振動方程式が以下のようになる.

R~=R0+αeβt+ikx\widetilde{R}=R_0+αe^{-βt+ikx}

β=ω2β=|ω^2 |

よって,この時のR~\widetilde{R}の振動が減衰する.つまり,振動は発生しない.減衰が最も早い時はkR0=0.697kR_0=0.697の時である.

次はω2>0ω^2>0の場合を考察する.波長を計算するためには方程式(2.2.17)を解かなければならない.よって次は下記の式を描画を行う.

f(kR0)=kR0(k2R21)I1(kR0)I0(kR0)f(kR_0)=kR_0(k^2R^2-1)\frac{I_1(kR_0)}{I_0(kR_0)}

図6

図6を観察すると,関数曲線が直線y=xy=xに接近していることがわかる.よってkR0kR_0が十分大きい時,f(R0)=kR0f(R_0)=kR_0として近似することが出来る.近似した方程式(2.3.17)の解は以下のようになる.

λ=2πk=2πσρv2(2.2.18)λ=\frac{2π}{k}=\frac{2πσ}{ρv^2}\tag{2.2.18}

3. 実験

図7

図7は高さを一定にして流量だけを変えて撮った実験写真である(高さ5cm).この時,液体が落ちるは一定であるため,流速が一定と考えられる.図を観察すると,液体表面の微小振動の波長は流量・断面半径と関係ないことがわかる.

次は実験数値を代入して式(2.2.18)を検証する.

v2=2ghv^2=2gh

λ=2πσρv2λ=\frac{2πσ}{ρv^2}

よって,λ=πσρgh=0.46mmλ=\frac{πσ}{ρgh}=0.46mmで,測定値と近いと考えられる.

4. おわりに

今回の研究では,自由落下流体柱の振動現象についての物理のモデルを解析し,Pythonを使い数値計算を行った,その結論は以下となる.

解析解:式(2.2.17)近似解:式(2.2.18)

今回は蛇口から出る水流を対象として分析を行ったが,この現象は他の流体でも起こっている.この現象の液体の表面張力と液体自身の運動量が原因となっているもので,表面エネルギーと質量がある流体はすべてこれと似た現象があると考えられる。しかし,今回の計算は液体の粘性が含まれていないため,オイルなどの粘性がある流体の現象を正確に計算するには別の解析が必要となる.

参考文献

[1] Wada, Yoshimasa. On the Steady Surface Ripples of a Cylindrical Flow. Journal of the Physical Society of Japan, 5(4), 259–262, (1950). doi:10.1143/jpsj.5.259

[2] Rayleigh, L. On The Instability Of Jets. Proceedings of the London Mathematical Society, s1 10(1), 4–13, (1878). doi:10.1112/plms/s1-10.

1

Comment