物性推算の基本について

状態方程式型 / 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式の推算精度の確認については、こちら(Internet Explorerでのみご覧いただけます)をご参照ください。

上記のリンク先ツールで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年に発表している補正式も合わせてご紹介します。

Claude F.Leibovici, ; "Variant and invariant properties from cubic equations of state; Fluid Phase Equilibria, 1993, Vol.84, no. 1, p.1-8

SRKのエンタルピー算出

理想気体のエンタルピー

理想気体のエンタルピーは温度の関数で表現されます。

PerryのhandbookNIST Webbookより、理想気体の比熱を求めるためのパラメーターが掲載されているので、積分することでエンタルピーを求めることができます。

  

積分してエンタルピーを求めるためには定数項を決める必要があります。定数項はシミュレーターによって異なりますが、例えばAVEVA社のPRO/II®、DYNSIM®では、0Cにおける液エンタルピーが0となるように決められています。(理想気体のエンタルピーは0Cにおけるエンタルピーが潜熱となります)※NIST Webbookではエンタルピーも記載されています

  

Ex.PRO/II®の理想気体エンタルピー係数の確認画面(係数は消しています。温度範囲外になると、外挿されるため注意が必要です)

GoaslSeek GoaslSeek

実在気体のエンタルピー

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

なお、偏倚の式導出は以下の書籍が参考になります。

Integrated Design and Simulation of Chemical Processes

$$\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のフガシティー計算

成分iのフガシティーは以下の式より求めることができます。

$$ln \Phi_i = \frac{1}{RT}\int_V^\infty\left[ \left( \frac{\partial P}{\partial n_i} \right)_{T,V,n_j} - \frac{RT}{V}\right]dV-ln(z)$$

SRKの純成分蒸気圧算出

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

平衡定数の初期値として、よく以下のWilson式が使われます。(臨界温度:Tc、臨界圧力:Pc、標準沸点:NBPが必要です)

準備中。

Wilson式の出典は以下

Wilson, G., “A modified Redlich-Kwong equation of state applicable to general physical data calculations”; 65th AIChE National meeting,1968.

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のページでも実例を掲載していこうと思います。

アルファ関数について

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

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式の出典は以下

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.

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

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}+\frac{c_{ij}}{T^2}+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

代表的な原子団寄与法

原子団寄与法による推算について:原子団寄与法による各種物性(標準沸点、臨界物性等)の求め方、及び比較結果をまとめていきます。

以下の本が参考になります。

Jobakの推算法

Lydersenの原子団寄与法に基づいた推算法

標準沸点、臨界温度、臨界圧力の推算

$$T_{b}(K) = 198 + \Sigma (vck)$$ $$T_{c}(K) = T_{b} \left[ 0.584 + 0.965{\sum_{k} N_{k}(tck)} - {\sum_{k} N_{k}(tck)}^2 \right]^{-1}$$ $$P_{c}(bar) = \left[ 0.113 + 0.0032 N_{atoms} - {\sum_{k} N_{k}(pck)} \right]^{-2}$$

Wilson-Jaspersonの推算法

標準沸点、臨界温度、臨界圧力の推算

$$T_{b}(K) = 204.359 ln \left[ \sum_{k} N_{k}(tblk) + W \sum_{j} M_{j}(tb2j) \right]$$ $$T_{c}(K) = T_{b}\left[0.048721 - 0.019846N_r + {\sum_{k} N_{k}(\Delta tck)} + \sum_{j} M_{j}(\Delta tcj)\right]^{0.2}$$ $$P_{c}(bar) = 0.0186233 \cdot T_c / \left[ -0.96601 + exp(Y) \right] $$ $$Y = -0.0092295-0.0290403 \cdot N_r + 0.041 \left[ \sum_{k} N_{k}(\Delta pck) + \sum_{j} M_{j}(pcj) \right]$$

Marrero-Pradilloの推算法

標準沸点、臨界温度、臨界圧力の推算

$$exp(T_b/T_{b0}) = \sum_{i} N_{i}T_{b1i} + \sum_{j} M_{j}T_{b2j} + \sum_{k} M_{k}T_{b3k} $$ $$exp(T_c/T_{c0}) = \sum_{i} N_{i}T_{c1i} + \sum_{j} M_{j}T_{c2j} + \sum_{k} M_{k}T_{c3k} $$ $$(P_c-P_{c1})^{-0.5} = \sum_{i} N_{i}P_{c1i} + \sum_{j} M_{j}P_{c2j} + \sum_{k} M_{k}P_{c3k} $$

出典

J.Marrero,R.Gani,"Group-contribution based estimation of pure component properties",Fluid Phase Equilibria,2001, p.183–p.208

Constantinou-Ganiの推算法

標準沸点、臨界温度、臨界圧力の推算

