ギブスの自由エネルギーと圧力pのグラフ

ファンデルワールスの状態方程式

上記をP-V線図にプロットすると以下になる。

ギブスの自由エネルギーgA-gMは体積Vを用いて以下のように表される。

いま、ギブスの自由エネルギーの最下点（極小点）gMのエネルギーを仮にゼロとすると（任意の定数なのでゼロとおいても差し支えない）、左辺gA-gMはgAとなる。
この方程式(7.12)式を、圧力pが与えられたときにVについて解くことを考える。
p>pNおよびp

clc; clear; close all;
a = 3.59;
b= 0.0427;
R = 0.08206;
T = 300;
v0 = 0.1;
vf = 0.2;
v = linspace(v0,vf,1000);%m3
p  = R*T./(v-b) - a./v.^2;
semilogx(v,p,'b');
ylim([68 72]);
xlabel('v[L]');
ylabel('p[atm]')
title('ファンデルワールスの状態方程式');
hold on;

TF = islocalmin(p);
plot(v(TF),p(TF),'r*');
minV = v(find(TF==1));
minP = p(find(TF==1));
disp(sprintf('Minimum: (v,p)=(%f,%f)',minV,minP));

TF = islocalmax(p);
plot(v(TF),p(TF),'r*');
maxV = v(find(TF==1));
maxP = p(find(TF==1));
disp(sprintf('Maximum: (v,p)=(%f,%f)',maxV,maxP));

big;

PM = 69.1;
PN = 69.8;
Ps = linspace(PM,PN,20);
N = length(Ps);
vss = zeros(N,3);
syms v;
g = -2*a/v - R*T*log(v-b) + R*T*b/(v-b);
g1s = [];
g2s = [];
g3s = [];
for idx=1:N
    p = Ps(idx);
    eqn = R*T./(v-b) - a./v.^2 - p == 0;    
    vs = vpasolve(eqn);%解が3つ求まる
    vss(idx,:)  = vs;
    v1 = vs(1);
    g1 = subs(g,v1);
    g1s(end+1) = g1;
    v2 = vs(2);
    g2 = subs(g,v2);
    g2s(end+1) = g2;
    v3 = vs(3);
    g3 = subs(g,v3);
    g3s(end+1) = g3;
end

figure();
plot(Ps,g1s,'k');
hold on;
plot(Ps,g2s,'k');
plot(Ps,g3s,'k');
extension = 0.3;
PM_lower = linspace(PM-extension,PM,10);
gs_lower = [];
offset = 4;
for idx=1:length(PM_lower)-offset
    p = PM_lower(idx);
    eqn = R*T./(v-b) - a./v.^2 - p == 0;    
    v_ans = vpasolve(eqn);%解が1つ求まる
    v_ans = v_ans(imag(v_ans)==0);
    g_single = subs(g,v_ans);
    gs_lower(end+1) = g_single;
end
plot([PM_lower(1:end-offset),Ps(1)],[gs_lower,g3s(1)],'k');

PN_upper = linspace(PN,PN+extension,10);
gs_upper = []
for idx=offset:length(PN_upper)
    p = PN_upper(idx);
    eqn = R*T./(v-b) - a./v.^2 - p == 0;    
    v_ans = vpasolve(eqn);
    v_ans = v_ans(imag(v_ans)==0);
    g_single = subs(g,v_ans);
    gs_upper(end+1) = g_single;
end

plot([Ps(end),PN_upper(offset:end)],[g1s(end),gs_upper],'k');
ylabel("Gibbs free energy");
xlabel('p[atm]');
xlim([68.7,70.2]);
title(['ファンデルワールスの状態方程式の',newline,'体積の3重解区間における',newline,'ギブスの自由エネルギー']);

big;

できた

上記のグラフは以下の図に一致する