物性推算の基本について

状態方程式型 / Equation of State

状態方程式型の物性推算法:Soave-Redlich-Kwong(SRK)、Peng-Robinson(PR)、Benedict-Webb-Rubin-Starling(BWRS)を用いた各種物性の求め方の基本をまとめていきます。

SRK式:密度、エンタルピー、エントロピー、蒸気圧の算出

PR式:密度、エンタルピー、エントロピー、蒸気圧の算出

アルファ関数について

状態方程式の混合則修正式

BWRS式

活量係数型 /Liquid Activity

活量係数型の物性推算法:Wilson、Non-random two-liquid(NRTL)を用いた各種物性の求め方の基本をまとめていきます。

活量係数型相平衡の基本

Wilson式

NRTL式

Binary Parameterの回帰計算について

Vapor Fugacityの補正について

石油成分 / Assay

石油成分の物性推算法についてまとめていきます。

単成分の物性推算

蒸留曲線からの物性推算

Soave-Redlich-Kwong(SRK)式

SRK式

$$p = {\frac{RT}{v-b} - \frac{\alpha(T)a}{v(v+b)}} \tag{1.1}$$ $$\alpha(T) = [1+m(1-\sqrt{Tr})]^2 \tag{1.2}$$ $$m = 0.48+1.574\omega-0.176\omega^2 \tag{1.3}$$

混合則(i成分とj成分の相互作用パラメーター)

$$a = \sum_{i} \sum_{j}x_ix_j(a_ia_j)^{0.5}(1-k_{ij}) \tag{1.4}$$ $$b = \sum_{i} x_ib_i \tag{1.5}$$

$$a_i = 0.42748\frac{R^2T^2_{ci}}{P_{ci}} \tag{1.6}$$ $$b_i = 0.08664\frac{RT_{ci}}{P_{ci}} \tag{1.7}$$

※ωはPitzerの偏心因子。

※以下Soave氏の論文ではaiの係数(Redlich Kwong定数)は0.42747となっておりますが、臨界物性から求めますと、正確には0.42748となります。

Soave,G.;"Equilibrium constants from a modified Redkh-Kwong equation of state"; Chemical Engineering Science, 1972, Vol.27, p.1197-1203.

Redlich Kwong定数の導出はこちらをご参照ください。

本サイトでは、aiは0.42748として計算しております。

SRKの液密度修正推算式例

状態方程式型の推算法では、液密度の推算精度がよくありません。

このため、液密度の推算はAPI法、RACKET法、COSTALD法等がよく用いられますが、ここではSRK式から求まるVに対する修正式としてSRK-Peneloux式をご紹介します。

SRK-Peneloux式は以下になります。

$$V_{correct}^L = \sum_{i=1}^n V_{SRK}^L-0.40768(0.29441-Z_{RA})\frac{RT_{ci}}{P_{ci}}x_i \tag{1.8}$$ $$Z_{RA}=0.29056-0.08775\omega_i \tag{1.9}$$

Peneloux式の推算精度の確認については、こちらをご参照ください。

上記のリンク先ツールでCalculateボタンを押していただくと、各物性推算法の密度を求めることができます。

※silverlightで作成しております。chrmoeではバージョンによって表示されないので、その場合Internet Exploereでご確認ください

・参考としてモル濃度:Methane:0.3[-]、Ethane:0.3[-]、Proapne:0.4[-]、圧力:1000[kPa]、温度:100[K]の密度計算を確認します。

REFPROPで上記の密度を求めると、642.47[kg/m3]となりますが、SRK-Peneloux式で比較的精度良い結果が得られることが確認できます。

※SRK-Peneloux式の出展

Peneloux A, Rauzy E and Freze R, ; "A consistent correction for Redlich-Kwong-Soave volumes"; Fluid Phase Equilibria,1982, Vol.8, no.1, p.7-23

また、もう一つの密度の補正計算式として、Claude氏が1993年に発表している補正式も合わせてご紹介します。

SRKのエンタルピー算出▼


理想気体のエンタルピー

理想気体のエンタルピーは温度の関数で表現することができます。

Perryのhandbookでは理想気体の比熱が温度の関数として記載されています。

  

理想気体のエンタルピーと実在気体のエンタルピーの偏倚は以下の式で表されます。

