Phase-field simulation of abnormal grain growth due to inverse pinning

Available online at www.sciencedirect.com

Acta Materialia 55 (2007) 6881–6894 www.elsevier.com/locate/actamat

Phase-field simulation of abnormal grain growth due to inverse pinning Yoshihiro Suwa a

a,*

, Yoshiyuki Saito b, Hidehiro Onodera

c

Computational Materials Science Center, National Institute for Materials Science, Tsukuba 305-0047, Japan b Department of Materials Science and Engineering, Waseda University, Tokyo 169-8555, Japan c Materials Engineering Laboratory, National Institute for Materials Science, Tsukuba 305-0047, Japan Received 13 April 2007; received in revised form 30 August 2007; accepted 30 August 2007 Available online 18 October 2007

Abstract The possibility of abnormal grain growth due to inverse pinning was verified using phase-field simulations. In bicrystalline systems with circular precipitates, the perfect wetting condition is required for the long-distance migration of the interface between the matrix grains. If the distance between precipitates that are perpendicular to the interface exceeds a critical value, the migration is not observed irrespective of the wetting condition. In polycrystalline systems, abnormal grain growth occurs with the aid of the driving force for grain growth even though llim exceeded the critical value, where llim is the minimum distance between precipitates. Furthermore, the perfect wetting condition is not required for the abnormal grain growth in the polycrystalline systems. These facts enlarge the possibility of inverse pinning in real alloy systems.  2007 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Computer simulation; Abnormal grain growth; Pinning; Phase-field modeling

1. Introduction The grain growth processes that occur in polycrystalline materials during annealing are usually classified into two types. The first type is a self-similar coarsening process, which is called normal grain growth. In normal grain growth, an invariant distribution of scaled grain sizes develops, and the grain growth obeys a power-law kinetics with a characteristic exponent of 1/2 [1–20]. The second type is called abnormal grain growth and is characterized by the coarsening of a few grains at the expense of the surrounding matrix. In such conditions, the self-similar aspect of the coarsening process is lost. Different mechanisms have been proposed for abnormal grain growth, although the actual physical mechanism responsible for this phenomenon remains largely unknown. Defect-induced strains can induce isolated grain growth [21], as well as the same capillary forces responsible for coarsening when anisotropy of grain boundary energies *

Corresponding author. E-mail address: [email protected] (Y. Suwa).

or mobilities exists [22,23]. The conditions for abnormal growth due to variable interfacial energies or mobilities were recently examined [24,25] using a mean-field treatment of the matrix grains. For the case of a single grain with boundary properties that differ from those of the surrounding matrix, it was found that higher boundary mobility generally promotes abnormal growth whereas higher boundary energy constrains it. The detailed behavior can be quite complex, depending on the ranges of the model parameters selected. This includes abnormal growth only up to a limiting grain size, or lower bounds in the initial relative size of the grain for abnormal growth to occur. Abnormal grain growth has also been shown to occur when grain boundaries pin due to, for example, existing precipitate phases or other defects. Simplified models have been proposed that introduce grain boundary drag forces that lead to ultimate pinning (Zener pinning) [26–28], while the role of thermal fluctuations in overcoming pinning has been analyzed by Monte Carlo simulation [29]. Dispersed precipitates have been generally considered to be the inhibitors for boundary migration. However, it is assumed that, when anisotropic precipitates are present in

1359-6454/$30.00  2007 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2007.08.045

6882

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

the matrix and the boundary energy of these precipitates in a grain of type I is smaller than that in the adjacent grain of type II (for simplicity, we consider bicrystalline sample), the grain boundary between two grains is induced to migrate by these precipitates. Nishizawa called this phenomenon ‘‘inverse pinning’’ [30]. Nishizawa suggested that a typical example of inverse pinning was observed in steels when AlN precipitates were present. The AlN is an anisotropic ionic compound with a hexagonal B4 structure that precipitates in steel in the form of platelets or fine needles along specific planes of the parent metallic crystals. Consequently, the selective growth of a favorably oriented grain of type I occurs due to inverse pinning if it has the same orientation as the parent crystal of the group of AlN. In addition, Miodownik [31] suggests that it is very much an open question whether pinning force anisotropy (i.e. the boundary energy is a strongly anisotropic function of boundary plane and/or misorientation) is the cause of abnormal growth or preferred orientations in some systems, leading to strong textures and the attendant misorientation distribution function of nanocrystalline nickel alloys and Fe–Si systems [32]. The main purpose of this study is to verify the possibility of abnormal grain growth due to inverse pinning by using phase-field simulations. All calculations will be performed on a two-dimensional (2D) lattice with a periodic boundary condition. For simplicity, it is assumed that there exist only two groups of orientations (A) and (B) for matrix grains and one for precipitates (P). The system with a single texture component (A) in a matrix of randomly orientated grains (B) is considered. Further, we assume that the precipitates (P) have a special orientation relationship with A, where the difference in the boundary properties of A/ A, A/B, and B/B boundaries is not introduced to emphasize the effects of inverse pinning. We evaluate whether a preferred growth of A grains occurs or not, provided that only interfaces BP have a higher boundary energy and the other interfaces have equivalent and smaller boundary energies. We will perform numerical simulation for bicrystalline systems in Sections 3.1–3.3 and for polycrystalline systems in Sections 3.4 and 3.5. A numbers of simulations with finely dispersed precipitates have been performed using Monte Carlo Potts models [29,33–41], the vertex model [42,43] and the phase-field model [41,44–46]. However, for both experiments and numerical simulations, there have been no systematic investigations of anisotropic pinning. In experiments, it is difficult to change the volume fraction of precipitates and their distribution without changing other variables, such as composition. In numerical simulation, in order to treat the anisotropy in precipitate–matrix interfaces precisely, a higher spatial resolution is required that requires huge computational costs. Systematic investigations have been feasible, provided the volume fraction of precipitates is relatively high [47,48]. However, the computational costs for obtaining statistically meaningful results become higher with decreas-

