Octaveの精義

ねがてぃぶろぐでは、数値計算にScilabを利用しています。
数値計算といえばMATLABがもっとも有名ですが、高価です。ScilabとOctaveはどちらもMATLABに良く似た数値演算ソフトで、無償で利用することが出来ます。

ScilabにせよOctaveにせよ、超有名なMATLABのクローンであるという立ち位置なので勉強するための書籍の数もMATLABに劣ります。特に初心者のうちは、サンプルとして紹介されているものをコピーして使うのが上達の近道であると私は考えているので、参考に出来るものの量が少ないことは深刻な問題です。

そんなわけで今回は、Octaveの参考書籍ではありますがOctaveの精義を紹介します。




常微分方程式


Scilabのプログラミングをするにあたって、サンプルコードを見つけにくいのが常微分方程式です。
常微分方程式の解法の解説サイトとして、私が一番読みやすいと感じるのはSCILABで微分方程式を解く。です。

しかしながら、常微分方程式といえば、現実的な問題に対する数学モデルです。そういった意味で、SCILABで微分方程式を解く。の説明は、xとかyとか抽象的なことを言われても良くわからんと思ってしまいます(最初は)。

そこで、少しでも計算のイメージが沸くであろう基本的な力学の計算をScilabで霧雨粒の落下運動Scilabで大きな雨粒の落下運動で行いました。これらは微分方程式による物理現象のモデル化(PDF)で紹介されている事例の追試です。


002_20130702082152.png

Fig.2: 霧雨粒の落下運動のシミュレーション。微分方程式による物理現象のモデル化(PDF)で紹介されているスクリプトのScilab移植版。


実を言うと今回紹介するOctaveの精義の著者は(おそらく)微分方程式による物理現象のモデル化(PDF)の著者と同じ方です。いくつか共通するソースコードが見つかりました。

Octaveの精義のはじめにでは

著者はOctaveを主として常微分方程式の計算とその結果のグラフ化に使っている


と書いています。また

原則的にはスクリプトは必ず示すことにする.なぜなら「百聞は一見に如かず」のとおり,私の経験では,スクリプトの現物が役立つことが多かったからです.


とも有ります。

実際に内容はその通りで、基本的な説明や様々な分野のための数値計算の教科書として一通りの応用事例も含んでいるものの、常微分方程式の説明に詳しく、サンプルのソースコードが常に書かれています。

ねがてぃぶろぐでは、今後も微分方程式の数値解をScilabを利用して計算していくつもりですが、手っ取り早く常微分方程式を解くコードが欲しいと考えるなら、Octaveの精義を参考にOctaveのスクリプトを書くほうが簡単かもしれません。

蛇足


共通するスクリプトを見て著者が同じことに気が付いたといえば、Scilab入門書籍2冊で紹介したScilab入門―電気電子工学で学ぶ数値計算ツールは、Scilab つかいませんかと同じ方が書いているようです。ネットは狭いですね。



蛇足ついでにもうひとつ。
現実の物理現象を微分方程式を含んだ数学モデルにするという観点から、微分方程式で数学モデルを作ろう道具としての微分方程式―「みようみまね」で使ってみよう (ブルーバックス)も読みたいと思っているのですが、なかなか時間が取れません。いつかレビューしたいと思うのですが。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 常微分方程式 

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 

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  

Scilabで常微分方程式:摩擦のある運動

今回は理工系学生の遊歩道で公開されているScilabを使ってみよう2 - 摩擦のある振動のプログラムを走らせます。自分ではプログラムしません。
したがって、内容は全く同じなのですが、元記事では数式をまとめたPDFがリンク切れしているので、その点だけ補完します。


001_20130604224856.png

Fig.1: 摩擦のある床の上でのバネにつながれた物体の運動


摩擦が働く床上に物体があり、ばねで壁につながれています。摩擦の大きさは速さに比例すると考えましょう。このとき、物体はどのように振動するのかを考えてみましょう。



常微分方程式の解法


常微分方程式の解は、常微分方程式ソルバ(ordinary differential equation solver; odeソルバ)を使ってとくことができます。odeソルバの中身はルンゲ・クッタ法(runge-kutta)という有名な方法です。

Scilabでは、これがode関数として実装されています。実際の問題から、odeを使うために私たちがしなければならないことは、解くべき方程式を