なお、偏倚の式導出はこちらをご参照ください。

$$\left[h-h^{ideal}\right]_{T,P}=\int_v^\infty\left[P-T\left(\frac{\partial P}{\partial T}\right)_v\right]dv+Pv-RT \tag{1.10}$$

ここで、SRKであれば、圧力は1.1式で表せるため、\(\left(\frac{\partial P}{\partial T}\right)_v\)は、以下となります。

$$\left(\frac{\partial P}{\partial T}\right)_v=\frac{R}{v-b}-\frac{\partial a (T)}{\partial T}\frac{1}{v\left(v+b\right)}$$

よって、1.10式は、以下となります。

$$\left[h-h^{ideal}\right]_{T,P}=\int_v^\infty\left[P-T\frac{R}{v-b}-\frac{\partial a (T)}{\partial T}\frac{1}{v\left(v+b\right)}\right]dv+Pv-RT \tag{1.11}$$ $$\left[h-h^{ideal}\right]_{T,P}=\frac{1}{b}\int_v^\infty\left[\left(\frac{\partial a (T)}{\partial T}T-a (T)\right) \left(\frac{1}{v}-\frac{a}{v+b}\right)\right]dv+Pv-RT \tag{1.12}$$ $$\left[h-h^{ideal}\right]_{T,P}=RT\left(\frac{1}{bRT}\left(\frac{\partial a (T)}{\partial T}\right)\right)\int_v^\infty\left[ln\frac{v}{v+b}\right]dv+Pv-RT \tag{1.13}$$

ここで、a(T)の微分を求めておきます。SRK(1972)のα関数は(1.2)式(以下再掲します)の形ですので、これを代入し、微分していきます。

$$\alpha(T) = [1+m(1-\sqrt{Tr})]^2 \tag{1.2}$$ $$\frac{\partial a(T)}{\partial T} = a\frac{\partial \alpha(T)}{\partial T} \tag{1.14}$$ $$\frac{\partial a(T)}{\partial T} = 2a\sqrt{\alpha (T)}\frac{\partial}{\partial T}\left(1+m\left(1-\sqrt{\frac{T_c}{T}}\right)\right) \tag{1.15}$$ $$\frac{\partial a(T)}{\partial T} = 2a\sqrt{\alpha (T)} \cdot \frac{1}{2} \cdot \left(\frac{-m}{\sqrt{T \cdot T_c}}\right) \tag{1.16}$$ $$\frac{\partial a(T)}{\partial T} = -ma\sqrt{\frac{\alpha (T)}{T \cdot T_c}} \tag{1.17}$$

a(T)の微分\(\frac{\partial a(T)}{\partial T}\)が求まったので、式1.17を式1.13に代入します。

$$\left[h-h^{ideal}\right]_{T,P}=RT\left(\frac{1}{bRT}\left(-maT\sqrt{\frac{\alpha (T)}{T \cdot T_c}}\right)\right)\int_v^\infty\left[ln\frac{v}{v+b}\right]dv+Pv-RT \tag{1.18}$$

これでPVTだけの式になったので、P、TからVを求めてエンタルピーを求めることができます。

上式のままだと計算しづらいので以下少し整理します。

準備中

SRKのエントロピー算出

準備中。

SRKの純成分蒸気圧算出

状態方程式法である温度における純成分の蒸気圧を求めるには、液相、気相のフガシティーが等しくなる圧力を繰り返し計算で求める必要があります。

以下蒸気圧の初期値として、Tc、Pc、NBPから求めていきます。

Peng-Robinson(PR)式

PR式(1976年)

$$p = {\frac{RT}{v-b} - \frac{\alpha(T)a}{v(v+b)+b(v-b)}} \tag{2.1}$$ $$\alpha(T) = {1+m(1-\sqrt{Tr})}^2 \tag{2.2}$$ $$m = 0.37464+1.54226\omega-0.26992\omega^2 \tag{2.3}$$ $$a = \sum_{i} \sum_{j}x_ix_j(a_ia_j)^{0.5}(1-k_{ij}) \tag{2.4}$$ $$b = \sum_{i} x_ib_i \tag{2.5}$$ $$a_i = 0.45724\frac{R^2T^2_{ci}}{P_{ci}} \tag{2.6}$$ $$b_i = 0.0778\frac{RT_{ci}}{P_{ci}} \tag{2.7}$$

