Numerical analysis of phase decomposition in A-B binary alloys using Cahn-Hilliard equations

Numerical analysis of phase decomposition in A-B binary alloys using Cahn-Hilliard equations

Susana Lezama-Alvarez^I; Erika O. Avila-Davila^II; Victor M. Lopez-Hirata^{I, *}; Jorge L. Gonzalez-Velazquez^I

^IInstituto Politécnico Nacional - IPN, Escuela Superior de Ingeniería Química e Industrias Extractivas - ESIQIE, Apartado Postal 118-556, México, D.F. 07051

^IIInstituto Tecnológico de Pachuca - ITP, Pachuca de Soto Hgo. C.P. 42080, Mexico

ABSTRACT

The analysis of phase decomposition was carried out using the nonlinear and linear Cahn-Hilliard equations in a hypothetical A-B alloy system with a miscibility gap. These equations were solved by the explicit finite difference method assuming a regular solution model. The supersaturated solid solution and decomposed phases were considered to have an fcc structure. Different aging temperatures and thermodynamic interaction parameters Ω_A-B were used to simulate different alloy systems. The numerical simulation results showed that the growth kinetics of phase decomposition in the alloy with 30at.% A was slower than that of 50 at.% A. Additionally, the start time and modulation wavelength of phase decomposition are strongly affected by the thermodynamic interaction parameter Ω_A-B value. The numerical simulation results showed that the growth kinetics of phase decomposition with the linear equation is slower than that with the nonlinear one.

Keywords: A-B binary alloys, phase decomposition, linear and nonlinear Cahn-Hilliard equations, microstructural simulation

1. Introduction

The microstructure evolution during the heat treatments of industrial alloys plays an important role in the mechanical properties of final products. Nowadays, the use of numerical methods has become a good alternative for the modelation and simulation of microstructural evolution during the heating of alloys either as a part of the heat treatment or during the service-operation of industrial components. Since these methods have permitted to simulate different heat treating or operating conditions, which are sometimes difficult to obtain in a practical way for instance, very prolonged times¹. Additionally, the use of numerical methods can be used to understand in more detail the mechanism and growth kinetics of different phase transformations, which occur during the heating of alloys^2,3.

Several numerical methods⁴ have been used to analyze the phase transformations in alloys. One of these methods is the called phase-field method which usually is based on a solution of the nonlinear Cahn-Hilliard equation using mainly thermodynamic and atomic diffusion data⁵. The phase-field method has been used to simulate phase transformations in different alloy systems^2,3,6,7.

One of the typical applications of the phase-field method is the microstructural simulation of the phase decomposition in alloy systems with a miscibility gap⁸. This kind of alloy systems is important for different industrial alloys such as, Cu-Ni, Al-Zn, Fe-Cr, Cu-Ni-Fe, Cu-Ni-Sn, Cu-Ni-Cr, Fe-Cr-Co, etc. Thus, it is necessary to study the effect of the different parameters of Cahn-Hilliard equations on the growth kinetics and the microstructural evolution of this type of spinodally-decomposed alloys.

Additionally, it has been pointed out^9,10 that the linear Cahn-Hilliard equation is suitable to analyze the early stage of phase decomposition during aging of these alloy systems. It is also interesting to analyze the application of the linear equation to the microstructural simulation for this kind of phase transformation and for the comparison to the nonlinear equation. It is important to mention that there is no such comparison work reported in the literature. Furthermore, most of the simulation works^11,12 reported in the literature, using the phase field method based on the nonlinear Cahn-Hilliard equation for hypothetic binary A-B alloys, use mathematical functions to model the spinodal curve in a free energy vs. composition diagram, but these ones do not represent thermodynamic parameters which can be associated with a specific alloy system. Therefore, it is necessary to use free energy functions based on thermodynamic solution models in order to approach the microstructure simulation to that taking place in real alloy systems¹³.

Thus, the goal of this work is to analyze the phase decomposition process in a binary A-B alloy system using both the linear and nonlinear Cahn-Hilliard equations, in order to understand the effect of the equation parameters on the microstructural evolution and growth kinetics for this transformation.

2. Numerical Method

2.1. Cahn-Hilliard equations

