Scilabで数値積分:地球深部の密度と圧力

Scilabで数値積分: 固体の比熱ではintegrate関数を用いてデバイモデルから結晶固体の比熱を計算しました。
またScilabで乱数による関数の積分ではモンテカルロ法を用いて関数の積分を行いました。

これらはいずれも既知の関数f(x)について具体的な積分範囲、たとえばa≦x≦bでの積分値

\int_a^b f(x)\mathrm{d}x

を求める手法でした。

これらの方法では関数f(x)の形が式として分かっている必要があります。
しかしながら、離散的な数値データに対する積分値が計算したい場合もあります。

001_201311020335151bb.png

Fig.1: 地球の内部構造(世界初!地球中心部の超高圧高温状態を実現 ~ようやく手が届いた地球コア~より)


今回は地球内部の密度を表す数値表から、数値積分を使って重力加速度と圧力の深さ依存性を計算します。



標準地球モデル:PREM


地球内部の構造は、地震波観測からfig.1のような層構造をしていることが知られています。
もっとも有名な地球内部構造のモデルはDziewonski and Anderson (1981)によるPREM(Preliminary Reference Earth Model)と呼ばれるものです。

002_20131102033515160.png

Fig.2: 地球内部の地震波速度分布と密度分布。横軸は地球の中心からの距離(半径)。


PREMにはP波速度、S波速度、密度の深さごとの値、及びこれら3つの数値から計算できる重力加速度、圧力などが表として掲載されています。(参考:PREMの表地球を貫通するトンネル内の落下運動)

今回はPREMの密度から重力加速度と圧力を計算し、PREMの論文の中の値と比較して計算結果のチェックを行います。

重力加速度の計算


高校物理では重力加速度は約9.8[m/s]として計算をすると思います。
重力加速度はは大きな質量を持つ地球から万有引力によって引っ張られる力に起因しているため、地球内部では重力加速度の値が地球表面の値とは異なります。

地球の中心からr[m]での密度をρ(r)[g/m3]とするとr[m]より深い部分の質量M(r)[g]は

M(r) = 4 \pi \int_0^r u^2 \rho(u) \mathrm{d}u

重力加速度は、万有引力の法則より

g(r) = G\frac{M(r)}{r} = \frac{4\pi G}{r}\int_{0}^{r}u^2\rho(u)\mathrm{d}u

ここでGは万有引力定数G=6.67259×10-14[m3/s2/g]です。
PREMの数値データはPREM.txtに保存しfscanfMat等を用いてScilabから読み出しを行います。離散データの数値積分にはinttrapを用いました。

重力加速度計算の積分のところだけ抜き出すと、以下のようになります。

// *** 重力加速度の計算 ***
Gravity_num = zeros(Radius);
for i = 2:length(Gravity_num)
Gravity_num(i) = 4 * %pi * grav / (Radius(i) ^ 2) * inttrap(Radius(1:i), Radius(1:i) .^ 2 .* RHO(1:i));
end


圧力の計算


地球内部では非常に長い地質学的時間のなかで静水圧平衡が成り立っていると考えることが出来るため、圧力の計算が出来ます。

\frac{\mathrm{d}P}{\mathrm{d}r} = - \rho(r)g(r)

これを素直に積分すると以下の式が得られます。

\int\mathrm{d}P = - \int\rho(r)g(r)\mathrm{d}r

文字を置き換えて

\int\mathrm{d}Q = - \int\rho(u)g(u)\mathrm{d}u

地球表面rsでの圧力Ps(つまり大気圧)から地球深部rの圧力Pまで積分すると

\int_{P_s}^{P}\mathrm{d}Q = - \int_{r_s}^r\rho(u)g(u)\mathrm{d}u

大気圧はPs=1013.25[hPa]ですが、これは地球内部の圧力である数[GPa](数万気圧)と比べて充分小さいので左辺はほぼ地球深部の圧力Pと考えることが出来ます。

\int_{P_s}^{P}\mathrm{d}Q = P - P_s \simeq P

