Microstructural Evolution in Elastically-stressed Solids: A Phase-field Simulation

Simulation of microstructures under different processing conditions is important for fine- tuning the processing window as well as to understand the mechanisms. Phase field simulation has gained importance for problems with diffuse interfaces. Since in this simulation, thermodynamic driving forces (chemical as well as non-chemical) and kinetic constraints have been naturally incorporated, it has the potential to simulate microstructures under different processing and service conditions. In this paper, DMRL’s initiatives on using phase field simulations to understand microstructural evolution in both the phase separating and precipitating model systems have been presented. The influence of misfit stresses on the morphology of microstructures has been described. Output from actual thermodynamic calculations can be combined with these simulations to study systems of technological importance.

Technologically important materials such as Ni-base superalloys exhibit essentially two-phase microstructures. In Ni-base superalloys used for high temperature applications such as aeroengine turbine blades, the second phase (γ’) is precipitated through a solid solid transformation from the parent matrix phase (γ) through a multi-step solutionising and multi-step ageing process1-3. While the matrix has a disordered fcc structure, the precipitate has an ordered fcc (L12) structure and more often than not the compositions of the two phases are different. This difference in composition and structural arrangement results in a difference in the lattice parameters of the two phases. However, to keep the continuity of the lattice intact (i.e., to maintain coherency), the interface is strained (with the stresses within elastic limit) and such systems are called elastically-stressed solids.

Coherency stresses have a direct role in deciding the morphology of the precipitate phase4-7. In a two-phase system, the morphology of the second phase (minority phase) is decided by the sum of interfacial energy and the elastic energy8,9. In systems with isotropic interfacial energy, if the precipitate phase is completely coherent with the matrix (i.e., no lattice parameter misfit), then spherical shaped precipitates are expected.

In various misfit-controlled Ni-X alloys, it was found that the alloys with lower misfit had spherical precipitates whereas the alloys with higher misfit had cuboidal or elongated precipitates10. While the interfacial energy (which is more or less isotropic) drives a spherical shape, the elastic energy drives elongated shape and the precipitates acquire a shape for which the total energy is a minimum. The shapes corresponding to minimum energy are equilibrium shapes and the rate at which this shape is achieved is affected by kinetics of the transformation.

Phase-field simulation11-20 has emerged as a successful tool in predicting the morphology of microstructures considering both thermodynamical driving force as well as kinetic constraints. Its application spans a wide area of materials simulation, ranging from solidification microstructure, grain growth in single-phase and multi-phase materials, precipitation kinetics, effect of magnetic field on microstructures, electromigration, etc. An overview by Chen19 and the references therein have enormous information. Since tracking of interfaces is not required in phase-field simulations, it is more attractive and convenient to deal with problems involving interfaces in contrast to sharp interface model where it is very difficult to track interfaces11. Moreover, in these simulations, no assumption is made with regards to the thickness of interfaces, and hence, these are well-suited for studying this class of problems21. This method has the potential to predict the microstructural changes during the entire thermal treatment starting from dendrites arising during solidification, solutionising, precipitation during high temperature ageing, coarsening of precipitates and changes in the microstructures during service, especially arising due to externally applied stresses at high temperatures and dislocation interactions.

This study demonstrates the initiative in establishing expertise in phase-field simulation of microstructural evolution during solid → solid transformations. The paper covers essential mathematical background of phase-field simulation including the micromechanics of defects (specifically addressing solid-state precipitates). The application of phase-field simulation to microstructural evolution in a spinodally decomposing (2-D as well as 3-D) systems and a precipitating (2-D model Ni-Al) system has been demonstrated. The effect of elastic stresses on the morphological evolution has been shown qualitatively.

In a phase-field simulation, the energy of the system is expressed in terms of appropriate field variables and their derivatives11-20. There are two types of field variables: conserved and non-conserved. Composition is a conserved variable since the overall composition in a closed system has to be constant, whereas grain orientations (in grain growth problems), order-parameters (in solidification or precipitation problems) are non-conserved variables. The spatial derivatives of these variables account for changes wrt adjacent domain, and hence, represent interfacial contribution. Apart from these, long-range contribution such as elastic energy, magnetic energy (which themselves can be dependent on the field variables) can be added to the overall energy of the system. In a typical phase-field problem, the energy is expressed as a sum of chemical and non-chemical energies19. For example, in systems such as Ni-base superalloys, the non-chemical free energy is the elastic energy (Fel), and hence, the total energy, F, is expressed as

$F={F}_{ch}+{F}_{el}$ (1)

where, Fch is the chemical free energy given as

${F}_{ch}=\mathrm{\int }\left[\begin{array}{l}f\left({c}_{1},{c}_{2},\cdots {c}_{n},{\text{η}}_{1},{\text{η}}_{2},\cdots ,{\text{η}}_{m}\right)\\ +\sum _{i=1}^{n}{\text{α}}_{i}{\left(\nabla {c}_{i}\right)}^{2}+\sum _{i=1}^{m}{\text{β}}_{i}{\left(\nabla {\text{η}}_{i}\right)}^{2}\end{array}\right]\text{\hspace{0.17em}}d\mathbit{r}$ (2)

In Eqn (2), f is the local free energy density expressed as a function of conserved variables, c1, c2, …, cn (compositions in a multi-component alloy) and non-conserved variables, η1, η2,… ηm (order parameters indicating whether the location is matrix or precipitate and if precipitate, to which of the different possible orientations it belongs to). αi and βi are respectively the gradient energy coefficients for the conserved and the non-conserved variables. The functional form of local free energy functional is the key that would differentiate between different phase-field models.

