vdW数値計算(蒸気圧曲線)

むかーしむかしのノートを記します.

超不完全なので随時更新するかもしれないし、二度と更新しないかもしれません.

---

vdW状態方程式

(p+\dfrac{an^2}{V^2})(V-nb)=nRT

p=\dfrac{nRT}{V-nb}-\dfrac{an^2}{V^2}

Tを変えると臨界点C(p_c,V_c,T_c) がある.

f:id:Atsumu:20200611141022p:plain

T=T_cの変曲点で(p_c,V_c,T_c)を決めると、

\dfrac{\partial p}{\partial V}=\dfrac{-nRT}{(V-nb)^2}+\dfrac{2an^2}{V^3}

\dfrac{-nRT_c}{(V_c-nb)^2}+\dfrac{2an^2}{V_c^3}=0

\dfrac{nRT_c}{(V_c-nb)^2}=\dfrac{2an^2}{V_c^3}...(1)

 

\dfrac{\partial^2 p}{\partial V^2}=\dfrac{2nRT}{(V-nb)^3}-\dfrac{6an^2}{V^4}

\dfrac{\partial^2 p}{\partial V^2}=0

\dfrac{2nRT_c}{(V_c-nb)^3}=\dfrac{6an^2}{V_c^4}...(2)

(1),(2)の辺々割って整理すると、V_c=3nb...(3)

(1)に代入して整理すると、T_c=\dfrac{8a}{27Rb}...(4)

 V_c,T_cをvdW状態方程式に代入して、p_c=\dfrac{8a}{27Rb}...(5)

以上より、(p_c,V_c,T_c)=(\dfrac{8a}{27Rb},3nb,\dfrac{8a}{27Rb}) 

臨界点C を基準に各状態量を表すと、

p=p_c \hat{p},V=V_c \hat{V},T=T_c \hat{T}

状態方程式に代入して、

 (p_c \hat{p}+\dfrac{an^2}{(V_c \hat{V})^2})(V_c \hat{V}-nb)=nR\dfrac{8a}{27bR}\hat{T}

(3)~(5)を 代入して、整理すると無次元化したvdW状態方程式

\hat{p}\hat{V}=\dfrac{8\hat{T}\hat{V}}{3\hat{V}-1}-\dfrac{3}{\hat{V}}

この圧力の表式を用いて計算してみる.

---

化学ポテンシャルの導出は長いので、また現実逃避したいときに書きます.

結果だけ書くと、

新たな基準点O(p_o,V_o,T_o)(V_o\gg nb)を設けて、

z=\dfrac{3}{2},T_o=V_o=1,S_o=0の条件下で、(多分簡単のため.意味があるかもしれないので、いつか考えます.)

---

μ=z\hat{t}(1-\ln{\hat{t}})-\hat{t}\ln{\dfrac{3\hat{V}-1}{3}}+\dfrac{3\hat{V}\hat{T}}{3\hat{V}-1}-\dfrac{9}{4\hat{V}}

 

これらをmathematicaに定義すると、

f:id:Atsumu:20200609172145p:plain

 試しに、t=0.9(一定)で、V=\dfrac{1}{3}\rightarrow10まで書き、

FindMaximumで極大点、FindMinimumで極小点を求めます.

ここで求めた体積V_M,V_mは次の式を満たすVを探すヒントに使います.

p(V_1)=p(V_2)かつμ(V_1)=μ(V_2)が成り立つ(V_1,V_2)は、上で求めたV_M,V_mに対して、V_1\leq V_m,V_M\leq V_1にあります.

FindRootでp(V_1)=p(V_2),μ(V_1)=μ(V_2)連立方程式を解いて、(V_l,V_g)=(V_1,V_2)と置いておきます.

次に(V_l,V_g)のときのp,μをそれぞれ計算して、

最後の行で、数値をリストにまとめます.

f:id:Atsumu:20200609172325p:plain

 そんなことを延々と続けて、大元のリストを作ります.

f:id:Atsumu:20200609172353p:plain

 次に、大元のリストから{T,P}のリストを作り、Interpolationでプロット間を補完します.

蒸気圧曲線です.(T,P)=(1,1)で臨界点です.

f:id:Atsumu:20200609172429p:plain

ちなみに、クラウジウス-クラペイロンの式から

---

導出いつかやるかも.

---

P(T)=P_o\exp{(-\dfrac{Q}{RT})}です.