$$T_{b}(K) = 204.359 ln \left[ \sum_{k} N_{k}(tblk) + W \sum_{j} M_{j}(tb2j) \right]$$ $$T_{c}(K) = 181.128 ln \left[ {\sum_{k} N_{k}(tclk)} + W \sum_{j} M_{j}(tc2j)\right]$$ $$P_{c}(bar) = \left[ {\sum_{k} N_{k}(pclk)} + W \sum_{j} M_{j}(pc2j) + 0.10022 \right]^{-2} + 1.3705 $$

各推算法の精度

構造異性体を考慮できるのはMarrero-Pradilloの推算法だけですが、以下各推算法について、o-xylenの各物性について、DIPPRのデータと比較してみます。

Jobakの推算法

標準沸点の推算用係数

Nk tbk Nk(tbk)
CH 4 26.73 106.92
CH3 2 23.58 47.16
C 2 31.01 62.02

臨界温度、臨界圧力の推算用係数

Nk tck pck Nk(pck) Nk(pck)
CH 4 0.0082 0.0011 0.0328 0.0044
CH3 2 0.0141 -0.0012 0.0282 -0.0024
C 2 0.0143 0.0008 0.0286 0.0016

推算結果

Jobak DIPPR Error %
Tb[K] 414.10 417.58 0.84
Tc[K] 625.1 630.3 0.83
Pc[K] 35.86 37.32 4.08
Wilson-Jaspersonの推算法

標準沸点の推算用係数

Nk tbk Nk(tbk)
ACH 4 0.9297 3.7188
ACCH3 2 1.9669 3.9338

臨界温度、臨界圧力の推算用係数

Nk tck pck Nk(pck) Nk(pck)
C 8 0.008532 0.72983 0.068256 5.83864
H 10 0.002793 0.1266 0.02793 1.266

推算結果

Wilson DIPPR Error %
Tb[K] 415.90 417.58 -0.40
Tc[K] 630.8 630.3 0.08
Pc[K] 36.49 37.32 -2.27
Constantinou-Ganiの推算法

標準沸点、臨界温度、臨界圧力の推算

標準沸点の推算用係数

Nk tbk Nk(tbk)
ACH 4 0.9297 3.7188
ACCH3 2 1.9669 3.9338

臨界温度、臨界圧力の推算用係数

Nk tck pck Nk(pck) Nk(pck)
ACH 4 3.7337 0.0075 14.9348 0.03
ACCH3 2 8.213 0.0194 16.426 0.0388

推算結果

CG DIPPR Error %
Tb[K] 415.90 417.58 0.41
Tc[K] 624.1 630.3 1.00
Pc[K] 36.37 37.32 2.60
Marrero-Pradilloの推算法

標準沸点の推算用係数

Nk tbk Nk(tbk)
ACH 4 0.8365 3.346
ACCH3 2 1.5653 3.1306

臨界温度、臨界圧力の推算用係数

Marrero-Pradilloの推算法では、構造異性体も考慮することが可能です。結合箇所に応じて第2項までパラメーターが用意されています。以下のテーブルは第1項、第2項の計算結果を積算した値です。

Nk tck pck Nk(pck) Nk(pck)
aCH 4 2.0337 0.00726 8.1348 0.02904
aC-CH3 2 3.4611 0.020907 6.9222 0.041814

推算結果

MP DIPPR Error %
Tb[K] 410.22 417.58 1.76
Tc[K] 622.2 630.3 1.29
Pc[K] 36.72 37.32 1.61

代表的な輸送物性の推算

Chapman-Enskogの式を用いた物性推算法について、以下にまとめています。

以下粘度、熱電伝度、拡散係数の推算において、主にChapman-Enskogの式を用いますが、その際Lennard-Jonesパラメーターが必要となります。パラメーターが得られない場合の推算式を以下紹介します。

なお、The Properties of Gases and Liquidsでは注意点として、極性物質についてはLennard-Jonesパラメーターのみでは不十分であることが述べられており、Brokawの補正式が紹介されていますので、極性物質の推算を行う場合は、じ。

$$ \sigma \left(\frac{p_c}{T_c} \right)^{1/3} = E-10(2.3551-0.087) $$ $$ \frac{\epsilon / k_B}{T_c} = 0.7915 + 0.16930 $$

ここで、\(p_c\):臨界圧力[atm]、\(p_c\):臨界温度[K]、\( \omega \):偏心因子[-]、\(\sigma\):衝突直径[m]です。

上式の参考書籍は以下になります。なお、以下の参考書籍のAppendixに主な気体のLennard-Jonesパラメーターが記載されています。

例えば上記の書籍に記載されているH2のパラメーターについて精度を確認してみます。

NISTのWebbookから、Tc : 33.18[K]、Pc[bar]を使って\(\sigma\)を求めると、3.113e-10。

一方で上記書籍の数値は\(\sigma\)=2.827e-10ですので誤差10%程度です。