2.1 Field Evolution

In a phase-field simulation, the microstructural evolution is obtained by solving Cahn-Hilliard and Allen-Cahn equations for the conserved and non-conserved variables, respectively19,21,22. These equations describe the time evolution of the field variables in terms of variational derivatives of the free energy functional wrt the field variables. Mathematically, these are stated as

$\frac{\partial {c}_{i}}{\partial t}=M{\nabla }^{2}\frac{\delta F}{\delta {c}_{i}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}i=1\cdots n$ (3)

$\frac{\partial {\text{η}}_{i}}{\partial t}=-L\frac{\delta F}{\delta {\text{η}}_{i}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}i=1\cdots m$ (4)

It is to be noted that both ci and ηi are functions of position and time. The coefficients M and L are the kinetic coefficients which would decide how fast the microstructures would evolve. Though, in general, these coefficients can be functions of field variables as well as can have directional dependence, here these are considered as isotropic scalar values and independent of concentration and order parameters.

2.2 Numerical Solution

It is convenient to solve the field kinetic equations in Fourier space. Assuming periodic boundary conditions (i.e., if we solve a 2-D problem, the problem space is repeated infinitely in both x and y directions filling the entire 2-D space), Fourier transformation (wrt position coordinates) of Eqns (3) and (4) will convert these into a set of algebraic equations19, as :

$\frac{{\stackrel{˜}{c}}_{i}^{t+\mathrm{\Delta }t}-{\stackrel{˜}{c}}_{i}^{t}}{\mathrm{\Delta }t}=-M{k}^{2}\left({\stackrel{˜}{g}}^{t}+2{k}^{2}{\text{α}}_{i}{\stackrel{˜}{c}}_{i}^{t+\mathrm{\Delta }t}\right)$ (5)

$\frac{{\stackrel{˜}{\eta }}_{i}^{t+\mathrm{\Delta }t}-{\stackrel{˜}{\eta }}_{i}^{t}}{\mathrm{\Delta }t}=-L\left({\stackrel{˜}{h}}^{t}+2{k}^{2}\beta {\stackrel{˜}{\eta }}_{\text{i}}^{t+\mathrm{\Delta }t}\right)$ (6)

In Eqns (5) and (6), the ~ above the quantities represent Fourier transforms, g and h are the derivatives of f wrt c and η, respectively, and k is the length of Fourier vector. The superscripts, t and tt represent the time at which the quantities are evaluated. In the above formalism, semi-implicit method has been used, wherein c and η occurring on the RHS of the equations are taken at time tt, whereas the nonlinear terms g and h are taken at time t.

2.3 Free Energy Functional

Two types of free energy functionals have been considered: (i) a spinodally decomposing system, and (ii) precipitating systems. For a binary spinodally decomposing system, simplest form of free energy is a double-well potential21. The double-well potential is mathematically represented as

$f={c}^{2}{\left(1-c\right)}^{2}$ (7)

where c is the concentration of the solute phase. This potential mimics a system with free energy minima at concentrations c = 0 and c = 1, as shown in Fig. 1. The physical implication is that in an A-B type of phase separating system, the two phases will have the compositions of that of pure A and pure B.

In case of precipitating system such as γ-γ’, the free energy functional is much more complicated, and involves a set of non-conserved variables (η1, η2, η3) apart from the conserved variable, c. A typical functional is given as20

$\begin{array}{l}f\left(c,{\text{η}}_{1},{\text{η}}_{2},{\text{η}}_{3}\right)={A}_{1}{\left(c-{c}_{1}\right)}^{2}+{A}_{2}\left({c}_{2}-c\right)+\left({\text{η}}_{1}^{2}+{\text{η}}_{2}^{2}+{\text{η}}_{3}^{2}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-{A}_{3}{\text{η}}_{1}{\text{η}}_{2}{\text{η}}_{3}+{A}_{4}\left({\text{η}}_{1}^{4}+{\text{η}}_{2}^{4}+{\text{η}}_{3}^{4}\right)\end{array}$ (8)

where A1, A2, A3, A4, c1, c2 are constants for a particular system of interest. In this free energy functional, the disorderd γ region is represented by the order parameter set values (η1, η2, η3) = (0, 0, 0), whereas the four different variants of the ordered γ’ phase are represented by (1, 1, 1) η0, (1, -1, -1) η0, (-1, 1, -1) η0, (-1, -1, 1) η0, where η0 is the equilibrium value of the order parameter. Considering A1, A2, A3, A4, c1, c2 as 40.0, 17.0, 46.8, 15.0, 5.0, 0.05, 0.22 respectively and M = 4.0 and L = 4.0, the calculated free energy20 curve is shown in Fig. 2.

2.4 Elastic Energy

As noted earlier, two-phase Ni-base superalloys are under the influence of elastic stresses arising out of coherency strains (lattice misfit). The total elastic energy of the stressed system is given by22

${F}_{el}=\frac{1}{2}\underset{V}{\int }{\sigma }_{ij}^{el}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\epsilon }_{ij}^{el}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}d\mathbit{r}$ (9)

