van der Waals interaction and figures in Schroeder Sec 5.3

Remove["Global`*"]

Note that I use “Nt” instead of “N” because “N” is a protected name in Mathematica

$Assumptions={P>0,V>0,T>0,Nt>0,k>0,a>0,b>0};

Find the critical point from the equation of state

Write down the van der Waals equation of state and solve for the pressure

vdw=(P+a Nt^2/V^2)(V-Nt b)==Nt k T;
solP=Solve[vdw,P];
pressure=P/.solP[[1]];

Solve for the first and second derivatives of P to be zero

sol=Solve[{D[pressure,V]==0,D[pressure,{V,2}]==0}, {V,T}];
Vc=V/.sol[[1]]
Tc=T/.sol[[1]]
Pc=pressure/.sol[[1]]

Derive the scaled version of the equation of state

Simplify[vdw/.{P->p Pc,V->v Vc,T->t Tc}]

Plot Figure 5.20

p=Simplify[(pressure/Pc)/.{V->v Vc,T->t Tc}]

Plot[p/.t->{0.8,0.9,1.0,1.1,1.2},{v,0.4,4},
PlotRange->{{0,4},{-0.5,2.5}},PlotStyle->Black]

Plot the right side of Fig 5.21

Define the PV curve for T=0.9Tc

p0p9=p/.t->0.9;

Find the point labeled 1

p1=0.5;
v1=v/.NSolve[p0p9==p1][[1]];

Find the points labeled 2,4,6. This is just a guess! I am estimating the value of p where the flat line is drawn in the figure. In fact, this is supposed to be the “vapor pressure” where the phase transition takes place. The method to find the exact value is described later in the section, and is carried through at the end of this notebook.

pFlat=0.6;
solFlat=NSolve[p0p9==pFlat,v];

p2=pFlat;
v2=v/.solFlat[[1]];
p4=pFlat;
v4=v/.solFlat[[2]];
p6=pFlat;
v6=v/.solFlat[[3]];

Find the points labeled 3 and 5

solExtr=NSolve[D[p0p9,v]==0,v];
v3=v/.solExtr[[1]];
p3=p0p9/.v->v3;
v5=v/.solExtr[[2]];
p5=p0p9/.v->v5;

Pick a point for 7 that is just smaller than 6

v7=0.94 v6;
p7=p0p9/.v->v7;

Make the seven points into a list and display them

vs={v1,v2,v3,v4,v5,v6,v7}
ps={p1,p2,p3,p4,p5,p6,p7}

Now plot the right side of Figure 5.21

Show[
Plot[p0p9,{v,0.4,4},
PlotRange->{{0,4.0},{0.3,1.0}},PlotStyle->Black],ListPlot[Table[Callout[{vs[[i]],ps[[i]]},i],{i,1,7}],
PlotMarkers->{Automatic,10},PlotStyle->Black]
]

Determine the Gibbs Free Energy and plot the left side of Fig 5.21

Follow (5.54) and (5.55) to get (5.56)

dGdV=V D[pressure,V];
G=Integrate[dGdV,V]

Simplify it down to get rid of the constants

FullSimplify[G/(Nt k T)/.{T->t Tc,V->v Vc}]

g=Simplify[%+Log[b Nt]]

Make the plot in Fig 5.21

g0p9=g/.t->0.9;
gs=g0p9/.v->vs;ListLinePlot[Table[Callout[{ps[[i]],gs[[i]]},i],{i,1,7}],
PlotRange->{{0.4,0.9},{-2.9,-2.6}},
PlotMarkers->{Automatic,10},PlotStyle->Black]

Points 2 and 6 are not at the same spot, because we just guessed at a value for the vapor pressure. Now let’s try to get it right.

Calculation of the vapor pressure for T=0.9 Tc

Above we just guessed the value pFlat for the vapor pressure. The fact that our points 2 and 6 aren’t at the same point in G vs P means that we didn’t guess quite right. So let’s do it right, using the Maxwell Construction and then check it using the plot of G vs P.

Calculate the integral of pressure vs volume (in scaled variables)

p0p9Int=Integrate[p0p9,v]

Use this function to determine the integral of the “lobes” for different choices of pressure

vp=Range[0.5,0.725,0.01];
lvp=Length[vp];
vMin=Table[Min[v/.NSolve[p0p9==vp[[i]]]],{i,1,lvp}];
vMax=Table[Max[v/.NSolve[p0p9==vp[[i]]]],{i,1,lvp}];
pint=Table[
Re[(p0p9Int/.v->vMax[[i]])-(p0p9Int/.v->vMin[[i]])
-(vMax[[i]]-vMin[[i]])vp[[i]]],{i,1,lvp}];
points=Transpose[{vp,pint}];
ListLinePlot[points,PlotRange->{{0.45,0.75},{-0.5,0.5}}]

The sum of the lobes is zero near 0.65, not the 0.60 we originally chose. We can get very close to the right vapor pressure by making an interpolation function from the points and finding where it crosses zero.

ifun = Interpolation[points];
vpVal=x/.FindRoot[ifun[x],{x,0.65}]

Draw the curve that is G vs P, and show that it looks like Fig 5.21, but show the vapor pressure

Show[
ParametricPlot[{p0p9,g0p9},{v,0.575,3.5},
PlotRange->{{0.4,0.875},{-2.875,-2.625}},PlotStyle->Black],
ListLinePlot[{{vpVal,-2.85},{vpVal,-2.685}},PlotStyle->Red]
]

It looks like we got it right. Finally, plot Fig 5.21 again but with the right choice

solFlat=NSolve[p0p9==vpVal,v];

p2=vpVal;
v2=v/.solFlat[[1]];
p4=vpVal;
v4=v/.solFlat[[2]];
p6=vpVal;
v6=v/.solFlat[[3]];

vs={v1,v2,v3,v4,v5,v6,v7};
ps={p1,p2,p3,p4,p5,p6,p7};

Show[
Plot[p0p9,{v,0.4,4},
PlotRange->{{0,4.0},{0.3,1.0}},PlotStyle->Black],ListPlot[Table[Callout[{vs[[i]],ps[[i]]},i],{i,1,7}],
PlotMarkers->{Automatic,10},PlotStyle->Black]
]

gs=g0p9/.v->vs;ListLinePlot[Table[Callout[{ps[[i]],gs[[i]]},i],{i,1,7}],
PlotRange->{{0.4,0.9},{-2.9,-2.6}},
PlotMarkers->{Automatic,10},PlotStyle->Black]