3.029 Spring 2022
Lecture 15 - 03/28/2022

Chemical Potential Common Tangent

Last lecture, we introduced partial molar quantities, and the tangent-intercept rule

In particular we saw how the chemical potential can be obtained graphically from molar free energy plots

In this lecture, we will investigate miscibility gaps using the common-tangent construction

Solution Models Overview

Recall the molar GIbbs free energy of mixing for an ideal solution, is given simply by the configurational entropy

molarGibbsEnergy["Ideal"][T_,R_:8.314][xB_]=R T ((1-xB) Log[1-xB]+xB Log[xB]);

Similarly for the regular solution, we introduce an enthalpic contribution due to inter-particle interactions

Where we’ve introduced

molarGibbsEnergy["Regular"][ω_,T_,R_:8.314][xB_]=(1-xB) xB ω+R T ((1-xB) Log[1-xB]+xB Log[xB]);

Finally, for a more general solution model - we additionally introduce activity coefficients in the configurational entropy

This allows us to introduce asymmetric molar free energies w/ different chemical potentials of the pure compounds

molarGibbsEnergy["Regular-Activity"][γA_,γB_,ω_,T_,R_:8.314][xB_]=(1-xB) xB ω+R T ((1-xB) Log[γA(1-xB)]+xB Log[γB xB]);

Manipulate[
Plot[{
molarGibbsEnergy["Ideal"][temperature][xB],
molarGibbsEnergy["Regular"][ω,temperature][xB],
molarGibbsEnergy["Regular-Activity"][γA,γB,ω,temperature][xB]
},{xB,0,1},PlotStyle->Thread[{{Red,Blue,Black},Thick}],Frame->True,LabelStyle->Directive[Black,16],ImageSize->550,PlotLegends->LineLegend[Automatic,{"Ideal Solution Model","Regular Solution Model","Regular-Activity Solution Model"}],FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}],
{{temperature,1},.001,1},
Delimiter,
{{ω,1},-5,5},
Delimiter,
{{γA,0.8},0.5,1.5},{{γB,1.2},0.5,1.5},
Paneled->False,SaveDefinitions->True
]

Common Tangent Construction

When the Enthalpy of mixing is positive, in  the regular and general solution models above, there’s a competition between molecular interactions and the configurational entropy at low temperatures, manifesting in the loss of convexity

As we’ve seen before, a loss of convexity in thermodynamic potentials suggests a phase transition

In particular, this introduces a range of compositions where the free energy of a single homogeneous solution is higher than the weighted-average of the mechanical mixture of two solutions with the same average composition

As such, the homogenous solution separates into two homogenous solutions with different compositions each to reach equilibrium

At a single temperature, these limiting compositions can be computed with the common tangent construction

Similar to L05 and P02, we will compute the common tangent construction numerically using FindRoot and the appropriate equations given below

where  are the two limiting compositions of the two-phase region on the A-rich and B-rich sides respectively

chemicalPotentialA["Regular-Activity"][γA_,γB_,ω_,T_,R_:8.314][xB_]=molarGibbsEnergy["Regular-Activity"][γA,γB,ω,T,R][xB]-xB molarGibbsEnergy["Regular-Activity"][γA,γB,ω,T,R]'[xB];

First, we want to check if the function indeed loses convexity

We’ll use a built-in function we introduced in L06

convexQ[modelName_][params__]:=FunctionConvexity[{molarGibbsEnergy[modelName][params][xB],1>xB>0},xB]===1

E.g. for negative and positive interaction parameters at a fixed temperature

convexQ["Regular-Activity"][0.8,1.1,-1,0.25]
convexQ["Regular-Activity"][0.8,1.1,5,0.25]

And for various temperatures at fixed positive interaction parameters

convexQ["Regular-Activity"][0.8,1.1,5,0.25]
convexQ["Regular-Activity"][0.8,1.1,5,0.35]

This suggest that, for positive interaction parameters, there exists a critical temperature at which the Gibbs free energy loses convexity

Minimize[{molarGibbsEnergy["Regular-Activity"][0.8,1.1,5,T]''[xB],1>=xB>=0,T>0},xB,Reals]
Solve[%[[1,1,1,1]] ==0,T]

We’ll use a similar lookup-table approach for initial guesses

Notice we now have a multidimensional coordinate space (γA,γB,ω,T)

In-fact, let’s also use the same mechanism to store our results

Clear[commonTangentConstruction,commonTangentResults]

commonTangentResults["Regular-Activity"]=<|
{0.8,1.1,5,0.25}->
<|
"compositions"->{0.15,0.85},
"function"->molarGibbsEnergy["Regular-Activity"][0.8,1.1,5,0.25]|>
|>;