where, ${\epsilon }_{ij}^{el}\left(\mathbit{r}\right)={\epsilon }_{ij}\left(\mathbit{r}\right)-{\epsilon }_{ij}^{0}\left(\mathbit{r}\right)$ is the elastic strain (difference between the total strain and the misfit strain). The misfit strain itself is a function of composition and it can be assumed to vary linearly with composition (Vegard’s law). i.e., ${\epsilon }_{ij}^{0}\left(\mathbit{r}\right)={\epsilon }_{ij}^{0}\text{\hspace{0.17em}}\text{\hspace{0.17em}}c\left(\mathbit{r}\right)$ . Assuming homogeneous elasticity (i.e., elastic constants are independent of composition, and hence, position), the elastic stress is given by ${\sigma }_{ij}^{0}\left(\mathbit{r}\right)={C}_{ijkl}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\epsilon }_{kl}^{0}\left(\mathbit{r}\right)$ , where Cijkl are the elastic constants of the homogeneous system. In all these expressions, summation over repeated indices (where the indices vary from 1 to 3) is understood without explicitly stating. Though the Eqn (9) appears compact, evaluation of elastic strain is not straightforward. Considering, mechanical equilibrium of the system (Appendix 1), the final result (which is easily programmable) is given as20,22

${F}_{el}=\frac{1}{2}\underset{V}{\int }B\left(n\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{|\stackrel{˜}{c}\left(\mathbit{r}\right)|}^{2}\frac{dk}{{\left(2\text{π}\right)}^{3}}$ (10)

where, $B\left(n\right)={C}_{ijkl}{\epsilon }_{ij}^{0}{\epsilon }_{kl}^{0}-{n}_{j}{\sigma }_{ij}^{0}{\text{Ω}}_{ki}\left(n\right){\sigma }_{kl}^{0}{n}_{l}$ and ${\text{Ω}}_{ik}^{\text{-}1}\left(n\right)={C}_{ijkl}{n}_{j}{n}_{l}$ is inverse Green function of anisotropic elasticity. Here, $n=k/|k|$ is the unit vector in Fourier space with k being the reciprocal vector. $\stackrel{˜}{c}\left(k\right)$ is the Fourier transform of the concentration field given as

$\stackrel{˜}{c}\left(k\right)=\underset{V}{\int }c\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{e}^{-ik\cdot \mathbit{r}}\text{\hspace{0.17em}}d\mathbit{r}$

The variational derivative of the elastic energy wrt to composition (all quantities in Fourier space) is given as $B\left(\mathbit{n}\right)\text{\hspace{0.17em}}B\left(\mathbit{n}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{˜}{c}\left(\mathbit{k}\right)$ . Incorporating the elastic energy, the semi-implicit equation for compositional evolution Eqn (5) is modified as

$\frac{{\stackrel{˜}{c}}_{i}^{t+\mathrm{\Delta }t}-{\stackrel{˜}{c}}_{i}^{t}}{\mathrm{\Delta }t}=-M{k}^{2}\left({\stackrel{˜}{g}}^{t}+\left[2{k}^{2}{\text{α}}_{i}+B\left(n\right)\right]{\stackrel{˜}{c}}_{i}^{t+\mathrm{\Delta }t}\right)$ (16)

The equation for order parameter evolution remains unchanged. Incorporating the elastic energy and the free energy functional of the precipitating system, a program in C has been developed by the author.

Here, application of the developed code has been demonstrated through simulations of 2-D and 3-D spinodally decomposing systems (using the free energy in Eqn (9)) and of 2-D precipitating systems (using the free energy in Eqn (10)). For the 2-D simulations of the spinodal decomposition, binary A-B systems with 0.25B and 0.5B have been considered. The effect of elastic (both isotropic and anisotropic) stresses on the microstructural evolution has been dealt with. For the 3-D spinodal decomposition, binary A-B system with 0.3B and 0.5B, both in the absence and presence of elastic stresses have been considered. Precipitating system has been simulated for 2-D model Ni-Al systems with 15 per cent and 18 per cent Al, in the presence of elastic stresses. These demonstrations serve a pedagogical role in understanding the individual effect of different factors that affect microstructural morphology.

3.1 Microstructural Evolution in Phase Separating System

3.1.1 Two Dimensional Systems

The microstructural evolution in a spinodally decomposing system as a function of time has been studied in the absence as well as presence of misfit stresses. The system size in each case is 1024 x 1024. Starting configurations for these simulations are generated by introducing a random fluctuation (to a maximum extent of 1 per cent of desired composition) at every point of the simulation system. This has been done to mimic the compositional fluctuations that trigger spinodal decomposition. Microstructural evolution has been tracked using free energy functional in Eqn (9) and using Eqn (8) (evolution of composition field as a function of time).

System with 25 per cent B: Figure 3 shows the microstructural evolution for a binary A-B alloy of composition c = 0.25B as a function of time. The gray areas represent regions that are richer in B. The segregation of second phase starts appearing around t = 1000 [Fig. 3(a)]. Well developed particles can be seen in Fig. 3(b) corresponding to a time of t = 2000. As the system evolves further, coarsening as well as coalescence of the second phase particles is evident from the decrease in the number density and increase in the size of particles. It can be seen that most of the particles have circular morphology due to the isotropic nature of the interfacial energy. Some particles are seen to have an elongated morphology, which is due to coalescence; eventually these particles also acquire circular morphology through redistribution of the solute atoms.

System with 50 per cent B: Figure 4 shows the microstructural evolution as a function of time in a system with composition c = 0.5B. Even from early times of evolution, the microstructure is an interconnected structure of A-rich (dark) and B-rich (bright) phases. The microstructure evolves in a self-similar fashion with increasing time. As time progresses, the phases can be seen to coarsen. Isolated particles with convex shapes are seen to disappear to the benefit of growth of lobe-like features of adjacent regions due to curvature effect or Ostwald ripening24. The chemical potential of solute around a particle scales as inverse of its radius of curvature. Hence, the chemical potential of solute around a smaller particle is higher compared to that around a bigger particle. As diffusion occurs down the chemical potential gradient, solute atoms diffuse from smaller particles to their bigger counterparts. This leads to eventual disappearance of smaller particles and coarsening of the larger ones.

Effect of isotropic elastic stress: Figures 5 and 6 show the evolution of microstructures in systems with composition 0.25B and 0.5B, respectively, but now with isotropic elastic stress arising due to lattice parameter misfit strain of 0.01 (i.e., lattice parameter of B differs from that of A by 1 per cent). The shear modulus and Poisson’s ratio considered are 400 (non-dimensional units) and 0.3, respectively. Comparing Fig. 3 (corresponding to t = 2000) with Fig. 5 (corresponding to t = 4000), it is clear that the nucleation is delayed in the later case. In contrast, for a composition of c = 0.5 B, there is no appreciable delay in the initiation of phase separation. This is easy to understand as this composition lies well within the coherent spinodal whereas a composition of c = 0.25B is presumably near the edge of coherent spinal (the edge of chemical spinodal being at c = 0.21B). It is to be noted that the phase boundaries of coherent spinodal (due to coherency stresses) occur within the chemical spinodal, and hence, the coexisting two-phase filed is narrowed24. The incorporation of elastic energy provides a barrier to the spontaneously evolving system, and hence, demands more fluctuations to induce phase separation. However, once the phase separation sets in, further growth is a spontaneous process. Accordingly, during late stages of evolution, there is no appreciable difference in the morphologies of microstructures of systems with and without elastic stresses.

Effect of anisotropic elastic stress: Figures 7 and 8 show the microstructural evolution in systems with composition 0.25B and 0.5B, respectively, in the presence of cubic elastic misfit (effective shear modulus of 400, effective Poisson’s ratio of 0.3 and Zener anisotropy factor of 3.0 were chosen). From Fig. 7, the effect of elastic stress is evident even during the early stages of phase separation. Since elastic stresses are of long range in nature, the second phase particles are aligned along the x- and y- axes (which are softer directions for the particular set of elastic constants chosen). During late stages, bigger second phase particles are elongated along the x- and y- axes. In the case of precipitates with cubic misfit, elastic energy is minimum for needle shaped particles. When the interfacial energy is isotropic, it is minimum for circular shapes. In the presence of both elastic and interfacial energies, for small sizes, interfacial energy is dominant, and hence, the particles are circular in nature, and for large sizes, elastic energy contribution is more, and hence, the shapes are elongated with edges along the x and y axes, with curved corners. In the case of c = 0.5B, the interconnected microstructures are elongated along the x- and y- axes.

3.1.2 Three-dimensional Systems

Three dimensional (3-D) microstructural evolution in binary phase separating systems has been presented in Figs 9-12, for different cases similar to those presented for 2-D systems. OpenDX program25 has been used to render these plots. A system size of 128 x 128 x 128 has been chosen, since study of higher sizes requires more memory and is time-consuming as well. However, it should be noted that for a qualitative study, this system size is adequate to discern different microstructural features. For studying the influence of coherency stresses, similar parameters as in 2-D cases have been chosen, except that the effective shear modulus is considered as 200. In the concentration maps, red corresponds to pure B and blue corresponds to 50 per cent B. Isoconcentration maps corresponding to 50 per cent B have also been shown to clearly delineate boundaries of solute-rich regions (in these maps, green and grey correspond to outer and inner surfaces, respectively). The results are qualitatively similar to those explained in the two dimensional cases. The complexity of interconnected microstructures in 3-D cases is quite evident from these figures. Two-dimensional sections of these microstructures can however produce artifacts such as circular and elongated morphologies depending on the orientation of the cutting plane. In other words, it is important to use 3-D visualisation techniques for analysing 3-D simulation results.

3.2 Microstructures in Model Ni-Al System

Microstructural evolution in a model Ni-Al system has been studied using free energy expression given in Eqn (10). While in the case of phase separating system, microstructural evolution is tracked using Eqn. (7) (corresponding to composition), in this case, equations corresponding to order parameters [(Eqn. (8)] have also been solved. In this system, since there is a barrier for forming nuclei of second phase, sustained noise in composition as well as order parameters (mimicking fluctuations) has been introduced till some incubation time. Care has been taken to keep the overall composition at the initially desired level, since it is a conserved variable. For times greater than the incubation time, the incorporation of noise has been stopped. In some literature, noise is incorporated adhering to fluctuation-dissipation theorem26.

In Fig. 13, microstructural evolution corresponding to a composition of 15 per cent Al has been shown as a function of time. The incubation time in this case been chosen as 10 units. After t = 15, the precipitates (γ’) are well isolated and mostly have a circular morphology, due to the dominant nature of (isotropic) interfacial energy. Solute depletion in matrix near the precipitate-matrix interface is evident. After t = 50, the precipitates have square morphology, driven by the elastic energy arising due to misfit. Further, alignment and elongation of precipitates, both driven by long-range elastic energy, can be seen after t = 400. The four different colours of the precipitates correspond to the four possible variants of γ’ precipitates. Comparing the figures corresponding to t = 50 with t = 400, it can be seen that two particles (at top left corner) have coalesced to form a single particle. Since these two particles belong to the same variant, there is no antiphase boundary separating them, and hence, coalescence has taken place. Other particles, even if they are spatially located nearby, cannot undergo coalescence because they belong to different variants. Such particles can grow only through coarsening.

Microstructural evolution in model Ni-Al system containing 18 per cent Al is shown in Fig. 14. Here, because of large supersaturation of solute Al atoms, it was suffi cient to provide an incubation period of 2 time units, at which point of time, enough nuclei had formed. The circular precipitates at early stages (t = 7) grow in size by acquiring solute atoms from the matrix. The coalescence of precipitates belonging to the same variant leads to precipitates having non-circular morphology (Figs. 14(b) and 14(c)). Towards the late stage (t > 500), most of the precipitates have straight edges parallel to the x- and y- axes in accordance with the cubic anisotropic condition. Further, due to the long-range nature of elastic energy, the precipitates are aligned along the x- and y- axes. The microstructure mimics a typical fully aged microstructure of nickel base superalloy, wherein cuboidal precipitates are embedded in a disordered matrix.

Microstructural simulation using phase-field method is fast gaining importance. In this paper, DMRL’s initiatives in using phase-field simulations have been described. Detailed description of phase-field simulation methodology has been given. Results of case studies concerning phase separating system and precipitating system have been presented. Influence of elastic stresses arising due to difference in lattice parameter between two phases has been shown. Delay in initiation of phase separation due to the barrier imposed by elastic stresses has been shown. With the incorporation of elastic anisotropy, alignment and elongation of precipitates have been observed. Simulated microstructures for model Ni-Al system resemble experimentally observed γ-γ’ microstructures.

This study forms a part of DMRL’s effort on multiscale simulation of microstructural evolution. This will be further taken forward by incorporating input from thermodynamic calculations of real systems (free energies calculated from ThermoCalc software which is available in DMRL will be used for this). The presently developed code can handle only homogeneous elasticity problem (i.e., elastic constants of the precipitate and matrix are equal). Extension of this code to handle elastic inhomogeneity and effect of externally applied stresses is underway. The present code for two-dimensional systems for the precipitating system will be extended to three-dimensional systems.

The author thanks Director, DMRL for allowing him to carry out this work and Dr M. Vijayakumar, Sc G, for his support and encouragement. The author also thanks Prof T.A. Abinandanan, IISc, Bengaluru, for providing the kernel of the phase field simulation code for simulating microstructural evolution of incoherent miscibility gap system and for detailed discussions on the fundamentals of micromechanics of defects.

1. Das, N.; Singh, S.; Hazari, N.; Chatterjee, D. & Praveen, V.V.N.S.S.C. Indigenous cast superalloys and investment casting technology for gas turbine components. Metals, Mater. Proce., 2007, 19, 189-202.

2. Reed, Roger C. The superalloys fundamentals and applications. Cambridge University Press, Cambridge, 2006.

3. Durant-Charre, Madeline. The microstructure of superalloys. Gordon Breach Science Publishers, Amsterdam, 1997.

4. Johnson, W.C. Influence of elastic stress on phase transformations, In Lectures on the theory of phase transformations, edited by H.I. Aaronson. TMS, Warrendale, Pennsylvania, 1999. pp. 35-134.

5. Doi, M. Elasticity effects on the microstructure of alloys containing coherent precipitates. Prog. Mater. Sci., 1996, 40, 79-180.

6. Fratzl, P; Penrose, O & Lebowitz, J.L. Modelling of phase separation in alloys with coherent elastic misfit. J. Stat. Phy., 1999, 95(5/6), 1429-503.

7. Voorhees, P.W. & Johnson, W.C. The thermodynamics of elastically stressed crystals. In Solid State Physics: Advances in Research and Applications. Vol 59. Elsevier Academic Press, 2004. pp. 1-201.

8. Johnson, W.C. & Cahn, J.W. Elastically induced shape bifurcations of inclusions. Acta Metallurgica, 1984. 32, 1925-933.

9. Sankarasubramanian, R.; Jog, C.S. & Abinandanan, T.A. Symmetry-breaking transitions in equilibrium shapes of coherent precipitates: Effect of elastic anisotropy and inhomogeneity. Metall. Mater. Trans. A, 2002, 33A, 1083-1090.

10. Li, F.; Prikhodko, S.V.; Ardell, A.J. & Kim, D. Morphological evolution of γ’-type particles in Ni-base alloys: shape characterisation. In Twelth Proceedings of the International Conference on Solid-Solid Transformations. The Japan Institute of Metals, 1999. pp. 545-52.

11. Leo, P.H.; Lowengrub, J.S. & Jou, H.J. A diffuse interface model for microstructural evolution in elastically stressed solids. Acta Materialia, 1998, 46, 2113-130.

12. Zhu, J.; Chen, L.Q.; & Shen, J. Morphological evolution during phase separation and coarsening with strong inhomogeneous elasticity. Modell. Simul. Mat. Sci. Engg., 2001, 9, 499-511.

13. Hu, S.Y.; & Chen, L.Q.; A phase-field model for evolving microstructures with strong elastic inhomogeneity. Acta Materialia, 2001, 49, 1879-890.

14. Jin, Y.M.; & Wang, Y.U. & Khachaturyan, A.G. Three-dimensional phase field microelasticity theory and modeling of multiple cracks and voids. App. Phy. Lett., 2001, 79, 3071-073.

15. Wang, Y.U.; Jin, Y.M. & Khachaturyan, A.G. Phase field microelasticity theory and simulation of multiple voids and cracks in single crystals and polycrystals under applied stress. J. Appl. Phy., 2001, 91, 6435-451.

16. Wang, Y.U.; Jin, Y.M. & Khachaturyan, A.G. Three-dimensional phase field microelasticity theory of a complex elastically inhomogeneous solid. Appl. Phy. Lett., 2002, 80, 4513-515.

17. Wang, Y.U.; Jin, Y.M. & Khachaturyan, A.G. Phase field microelasticity theory and modeling of elastically and structurally inhomogeneous solid. J. Appl. Phy., 2002, 92, 1351-360.

18. Gururajan, M.P. Elastic inhomogeneity effects on microstructures – A phase field study. Indian Institute of Science, Bengaluru, 2006. PhD Thesis.

19. Chen, L.Q. Phase field models for microstructure evolution. Ann. Rev. Mater. Res., 2002, 32, 113-40.

20. Li, D.Y. & Chen, L.Q., Shape evolution and splitting of coherent particles under applied stresses. Acta Materialia, 1999, 47, 247-57.

21. Cahn, J.W. & Hilliard, J.E. Free energy of a nonuniform system I, Interfacial free energy. J. Chem. Phy., 1958, 28, 258-67.

22. Allen, S.M. & Cahn, J.W. Mechanisms of phase transformations within the miscibility gap of Fe-rich Fe-Al alloys. Acta Metallurgica, 1976, 24, 425-37.

23. Khachaturyan, A.G. Theory of structural transformations in solids. John Wiley & Sons, Chicago, 1983.

24. Porter, D.A., & Easterling, K.E. Phase transformations in metals and alloys. Nelson Thornes Ltd, Cheltenham, UK, 2001.

25. http://www.opendx.org/ (Assessed on in 2007)

26. Lifshitz, E.M. & Pitaevskii, L.P. Statistical physics, Pt-I: Landau and Lifshitz course on theoretical physics. Pergamon Press, Oxford, 1980, 5.

 Dr R. Sankarasubramanian completed his integrated ME and PhD from the Department of Metallurgy, Indian Institute of Science, Bengaluru, in 1994 and 2000, respectively. Presently, he is working as Scientist D at the Defence Metallurgical Research Laboratory (DMRL), Hyderabadand his research interest is multi-scale materials modelling and simulation.

# Appendix-1

To find out the elastic energy of coherent system, we start from the equation of mechanical equilibrium.

${\sigma }_{ij,j}\left(\mathbit{r}\right)=0,$ (A.1)

where, σij (r) are stress components and the subscript, j represents derivative of the quantity along the direction rj. Here, Einstein’s summation notation over repeated indices (all indices ranging from 1 to 3) has been followed. Assuming, the coherency strains to be in the linear elastic regime, the stress is related to strain through constitutive relation as

${\sigma }_{ij}={C}_{ijkl}\left({u}_{k,l}\left(\mathbit{r}\right)-{\epsilon }_{kl}^{0}\left(\mathbit{r}\right)\right),$ (A.2)

where, Cijkl are the elastic constants, uk,l (r) are the displacement gradients and ${\text{ε}}_{kl}^{0}\left(\mathbit{r}\right)$ are the components of the transformation strain or eigen strain. Following Vegard’s law, the transformation strain can be assumed to scale linearly with compostion, c(r), and is given as

${\text{ε}}_{kl}^{0}\left(\mathbit{r}\right)={\text{ε}}_{kl}^{0}\text{\hspace{0.17em}}\text{\hspace{0.17em}}c\left(\mathbit{r}\right)$ (A.3)

Substituting Eqn (2) in Eqn (1),

$\begin{array}{l}{C}_{ijkl}{\left({u}_{k,l}\left(\mathbit{r}\right)-{\text{ε}}_{kl}^{0}\left(\mathbit{r}\right)\right)}_{,j}=0\\ \text{or}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{C}_{ijkl}{u}_{k,lj}\left(\mathbit{r}\right)={C}_{ijkl}{\text{ε}}_{kl,j}^{0}\left(\mathbit{r}\right)\end{array}$ (A.4)

Using Fourier transform of the displacement field, the above equation can be simplified. Fourier transform is given by

$F\left({u}_{i}\left(\mathbit{r}\right)\right)={\stackrel{˜}{u}}_{i}\left(k\right)=\underset{-\infty }{\overset{\infty }{\int }}{u}_{i}\left(\mathbit{r}\right){e}^{-ik\cdot \mathbit{r}}d\mathbit{r}$ (A.5)

where, k is the Fourier vector and the ~ indicates Fourier transform.

It is easy to see that

$F\left({u}_{k,lj}\left(\mathbit{r}\right)\right)=\left(i{k}_{l}\right)\left(i{k}_{j}\right){\stackrel{˜}{u}}_{k}\left(k\right)=-{k}_{l}{k}_{j}{\stackrel{˜}{u}}_{k}\left(k\right)$ (A.6)

$\mathrm{and} F\left({\epsilon }_{kl,j}^{0}\left(\mathbit{r}\right)\right)=i{k}_{j}{\stackrel{˜}{\epsilon }}_{kl}^{0}\left(k\right)=i{k}_{j}{\epsilon }_{kl}^{0}\stackrel{˜}{c}\left(k\right)$ (A.7)

where, c(k) is the Fourier transform of the composition c(r). Applying Fourier transform to equation (A.4) and using equations (A.6) and (A.7),

$-{k}_{l}{k}_{j}{C}_{ijkl}{\stackrel{˜}{u}}_{k}\left(k\right)=i{k}_{j}{C}_{ijkl}\text{\hspace{0.17em}}{\epsilon }_{kl}^{0}\stackrel{˜}{c}\left(k\right)$

Writing in terms of unit vector in Fourier space ${n}_{i}={k}_{i}/|k|$

$-{|k|}^{2}{n}_{l}{n}_{j}{C}_{ijkl}{\stackrel{˜}{u}}_{k}\left(k\right)=i{n}_{j}|k|{C}_{ijkl}{\text{ε}}_{kl}^{0}\stackrel{˜}{c}\left(\mathbit{r}\right)$ (A.8)

Writing ${n}_{l}{n}_{j}{C}_{ijkl}={\text{Ω}}_{ik}^{-1}\left(n\right)$ and ${C}_{ijkl}{\text{ε}}_{kl}^{0}={\text{σ}}_{kl}^{0}$ , the displacement field can be written from Eqn (A.8) as:

${\stackrel{˜}{u}}_{k}\left(k\right)=-i{n}_{j}{\sigma }_{ij}^{0}{\text{Ω}}_{ik}\left(n\right)\stackrel{˜}{c}\left(k\right){|k|}^{-1}$ (A.9)

where, Ωik (n) is the inverse of ${\text{Ω}}_{ik}^{-1}\left(n\right)$ . Displacement gradient can be obtained by taking the derivative of Eqn (A.9) as:

${\stackrel{˜}{u}}_{k,l}\left(k\right)={\stackrel{˜}{u}}_{k}\left(k\right)i{k}_{l}={n}_{j}{n}_{l}{\text{σ}}_{ij}^{0}{\text{Ω}}_{ik}\left(n\right)\stackrel{˜}{c}\left(k\right)$ (A.10)

Strain is given by

${\text{ε}}_{ij}=\frac{1}{2}\left({u}_{i,j}+{u}_{j,i}\right)$ (A.11)

The elastic energy, Fel, is given by,

${F}_{el}=\frac{1}{2}\underset{V}{\int }{C}_{ijkl}\left({\text{ε}}_{ij}\left(\mathbit{r}\right)-{\text{ε}}_{ij}^{0}\left(\mathbit{r}\right)\right)\left({\text{ε}}_{kl}\left(\mathbit{r}\right)-{\text{ε}}_{kl}^{0}\left(\mathbit{r}\right)\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}d\mathbit{r}$ (A.12)

where, the integration is performed over the volume, V, of the system. Exploiting the symmetries of the components Cijkl,

$\begin{array}{l}{F}_{el}=\frac{1}{2}\underset{V}{\int }{C}_{ijkl}{\text{ε}}_{ij}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\text{ε}}_{kl}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}d\mathbit{r}-\underset{V}{\int }{C}_{ijkl}{\text{ε}}_{ij}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\text{ε}}_{kl}^{0}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}d\mathbit{r}\text{\hspace{0.17em}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}\underset{V}{\int }{C}_{ijkl}{\text{ε}}_{ij}^{0}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\text{ε}}_{kl}^{0}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}d\mathbit{r}\end{array}$ (A.13)