したがって圧力は以下のような積分で計算することが出来ます。

P = \int_{r}^{r_s}\rho(u)g(u)\mathrm{d}u

この積分を行う部分をScilabのプログラムから抜粋すると以下のようになります。

// *** 圧力の計算 ***
Pressure_num = zeros(Radius);
for i = 1:length(Pressure_num)-1
Pressure_num(i) = inttrap(Radius(i:$), RHO(i:$) .* Gravity_num(i:$));
end


計算結果


重力加速度と圧力の計算結果をFig.3に示します。

003_20131102033515738.png

Fig.3: 地球内部の重力加速度(上)と圧力(下)


赤の実線で示されたのがPREMの表にある値で、青の破線で示したのが今回積分から計算された値です。多少の誤差が発生しているようですが、おおよそ上手く計算できています。

全ソースコードはDeepEarth_sce.txtです。

clear;

// *** 物理定数 ***
grav = 6.67259e-14; // 万有引力定数 (m^3/s^2 g)

// *** PREMのテーブルを読み出し ***
X = fscanfMat('PREM.txt');
Radius = 1E3 * X(:,1); // 半径 (m)
Vp = X(:,2); // P波速度 (m/s)
Vs = X(:,3); // S波速度 (m/s)
RHO = 1E3 * X(:,4); // 密度 (g/m^3)
Ks = 1E12 * X(:,5); // 断熱体積弾性率(Pa)
Mu = 1E12 * X(:,6); // 剛性率(Pa)
Nu = 1E12 * X(:,7); //
Pressure = 1E12 * X(:,8); // 圧力(Pa)
Gravity = X(:,9); // 重力加速度 (m/s^2)

// *** 重力加速度の計算 ***
Gravity_num = zeros(Radius);
for i = 2:length(Gravity_num)
Gravity_num(i) = 4 * %pi * grav / (Radius(i) ^ 2) * inttrap(Radius(1:i), Radius(1:i) .^ 2 .* RHO(1:i));
end

// *** 圧力の計算 ***
Pressure_num = zeros(Radius);
for i = 1:length(Pressure_num)-1
Pressure_num(i) = inttrap(Radius(i:$), RHO(i:$) .* Gravity_num(i:$));
end

// *** 地震波速度と密度のプロット ***
scf(0)
// 地震波速度
xsetech([0,0,1,0.45]);
plot(1E-3 * Radius, 1E-3 * Vp,'-r');
plot(1E-3 * Radius, 1E-3 * Vs,'-b');
legend(["Vp","Vs"],1);
ylabel("Velocity (km/s)");
xlabel("Radius (km)");
// 密度
xsetech([0,0.45,1,0.45]);
plot(1E-3 * Radius, 1E-3 * RHO,'-r');
ylabel("Density (kg/s^3)");
xlabel("Radius (km)");

// *** 重力加速度と圧力のプロット ***
scf(1)
// 重力加速度
xsetech([0,0,1,0.45]);
plot(1E-3 * Radius,Gravity,'-r');
plot(1E-3 * Radius,Gravity_num,'--b');
legend(["Table","Numerical"],4);
ylabel("Acceleration of gravity (m/s^2)");
xlabel("Radius (km)");
// 圧力
xsetech([0,0.45,1,0.45]);
plot(1E-3 * Radius,1E-12 * Pressure,'-r');
plot(1E-3 * Radius,1E-12 * Pressure_num,'--b');
legend(["Table","Numerical"],1);
ylabel("Pressure (GPa)");
xlabel("Radius (km)");


台形積分とスプライン積分


今回使用したinttrapは、データ店の間を直線で補完する(要するに台形の面積を計算する)台形積分です。
もう少し高度な計算をするスプライン積分もintsplinを使って計算することが出来ます。
入力するパラメータも同じなので、多くの場合そのまま置き換えることが出来るでしょう。

ただし今回のPREMのようにデータが不連続となる点で複数の値を持っているような表の積分を行うときにはintsplinでは下記のようなエラーが出ます。

