316-2: Microstructural Dynamics II

D. Joester, K.R. Shull, P.W. VoorheesDepartment of Materials Science and EngineeringNorthwestern University

Table of Contents

Catalog Description (316-1,2)

Principles underlying development of microstructures. Defects, diffusion, phase transformations, nucleation and growth, thermal and mechanical treatment of materials. Lectures, laboratory. Prerequisite: 315 or equivalent.

Course Outcomes

At the conclusion of 316-2 students will be able to:
  1. Predict nucleation rates from thermodynamic data
  2. Describe where precipitates are likely to form in a multicomponent material
  3. Design processing histories to obtain a desired microstructure
  4. Correctly use and interpret TTT diagrams

3 Introduction

Much of this class is concerned with the appearance and growth of a new phase, referred to as β in Figure 3.1, from a matrix phase, α . In many cases the process can be divided into the following two phases:
  1. A nucleation phase, where individual very small (typically in the nanometer range) β precipitates are observed.
  2. A growth phase where the β precipitates grown in size.
image: 4_home_ken_Mydocs_MSEcore_316-2_figures_drawing31.svg Figure 3.1: Schematic Representation of Nucleation and Growth.
The process is actually much more interesting than one might think by looking at the simple example shown in Figure 3.1. Consider, for example, the images of snowflakes shown in figure 3.2.
image: 5_home_ken_Mydocs_MSEcore_316-2_figures_snowflake1.png   image: 6_home_ken_Mydocs_MSEcore_316-2_figures_snowflake2.png   image: 7_home_ken_Mydocs_MSEcore_316-2_figures_snowflake3.png
image: 8_home_ken_Mydocs_MSEcore_316-2_figures_snowflake4.png   image: 9_home_ken_Mydocs_MSEcore_316-2_figures_snowflake5.png   image: 8_home_ken_Mydocs_MSEcore_316-2_figures_snowflake4.png
Figure 3.2: Snowflake images .
image: 10_home_ken_Mydocs_MSEcore_316-2_figures_snowmorphology.png
Figure 3.3: Morphology diagram for snowflakes .
If you are not motivated by snowflakes, plenty of modern technological examples exist that illustrate the concepts developed in this course. Silicon nanowires illustrated below are one excellent example.
image: 11_home_ken_Mydocs_MSEcore_316-2_figures_nanowire_forest.png   image: 12_home_ken_Mydocs_MSEcore_316-2_figures_nanowire_falsecolor.png
Figure 3.4: Two examples of silicon nanowires .
image: 13_home_ken_Mydocs_MSEcore_316-2_figures_eutectic_alpha_ppt.svg
Figure 3.5: Generic eutectic phase diagram. A solid phase appears from the liquid as the temperature is cooled along the path indicated by the arrow.
The simplest example of these concepts is the freezing of pure water. The equilibrium phase at a given temperature is the one with the lowest molar free energy, G m . A schematic representation of these free energies is shown in Figure 3.6. The free energies cross at 0 C, with the liquid having a lower free energy at higher temperatures and the solid having a lower free energy at lower temperatures.
image: 14_home_ken_Mydocs_MSEcore_316-2_figures_liquid_and_solid_free_energy.svg
Figure 3.6: Schematic representation of the molar free energies of the solid and liquid phases of water.
Of course if we cool water below 0 °C, ice is not formed immediately - the solidification process takes time. Under appropriate conditions, it can be possible to maintain a liquid below its equilibrium melting point indefinitely, without ever seeing a solid. The ability to supercool a liquid is related to the fact that the interface between the solid and the liquid contributes a positive contribution to the overall free energy of the system.
Thermodynamics can be used to calculate the properties of equilibrium phases, including unstable equilibrium phases like the critical nucleus. The difference between a stable equilibrium (absolute minimum in free energy), a metastable equilibrium (local minimum in the free energy) and an unstable equilibrium (local maximum in the free energy) is illustrated in Figure 3.7.
image: 15_home_ken_Mydocs_MSEcore_316-2_figures_drawing1.svg
Figure 3.7: Energy as a function of nucleus size, illustrating stable, metastable and unstable states.
We need to compare the appropriate thermodynamic potentials for the cases where a critical nucleus is present and where it is absent. As illustrated schematically in Figure 3.8, we assume that the critical nucleus is spherical, with a radius of R * .
image: 16_home_ken_Mydocs_MSEcore_316-2_figures_drawing32.svg
Figure 3.8: Initial and final states for the calculation of the work to form a critical nucleus.
Because of Laplace pressure, the pressure inside the critical nucleus is not equal to the pressure outside the critical nucleus. The correct thermodynamic potential to use is called the omega potential. The omega potential is relevant when a region of fixed volume is transformed from one phase to the other, in conditions where atoms are freely exchanged between the two phases. We consider this situation in more detail in the following section.

Homogeneous Nucleation in Unary Systems

Consider solidification of a pure compound from its melt as an example for a phase transformation L β in a unary system. Let's assume that the liquid phase L is homogeneous, i.e. there are no foreign objects or other surfaces in the system (such as container walls). We would like to understand the formation of a very small volume V β * of the solid phase, i.e. the nucleus, from the liquid. Note that we use β * to indicate the nucleus (Figure 4). This is to remind ourselves that the composition of the nucleus is not necessarily identical to that of the final phase β at equilibrium (this is for the general case; in a unary system the composition cannot change). We will further assume that the nucleus shall be spherical with radius R * .
image: 17_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_02_Fig1.svg
Figure 4.1: Homogeneous nucleation in a unary system. Initially, the system only consists of the parent phase, here the liquid phase L. After nucleation, one infinitesimally small, spherical nucleus β * is present also, surrounded by the parent phase L.
Before the formation of the nucleus, our system only contains L, at temperature T L and pressure P L . After, we also have to consider the temperature T β * and the pressure P β * of the nucleus.
The solidification shall be isothermal. Therefore, T L T β * (4.1)
What about the pressure? There is an interface between the nucleus and the liquid phase. The creation of this interface requires energy. As a consequence, the pressure in the nucleus is higher than in the liquid.
P β * > P L (4.2)
As an analogy, consider inflating a rubber balloon against atmospheric pressure. Because of the energy stored in the stretched rubber membrane, the pressure on the inside of the balloon is higher than that on the outside.
We therefore need to determine the pressure differential
Δ P= P β * - P L (4.3)
We will use the grand canonical (aka Landau, Omega) potential in the derivation of the free energy of the nucleus.
In our system, the grand canonical potential of the final state can be written as the sum Ω = Ω bulk β * + Ω exc + Ω L , (4.4) where Ω bulk β * is contribution from the bulk of the nucleus, Ω exc is the surface excess free energy, i.e. the contribution from the interface, and Ω L is the contribution from the liquid phase.
We can write the following: Ω exc = γ L β * A L β * , (4.5) where γ L β * is the interfacial free energy, with units [ γ ] = J m 2 , of the interface between the nucleus and the liquid phase, and A L β * is the interfacial area.
Recall that: Ω =-PV (4.6)
We can therefore write Eq. 4.4 as Ω =- P β * V β * - P L V L + γ L β * A L β * (4.7)
Generally speaking, we can see in Eq. 4.7 that the free energy of a small cluster (embryo) of the newly formed phase consists of negative terms that favor the formation, and positive terms that oppose the formation. The negative terms scale with the volume of the embryo, whereas the positive ones scale with the surface area. Therefore, if the surface-to-volume ratio is large, the embryo is unstable with respect to L. This means that it is much more likely to shrink than to grow. Conversely, if the surface-to-volume ratio is small, the embryo is stable with respect to L, and it is more likely to grow than to shrink. There is a critical size, at which the embryo is in unstable equilibrium, meaning that both shrinking and growing will reduce the free energy. This critical embryo is called the nucleus. In Figure 4.2, the free energy of an embryo of liquid water formed from water vapor is plotted against the radius. Note that this represents the condensation of a vapor rather than the solidification of a melt. Regardless, the thermodynamics are the same.
image: 18_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_02_Fig2.svg
Figure 4.2: Free energy change as a function of the radius of an embryo of the new phase (here: spherical embryo of liquid water condensing from water vapor at a supercooling of 10 K). Free energy is released as the volume of the embryo increases. This component scales with - R 3 (red line). However, energy is required to create the interface between the embryo and the parent phase. This term scales with + R 2 (green line). The latter term dominates at small R . Consequently, the sum of the two terms (blue line) goes through a maximum at the critical radius, R * .
Because the nucleus is in unstable equilibrium, . d Ω dR | R * =0 (4.8)
Performing this differentiation on Eq. 4.7, - P β * . d V β * dR | R * - P L . d V L dR | R * + γ L β * . d A L β * dR | R * + . d γ L β * dR | R * A L β * =0 (4.9)
We now make the very important assumption that the interfacial free energy γ does not change with the radius. In other words, that the interfacial free energy of the infinitesimally small nucleus is the same as that of a macroscopic interface. This assumption, known as the “capillarity assumption”, is convenient, and very likely wrong. How wrong, we do not know. So, with . d γ L β * dR | R * =0 , we write
- P β * . d V β * dR | R * - P L . d V L dR | R * + γ L β * . d A L β * dR | R * =0 (4.10)
The total volume of our system is constant. Therefore, any increase in the volume of the nucleus must equal the decrease in volume of the liquid (Figure 4.10).
image: 19_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_02_Fig3.svg
Figure 4.3: Volume conservation during nucleation. As the radius of the nucleus increases by dR , its volume increases by dV dR . Because the total volume of the system is constant, the volume of the parent phase L shrinks by the same amount.
This allows us to rewrite Eq. 4.10, - P β * . d V β * dR | R * + P L . d V β * dR | R * + γ L β * . d A L β * dR | R * =0 (4.11) -( P L - P β * ) . d V β * dR | R * = γ L β * . d A L β * dR | R * (4.12) Δ P = γ L β * . d A L β * dR | R * . d V β * dR | R * (4.13)
Using the assumption that the nucleus is spherical, we can write V β * ( R ) = 4 3 π R 3 . d V β * dR | R * =4 π R * 2 (4.14) A L β * ( R ) =4 π R 2 . d A L β * dR | R * =8 π R * (4.15)
With this,
Δ P= γ L β * 8 π R * 4 π R 2 = 2 γ L β * R * (4.16)
This pressure differential is called the Laplace pressure for the spherical nucleus (see also course notes for MSE 316-1). Note that for the derivation, we need not have assumed that the nucleus is spherical until solving Eqs. 4.15 and 4.15. For nuclei of different geometries, one could define a convenient critical length scale, L * , instead. The Laplace pressure for such a nucleus would then follow from evaluating the differentials in Eqs. 4.15 and 4.15 at L * .
Another way to evaluate the Laplace pressure for any interface is by considering the Laplace theorem Δ P= γ ( 1 R 1 + 1 R 2 ) (4.17) , where R 1 and R 2 are the principle radii of curvature.
Returning to the description of the nucleus, we can determine the reversible work of nucleation, W r * , as the difference between the free energy of the system at R= R * and the free energy at R=0 , where the volume and area of the nucleus vanish. W r * = Ω ( R= R * ) - Ω ( R=0 ) (4.18) =- P β * V β * + γ A- P L V L -( - P L V 0 L ) , (4.19) where V 0 L = V L + V β * is the volume of L at before a nucleus forms, i.e. at R=0 . Plugging in, and collecting terms in two steps W r * =- P β * V β * + γ A- P L ( V L -( V L + V β * ) ) (4.20) =-( P β * - P L ) V β * + γ A (4.21) =- Δ P V β * + γ A (4.22)
Assuming again that the nucleus is spherical, W r * =- Δ P( 4 3 π R *3 ) + γ ( 4 π R *2 ) (4.23)
Using Eq. 4.16, we write the critical radius R * as a function of the Laplace pressure R * = 2 γ Δ P (4.24)
substitute into Eq. 4.23, and simplify W r * =- Δ P( 4 3 π ( 2 γ Δ P ) 3 ) + γ ( 4 π ( 2 γ Δ P ) 2 ) (4.25) =- 32 π 3 γ 3 Δ P 2 + 48 π 3 γ 3 Δ P 2 (4.26) W r * = 16 π 3 γ 3 Δ P 2 (4.27)
With this, we have determined the critical radius and reversible work of the spherical nucleus as a function of the Laplace pressure. It would be more convenient if we could replace this pressure, which is difficult to measure, by some other expression that we can evaluate. We will do this in the next section.
Recall that the nucleus is in unstable equilibrium with the parent phase, i.e. the phase that it formed from. At some temperature T , we can therefore write the following equality G m β * ( T, P β * ) = G m L ( T, P L ) (4.28)
Let's consider a phase transformation that occurs as the temperature is lowered, such as a solidification from the melt. Let Δ T= T m -T be the supercooling below the melting point T m . We expect that the solidification is spontaneous if Δ T>0 . Recall also that Δ P= P β * - P L .
The molar Gibbs free energies in Eq. 4.28 are energy surfaces of two variables. To simplify things, let's approximate G m as a linear function in both T and P . We do this by writing the Taylor expansion of G m and terminating it after the 1 st order terms. Because we know something about the equilibrium state, let's develop the Taylor expansion around the equilibrium point ( T m , P L ) .
G m ( T,P ) G m ( T m , P L ) + . G m T | T m ( T- T m ) + . G m P | P L ( P- P L ) (4.29)
Recall that G m T =- S m (4.30) G m P = V m (4.31)
Let's now write the l.h.s. and r.h.s. of Eq 4.28 in terms of the Taylor expansion given in Eq. 4.29, and substitute the thermodynamic relationships given in Eqs. 4.31 and 4.31: G m β * ( T, P β * ) = G m β * ( T m , P L ) - S m β * ( T- T m ) + V m β * ( P β * - P L ) (4.32) G m L ( T, P L ) = G m L ( T m , P L ) - S m L ( T- T m ) + V m L ( P L - P L ) (4.33)
Note that by using Eq. 4.31, we assume that the molar volume V m is a constant, and does not change with changing pressure. This assumption of an incompressible nucleus is oftentimes reasonable in condensed matter, but can fail miserably when considering vapor phases.
Note that the third term on the r.h.s. of Eq. 4.33 vanishes. Further, because ( T m , P L ) is the point of equilibrium, the first term on the r.h.s. of Eq. 4.33 is equal to the first term on the r.h.s. of Eq. 4.33. Finally, with the definitions for Δ P and Δ T , we rewrite Eq. 4.28 as:
S m β * Δ T+ V m β * Δ P= S m L Δ T (4.34) V m β * Δ P=( S m L - S m β * ) Δ T (4.35)
Note that S m L - S m β * = Δ S m (4.36) where Δ S f,m is the molar standard entropy of fusion (melting), i.e. the entropy gain associated with the opposite of the phase transformation that we are considering. This is simply the way the signs work out. When considering other phase transformations, for example melting, pay particular attention to getting the signs right in this step.
We substitute Eq. 4.36 into Eq. 4.35 and rearrange: Δ P= Δ S m V m β * Δ T (4.37)
Recall that at equilibrium, Δ G m ( T m ) =0 (4.38) Δ H m - T m Δ S m =0 (4.39) Δ S m = Δ H m T m (4.40)
This allows us to rewrite Eq. 4.37: Δ P= Δ H f,m V m β * T m Δ T (4.41)
Equations 4.37 and 4.41 therefore allow us to estimate the Laplace pressure for any unary system for which we know the molar standard entropy of fusion, or the molar standard enthalpy of fusion and the melting temperature. These values are available for many pure phases. We can now express the critics radius and the reversible work of nucleation for a spherical, incompressible nucleus, using the capillarity approximation, as follows:
R * = 2 γ Δ P (4.42) = 2 V m β * Δ S m γ Δ T (4.43) = 2 V m β * T m Δ H f,m γ Δ T (4.44) W r * = 16 π 3 γ 3 Δ P 2 (4.45) = 16 π 3 ( V m β * Δ S f,m ) 2 γ 3 Δ T 2 (4.46) = 16 π 3 ( V m β * T m Δ H f,m ) 2 γ 3 Δ T 2 (4.47)
How do the Laplace pressure, critical radius, and the reversible work of nucleation depend on supercooling? Consider their graphs in Figure 4.4.
image: 20_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_03_Fig1.svg
Figure 4.4: Scaling of the Laplace pressure Δ P , the critical radius R * , and the reversible work of nucleation W r * , with supercooling Δ T . Graphs shown here correspond to the formation of a nucleus of liquid water from water vapor. A) Note the linear increase of the Laplace pressure with supercooling. At Δ T=0 , i.e. T= T b , Δ P=0 . This means that the interface is flat (infinite radius of curvature). B) The critical radius is inversely proportional to Δ T . From this follows that small radii are expected at high supercooling, i.e. for temperatures far below the equilibrium temperature (here: T b ). At low supercooling, the critical radius increases quickly and approaches infinity as Δ T0 . An infinite radius is equivalent to a flat interface, but also means that nucleus is the size of the entire system. C) The reversible work is inversely proportional to Δ T 2 and shows the same qualitative behavior as the critical nucleus. Because of the higher exponent, the effects are more pronounced. Close to the equilibrium temperature, W r * , which is an energetic barrier, becomes very large. Nucleation is therefore unlikely/very slow. The barrier quickly drops with increasing supercooling. In (D), the reversible work is expressed in units of kT .
As an example, consider the condensation of water vapor at Δ T=10 K . In this case, the parent phase is the vapor phase (vap), and the nucleus is the liquid phase L. In Eq. 4.47 and Eq. 4.47, we further have to replace the entropy of fusion by the entropy of vaporisation. Keep in mind that we also replace the melting point by the boiling point.
R * = 2 V m L Δ S v ap,m γ Δ T (4.48) W r * = 16 π 3 ( V m L Δ S v ap,m ) 2 γ 3 Δ T 2 (4.49)
We find the following values in the literature: γ =72 m J m 2 , Δ S vap,m =118.19 J m olK , T b =373.16 K . The molar volume of the nucleus at the boiling point can be determined from the density of water at 372.16 K, ρ 3 73.16K H 2 O =0.95805 g c m 3 , and the molecular weight of water, M w =18 g m ol . We find the molar volume as: V m L = M w ρ 3 73.16K H 2 O = 18 g m ol 0.95805 g c m 3 = 181 0 -6 0.95805 m 3 m ol 1.8791 0 -5 m 3 m ol (4.50)
Plugging in, R * = 2 V m L Δ S v ap,m γ Δ T (4.51) = 21.8791 0 -5 m 3 m ol 118.19 J m olK 0.072 J m 2 10 K (4.52) =2.28911 0 -9 m (4.53) 2.3 n m (4.54) (4.55) W r * = 16 π 3 ( V m L Δ S v ap,m ) 2 γ 3 Δ T 2 (4.56) = 16 π 3 ( 1.8791 0 -5 m 3 m ol 118.19 J molK ) 2 ( 0.072 J m 2 ) 3 ( 10 K ) 2 (4.57) =2.19541 0 -17 J (4.58) 2.201 0 -17 J (4.59)
Let's also calculate the Laplace pressure for this case. Plugging into Eq. 4.37: Δ P = Δ S f,m V m β * Δ T (4.60) = 118.19 J m olK 1.8791 0 -5 m 3 m ol 10 K (4.61) =6.29001 0 7 J m 3 (4.62) 6.291 0 7 N m m 3 =62.9 M Pa (4.63)
Let's run a quick sanity check on the assumption that the nucleus is incompressible. The bulk modulus of water is K=2.15 GPa . Using the definition of K and assuming linear bulk elasticity
K =- dP dV V 0 (4.64) K =- Δ P Δ V V 0 (4.65) Δ V V 0 =- Δ P K (4.66) (4.67)
Plugging in, we find the volumetric strain Δ V V 0 =- 62.9 M Pa 2.15 G Pa (4.68) =-0.0293 (4.69) -3% (4.70)
Therefore, we would expect that assuming that the nucleus is incompressible, and ignoring the dependency of V m on pressure will introduce only a small error into our calculation. Obviously, if the bulk modulus is higher, or the Laplace pressure lower, the error will have less of an impact.

Kinetics of Homogeneous Nucleation in Unary Systems

Microscopically, one of the assumptions of classical nucleation theory is that embryos of the new phase grow by addition of one monomer unit at a time, and shrink by dissociation of one monomer unit. In other words, embryos do not interact with each other. For an embryo with n monomer units, we can define forward and backward rates for these processes, as shown in Figure 1. At any given time t in our system, we can then count the number of embryos of size n and divide by the total volume to get the number density ρ ( n,t ) .
image: 21_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_04_Fig1.svg
Figure 5.1: Embryos are assumed to grow by addition, and shrink by dissociation of one monomer unit. The monomer unit could be a metal atom, a water molecule, or a protein. This is reasonable as long as the density of embryos is so low that they do not encounter one another. For each embryo comprised of n monomer units, a forward (growth) rate f( n ) and a backward (shrinkage) rate b( n ) is defined.
The number of embryos of size n that are formed per unit time and volume, can the be expressed as J( n ) = ρ ( n-1 ) A( n-1 ) β ( n-1 ) - ρ ( n ) A( n ) α ( n ) , (5.1) where A is the surface area of the embryo, β is the total flux of monomers towards the embryo, and α is the evaporative flux, away from the embryo. These fluxes depend on the forward and backward rates.
Becker and Döring used this approach to determine the change in the number density of clusters over time and found that
ρ ( n,t ) ( t ) =J( n ) -J( n+1 ) (5.2)
While a detailed analysis of this problem is beyond the scope of this class, a useful outcome is that for many systems relevant to us, it can be shown that a steady state is reached within 1 μ s . Without going into detail, we will use the following formula for the steady state nucleation current N ˙ V =Z ν A nuc ρ e - W r * kT , (5.3) [ N ˙ V ] = 1 sc m 3 (5.4) where Z is the dimensionless Zeldovich factor, i.e. the probability that a nucleus will become supercritical, ν is the impingement rate of monomers on the surface of the nucleus ( [ ν ] = 1 s c m 2 ), A nuc is the surface area of the nucleus ( [ A nuc ] = c m 2 ), and ρ is the number density of monomers in the parent phase ( [ ρ ] = 1 c m 3 ). Because of the strong dependence of the reversible work of nucleation, and therefore the entire exponent, on temperature, we can frequently approximate the prefactor as a constant independent of temperature. In many unary (!!) condensed matter systems, Z ν A nuc 1 0 11 1 s and ρ 1 0 22 1 c m 3 . We can therefore use the following approximation in many cases: N ˙ V 1 0 33 1 s c m 3 e - W r * kT (5.5)
A technically useful nucleation current, i.e. a nucleation current that leads to a phase transformation in a “reasonable” amount of time, is N ˙ V 1 1 s c m 3 . Using eq. 5.5, what is the maximum value that the reversible work of nucleation can take? N ˙ V 1 0 33 1 s c m 3 e - W r * kT 1 1 s c m 3 (5.6) e - W r * kT 1 0 -33 (5.7) W r * kT log 10 e33 (5.8) W r * kT 33 log 10 e (5.9) W r * 76kT (5.10)
This means that the barrier to nucleation, i.e. the reversible work, in a unary system, generally needs to be smaller or equal to 76kT (at the relevant temperature) in order to occur at a rate that is 'reasonable', i.e. of some technical use. We can use this value to quickly estimate whether we can expect to observe nucleation or not. Consider for example the condensation of water from its vapor for Δ T20 K will W r * <76kT. Please note that this quick-and-dirty approach is only valid for unary systems.
It is instructive to consider the nucleation current N ˙ V as a function of Δ T . We can attempt to do so analytically, by considering all temperature-dependent variables in eq. 5.5, or numerically, e.g. using realistic values for all materials properties.For the analytical approach, write N ˙ V in a way that makes the temperature dependence more apparent: N ˙ V 1 0 33 1 s c m 3 exp ( - W r * kT ) (5.11) N ˙ V exp ( - const Δ T 2 T ) (5.12)
Note that we have ignored small effects, such as the dependence of the molar volume (of a solid or liquid!) on temperature. Using the definition of Δ T= T e q -T , N ˙ V exp ( - const Δ T 2 ( T e q - Δ T ) ) (5.13)
We are interested in the interval from Δ T=0( T= T e q ) to T=0 K ( Δ T= T e q ). By sketching out the terms in the denominator of the exponent, we can quickly get an idea of the overall shape of the graph of N ˙ V (Figure 5.2).
image: 22_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_04_Fig2.svg
Figure 5.2: Sketching the graph of N ˙ V as a function of Δ T , Eq. 5.13. a) The denominator of the fraction in the exponent consists of two terms, T e q - Δ T (green) and Δ T 2 (blue). The denominator is the product of these terms (red). The product goes to zero at both side of the Δ T interval and has a maximum in between. Note that the each of the functions was sketched on a different y -axis. b) Next, we sketch the negative reciprocal of the product in (a). Note that the graph goes to negative infinity on either side of the Δ T interval and has a maximum at a negative value in between. c) Finally, we sketch Eq. 5.13, and find that N ˙ V goes to zero at either end of the interval, and goes through a maximum in between.
Using an appropriate graphing software, we can explore a specific example, for example the condensation of liquid water from vapor at ambient pressure (Figure 5.3).
image: 23_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_04_Fig3.svg
Figure 5.3: Nucleation of liquid water from its vapor. Nucleation current N ˙ V is plotted on the x -axis, supercooling Δ T on the left y -axis, and the absolute temperature T on the right y -axis. Note that in b) and d) the very strong dependence of N ˙ V on temperature was made more apparent by plotting on a log scale. a) and b) show the entire physically relevant temperature range. c) and d) show the range in which N ˙ V passes through 1 s -1 c m -3 . Note that N ˙ V 0 for Δ T0 ; this is because the reversible work of nucleation, W r * Δ T -2 , goes to infinity as T T e q . Also, N ˙ V 0 for T0 , because even a small barrier becomes insurmountable as the thermal energy of the system goes to zero. Finally, there is a pronounced maximum in N ˙ V , here at Δ T250 K. At this temperature ( T-15 0 C ), we would also have to consider direct nucleation of solid ice. This plot was created using Lecture_04_Fig3.m.

Heterogeneous Nucleation in Unary Systems

In this section, we will consider the impact of heterogeneity in the system. Specifically, we will assume the presence of a flat surface of some substrate σ . As before, we will look at the formation of a nucleus (phase: β * ) from a parent phase L. The major difference to the homogeneous case is that the shape of a nucleus formed at the interface between σ and L is generally not spherical. In addition, when determining the excess excess free energy of the nucleus, we will also have to consider interfacial free energies other than γ L β * . Let's first consider the shape of the nucleus.
image: 24_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_05_Fig1.svg
Figure 6.1: Force balance at the triple line phase boundary for a small (axisymmetric) volume of the phase β * at the interface between σ L. a) General case in which both horizontal and vertical forces are considered. b) Assuming an infinitely stiff substrate σ , vertical forces are ignored and β * takes the shape of a spherical cap with contact angle θ .
Recall from MSE316-1, that the equilibrium shape at a triple line phase boundary ( L - σ - β * ) results from the balance of the vector surface tension of the L - σ , σ - β * , and L - β * interfaces (Figure 1A). If we assume that σ is infinitely stiff, we can ignore the vertical force balance (Figure 1B). In the latter case, β * takes the shape of a spherical cap with radius R * and contact angle θ .
We find θ by considering the horizontal force balance only γ L σ = γ σ β * + γ L β * cos θ (6.1) cos θ = γ L σ - γ σ β * γ L β * (6.2)
image: 25_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_05_Fig2.svg
Figure 6.2: 3D rendering of a nucleus of the phase β * in the shape of a spherical cap, at the interface between substrate σ and the parent phase L.
To find the reversible work of heterogeneous nucleation, we consider the change in the grand canonical potential as we go from the initial to the final state W r *, h et = Ω f - Ω i (6.3) =- Δ P V β * + γ L β * A L β * + γ σ β * A σ β * - γ L σ A L σ , (6.4) where V β * is the volume of the nucleus, A L β * is the interfacial area between the nucleus and L, and A σ β * is the area of the interface between the substrate and the nucleus. While both these areas are created during nucleation, the area of the interface between the substrate and L, A L σ , is removed in the process. Therefore, it enters into Eq. (3) as a negative term. Note that A σ β * = A L σ (Figure 2).
By writing the volume and area terms in Eq. (3) as functions of θ , one can show that W r *, h et = W r * S( θ ) , (6.5) where W r * is the reversible work of nucleation for the homogeneous case and the structure factor S( θ ) is defined as S( θ ) = 1 4 ( 2+ cos θ ) ( 1- cos θ ) 2 (6.6) Inspection of Eq. (6) and its graph (Fig. 3A,B) reveals the following properties
image: 26_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_05_Fig3_alt.svg
Figure 6.3: Plot of the structure factor for the spherical cap, S( θ ) , vs. the contact angle θ using normal (a) and logarithmic (b) scaling of the y -axis.
Note that the critical radius R * is the same in both homogeneous and heterogeneous nucleation and can be calculated as discussed earlier. However, the volume and surface area of the spherical cap are generally smaller. R * , het = R * , hom = R * (6.7)
Using an approach similar to the one discussed in the homogeneous case, the heterogeneous nucleation current can be shown to be N ˙ V h et =Z ν A n uc ρ e - W r * kT S( θ ) (6.8)
However, shape of the nucleus and the substrate impact the terms of the prefactor. We can frequently estimate N ˙ V h et in unary condensed matter systems as N ˙ V h et 1 0 26 1 s c m 2 A t ot σ V t ot e - W r * kT S( θ ) , (6.9) where A t ot σ is the total area of the nucleating surface, and V t ot is the total volume of the system. In a simple case, this area could simply be the total area of the container walls, and the volume that of the container (and therefore, L before the transformation). However, nucleators could also be finely dispersed in the volume itself. Finally, in most real systems we would have to at least consider that there may be a distribution of nucleators of different potency (i.e. differing contact angle) on the container walls and dispersed in the bulk.
image: 27_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_05_Fig4.svg
Figure 6.4: Heterogeneous nucleation in a v-shaped scratch. The angle enclosed by the sidewalls of the scratch and the xy plane is given by α . Note that the nucleus here takes the shape of a segment of a cylinder with radius R * .
The structure factor also depends strongly on the assumptions we make regarding the shape of the substrate. For example, if we consider a nucleus forming in a v-shaped scratch in the surface (Fig. 4), we find that one can write the reversible work of nucleation as W r *, h et = W r * f( α , θ ) , (6.10) where α <9 0 is the angle between the side walls and the horizontal plane, and θ is the contact angle. While f( α , θ ) looks too intimidating to be shown here, its behavior is really very interesting: α 0 f( α , θ ) S( θ ) (6.11) α >0 f( α , θ ) <S( θ ) (6.12) θ - α +0 f( α , θ ) 0 (6.13) θ < α f( α , θ ) =0 (6.14)
This means that any scratch in a surface is a better nucleator than the smooth surface. Furthermore, if there are many scratches with a wide distribution of α in the substrate, and the contact angle θ is smaller than 9 0 , there will be at least a few for which the structure factor becomes vanishingly small. As a consequence, the reversible work of nucleation will become very small, and nucleation will occur quickly at these sites. As only one nucleus is in principle sufficient to initiate phase transformation (if growth is reasonably fast, which is not necessarily the case in solid state transformations), the most powerful nucleators can play an important role. It is therefore generally true that rough or porous substrates are more efficient nucleators.
However, if scratches or pores set in a smooth surface get very small, nuclei may not be able to grow out of them and into the bulk. Why? Hint: Consider the radius of curvature of the newly formed phase as it grows out of the scratch and compare it to the critical radius.

