Phase Separations

Reference

Title

Phase Separations

Author

Thomas G. Anderson (tga@stanford.edu)

Summary

This package computes the compositions of coexisting phases in a binary mixture.

Copyright

© Copyright 2002 Thomas G. Anderson

Permission is granted to distribute this file for any purpose except for inclusion in commercial software or program collections. This copyright notice must remain intact.

Notebook Version

1.0

Mathematica Version

4.1.5

History

1.0	21 January 2002	First public version

Keywords

phase separation, phase diagram, binary mixture, miscibility

Sources

Package template adapted from Maeder, R. "Programming in Mathematica," 3d Ed.

Warnings

none

Limitations

BinaryPhaseSeparations uses a finite sample of points in the free energy function to determine the regions of phase separation, so it is possible that small two-phase regions will be missed. The points are chosen by the adaptive sampling algorithm used in the built-in Plot function, so only features resolvable by Plot will be detected by BinaryPhaseSeparations. In other words, if you can see a region of phase separation by plotting the free energy function, it can be found.

Discussion

In chemical thermodynamics, a mixture of two or more components is unstable with respect to phase separation if the second derivative of the free energy with respect to composition is negative. At equilibrium, the chemical potential of the components is the same between phases. Geometrically speaking, this means that the free energy curve has a common tangent at the compositions of the coexisting phases. Conversely, the composition of coexisting phases may be found by finding a line that is tangent to the free energy curve at two points; this is the method of double tangent construction.

This package implements the method of double tangent construction using an adaptive algorithm to sample the free energy function find the approximate composition of coexisting phases, followed by refinement to arbitrary precision. Multiple distinct regions of phase separation are found automatically by a simple recursive algorithm.

Briefly, the algorithm works as follows. First, the free energy function is sampled using the built-in Plot function (with the graphic output suppressed). These points are used to numerically compute the curvature of the function. If any points of negative curvature are found, the algorithm proceeds to the next step. The composition having the smallest (most negative) curvature is identified as (Subscript[x, 0]). The point of minimum free energy is taken as the first approximation of one phase, (Subscript[x, 1]). The first approximation of the second phase, (Subscript[x, 2]), is taken to be that point on the other side of (Subscript[x, 0]) for which the slope of the line between it and ((Subscript[x, 1]), (Subscript[y, 1])) is closest to zero. Then the set of points is skewed so that (Subscript[x, 1]) and (Subscript[x, 2]) have the same y coordinate. The process is repeated until the result no longer changes. To refine the values of (Subscript[x, 1]) and (Subscript[x, 2]), the free energy of the two phases present at (Subscript[x, 0]), (Subscript[y, sep]) = [((Subscript[x, 0]) - (Subscript[x, 1])) 
                (Subscript[y, 2]) + ((Subscript[x, 2]) - (Subscript[x, 0])) (Subscript[y, 1])]/((Subscript[x, 2]) - (Subscript[x, 1])), is minimized with respect to (Subscript[x, 1]) and (Subscript[x, 2]).

Requirements

Requires the Standard Packages NumericalMath`NLimit` and Utilities`FilterOptions`.

Interface

This part declares the publicly visible functions, options, and values.

Set up the package context, including public imports

\(\(BeginPackage["PhaseSeparations`", {"NumericalMath`NLimit`"}];\)\)

Usage messages for the exported functions and the context itself

PhaseSeparations::usage = "PhaseSeparations is a package that \
computes the compositions of coexisting phases in multicomponent \
mixtures.";

BinaryPhaseSeparations::usage = "BinaryPhaseSeparations[f, x] \
calculates the composition intervals for phase separation in a binary mixture \
with free energy function f[x] over the composition range x = 0 to 1. \
BinaryPhases[f, {x, xmin, xmax}] restricts the composition range to {xmin, \
xmax}.";

DrawPhaseDiagram::usage = "DrawPhaseDiagram[{{int1, p1}, {int2, \
p2}, ...}] draws a phase diagram by interpolating between the given intervals \
of phase coexistence inti at the values of the potential pi.";