!--error 999 
splin: 入力引数 #1 の値が間違っています: (厳密に)昇順ではありません, もしくは +-inf を検出しました.


このような場合は、今回のように台形積分inttrapを使います。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 数値積分  マントル 

Scilabで乱数による関数の積分

Scilabには関数の積分の計算を行うintegrateがあるため簡単な関数の積分には必要が無いのですが、乱数を使った積分の方法について書きます。
[0,1]の範囲で一様分布な乱数列を{Xk}k=1,2...nとすると関数f(x)の積分は以下のようにして計算することが出来ます。

\begin{eqnarray}\int_a^b f(x) \mathrm{d}x & = & \int_0^1 f(a+(b-a)r)(b-a)\mathrm{d}r \\ & = & \lim_{n \to \infty} \frac{1}{n}\sum_{k=1}^{n}f(a+(b-a)X_k)(b-a) \end{eqnarray}


乱数による積分の考え方


[0,1]の範囲の一様分布な乱数列{Xk}k=1,2...nはScilabを使うと簡単に得られます。(参考:Scilabで乱数の生成,Scilabでコイン投げ)

このとき関数f(x)の積分は、f(Xk)の平均から求めることが出来るでしょう。

\int_{0}^{1}f(x)\mathrm{d}x = \lim_{n \to \infty} \frac{1}{n} \sum_{k=1}^{n}f(X_k)

例として、以下に挙げる具体的な関数f(x)に対して積分の計算をします。

f(x) = \exp(-x^2)

関数f(x)のグラフはFig.1のような形になります。積分は0≦x≦1の範囲の面積を求めることと同じです。


001_20130926022100898.png
Fig.1: 関数f(x)=exp(-x2)のグラフ


Scilabで楽しむ確率論(PDF)では、本文中に書いてある数式と実際にScilabで計算している数式が異なっています。本エントリではプログラムの方にあわせることにします。

モンテカルロ法以外の計算法


実際のところ、この程度の関数なら乱数を用いたモンテカルロ法を使うよりも、素直に数値積分をするほうがはるかに良い結果が得られます。

Scilabではintegrateをつかって簡単に数値積分が出来ます。(参考:Scilabで数値積分: 固体の比熱)
また、もっと基本的な方法として、解析的に原始関数を決められないか考えるという方法もあります。Maximaを使って不定積分を行うと以下のような結果が得られました。

\int f(x)\mathrm{d}x = \frac{\sqrt{\pi} \mathrm{erf}(x)}{2} + C(C: 積分定数)

まあこの場合は誤差関数erfが入っているので数値積分から逃れられてはいないというのも事実ではありますが。

数値積分、解析的な積分、乱数を使った積分のそれぞれの方法で積分値を計算するScilabのプログラムがfunc_sce.txtです。
数値積分と解析的な積分の両方の方法で、積分の値は0.7468241となりました。モンテカルロ法で求めた値もおおよそ同じくらいになります。

clear;

// *** 関数の定義 ***
// 関数f(x)
deff('y = f(x)','y = exp(- x .^ 2)');
// 原始関数F(x)
deff('y = F(x)','y = sqrt(%pi) .* erf(x) ./ 2');

X = linspace(0,5);

// *** グラフの描画 ***
plot(X,f(X));
xlabel("x");
ylabel("f(x)");

// *** 積分の計算 ***
// 解析的な積分
F(1) - F(0)
// 数値積分
integrate('f(x)','x',0,1)
// モンテカルロ積分
n = 1000;
mean(f(rand(1,n)))


ヒストグラム


このモンテカルロ法による積分を複数回行ったときのヒストグラムをFig.2へ示します。
確かにintegrateで求めた値を中心にばらついていることが分かります。


002_20130926022059bb7.png
Fig.2: 積分結果のばらつきを示すヒストグラム


積分範囲の変更


次に積分範囲がa≦x≦bの場合へ拡張します。方法は単純モンテカルロ積分 (2)の通りです。