Nucleation in Binary Systems

In the next couple of sections, we will take a look at nucleation in binary systems. Consider a typical precipitation in the solid state (Figure 1A). In this eutectic system with components 1 and 2, and the terminal solutions α (rich in component 2) and β (rich in component 1), we start with α at some initial mole fraction X 1 i and cool through the solvus line into the two-phase field, then hold at some temperature T . We are therefore considering the phase transformation α α + β . For now, let's pretend the system is homogeneous. From the phase diagram, we find the undercooling Δ T= T e q -T (here, the relevant equilibrium temperature is given by the solvus line at X 1 i ). Drawing the tie-line at T , we find the equilibrium mole fraction of component 1 in α , X 1 ,eq α , and the equilibrium mole fraction of component 1 in β , X 1 ,eq β . The supersaturation of undercooled α at T is Δ X= X 1 i - X 1 ,eq α .
Next, let's first look at the driving force for phase transformation. For this, consider the GX diagram corresponding to temperature T that underlies the TX (phase) diagram (Figure 1B). The free energy of the α phase as a function of the mole fraction of component 1 at temperature T and pressure P α is given by G m α ( T, P α , X 1 ) , that of β by G m β ( T, P α , X 1 ) . The initial state of the system is that of a homogenous, but supersaturated/undercooled α with G m α ( T, P α , X 1 i ) (point I). To reduce the free energy, β precipitates enriched in component 1 form from the matrix α , depleting the matrix of component 1 in the process. Equilibrium is reached where X 1 = X 1 i intersects the common tangent (point F). In this state, the mole fraction of component 1 in α is that of the point of tangency, X 1 ,eq α , and that in β is X 1 ,eq β . The free energy of the system in this state can be calculated by calculating the weighted sum X α G m α ( T, P α , X 1 ,eq α ) + X β G m β ( T, P α , X 1 ,eq β ) , using the lever rule to determine X α and X β . The change in free energy associated with the phase transformation, Δ G m α α + β , is identical to the negative distance FI ¯ . Note, however, that we consider here both phases at the same pressure P α . This means that the Laplace pressure and the curvature of the interface is zero. This corresponds to complete phase separation as shown in Figure 1C.
image: 28_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_06_07_Fig1.svg
Figure 7.1: Phase transformation in a binary system. a) Generic binary eutectic phase ( TX ) diagram. b) GX diagram for the system at temperature T in (a). The free energy change Δ G m α α + β corresponds to the energy per mole of α converted, i.e. the overall driving force. c) Schematic drawing that corresponds to the initial and final equilibrium states referenced in (a) and (b). Note that the α - β interface has no curvature, therefore Δ P=0 !

Instead of looking at the free energy change for the overall phase transformation, we need to consider the free energy change associated with the formation of the nucleus from the parent phase (Figure 2). Let's assume the nucleus is spherical with radius R * and incompressible.
image: 29_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_06_07_Fig2.svg
Figure 7.2: Schematic drawing showing the initial state of the homogeneous α at ( T, P α , X 1 i ) , and the unstable equilibrium of one very small, spherical nucleus β * with radius R * at ( T, P β * , X 1 ) with the surrounding matrix α . Note that the transformation is isothermal and the nucleus so small that its formation does not change the composition of α .
Recall that if we assume incompressibility, G m P = V m (7.1) Δ G m Δ P = V m (7.2) G m ( T, P f , X 1 ) - G m ( T, P i , X 1 ) ( P f - P i ) = V m (7.3) G m ( T, P f , X 1 ) = G m ( T, P i , X 1 ) + V m Δ P (7.4)
This means in the GX diagram, increasing the pressure from an initial value P i to a final value P f for any given phase simply shifts the free energy to higher energy by the amount V m Δ P , where Δ P= P f - P i .
For the nucleus, this means that G m β ( T, P β * , X 1 ) = G m β ( T, P α , X 1 ) + V m β Δ P (7.5)
V m β Δ P0 is therefore the contribution of the interface to the molar free energy of the nucleus. In unstable equilibrium, this must be exactly balanced by the free energy released when a solid of the composition of the nucleus is formed from the matrix. One way to think about this scenario is to consider the corresponding GX diagram (Fig. 3).
The major difference between Fig. 3 and Fig. 1B is that we now also consider the unstable equilibrium of the matrix at composition X 1 i with the nucleus. We know that there must be a common tangent between the free energy curves for the matrix and the nucleus. One of the points of tangency is G m α ( P α , X 1 i ) and the slope of the tangent is equal to the slope of G m α at X 1 i . A priori, we do not know what value P β * takes, but really all we need to do is shift G m β ( P α , X 1 ) vertically until there is a point of tangency to the tangent in G m α . The vertical shift is then exactly equal to V m β Δ P and the point of tangency indicates the composition of the nucleus, X 1 .
image: 30_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_06_07_Fig3.svg
Figure 7.3: Schematic GX diagram for a binary system at temperature T similar to the one in Fig. 1A. The molar free energy for α is given for P α (purple line). For β , molar free energies at P α (green dashed line) and P β * (green line) are given; the latter is shifted to higher free energy by V m β Δ P . The first of two common tangents is for the final equilibrium between matrix α at composition X 1 ,eq α and pressure P α , and precipitate β at composition X 1 ,eq β and pressure P α (cyan dashed line), similar to the common tangent in Figure 1B. The second is for the unstable equilibrium between the matrix α with composition X 1 i and pressure P α , and nucleus with composition X 1 and pressure P β * (red dashed line).
By inspection, we find that G m β ( P α , X 1 ) + V m Δ P- G m α ( P α , X 1 i ) = . G m α X 1 | X 1 i ( X 1 - X 1 i ) (7.6) - V m Δ P= G m β ( P α , X 1 ) - G m α ( P α , X 1 i ) - . G m α X 1 | X 1 i ( X 1 - X 1 i ) (7.7)
Now let's see whether we can show that the r.h.s. of Eq. (7.7) is indeed the molar free energy change for the formation of a nucleus, i.e. the Gibbs free energy change associated with the formation of an infinitesimal amount of material with a composition X 1 that is different from the initial composition X 1 i . This amount shall be so small that the mole fractions of the componenta 1 and 2 in the matrix shall not change. We can then calculate the Gibbs free energy change by
Note that the free energy required to remove material is not equal to the free energy gained because the system is closed and the total number of atoms is fixed.
Step 1. The free energies of the initial and final states of the matrix α are then G α , i = G α ( n 1 i , n 2 i ) (7.8) G α , f = G α ( n 1 f , n 2 f ) (7.9) , where the number of atoms of component j in the initial state is given by n j i and that in the final state given by n j f .
For a very small number of atoms n j = n j f - n j i removed, we can expand G α , f around G α , i to first order only. At constant T and P α ,
Using the definition of the chemical potential μ j = G n j , the free energy change can be written as Δ G r = G α , f - G α , i (7.10) = μ 1 α ( X 1 i ) ( n 1 f - n 1 i ) + μ 2 α ( X 1 i ) ( n 2 f - n 2 i ) (7.11) =- μ 1 α ( X 1 i ) n 1 - μ 2 α ( X 1 i ) n 2 (7.12)
To find the molar free energy change associated with the removal from the matrix, we divide by the total number of moles of removed, n tot = n 1 + n 2 . Δ G m r =- μ 1 α ( X 1 i ) X 1 - μ 2 α ( X 1 i ) X 2 (7.13)
Step 2. The free energy of formation for an infinitesimal amount of material of composition X 1 is simply its molar free energy. G m β ( X 1 ) = μ 1 β ( X 1 ) X 1 + μ 2 β ( X 1 ) X 2 (7.14)
Step 3. The total free energy change for removal and formation of material of composition X 1 from a matrix of composition X 1 i is then the sum of Eq. (7.13) and Eq. (7.14). Δ G m α β * = Δ G m r + G m β ( X 1 ) (7.15) =[ μ 1 β ( X 1 ) - μ 1 α ( X 1 i ) ] X 1 +[ μ 2 β ( X 1 ) - μ 2 α ( X 1 i ) ] X 2 (7.16) This is a very general result that can be used even if the difference in mole fractions is large, as long as only a small amount of β * is made.
But for our purposes, we just substitute Eq. (7.13) in Eq. (7.16) Δ G m α β * = G m β ( X 1 ) - μ 1 α ( X 1 i ) X 1 - μ 2 α ( X 1 i ) X 2 (7.17)
Separately, we express G m α ( X 1 i ) in terms of the chemical potentials, rearrange to 0=- G m α ( X 1 i ) + μ 1 α ( X 1 i ) X 1 i + μ 2 α ( X 1 i ) X 2 i (7.18) , and then add the l.h.s. of Eq. (7.17) to the l.h.s of Eq. (7.18) and the r.h.s. of Eq. (7.17) to the r.h.s of Eq. (7.18)
Δ G m α β * = G m β ( X 1 ) - G m α ( X 1 i ) + μ 1 α ( X 1 i ) ( X 1 i - X 1 ) + μ 2 α ( X 1 i ) ( X 2 i - X 2 ) (7.19)
Substituting X 2 =1- X 1 Δ G m α β * = G m β ( X 1 ) - G m α ( X 1 i ) +[ μ 1 α ( X 1 i ) - μ 2 α ( X 1 i ) ] ( X 1 i - X 1 ) (7.20)
Finally, with we can express the change in free energy in terms of the Gibbs free energies
As we fixed the pressure at P α , this is identical to
Finally, comparison of Eq. (0.1) and Eq. (7.7) reveals that the right hand sides are identical. This means that the left hand side must also be identical
One can further show that for small Δ X= X 1 i - X 1 ,eq α
Using the definition of the supersaturation and substituting into Eq. (0.1) , where C is (nearly) constant
With this, we can calculate the critical radius for the spherical, incompressible nucleus, and the reversible work of nucleation for the binary case. R * = 2 V m C γ Δ X (7.21) W r * = 16 π 3 ( V m C ) 2 γ 3 Δ X 2 (7.22)
These expressions resemble those we found for unary systems. Note how the supersaturation Δ X replaces the undercooling Δ T , and C takes the places of Δ S m . Because Δ T and Δ X are not independent variables, they can be used interchangeably when discussing effects qualitatively.

Strain Effects in Nucleation

Some of the most interesting phase transformations in materials science occur in the solid state. In this case, we need to consider the effect of strain at the interface between the nucleus and the matrix. Recall that strain is unique to solid-solid, coherent or semi-coherent interfaces. Strain is caused by lattice mismatch across the interface and results in an elastic deformation of bonds. We can define an elastic strain energy density, W V e l , that is equal to the total strain energy divided by the volume of the nucleus W V e l = W t ot e l V β [ W V e l ] = J m 3 = Pa (8.1)
We can therefore write down the grand canonical potential of the system Ω =- P β V β - P α V α + γ α β A α β + W V e l V β (8.2)
In unstable equilibrium, for a spherical nucleus with critical radius R *
Using
, and collecting terms
For the spherical nucleus, we substitute volume and area, differentiate, and rearrange R * = 2 γ α β Δ P- W V e l (8.3) It is equally straightforward to show that W r * = 16 π 3 γ α β 3 ( Δ P- W V e l ) 2 (8.4)
Therefore, both R * and W r * increase with increasing W V e l . In other words, strain makes it more difficult for nucleation to occur. Note that for W V e l Δ P , R * and W r * , similar to what happens in strain-free systems when Δ T (unary case) or Δ X (binary case) approach zero, i.e. when the state of the system approaches the solvus line. Strain can therefore be thought as shifting the solvus to lower temperature. For the unary case We introduce the effective undercooling/supersaturation u nary Δ P- W V e l = C u Δ T e ff (8.5) C u Δ T- W V e l = C u Δ T e ff (8.6) Δ T e ff = Δ T- W V e l C u (8.7) b inary Δ P- W V e l = C b Δ X e ff (8.8) Δ X e ff = Δ X- W V e l C b (8.9) This allows us to define define the shift of the equilibrium temperature in the unary case, and of the solvus line in the binary case (Figure 1). u nary Δ Δ T= W V e l C u b inary Δ Δ X= W V e l C b (8.10)
image: 31_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_08_09_Fig1.svg
Figure 8.1: Effect of strain on the solvus in a binary system. a) Schematic GX diagram. Elastic strain energy shifts G m β to higher energy. Note that this shift is addition to the shift Δ P V m that results from any curvature of the interface. As a consequence, a new equilibrium is established and the equilibrium compositions shift to new values, X 1 ,eq α ,str and X 1 ,eq β ,str . b) In the phase diagram, the shift of the equilibrium in the presence of strain is equivalent to a shift of the solvus by Δ Δ X= X 1 ,eq α ,str - X 1 ,eq α . As a consequence the effective supersaturation is reduced.
One could conclude that the reversible work of nucleation for coherent, strained nuclei should always be larger than that for incoherent nuclei, where there is no strain. However, the interfacial free energy of the coherent interface, γ C , can be considerably smaller than that of the totally incoherent interface, γ i . We therefore have to compare W r * , c = 16 π 3 γ C 3 ( C b Δ X- W V e l ) 2 (8.11) W r * , i = 16 π 3 γ i 3 ( C b Δ X ) 2 (8.12)
Comparing the graphs of Eq. (8.12) and Eq. (8.12), it is apparent that the asymptote of W r * , c is shifted by Δ Δ X= W V e l C b towards greater Δ X (Figure 2). For small Δ X , therefore W r * , i < W r * , c , and we expect incoherent nuclei to form faster. However, as Δ X is small, the nucleation current will be small. For large Δ X , we can neglect the contribution of the strain energy density to the denominator in Eq. (8.12). With γ i > γ C , we find that the barrier W r * , i > W r * , c ; therefore coherent nuclei form more readily.
image: 32_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_08_09_Fig2.svg
Figure 8.2: Plot of the reversible work of nucleation for incoherent (blue line) and coherent nuclei (red line), using γ i =0.5 Jm -2 ; γ C =0.2 Jm -2 ; W V e l =2 GPa ; C b =50 GPa .
Q: Derive an expression for the supersaturation at which there is a crossover between W r * , c and W r * , i . Discuss the dependence of Δ X xo on the ratio of the interfacial free energies and other relevant properties.
Q: Replot Figure 2 with the W r * expressed in units of kT . Assume that the solvus temperature at Δ X=0 is 473 K, and that the solvus is linear with a slope of 100 K. Discuss the expected nucleation rate for coherent and incoherent nuclei as a function of Δ X .
If we assume that the nucleus is not only spherical and incompressible, and that the capillarity assumption holds, but also that the interface between nucleus ( β ) and matrix ( α ) is coherent and that both are elastically isotropic, the elastic strain energy density has a simple form W V e l =18 ε 2 μ α K β 3 K β +4 μ α , (8.13) where K β is the bulk modulus of the nucleus, μ α is the shear modulus of the matrix, and ε is the misfit parameter.
We can then consider two limiting cases,
Finally, if we assume that matrix and nucleus both have cubic structure ε = 1 3 V m β - V m α V m α (8.16) ε a β - a α a α , (8.17) where V m is the molar volume and a the lattice parameter.

Kinetics of Nucleation in the Solid State

The nucleation current for homogeneous nucleation in unary systems can be written in the form of an Arrhenius law with several prefactors that, to first order, we can approximate by a constant. N ˙ V =Z ν A n uc ρ e - W r * kT 1 0 33 1 s c m 3 e - W r * kT (9.1) Recall that ρ is the number density of monomer units and the product Z ν A n uc can be interpreted as the frequency at which monomer units attach successfully to the nucleus. In unary systems, we can get away with ignoring the temperature dependence of this frequency. However, in binary and multicomponent systems, where thermally activated diffusion becomes much more important, we write instead N ˙ V = ω e - Δ G d kT ρ e - W r * kT = ω ρ e - Δ G d + W r * kT , (9.2) where ω is again an attachment frequency with units of [s -1 ], and Δ G d is the free energy of activation for the relevant diffusive processes. We can directly see that all other things remaining equal, the nucleation rate decreases if Δ G d increases. We therefore expect nucleation to be more rapid where diffusion is “easy”.
Previously, we found that the nucleation rate in unary systems increases with increasing supercooling, because the reversible work of nucleation decreases ( W r * Δ T -2 ), but eventually decreases again at very low temperatures because it becomes impossible for the system to overcome even a small (but finite) barrier (section 4, Figures 2 and 3). In binary systems, the reversible work of nucleation instead is proportional to Δ X -2 (section XYZ). Δ T is of course a function of Δ X , and for the sake of the following discussion let's pretend this function is linear. This is equivalent of assuming that the solvus is a straight line. If so, then we only need to consider the impact of the additional exponential term in binary systems (Figure 1). The term e - Δ G d kT goes to unity at high temperatures and decreases rapidly as the temperature decreases. Generally speaking, this reduces the nucleation rate across the board, but especially so at low temperatures. It also shifts the maximum rate to slightly higher temperature, i.e. lower supercooling.
image: 33_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_10_11_Fig1.svg
Figure 9.1: Dependence of the nucleation current N ˙ V (black line) on temperature T . N ˙ V is proportional to the product of two exponential terms. One is dependent on the reversible work of nucleation (blue dashed line), which itself is dependent on temperature, and approaches zero for temperatures approaching the equilibrium temperature (where the barrier is infinite) and absolute zero (where the thermal energy is insufficient to overcome even a small, but finite barrier). This term goes through a maximum between at an intermediate temperature. The other term represents diffusive processes (red dashed line), which increase with increasing temperature but are strongly attenuated at low temperature. The product of the two terms, and therefore the nucleation current, has the same general shape as the first term. However, the maximum is shifted towards higher T .
It is often useful to determine in which of two similar systems nucleation would occur at a higher rate. Consider for example the precipitation α α + β from a binary mixture with slightly different initial compositions, X 0 1 and X 0 2 (Figure 2A). The solvus line indicates the equilibrium temperature, below which nucleation of β may occur. From X 0 1 > X 0 2 follows that T e q 1 > T e q 2 . Let's first look at the case where we compare rates at the same absolute temperature. Clearly, Δ T 1 > Δ T 2 . As a consequence, the reversible work of nucleation is lower, and the nucleation current is therefore greater for the mixture with the higher initial mole fraction for most of the temperature range that is relevant for processing (Figure 2B). At very low temperatures, where the nucleation rate is attenuated by 'freezing out' of the thermally activated processes, both rates become nearly identical.
image: 34_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_10_11_Fig2.svg
Figure 9.2: A. hypothetical binary phase diagram for precipitation in the solid state. B. Comparison of the homogeneous nucleation current of two binary mixtures with initial compositions X 0 1 and X 0 2 as a function of the temperature. C. Comparison of the homogeneous nucleation current in the same system as (B), but as a function of the supercooling Δ T .
If we, however, consider the rates at identical undercooling Δ T , the reversible work of nucleation is the same for both compositions. The barrier to diffusion being independent of temperature, one might conclude that the rates should be identical. However, the absolute temperature does go into the rate in form of the factor of (kT ) -1 in the exponent. In the example shown in Figure 2C, this results in a higher nucleation rate for the mixture with the lower initial mole fraction at low supercooling, and a crossover near the maximum.
Clearly, the assumption that the solvus is a linear function breaks down even for small supercooling. In many terminal solid solutions, the solvus approaches the the T -axis. How does this affect the argument made above?
Now consider a system at a fixed initial composition. How would strain impact the graph of the nucleation rate vs. temperature?
Recall that we previously found that heterogeneous nucleation is frequently much faster than homogeneous nucleation. However, it is important to remember that both can and will occur in any given system. In complex systems, such as solid materials, there are many internal features that can serve as heterogeneous nucleators. Before we can address kinetics, we should develop at least a qualitative picture of these. In principle, any inhomogeneity that lowers the free energy of a nucleus can act as a nucleator.
Roughly in order of increasing free energy of formation of such defects, there are
To better understand how defects can act as nucleators, consider the lattice strain near a vacancy, substitutional defect, or an edge dislocation (Figure 3).
image: 35_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_10_11_Fig3.svg
Figure 9.3: Lattice strain in the vicinity of defects. A. Vacancies result in tensile (blue) strains. B. Depending on the radius, substitutional defects can result in either compressive (here: red) or tensive strain. C. Edge dislocations are associated with both tensile (blue) and compressive (red) strain. A nucleus with a molar volume that is higher than that of the matrix would preferentially form in volumes with tensile strain, whereas a nucleus with a smaller molar volume would find more favorable conditions in a volume with compressive strain.
In case of the vacancy, the absence of an atom creates a region of tensile strain in the lattice. For a substitutional defect with an atomic radius that is larger than the matrix, a compressive strain results. In an edge dislocation, the insertion of a half plane results in regions with compressive strain and tensile strain. If we assume that a phase β nucleates from α , and that the molar volume of the precipitate is greater than that of the matrix, i.e. V m β > V m α , then we would expect that the free energy required to form a nucleus is lower in a volume where there is tensile strain. Formation of the nucleus would not only reduce the existing tensile strain, but also result in a lower final compressive strain. We can express this in very general terms as W r * =( - Δ P+ W V e l ) V m β * + γ A- Δ G d efect , (9.3) where Δ G d efect represents the free energy gained from the interaction of the defect and the nucleus. It can be rather tricky to accurately quantify this for any given system.
We have previously treated nucleation at free surfaces, the simplest case of heterogeneous nucleation on an infinitely stiff, flat surface, and found that the nucleus takes the shape of a spherical cap (section 5; Figure 4A). A conceptually similar approach can be used to model nucleation at grain boundaries, i.e. interfaces between single crystalline domains of the same phase, and interphase boundaries, i.e. interfaces between single crystalline domains of different phases.
image: 36_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_10_11_Fig4.svg
Figure 9.4: Rendering of three-dimensional models of nuclei formed at a free surface (A) and at grain boundaries (B-C). A. At a free surface (gray plane), the nucleus takes the shape of a spherical cap with contact angle θ and radius R * . B. At an idealized grain boundary (gray plane) represented by one value for the surface free energy, the nucleus has a lentil shape that consists of two spherical caps, again with contact angle θ and radius R * . C. At an idealized grain edge, i.e. the line of intersection between three grain boundaries, the nucleus can be thought of as the volume of intersection of three spheres of radius R * . D. At an idealized grain corner, i.e. the point of intersection between four grain edges, the nucleus takes the shape of volume of intersection of four spheres of radius R * .
Let's consider the simplest case, a planar interface between two grains of matrix α (Figure 4B). For simplicity, let's further assume that the interface between a nucleus of β and the matrix is totally incoherent with a surface free energy γ α β , and independent of the nature of the α - α interface (with surface free energy γ α α ) that has five degrees of freedom. Under these conditions we expect that the β nucleus will take a lentil shape that consists of two equal halves that are each a spherical cap. Considering the force balance, we find for the contact angle
cos θ = γ α α 2 γ α β (9.4)
We can write the reversible work of nucleation by considering the volume of the nucleus and the interfacial areas generated and destroyed:
W r ,gb * =- Δ P V m β +2 γ α β A α β - γ α α A α α (9.5)
One can show that eq. 9.5 can be written as the product of the reversible work of homogeneous nucleation and a structure factor that depends on θ .
W r ,gb * = W r ,hom * S g b ( θ ) (9.6)
A similar argument can be made for nucleation at two-dimensional grain edges (Figure 4C), where three planar grain boundaries meet (all dihedral angles shall be 12 0 ), or at tetrahedral grain corners (Figure 4D) . We can now compare the reversible work at different (idealized) grain boundaries (Figure 4B) and at free surfaces (Figure 4A) that we discussed previously:
W r ,het * = W r ,hom * S( θ ) (9.7) W r ,gb * = W r ,hom * S g b ( θ ) (9.8) W r ,ge * = W r ,hom * S g e ( θ ) (9.9) W r ,gc * = W r ,hom * S g c ( θ ) (9.10)
, where S( θ ) = 1 4 ( cos θ +2 ) (1- cos θ ) 2 . Without going into detail of the derivation and final form of the structure factors at grain boundaries, edges, and corners, a comparison of the graphs is very informative (Figure 5). Note that the contact angle for the free surface depends on three different surface energies, whereas the for all other cases, it only depends on γ α α and γ α β * . It is therefore not meaningful to plot the structure factor S( θ ) on the same axes as the other three.
image: 37_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_10_11_Fig5.svg
Figure 9.5: Graph of the structure factors for an idealized grain boundary (red), grain edge (green), and grain corner (blue) plotted on a linear (A) and logarithmic (B) vertical axis against the contact angle on the horizontal axis. Note that the structure factors all converge on 1 for θ 9 0 . At this contact angle, the nucleus is spherical and thus identical to the nucleus formed by homogeneous nucleation. For all angles smaller than 9 0 , the reversible work of nucleation at a grain corner is smaller than that at grain boundaries, and that at grain corners is smaller still. Note that structure factors for grain edges and corners have vertical asymptotes, meaning the go to zero for an angle greater than 0 . This can be seen better in (B). Why is that and what are the exact angles?
Imagine an idealized polycrystal of α , where all grains are identical in size, and all grain boundaries, edges, and corners can be described using the assumptions we made above. What is the grain boundary area, grain edge length, and number of grain corners (per unit volume) as a function of the grain size? Where is nucleation β going to occur first? Is nucleation only going to occur there? Briefly describe the sequence of nucleation events that you would expect in such a material.

Solidification in Unary Systems

