Pressure-Temperature Phase Diagram for Water

clausius[T_]:=0.101325*Exp[-5268.134*((1/T)-(1/373))]

ice={{252.,203.5357},{253.,195.8644},{254.,188.1245},{255.,180.3019},{256.,172.3814},{257.,164.3468},{258.,156.1811},{259.,147.866},{260.,139.3821},{261.,130.7087},{262.,121.8238},{263.,112.7041},{264.,103.3243},{265.,93.6579},{266.,83.6765},{267.,73.3496},{268.,62.645},{269.,51.5281},{270.,39.962},{271.,27.90750},{272.,15.32260},{273.,2.1627},{273.1,1.},{273.16,0.0006117}};

fice=Quiet@Interpolation[ice,InterpolationOrder->1];

Manipulate[
Module[{T,P,Pc,Tc,R,data,pressure,diagram,TempAtHeat,Ti,cpIce,cpWater,cpVapor,heatSublime,heatMelt,heatVapor,heatTriple,Hm,Hv,Hv0,Hs,νLiq,νIce,νVap,title,comp,phases},
Module[{g00,g01,g02,g03,g04,t1,t2,r1,r20,r21,r22,T2,π0,πp,τ,r2p,g0p,g,ρ,Tt},
g00=-0.632020233335886*10^6;
g01=0.655022213658955;
g02=-0.189369929326131*^-7;
g03=0.339746123271053*^-14;
g04=-0.556464869058991*^-21;
t1=0.368017112855051*^-1+(0.510878114959572*^-1)I;
t2=0.337315741065416+0.335449415919309*I;
r1=0.447050716285388*^2+(0.656876847463481*^2)I;
r20=-(0.725974574329220*^2)-(0.781008427112870*^2)I;
r21=-(0.557107698030123*^-4)+(0.464578634580806*^-4)I;
r22=(0.234801409215913*^-10)-(0.285651142904972*^-10)I;
Tt=273.16;
π0=(101325/611.6571);
πp[p_]:=(p/611.6571);
τ[T_]:=(T/Tt);
r2p[p_]:=r21*(1/611.6571)(((πp[p]-π0))^0)+r22*(2/611.6571)(((πp[p]-π0))^1);
g0p[p_]:=g01*(1/611.6571)*(((πp[p]-π0))^0)+g02*(2/611.6571)*(((πp[p]-π0))^1)+g03*(3/611.6571)*(((πp[p]-π0))^2)+g04*(4/611.6571)*(((πp[p]-π0))^3);
g[T_,p_]:=g0p[p]+Tt*Re[r2p[p]*((t2-τ[T])*Log[t2-τ[T]]+(t2+τ[T])*Log[t2+τ[T]]-2*t2*Log[t2]-(((τ[T])^2)/t2))](*specific volume of ice, (m^3)/kg*);
ρ[T_,p_]:=(((g[T,p]*.001*1*(10^6)))^(-1))(*density of ice, g/(cm^3)*);
νIce[T_,p_]:=g[T,p]*1000*0.01801528(*molar specific volume of ice, L/mol*)];
Module[{IJ,R2,πp,τ,γπ,νmassLiq},
IJ={{1,0,-2,0.14632971213167},{2,0,-1,-0.84548187169114},{3,0,0,-3.7563603672040},{4,0,1,3.3855169168385},{5,0,2,-0.95791963387872},{6,0,3,0.15772038513228},{7,0,4,-0.016616417199501},{8,0,5,0.81214629983568*^-3},{9,1,-9,0.28319080123804*^-3},{10,1,-7,-0.60706301565874*^-3},{11,1,-1,-0.18990068218419*^-1},{12,1,0,-0.32529748770505*^-1},{13,1,1,-0.21841717175414*^-1},{14,1,3,-0.52838357969930*^-4},{15,2,-3,-0.47184321073267*^-3},{16,2,0,-0.30001780793026*^-3},{17,2,1,0.47661393906987*^-4},{18,2,3,-0.44141845330846*^-5},{19,2,17,-0.72694996297594*^-15},{20,3,-4,-0.31679644845054*^-4},{21,3,0,-0.28270797985312*^-5},{22,3,6,-0.85205128120103*^-9},{23,4,-5,-0.22425281908000*^-5},{24,4,-2,-0.65171222895601*^-6},{25,4,10,-0.14341729937924*^-12},{26,5,-8,-0.40516996860117*^-6},{27,8,-11,-0.12734301741641*^-8},{28,8,-6,-0.17424871230634*^-9},{29,21,-29,-0.68762131295531*^-18},{30,23,-31,0.14478307828521*^-19},{31,29,-38,0.26335781662795*^-22},{32,30,-39,-0.11947622640071*^-22},{33,31,-40,0.18228094581404*^-23},{34,32,-41,-0.93537087292458*^-25}};

R2=0.461526;
πp[p_]:=(p/16.53);
τ[T_]:=(1386/T);

γπ[T_,p_]:=Sum[-IJ[[i]][[4]]*IJ[[i]][[2]]*(((7.1-πp[p]))^(IJ[[i]][[2]]-1))*(((τ[T]-1.222))^(IJ[[i]][[3]])),{i,1,34}];

νmassLiq[T_,p_]:=((R2*T)/p)*πp[p]*γπ[T,p]; (*specific volume of water, L/kg*)
νLiq[T_,p_]:=νmassLiq[T,p]*0.01801528(*specific volume of water, L/mol*)

];
Module[{Jn,IJ,γoπ,γrπ,πp,τ,νmassVap,R2},
R2=0.461526;

Jn={{0,-9.6927686500217},{1,10.086655968018},{-5,-0.56087911283020*^-2},{-4,0.71452738081455*^-1},{-3,-0.40710498223928},{-2,1.4240819171444},{-1,-4.3839511319450},{2,-0.28408632460772},{3,0.021268463753307}};

IJ={{1,1,0,-0.17731742473213*^-2},{2,1,1,-0.17834862292358*^-1},{3,1,2,-0.45996013696365*^-1},{4,1,3,-0.57581259083432*^-1},{5,1,6,-0.50325278727930*^-1},{6,2,1,-0.33032641670203*^-4},{7,2,2,-0.18948987516315*^-3},{8,2,4,-0.39392777243355*^-2},{9,2,7,-0.43797295650573*^-1},{10,2,36,-0.26674547914087*^-4},{11,3,0,0.20481737692309*^-7},{12,3,1,0.43870667284435*^-6},{13,3,3,-0.32277677238570*^-4},{14,3,6,-0.15033924542148*^-2},{15,3,35,-0.40668253562649*^-1},{16,4,1,-0.78847309559367*^-9},{17,4,2,0.12790717852285*^-7},{18,4,3,0.48225372718507*^-6},{19,5,7,0.22922076337661*^-5},{20,6,3,-0.16714766451061*^-10},{21,6,16,-0.21171472321355*^-2},{22,6,35,-0.23895741934104*^2},{23,7,0,-0.59059564324270*^-17},{24,7,11,-0.12621808899101*^-5},{25,7,25,-0.38946842435739*^-1},{26,8,8,0.11256211360459*^-10},{27,8,36,-0.82311340897998*^1},{28,9,13,0.19809712802088*^-7},{29,10,4,0.10406965210174*^-18},{30,10,10,-0.10234747095929*^-12},{31,10,14,-0.10018179379511*^-8},{32,16,29,-0.80882908646985*^-10},{33,16,50,0.10693031879409},{34,18,57,-0.33662250574171},{35,20,20,0.89185845355421*^-24},{36,20,35,0.30629316876232*^-12},{37,20,48,-0.42002467698208*^-5},{38,21,21,-0.59056029685639*^-25},{39,22,53,0.37826947613457*^-5},{40,23,39,-0.12768608934681*^-14},{41,24,26,0.73087610595061*^-28},{42,24,40,0.55414715350778*^-16},{43,24,58,-0.94369707241210*^-6}};

πp[p_]:=(p/1);
τ[T_]:=(540/T);
γoπ[p_]:=(1/(πp[p]));
γrπ[T_,p_]:=Sum[IJ[[i]][[4]]*IJ[[i]][[2]]*(πp[p]^(IJ[[i]][[2]]-1))*(τ[T]-0.5)^IJ[[i]][[3]],{i,1,43}];

νmassVap[T_,p_]:=((R2*T)/p)*πp[p]*(γoπ[p]+γrπ[T,p]); (*specific volume of vapor, L/kg*)
νVap[T_,p_]:=νmassVap[T,p]*0.01801528(*specific volume of water vapor, L/mol*)];

Ti=Which[thermoBaric==1,
Which[location==1,250,
location==2,273.07,
location==3,325,
location==4,273.16],
thermoBaric==2,
Which[
location==1,200,
location==2,200,
location==3,303,
location==4,200]];

cpWater=75.3;(*liquid water, J/molK*)
cpIce=31.01;(*used approximation cp=-35.3622+7.68712(T^0.4) J/molK and averaged between 175K-265K.*)
cpVapor=33.75;(*J/molK, approximately constant between 273-373K at 0.1 MPa*)

Hm=6020;(*water standard heat of fusion, J/mol*)
Hv=40650;(*water standard heat of vaporization, J/mol*)
Hv0=44480;(*approximate heat of vaporization at 0 C (273K)*)
Hs=51000;(*approximate heat of sublimation @250K, J/mol*)

heatSublime=(250-Ti)*cpIce;(*amount of heat to get from initial temperature to sublimation point @250K*)
heatMelt=(273-Ti)*cpIce;(*amount of heat to get from initial temperature to melting point @273K*)
heatVapor=(373-Ti)*cpWater;
heatTriple=(273.16-Ti)*cpIce;

(* data = {temperature,%solid,%liquid,%vapor} *)
data=Which[thermoBaric==1,
Which[location==1,
Which[pr1>-9.4,{250,1,0,0,1},-9.4≥pr1≥-21.4,{250,((-21.4-pr1)/(-12)),0,((pr1+9.4)/(-12)),2},pr1<-21.4,{250,0,0,1,1}],
location==2,
Which[pr2>0.25,{273.07,0,1,0,1},0.25≥pr2≥-1,{273.07,((pr2-.25)/(-1.25)),((-1-pr2)/(-1.25)),0,2},pr2<-1,{273.07,1,0,0,1}],
location==3,
Which[pr3>-4.4,{325,0,1,0,1},-4.4≥pr3≥-7.8,{325,0,((-7.8-pr3)/(-3.4)),((pr3+4.4)/(-3.4)),2},pr3<-7.8,{325,0,0,1,1,1}],
location==4,
Which[pr4>-7.4,{273.16,0,1,0,1},-7.4≥pr4≥-10.8,{273.16,2*(((pr4+7.4)/(-3.4)))*(1-((pr4+7.4)/(-3.4))),1-2*(((pr4+7.4)/(-3.4)))+(((((pr4+7.4)/(-3.4))))^2),(((((pr4+7.4)/(-3.4))))^2),3},pr4<-10.8,{273.16,0,0,1,1}]],
thermoBaric==2,
Which[location==1,
Which[
1000*heat1≤heatSublime,{Ti+((1000*heat1)/cpIce),1,0,0,1},
heatSublime<1000*heat1<heatSublime+Hs,{250,((heatSublime+Hs-1000*heat1)/Hs),0,1-((heatSublime+Hs-1000*heat1)/Hs),2},
heatSublime+Hs≤1000*heat1,{250+(((1000*heat1-heatSublime-Hs))/cpVapor),0,0,1,1}],
location==2,Which[
1000*heat2≤heatMelt,{Ti+((1000*heat2)/cpIce),1,0,0,1},
heatMelt<1000*heat2<heatMelt+Hm,{273,((heatMelt+Hm-1000*heat2)/Hm),1-((heatMelt+Hm-1000*heat2)/Hm),0,2},
heatMelt+Hm≤1000*heat2,{273+(((1000*heat2-heatMelt-Hm))/cpWater),0,1,0,1}],
location==3,Which[
1000*heat3≤heatVapor,{Ti+((1000*heat3)/cpWater),0,1,0,1},
heatVapor<1000*heat3<heatVapor+Hv,{373,0,((heatVapor+Hv-1000*heat3)/Hv),1-((heatVapor+Hv-1000*heat3)/Hv),2},
heatVapor+Hv≤1000*heat3,{373+(((1000*heat3-heatVapor-Hv))/cpVapor),0,0,1,1}],
location==4,Which[
1000*heat4≤heatTriple,{Ti+((1000*heat4)/cpIce),1,0,0,1},
heatTriple<1000*heat4<heatTriple+Hs,{273.16,1-2*(((1000*heat4-heatTriple)/Hs))+(((((1000*heat4-heatTriple)/Hs)))^2),2*(((1000*heat4-heatTriple)/Hs))*(1-((1000*heat4-heatTriple)/Hs)),(((((1000*heat4-heatTriple)/Hs)))^2),3},
heatTriple+Hs≤1000*heat4,{273.16+(((1000*heat4-heatTriple-Hs))/cpVapor),0,0,1,1}]]];

pressure=Which[thermoBaric==1,
Which[location==1,Which[pr1>-9.4,pr1,-9.4≥pr1≥-21.4,-9.4,pr1<-21.4,-9.4+(((pr1+21.4))/5)],
location==2,Which[pr2>0.25,pr2,0.25≥pr2≥-1,0.25,pr2<-1,-1+pr2],
location==3,Which[pr3>-4.4,pr3,-4.4≥pr3≥-7.8,-4.4,pr3<-7.8,-4.4-((-7.8-pr3)/1.5)],
location==4,Which[pr4>-7.4,pr4,-7.4≥pr4≥-10.8,-7.4,pr4<-10.8,-7.4-((-10.8-pr4)/1.5)]],
thermoBaric==2,
Which[location==1,Log[1*^-4],
location==2,Log[0.1],
location==3,Log[0.1],
location==4,Log[0.00061]]];

Pc=22.12(*critical pressure, MPa*);
Tc=647.096(*critical temperature, K*);
R=8.314*^-6(*gas constant, m^3 MPa / (K * mol)*);

title=If[gibbsrule==True,Style[Text[Row@{Style["F",Italic]," = ",Style["C",Italic]," - ",Style["P",Italic]," + 2 = 1 - ",data[[5]]," + 2 = ",Plus[1,-data[[5]],2]}],17],Style[Text[" "],17]];

diagram=Show[
ListLogPlot[ice,Joined->True,PlotStyle->{Thick,Black},PlotRange->{{180,670},{1*^-7,100}},Frame->True,FrameLabel->{"temperature (K)","pressure (MPa)"},LabelStyle->{17,Black},ImageSize->1.05{550,395},AspectRatio->Full],
LogPlot[clausius[T],{T,200,273.15},PlotStyle->{Thick,Black}],
LogPlot[clausius[T],{T,273.15,647.096},PlotStyle->{Thick,Black}],
Graphics@{Text[Style[#1,18,#3],{#2[[1]],Log@#2[[2]]}]&@@@{{"solid",{220,0.02},{Opacity[0.9],Darker[Green]}},{"liquid",{350,10},{Opacity[0.7],Darker[Red]}},{"vapor",{550,0.1},Lighter[Blue,0.25]}},
PointSize@0.02,Point@{data[[1]],pressure},GrayLevel@0.55,PointSize@0.015,Point[{273.16,Log[0.0006]}],Point[{647.1,Log@40}]}];

Text1=Which[
pr1>-9.4,Text[ScientificForm[νIce[250,Exp[pr1]*(10^6)],4]],
-9.4≥pr1≥-21.4,Text[ScientificForm[(((pr1+9.4))/(-12))*νVap[250,Exp[-9.4]]+(((pr1+21.4))/12)*νIce[250,Exp[-9.4]*(10^6)],4]],
pr1<-21.4,Text[ScientificForm[νVap[250,Exp[-9.4+(((pr1+21.4))/5)]],4]]];
Text2=Which[
pr2>0.25,Text[ScientificForm[νLiq[273.07,Exp[pr2]],4]],
0.25≥pr2≥-1,Text[ScientificForm[(((pr2-0.25))/(-1.25))*νIce[273.07,Exp[0.25]*(10^6)]+(((pr2+1))/1.25)*νLiq[273.07,Exp[0.25]],4]],
pr2<-1,Text[ScientificForm[νIce[273.07,Exp[1.25+pr2]*(10^6)],4]]];
Text3=Which[
pr3>-4.4,Text[ScientificForm[νLiq[325,Exp[pr3]],4]],
-4.4≥pr3≥-7.8,Text[ScientificForm[(((pr3+4.4))/(-3.4))*νVap[325,Exp[-4.4]]+(((pr3+7.8))/3.4)*νLiq[325,Exp[-4.4]],4]],
pr3<-7.8,Text[ScientificForm[νVap[325,Exp[-4.4-((-7.8-pr3)/1.5)]],4]]];
Text4=Which[
pr4>-7.4,Text[ScientificForm[νLiq[273.16,Exp[pr4]],4]],
-7.4≥pr4≥-10.8,Text[ScientificForm[(((((pr4+7.4)/(-3.4))))^2)*νVap[273.16,Exp[-7.4]]+(1-(((((pr4+7.4)/(-3.4))))^2))*νLiq[273.16,Exp[-7.4]],4]],
pr4<-10.8,Text[ScientificForm[νVap[273.16,Exp[-7.4-((-10.8-pr4)/1.5)]],4]]];

Show[diagram,Epilog->{
Inset[
Labeled[
RectangleChart[{{0.2,data[[2]]},{0.2,data[[3]]},{0.2,data[[4]]}},
ChartLabels->Placed[{"solid","liquid","vapor"},Above],
LabelStyle->{17,Black},
If[thermoBaric==1&&location==2,ImageSize->{280,205},ImageSize->{280,220}],
PlotRange->{0,1.2},
Frame->True,
ChartStyle->{Opacity[0.4,Darker[Green]],Opacity[0.4,Darker[Red]],Lighter[Blue,0.65]},
BarSpacing->0.02,
AspectRatio->Full,
ImagePadding->30,
PlotRangeClipping->False,
ImageMargins->{{0,0},{0,5}},
Background->Directive[Opacity[0]]],

{"relative amount",Spacer[{0,0}]},{Left,Bottom},RotateLabel->True,LabelStyle->{Black,17},ImageMargins->{{0,0},{-7,-10}}],
If[thermoBaric==1&&location==2,{515,Log@0.000015},{510,Log@0.00004}],Background->Directive[Opacity[0]]],

If[thermoBaric==1&&location==2,
{
Inset[Framed[LogPlot[fice[x],{x,273,274},PlotRange->{{273,273.18},{0.01,6}},ImageSize->{220,140},TicksStyle->{14,Black},Ticks->{{{273.00,""}},{{0.1,""},{1,""},{5,""}}},Epilog->{PointSize[0.05],Point@{273.07,pressure}}]],
{440,-0.1},Background->White],
Line[{{277,0.25},{322,4.7}}],
Line[{{277,0.25},{322,-4.8}}],
Style[Text["5",{354,3.7}],14],
Style[Text["1",{354,2}],14],,
Style[Text["0.1",{359,-0.6}],14],
Style[Text["273.00",{347,-4}],14],
Style[Text["273.16",{515,-4}],14],
Rotate[Style[Text["pressure (MPa)",{332,0.5}],15],π/2],
Style[Text["temperature (K)",{430,-4}],15]},
Null
]
},PlotLabel->title,PlotRangeClipping->False]
],
Column[{Row@{Spacer[350],Control[{{gibbsrule,False,"display Gibbs phase rule"},Checkbox}]},
Row@{Control[{Text1,None}],Control[{Text2,None}],Control[{Text3,None}],Control[{Text4,None}],
Control[{{thermoBaric,2,""},{1:>"isothermal",2:>"isobaric"},ControlType->PopupMenu}],Control[{{location,2,""},{1:>"sublimation",2:>"melting",3:>"vaporization",4:>"triple point"},ControlType->PopupMenu}],
Spacer[10],
PaneSelector[{
1->PaneSelector[{
1->Row@{Control[{{pr1,2,"specific volume (L/mol)"},2,-34,-0.1,ImageSize->Small}],Spacer[10],Dynamic@Text1},
2->Row@{Control[{{pr2,1.5,"specific volume (L/mol)"},1.5,-2.25,-0.1,ImageSize->Small}],Spacer[10],Dynamic@Text2},
3->Row@{Control[{{pr3,-1,"specific volume (L/mol)"},-1,-11.2,-0.1,ImageSize->Small}],Spacer[10],Dynamic@Text3},4->Row@{Control[{{pr4,-2,"specific volume (L/mol)"},-2,-14.2,-0.1,ImageSize->Small}],Spacer[10],Dynamic@Text4}},Dynamic@location],
2->PaneSelector[{
1->Control[{{heat1,0,"heat (kJ)"},0,56,Slider,Appearance->"Labeled"}],
2->Control[{{heat2,0,"heat (kJ)"},0,14,Slider,Appearance->"Labeled"}],
3->Control[{{heat3,0,"heat (kJ)"},0,50,Slider,Appearance->"Labeled"}],
4->Control[{{heat4,0,"heat (kJ)"},0,56,Slider,Appearance->"Labeled"}]},Dynamic@location]
},Dynamic@thermoBaric]
}}],
SaveDefinitions->True]

Phase behavior for water is represented on a log pressure versus temperature diagram. Use the first drop-down menu to select for either an isothermal or isobaric process. Select the particular phase transition (sublimation, melting, vaporization, triple point) from the second drop-down menu. For isobaric transitions, you can add heat with the slider. For isothermal transitions, the specific volume is changed. The relative amounts of each phase present are shown in the bar graph. A log pressure-temperature graph with an expanded temperature scale is inserted for isothermal melting because of the steepness of the phase transition line. At the triple point, the three phases are present simultaneously. Select "display Gibbs phase rule" to show the phase rule applied at the location of the black dot.

The Gibbs phase rule, which determines how many intensive variables (e.g., temperature, pressure) can be varied independently, is:

F=C-P+2,

where F is the number of degrees of freedom, C is the number of components and P is the number of phases.

In a single-phase region, heat may be added isobarically to increase the temperature:

ΔH=n (Subscript[c, p]) ((Subscript[T, 2])-(Subscript[T, 1]))=Q,

where Q is the amount of heat added, (Subscript[c, p]) (kJ/kmol K) is the constant-pressure heat capacity, n is the number of moles, and (Subscript[T, 1]) and (Subscript[T, 2]) (K) are the initial and final temperatures. For isothermal processes, changes in the specific volume (L/mol) accompany phase changes; heat must be added or removed for these phase changes to take place.

In a two-phase region, the amount of each phase can be determined from:

((Subscript[P, H])/((Subscript[P, H]) (Subscript[P, L])))=((H-(Subscript[H, L]))/((Subscript[H, H])-(Subscript[H, L]))),

where (Subscript[P, H]) (mol) is the amount of the higher-energy phase, (Subscript[P, L]) (mol) is the amount of the lower-energy phase, H (kJ/mol) is the specific enthalpy of the system, (Subscript[H, H]) (kJ/mol) is the specific enthalpy of the saturated high-energy phase and (Subscript[H, L]) (kJ/mol) is the specific enthalpy of the saturated low-energy phase. Note that H, (Subscript[H, H]) and (Subscript[H, L]) must be evaluated at the same temperature and pressure.

Within the triple point, the amounts of each phase composition is complex [1, 2].

References

[1] Revised Release on the Equation of State 2006 for (Subscript[H, 2])O Ice Ih, Report Number IAPWS R10-06(2009), Doorwerth, Netherlands: The International Association for the Properties of Water and Steam. (Aug 3, 2018) www.iapws.org/relguide/Ice-2009.html.

[2] Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, Report Number IAPWS R7-97(2012), Lucerne, Switzerland: The International Association for the Properties of Water and Steam. (Aug 3, 2018) www.iapws.org/relguide/IF97-Rev.html.

[3] J. R. Elliott and C. T. Lira, Introductory Chemical Engineering Thermodynamics, 2nd ed., Upper Saddle River, NJ: Prentice Hall, 2012.

Contributed by: Neil C. Hendren

Additional contributions by: John L. Falconer and Rachael L. Baumann

(University of Colorado Boulder, Department of Chemical and Biological Engineering)