a \leq x \leq b
0 \leq x - a \leq b - a
0 \leq \frac{x - a}{b - a} \leq 1

\frac{x-a}{b-a} = r

とおくと0≦r≦1となるので、置換積分を考えればよいことになります。

x = a+(b-a)r
\frac{\mathrm{d}x}{\mathrm{d}r} = b-a
\mathrm{d}x=(b-a)\mathrm{d}r

したがって

\begin{eqnarray}\int_a^b f(x) \mathrm{d}x & = & \int_0^1 f(a+(b-a)r)(b-a)\mathrm{d}r \\ & = & \lim_{n \to \infty} \frac{1}{n}\sum_{k=1}^{n}f(a+(b-a)X_k)(b-a) \end{eqnarray}

となります。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 数値積分 乱数 モンテカルロ解析 Maxima 

Scilabで二酸化炭素の状態方程式 その2

Scilabで二酸化炭素の状態方程式 その1では、Scilabを用いてファンデルワールスの状態方程式を解くプログラムを作成しました。

\left(P + \frac{a^2}{V} \right) (V-b)=RT

(P: 圧力, V: モル体積, R: 気体定数, T: 絶対温度, a, b: ファンデルワールス定数)

実はファンデルワールスの状態方程式は気体の体積だけでなく、液体や気体・液体共存のときのP,V,Tも計算することが出来ます。
この計算を行う際にMaxwellの規則を利用します(参考:ファン・デル・ワールスの状態方程式 (クラウジウス=クラペイロンの式、ジュール=トムソン効果)3.液体・気体共存領域)。Maxwellの規則ではグラフの面積が等しくなる条件を用いますが、こういった条件を探すのは、解析的には難しく、数値計算の出番となります。

004_20130811163805e0f.png

Fig.1: T = 250 (K)のときの二酸化炭素の圧力Pとモル体積Vの関係。グラフ中の赤の水平線の部分が液体・気体の共存領域、それよりも右側が気体の安定領域で左側が液体の安定領域。赤の水平線は、赤の水平線と赤の点線で囲まれた二つの領域の面積が等しくなる条件から引かれる。(Maxwellの規則)



二酸化炭素の相図


常温・常圧の二酸化炭素は気体ですが、固体の二酸化炭素といえばドライアイスです。ドライアイスは室温においておくと液体にならずに気体になりますが、5.1×105(Pa)つまり5気圧以上では水などと同様に液体になります。
以下に示すのは、ファン・デル・ワールスの状態方程式 (クラウジウス=クラペイロンの式、ジュール=トムソン効果)から引用した二酸化炭素の相図です。相図(Phase diagram, 状態図とも)は物質の状態(固体,液体,気体)と温度・圧力・化学組成の関係性を表した図です。



液体・気体共存領域


ファンデルワールスの状態方程式の1つの特徴として気体と液体の相転移を表すことができる点が挙げられます。
このときの液体、気体のそれぞれのモル体積は、Fig.1に示したとおりファンデルワールスの状態方程式と水平線の一番左の交点と一番右の交点として表されます。
またこの水平線は、水平線と状態方程式の作る二つの領域の面積が等しくなる条件から引くことが出来ます。
これをMaxwellの規則と呼び、数式で表すと以下のようになります。

\int_{V_l}^{V_g}P_{vdw}\mathrm{d}V - \int_{V_l}^{V_g}P_k \mathrm{d}V = 0

ただしPvdwがファンデルワールスの状態方程式で計算した圧力、Pkが水平線の圧力です。これらの交点のうち一番左の物が液体の体積Vlで一番右のものが気体の体積Vgです。

数値計算


色々な教科書に水平線Pkは面積が等しくなる条件から決定することが出来ると簡単に書かれていますが、Pkを決めないと交点Vg,Vlが決まらないため、Pkを決めるためには数値計算が必要です。

この計算にはScilabで金属の化学ポテンシャルScilabでおもりの吊り下げで利用した非線形方程式ソルバfsolveが使えます。