また、\(\epsilon / k_B \)については、推測値:31.9に対し、上記書籍記載の数値は59.7ですので、こちらは90%近い誤差が生じています。

Lennard-Jonesパラメーターについては極力文献値等を探した方が良さそうです。

上述の参考書籍には以下4つの前提条件が記載されています。以降のChapman-Enskogの式を採用する場合場合、念頭に置いてください。

1.気体が十分に希薄で、2分子間の衝突のみが想定できること

2.衝突時の分子運動は古典力学で表現できること

3.全ての衝突が弾性衝突であること

4.分子間の相互作用は分子間の中心のみで考えられること

低圧域の粘度推算:Chapman-Enskogの式(1916)

$$ \mu = 2.6693E^{-26} \left( \frac{(MT)^{1/2}}{\sigma^2 \Omega^*} \right)$$

ここで、\( \mu \):粘度[Pa/sec]、M:分子量、\(\sigma\):衝突直径[m]、 \(\Omega^*\):還元衝突積分(reduced collision integral)

ここで、還元衝突積分は無次元温度:\(T^* = T/ (\epsilon / k_B) \)の関数となっている。ただし、\(k_B\):ボルツマン定数(1.3806504×E-23[J/K])、\( \epsilon \):特徴的な結合エネルギー(characteristic energy of interaction)。無次元温度の適用範囲は、0.3~100。

例題14.1で空気の粘度が計算されており、313.15Kの1点データとはいえ、誤差0.4%と精度が良いことが示されています。

衝突積分は上記の参考書籍のAppendixより、下式で表されます。

$$\Omega^* = \frac{1.16145}{(T^*)^{0.14874}} + \frac{0.52487}{exp(0.77320 \cdot T^*)} + \frac{2.16178}{exp(2.43787 \cdot T^*)}$$

また、0.4 <\( T^* \) <1.4においては、以下の簡易式も紹介されています。

$$\Omega^* = \frac{1.1604}{(T^*)^{0.5}} $$

H2しか確認していませんが、無次元温度が上記の範囲内であれば十分な精度がありそうです。

参考

Chung TH, Ajlan M, Lee LL, Starling KE, "Generalized multiparameter correlation for nonpolar and polar fluid transport properties.", Ind Eng Chem, 1988, 27, 671-9.

粘度の混合測

熱伝導度推算

Chapman-Enskogの式

$$ k = 8.3224E^{-22} \left( \frac{(T/M)^{1/2}}{\sigma^2 \Omega^*} \right)$$

ここで、熱電導度:k[W/m-K]。その他のパラメーターについては粘度の項をご参照ください。

熱伝導度の混合測

相互拡散係数推算

$$ D = D_{AB} = D_{BA} = 1.8829E^{-22} \frac{ \{T^3 \left[ (1/M_A) + (1/M_B) \right] \}^{1/2} }{p \sigma_{AB}^2 \Omega_D}$$

ここで、衝突積分は上記の参考書籍のAppendixより、下式で表されます。

$$\Omega^* = \frac{1.06036}{(T^*)^{0.15610}} + \frac{0.19300}{exp(0.47635 \cdot T^*)} + \frac{1.03587}{exp(1.529 \cdot T^*)} + \frac{1.76474}{exp(3.89411 \cdot T^*)}$$

拡散係数については、上記の‘Chapman-Enskogの式以外にも、例えばThe Properties of Gases and Liquidsでは、Wilke-LeeやFullerの方法が紹介されています。

こちらの本では、上述の3つの方法の精度について、CO2、N2の実験結果との比較を行っていますが、いずれの方法も大きな差異はなく、いずれもそれなりの精度がありそうです。(ただ、注意点として、あくまで上記の例はデータバンクのある成分同士の比較であって、アミンではどの程度の精度が出るか分からないと記載されています。)

液相の拡散係数について、The Properties of Gases and LiquidsではWile-ChangTyn-CalusHayduk-MinhasNakanishiの4つの式が紹介されています。

$$ D_{AB} = \frac{7.4e^{-8}(\phi M_B)^{1/2}T}{\eta_B V_A^{0.6}} $$

ここで、\( D_{AB}\)はB溶液中における溶質Aの相互拡散係数[cm/s]、\(M_B\)は溶液Bの分子量[g/mol]、\( \eta_B \)は溶液Bの粘度[cP]、\(V_A\)は溶質Aの標準沸点における分子容量[cm3/mol]、\(\phi \)は溶液Bの会合係数になります。

会合係数は、溶液が水の場合2.6、メタノールで1.9、エタノールで1.5、会合を起こさない場合は1.0です。

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

以下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:Internet Explorerのみ対応)

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

PetroTool

蒸留曲線からの物性推算

準備中。

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

準備中。

注意点1

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

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

注意点2

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

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

Page Top