The linear Cahn-Hilliard equation ^9,10 is the basis for the theory of the spinodal decomposition in alloys, developed by Cahn and Hilliard. This theory has been used to analyze the phase decomposition in numerous alloys⁸. The most common mathematical expression is the following:

where c is the concentration of either A or B elements as a function of a position vector and time t, M is the atomic mobility, f the local free energy, and K the gradient energy coefficient.

In contrast, the nonlinear Cahn-Hilliard equation has been the base for numerical simulation of different metallurgical phenomena⁵ such as, solidification, recrystallization, phase decomposition, etc. This equation is expressed as follows:

The same type of parameters of the linear equation is also involved in the nonlinear one. The linearization of equation, Equation 1, can be obtained from Equation 2, if c is assumed to be only slightly different from its average value⁹.

Both equations are partial differential equations and therefore, they can be solved using the finite difference method⁴.

2.2. Thermodynamical parameters

As stated above, one of the important parameters of the nonlinear Cahn-Hilliard equation corresponds to the local free energy f which can be defined in a simple way using the strict regular solution model¹³ for a binary alloy system as follows:

where R is the gas constant, T is the absolute temperature. X_A and X_B are the mole fractions of A and B, respectively. f_A and f_B are the molar free energy of pure element A and B, respectively, and Ω_A-B is the interaction parameter¹³ between the A and B atoms. Ω_A-B >0 when A and B atoms are repulsive. On the contrary, Ω_A-B <0 when A and B atoms are attractive. Figure 1 shows the plots of free energy f versus composition c for the A-B binary system at different temperatures considering the value of Ω_A-B as a multiple value of gas constant R, 1000 R, 1500R and 2000R. The values of f_A and f_B were assumed to be equal to R (T-T_A) and R (T-T_B), respectively. T_A and T_B correspond to the melting point of A and B, respectively, shown in Table 1.

Figure 1.
Plot of free energy f, vs. composition c.

Table 1.
Values of parameters for simulation.

The f vs. c curve at high temperatures is convex which indicates that the single phase is stable over all the entire composition range. At the lowest temperature, the f vs. c curve shows two minima and one maximum. These points are defined where df/dc=0 and all the minima, corresponding to each temperature, represent the miscibility gap or binodal curve shown in Figure 2. The inflection points, where the f vs. c curve changes from convex to concave downward, are defined by the second derivative df ²/dc²=0. All the inflection points for each temperature form the spinodal curve, also shown in Figure 2. The spinodal curve is parabolic and it is inscribed in the binodal curve at its vertex (c=0.5). Therefore, the critical temperature¹³ of the two-phase decomposition T_C can be expressed by:

Figure 2.
Equilibrium phase diagram for Ω_A-B= 16628 J/mol.

2.3. Diffusion parameters

The diffusion parameter of both Cahn-Hilliard equations is the atomic mobility M, which is usually defined¹⁰ as a function of the interdiffusion coefficient D by the following expression:

D can be defined⁸ for a binary alloy system as follows:

where D_A and D_B are the diffusion coefficients of A and B, respectively, and they can be expressed by an Arrhenius equation as follows⁸:

where D₀ is the frequency factor and it depends on the crystalline structure and Q the activation energy for atomic diffusion and it is a function of the melting point⁸. These parameters are also shown in Table 1.

2.4. Other parameters

The gradient energy coefficient K can be defined as proposed by Hilliard⁹:

where H^M is the heat of mixing per volume unit and r₀ is the nearest-neighbor distance. The heat of mixing H^M was determined according to the following equation¹⁰:

The nearest-neighbor distance r_o was estimated considering that the solution treated and decomposed phases had an fcc crystalline structure. The lattice parameter a is shown in Table 1.

Additionally, the elastic-strain energy due to the coherency between the decomposed phases was neglected since it was considered a similar lattice parameter of the decomposed phases in order to have a very low lattice misfit, which is related directly to the coherency elastic-stain energy¹⁰.

2.5. Computer program

Both Can-Hilliard equations, Equations 1 and 2, were solved numerically in two dimensions using the explicit finite difference method with 101 × 101 points-square grit with a mesh size of 0.25 nm and a time-step size of 1 s. The computer program was coded in FORTRAN 95. Table 1 shows the numerical values of variables used for this simulation. The selected chemical compositions were 30 and 50 at. % of a hypothetic element A, which correspond to the asymmetrical and symmetrical positions in the miscibility gap of Figure 2. The aging temperatures and times were 400-900 K and 1000-100000000 s, respectively.