まずPkとPvdw(V)が与えられているとして3つの交点の値を求めます。
交点の値は、ファンデルワールスの状態方程式を変形して得られる以下の三次方程式の解です。

P_k V^3 - (b P_k + RT) V^2 + aV - ab = 0

Scilabでは多項式の解はrootsを用いて簡単に計算できます。解は複素数の要素を持つ縦ベクトルになります。Scilab言語には(C言語のような)型の宣言が無いのでミスを犯しやすいのですが、解の全ての要素が実部しか持たなかったとしても、変数の型としては複素数なので、realをつかって実数型へ変換する必要があります。

数値積分にはScilabで数値積分: 固体の比熱Scilabで関数フィッティング: 金属の電気抵抗で使ったintegrateを利用します。

最終的なプログラムはCO2-subcritical_sce.txtです。

clear;

// *** 入力パラメータ ***
// 物理定数
r = 8.31 // 気体定数 (J/K/mol)

// 二酸化炭素のファンデルワールス状態方程式
a = 3.65E-1; // Pa m^6 / mol^2
b = 4.28E-5; // m^3 / mol
// 温度 (K)
t = 250;

// *** ファンデルワールスの状態方程式 **
function p = pvdw(v)
p = r .* t ./ (v - b) - a ./ v ./ v
endfunction

// *** 解くべき方程式 ***
function e = func(pk)
Vx = roots([pk, -1 * (b * pk + r * t), a, - a * b]);
vg = max(real(Vx));
vl = min(real(Vx));
s1 = integrate('pvdw(v)','v',vl,vg);
s2 = (vg - vl) * pk;
e = s1 - s2;
endfunction

// *** 非線形方程式の数値解 ***
pk0 = %eps;
pk = fsolve(pk0,func);

Vx = roots([pk, -1 * (b * pk + r * t), a, - a * b]);
vg = max(real(Vx));
vl = min(real(Vx));

// *** 圧力の計算とプロット ***
// モル体積ベクトル
Vl = linspace(b,vl,1000) + %eps; // (m^3 / mol)
Vm = linspace(vl,vg,1000); // (m^3 / mol)
Vg = linspace(vg,1e-3,1000); // (m^3 / mol)
// 温度
plot(Vl,pvdw(Vl),'-r');
plot(Vm,pk,'-r');
plot(Vm,pvdw(Vm),'--r');
plot(Vg,pvdw(Vg),'-r');
// *** グラフの装飾 ***
xlabel("Molar volume (m^3/mol)");
ylabel("Pressure (Pa)");
plot([b,b],[-1E8,2E7],'--k');
plot([0,0.001],[0,0],'--k');
zoom_rect([0,-2E6,1E-3,2E7]);


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 非線形方程式ソルバ 数値積分 状態方程式 熱力学 

Scilabで単振り子 その1 解析解との比較

これまで微分方程式による物理現象のモデル化(PDF)から運動学の章で紹介されている微分方程式をScilabを用いて計算してきました。
今回から4回は単振り子の微分方程式をScilabで計算します。

001_20130720182725.png

Fig.1: 単振り子の問題設定



単振り子の微分方程式と解析解


導出は次回に回しますが、単振り子の運動を表す微分方程式は、以下のように書くことが出来ます。

\frac{\mathrm{d}\theta}{\mathrm{d}t} = q
\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2} = - \frac{g}{L}\sin\theta

(θ: 振れ角, t: 時間, q: 角速度, g: 重力加速度, L: 糸の長さ)

この微分方程式は、解析的に解を求めることが出来ます。
例えば、初角度θ0から初角速度q0=0で振動を開始したとすると下記のようになります。

t(\theta) = - \sqrt{\frac{L}{2g}}\int^{\theta}_{\theta_0}\frac{\mathrm{d}\theta}{\sqrt{\cos\theta - \cos\theta_0}}

しかしながら、この積分は初等関数で表せる形に出来ないので、結局は数値積分をせざるを得ません。

今回は、これまでと同様にScilabの常微分方程式ソルバodeを用いて振りこの運動をシミュレーションするとともに、解析解に対する数値積分も行い、これら二つを比較します。