BinaryPhaseDiagram::usage = "BinaryPhaseDiagram[f, x, {p, pmin, \
pmax}] calculates the phase diagram for a binary mixture with a free energy \
function f[x, p] over the composition range x = 0 to 1 and p = pmin to pmax, \
where p is a potential such as temperature. BinaryPhaseDiagram[f, {x, xmin, \
xmax}, {p, pmin, pmax}] restricts the composition range to {xmin, \
xmax}.";

SepPoints::usage = "SepPoints is an option for \
BinaryPhaseSeparations that specifies how many sample points to \
use.";

SepDivision::usage = "SepDivision is an option for \
BinaryPhaseSeparations which specifies the maximum amount of subdivision to \
be used in attempting to approximate the free energy curve.";

\(\(VerticalPoints::usage = "VerticalPoints is an option for BinaryPhaseDiagram that specifies how many values to use on the y axis.";\)\)

VerticalDivision::usage = "VerticalDivision is an option for \
BinaryPhaseDiagram which specifies the maximum amount of subdivision to be \
used on the y axis.";

Error messages for the exported objects

BinaryPhaseSeparations::noref = "The separation interval `1` could \
not be refined.";

Implementation

Begin the private context (implementation part)

Begin["PhaseSeparations`Private`"];

Read in any hidden imports

Needs["Utilities`FilterOptions`"]

Unprotect any system functions for which definitions will be made

n/a

Definition of auxiliary functions and local (static) variables

skew

Skews the curve given by the points vals about the point (Subscript[x, 0]) by addition of a linear function. The slope of the added function is determined by finding the minimum point on the curve, (Subscript[P, 1]), and then finding that point (Subscript[P, 2]) on the other side of (Subscript[x, 0]) such that the slope of the line (Subscript[P, 1]) (Subscript[P, 2]) is closest to zero.

To avoid unnecessary computation, points on the side of (Subscript[P, 1])away from (Subscript[x, 0]) are discarded, as is the point (Subscript[x, 0]).