Cf)PR式(1978年)

蒸気圧の精度を向上させた1978年の式を以下ご紹介します。

$$p = {\frac{RT}{v-b} - \frac{\alpha(T)a}{v(v+b)+b(v-b)}} \tag{2.8}$$ $$\alpha(T) = {1+m(1-\sqrt{Tr})}^2 \tag{2.9}$$ $$m = 0.37464+1.54226\omega-0.26992\omega^2 (\omega \leqq 0.491) \tag{2.10}$$ $$m = 0.379642+1.48503\omega-0.164423\omega^2 + 0.016666\omega^3 (\omega \gt 0.491) \tag{2.11}$$ $$a = \sum_{i} \sum_{j}x_ix_j(a_ia_j)^{0.5}(1-k_{ij}) \tag{2.12}$$ $$b = \sum_{i} x_ib_i \tag{2.13}$$ $$a_i = 0.457235529\frac{R^2T^2_{ci}}{P_{ci}} \left(1+m_i \cdot \left(1- \sqrt{ \frac{T}{T_{ci}} } \right)\right)^2 \tag{2.14}$$ $$b_i = 0.0777960739\frac{RT_{ci}}{P_{ci}} \tag{2.15}$$

※ここで、\(a_i\)、\(b_i\)の係数は、それぞれ以下の算式より求められています。

$$X = \frac{-1+\sqrt[3]{6\sqrt{2}+8}-\sqrt[3]{6\sqrt{2}-8}}{3} \approx 0.253076587 \tag{2.16}$$

\(a_i\)の係数

$$ \frac{8 \cdot \left(5X+1\right)}{49-37X} \tag{2.17}$$

\(b_i\)の係数

$$ \frac{X}{X+3} \tag{2.18}$$

※PR(1978)式の出展

Peng, D. Y, and Robinson, D. B. ; "The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs"  Gas Processors Association, 1978,

Cf)PR式(1976年)とPR式(1978年)の差異について

偏心因子の大きい成分(アルカンであれば、カーボン数13以上)について、1976式と1978式でどの程度蒸気圧の精度に差異が生じるのか確認してみます。

準備中

PRの液密度修正推算式例

状態方程式型の推算法では、液密度の推算精度がよくありません。

このため、液密度の推算はAPI法RACKET法COSTALD法等を用いますが、ここではPR式から求まるVに対する修正式としてPR-Peneloux式をご紹介します。

Peneloux式のvolume補正係をvolume shift parameter推算精度の確認については、こちらをご参照ください。

$$V_{correct}^L = \sum_{i=1}^n V_{PR}^L-0.11537(0.44064-Z_{RA})\frac{RT_{ci}}{P_{ci}}x_i \tag{2.8}$$ $$Z_{RA}=0.29056-0.08775\omega_i \tag{2.9}$$

Peneloux式の推算精度の確認については、こちらをご参照ください。

PRのエンタルピー算出

準備中。

PRのエントロピー算出▼

準備中。

状態方程式法の解法について

状態方程式法は式変形によってz、あるいはvの3次方程式に変形することが可能です。

具体的な式変形の流れは以下をご参照ください。

例としてSRK式を3次式の形に変形していきます。(1.1)式からスタートし、できるだけ式を省略せずに以下記載します。

$$p = {\frac{RT}{v-b} - \frac{\alpha(T)a}{v(v+b)}} \tag{4.1}$$

両辺に(v-b)v(v+b)をかけます。

$$p(v-b)v(v+b) =RTv(v+b) - \alpha(T)a(v-b) \tag{4.2}$$ $$pv^3-pb^2v-RTv^2-bRTv + \alpha(T)av-\alpha(T)ab =0 \tag{4.3}$$ $$pv^3-RTv^2+\left(-pb^2-bRT + \alpha(T)a\right)v-\alpha(T)ab =0\tag{4.4}$$

両辺に\(\frac{p^2}{R^3T^3}\)をかけます。さらに、\(\frac{pv}{RT}=z\)を使って、式変形します。

$$z^3-z^2+\left(\frac{-p^2b^2}{R^2T^2}-\frac{bp}{RT} + \frac{\alpha(T)ap}{R^2T^2}\right)z-\frac{\alpha(T)aP}{R^2T^2}\frac{bP}{RT} =0\tag{4.5}$$