In this section, we will look at solidification of pure phases, for instance the solidification of a metal or polymer melt, or the freezing of water. We are particularly interested in what happens at the interface between the solid and the liquid. We differentiate two general cases, diffuse interfaces and (atomically) flat interfaces (Figure 1).
image: 38_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_12_Fig1.svg
Figure 10.1: Schematic drawing of idealized structure of atomic solid-liquid interfaces. A. Diffuse, atomically rough interface of crystalline solid (black) and liquid (blue). B. Atomically flat interface.
Diffuse interfaces result when atoms or molecules in the liquid that attach to the surface of the solid have a high probability of sticking to it. This results in an atomically rough surface. Movement of the interface is then controlled by how quickly atoms are transported to the interface by diffusion. A useful rule of thumb is that if the entropy of fusion is small, Δ S f 10 J K , which is true for many metals, diffuse interfaces form, and the process of solidification is said to be diffusion-controlled. Atomically flat interfaces form when the entropy of fusion is larger. In this case, atoms will stick only rarely, and if they stick will migrate to ledges and kink sites to incorporate into the lattice. This process results in rather flat surfaces and is said to be interface controlled. For the remainder of the chapter, we will only consider diffusion-controlled processes.
During solidification, the solid-liquid interface moves with a velocity v ( [ v ] = m s ). As some liquid is converted into solid, heat is released. The amount of heat released per unit volume is the latent heat of solidification, L V ( [ L V ] = J m 3 ) , which is related to the standard enthalpy of fusion Δ H f .
In a unidirectional solidification, the amount of heat released per unit area of interface and unit time is q i =v L V , ( [ q i ] = Jm m 3 s = W m 2 ). This heat will increase the temperature of the interface unless it is transported away. Heat transport depends on the thermal conductivity k ( [ k ] = W mK ) and the temperature gradient dT dx :
q=-k dT dx (10.1)
For unidirectional solidification in the positive x-axis direction, we can then write k s d T s dx = k L d T L dx +v L V (10.2) , where the subscript s indicates heat transport in the solid and the subscript L indicates heat transport in the liquid.
image: 39_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_12_Fig2.svg
Figure 10.2: Unidirectional solidification of a cold solid into a hot liquid. A. One dimensional temperature profile across the interface. B. Two dimensional contour map of the same scenario as in (A). Note that contour lines correspond to isotherms. C. A small fluctuation in growth speed results in a protrusion (bump) in the interface. Close to the bump, isotherm shape is affected, but further away isotherms remain the same. D. Comparison of 1D temperature profiles across the flat part of the interface (red/black) and across the bump (green). Note decrease in interface velocity at the bump compared to flat interface.
Consider the scenario in Figure 2A. Here, the solid is cold and the liquid is hot. There is a steep temperature gradient in the solid, and heat is moved away from the interface. The temperature gradient in the liquid is shallower, but heat is transported towards the interface. Under these conditions, we can rewrite Eq. (2) to find the interface velocity
v= k s d T s dx - k L d T L dx L V (10.3)
If the heat conductivities in the solid and the liquid are similar ( k s k L ), then we can see that from d T s dx > d T L dx follows that the interface velocity is positive, meaning that the solid grows at the expense of the liquid.
Now let's consider the growth front in two dimensions (Figure 2B). An important question is whether the interface will keep its flat (but not atomically flat) shape during the solidification. A priori, the amount of heat removed through the solid and the amount of heat delivered through the liquid are the same all along the interface. We therefore expect that the velocity is constant everywhere. However, let's assume that by some fluctuation, a small protrusion or bump forms on the interface (Figure 2C). This will result in a small change of the isotherm shape in the vicinity of the bump, but we expect that the temperature contours far away from the bump will be identical to the case where there is no bump. Inspecting Figure 2C or the one-dimension temperature profile in Figure 2D, we can see that the temperature gradient in the liquid ahead of the bump becomes steeper compared to the flat interface. At the same time, the temperature gradient in solid becomes shallower. This means that more heat arrives and less heat is removed at the interface. As a consequence, the numerator in Eq. (3) and therefore the velocity of the interface of the bump decreases. This means that the flat interface will catch up with the bump, and the flat interface shape will be restored. We call this interface stable.

Now consider the case where the solid is hot and the liquid is supercooled (Figure 3). This scenario describes a small piece of solid phase suspended in the liquid phase. There is no temperature gradient in the solid, and heat is transported away from the interface into the supercooled liquid (Figure 3A). As before, let's assume the growth velocity in the x -axis direction is uniform and add another spatial dimension (Figure 3B). If by some fluctuation, a small bump forms where the interface juts into the liquid (Figure 3C), we now find the temperature gradient ahead of the bump is steeper than ahead of the flat interface (Figure 3C&D). As a consequence, the numerator in Eq. (3) and therefore the interface velocity of the bump increases, further increasing the steepness of the gradient, resulting in further increase of the growth speed. As a consequence, the interface is not stable against small perturbations.
image: 40_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_12_Fig3.svg
Figure 10.3: Unidirectional solidification of a hot solid into a supercooled liquid. A. One dimensional temperature profile across the interface. B. Two dimensional contour map of the same scenario as in (A). Note that contour lines correspond to isotherms. C. A small fluctuation in growth speed results in a protrusion (bump) in the interface. Close to the bump, isotherm shape is affected, but further away isotherms remain the same. D. Comparison of 1D temperature profiles across the flat part of the interface (red/black) and across the bump (green). Note that the interface velocity at the bump increases against the velocity of the flat interface.
A common consequence of growth conditions that result in instable interfaces is the development of a dendritic microstructure (Figure 4). Accelerated growth of the initial protrusion results in a finger-like protrusion. The side-walls of this primary arm are also affected by interface instability and secondary arms can form. This process remains poorly understood. Nevertheless, secondary arm spacing in some systems provides information on the thermal history of a material. Secondary arms can branch again, resulting in tertiary arms. Of course, branching can and will occur in all three spatial dimensions. Very generally speaking, primary dendrites are oriented towards the steepest temperature gradient. However, as dendrites are single crystals, orientation-dependent changes in growth speed affect the orientation of primary, secondary, and tertiary arms.
image: 41_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_12_Fig4.svg
Figure 10.4: Development of dendritic microstructure.

Solidification in Binary Systems

In this section, we will extend our analysis of solidification to binary systems. The concepts discussed are general, and apply to solidification of molten alloys, freezing of biological specimens, of the formation of microstructures from polymer melts, as long as growth is diffusion controlled, i.e. occurs with a diffuse interface.
Before we begin, let's consider an idealized binary phase diagram (Figure 1) in which the coexistence lines, and in particular liquidus T L ( X ) and solidus T s ( X ) lines are all linear functions.
image: 42_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_13_14_Fig1.svg
Figure 11.1: An idealized binary phase diagram where coexistence lines are linear functions of the composition X . Note that we drop the index for the component for convenience.
We can then write T L ( X ) = T m α + b L X (11.1) T s ( X ) = T m α + b s X , (11.2) where b L = d T L dx is the slope of the liquidus and b s = d T s dx is the slope of the solidus.
For any given isotherm (tie line) in the α +L 2-phase field, the intercept of the tie line with the solidus shall be at X s and the intercept with the liquidus at X L .
Therefore, T= T L ( X L ) = T s ( X s ) (11.3) T m α + b L X L = T m α + b s X s , (11.4) b L X L =+ b s X s (11.5) X s X L = b L b s k (11.6) where k is a constant that will become useful.
Consider now the solidification of L with initial composition X 0 (Figure 2A). Solidification will begin at T L ( X 0 ) and the first formed solid will have composition k X 0 . The last liquid to solidify will have composition X 0 k and will solidify at T L ( X 0 k ) . At some temperature T , where T L ( X 0 ) >T> T L ( X 0 k ) , the equilibrium composition of the solid is X s ( T ) and that of the liquid is X L ( T ) , where X s ( T ) X L ( T ) =k .
image: 43_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_13_14_Fig2.svg
Figure 11.2: Solidification from the melt. A. Idealized phase diagram in which liquidus and solidus are linear functions of composition. Solidification begins at T L ( X 0 ) and ends at T L ( X 0 k -1 ) . At any intermediate temperature T , the equilibrium composition of the solid is X s ( T ) and that of the liquid is X L ( T ) . B. Concentration profile across the interface in the direction of solidification, at some intermediate temperature T .
If we assume that the system is at equilibrium at any time, we can draw the concentration profile we expect to see for a system that solidifies unidirectionally (Figure 2B). For unidirectional solidification we can define the degree of solidification, f s , as the volume fraction of solid that has formed. This is equivalent to the x -position of the interface, x i , divided by the total length in the direction of solidification, x t ot . f s = V s V s + V L = x i x t ot (11.7) (11.8) If we assume that the molar volume of the liquid and the solid are identical, it is straightforward to determine f s for a given temperature using the lever rule. At some temperature T between T L ( X 0 ) and T L ( X 0 k ) , the composition of the solid will be X s ( T ) and that of the liquid X L ( T ) .
Let's look at an example, where T m α =1000 K , b L =-400 K , b s =-800 K , and thus k= 1 2 . With X 0 =0.15 , we find the temperature at which the first solid forms to be T L ( 0.15 ) =940 K , and the temperature which the entire system is solid T L ( 20.15 ) =880 K . The composition of the first-formed solid is X s (940 K )=0.075 and that of the last bit of liquid is X L (880 K )=0.3 (Figure 3A). Using the lever rule, we can determine f s for several temperatures between the onset and end of solidification (Figure 3B). Using equations (1) and (2), or by looking at the phase diagram, we can determine the composition of the solid and the liquid at these intermediate temperatures (Figure 3C).
image: 44_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_13_14_Fig3.svg
Figure 11.3: An example of a binary solidification with perfect mixing in the solid and liquid phase. A. Idealized phase diagram where k= 1 2 . B. Plot of the temperature of the interface as a function of the fractional degree of solidification ( f s ), calculated using the lever rule. C. Plot of the equilibrium composition of the solid and the liquid as a function of the fractional degree of solidification.
This then allows us to draw concentration profiles in the direction of the solidification (Figure 4). It is a good exercise to find an analytical solution for the dependence of f s on T and of the composition of the liquid and solid phases on f s . This will allow you to reproduce the plots in Figures 3 and 4 using Matlab.
image: 45_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_13_14_Fig4.svg
Figure 11.4: Plots of the concentration profiles in the direction of solidification for the system in Figure 3, at the very beginning and end of the solidification process, and for several intermediate stages. Note that T and f s are interdependent as shown in Figure 3B.

The assumption that the system goes to equilibrium at each temperature is unrealistic. If a solidification occurs at a reasonable rate, diffusion in the solid is not likely to have a major impact. Therefore, let's assume there is no diffusion in the solid, but perfect mixing, by convection and diffusion, in the liquid. The Scheil equations, also known as the non-equilibrium lever rule, then describe the composition of the solid and the liquid as a function of the fractional degree of solidification. X s =k X 0 (1- f s ) k-1 (11.9) X L = X 0 (1- f s ) k-1 (11.10)
For a derivation of the Scheil equations, see P&E. The composition of the solid and the liquid phases for the system in Figure 3A are shown in Figure 5. Note that while the composition of the liquid changes with time, the composition of the solid changes in space. The first formed solid, which precipitates at the highest temperature, has the lowest mole fraction of the minority component. As the phase transformation progresses, and the temperature drops, the rejected component accumulates in the liquid phase, and its concentration in the solid slowly increases. Unlike the perfectly mixed sample discussed previously, the mole fraction of the minority component in the liquid can go to very high values. Indeed, the model predicts that for k<1 , X L as f s 1 . This is clearly unphysical. Going back to the phase diagram in Figure 1, what would you expect to happen instead?
image: 46_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_13_14_Fig5.svg
Figure 11.5: Unidirectional solidification with no mixing in the solid and perfect mixing in the liquid is described by the Scheil equations. Here, the composition of the solid and the liquid as a function of the fractional degree of solidification are plotted for the system described by Figure 3A.
Using the Scheil equations, or Figure 5, we can predict the shape of concentration profiles in the direction of solidification at different values for f s (Figure 6). Note that the prediction is unphysical for high f s .
image: 47_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_13_14_Fig6.svg
Figure 11.6: Concentration profiles in the direction of solidification for the system in Figure 5. Note that the resulting solid is graded in terms of its composition. However, the average concentration of the minority component in the solid, X s , approaches X 0 for high f s .
Finally, let's consider the scenario that there is no mixing in the solid, and only diffusion, but no convection in the liquid phase. As before, the first solid formed has composition k X 0 . (Figure 7A). Solute rejected from the solid enters the liquid. It is transported away from the interface by diffusion. At the interface, the mole fraction of the solute is higher than the initial value. This is referred to a solute pile-up. As a consequence of the pile-up, the temperature has to drop before more solid can form. Once it does, the solid that precipitates will have a higher mole fraction of solute. However, some solute is still rejected, driving up the mole fraction of the solute in the liquid further. One can show that after an initial transient, a steady state is established (Figure 7B). At steady state, the mole fraction of the solute in the solid is X 0 , and that in the liquid is X 0 k -1 . The temperature is the liquidus temperature T L ( X 0 k -1 ) .
image: 48_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_13_14_Fig7.svg
Figure 11.7: Concentration profile for a unidirectional solidification with no mixing in the solid and diffusive transport only in the liquid. A. As solute is rejected from the first-precipitated solid, it piles up on the liquid side of the interface. As a result, the temperature needs to drop for further soldification to occur, and the mole fraction of the solute in the solid increases during the initial transient. B. After some time, a steady state is reached where the temperature of the interface is constant, the composition of the solid is X 0 , and the composition of the liquid in local equilibrium with the solid is X 0 k -1 .
In steady state, the mole fraction of the solute on the liquid side of the interface at steady state is decays from X 0 k -1 right at the interface to X 0 far from the interface. One can show that the the composition of the liquid as a function can be written as X L ( x ) = X 0 [ 1+ 1-k k e - v D x ] (11.11) = X 0 +[ X 0 k - X 0 ] e - v D x (11.12) , where v is the interface velocity, D is the diffusivity of the solute in the liquid, and x is the distance from the interface.
Having an analytical solution for the local mole fraction of the solute, or, in other words, knowing the shape of the pile-up, allows us to predict the local liquidus temperature, T L ( x ) . The local liquidus temperature simply tells us at what temperature we would expect solidification to be possible given the local mole fraction of the solute. Using Eq. (1), we write
T L ( X L ( x ) ) = T m α + b L X L ( x ) (11.13) (11.14)
Using Eq. (12), T L ( x ) = T m α + b L X 0 + b L [ X 0 k - X 0 ] e - v D x (11.15) = T m α + b L X 0 +[ b L X 0 k - b L X 0 ] e - v D x (11.16) Note that T m α + b L X 0 = T L ( X 0 ) . T L ( x ) = T L ( X 0 ) +[ b L X 0 k - b L X 0 ] e - v D x (11.17) Next, let's expand the square brackets by adding T m α - T m α =0 T L ( x ) = T L ( X 0 ) +[ T m α + b L X 0 k -( T m α + b L X 0 ) ] e - v D x (11.18) Using T m α + b L X 0 k = T L ( X 0 k ) and, once more, T m α + b L X 0 = T L ( X 0 ) , we rewrite as T L ( x ) = T L ( X 0 ) +[ T L ( X 0 k ) - T L ( X 0 ) ] e - v D x (11.19)
From Eq. 19, we see that the exponential decay in solute mole fraction with in creasing distance x from the interface (Figure 8A) causes a corresponding increase in the local liquidus temperature (Figure 8B). At the interface, x=0 , and T L ( x ) = T L ( X 0 k ) . Far from the interface, x and T L ( x ) T L ( X 0 ) . If the local temperature T( x ) is below the local liquidus temperature, the system is locally supercooled. Because this supercooling is a consequence of the local concentration of the constituents (components) of the binary, this phenomenon is referred to as constitutional supercooling.
Constitutional supercooling is a requirement for interface instability in binary system. Consider a system where the actual temperature in the solid and liquid is given by T( x ) . If the thermal conductivities are approximately equal, we expect that the temperature gradient in the solid is a bit steeper than in the liquid, resulting in a positive interface velocity. A protrusion that forms on the interface may thus jut into liquid that has a lower local mole fraction of solute than at the flat interface, a higher local liquidus temperature, and thus a higher supercooling. The protrusion could then grow more rapidly, resulting in an unstable interface. The only way to avoid constitutional supercooling is my increasing the temperature gradient in the liquid such that it is steeper than the slope of the local liquidus at the interface.
image: 49_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_13_14_Fig8.svg
Figure 11.8: Constitutional undercooling in binary systems. A. Plot of local composition X L ( x ) of the liquid ahead of the interface (red line) in unidirectional solidification in steady state. B. Plot of the local liquidus temperature T L ( x ) (blue line) for the same system as in (A). If the actual temperature in the solid and liquid is described by T( x ) (green line), the local actual temperature is below the local liquidus temperature in the shaded region. The system is said to be constitutionally supercooled in this region. Constitutional supercooling can only be avoided if the gradient of the temperature at the interface exceeds a critical value (green dashed line).
We can determine this slope from Eq. 19: . dT dx | crit . dT( x ) dx | x=0 = T L ( X 0 ) - T L ( X 0 k ) D v (11.20)
Note that the numerator is a characteristic of the phase diagram and initial composition. The denominator is the characteristic thickness of the diffusion layer on the liquid side of the interface, i.e. the length x over which the exponential term falls from 1 to e -1 .

Precipitate Growth

In this section, we will look at growth of a small precipitate of β phase that shall be rich in component 2, from a matrix α that is rich in component 1. We know that if growth is isothermal at some temperature T< T solvus , and the initial concentration of component 2 in α is C 0 , then we can find the equilibrium concentration of component 2 in the matrix, C e q α , and that in the precipitate C e q β using a tie line in the phase diagram (Figure 1). We shall define the supersaturation Δ C C 0 - C e q α and the undercooling as Δ T T solvus ( C 0 ) -T .
image: 50_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_16_Fig1.svg
Figure 12.1: Phase diagram for β precipitate growing in α matrix with initial concentration C at some temperature T .
Many precipitates are plate, disk, or needle shaped, with coherent interfaces and incoherent interfaces. Because of strain at the the coherent interface, the local supersaturation is lower and growth is typically slower than at incoherent interfaces (Figure 2a).
image: 51_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_16_Fig2.svg
Figure 12.2: A. Schematic drawing of a precipitate growing only in the direction of the incoherent interface, i.e. the x -direction. The curvature of the interface is ignored. The origin is in the center of the precipitate, and the position of the interface at time t is given by h( t ) . B. Concentration profile for component 2 in the x -direction. Note that mass conservation dictates that the shaded areas A 1 and A 2 are identical.
Consider the concentration of component 2 along a hypothetical line from the midpoint of the precipitate in the growth direction. A plot of this concentration C 2 versus the distance x is called a concentration profile (Figure 2B). From x=0 right up to the interface at x=h( t ) , we are inside the precipitate and C 2 = C e q β . Just to the right of the interface, the matrix is in local equilibrium with the precipitate, and C 2 = C e q α . At large x , the concentration of component two should approach the initial concentration C 2 C 0 . Thus, there is a concentration gradient. Units of component 2 diffuse down this gradient until the local concentration at the interface is larger than the equilibrium value. This allows the precipitate to grow and return to local equilibrium. At any given time, mass conservation dictates that the shaded areas in Figure 2B are the same. Therefore, h( t ) ( C e q β - C 0 ) = h( t ) ( C 0 -C( x ) ) x (12.1)
To simplify this a little bit, let's assume we are in quasi steady state. This means that the concentration gradient on the right side of the interface is a linear function of the distance (Figure 3A).
image: 52_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_16_Fig3.svg
Figure 12.3: A. Concentration profile in quasi-steady state. B. To move the interface by an infinitesimal amount dh to the right, an amount of component 2 corresponding to the shaded area has to diffuse down the concentration gradient.
Therefore, h( t ) ( C e q β - C 0 ) = 1 2 ( C 0 - C e q α ) L= 1 2 Δ CL (12.2) L = 2( C e q β - C 0 ) h( t ) Δ C (12.3)
To grow the precipitate by an infinitesimal volume element dV=Adh , where A is the cross sectional area and dh the infinitesimal change in thickness, we have to convert an volume element of the matrix α at C e q α to β at C e q β . This requires that the following number of units of component 2 enter the volume of α from the right. d N 2 =-( C e q β - C e q α ) Adh (12.4)
Component 2 units arrive at the interface by diffusion down the concentration gradient. In the general case, this number can be calculated if we know the flux, the cross sectional area, and the infinitesimal amount of time dt that it takes for the interface to move. d N 2 = J 2 Adt=-D C x Adt (12.5)
Therefore, -( C e q β - C e q α ) Adh=-D C x Adt (12.6) dh dt = D C x ( C e q β - C e q α ) (12.7)
Assuming steady state, we can substitute C x = Δ C L dh dt = D Δ C ( C e q β - C e q α ) L (12.8)
Using Eq. (3), hdh dt = D Δ C 2( C e q β - C e q α ) ( C e q β - C 0 ) (12.9) hdh= D Δ C 2( C e q β - C e q α ) ( C e q β - C 0 ) dt (12.10) hdh = D Δ C 2( C e q β - C e q α ) ( C e q β - C 0 ) dt (12.11) 1 2 h 2 = Δ CDt 2( C e q β - C e q α ) ( C e q β - C 0 ) (12.12) h= Δ C Dt ( C e q β - C e q α ) ( C e q β - C 0 ) (12.13)
Most of the time, growth will occur at small supersaturation. Assuming Δ C0 (12.14) C 0 - C e q α 0 (12.15) C 0 C e q α (12.16) , we can replace the term ( C e q β - C 0 ) in the denominator of Eq. 13 by ( C e q β - C e q α ) . h( t ) Δ C Dt ( C e q β - C e q α ) ( C e q β - C e q α ) (12.17) Δ C Dt ( C e q β - C e q α ) (12.18)
Using this approximation, we can write the interface velocity v( t ) = dh( t ) dt Δ C 2( C e q β - C e q α ) D t (12.19)
With the final assumption that the molar volume V m is independent of the concentration of the components C i , we can replace C i with X i .
h( t ) Δ X Dt ( X e q β - X e q α ) (12.20) v( t ) Δ X 2( X e q β - X e q α ) D t (12.21)
Note that the thickness of the interface h( t ) Δ X Dt , meaning that growth is faster at higher supersaturation. However, for the thickness to grow to twice its value, the time required increases by a factor of 4! For the interface velocity, v( t ) Δ X D t , meaning again that the velocity increases directly proportional to the supersaturation. However, the growth speed decreases over time! This is because the concentration gradient in the depletion zone gets increasingly shallow. Finally, through the diffusivity, temperature has a very strong impact on interface position and velocity.
Q. Sketch a plot of the temperature (y-axis) vs. the interface velocity (x-axis). Which variables depend (strongly) on T ? What is the (approximate) proportionality of the individual terms?

Growth and Curvature

In this section, we will consider the effect of curvature at the interface of a precipitate and the matrix. Curvature effects are universal and important both in condensed matter and the gas phase. We will here only consider binary systems. We are interested in the equilibrium of a spherical particle with radius r of phase β in a matrix α (Figure 1a). As before, we can describe this situation in a GX diagram (Figure 1b).
image: 53_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_17_18_Fig1.svg
Figure 13.1: A. Schematic drawing of a β phase precipitate with radius r in equilibrium with the surrounding α matrix. The curved interface gives rise to a Laplace pressure Δ P= P β - P α . B. GX diagram of α in equilibrium with a β -precipitate with a flat interface ( r= , blue common tangent) and in equilibrium with a β -precipitate with a curved interface of radius r (red common tangent). The points of tangency reveal the mole fraction of component 2 in the matrix and the precipitate for the two equilibria.
While the molar free energy of the α phase is described by G m α , that of β depends on the radius of curvature. If β is incompressible, the molar free energy of a spherical particle with radius r is shifted towards higher free energy by Δ P V m = 2 γ V m r (13.1)
Using the common tangent approach, we can determine the equilibrium mole fraction of component 2 in the two phases, X r α , and X r β . Note that for r , the system goes to its final equilibrium. Therefore, X α = X e q α and X β = X e q β and the nominal supersaturation is Δ X= X r α - X α .
Recall that for incompressible systems, the free energy of formation of a spherical particle of phase β from a matrix α can be expressed as Δ G m α α + β =- . 2 G m X 2 2 | X α Δ X( X β - X α ) (13.2)
Furthermore, we found that Δ G m α α + β =- Δ P V m (13.3) 2 γ V m r = . 2 G m X 2 2 | X α Δ X( X β - X α ) (13.4)
How can we evaluate the second derivative? Rewrite 2 G m X 2 2 = X 2 [ G m X 2 ] (13.5)
The first derivative of the molar free energy is the difference in chemical potentials 2 G m X 2 2 = X 2 [ μ 2 - μ 1 ] (13.6) = X 2 [ RT ln a 2 -RT ln a 1 ] (13.7) , which can be expressed in terms of the activities of components 1 and 2. Now if we consider ideal dilute solutions, we can use Rault's law for the solvent lim X 2 0 a 1 = X 1 =1- X 2 (13.8)
, and Henry's law for the solute lim X 2 0 a 2 = Γ X 2 (13.9) , where Γ is the activity coefficient.
Substituting (8) and (9) into (7), 2 G m X 2 2 = X 2 [ RT ln ( Γ X 2 ) -RT ln ( 1- X 2 ) ] (13.10) =RT X 2 [ ln Γ + ln X 2 - ln ( 1- X 2 ) ] (13.11)
Using lim x0 ln ( 1-x ) =-x (why? write down Taylor expansion and you will know), and differentiating 2 G m X 2 2 =RT[ 0+ 1 X 2 +1 ] (13.12)
In the limit of small X 2 , the second term in the sum is much larger that the third, and we can write . 2 G m X 2 2 | X α RT X α (13.13)
With the final assumption that X r β X β 1 and X α 0 , the difference X β - X α 1 . Substituting into Eq. (4) and rearranging, we arrive at the Gibbs-Thomson equation. 2 γ V m r = X α Δ X (13.14) Δ X= 2 γ V m RTr X α (13.15) X r α = 2 γ V m RTr X α + X α (13.16) X r α X α = 2 γ V m RTr +1 (13.17)
The Gibbs-Thomson equation is a very important result. It allows us to predict the equilibrium mole fraction of components in a matrix that is in equilibrium with a precipitate of a given radius (of curvature). From Eqs. 15 and 17, we see that X r α X α (13.18) X r α 1 r , Δ X 1 r (13.19) r X r α X α (13.20)
In other words, the mole fraction of component 2 in the matrix that is in equilibrium with a small precipitate (small radius of curvature) is higher than that in a matrix that is in equilibrium with a large precipitate (large radius of curvature). We can also say that the matrix in equilibrium with the smaller precipitate has higher nominal supersaturation. Note, however, that the nominal supersaturation is calculated with respect to a precipitate with infinite radius.
Sometimes it is more useful to consider the effective supersaturation of the matrix with respect to a precipitate of a given radius. Let's assume that the mole fraction of component 2 in the matrix is X 0 α , and that at this mole fraction precipitates of radius r * are in equilibrium, i.e. X r * α = X 0 α . The effective supersaturation of the matrix for a particle with radius r is then Δ X e ff = X r * α - X r α (13.21) =[ 2 γ V m RT r * X α + X α ] -[ 2 γ V m RTr X α + X α ] (13.22) = 2 γ V m RT X α [ 1 r * - 1 r ] (13.23) = 2 γ V m RT X α [ 1 r * - 1 r ] r * r * , (13.24) = 2 γ V m RT r * X α [ 1- r * r ] (13.25) = Δ X 0 [ 1- r * r ] (13.26) where Δ X 0 = X r * α - X α = 2 γ V m RT r * X α (13.27) is the nominal supersaturation of the matrix in equilibrium with a ppt of radius r * .
image: 54_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_17_18_Fig2.svg
Figure 13.2: Schematic drawing of three β precipitates with different radii in the α matrix. Only the particle with radius r * is in equilibrium with the matrix with mole fraction X 0 = X r * α . Consequently, Δ X e ff =0 for this precipitate. Eq. 26 tells us that Δ X e ff <0 for the small precipitate, and Δ X e ff >0 for the large precipitate. A negative effective supersaturation results in dissolution of the precipitate until local equilibrium is established, and positive effective supersaturation will lead to growth.
Inspection of Eq. 26 reveals that the effective supersaturation is positive for r> r * , equal to zero for r= r * , and negative for r< r * . This means that only particles with a radius that is larger than r * will grow, thereby take up component 2 from the supersaturated matrix, and thus lower the local supersaturation until equilibrium is reached. Ppt with radii smaller than r * on the other hand will shrink, releasing component 2 into the matrix and reducing the local undersaturation in the matrix. This is sometimes expressed as saying that the precipitate dissolves in the matrix. Similarly, we can say that the Gibbs-Thomson equation predicts that the solubility of small precipitates is higher than that of large precipitates.
Only ppt that are exactly at radius r will neither shrink nor grow. Note that the equilibrium that the latter precipitates are in is unstable, meaning that if they grow by even a tiny amount in a random fluctuation, they will keep growing. This is exactly the same situation that we discussed for the nucleus. Note that for nucleation, we typically deal with high supersaturation and consider the formation of one nucleus while the matrix composition remains constant. In growth, we consider that many precipitates are already present, and the supersaturation is much lower.
Q: Plot the ratio of Δ X e ff over Δ X 0 vs. r r * .
In the previous paragraphs, we considered precipitates of different radii in a matrix with a fixed mole fraction of component 2, essentially treating them as independent systems. What happens if we instead assume that the matrix shall have an average (but not constant) composition X 0 = X r * α , and that there are precipitates of different sizes distributed in the matrix? Consider the hypothetical situation in Figure 3, where three precipitates, one with r 1 < r * , one with r 2 = r * , and one with r 3 > r * , sit next to each other.
image: 55_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_17_18_Fig3.svg
Figure 13.3: A. Schematic drawing of three β precipitates with different radii in the α matrix. Each precipitate is in local equilibrium, meaning that the matrix just outside the interface has the composition X r i α . The average composition of α is about X r 2 α . B. Concentration profile along the dashed line in (A). Because X r 1 α > X r 2 α , there is a concentration gradient from the interface of the small precipitate to the surrounding matrix that results in a diffusive flux of component 2 away from the interface. This results in dissolution of the precipitate and the interface move inward. For the large precipitate, X r 3 α < X r 2 α . As a consequence, component 2 diffuses toward the interface, and the precipitate grows in order to maintain equilibrium.
For each of the precipitates, the composition of the matrix in local equilibrium is set by the Gibbs-Thomson equation. This means that the mole fraction of component 2 in the matrix close to the smallest precipitate is higher than the average composition, that in matrix in local equilibrium with the intermediate precipitate is exactly equal to the average composition, and that in the matrix in local equilibrium with the large precipitate is lower than the average value.
image: 56_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_17_18_Fig4.svg
Figure 13.4: A. Schematic drawing of spherical particle, defining the abscissa ( ρ ) for the radial concentration profiles in B and C. B. Plot of the radial concentration profile for a spherical precipitate in a supersaturated matrix. If the particle is in local equilibrium, the concentration of the components in the matrix just outside the particle are given by the Gibbs-Thomson equation. C. Assuming quasi steady state, the concentration gradient becomes linear over the distance L .
We now know how to determine whether a precipitate shrinks, remains the same, or grows, and what the effective supersaturation is. The next step is to find the interface velocity. In order to account for curvature, we need to consider growth in three dimensional space. For a spherical precipitate, it is sufficient to determine the radial growth velocity. To do so, consider the concentration profile from the center of the precipitate in the radial direction (Figure 4A). While the profile looks very similar to the disk, plate, or needle shaped precipitates we considered previously, the shaded areas are not identical in size. Based on the shape of the precipitate, conservation of mass, and using a quasi-steady-state approximation (Figure 4B), however, it is possible to determine the size of the depletion zone L=kr (13.28) , where k is a geometric factor that takes the value k=1 for a spherical precipitate, and r( t ) is the radius of the precipitate at time t . L is therefore dependent on the radius of the precipitate and will increase (or decrease) as it grows (or shrinks). In this case, the radial interfacial velocity takes a form that is very similar to the growth velocity we derived for the unidirectional movement of a flat interface. v= dr( t ) dt = D C β - C r α d C 2 d ρ D C β - C r α Δ C 2 Δ ρ (13.29)
Using Figure 4B to determine the slope of the concentration gradient, v D C β - C r α C α - C r α L = D C β - C r α Δ C e ff kr( t ) (13.30) , where Δ C e ff is the effective supersaturation of the matrix with respect to a precipitate with radius r (from here on, we use r instead of r( t ) even though r remains a function of time). If we assume that the molar volume is independent of the concentration, we can replace concentration by mole fraction and write
v D X β - X r α Δ X e ff kr (13.31)
Using Gibbs-Thomson, we can express the effective supersaturation as a function of the radius and the nominal supersaturation.
v D Δ X k( X β - X r α ) 1 r ( 1- r * r ) (13.32)
As before, r * is the radius of a precipitate that is in equilibrium with a matrix that has the average composition X r * α X α . Assuming that X r α X α , and for a spherical precipitate where k=1
v D Δ X ( X β - X α ) 1 r ( 1- r * r ) (13.33)
Inspection of Eq. and comparison with the growth velocity of a precipitate with a flat interface (see section 11), v Δ X 2( X e q β - X e q α ) D t (13.34) reveals that in both cases the velocity v Δ X . For the curved interface, vD , whereas for the flat interface v D . In either case, the temperature, though the diffusivity, has a strong impact on the growth velocity (what other variables are dependent on T?). Finally, for the flat interface, v t - 1 2 . For the curved interface, the velocity also depends on r , which of course depends on time as well.
v 1 r ( 1- r * r ) (13.35)
Inspection of Eq (35) reveals that the growth velocity is negative for r< r * , zero for r= r * , and positive for r> r * , as expected from the discussion of the effective supersaturation.
Q. Sketch the dependence of growth velocity v on r and show that the velocity has a global maximum at r=2 r * (Figure 5).
image: 57_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_17_18_Fig5.svg
Figure 13.5: Dependence of growth velocity of a curved interface on the radius of curvature.

