Phase-field modeling of bimodal particle size ...
Acta Materialia 51 (2003) 1123–1132 www.actamat-journals.com
Phase-field modeling of bimodal particle size distributions during continuous cooling Y.H. Wen a,∗, J.P. Simmons b, C. Shen c, C. Woodward a, Y. Wang c b
a UES, Inc., Dayton, OH 45432, USA Air Force Research Laboratory, AFRL/MLLM, Wright-Patterson AFB, OH 45433, USA c The Ohio State University, 2041 College Road, Columbus, OH 43210, USA
Received 11 June 2002; accepted 21 October 2002
Abstract Microstructures in Nickel-base alloys typically contain a two-phase mixture of γ/γ⬘. The microstructure having a bimodal size distribution of γ⬘ is of particular interest because it has important property consequences [1]. In this paper, the phase-field method with an explicit nucleation algorithm is employed to investigate the microstructural development during a continuous cooling with various cooling rates. It is demonstrated that bimodal particle size distributions can be achieved at an intermediate cooling rate due to a coupling between diffusion and undercooling, in which the system experiences two peaks of well-isolated nucleation events. It is suggested that this is caused by soft impingement, followed by a renewal of driving force for nucleation, followed by a subsequent soft impingement. Under very high cooling rates, the microstructure becomes unimodal, because undercooling always outruns diffusion and the microstructure never reaches soft impingement. 2002 Acta Materialia Inc. Published by Elsevier Science Ltd. All rights reserved. Keywords: Phase-field approach; Bimodal particle size distributions; Continuous cooling; Ni-base alloy
1. Introduction Alloy properties can be tailored by modifying processing conditions that affect microstructure formations. Alloy development cycles typically require between 10 and 20 yr for completion, while engine design cycles are typically 3–5 yr and are under pressure to be shortened. Unless the alloy development cycle can be shortened, engine devel-
Corresponding author. Tel.: +1-37-255-3514; fax: +1-937255-0445. E-mail address: [email protected] (Y.H. Wen). ∗
opment will use off-the-shelf alloys and alloy development will no longer play a direct role in engine design. One obvious choice for shortening this cycle time is to employ simulation methods towards modeling alloy behavior. This work explores some of the details of the phase-field model [2,3,4] when it is adapted to simulate microstructural evolution under more realistic heat treating processes. In particular, we investigate the formation of bimodal distributions that have been reported to form in superalloy systems under continuous cooling conditions [1,5]. Microstructural development under non-isothermal conditions involves nucleation, growth,
1359-6454/03/$30.00 2002 Acta Materialia Inc. Published by Elsevier Science Ltd. All rights reserved. doi:10.1016/S1359-6454(02)00516-5
1124
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
and coarsening (for a general review, see [6]), in which the three mechanisms are often overlapping and even competing [7]. Concurrent nucleation and growth was first treated theoretically for near-critical liquids by Langer and Schwartz [8] based on a mean field model. This treatment was extended to handle the precipitation in alloys with rather high supersaturation by Wendt and Haasen [9] and later by Kampmann and Wagner to account for non-linear contributions to capillarity and to eliminate assumptions as to the shape of the particle size distribution [7,10]. A similar numerical model was formulated to predict the formation of γ⬘ precipitation in Ni-base superalloys during the cooling process [1]. While the Langer–Schwartz approach can successfully predict particle size distributions during a concomitant nucleation, growth, and coarsening process, it can not provide additional information such as interparticle spacings. The phase-field method has received some attention in the literature recently because of its ability to simulate realistic microstructures as well as to incorporate the effects of elastic interactions, overlapping diffusion fields, interparticle correlation effects, and other fine scale properties. Computer simulations based on phase-field method has been proved to be very useful in studying the formation and dynamic evolution of complex microstructures. For example, some of the general features of the morphologic patterns formed during a hexagonal-to-orthorhombic transformation have been successfully predicted [11,12]. It was shown that the long-range elastic interactions arising from the lattice accommodation among different orientation domains of the orthorhombic phase dominate the domain morphologies and the kinetics of domain coarsening [13]. The simulation of the microstructural development during precipitation of a coherent orthorhombic phase from an a2 matrix in a Ti– Al–Nb system demonstrated how volume fractions of the precipitate could dramatically influence their morphology and mutual arrangement. The complex lamellar structure in γTiAl system was also simulated using a 3D phase-field model without any prior assumption as how the microstructure may develop [14]. This work uses the explicit nucleation phasefield method, as developed by Simmons, Shen, and
Wang [15] to explore the formation of bimodal distributions under continuous cooling conditions. The model is described in Section 2. Simulation results are presented in Section 3, demonstrating that bimodal particle size distribution could be achieved at an intermediate cooling rate. In Section 4, we compare our simulation results with experimental ones and discuss the formation of bimodal particle size distribution in terms of the competition of nucleation and growth during a continuous cooling process. Brief concluding remarks are presented in Section 5. 2. Computer model 2.1. The phase-field method An arbitrary multi-phase and/or multi-domain microstructure is described by a few mesoscopic field parameters. To describe the L12 ordered structure of γ⬘ phase, we need to employ three long range order (lro) parameters, h1(r,t), h2(r,t), and h3(r,t). To describe the concentration profile, we need another field parameter, c(r,t), where r represents spatial coordinate and t is time [16]. The spatio-temporal evolution of these four phase-field parameters describes the microstructural evolution of the two-phase mixture. The spatio-temporal evolution of the lro parameters can be obtained by solving the time-dependent Ginzburgh–Landau equation while the temporal evolution of the concentration field can be described by the non-linear Cahn–Hilliard diffusion equation, i.e. dF ∂hp ⫽ ⫺L ; p ⫽ 1,2,3. ∂t dhp ∂c dF ⫽ Mⵜ2 , ∂t dc
(1)
where L and M are respectively structural relaxation and diffusion mobilities, F is the total free energy of the system. In the mean field approximation, M is the atomic mobility and can be related to diffusion coefficient D through M ⫽ c0(1⫺ c0)D / RT in a dilute solution, where c0 represents the mean composition of the alloy, R is the gas constant and T is absolute temperature.
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
To numerically solve Eq. (1), one needs to express the total free energy F as a function of the aforementioned four phase-field parameters. For the coherent precipitation under consideration, the total free energy should consist of two terms: the elastic strain energy (Fel) and the chemical free energy (Fch), namely F ⫽ F el ⫹ F ch. In the following, we describe the formulation of Fel and Fch in terms of the four chosen phase-field parameters. Elastic strain energy arises due to the misfit between γ and γ⬘ phases. The misfit is considered to be dependent solely on the composition. As a result, the elastic strain energy can be formulated in terms of the concentration profile c(r, t). The total elastic strain energy can be divided into two terms [17], i.e., Fel ⫽ E0 ⫹ Erelax with E0 ⫽ c11 ⫹ c12 V ⬍ ⌽(n) ⬎ 苸20c¯ (1⫺c¯ ) (c ⫹ c12) 3⫺ 2 11 c11 1 ’ dk and Erelax ⫽ ⫺ B(n)|c˜ (k)|2 [16]. In these 2 (2p)3 expressions, V is the volume of the whole system, cij are elastic constants of the system, n ⫽ k / k is a unit vector in reciprocal space and ni is its ith component, 苸0 ⫽ da(c) / (ao dc) is the concentration coefficient of crystal lattice parameter, a(c) is the lattice parameter of a solid solution with concentration c and a0 is the lattice parameter of pure solvent. c¯ is the total atomic fraction of the solute
冉
冕
冊
atoms, c˜ (k) is the Fourier transform of c(r),
冕
’
indi-
cates that the point n ⫽ 0 is excluded in the integration, B(n) ⫽ [(c11 ⫹ 2c12)2苸20 / c11]b(n) with b(n) ⫽ ⌽(n)⫺ ⬍ ⌽(n) ⬎ , where ⬍⬎ denotes an average over all directions of the unit vector n, and for the 2-D case ⌽(n) ⫽ c11(1 ⫹ 2zn21n22) / [c11 ⫹ z(c11 ⫹ c12)n21n22], where z ⫽ (c11⫺c12⫺2c44) / c44 characterizes the elastic anisotropy of the system. Assuming that the lro parameters and the concentration profile change appreciably on a length scale which is much larger than the interatomic distance, the non-equilibrium chemical free energy as a functional of the field parameters can be approximated by the Ginzburg–Landau coarse-grained free energy functional [18], which contains a local specific free energy, f(c,h1,h2,h3), and gradient energy terms, e.g.
Fch ⫽
1125
冕冋
1 1 c(ⵜc)2 ⫹ h 2 2
V
冘
(ⵜhp)2
(2)
p
册
⫹ f(c,h1,h2,h3)dV
where c and h are gradient energy coefficients and, in principle, can be fitted to the data from experiment or atomistic calculations. The gradient terms provide an energy penalty to non-homogeneities in composition and lro parameters which take place mainly at interfaces. The local specific free energy f(c,h1,h2,h3) defines the basic thermodynamic properties of the system. It can be approximated by a Landau-type polynomial expansion. In principle, the polynomial may include any terms allowed by the symmetry operations with respect to the parent phase [19]. In this work, we use the following formulation: f(c,h1,h2,h3) ⫽ ⫹
a1(T) [c⫺c1(T)]2 2
a2(T) [c2(T)⫺c](h21 ⫹ h22 ⫹ h23) 2
(3)
a3(T) a4(T) 4 ⫺ h1h2h3 ⫹ (h1 ⫹ h42 ⫹ h43) 3 4 which is the same as used by Wang et al.[16]. In their work, all phenomenological constants, i.e. (c1,c2,a1,a2,a3,a4), were fixed in value for isothermal condition. In this work, however, their values are updated according to the temperature history, which will be further described in Section 2.4. 2.2. Dimensionless phase-field kinetic equations For the convenience of numerical solution, Eq. (1) was transformed into a dimensionless form through the introduction of a reduced time, defined as t ⫽ LtQ ch, and reduced spatial coordinates r∗ with r∗i ⫽ xi / l, where Θch is a positive constant introduced as a characteristic local chemical free energy and l is the length unit of the computational grid size [16]. Substituting the energy formulation described in the above section into Eq. (1) and using these reduced parameters, the dimensionless form of Eq. (1) reads,
1126
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
冉
冊
∂f∗ ∂hp ; p ⫽ 1,2,3. ⫽ ⫺ ⫺j1ⵜ∗2hp ⫹ ∂t ∂hp
冉
(4)
冊
∂f∗ ∂c ⫽ cⵜ∗2 ⫺j2ⵜ∗2c ⫹ ⫺m{b(n)c˜ (k∗)}r∗ , ∂t ∂c where {b(n)c˜ (k∗)}r∗ is the inverse Fourier transform of b(n)c˜ (k∗). c ⫽ M / Ll2 characterizes relative mobility for diffusion and interface migration, j1 ⫽ h / ⌰chl2 and j2 ⫽ c / ⌰chl2 characterize the relative strength of gradient energies with respect to the chemical driving force, f∗ ⫽ f / ⌰ch is a rescaled local chemical free energy, and m ⫽ (c11 ⫹ 2c12)2苸20 / (c11⌰ch) characterizes the strength of elastic strain energy relative to the chemical driving force in governing the microstructural evolution. The microstructural evolution is governed by a set of dimensionless parameters, namely c, ϕ1, ϕ2, m, and f∗. In Section 2.4, we will describe how their values are selected. 2.3. Nucleation model By solving Eq. (4), the growth and coarsening of precipitates are fully determined. A missing component of the model is the description of nucleation of γ⬘ from γ phase. Conventionally, the Langevin noise term [2] is used to ‘punch out’ nuclei at the beginning of the simulation and that it is switched off at some point and the transformation is allowed to proceed [16]. To simulate nucleation with a noise term will require that we set an extremely small value on the time stepsize in order to catch the (rare event) nucleation of a particle. With the stepsize that small, it is not possible to simulate growth due to computational limitations. Simmons et al. [15] developed a hybrid model recently, in which stochastic nucleation events during isothermal processing are explicitly incorporated to a phase-field model. Nuclei were introduced into individual cells randomly, but with the probability related to the corresponding nucleation rate at those cells. With this treatment of nucleation process, the random noise terms are no longer needed in the phase-field kinetic equations. The nucleation rate (J∗) was related to the supersaturation (⌬c) through J∗=1e-2/⌬c for a 2D formulation, where 1 and 2 were assumed to be constants for isothermal application [15].
To account for the effect of undercooling during a cooling process, we replaced ⌬c by (⌬c ⫹ r⌬T), where r is a phenomenological constant characterizing relative contribution of undercooling with respect to supersaturation for determining the nucleation rate. Furthermore, the prefactor 1 in classic theory is considered to be proportional to the rate of atomic attachment to the critical nuclei, which is temperature sensitive. To take this into account, 1 is replaced by 1χ, and the nucleation rate is evaluated using the expression J∗ ⫽ 1ce-[2/(⌬c+r⌬T)],
(5)
where 1, 2 and ρ are assumed to be constants and they take the values of 1020 cell⫺1, 3.0, and 10⫺4 K⫺1, respectively, for all simulations presented in this work. 2.4. Values for the dimensionless parameters In order to attempt to be close to realistic parameters, the Ni–Al binary alloy was used as a guide in choosing values for some phenomenological constants, but it should be mentioned that these results will not compare quantitatively with Ni–Al. Work is underway in our lab. calibrating the phasefield model to commercial superalloys. Similar to the variation of diffusion coefficient D with temperature, we assume that c decreases with temperature following an exponential law, i.e. c(T) ⫽ c 0 exp( ⫺ Q / RT), where c0 and Q are two constants to be determined. c was set to 0.064 for isothermal modeling at 1273 K to get reasonable microstructure [15]. We retained that value for this work. Ardell measured the diffusion coefficients at 625 and 715 °C for Ni–Al binary alloy and they are, respectively, 9.20 × 10⫺16 and 2.79 × 10⫺14 cm2 s⫺1 [20]. We require that the ratio of c at the two temperatures equals the ratio of their corresponding diffusion coefficients. As a result, c0 and Q are determined to be 1.91 × 1010 and 2.80 × 105 J mol⫺1, respectively. The explicit form of c as a function of temperature reads c(T) ⫽ 1.91 × 1010exp(⫺2.80 × 105 / RT).
(6)
The interfacial energy coefficients h and c can generally be related to the interatomic interaction energies [16,18]. An alternative way to determine
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
1127
them is by fitting the specific APB and interfacial energy calculated from the computer simulation to the corresponding experimental values [21,11]. These two coefficients are generally temperature dependent. As a result, ϕ1 and ϕ2 should vary with temperature. Based on Ardell’s experimental measurement on Ni–Al binary alloy, however, the temperature dependence of the specific interfacial energy is quite weak [20]. For example, its values are 14.4 and 14.2 erg cm⫺2 at 625 and 715 °C, respectively. So in this work, we ignored the temperature dependency of ϕ1 and ϕ2, they are assumed to be 2.5 and ⫺3, respectively. Wang et al. discussed why f2 should take a negative value [16]. m is assumed to be 2100 to fit microstructures. The rescaled local chemical free energy f∗ can be expressed as f∗ ⫽
a1∗(T) a∗2 (T) [c⫺c1(T)]2 ⫹ [c2(T)⫺c] 2 2
a3∗(T) a∗4 (T) (h21 ⫹ h22 ⫹ h23)⫺ h1h2h3 ⫹ 3 4
Fig. 1.
Schematic illustration of the heat treatment scheme.
(7)
(h41 ⫹ h42 ⫹ h43). The local chemical free energy at equilibrium corresponds to the order parameter that minimizes the free energy expression. If the phase is disordered, h 1 ⫽ h 2 ⫽ h 3 ⫽ 0. If the phase is ordered, the h assumed a non-zero value at equlibrium. The free energy expression for the disordered phase g ) can be found by substituting 0 in for h in Eq. (fch (7). That for the ordered phase (fg⬘ ch) can be found by solving for the h that minimizes the total free energy and substituting this value back into Eq. (7). The common tangent of these two free energy curves determines the equilibrium compositions of γ and γ⬘. It is required that these two values equal to the equilibrium composition of γ and γ⬘ (c1 and c2), respectively. All the phenomenological constants, i.e., (a∗1 ,a∗2 ,a∗3 ,a∗4 ,c1,c2), are adjusted accordingly to fit the phase boundary of Ni–Al binary alloy [22]. 2.5. Simulations A relevant portion of the Ni–Al phase diagram is illustrated in Fig. 1 for describing our heat-treatment scheme. All simulations started with homo-
geneous γ phase which corresponds to experimentally holding at point A for a long time. The sample was then quenched to point B (1000 °C), which is 200 °C below the equilibrium temperature for the γPγ⬘ phase transformation for the chosen mean composition, and followed by a continuous cooling at a fixed cooling rate to 550 °C. All simulations presented in this work were performed in 2-D on a 512 × 512 square mesh. Periodic boundary conditions were applied.
3. Results In the figures to show simulated results, the shades of grey represent the value of composition c(r∗, t), e.g. the higher the value, the brighter the shade. Therefore, the white and black areas represent γ⬘ and γ phase, respectively. We begin with a simulation with a cooling rate of 1 K t⫺1. The microstructural evolution was accompanied with nucleations up to t ⫽ 25, which will be referred to as the first peak of nucleation. This was followed by a long period of time where nucleation was completely shut down.
1128
Fig. 2.
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
Microstructural evolution with the cooling rate at 1 K t-1 (a) t ⫽ 20; (b) t ⫽ 40; (c) t ⫽ 150; (d) t ⫽ 450.
A new peak of nucleation started at t ⫽ 250 near the end of the simulation. This indicates that the system experienced two peaks of well-separated nucleations. In Fig. 2, we show the microstructures at three representative stages, i.e., near the end of the first peak of nucleation [Fig. 2(a)], during the shut-off of nucleation [Fig. 2(b,c)], and at the end of the simulation [Fig. 2(d)]. While most of those particles formed during the first peak of nucleation grow with time there are some exceptions. In Fig. 2(b), a few such particles are highlighted. They initially grow [compare their sizes in Fig. 2(a,b)] and dissolve later on [compare their sizes in Fig. 2(b,c)]. This competitive precipitate dissolution during precipitation was observed experimentally in the Ni–Ni3Al system (see [6] and references therein). During the second peak of nucleation, some smaller particles are introduced [Fig. 2(d)]. The bimodal dual size of γ⬘ particles is a result of the two peaks of well-separated nucleation with those bigger particles being formed during the first peak of nucleation. The mean nucleation rate initially increased and soon it decreased sharply for more than ten orders of magnitude by 900 °C which corresponds to t ⫽ 100, indicating an actual shut-off of the nucleation. This sharply decrease was followed by a increase and a second peak of nucleation occurred at a later stage. The influence of fast cooling on microstructures was investigated and the final microstructures are shown in Fig. 4 for three representative cases with cooling rates of 5, 25, and 100 K t⫺1 respectively. For the case under cooling rate of 5 K t⫺1 [Fig. 4(a)], the nucleation was shut down only briefly between t ⫽ 20 and t ⫽ 30. As a result, the
Fig. 3. Plot of the evolution of mean nucleation rate ⬍NR⬎ per phase-field cell in logarithm with temperature in °C. The corresponding cooling rate is 1 K t⫺1.
size difference between the larger particles and the smaller ones are much less significant as compared to that shown in Fig. 2(d). When the cooling rate was further increased to 25 and 100 K t⫺1, both cases had but a single peak of nucleation. The nucleation was active from the beginning to a complete shut-off when the depletion zone overlaps completely. Accordingly, the sizes of the particles
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
Fig. 4.
1129
Simulated final microstructures under higher cooling rates: (a) 5 K t⫺1; (b) 25 K t⫺1; (c) 100 K t⫺1.
are quite uniform. The final microstructures are featured with fine particles of limited growth and coarsening due to the successive nucleation and the largely handicapped kinetics associated with the fast cooling. Comparing the microstructures shown in Fig. 2(d) and those in Fig. 4, there is a clear tendency for the particle density to increase and the overall particle size to decrease with increasing cooling rate. A plot of the the mean nucleation rate against temperature for the case at cooling rate of 100 K t⫺1 (Fig. 5) demonstrates that the mean nucleation rate keeps increasing from the very beginning until very later stage, which is very different from what was shown in Fig. 3 for a case with an intermediate cooling rate.
4. Discussion Gabb et al. formulated a kinetic model for describing γ⬘ formation in a Ni-base disk superalloy [1]. The model gave reasonable predictions of primary cooling γ⬘ size over a very broad range of cooling conditions. They found that the volume fraction of primary cooling γ⬘ decreased with higher cooling rate in the web and rim of the disk, allowing a larger volume fraction of tertiary γ⬘ precipitation which lead to the bimodal size distributions. Babu et al. also characterized the microstructures under various cooling rates including 0.17, 1, 10, 75 K s⫺1, and water-quenching (refer to Fig. 3 in their paper [5]). They found that: (1) particle density increases and the particle size decreases
Fig. 5. Plot of the evolution of mean nucleation rate ⬍NR⬎ per phase-field cell in logarithm with temperature in °C. The corresponding cooling rate is 100 K t⫺1.
with increasing cooling rate; (2) morphological change of precipitates from cuboidal to spherical with increasing cooling rate; (3) bimodal microstructure can be found at cooling rate of 1 K s⫺1, an intermediate cooling rate in their experiments. These prior works indicate that the cooling rate is an important factor in the formation of bimodal
1130
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
particle size distribution. Our simulation results are in line with these observations. Moreover, our simulation results demonstrate that bimodal particle size distribution is due to two peaks of nucleations. The first peak of nucleation is shut down by soft impingement. With further cooling, the second peak of nucleation shows up, bringing in those smaller γ⬘ particles. The extended period of shut-off of nucleation in-between the two peaks of nucleation results in the large size difference between the γ⬘ particles nucleated in the two peaks of nucleation. Nucleation shuts off with soft impingement. This will occur when the interparticle spacing becomes smaller than the diffusion distance, √Dt [6]. We check here if our simulation results are consistent with the √Dt criterion. A simple isothermal case at 1000 °C is used for this purpose, which was realized by setting the cooling rate at zero. The simulated microstructural evolution is shown in Fig. 6. The simulation revealed that nucleations were active only at the early stage of the evolution. The last nucleation event happened at t ⫽ 35 ± 10, where the standard deviation was obtained from multiple runs, and the microstructural evolution thereafter is controlled by growth and coarsening of those previously formed nuclei. The nucleation rate [Fig. 6(c)] started to decrease from the very beginning and by t ⫽ 100 the mean nucleation rate decreased by about 15 orders of magnitude.
To assess √Dt at a given reduced time t, we need to convert the dimensionless t to a real world time t and assess the diffusion coefficient D. In our formulation, t was defined as: t ⫽ Lt⌰ch.
(8)
L can be estimated through the expression c ⫽ M / (Ll2). c was set to 0.064 at 1000 °C. The mobility M at 1000 °C is estimated to be 1.45 × 10⫺19 m2 mol J⫺1 s⫺1 using the data provided in [23,24]. Using these two numbers, L is estimated to be 2.26 × 10⫺18 / l 2 mol J⫺1 s⫺1, where l is in the unit of meter. ⌰ch can be estimated through the expression m ⫽ (c11 ⫹ 2c 12)2苸20 / (c11⌰ch), where m ⫽ 2100 is an input data in our simulation, cij are elastic constants with (c 11,c 12,c 44) ⫽ (2.31, 1.49, 1.17) × 1012 erg cm⫺3 [25], 苸0 ⫽ (aγ⬘⫺aγ) / aγ(cγ⬘⫺cγ) ⫽ 0.049 assuming Vegard’s law where (aγ⬘⫺aγ) / aγ⯝0.0056 is the γ/γ⬘ lattice misfit [26] and cγ⬘ ⫽ 0.24 and cγ ⫽ 0.125 are equilibrium composition of γ⬘ and γ, respectively. Using these numbers, ⌰ch is estimated to be 1.39 × 106 J m⫺3. Substituting values of L and ⌰ch into Eq. (8) we get t ⫽ 2.58 × 1016 l 2 t s. The diffusion coefficient D was related to the ∂2f mobility M through D ⫽ 2M ⫽ a1(T)M[27] ∂c where f is the dimensional free energy as defined in Eq. (3). Using the input value of a1(1000 °C) in our model and the value of M at 1000 °C one has D ⫽ 3.63 × 10⫺15 m2 s⫺1 at 1000 °C. Using this
Fig. 6. Isothermal treatment at 1000 °C. (a) Microstructures shortly after the shup-off of nucleation event (t ⫽ 50). The line on the upper-left corner shows the length scale of evaluated √Dt, see text for more detail. (b) Final microstructures at t ⫽ 450. (c) Mean nucleation rate ⬍NR⬎ per phase-field cell in logarithm with reduced time.
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
assessed D and the converted t, the estimated √Dt at t ⫽ 35 is 57l. This length scale is shown on the upper-left corner of Fig. 6(a). It is about half of the average interparticle distances. Therefore, our simulation result is consistent with the √Dt rule. The √Dt rule indicates that there is a typical length-scale for a given temperature. Nucleation will shut down by soft impingement when interparticle distances are below this length-scale. Due to the fact that D changes exponentially with temperature, this length-scale √Dt changes with temperature accordingly. Under lower temperature, the interparticle distances are smaller because of a smaller D to prevent further nucleation and vice versa. To understand microstructures under continuous cooling using the √Dt rule, we need to understand the role of the cooling rate first. If the cooling rate is extremely slow, the microstructure should be similar to the isothermal case due to the fact that D is changing very slowly. If the cooling rate is very high, the length-scale of shut-off of nucleation is constantly decreasing—so fast that nucleation will never stop in the temperature range of our simulation. For the case with an intermediate cooling rate, the situation is more complicated. The nucleation events can be totally shut down after the first peak of nucleation because of the soft impingement at higher temperature (with a higher interparticle distances). With further cooling, our simulations suggest that superstaturation increases faster than diffusion depletion at lower temperature, which makes the second peak of nucleation possible as we simulated under the cooling rate of 1 K t⫺1.
5. Concluding remarks The phase-field method is used to investigate the formation of bimodal particle size distributions under continuous cooling conditions. It is demonstrated that bimodal particle size distributions could be achieved at an intermediate cooling rate. At an intermediate cooling rate, the system develops a bimodal microstructure through two peaks of well-isolated nucleation events. The first
1131
peak of nucleation events is shut down by soft impingement. With further cooling superstaturation increases faster than diffusion depletion leading to the second peak of nucleation events.
Acknowledgments YHW acknowledges support from the Air Force Research Laboratory through contract #F3361501-5214 with UES Inc. CS and YZW acknowledge the support from NSF under grants DMR-0080766 and DMR9905725.
References [1] Gabb TP et al. In: Superalloys 2000. Warrendale, PA: TMS; 2000. p. 405. [2] Gunton JD, Miguel MS, Sahni PS. In: Phase transitions and critical phenomena. New York: Academic Press; 1983. p. 267. [3] Elder K. Comput. Phys. 1993;7:27. [4] Chen LQ, Wang YZ. JOM 1996;48:13. [5] Babu SS, Miller MK, Vitek JM, David SA. Acta Mater. 2001;49:4149. [6] Doherty RD. In: Physical metallurgy. North-Holland, Amsterdam: Elsevier; 1996. p. 1364. [7] Wagner R, Kampmann R. In: Materials science and technology: a comprehensive treatment. New York: VCH; 1991. p. 213. [8] Langer JS, Schwartz AJ. Phys. Rev. 1980;21:948. [9] Wendt H, Haasen P. Acta Metall. 1983;31:1649. [10] Kampmann R, Wagner R. In: Decomposition of alloys: the early stage. New York: Pergamon Press; 1983. p. 91. [11] Wen YH, Wang Y, Chen LQ. Phil. Mag. A 2000;80:1967. [12] Wen YH, Wang Y, Chen LQ. Acta Metall. 1999;47:4375. [13] Wen YH, Wang Y, Chen LQ. Acta Metall. 2002;50:13. [14] Wen YH, Chen LQ, Hazzledine PM, Wang Y. Acta Metall. 2001;49:2341. [15] Simmons JP, Shen C, Wang Y. Scripta Mater. 2000;43:935. [16] Wang Y, Banerjee D, Su CC, Khachaturyan AG. Acta Mater. 1998;46:2983. [17] Khachaturyan AG. Theory of structural transformations in solids. New York: John Wiley & Sons, 1983. [18] Cahn JW, Hilliard JE. J. Chem. Phys. 1958;28:258. [19] Izyumov YuA, Syromyatnikov VN. Phase transitions and crystal symmetry, fundamental theories of physics. New York: Kluwer, 1990. [20] Ardell AJ. Acta Metall. 1968;16:511. [21] Wang Y, Khachaturyan A. Phil. Mag. A 1995;72:1161. [22] Okamoto H. J. Phase Equil. 1993;14:257.
1132
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
[23] Engstom A, Argen J. Z. Metallkd 1996;87:92. [24] Campbell CE, Boettinger WJ, Kattner UR. Acta Mater. 2002;50:775. [25] Pottenbohm H, Neitze G, Nembach E. Mater. Sci. Eng. 1983;60:189.
[26] Miyazaki T, Imamura M, Kozaki T. Mater Sci Eng 1982;54:9. [27] Hilliard JE. In: Phase transformation. Materials Park, OH: ASM International; 1970. p. 497.
Phase-field modeling of bimodal particle size distributions during continuous cooling Y.H. Wen a,∗, J.P. Simmons b, C. Shen c, C. Woodward a, Y. Wang c b
a UES, Inc., Dayton, OH 45432, USA Air Force Research Laboratory, AFRL/MLLM, Wright-Patterson AFB, OH 45433, USA c The Ohio State University, 2041 College Road, Columbus, OH 43210, USA
Received 11 June 2002; accepted 21 October 2002
Abstract Microstructures in Nickel-base alloys typically contain a two-phase mixture of γ/γ⬘. The microstructure having a bimodal size distribution of γ⬘ is of particular interest because it has important property consequences [1]. In this paper, the phase-field method with an explicit nucleation algorithm is employed to investigate the microstructural development during a continuous cooling with various cooling rates. It is demonstrated that bimodal particle size distributions can be achieved at an intermediate cooling rate due to a coupling between diffusion and undercooling, in which the system experiences two peaks of well-isolated nucleation events. It is suggested that this is caused by soft impingement, followed by a renewal of driving force for nucleation, followed by a subsequent soft impingement. Under very high cooling rates, the microstructure becomes unimodal, because undercooling always outruns diffusion and the microstructure never reaches soft impingement. 2002 Acta Materialia Inc. Published by Elsevier Science Ltd. All rights reserved. Keywords: Phase-field approach; Bimodal particle size distributions; Continuous cooling; Ni-base alloy
1. Introduction Alloy properties can be tailored by modifying processing conditions that affect microstructure formations. Alloy development cycles typically require between 10 and 20 yr for completion, while engine design cycles are typically 3–5 yr and are under pressure to be shortened. Unless the alloy development cycle can be shortened, engine devel-
Corresponding author. Tel.: +1-37-255-3514; fax: +1-937255-0445. E-mail address: [email protected] (Y.H. Wen). ∗
opment will use off-the-shelf alloys and alloy development will no longer play a direct role in engine design. One obvious choice for shortening this cycle time is to employ simulation methods towards modeling alloy behavior. This work explores some of the details of the phase-field model [2,3,4] when it is adapted to simulate microstructural evolution under more realistic heat treating processes. In particular, we investigate the formation of bimodal distributions that have been reported to form in superalloy systems under continuous cooling conditions [1,5]. Microstructural development under non-isothermal conditions involves nucleation, growth,
1359-6454/03/$30.00 2002 Acta Materialia Inc. Published by Elsevier Science Ltd. All rights reserved. doi:10.1016/S1359-6454(02)00516-5
1124
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
and coarsening (for a general review, see [6]), in which the three mechanisms are often overlapping and even competing [7]. Concurrent nucleation and growth was first treated theoretically for near-critical liquids by Langer and Schwartz [8] based on a mean field model. This treatment was extended to handle the precipitation in alloys with rather high supersaturation by Wendt and Haasen [9] and later by Kampmann and Wagner to account for non-linear contributions to capillarity and to eliminate assumptions as to the shape of the particle size distribution [7,10]. A similar numerical model was formulated to predict the formation of γ⬘ precipitation in Ni-base superalloys during the cooling process [1]. While the Langer–Schwartz approach can successfully predict particle size distributions during a concomitant nucleation, growth, and coarsening process, it can not provide additional information such as interparticle spacings. The phase-field method has received some attention in the literature recently because of its ability to simulate realistic microstructures as well as to incorporate the effects of elastic interactions, overlapping diffusion fields, interparticle correlation effects, and other fine scale properties. Computer simulations based on phase-field method has been proved to be very useful in studying the formation and dynamic evolution of complex microstructures. For example, some of the general features of the morphologic patterns formed during a hexagonal-to-orthorhombic transformation have been successfully predicted [11,12]. It was shown that the long-range elastic interactions arising from the lattice accommodation among different orientation domains of the orthorhombic phase dominate the domain morphologies and the kinetics of domain coarsening [13]. The simulation of the microstructural development during precipitation of a coherent orthorhombic phase from an a2 matrix in a Ti– Al–Nb system demonstrated how volume fractions of the precipitate could dramatically influence their morphology and mutual arrangement. The complex lamellar structure in γTiAl system was also simulated using a 3D phase-field model without any prior assumption as how the microstructure may develop [14]. This work uses the explicit nucleation phasefield method, as developed by Simmons, Shen, and
Wang [15] to explore the formation of bimodal distributions under continuous cooling conditions. The model is described in Section 2. Simulation results are presented in Section 3, demonstrating that bimodal particle size distribution could be achieved at an intermediate cooling rate. In Section 4, we compare our simulation results with experimental ones and discuss the formation of bimodal particle size distribution in terms of the competition of nucleation and growth during a continuous cooling process. Brief concluding remarks are presented in Section 5. 2. Computer model 2.1. The phase-field method An arbitrary multi-phase and/or multi-domain microstructure is described by a few mesoscopic field parameters. To describe the L12 ordered structure of γ⬘ phase, we need to employ three long range order (lro) parameters, h1(r,t), h2(r,t), and h3(r,t). To describe the concentration profile, we need another field parameter, c(r,t), where r represents spatial coordinate and t is time [16]. The spatio-temporal evolution of these four phase-field parameters describes the microstructural evolution of the two-phase mixture. The spatio-temporal evolution of the lro parameters can be obtained by solving the time-dependent Ginzburgh–Landau equation while the temporal evolution of the concentration field can be described by the non-linear Cahn–Hilliard diffusion equation, i.e. dF ∂hp ⫽ ⫺L ; p ⫽ 1,2,3. ∂t dhp ∂c dF ⫽ Mⵜ2 , ∂t dc
(1)
where L and M are respectively structural relaxation and diffusion mobilities, F is the total free energy of the system. In the mean field approximation, M is the atomic mobility and can be related to diffusion coefficient D through M ⫽ c0(1⫺ c0)D / RT in a dilute solution, where c0 represents the mean composition of the alloy, R is the gas constant and T is absolute temperature.
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
To numerically solve Eq. (1), one needs to express the total free energy F as a function of the aforementioned four phase-field parameters. For the coherent precipitation under consideration, the total free energy should consist of two terms: the elastic strain energy (Fel) and the chemical free energy (Fch), namely F ⫽ F el ⫹ F ch. In the following, we describe the formulation of Fel and Fch in terms of the four chosen phase-field parameters. Elastic strain energy arises due to the misfit between γ and γ⬘ phases. The misfit is considered to be dependent solely on the composition. As a result, the elastic strain energy can be formulated in terms of the concentration profile c(r, t). The total elastic strain energy can be divided into two terms [17], i.e., Fel ⫽ E0 ⫹ Erelax with E0 ⫽ c11 ⫹ c12 V ⬍ ⌽(n) ⬎ 苸20c¯ (1⫺c¯ ) (c ⫹ c12) 3⫺ 2 11 c11 1 ’ dk and Erelax ⫽ ⫺ B(n)|c˜ (k)|2 [16]. In these 2 (2p)3 expressions, V is the volume of the whole system, cij are elastic constants of the system, n ⫽ k / k is a unit vector in reciprocal space and ni is its ith component, 苸0 ⫽ da(c) / (ao dc) is the concentration coefficient of crystal lattice parameter, a(c) is the lattice parameter of a solid solution with concentration c and a0 is the lattice parameter of pure solvent. c¯ is the total atomic fraction of the solute
冉
冕
冊
atoms, c˜ (k) is the Fourier transform of c(r),
冕
’
indi-
cates that the point n ⫽ 0 is excluded in the integration, B(n) ⫽ [(c11 ⫹ 2c12)2苸20 / c11]b(n) with b(n) ⫽ ⌽(n)⫺ ⬍ ⌽(n) ⬎ , where ⬍⬎ denotes an average over all directions of the unit vector n, and for the 2-D case ⌽(n) ⫽ c11(1 ⫹ 2zn21n22) / [c11 ⫹ z(c11 ⫹ c12)n21n22], where z ⫽ (c11⫺c12⫺2c44) / c44 characterizes the elastic anisotropy of the system. Assuming that the lro parameters and the concentration profile change appreciably on a length scale which is much larger than the interatomic distance, the non-equilibrium chemical free energy as a functional of the field parameters can be approximated by the Ginzburg–Landau coarse-grained free energy functional [18], which contains a local specific free energy, f(c,h1,h2,h3), and gradient energy terms, e.g.
Fch ⫽
1125
冕冋
1 1 c(ⵜc)2 ⫹ h 2 2
V
冘
(ⵜhp)2
(2)
p
册
⫹ f(c,h1,h2,h3)dV
where c and h are gradient energy coefficients and, in principle, can be fitted to the data from experiment or atomistic calculations. The gradient terms provide an energy penalty to non-homogeneities in composition and lro parameters which take place mainly at interfaces. The local specific free energy f(c,h1,h2,h3) defines the basic thermodynamic properties of the system. It can be approximated by a Landau-type polynomial expansion. In principle, the polynomial may include any terms allowed by the symmetry operations with respect to the parent phase [19]. In this work, we use the following formulation: f(c,h1,h2,h3) ⫽ ⫹
a1(T) [c⫺c1(T)]2 2
a2(T) [c2(T)⫺c](h21 ⫹ h22 ⫹ h23) 2
(3)
a3(T) a4(T) 4 ⫺ h1h2h3 ⫹ (h1 ⫹ h42 ⫹ h43) 3 4 which is the same as used by Wang et al.[16]. In their work, all phenomenological constants, i.e. (c1,c2,a1,a2,a3,a4), were fixed in value for isothermal condition. In this work, however, their values are updated according to the temperature history, which will be further described in Section 2.4. 2.2. Dimensionless phase-field kinetic equations For the convenience of numerical solution, Eq. (1) was transformed into a dimensionless form through the introduction of a reduced time, defined as t ⫽ LtQ ch, and reduced spatial coordinates r∗ with r∗i ⫽ xi / l, where Θch is a positive constant introduced as a characteristic local chemical free energy and l is the length unit of the computational grid size [16]. Substituting the energy formulation described in the above section into Eq. (1) and using these reduced parameters, the dimensionless form of Eq. (1) reads,
1126
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
冉
冊
∂f∗ ∂hp ; p ⫽ 1,2,3. ⫽ ⫺ ⫺j1ⵜ∗2hp ⫹ ∂t ∂hp
冉
(4)
冊
∂f∗ ∂c ⫽ cⵜ∗2 ⫺j2ⵜ∗2c ⫹ ⫺m{b(n)c˜ (k∗)}r∗ , ∂t ∂c where {b(n)c˜ (k∗)}r∗ is the inverse Fourier transform of b(n)c˜ (k∗). c ⫽ M / Ll2 characterizes relative mobility for diffusion and interface migration, j1 ⫽ h / ⌰chl2 and j2 ⫽ c / ⌰chl2 characterize the relative strength of gradient energies with respect to the chemical driving force, f∗ ⫽ f / ⌰ch is a rescaled local chemical free energy, and m ⫽ (c11 ⫹ 2c12)2苸20 / (c11⌰ch) characterizes the strength of elastic strain energy relative to the chemical driving force in governing the microstructural evolution. The microstructural evolution is governed by a set of dimensionless parameters, namely c, ϕ1, ϕ2, m, and f∗. In Section 2.4, we will describe how their values are selected. 2.3. Nucleation model By solving Eq. (4), the growth and coarsening of precipitates are fully determined. A missing component of the model is the description of nucleation of γ⬘ from γ phase. Conventionally, the Langevin noise term [2] is used to ‘punch out’ nuclei at the beginning of the simulation and that it is switched off at some point and the transformation is allowed to proceed [16]. To simulate nucleation with a noise term will require that we set an extremely small value on the time stepsize in order to catch the (rare event) nucleation of a particle. With the stepsize that small, it is not possible to simulate growth due to computational limitations. Simmons et al. [15] developed a hybrid model recently, in which stochastic nucleation events during isothermal processing are explicitly incorporated to a phase-field model. Nuclei were introduced into individual cells randomly, but with the probability related to the corresponding nucleation rate at those cells. With this treatment of nucleation process, the random noise terms are no longer needed in the phase-field kinetic equations. The nucleation rate (J∗) was related to the supersaturation (⌬c) through J∗=1e-2/⌬c for a 2D formulation, where 1 and 2 were assumed to be constants for isothermal application [15].
To account for the effect of undercooling during a cooling process, we replaced ⌬c by (⌬c ⫹ r⌬T), where r is a phenomenological constant characterizing relative contribution of undercooling with respect to supersaturation for determining the nucleation rate. Furthermore, the prefactor 1 in classic theory is considered to be proportional to the rate of atomic attachment to the critical nuclei, which is temperature sensitive. To take this into account, 1 is replaced by 1χ, and the nucleation rate is evaluated using the expression J∗ ⫽ 1ce-[2/(⌬c+r⌬T)],
(5)
where 1, 2 and ρ are assumed to be constants and they take the values of 1020 cell⫺1, 3.0, and 10⫺4 K⫺1, respectively, for all simulations presented in this work. 2.4. Values for the dimensionless parameters In order to attempt to be close to realistic parameters, the Ni–Al binary alloy was used as a guide in choosing values for some phenomenological constants, but it should be mentioned that these results will not compare quantitatively with Ni–Al. Work is underway in our lab. calibrating the phasefield model to commercial superalloys. Similar to the variation of diffusion coefficient D with temperature, we assume that c decreases with temperature following an exponential law, i.e. c(T) ⫽ c 0 exp( ⫺ Q / RT), where c0 and Q are two constants to be determined. c was set to 0.064 for isothermal modeling at 1273 K to get reasonable microstructure [15]. We retained that value for this work. Ardell measured the diffusion coefficients at 625 and 715 °C for Ni–Al binary alloy and they are, respectively, 9.20 × 10⫺16 and 2.79 × 10⫺14 cm2 s⫺1 [20]. We require that the ratio of c at the two temperatures equals the ratio of their corresponding diffusion coefficients. As a result, c0 and Q are determined to be 1.91 × 1010 and 2.80 × 105 J mol⫺1, respectively. The explicit form of c as a function of temperature reads c(T) ⫽ 1.91 × 1010exp(⫺2.80 × 105 / RT).
(6)
The interfacial energy coefficients h and c can generally be related to the interatomic interaction energies [16,18]. An alternative way to determine
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
1127
them is by fitting the specific APB and interfacial energy calculated from the computer simulation to the corresponding experimental values [21,11]. These two coefficients are generally temperature dependent. As a result, ϕ1 and ϕ2 should vary with temperature. Based on Ardell’s experimental measurement on Ni–Al binary alloy, however, the temperature dependence of the specific interfacial energy is quite weak [20]. For example, its values are 14.4 and 14.2 erg cm⫺2 at 625 and 715 °C, respectively. So in this work, we ignored the temperature dependency of ϕ1 and ϕ2, they are assumed to be 2.5 and ⫺3, respectively. Wang et al. discussed why f2 should take a negative value [16]. m is assumed to be 2100 to fit microstructures. The rescaled local chemical free energy f∗ can be expressed as f∗ ⫽
a1∗(T) a∗2 (T) [c⫺c1(T)]2 ⫹ [c2(T)⫺c] 2 2
a3∗(T) a∗4 (T) (h21 ⫹ h22 ⫹ h23)⫺ h1h2h3 ⫹ 3 4
Fig. 1.
Schematic illustration of the heat treatment scheme.
(7)
(h41 ⫹ h42 ⫹ h43). The local chemical free energy at equilibrium corresponds to the order parameter that minimizes the free energy expression. If the phase is disordered, h 1 ⫽ h 2 ⫽ h 3 ⫽ 0. If the phase is ordered, the h assumed a non-zero value at equlibrium. The free energy expression for the disordered phase g ) can be found by substituting 0 in for h in Eq. (fch (7). That for the ordered phase (fg⬘ ch) can be found by solving for the h that minimizes the total free energy and substituting this value back into Eq. (7). The common tangent of these two free energy curves determines the equilibrium compositions of γ and γ⬘. It is required that these two values equal to the equilibrium composition of γ and γ⬘ (c1 and c2), respectively. All the phenomenological constants, i.e., (a∗1 ,a∗2 ,a∗3 ,a∗4 ,c1,c2), are adjusted accordingly to fit the phase boundary of Ni–Al binary alloy [22]. 2.5. Simulations A relevant portion of the Ni–Al phase diagram is illustrated in Fig. 1 for describing our heat-treatment scheme. All simulations started with homo-
geneous γ phase which corresponds to experimentally holding at point A for a long time. The sample was then quenched to point B (1000 °C), which is 200 °C below the equilibrium temperature for the γPγ⬘ phase transformation for the chosen mean composition, and followed by a continuous cooling at a fixed cooling rate to 550 °C. All simulations presented in this work were performed in 2-D on a 512 × 512 square mesh. Periodic boundary conditions were applied.
3. Results In the figures to show simulated results, the shades of grey represent the value of composition c(r∗, t), e.g. the higher the value, the brighter the shade. Therefore, the white and black areas represent γ⬘ and γ phase, respectively. We begin with a simulation with a cooling rate of 1 K t⫺1. The microstructural evolution was accompanied with nucleations up to t ⫽ 25, which will be referred to as the first peak of nucleation. This was followed by a long period of time where nucleation was completely shut down.
1128
Fig. 2.
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
Microstructural evolution with the cooling rate at 1 K t-1 (a) t ⫽ 20; (b) t ⫽ 40; (c) t ⫽ 150; (d) t ⫽ 450.
A new peak of nucleation started at t ⫽ 250 near the end of the simulation. This indicates that the system experienced two peaks of well-separated nucleations. In Fig. 2, we show the microstructures at three representative stages, i.e., near the end of the first peak of nucleation [Fig. 2(a)], during the shut-off of nucleation [Fig. 2(b,c)], and at the end of the simulation [Fig. 2(d)]. While most of those particles formed during the first peak of nucleation grow with time there are some exceptions. In Fig. 2(b), a few such particles are highlighted. They initially grow [compare their sizes in Fig. 2(a,b)] and dissolve later on [compare their sizes in Fig. 2(b,c)]. This competitive precipitate dissolution during precipitation was observed experimentally in the Ni–Ni3Al system (see [6] and references therein). During the second peak of nucleation, some smaller particles are introduced [Fig. 2(d)]. The bimodal dual size of γ⬘ particles is a result of the two peaks of well-separated nucleation with those bigger particles being formed during the first peak of nucleation. The mean nucleation rate initially increased and soon it decreased sharply for more than ten orders of magnitude by 900 °C which corresponds to t ⫽ 100, indicating an actual shut-off of the nucleation. This sharply decrease was followed by a increase and a second peak of nucleation occurred at a later stage. The influence of fast cooling on microstructures was investigated and the final microstructures are shown in Fig. 4 for three representative cases with cooling rates of 5, 25, and 100 K t⫺1 respectively. For the case under cooling rate of 5 K t⫺1 [Fig. 4(a)], the nucleation was shut down only briefly between t ⫽ 20 and t ⫽ 30. As a result, the
Fig. 3. Plot of the evolution of mean nucleation rate ⬍NR⬎ per phase-field cell in logarithm with temperature in °C. The corresponding cooling rate is 1 K t⫺1.
size difference between the larger particles and the smaller ones are much less significant as compared to that shown in Fig. 2(d). When the cooling rate was further increased to 25 and 100 K t⫺1, both cases had but a single peak of nucleation. The nucleation was active from the beginning to a complete shut-off when the depletion zone overlaps completely. Accordingly, the sizes of the particles
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
Fig. 4.
1129
Simulated final microstructures under higher cooling rates: (a) 5 K t⫺1; (b) 25 K t⫺1; (c) 100 K t⫺1.
are quite uniform. The final microstructures are featured with fine particles of limited growth and coarsening due to the successive nucleation and the largely handicapped kinetics associated with the fast cooling. Comparing the microstructures shown in Fig. 2(d) and those in Fig. 4, there is a clear tendency for the particle density to increase and the overall particle size to decrease with increasing cooling rate. A plot of the the mean nucleation rate against temperature for the case at cooling rate of 100 K t⫺1 (Fig. 5) demonstrates that the mean nucleation rate keeps increasing from the very beginning until very later stage, which is very different from what was shown in Fig. 3 for a case with an intermediate cooling rate.
4. Discussion Gabb et al. formulated a kinetic model for describing γ⬘ formation in a Ni-base disk superalloy [1]. The model gave reasonable predictions of primary cooling γ⬘ size over a very broad range of cooling conditions. They found that the volume fraction of primary cooling γ⬘ decreased with higher cooling rate in the web and rim of the disk, allowing a larger volume fraction of tertiary γ⬘ precipitation which lead to the bimodal size distributions. Babu et al. also characterized the microstructures under various cooling rates including 0.17, 1, 10, 75 K s⫺1, and water-quenching (refer to Fig. 3 in their paper [5]). They found that: (1) particle density increases and the particle size decreases
Fig. 5. Plot of the evolution of mean nucleation rate ⬍NR⬎ per phase-field cell in logarithm with temperature in °C. The corresponding cooling rate is 100 K t⫺1.
with increasing cooling rate; (2) morphological change of precipitates from cuboidal to spherical with increasing cooling rate; (3) bimodal microstructure can be found at cooling rate of 1 K s⫺1, an intermediate cooling rate in their experiments. These prior works indicate that the cooling rate is an important factor in the formation of bimodal
1130
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
particle size distribution. Our simulation results are in line with these observations. Moreover, our simulation results demonstrate that bimodal particle size distribution is due to two peaks of nucleations. The first peak of nucleation is shut down by soft impingement. With further cooling, the second peak of nucleation shows up, bringing in those smaller γ⬘ particles. The extended period of shut-off of nucleation in-between the two peaks of nucleation results in the large size difference between the γ⬘ particles nucleated in the two peaks of nucleation. Nucleation shuts off with soft impingement. This will occur when the interparticle spacing becomes smaller than the diffusion distance, √Dt [6]. We check here if our simulation results are consistent with the √Dt criterion. A simple isothermal case at 1000 °C is used for this purpose, which was realized by setting the cooling rate at zero. The simulated microstructural evolution is shown in Fig. 6. The simulation revealed that nucleations were active only at the early stage of the evolution. The last nucleation event happened at t ⫽ 35 ± 10, where the standard deviation was obtained from multiple runs, and the microstructural evolution thereafter is controlled by growth and coarsening of those previously formed nuclei. The nucleation rate [Fig. 6(c)] started to decrease from the very beginning and by t ⫽ 100 the mean nucleation rate decreased by about 15 orders of magnitude.
To assess √Dt at a given reduced time t, we need to convert the dimensionless t to a real world time t and assess the diffusion coefficient D. In our formulation, t was defined as: t ⫽ Lt⌰ch.
(8)
L can be estimated through the expression c ⫽ M / (Ll2). c was set to 0.064 at 1000 °C. The mobility M at 1000 °C is estimated to be 1.45 × 10⫺19 m2 mol J⫺1 s⫺1 using the data provided in [23,24]. Using these two numbers, L is estimated to be 2.26 × 10⫺18 / l 2 mol J⫺1 s⫺1, where l is in the unit of meter. ⌰ch can be estimated through the expression m ⫽ (c11 ⫹ 2c 12)2苸20 / (c11⌰ch), where m ⫽ 2100 is an input data in our simulation, cij are elastic constants with (c 11,c 12,c 44) ⫽ (2.31, 1.49, 1.17) × 1012 erg cm⫺3 [25], 苸0 ⫽ (aγ⬘⫺aγ) / aγ(cγ⬘⫺cγ) ⫽ 0.049 assuming Vegard’s law where (aγ⬘⫺aγ) / aγ⯝0.0056 is the γ/γ⬘ lattice misfit [26] and cγ⬘ ⫽ 0.24 and cγ ⫽ 0.125 are equilibrium composition of γ⬘ and γ, respectively. Using these numbers, ⌰ch is estimated to be 1.39 × 106 J m⫺3. Substituting values of L and ⌰ch into Eq. (8) we get t ⫽ 2.58 × 1016 l 2 t s. The diffusion coefficient D was related to the ∂2f mobility M through D ⫽ 2M ⫽ a1(T)M[27] ∂c where f is the dimensional free energy as defined in Eq. (3). Using the input value of a1(1000 °C) in our model and the value of M at 1000 °C one has D ⫽ 3.63 × 10⫺15 m2 s⫺1 at 1000 °C. Using this
Fig. 6. Isothermal treatment at 1000 °C. (a) Microstructures shortly after the shup-off of nucleation event (t ⫽ 50). The line on the upper-left corner shows the length scale of evaluated √Dt, see text for more detail. (b) Final microstructures at t ⫽ 450. (c) Mean nucleation rate ⬍NR⬎ per phase-field cell in logarithm with reduced time.
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
assessed D and the converted t, the estimated √Dt at t ⫽ 35 is 57l. This length scale is shown on the upper-left corner of Fig. 6(a). It is about half of the average interparticle distances. Therefore, our simulation result is consistent with the √Dt rule. The √Dt rule indicates that there is a typical length-scale for a given temperature. Nucleation will shut down by soft impingement when interparticle distances are below this length-scale. Due to the fact that D changes exponentially with temperature, this length-scale √Dt changes with temperature accordingly. Under lower temperature, the interparticle distances are smaller because of a smaller D to prevent further nucleation and vice versa. To understand microstructures under continuous cooling using the √Dt rule, we need to understand the role of the cooling rate first. If the cooling rate is extremely slow, the microstructure should be similar to the isothermal case due to the fact that D is changing very slowly. If the cooling rate is very high, the length-scale of shut-off of nucleation is constantly decreasing—so fast that nucleation will never stop in the temperature range of our simulation. For the case with an intermediate cooling rate, the situation is more complicated. The nucleation events can be totally shut down after the first peak of nucleation because of the soft impingement at higher temperature (with a higher interparticle distances). With further cooling, our simulations suggest that superstaturation increases faster than diffusion depletion at lower temperature, which makes the second peak of nucleation possible as we simulated under the cooling rate of 1 K t⫺1.
5. Concluding remarks The phase-field method is used to investigate the formation of bimodal particle size distributions under continuous cooling conditions. It is demonstrated that bimodal particle size distributions could be achieved at an intermediate cooling rate. At an intermediate cooling rate, the system develops a bimodal microstructure through two peaks of well-isolated nucleation events. The first
1131
peak of nucleation events is shut down by soft impingement. With further cooling superstaturation increases faster than diffusion depletion leading to the second peak of nucleation events.
Acknowledgments YHW acknowledges support from the Air Force Research Laboratory through contract #F3361501-5214 with UES Inc. CS and YZW acknowledge the support from NSF under grants DMR-0080766 and DMR9905725.
References [1] Gabb TP et al. In: Superalloys 2000. Warrendale, PA: TMS; 2000. p. 405. [2] Gunton JD, Miguel MS, Sahni PS. In: Phase transitions and critical phenomena. New York: Academic Press; 1983. p. 267. [3] Elder K. Comput. Phys. 1993;7:27. [4] Chen LQ, Wang YZ. JOM 1996;48:13. [5] Babu SS, Miller MK, Vitek JM, David SA. Acta Mater. 2001;49:4149. [6] Doherty RD. In: Physical metallurgy. North-Holland, Amsterdam: Elsevier; 1996. p. 1364. [7] Wagner R, Kampmann R. In: Materials science and technology: a comprehensive treatment. New York: VCH; 1991. p. 213. [8] Langer JS, Schwartz AJ. Phys. Rev. 1980;21:948. [9] Wendt H, Haasen P. Acta Metall. 1983;31:1649. [10] Kampmann R, Wagner R. In: Decomposition of alloys: the early stage. New York: Pergamon Press; 1983. p. 91. [11] Wen YH, Wang Y, Chen LQ. Phil. Mag. A 2000;80:1967. [12] Wen YH, Wang Y, Chen LQ. Acta Metall. 1999;47:4375. [13] Wen YH, Wang Y, Chen LQ. Acta Metall. 2002;50:13. [14] Wen YH, Chen LQ, Hazzledine PM, Wang Y. Acta Metall. 2001;49:2341. [15] Simmons JP, Shen C, Wang Y. Scripta Mater. 2000;43:935. [16] Wang Y, Banerjee D, Su CC, Khachaturyan AG. Acta Mater. 1998;46:2983. [17] Khachaturyan AG. Theory of structural transformations in solids. New York: John Wiley & Sons, 1983. [18] Cahn JW, Hilliard JE. J. Chem. Phys. 1958;28:258. [19] Izyumov YuA, Syromyatnikov VN. Phase transitions and crystal symmetry, fundamental theories of physics. New York: Kluwer, 1990. [20] Ardell AJ. Acta Metall. 1968;16:511. [21] Wang Y, Khachaturyan A. Phil. Mag. A 1995;72:1161. [22] Okamoto H. J. Phase Equil. 1993;14:257.
1132
Y.H. Wen et al. / Acta Materialia 51 (2003) 1123–1132
[23] Engstom A, Argen J. Z. Metallkd 1996;87:92. [24] Campbell CE, Boettinger WJ, Kattner UR. Acta Mater. 2002;50:775. [25] Pottenbohm H, Neitze G, Nembach E. Mater. Sci. Eng. 1983;60:189.
[26] Miyazaki T, Imamura M, Kozaki T. Mater Sci Eng 1982;54:9. [27] Hilliard JE. In: Phase transformation. Materials Park, OH: ASM International; 1970. p. 497.