\begin{pmatrix}y^{'}\\y^{''}\\\vdots\\y^{(n)}\end{pmatrix}


の形に変形するだけです。

摩擦ある床の上でのバネにつながれた物体の運動


今回は理工系学生の遊歩道で公開されているScilabを使ってみよう2 - 摩擦のある振動の問題の通り、摩擦のある床の上に置かれた、バネにつながれた物体の運動(位置と速度)が時間とともにどのように変化するのかを計算します。

物体に働くバネの力は、バネ定数をkとおくとフックの法則より

F_k = -kx


また、摩擦力は単純に速度に比例するとして、その比例定数をRとおくと

F_R = -Rv

となります。

これらを F = ma に代入すると、運動方程式は

-kx-Rv = ma


となります。

さてここから冒頭の行列の形へ変形します。
運動方程式は、位置xと時間tの関数なので、目標は以下の形に変形することです。

texclip20130525031606.png


位置の時間微分は速度なので

\dot{x}=v


位置の二回微分、つまり速度の微分は加速度なので、運動方程式から

\ddot{x}=\dot{v}=a=-\frac{k}{m}x-\frac{R}{m}\dot{x}


従って求める式は

texclip20130525032106.png


つまり

texclip20130525032349.png


です。

計算の結果に関して、位置と速度をプロットしました。


002_20130605114941.png

Fig.2: 計算結果、上が位置、下が速度の時間変化。


関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 常微分方程式 ode  

第23回CMD®ワークショップ参加報告

9月2日(月)から6日(金)にかけて大阪大学で開催された第23回コンピュテーショナル・マテリアルズ・デザイン(CMD®)ワークショップに参加させていただきました。


001_201309080457131b0.jpg

Fig.1: 大阪大学ナノサイエンスデザイン教育研究センター



私が受講したのは、最も初歩的なビギナーコースで、以下の3つのコードの使い方を教えていただきました。


私の一番の目的はAkaiKKR(Machikaneyama)の使い方を教えていただくことでした。CygwinでAkaiKKR(Machikaneyama)などでいくつかの予習をして受講したので、これまで疑問だった点の質問をすることも出来ました。

ビギナーコースの参加者は、大阪大学やそれ以外の大学の学部生から企業で働いている社会人の方までさまざまでした。

特に社会人の方は、むしろ計算を本職にしている人は少ないように感じました。
自分は実験屋さんで、計算屋さんと話をする機会があるのだけれど、彼らとのコミュニケーションをスムーズにする
ために第一原理計算の背景を知っておきたい、こういう人たちが多かった印象です。

実際にコンピュテーショナル・マテリアルズ・デザイン(CMD®)ワークショップに参加するに当たって気になるのは、参加者の知識・技術レベルはどの程度のものが要求されるのかということだと思います。公式のFAQには以下の様にあります。

どのコースを選択したら良いのか?基準は?

ビギナーコースは大学院生の実習でもあることから、学部卒の基本的な知識があれば、UNIXの経験はなくても参加は可能です。コースの中にUNIX講座を設けています。ただし、もちろん自身である程度習得されているほうが望ましいです。


学部卒の基本的な知識というのは、恐らく量子力学と固体物理の基本的な部分だと思います。
UNIXの知識は無くても参加できるとあり、実際にコマンドの使い方が分からない受講生の方がチューターに逐一質問するという光景は頻繁に見られました。

しかしながら、コンピュータの知識も量子力学や固体物理の知識も両方とも自信がないとなると、受講は二重苦になってしまうため、あらかじめコマンドラインを触ってみるという事をやっておいた方が良いと思います。
CygwinでAkaiKKR(Machikaneyama)に書いたとおりCygwinを入れておけばAkaiKKR(Machikaneyama)の計算も自分で出来るのでオススメです。

実際に使うコマンドは多くなく
  • cd
  • ls
  • cp
  • mkdir
とかの一般に良く使われるやつです。

実習で使うPCは、大阪大学のものを貸していただくことも、自分でノートPCを持ち込むことも出来ました。
とはいえ実際に計算を行うのは大阪大学のサーバーなので、PCに要求される機能はコマンドラインからsshでサーバーに接続することが出来るという点だけなのだと思います。私は自信がなかったので当日は大阪大学のものをお借りしましたが、家に帰ってからCygwinのsshで大阪大学のサーバーへ接続してみたところ、何の問題も無く成功しました。

コンピュテーショナル・マテリアルズ・デザイン(CMD®)ワークショップへの参加は、私にとって非常に有意義なものでした。習った手法を生かした計算は、また別のエントリに書こうと思います。

次回の開催は以下の通りです。

日時: 平成26年2月24日(月)~28日(金)
場所: 理化学研究所計算科学研究機構(神戸市中央区港島町7-1-26)

まだウエブページ上での告知はなされていませんが、申し込みの締め切りは平成26年1月19日との事です。

このエントリを読んで興味を持った方は、是非、参加をすることをオススメします。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: CMD 計算機マテリアルデザイン AkaiKKR Machikaneyama STATE-Senri ABCAP 第一原理計算 

AkaiKKRで鉄の安定相と格子定数

AkaiKKRでダイヤモンド型構造半導体AkaiKKRでニッケル・鉄・コバルトでは、色々な結晶構造を持つ固体の電子構造の計算を行いましたが、その格子定数はあらかじめ分かっているものとして入力ファイルを作成しました。

今回は、全エネルギー最小化という条件から、結晶構造や格子定数を第一原理的に決定することを目的にα, γ, δ-Feの計算をAkaiKKRで行い、計算結果をMurnaghanの状態方程式へフィッティングしました。

その結果α-Feが最も安定で、その格子定数は a=2.81(Å) となる事がわかりました。これは、実測値である a=2.87 (Å) と極めて近い値です。


第一原理計算と格子定数


『第一原理計算とは』といった説明は、調べてみると色々な方が説明されていますが、私にとって直感的に分かりやすかったのは第一原理計算 (基礎編)の下記のものです。

第一原理(英語はfirst-principlesまたはab initio(英語ではない))とは、なんら実験データや経験パラメーターを使わないで理論計算をする方法の総称。でも、この業界では電子状態計算、平たく言うとシュレディンガー方程式を解く計算のことを暗に指すことが多い。


第一原理計算というために大切なことは、実験データを見てパラメータをいじるということをしていないということです。

とは言うものの実際に第一原理計算パッケージを使うときには、これまでAkaiKKRを使ったエントリにも当てはまることですが、結晶構造やその格子定数は計算の前に入力する必要があります。

例えば鉄の場合について考えます。

鉄は常温では体心立方(body-centered cubic; bcc)構造となり(フェライト; α-Fe)強磁性を持ちます。温度を上げていくと911℃で面心立方(face-centered cubic; fcc)構造に相転移します(オーステナイト; γ-Fe)。さらに温度を上げていくと1392℃で再びbcc構造へ戻りますが(デルタフェライト; δ-Fe)、このときは既にキュリー温度を超えているため常磁性です。

AkaiKKRの入力ファイルを作成するとき、常温(というよりは絶対零度)でα, γ, δのどの鉄が安定なのか、また、その格子定数がいくつなのかは実験データを見ないと分からないことになります。

第一原理計算が実験データを使わないとは言っても、入力する結晶構造と格子定数だけは多くの場合大目に見てくれるというのが暗黙の了解の様ですが、厳しい見方をするなら第一原理でないということも出来ます。

この問題は、いろいろな結晶構造でいろいろな格子定数を使って計算をした後、全エネルギーを比較することによって解決することが出来ます。
実例として、アドバンスソフト社が発売しているAdvance/PHASEで行われた計算例が公開されています。

basis_spin_01.png
Fig.1: Advance/PHASEによるFeの全エネルギーと体積の関係


今回のエントリではこの計算をAkaiKKRで行います。

specx.fの設定


今回の計算では、specx.fを下記の設定にしてmakeしました。
     & (natmmx=4, ncmpmx=4, msizmx=198, mxlmx=3, nk1x=2200, nk3x=2688,
& msex=201, ngmx=15, nrpmx=650, ngpmx=650, npmx=350, msr=400)


入力ファイル


磁性を考慮した(magtyp=mag)bcc鉄(alpha_in.txt)、及び考慮しない(magtyp=nmag)bcc鉄(delta_in.txt)とfcc鉄(gamma_in.txt)の3つの入力ファイルを用意し、それぞれ1原子あたりの体積が8-15(Å3)となるようにしました。

Murnaghanの状態方程式へのフィッティング


計算結果の体積とエネルギーの対応関係(Energy_dat.txt)から、最もエネルギーが低いときの体積を求めるため、gnuplotを用いてMurnaghanの状態方程式へ最小二乗法フィッティングを行います。

Murnaghanの状態方程式は下記の式で表されます。

E_{tot}(V) = \frac{BV}{B^{'}(B^{'}-1)}\left[B^{'}\left(1-\frac{V_0}{V}\right)+\left(\frac{V_0}{V}\right)^{B^{'}}-1\right]+E_0

