- Email: [email protected]
Anisotropic grain growth with pore drag under applied loads
Materials Science and Engineering A 412 (2005) 271–278
Anisotropic grain growth with pore drag under applied loads X.N. Jing a,b , J.H. Zhao a , G. Subhash b,∗ , X.-L. Gao c a
CAS Key Laboratory of Mechanical Behavior and Design of Materials, University of Science and Technology of China, Hefei, Anhui 230027, PR China b Department of Mechanical Engineering and Engineering Mechanics, Michigan Technological University, 1400 Townsend Drive, Houghton, MI 49931-1295, USA c Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843-3123, USA Accepted 31 August 2005
Abstract In the final stage of sintering of ceramics, residual pores co-evolve with grain boundaries because of incomplete densification. Their interactions coupled with external loads are critical to the microstructural evolution of structural ceramics. A modified two-dimensional (2D) diffuse-interface phase field model, which differs from the boundary-tracking methods, is utilized to investigate the effects of stochastically distributed pore drag on grain growth kinetics and morphological evolution process of ceramics under applied loads. Contributions from both the boundary energy and elastic strain energy caused by pore drag forces and applied loading are incorporated in the modified phase field model to describe the isotropic or cubically anisotropic behaviors of polycrystalline materials. The temporal evolution of the spatially dependent grain orientation variables is determined by numerically solving non-linear Ginzburg–Landau equations using a semi-implicit Fourier-spectral method. Numerical results show that the anisotropic strain energy dominates the non-self-similar growth manner, leads to ordered grain morphologies and changes the growth rate. © 2005 Elsevier B.V. All rights reserved. Keywords: Phase field model; Grain growth; Pore drag; Morphological evolution
1. Introduction Mechanical properties of polycrystalline materials (e.g., ceramics) depend strongly on the size, distribution and morphology of grains. Thus, how to control grain growth during processing is a critical issue for the development of polycrystalline materials [1–3]. Grain growth can be viewed as a process of grain boundary migration to decrease the total grain boundary areas and the total free energy of the material system, both of which are driven by mean curvatures of grain boundaries [1,2]. In a single-phase material, the only process which occurs during grain growth is local atomic re-arrangement. However, in a multi-phase material, the long-range diffusion dominates the migration of grain boundary. Hence, the kinetics of grain growth is strongly affected by the presence or absence of solute or impurity migration and segregation at grain boundaries. The pinning effect of pores is of great practical importance in the sintering of high-quality ceramics for structural applications, where den-
∗
Corresponding author. Tel.: +1 906 487 3161; fax: +1 906 487 2822. E-mail address: [email protected] (G. Subhash).
0921-5093/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2005.08.220
sification and small grain sizes are often required to obtain good strength and toughness. In the last two decades, many models have been proposed to predict the time dependence of average grain size and size distribution in polycrystalline materials [3–7]. Due to the complexity of topologies and coupled multiple interactions, the analytical modeling of microstructure evolution is often very difficult. As a result, computer simulation plays a key role in exploring the details of grain growth and validating the analytical models. The sharp-interface approaches, such as Vortex model [8], Potts models [9–11] and Cellular Automata method [12], have been used in simulations of grain growth. Recently, simulating grain growth using a continuum-based diffuse-interface phase field model has been developed by several researchers [13–19]. Chen and co-workers simulated grain growth in single-phase materials [13,14] and studied the microstructure evolution in volume-conserved two-phase material systems with finite interface thickness [15] by considering grains of different crystallographic orientations represented by a set of non-conserved order parameter fields. In their model, the effects of solute or precipitate drag on grain growth rate and size distribution are taken into account without ad hoc assumptions
272
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
about grain shapes. An independent phase field model was developed in [16–19], which adopted two order parameters to describe a grain structure; one denoting the crystalline order and the other describing predominantly the local orientation of the crystal. In this model, the grain rotation can be simulated by the relaxation of the local orientation order parameter. During the final-stage sintering of ceramics, most residual pores are isolated, residing on grain boundaries, three-grain junctions or four-grain junctions due to incomplete sintering process [1,2]. The volume fraction of pores decreases as atoms elsewhere diffuse along the grain boundaries into pores, and grains grow as atoms migrate between grains. Yu and Suo [20] developed an axisymmetric model to study the steady and transient motion of pore–grain boundary separation using the finite element method. Wang and Li [21] further considered the shrinkage of grain boundary pores under pressure. However, it appears that the effects of pore shrinking and pinning on the grain growth process during the final-stage sintering have not been fully discussed. To progress from this situation, the effects of the Zener pinning of randomly non-conservative pores moving with grain boundaries on grain growth rate and anisotropic growth mechanism are studied from the energy change point of view using a diffuse-interface phase field model reported in [13,14]. 2. Phase field model description In the diffuse-interface field model developed by Fan and Chen [13,14], the microstructure of a polycrystalline material specimen is described by using a set of p continuous nonconserved long-range order parameters {ηi (r,t)} (i = 1, 2, . . ., p) defined at a given time t at each position r within the specimen, i.e. η1 (r, t), η2 (r, t), . . . , ηp (r, t),
(1)
where p is the number of possible orientations of a grain in space. In real materials, the number of orientations is infinite (i.e., p → ∞), but a finite value of p (p > 32) is often sufficient for reasonably accurate simulations of grain growth. ηi (r,t) (i = 1, 2, . . ., p) serve as orientation field variables to distinguish different crystallographic orientations of grains and are taken to be continuous in space. When the total free energy is minimized, only one field variable within each grain takes on a value of 1 or −1, and all other field variables are zero-valued. Across a grain boundary between two adjacent grains, ηi changes continuously between 1 and 0 or −1 and 0. According to the diffuse-interface theories [22], the total free energy F of the material specimen can be described as the function of a set of field variables: p κi 2 F= f0 (η1 (r), η2 (r), . . . , ηp (r)) + dV, [∇ηi (r)] 2 V i=1 (2) where f0 is the local free energy density as a function of ηi (r), ∇ηi (r) is the gradient of ηi (r), κi are the gradient energy coefficients which are positive constants, and V is the volume of the specimen. The gradient term is related to the boundary energy.
For a planar grain boundary between a grain of orientation i and another grain of orientation j, the grain boundary energy γ gb is given by [13]: ∞ κi dηi 2 κj dηj 2 γgb = f (ηi , ηj ) + dx, (3) + 2 dx 2 dx −∞ where f(ηi ,ηj ) = f0 (0, . . ., ηi , ηj , . . ., 0) − fmin , with fmin being the minima of f0 . In order to model the temporal and spatial evolution of the material microstructure, an approximate Landau type free energy density is proposed as [13,14]: p 2 p p p A 1 2 A2 2 A2 2 2 f0 = − ηi + ηi + A3 − ηi η j , 2 4 2 i=1 i=1 i=1 j=i (4) where A1 , A2 , A3 are phenomenological constants whose values are determinable from grain boundary energy and grain boundary width. This free energy density has 2p potential wells in the p-field space, which represent the equilibrium free energies of crystalline grains in 2p different orientations. The kinetic equations that describe microstructure evolution with time and space are given by the time-dependent Ginzburg–Landau equations (e.g. [23]) as ∂ηi (r, t) δF = −Li ∂t δηi (r, t)
(i = 1, 2, . . . , p),
(5)
where Li are the kinetic coefficients related to the grain boundary mobility, ∂ stands for partial differentiation, and δ denotes the first variation. During the final stage of sintering of ceramics, pores and grains co-evolve to reduce the free energy of the sintered material system. With different migration velocity ratios between pores and grains, pores can keep up with grain boundary motion or get trapped inside the grains. But from experimental observations [1,2], most residual pores distribute on grain boundaries and move with boundary migration because the free energy is lower for the pore occupying the grain boundaries to reduce the
Fig. 1. Schematic of pore–boundary interaction.
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
total grain boundary area. Mobile pores exert a pinning force on the grain boundary migration, thereby modifying the growth rate, as illustrated in Fig. 1. Under the assumptions of flat grain boundary randomly distributed with two-phase particles, Zener proposed that the pinning pressure exerted by a particle distribution on a unit area of grain boundary is proportional to the volume fraction of pores [24]. Yu and Suo [20] assumed that the pinning pressure on the grain boundary migration caused by a mobile pore is to be independent of the pore radius in two dimensions. Therefore, in this study, a constant pinning pressure by randomly distributed pores (i.e., interaction energy per unit volume pore on a planar grain boundary) σ * , is logically used. The elastic strain energy induced by local pinning forces and
273
applied loads on the specimen boundary is incorporated in the phase field model to simulate the stress-modified grain growth process. Rather than modeling the interaction of discrete pores with grain boundaries in [20,21], the whole pinning interaction exerted by randomly residual pores during the final-stage sintering is taken into account by the mean-field method neglecting localized change of each pore. Different from volume-conserved two-phase material systems in [15], grains grow with decreasing relative pore densities. As mentioned above, most pores prefer to reside on grain boundaries from the experimental observations. Hence, the volume distribution assumed to of pores is reasonably be proportional to (1 − a np=1 η2p ) where np=1 η2p approaches
Fig. 2. Topological evolution of isotropic grain growth under hydrostatic pressure.
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
274
0 inside grains, but 1 at grain boundaries, and a is a constant related to specific material system, which can be determined in SEM testing figures. Thus, the Zener pinning stress on grain boundaries exerted by pores, which measures the interaction between pores and grains per volume, can be represented as n σij0 (r) = σij∗ 1 − a η2p (r) , (6) p=1
where σij0 (r) are the components of the Zener pinning stress, σij∗ = σ ∗ δij , δij are the components of the second-order identity tensor, and n is the number of orientation variables. In our simulation, a = 1 is adopted because of the lack of experimental data. The equivalent stress-free strain ε0ij (r) associated with the Zener pinning stress σij0 (r) can be written as n ε0ij (r) = ε∗ij 1 − η2p (r) , (7) p=1
where the strain ε∗ij is related to σij∗ through σij∗ = Cijkl ε∗kl = Cijkl ε∗ δkl , with Cijkl being the elasticity tensor and ε* being the equivalent pinning strain. According to Khachaturyan’s elastic phase transformation theory [25], the elastic strain energy of a deformed material can be expressed as 1 el Eel = Cijkl εel ij (r)εkl (r) dV 2 V 1 = Cijkl ε¯ ij + δεij (r) − ε0ij (r) 2 V × ε¯ kl + δεkl (r) − ε0kl (r) dV, (8)
where (· · ·) represents the volume average of (· · ·), V the total −1 volume of the specimen, Sijkl (=Cijkl ) the components of the elasa tic compliance tensor, σij the components of the applied stress, η2q (r) = η2p (r) exp(−ig × r) d3 r the Fourier transform of g ∗ η2p (r), and η2p (r) is the complex conjugate of η2p (r) , g
−1 (n) = nj Cijkl nl = Ωik
p=1
−1 −
n
g
Cijkl =
2µν δij δkl + µ(δik δjl + δil δjk ), 1 − 2ν
Gik (g) =
1 2
(10)
g i gk g2 δik − , µ 2µ(1 − ν)
(11)
where µ is the shear modulus and ν is the Poisson’s ratio. For elastically cubic materials (e.g. [26]), it gives: Cijkl = C12 δij δkl + C44 (δik δjl + δil δjk ) + B
3
δip δjp δkp δlp ,
p=1
(12) Gik (g) =
C44 − ×
δik + Bgi2
g2
(C44
(C12 + C44 )gi gk + Bgi2 )(C44 g2 + Bgk2 )
g2
1 (1 + (C12 + C44 ) p f (gp ))
(13)
where C12 , C44 , B = C11 − C12 − 2C44 are elastic constants, and f (gp ) =
gp2 . C44 g2 +Bgp2
The rules in summation of dummy index
are not applied in Eq. (13).
q=1
η2p (r) 1 −
p=1
−
1 1 gj Cijkl gl = 2 G−1 (n), 2 g g ik
where nj = gj , g2 = gk gk , and Gik (g) is the Green function tensor in Fourier space [25]. For isotropic elastic materials, it is known that (e.g. [26]):
εel ij (r)
is the elastic strain, ε¯ ij the applied homogeneous where strain, δεkl (r) the heterogeneous strain, and ε0kl (r) is the stressfree strain. For an elastically inhomogeneous body under traction boundary conditions, Eq. (8) can be evaluated to give [25]: n n V Eel = Cijkl ε∗ij ε∗kl 1 − η2p (r) 1 − η2q (r) 2
g
Bpq (n) ≡ nj σji (p)Ωik (n)σkl (q)nl , with Ωik (n) being the unit Green function tensor that is inverse to the tensor
n
η2q (r)
q=1
n
∗
1 B (n)1 − η2p (r) 3 pq (2π) V p=1
× 1 −
n
g
η2q (r) dV −
q=1
− Vσija ε∗ij 1 −
g n p=1
V a Sijkl σija σkl 2
η2p (r),
(9)
Fig. 3. Temporal evolution of the mean grain size (with and without hydrostatic pressure).
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
It is known from the theory of interface migration [27] that the boundary migration speed v can be expressed as ¯ v = Mµ = MΩγgb k,
(14)
where M is the grain boundary mobility, Ω the atomic volume, and k¯ is the mean curvature of a grain boundary. In the current study, Eq. (14) is modified to include the effects of the stress σ * induced by the pore Zener pinning and external loading, which gives v = Mµ = MΩ(γgb k¯ + σ ∗ ).
(15)
275
From this equation, the boundary motion is determined by both the local mean curvature and the elastic stress field. As a result, the strain energy density Eel resulting from pore pinning forces and applied loads is incorporated in the total free energy (see Eq. (2)). Normally, the presence of pores not only results in a change of the total free energy, but also modifies the grain boundary mobility. But many numerical results (e.g. [28–30]) verified that the change of boundary mobility plays little role in the stochastic evolution of the microstructure compared to the effects of energy change. Therefore, the constant boundary mobility (the Li terms in Eq. (5)) is used here.
Fig. 4. Topological evolution of cubically anisotropic grain growth under hydrostatic pressure (B > 0).
276
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
The modified phase field model described above is implemented in our computational program for simulating the timedependent grain growth process of structural ceramics, which is discussed in the next section. 3. Simulation results To facilitate the computation, the parameters involved in the modified phase field model are made non-dimensional. The length scale and the time scale are set to be b and τ, respectively, so that xi∗ = xbi and t ∗ = τt give non-dimensional length and time. In the Landau free energy density polynomial given in
Eq. (4), Ai (i = 1, 2, 3) can be replaced by dimensionless paramAi eter A∗i = f , where f is defined near Eq. (3). Also, ki and Li introduced in Eqs. (2) and (5), respectively, are replaced by dimensionless parameters κi∗ = 2.0 and L∗i = 1.0. Following [13], the values of the dimensionless constants in our simulations are taken to be: A∗1 = 1.0, A∗2 = 1.0, A∗3 = 1.0. Consider a 2D specimen that can be represented by periodical arrays of a unit cell whose shape is a square with the side length of L. Periodic boundary conditions are applied to the cell. The cell is divided into N × N sub-cells, with the distance between L neighboring grids being d = N . The value of N (an integer) can be chosen in such a way that the value of d equals b, which is
Fig. 5. Topological evolution of cubically anisotropic grain growth under hydrostatic pressure (B < 0).
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
taken to be b = ηeq
κi f ,
277
where ηeq is the equilibrium value of
orientation variables. And the intrinsic time can be defined by 1 . τ = Lf Extending from Eq. (5), the kinetic equations including the elastic strain energy shown in Eq. (9) are given ∂f0 δEel ∂ηi (r, t) 2 = −L − κi ∇ ηi (r, t) + . (16) ∂t ∂ηi (r, t) δηi (r, t) In order to improve computation efficiency, the semi-implicit Fourier transform scheme [31] is applied to solve the above temporal kinetic equations, which yields η˜ n+1 (k, t) − η˜ ni (k, t) i t ∗ ∂f˜ 0 k12 k22 n ∗ ∗ 2 n+1 = −Li + κi k η˜ i (k, t) + ω 4 η˜ i (k, t) , ∂˜ηni (k, t) k (17) ˜ el , represent, respectively, the Fourier where η˜ i , f˜ 0 and E transform of ηi , f0and Eel , k = (k1 , k2 ) is a vector in the Fourier space, k =
k12 + k22 is the magnitude of k, and ω =
2µν (C11 +2C12 )2 ε∗2 B. For isotropic elastic materials, C12 = 1−2ν , (C12 +2C44 )2 f C44 = µ and C11 = C12 + 2C44 , which lead to a simpler expres-
sion for ω. For illustration purpose, 256 × 256 points are adopted with the time step t* = 0.1, the spatial step x* = 1 and the magnitude of ω, |ω| = 0.5. To visualize the microstructural features using orientation field variables, define [13]: φ(r) =
p
[ηi (r)]2 .
(18)
i=1
The value of φ(r) is illustrated by the gray level of figures. The higher values close to 1.0 within grains are marked as white, and the lower values at grain boundaries marked as black. First, the isotropic elastic material discussed in Section 2 is considered. The associated elastic constants and Green’s function are given in Eqs. (10) and (11), respectively. When the specimen is subjected to a hydrostatic pressure, the time-dependent evolution process of the microstructure is shown in Fig. 2. The initial microstructure is formed from the liquid phase by a small perturbation on all field variables, and crystallization occurs by minimizing the bulk free energy. The microstructure evolves in a way determined by the reduction of the excess chemical potential and strain energy. It is shown that the smaller grains get absorbed into the larger grains so that the number of grain sides and the grain sizes increase gradually. Accordingly, compared to the free grain growth without stress effects obtained in [13,14], the applied hydrostatic pressure and the intrinsic pinning force do not affect the geometric shape of grains, and hence the grains evolve in a self-similar manner into equi-axial shapes. Furthermore, the growth rates in the two aforementioned cases are compared to each other in Fig. 3. The pinning force by mobile pores at grain boundaries does not change diffusion paths of atoms and migration ways of grain boundaries, but reduces the
Fig. 6. Temporal evolution of mean grain size for B > 0 and B < 0.
migration rate of grain boundary. As a result, the average grain size decreases in the stress-induced growth condition because of constraint effects of pores on grain boundary, which is beneficial to control the abnormal grain growth and enhancement of sintered materials. The effects of anisotropy are now considered for elastically cubic materials (e.g., cubic crystals) discussed in Section 2. Cubic crystals are widely used, and perhaps the simplest, anisotropic materials. The elastic constants and the Green’s function for cubic crystals are given in Eqs. (12) and (13), respectively. When the anisotropic factor B > 0 (i.e., C12 + 2C44 < C11 ), the elastic strain energy density is smaller along the diagonal direction. When B < 0 (i.e., C12 + 2C44 < C11 ), the elastic strain energy density is smaller along the axial direction. Figs. 4 and 5 show the temporal topological evolution features with B > 0 and B < 0, respectively. The anisotropic behavior of the material destroys the self-similar growth mechanism and encourages the orientation growth along specific directions. When B > 0, grains elongate along 45◦ direction because of the low strain energy density along this direction where the elastic modulus C12 + 2C44 is smaller; conversely, when B < 0, grains evolve in columnar shape, thereby shortening along the diagonal direction constrained by the high strain energy density. It can therefore be concluded that the direction of the smaller elastic strain energy density is favorable for microstructure evolution. For the same hydrostatic pressure applied, it is observed from Figs. 3 and 6 that the growth rate of elastically cubic materials is lower than that of isotropic materials, which may be caused by the enhanced diffusion barrier along specific orientations or by the change of grain boundary migration paths. In Fig. 6, the average size of diagonal grains is seen to be bigger than that of columnar grains. The total free energy can be reduced not only by simply decreasing the overall grain boundary area but also by changing the orientations of grains. 4. Conclusions In the final stage of sintering of ceramics, residual pores migrate along the grain boundaries. In an attempt to understand and predict the interplay between the grain boundary migration and the pore drag coupled with an applied hydrostatic pressure
278
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
as well as their effects on grain growth kinetics and morphological evolution, a modified phase field model is proposed and utilized to simulate the microstructural evolution of two crystalline materials—isotropic and cubic. Our simulation results reveal that the motion of grain boundaries is retarded by the pore pinning forces. For isotropic elastic materials, grains grow in a self-similar manner, and for elastically cubic materials, grains grow in preferential directions. By the minimization principle of the total free energy, grains always grow along the direction of the smaller elastic strain energy. Furthermore, the growth rate slows down because of the Zener pinning effects and the applied hydrostatic pressure. The different grain growth rate may be due to two factors: one is anisotropic modification of strain energy to potential barriers for grain boundary migration; and the other is related to the diverse diffusivity and migration mechanisms. Acknowledgements The work reported here is funded by a grant from NSF (Grant # CMS-0324461, with Dr. Ken Chong as the program manager) and a grant from NSFC of China (Grant No. 50072025). These supports are gratefully acknowledged. References [1] M.N. Rahaman, Ceramic Processing and Sintering, Marcel Dekker, New York, 2003. [2] R.M. German, Sintering Theory and Practice, John Wiley, New York, 1996. [3] H.V. Atkinson, Acta Metall. 36 (1988) 469–491. [4] W.W. Mullins, J. Appl. Phys. 27 (1956) 900–904.
[5] G. Abbruzzese, I. Heckelmann, K. Lucke, Acta Metall. Mater. 40 (1992) 519–532. [6] G. Abbruzzese, K. Lucke, Mater. Sci. Forum 204 (1996) 55–70. [7] P.R. Rios, K. Lucke, Scripta Mater. 44 (2001) 2471–2475. [8] C.V. Thompson, H.J. Frost, F. Spaepen, Acta Metall. 35 (1987) 887–890. [9] M.P. Anderson, D.J. Srolovitz, G.S. Grest, P.S. Sahni, Acta Metall. 32 (1984) 783–791. [10] D.J. Srolovitz, M.P. Anderson, P.S. Grest, P.S. Sahni, Acta Metall. 32 (1984) 1429–1438. [11] V. Tikare, J.D. Cawley, Acta Metall. 46 (1998) 1333–1342. [12] J. Geiger, A. Roosz, P. Barkoczy, Acta Metall. 49 (2001) 623–629. [13] D. Fan, L.Q. Chen, Acta Mater. 45 (1997) 611–622. [14] D.N. Fan, L.Q. Chen, S.P. Chen, Mater. Sci. Eng. A 238 (1997) 78–84. [15] D. Fan, L.Q. Chen, Acta Mater. 45 (1997) 4145–4154. [16] R. Kobayashi, J.A. Warren, W.C. Carter, Physica D 119 (1998) 415–423. [17] R. Kobayashi, J.A. Warren, W.C. Carter, Physica D 140 (2000) 141–150. [18] J.A. Warren, R. Kobayashi, W.C. Carter, J. Cryst. Growth 211 (2000) 18–20. [19] A.E. Lobkovsky, J.A. Warren, J. Cryst. Growth 225 (2001) 282–288. [20] H.H. Yu, Z. Suo, J. Mech. Phys. Solids 47 (1999) 1131–1155. [21] H. Wang, Z.H. Li, Metall. Mater. Trans. A 34 (2003) 1493–1500. [22] J.W. Cahn, J.E. Hilliard, J. Chem. Phys. 28 (1958) 258–267. [23] S.M. Allen, J.W. Cahn, Acta Metall. 27 (1979) 1085–1095. [24] N. Moelans, B. Blanpain, P. Wollants, Acta Mater. 53 (2005) 1171–1181. [25] A.G. Khachaturyan, Theory of Structural Transformation in Solids, John Wiley, New York, 1983. [26] T. Mura, Micromechanics of Defects in Solids, Kluwer Academic Publishers, Dordrecht, 1987. [27] Z. Suo, Adv. Appl. Mech. 33 (1997) 193–294. [28] A. Kazaryan, Y. Wang, S.A. Dregia, et al., Acta Metall. 50 (2002) 2491–2502. [29] A. Kazaryan, Y. Wang, S.A. Dregia, et al., Phys. Rev. B 61 (2000) 14275–14278. [30] M. Upamanyu, G.N. Hassold, A. Kazaryan, et al., Interf. Sci. 10 (2002) 201–216. [31] L.Q. Chen, J. Shen, Comput. Phys. Commun. 108 (1998) 147–158.
Anisotropic grain growth with pore drag under applied loads X.N. Jing a,b , J.H. Zhao a , G. Subhash b,∗ , X.-L. Gao c a
CAS Key Laboratory of Mechanical Behavior and Design of Materials, University of Science and Technology of China, Hefei, Anhui 230027, PR China b Department of Mechanical Engineering and Engineering Mechanics, Michigan Technological University, 1400 Townsend Drive, Houghton, MI 49931-1295, USA c Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843-3123, USA Accepted 31 August 2005
Abstract In the final stage of sintering of ceramics, residual pores co-evolve with grain boundaries because of incomplete densification. Their interactions coupled with external loads are critical to the microstructural evolution of structural ceramics. A modified two-dimensional (2D) diffuse-interface phase field model, which differs from the boundary-tracking methods, is utilized to investigate the effects of stochastically distributed pore drag on grain growth kinetics and morphological evolution process of ceramics under applied loads. Contributions from both the boundary energy and elastic strain energy caused by pore drag forces and applied loading are incorporated in the modified phase field model to describe the isotropic or cubically anisotropic behaviors of polycrystalline materials. The temporal evolution of the spatially dependent grain orientation variables is determined by numerically solving non-linear Ginzburg–Landau equations using a semi-implicit Fourier-spectral method. Numerical results show that the anisotropic strain energy dominates the non-self-similar growth manner, leads to ordered grain morphologies and changes the growth rate. © 2005 Elsevier B.V. All rights reserved. Keywords: Phase field model; Grain growth; Pore drag; Morphological evolution
1. Introduction Mechanical properties of polycrystalline materials (e.g., ceramics) depend strongly on the size, distribution and morphology of grains. Thus, how to control grain growth during processing is a critical issue for the development of polycrystalline materials [1–3]. Grain growth can be viewed as a process of grain boundary migration to decrease the total grain boundary areas and the total free energy of the material system, both of which are driven by mean curvatures of grain boundaries [1,2]. In a single-phase material, the only process which occurs during grain growth is local atomic re-arrangement. However, in a multi-phase material, the long-range diffusion dominates the migration of grain boundary. Hence, the kinetics of grain growth is strongly affected by the presence or absence of solute or impurity migration and segregation at grain boundaries. The pinning effect of pores is of great practical importance in the sintering of high-quality ceramics for structural applications, where den-
∗
Corresponding author. Tel.: +1 906 487 3161; fax: +1 906 487 2822. E-mail address: [email protected] (G. Subhash).
0921-5093/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2005.08.220
sification and small grain sizes are often required to obtain good strength and toughness. In the last two decades, many models have been proposed to predict the time dependence of average grain size and size distribution in polycrystalline materials [3–7]. Due to the complexity of topologies and coupled multiple interactions, the analytical modeling of microstructure evolution is often very difficult. As a result, computer simulation plays a key role in exploring the details of grain growth and validating the analytical models. The sharp-interface approaches, such as Vortex model [8], Potts models [9–11] and Cellular Automata method [12], have been used in simulations of grain growth. Recently, simulating grain growth using a continuum-based diffuse-interface phase field model has been developed by several researchers [13–19]. Chen and co-workers simulated grain growth in single-phase materials [13,14] and studied the microstructure evolution in volume-conserved two-phase material systems with finite interface thickness [15] by considering grains of different crystallographic orientations represented by a set of non-conserved order parameter fields. In their model, the effects of solute or precipitate drag on grain growth rate and size distribution are taken into account without ad hoc assumptions
272
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
about grain shapes. An independent phase field model was developed in [16–19], which adopted two order parameters to describe a grain structure; one denoting the crystalline order and the other describing predominantly the local orientation of the crystal. In this model, the grain rotation can be simulated by the relaxation of the local orientation order parameter. During the final-stage sintering of ceramics, most residual pores are isolated, residing on grain boundaries, three-grain junctions or four-grain junctions due to incomplete sintering process [1,2]. The volume fraction of pores decreases as atoms elsewhere diffuse along the grain boundaries into pores, and grains grow as atoms migrate between grains. Yu and Suo [20] developed an axisymmetric model to study the steady and transient motion of pore–grain boundary separation using the finite element method. Wang and Li [21] further considered the shrinkage of grain boundary pores under pressure. However, it appears that the effects of pore shrinking and pinning on the grain growth process during the final-stage sintering have not been fully discussed. To progress from this situation, the effects of the Zener pinning of randomly non-conservative pores moving with grain boundaries on grain growth rate and anisotropic growth mechanism are studied from the energy change point of view using a diffuse-interface phase field model reported in [13,14]. 2. Phase field model description In the diffuse-interface field model developed by Fan and Chen [13,14], the microstructure of a polycrystalline material specimen is described by using a set of p continuous nonconserved long-range order parameters {ηi (r,t)} (i = 1, 2, . . ., p) defined at a given time t at each position r within the specimen, i.e. η1 (r, t), η2 (r, t), . . . , ηp (r, t),
(1)
where p is the number of possible orientations of a grain in space. In real materials, the number of orientations is infinite (i.e., p → ∞), but a finite value of p (p > 32) is often sufficient for reasonably accurate simulations of grain growth. ηi (r,t) (i = 1, 2, . . ., p) serve as orientation field variables to distinguish different crystallographic orientations of grains and are taken to be continuous in space. When the total free energy is minimized, only one field variable within each grain takes on a value of 1 or −1, and all other field variables are zero-valued. Across a grain boundary between two adjacent grains, ηi changes continuously between 1 and 0 or −1 and 0. According to the diffuse-interface theories [22], the total free energy F of the material specimen can be described as the function of a set of field variables: p κi 2 F= f0 (η1 (r), η2 (r), . . . , ηp (r)) + dV, [∇ηi (r)] 2 V i=1 (2) where f0 is the local free energy density as a function of ηi (r), ∇ηi (r) is the gradient of ηi (r), κi are the gradient energy coefficients which are positive constants, and V is the volume of the specimen. The gradient term is related to the boundary energy.
For a planar grain boundary between a grain of orientation i and another grain of orientation j, the grain boundary energy γ gb is given by [13]: ∞ κi dηi 2 κj dηj 2 γgb = f (ηi , ηj ) + dx, (3) + 2 dx 2 dx −∞ where f(ηi ,ηj ) = f0 (0, . . ., ηi , ηj , . . ., 0) − fmin , with fmin being the minima of f0 . In order to model the temporal and spatial evolution of the material microstructure, an approximate Landau type free energy density is proposed as [13,14]: p 2 p p p A 1 2 A2 2 A2 2 2 f0 = − ηi + ηi + A3 − ηi η j , 2 4 2 i=1 i=1 i=1 j=i (4) where A1 , A2 , A3 are phenomenological constants whose values are determinable from grain boundary energy and grain boundary width. This free energy density has 2p potential wells in the p-field space, which represent the equilibrium free energies of crystalline grains in 2p different orientations. The kinetic equations that describe microstructure evolution with time and space are given by the time-dependent Ginzburg–Landau equations (e.g. [23]) as ∂ηi (r, t) δF = −Li ∂t δηi (r, t)
(i = 1, 2, . . . , p),
(5)
where Li are the kinetic coefficients related to the grain boundary mobility, ∂ stands for partial differentiation, and δ denotes the first variation. During the final stage of sintering of ceramics, pores and grains co-evolve to reduce the free energy of the sintered material system. With different migration velocity ratios between pores and grains, pores can keep up with grain boundary motion or get trapped inside the grains. But from experimental observations [1,2], most residual pores distribute on grain boundaries and move with boundary migration because the free energy is lower for the pore occupying the grain boundaries to reduce the
Fig. 1. Schematic of pore–boundary interaction.
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
total grain boundary area. Mobile pores exert a pinning force on the grain boundary migration, thereby modifying the growth rate, as illustrated in Fig. 1. Under the assumptions of flat grain boundary randomly distributed with two-phase particles, Zener proposed that the pinning pressure exerted by a particle distribution on a unit area of grain boundary is proportional to the volume fraction of pores [24]. Yu and Suo [20] assumed that the pinning pressure on the grain boundary migration caused by a mobile pore is to be independent of the pore radius in two dimensions. Therefore, in this study, a constant pinning pressure by randomly distributed pores (i.e., interaction energy per unit volume pore on a planar grain boundary) σ * , is logically used. The elastic strain energy induced by local pinning forces and
273
applied loads on the specimen boundary is incorporated in the phase field model to simulate the stress-modified grain growth process. Rather than modeling the interaction of discrete pores with grain boundaries in [20,21], the whole pinning interaction exerted by randomly residual pores during the final-stage sintering is taken into account by the mean-field method neglecting localized change of each pore. Different from volume-conserved two-phase material systems in [15], grains grow with decreasing relative pore densities. As mentioned above, most pores prefer to reside on grain boundaries from the experimental observations. Hence, the volume distribution assumed to of pores is reasonably be proportional to (1 − a np=1 η2p ) where np=1 η2p approaches
Fig. 2. Topological evolution of isotropic grain growth under hydrostatic pressure.
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
274
0 inside grains, but 1 at grain boundaries, and a is a constant related to specific material system, which can be determined in SEM testing figures. Thus, the Zener pinning stress on grain boundaries exerted by pores, which measures the interaction between pores and grains per volume, can be represented as n σij0 (r) = σij∗ 1 − a η2p (r) , (6) p=1
where σij0 (r) are the components of the Zener pinning stress, σij∗ = σ ∗ δij , δij are the components of the second-order identity tensor, and n is the number of orientation variables. In our simulation, a = 1 is adopted because of the lack of experimental data. The equivalent stress-free strain ε0ij (r) associated with the Zener pinning stress σij0 (r) can be written as n ε0ij (r) = ε∗ij 1 − η2p (r) , (7) p=1
where the strain ε∗ij is related to σij∗ through σij∗ = Cijkl ε∗kl = Cijkl ε∗ δkl , with Cijkl being the elasticity tensor and ε* being the equivalent pinning strain. According to Khachaturyan’s elastic phase transformation theory [25], the elastic strain energy of a deformed material can be expressed as 1 el Eel = Cijkl εel ij (r)εkl (r) dV 2 V 1 = Cijkl ε¯ ij + δεij (r) − ε0ij (r) 2 V × ε¯ kl + δεkl (r) − ε0kl (r) dV, (8)
where (· · ·) represents the volume average of (· · ·), V the total −1 volume of the specimen, Sijkl (=Cijkl ) the components of the elasa tic compliance tensor, σij the components of the applied stress, η2q (r) = η2p (r) exp(−ig × r) d3 r the Fourier transform of g ∗ η2p (r), and η2p (r) is the complex conjugate of η2p (r) , g
−1 (n) = nj Cijkl nl = Ωik
p=1
−1 −
n
g
Cijkl =
2µν δij δkl + µ(δik δjl + δil δjk ), 1 − 2ν
Gik (g) =
1 2
(10)
g i gk g2 δik − , µ 2µ(1 − ν)
(11)
where µ is the shear modulus and ν is the Poisson’s ratio. For elastically cubic materials (e.g. [26]), it gives: Cijkl = C12 δij δkl + C44 (δik δjl + δil δjk ) + B
3
δip δjp δkp δlp ,
p=1
(12) Gik (g) =
C44 − ×
δik + Bgi2
g2
(C44
(C12 + C44 )gi gk + Bgi2 )(C44 g2 + Bgk2 )
g2
1 (1 + (C12 + C44 ) p f (gp ))
(13)
where C12 , C44 , B = C11 − C12 − 2C44 are elastic constants, and f (gp ) =
gp2 . C44 g2 +Bgp2
The rules in summation of dummy index
are not applied in Eq. (13).
q=1
η2p (r) 1 −
p=1
−
1 1 gj Cijkl gl = 2 G−1 (n), 2 g g ik
where nj = gj , g2 = gk gk , and Gik (g) is the Green function tensor in Fourier space [25]. For isotropic elastic materials, it is known that (e.g. [26]):
εel ij (r)
is the elastic strain, ε¯ ij the applied homogeneous where strain, δεkl (r) the heterogeneous strain, and ε0kl (r) is the stressfree strain. For an elastically inhomogeneous body under traction boundary conditions, Eq. (8) can be evaluated to give [25]: n n V Eel = Cijkl ε∗ij ε∗kl 1 − η2p (r) 1 − η2q (r) 2
g
Bpq (n) ≡ nj σji (p)Ωik (n)σkl (q)nl , with Ωik (n) being the unit Green function tensor that is inverse to the tensor
n
η2q (r)
q=1
n
∗
1 B (n)1 − η2p (r) 3 pq (2π) V p=1
× 1 −
n
g
η2q (r) dV −
q=1
− Vσija ε∗ij 1 −
g n p=1
V a Sijkl σija σkl 2
η2p (r),
(9)
Fig. 3. Temporal evolution of the mean grain size (with and without hydrostatic pressure).
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
It is known from the theory of interface migration [27] that the boundary migration speed v can be expressed as ¯ v = Mµ = MΩγgb k,
(14)
where M is the grain boundary mobility, Ω the atomic volume, and k¯ is the mean curvature of a grain boundary. In the current study, Eq. (14) is modified to include the effects of the stress σ * induced by the pore Zener pinning and external loading, which gives v = Mµ = MΩ(γgb k¯ + σ ∗ ).
(15)
275
From this equation, the boundary motion is determined by both the local mean curvature and the elastic stress field. As a result, the strain energy density Eel resulting from pore pinning forces and applied loads is incorporated in the total free energy (see Eq. (2)). Normally, the presence of pores not only results in a change of the total free energy, but also modifies the grain boundary mobility. But many numerical results (e.g. [28–30]) verified that the change of boundary mobility plays little role in the stochastic evolution of the microstructure compared to the effects of energy change. Therefore, the constant boundary mobility (the Li terms in Eq. (5)) is used here.
Fig. 4. Topological evolution of cubically anisotropic grain growth under hydrostatic pressure (B > 0).
276
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
The modified phase field model described above is implemented in our computational program for simulating the timedependent grain growth process of structural ceramics, which is discussed in the next section. 3. Simulation results To facilitate the computation, the parameters involved in the modified phase field model are made non-dimensional. The length scale and the time scale are set to be b and τ, respectively, so that xi∗ = xbi and t ∗ = τt give non-dimensional length and time. In the Landau free energy density polynomial given in
Eq. (4), Ai (i = 1, 2, 3) can be replaced by dimensionless paramAi eter A∗i = f , where f is defined near Eq. (3). Also, ki and Li introduced in Eqs. (2) and (5), respectively, are replaced by dimensionless parameters κi∗ = 2.0 and L∗i = 1.0. Following [13], the values of the dimensionless constants in our simulations are taken to be: A∗1 = 1.0, A∗2 = 1.0, A∗3 = 1.0. Consider a 2D specimen that can be represented by periodical arrays of a unit cell whose shape is a square with the side length of L. Periodic boundary conditions are applied to the cell. The cell is divided into N × N sub-cells, with the distance between L neighboring grids being d = N . The value of N (an integer) can be chosen in such a way that the value of d equals b, which is
Fig. 5. Topological evolution of cubically anisotropic grain growth under hydrostatic pressure (B < 0).
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
taken to be b = ηeq
κi f ,
277
where ηeq is the equilibrium value of
orientation variables. And the intrinsic time can be defined by 1 . τ = Lf Extending from Eq. (5), the kinetic equations including the elastic strain energy shown in Eq. (9) are given ∂f0 δEel ∂ηi (r, t) 2 = −L − κi ∇ ηi (r, t) + . (16) ∂t ∂ηi (r, t) δηi (r, t) In order to improve computation efficiency, the semi-implicit Fourier transform scheme [31] is applied to solve the above temporal kinetic equations, which yields η˜ n+1 (k, t) − η˜ ni (k, t) i t ∗ ∂f˜ 0 k12 k22 n ∗ ∗ 2 n+1 = −Li + κi k η˜ i (k, t) + ω 4 η˜ i (k, t) , ∂˜ηni (k, t) k (17) ˜ el , represent, respectively, the Fourier where η˜ i , f˜ 0 and E transform of ηi , f0and Eel , k = (k1 , k2 ) is a vector in the Fourier space, k =
k12 + k22 is the magnitude of k, and ω =
2µν (C11 +2C12 )2 ε∗2 B. For isotropic elastic materials, C12 = 1−2ν , (C12 +2C44 )2 f C44 = µ and C11 = C12 + 2C44 , which lead to a simpler expres-
sion for ω. For illustration purpose, 256 × 256 points are adopted with the time step t* = 0.1, the spatial step x* = 1 and the magnitude of ω, |ω| = 0.5. To visualize the microstructural features using orientation field variables, define [13]: φ(r) =
p
[ηi (r)]2 .
(18)
i=1
The value of φ(r) is illustrated by the gray level of figures. The higher values close to 1.0 within grains are marked as white, and the lower values at grain boundaries marked as black. First, the isotropic elastic material discussed in Section 2 is considered. The associated elastic constants and Green’s function are given in Eqs. (10) and (11), respectively. When the specimen is subjected to a hydrostatic pressure, the time-dependent evolution process of the microstructure is shown in Fig. 2. The initial microstructure is formed from the liquid phase by a small perturbation on all field variables, and crystallization occurs by minimizing the bulk free energy. The microstructure evolves in a way determined by the reduction of the excess chemical potential and strain energy. It is shown that the smaller grains get absorbed into the larger grains so that the number of grain sides and the grain sizes increase gradually. Accordingly, compared to the free grain growth without stress effects obtained in [13,14], the applied hydrostatic pressure and the intrinsic pinning force do not affect the geometric shape of grains, and hence the grains evolve in a self-similar manner into equi-axial shapes. Furthermore, the growth rates in the two aforementioned cases are compared to each other in Fig. 3. The pinning force by mobile pores at grain boundaries does not change diffusion paths of atoms and migration ways of grain boundaries, but reduces the
Fig. 6. Temporal evolution of mean grain size for B > 0 and B < 0.
migration rate of grain boundary. As a result, the average grain size decreases in the stress-induced growth condition because of constraint effects of pores on grain boundary, which is beneficial to control the abnormal grain growth and enhancement of sintered materials. The effects of anisotropy are now considered for elastically cubic materials (e.g., cubic crystals) discussed in Section 2. Cubic crystals are widely used, and perhaps the simplest, anisotropic materials. The elastic constants and the Green’s function for cubic crystals are given in Eqs. (12) and (13), respectively. When the anisotropic factor B > 0 (i.e., C12 + 2C44 < C11 ), the elastic strain energy density is smaller along the diagonal direction. When B < 0 (i.e., C12 + 2C44 < C11 ), the elastic strain energy density is smaller along the axial direction. Figs. 4 and 5 show the temporal topological evolution features with B > 0 and B < 0, respectively. The anisotropic behavior of the material destroys the self-similar growth mechanism and encourages the orientation growth along specific directions. When B > 0, grains elongate along 45◦ direction because of the low strain energy density along this direction where the elastic modulus C12 + 2C44 is smaller; conversely, when B < 0, grains evolve in columnar shape, thereby shortening along the diagonal direction constrained by the high strain energy density. It can therefore be concluded that the direction of the smaller elastic strain energy density is favorable for microstructure evolution. For the same hydrostatic pressure applied, it is observed from Figs. 3 and 6 that the growth rate of elastically cubic materials is lower than that of isotropic materials, which may be caused by the enhanced diffusion barrier along specific orientations or by the change of grain boundary migration paths. In Fig. 6, the average size of diagonal grains is seen to be bigger than that of columnar grains. The total free energy can be reduced not only by simply decreasing the overall grain boundary area but also by changing the orientations of grains. 4. Conclusions In the final stage of sintering of ceramics, residual pores migrate along the grain boundaries. In an attempt to understand and predict the interplay between the grain boundary migration and the pore drag coupled with an applied hydrostatic pressure
278
X.N. Jing et al. / Materials Science and Engineering A 412 (2005) 271–278
as well as their effects on grain growth kinetics and morphological evolution, a modified phase field model is proposed and utilized to simulate the microstructural evolution of two crystalline materials—isotropic and cubic. Our simulation results reveal that the motion of grain boundaries is retarded by the pore pinning forces. For isotropic elastic materials, grains grow in a self-similar manner, and for elastically cubic materials, grains grow in preferential directions. By the minimization principle of the total free energy, grains always grow along the direction of the smaller elastic strain energy. Furthermore, the growth rate slows down because of the Zener pinning effects and the applied hydrostatic pressure. The different grain growth rate may be due to two factors: one is anisotropic modification of strain energy to potential barriers for grain boundary migration; and the other is related to the diverse diffusivity and migration mechanisms. Acknowledgements The work reported here is funded by a grant from NSF (Grant # CMS-0324461, with Dr. Ken Chong as the program manager) and a grant from NSFC of China (Grant No. 50072025). These supports are gratefully acknowledged. References [1] M.N. Rahaman, Ceramic Processing and Sintering, Marcel Dekker, New York, 2003. [2] R.M. German, Sintering Theory and Practice, John Wiley, New York, 1996. [3] H.V. Atkinson, Acta Metall. 36 (1988) 469–491. [4] W.W. Mullins, J. Appl. Phys. 27 (1956) 900–904.
[5] G. Abbruzzese, I. Heckelmann, K. Lucke, Acta Metall. Mater. 40 (1992) 519–532. [6] G. Abbruzzese, K. Lucke, Mater. Sci. Forum 204 (1996) 55–70. [7] P.R. Rios, K. Lucke, Scripta Mater. 44 (2001) 2471–2475. [8] C.V. Thompson, H.J. Frost, F. Spaepen, Acta Metall. 35 (1987) 887–890. [9] M.P. Anderson, D.J. Srolovitz, G.S. Grest, P.S. Sahni, Acta Metall. 32 (1984) 783–791. [10] D.J. Srolovitz, M.P. Anderson, P.S. Grest, P.S. Sahni, Acta Metall. 32 (1984) 1429–1438. [11] V. Tikare, J.D. Cawley, Acta Metall. 46 (1998) 1333–1342. [12] J. Geiger, A. Roosz, P. Barkoczy, Acta Metall. 49 (2001) 623–629. [13] D. Fan, L.Q. Chen, Acta Mater. 45 (1997) 611–622. [14] D.N. Fan, L.Q. Chen, S.P. Chen, Mater. Sci. Eng. A 238 (1997) 78–84. [15] D. Fan, L.Q. Chen, Acta Mater. 45 (1997) 4145–4154. [16] R. Kobayashi, J.A. Warren, W.C. Carter, Physica D 119 (1998) 415–423. [17] R. Kobayashi, J.A. Warren, W.C. Carter, Physica D 140 (2000) 141–150. [18] J.A. Warren, R. Kobayashi, W.C. Carter, J. Cryst. Growth 211 (2000) 18–20. [19] A.E. Lobkovsky, J.A. Warren, J. Cryst. Growth 225 (2001) 282–288. [20] H.H. Yu, Z. Suo, J. Mech. Phys. Solids 47 (1999) 1131–1155. [21] H. Wang, Z.H. Li, Metall. Mater. Trans. A 34 (2003) 1493–1500. [22] J.W. Cahn, J.E. Hilliard, J. Chem. Phys. 28 (1958) 258–267. [23] S.M. Allen, J.W. Cahn, Acta Metall. 27 (1979) 1085–1095. [24] N. Moelans, B. Blanpain, P. Wollants, Acta Mater. 53 (2005) 1171–1181. [25] A.G. Khachaturyan, Theory of Structural Transformation in Solids, John Wiley, New York, 1983. [26] T. Mura, Micromechanics of Defects in Solids, Kluwer Academic Publishers, Dordrecht, 1987. [27] Z. Suo, Adv. Appl. Mech. 33 (1997) 193–294. [28] A. Kazaryan, Y. Wang, S.A. Dregia, et al., Acta Metall. 50 (2002) 2491–2502. [29] A. Kazaryan, Y. Wang, S.A. Dregia, et al., Phys. Rev. B 61 (2000) 14275–14278. [30] M. Upamanyu, G.N. Hassold, A. Kazaryan, et al., Interf. Sci. 10 (2002) 201–216. [31] L.Q. Chen, J. Shen, Comput. Phys. Commun. 108 (1998) 147–158.