ここで、\(\frac{\alpha(T)aP}{R^2T^2}=A\)、\(\frac{bp}{RT}=B\)とすると、

$$z^3-z^2+\left(-B^2-B+A\right)z-AB =0\tag{4.6}$$

大分スッキリした形になりました。A、Bはa,b,R,Tから求められますので、後はこの3次方程式を解くだけです。

3次方程式の解法については、非線形方程式の解法をご参照ください。

ゆくゆくVBAのページでも実例を掲載していこうと思います。

注意点

本サイトではSI単位系をベースとしております。

異なる単位系で計算される場合は単位換算を忘れないようにご注意ください。

注意点2

臨界物性や相互作用パラメーターについては当サイトには載せておりません。必要なデータについてはDIPPR等ご参照ください

アルファ関数について

主なアルファ関数を以下ご紹介します。

RK(1949)

$$\alpha (T)=\frac{1}{\sqrt{T}}$$

Wilson(1966)

$$\alpha (T)=\left[T_r+\left(1.57+1.62\omega \right)\left(1-T_r\right)\right]$$

SRK(1972)

$$\alpha (T)=\left[1+m\left(1-T_r^{0.5}\right)\right]$$

※mについては、こちらをご確認ください。

Soave(1993)

Lee-Kesler(1975)式で各Tr、ωで蒸気圧が等しくなるように係数を決めた式

※Lee-Kesler方のPitzer式は、NBP~超臨界域の推算精度は良いが、対臨界温度が低い領域での推算精度が著しく低下する、、、そうです。一度確認してみないと実感がわかないですが。

$$\alpha (T)=1+m\left(1-T_r\right)+n\left(1-T_r^{0.5}\right)^2$$ $$m=0.484+1.515\omega-0.44\omega^2$$ $$n=2.756m-0.700$$

Twu et al(1995):以下の論文によります

$$\alpha (T)=T_r^{N\left(M-1\right)}e^L\left(1-T_r^{NM}\right)$$

※α関数確認のソースは以下論文になります。(Part1はPeng-Robinson用)

Chorng H. Twu, John E. Coon, and John R. Cunningham; "A NEW GENERALIZED ALPHA FUNCTION FOR A CUBIC EQUATION OF STATE PART 2. REDLICH-KWONG EQUATION"; Fluid Phase Equilibria,1995, p.61-69

SRKの修正式としてTwuのα関数を用いた場合にどの程度精度が変わるのか、以下確認してみました。

上述の論文で、L,M,Nのパラメーターリストとして、Argon、Benzene、Cyclohexaneや主要なハイドロカーボンのデータがあります。

まずはmethaneについて確認してみます。

methaneのパラメーターとして次の論文中の数値を使用します。 Tc:190.58[K] 46.04[bar] L:0.106750[-] M:0.920161[-] N:3.09674[-]

準備中。

状態方程式法の各種混合則(Mixing Rule)について

SRK、PRの混合則は、非理想性の高い成分の組み合わせには適していません。

主な修正混合則を以下ご紹介します。

Panagiotopoulos-Reid Mixing Rule (1985)

$$a=\sum_i \sum_j y_i y_j a_{ij}$$ $$a_{ij}=\sqrt{a_i a_j}\left[1-k_{ij}+(k_{ij}-k_{ji})x_i\right]$$

Benedict-Webb-Rubin-Starling(BWRS)式

BWRS式(1972年)

Virial型の状態方程式の例として、BWRS式をご紹介します。

BWRS式は、軽炭化水素の物性推算に適した推算法で、重質炭化水素や、極性の強い物質の物性推算には適していません。

$$P = \rho RT+ \left( B_0RT-A_0-\frac{C_0}{T^2}+\frac{D_0}{T^3}-\frac{E_0}{T^4} \right) \rho ^2 + \left( bRT-a-\frac{d}{T} \right) \rho ^3 $$ $$+ \alpha \left( a + \frac{d}{T}\right) \rho^6 + \frac{c\rho^3}{T^2} \left(1+\gamma \rho^2 \right) exp \left(-\gamma \rho^2\right) \tag{5.1}$$