ここでEtotが全エネルギー、Bが体積弾性率、B'が体積弾性率の圧力微分、Vが体積です。
フィッティングとグラフの描画を行うgnuplotのスクリプトはEnergy_plt.txtです。

結果


結果をFig.2に示します。Advance/PHASEで行われた計算結果と比較すると(定量的なものはともかく)強磁性bcc構造のα-Feが1原子辺りの体積が11(Å3)のあたりで最もエネルギーが低く安定であることが再現できています。


001_20130309181913.png
Fig.2: AkaiKKRによるFeの全エネルギーと体積の関係


実際に測定された鉄の格子定数、Advance/PHASEによる計算結果と今回AkaiKKRで計算したものをTable.1にまとめました。

実測PHASEAkaiKKR
a (Å)2.872.8452.81
B (GPa)168177.237190
table.1: 強磁性bcc鉄の格子定数と体積弾性率


実測とAdvance/PHASEの値は第一原理シミュレータ入門 PHASE&CIAOから引用しました。

ひょっとしたらたまたまなのかも知れませんが、格子定数に関しては非常によく一致しています。

Appendix: β-Fe?


今回のエントリでは「強磁性bcc鉄=α-Fe」「常磁性bcc鉄=δ-Fe」という様な書き方をしてしまいましたが、これは厳密には間違いです。鉄のキュリー温度は770℃なので、常温付近では常磁性だったbcc鉄はfcc(γ-Fe)へ相転移する911℃よりも低温で常磁性になります。
この770-911℃の間で安定な常磁性bcc鉄のことを以前はβ-Feと呼んでいたようですが、現在ではα-Feの一部として扱われています。