ing volume fraction. In this study, for realizing the simulation with the anisotropy in precipitate–matrix interfaces, we utilize an efficient algorithm called the active parameter tracking (APT) algorithm [49–51] coupled with parallel coding techniques. This enables us to perform large-scale calculations, such as 20482 grid points. 2. Method 2.1. Phase-field model To represent the temporal evolution of polycrystalline material with finely dispersed precipitates we utilized the multi phase-field (MPF) model proposed by Steinbach and Pezzola [52]. In this paper, a set of continuous field variables, /1(r, t), /2(r, t), . . ., /N(r, t), are defined to distinguish the types of constituents (matrix grains and precipitates) and their orientation, where /i(r, t) represents the existence ratio of each constituent at a position r and time t. As described in Section 2.2, in order to avoid the coalescence of constituents that have the same field number i, we apply a unique number to each constituent (i.e. N is assumed to be the total number of constituents). We derive governing equations according to Ref. [53] as follows: The sum of phase fields at any position in the system is conserved. N X

/i ðr; tÞ ¼ 1

ð1Þ

i¼1

The free energy functional of a system of volume V is given by !# Z " N X P T f þ f þ kL /i  1 dV ð2Þ F ¼ V

i

where kL is the Lagrange multiplier accounting for the constraint Eq. (1) and fP can be defined by " # N X N X 2ij P f ¼  r/i  r/j þ xij /i /j ð3Þ 2 j>i i as in Steinbach and Pezzolla [52], where ij is the gradient energy coefficient and xij is the height of the double-well potential. The parabolic double-well potential /i/j is defined in the boundary region only where 0 < /i < 1 and 0 < /j < 1. For the thermodynamic potential f T, we assume a rule of mixture fT ¼

N X

/i f i ðci Þ

ð4Þ

i

where fi(ci) is the free energy density of phase i with composition ci. Any point (x, y) in the system is assumed to be a mixture of N phases, with the fraction of /i for phase i. However, the compositions of the phases coexisting at a given point are not independent of each other, but are constrained by the following condition [54]:

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

fc11 ½c1 ðr; tÞ ¼ fc22 ½c2 ðr; tÞ ¼    ¼ fcNN ½cN ðr; tÞ  fc ðr; tÞ

ð5Þ

where subscripts ci with the free energy densities denote the derivatives with respect to ci. This condition ensures that the coexisting phases have an equal chemical potential, which is defined by the difference in the chemical potentials between a solute and a solvent. With the imposed equal chemical potential condition Eq. (5), the extra chemical double-well potential at the boundary region can be removed, resulting in the relaxation of the limitation on boundary width. The average composition of a mixture is also given by the mixture rule N X cðr; tÞ ¼ /i ci ð6Þ i

With the restrictions of Eqs. (1), (5) and (6), the compositions ci (i = 1, 2, . . ., N) of each phase can be expressed as functions of c and /i. We define a step function si(r, t) = 1 if /i > 0 and si(r, t) = 0 otherwise. Then, the number of phases coexisting in a given point can be written as nðr; tÞ ¼

N X

si ðr; tÞ

ð7Þ

i¼1

Following Ref. [52], we introduce the interface field wij as wij ¼ /i  /j

ð8Þ

which is defined only within the boundary region where si 6¼ 0 and sj 6¼ 0. We can then write an identity ! N 1 X /i ¼ sij wij þ 1 ð9Þ n j6¼i with the notation sij = sisj. It then follows that   N N owij 1 X o/i 1 X o/i o/j ¼ ¼  sij sij n j6¼i n j6¼i ot ot ot ot

potential condition (Eq. (5)), the second term on the right-hand side of Eq. (12) becomes N N X X dF T ocj ocj ¼ f i ðci Þ þ /j fcjj ðcj Þ ¼ f i ðci Þ þ fc /j d/i o/i o/i j j

ð14Þ Differentiating Eq. (6) with /i yields ci þ

N X j

/j

ocj ¼0 o/i

ð15Þ

and we obtain dF T ¼ f i ðci Þ  ci fc d/i

ð16Þ

Using Eqs. (10)–(16), we find the phase-field equation " # N o/i 2 X dF dF ¼ sij M ij  ð17Þ N j6¼i d/i d/j ot where

" # N X 2ij 2 dF r /j þ xij /j þ f i ðci Þ  ci fc ¼ d/i 2 j6¼i

ð18Þ

with xij = xji and ij = ji. It should be emphasized that Eq. (18) was derived under the condition of the equal chemical potential among the coexisting phases, given by Eq. (5). The diffusion equation, where the mass should be conserved, is written as oc dF ¼ rM c r ot dc

ð19Þ

where Mc is the mobility of the diffusion field. The derivative can be modified into ð10Þ

Note that the change in the phase field takes place only in the boundary region where sij is not vanishing. Eq. (10) indicates that the phase-field change at a given point can be decomposed into a sum of interface field changes at the point [52]. We then make the relaxation ansatz " # o/i o/j dF dF  ¼ 2M ij  ð11Þ d/i d/j ot ot With definitions (3) and (4), the functional derivatives in Eq. (11) are given by dF dF P dF T ¼ þ  kL : d/i d/i d/i

6883

ð12Þ

The first term on the right-hand side of Eq. (12) is given by " # N 2ij 2 dF P X r /j þ xij /j ¼ ð13Þ d/i 2 j6¼i Eqs. (7)–(13) are the same as in Ref. [52]. With the thermodynamic potential (Eq. (4)) and the equal chemical

N N X dF X oci oci ¼ ¼ fc ¼ fc /i fci ðci Þ /i ð20Þ dc oc oc i i P where we used Eqs. (2) and (5) and /i oci =oc ¼ 1 from Eq. (6). If we put Mc = D/fcc, where D is the diffusivity dependent on the phase fields, the diffusivity in a bulk phase can be kept constant and the diffusion equation (19) becomes N oc D D X ¼ r rfc ¼ rDrc þ r  fc/ r/i ot fcc fcc i

ð21Þ

With the relationship fc/i =fcc ¼ ci , which can be easily derived from Eqs. (5) and (6), we find N X oc ¼ r  Drc  r  D ci r/i ot i

ð22Þ

By using Eq. (6) again, this equation may be written in an equivalent, but more compact and convenient, form: N X oc ¼rD /i rci ot i

ð23Þ

6884

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