Coarsening

In this section, we will consider what happens after nucleation and some growth have reduced the supersaturation to the point that the nucleation current has become negligible. Let's assume that we have a binary system with β precipitates (rich in component 2) in an α matrix (rich in component 1), and that prior nucleation and growth have resulted in a number of spherical ppt with radius r i (Figure 1).
image: 58_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_19_Fig1.svg
Figure 14.1: A. Schematic depiction of N=100 spherical β precipitates with radius r i in α matrix. B. Histogram of particle radii in (A).
For a system with N precipitates, the average precipitate radius is then r= 1 N 1 N r i (14.1) Let's further assume that all ppt shall be in local equilibrium with the matrix. With Gibbs-Thomson, c r i α c α = 2 γ V m β RT r i +1 (14.2) c r i α = 2 γ V m β c α RT r i + c α (14.3)
The average concentration of component 2 in the matrix is then c c r α = 2 γ V m β c α RTr + c α (14.4)
Now let's consider the growth velocity of a precipitate with radius r i v= dr dt =- J 2 c β - c r i α (14.5) , where J 2 is the flux of component two towards the interface. Using Fick's first law, v= D c β - c r i α c 2 x (14.6) , where D is the diffusivity of component 2 in the matrix α . Assuming quasi-steady state, i.e. a linear concentration gradient from the interface to the bulk matrix, v= D c β - c r i α Δ c 2 Δ x = D k r i c - c r i α c β - c r i α (14.7) , where k=1 for spherical particles.This is known as the Zener Growth Law.
Using Eq. 4 and 7, v= D r i ( c β - c r i α ) 2 γ V m β c α RT ( 1 r - 1 r i ) (14.8)
If we further assume that the supersaturation is small, c r i α c α , and therefore v= 2 γ D V m β c α RT( c β - c α ) 1 r i 2 ( r i r -1 ) (14.9)
Inspecting Eq. 9, we find that the growth velocity is positive for r i >r , zero for r i =r , and negative for r i <r ! Only precipitates that are larger than average will grow, those that are smaller will shrink. The growth velocity has a maximum for precipitates that are twice the average size, and decreases again for larger precipitates (show that this is true).
Based on this, and a similar growth law for interface-controlled growth, Lifshitz, Slyosow, and Wagner (LSW) derived a theory that describes the time evolution of the average precipitate radius, using the following assumptions:
Under these conditions, Lifshitz and Slyosow found for diffusion controlled growth r( t ) 3 =r( t=0 ) 3 + 8 γ V m β D c α 9RT( c β - c α ) t (14.10) r( t ) t 1 3 (14.11)
image: 59_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_19_Fig2.svg
Figure 14.2: Diffusion-limited growth of cementite precipitates in α -Fe matrix during coarsening. Log-log-plot of the ratio of the average precipitate radius at time t over the average radius at t=0 against time at different temperatures.
For interface controlled growth, Wagner found r( t ) 2 = 64 γ V m β c α 81RT( c β - c α ) k s t (14.12) r( t ) t 1 2 (14.13)
Note that the diffusivity D only appears in the former, the attachment rate k s only in the latter equations. In either case, the rate by which the average radius grows decreases with time.
Note that these equations only describe how the average precipitate radius evolves with time. One might ask what the particle size distribution is. We will not discuss this in detail, but one interesting result is that for diffusion controlled processes, one can show that the fraction Ψ ( ρ ,t ) of particles with reduced radius ρ = r r is independent of time (Figure 3).
image: 60_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_19_Fig3.svg
Figure 14.3: Precipitate size distribution is independent of time. There is a maximum at ρ 1.13 , meaning that the most frequently observed precipitates are slightly larger than the average. The probability for larger radii drops off quickly and goes to zero at ρ =1.5 . Note that this size distribution is a prediction based on theory; actual size distributions can and do differ.
The LSW theory of coarsening can be applied to many different materials, including metals, ceramics, polymers, and biomaterials, and holds for liquids, solids, and with some limitations, gases.
Q: Discuss how coarsening in diffusion controlled systems is affected by temperature

Kinetics of Phase Transformations

In this section, we will connect microscopic events, such as nucleation and growth of small volumes of a new phase, to the macroscopic rate of transformation. Let's consider a simple case of α β , meaning that the entire volume of α is transformed. Examples for such a transformation would be solidification from the melt in a unary system, e.g. freezing of water, or the transformation of γ -Fe (austenite, FCC) to α -Fe (ferrite, BCC) in pure iron. This phase transformation shall occur isothermally, i.e. at a constant undercooling. Let's further assume that nucleation is random in the entire volume and occurs at a steady state rate N ˙ V . Once a spherical nucleus has formed, it shall grow with a constant and isotropic growth velocity v .
image: 61_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_20_Figure1.svg
Figure 15.1: Simulation of a α β phase transformation in 2D with constant nucleation and growth rate. A. At time τ , four nuclei are present. B. At time τ +d τ , the four nuclei formed in the previous interval have grown into particles (gray circles), and four new nuclei have formed randomly in the matrix. Growth and nucleation continue in this fashion until growing β -particles impinge on each other (E). Areas highlighted in red indicate the difference between the extended area, i.e. the area that particles would take if their growth would continue unimpeded, and the actual area of β formed.
Given these conditions, we can understand the transformation of α into β by considering what is going on in short time intervals (Figure 1). During very first interval (not shown in Fig 1), a number of nuclei appears in the matrix. In the next interval, the nuclei grow in radius, and a new batch of nuclei forms. In the following interval, again new nuclei form and the existing β -particles grow. This can continue for some time, but eventually, neighboring particles will impinge on each other and will no longer be able to grow evenly in all directions (Figure 1E).
We would like to determine the volume-based, fractional conversion f( t ) = V β ( t ) V t ot (15.1)
However, because of the problem that we eventually run out of α phase for the β particles to grow into, the time dependence of V β ( t ) is not immediately obvious. Ignoring the problem for now, lets define the extended volume of β , i.e. the volume β would take if there was an infinite amount of α such that β particles never impinge on each other. If there are n β -particles, then the extended volume is simply the sum over all of them
V ex β ( t ) = i n V i β (15.2)
The extended fractional conversion is then f ex ( t ) = V ex β ( t ) V t ot = V ex β ( t ) V α ( t ) + V β ( t ) (15.3) , where V α ( t ) = V t ot - V β ( t ) is the remaining volume that β particles can grow into.
We can write the available amount of α also in terms of the available volume fraction 1-f( t ) =1- V β ( t ) V t ot (15.4)
To find the fractional conversion we simply multiply the extended fractional conversion with ( 1-f ) to account for the reduction in the available α phase.
df=d f ex ( 1-f ) (15.5) df 1-f =d f ex (15.6)
Using df=-d( 1-f ) , - d( 1-f ) 1-f =d f ex (15.7)
With f( 0 ) =0 , we can integrate - d( 1-f ) 1-f = f ex (15.8) ln ( 1-f ) =- f ex (15.9) f=1- e - f ex (15.10)
This is a useful result because it allows us to determine the fractional conversion in a straightforward manner if we can find an expression for the extended fractional conversion. Let's look at the specific example outline above, where the nucleation rate and the growth rate are both constant.
In any time interval d τ , we therefore expect dN= N ˙ V Vd τ (15.11) nuclei to form in the volume V. Any nucleus formed at time τ and observed at time t will have grown by dV= 4 π 3 ( r t 3 - r τ 3 ) (15.12)
Because r t r τ , and with r t =v( t- τ ) , dV 4 π 3 ( v( t- τ ) ) 3 (15.13)
The contribution to the extended volume of all nuclei formed during a time d τ at some time τ is simply the product: d V ex β =dVdN= 4 π 3 ( v( t- τ ) ) 3 N ˙ V Vd τ (15.14) d V ex β V = 4 π 3 N ˙ V v 3 (t- τ ) 3 d τ (15.15)
Integrating up until the time of observation, d V ex β V = 4 π 3 N ˙ V v 3 0 t (t- τ ) 3 τ (15.16) f ex ( t ) = π 3 N ˙ V v 3 t 4 (15.17)
Combining Eqs. 10 and 17, we can write the fractional conversion f( t ) =1- e -k t 4 (15.18) , where k= π 3 N ˙ V v 3 (Figure 2A).
image: 62_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_20_Figure2.svg
Figure 15.2: A. Plot of the fractional conversion vs. logarithmic time for the conversion of α to β described by Eq. (18), with N ˙ V =1 m -3 s -1 , v=100 μ m s -1 , k=1 0 -12 s -4 , and n=4 . B. Plotting lg ln (1-f ) -1 against lg t linearizes the JMAK equation. A linear function (red dashed line) fitted to data points indicated by blue circles in (A) has slope n=4 and intercept lg k=-12 .
This is but one example of the general form of the Johnson-Mehl-Avrami-Kolmogorov (JMAK) equation that describes macroscopic phase transformations. f( t ) =1- e -k t n , (15.19) where k and n are generally fit parameters. Note that it is possible to derive a JMAK equation with analytical expressions for k and n based on an arbitrarily complex model for the nucleation and growth processes in the sample. Where analytical solutions fail, it may be possible to use numerical methods. However, these solutions are not unique - more than one model could result in a given combination of k and n . Fitting experimental data to the JMAK equation therefore DOES NOT give us information on the mechanism of phase transformation.
Fitting the JMAK equation used to require linearizing it for regression. While this is no longer necessary, it is still useful. To linearize, we rewrite Eq. 19 (Figure 2B). ln ( 1-f ) =-k t n (15.20) lg ln 1 ( 1-f ) =n lg t+ lg k (15.21)

Spinodal Decomposition

In the previous section, we realized that if the system is in the spinodal regime, the free energy of unmixing is negative. However, for there to be a diffusive flux, there needs to be a concentration gradient, even if uphill diffusion is possible. In principle then, as long as the concentration profile is completely flat, the flux should be zero and no phase separation should occur. However, components will still undergo thermally activated diffusion, and at any one moment in time, a one-dimensional concentration profile will show very small, random deviations from the average concentration (Fig 1 A).
r image: 63_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_21_23_Fig1.svg
Figure 16.1: A. Due to thermal fluctuations, the concentration profile of any binary system (here: one dimensional) is never completely flat. B. As any concentration profile can be thought of as a superposition of sine and cosine functions, we will consider the free energy change associated with changing the concentration profile from a completely flat line to a sine wave with wavelength λ and amplitude δ .
To understand how these local fluctuations affect the stability of the system, we model a concentration fluctuation by a sine function with wavelength λ and amplitude δ (Fig. 1B). C- C 0 = δ sin 2 π z λ (16.1)
Next, we determine the change in the Gibbs free energy per unit volume that is associated with the change from a flat concentration profile to the sine-shaped fluctuation. To do that, we simply average the concentration-dependent Gibbs free energy over one period and subtract the Gibbs free energy at C 0 . Δ G V = 1 λ 0 λ G V ( C ) z - G V ( C 0 ) (16.2)
To evaluate the integral, write the Taylor expansion of the Gibbs free energy around the average concentration and terminate after the second order term G V ( C ) = G V ( C 0 ) + . G V C | C 0 ( C- C 0 ) + 1 2 . 2 G V C 2 | C 0 (C- C 0 ) 2 (16.3)
Substituting Eq (3) into (2), Δ G V = 1 λ 0 λ [ G V ( C 0 ) + . G V C | C 0 ( C- C 0 ) + 1 2 . 2 G V C 2 | C 0 (C- C 0 ) 2 ] z - G V ( C 0 ) (16.4) = 1 λ . G V C | C 0 0 λ ( C- C 0 ) z + 1 2 λ . 2 G V C 2 | C 0 0 λ (C- C 0 ) 2 z (16.5)
Using Eq (1), Δ G V = 1 λ . G V C | C 0 0 λ δ sin 2 π z λ z + 1 2 λ . 2 G V C 2 | C 0 0 λ δ 2 sin 2 2 π z λ z (16.6)
The first term in Eq. (6) is equal to zero, as we are integrating over an entire period. Therefore, Δ G V = δ 2 2 λ . 2 G V C 2 | C 0 0 λ sin 2 2 π z λ z (16.7) = δ 2 2 λ . 2 G V C 2 | C 0 λ 2 (16.8) = δ 2 4 . 2 G V C 2 | C 0 (16.9)
This is an interesting result that confirms our qualitative assessment. As long as the second derivative of the Gibbs free energy is negative, any fluctuation with an amplitude δ greater than zero lowers the free energy of the system. The wavelength λ doesn't even occur in Eq. 9. Is this reasonable? A fluctuation with a short wavelength is equivalent to a rapid concentration change over a small distance, i.e. a steep concentration gradient (Fig. 2).
image: 64_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_21_23_Fig2.svg
Figure 16.2: As wavelength lambda decreases, increasingly steep concentration gradients (dashed lines) result.
It is reasonable to assume that there is an energetic price for establishing a steep gradient. We will therefore add a term to Eq. 2 that depends on the magnitude of the concentration gradient. Δ G V = 1 λ 0 λ [ G V ( C ) + κ ( C z ) 2 ] z - G V ( C 0 ) , (16.10) where κ >0 is a proportionality constant that expresses how strongly the free energy depends on the magnitude of the concentration gradient.
Evaluating the integral is straightforward Δ G V = δ 2 4 . 2 G V C 2 | C 0 + κ λ 0 λ ( C z ) 2 z (16.11) = δ 2 4 . 2 G V C 2 | C 0 + κ λ 0 λ ( z [ δ sin 2 π z λ ] ) 2 z (16.12) = δ 2 4 . 2 G V C 2 | C 0 + κ δ 2 λ 0 λ ( 2 π λ cos 2 π z λ ) 2 z (16.13) = δ 2 4 . 2 G V C 2 | C 0 + 4 π 2 κ δ 2 λ 3 0 λ cos 2 2 π z λ z (16.14) = δ 2 4 . 2 G V C 2 | C 0 + 4 π 2 κ δ 2 λ 3 λ 2 (16.15) = δ 2 4 . 2 G V C 2 | C 0 + 2 π 2 κ δ 2 λ 2 (16.16) = δ 2 2 [ . 2 G V 2 C 2 | C 0 + 4 π 2 κ λ 2 ] (16.17)
Now if the second derivative of the Gibbs free energy is negative, and with κ >0 , the two terms in the parentheses have opposite sign, indicating that Δ G V changes sign for some critical value of λ . To find λ C , we write Δ G V ( λ C ) =0 (16.18) δ 2 2 [ . 2 G V 2 C 2 | C 0 + 4 π 2 κ λ C 2 ] =0 (16.19) . 2 G V 2 C 2 | C 0 =- 4 π 2 κ λ C 2 (16.20) λ C = π - 8 κ . 2 G V C 2 | C 0 (16.21)
Next, we determine for which λ spinodal decomposition is spontaneous: Δ G V ( λ C ) <0 (16.22) λ > π - 8 κ . 2 G V C 2 | C 0 (16.23) λ > λ C (16.24)
This means that only fluctuations with a wavelength greater than the critical wavelength λ C will form spontaneously. A rule of thumb is that in metals, λ C is on the order of 50 nm. In polymers, it is on the order of 1 μ m.
Note that in spindle decomposition, there are no sharp phase boundaries or interfaces. Instead, there are gradual changes in the composition.
A continuous Gibbs free energy curve implies that only the composition, but not the structure of the material changes, meaning that the system is coherent. As a consequence, we expect that there will be strain as differences in composition will result in local variation of the lattice parameter. Recall that we used the elastic strain energy density W V e l when we looked at the effect of strain on nucleation. For the elastic strain energy in spinodal decomposition we use W V e l 6 ε 2 μ (16.25)
The misfit parameter is the relative change in lattice parameter ε = Δ a a (16.26) where Δ a= a C Δ C (16.27) With η 1 a a C , (16.28) we can write the misfit parameter as a function of Δ C=C- C ε = η ( C- C ) , (16.29) and finally the elastic strain energy density W V e l 6 η 2 μ ( C- C ) 2 (16.30)
We can then account for strain by adding a term to Eq. (10) Δ G V = 1 λ 0 λ [ G V ( C ) + κ ( C z ) 2 +6 η 2 μ ( C- C ) 2 ] z - G V ( C 0 ) (16.31)
Integration is straightforward and we find Δ G V = δ 2 2 [ . 2 G V 2 C 2 | C 0 +6 η 2 μ + 4 π 2 κ λ 2 ] (16.32) λ C = π - 8 κ . 2 G V C 2 | C 0 +6 η 2 μ (16.33) Inspection reveals that the (positive) strain energy contribution shifts the free energy change to more positive values, resisting spinodal decomposition. This is equivalent to shifting the binodal and spinodal lines to lower temperature in the phase diagram. The strain energy term further makes the denominator in Eq (33) less negative. The critical wavelength λ C therefore increases with increasing strain energy.
Take a moment to reflect on the similarities and differences of nucleation and spinodal decomposition. In nucleation, a small amount of a new phase forms that has a much higher concentration of one component. This is equivalent to a very high amplitude fluctuation. Because the second derivative of the free energy is positive, diffusion against the concentration gradient increases the free energy of the system and the process is not spontaneous. There is an interface between nucleus and matrix and the formation of the nucleus is opposed by the interfacial free energy. The critical radius indicates the size of a particle that is in unstable equilibrium, meaning that both growing and shrinking will lower its free energy. As soon as a particle becomes supercritical, it is likely to keep growing. If nucleus and matrix are coherent, then strain will oppose the formation of a nucleus, even though the interfacial free energy may be considerably reduced. One way to think about the effect of strain is that it shifts the coexistence line, e.g. the solvus, to lower temperatures.
image: 65_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_21_23_Fig3.svg
Figure 16.3: GX plots for a regular solution model with critical temperature of 1000K without strain (A) and in the presence of strain (B). The effect of strain is to lower the critical temperature, here to a value of 800K. Minima (open circle) and inflections points (x) are indicated on each curve. C. In this phase diagram, the binodal line was constructed from the minima of the GX plot in A (blue line) and B (red line) and the spinodal line was constructed from the inflection points (dashed blue and red lines). Note that while strain shifts the binodal and spinodal lines to smaller temperatures, the effect is strongest for X= 1 2 , but is attenuated as X increases or decreases. The overall effect is one of vertical compression.
In spinodal decomposition, the negative curvature of the Gibbs free energy allows mass transport against the concentration gradient. Thermal fluctuations in the composition can therefore grow in amplitude. There is no sharp interface between regions of higher and lower concentration, but creating a concentration gradient is associated with an increase in free energy that resists the phase transformation, similar to the effect of the interfacial free energy in nucleation. As a consequence, only concentration fluctuations with a wavelength greater than the critical wavelength, i.e. a sufficiently shallow gradient, will lower the free energy of the system. Strain shifts the free energy to more positive values and therefore resists the phase transformation. Similar to the case of nucleation, this strain can be thought to shift the coexistence line, i.e. the spinodal, to lower temperatures. In case of a regular solid solution, the shift of the spinodal can be expressed as a reduction of the critical temperature (Figure 3).
In the following we'll take a quick look at the time evolution of spinodal decomposition. Recall that we can write the chemical diffusivity as D=M 2 G V C 2 , (16.34) where M>0 is the mobility.
Let's rewrite Fick's first and second law using this definition J =-D C z =-M 2 G V C 2 C z =-M z G V C (16.35) C t =- J z =M 2 z 2 G V C (16.36)
To write out the Fick's second law, we therefore need to find G V C . Recall that we can write the free energy of the system with a fluctuation with wavelength lambda as G V = 1 λ 0 λ [ G V ( C ) + κ ( C z ) 2 + W V e l ( C ) ] z (16.37)
One can show that the partial derivative of the Gibbs free energy per unit volume with respect to concentration is G V C = . 2 G V C 2 | C 0 ( C- C ) -2 κ 2 C z 2 + s trainterm (16.38)
Plugging into Eq. (36), we find C t =M . 2 G V C 2 | C 0 2 C z 2 -2 κ 4 C z 4 + s trainterm (16.39)
Eq. (36) is known as Cahn's Diffusion Equation. The time evolution of spinodal decomposition is described by a solution to this differential equation:
C( z,t ) - C = A m ( λ ) e R( λ ) t cos 2 π z λ , (16.40) where A m ( λ ) is the amplitude of the fluctuation with wavelength λ at t=0 , and R( λ ) =-M ( 2 π λ ) 2 [ . 2 G V C 2 | C 0 +2 κ ( 2 π λ ) 2 + T s train ] , (16.41) where T s train is a positive term that reflects coherency strain.
Note that in (40) and (41) we treat all spatial frequencies λ . We can think of any initial concentration profile as a sum of cosines, i.e. a Fourier series. At t=0 , their amplitudes (Fourier coefficients) are described by A m ( λ ) . These amplitudes evolve over time, as described by e R( λ ) t . Let's take a closer look at R( λ ) .
To do so, let's focus on the important variable and collect everything else in positive constants: R( λ ) =- A 1 λ 2 [ . 2 G V C 2 | C 0 + T s train + A 2 λ 2 ] (16.42) =- A 1 λ 2 [ . 2 G V C 2 | C 0 + T s train ] - A 1 A 2 λ 4 , (16.43) where A 1 =M(2 π ) 2 and A 2 =2 κ (2 π ) 2 .
Inspection of Eq. (43) reveals that if the sum of the terms in parentheses is negative, then R is the sum of a positive term that is proportional to λ -2 and a negative term proportional to λ -4 (Figure). This immediately tells us that R( λ ) will change sign for some critical value of λ , that R( λ ) <0 for λ < λ C , and that R( λ ) >0 for λ > λ C . Furthermore, R has a maximum at some λ max > λ C and asymptotically approaches the x -axis for large λ .
To find λ C , we write R( λ C ) =0 (16.44) - A 1 λ C 2 [ . 2 G V C 2 | C 0 + T s train ] - A 1 A 2 λ C 4 =0 (16.45)
Because λ C >0 and A 1 >0 we can multiply both sides with with λ C 4 - A 1 λ C 2 [ . 2 G V C 2 | C 0 + T s train ] - A 1 A 2 =0 (16.46) λ C 2 = - A 2 . 2 G V C 2 | C 0 + T s train (16.47) λ C = π -8 κ . 2 G V C 2 | C 0 + T s train (16.48)
Not surprisingly, we find that the critical wavelength is the same as before. We can express R in terms of λ C to better understand its behavior R( λ ) =2 κ M ( 2 π λ ) 4 [ λ 2 λ C 2 -1 ] (16.49)
The significance of the sign of R( λ ) is that it affects how the amplitude of the fluctuation with wavelength λ evolves over time. We can differentiate three cases:
This means that the amplitude of fluctuations with a wavelength smaller than λ C will be dampened with time, that of fluctuations with a wavelength equal to λ C will remain the same, and that of fluctuations with a wavelength greater than λ C will be amplified. Note that the derivation given here is only valid for short times and small amplitudes. This can be seen from the fact that the the amplitude for the amplified wavelengths goes to infinity with time, which is unphysical. Also note that while Cahn and Hilliard assumed that the mobility is independent of the concentration, this is not true and becomes important as spinodal decomposition progresses (see for example Nauman and He, Chemical Engineering Science 2001, 56, 1999-2018).
Recall that after R( λ ) passes through its maximum it approaches zero for large λ , meaning that fluctuations with a wavelength near the maximum will be amplified the most and fluctuations with very long wavelengths will be amplified very little. We find λ max in straightforward fashion . dR( λ ) d λ | λ max =0 (16.50) λ max = 2 λ C (16.51)
As an example, let's assume that a binary solid solution is quenched into the spinodal regime and that at t=0 . Initially (column 1 in Fig. 4) periodic concentration fluctuations with λ λ C = 0.2, 0.5, 0.75, 1, ( 2) , and 2.5 (rows 1-6) shall all have the same amplitude. Adding these fluctuations up gives the concentration profile in the one dimensional sample (bottom row). With increasing time (columns 15 ), the amplitudes are modulated as described by Eq. (40). For fluctuations with wavelengths smaller than λ C , the amplitude is rapidly attenuated. This can be seen by comparing the amplitude of the fluctuations in each of rows 1 through 3 across the columns, from left to right. For λ = λ C , the amplitude stays the same (row 4), and for λ > λ C the amplitude increases. The fluctuation with λ = 2 λ C is amplified most rapidly (cf. row 5 vs row 6).
image: 66_home_ken_Mydocs_MSEcore_316-2_figures_Lecture_21_23_Figure4.svg
Figure 16.4: Evolution of the amplitude of 1D concentration fluctuations with time in spinodal decomposition. Note how rapid attenuation of low wavelength (high frequency) fluctuations creates smooth profiles. The Matlab code used to generate this plot is given ... A Matlab code to simulate spinodal decomposition in 2D is available on Canvas (CahnHilliard_DJ.m).

Case Study: Precipitation of Cu-rich Precipitates from Al:

Al is soft - can be precipitation hardened by adding Cu
Pure Al α phase ( κ in above drawing): FCC, a = 4.04Å, valence = 3+
Pure Cu: FCC, a = 3.62Å, valence =1+
Equilibrium precipitate when Al is supersaturated with Cu: θ - CuAl 2 - tetragonal
 
Al-rich region of phase diagram:
image: 67_home_ken_Mydocs_MSEcore_316-2_figures_drawing47.svg
Figure 17.1: Cu-Al phase diagram.
 
