Scilabで大きな雨粒の落下運動

Scilabで霧雨粒の落下運動では、雨滴の落下速度を参考に微分方程式による物理現象のモデル化(PDF)の例題2の再現をScilabで行いました。

今回は、同じPDFの問題16にあたる速度の2乗に比例する空気抵抗のある場合の計算を行います。

001_20130704123801.jpg

Fig.1: 落下する雨粒の受ける空気抵抗は、雨粒の大きさによって式の形が変わる。粒径の大きな雨粒は、速度の2乗に比例した空気抵抗を受けて落下する。写真は (c) 写真素材 素材三昧



速度の2乗に比例する抵抗


Scilabで霧雨粒の落下運動では、粒径の小さい0.01mmの霧雨粒の落下運動のシミュレーションをScilabを用いて行いました。
粒径の小さい霧雨粒は、落下速度に比例する空気抵抗が働くと考えましたが、粒径が大きくなると速度の2乗に比例する空気抵抗が働くことになります。

前回と同様にして、速度の2乗に比例する運動方程式を立てると

m \frac{\mathrm{d}v}{\mathrm{d}t} = mg - k v^2

(m: 雨粒の質量, v: 雨粒の落下速度, t: 時間, g: 重力加速度, k: 比例定数)

となります。Scilabのode関数を利用するためには dv/dt=... の形にする必要があるので両辺を質量mで割って

\frac{\mathrm{d}v}{\mathrm{d}t} = g - \frac{k}{m}v^2

とします。

この式の解析解も既に知られていて

v(t) = \sqrt{\frac{mg}{k}} \tanh \left(t \sqrt{\frac{kg}{m}}\right)

です。

比例定数kの見積もり


解析解から得られる終端速度は

v_{\infty} = \lim_{t \to \infty} v(t) = \sqrt{\frac{mg}{k}}

です。

前回はストークスの式から終端速度を求めましたが、(後述する)レイノルズ数Reの値によってこの式は形を変えます。(参考:終端速度)

(i) Re < 2のとき

v_{\infty} = \frac{2 r^2 (\rho_{\mathrm{H_2 O}- \rho_{\mathrm{air}}}) g}{9 \mu}


(ii) 2 < Re < 500のとき

v_{\infty} = 2 r \left\{\frac{4}{255} \frac{(\rho_{\mathrm{H_2 O}} - \rho_{\mathrm{air}})^2 g^2}{\rho_{\mathrm{air}}\mu} \right\}^{\frac{1}{3}}

(iii) 500 < Re < 105のとき

v_{\infty} = \left\{\frac{8r}{3 \times 0.44}\frac{(\rho_{\mathrm{H_2 O}} - \rho_{\mathrm{air}})g}{\rho_{\mathrm{air}}} \right\}^{\frac{1}{2}}

(i)のときがストークスの式です。今回は雨粒の半径が1mm程度の計算を行うために一番下の(iii)の式を利用します。雨粒の大きさとレイノルズ数の関係はAppendixを見てください。

これらの終端速度が一致するはずなので

\sqrt{\frac{mg}{k}} = \left\{\frac{8r}{3 \times 0.44}\frac{(\rho_{\mathrm{H_2 O}} - \rho_{\mathrm{air}})g}{\rho_{\mathrm{air}}} \right\}^{\frac{1}{2}}

更に

\rho_{\mathrm{H_2 O}} - \rho_{\mathrm{air}} \simeq \rho_{\mathrm{H_2 O}}
m = \rho_{\mathrm{H_2 O}} \times \frac{4}{3} \pi r^3

を用いて整理すると

k = 0.22 \pi \rho_{\mathrm{air}} r^2

となります。雨滴の落下速度では

k = 0.235 \pi \frac{\mu}{\nu} r^2

とかかれており、これはν= μ/ρを用いて(レイノルズ数)

k = 0.235 \pi \rho_{\mathrm{air}} r^2

となり(先頭の係数がちょっとちがいますが)同じ形になります。

数値シミュレーション


説明記号
水の密度ρH2O1.0 × 106 (g/m3)
空気の密度ρair1.3 × 103(g/m3)
空気の分子粘性係数μ1.8 × 10-2 (g s/m)
重力加速度g9.8 (m/s2)
table.1: 物理定数


これを踏まえたScilabソースコードを以下に示します。