強磁性bcc鉄(α-Fe)
 |
 | 770℃
 ↓
常磁性bcc鉄(これもα-Fe,以前はβ-Feと呼ばれていた。)
 |
 | 911℃
 ↓
fcc鉄(γ-Fe)
 |
 | 1392℃
 ↓
常磁性bcc鉄(δ-Fe)

AkaiKKRでは絶対零度の計算しか出来ないため『もしも絶対零度のとき常磁性bcc構造が安定だったとしたらどうなるか』という計算を(本来高温の相である)δ-Feと呼ぶのはおかしいのですが、表記をややこしくしないためこう書きました。

なお、第一原理シミュレータ入門 PHASE&CIAOの表3.4にはα,β,γと書いてありますが、これは単なるα,γ,δの誤植だと思われます。

関連エントリ




参考URL




付録


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


参考文献/使用機器





フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR 安定相 状態方程式 

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

LTspiceAkaiKKRScilabmachikaneyamaKKRPSoCOPアンプPICCPA強磁性モンテカルロ解析常微分方程式トランジスタode状態密度インターフェースDOSPDS5022ecaljスイッチング回路定電流半導体シェルスクリプトレベルシフト乱数HP6632A温度解析ブレッドボードI2CR6452A分散関係トランジスタ技術可変抵抗確率論数値積分反強磁性セミナー非線形方程式ソルバ絶縁バンドギャップ熱設計偏微分方程式バンド構造GW近似カオス三端子レギュレータLEDフォトカプラシュミットトリガISO-I2CA/DコンバータLM358USBカレントミラーTL431マフィンティン半径PC817C数値微分アナログスイッチ発振回路サーボ直流動作点解析74HC40532ちゃんねる標準ロジックチョッパアンプLDAアセンブラFFTbzqltyイジング模型ブラべ格子開発環境補間量子力学電子負荷BSchパラメトリック解析単振り子基本並進ベクトル熱伝導繰り返しGGAMaximaTLP621ewidthSMP相対論抵抗位相図ランダムウォークスピン軌道相互作用六方最密充填構造不規則合金FETコバルト失敗談QSGWcygwinスレーターポーリング曲線スイッチト・キャパシタラプラス方程式gfortranキュリー温度状態方程式条件分岐格子比熱TLP552LM555TLP521三角波NE555過渡解析FXA-7020ZRWriter509テスタ詰め回路MCUマントルダイヤモンドQNAPデータロガーガイガー管自動計測UPS井戸型ポテンシャルawk第一原理計算仮想結晶近似ブラウン運動差し込みグラフ平均場近似fsolve起電力熱力学OpenMPスーパーセル固有値問題最適化最小値VCAシュレディンガー方程式VESTAubuntu最大値面心立方構造PGAOPA2277L10構造非線型方程式ソルバ2SC1815fccフェルミ面等高線ジバニャン方程式ヒストグラム確率論マテリアルデザイン正規分布結晶磁気異方性interp1フィルタ初期値ウィグナーザイツ胞c/aルチル構造岩塩構造スワップ領域リジッドバンド模型edeltBaOウルツ鉱構造重積分SIC二相共存ZnOquantumESPRESSOCapSensegnuplotmultiplot全エネルギー固定スピンモーメントFSM合金ノコギリ波フォノンデバイ模型ハーフメタル半金属TeXifortTS-110不規則局所モーメントTS-112等価回路モデルパラメータ・モデルヒストグラムExcel円周率GimpトラックボールPC直流解析入出力文字列マンデルブロ集合キーボードフラクタル化学反応三次元Realforce縮退日本語最小二乗法関数フィッティング疎行列シンボル線種ナイキスト線図陰解法負帰還安定性熱拡散方程式EAGLECrank-Nicolson法連立一次方程式P-10クーロン散乱Ubuntu境界条件MBEHiLAPW軸ラベルトランスCK1026MAS830L凡例PIC16F785LMC662AACircuit両対数グラフ片対数グラフグラフの分割specx.f

最新コメント
リンク

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