BCC Alloys in the Bragg-Williams Approximation

The Model

Samuel M. Allen	April 19, 2006

Clear["Global`*"]
<< Graphics`Legend`

Define values for Agogadro and Boltzmann constants

N0 = 6.02214×(10^23);
kB = 1.3806503×(10^(-23));

Bragg-Williams model for b.c.c. alloys with first-neighbor pairwise interactions.

Order parameters are the atom fraction of B, XB, and the long-range order parameter, η (see KoM Fig. 17.4 and Eqs. 17.7–8 for definitions).

Enter expressions for internal energy,entropy, and free energy as in KoM Eqs. 17.13–14. In the formuation below ΔF has units of energy/mole

ΔU[
          XB_, η_] := -4 N0 W (XB (1 - 
                    XB) + (η^2));
ΔS [
          XB_, η_] := -N0 kB ((XB + η) 
                Log[XB + η] + (XB - η) 
                Log[XB - η] + (1 - XB + η) 
                Log[1 - XB + η] + (1 - XB - η) 
                Log[1 - XB - η]);
ΔF[XB_, η_, T_] := ΔU[XB, η] - 
        T ΔS[XB, η]

Look first at the possibility for order-disorder transitions by calculating the order-disorder transition temperature at XB = 0.5:

SecDerivEta = (∂_(η, η))ΔF[
          XB, η, T]
Fetaeta = Limit[SecDerivEta, η -> 0]
Fetaeta /. XB -> 0.5