clear;

// *** 入力パラメータ ***
r = 1E-3; // 雨粒の半径 (m)

// 物理定数
rho_h2o = 1.0E6; // 水の密度 (g/m^3)
rho_air = 1.3E3; // 空気の密度 (g/m^3)
mu = 1.8E-2; // 空気の分子粘性係数 (g s/m)
//nu = mu / rho_air; // 空気の動粘性係数
g = 9.8; // 重力加速度 (m/s^2)

// *** 運動方程式 ***
m = rho_h2o * 4 * %pi * r ^ 3 / 3; // 雨粒の重さ (g)
//k = 0.235 * mu * %pi * r * r / nu; // 空気抵抗の係数 (g/s)
k = 0.22 * %pi * rho_air * r ^ 2; // 空気抵抗の係数 (g/s)
v0 = 0; // 初期速度 (m/s)
// 解くべき方程式の定義
deff("vdot = df(t,v)", "vdot = g - (k / m) * v ^ 2");

// *** 常微分方程式の計算 ***
// 時間ベクトル(s)
T = linspace(0, 5, 101);
// 常微分方程式の数値解
V = ode(v0, 0, T, df);

// *** 速度の描画 ***
// 数値計算(赤丸)
plot(T, V, 'or');
// 解析解(緑破線)
plot(T, sqrt(m * g / k) * tanh(sqrt(k * g / m) .* T), '--g');

// *** グラフの体裁 ***
// 軸ラベル
xlabel("Time (s)");
ylabel("Velocity (m/s)")
// 凡例
legend(["Numerical";"Analytical"],4)


計算結果


以下にScilabによる計算結果を示します。

002_20130704124843.png

Fig.2: 大きな雨粒の落下運動。横軸が時間、縦軸が速度。赤丸で示したのが数値計算の結果で、緑破線が解析解の結果。半径1mm程度の雨粒は落下を開始してから数秒ほどで終端速度へ達し、その速度は毎秒7メートル程度である。


Appendix: レイノルズ数


終端速度の場合わけは無次元パラメータであるレイノルズ数Reによって決まります。

Re = \frac{\rho_{\mathrm{air}}vL}{\mu}

air: 空気の密度, v: 平均速度, μ: 空気の粘性係数, L: 特性長さ)

特性長さというのがどのように決めるべき値なのか難しいのですが、レイノルズ数によると流れの中の球に関し、特性長さは球の直径であり、特性速度は、球からすこし離れた場所の球の運動が流体の検査体を乱さないような場所にある流体と球との相対速度である。との事なので雨粒の直径を用い、平均速度は終端速度としました。

これをまとめると

Re = \frac{2r \rho_{\mathrm{air}}v_{\infty}}{\mu}

となります。

Scilabを用いて雨粒の半径と(i)-(iii)それぞれの式から計算した終端速度、及びレイノルズ数を計算しました。

003_20130704124843.png

Fig.3: 雨粒の半径と終端速度(上)、雨粒の半径とレイノルズ数(下)。場合わけの式により(i)が赤、(ii)が緑、(iii)が青でプロットしてある。レイノルズ数のグラフにはどの式を使うかのしきい値となるRe = 2, 500, 105が比較用として破線でプロットしてある。


clear;

// *** 物理定数 ***
rho_h2o = 1.0E6; // 水の密度 (g/m^3)
rho_air = 1.3E3; // 空気の密度 (g/m^3)
mu = 1.8E-2; // 空気の分子粘性係数 (g s/m)
g = 9.8; // 重力加速度 (m/s^2)

// *** 終端速度の関数 ***
// Re < 1
function v = v1(r)
v = 2 .* r .^ 2 * (rho_h2o - rho_air) .* g ./ 9 ./ mu
endfunction
// 2 < Re < 500
function v = v2(r)
v = 2 * r * (4 * (rho_h2o - rho_air) ^ 2 * g ^ 2 / 255 / rho_air / mu) ^ (1 / 3)
endfunction
// 500 < Re < 10E5
function v = v3(r)
v = sqrt(8 * r * (rho_h2o - rho_air) * g / (3 * 0.44 * rho_air))
endfunction

// *** 雨粒の半径 ***
R1 = logspace(-7,-3);
R2 = logspace(-6,-2);
R3 = logspace(-5,-1);
R = logspace(-7,-1);