Now, we can use Parseval’s theorem which is an identity relating integrals over real space to that over Fourier space.

$\int f\left(\mathbit{r}\right){g}^{*}\left(\mathbit{r}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}d\mathbit{r}=\frac{1}{{\left(2\text{π)}}^{\text{3}}}\int \stackrel{˜}{f}\left(k\right)\text{\hspace{0.17em}}{\stackrel{˜}{g}}^{*}\left(k\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}dk$ (A.14)

where, the superscript * represents complex conjugate. Applying Parseval’s theorem to Eqn (A.13) and making use of the fact that the strains in that equation are real valued functions (and therefore the complex conjugates are the strains themselves), one gets

$\begin{array}{l}{F}_{el}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\stackrel{˜}{\epsilon }}_{ij}\left(k\right){\left[{\stackrel{˜}{\epsilon }}_{kl}\left(k\right)\right]}^{*}\frac{dk}{{\left(2\pi \right)}^{3}}-\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\stackrel{˜}{\epsilon }}_{ij}\left(k\right){\left[{\stackrel{˜}{\epsilon }}_{kl}^{0}\left(k\right)\right]}^{*}\frac{dk}{{\left(2\pi \right)}^{3}}\\ +\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\stackrel{˜}{\epsilon }}_{ij}^{0}\left(k\right){\left[{\stackrel{˜}{\epsilon }}_{kl}\left(k\right)\right]}^{*}\frac{dk}{{\left(2\pi \right)}^{3}}\end{array}$ (A.15)

