Diffusional Growth of a Spherical Precipitate from Supersaturated Solution

Samuel M. Allen	April 28, 2006

Clear["Global`*"]

The Problem

Let's model the growth kinetics in an Ni– 48 at.% Cr alloy that has been solution treated by annealling at 1573 K and quenched to 1273 K where it nucleates a Cr-rich particle.The Ni–Cr equilibrium diagram appears below.

\(At\ 1573\ K\ the\ alloy\ composition\ is\ \((c\^α∞ = 0.48)\)\)
\(At\ 1273\ K\ the\ Ni\ phase\ composition\ is\ \((c\^αβ = 0.44)\)\)
\(At\ 1273\ K\ the\ Cr\ phase\ composition\ is\ \((c\^βα = 0.93)\)\)

Assume growth of a single spherical β particle in an infinite matrix and follow the formulation in KoM Section 20.2.1, page 513.

Show[Import["CrNi.tiff"]]

cα∞ = 0.48;
cαβ = 0.44;
cβα = 0.93;

The Solution

The solution to this growth problem is governed by the supersaturation parameter \(\(cα∞ - cαβ\)\/\(cβα - cαβ\)\):

SaturationPar = ((cα∞ - cαβ)/(cβα - cαβ))

Solve Eq. 20.45 for ηR

Sol = 
      FindRoot[η(R^3) ==  (SaturationPar/2) 
            Exp[-η(R^2)]/((1/ηR) Exp[-η(R^2)] - (Sqrt[π]) Erfc[ηR]), {ηR, 0.5}]

EtaR = Sol[[1, 2]]

Particle Size and Growth Rate

The particle growth is parabolic, given by KoM Eq. 20.40

R[t_] := EtaR (Sqrt[4 Dα t])

Let the diffusivity Dα = (10^(-16)) (m^2)/s.
Calculate the particle radius and growth rate of a particle that has been growing for 100 s:

tGrowth = 100;
Dα = (10^(-16));
R[tGrowth]

Vel[t_] := (∂_t) R[t]

Vel[t] /. t -> tGrowth

So at 1273 K a particle grows to about 50 nm in radius in 100 s, and the growth rate at that time is about 2.5 Ångstrom/s.

Concentration profile in Matrix

The concentration is given as a function of the scaling parameter, η, by KoM Eq. 20.43

cOfEta[η_] := 
      cα∞ - (cα∞ - 
              cαβ) (((1/η) Exp[-(η^2)] - (Sqrt[π]) Erfc[η])/((1/EtaR) Exp[-(EtaR^2)] - (Sqrt[π]) Erfc[EtaR]))

Plot[cOfEta[η], {η,  .001, 1}, PlotRange -> {0, 1}, 
      AxesLabel -> {"η", "c"}, 
      PlotLabel -> "Composition vs. Scaling Variable"]

The c(η) plot is a "master plot" that characterizes all diffusional growth at a given supersaturation parameter. From it, one can compute c(r) for any growth time. An example follows for a particle that has grown for 100 s.

cOfR[R_, t_] := 
      cOfEta[η] /. η -> r/(Sqrt[4 Dα t])

P1 = 
      Plot[cOfR[r, tGrowth], {r, R[tGrowth], 0.5×(10^(-6))}, 
        PlotRange -> {0.4, 0.5}, AxesLabel -> {"r", "c"}, 
        PlotLabel -> "Composition vs. Radius (m)"]

Comparison with Solution via Invariant-Field Approximation

Recall that the invariant-field approximation involves assuming that the composition profile is given by the solution to Laplace's equation (∇^2) c = 0, with the boundary concentrations being c(∞) = cα∞ and c(RInterface) = cαβ. Hence

cIFA[r_] := 
      cα∞ - (((cα∞ - cαβ) R[tGrowth])/r)

P2 = 
      Plot[cIFA[r], {r, R[tGrowth], 0.5×(10^(-6))}, 
        PlotRange -> {0.4, 0.5}, AxesLabel -> {"r", "c"}, 
        PlotLabel -> "Composition vs. Radius (m)", 
        PlotStyle -> Dashing[{ .01,  .01}]]