For the purpose of numerical simulation, the set of phasefield equations (Eq. (17)) and the diffusion equation (23) must be solved numerically by discretizing them in space and time. The second-order central difference method and the simple explicit Euler equation are used for the discretization with respect to space and time, respectively. 2.2. Active parameter tracking algorithm In MPF models, the momory consumption increases linearly with the number of field variables, N. Furthermore, the computational cost is mainly proportional to N2. In the standard approach, we have to reduce the value of N to reduce the computational costs. However, there are two immediate drawbacks in the use of a limited number of field variables. Since each field variable represents a particular grain orientation, intermediate values of orientations are excluded. This may be viewed as restricting the allowed grain orientations to a finite discrete set of angles (with a loss of rotational invariance of the free energy) or each order parameter as representing a range of orientations. Furthermore, another undesirable consequence of limiting the number of field variables is the coalescence of grains during grain growth. Coalescence is the situation in which two grains, which have the same field variables, come into contact and instantaneously form a single large grain. This leads to incorrect growth rates. Krill and Chen [20] modified Fan and Chen’s phase-field algorithm in order to perform 3D simulations with a small number of field variables by applying the dynamic grainorientation reassignment scheme. However, their calculations were limited to isotropic grain growth only. Suwa et al. extended Krill and Chen’s scheme to an anisotropic system [55] and a system containing finely dispersed second-phase particles [46] by introducing the distance between neighboring grains and utilizing parallel coding techniques. In spite of their efforts, the available number of grain orientations is still limited. More recently, Vedantam and Patnaik [49] devised an efficient new algorithm called the active parameter tracking (APT) algorithm for solving MPF equations numerically. Gruber et al. [50] and Kim et al. [51] also devised essentially the same algorithm. The APT algorithm does not engender any additional computational expense associated with increasing the value of N. In the standard approach, the values of all the field variables at every grid point are stored and evolved at every iteration step. However, at any given point away from the grain boundaries, only one field variable is active (has a nonzero value). And near a grain boundary a small number of field variables corresponding to the adjacent grains are active. Even at the junction of the boundaries between several adjacent grains, the number of active field variables is small. Therefore, in the APT algorithm, we evolve only the active field variables instead of evolving all the field variables at every grid point. In this paper, we utilized the APT algorithm for the temporal evolution of {/i(r, t)} (i = 1 . . . N). In order to further

expedite the calculations, we coupled the APT algorithm with parallel coding techniques. The parallelization of the APT algorithm is quite easy compared with that of the dynamic grain-orientation reassignment, when we employ the second-order central difference method and the simple explicit Euler scheme for discretization with respect to space and time. The number of active field variables is used as the value of n(r, t) because only the active variables are maintained at each grids point. The total amount of memory consumption is approximately proportional to the maximum number of n(r, t), nmax. In order to reduce the value of nmax we introduce the threshold value, /th, for {/i(r,t)} (i = 1 . . . N). Further, if the value of /i(r, t) becomes smaller than that of /th, /i is forced to be inactive at position r and time t. The value of /th is set to be 1.0 · 1040. The use of the parabolic double-well potential in Eq. (3) instead of a quadratic one enables us to represent grain boundaries precisely in space with smaller grid points and also decrease the value of nmax. The value of nmax is assumed to be 12. 2.3. Model parameters and simulation procedure All calculations are performed on a 2D lattice with a periodic boundary condition. For simplicity, it is assumed that there exist only two groups of orientations (A (k = 1) and B (k = 2)) for matrix grains and one for precipitates (P (k = 3)). If two constituents have the same value of k, they belong to the same class of constituents (A, B or C), and all constituents (i = 1 . . . N) belong to one of the constituent classes (k = 1, 2, 3) (hereafter we use the word ‘‘phase’’ instead of ‘‘class of constituent’’ for simplicity). As mentioned in the introduction, we consider the system with a single texture component (A) in a matrix of randomly orientated grains (B). Further, we assume that the precipitates (P) have a special orientation relationship with A, where the difference in the properties of the A/A, A/B and B/B boundaries is not introduced to emphasize the effects of inverse pinning. The target materials are not specified in this paper. However, in order to estimate the driving forces in an actual unit (J mol1), the physical properties are assumed as follows: The diffusivity is set to be D = 1.0 · 1014 m2 s1 and phase-field mobilities are set to be {Mkl} (k, l = 1, 2, 3) = M/ = 4.0 · 108 m3 J1 s1. They are assumed to be constant values. Physical grain boundary mobility is calculated to be Mp = 2.27 · 1014 m4 J1 s1 [51] by using other physical properties described in this section. Further, the solute drag effect between matrix grains is not taken into account. Because the coarsening of precipitates is not the main subject of this paper, the value of D is changed to zero at t = 25 s to suppress a decrease in the area fraction of precipitates suffering from the Gibbs– Thomson effect. In order to simplify the calculation of the compositions ci (i = 1, 2, . . ., N) with the restrictions of Eq. (5), the

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

6885

chemical free energy density of each phase is approximated as parabolic functions [56]: f i ¼ K k ðci  ak Þ

2

ðk ¼ 1; 2; 3Þ

ði ¼ 1; 2; . . . ; N Þ

ð24Þ

where the coefficient, Kk, which determines the stability of each phase, is set to be Kk = 1.0 · 103 (k = 1, 2, 3) J mol1; the equilibrium composition of each phase is ak = 0.1 (k = 1, 2) and a3 = 0.9, respectively. The maximum value of boundary energy is rmax = 1.0 J m2 and the molar volume is 7.0 · 106 m3 mol1. The lattice step size Dx is set to be 1.0 · 107 m and the boundary width is assumed to be 7 · Dx. Hereafter, length and area are expressed in m and m2, respectively. Then, the parameters kl and xkl are calculated as kl = 4/p(3.5Dxrkl)0.5 and xkl = 2rkl/3.5Dx [54], where rkl is the boundary energy between phases k and l. A step for time integration, Dt, of 2.5 · 102 s is employed. The system size of (5.12 · 105)2 is employed in Sections 3.1–3.3 and that of (2.048 · 104)2 is employed in Sections 3.4 and 3.5. All calculations in this paper have been performed on the Numerical Materials Simulator (HITACHI SR11000) at the National Institute for Materials Science (NIMS) with two nodes (number of central processing units: 16 · 2 = 32, memory for application codes: 24 · 2 = 48 GB). The longest simulation time was about 2 days for 6.0 · 105 steps (1.5 · 104 s). 3. Results and discussion 3.1. Bicrystalline system with plate-shaped precipitates Nishizawa [30] proposed the driving force due to inverse pinning for platelet precipitates to be DGinv-p  ð2f p V m =dp ÞðrBP  rAP Þ ðJ mol1 Þ

