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 非線形方程式ソルバ 数値積分 状態方程式 熱力学 

comment

Secret

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

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性OPアンプPICCPA常微分方程式モンテカルロ解析ecaljodeトランジスタ状態密度インターフェースDOS定電流スイッチング回路PDS5022半導体シェルスクリプトレベルシフト乱数HP6632AR6452AI2C可変抵抗分散関係トランジスタ技術ブレッドボード温度解析反強磁性確率論バンドギャップセミナー数値積分熱設計非線形方程式ソルババンド構造絶縁偏微分方程式ISO-I2CLM358フォトカプラ三端子レギュレータカオスLEDシュミットトリガGW近似A/Dコンバータ発振回路PC817C直流動作点解析USBマフィンティン半径数値微分アナログスイッチTL43174HC4053カレントミラーサーボ量子力学単振り子チョッパアンプ補間2ちゃんねる開発環境bzqltyFFT電子負荷LDAイジング模型BSch基本並進ベクトルブラべ格子パラメトリック解析標準ロジックアセンブラ繰り返し六方最密充填構造SMPコバルトewidthFET仮想結晶近似QSGW不規則合金VCAMaximaGGA熱伝導cygwinスレーターポーリング曲線キュリー温度スイッチト・キャパシタ失敗談ランダムウォークgfortran抵抗相対論位相図スピン軌道相互作用VESTA状態方程式TLP621ラプラス方程式TLP552条件分岐NE555LM555TLP521マントル詰め回路MCUテスタFXA-7020ZR三角波過渡解析ガイガー管自動計測QNAPUPSWriter509ダイヤモンドデータロガー格子比熱熱力学awkブラウン運動起電力スーパーセル差し込みグラフ第一原理計算フェルミ面fsolve最大値xcrysden最小値最適化ubuntu平均場近似OpenMP井戸型ポテンシャルシュレディンガー方程式固有値問題2SC1815結晶磁気異方性OPA2277非線型方程式ソルバTeXgnuplot固定スピンモーメントFSMPGAc/a全エネルギーfccフラクタルマンデルブロ集合正規分布縮退初期値interp1multiplotフィルタ面心立方構造ウィグナーザイツ胞L10構造半金属二相共存SICZnOウルツ鉱構造BaO重積分クーロン散乱磁気モーメント電荷密度化学反応CIF岩塩構造CapSenseノコギリ波デバイ模型ハーフメタルキーボードフォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt合金等高線線種凡例シンボルトラックボールPC軸ラベルグラフの分割トランス文字列CK1026MAS830L直流解析Excel不規則局所モーメントパラメータ・モデル入出力日本語最小二乗法等価回路モデルヒストグラムGimp円周率TS-110TS-112PIC16F785LMC662三次元specx.fifortUbuntu疎行列不純物問題Realforceジバニャン方程式ヒストグラム確率論マテリアルデザインP-10境界条件連立一次方程式熱拡散方程式AACircuitHiLAPW両対数グラフ片対数グラフ陰解法MBEナイキスト線図負帰還安定性Crank-Nicolson法EAGLE関数フィッティング

最新コメント
リンク

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