or Fel = I1 + I2 + I3. We will show simplification of each of the above terms.

${I}_{1}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\stackrel{˜}{\epsilon }}_{ij}\left(k\right) {\stackrel{˜}{\epsilon }}_{kl}^{*}\left(k\right)\frac{dk}{{\left(2\pi \right)}^{3}}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\stackrel{˜}{u}}_{i,j}\left(k\right){\stackrel{˜}{u}}_{k,l}^{*}\left(k\right)$ (A.16)

While writing the above, the symmetry property that Cijkl = Cklij, has been exploited. Substituting the expression for ũi, j from Eqn (A.10) gives

${I}_{1}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{n}_{j}{\text{σ}}_{mn}^{0}{n}_{n}{\text{Ω}}_{mi}\left(n\right)\stackrel{˜}{c}\left(k\right){n}_{l}{\sigma }_{pq}^{0}{n}_{q}{\text{Ω}}_{pk}\left(n\right){\stackrel{˜}{c}}^{*}\left(k\right)\frac{dk}{{\left(2\pi \right)}^{3}}$ (A.17)

Using the relation ${C}_{ijkl}{n}_{j}{n}_{l}={\text{Ω}}_{ik}^{-1}\left(n\right)$ and the identity ${\text{Ω}}_{ik}^{-1}\left(n\right){\text{Ω}}_{mi}\left(n\right)={\delta }_{km}$ , Eqn (A.17) can be simplified as