ð25Þ

where fp is the volume fraction of the precipitates and dp is the thickness of the precipitate. If we adjust the value of fp by changing the distance between precipitates, d, with dp kept constant, Eq. (25) can be rewritten as DGinv-p  ð2V m =dÞðrBP  rAP Þ

ðJ mol1 Þ

ð26Þ

where the condition dp  d is assumed. In this section we investigate the effect of the value of (rBP  rAP) and d on the driving force of inverse pinning for plate-shaped precipitates. First, in order to investigate the effect of boundary energy difference, (rBP  rAP), we have performed the simulation under the following boundary energy conditions: CASE1 (rBP = rPB = 1.0 J m2 and rkl = 0.6 J m2 (k, l = 1, 2, 3)) and CASE2 (rBP = rPB = 1.0 J m2 and rkl = 0.4 J m2 (k, l = 1, 2, 3)) with the value of d kept constant at d = 6.4 · 106. The temporal evolutions of the calculated microstructure are shown in Fig. 1. Since the simulations are performed in 2D, the precipitates are considered to be the plates having infinite thickness in the third dimension. In this figure, the black boundary represents the high-energy

Fig. 1. Temporal evolutions of simulated microstructure for bicrystalline system with plate-shaped precipitates. In order to investigate the effect of boundary energy difference, (rBP  rAP), we performed the simulation under the following conditions: CASE1 (rBP = rPB = 1.0 J m2 and rkl = 0.6 J m2 (k, l = 1, 2, 3)) and CASE2 (rBP = rPB = 1.0 J m2 and rkl = 0.4 J m2 (k, l = 1, 2, 3)), with the value of d kept constant at d = 6.4 · 106 [m].

boundary and white one represents the low-energy boundary. Because the boundary energy is normalized at each simulation run, the darkness of a boundary does not reflect the absolute value of boundary energy. This figure shows the continuous migration of interface AB in both the conditions. This indicates that the condition S > 0 is not required for its migration, where the spreading coefficient, S, is defined as S = rBP  rAB  rAP and the condition S > 0 implies the perfect wetting condition [57,58]. Next, we have investigated the effects of distance, d, with CASE1. The position of interface AB is determined by the area of A. Then, the mean velocity of interface AB, ÆvABæ, is calculated by using the temporal evolution of its position. For the purpose of image analysis, a function O(r, t) is defined as follows: when /q(r, t) has the maximum value among all the field variables at a grid point, r, O(r, t) is set to be the constituent number, q, and the value of O(r, t) is assigned to all the grid points in the system. Because the perfect wetting of precipitates by grain A makes the evaluation of area difficult, we have selected CASE1. As mentioned above, the value of ÆvABæ is calculated by image analysis. Further, the driving force due to inverse pinning is estimated from the relation hvAB i ¼ M p;AB DGinv-p =V m

ðm s1 Þ;

ð27Þ

where Mp,AB is the physical mobility of interface AB. The values of DGinv-p for various values of d calculated from simulation results are plotted in Fig. 2. For comparison, the values estimated using Eq. (26) are also plotted. This figure indicates that the simulated values are in good agreement with the estimated values if the value of d is sufficiently large. The increase in discrepancy with d is

6886

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

Fig. 2. Driving force, DGinv-p, for various values of d calculated from simulation results (solid circles). For comparison, the values estimated by Eq. (26) are also plotted (open circles and solid line).

attributed to two main reasons. One reason is the assumption dp  d in Eq. (26), the other is that the kinetics of triple junction among grains A and B and precipitate P are not taken into account in Eqs. (25) and (26) [57,58]. 3.2. Bicrystalline system with regularly arranged circular precipitates In the preceding section, we compared the simulated driving force of inverse pinning to that estimated from Eq. (26) by assuming plate-shaped precipitates. In this and the next section, bicrystalline systems containing many dispersed precipitates are considered. For simplicity, the shape of the precipitates is assumed to be circular. Since the simulations are performed in 2D, the precipitates are considered to be circular cylinders having infinite thickness in the third dimension. The configurations of each constituent are described in Fig. 3 and the detailed information about the configurations is illustrated schematically in Fig. 4. We shall consider systems with regularly arranged

circular precipitates in this section. In Figs. 3 and 4, the lengths a and b indicate the center-to-center spacing between the precipitates that are parallel to AB and the precipitates that are perpendicular to AB, respectively; r is the radius of precipitates and lint is the boundary width. First, in order to investigate the effect of boundary energy difference, (rBP  rAP), we performed the simulation under the following conditions: CASE1 and CASE2, with the values of r, a and b kept constant at r = 4.4 · 107, a = 6.4 · 106 and b = 1.6 · 106. Fig. 5 shows that the continuous migration of interface AB is not observed with CASE1. In the preceding section, the perfect wetting condition, S > 0, is not required for the continuous migration of interface AB. However, in the case of the bicrystalline system with finely dispersed precipitates, the condition S > 0 is required for the migration across the row of precipitates because all BP interfaces should divide into AB and AP interfaces for the continuous migration of interface AB. Furthermore, we have confirmed the critical distance blim, and if b exceeds blim, the migration across the row of precipitates cannot be observed irrespective of the boundary energy condition, where blim is estimated as blim = 2(r + lint). It is worth noting that blim is determined not only by physical requirements but also by model requirements (i.e. the diffuseness of an interface in the phase-field model). Therefore, it is highly dependent on the boundary width, lint. The value of blim can be estimated as blim = 2r, provided that the interface has zero thickness and there is no long-range interaction between grain boundary and precipitate. In this section we utilized circular precipitates for simplicity; however, we would have a critical distance, blim, even if different shaped precipitates were utilized.

Fig. 3. Configurations of each constituent for the bicrystalline systems with circular precipitates. To save space, only the left parts of a complete system are presented for each configuration. The lengths a and b indicate the center-to-center spacing between the precipitates that are parallel to AB and the precipitates that are perpendicular to AB, respectively; r is the radius of precipitates and lint is the boundary width.

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