// *** 終端速度のプロット ***
xsetech([0,0,1,0.45],[1E-7,1E-6,1E-1,1E3],"ll")
plot(R1,v1(R1),"-r");
plot(R2,v2(R2),"-g");
plot(R3,v3(R3),"-b");
xlabel("Radius (m)")
ylabel("Terminate velocity (m/s)")

// *** レイノルズ数のプロット ***
xsetech([0,0.5,1,0.45],[1E-7,1,1E-1,1E7],"ll")
plot(R,2,"--k");
plot(R,500,"--k");
plot(R,1E5,"--k");
plot(R1,2 .* R1 .* rho_air .* v1(R1) ./ mu,"-r");
plot(R2,2 .* R2 .* rho_air .* v2(R2) ./ mu,"-g");
plot(R3,2 .* R3 .* rho_air .* v3(R3) ./ mu,"-b");
xlabel("Radius (m)")
ylabel("Reynolds number"


関連エントリ




参考URL




付録


このエントリで使用したScilab用ソースコードを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。


tag: Scilab 常微分方程式 ode 

comment

Secret

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAPICOPアンプecalj常微分方程式モンテカルロ解析状態密度トランジスタodeDOSインターフェース定電流スイッチング回路PDS5022半導体シェルスクリプト分散関係レベルシフト乱数HP6632AR6452A可変抵抗トランジスタ技術温度解析ブレッドボードI2C反強磁性確率論数値積分セミナーバンドギャップバンド構造偏微分方程式非線形方程式ソルバ熱設計絶縁ISO-I2Cカオス三端子レギュレータLM358GW近似マフィンティン半径A/DコンバータフォトカプラシュミットトリガLEDPC817C発振回路数値微分直流動作点解析サーボカレントミラーTL431アナログスイッチUSB74HC4053bzqltyVESTA補間電子負荷アセンブライジング模型BSch量子力学単振り子2ちゃんねるチョッパアンプLDA開発環境基本並進ベクトルFFT標準ロジックブラべ格子パラメトリック解析抵抗SMPMaxima失敗談ラプラス方程式繰り返し位相図スイッチト・キャパシタ熱伝導状態方程式キュリー温度gfortranコバルトTLP621不規則合金Quantum_ESPRESSO六方最密充填構造ランダムウォーク相対論ewidthスピン軌道相互作用FETQSGWVCAcygwinスレーターポーリング曲線GGA仮想結晶近似PWscfシュレディンガー方程式LM555ハーフメタル固有値問題NE555最小値ガイガー管QNAPUPS自動計測ダイヤモンドマントルTLP552格子比熱最適化MCU井戸型ポテンシャル最大値xcrysdenCIF条件分岐詰め回路フェルミ面差し込みグラフスーパーセルfsolveブラウン運動awk過渡解析起電力三角波第一原理計算FXA-7020ZRWriter509Ubuntuテスタ熱力学データロガーTLP521OpenMPubuntu平均場近似MAS830LトランスCK1026PIC16F785PGA2SC1815EAGLEノコギリ波負帰還安定性ナイキスト線図MBEOPA2277P-10フィルタCapSenseAACircuitLMC662文字列固定スピンモーメントFSMTeX結晶磁気異方性全エネルギーc/a合金multiplotgnuplot非線型方程式ソルバL10構造正規分布等高線ジバニャン方程式初期値interp1fcc面心立方構造ウィグナーザイツ胞半金属デバイ模型電荷密度重積分SIC二相共存磁気モーメント不純物問題PWgui擬ポテンシャルゼーベック係数ZnOウルツ鉱構造edeltquantumESPRESSOフォノンリジッドバンド模型スワップ領域BaO岩塩構造ルチル構造ヒストグラム確率論マテリアルデザインフラクタルマンデルブロ集合キーボードRealforceクーロン散乱三次元疎行列縮退化学反応関数フィッティング最小二乗法Excel直流解析PCTS-110TS-112日本語パラメータ・モデル等価回路モデルcif2cell入出力陰解法熱拡散方程式HiLAPW両対数グラフCrank-Nicolson法連立一次方程式specx.fifort境界条件片対数グラフグラフの分割円周率ヒストグラム不規則局所モーメントGimpシンボル軸ラベル凡例線種トラックボール

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