プログラミング


常微分方程式ソルバはode関数を使います。(参考:常微分方程式のタグが付いたエントリ)

解析解のほうの数値積分はintegrate関数を用いました。(参考:Scilabで数値積分: 固体の比熱,integrate - 求積法による積分)

// *** 解析解とソルバ解の共通部分 ***
// 初期振幅の入力
g = 9.8; // 重力加速度
l = g / (2 * %pi) ^ 2; // 糸の長さ
// 初期条件
th0 = input("th0 = "); // 初角度
q0 = 0; // 初角速度
X0 = [th0; q0];

// *** 常微分方程式ソルバによる解 ***
// 解くべき微分方程式の定義
function dx = pend(t,x)
dx(1) = x(2);
dx(2) = - g / l * sin(x(1));
endfunction
T = linspace(0,2,200); // 時間ベクトル
TH = ode(X0, 0, T, pend); // 常微分方程式ソルバ
plot(T,TH(1,:),'-r'); // プロット

// *** 解析解 ***
THanaly = linspace(- th0, th0, 20);
Tanaly = - sqrt(l / (2 * g)) * integrate('1 ./ sqrt(cos(th) - cos(th0))','th',th0,THanaly);
plot(Tanaly,THanaly,'og');

// *** グラフの体裁 ***
legend("o.d.e","integration",1);
xlabel("t[s]");
ylabel("$\theta \mathrm{[rad]}$");


計算結果


以下に、初角速度q0=0で初角度θ0を0.1ラジアン、2.9ラジアンの2通りの値で計算した結果を示します。

002_20130720182725.png

003_20130720182724.png

Fig.2-3: 厳密な単振り子の解。sinθ≒θの近似が成り立つときの周期が1秒となるように糸の長さLを調整している。振幅(初角度)が小さいとき(θ0=0.1; グラフ上)は周期が1秒となっているが、振幅(初角度)が大きいとき(θ0=2.9; グラフ下)は周期が1秒を大きく超えている。


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 常微分方程式 ode 数値積分 単振り子 

Scilabで金属の化学ポテンシャル

Scilabの非線形方程式ソルバfsolveを用いて、積分を含む方程式である金属の化学ポテンシャルを計算しました。
対象とする系は、自由電子近似が良く成り立つはずのアルミニウムとし、数値計算による化学ポテンシャルの計算結果とゾンマーフェルト展開を用いた近似値の比較を行いました。

その結果、この系に関しては、数値計算とゾンマーフェルト展開を用いた計算結果がアルミニウムの融点以上までほとんど同じ結果を示すことがわかりました。

001_20130609181459.pngFig.1: アルミニウムの化学ポテンシャルの温度依存性、ゾンマーフェルト展開を用いた近似(青実線)とfsolveを利用した数値計算(赤点線)。インセットは自由電子近似による電子の状態密度(青)、低温時(10K,緑)と高温時(赤)の電子の分布は絶対零度の状態密度に各温度のフェルミ分布関数を掛けたもの、エネルギーの原点はフェルミ準位。



金属電子論と数値計算


Scilabで差し込みグラフ:金属の比熱では実験的に求められたアルミニウムの電子比熱係数γ=1.35(mJ/mol/K^2)という値をもとに電子比熱の温度依存性をグラフにプロットしました。

これに対して、次は理論的にアルミニウムの電子比熱を予想してみます。
差し当たり今回は、電子比熱を計算する前段階として、アルミニウムの化学ポテンシャルを自由電子近似と数値計算から求めます。

金属の電子構造を計算するに当たり、次の2つの近似がよく使われます。
  • 自由電子近似(free electron approximation)
  • ゾンマーフェルト展開(Sommerfeld expansion)

自由電子近似はナトリウムやアルミニウムといった比較的単純な金属に対しては良い近似になることが知られています。また、ゾンマーフェルト展開は温度が充分低い条件で成り立ちます。