commonTangentConstruction[modelName_][params__]:=Module[{
reversedGuesses,guessesX,sol,pw,initialGuess={0.8,1.1,5,0.25}},

If[{params}!=initialGuess&&KeyExistsQ[commonTangentResults[modelName],{params}],
Return[commonTangentResults[modelName][{params}]]];

If[convexQ[modelName][params],
AssociateTo[
commonTangentResults[modelName],{params}-><|
"compositions"->{-1,2},
"function"->molarGibbsEnergy[modelName][params]|>];
Return[commonTangentResults[modelName][{params}]]];

reversedGuesses=AssociationMap[Reverse,commonTangentResults[modelName]];
guessesX=First[Nearest[reversedGuesses,{params}]]["compositions"];

sol=Quiet@FindRoot[{

chemicalPotentialA[modelName][params][x1]==chemicalPotentialA[modelName][params][x2],

molarGibbsEnergy[modelName][params]'[x1]==molarGibbsEnergy[modelName][params]'[x2]

},
{x1,guessesX[[1]]},{x2,guessesX[[2]]}];

pw=Piecewise[{{
molarGibbsEnergy[modelName][params][x1]+molarGibbsEnergy[modelName][params]'[x1] (#-x1),x2>#>x1}
},
molarGibbsEnergy[modelName][params][#]]&/.sol;

AssociateTo[
commonTangentResults[modelName],
{params}-><|
"compositions"->{x1,x2}/.sol,
"function"->pw|>];

commonTangentResults[modelName][{params}]
]

Let’s see how this looks interactively!

Manipulate[
Block[{asc,pts,func},
asc=commonTangentConstruction["Regular-Activity"][γA,γB,ω,temperature];
{pts,func}=Values[asc];
Plot[{func[xB],molarGibbsEnergy["Regular-Activity"][γA,γB,ω,temperature][xB]},{xB,0,1},Filling->{1->{2}},GridLines->{pts,None},GridLinesStyle->Directive[Gray, Dashed],PlotStyle->{{Red,Thick},{Black,Thick}},Frame->True,LabelStyle->Directive[Black,16],ImageSize->550,PlotLegends->LineLegend[Automatic,{"Common Tangent Construction","Regular-Activity Solution Model"}],FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}]],
{{temperature,0.25},.001,1},
Delimiter,
{{ω,5},-5,5},
Delimiter,
{{γA,0.8},0.5,1.5},{{γB,1.2},0.5,1.5},
Paneled->False
]

Miscibility Gap

Note we got different limiting concentrations for each temperature

Repeating the construction over a range of temperatures, we can construct the miscibility gap of the two-phase region

i.e. concentrations for which the homogeneous solution phase-separates

we will connect these using tie-lines to illustrate the coexistence between the two phases

Monitor[
Do[commonTangentConstruction["Regular-Activity"][0.8,1.1,5,T],{T,Join[Range[0.25,0.3,0.005],Range[0.25- 0.005,0.15,-0.005]]}],
T]

interpolateTielines[selectedCompositions_][minT_,maxT_,ΔT_]:=Block[{Ts,newTs,minComp,maxComp,minInt,maxInt,newComps},
Ts=Keys[selectedCompositions][[All,-1]];
newTs=Range[minT,maxT,ΔT];

{minComp,maxComp}=((Values[selectedCompositions[[All,"compositions"]]])^\[Transpose]);
minInt=Interpolation[Thread[{Ts,minComp}]];
maxInt=Interpolation[Thread[{Ts,maxComp}]];

newComps=((Through[{minInt,maxInt}[newTs]])^\[Transpose]);
MapThread[Line[Thread[{#1,#2}]]&,{newComps,newTs}]
]

miscibilityGap=With[{
tieLines=interpolateTielines[KeySelect[commonTangentResults["Regular-Activity"],Take[#,3]=={0.8,1.1,5}&]][0.15,0.3,0.0025]},
Graphics[tieLines,AspectRatio->1/2,FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Temperature, T"},Frame->True,PlotRange->{{0,1},{0.15,0.35}},PlotRangeClipping->True,ImageSize->500,BaseStyle->18]]

Lever Rule

The miscibility gap indicates the unstable range of compositions the two phases separate into

The composition of each phase is determined by the limits of the miscibility gap

In order to know the fraction of each phase in the two phase region, the lever rule is used

where the lever lengths are determined by the compositions as

For example, for T=0.25, (Subscript[x, B])=0.75

ClearAll[MPLColorMap]
<<"http://pastebin.com/raw/pFsb4ZBS";

phaseFractions[modelName_][params__][xB_]:=Block[{commonTangent,comps},
commonTangent=commonTangentConstruction[modelName][params];
comps=commonTangent["compositions"];

If[Not[#1<xB<#2&@@comps],Return[{{xB,1}}]];
If[comps=={-1,2},Return[{{xB,1}}]];

{{#1,((#2-xB)/(#2-#1))},{#2,((xB-#1)/(#2-#1))}}&@@comps
]

phaseFractionGraphics[modelName_][params__][T_,xB_]:=Block[{ϕ=phaseFractions[modelName][params][xB],
cf=MPLColorMap["Magma"]},
If[Length[ϕ]==1,
Graphics[{cf[xB],Disk[{xB,T+0.075},Offset[20]],Text["100%",{xB,T+0.11}]}],
Graphics[{
{Arrowheads[{0.025,0.025}],Arrow[{{xB,T+0.095},{ϕ[[1,1]],T+0.095}}],Arrow[{{xB,T+0.095},{ϕ[[2,1]],T+0.095}}],
Thick,Line[Thread[{ϕ[[All,1]],T+0.075}],VertexColors->cf/@ϕ[[All,1]]]
},
{cf[#1],Disk[{#1,T+0.075},Offset[5 (Sqrt[1+15 #2])]],Text[StringTemplate["`1`%"][Round[#2 100]],{#1,T+0.11}]}&@@@ϕ
}]
]
]

With[{T=0.25,xB=0.75},
Show[miscibilityGap,phaseFractionGraphics["Regular-Activity"][0.8,1.1,5,T][T,xB],GridLines->{{xB},{T}},PlotRange->{{-0.05,1.05},{0.15,0.425}},GridLinesStyle->Directive[Red,Thick,Dashed]]]

Let’s visualize this interactively

Manipulate[
Show[miscibilityGap,phaseFractionGraphics["Regular-Activity"][0.8,1.1,5,T][T,xB],GridLines->{{xB},{T}},PlotRange->{{-0.05,1.05},{0.15,0.425}},GridLinesStyle->Directive[Red,Thick,Dashed]],{{T,0.275,"Temperature"},0.15,0.305},{{xB,0.5,"Composition"},0,1},Paneled->False]

Convex Hull Construction

In a binary compound, there are more than two competing phases and finding all common tangent constructions between these phases using the previous technique can become cumbersome

Instead, we’ll demonstrate an equivalent technique using a computational geometry technique called the convex hull

In-fact, we’ve already seen this when we were exploring the Legendre transform using supporting hyperplanes

Let’s start by defining a binary compound with three phases given by

Clear[phase]
phase[name_][T_][xs_List]:=phase[name][T]/@xs

phase["A"][T_][xB_]=molarGibbsEnergy["Regular-Activity"][0.9,1.1,3,T][xB];
phase["B"][T_][xB_]=molarGibbsEnergy["Regular-Activity"][1.2,0.9,3,T][xB];
phase["C"][T_][xB_]=molarGibbsEnergy["Regular-Activity"][0.9,0.9,4,T][xB];

Manipulate[Plot[{
phase["A"][T][xB],
phase["B"][T][xB],
phase["C"][T][xB]
},{xB,0,1},PlotStyle->Thread[{{Red,Blue,Black},Thick}],Frame->True,LabelStyle->Directive[Black,16],ImageSize->550,PlotLegends->LineLegend[Automatic,{"Phase A","Phase B","Phase C"}],FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}],{{T,0.3,"Temperature"},0.1,1},Paneled->False,SaveDefinitions->True]

Note our x Log [x] function is Indeterminate at x=0 | x =1

phase["A"][T][0]

As such, we also define the following limits

phase["A"][T_][0]=phase["A"][T_][0.0]=Limit[phase["A"][T][x],x->0];
phase["A"][T_][1]=phase["A"][T_][1.0]=Limit[phase["A"][T][x],x->1];

phase["B"][T_][0]=phase["B"][T_][0.0]=Limit[phase["B"][T][x],x->0];
phase["B"][T_][1]=phase["B"][T_][1.0]=Limit[phase["B"][T][x],x->1];

phase["C"][T_][0]=phase["C"][T_][0.0]=Limit[phase["C"][T][x],x->0];
phase["C"][T_][1]=phase["C"][T_][1.0]=Limit[phase["C"][T][x],x->1];

We are interested in computing the convex-hull of the epigraph of the functions

and more specifically - the convex hull of the epigraph on the minimum of the functions

We’ll start by sampling the functions at evenly sampled concentration increments and returning the minimum values

minimumPoints[functionList_, numPts_] :=Thread[{Subdivide[numPts],Min/@Thread[Through[functionList[Subdivide[numPts]]]]}]

Plot[{
phase["A"][0.3][xB],
phase["B"][0.3][xB],
phase["C"][0.3][xB]
},{xB,0,1},PlotStyle->LightGray,Frame->True,LabelStyle->Directive[Black,16],ImageSize->550,Epilog->{PointSize[Large],Opacity[0.15],RGBColor[0,1/3,0],Point[minimumPoints[{phase["A"][0.3],phase["B"][0.3],phase["C"][0.3]},1000]]},FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}]

Next, we “convexify” the function by only returning the points on the convex envelope or ‘hull’

getConvexHullVertices[candidatePoints_] := Module[{
maxY = Max[Last /@ candidatePoints],
amendedPoints
}, 
amendedPoints = Join[{{0., 1 + maxY}}, candidatePoints, {{1., 1 + maxY}}]; 
MeshCoordinates[ConvexHullMesh[amendedPoints]][[2;;-2]]
]

With[{ch=getConvexHullVertices[minimumPoints[{phase["A"][0.3],phase["B"][0.3],phase["C"][0.3]},1000]]},
Plot[{
phase["A"][0.3][xB],
phase["B"][0.3][xB],
phase["C"][0.3][xB]
},{xB,0,1},PlotStyle->LightGray,Frame->True,LabelStyle->Directive[Black,16],ImageSize->550,Epilog->{PointSize[Large],Opacity[0.15],RGBColor[0,1/3,0],Point[ch]},FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}]]

As we can see, the function only returns the points on the epigraph - without connecting the common tangent points

Each “set” of points on the epigraph corresponds to a homogeneous single-phase region b/w two phases
E.g. the left-most points correspond to a single-phase region b/w Phase-A and Phase-B

We can use this to split the epigraph according to these single phase-regions

With[{ch=Split[getConvexHullVertices[minimumPoints[{phase["A"][0.3],phase["B"][0.3],phase["C"][0.3]},1000]],First[#2]-First[#1]<1.001/1000&]},
Plot[{
phase["A"][0.3][xB],
phase["B"][0.3][xB],
phase["C"][0.3][xB]
},{xB,0,1},PlotStyle->LightGray,Frame->True,LabelStyle->Directive[Black,16],ImageSize->550,Epilog->{PointSize[Large],Opacity[0.15],RGBColor[0,1/3,0],Point[ch[[1]]],Purple,Point[ch[[2]]],Cyan,Point[ch[[3]]]},FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}]]

In-fact, we can color the single-phase region points according to their composition

and connect the various miscibility gaps with dotted lines

ClearAll[MPLColorMap]
<<"http://pastebin.com/raw/pFsb4ZBS";

freeEnergyLineGraphicAndMiscibilityGaps[functionList_,numPts_] := 
Module[{convexHullPoints = getConvexHullVertices[minimumPoints[functionList, numPts]],
singlePhasePoints,miscibilityGaps,cf=MPLColorMap["Magma"],cols},
singlePhasePoints=Split[convexHullPoints,First[#2]-First[#1]<1.001/numPts&];
miscibilityGaps={Last[#1],First[#2]}&@@@Partition[singlePhasePoints,2,1];
cols=Map[cf[First[#]]&,singlePhasePoints,{2}];

{Graphics[{Thick,MapThread[Line[#1,VertexColors->#2]&,{singlePhasePoints,cols}],Black,Dotted,Line/@miscibilityGaps}],miscibilityGaps}
]

Show[
Plot[{
phase["A"][0.3][xB],
phase["B"][0.3][xB],
phase["C"][0.3][xB]
},{xB,0,1},PlotStyle->LightGray,Frame->True,LabelStyle->Directive[Black,16],ImageSize->550,FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}],
First@freeEnergyLineGraphicAndMiscibilityGaps[{phase["A"][0.3],phase["B"][0.3],phase["C"][0.3]},1000]]

Let’s see how this looks as a function of temperature interactively

Manipulate[
Show[
Plot[{
phase["A"][T][xB],
phase["B"][T][xB],
phase["C"][T][xB]
},{xB,0,1},PlotStyle->LightGray,Frame->True,LabelStyle->Directive[Black,16],ImageSize->550,FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}],
First@freeEnergyLineGraphicAndMiscibilityGaps[{phase["A"][T],phase["B"][T],phase["C"][T]},1000]],
{{T,0.3,"Temperature"},0.1,0.75},Paneled->False]

Binary Phase Diagram

Similar to the miscibility gap above, we can put all this together to plot a binary phase diagram as a function of temperature

backgroundGraphic[Tmin_,Tmax_]:=backgroundGraphic[Tmin,Tmax]=DensityPlot[x,{x,0,1},{y,Tmin,Tmax},ColorFunction->MPLColorMap["Magma"]][[1]]

miscibilityGaps[functionList_,temperature_,numPts_] := 
Module[{convexHullPoints = getConvexHullVertices[minimumPoints[Through[functionList[temperature]], numPts]],singlePhasePoints,miscibilityGaps,cf=MPLColorMap["Magma"],cols},
singlePhasePoints=Split[convexHullPoints,First[#2]-First[#1]<1.001/numPts&];
miscibilityGaps={Last[#1],First[#2]}&@@@Partition[singlePhasePoints,2,1];
Line/@Map[Thread[{#,temperature}]&,miscibilityGaps[[All,All,1]],{2}]
]

Graphics[{backgroundGraphic[0.1,0.6],White,miscibilityGaps[{phase["A"],phase["B"],phase["C"]},#,1000]&/@Range[0.1,0.6,0.005]},PlotRange->{{0,1},{0.1,0.6}},FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Temperature, T"},Frame->True,PlotRangeClipping->True,AspectRatio->1/GoldenRatio,ImageSize->500,LabelStyle->Directive[Black,16]]

Let’s put all of this together in an interactive phase diagram explorer!

showTieLines[{Tmin_,Tmax_,ΔT_}][tieLines_,temperature_]:=Graphics[{
backgroundGraphic[Tmin,Tmax],
{White,Values[tieLines]},
{LightGray,Line[{{0,temperature},{1,temperature}}]}
},PlotRange->{{0,1},{Tmin,Tmax}},FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Temperature, T"},Frame->True,PlotRangeClipping->True,ImageSize->500,LabelStyle->Directive[Black,16],ImagePadding->{{75,15},{50,15}},AspectRatio->1/GoldenRatio]

showFreeEnergyGraphicAndReturnMiscibilityGaps[functionList_,temperature_,numPts_]:=Block[{fs=Through[functionList[temperature]],graphic,tieline},
{graphic,tieline}=freeEnergyLineGraphicAndMiscibilityGaps[fs,numPts];
graphic=Show[Plot[Through[fs[xB]],{xB,0,1},PlotStyle->LightGray,Frame->True,LabelStyle->Directive[Black,16],ImageSize->500,ImagePadding->{{75,15},{50,15}},FrameLabel->{"Mole Fraction B, (Subscript[x, B])","Molar Gibbs Energy, (UnderBar[G])"}],graphic];
{graphic,tieline}
]

phaseDiagramExplorer[functionList_,{Tmin_,Tmax_,ΔT_}, numPts_]:= 
DynamicModule[{
tieLines = <||>,allTemperatures,meanTemperature,showFreeEnergyAndPhaseDiagramGraphic
},

allTemperatures = Range[Tmin, Tmax, ΔT];
meanTemperature = allTemperatures[[Round[Length[allTemperatures]/2]]];

showFreeEnergyAndPhaseDiagramGraphic[temperature_]:=Block[{graphic,tieline},
If[KeyExistsQ[tieLines,temperature],

Column[{
showFreeEnergyGraphicAndReturnMiscibilityGaps[functionList,temperature,numPts]//First,
showTieLines[{Tmin,Tmax,ΔT}][tieLines,temperature]}
],

{graphic,tieline}=showFreeEnergyGraphicAndReturnMiscibilityGaps[functionList,temperature,numPts];
AssociateTo[tieLines,temperature->Line/@Map[Thread[{#,temperature}]&,tieline[[All,All,1]],{2}]];
Column[{graphic,showTieLines[{Tmin,Tmax,ΔT}][tieLines,temperature]}]
]
];

Manipulate[
If[computeAll && Length[tieLines] < Length[allTemperatures], 
 temperature= First[Sort[Complement[allTemperatures, Keys[tieLines]]]]];

If[Length[tieLines]≥ Length[allTemperatures], computeAll=False];

showFreeEnergyAndPhaseDiagramGraphic[temperature],

{{temperature,meanTemperature,Style["Temperature",FontSize->12]},allTemperatures,Slider},
Delimiter,
{{computeAll,False,Style["Compute all tie-lines",FontSize->12]},{False, True}},Paneled->False
]
]

phaseDiagramExplorer[{phase["A"],phase["B"],phase["C"]},{0.1,0.6,0.005},1000]