skew[pts_, x0_] := 
      Module[{pmin, lpts, rpts, slope}, 
pmin = 
          Reverse[First[Sort[Reverse /@ pts]]]; 
If[
          pmin[[1]] < 
            x0, 
 (* 
            minimum is to left of x0 *) 
(
lpts = 
              Select[pts, (pmin[[1]] ≤ #[[1]] < 
                      x0) &]; 
rpts = 
              Select[pts, (x0 < #[[1]]) &]; 
slope = 
              Min[(Divide @@ Reverse[# - pmin]) & /@ 
                  rpts];
), 
 (* 
            minimum is to right of x0 *) 
(
lpts = 
              Select[pts, (#[[1]] < 
                      x0) &]; 
rpts = 
              Select[pts, (x0 < #[[1]] ≤ 
                      pmin[[1]]) &]; 
slope = 
              Max[(Divide @@ Reverse[# - pmin]) & /@ 
                  lpts];
)
]; 
{#[[1]], #[[2]] - 
                slope  #[[1]]} & /@ 
          Join[lpts, rpts]
]

concavity

Determines the concavity for a curve running through three points; this expression corresponds to the coefficient of the quadratic term of a parabola fit to the points.

concavity[{x1_, y1_}, {x2_, y2_}, {x3_, 
          y3_}] := (x1 (y3 - y2) + x2 (y1 - y3) + 
            x3 (y2 - y1))/((x1 - x2) (x1 - x3) (x2 - 
                x3))

findSeparations

This function recursively determines the intervals of phase separation over the given composition range. First, the point of greatest instability (negative concavity) is identified and its phases calculated. Then the function calls itself for the composition ranges on either side of this separation. If no unstable point is found, returns an empty list.

findSeparations[pts_ /; Length[pts] < 3] = {};

findSeparations[pts_] := 
      Module[{xc, xcmin, skewpts, x1, x2}, 
xc = 
          Map[{#[[2, 1]], 
                concavity @@ #} &, 
            Partition[pts, 3, 1]]; 
xcmin = 
          Reverse[First[Sort[Reverse /@ xc]]]; 
If[
          xcmin[[2]] ≥ 0, 
          Return[{}]]; 

skewpts = 
          FixedPoint[
            skew[#, xcmin[[1]]] &, pts, 
            10]; 
{x1, 
            x2} = Sort[Reverse /@ skewpts][[{1, 2}, 
            2]]; 

{findSeparations[
            Select[pts, (#[[1]] < 
                    x1) &]], 
Interval[{x1, 
              x2}], 
findSeparations[
            Select[pts, (#[[1]] > 
                    x2) &]]}
]

refineSeparation

Refines a calculated phase separation by minimizing the total free energy of the two phases. For a point (Subscript[x, 0]) inside the region of phase separation, {(Subscript[x, 1]), (Subscript[x, 2])}, the total free energy of the coexisting phases is, by the lever rule,

f ((Subscript[x, 0])) = ((((Subscript[x, 2]) - (Subscript[x, 0])))/(((Subscript[x, 2]) - (Subscript[x, 1])))) 
          f ((Subscript[x, 1])) + ((((Subscript[x, 0]) - (Subscript[x, 1])))/(((Subscript[x, 2]) - (Subscript[x, 1])))) 
          f ((Subscript[x, 2]))

Attributes[refineSeparation] = {HoldFirst};

refineSeparation[f_, x_Symbol, Interval[{x1o_, x2o_}], opts___] := 
      Module[{x0, 
          fm}, 
x0 = (x1o + x2o)/
            2; 
fm = 
          FindMinimum[
If[! (0 < x1 < x0 < x2 < 
                    1), $MaxMachineNumber, 
ReleaseHold[((x2 - x0) (Hold[f] /. 
                            x -> x1) + (x0 - x1) (Hold[f] /. 
                            x -> x2))/(x2 - 
                      x1)]], 
{x1, {x1o, 
                x1o + (x2o - x1o)/100}}, 
{x2, {x2o, 
                x2o - (x2o - x1o)/
                    100}}, 
opts
]; 
If[
          Head[fm] === 
            FindMinimum, (Message[BinaryPhaseSeparations::noref, {x1, x2}]; 
            Interval[{x1o, x2o}]), 
Interval[{x1, x2} /. 
              fm[[2]]]]
]

connectIntervals

For a pair of binary phase separation intervals at consecutive potential values, connects the endpoints of intervals that overlap in composition.

connectIntervals[{i1_Interval, p1_}, {i2_Interval, p2_}] := 
      Outer[
If[
            IntervalIntersection[#1, #2] === 
              Interval[], {}, {Line[{{#1[[1, 
                      1]], 
                    p1}, {#2[[1, 1]], 
                    p2}}], Line[{{#1[[1, 
                      2]], 
                    p1}, {#2[[1, 2]], 
                    p2}}]}] &, 
List @@ Split[i1], 
        List @@ Split[i2]
]

joinLineSegments

List of rules for joining together Lines that have common endpoints.

joinLineSegments = {{l1___, Line[{p1__, p2_}], l2___, 
              Line[{p2_, p3__}], l3___} :> {Line[{p1, p2, p3}], 
              l1, l2, l3}, {l1___, Line[{p2_, p3__}], l2___, 
              Line[{p1__, p2_}], l3___} :> {Line[{p1, p2, p3}], 
              l1, l2, l3}};

criticalPoints

For a pair of binary phase separation intervals at consecutive potential values, identifies separation regions that are present at one potential but not the other, and interpolates their coordinates to obtain miscibility critical points.

criticalPoints[{i1_Interval, p1_}, {i2_Interval, p2_}] := 
      Join[Select[
            List @@ Split[
                i1], (IntervalIntersection[#, i2] === Interval[]) &], 
          Select[List @@ 
              Split[i2], (IntervalIntersection[#, i1] === 
                  Interval[]) &]] /. 
        Interval[{x1_, x2_}] :> 
          Point[{(x1 + x2)/2, (p1 + p2)/2}]

Definition of the exported functions

BinaryPhaseSeparations

Calculates the intervals of two-phase coexistence for a given free energy function. In the event that a free energy function is not well-behaved near the edges of the composition range, we allow explicit limits on the composition to be provided as optional arguments.

Attributes[BinaryPhaseSeparations] = {HoldFirst, 
          ReadProtected};

Options[BinaryPhaseSeparations] = 
        Join[Options[FindMinimum], {SepPoints -> 50, 
            SepDivision -> 30. }];

BinaryPhaseSeparations[f_, x_Symbol, opts___?OptionQ] := 
      BinaryPhaseSeparations[f, {x, 0, 1}, opts]

BinaryPhaseSeparations[f_, {x_Symbol, xmin_?NumericQ, xmax_?NumericQ}, 
          opts___?OptionQ] /; (0 ≤ xmin < xmax ≤ 
            1) := Module[{optsfm, sp, sd, pts}, 
optsfm = 
          FilterOptions[FindMinimum, opts, 
            Sequence @@ 
              Options[BinaryPhaseSeparations]]; 
{sp, 
            sd} = {SepPoints, SepDivision} /. 
            Flatten[{opts, 
                Options[
                  BinaryPhaseSeparations]}]; 
pts = Plot[
              f, {x, xmin, xmax}, DisplayFunction -> Identity, 
              PlotPoints -> sp, 
              PlotDivision -> sd][[1, 1, 1, 
            1]]; 
IntervalUnion @@ (refineSeparation[f, x, #, Evaluate[optsfm]] & /@ 
              Flatten[findSeparations[pts]])
]

DrawPhaseDiagram

Given a list of phase coexistence regions {{(Subscript[int, 1]),(Subscript[p, 1])},{(Subscript[int, 2]),(Subscript[p, 2])},...} where the (Subscript[int, i]) are the separation Intervals at the potentials (Subscript[p, i]), "connects the dots" to draw the phase diagram over the range of potentials given.

Attributes[DrawPhaseDiagram] = {ReadProtected};

Options[DrawPhaseDiagram] = 
        Join[{PlotStyle -> Automatic, 
            CriticalPointStyle -> PointSize[0.02]}, 
          Options[Graphics] /. {(Frame -> _) -> (Frame -> {True, True, False, False})}];

DrawPhaseDiagram[seps : {{_Interval, _?NumericQ} .. }, 
        opts___?OptionQ] := 
      Module[{optsgr, ps, cs, splt, lins, cpts}, 
optsgr = 
          FilterOptions[Graphics, opts, 
            Sequence @@ Options[DrawPhaseDiagram]]; 
{ps, 
            cs} = {PlotStyle, CriticalPointStyle} /. 
              Flatten[{opts, Options[DrawPhaseDiagram]}] /. 
            Automatic -> 
              GrayLevel[0]; 

lins = 
          Flatten[Table[
                connectIntervals[
                  seps[[i]], 
                  seps[[i + 1]]], {i, 1, 
                  Length[seps] - 1}]] //. 
            joinLineSegments; 
splt = 
          Split[seps, (Length[#1[[1]]] == Length[#2[[1]]]) &]; 
        cpts = Table[
            criticalPoints[
              splt[[i, -1]], 
              splt[[i + 1, 1]]], {i, 1, 
              Length[splt] - 1}]; 

Show[
          Graphics[{{ps, lins}, {cs, cpts}}], 
          optsgr]
]

BinaryPhaseDiagram

Attributes[BinaryPhaseDiagram] = {HoldFirst, 
          ReadProtected};

Options[BinaryPhaseDiagram] = 
        Join[Options[BinaryPhaseSeparations], {VerticalPoints -> 25, 
            VerticalDivision -> 50. }, 
          Options[DrawPhaseDiagram]];

BinaryPhaseDiagram[f_, x_Symbol, p_List, opts___Rule] := 
      BinaryPhaseDiagram[f, {x, 0, 1}, p, opts]

BinaryPhaseDiagram[
        f_, {x_Symbol, xmin_?NumericQ, 
            xmax_?NumericQ} /; (0 ≤ xmin < xmax ≤ 
              1), {p_Symbol, pmin_?NumericQ, 
            pmax_?NumericQ} /; (pmin < pmax), opts___Rule] := 
      Module[{optsbps, optsdpd, pr, vp, vd, ps, seps = {}, 
          srt}, 
optsbps = 
          FilterOptions[BinaryPhaseSeparations, opts, 
            Sequence @@ 
              Options[BinaryPhaseDiagram]]; 
optsdpd = 
          FilterOptions[DrawPhaseDiagram, opts, 
            Sequence @@ 
              Options[BinaryPhaseDiagram]]; 
pr = PlotRange /. {optsdpd} /. ({Automatic, 
                      pr_} :> {{xmin - (xmax - xmin)/40. , 
                        xmax + (xmax - xmin)/40}, pr}) /. ({xr_, 
                    Automatic} :> {xr, {pmin - (pmax - pmin)/
                          40. , pmax + (pmax - pmin)/
                          40. }}) /. (Automatic -> {{xmin - (xmax - xmin)/40. , 
                    xmax + (xmax - xmin)/40. }, {pmin - (pmax - pmin)/
                        40. , pmax + (pmax - pmin)/
                        40. }}); 
{vp, 
            vd} = {VerticalPoints, VerticalDivision} /. 
            Flatten[{opts, 
                Options[
                  BinaryPhaseDiagram]}]; 

Plot[
(
ps = 
              BinaryPhaseSeparations[With[{p = y}, f], x, 
                optsbps]; 
seps = 
              Join[seps, {{ps, y}}]; 
If[ps === Interval[], 
              0, Min[(#[[2]] - #[[1]]) & /@ (List @@ 
                      ps)]]
), 
{y, 
            pmin, pmax}, DisplayFunction -> Identity, 
          PlotPoints -> vp, 
          PlotDivision -> 
            vd]; 

srt = 
          Sort[Union[
              seps], (#1[[2]] < #2[[2]]) &]; 

{srt, 
          DrawPhaseDiagram[srt, PlotRange -> pr, 
            optsdpd]}
]

Definitions for system functions

n/a

Restore protection of system symbols

n/a

End the private context

End[];

Epilog

This section protects exported symbols and ends the package.

Protect exported symbols

Protect[Evaluate[$Context <> "*"]]

End the package context

EndPackage[]

Examples, Tests

Disable spelling error messages and set options for Plot.

Off[General::spell, General::spell1]

SetOptions[Plot, Axes -> False, 
        Frame -> {True, True, False, False}];

This reads in the package.

<< PhaseSeparations`

1. Regular solution

f1[x_, t_] := 
      t (x Log[x] + (1 - x) Log[1 - x]) + 
        2  x (1 - x)

Plot[f1[x, 0.8], {x, 0, 1}]

BinaryPhaseSeparations[f1[x, 0.8], x]

A phase diagram may be constructed manually, using BinaryPhaseSeparations and DrawPhaseDiagram...

seps = 
      Table[{BinaryPhaseSeparations[f1[x, p], x], p}, {p, 0.6, 1.1, 
          0.05}]

DrawPhaseDiagram[seps]

... or automatically, using BinaryPhaseDiagram:

BinaryPhaseDiagram[f1[x, t], x, {t, 0.6, 1.1}] // Short

2. Multiple critical points & three-phase lines

More complicated free energy functions can produce phase diagrams with more than one critical point and/or splitting of a single phase coexistence region into two.

f2[x_, t_] := 
      t (x Log[x] + (1 - x) Log[1 - x]) + 2  x^2 (1 - x) - 
        500 (x (1 - x))^6

BinaryPhaseDiagram[f2[x, t], x, {t, 0.1, 2. }];

The near-critical regions of the phase diagrams may be better resolved by increasing the value of the SepPoints option. Plotting of the critical points may be suppressed by setting their style to a PointSize of zero.

BinaryPhaseDiagram[f2[x, t], x, {t, 0.1, 2. }, SepPoints -> 200, 
        CriticalPointStyle -> PointSize[0]];