θ cannot form coherent interfaces with α , γ for this interface is high:
W r * γ 3 , v exp ( - W r * / k B T ) , nucleation rate of θ is very slow
Other non-equilibrium precipitates - Guinier-Preston (GP) zones, θ " , θ '
Al: FCC, a = 4.04Å
image: 68_home_ken_Mydocs_MSEcore_316-2_figures_drawing48.svg
Figure 17.2: Structure of a GP zone
GP Zones: FCC, totally coherent, very small disks: 2 atomic layers thick, 50 atoms in diameter
 
θ " : tetragonal, FCC distorted in one direction
image: 69_home_ken_Mydocs_MSEcore_316-2_figures_drawing49.svg
Figure 17.3: Structure of the θ " phase.
all interfaces coherent with α phase
 
θ ' : tegragonal, only one interface coherent with α phase:
image: 70_home_ken_Mydocs_MSEcore_316-2_figures_drawing50.svg
Figure 17.4: Structure of the θ ' phase.
 
θ phase: totally incoherent
image: 71_home_ken_Mydocs_MSEcore_316-2_figures_drawing51.svg
Figure 17.5: Structure of the θ phase.
 
Rank the values of γ , nucleation rate and thermodynamic stability:
precipitate
γ
v
G m
GP Zone
1
4
4
θ "
2
3
3
θ '
3
2
2
θ
4
1
1
Table 1: Order of the surface energy, nucleation rate and free energy for different Cu-rich precipitate phases in Al, ranked from lowest (1) to highest (4).
 
Most unstable phase appears first because of higher nucleation rate:
Map time dependence with TTT diagram:
 
General time dependence for a transformation:
image: 72_home_ken_Mydocs_MSEcore_316-2_figures_drawing52.svg
Figure 17.6: Generic TTT diagram.
 
Unstable compounds only appear at larger undercoolings - different solvus for each phase
image: 73_home_ken_Mydocs_MSEcore_316-2_figures_drawing53.svg
Figure 17.7: TTT diagram for precipitate formation in the Cu-Al system.
 
At a given undercooling, each phase is in eq'm with a different concentration:
image: 74_home_ken_Mydocs_MSEcore_316-2_figures_drawing54.svg
Figure 17.8: Schematic representation of the solvus lines for different precipitate phases in the Al-Cu system.
 
Each concentration is affected by a local curvature:
X r = X e ( 1+ 2 γ V m RTr ) (17.1)
 
Once an equilibrium phase is formed, it grows at the expense of a non-eq'm phase because of the lower concentration of Cu in eq'm with it:
image: 75_home_ken_Mydocs_MSEcore_316-2_figures_drawing55.svg
Figure 17.9: Cu flux from a more soluble to a less-soluble precipitate phase.
Matrix concentration of B decreases from X 0 to X GB to X θ " to X θ ' to X θ as different precipitates appear and grow to an appreciable size.
 
Types of precipitation:
  1. All nucleation occurs at the beginning: site saturation (heterogeneous nucleation sites used up initially)
    image: 76_home_ken_Mydocs_MSEcore_316-2_figures_drawing56.svg
    Figure 17.10: Site saturation.
  2. Constant nucleation rate throughout transformation:
    image: 77_home_ken_Mydocs_MSEcore_316-2_figures_drawing57.svg
    Figure 17.11: Constant nucleation rate throughout precipitation.
  3. Cellular precipitation: entire parent phase is consumed (recrystallization, formation of pearlite)
    image: 78_home_ken_Mydocs_MSEcore_316-2_figures_drawing58.svg
    Figure 17.12: Cellular precipitation.
 
What is time dependence of transformation?
f( t ) = fractiontransformedattimet finafractiontransformed (depends on temperature)
image: 79_home_ken_Mydocs_MSEcore_316-2_figures_drawing59.svg
Figure 17.13: Time dependence of transformation
Ex: derivation of f( t ) for constant v
 
Vol = volume of a single cell
v = growth velocity
r = cell radius
 
Vo( t ) = 4 3 π r 3 = 4 3 π (vt ) 3 (nucleation occurs at t=0 )
 
If nucleation does not occur until t= τ :
Vol( t, τ ) = 4 3 π v 3 (t- τ ) 3 (17.2)
v d τ = number of nuclei formed during time increment d τ .
f( t ) = v 0 t Vol( t, τ ) d τ = v 4 3 π v 3 0 t (t- τ ) 3 τ = v π 3 v 3 t 4 (17.3)
for long t , f( t ) =1 :
detailed solution gives:
f( t ) =1- exp ( - π 3 v v 3 t 4 ) Johnson-Mehl-Avrami equation
(in agreement with previous expression at small t , since 1- exp ( -x ) x for small x
This is a specific example of the following more general expression:
f( t ) =1- exp ( -k t n ) (17.4)
n is found to vary between 1 and 4
temperature dependence is in k :
This is the simplest equation that has the basic behavior observed experimentally:
image: 80_home_ken_Mydocs_MSEcore_316-2_figures_drawing60.svg
Figure 17.14: Behavior of the Johnson-Mehl Avrami equation
k is low at high temp. because N v is small (low undercooling)
k is low at low temp. because v is small (slow diffusion)

Case Study: Mineralization from Solution

Basic Carbonate forming reaction: (all concentrations in moles/liter)
CaC O 3 C a 2+ +C O 3 2- [ C a 2+ ] [ C O 3 2- ] =6x1 0 -9 (18.1)
[ C a 2+ ] [ C O 3 2- ] >6x1 0 -9  for aragonite, vaterite
C O 3 2- - concentration affected by pH
HC O 3 - H + +C O 3 2- [ H + ] [ C O 3 2- ] [ HC O 3 - ] =5.61x1 0 -11 (18.2)
H 2 C O 3 H + +HC O 3 - [ H + ] [ HC O 3 - ] [ H 2 C O 3 ] =1.5x1 0 -4 (18.3)
Dissolved C O 2 in equilibrium with carbonic acid
H 2 O+C O 2 H 2 C O 3 [ H 2 C O 3 ] [ C O 2 ] =1.70x1 0 -3 (18.4)
Dissolved C O 2 in equilibrium with atmospheric C O 2
dissolvedC O 2 atmosphericC O 2 P C O 2 [ C O 2 ] =29.76atm/mol/l (18.5)
Basic water equilibrium
H 2 O H + +O H - [ H + ] [ O H - ] =1 0 -14 (18.6)
Charge neutrality
2[ C a 2+ ] +[ H + ] =[ HC O 3 - ] +2[ C O 3 2- ] +[ O H - ] +[ C l - ] (18.7)
Other equilibria. If Chloride ions are present:
HCl H + +C l - [ H + ] [ C l - ] HCl =1 0 -8 (18.8)
So we have 8 equations and 10 unknowns:
[ C a 2+ ] , [ C O 3 2- ] , [ H + ] , [ HC O 3 - ] , [ H 2 C O 3 ] , [ C O 2 ] , PC O 2 , [ O H - ] , [ HC ] , [ C - ]
So we need to specify 2 quantities, generally [ C a 2+ ] and P CO2 to completely specify the problem. Note: pH=-og[ H + ]
Increasing [ C O 2 ] increases carbonic acid concentration (Eq.(4)), decreases pH , increases [ H + ]

Case Study: Nanoparticle Nucleation and Growth

image: 81_home_ken_Mydocs_MSEcore_316-2_figures_nanoparticle_pictures.png
Figure 19.1: Metallic Nanoparticles.[5]

19.1 Motivation

image: 82_home_ken_Mydocs_MSEcore_316-2_figures_nanoparticle_colors.png
Figure 19.2: Size dependent optical properties of metallic nanoparticles.

19.2 Wulff Constuction

image: 83_home_ken_Mydocs_MSEcore_316-2_figures_wulff_construction.png
Figure 19.3: The Wulff Construction.
  1. Create polar plot of surface energy (red)
    1. Find intersections with lines drawn from the origin (brown)
    2. Draw perpendicular planes (black)
    3. Eq'm shape defined by innermost planes (blue)

19.3 Kinetic Control of Particle Shape

image: 84_home_ken_Mydocs_MSEcore_316-2_figures_particle_shapes.png
Figure 19.4: Dependence of particle shape on the growth rate ratio.
image: 85_home_ken_Mydocs_MSEcore_316-2_figures_particle_shape_111_100.png
Figure 19.5: Crystal facets for a cubic material.

19.4 Controlling Growth Rates


image: 86_home_ken_Mydocs_MSEcore_316-2_figures_polyvinylpyrrolidone.png
Figure 19.6: Chemical structure of poly(vinyl pyrrolidone).
image: 87_home_ken_Mydocs_MSEcore_316-2_figures_silver_cubes.png
Figure 19.7: Cubic particles form rapid growth in 111 directions.
Nucleation Rate and Size Distribution
image: 88_home_ken_Mydocs_MSEcore_316-2_figures_drawing61.svg
Figure 19.8: Time dependent reactant concentration profiles leading to two types of particle size distributions.
image: 89_home_ken_Mydocs_MSEcore_316-2_figures_optical_growth_determination.png
Figure 19.9: Optical Determination of Nucleation and Growth.

20 Interface Stability

In this section we are concerned with factors that cause the interface between a solid and a liquid to propagate in a non planar fashion, resulting in a dendritic microstructure of the sort shown below in Figure. 20.1.
image: 90_home_ken_Mydocs_MSEcore_316-2_figures_voorhees_dendrite.png
Figure 20.1: Image of a growth highlighting the curvature.
image: 91_home_ken_Mydocs_MSEcore_316-2_figures_drawing38.svg
Figure 20.2: Solidification in a cold mold.
Heat flow governed by thermal conductivity: q=-k dT dx
q = heat flux, k = thermal conductivity
Growth into supercooled liquid:
image: 92_home_ken_Mydocs_MSEcore_316-2_figures_drawing13.svg
Figure 20.3: Interface instability due to solidification into a supercooled liquid.
Growth by removal of heat from the wall (more common)
image: 93_home_ken_Mydocs_MSEcore_316-2_figures_drawing81.svg
Figure 20.4: Solidification against a cold surface.
Alloys: similar, except formation of dendrites can be controlled by movement of atoms; dendrites can form even when T increases into liquid
Solidify a material with overall composition C 0 : higher impurity concentration in liquid phase
image: 94_home_ken_Mydocs_MSEcore_316-2_figures_drawing39.svg
Figure 20.5: Phase diagram for immiscible solid phases.
Composition profile:
image: 95_home_ken_Mydocs_MSEcore_316-2_figures_drawing14.svg
Figure 20.6: Impurity concentration profile at a solid/liquid interface.
V = velocity of solid/liquid interface
D = diffusion of coefficient of impruity in liquid
How does this impurity buildup affect stability of the interface (tendency to form dendrites)
Use concentration profile to define a liquidus temperature that depends on position:
image: 96_home_ken_Mydocs_MSEcore_316-2_figures_drawing15.svg
Figure 20.7: Significance of local liquid temperature.
image: 97_home_ken_Mydocs_MSEcore_316-2_figures_Lect3_Fig7.png
Figure 20.8: Mapping local concentrations to local liquidus temperature.
So we can convert concentration profile to liquidus temperature as a function of distance:
image: 98_home_ken_Mydocs_MSEcore_316-2_figures_drawing16.svg
Figure 20.9: Plot of the local liquidus temperature.
If the temperature is less than the local T , solid will form, if T> T , no solid forms
Stability criterion determined by comparing actual temperature profiles
image: 99_home_ken_Mydocs_MSEcore_316-2_figures_drawing17.svg
Figure 20.10: Graphical criterion for interface stability.
Case 1: Actual T always above T = stable interface
Case 2: Actual T below T , solid protrusion can grow:
image: 100_home_ken_Mydocs_MSEcore_316-2_figures_drawing18.svg
Figure 20.11: Growth of an unstable interface.
Note: interface is always stable for a sufficiently large temperature gradient into the sample.

Eutectic Solidification

Thermodynamics review - melting point depression
Thermodynamics: Equation for liquidus line
X β ( T ) = liquidus comp. for α side of diagram
X S β ( T ) = solidus comp. for α side of diagram
G B = partial molar free energy of liquid B
G B S = partial molar free energy of solid B
T M A , T M B - melting points of pure components
 
At eq'm:
G B L ( X L β ( T ) ) = G B S ( X S β ( T ) ) (21.1)
If X S β ( T ) 1 (nearly pure solid phase):
G B S ( X S β ( T ) ) = G B L X S β ( T ) - Δ H B ( 1-T/ T m B ) (21.2)
Combine (1) and (2):
G B ( X β ) - G B ( X S β ) =- Δ H B ( 1- T T m B ) (21.3)
Assume liquids form an ideal solution:
G B L ( X L β ) - G B L ( X S β ) =RT ln X L β (21.4)
Combine (3) and (4):
RT ln X β =- Δ H B ( 1- T T m B ) (21.5)
Melting point reduction due to increased entropy of liquid phase
Assumptions:
  1. Ideal mixing in liquid (no enthalpy of mixing)
  2. Solid phases are nearly pure
Optimum phase size for eutectic microstructure:
(draw alternating α , β , show wavelength, λ )
Overall enthalpy of melting of solid at composition X e :
Δ H e = X B Δ H B +( 1- X B ) Δ H A (weighted average of component enthalpies)
 
Two contributions to free energy of solidification:
  1. Bulk term:
    Δ G buk V = - Δ H e Δ T V m T e , Δ T = undercooling ( T e -T )
  2. Interfacial term:
    within one wavelength - total interfacial area = 2A ( A = cross section) - 2 α / β interfaces per wavelength)
    volume = λ A
    Δ G int V = 2 γ α β λ
    Δ G tota V = 2 γ α β λ - Δ H e Δ T V m T e
    Δ G tota <0 for λ > λ * , λ * = 2 γ α γ V m T e Δ H Δ T
    compare to r * for nucleation of a solid phase: r * = -2 γ ( Δ G buk /V ) = 2 γ V m T e Δ H Δ T
Alternate Approach: growth of curved lamellae
image: 101_home_ken_Mydocs_MSEcore_316-2_figures_drawing62.svg
Figure 21.1:
R α R β λ , details depend on relative values of γ α β , γ a , γ β .
Supercooled eutectic alloy (composition X B = X E ) - consider α phase
image: 102_home_ken_Mydocs_MSEcore_316-2_figures_drawing63.svg
Figure 21.2:
Conc. gradient drives B away from α / interface:
image: 103_home_ken_Mydocs_MSEcore_316-2_figures_drawing64.svg
Figure 21.3:
Now consider β phase
Supercooled eutectic alloy (composition X B = X E ) - consider α phase
image: 104_home_ken_Mydocs_MSEcore_316-2_figures_drawing65.svg
Figure 21.4:
Conc. gradient drives B toward from β / interface:
image: 105_home_ken_Mydocs_MSEcore_316-2_figures_drawing66.svg
Figure 21.5:
Simultaneous growth of α and β minimizes required diffusion distance:
In liquid: B is enhanced near α , depleted near β .
image: 106_home_ken_Mydocs_MSEcore_316-2_figures_drawing67.svg
Figure 21.6:
Concentration gradient in liquid: Δ X λ λ
Diffusive flux in liquid: D Δ X λ
What determines Δ X λ ?
Liquidus lines for finite λ are suppressed in comparison to eq'm ( λ = ) values.
Eutectic point of phase diagram:
image: 107_home_ken_Mydocs_MSEcore_316-2_figures_drawing68.svg
Figure 21.7:
  1. Assume liquidus lines are linear in this regime - Must have Δ T 0 linearly related to Δ T λ .
  2. Must have Δ X λ =0 for λ = λ *
  3. Must have Δ X λ = Δ X 0 for λ =
  4. λ * 1 Δ T 0
All this gives: Δ X λ = Δ X 0 ( 1- λ * λ )
If growth is diffusion limited:
V Δ X 0 D λ ( 1- λ * λ ) (21.6)
Maximum growth velocity for λ =2 λ * .

Eutectoid Reactions

Eutectic transformation:
V Δ X 0 D λ ( 1- λ * λ )
Δ X 0 Δ T 0
Maximum V for λ =2 λ * .
image: 108_home_ken_Mydocs_MSEcore_316-2_figures_drawing69.svg
Figure 22.1:
Eutectoid transformation: High temp. phase is a solid
image: 109_home_ken_Mydocs_MSEcore_316-2_figures_drawing70.svg
Figure 22.2:
γ α +F e 3 C
Lamellar microstructure of α , F e 3 C = pearlite
From lever rule: mostly ferrite
Nucleation occurs preferentially at γ grain boundaries:
image: 110_home_ken_Mydocs_MSEcore_316-2_figures_drawing71.svg
Figure 22.3:
Nucleation:
  1. F e 3 C nucleates first
    image: 111_home_ken_Mydocs_MSEcore_316-2_figures_drawing72.svg
    Figure 22.4:
  2. region near F e 3 C ppt. depleted in C
  3. α nucleates
    image: 112_home_ken_Mydocs_MSEcore_316-2_figures_drawing73.svg
    Figure 22.5:
    Incoherent interfaces of α / γ 2 and F e 3 C/ γ 2
    Coherent interfaces for α / γ 1 and F e 2 C/ γ 1 .
    F e 3 C is orthorhombic: abc , α = β = γ =90 º - γ is FCC
Preferred orientation for coherent F e 3 C γ interfaces:
(100 ) F e 3 C (111 ) γ ; (010 ) F e 3 C (110 ) γ ; (001 ) F e 3 C (112 ) γ
Growth rate limited by C diffusion in γ phase (like eutectic case).
image: 113_home_ken_Mydocs_MSEcore_316-2_figures_drawing74.svg
Figure 22.6:
Data: Growth rate (Hull, 1942)
image: 114_home_ken_Mydocs_MSEcore_316-2_figures_drawing75.svg
Figure 22.7:
Wavelength (Pellissier, 1942)
image: 97_home_ken_Mydocs_MSEcore_316-2_figures_Lect3_Fig7.png
Figure 22.8:

Review Questions

Matlab Example Scripts

24.1 Free energy of water droplets formed from water vapor at 10 K (Figure 4.2)

%% Lecture 2 Figure 2
% free energy of water droplets formed from water vapor at 10 K
% supercooling, as a function of the radius
% all data from Wikipedia

close all;
clear all;

Mw          = 18;% [g/mol]
DeltaS0_f   = 22;% [J/mol/K]
T_m         = 273.15;% [K]
DeltaS0_vap = 118.89;% [J/mol/K]
T_b         = 373.15;% [K]

rho_363     = 0.96531;% [g/cm^3]
Vm_363      = Mw/rho_363*1E-6;% [m^3/mol] 

gamma       = 0.072;% [J/m^2] this is for the liquid-vapor interface

P_vap       = 10^5;% [Pa]
DeltaT      = 10;  % [K] supercooling

% calculate the Laplace pressure
DeltaP = DeltaS0_vap/Vm_363*DeltaT;

% calculate the critical radius
Rstar = 2*gamma/DeltaP;% [m]

% calculate the reversivle work of nucleation
Wstar = 16/3*pi*gamma^3/DeltaP^2; %[J]

% define radii for plot of the free energy of a small spherical particle of the new
% phase
R = logspace(-10,-8,100);
% calculate free energy difference (Delta Omega) at the radii
omega = -DeltaP * 4/3*pi*R.^3 + gamma * 4*pi*R.^2;% [J]

% plot
figure;
% plot omega
plot(R*1E9,omega,'b'); hold on;
% plot bulk term only
plot(R*1E9,-DeltaP * 4/3*pi*R.^3,'r');
% plot excess (surface) term only
plot(R*1E9,gamma * 4*pi*R.^2,'g');
% draw horizontal line at 
plot(R*1E9,zeros(size(R)),'k--');
% draw vertical line at critical radius
plot([Rstar*1E9, Rstar*1E9],[Wstar 0],'k:');
% draw horizontal line at maximum
plot([0, Rstar*1E9],[Wstar Wstar],'k:');


set(gca,'TickDir','out');
xlim([10^-1,4]);
ylim([-.5*10^-17,.5*10^-17]);
xlabel('\rightarrow R [nm]');
ylabel('\Delta\Omega(R) [J] for one nucleus');

24.2 Laplace pressure, critical radius, and reversible work of nucleation for condensation of water from vapor, as a function of temperature (Figure 4.4)

%% Lecture 2 Figure 2
% Laplace pressure, critical radius, and reversible work of nucleation for
% condensation of water from vapor, as a function of temperature
% all data from Wikipedia

close all;
clear all;

k           = 1.38E-23; %J/K Boltzmann constant
Mw          = 18;% [g/mol]
DeltaS0_f   = 22;% [J/mol/K]
T_m         = 273.15;% [K]
DeltaS0_vap = 118.89;% [J/mol/K]
T_b         = 373.15;% [K]
rho_373     = 0.95805;% [g/cm^3]
Vm_373      = Mw/rho_373*1E-6;% [m^3/mol] 
gamma       = 0.072;% [J/m^2] this is for the liquid-vapor interface
P_vap       = 10^5;% [Pa]


DeltaP = @(DT) DeltaS0_vap/Vm_373*DT;
Rstar = @(DT) 2*gamma./DeltaP(DT);% [m]
Wstar = @(DT) 16/3*pi*gamma^3./DeltaP(DT).^2; %[J]

DeltaT      = [0.5:0.1:30];% [K] supercooling
T           = T_b-DeltaT;

figure;
subplot(1,4,1);
plot(DeltaT,DeltaP(DeltaT)/1E6,'b'); hold on;
set(gca,'TickDir','out');
axis square;
% xlim([10^-1,4]);
% ylim([-.5*10^-17,.5*10^-17]);
xlabel('\rightarrow \DeltaT [K]');
ylabel('\DeltaP [MPa]');

subplot(1,4,2);
plot(DeltaT,Rstar(DeltaT)*1e9,'b'); hold on;
set(gca,'TickDir','out');
axis square;

% xlim([10^-1,4]);
% ylim([-.5*10^-17,.5*10^-17]);
xlabel('\rightarrow \DeltaT [K]');
ylabel('R^* [nm]');

subplot(1,4,3);
plot(DeltaT,Wstar(DeltaT),'b'); hold on;
set(gca,'TickDir','out');
axis square;

% xlim([10^-1,4]);
% ylim([-.5*10^-17,.5*10^-17]);
xlabel('\rightarrow \DeltaT [K]');
ylabel('W^*_r [J]');

subplot(1,4,4);
semilogy(DeltaT,Wstar(DeltaT)/k./T,'b'); hold on;
set(gca,'TickDir','out');
axis square;
xlabel('\rightarrow \DeltaT [K]');
ylabel('W^*_r [kT]');

24.3 General Nucleation Rate (Figure 5.2)

%% Lecture 4 Figure 2 how to
close all
clear all
clc

Teq = 373;

f1 = @(DT) DT.^2
f2 = @(DT) Teq-DT;

DT = linspace(0,Teq,1000);

T1n = f1(DT)/max(f1(DT));
T2n = f2(DT)/max(f2(DT));
figure;

subplot(131);
plot(DT, T1n,'b',DT, T2n,'g',DT,T1n.*T2n,'r');
set(gca,'TickDir','out');
xlim([0,Teq]);
T3 = -1./(T1n.*T2n);

subplot(132);
plot(DT,T3,'r');
set(gca,'TickDir','out');

xlim([0,Teq]);
ylim([-1000,0]);

subplot(133);
plot(DT,exp(T3),'r');
set(gca,'TickDir','out');

xlim([0,Teq]);

%% account for diffusive term
figure;

max_y3=max(exp(T3));
y4=exp(T3).*exp(-200./f2(DT));
max_y4=max(y4);
y4_n = y4/max(y4)*max_y3;
plot(DT,exp(T3),'b');hold on;

%plot(DT,exp(-200./f2(DT)),'g');
plot(DT,y4_n,'r');

set(gca,'TickDir','out');

24.4 Nucleation Rate for Water Condensation (Figure 5.3)

%% Lecture 4 Figure 3
% Nucleation current for condensation of water from vapor, as a function of temperature
% all data from Wikipedia

% Derk Joester
% 9/23/2017

close all;
clear all;

k           = 1.38E-23;        % [J/K]     Boltzmann constant
Mw          = 18;              % [g/mol]   molecular weight of water
DeltaS0_f   = 22;              % [J/mol/K] standard molar entropy of fusion
T_m         = 273.15;          % [K]       melting point
DeltaS0_vap = 118.89;          % [J/mol/K] standard molar entropy of vaporization
T_b         = 373.15;          % [K]       boiling point
rho_373     = 0.95835;         % [g/cm^3]  density of liquid water at T_m
Vm_373      = Mw/rho_373*1E-6; % [m^3/mol] molar volume of liquid water at T_m
gamma       = 0.072;           % [J/m^2]   surface energy for water liquid-vapor interface
P_vap       = 10^5;            % [Pa]      atmospheric pressure
A           = 10^33;           % [s^-1 cm^-3] kinetic prefactor

DeltaP = @(DT) DeltaS0_vap/Vm_373*DT;         % [Pa]
Rstar = @(DT) 2*gamma./DeltaP(DT);            % [m]
Wstar = @(DT) 16/3*pi*gamma^3./DeltaP(DT).^2; % [J]
NdotV = @(DT) A*exp(-Wstar(DT)/k./(T_b-DT));  % [s^-1 cm^-3]

DeltaT      = [1:0.1:T_b];% [K] supercooling 
T           = T_b-DeltaT; % [K] actual temperature

% find DeltaT at which NdotV(x)==1, choose positive solution closest to T_b
syms x
S_DT = vpasolve(NdotV(x)==1,x,[1 50]);

figure;
subplot(2,2,1);

yyaxis left
plot(NdotV(DeltaT),DeltaT,'b');
set(gca,'YDir','reverse');
set(gca,'TickDir','out');
axis square;
ylim([0,T_b]);
ylabel('\leftarrow \DeltaT [K]');
xlabel('N_V [s^-^1 cm^-^3]');

yyaxis right
plot(NdotV(DeltaT),T,'b');
ylim([0,T_b]);
ylabel('\rightarrow T [K]');

subplot(2,2,2);
yyaxis left
semilogx(NdotV(DeltaT),DeltaT,'b');hold on;
semilogx([1,1],DeltaT([1,end]),'k--');
set(gca,'YDir','reverse');
set(gca,'TickDir','out');
set(gca,'XTick',logspace(-50,50,11));
axis square;
ylim([0,T_b]);
xlim([1e-50,1e+50]);
ylabel('\leftarrow \DeltaT [K]');
xlabel('N_V [s^-^1 cm^-^3]');

yyaxis right
semilogx(NdotV(DeltaT),T,'b');
ylim([0,T_b]);
xlim([1e-50,1e+50]);
ylabel('\rightarrow T [K]');
set(gca,'XTick',logspace(-50,50,11));


subplot(2,2,3);

yyaxis left
plot(NdotV(DeltaT),DeltaT,'b');
set(gca,'YDir','reverse');
set(gca,'TickDir','out');
axis square;
ylim([5,30]);
ylabel('\leftarrow \DeltaT [K]');
xlabel('N_V [s^-^1 cm^-^3]');

yyaxis right
plot(NdotV(DeltaT),T,'b');
ylim([T_b-30,T_b-5]);
ylabel('\rightarrow T [K]');

subplot(2,2,4);
yyaxis left
semilogx(NdotV(DeltaT),DeltaT,'b');hold on;
plot([1,1],[S_DT,30],'k--');
plot([1e-50,1],[S_DT,S_DT],'k--');

set(gca,'YDir','reverse');
set(gca,'TickDir','out');
set(gca,'XTick',logspace(-50,50,11));
axis square;
ylim([5,30]);
xlim([1e-50,1e+50]);
ylabel('\leftarrow \DeltaT [K]');
xlabel('N_V [s^-^1 cm^-^3]');

yyaxis right
semilogx(NdotV(DeltaT),T,'b');
ylim([T_b-30,T_b-5]);
xlim([1e-50,1e+50]);
ylabel('\rightarrow T [K]');
set(gca,'XTick',logspace(-50,50,11));

24.5 Figure TOC entry

%% Lecture 5 Figure 3
% Plots of the structure factor S(theta) for spherical cap geometry

close all
clear all

S = @(x) 0.25*(2+cos(x)).*(1-cos(x)).^2;
theta = 0:0.1:180;

figure;

subplot(121);
plot(theta,S(deg2rad(theta)),'b');
set(gca,'TickDir','out');
xlabel('\theta [?]');
ylabel('S(\theta)');
xlim([0,180]);
ylim([1e-10,1]);
grid on

subplot(122);
semilogy(theta,S(deg2rad(theta)),'b');
set(gca,'TickDir','out');
xlabel('\theta [?]');
ylabel('S(\theta)');
xlim([0,180]);
ylim([1e-10,1]);
grid on

%% some solutions
syms x
Sv = [1e-5, 1e-3, 1e-2, 1e-1, .5];

for i=1:length(Sv)
    thv(i) = rad2deg(double((vpasolve(S(x)==Sv(i),x,[0 pi]))));
end

24.6 Figure TOC entry

% Draw spherical caps 
clear all
close all

S = @(x) 0.25*(2+cos(x)).*(1-cos(x)).^2;
syms x
Sv = [1e-5, 1e-3, 1e-2, 1e-1, .5, 0.9, 0.99, 0.999];

