Scilabで霧雨粒の落下運動

Scilabで常微分方程式:摩擦のある運動ではScilabで常微分方程式を解く例として、Scilabを使ってみよう2 - 摩擦のある振動:理工系学生の遊歩道の計算を追試しました。

今回は、もう少し簡単な例として雨粒の落下運動のシミュレーションを行います。具体的にはOctaveで書かれた微分方程式による物理現象のモデル化(PDF)の例題2(雨粒の落下運動)をScilabで書き直します。


001_20130702073402.jpg

Fig.1: 雨粒は、重力による加速と空気抵抗による減速がつりあって地表付近ではほぼ一定の速度で落下してくる。画像はフリー画像素材を使わせていただいています。(c) Jonathan Kos-Read


ただし元PDFでは、雨粒の重さが1kgもあったりと数値が余りにも非現実的なので雨粒の落下速度を参考に半径0.01mmの霧雨粒の運動をシミュレーションします。



速度に比例する抵抗


速度に比例する抵抗が働く物体の運動方程式は

m \frac{\mathrm{d}v}{\mathrm{d}t} = mg - kv

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

となります。この両辺を質量mで割って

\frac{\mathrm{d}v}{\mathrm{d}t} = g - \frac{k}{m}v<br />

と dv/dt=... の形にすればScilabのode関数で解けるようになります。(参考:Scilabで常微分方程式:摩擦のある運動)

またこの式の解は、解析的な解が知られていて以下のように書き表せます。

v(t) = \frac{mg}{k}\left\{1 - \exp\left(- \frac{k}{m}t\right) \right\}

比例係数kの見積もり


運動方程式の中の空気抵抗の係数kを雨滴の落下速度に倣って決定します。

前述の解析解を t → ∞ とすると雨粒の速度が一定になった後の速度(終端速度)が得られます。

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


またストークスの式から流体中を落下する物体の終端速度は、流体と物体の間の密度差や流体の粘性率などから

v_{\infty} = \frac{2 r^2 \Delta \rho g}{9 \mu}

(r: 雨粒の半径, Δρ: 水と空気の密度差, g: 重力加速度, μ: 空気の分子粘性係数)

空気の密度は、水の密度に比べて無視できるほど小さいので空気と水の密度差は実質的に水の密度と考えることができます。(参考:水・空気の物性)

\Delta \rho \simeq \rho_{\mathrm{H_2 O}}

解析解から得られた終端速度とストークスの式から得られた終端速度は等しいはずなので

\frac{mg}{k} = \frac{2 r^2 \rho_{\mathrm{H_2 O}} g}{9 \mu}

ここで、水滴の重さは

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

したがって

k = 6 \pi \mu r

と求まります。

数値シミュレーション


入力する物理定数はTable.1のようにしました。
雨滴の落下速度では、水の密度が103 (kg/m-3)と書かれていますが、これは103(kg/m3)の誤植でしょう。

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


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

clear;

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

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

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

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

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

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


計算結果


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


002_20130702082152.png

Fig.2: 霧雨粒の落下運動。横軸が時間、縦軸が速度。赤丸で示したのが数値計算の結果で、緑破線が解析解の結果。半径0.01mm程度の雨粒は落下を開始してから0.01秒以内にほとんど終端速度へ達する。


雨粒の重さmや摩擦係数kといったパラメータの違いが有る為、横軸・縦軸のスケールは変わっていますが元PDFの結果と相似形のグラフが描けました。
またode関数を使った数値計算の結果は、解析解をよく再現することが分かります。

秒速1センチなんだって。半径0.01mmの霧雨粒の落ちるスピード。

関連エントリ




参考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シンボル軸ラベル凡例線種トラックボール

最新コメント
リンク

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