ゾンマーフェルト展開は、複雑な複雑な積分の計算を、比較的計算しやすい解析的な式に変形するための級数展開なのですが、Scilabなどで数値計算を行うことでこの近似に頼らずに電子比熱が計算できます。
今回は「自由電子近似+ゾンマーフェルト展開」と「自由電子近似+数値計算」の両方を行い比較します。

金属の化学ポテンシャル


金属の電子比熱は、以下の式で定義されます。

c_e = \frac{\mathrm{d} u_e}{\mathrm{d} T}

ここで、電子系の内部エネルギーは

u_e = \int^{\infty}_{-\infty}\epsilon f(\epsilon) D(\epsilon)\rm{d}\epsilon

さらにD(ε)は電子の状態密度(density of states)でf(ε)はフェルミ・ディラック分布関数です。

f(\epsilon) = \frac{1}{\exp{\left(\frac{\epsilon-\mu}{k_B - T}\right)}+1}

まずD(ε)ですが、これは金属の個性を現す重要なパラメータなのですが、自由電子近似と金属の体積、価電子数から求めることが出来ます(後述)。
一方で、フェルミ・ディラック分布関数内の化学ポテンシャルμは、温度の関数となり少々厄介です。本エントリの主題は、この化学ポテンシャルをScilabの非線形方程式ソルバーを用いて数値的に計算し、ゾンマーフェルト展開によって近似的に求められた化学ポテンシャルと比較を行うことです。

自由電子近似


詳細には立ち入りませんが、自由電子近似による金属の状態密度とフェルミエネルギー(絶対零度における化学ポテンシャル)は、リュードベリ原子単位系では以下の様に表されます。

D(\epsilon) = \frac{V}{2 \pi^2} \epsilon^{\frac{1}{2}}

\epsilon_F = \left( \frac{3 \pi^2 n}{V} \right)^{\frac{2}{3}}

ここでVは単位物質量あたりの体積です。具体的には格子体積やモル体積を用います。今回は差し当たり単位原子辺りの体積(atomic volume)を用いて計算します。nは電子の数で、今回は単位原子辺りの電子の数なので価電子数と一致します。(参考:金持 徹著固体電子論)

自由電子近似は、ナトリウムやアルミニウムでは良く成り立ちますが、遷移金属などには適用できません。

ゾンマーフェルト展開


これも詳細には立ち入らず結論だけ書きます。
ゾンマーフェルト展開の結果、内部エネルギーue、化学ポテンシャルμ、電子比熱ceは

u_e = u_{e0} + \frac{\pi^2}{6}k_B^2 D(\epsilon_F)T^2