for i=1:length(Sv)
    thv(i) = rad2deg(double((vpasolve(S(x)==Sv(i),x,[0 pi]))));
end;

h = @(theta,R) R*(1-cosd(theta));
a = @(theta,R) sqrt(h(theta,R).*(2*R-h(theta,R)));

Rstar = 1;

hv = h(thv,Rstar);
av = a(thv,Rstar);

center = [zeros(size(hv))',hv'-1];

for i=1:length(hv)
    subplot(3,3,i);
    hold on;
    clear intLbeta rotanv

    r1=[-av(i);0;0]-[center(i,:)';0]
    r2=[av(i);0;0]-[center(i,:)';0]

    rotvec = vrrotvec(r1,r2)
    if rotvec(3)==1
        rotvec(4)=2*pi-rotvec(4);
        rotvec(3)=-1;
    end;
    
    rotanv = linspace(0,rotvec(4),100);
    
    for j=1:length(rotanv)
      M=vrrotvec2mat([rotvec(1:3),rotanv(j)]);
      intLbeta(j,:)=M*r1+[center(i,:)';0];
    end;
    
    plot(intLbeta(:,1),intLbeta(:,2),'r');hold on
    plot([-3,3],[0,0],'k');
    plot([-3,-av(i)],[0,0],'b');
    plot([av(i),3],[0,0],'b');
    plot(center(i,1),center(i,2),'k+');

    if thv(i)<90
        xlim([-1.5*av(i),1.5*av(i)]);
        ylim([-1.5*av(i),1.5*av(i)]);
    else
        xlim([-2,2]);
        ylim([-2,2]);
    end;
    set(gca,'TickDir','out');
    axis equal
    titlestr = ['\theta = ',num2str(thv(i),'%3.1f'),'?; S(\theta) = ',num2str(Sv(i),'%3.2e')];
    title(titlestr);
    xlabel('\rightarrow x/R^*');
    ylabel('\rightarrow z/R^*');
end;

24.7 Figure TOC entry

% Lecture 8 Figure 2
close all
clear all
clc

k = 1.38E-23; % J/mol/K
T = 473; % K

gamma_incoherent = 0.5; % J/m^2
gamma_coherent   = 0.2; % J/m^2

C                = 5e10; % J/m^3 = Pa

W_V_el           = 2e9;  % J/m^3 = Pa 

DX1 = linspace(0,0.2,100);
DT1 = 100*DX1;

DX2 = linspace(W_V_el/C,0.2,100);
DT2 = 100*DX2;

W_i = 16*pi/3 * gamma_incoherent^3./(C*DX1).^2/k./(T-DT1);
W_c = 16*pi/3 * gamma_coherent^3./(C*DX2-W_V_el).^2/k./(T-DT2);

% figure;
% semilogy(DX1,W_i,'b',DX2,W_c,'r');
% set(gca,'TickDir','out');
% xlabel('\rightarrow \Delta X');
% ylabel('\rightarrow W^*_r [kT]');

figure;
plot(DX1,W_i,'b',DX2,W_c,'r');hold on;
ylim([0,200]);
plot([W_V_el/C,W_V_el/C],[0, 200],'k:');
set(gca,'TickDir','out');
xlabel('\rightarrow \Delta X');
ylabel('\rightarrow W^*_r [kT]');
legend('incoherent','coherent');
% figure;
% semilogy(DX1,W_i.*k.*(T-DT1),'b',DX2,W_c.*k.*(T-DT2),'r');
% set(gca,'TickDir','out');
% xlabel('\rightarrow \Delta X');
% ylabel('\rightarrow W^*_r [J]');

figure;
plot(DX1,W_i.*k.*(T-DT1),'b',DX2,W_c.*k.*(T-DT2),'r');hold on;
plot([W_V_el/C,W_V_el/C],[0, 5e-18],'k:');
ylim([0,5e-18]);
set(gca,'TickDir','out');
xlabel('\rightarrow \Delta X');
legend('incoherent','coherent');

24.8 Figure TOC entry

% Elastic Strain Energy Density

mu_alpha = 80E9; % Pa
K_beta   = 170E9;% Pa

Vm_alpha = 7.09E-6; % m^3/mol
Vm_beta  = linspace(Vm_alpha,Vm_alpha*1.2,100);

epsilon = 1/3*(Vm_beta-Vm_alpha)/Vm_alpha;

W_V_el = 18*epsilon.^2.*mu_alpha.*K_beta./(4*mu_alpha + 3*K_beta);

figure;
plot(epsilon, W_V_el/1E6);
xlabel('\epsilon');
ylabel('W_V^e^l [MPa]');

set(gca,'TickDir','out');

24.9 Figure TOC entry

%% Lecture 10 Figure 1
close all
clear all
clc

Teq = 1000; % [K] equilibrium temperature for phase transformation
DT = linspace(0,Teq,1000);

F1  = exp(-1e8*DT.^-2.*(Teq-DT).^-1);
F1n = F1/max(F1);

F2  = exp(-1e2*(Teq-DT).^-1);
F2n = F2/max(F2);

figure;
plot(DT,F1n,'b--');hold on
plot(DT,F2n,'r');hold on
yyaxis right
plot(DT,F1n.*F2,'b');hold on
set(gca,'TickDir','out');

%

F1  = @(Teq,DT) exp(-1e8*DT.^-2.*(Teq-DT).^-1);
F2  = @(Teq,DT) exp(-1e2*(Teq-DT).^-1);
Teq1 = 1000; % [K] equilibrium temperature for phase transformation at initial comosition X1
Teq2 = 900;  % [K] equilibrium temperature for phase transformation at initial comosition X2

DT1  = linspace(0,Teq1,1000);
T1   = Teq1-DT1;
DT2  = linspace(0,Teq2,1000);
T2   = Teq2-DT2;

N1 = F1(Teq1,DT1).*F2(Teq1,DT1);
N2 = F1(Teq1,DT1).*F2(Teq2,DT2);


figure;
plot(N1,T1,'b',N2,T2,'g');hold on;
ax = gca;
plot(ax,[0 ax.XLim(2)],[Teq1 Teq1],'b--');
plot(ax,[0 ax.XLim(2)],[Teq2 Teq2],'g--');
ax.TickDir = 'out';
ax.YLim(2) = Teq1+100;

%
figure;
plot(F1(Teq1,DT1), DT1, 'b--', N1,DT1,'b',F1(Teq2,DT2), DT2, 'g--', N2,DT2,'g', F2(Teq1,DT1), DT1,'r--');hold on;
ax = gca;
ax.TickDir = 'out';
ax.YDir = 'reverse';

24.10 Figure TOC entry

% shape of nucleus of a precipitate at a totally incoherent grain boundary,
% grain edge, and grain corner after Clemm and Fisher 1955
clear all
close all
clc
% case 1: planar grain boundary

gamma_aa =0.1;
gamma_ab =0.08;

theta = acos(gamma_aa/(2*gamma_ab));

azv = linspace(0,2*pi,36);
elv = linspace(pi/2-theta,pi/2,ceil((pi-theta)*18/pi));
r   = 1;

[AZ EL R] = meshgrid(azv,elv,r);
[X Y Z] = sph2cart(AZ, EL, R);
Z = Z-min(min(Z));

sphcap1 = surf2patch(X,Y,Z);
sphcap2 = surf2patch(X,Y,-Z);

% xy plane
xy.vertices = 1.25*[-r -r 0;-r r 0; r r 0; r -r 0];
xy.faces    = [1 2 3 4];

figure;
patch(sphcap1,'FaceColor','c','EdgeColor','k','AmbientStrength',0.5); hold  on;
patch(sphcap2,'FaceColor','c','EdgeColor','k','AmbientStrength',0.5);
patch('Faces',xy.faces,'Vertices',xy.vertices,'FaceColor','k','EdgeColor','none','FaceAlpha',0.2);

lighting gouraud
light('Position',[-2 0 4],'Style','local');
%light('Position',[0 0 -1],'Style','infinite')
%camlight(-60,90)
view(25,20);
axis equal
%axis off
grid on
xlim(1.25*[-r r]);
ylim(1.25*[-r r]);


24.11 Figure TOC entry

%% Lecture 10 Figure 1
% Plots of structure factors at a free interface, grain boundary, grain
% edge, and grain corner
% adapted from Clemm and Fisher, Acta Metall. 1955, 
close all
clear all

S_free = @(x) 0.25*(2+cos(x)).*(1-cos(x)).^2;

% Clemm and Fisher 
% derive a general form for the reversible work of nucleation for nucleus
% beta at the interface between two or more grains of alpha phase
% The ratio of this rev. work for the heterogeneous case over that for the
% homogeneous case is the structure factor with the general form

S = @(x,a,b,c) ((b-2*a.*cos(x)).^3)./(36*pi*c.^2);

a_gb = @(x) pi*(1-cos(x).^2);
b_gb = @(x) 4*pi*(1-cos(x));
c_gb = @(x) (2*pi/3)*(2+cos(x)).*(1-cos(x)).^2;

a_ge = @(x,beta) 3.*beta.*(1-cos(x).^2)-cos(x).*(3-4*cos(x).^2).^0.5;
b_ge = @(x,alpha, beta) 12*(pi/2-alpha-cos(x).*beta);
c_ge = @(x,alpha, beta) 2*(pi-2*alpha+cos(x).^2/3.*(3-4*cos(x).^2).^0.5-beta.*cos(x).*(3-cos(x).^2));
alpha = @(x) asin(1./(2*(1-cos(x).^2).^.5));
beta  = @(x) acos(cos(x)./(3*(1-cos(x).^2)).^.5);

a_gc  = @(x,phi,K) 3*(2*phi.*(1-cos(x).^2)-K.*((1-cos(x).^2-K.^2/4).^.5-K.^2/sqrt(8)));
b_gc  = @(x,phi,delta) 24*(pi/3-cos(x).*phi-delta);
c_gc  = @(x,phi,K,delta) 2*(4*(pi/3-delta)+cos(x).*K.*((1-cos(x).^2-K.^2/4).^.5-K.^2/sqrt(8))-2*cos(x).*phi.*(3-cos(x).^2));
K     = @(x) 4/3*(3/2-2*cos(x).^2).^.5-2/3*cos(x);
phi   = @(x,K) asin(K./(2*(1-cos(x).^2).^.5))
delta = @(x,K) acos((sqrt(2)-cos(x).*(3-K.^2).^.5)./(K.*(1-cos(x).^2).^.5));

theta_free = deg2rad(0:0.1:180);

theta_gb = deg2rad(0:0.1:90);
S_gb = S(theta_gb,a_gb(theta_gb),b_gb(theta_gb),c_gb(theta_gb));

theta_ge = deg2rad(rad2deg(acos(sqrt(3)/2)):0.1:90);
S_ge = S(theta_ge,a_ge(theta_ge,beta(theta_ge)),b_ge(theta_ge,alpha(theta_ge),beta(theta_ge)),c_ge(theta_ge,alpha(theta_ge),beta(theta_ge)));

theta_gc = deg2rad(rad2deg(acos(0.817)):0.01:90);
S_gc = S(theta_gc,a_gc(theta_gc,phi(theta_gc,K(theta_gc)),K(theta_gc)),b_gc(theta_gc,phi(theta_gc,K(theta_gc)),delta(theta_gc,K(theta_gc))),c_gc(theta_gc,phi(theta_gc,K(theta_gc)),K(theta_gc),delta(theta_gc,K(theta_gc))));
%
figure;

subplot(121);
plot(rad2deg(theta_gc),S_gc,'b',rad2deg(theta_ge),S_ge,'g',rad2deg(theta_gb),S_gb,'r');hold on
%plot(rad2deg(theta_free),S_free(theta_free),'k')
set(gca,'TickDir','out');
xlabel('\theta [\circ]');
ylabel('S(\theta)');
%xlim([0,90]);
ylim([0,1]);
grid on
legend('S_g_c(\theta)','S_g_e(\theta)','S_g_b(\theta)');

subplot(122);
semilogy(rad2deg(theta_gc),S_gc,'b',rad2deg(theta_ge),S_ge,'g',rad2deg(theta_gb),S_gb,'r');hold on
%semilogy(rad2deg(theta_free),S_free(theta_free),'k')
set(gca,'TickDir','out');
xlabel('\theta [\circ]');
ylabel('S(\theta)');
xlim([0,90]);
ylim([1e-6,1]);
grid on
legend('S_g_c(\theta)','S_g_e(\theta)','S_g_b(\theta)');

24.12 Figure TOC entry

%% grain corner precipitates

% this seems to work pretty well at larger d, i.e. displacement of the
% sphere along the axis, but runs into trouble with the surface mesh at
% smaller values. The code could probably be simplified by making better
% use of symmetry properties.

clear all
close all
clc

% Origin
O = [0,0,0];
% vertices of a tetrahedron
T1 = [1,0,-1/sqrt(sym(2))];
T1 = T1/norm(T1);
T2 = [-1,0,-1/sqrt(sym(2))];
T2 = T2/norm(T2);
T3 = [0,1,1/sqrt(sym(2))];
T3 = T3/norm(T3);
T4 = [0,-1,1/sqrt(sym(2))];
T4 = T4/norm(T4);

% normal vectors to planes

nP24 = cross(T2,T4)/norm(cross(T2,T4));
nP23 = cross(T2,T3)/norm(cross(T2,T3));
nP34 = cross(T3,T4)/norm(cross(T3,T4));

% 
T.vertices = double([T1;T2;T3;T4]);
T.faces    = [[1,2,3];[1 2 4];[1 3 4];[2 3 4]];



% define sphere S1, offset in the positive T1 direction
d=0.6;R=1;
C1 = d*R*T1;

% circle of intersection between S1 and the 3 planes OT2T3, OT3T4, and OT2T4
d23 = dot(C1-T2,nP23);
C23 = (C1-d23*nP23)';
r23 = sqrt(R^2-d23^2);

syms x
assume(x>0);
S = solve(norm(T2*x-C1)-R);

I2 = (T2*S)';
% rotation by 2pi/3 around T1 gives I3
I3 = vrrotvec2mat(double([T1,2*pi/3]))*I2;
% rotation by -2pi/3 around T1 gives I4
I4 = vrrotvec2mat(double([T1,-2*pi/3]))*I2;
% rotation by 2pi/3 around T4 gives I1
I1 = vrrotvec2mat(double([T4,2*pi/3]))*I2;

%plot3(C1(1),C1(2),C1(3),'rx'); % center of sphere
%plot3(C23(1),C23(2),C23(3),'gx');  % center of circle of intersection
% plot3(I2(1),I2(2),I2(3),'go'); % intercept of sphere with T2 axis
% plot3(I3(1),I3(2),I3(3),'bo'); % intercept of sphere with T3 axis
% plot3(I4(1),I4(2),I4(3),'mo'); % intercept of sphere with T4 axis
% plot3(I1(1),I1(2),I1(3),'ro'); % intercept of sphere with T1 axis

% determine points on arc from I2 to I3
alpha = acos(double(dot(I2-C23,I3-C23)/r23^2)); 
phi = linspace(0,alpha,20);
for i=1:length(phi)
    arc23{i}=double(C23+vrrotvec2mat(double([nP23,phi(i)]))*(I2-C23));
    arc34{i}=vrrotvec2mat(double([T1,2*pi/3]))*arc23{i};
    arc24{i}=vrrotvec2mat(double([T1,-2*pi/3]))*arc23{i};
    arc12{i}=vrrotvec2mat(double([T2,2*pi/3]))*arc23{i};
    arc13{i}=vrrotvec2mat(double([T1,2*pi/3]))*arc12{i};
    arc14{i}=vrrotvec2mat(double([T1,-2*pi/3]))*arc12{i};
end;
arc23mat = cell2mat(arc23);
arc34mat = cell2mat(arc34);
arc24mat = cell2mat(arc24);
arc12mat = cell2mat(arc12);
arc13mat = cell2mat(arc13);
arc14mat = cell2mat(arc14);
plot3(arc23mat(1,:),arc23mat(2,:),arc23mat(3,:),'g');hold on;
plot3(arc34mat(1,:),arc34mat(2,:),arc34mat(3,:),'b');
plot3(arc24mat(1,:),arc24mat(2,:),arc24mat(3,:),'m');
plot3(arc12mat(1,:),arc12mat(2,:),arc12mat(3,:),'r');
plot3(arc13mat(1,:),arc13mat(2,:),arc13mat(3,:),'c');
plot3(arc14mat(1,:),arc14mat(2,:),arc14mat(3,:),'y');

% draw hoops of sphere segment

% iterate over points on arc23

for i=1:1:length(arc24)
    di = double(dot(C1-arc24{i}',nP23));
    Ci = double((C1-di*nP23))';
    ri = double(sqrt(R^2-di^2));
%     plot3(arc24{i}(1),arc24{i}(2),arc24{i}(3),'kx');
%     plot3(arc34{end-i+1}(1),arc34{end-i+1}(2),arc34{end-i+1}(3),'kx');
    
    phi =linspace(0,double(acos(dot(arc24{i}-Ci,arc34{end-i+1}-Ci)/(norm(arc24{i}-Ci)*norm(arc34{end-i+1}-Ci)))),10);
%     plot3(Ci(1),Ci(2),Ci(3),'kx');
    for j=1:length(phi)
        arci(:,j) = Ci+double(vrrotvec2mat(double([nP23,phi(j)]))*(arc24{i}-Ci));
    end;
    h1(i)=plot3(arci(1,:),arci(2,:),arci(3,:),'k');
end;

ax = gca;
h2 = copyobj(h1,ax);
rotate(h2,T4,120);
h3 = copyobj(h2,ax);
rotate(h3,T4,120);
h4 = copyobj(h1,ax);
rotate(h1,T3,-120);

% convert to patch
XData=[];YData=[];ZData=[];
for i=1:length(ax.Children)
    XData = [XData,ax.Children(i).XData];
    YData = [YData,ax.Children(i).YData];
    ZData = [ZData,ax.Children(i).ZData];
end;

DT = delaunayTriangulation([XData,0]',[YData,0]',[ZData,0]');
[FBtri,FBpoints]=freeBoundary(DT);
trisurf(FBtri,FBpoints(:,1),FBpoints(:,2),FBpoints(:,3),'Facecolor','c','FaceAlpha',0.8,'EdgeColor','none');

% draw tetrahedron
%patch('Faces',T.faces,'Vertices',T.vertices,'FaceColor','r','EdgeColor','k','FaceAlpha',0.2);hold on;
% draw grain corner
plot3([0,1.5*T1(1)],[0,1.5*T1(2)],[0,1.5*T1(3)],'r');hold on;
plot3([0,1.5*T2(1)],[0,1.5*T2(2)],[0,1.5*T2(3)],'g');
plot3([0,1.5*T3(1)],[0,1.5*T3(2)],[0,1.5*T3(3)],'b');
plot3([0,1.5*T4(1)],[0,1.5*T4(2)],[0,1.5*T4(3)],'m');
% draw facets
F.vertices = double(1.5*[O;T1;T2;T3;T4]);
F.faces = [[1 2 3];[1 2 4];[1 2 5];[1 3 4];[1 3 5];[1 4 5]];
patch('Faces',F.faces,'Vertices',F.vertices,'FaceColor','k','EdgeColor','none','FaceAlpha',0.2);

axis equal
grid on
lighting gouraud
camlight

24.13 Figure TOC entry

%% Drawing Nuclei at Grain Edges and Grain Corners
clear all
close all

syms x

a1 = [1;0;0];
a2 = [-sym(1/2);sqrt(sym(3))/2;0];
a3 = [-sym(1/2);-sqrt(sym(3))/2;0];
z  = [0;0;1];
na1 = cross(a1,z);
na2 = cross(a2,z);
na3 = cross(a3,z);


R = 1;
d = 0.5;

% plot axes
% plot3([0,a1(1)],[0,a1(2)],[0,a1(3)],'r');hold on;
% plot3([0,a2(1)],[0,a2(2)],[0,a2(3)],'b');hold on;
% plot3([0,a3(1)],[0,a3(2)],[0,a3(3)],'g');hold on;
% plot3([0,0],[0,0],[-1,1],'k');hold on;



% find origin of sphere
O1 = d*R*a1;O2 = d*R*a2;O3 = d*R*a3;
% plot3([O1(1)],[O1(2)],[O1(3)],'bo');
% plot3([O2(1)],[O2(2)],[O2(3)],'go');
% plot3([O3(1)],[O3(2)],[O3(3)],'ro');
% draw wire model of sphere
azv = [0:pi/15:2*pi];
elv = [-pi/2:pi/15:pi/2];
[AZ EL RA] = meshgrid(azv,elv,R);
[X Y Z] = sph2cart(AZ, EL, RA);
X = X+O1(1);Y = Y+O1(2);Z = Z+O1(3);
sph1 = surf2patch(X,Y,Z);
%patch(sph1,'FaceColor','none','EdgeColor','k');

% find origin of circle of intersection between sphere and the a2-z plane
S=solve(dot(x*a2,O1-x*a2));
O12 = S(find(not(S==0)))*a2;
% plot3(O12(1),O12(2),O12(3),'bx');

% find origins on a3-z and a1-z plane by rotation
O23 = vrrotvec2mat([z;2*pi/3])*O12;
% plot3(O23(1),O23(2),O23(3),'gx');
O31 = vrrotvec2mat([z;2*pi/3])*O23;
% plot3(O31(1),O31(2),O31(3),'rx');

% find radius of circle of intersection between sphere and the a2-z plane
R12 = sqrt(R^2-norm(O12-O1)^2);

% find circle segments on positive half planes
S=solve(norm(x*z-O12)-R12);
% plot3(S*z(1),S*z(2),S*z(3),'bo');
% plot3(S*z(1),S*z(2),-S*z(3),'bo');
alpha = acos(norm(O12)/R12);
phi = [-alpha:pi/60:alpha];
C12=zeros(3,length(phi));
C23=zeros(3,length(phi));
C31=zeros(3,length(phi));

for i=1:length(phi)
    C12(:,i) = vrrotvec2mat(double([na2;phi(i)]))*(R12*a2);
    C23(:,i) = vrrotvec2mat([z;2*pi/3])*vrrotvec2mat(double([na2;phi(i)]))*(R12*a2);
    C31(:,i) = vrrotvec2mat([z;-2*pi/3])*vrrotvec2mat(double([na2;phi(i)]))*(R12*a2);
end;
C12X = O12(1)+C12(1,:);
C12Y = O12(2)+C12(2,:);
C12Z = O12(3)+C12(3,:);

C23X = O23(1)+C23(1,:);
C23Y = O23(2)+C23(2,:);
C23Z = O23(3)+C23(3,:);

C31X = O31(1)+C31(1,:);
C31Y = O31(2)+C31(2,:);
C31Z = O31(3)+C31(3,:);

plot3(C12X,C12Y,C12Z,'b');
plot3(C23X,C23Y,C23Z,'g');
plot3(C31X,C31Y,C31Z,'r');

% draw hoops for surface 1
[az12, el12] = cart2sph(double(C12X-O1(1)),double(C12Y-O1(2)),double(C12Z-O1(3)));
[az23, el23] = cart2sph(double(C23X-O1(1)),double(C23Y-O1(2)),double(C23Z-O1(3)));

N = 10;
[Xhoop_old,Yhoop_old,Zhoop_old] = sph2cart(az12(1)*ones(1,N),el12(1)*ones(1,N),R);
for i=2:2:length(el12)
    azhoop=linspace(az12(i),az23(i)+2*pi,N);
    [Xhoop,Yhoop,Zhoop] = sph2cart(azhoop,el12(i)*ones(size(azhoop)),1);
    plot3(Xhoop+O1(1),Yhoop+O1(2),Zhoop+O1(3),'k');  
    
    for j=1:length(Xhoop)
        plot3([Xhoop_old(j) Xhoop(j)]+O1(1),[Yhoop_old(j) Yhoop(j)]+O1(2),[Zhoop_old(j) Zhoop(j)]+O1(3),'k');hold on;
    end;
    
    Xhoop_old=Xhoop;
    Yhoop_old=Yhoop;
    Zhoop_old=Zhoop;    
end;
[Xhoop,Yhoop,Zhoop] = sph2cart(az12(1)*ones(1,N),el12(end)*ones(1,N),R);
for j=1:length(Xhoop)
    plot3([Xhoop_old(j) Xhoop(j)]+O1(1),[Yhoop_old(j) Yhoop(j)]+O1(2),[Zhoop_old(j) Zhoop(j)]+O1(3),'k');
end;


% draw hoops for surface 2
[az23, el23] = cart2sph(double(C23X-O2(1)),double(C23Y-O2(2)),double(C23Z-O2(3)));
[az31, el31] = cart2sph(double(C31X-O2(1)),double(C31Y-O2(2)),double(C31Z-O2(3)));

N = 10;
[Xhoop_old,Yhoop_old,Zhoop_old] = sph2cart(az23(1)*ones(1,N),el23(1)*ones(1,N),R);
for i=2:2:length(el23)
    azhoop=linspace(az23(i),az31(i),N);
    [Xhoop,Yhoop,Zhoop] = sph2cart(azhoop,el23(i)*ones(size(azhoop)),1);
    plot3(Xhoop+O2(1),Yhoop+O2(2),Zhoop+O2(3),'k');  
    
    for j=1:length(Xhoop)
        plot3([Xhoop_old(j) Xhoop(j)]+O2(1),[Yhoop_old(j) Yhoop(j)]+O2(2),[Zhoop_old(j) Zhoop(j)]+O2(3),'k');
    end;
    
    Xhoop_old=Xhoop;
    Yhoop_old=Yhoop;
    Zhoop_old=Zhoop;    
end;
[Xhoop,Yhoop,Zhoop] = sph2cart(az23(1)*ones(1,N),el23(end)*ones(1,N),R);
for j=1:length(Xhoop)
    plot3([Xhoop_old(j) Xhoop(j)]+O2(1),[Yhoop_old(j) Yhoop(j)]+O2(2),[Zhoop_old(j) Zhoop(j)]+O2(3),'k');
end;

% draw hoops for surface 3
[az31, el31] = cart2sph(double(C31X-O3(1)),double(C31Y-O3(2)),double(C31Z-O3(3)));
[az12, el12] = cart2sph(double(C12X-O3(1)),double(C12Y-O3(2)),double(C12Z-O3(3)));

N = 10;
[Xhoop_old,Yhoop_old,Zhoop_old] = sph2cart(az31(1)*ones(1,N),el31(1)*ones(1,N),R);
for i=2:2:length(el31)
    azhoop=linspace(az31(i),az12(i),N);
    [Xhoop,Yhoop,Zhoop] = sph2cart(azhoop,el31(i)*ones(size(azhoop)),1);
    plot3(Xhoop+O3(1),Yhoop+O3(2),Zhoop+O3(3),'k');  
    
    for j=1:length(Xhoop)
        plot3([Xhoop_old(j) Xhoop(j)]+O3(1),[Yhoop_old(j) Yhoop(j)]+O3(2),[Zhoop_old(j) Zhoop(j)]+O3(3),'k');
    end;
    
    Xhoop_old=Xhoop;
    Yhoop_old=Yhoop;
    Zhoop_old=Zhoop;    
end;
[Xhoop,Yhoop,Zhoop] = sph2cart(az31(1)*ones(1,N),el31(end)*ones(1,N),R);
for j=1:length(Xhoop)
    plot3([Xhoop_old(j) Xhoop(j)]+O3(1),[Yhoop_old(j) Yhoop(j)]+O3(2),[Zhoop_old(j) Zhoop(j)]+O3(3),'k');
end;

ax = gca;
% convert to patch?

XData=[];YData=[];ZData=[];
for i=1:length(ax.Children)
    XData = [XData,ax.Children(i).XData];
    YData = [YData,ax.Children(i).YData];
    ZData = [ZData,ax.Children(i).ZData];
end;
%
DT = delaunayTriangulation([XData,0]',[YData,0]',[ZData,0]');
[FBtri,FBpoints]=freeBoundary(DT);
trisurf(FBtri,FBpoints(:,1),FBpoints(:,2),FBpoints(:,3),'Facecolor','c','FaceAlpha',0.8,'EdgeColor','none');
lighting gouraud
camlight
axis equal

%% plot half planes

% a1-z plane
p1.vertices = R*[0 0 1;1 0 1; 1 0 -1; 0 0 -1];
p2.vertices = double(R*[0 0 1;a2(1) a2(2) 1; a2(1) a2(2) -1; 0 0 -1]);
p3.vertices = double(R*[0 0 1;a3(1) a3(2) 1; a3(1) a3(2) -1; 0 0 -1]);
p.faces    = [1 2 3 4];

patch('Faces',p.faces,'Vertices',p1.vertices,'FaceColor','k','EdgeColor','none','FaceAlpha',0.2);
patch('Faces',p.faces,'Vertices',p2.vertices,'FaceColor','k','EdgeColor','none','FaceAlpha',0.2);
patch('Faces',p.faces,'Vertices',p3.vertices,'FaceColor','k','EdgeColor','none','FaceAlpha',0.2);

lighting gouraud
light('Position',[-2 0 4],'Style','local');
view(25,20);

24.14 Figure TOC entry

%% solidification in binary systems
clear all 
close all
clc

R = 8.314; % J mol^-1 K^-1
k = 0.5;
X_0 = 0.05;
p = '/Users/derk/Box Sync/Teaching/MSE 316-2 Fall 2018/Notes (F18)/Tex Files and Figures/';

%% perfect mixing in liquid and solid
Tm  = 1000; %K
b_L = -400; %K
b_s = -800; %K
k = b_L/b_s;

X_0 = 0.15;

T_s = @(X) Tm + b_s*X; 
T_L = @(X) Tm + b_L*X;

X_L = @(T) (T - Tm)/b_L;
X_s = @(T) k*X_L(T);

f_s = @(T) (X_L(T)-X_0)./(X_L(T)*(1-k)); % assuming molar volumes indentical 

T_i = T_L(X_0);
T_f = T_s(X_0);

% plot PD
close all
figure;
X = [0:0.01:0.5];
plot(X,T_L(X),'b');hold on;
plot(X,T_s(X),'g');
ax=gca;
plot([X_0,X_0],ax.YLim,'k:');
plot([0,X_0],[T_i,T_i],'k:');
plot([0,X_0/k],[T_f,T_f],'k:');
legend('liquidus','solidus');
ax=gca;
ax.TickDir='out';
saveas(gcf,[p,'Lecture13_Figure3A'],'epsc');
% plot f_s vs T
figure;
T = T_i:-1:T_f;
plot(f_s(T),T);
ylabel('\rightarrow T [K]');
xlabel('\rightarrow f_s');
ax=gca;
ax.TickDir='out';
saveas(gcf,[p,'Lecture13_Figure3B'],'epsc');
% plot X_L and X_s vs f_s
figure;
plot(f_s(T),X_L(T),'b');hold on;
plot(f_s(T),X_s(T),'g');
ylabel('\rightarrow X');
xlabel('\rightarrow f_s');
legend('liquid','solid');
ax=gca;
ax.TickDir='out';
saveas(gcf,[p,'Lecture13_Figure3C'],'epsc');

% plot snapshots of solidification
N = 3; % number of intermediate temperatures
figure;
Tv = [T_i:(T_f-T_i)/(N+1):T_f];

for i=1:N+2
    subplot(N+2,1,i);
    plot([f_s(Tv(i)),1],[X_L(Tv(i)),X_L(Tv(i))],'bo-');hold on;
    plot([0,f_s(Tv(i))],[X_s(Tv(i)),X_s(Tv(i))],'go-');
    plot([0,1],[X_0,X_0],'k:');
    plot([f_s(Tv(i)),f_s(Tv(i))],[0,X_0/k+0.05],'r');
    ax=gca;
    ax.YLim = [0,X_0/k+0.05];
    ax.TickDir='out';
    ylabel('\rightarrow X');
end
xlabel('\rightarrow f_s');
legend('liquid','solid','X_0','interface');
saveas(gcf,[p,'Lecture13_Figure4'],'epsc');

%% no diffusion in solid, perfect mixing in liquid
close all
% Scheil equation
X_s = @(f_s) k*X_0*(1-f_s).^(k-1);
X_L = @(f_s) X_0*(1-f_s).^(k-1);

f = [0:0.01:1];
figure;
plot(f,X_s(f));hold on
plot(f,X_L(f));hold on
plot([0,1],[X_0,X_0],'k:'); % initial composition
ax=gca;
ax.TickDir='out';
xlabel('\rightarrow fractional progress of solidification f_{\rm{s}}');
ylabel('\rightarrow X');
legend('solid','liquid','initial');
saveas(gcf,[p,'Lecture13_Figure5'],'epsc');

% plot snapshots of solidification
N = 3; % number of intermediate temperatures
f_sv = [0:1/(N+1):1];

figure;
for i=1:N+2
    subplot(N+2,1,i);
    plot([f_sv(i),1],[X_L(f_sv(i)),X_L(f_sv(i))],'bo-'); hold on;
    f_s = 0:0.01:f_sv(i);
    plot(f_s,X_s(f_s),'g');
    plot([0,1],[X_0,X_0],'b:');
    X_S_av = mean(X_s(f_s));
    plot(f_s([1,end]),[X_S_av,X_S_av],'m:');
    plot([f_sv(i),f_sv(i)],[0,X_L(0.9)],'r');
    ax=gca;
    ax.YLim = [0,X_L(0.9)];
    ax.TickDir='out';
    ylabel('\rightarrow X');
end
legend('liquid','solid','X_0','mean(X_s)','interface');
xlabel('\rightarrow fractional progress of solidification f_{\rm{s}}');
saveas(gcf,[p,'Lecture13_Figure6'],'epsc');

%% no diffusion in solid, only diffusion in liquid
close all;

D = 1E-9; % m^2 s^-1
x = [0:1e-8:1e-5];%m
v = 1e-3; % m/s

X_L = X_0*(1+(1-k)/k*exp(-x*v/D));
figure;
plot(x,X_L);hold on;
plot([0,0],[0,X_L(1)+0.05],'r');
plot([-2*x(end),x(end)],[X_0,X_0],'b:');
plot([-2*x(end),x(end)],[k*X_0,k*X_0],'b:');
plot([-2*x(end),x(end)],[X_0/k,X_0/k],'b:');
ax=gca;
ax.TickDir='out';
ax.YLim = [0,X_L(1)+0.05];
ax.XLim = [-2*x(end),x(end)];
saveas(gcf,[p,'Lecture13_Figure7'],'epsc');

%% local liquidus temperature
figure;
% plot local liquidus
plot(x,T_L(X_L)); hold on;
% indicate interface
plot([0,0],[T_L(X_0/k)-100,T_L(X_0)+50],'r');
% plot T_L(X_0) isotherm
plot([-2*x(end),x(end)],[T_L(X_0),T_L(X_0)],'b:');
% plot T_L(X_0/k) isotherm
plot([-2*x(end),x(end)],[T_L(X_0/k),T_L(X_0/k)],'b:');
% plot critical T-gradient
fun_gradT_crit = @(x) T_L(X_0/k)+(T_L(X_0)-T_L(X_0/k))*v/D*x;
xv = [0:1e-6:D/v];
plot(xv,fun_gradT_crit(xv),'m');
ax=gca;
ax.TickDir='out';
ax.YLim = [T_L(X_0/k)-100,T_L(X_0)+50]
ax.XLim = [-2*x(end),x(end)];
saveas(gcf,[p,'Lecture13_Figure8'],'epsc');

24.15 Figure TOC entry

%% Lecture 11 Figure 1

% precipitate growth: interface position and velocity
k = 1.38e-23; %J/K
R = 8.314; %J/mol/K
Qd= 136000; %J/mol
D0= 6.5e-5; %m^2/s

DeltaC = 0.05;
C_beta = 0.98;
C_alpha= 0.02;

D = @(T)   D0 * exp(-Qd/R./T);
h = @(T,t) DeltaC/(C_beta-C_alpha)*sqrt(D(T).*t);
v = @(T,t) DeltaC/2/(C_beta-C_alpha)*sqrt(D(T)./t);

DeltaC = 0.05;
C_beta = 0.95;
C_alpha= 0.05;

Tv = 0:1:700;

close all
figure
plot(sqrt(D(Tv)),Tv)

24.16 Figure TOC entry

%% Lecure 12: Precipitate Growth
clear all 
close all
clc
% 1D growth rate of an incoherent precipitate as a function of temperature
% example: carbon diffusing in ferrite
D0 = 6.2E-7; % m^2s^-1 C in ferrite
Qd = 80000;  % J mol^-1
R  = 8.314; % J mole^-1 K^-1
D = @(T) D0*exp(-Qd/R./T);

T_eu   = 727+273; %K
Xs_max = 0.022;
X_0 = 0.02;

Ts   = @(X) T_eu/Xs_max*X;
Xs   = @(T) Xs_max/T_eu*T;

Tv = linspace(0,Ts(X_0),1000);

Delta_X = X_0 - Xs(Tv);

yyaxis left
plot(Tv,Delta_X/max(Delta_X),'b');hold on;

yyaxis right
plot(Tv,D(Tv)/max(D(Tv)),'g');


plot(Tv,Delta_X.*D(Tv)/max(Delta_X.*D(Tv)),'r')

24.17 Figure TOC entry

%% Lecture 16,17 Figure 5
clear all
% growth velocity of a curved interface
R  = 8.314;   % J mol^-1 K^-1
Na = 6.023E23;% mol^-1
kb = R/Na;    % J K^-1

D = @(D_0,Qd,T) D_0*exp(-Qd/R/T);
v = @(D,DX_0, r, r_star, X_ppt, X_matrix_eq) (D .* DX_0)./(X_ppt - X_matrix_eq).* r.^-1 .* (1 - r_star./r);
w2X = @(w1,M1,M2) w1/M1/(w1/M1+(1-w1)/M2);  
X2w = @(X1,M1,M2) X1*M1/(X1*M1+(1-X1)*M2);
r_star = @(X_matrix_0,X_matrix_eq,gamma,Vm_ppt,T) 2*gamma*Vm_ppt/R/T/(X_matrix_0/X_matrix_eq-1);

% for visualization of terms going into v
vT1 = @(D,DX_0, r, r_star, X_ppt, X_matrix_eq) (D .* DX_0)./(X_ppt - X_matrix_eq).* r.^-1;
vT2 = @(D,DX_0, r, r_star, X_ppt, X_matrix_eq) -(D .* DX_0)./(X_ppt - X_matrix_eq).* r_star./r.^2;

% example 1
% Cu in Al, solid state
D_0_Cu         = 6.5E-5; % m^2 s^-1 from Callister
Qd_Cu          = 136000; % J mol^-1
M_Al           = 27;     % g mol^-1
M_Cu           = 63.55;  % g mol^-1
Vm_theta       = 9e-6;   % m^3 mol^-1 from Y. Ocak et al. / Thin Solid Films 518 (2010) 4322?4327 and ref therein
gamma_theta_Al = 88e-3;  % J m^2 from Y. Ocak et al. / Thin Solid Films 518 (2010) 4322?4327 and ref therein

% calcuate diffusivity at 400?C
T      = 673;    % K
D_Cu = D(D_0_Cu,Qd_Cu,T);

% composition of alpha Al and theta CuAl2 at 400?C, from phase diagram
w_alpha_eq = 0.03; % (w/w)
w_theta = 0.5369; % (w/w)
X_alpha_eq = w2X(w_alpha_eq,M_Cu,M_Al);
X_theta = w2X(w_theta,M_Cu,M_Al);

% set up starting condition
w_0 = 0.0305; 
X_0 = w2X(w_0,M_Cu,M_Al);
DeltaX0 = X_0-X_alpha_eq;

r_crit = r_star(X_0,X_alpha_eq,gamma_theta_Al,Vm_theta,T); % m
r = logspace(floor(log10(r_crit/10)),ceil(log10(r_crit*10)),100); % m
v_theta = v(D_Cu,DeltaX0, r, r_crit, X_theta, X_alpha_eq);
v_max = v(D_Cu,DeltaX0, 2*r_crit, r_crit, X_theta, X_alpha_eq);

v_theta_T1 = vT1(D_Cu,DeltaX0, r, r_crit, X_theta, X_alpha_eq);
v_theta_T2 = vT2(D_Cu,DeltaX0, r, r_crit, X_theta, X_alpha_eq);


close all
figure;
subplot(121);
plot(r*1e9,v_theta*1e9);hold on;
plot(r*1e9,v_theta_T1*1e9,'g--',r*1e9,v_theta_T2*1e9,'r--');

set(gca,'TickDir','out','XMinorTick','on','XMinorGrid','on');
xlim([0,15*r_crit*1e9]);
ylim([-2*v_max*1e9,4*v_max*1e9]);
grid on
xlabel('\rightarrow r [nm]');
ylabel('\rightarrow v [nm/s]');

subplot(122);
plot(r/r_crit,v_theta*1e9);hold on;
plot(r/r_crit,v_theta_T1*1e9,'g--',r/r_crit,v_theta_T2*1e9,'r--');

set(gca,'TickDir','out','XMinorTick','on','XMinorGrid','on');
ylim([-2*v_max*1e9,4*v_max*1e9]);
xlim([0,15]);
grid on
xlabel('\rightarrow r/r^*');
ylabel('\rightarrow v [nm/s]');
% example 2
% protein in water / polymer in solvent

24.18 Figure TOC entry

%% Coarsening
clear all
close all

imx = 200;
beta = zeros(imx,imx);

centers = rand(10,2)*imx;

[V C] =voronoin(centers);
voronoi(centers(:,1),centers(:,2)); hold on;
syms x1 x2 m
X =[x1, x2];
%%
j = 3
c = C{j}
%for i=1:length(C{j}) 

i=1;
    V1V2 = V(c(2),:)-V(c(1),:)
    plot(V(c(2),1),V(c(2),2),'go',V(c(1),1),V(c(1),2),'go',centers(j,1),centers(j,2),'r+');
    
    OX   = X-centers(j,:);
    
    S=solve([dot(V1V2,OX)==0, X == V(c(1),:) + m*V1V2],[x1 x2 m])
    
    double([S.x1, S.x2, S.m])
    double(V(c(1),:) + S.m*V1V2)
    
    plot(double(V(c(1),1) + S.m*V1V2),double(V(c(1),2) + S.m*V1V2),'ro');
    
    %c=circshift(c,[0 -1]);

24.19 Figure TOC entry

%% Spinodal decomposition in one dimension

% initial system described by fluctuations with wavelengths lambda_i
% vector vl gives the wavelength of these fluctuations as multiples of
% lambda_c, the critical wavelength

lambda_c = 50e-9; % m
vl = [0.2 0.5, 0.75, 1, sqrt(2), 2.5];
lambda = lambda_c*vl;

t = [0,10,50,200,1000]; %

% the matrix Am(lambda,t) gives the amplitude of the fluctuations w/wavelengths given
% in vl at time t

Am = zeros(length(vl),length(t));
Am(:,1) = ones(size(vl)); % initial Amplitudes are all 1

% calculate the solution to Cahn's Diffusion equation 
% here, we don't use an explicit example, so kappa and M are choses purely
% for convenience of visualization
kappa =1;  % Jm^5mol^-2
M = 1e-35; % mol^2 J^-1m^-1s^-1

R = @(l,l_c) 2*kappa*M*(2*pi./l).^4.*(l.^2/l_c^2-1);

for i=2:length(t)
    Am(:,i)=exp(R(lambda,lambda_c)*t(i));
end;
Amax = max(max(Am));

close all
z = linspace(0,5*lambda_c,1000);
C = zeros(length(vl),length(t),length(z));
for j=1:length(t)
    for i = 1:length(vl)
        subplot(length(vl)+1,length(t),(i-1)*length(t)+j);
        C(i,j,:) = Am(i,j)*cos(2*pi*z/lambda(i));
        plot(z/lambda_c,squeeze(C(i,j,:)));hold on
        ylim([-Amax,Amax]);
        xlim([z(1),z(end)]/lambda_c);
        set(gca,'TickDir','out','XMinorGrid','on','YMinorGrid','on','XMinorTick','on','YMinorTick','on');
        grid on
    end;
    subplot(length(vl)+1,length(t),i*length(t)+j);
    plot(z/lambda_c,sum(squeeze(C(:,j,:)),1));hold on
    %ylim([-Amax,Amax]);
    xlim([z(1),z(end)]/lambda_c);
    xlabel('\rightarrow z/\lambda_c');
    ylabel('\rightarrow C(z,t)-C_0');
    set(gca,'TickDir','out','XMinorGrid','on','YMinorGrid','on','XMinorTick','on','YMinorTick','on');
    grid on
end;

24.20 Figure TOC entry

%% Spinodal Decomposition Figures
clear all
close all
% Figure 1

z = 0:1:100;
c = 0.1*rand(size(z));
c = c-mean(c);

lambda = 100;
delta  = 0.1;

figure;
subplot(1,2,1);
plot(z,c);
ylim([-2*delta,2*delta]);
set(gca,'TickDir','out','YMinorTick','off','YTick',[-2:1:2]*delta,'YTickLabel',{'','','C_\circ','',''});
xlabel('\rightarrow z');
ylabel('\rightarrow C(z)-C_\circ');
grid on

subplot(1,2,2);
plot(z/lambda,delta*sin(2*pi*z/lambda));
ylim([-2*delta,2*delta]);
set(gca,'TickDir','out','YMinorTick','off','XTick',[0:.25:1],'XTickLabel',{'0','\lambda/4','\lambda/2','3\lambda/4','\lambda'},'YTick',[-2:1:2]*delta,'YTickLabel',{'-2\delta','-\delta','C_\circ','\delta','2\delta'});
xlabel('\rightarrow z');
ylabel('\rightarrow C(z)-C_\circ');
grid on

%% Figure 2
clear all
close all

f = @(delta, lambda, z) delta*sin(2*pi*z/lambda);
syms delta lambda z
df= matlabFunction(diff(f(delta, lambda, z),z),'Vars',[delta, lambda, z]);
clear delta lambda z

lambda = [20,50,100];
delta = 1;

for i=1:length(lambda)
    z = 0:0.01:lambda(i);
    plot(z,f(delta,lambda(i),z));hold on; % plot concentration fluctuation
    slope = df(delta, lambda(i), lambda(i)/2); % find slope at inflection point
    dz = -delta/slope;
    plot([lambda(i)/2-dz, lambda(i)/2+dz],[delta, -delta]);
    ylim([-2*delta,2*delta]);
    set(gca,'TickDir','out','YMinorTick','off','YTick',[-2:1:2]*delta,'YTickLabel',{'-2\delta','-\delta','C_\circ','\delta','2\delta'});
    xlabel('\rightarrow z');
    ylabel('\rightarrow C(z)');
    grid on
end;
%% Figure 3
clear all
close all


G = @(X2,T,Tc,R) R*T.*(1-X2).*log(1-X2)+R*T.*(X2).*log(X2)+2*R*Tc.*(1-X2).*X2;
syms X2 T Tc R
dG = matlabFunction(diff(G(X2,T,Tc,R),X2),'Vars',[X2,T,Tc,R]);
d2G = matlabFunction(diff(dG(X2,T,Tc,R),X2),'Vars',[X2,T,Tc,R]);
Tsp = matlabFunction(solve(d2G(X2,T,Tc,R),T),'Vars',[X2,Tc]);
Xsp = matlabFunction(solve(d2G(X2,T,Tc,R),X2),'Vars',[T,Tc]);

R = 8.314; % J mol^-1 K^-1
X     = [0:0.001:1];
Tcrit = [1000,800]; % K
c = ['b','r'];
s = {'-','-'};
cm = colormap(jet(64));
Tmin = 0;
Tmax = max(Tcrit);

figure;
for j=1:length(Tcrit)
    subplot(2,1,j);
    Tact  = [0,100:50:Tcrit(j)-200,Tcrit(j)-190:10:Tcrit(j)-20,Tcrit(j)-18:2:Tcrit(j)-2,Tcrit(j)];  % K for Fig 3b
    %Tact  = [0,100:100:Tcrit(j)-2,Tcrit(j)];  % K for Fig 3a
    i=1;col = cm(1+floor((Tact(i)-Tmin)/Tmax*(length(cm)-1)),:);
    Xspin(:,i) = [1,0];
    Gspin(:,i) = G(Xspin(:,1),Tact(1),Tcrit(j),R);
    Xbin(:,i) = [1,0];
    Gbin(:,i) = G(Xbin(:,1),Tact(1),Tcrit(j),R);
    plot(Xspin(:,i),Gspin(:,i),'MarkerEdgeColor',col,'Marker','o','LineStyle','none');hold on;
    plot(Xbin(:,i),Gbin(:,i),'MarkerEdgeColor',col,'Marker','x','LineStyle','none');
    plot(X,G(X,Tact(i),Tcrit(j),R),'Color',col,'LineStyle',s{j});hold on;
    colormap(jet(64));
    title(['T_c = ',num2str(Tact(end)),' K']);
    for i=2:length(Tact)-1;
        col = cm(1+floor((Tact(i)-Tmin)/Tmax*(length(cm)-1)),:);
        plot(X,G(X,Tact(i),Tcrit(j),R),'Color',col,'LineStyle',s{j});hold on;
        
        Xspin(:,i) = Xsp(Tact(i),Tcrit(j));
        Gspin(:,i) = G(Xspin(:,i),Tact(i),Tcrit(j),R);
        
        Xbin(:,i) = [fzero(@(X) dG(X,Tact(i), Tcrit(j), R),[0.5+1e-10 1-1e-10]), fzero(@(X) dG(X,Tact(i), Tcrit(j), R),[1e-10 0.5-1e-10])];
        Gbin(:,i) = G(Xbin(:,i),Tact(i),Tcrit(j),R);
        
        plot(Xspin(:,i),Gspin(:,i),'MarkerEdgeColor',col,'Marker','o','LineStyle','none');hold on;
        plot(Xbin(:,i),Gbin(:,i),'MarkerEdgeColor',col,'Marker','x','LineStyle','none');
    end;
    i=i+1;
    col = cm(1+floor((Tact(i)-Tmin)/Tmax*(length(cm)-1)),:);
    
    plot(X,G(X,Tact(i),Tcrit(j),R),'Color',col,'LineStyle',s{j});hold on;
    
    Xspin(:,i) = [.5 .5];
    Gspin(:,i) = G(Xspin(:,i),Tact(i),Tcrit(j),R);  
    plot(Xspin(:,i),Gspin(:,i),'MarkerEdgeColor',col,'Marker','o','LineStyle','none');
 
    Xbin(:,i) =  [.5 .5];
    Gbin(:,i) = G(Xbin(:,i),Tact(i),Tcrit(j),R);
    plot(Xbin(:,i),Gbin(:,i),'MarkerEdgeColor',col,'Marker','x','LineStyle','none');
    
    data.spinodal.Tc = Tcrit(j); 
    data.spinodal.T = Tact;
    data.spinodal.X = Xspin;
    data.spinodal.G = Gspin;
    data.binodal.Tc = Tcrit(j); 
    data.binodal.T = Tact;
    data.binodal.X = Xbin;
    data.binodal.G = Gbin;
    GXTX{j}=data;
    clear Xspin Xbin Gspin Gbin
    set(gca,'TickDir','out');
    xlabel('\rightarrow X');
    ylabel('\rightarrow G');
    grid on
    colorbar
end;
 
%
figure;
for j=1:length(GXTX)
    plot(GXTX{j}.binodal.X(1,:),GXTX{j}.binodal.T,[c(j),'-'],GXTX{j}.binodal.X(2,:),GXTX{j}.binodal.T,[c(j),'-']);hold on;
    plot(GXTX{j}.spinodal.X(1,:),GXTX{j}.spinodal.T,[c(j),':'],GXTX{j}.spinodal.X(2,:),GXTX{j}.spinodal.T,[c(j),':']);hold on;
end;
set(gca,'TickDir','out');
xlabel('\rightarrow X');
ylabel('\rightarrow T');
grid on    


25 Lab Examples

In this section we take a closer look at some of the data from the laboratories.

25.1 Coarsening

image: 115_home_ken_Mydocs_MSEcore_316-2_figures_lab2_fig1.png
Figure 25.1: Phase diagram for the N H 4 N O 3 /NaN O 3 system.
image: 116_home_ken_Mydocs_MSEcore_316-2_figures_lab_dendrite_image.png
Figure 25.2: Observed Dendritic Morphology
Note: λ 1 (primary arm spacing) 10 λ 2 (secondary arm spacing)

25.2 Secondary Arm Spacing

image: 117_home_ken_Mydocs_MSEcore_316-2_figures_lab_dendrite_spacing.svg
Figure 25.3: Measured Values of the Secondary Arm Spacing

25.3 Surface to Volume Ratio

image: 118_home_ken_Mydocs_MSEcore_316-2_figures_lab_dendrite_svplot.svg
Figure 25.4: Measured surface to volume ratio (units missing)

25.4 Temperature Dependence: K

image: 119_home_ken_Mydocs_MSEcore_316-2_figures_lab_dendrite_svplotv2.svg
Figure 25.5: Temperature Dependence of the Coarsening Constant
%some initialization stuff
clear all
close all
set(0, 'DefaultAxesFontSize', 16);
set(0, 'DefaultTextFontSize', 16);
set(0, 'DefaultLineLineWidth', 2);
set(0, 'DefaultFigurePaperposition',[0 0 6 4])
set(0, 'DefaultFigurePaperSize',[6 4])

%input the data
t=60*[0 6 12 18 24 36 48 60];
l=[7.94 13.3 14.3 16.7 16.6 25 28.6 28.6];
sv=[0.364 0.32 0.258 0.222 0.196 0.178 0.16 0.16];

% create the first plot
plot(t, l.^3, '+b');
hold on
spacingfit=polyfit(t, l.^3, 1); % generates a linear curve fit
plot(t,polyval(spacingfit,t),'b-'); % plots the fit
hold off
xlabel('t (s)')
ylabel('\lambda_{2}^{3} (\mum^{3})')
text(200,2.2e4,'\lambda_{2}^{3}=A+Kt')
ktext=num2str(1e-18*spacingfit(1),'%5.1e');
text(200,1.8e4,['K =' ktext ' m^3/s'])
print(gcf,'-depsc2','../316-2_figures/lab_dendrite_spacing.eps')

% now make the sv plot
figure % creates a new figure window
plot(t, sv.^-3, 'b+');
hold on
svfit=polyfit(t, sv.^(-3),1);
plot(t,polyval(svfit,t))
xlabel('t (s)')
ylabel('S_{v}^{-3} (unspecified units)')
text(200,270,'S_{v}^{-3}=A+Kt')
ktext=num2str(1e-18*svfit(1),'%5.1e');
text(200,220,['K =' ktext])
print(gcf,'-depsc2','../316-2_figures/lab_dendrite_svplot.eps')

% now include another dataset (from Monday PM group)
figure
t=[0 80 160 240 320]; % time in seconds
sv1311=[0.258 0.157 0.10 0.099 0.075];
sv1338=[0.144 0.113 0.096 0.085 0.0756];
plot(t, sv1311.^-3, 'b+', t, sv1338.^-3, 'ro');
svfit1=polyfit(t, sv1311.^(-3),1);
svfit2=polyfit(t, sv1338.^(-3),1);
hold on
plot(t,polyval(svfit1,t),'b-', t, polyval(svfit2,t),'r-')
xlabel('t (s)')
ylabel('S_{v}^{-3} (\mum^{3})')
% now generate legend info
ktext1=num2str(1e-18*svfit1(1),'%5.1e');
ktext2=num2str(1e-18*svfit2(1),'%5.1e');
legendtext{1}=['131.1 ^{\circ}C, K=' ktext1 ' m^{3}/s'];
legendtext{2}=['133.8 ^{\circ}C, K=' ktext2 ' m^{3}/s'];
legend(legendtext,'location','best')
print(gcf,'-depsc2','../316-2_figures/lab_dendrite_svplotv2.eps')

25.5 LSW Theory

K= 8 γ V m D C e 9RT( C β - C e )

25.6 Estimating D

λ 1 D V
DV λ 1 10 λ 2 V
T ( C)
V ( μ m/s)
λ 2 ( μ m)
10 λ 2 V( m 2 s )
130.4
150
9.3
1.4x10 -9
130.9
130
11
1.4x10 -9
133.6
86
13
1.1x10 -9
135
45
16
0.7x10 -9

25.7 Estimating the Expected Coarsening Constant

K γ V m D RT
V m =( 80g mol ) ( c m 3 1.7 g ) =47 c m 3 /mol=4.7x1 0 -5 m 3 /mol
D1 0 -9 m 2 /s
γ 0.1 J/ m 2
K ( 0.1 J/ m 2 ) ( 4.7x1 0 -5 m 3 /mo ) ( 1 0 -9 m 2 /s ) ( 8.314 J.mo-K ) ( 400 K ) =1.4x1 0 -18 m 3 /s

25.8 Heat treatment of Al alloys

image: 120_home_ken_Mydocs_MSEcore_316-2_figures_Cu_Al_age-130.PNG
Figure 25.6: Heat Treatment of Al Alloys (from P&E)
image: 121_home_ken_Mydocs_MSEcore_316-2_figures_Cu_Al_age-190.PNG
Figure 25.7: (P&E 5.5.4)

26 Appendix: Thermodynamic Data for Cu

%% Phase transformations in elemental copper
% This script uses NIST thermochemical data to 
%       - create plots of thermodynamic functions of state
%       - calculate phase transformation temperatures, and the changes in
%       standard enthalpy, entropy, and Gibbs free energy associated with
%       the phase transformation
%       - compare the change in Gibbs free energy calculated from NIST data
%       to the linearized Gibbs free energy estimate near the phase
%       transformation temperature. 

close all;  clear all; filename = 'Cu_Shomate_';
%% NIST Data
% The NIST Chemistry webbook (URL: http://webbook.nist.gov/chemistry/)
% gives Shomate coefficients that can be used to calculate the heat
% capacity, standard enthalpy, absolute entropy, and standard Gibbs free
% energy as a function of temperature. The Shomate coefficients are define
% only over a certain range of T. Extrapolation beyond this range is
% dangerous.

% NIST Chemistry webbok data for iron can be found here:
% http://webbook.nist.gov/cgi/cbook.cgi?Name=copper&Units=SI
% This is actually a fairly straightforward case, as there only two phase
% transformations: solid (alpha, hcp) -> L and L -> gas

Cu.phase={'s','L','g'};
Cu.Tmin = [298, 1358, 2843.261];
Cu.Tmax = [1358, 2843, 6000];

S=[[17.72891,  32.84450,   -80.48635];
   [28.09870,  -0.000084,   49.35865];
   [-31.25289,  0.000032,   -7.578061];
   [13.97243,  -0.000004,    0.404960];
   [0.068611,  -0.000028,  133.3382];
   [-6.056591, -1.804901,  519.9331];
   [47.89592,  73.92310,   193.5351];
   [0,         11.85730,   337.6003];];

%% Calculating Thermodynamic Functions of State from Shomate Coefficients
% We use here anonymous functions to calculate the thermodynamic functions
% of state. Input arguments are a column vector of exactly 8 Shomate
% coefficients (A-H) and a row vector with temperature values where t = T/1000 [K]. Please note that H0
% is really H0, not H0-H0_298! Cp0 and S0 are given in units of J mol^-1 K^-1, H0
% and G0 in units of kJ mol^-1.

Cp0 = @(Shomate,t) Shomate*[0*t+1 ; t     ; t.^2  ; t.^3  ;  t.^-2    ; 0*t  ; 0*t  ; 0*t]; %[J mol^-1 K^-1]

H0  = @(Shomate,t) Shomate*[t     ; t.^2/2; t.^3/3; t.^4/4; -t.^-1    ; 0*t+1; 0*t  ; 0*t]; %[kJ mol^-1]

S0  = @(Shomate,t) Shomate*[log(t); t     ; t.^2/2; t.^3/3; -0.5*t.^-2; 0*t  ; 0*t+1; 0*t]; %[J mol^-1 K^-1]
 
G0 = @(Shomate, T) H0(Shomate, T/1000)-T/1000.*S0(Shomate, T/1000); % [kJ mol^-1]

%% Calculate phase transformation temperatures
% We define the Gibbs free energy change for four different
% transformations in the system as anonymous functions. To solve for the
% phase transformation tempreature at which the Gibbs free energy change is
% zero, we find the root of the anonymous functions using a guess based on
% the phase transformation temperatures given in phase diagrams.

% vaporization: L->gas
DeltaG_Lg= @(T) G0(S(:,3)',T)-G0(S(:,2)',T);
T_Lg = fzero(DeltaG_Lg,3130);

% melting/fusion: s->L
DeltaG_sL= @(T) G0(S(:,2)',T)-G0(S(:,1)',T);
T_sL = fzero(DeltaG_sL,1800);


%% Plot Functions of State

% T-ranges
T1=linspace(298,T_sL,1000);              
T2=linspace(T_sL,T_Lg,1000);             
T3=linspace(T_Lg,6000,1000);

h1=figure;
set(gcf,'DefaultLineLineWidth',1.5)

c=colormap(lines(3));
plot(T1,Cp0(S(:,1)',T1/1000),'Color',c(1,:));hold on;
plot(T2,Cp0(S(:,2)',T2/1000),'Color',c(2,:));hold on;
plot(T3,Cp0(S(:,3)',T3/1000),'Color',c(3,:));hold on;

title('Heat capacity (specific heat) of elemental Cu');
ylabel('C_p^0 [J mol^-^1 K^-^1]');
xlabel('T [K]');
legend(Cu.phase);
set(gca,'TickDir','out');
grid on;

h2=figure;
set(gcf,'DefaultLineLineWidth',1.5)
plot(T1,H0(S(:,1)',T1/1000),'Color',c(1,:));hold on;
plot(T2,H0(S(:,2)',T2/1000),'Color',c(2,:));hold on;
plot(T3,H0(S(:,3)',T3/1000),'Color',c(3,:));hold on;

title('Standard Enthalpy of elemental Cu');
ylabel('H^0-H_2_9_8 [kJ mol^-^1]');
xlabel('T [K]');
legend(Cu.phase);
set(gca,'TickDir','out');
grid on;

h3=figure;
set(gcf,'DefaultLineLineWidth',1.5)
plot(T1,S0(S(:,1)',T1/1000),'Color',c(1,:));hold on;
plot(T2,S0(S(:,2)',T2/1000),'Color',c(2,:));hold on;
plot(T3,S0(S(:,3)',T3/1000),'Color',c(3,:));hold on;

title('Absolute Entropy of elemental Cu');
ylabel('S^0 [J mol^-^1 K^-^1]');
xlabel('T [K]');
legend(Cu.phase);
set(gca,'TickDir','out');
grid on;

h4=figure;
set(gcf,'DefaultLineLineWidth',1.5)
plot(T1,G0(S(:,1)',T1),'Color',c(1,:));hold on;
plot(T2,G0(S(:,2)',T2),'Color',c(2,:));hold on;
plot(T3,G0(S(:,3)',T3),'Color',c(3,:));hold on;

title('Gibbs Free Energy of elemental Cu');
ylabel('G^0 [kJ mol^-^1]');
xlabel('T [K]');
legend(Cu.phase);
set(gca,'TickDir','out');
grid on;


%% Plot Close Ups of Phase Transformations
h5=figure;
set(gcf,'DefaultLineLineWidth',1.5)
% solid-liquid transition
DT = 2;
T_tr = T_sL;
p1=1;
p2=2;

Tr_left = (T_tr-DT):0.1:T_tr;
Tr_right = T_tr:0.1:(T_tr+DT);
Tr = (T_tr-DT):0.1:(T_tr+DT);

subplot(4,2,1);
ycent=mean([H0(S(:,p1)',Tr_left(end)/1000), H0(S(:,p2)',Tr_left(end)/1000)]);
yrange=abs(H0(S(:,p1)',Tr_left(end)/1000)-ycent);
DeltaH0(1) = H0(S(:,p2)',T_tr/1000)-H0(S(:,p1)',T_tr/1000);

plot(Tr_left,H0(S(:,p1)',Tr_left/1000),'Color',c(p1,:));hold on;
plot(Tr_right,H0(S(:,p1)',Tr_right/1000),'Color',c(p1,:),'LineStyle',':');
plot(Tr_left,H0(S(:,p2)',Tr_left/1000),'Color',c(p2,:),'LineStyle',':');hold on;
plot(Tr_right,H0(S(:,p2)',Tr_right/1000),'Color',c(p2,:));
plot([T_tr T_tr],[0 600],'k:');
plot([T_tr T_tr],[H0(S(:,p1)',Tr_left(end)/1000), H0(S(:,p2)',Tr_left(end)/1000)],'m','LineWidth',2);
xlim([Tr(1) Tr(end)]);
ylim(ycent+[-yrange yrange]*1.5);
ylabel('H^0-H^0_2_9_8 [kJ/mol]');
xlabel('T [K]');
set(gca,'TickDir','out');

subplot(4,2,3);
ycent=mean([S0(S(:,p1)',Tr_left(end)/1000), S0(S(:,p2)',Tr_left(end)/1000)]);
yrange=abs(S0(S(:,p1)',Tr_left(end)/1000)-ycent);
DeltaS0(1) = S0(S(:,p2)',T_tr/1000)-S0(S(:,p1)',T_tr/1000);

plot(Tr_left,S0(S(:,p1)',Tr_left/1000),'Color',c(p1,:));hold on;
plot(Tr_right,S0(S(:,p1)',Tr_right/1000),'Color',c(p1,:),'LineStyle',':');
plot(Tr_left,S0(S(:,p2)',Tr_left/1000),'Color',c(p2,:),'LineStyle',':');hold on;
plot(Tr_right,S0(S(:,p2)',Tr_right/1000),'Color',c(p2,:));
plot([T_tr T_tr],[0 600],'k:');
plot([T_tr T_tr],[S0(S(:,p1)',Tr_left(end)/1000), S0(S(:,p2)',Tr_left(end)/1000)],'m','LineWidth',2);
xlim([Tr(1) Tr(end)]);
ylim(ycent+[-yrange yrange]*1.5);
ylabel('S^0 [J/(mol K)]');
xlabel('T [K]');
set(gca,'TickDir','out');

subplot(4,2,5);
plot(Tr_left,G0(S(:,p1)',Tr_left),'Color',c(p1,:));hold on;
plot(Tr_right,G0(S(:,p1)',Tr_right),'Color',c(p1,:),'LineStyle',':');
plot(Tr_left,G0(S(:,p2)',Tr_left),'Color',c(p2,:),'LineStyle',':');hold on;
plot(Tr_right,G0(S(:,p2)',Tr_right),'Color',c(p2,:));
plot([T_tr T_tr],[0 -600],'k:');
xlim([Tr(1) Tr(end)]);
ylim(G0(S(:,p1)',Tr_left(end))+[-.2 .2]);
ylabel('G^0 [kJ/mol]');
xlabel('T [K]');
set(gca,'TickDir','out');

subplot(4,2,7);
y1=G0(S(:,p2)',Tr(1))-G0(S(:,p1)',Tr(1));
y2=G0(S(:,p2)',Tr(end))-G0(S(:,p1)',Tr(end));

plot(Tr,G0(S(:,p2)',Tr)-G0(S(:,p1)',Tr),'b');hold on; % Note: could be replaced with DeltaG function
plot(Tr,-DeltaS0(1)/1000*(Tr-T_tr),'r:');

plot([Tr(1) Tr(end)],[0,0],'k:');
plot([T_tr T_tr],[y2 y1],'k:');
xlim([Tr(1) Tr(end)]);
ylim([y2 y1]);

ylabel('\DeltaG [kJ/mol]');
xlabel('T [K]');
set(gca,'TickDir','out');

% liguid-gas transition
DT = 2;
T_tr = T_Lg;
p1=2;
p2=3;

Tr_left = (T_tr-DT):0.1:T_tr;
Tr_right = T_tr:0.1:(T_tr+DT);
Tr = (T_tr-DT):0.1:(T_tr+DT);

subplot(4,2,2);
ycent=mean([H0(S(:,p1)',Tr_left(end)/1000), H0(S(:,p2)',Tr_left(end)/1000)]);
yrange=abs(H0(S(:,p1)',Tr_left(end)/1000)-ycent);
DeltaH0(2) = H0(S(:,p2)',T_tr/1000)-H0(S(:,p1)',T_tr/1000);

plot(Tr_left,H0(S(:,p1)',Tr_left/1000),'Color',c(p1,:));hold on;
plot(Tr_right,H0(S(:,p1)',Tr_right/1000),'Color',c(p1,:),'LineStyle',':');
plot(Tr_left,H0(S(:,p2)',Tr_left/1000),'Color',c(p2,:),'LineStyle',':');hold on;
plot(Tr_right,H0(S(:,p2)',Tr_right/1000),'Color',c(p2,:));
plot([T_tr T_tr],[0 600],'k:');
plot([T_tr T_tr],[H0(S(:,p1)',Tr_left(end)/1000), H0(S(:,p2)',Tr_left(end)/1000)],'m','LineWidth',2);
xlim([Tr(1) Tr(end)]);
ylim(ycent+[-yrange yrange]*1.5);
ylabel('H^0-H^0_2_9_8 [kJ/mol]');
xlabel('T [K]');
set(gca,'TickDir','out');

subplot(4,2,4);
ycent=mean([S0(S(:,p1)',Tr_left(end)/1000), S0(S(:,p2)',Tr_left(end)/1000)]);
yrange=abs(S0(S(:,p1)',Tr_left(end)/1000)-ycent);
DeltaS0(2) = S0(S(:,p2)',T_tr/1000)-S0(S(:,p1)',T_tr/1000);

plot(Tr_left,S0(S(:,p1)',Tr_left/1000),'Color',c(p1,:));hold on;
plot(Tr_right,S0(S(:,p1)',Tr_right/1000),'Color',c(p1,:),'LineStyle',':');
plot(Tr_left,S0(S(:,p2)',Tr_left/1000),'Color',c(p2,:),'LineStyle',':');hold on;
plot(Tr_right,S0(S(:,p2)',Tr_right/1000),'Color',c(p2,:));
plot([T_tr T_tr],[0 600],'k:');
plot([T_tr T_tr],[S0(S(:,p1)',Tr_left(end)/1000), S0(S(:,p2)',Tr_left(end)/1000)],'m','LineWidth',2);
xlim([Tr(1) Tr(end)]);
ylim(ycent+[-yrange yrange]*1.5);
ylabel('S^0 [J/(mol K)]');
xlabel('T [K]');
set(gca,'TickDir','out');

subplot(4,2,6);
plot(Tr_left,G0(S(:,p1)',Tr_left),'Color',c(p1,:));hold on;
plot(Tr_right,G0(S(:,p1)',Tr_right),'Color',c(p1,:),'LineStyle',':');
plot(Tr_left,G0(S(:,p2)',Tr_left),'Color',c(p2,:),'LineStyle',':');hold on;
plot(Tr_right,G0(S(:,p2)',Tr_right),'Color',c(p2,:));
plot([T_tr T_tr],[0 -600],'k:');
xlim([Tr(1) Tr(end)]);
ylim(G0(S(:,p1)',Tr_left(end))+[-.2 .2]);
ylabel('G^0 [kJ/mol]');
xlabel('T [K]');
set(gca,'TickDir','out');

subplot(4,2,8);
y1=G0(S(:,p2)',Tr(1))-G0(S(:,p1)',Tr(1));
y2=G0(S(:,p2)',Tr(end))-G0(S(:,p1)',Tr(end));

plot(Tr,G0(S(:,p2)',Tr)-G0(S(:,p1)',Tr),'b');hold on;
plot(Tr,-DeltaS0(2)/1000*(Tr-T_tr),'r:');

plot([Tr(1) Tr(end)],[0,0],'k:');
plot([T_tr T_tr],[y2 y1],'k:');
xlim([Tr(1) Tr(end)]);
ylim([y2 y1]);

ylabel('\DeltaG [kJ/mol]');
xlabel('T [K]');
set(gca,'TickDir','out');

%% Output data
clc;
fprintf('Phase transformation:\n');
fprintf('s->L:\n');
fprintf('T_tr    = %4.2f K.\n',T_sL);
fprintf('DeltaH0 = %4.2f kJ mol^-1.\n',DeltaH0(1));
fprintf('DeltaS0 = %4.2f J mol^-1 K^-1.\n\n',DeltaS0(1));

fprintf('L->gas:\n');
fprintf('T_tr    = %4.2f K.\n',T_Lg);
fprintf('DeltaH0 = %4.2f kJ mol^-1.\n',DeltaH0(2));
fprintf('DeltaS0 = %4.2f J mol^-1 K^-1.\n\n',DeltaS0(2));

%% How good is the linear approximation of DeltaG near T_tr?
% careful here: extrapolation beyond the range over which Shomate
% coefficients are defined may not represent reality very well, or at all.
h6=figure;
set(gcf,'DefaultLineLineWidth',1.5)
DT = 50;

T_tr = T_Lg;
p1=2;
p2=3;
i=2;
Tr = (T_tr-DT):0.1:(T_tr+DT);

subplot(223);
y1=G0(S(:,p2)',Tr(1))-G0(S(:,p1)',Tr(1));
y2=G0(S(:,p2)',Tr(end))-G0(S(:,p1)',Tr(end));

plot(Tr,G0(S(:,p2)',Tr)-G0(S(:,p1)',Tr),'b');hold on;
plot(Tr,-DeltaS0(i)/1000*(Tr-T_tr),'r:');
plot([Tr(1) Tr(end)],[0,0],'k:');
plot([T_tr T_tr],[y2 y1],'k:');
xlim([Tr(1) Tr(end)]);
ylim([y2 y1]);

ylabel('\DeltaG [kJ/mol]');
xlabel('T [K]');
set(gca,'TickDir','out');grid on;

subplot(224);
plot(Tr,100*(G0(S(:,p2)',Tr)-G0(S(:,p1)',Tr)+DeltaS0(i)/1000*(Tr-T_tr))./(-DeltaS0(i)/1000*(Tr-T_tr)),'b');hold on;
grid on;
%plot([T_tr T_tr],[y2 y1],'k:');
xlim([Tr(1) Tr(end)]);
%ylim([y2 y1]);
ylabel('(\DeltaG-\DeltaG_e_s_t)/\DeltaG_e_s_t [%]');
xlabel('T [K]');
set(gca,'TickDir','out');

%
T_tr = T_sL;
p1=1;
p2=2;
i=1;
Tr = (T_tr-DT):0.1:(T_tr+DT);
subplot(221);
y1=G0(S(:,p2)',Tr(1))-G0(S(:,p1)',Tr(1));
y2=G0(S(:,p2)',Tr(end))-G0(S(:,p1)',Tr(end));

plot(Tr,G0(S(:,p2)',Tr)-G0(S(:,p1)',Tr),'b');hold on;
plot(Tr,-DeltaS0(i)/1000*(Tr-T_tr),'r:');
plot([Tr(1) Tr(end)],[0,0],'k:');
plot([T_tr T_tr],[y2 y1],'k:');
xlim([Tr(1) Tr(end)]);
ylim([y2 y1]);

ylabel('\DeltaG [kJ/mol]');
xlabel('T [K]');
set(gca,'TickDir','out');grid on;

subplot(222);
plot(Tr,100*(G0(S(:,p2)',Tr)-G0(S(:,p1)',Tr)+DeltaS0(i)/1000*(Tr-T_tr))./(-DeltaS0(i)/1000*(Tr-T_tr)),'b');hold on;
grid on;
%plot([T_tr T_tr],[y2 y1],'k:');
xlim([Tr(1) Tr(end)]);
%ylim([y2 y1]);
ylabel('(\DeltaG-\DeltaG_e_s_t)/\DeltaG_e_s_t [%]');
xlabel('T [K]');
set(gca,'TickDir','out');

%% save figures h1=h6 as .eps files

print(h1,[filename,'Cp'],'-deps','-cmyk','-opengl');
print(h2,[filename,'S'],'-deps','-cmyk','-opengl');
print(h3,[filename,'H'],'-deps','-cmyk','-opengl');
print(h4,[filename,'G'],'-deps','-cmyk','-opengl');
print(h5,[filename,'PTs'],'-deps','-cmyk','-opengl');
print(h6,[filename,'Linearity'],'-deps','-cmyk','-opengl');

316-2 Problems

27.1 General

  1. Write a paragraph discussing the relevance of phase transformations in your daily life.

27.2 Laplace Pressure Derivation

  1. Derive the expression for the Laplace pressure inside a long cylinder of radius R .

27.3 Homogeneous Nucleation

  1. Consider the following data for nickel:
    Melting point
    1452 C
    Molar entropy of solid at T m
    56.07 J/K
    Molar entropy of liquid at T m
    66.27 J/K
    Solid density
    8.9 g/c m 3
    Molar mass
    58.7
    In their classic experiment Turnbull and Cech studied the undercooling of small droplets for a number of different metals . Assuming that nucleation in the droplets occurs homogeneously and using the data given below calculate the following at 1100 °C and 1200 °C:
    1. The molar volume of nickel.
    2. The work of nucleation ( W R * ).
    3. The dimensionless ratio, W R * / k B T .
    4. The radius of the critical nucleus.
    5. The pressure of the critical nucleus in pascals (assume the surrounding liquid is at atmospheric pressure).
    6. The molar enthalpy of melting at T m .
    7. Suppose a Ni droplet with a volume of about 100 μ m 3 is solidified. Approximate the temperature to which the droplet must be cooled in order for solidification to occur by homogeneous nucleation.
  2. Import the file labeled ElementData.mat that includes the required data for various elements on the periodic table into Matlab and:
    1. Derive the expressions for Δ P , R * , W R * , and W R * / k B T in terms of T m , Δ T , V m , Δ S f , and γ .
    1. Plot V m , Δ S f , γ , Δ P , R * , W R * , and W R * / k B T using Δ T=100K versus atomic number ( Z ) and label all axes including units and each data point with the chemical symbol corresponding to the element. Hint: You should only consider those elements for which the values of γ are included in the ElementData.mat file. For both W R * and W R * / k B T plot the y axis on a log scale. Also, in order to label the data points with the chemical symbol you will need to use the text(x, y, 'string') function. You may want to use subplots.
      The ElementData.mat file has the following format:
      ElementData = 
      Name: {118x1 cell}
      Symbol: {118x1 cell} 
      DeltaH0f: [118x1 double] 
      Tm: [118x1 double] 
      Z: [118x1 double] 
      Aw: [118x1 double] 
      rho: [118x1 double] 
      gamma: [118x1 double] 
      Vm: [118x1 double] 
      DeltaS0f: [118x1 double] 
      Structure: {118x1 cell} 
      Units: {1x8 cell}
    2. Discuss the plots from part (b) with respect to trends in the periodic table, which variables are really important, outliers, and rules of thumb i.e. typical range of values or average value. Does homogeneous nucleation ever really happen?
    3. Now replot the data for both Δ T=352K and Δ T=252K and compare the R * and W R * values obtained for Ni to those you calculated in question 3..
  3. Derive expressions for R * and W R * for a cuboidal nucleus.
  4. In the derivations for nucleation in this course we assume that the nucleus is incompressible. Show that this is a valid assumption for solidification of Ni with γ =2.38J/ m 2 and R * =1nm . Hint: Assume that the material is linearly elastic and isotropic. Therefore, you can calculate the bulk modulus using a simple relationship which is a function of Young's modulus and Poisson's ratio. Please cite your source for the values of E and ν that you use.

27.4 Surface and Interface Effects

  1. The surface free energy of solid gold at its melting point (1063ºC) is 1.400J/ m 2 . The surface energy of liquid gold at this temperature is 1.128J/ m 2 , and the interfacial energy for the gold solid/liquid interface is = 0.132J/ m 2 . The latent heat of fusion for gold is 1.2x1 0 9 J/ m 3 .
    1. What is the contact angle for liquid gold on a solid gold surface at 1063ºC ?
    1. Is there thermodynamic barrier for the melting of a gold surface?
    2. Suppose a thin liquid gold layer of thickness δ exists at the surface of gold at 1058 C (5 below the equilibrium melting point). By comparing to the free energy of a gold surface that does not have this liquid layer, estimate the maximum thickness of the liquid layer that will be thermodynamically stable at this temperature.
    3. Very small gold particles have melting points that differ from the melting point of bulk gold. From the analysis given above, do you expect the melting point of a particle with a diameter of 2 nm to be higher or lower than the melting point of bulk gold? Give a brief explanation for your answer.
  2. Suppose precipitates form at grain boundaries within the matrix phase, with geometries that look like the following:
    image: 122_home_ken_Mydocs_MSEcore_316-2_figures_drawing77.svg
    What is the ratio of the grain boundary free energy to the interfacial energy between the precipitate and the matrix phase?
  3. Water beads up on a freshly waxed car to form droplets with a contract angle of 80 . What is the interfacial free energy for the wax/water interface, if the surface energy of the wax is 0.025 J/ m 2 ? (Note: you'll need to look up the surface energy of water to do this problem).
  4. An oil droplet ( δ phase) is placed on the water surface (phase β ) in contact with air (phase α ). The schematic of the cross section of the droplet is as describe in class (and repeated below). The surface free energy of water (against air) is 0.072 J/ m 2 . If the measured values of θ 1 and θ 2 in the figure below are 37 and 23 , respectively, what are the values of the oil surface energy and the oil/water interfacial energy.
    image: 123_home_ken_Mydocs_MSEcore_316-2_figures_drawing3.svg

27.5 Heterogeneous Nucleation

  1. Derive the structure factor, S( θ ) .
  2. Suppose that nucleation of a solid, single component metal occurs heterogeneously at a wall. Based on the values given for Ni in problem 3., what contact angle for the critical nucleus must be obtained in order to increase the minimum temperature required for solidification by 50°C?

27.6 Nucleation in a Binary System

  1. 3. Consider the formation of a nucleus β * with composition X β * from metastable α with composition X 0 α .At temperature T , the composition of stable α is X eq α , that of stable β is X eq β (all X refer to X 1 ). In class we derived an expression for the molar Gibbs free energy of formation for the nucleus:
    Δ G m α β * = G m β ( X β * ) - G m α ( X 0 α ) - G α X | X 0 α ( X β * - X 0 α ) (27.1)
    Show that for X 0 α - X eq α 0 and X β * - X eq β 0 , Eq. (1) can be rewritten in the following form:
    Δ G m α β * =- δ 2 G α δ X 2 | X 0 α ( X 0 α - X eq α ) ( X eq β - X eq α )
    Hint: Express G m β ( X β * ) in terms of G m α . Approximate all terms at non-equilibrium compositions as Taylor expansions around suitable equilibrium values.
  2. In class we used the definition of the misfit parameter for a β nucleus in an α matrix as
    ϵ = 1 3 ( V m β - V m α V m α )
    i.e. one third of the volume strain. Show that for cubic systems, the misfit parameter can be approximated as
    ϵ cubic = a β - a α a α
    where a is the lattice parameter. Hint: Write Δ V in terms of ϵ cubic and look at the behavior as ϵ cubic 0 .
  3. A coherent precipitate nucleates much more easily than does an incoherent particle of the same precipitate. To illustrate this:
    1. What is the ratio of W R * for the two types of precipitate if γ coherent =30 ergs/c m 2 and γ incoherent =300 ergs/c m 2 ? Assume that the precipitate is unstrained.
    1. If the chemical driving force ( Δ G v ) is given by -50 Δ T/ T e cal/c m 3 , T e =1000 K , the misfit strain is 0.001 for the coherent precipitate and zero for the incoherent precipitate, at what Δ T are the W R * 's for the two equal? Assume a shear modulus of the matrix of 5.46x1 0 10 Pa and bulk modulus of the precipitate of 15x1 0 10 Pa .
    2. Repeat the previous calculation using a misfit strains of 0.01 and 0.1.
    3. If the number of nuclei formed per cubic centimeter per second is given by N=1 0 27 \ exp( - W R * /kT ) , what is the rate of coherent nucleation at Δ T=25K and 250K with a misfit of 0.01? What is it for incoherent nucleation at these same values of Δ T ?
  4. Consider the following Al-Cu phase diagram:
    image: 124_home_ken_Mydocs_MSEcore_316-2_figures_Al-Cu_phase_diagram.jpg
    Suppose that a dispersion of roughly spherical θ precipitates is formed at 300 °C. Estimate the precipitate radius for which Cu solubility in the α phase (the Al-rich phase) will be increased by 25% in comparison to a flat α / θ interface. Assume an interfacial free energy for the α / θ interface of 0.3 J/ m 2 and a molar volume for the α and β phases of 7 c m 3 .
  5. Consider the Co-Cu phase diagram shown below:
    image: 125_home_ken_Mydocs_MSEcore_316-2_figures_Co-Cu_phase_diagram.jpg
    1. Plot the equilibrium activity of Cobalt as a function of composition across the entire phase diagram at 900ºC.
    2. Suppose the interfacial free energy for the Cu/Co interface is 300 mJ/ m 2 . Develop an expression for r * , the critical radius for a cobalt precipitate, as function of the atomic % cobalt in the alloy.
    3. Calculate W r * for a Copper rich alloy at 900ºC with a cobalt composition that exceeds the equilibrium composition by a factor of 1.15.

27.7 Spinodal Decomposition

  1. A and B form a regular solution with a positive heat of mixing so that the A-B phase diagram contains a miscibility gap.
    1. Starting from G= X A G A + X B G B + varOmega X A X B +RT( X A ln X A + X B ln X B ) , derive an equation for d 2 G/d X B 2 , assuming G A = G B =0 .
    1. Use the above equation to calculate the temperature at the top of the miscibility gap T c in terms of varOmega .
    2. Using MATLAB plot the miscibility gap for this system.
    3. On the same diagram plot the chemical spinodal.
  2. For a homogeneous alloy of composition X 0 decomposes into two parts, one with composition X 0 + Δ X and the other with composition X 0 - Δ X , show that the total chemical free energy will change by an amount Δ G c given by
    Δ G c = 1 2 d 2 G d X 2 ( Δ X ) 2
    Hint: Express G( X 0 + Δ X ) and G( X 0 - Δ X ) as Taylor series.
  3. Describe the effect of each of the following, and briefly explain your answer.
    1. The effect of coherent strains on the characteristic wavelength of the two-phase structure formed by spinodal decomposition.
    2. The effect of a reduction of the surface free energy on the nucleation rate.
    3. The effect of a decrease in the contact angle of a precipitate on its heterogeneous nucleation rate.
    4. Can a diffusion coefficient ever be negative? If so, when is this the case?

27.8 Constitutional Undercooling and the 'Mushy Zone'

  1. In our classroom discussion of interface stability, we considered the case where impurities decrease the melting point. Suppose that the impurities increase the melting point, so that the phase diagram looks like this:
    image: 126_home_ken_Mydocs_MSEcore_316-2_figures_drawing85.svg
    Suppose sample with the composition indicated by the arrow is solidified, so that the front moves forward with a certain velocity.
    1. Sketch the behavior of the impurity concentration in the liquid phase just ahead of the solidification front. Reference any specific compositions to the corresponding compositions on the phase diagram.
    1. On a separate figure, sketch the liquidus temperature in the liquid phase just ahead of the solidification front. Reference