3. Results and Discussion

3.1. Concentration profiles

Figures 3 and 4 show the concentration profiles, the plot of the concentration of A element versus distance, for the B-50 and 30at.%A alloys, respectively, aged at 650 K for the case of Ω_A-B= 16628 J/mol (2000 R). The increase in the amplitude of the initial fluctuation with the aging time can be observed in both alloys. This fact confirms the occurrence of the phase by the spinodal decomposition mechanism in the A-B alloy system⁹. The phase separation reaction took place as follows:

Figure 3.
Concentration profile of the A-50at. % B alloy aged at 650 K for different times for Ω_A-B= 16628 J/mol.

Figure 4.
Concentration profile of the A-30at.%B alloy aged at 650 K for different times for Ω_A-B= 16628 J/mol.

The initial fluctuation has to get a minimum wavelength value in order to start the spinodal decomposition in both cases. This is in agreement with the theory of spinodal decomposition which states that there is a minimum wavelength in order to increase the modulation amplitude with time⁹. The evolution of the composition modulations was faster in the case of the symmetric alloy, 50 %at. A, than that of the asymmetric alloy, 30 at.% A. According to the spinodal decomposition theory of Cahn and Hilliard⁹, the amplitude of modulation, c-c_o, is given by

in which c_o is the average composition and A(β,t) is the amplitude of the Fourier component of the wave number β (β=2π/λ) at time t expressed in terms of the initial amplitude at time t=0 as follows:

Where R(β) is the amplification factor given b

Inside the spinodal ∂²f/∂c² <0 and R(β) >0 for all values of β. Thus, any modulation will grow, according to Equation 12. In order to obtain higher and faster fluctuation amplitudes, the value of the driving force for spinodal decomposition, ∂²f/∂c², must be also high, according to Equations 11 and 13. The variation of the second derivative with composition is shown in Figure 5 for the case of Ω_{A - B}= 16628 J/mol at T= 650 K. This figure shows that the symmetric alloy has a higher value of the driving force, ∂²f/∂c², than that of the asymmetric one and thus the kinetic behavior of phase decomposition kinetics is faster in the former alloy.

Figure 5.
Plot of second derivative of free energy f vs. composition c for Ω_A-B= 16628 J/mol at 650 K.

Figures 6 and 7 show the concentration profiles for the A-50at.%B alloy aged at 650 K for the case of Ω_A-B= 6628 J/mol using the nonlinear and linear Cahn-Hilliard equations, respectively. The increase in the amplitude of the initial fluctuation with aging time can be also observed in both cases. This increase confirms the occurrence of the phase separation by the spinodal decomposition mechanism in the A-B alloy system⁹.

Figure 6.
Concentration profile of the A-50at. % B alloy aged at 650 K for different times for Ω_A-B= 16628 J/mol with nonlinear equation.

Figure 7.
Concentration profile of the A-50at.%B alloy aged at 650 K for different times for Ω_A-B= 16628 J/mol with linear equation.

It is interesting to notice that the composition modulation is more sinusoidal and periodic in the case of the linear equation. In contrast, the coarsening process of the composition fluctuations is more evident with the nonlinear equation than in the linear one.

The initial fluctuation has to get a minimum wavelength value in order to start the spinodal decomposition in both cases. The simulation with the nonlinear equation requires a smaller minimum value of wavelength than that with the linear one. Besides, the modulation amplitude reached with the nonlinear equation is higher than that corresponding to the linear one. According to the spinodal decomposition theory of Cahn and Hilliard⁹, the minimum wavelength λ_min for spinodal decomposition is given by:

In order to obtain a lower value of minimum wavelength and higher fluctuation amplitude, the value of driving force for spinodal decomposition,

∂

²

f/

∂c

², must be also high, according to Equation 14. This fact suggests that the nonlinear equation involves a higher driving force than that of the linear equation. That is, the driving force is represented by the first derivative of the free energy with respect to the composition,

∂f/

∂c, in the nonlinear equation, Equation 2, and by the second derivative,

∂