$\begin{array}{l}{I}_{1}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{\delta }_{km}{\text{σ}}_{mn}^{0}{n}_{n}{\text{σ}}_{pq}^{0}{n}_{q}{\text{Ω}}_{pk}\left(n\right){|\stackrel{˜}{c}\left(k\right)|}^{2}\frac{dk}{{\left(2\pi \right)}^{3}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{n}_{n}{\text{σ}}_{kn}^{0}{\text{Ω}}_{pk}\left(n\right){\text{σ}}_{pq}^{0}{n}_{q}{|\stackrel{˜}{c}\left(k\right)|}^{2}\frac{dk}{{\left(2\pi \right)}^{3}}\end{array}$

Changing the dummy (or repeated) indices,

$\begin{array}{l}{I}_{1}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{n}_{j}{\sigma }_{ij}^{0}{\text{Ω}}_{ki}\left(n\right){\sigma }_{kl}^{0}{n}_{l}{|\stackrel{˜}{c}\left(k\right)|}^{2}\frac{dk}{{\left(2\pi \right)}^{3}}\\ {I}_{2} =-\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\stackrel{˜}{\epsilon }}_{ij}\left(k\right){\left[{\stackrel{˜}{\epsilon }}_{kl}^{0}\left(k\right)\right]}^{*}\frac{dk}{{\left(2\pi \right)}^{3}}=\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\stackrel{˜}{u}}_{i,j}\left(k\right){\text{ε}}_{kl}^{0}{\stackrel{˜}{c}}^{*}\left(k\right)\frac{dk}{{\left(2\pi \right)}^{3}}\\ =-\underset{-\infty }{\overset{\infty }{\int }}{\text{σ}}_{ij}^{0}{n}_{j}{\text{σ}}_{mn}^{0}{n}_{n}{\text{Ω}}_{mi}\left(n\right)\stackrel{˜}{c}\left(k\right){\stackrel{˜}{c}}^{*}\left(k\right)\frac{dk}{{\left(2\pi \right)}^{3}}\\ =-\underset{-\infty }{\overset{\infty }{\int }}{n}_{j}{\text{σ}}_{ij}^{0}{\text{Ω}}_{ki}\left(n\right){\text{σ}}_{kl}^{0}{n}_{l}{|\stackrel{˜}{c}\left(k\right)|}^{2}\frac{dk}{{\left(2\pi \right)}^{3}}\end{array}$ (A.18)