Cf)BWR式(1940年)

$$P = \rho RT+ \left( B_0RT-A_0-\frac{C_0}{T^2}\right) \rho ^2 + \left( bRT-a \right) \rho ^3 $$ $$+ a \alpha \rho^6 + \frac{c\rho^3}{T^2} \left(1+\gamma \rho^2 \right) exp \left(-\gamma \rho^2\right) \tag{5.2}$$

BWRS式における各パラメーター値が入手できない場合、偏心因子、臨界体積、臨界温度から求められます。

$$B_0=\frac{0.443690+0.115449\omega}{\rho_c}$$ $$A_0=\frac{\left(1.28438-0.920731\omega\right)RT_c}{\rho_c}$$ $$C_0=\frac{\left(0.356306+1.70871\omega\right)RT_c^3}{\rho_c}$$ $$D_0=\frac{\left(0.0307452+0.179433\omega\right)RT_c^4}{\rho_c}$$ $$E_0=\frac{\left(0.006450-0.022143\omega exp\left(-3.8\omega\right) \right)RT_c^5}{\rho_c}$$ $$b=\frac{0.528629-0.349261\omega}{\rho_c^2}$$ $$a=\frac{\left(0.484011+0.754130\omega \right)RT_c}{\rho_c^2}$$ $$d=\frac{\left(0.0732828+0.463492\omega \right)RT_c^2}{\rho_c^2}$$ $$c=\frac{\left(0.504087+1.32245\omega \right)RT_c^3}{\rho_c^2}$$ $$\alpha=\frac{0.0705233-0.044448\omega}{\rho_c^3}$$ $$\gamma=\frac{0.544979-0.270896\omega}{\rho_c^2}$$

※ただし、English単位系なので、注意。R=10.73[psi・ft3/(R・lb-mol)]

以下、ハロゲン化炭化水素の物性推算にも良いという記載があったので、試しに以下のペーパーに書かれているR11、R12、R13、R14、R22、R23、R113、R114について、BWRSでどの程度の精度が出るのか、REFPROPの計算結果と比較して確認してみました。

K. H. Aboul-Fotouh and K. E. Starling.; "USE OF A GENERALIZED MODIFIED BWR EQUATION OF STATE FOR HALOGENATED HYDROCARBON SATURATED THERMODYNAMIC PROPERTIES"; 1978,

詳細はこちらをご参照ください。

BWRS式は、SRKや、PRのように解析解を求めることができないので(数学に強い方なら解けるのかもしれませんが、、、)、収束計算で求めました。

BWRS式は、解から離れるとあっという間に発散してしまい、逐次代入法では全く収束しないため、Steffensen法で解いています。

ただ、Steffnesen法を用いても初期値が悪いとやはり収束しなかったので、やむをえず、一旦SRK式で密度を求め、その密度を初期値として与えて収束計算を行っています。

もっと上手い方法をご存知の方がいらっしゃいましたら、是非教えてください

使用したパラメーターは以下。

BWRSbase

検証結果(一部)は以下。

BWRScheck

REFPROPの適用範囲の制約もあり、暫定なので検証の圧力、温度がまばらですが、、、。

時間ができたらC#でBWRSの推算ツールも作成してみます。

相平衡の基本▼

準備中。

Wilson式▼

DECHEMAの式

$$ln\gamma_i=1-ln\sum_{j=1}x_j\lambda_{ij}-\sum_{k=1}\frac{\displaystyle x_k\lambda_{ki}}{\displaystyle \sum_{j=1}x_j\lambda_{kj}}$$ $$\lambda_{ij}=\frac{v_i^L}{v_j^L}exp\left[-\frac{\lambda_{ij}-\lambda_{ii}}{RT}\right]$$ $$A_{ij}=(\lambda_{ij}-\lambda_{ii})$$

ただし、R=1.98721[kcal/mol-K]

液液平衡計算は取り扱えない。

Wilson式の出典は以下

Excelを使った計算例はこちらをご参照ください。

NRTLVBA

Wilson,GM.; "Vapor-Liquid Equilibrium XI. A New Expression for the Excess Free Energy of Mixing"; J Amer Chem Soc, 1964, 86, p.127-130.

NRTL式▼

DECHEMAの式