Now superpose the two plots for a comparison

Show[P1, P2]

The invariant-field approximation is seen to underestimate the composition at a given radius, and hence overestimate the spatial extent of the diffusion field around the particle.

The particle growth rate is also underestimated in the invariant-field approximation:

VelIFA[
        t_] := ((Dα (∂_r) cIFA[r])/(cβα - cαβ))

VelIFA[t] /. {r -> R[tGrowth], t -> tGrowth}

The invariant-field approximation is valid in the limit of small supersaturation (in this example the supersaturation is about 8%). I've written this notebook so it is easy to adjust the supersaturation and explore just how small the supersaturation needs to be for a given error in the composition profile or the interface velocity.

Spherical Harmonics

Samuel M. Allen	May 2, 2006

Clear["Global`*"]
<< Graphics`Graphics3D`
<< Calculus`VectorAnalysis`
SetCoordinates[Spherical[r, θ, ϕ]];

Spherical harmonics arise from the solution to Laplace's equation in spherical coordinates. They are useful for describing shape perturbations of a sphere, and they take two parameters that relate to the number of "bumps" that distort the sperical shape.

Consider first a subset of spherical harmonics that describe the shape of a circle that has its radius perturbed by a periodic shperical harmonic function. Here the parameter n corresponds to the "wavelength" of the perturbation of the circle and the second parameter m is set to zero.

SHmzero[n_, θ_] := 
      SphericalHarmonicY[n, 0, θ, 0]

Here is a perturbation with 5 maxima

SHmzero[5, θ]

In the plots below the undistorted circle/sphere has unit radius. The variable "amp" sets the amplitude of the perturbation.

Amp =  .2;
n = 4;
ParametricPlot[{(1 + Amp SHmzero[n, θ]) 
          Cos[θ], (1 + Amp SHmzero[n, θ]) 
          Sin[θ]}, {θ, 0,  2  π}, 
      AspectRatio -> 1]

Now look at perturbations of a sphere. This routine plots an unperturbed sphere:

ParametricPlot3D[{Cos[θ] Cos[ϕ], 
        Sin[θ] Cos[ϕ],  Sin[ϕ]}, {θ, 0, 
        2  π}, {ϕ, -π/2, π/2}, 
      PlotPoints -> 50]

The following routine allows you to input the parameters n and m to the shperical harmonic function and plot the shape of a sphere that is modified by this type of perturbation.

First I plot n = 4 and m = 0, for which the horizontal mid-cross section should correspond to the perturbed circle illustrated above.

Amp = 0.2;
n = 4; m = 0; 
    ParametricPlot3D[{Abs[(1 + 
                Amp SphericalHarmonicY[n, m, θ, ϕ])] 
          Cos[θ] Cos[ϕ], 
        Abs[(1 + Amp SphericalHarmonicY[n, m, θ, ϕ])] 
          Sin[θ] Cos[ϕ], 
        Abs[(1 + Amp SphericalHarmonicY[n, m, θ, ϕ])] 
          Sin[ϕ]}, {θ, 0, 
        2  π}, {ϕ, -π/2, π/2}, 
      PlotPoints -> 50]

Now I show the n = 4 and m = 2 distortion, for comparison with the previous plot. Note the vertical "broadening" of the lobes in the 4, 2 harmonic function.

Amp = 0.2;
n = 4; m = 2; 
    ParametricPlot3D[{Abs[(1 + 
                Amp SphericalHarmonicY[n, m, θ, ϕ])] 
          Cos[θ] Cos[ϕ], 
        Abs[(1 + Amp SphericalHarmonicY[n, m, θ, ϕ])] 
          Sin[θ] Cos[ϕ], 
        Abs[(1 + Amp SphericalHarmonicY[n, m, θ, ϕ])] 
          Sin[ϕ]}, {θ, 0, 
        2  π}, {ϕ, -π/2, π/2}, 
      PlotPoints -> 50]