$\begin{array}{l}{I}_{3}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\stackrel{˜}{\epsilon }}_{ij}^{0}\left(k\right){\left[{\stackrel{˜}{\epsilon }}_{kl}^{0}\left(k\right)\right]}^{*}\frac{dk}{{\left(2\pi \right)}^{3}}\\ =\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\epsilon }_{ij}^{0}{\epsilon }_{kl}^{0}\stackrel{˜}{c}\left(k\right){\stackrel{˜}{c}}^{*}\left(k\right)\frac{dk}{{\left(2\pi \right)}^{3}}\end{array}$ (A.19)

$=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}{C}_{ijkl}{\epsilon }_{ij}^{0}{\epsilon }_{kl}^{0}{|\stackrel{˜}{c}\left(k\right)|}^{2}\frac{dk}{{\left(2\pi \right)}^{3}}$ (A.20)

Adding I1, I2 and I3, the elastic energy in Eqn (A.15) is simplified as

$\begin{array}{c}{F}_{el}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}\left[{C}_{ijkl}{\epsilon }_{ij}^{0}{\epsilon }_{kl}^{0}-{n}_{j}{\text{σ}}_{ij}^{0}{\text{Ω}}_{ki}\left(n\right){\text{σ}}_{kl}^{0}{n}_{l}\right]{|\stackrel{˜}{c}\left(k\right)|}^{2}\frac{dk}{{\left(2\pi \right)}^{3}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\frac{1}{2}\underset{-\infty }{\overset{\infty }{\int }}B\left(n\right){|\stackrel{˜}{c}\left(k\right)|}^{2}\frac{dk}{{\left(2\pi \right)}^{3}}\end{array}$ (A.21)

where,

$B\left(n\right)={C}_{ijkl}{\epsilon }_{ij}^{0}{\epsilon }_{kl}^{0}-{n}_{j}{\text{σ}}_{ij}^{0}{\text{Ω}}_{ki}\left(n\right){\text{σ}}_{kl}^{0}{n}_{l}$