We see that there will only be order-disorder transitions (i.e., a positive-temperature solution to 
	Limit(Subscript[((Subscript[F, η])η), η -> 0])=0 
when W>0. With this result, we can explore the simpler case of phase equilibria with miscibility gaps by considering only W<0 and setting η=0 in the free energy model. (In this way, we recover the regular solution model.)

Spinodal Decomposition

Calculation of Miscibility Gap and Chemical Spinodal in the Bragg-Williams Approximation

We'll need to make calculations and draw graphs so I begin by picking a suitable negative value for W

W = -1000×kB

Begin by evaluating the spinodal temperature at XB = 0.5

(∂_(XB, XB))ΔF[XB, 0, 
        T]
FXBXB = 
      Limit[(∂_(XB, XB))ΔF[XB, 0, T], 
        XB -> 0.5]

Solve[FXBXB == 0, T]

So the top of the miscibility gap (remember, this free-energy model is symmetric about XB = 0.5) is at T = 1000 K.

Plot ΔF vs. XB for T = 800, 900, 1000, 1100 K

Plot[{ΔF[XB, 0, 800], ΔF[XB, 0, 
          900], ΔF[XB, 0, 1000], ΔF[XB, 0, 
          1100]}, {XB, 0.0, 1}, AxesLabel -> {XB, ΔF}, 
      PlotStyle -> {{Hue[ .3]}, {Hue[ .5]}, {Hue[ .7]}, {Hue[ .9]}}, 
      PlotLegend -> {"800 K", "900 K", "1000 K", "1100 \
K"}, LegendPosition -> {1, -1}]

Here is a 3D representation of the free energy surface ΔF(XB,T):

Plot3D[ΔF[XB, 0, T], {XB, 0.01, 0.99}, {T, 300, 1500}, 
      AxesLabel -> {XB, T, ΔF}]

Mathematica's FindRoot function can be used to map out the miscibility gap. Note that the common tangent in this symmetric free energy model will have zero slope. FindRoot needs an initial guess and I start from 990 K and work downward along one branch of the curve.

First, initialize an array for storing pairs of values (T, XB) that will form the coexistence curve.

CoexCurvePts = {}

Here's the general procedure for calculating and extracting the value of a point on the coexistence curve:

T = 990
XBCoex = 
      FindRoot[(∂_XB) ΔF[XB, 0, T] == 0, {XB, 
          0.6}]
XBCoex[[1, 2]]

Now calculate a whole lot of points:

CoexCurvePts = {};
Do[
XBCoex = 
        FindRoot[(∂_XB) ΔF[XB, 0, T] == 0, {XB, 
            XBCoex[[1, 2]]}]; 
CoexCurvePts = 
        Append[CoexCurvePts, {XBCoex[[1, 2]], T}], {T, 990, 
        300, -10}]

BSide = 
      ListPlot[CoexCurvePts, 
        AxesLabel -> {"Atom Fraction,XB", "Temperature,K"}]

Now do the same for the A-rich compositions, adding the points to CoexCurvePts. Note I'm adding the point at the top of the miscibility gap as well. Also, I could have made use of the fact that the coexistence curve is symmetric about XB = 0.5, but I calculated the points instead.

CoexCurvePts = 
        Prepend[CoexCurvePts, {0.5, 1000}];
XBCoex[[1, 2]] = 0.4;
Do[
XBCoex = 
        FindRoot[(∂_XB) ΔF[XB, 0, T] == 0, {XB, 
            XBCoex[[1, 2]]}]; 
CoexCurvePts = 
        Prepend[CoexCurvePts, {XBCoex[[1, 2]], T}], {T, 990, 
        300, -10}]

Here is the calculated equilibrium diagram

CoexPlot = 
      ListPlot[CoexCurvePts, PlotJoined -> True,  
        AxesLabel -> {"Atom Fraction,XB", "Temperature,K"}, 
        PlotStyle -> {Thickness[0.005]}]

The same basic process can be used to calculate and plot the spinodals. I collect the two branches of the curve in separate arrays.

SpinodalCurvePtsR = {};
SpinodalCurvePtsL = {};
Do[
XBSpin = 
        Flatten[Solve[(∂_(XB, XB))ΔF[XB, 0, T] == 
              0, {XB}]]; 
SpinodalCurvePtsR = 
        Append[SpinodalCurvePtsR, {XBSpin[[2, 2]], 
            T}]; 
SpinodalCurvePtsL = 
        Append[SpinodalCurvePtsL, {XBSpin[[1, 2]], T}], {T, 990, 
        300, -10}]

Now add in the point at the top of the spinodal curve, join the two branches of the curve into a single ordered list, and plot the results.

SpinodalCurvePts = 
        Prepend[SpinodalCurvePtsR, {0.5, 1000}];
SpinodalCurvePts = 
        Join[Reverse[SpinodalCurvePtsL], 
          SpinodalCurvePts];
SpinodalPlot = 
      ListPlot[SpinodalCurvePts, 
        AxesLabel -> {"Atom Fraction,XB", "Temperature,K"}, 
        PlotJoined -> True]

Here is a combined equilibrium showing the coexistence curve and the spinodal superposed.

Show[CoexPlot, SpinodalPlot]

Diffuse Interfaces in the Bragg-Williams approximation

The gradient-energy coefficient Kc is related to W and ΔU and the lattice constant a of the crystal. Assuming a reasonable value for the lattice constant I calculate Kc for the system as

a = 0.29×(10^(-9));
VMolar  =  ((a^3)/2) N0;
VAtomic  =  ((a^3)/2);
ΔU50 = ΔU[0.5, 
          0];
Kc = (2/3) ΔU50 (( ((((Sqrt[3])/2  a))^2))/VMolar)

We will look at the energy of the interface between equilibrium phases at T = 900 K. We need the expression for the free energy density Δfhom(XB) as illustrated in KoM Fig. 18.1. First calculate the free energy density of the equilibrium phases. Note again that ΔF is symmetric about XB = 0.5 and the common tangent has zero slope.

T = 900;
αRoot = 
        FindRoot[(∂_XB) ΔF[XB, 0, T] == 
            0, {XB,  .2}];
XBEqα = αRoot[[1, 2]]
XBEqβ = 1 - XBEqα;
ΔFEq[T_] = ΔF[XBEqα, 0, 
        T]

Now calculating Δfhom(XB) on a per unit volume basis (units are thus Joule/(m^3))

ΔfHom[XB_, 
        T_] := ((ΔF[XB, 0, T] - ΔFEq[T])/VMolar)

\(Plot[ΔfHom[XB, T], {XB, XBEqα, XBEqβ}, 
      AxesLabel -> {"Atom Fraction,XB", "ΔfHom"}]\)

The excess free energy per unit area of the interface in J/(m^2) is

γ = 
      2  NIntegrate[(Sqrt[ Kc ΔfHom[XB, T]]), {XB, 
            XBEqα, XBEqβ}]

This is reasonable--about 20 mJ/(m^2).

Spinodal decomposition kinetics

Assume an Arrhenius law for the mobility Mc

Mc[T_] := A0 Exp[-(Q/(N0 kB T))]

Q =  200000;
A0 = (10^(-18));

Clear[FXBXB]
FXBXB[T_] := 
        Limit[(∂_(XB, XB))ΔF[XB, 0, T], 
          XB -> 0.5];

So the interdiffusivity at 900 K for a solution with XB = 0.5 will be

T = 900;
DTilde[T_] = Mc[T] FXBXB[T]

The amplification factor given by KoM Eq. 18.47 becomes

R[β_, 
        T_] :=  -Mc[T] (FXBXB[T]/VMolar  +  
            2 Kc (β^2)) (β^2)

Plot[R[β, T], {β, 0, (10^10)}, 
      AxesLabel -> {"Wavenumber, β", "Amplification Factor, \
R(β)"}]

Find kinetic and thermodynamic wavelengths at 900 K

SolMax = Solve[(∂_q) R[q] == 0, q]

λMax = (2  π/q) /. 
        SolMax[[3]]
RMax = R[q] /. SolMax[[3]]

1 / RMax would be a good estimate of the "characteristic time" that it would take a spinodal microstructure to develop.  In this case it would be about 200 s.

SolZero = FindRoot[R[p] == 0, {p, 4×(10^9)}]

λCrit = (2  π/p) /. SolZero

So at T = 900 K we find λcrit = 2.04 nm and  λmax = 2.88 nm.

Now plot the R(β,T) surface to get an idea of the temperature dependence of βMax and R(βMax)

Plot3D[R[β, T], {β, 0, 5×(10^9)}, {T, 800, 1000}, 
      PlotRange -> { .01, -0.01}, PlotPoints -> 50, 
      ViewPoint -> {-1.3, 2.4, 1. }, 
      AxesLabel -> {"β", "T", "R(β)"}]

Note the following trends:
As the undercooling below the spinodal temperature of 1000 K increases
	(1) The kinetic wavenumber βMax increases, so the kinetic wavelength decreases
	(2) There is a temperature at which the amplification rate R(βMax) is a maximum. In the plot above this occurs around 930 K.