$$ln\gamma_i=\frac{\displaystyle \sum_{j=1} \tau_{ji} G_{ji}x_j}{\displaystyle \sum_{k=1}G_{ki}x_k}+\sum_{j=1}\frac{\displaystyle x_j G_{ij}}{\displaystyle \sum_{k=1}G_{kj}x_k}\left(\tau_{ij}-\frac{\displaystyle \sum_{j=1} \tau_{kj}G_{kj}x_k}{\displaystyle \sum_{k=1}G_{kj}x_k}\right)$$ $$\tau_{ij}=\frac{g_{ij}-g_{jj}}{RT}$$ $$A_ij=(g_{ij}-g_{jj})$$

ただし、R=1.98721[kcal/mol-K]

上記の式形ですと、ある特定の温度におけるパラメーターで計算を行う場合はよいのですが、あるレンジ内で使えるパラメーターにしようと思うと、パラメーター数としては不十分な場合があります。

NRTLの場合、Simulatorrによって異なりますが、以下のような式形で表現されていることが多いかと思います。

$$\tau_{ij}=a_{ij}+\frac{b_{ij}}{T}+c_{ij}+d_{ij}ln(T)+e_{ij}T+f_ijTln(T)$$ $$G_ij=exp\left(-A{ji} \tau_{ij}\right)$$ $$A_{ji} = \alpha_{ji}+\beta_{ji}$$

この場合、DECHEMAの式の\(A_{ij}:[cal/mol]\)を、Simulatorの単位系に合わせる必要があります。

圧力一定の場合の計算例はこちらをご参照ください。

NRTL式の出典は以下

HRenon, H. and Prausnitz, J. M.; "Local Composition in Thermodynamic Excess Functions for Liquid Mixtures", AIChE J.,1968, vol.14, p.135-144.

Pythonを使った計算例はこちらをご参照ください。

NRTLplot

Excelを使った計算例はこちらをご参照ください。

NRTLVBA

Binary Parameterの回帰計算について▼

ExcelのSolverを使う方法、Pythonの最適化関数を使う方法をこちらに整理しています。

Vapor Fugacityの補正について▼

気相会合するような系ではVapor側のFugacityの補正が必要になります。

Hayden O'Conellの補正式をご紹介します。

活量係数型モデルの比較ツール▼

ツールでの比較

Wilson, NRTLの比較はこちらに計算ツールを作成しました。

※現状VLE(温度固定)のみです。

Wilsonのパラメーターは、以下の式形で搭載しています。

$$\lambda_{ij}=\frac{v_i^L}{v_j^L}exp\left[-\frac{\lambda_{ij}-\lambda_{ii}}{RT}\right]$$ $$A_{ij}=a_{ij}+\frac{b_{ij}}{T}$$

NRTLのパラメーターは、以下の式形で搭載しています。

$$\tau_{ij}=a_{ij}+\frac{b_{ij}}{T}+c_{ij}+d_{ij}ln(T)+e_{ij}T+f_ijTln(T)$$ $$G_ij=exp\left(-A{ji} \tau_{ij}\right)$$ $$A_{ji} = \alpha_{ji}+\beta_{ji}$$

ツールの使い方は、以下のように、計算温度、刻みを入力してCalcボタンを押すだけです。

※2016/10/8現在、各パラメーターのユーザー入力には対応しておりません。(水-ethanolのみ)

LACTtool

石油成分、単成分の物性推算▼

以下Twuの推算方法をご紹介します。

Chorng H. Twu ; "An internally consistent correlation for predicting the critical properties and molecular weights of petroleum and coal-tar liquids"  Fluid Phase Equilibria, 1984, 16, 137-150

Twuの推算法では、比重、標準沸点、分子量の内、2つ情報があれば、その他の物性(臨界物性、偏心因子、左記3つの物性の内の残り1つ)を推算することができます。

以下、主なハイドロカーボンに対して、比重、標準沸点から他の物性について、Excelで精度を確認してみます。

前準備

臨界温度

$$ T_c^0 = T_b * (0.533272 + 0.191017 \cdot 10 ^ {-3} \cdot Tb + 0.779681 \cdot 10 ^ {-7} \cdot Tb ^ 2 $$ $$ - 0.284376 \cdot 10 ^ {-10} \cdot Tb ^ 3 + 0.959468 \cdot 10 ^ {28} \cdot Tb ^ {-13}) ^ {-1}$$