²

f/

∂c

² in the linear equation, Equation 1. The variation of the first derivative and second derivative are shown in

Figures 8 and

5, respectively, for the case of Ω

_A-B= 16628 J/mol at T= 650 K. In general, the change in the first derivative with composition is higher between the inflexion points,

c= 0.2 and 0.8, than that of the second derivative. This fact means that the driving force for the nonlinear equation is higher than that in the linear one.

Figure 8.
Plot of first derivative of free energy f vs. composition c for Ω_A-B= 16628 J/mol.

3.2. Microstructural evolution

Figures 9 and 10 show the microstructural evolution predicted for the B-50 at.% A and 30at.%A alloys for the case of Ω_A-B=16628 J/mol at 650 K. The black and gray zones correspond to the A-rich and B-rich phases, respectively. The microstructure was determined using the two-dimension concentration results. The morphology of decomposed phases is irregular an interconnected for both alloys. This type of morphology is a characteristic at the early stages of the spinodal decomposition⁹. The interconnected microstructure is more clearly observed in the case of the alloy with 50 at.%A than at that corresponding to the other alloy composition. This fact can also be attributed to the higher driving force for the phase decomposition in the former alloy.

Figure 9.
Microstructural evolution for B-50at.%A alloy for Ω_A-B= 16628 J/mol at 650 K.

Figure 10.
Microstructural evolution for B-30at.%A alloy for Ω_A-B= 16628 J/mol at 650 K.

Additionally, the morphology of decomposed phases was irregular an interconnected for both equations in the case of Ω_A-B=16628 J/mol at 650 K, as expected in the early stages of spinodal decomposition⁹. Nevertheless, the size of decomposed phases predicted by the nonlinear equation is much smaller than that of the linear one.

3.3. Decomposition kinetics

To determine the modulation wavelength of the composition fluctuations, the concentration profiles of the present work were analyzed using the autocorrelation analysis¹². Figure 11 shows the growth kinetics of spinodal decomposition, plot of modulation wavelength λ versus aging time t, for the B-30 and 50 at.% A alloys for Ω_A-B=16628 J/mol at 650 K.

Figure 11.
Plots of modulation wavelength λ vs. time t for B-30 and 50at.%A alloys for Ω_A-B= 16628 J/mol at 650 K.

In both cases, the modulation wavelength remains almost constant with time, which is also a characteristic of the spinodal decomposition process⁹. The start of phase decomposition occurred first for the symmetric alloy than for the asymmetric one. This can also be attributed to the higher driving force decomposition of the former alloy.

According to the spinodal decomposition theory of Cahn and Hilliard⁹, the minimum wavelength λ_min for spinodal decomposition is given by Equation 14. It is evident that the modulation wavelength λ_min, is shorter if the alloy composition has either a higher driving force ∂²f/∂c² or a lower gradient energy coefficient K.

The modulation wavelengths λ_min , shown in Figure 11, are very close for the symmetric and asymmetric alloy compositions. This can be associated with very close values of K and ∂²f/∂c² for these compositions.

Figure 12 shows the plot of modulation wavelength λ versus aging time t for B-50at.%A alloy aged at 650 K for Ω_A-B values of 12471 and 16628 J/mol. The modulation wavelength λ_min is shorter as Ω_A-B increases. Besides, the start time is shorter for the phase decomposition with the highest value of Ω_A-B parameter, which can be attributed to a higher driving force. It is important to mention that the kinetic behavior of phase decomposition increases as the aging time increases in Figure 12. This fact can be related to the coarsening stage of the decomposed phases⁶.

Figure 12.
Plots of modulation wavelength λ vs. time t for B-50at.%A alloy at 650 K for Ω_A-B values of 12471 and 16628 J/mol.

The effect of the aging temperature on the modulation growth kinetics is shown in Figure 13 for B-50at.%A alloy with a Ω_A-B value of 16628 J/mol. The increase in the aging temperature caused a lower driving force, see Figure 1; nevertheless, the increase in the atomic mobility, due to the higher atomic diffusion, promoted that the phase decomposition started first at a higher aging temperature. The minimum modulation wavelength is very close for the two aging temperatures since they have the same interaction parameter.