\mu = \epsilon_F - \frac{\pi^2}{6}k_B^2 \frac{D^{'}(\epsilon_F)}{D(\epsilon_F)}T^2

c_e = \frac{\pi^2}{3}k_B^2 D(\epsilon_F)T

となります。(参考:無機化学特論I ゾンマーフェルト展開について)

ゾンマーフェルト展開は充分低温のみで適用できます。(参考:Sommerfeld 展開の導出- 武智公平)

自由電子近似+ゾンマーフェルト展開


ゾンマーフェルト展開に自由電子近似の結果を代入すると、金属の化学ポテンシャルや電子比熱の値が具体的に計算できます。

\mu = \epsilon_F \frac{\pi^2 k_B}{12 \epsilon_F}T^2

c_e = \frac{\pi^2 k_B^2 n}{2 \epsilon_F}T

押さえておかなくてはならない重要な点は、「自由電子近似」と「ゾンマーフェルト展開」は相互に全く関係の無い別の近似であり、適用できる条件もそれぞれ異なるということです。

数値計算


さて、本題です。
ここまでで「自由電子近似+ゾンマーフェルト展開」によって化学ポテンシャル(と電子比熱)を計算する準備が整いました。
ここからは、ゾンマーフェルト展開を使わずに化学ポテンシャルを計算することを考えます。

価電子数は温度に依存しないので

n = \int_{-\infty}^{\infty} D(\epsilon)f(\epsilon)\mathrm{d}\epsilon

をμに付いて解く事によって、化学ポテンシャルμをそれぞれの温度で計算することが出来ます。
この式は、解析的には解けないため、Scilabの非線形ソルバfsolveを利用します。

作成したソースコードはChemicalPotential_sce.txtです。

計算結果(自由電子近似+数値計算)


Fig.1に化学ポテンシャルの計算結果を示します。自由電子近似を用いたアルミニウムの場合、ゾンマーフェルト展開を用いた近似と数値積分を利用して求めた結果はほとんど変わらない結果を示します。(高温で若干ずれているのは、本当にゾンマーフェルト展開が妥当でなくなっていっているのか、数値計算の精度の問題かはちょっと自信がないです。計算するエネルギー範囲が狭いせいではなさそう・・・。)

アルミニウムの融点は660.32 ℃(933.47 K)なので、固体の間はゾンマーフェルト展開が破綻することは無いであろう事がわかります。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 非線形方程式ソルバ 電子比熱 数値積分 

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

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoCOPアンプCPAPIC強磁性モンテカルロ解析常微分方程式トランジスタode状態密度インターフェースecalj定電流スイッチング回路PDS5022DOS半導体乱数シェルスクリプトレベルシフトHP6632Aブレッドボード分散関係温度解析トランジスタ技術R6452A可変抵抗I2Cセミナー確率論反強磁性非線形方程式ソルバ絶縁偏微分方程式バンド構造熱設計数値積分バンドギャップカオスA/DコンバータフォトカプラシュミットトリガGW近似LEDLM358ISO-I2C三端子レギュレータ数値微分サーボ直流動作点解析カレントミラーマフィンティン半径TL431PC817C発振回路74HC4053USBアナログスイッチbzqltyFFTチョッパアンプ2ちゃんねる補間量子力学開発環境電子負荷標準ロジックパラメトリック解析アセンブラ基本並進ベクトルブラべ格子単振り子BSchLDAイジング模型繰り返しMaximaキュリー温度位相図状態方程式失敗談スピン軌道相互作用六方最密充填構造相対論FET抵抗コバルト不規則合金TLP621ewidthGGAQSGWgfortranランダムウォークラプラス方程式スイッチト・キャパシタcygwin熱伝導SMPスレーターポーリング曲線三角波格子比熱LM555条件分岐TLP552MCUNE555UPSTLP521QNAPマントルテスタFXA-7020ZR過渡解析詰め回路ガイガー管ダイヤモンド自動計測Writer509データロガー固有値問題VESTAスーパーセルOpenMP差し込みグラフ平均場近似起電力awk仮想結晶近似VCAubuntufsolveブラウン運動熱力学第一原理計算井戸型ポテンシャルシュレディンガー方程式面心立方構造fccウィグナーザイツ胞interp12SC1815L10構造非線型方程式ソルバFSMキーボードTeX結晶磁気異方性初期値OPA2277化学反応等高線ジバニャン方程式ヒストグラム確率論三次元フィルタRealforcePGAフェルミ面正規分布固定スピンモーメント全エネルギースワップ領域リジッドバンド模型edeltquantumESPRESSOルチル構造岩塩構造二相共存ZnOウルツ鉱構造BaOフォノンデバイ模型multiplotgnuplotc/aノコギリ波合金クーロン散乱ハーフメタル半金属CapSenseマンデルブロ集合マテリアルデザインSICGimpCK1026MAS830L円周率トランスPIC16F785凡例線種シンボルLMC662ヒストグラム不規則局所モーメント文字列疎行列TS-110TS-112Excel直流解析等価回路モデル入出力トラックボールPC軸ラベルAACircuitP-10フラクタル境界条件連立一次方程式Ubuntuifortパラメータ・モデルspecx.f関数フィッティング最小二乗法Crank-Nicolson法陰解法日本語EAGLEMBEグラフの分割負帰還安定性ナイキスト線図熱拡散方程式HiLAPW両対数グラフ片対数グラフ縮退

最新コメント
リンク

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