臨界圧力

$$ P_C^O = (3.83354 + 1.19629 \cdot \alpha ^ {0.5} + 34.8888 \cdot \alpha + 36.1952 * \alpha ^ 2 + 104.139 \cdot \alpha ^ 4) ^ 2$$

ただし、\(\alpha=1-T_b/T_c^0 \)。

臨界体積

$$ V_C^O = (1 - (0.419869 - 0.505839 \cdot \alpha - 1.56436 \cdot \alpha ^ 3 - 9481.7 \cdot \alpha ^ {14})) ^ {-8}$$

比重

$$ SG^O = 0.843593 - 0.128624 \cdot \alpha - 3.36159 \cdot \alpha ^ 3 - 13749.5 \cdot \alpha ^ {12}$$

分子量

$$T_b = \exp{(5.71419 + 2.71579 \cdot X - 0.28659 \cdot (\theta ^ 2) - 39.8544 / \theta - 0.122488 / (\theta ^ 2))} $$ $$ - 24.7522 \cdot \theta + 35.3155 \cdot (\theta ^ 2)$$

\(\theta\)を求めるには収束計算が必要ですので、ここではVBAでNewton法を使って求めます。

$$ \theta = ln (MW^0)$$ 前準備で求めたパラメーターから各物性の推算

臨界温度

$$ T_c = T_c^0 \cdot ((1 + 2 \cdot f_T) / (1 - 2 \cdot f_T)) ^ 2$$ $$ f_T = \Delta SG_T * (-0.362456 / (T_b ^ 0.5) + (0.0398285 - 0.948125 / (T_b ^ {0.5})) \cdot \Delta SG_T)$$ $$ \Delta SG_T = \exp(5(SG^0-SG))-1 $$

臨界圧力

$$ P_c = P_C^O \cdot (T_c / T_c^0) \cdot (V_C^O / V_c) \cdot ((1 + 2 \cdot f_p) / (1 - 2 \cdot f_p)) ^ 2 $$ $$ f_p = \Delta SG_P \cdot (2.53262 - 46.1955 / (T_b ^ {0.5}) - 0.00127885 * T_b + (-11.4277 + 252.14 / T_b ^ {0.5} $$ $$+ 0.00230535 * T_b) * \Delta SG_P) $$

臨界体積

$$ V_c = V_C^O * ((1 + 2 \cdot f_v) / (1 - 2 \cdot f_v)) ^ 2 $$ $$ f_v = \Delta SG_V \cdot (0.46659 / T_b ^ {0.5} + (-0.182421 + 3.01721 / T_b ^ {0.5}) \cdot \Delta SG_V) $$

分子量

$$ ln (MW) = Log(MW_0) \cdot ((1 + 2 \cdot f_M) / (1 - 2 \cdot f_M)) ^ 2 $$ $$ f_M = \Delta SG_M * (| x | + (-0.0175691 + 0.193168 / T_b ^ {0.5}) * \Delta SG_M) $$ $$ |x| = | 0.0123420 - 0.328086 / T_b ^ {0.5} | $$ $$ \Delta SG_M = \exp(5(SG^0-SG))-1 $$

MWの\(\theta\)の計算については、ブログを確認してください。エクセルでn-アルカンの物性を、Perryのengineering bookに記載されているSG、NBPから推算し、結果をやはりPerryの同著の数値と比較した結果が以下になります。

Twu

分子量の大きい成分の臨界体積と、一部臨界圧力のズレが大きくなっているものの、その他の物性については、誤差1%以下で精度としは問題無さそうです。

蒸気の推算は、比重、標準沸点からその他の物性を推算しましたが、もし入力した比重、標準沸点に誤差があった場合、最終的な計算結果にどの程度の影響を与えるでしょうか。

以下、比重、標準沸点にわざと数%の誤差を与えて、計算結果の誤差を確認してみます。

準備中。

簡単ですが、上記の推算ツールを作成しました。(2016/08/27)

使い方は以下のようにNBP(K)、SGを入力し、Calculateを押すだけです。

PetroTool

蒸留曲線からの物性推算▼

準備中。

石油成分、多成分の物性推算▼

準備中。

当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。

Page Top