Figure 13.
Plots of modulation wavelength λ vs. time t for B-50at.%A alloy at 650 K and 700 K for Ω_A-B= 16628 J/mol.

Figure 14 shows the growth kinetics of spinodal decomposition for the A-50 at.% B alloy with Ω_A-B values of 12471 J/mol (1500 R) and 16628 J/mol (2000 R) using both equations at 650 K. In both cases, the modulation wavelength remains almost constant with time, which is also a characteristic of the spinodal decomposition process⁹. It is evident that the modulation wavelength λ, is shorter for the nonlinear equation than for the linear one. In the case of the nonlinear equation, a higher value of Ω_A-B caused that the start of the phase decomposition growth kinetics was faster. This is attributed to the higher values of the first derivative as Ω_A-B increases. A higher first derivative value, or driving force, can be adopted as the reason for the shorter minimum wavelength λ_min for spinodal decomposition using the nonlinear equation. The higher values for minimum wavelength λ_min in the linear equation seem to be associated with the lower driving force, lower values of second derivative, as explained in section 3.1. Therefore, a longer aging time is needed to reach the longer wavelength λ_min for the start of the spinodal decomposition.

Figure 14.
Plots of modulation wavelength λ vs. time t for linear and nonlinear equations for Ω_A-B values of 12471 and 16628 J/mol.

3.4. Comparison with actual alloy systems

The use of low interaction Ω_A-B parameters causes that the T_c temperature of the miscibility gap to be located at low temperatures, see Equation 4. This promotes a low atomic diffusion which makes the phase decomposition process be very slow. For instance, Cu-Ni alloys¹² have a miscibility gap at low temperatures, lower than 573 K and its kinetic behavior of phase decomposition is very slow due to its low atomic diffusion at these temperatures. Furthermore, it has been reported that the phase decomposition takes longer aging times than several years. In contrast, the use of high Ω_A-B values makes the T_c temperature be higher which promotes the atomic diffusion and thus the appearance of the phase decomposition. For instance, the phase decomposition in Fe-Cr alloys occurs at higher temperatures⁶ and the phase decomposition process is much more rapid than that observed in Cu-Ni alloys. Table 2 shows a comparison of simulated results, the modulation wavelength and morphology of decomposed phases, for the present work using nonlinear equation with other simulated and experimental results for Cu-Ni and Fe-Cr alloys^14-17. The simulation results corresponding to Ω_A-B values of 8314 and 12471 J/mol have the T_c temperature of the miscibility gap located at temperatures between 500 and 750 K. This temperature range is close to that of the miscibility gap in Cu-Ni alloys^14-15. Furthermore, the minimum value of the modulation wavelength for the present work is in good agreement with the experimental and calculated values reported in the literature for Cu-Ni alloys. In contrast, the Ω_A-B and T_c values of 16628 J/mol are close to those in Fe-Cr alloys. The minimum value of modulation wavelength for the A-B alloy shows a good agreement with those reported in the literature for Fe-Cr alloys^16-17. Finally, the morphology prediction of the decomposed phases corresponds to that reported in the literature.

Table 2.
Comparison of present work simulated results with other works.

4. Conclusions

Both linear and nonlinear Cahn-Hilliard equations reproduced the basic characteristics of the morphology and growth kinetics of the spinodal decomposition in the early stages of aging for hypothetical A-B alloys. The increase in the thermodynamic interaction parameter Ω_A-B caused a shorter start time for the phase decomposition. The higher aging temperature promoted the faster phase decomposition due to the higher atomic diffusion. The kinetic behavior for the phase decomposition calculated by the linear equation is slower than that of the nonlinear one which can be attributed to the difference in the driving force parameter of each equation. The microstructure evolution and kinetics of phase decomposition of Cu-Ni and Fe-Cr alloys showed a good agreement with the present work results.

Acknowledgements

The authors wish to acknowledge the financial support from SIP-IPN-CONACYT 100584.

Received: August 15, 2012

Revised: January 31, 2013