6887

b

a

2r+lint

Grain A Grain B

lint

f p = π r 2 / ab

blim = 2(r + 0.5lint ) + lint Fig. 4. Precise description of the configuration. To save space, only a part around interface AB is presented. The lengths a and b indicate the centerto-center spacing between the precipitates that are parallel to AB and the precipitates that are perpendicular to AB, respectively; r is the radius of precipitates and lint is the boundary width.

Fig. 5. Temporal evolutions of simulated microstructure for bicrystalline system with circular precipitates. In order to investigate the effect of boundary energy difference, (rBP  rAP), we performed the simulation under the following conditions: CASE1 (rBP = rPB = 1.0 J m2 and rkl = 0.6 J m2 (k, l = 1, 2, 3)) and CASE2 (rBP = rPB = 1.0 J m2 and rkl = 0.4 J m2 (k, l = 1, 2, 3)), with the values of r,a and b kept constant at r = 4.4 · 107, a = 6.4 · 106 and b = 1.6 · 106.

Hereafter, b is kept constant at b = 16, which is smaller than blim for the various values of r used in the following calculations. In order to investigate the effect of the area fraction of precipitates, fp, we have performed the simulations with various values of a and r. The values of a = 3.2, 6.4, 12.8, 25.6 · 106 and r = 2.8, 4.4, 5.3, 6.6 · 107 are considered. The normalized velocities of interface AB with the simulation time for various values of a are shown in Fig. 6, where r is kept constant at r = 4.4 · 107. A time cycle for the oscillation of velocity, tv(a), corresponds to the time in which the interface moves by the distance b. Therefore, the mean velocity, ÆvABæ, is calculated as

Fig. 6. Normalized velocities of interface AB with the simulation time for various values of a. In these calculations, r is kept constant at r = 4.4 · 107. A time cycle for the oscillation of velocity, tv(a), corresponds to the time in which the interface moves by the distance b.

ÆvABæ = b/tv and the averaged driving force is estimated by using Eq. (27). In Fig. 7, the relation between fp and the averaged driving force is plotted for the various values of r. For comparison, the relation for plate-shaped precipitates is also plotted. Each relation is fitted to a relation of type DGinv-p = a Æ fb, where a and b are the fitting parameters. This figure indicates that the value of b is highly dependent on the precipitate shape and Eq. (26) is not applicable to the system containing dispersed circular precipitates. In the system with the plate-shaped precipitates, the morphology of interface AB does not change with d (i.e. fp). Therefore, its curvature is inversely proportional to d; further, b becomes 1, provided that d is sufficiently larger than dp. However, in the system containing the circular precipitates, the self-similarity of the interface is not maintained when a (i.e. fp) is varied as shown in Fig. 3. This is the main reason for the difference in the values of b. Further, b deviates from 2 when r = 2.8, 6.6 · 107. Although the reason for this is not completely clear, at least the value of r = 2.8 · 107 is not sufficiently large to obtain realistic grain boundary–precipitate interactions. On the other hand, for r = 6.6 · 107, the actual distance between precipitates becomes very small and this may affect the shape of the interface AB. Therefore, precipitates with radius r = 4.4 · 107 are used in the following sections.

Fig. 7. Relation between the driving force and fp for various values of r. For comparison, the relation for plate-shaped precipitates is also plotted. Each relation is fitted to a relation of type DGinv-p = a Æ fb, where a and b are the fitting parameters.

6888

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

Fig. 8. Temporal evolution of interface AB with Np = 300. This figure shows the continuous migration of interface AB.

3.3. Bicrystalline system with irregularly arranged circular precipitates In this section, we consider bicrystalline systems with irregularly arranged circular precipitates. We performed the simulation under CASE2 with r kept constant at r = 4.4 · 107. The precipitates were overlapped with a bicrystalline structure at the beginning of the simulation, such that they were separated from each other by a preset minimum distance, llim; llim was set to be llim = 2.4 · 106, which is larger than blim for the precipitates with r = 4.4 · 107. The total number of precipitates, Np, between 75 (fp = 0.020) and 300 (fp = 0.079) has been considered. In order to even out the data spread, three simulation runs were performed for each Np value, with each run having a different configuration of the precipitates. Fig. 8 depicts the temporal evolution of interface AB with Np = 300. The temporal evolution of the area fraction of A grain for various values of Np is shown in Fig. 9. This figure indicates that the continuous migration of interface AB is possible in the system with Np = 300, as shown in Fig. 8. Thus, the averaged driving force for the migration is estimated as 0.17 J mol1 by the linear approximation of the curve. In contrast, in the systems with Np = 150 and Np = 75, we could observe no continuous migration of interface AB. In order to discuss the limit of the interface migration, we define the distance, lmin(r), as follows: at each point, r, on interface AB, you can find the nearest precipitate in grain B. lmin(r) is defined as the distance between the center of the nearest precipitate and the point r. Further, interface AB can migrate until all its points simultaneously satisfy the conditions lmin(r) > r + lint and j(r) = 0, where j(r) is the mean curvature of interface AB at r. Note that the

Fig. 9. Temporal evolution of the area fraction of A grain, fA, for various values of Np (open symbols). For comparison, the temporal evolutions of the area fraction of precipitates, fp, for various values of Np (solid symbols) are also shown. The values of fp remain constant, as expected. The data points were averaged over the results of three separate runs at the indicated times.

migration speed at each point on the interface is proportional to j(r). Therefore, it is difficult to detect the migration of the interface with a sufficiently small j(r), even if j(r) is not zero. 3.4. Polycrystalline system with circular precipitates So far, we have considered bicrystalline systems only. In this and the next section, we shall consider polycrystalline systems with finely dispersed circular precipitates in order to investigate the effects of an additional driving force due to normal grain growth. All calculations have been performed on a 2D lattice of size (2.048 · 104)2. Initial configurations were constructed as follows: first, a polycrystalline microstructure was created by simulating normal grain growth by assigning one group of orientation for the entire system. The initial number of grains was 19,334. Twenty grains were randomly assigned to group

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