• 1. Janssens KGF, Raabe D, Kozeschnik E, Miodownik MA and Nestler B. Computational Materials Engineering. 1st ed. Elsevier Academic Press; 2007.
• 2. Avila-Davila EO, Melo-Maximo DV, Lopez-Hirata VM, Soriano-Vargas O, Saucedo-Muñoz ML and Gonzalez-Velazquez JL. Microstructural simulation in spinodally-decomposed Cu-70at.%Ni and Cu-46at.%Ni-4at.%Fe alloys. Materials Characterization 2009; 60:560-567. http://dx.doi.org/10.1016/j.matchar.2009.01.003
• 3. Montaño-Zuñiga IM, Sepulveda-Cervantes G, Lopez-Hirata VM, Rivas-Lopez DI and Gonzalez-Velazquez JL. Numerical simulation of recrystallization In bcc metals. Computational Materials Science 2010; 49:512-517. http://dx.doi.org/10.1016/j.commatsci.2010.05.042
• 4. Rapaz M, Bellete M and Deville M. Materials Modelling in Materials Science and Engineering 1st ed. Berlin: Springer-Verlag; 2010.
• 5. Chen LQ. Phase Field Modeling of Materials Microstructure. In: Xiao Guo Z, editor. Multiscale materials modeling. CRC Press LLC; 2007. p. 62-83. http://dx.doi.org/10.1533/9781845693374.62
• 6. Honjo M and Saito Y. Numerical simulation of phase separation in Fe-Cr binary and Fe-Cr-Mo ternary alloys with use of the Cahn-Hilliard equation. ISIJ International 2000; 40:914-919. http://dx.doi.org/10.2355/isijinternational.40.914
• 7. Kim SG and Kim WT. Phase field modeling of solidification. In: Yip S, editor. Handbook of Materials Modeling. Netherlands: Springer; 2005. p. 2106-2116.
• 8. Kostorz G. Phasetransformations in Materials 2nd ed. Wiley-VCH; 2001. PMid:11505305. http://dx.doi.org/10.1002/352760264X
• 9. Hilliard JE. Spinodal Decomposition. In: Aaronson HI, editor. Phase transformations. Metals Park Ohio ASM; 1970. p. 497-560.
• 10. Cahn JW. Spinodal decomposition. Transactions of the Metallurgical Society of AIME 1970; 242:89-103.
• 11. Avila-Davila EO, Lezama-Alvarez S, Saucedo-Muñoz ML, Lopez Hirata VM and Gonzalez-Velazquez JL. Numerical simulation of the spinodal decomposition in hypothetical A-B and A-B-C alloy systems. Revista de Metalurgia, 2012; 48:223-236.
• 12. Chen LQ. Computer Simulation of Spinodal Decomposition in Ternary Systems. Acta Metal Mater 1994; 42:3503-3513. http://dx.doi.org/10.1016/0956-7151(94)90482-0
• 13. Nishizawa T. Thermodynamics of Microstructure 1st ed. ASM International; 2008.
• 14. Lopez-Hirata VM, Sakurai T and Hirano K. A study of phase separation in Cu-Ni alloys by AP-FIM. Scripta Metallurgica et Materialia. 1992; 26:99-103. http://dx.doi.org/10.1016/0956-716X(92)90377-Q
• 15. Avila-Davila EO, Lopez Hirata VM, Saucedo-Muñoz ML and Gonzalez-Velazquez JL. Microstructural simulation of phase decomposition in Cu-Ni alloys. Journal of Alloys and Compounds 2008: 460:206-212. http://dx.doi.org/10.1016/j.jallcom.2007.05.070
• 16. Soriano-Vargas O Avila-Davila EO, Lopez Hirata VM, Saucedo-Muñoz ML, Cayetano-Castro N and Gonzalez-Velazquez JL. Effect of spinodal decomposition on The mechanical behavior of Fe-Cr alloys. Materials Science and Engineering: A. 2010; 527:2910-2914. http://dx.doi.org/10.1016/j.msea.2010.01.020
• 17. Lopez Hirata VM, Soriano-Vargas O, Dorantes-Rosales HJ and Saucedo-Muñoz ML. Phase decomposition in an Fe-40at.%Cr alloy after isothermal aging and its effect on hardening. Materials Characterization 2011; 62:789-792. http://dx.doi.org/10.1016/j.matchar.2011.05.012

Publication Dates

• Publication in this collection
  28 May 2013
• Date of issue
  Oct 2013

History

• Received
  15 Aug 2012
• Accepted
  31 Jan 2013