A; others were assigned to group B. Subsequently, precipitates were overlapped such that they separated from each other by a preset minimum distance, llim. Np was between 1250 (fp = 0.019) and 5000 (fp = 0.077). In Section 3.2, the migration of interface AB across the row of precipitates was not observed if the condition S > 0 was not satisfied. Therefore, we selected CASE2 at the beginning. In this case, the growth of B grains reaches its steady state during the early stages of simulation. Fig. 10 depicts the microstructural evolution of the polycrystalline system with Np = 5000, llim = 2.4 · 106 and r = 4.4 · 107. This figure indicates that the grains belonging to group A grow abnormally at the expense of those belonging to group B due to inverse pinning. On the other hand, the growth of B grains is prohibited, and the precipitates, which are located at the interface between A grains, also pin the interface. This is because, in these two cases, there is no additional driving force due to inverse pinning. In 2D systems, the grain boundaries have a tendency to become straight between the precipitates [44]. Thus, it is difficult for the grain boundaries to escape from the precipitates since all of the driving pressure for grain boundary movement is removed once a grain boundary is straight. This is the usual Zener pinning phenomenon in 2D systems. It is worth noting that an abnormally growing grain initiated by the differences in boundary properties has a limiting grain size even if one growing grain is embedded in the surrounding matrix [24,25]. On the other hand, an abnormally growing grain due to inverse pinning appears to have no limiting size, provided the number density of growing grains is sufficiently small. In Fig. 11, the distribution of n(r, t) is plotted against the simulation time. This figure indicates that the fraction of

6889

Fig. 11. Temporal evolution of the distribution of n(r, t). This figure indicates that the fraction of grid point, r, with n(r, t) = 1, increases with simulation time.

grid points, r with n(r, t) = 1, increases with the simulation time. This is attributed to the decrease in the number of triple junctions and the area of a boundary region. Therefore, after the system has reached its steady state (i.e. t > 10,000 s), we can observe no further evolution of the distribution. Further, the value of 8 for nmax might be sufficiently large and nmax should be microstructure-dependent (i.e. time-dependent) in order to perform a simulation with both efficiency and accuracy. Fig. 12 depicts the microstructural evolution of the polycrystalline system with Np = 1250, llim = 4.8 · 106 and r = 4.4 · 107. In this case, the growth of B grains does not reach a steady state during the early stages of simulation. Thus, almost all A grains are unable to grow because of their small number density. However, some A grains start to grow abnormally after the growth of B grains reaches its steady state. Note that, in the preceding section, we could not observe the continuous migration of interface AB if fp was not sufficiently large (i.e. fp 6 3.9). However,

Fig. 10. Microstructural evolution of the polycrystalline system with Np = 5000, llim = 2.4 · 106 and r = 4.4 · 107. This figure indicates that the grains belonging to group A grow abnormally at the expense of those belonging to group B.

6890

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

Fig. 12. The microstructural evolution of the polycrystalline system with Np = 1250, llim = 4.8 · 106 and r = 4.4 · 107. In this figure, the growth of B grains does not reach its steady state at an initial stage. Thus, almost all A grains are unable to grow because of their small number density. However, some A grains start to grow abnormally after the growth of B grains reaches its steady state.

in the polycrystalline systems, abnormal grain growth can occur with the aid of the driving force for grain growth once the growth of B grains reaches its steady state. The driving force for the grain growth may be proportional to the difference in the radius between a growing grain (belonging to group A) and the surrounding grains (belonging to group B) under an isotropic boundary energy condition [26]. If the mean radius of the surrounding B grains, ÆrBæ, remains constant and the radius of a growing grain is sufficiently larger than the value of ÆrBæ, the driving force is approximated to agrBBVm/ÆrBæ, where rBB is the boundary energy between B grains and ag is a dimensionless constant [24,25]. From our simulation results, the driving force is calculated as a value between 1 and 2 J mol1, depending on the value of ÆrBæ, where the ag is assumed to be unity. This calculated value is larger than the value due to inverse pinning obtained in the preceding section. This implies that the dominant driving force for the abnormal grain growth is not the driving force due to inverse pinning itself but, rather, the driving force for grain growth, at least at a later stage of the abnormal grain growth. The temporal evolutions of ÆrBæ for various numbers and configurations of precipitates are shown in Fig. 13. This figure indicates that the simulation time required for reaching the steady state becomes smaller as fp increases and the randomness of the configuration decreases. In the systems with 5000 precipitates, the decrease in the value of ÆrBæ is attributed to the significant decrease in the number of B grains, and B grains have been completely annihilated in the final configuration. Next, the temporal evolutions of fA for a various number and configurations of precipitates are shown in Fig. 14. This figure indicates that fA increases

Fig. 13. Temporal evolutions of the mean radius of B grains, ÆrBæ, for various number and configurations of precipitates. In systems with 5000 precipitates, the decrease in the value of ÆrBæ is attributed to the significant decrease in the number of B grains, and B grains have been completely annihilated in the final configuration.

at a faster speed as fp and the randomness of a configuration increase. In order to clarify the effect of configuration, the microstructural evolution of the polycrystalline system with Np = 5000, llim = 1.2 · 106 and r = 4.4 · 107 is shown in Fig. 15. The shape of abnormally grown A grains becomes irregular (far from circular shape) compared with those shown in Fig. 10. This fact can be interpreted as follows: the driving force for the abnormal growth of A grains is the reduction of the total boundary energy of the system. For the system with homogeneously distributed precipitates, the growing directions of A grains have little effect on the reduction speed of the total boundary energy. Therefore, the A grains maintain an almost circular shape. On the other hand, for the system with inhomogeneously

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

6891

this section enlarge the possibility of inverse pinning in real alloy systems. We also note that, in 3D systems, the condition S > 0 may not be required for the continuous migration of interface AB even in the bicrystalline systems, considering the geometry of the interaction between the interface and a spherical particle [44]. 3.5. Polycrystalline system with plate-shaped precipitates

Fig. 14. Temporal evolutions of the area fraction of A grains, fA, for various number and configurations of precipitates.

distributed precipitates, the directions have a significant effect on the reduction speed; the system selects the directions, at every time step, in which the total boundary energy can decrease most rapidly. This leads to irregularshaped A grains. Finally, we show the microstructural evolution of the polycrystalline system under CASE1 in Fig. 16, where Np = 5000, llim = 2.4 · 106 and r = 4.4 · 107. In this figure, we can observe at least the abnormal growth of A grains, even though the growth speed is noticeably slower than that observed in Fig. 10. In Section 3.2, the perfect wetting condition, S > 0, is required for the migration of interface AB across the row of precipitates. On the other hand, abnormal grain growth can occur with the aid of the driving force for grain growth even though the condition S > 0 is not satisfied. In real alloy systems, it appears difficult to satisfy the condition corresponding to S > 0. The facts validated in

In the preceding section, we considered the polycrystalline systems containing the circular precipitates. In this section, let us consider the polycrystalline systems with plateshaped precipitates in order to investigate the effects of the shape of precipitates on inverse pinning. Precipitates with an aspect ratio of 4 are utilized on the condition that the total number and area fraction of the precipitates and the center position of each precipitate are equal to those in the system containing the circular precipitates with the conditions Np = 5000 and llim = 2.4 · 106 (hereafter we refer to this case as ‘‘CIRC’’); moreover, CASE2 is also utilized. Two types of orientation configurations of the plate-shaped precipitates are considered: (1) randomly oriented precipitates; and (2) unidirectionally aligned precipitates. In the case of (1), the angle between the long axis of the each precipitate and the x-axis of the system is randomly assigned. On the other hand, in the case of (2), the angle between the long axis of all precipitates and the x-axis is equal to 3/8p. This value is selected to just equal fp in the case of CIRC; because the initial configuration is digitized at each finite difference grid point, fp is slightly dependent on the value of the angle. It should be noted that, in this section, only the difference in the shapes of the precipitates is considered; the boundary energy dependence on the inclination of the

Fig. 15. Microstructural evolution of the polycrystalline system with Np = 5000, llim = 1.2 · 106 and r = 4.4 · 107. The shape of abnormally grown A grains becomes irregular (far from circular shape) compared with those shown in Fig. 10.

6892

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

Fig. 16. Microstructural evolution of the polycrystalline system with the boundary energy condition CASE1 (rBP = rPB = 1.0 J m2 and rkl = 0.6 (k,l = 1,2,3)), where Np = 5000, llim = 2.4 · 106 and r = 4.4 · 107.

boundaries (i.e. the angle between the boundaries and the x-axis of the system) is not considered. Fig. 17 depicts the microstructural evolutions of the polycrystalline system with the orientation condition (1) shown in the upper column and that with the condition (2) shown in the lower column. In both the orientation conditions, the abnormal growth of A grains is observed, although an abnormally grown A grain with condition (2) is slightly elongated in the direction of the long axis of the precipitates.

The temporal evolutions of ÆrBæ for various orientation configurations of the precipitates are shown in Fig. 18; for comparison, that in the case of CIRC is also shown in this figure. This figure indicates that the plate-shaped precipitates more effectively pin the B grains. Although this can be attributed to the difference in the total lengths (area in 3D) of the precipitate–matrix interface for the plateshaped and circular precipitates [36], further investigation is needed to clarify the effects of the precipitates’ shape on the growth kinetics.

Fig. 17. Microstructural evolutions of the polycrystalline system with the plate-shaped precipitates. Orientation conditions (1) and (2) are shown in the upper and lower columns. In the first row of each column, only the upper left parts of complete structures are presented to show the exact shape and direction of each precipitate.

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894



Fig. 18. Temporal evolutions of the mean radius of B grains, ÆrBæ, for various orientation configurations of precipitates (i.e. (1) randomly oriented and (2) unidirectionally aligned precipitates). For comparison, that in the case of CIRC (the system containing the circular precipitates with the conditions Np = 5000 and llim = 2.4 · 106) is also shown.





Fig. 19. Temporal evolutions of the area fraction of A grains, fA, for various orientation configurations of precipitates (i.e. (1) randomly oriented and (2) unidirectionally aligned precipitates). For comparison, that in the case of CIRC (the system containing the circular precipitates with the conditions Np = 5000 and llim = 2.4 · 106) is also shown.

Finally, the temporal evolutions of fA for various orientation configurations of precipitates are shown in Fig. 19; for comparison, that in the case of CIRC is also shown. The system with the orientation condition (2) exhibits the smallest increasing speed of fA in the middle of the simulation run. This can be attributed to the randomness of the orientation of the precipitates, as is in the case described in the previous section. Meanwhile, the slowing down of the increase in speed during the later stage is due to the impingement among the abnormally grown A grains. At this stage, the choice of growing direction becomes limited; the effects of the randomness in the orientation of the precipitates may be smaller. Further, in all cases, the center positions of all precipitates are equal. Therefore, the time for reaching the steady state becomes almost equivalent in all cases. 4. Concluding remarks The possibility of inverse pinning in 2D systems has been investigated by using the multi phase-field (MPF) model. The following results were obtained.  In this study, to realize the simulation with the anisotropy in precipitate–matrix interfaces, we utilized an efficient algorithm, the active parameter tracking (APT)



6893

algorithm, coupled with parallel coding techniques. This enabled us to perform large-scale calculations such as 20482 grid points. We have validated the driving force for inverse pinning proposed by Nishizawa with the bicrystalline system containing plate-shaped precipitates. The simulated values are in good agreement with the estimated values if the distance between precipitates, d, is sufficiently large. The increase in discrepancy with decreasing d is attributed to two main reasons. First is the assumption dp  d in Eq. (26) and secondly, the kinetics of triple junction among grains A and B and precipitate P are not taken into account in Eqs. (25) and (26). Bicrystalline systems containing many dispersed precipitates have been considered. The perfect wetting condition, S > 0, is required for the migration of interface AB across the row of precipitates. Furthermore, we have confirmed the critical distance blim, and if b exceeds blim, the migration across the row of precipitates is not observed. The value of blim is determined not only by the physical requirements but also by the model requirement. Thus, it is highly dependent on the boundary width, lint. Polycrystalline systems with finely dispersed circular precipitates were considered. In these systems, abnormal grain growth can occur with the aid of the driving force for grain growth once the growth of B grains reaches its steady state. Furthermore, the perfect wetting condition, S > 0, is not required for the abnormal grain growth in these systems. These facts enlarge the possibility of inverse pinning in real alloy systems. Finally, polycrystalline systems with finely dispersed plate-shaped precipitates were considered. The simulation results indicate that the plate-shaped precipitates more effectively pin the B grains, and the system with the unidirectionally aligned precipitates shows the smallest increasing speed of the area fraction of A grains in the middle of the simulation run.

Acknowledgements This study was partially supported by the Next Generation Supercomputer Project, Nanoscience Program, MEXT, Japan. It was also supported by the CREST, Japan Science and Technology Agency. References [1] Humphreys FJ, Hatherly M. Recrystallization and related annealing phenomena. Oxford: Pergamon Press; 1996 [chapter 9]. [2] Burke J. Trans Metall Soc AIME 1949;180:73. [3] Smith C. In: Metals interfaces. Cleveland, OH: American Society for Metals; 1952. p. 5–6. [4] Burke J, Turnbull D. Prog Metal Phys 1952;3:220. [5] Mullins W. J Appl Phys 1956;27:900. [6] Mullins W. Acta Mater 1998;46:6219. [7] Frost HJ, Thompson CV. Curr Op Solid State Mater Sci 1996;1:361.

6894

Y. Suwa et al. / Acta Materialia 55 (2007) 6881–6894

[8] Fan D, Chen L-Q. Acta Mater 1997;45:611. [9] Grest GS, Anderson MP, Srolovitz DJ. Phys Rev B 1988;38: 4752. [10] Marthinsen K, Hunderi O, Ryum N. Acta Mater 1996;44:1681. [11] Anderson MP, Grest GS, Srolovitz DJ. Phil Mag B 1989;59:293. [12] Glazier JA. Phys Rev Lett 1993;70:2170. [13] Saito Y. ISIJ Int 1998;38:559. [14] Miyake A. Contrib Mineral Petrol 1998;130:121. [15] Song XY, Liu GQ. Univ Sci Tech Beijing 1998;5:129. [16] Fuchizaki K, Kusaba T, Kawasaki K. Phil Mag B 1995;75:333. [17] Weygand D, Bre´chet Y, Le´pinoux J, Gust W. Phil Mag B 1999;79:703. [18] Xue X, Righetti F, Telley H, Liebling TM. Phil Mag B 1995;75:576. [19] Wakai F, Enomoto N, Ogawa H. Acta Mater 2000;48:1297. [20] Krill CE, Chen L-Q. Acta Mater 2002;50:3057. [21] Rollett A, Srolovitz D, Anderson M, Doherty R. Acta Metall 1992;40:3475. [22] Rollett A, Srolovitz D, Anderson M. Acta Metall 1989;37:1227. [23] Grest G, Anderson M, Srolovitz D, Rollett A. Scripta Metall 1990;24:661. [24] Humphreys FJ. Acta Mater 1997;45:4231. [25] Rollett A, Mullins W. Scripta Mater 1997;36:975. [26] Hillert M. Acta Metall 1965;13:227. [27] Anderson I, Gront O, Ryum N. Acta Metall Mater 1995;43:2689. [28] Humphreys FJ. Acta Mater 1997;45:5031. [29] Messina R, Soucail M, Kubin L. Mater Sci Eng A 2001;308:258. [30] Nishizawa T. ISIJ Int 2000;40:1269. [31] Miodownik MA. Scripta Mater 2006;54:993. [32] Ushigami Y, Kumano T, Haratani T, Nakamura S, Takebayashi S, Kubota T. Mater Sci Forum 2004;467–470:853. [33] Srolovitz DJ, Anderson MP, Grest GS, Sahni PS. Acta Metall 1984;32:1429.

[34] Anderson MP, Grest GS, Doherty RD, Li K, Srolovitz DJ. Scripta Metall 1989;23:753. [35] Hassold GN, Holm EA, Srolovitz DJ. Scripta Metall 1990;24:101. [36] Kad BK, Hazzledine PM. Mater Sci Eng A 1997;238:70. [37] Gao J, Thompson RG, Patterson BR. Acta Mater 1997;45:3653. [38] Soucail M, Messina R, Cosnuau A, Kubin LP. Mater Sci Eng A 1999;271:1. [39] Miodownik M, Martin JW, Cerezo A. Phil Mag A 1999;79:203. [40] Miodownik M, Holm EA, Hassold GN. Scripta Mater 2000;42:1173. [41] Harun A, Holm EA, Clode MP, Miodownik MA. Acta Mater 2006;54:3261. [42] Weygand D, Bre´chet Y. Le´pinoux J. Acta Mater 1999;47:961. [43] Riege SP, Thompson CV, Frost HJ. Acta Mater 1999;47:1879. [44] Moelans N, Blanpain B, Wollants P. Acta Mater 2005;53:1771. [45] Moelans N, Blanpain B, Wollants P. Acta Mater 2006;54:1175. [46] Suwa Y, Saito Y, Onodera H. Scripta Mater 2006;55:407. [47] Fan D, Chen L-Q, Chen S-PP. J Am Cream Soc 1998;81:526. [48] Ohnuma I, Ishida K, Nishizawa T. Phil Mag A 1999;79:1131. [49] Vedantam S, Patnaik BSV. Phys Rev E 2006;73:016703. [50] Gruber J, Ma N, Wang Y, Rollett AD, Rohrer GS. Model Simul Mater Sci Eng 2006;14:1189. [51] Kim SG, Kim DI, Kim WT, Park YB. Phys Rev E 2006;74:061605. [52] Steinbach I, Pezzolla F. Physica D 1999;134:385. [53] Kim SG, Kim WT, Suzuki T, Ode M. J Cryst Growth 2004;261:135. [54] Kim SG, Kim WT, Suzuki T. Phys Rev E 1999;60:7186. [55] Suwa Y, Saito Y, Onodera H. Comput Mater Sci 2007;40:40. [56] Zhu JZ, Wang T, Ardell AJ, Zhou SH, Liu ZK, Chen L-Q. Acta Mater 2004;52:2837. [57] de Gennes PG. Rev Mod Phys 1985;57:827. [58] de Gennes PG, Brochard-Wyart F, Que´re´ D, Gouttes, bulles perles et ondes; translated by Okurmura K, Yoshioka Syoten; 2003 [in Japanese].