- Email: [email protected]
Numerical simulation of diffusion induced phase separation and coarsening in binary alloys
Computational Materials Science 50 (2011) 1359–1364
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Numerical simulation of diffusion induced phase separation and coarsening in binary alloys Denis Anders *, Kerstin Weinberg Lehrstuhl für Festkörpermechanik, Institut für Mechanik und Regelungstechnik, FB11, Universität Siegen, Paul-Bonatz-Straße 9-11, 57068 Siegen, Germany
a r t i c l e Article history: Received 1 October Received in revised Accepted 18 March Available online 18
i n f o 2009 form 11 March 2010 2010 April 2010
Keywords: Binary alloys Solder Spinodal decomposition Phase field Damage
a b s t r a c t Solder materials are solid mixtures and show a variety of micro-structural changes, for example, separation of phases and phase coarsening (Ostwald ripening). The micro-structural changes influence strength and life expectation of these materials. They may be of significant effect, in particular, in very small components such as solder joints in microelectronic packages. In order to analyze the micro-morphological evolution with a diffusion theory of heterogeneous solid mixtures we apply an extended Cahn–Hilliard phase field model. The underlying equations involve a bipotential operator and, accordingly, fourth-order spatial derivatives. Thus, the variational formulation of the Cahn–Hilliard problem requires approximation functions which are piecewise smooth and globally C 1 -continuous; standard finite element techniques cannot be applied here. In our contribution we fulfil the continuity requirement by means of rational B-spline finite element basis functions. To illustrate the versatility of this approach numerical simulations of phase decomposition and coarsening controlled by diffusion and by mechanical loading are discussed and compared with experimental results. Studies of error and convergence are provided. Ó 2010 Elsevier B.V. All rights reserved.
1. The problem of phase decomposition and coarsening Solder joints in microelectronic circuit units are required to provide mechanical as well as electrical connections between different components. The life expectation of these joints determines the reliability of the whole structure. The overall material properties of the solder are affected by the decomposition of the alloy into phases of different composition. Experimental observations show that separation and coarsening of phases generate domains, which are susceptible to mechanical failure. Voids, crack initiation and crack propagation are prevalent observed along phase boundaries. The distribution of components in an alloy can be described by their mass concentration, cI ; cII , etc. In a binary alloy it simply holds cII ¼ 1 cI ¼: c. In order to investigate the evolution of the concentration field cðx; tÞ in a representative domain X within time t we employ a phase field model. Originally introduced by Cahn and Hilliard for diffusion driven phase separation [4–7], this approach has been extended by several contributors to account for additional fields, e.g., thermo-mechanical stresses [8–10,12]. In our context the problem reads: Find c : X ð0; tÞ ! R such that
* Corresponding author. E-mail addresses: [email protected] (D. Anders), [email protected] (K. Weinberg). 0927-0256/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2010.03.030
@c ¼ r j ¼ r ðMrlÞ @t @ ðF th þ F el Þ kDc in X ð0; tÞ ¼ r Mr @c
ð1Þ
holds. Here j is the flux of concentrations driven by the negative product of mobility M and the gradient of the chemical potential l. The initial concentration is given,
cðx; 0Þ ¼ c0 ðxÞ in X;
ð2Þ
and, in order to account for the conservation of mass, homogeneous Neumann boundary conditions are prescribed.
@FðcÞ kDc n ¼ 0 in @ X ð0; tÞ @c rc n ¼ 0 in @ X ð0; tÞ j n ¼ Mrl n ¼ Mr
ð3Þ ð4Þ
Here and above the free Helmholtz energy density F is composed of thermal and elastic components FðcÞ ¼ F th þ F el . The mobility tensor 2 M is assumed to be isotropic and constant. The parameter k ¼ jl is related to surface energy density j and length l of the transition regions between the domains of each phase. Vector n denotes the unit outward normal to @ X. The evolution of the concentration field can roughly be divided in two steps. Starting from an almost homogeneous solution at first phases separate. The driving quantity is here the minimization of the free Helmholtz energy density. In the strain free case
1360
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
ðF el ¼ 0Þ it can be obtained from the experimentally measured Gibbs free energy of the specific solution.
FðcÞ ¼ F th ðcÞ ¼
h ½c lnðcÞ þ ð1 cÞ lnð1 cÞ þ hs cð1 cÞ 2
ð5Þ
solution ch and the discrete test functions uh can be depicted as a linear combination of basis functions N i in Vh ,
uh ¼
n X
ai N i
and ch ¼
n X
i¼1
Here h and hs denote specific material parameters with h ¼ NkT; hs ¼ NkT s where N is Avogadro’s number, k Boltzmann’s constant, T the absolute temperature and T s the critical temperature (solidification temperature). For T > T s Eq. (5) is a convex function. A separation into phases is energetically disadvantageous in comparison to the homogeneous mixture. Thus, the mixture remains stable to all fluctuations. For T < T s Eq. (5) forms a double-well potential, i.e., it has two relative minima and a concave region (spinodal region) in between. The homogeneous mixture is here unstable and will decompose into two phases with concentration ca (a-phase) and cb (b-phase). It holds
@F Fðcb Þ Fðca Þ @F ðca ; TÞ ¼ ðcb ; TÞ: ¼ @c cb ca @c
ð6Þ
After separation the phases start to grow. This coarsening process is driven by the minimization of the interfacial energy density in such manner that the size of the emerged phase regions increases at the expense of their number. The interfacial energy contribution is described by the bipotential operator in Eq. (1). The diffusion based Cahn–Hilliard model is extended by an elastic energy density [9,10]
F el ðc; uÞ ¼
1 ððuÞ th ðcÞÞ : CðcÞ : ððuÞ th ðcÞÞ; 2
ð7Þ
where u ¼ uðx; tÞ denotes the displacement field. Its symmetric gradient is the linearized strain tensor ðuÞ ¼ 12 ðru þ ruT Þ. The elastic properties enter Eq. (7) via elasticity tensor CðcÞ and thermal strain tensor th ðcÞ ¼ ath ðcÞDT; both depending on concentration field cðx; tÞ.
bj N j ;
j¼1
where n is the dimension of Vh . Employing Galerkin orthogonality the discrete problem reads: Find ch 2 Vh such that
Z
Z Z 2 @ch @ F th ðch Þ rch rNi dX Ni dX þ Mk Dch DN i dX þ M @c2 X @t X X Z @F el ðch Þ þM r rNi dX ¼ 0 for i ¼ 1; . . . ; n @c X
ð9Þ
Please note that Eq. (9) is non-linear of ch . The required continuity ch 2 Vh H2 ðXÞ can only be fulfilled by basis functions N j which are piecewise smooth and globally at least C 1 -continuous. Standard finite element basis functions provide C 0 continuity only and cannot be used in this context. Extending the finite element space by Hermitian polynomials is elaborate and restricted to rectangular domains [11]. A mixed finite element approach, suggested, e.g., in [3], is also numerically expensive because it introduces additional unknowns to the primal degrees of freedom. Moreover, the formulation of boundary conditions is not completely clear here. In virtue of these difficulties we follow an approach of [1,2] and apply rational B-splines as finite element basis functions. These basis functions offer high-order accuracy, robustness, and, most importantly, C 1 -continuity. Univariate basic splines (B-splines) are defined recursively for a given order p ¼ 1; 2; 3; . . . and a knot vector n ¼ fni g; i ¼ 1; 2; 3; . . . by the following formula:
Ni;p ðxÞ :¼
niþpþ1 x x ni Ni;p1 ðxÞ þ Niþ1;p1 ðxÞ: niþp ni niþpþ1 niþ1
ð10Þ
Stop function for this recursion is the jump function:
1 if ni 6 x < niþ1
2. Discretization of the problem and numerical solution
Ni;0 ðxÞ :¼
For numerical approximation we derive the primal variational form of the Cahn–Hilliard model. Multiplying Eq. (1) with a test function u, applying Green’s formula (integration by parts) and making use of the homogeneous boundary conditions (3) and (4) give, together with Eq. (2), the problem: Find c 2 H2 ðXÞ such that
For a uniform knot vector n ¼ f0; 1; 2; 3; . . .g and order p ¼ 1 the Bspline functions coincide with piecewise linear finite element basis functions which are C 0 -continuous, see Fig. 1. For higher order, p P 2, shape and properties of B-spline functions differ significantly from standard higher order finite element basis functions. The quadratic B-Spline function of Fig. 1 is globally smooth and C 1 -continuous. Therefore, it will be employed here as basis function in our spatial discretization. To analyze two- and three-dimensional domains multivariate B-spline functions can be derived by tensor products of univariate functions. In two dimensions it holds, see Fig. 2:
Z X
Z Z @c @ðF th þ F el Þ u dX ¼ M r ru dX k DcDu dX @t @c X X
8 u 2 H2 ðXÞ:
ð8Þ
2
Here H ðXÞ denotes the Sobolev space of square integrable functions with square integrable derivatives of first and second order.
0
otherwise:
:
Nij;p ðx; yÞ :¼ Ni;p ðxÞ Nj;p ðyÞ
2.1. Spatial discretization For spatial discretization of Eq. (8) we employ the Galerkin method. To this end the space of test/ansatz functions H2 ðXÞ is approximated by a finite dimensional subspace Vh . The discrete
ð11Þ
ð12Þ
The support of each B-spline Ni;p ðxÞ is compact and contained in the interval ½ni ; niþpþ1 . Hence, the numerical integration required to solve Eq. (9) can be performed with standard algorithms. Additionally, B-spline functions establish a partition of unity,
Fig. 1. Basic splines N i;p of order p = 0, 1, 2 for uniform knot vector n ¼ f0; 1; 2; 3; . . .g.
1361
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
Fig. 2. Two-dimensional B-spline of order p = 2.
Fig. 3. Evolution of concentration in the test example.
n X
Ni;p ðxÞ ¼ 1 8x 2 ½n1 ; nnþpþ1 :
ð13Þ
with fgnþ1 ¼ 12 ðfgn þ fgnþ1 Þ denoting a midpoint evaluation. 2
i¼1
Each spline is a non-negative function, i.e., N i;p ðxÞ P 0, for all x. One of the most important features in our context is ability of B-spline functions to map jumps. The phase boundaries in the Cahn–Hilliard model are typically sharp layers which require a high spatial resolution. Hermitian and Lagrangian interpolation polynomials (used as standard finite element basis) tend to oscillate while representing jumps in a solution function. 2.2. Temporal discretization The temporal discretization employed here is straightforward. The first order time derivative is approximated by finite differences
@c cnþ1 cn ¼ @t Dt
ð14Þ
with (uniform) time step Dt ¼ t nþ1 t n . For time integration we employ an implicit Crank-Nicholson scheme which delivers a second order accuracy [11]. Consequently the fully discrete problem reads:
Z
Z cnþ1 cn Ni dX þ Mk Dchnþ1 DNi dX 2 Dt X X Z 2 h h þM @ c F th ðc Þ nþ1 rcnþ1 rNi dX 2 2 ZX h þ M ½rð@ c F el ðc ÞÞnþ1 rNi dX ¼ 0 for i ¼ 1; . . . ; n ; X
2
2.3. Error analysis and convergence In order to analyze the convergence of our numerical approach we studied several examples in detail; a one-dimensional model problem is illustrated in Fig. 3. Because analytical solutions of problem (1) are not known we refer here to a numerical reference solution cref obtained on a very fine mesh with uniform mesh size h ¼ 29 . An initial concentration field (without random noise) is prescribed and its decomposition into phases ca ¼ 0:08 and cb ¼ 0:92 is observed. We evaluate the error between the reference solution and respective solutions on coarser meshes after 1000 time steps of constant step size Dt ¼ 2 105 . Table 1 shows the error measured in the L2 - and the L1 -norm and the corresponding concentrations of a- and b-phase. The results indicate second order convergence (Fig. 3). That actually coincides with the a priori error estimate of finite element
Table 1 Evaluation of the numerical error for a test example. h
kcref ch kL1 ðXÞ 5
ð15Þ
2 26 27 28
2
9.243210 2.1718102 1.1571102 2.2706103
kcref ch kL2 ðXÞ 2
1.709810 4.8621103 1.9447103 4.0361104
ca
cb
0.0292 0.0646 0.0738 0.0746
0.9560 0.9252 0.9195 0.9190
1362
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
approximations applying an implicit Crank-Nicholson discretization of time [13,14]: 2
kcðnDtÞ chn kL2 ðXÞ 6 cðkc0 ch0 kL2 ðXÞ þ h þ Dt 2 Þ;
the concentration profile is randomly perturbed over the domain but fulfills the boundary conditions. The concentration within the phases is ca ¼ cI and cb ¼ cII , respectively.
ð16Þ
where c is independent of mesh size h and time step size Dt. 3. Examples At first we study a model problem at temperature T ¼ 23 T s . Eq. (1) is formulated in dimensionless form with mobility M = I, phase boundary constant k ¼ 2:5 105 and process time t 2 ð0; tÞ. The domain X ¼ ð0; 1Þ2 is discretized with 33 finite elements in every direction. The fractions of both components are assumed to be equal. Consequently, the initial concentration is set c0 ðxÞ ¼ 0:5;
3.1. Short-term investigations In order to illustrate the influence of the interfacial energy on spinodal decomposition we compare computations with two different values of parameter k in Eq. (1), k1 ¼ 4k and k2 ¼ k. The specimen is strain free, F el ¼ 0. The spinodal decomposition of the alloy at different times is displayed in Fig. 4. Here and below the concentration ca and cb are displayed in blue and red, respectively, yellow marks regions of intermediate concentration. The influence of parameter k is twofold. On one hand, it determines the transition region between the
Fig. 4. Decomposition with different interfacial energy, k1 > k2 .
Fig. 5. Comparison between simulation and experimental results.
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
two phases, i.e., the width of phase boundaries. On the other hand, the magnitude of surface tension determines the speed of the decomposition process. The lower the interfacial energy the more the process is dominated by the free Helmholtz energy density FðcÞ. As can be seen from Fig. 4 the mixture decomposes much slower for a higher value of k (upper row). 3.2. Long-term investigations Here we simulate the full process of phase decomposition and coarsening in a strain free case. As outlined in Section 2 the diffusion process shows two stages. Initially we observe a (relatively) fast spinodal decomposition, see t ¼ 0; . . . ; 0:05 in Fig. 6. Consecutively follows a (very) slow process of phase coarsening up to full phase separation.
1363
By now both components have the same concentration, cI ¼ cII , and we compute lamellae-like shapes of separated phases therefore. At next we assume component I with a mass fraction of 35% to be embedded in component II. Hence, the initial concentration is set c0 ðxÞ ¼ 0:65 with a small random fluctuation, see Fig. 7. During phase separation we observe now spherical islands of the aphase emerging in the second-component matrix phase. Long-term simulations of phase separation are a challenge for numerical approximation. On one hand a large number of time steps is necessary to cover the full process, on the other hand the size of time step Dt is limited. Moreover, to compute the non-linear problem in every time step a Newton–Raphson iteration is required. One way to improve the numerical performance is, of course, a time-adaptive numerical scheme. However, a significant role plays also the size of the system of equations in the spatial dis-
Fig. 6. Long-term investigation of a specimen with equal phase fractions.
Fig. 7. Long-term investigation of a specimen with unequal phase fractions.
1364
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
Table 2 Material specific parameters for an eutectic Ag-Cu alloy. g1
GJ m3
5.20
g2
GJ m3
7.27
g3
mol m3
1.11105
v1
GJ m3
2.97
cretization. A comparison of the necessary degrees of freedom nicely illustrates the advantage of our approach. Our 33 33 element mesh has 1225 unknown coefficients for the concentration field. The same mesh with Hermitian polynomials would need more than 18,000 degrees of freedom. A mixed formulation with piecewise linear basis functions still requires 2312 degrees of freedom, not yet guaranteed a fine resolution of phase boundaries. 4. Comparison of simulation and experiment Binary solder alloy is commonly utilized with eutectic concentration of components. The solid mixture then decomposes into phases a and b of specific concentration 0 < ca;b < 1 and not into its components. The application of our model problem to such practical problems is straightforward. However, it requires the material data and Helmholtz free energy function for the specific solder alloy. In the following we analyze a brazing Ag–Cu solder. The aging experiments of the eutectic alloy (71 wt.% Ag, 29 wt.% Cu) have been performed at high temperature, T = 1000 K. The specimen is assumed to be free of mechanical stresses, see [12] for more details. Because of the long time-scale the experiments focus on phase separation and the early stages of coarsening only. In order to simulate the aging process by our numerical approach the thermal component of the free Helmholtz energy density F th , is modelled with a Margules ansatz,
F th ðcÞ ¼ g 1 c þ g 2 ð1 cÞ þ g 3 RT½c lnðcÞ þ ð1 cÞ lnð1 cÞ þ v1 c2 ð1 cÞ þ v2 cð1 cÞ2 ; where R ¼ 8:314 ðJ=mol KÞ is the gas constant. The calculation of the parameters g 1 ; g 2 ; g 3 ; v1 ; v2 from values of a chemical database has been performed in [12] (Table 2). Parameter k is determined
v2
GJ m3
3.01
k ð NÞ
ca
cb
3.431010
0.063
0.945
from the higher gradient coefficients calculated in [12]. Assuming a phase boundary width of about 20 nm the value corresponds well to the known values of surface tension of Ag–Cu alloys. The size of the computational domain is about 1/3 of the experimental one displayed in Fig. 5. The micrographs in Fig. 5 (second row) clearly show the different phases and the observed coarsening process on an area of about 10 lm 10 lm. The light areas denote the Ag-rich a-phase ðca ¼ 0:063Þ and the dark areas represent the Cu-rich b-phase ðcb ¼ 0:945Þ. In the pictures generated from our simulation, see Fig. 5 (first row), the blue areas depict the a-phase and the red areas represent the b-phase. The state shortly after solidification with its fine lamellae-like topology is well reproduced by our simulation results. The phase separation and Ostwald ripening (in early stages) is portrayed in accordance to the correlating time-scale as well. In conclusion, our numerical model shows a good agreement with the experimental results. References [1] J. Cottrell, T.J.R. Hughes, Y. Bazilevs, Wiley, 2009. [2] H. Gomez, V.M. Calo, Y. Bazilevs, T.J.R. Hughes, Comput. Methods Appl. Mech. Eng. 197 (2008) 4333–4352. [3] R.L.J.M. Ubachs, P.J.G. Schreurs, M.G.D. Geers, J. Mech. Phys. Solids 52 (2004) 1763–1792. [4] J.W. Cahn, J.E. Hilliard, J. Chem. Phys. 28 (1958) 258–267. [5] J.W. Cahn, J. Chem. Phys. 30 (1959) 1121–1124. [6] J.W. Cahn, J.E. Hilliard, J. Chem. Phys. 31 (1959) 688–699. [7] J.W. Cahn, Acta Metall. 9 (1961) 795–801. [8] N. Moelans, B. Blanpain, P. Wollants, CALPHAD 32 (2008) 268–294. [9] H. Garcke, Proc. R. Soc. Edinb., Sec. A, Math. 133 (2) (2003) 307–331. [10] A. Fratzl, O. Penrose, J.L. Lebowitz, J. Stat. Phys. 95 (1999) 1429–1503. [11] O.C. Zienkiewicz, R.L. Taylor, McGraw–Hill Book Company, 2003. [12] T. Böhme, Dissertationsschrift, Technische Universität Berlin, 2008. [13] C.M. Elliott, Math. models for phase change problems, in: J.F. Rodrigues, (Ed.), Int. Series Num. Math., vol. 88, 1989, pp. 35–73. [14] C.M. Elliott, S. Larsson, Math. Comput. 58 (1992) 603–630.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Numerical simulation of diffusion induced phase separation and coarsening in binary alloys Denis Anders *, Kerstin Weinberg Lehrstuhl für Festkörpermechanik, Institut für Mechanik und Regelungstechnik, FB11, Universität Siegen, Paul-Bonatz-Straße 9-11, 57068 Siegen, Germany
a r t i c l e Article history: Received 1 October Received in revised Accepted 18 March Available online 18
i n f o 2009 form 11 March 2010 2010 April 2010
Keywords: Binary alloys Solder Spinodal decomposition Phase field Damage
a b s t r a c t Solder materials are solid mixtures and show a variety of micro-structural changes, for example, separation of phases and phase coarsening (Ostwald ripening). The micro-structural changes influence strength and life expectation of these materials. They may be of significant effect, in particular, in very small components such as solder joints in microelectronic packages. In order to analyze the micro-morphological evolution with a diffusion theory of heterogeneous solid mixtures we apply an extended Cahn–Hilliard phase field model. The underlying equations involve a bipotential operator and, accordingly, fourth-order spatial derivatives. Thus, the variational formulation of the Cahn–Hilliard problem requires approximation functions which are piecewise smooth and globally C 1 -continuous; standard finite element techniques cannot be applied here. In our contribution we fulfil the continuity requirement by means of rational B-spline finite element basis functions. To illustrate the versatility of this approach numerical simulations of phase decomposition and coarsening controlled by diffusion and by mechanical loading are discussed and compared with experimental results. Studies of error and convergence are provided. Ó 2010 Elsevier B.V. All rights reserved.
1. The problem of phase decomposition and coarsening Solder joints in microelectronic circuit units are required to provide mechanical as well as electrical connections between different components. The life expectation of these joints determines the reliability of the whole structure. The overall material properties of the solder are affected by the decomposition of the alloy into phases of different composition. Experimental observations show that separation and coarsening of phases generate domains, which are susceptible to mechanical failure. Voids, crack initiation and crack propagation are prevalent observed along phase boundaries. The distribution of components in an alloy can be described by their mass concentration, cI ; cII , etc. In a binary alloy it simply holds cII ¼ 1 cI ¼: c. In order to investigate the evolution of the concentration field cðx; tÞ in a representative domain X within time t we employ a phase field model. Originally introduced by Cahn and Hilliard for diffusion driven phase separation [4–7], this approach has been extended by several contributors to account for additional fields, e.g., thermo-mechanical stresses [8–10,12]. In our context the problem reads: Find c : X ð0; tÞ ! R such that
* Corresponding author. E-mail addresses: [email protected] (D. Anders), [email protected] (K. Weinberg). 0927-0256/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2010.03.030
@c ¼ r j ¼ r ðMrlÞ @t @ ðF th þ F el Þ kDc in X ð0; tÞ ¼ r Mr @c
ð1Þ
holds. Here j is the flux of concentrations driven by the negative product of mobility M and the gradient of the chemical potential l. The initial concentration is given,
cðx; 0Þ ¼ c0 ðxÞ in X;
ð2Þ
and, in order to account for the conservation of mass, homogeneous Neumann boundary conditions are prescribed.
@FðcÞ kDc n ¼ 0 in @ X ð0; tÞ @c rc n ¼ 0 in @ X ð0; tÞ j n ¼ Mrl n ¼ Mr
ð3Þ ð4Þ
Here and above the free Helmholtz energy density F is composed of thermal and elastic components FðcÞ ¼ F th þ F el . The mobility tensor 2 M is assumed to be isotropic and constant. The parameter k ¼ jl is related to surface energy density j and length l of the transition regions between the domains of each phase. Vector n denotes the unit outward normal to @ X. The evolution of the concentration field can roughly be divided in two steps. Starting from an almost homogeneous solution at first phases separate. The driving quantity is here the minimization of the free Helmholtz energy density. In the strain free case
1360
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
ðF el ¼ 0Þ it can be obtained from the experimentally measured Gibbs free energy of the specific solution.
FðcÞ ¼ F th ðcÞ ¼
h ½c lnðcÞ þ ð1 cÞ lnð1 cÞ þ hs cð1 cÞ 2
ð5Þ
solution ch and the discrete test functions uh can be depicted as a linear combination of basis functions N i in Vh ,
uh ¼
n X
ai N i
and ch ¼
n X
i¼1
Here h and hs denote specific material parameters with h ¼ NkT; hs ¼ NkT s where N is Avogadro’s number, k Boltzmann’s constant, T the absolute temperature and T s the critical temperature (solidification temperature). For T > T s Eq. (5) is a convex function. A separation into phases is energetically disadvantageous in comparison to the homogeneous mixture. Thus, the mixture remains stable to all fluctuations. For T < T s Eq. (5) forms a double-well potential, i.e., it has two relative minima and a concave region (spinodal region) in between. The homogeneous mixture is here unstable and will decompose into two phases with concentration ca (a-phase) and cb (b-phase). It holds
@F Fðcb Þ Fðca Þ @F ðca ; TÞ ¼ ðcb ; TÞ: ¼ @c cb ca @c
ð6Þ
After separation the phases start to grow. This coarsening process is driven by the minimization of the interfacial energy density in such manner that the size of the emerged phase regions increases at the expense of their number. The interfacial energy contribution is described by the bipotential operator in Eq. (1). The diffusion based Cahn–Hilliard model is extended by an elastic energy density [9,10]
F el ðc; uÞ ¼
1 ððuÞ th ðcÞÞ : CðcÞ : ððuÞ th ðcÞÞ; 2
ð7Þ
where u ¼ uðx; tÞ denotes the displacement field. Its symmetric gradient is the linearized strain tensor ðuÞ ¼ 12 ðru þ ruT Þ. The elastic properties enter Eq. (7) via elasticity tensor CðcÞ and thermal strain tensor th ðcÞ ¼ ath ðcÞDT; both depending on concentration field cðx; tÞ.
bj N j ;
j¼1
where n is the dimension of Vh . Employing Galerkin orthogonality the discrete problem reads: Find ch 2 Vh such that
Z
Z Z 2 @ch @ F th ðch Þ rch rNi dX Ni dX þ Mk Dch DN i dX þ M @c2 X @t X X Z @F el ðch Þ þM r rNi dX ¼ 0 for i ¼ 1; . . . ; n @c X
ð9Þ
Please note that Eq. (9) is non-linear of ch . The required continuity ch 2 Vh H2 ðXÞ can only be fulfilled by basis functions N j which are piecewise smooth and globally at least C 1 -continuous. Standard finite element basis functions provide C 0 continuity only and cannot be used in this context. Extending the finite element space by Hermitian polynomials is elaborate and restricted to rectangular domains [11]. A mixed finite element approach, suggested, e.g., in [3], is also numerically expensive because it introduces additional unknowns to the primal degrees of freedom. Moreover, the formulation of boundary conditions is not completely clear here. In virtue of these difficulties we follow an approach of [1,2] and apply rational B-splines as finite element basis functions. These basis functions offer high-order accuracy, robustness, and, most importantly, C 1 -continuity. Univariate basic splines (B-splines) are defined recursively for a given order p ¼ 1; 2; 3; . . . and a knot vector n ¼ fni g; i ¼ 1; 2; 3; . . . by the following formula:
Ni;p ðxÞ :¼
niþpþ1 x x ni Ni;p1 ðxÞ þ Niþ1;p1 ðxÞ: niþp ni niþpþ1 niþ1
ð10Þ
Stop function for this recursion is the jump function:
1 if ni 6 x < niþ1
2. Discretization of the problem and numerical solution
Ni;0 ðxÞ :¼
For numerical approximation we derive the primal variational form of the Cahn–Hilliard model. Multiplying Eq. (1) with a test function u, applying Green’s formula (integration by parts) and making use of the homogeneous boundary conditions (3) and (4) give, together with Eq. (2), the problem: Find c 2 H2 ðXÞ such that
For a uniform knot vector n ¼ f0; 1; 2; 3; . . .g and order p ¼ 1 the Bspline functions coincide with piecewise linear finite element basis functions which are C 0 -continuous, see Fig. 1. For higher order, p P 2, shape and properties of B-spline functions differ significantly from standard higher order finite element basis functions. The quadratic B-Spline function of Fig. 1 is globally smooth and C 1 -continuous. Therefore, it will be employed here as basis function in our spatial discretization. To analyze two- and three-dimensional domains multivariate B-spline functions can be derived by tensor products of univariate functions. In two dimensions it holds, see Fig. 2:
Z X
Z Z @c @ðF th þ F el Þ u dX ¼ M r ru dX k DcDu dX @t @c X X
8 u 2 H2 ðXÞ:
ð8Þ
2
Here H ðXÞ denotes the Sobolev space of square integrable functions with square integrable derivatives of first and second order.
0
otherwise:
:
Nij;p ðx; yÞ :¼ Ni;p ðxÞ Nj;p ðyÞ
2.1. Spatial discretization For spatial discretization of Eq. (8) we employ the Galerkin method. To this end the space of test/ansatz functions H2 ðXÞ is approximated by a finite dimensional subspace Vh . The discrete
ð11Þ
ð12Þ
The support of each B-spline Ni;p ðxÞ is compact and contained in the interval ½ni ; niþpþ1 . Hence, the numerical integration required to solve Eq. (9) can be performed with standard algorithms. Additionally, B-spline functions establish a partition of unity,
Fig. 1. Basic splines N i;p of order p = 0, 1, 2 for uniform knot vector n ¼ f0; 1; 2; 3; . . .g.
1361
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
Fig. 2. Two-dimensional B-spline of order p = 2.
Fig. 3. Evolution of concentration in the test example.
n X
Ni;p ðxÞ ¼ 1 8x 2 ½n1 ; nnþpþ1 :
ð13Þ
with fgnþ1 ¼ 12 ðfgn þ fgnþ1 Þ denoting a midpoint evaluation. 2
i¼1
Each spline is a non-negative function, i.e., N i;p ðxÞ P 0, for all x. One of the most important features in our context is ability of B-spline functions to map jumps. The phase boundaries in the Cahn–Hilliard model are typically sharp layers which require a high spatial resolution. Hermitian and Lagrangian interpolation polynomials (used as standard finite element basis) tend to oscillate while representing jumps in a solution function. 2.2. Temporal discretization The temporal discretization employed here is straightforward. The first order time derivative is approximated by finite differences
@c cnþ1 cn ¼ @t Dt
ð14Þ
with (uniform) time step Dt ¼ t nþ1 t n . For time integration we employ an implicit Crank-Nicholson scheme which delivers a second order accuracy [11]. Consequently the fully discrete problem reads:
Z
Z cnþ1 cn Ni dX þ Mk Dchnþ1 DNi dX 2 Dt X X Z 2 h h þM @ c F th ðc Þ nþ1 rcnþ1 rNi dX 2 2 ZX h þ M ½rð@ c F el ðc ÞÞnþ1 rNi dX ¼ 0 for i ¼ 1; . . . ; n ; X
2
2.3. Error analysis and convergence In order to analyze the convergence of our numerical approach we studied several examples in detail; a one-dimensional model problem is illustrated in Fig. 3. Because analytical solutions of problem (1) are not known we refer here to a numerical reference solution cref obtained on a very fine mesh with uniform mesh size h ¼ 29 . An initial concentration field (without random noise) is prescribed and its decomposition into phases ca ¼ 0:08 and cb ¼ 0:92 is observed. We evaluate the error between the reference solution and respective solutions on coarser meshes after 1000 time steps of constant step size Dt ¼ 2 105 . Table 1 shows the error measured in the L2 - and the L1 -norm and the corresponding concentrations of a- and b-phase. The results indicate second order convergence (Fig. 3). That actually coincides with the a priori error estimate of finite element
Table 1 Evaluation of the numerical error for a test example. h
kcref ch kL1 ðXÞ 5
ð15Þ
2 26 27 28
2
9.243210 2.1718102 1.1571102 2.2706103
kcref ch kL2 ðXÞ 2
1.709810 4.8621103 1.9447103 4.0361104
ca
cb
0.0292 0.0646 0.0738 0.0746
0.9560 0.9252 0.9195 0.9190
1362
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
approximations applying an implicit Crank-Nicholson discretization of time [13,14]: 2
kcðnDtÞ chn kL2 ðXÞ 6 cðkc0 ch0 kL2 ðXÞ þ h þ Dt 2 Þ;
the concentration profile is randomly perturbed over the domain but fulfills the boundary conditions. The concentration within the phases is ca ¼ cI and cb ¼ cII , respectively.
ð16Þ
where c is independent of mesh size h and time step size Dt. 3. Examples At first we study a model problem at temperature T ¼ 23 T s . Eq. (1) is formulated in dimensionless form with mobility M = I, phase boundary constant k ¼ 2:5 105 and process time t 2 ð0; tÞ. The domain X ¼ ð0; 1Þ2 is discretized with 33 finite elements in every direction. The fractions of both components are assumed to be equal. Consequently, the initial concentration is set c0 ðxÞ ¼ 0:5;
3.1. Short-term investigations In order to illustrate the influence of the interfacial energy on spinodal decomposition we compare computations with two different values of parameter k in Eq. (1), k1 ¼ 4k and k2 ¼ k. The specimen is strain free, F el ¼ 0. The spinodal decomposition of the alloy at different times is displayed in Fig. 4. Here and below the concentration ca and cb are displayed in blue and red, respectively, yellow marks regions of intermediate concentration. The influence of parameter k is twofold. On one hand, it determines the transition region between the
Fig. 4. Decomposition with different interfacial energy, k1 > k2 .
Fig. 5. Comparison between simulation and experimental results.
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
two phases, i.e., the width of phase boundaries. On the other hand, the magnitude of surface tension determines the speed of the decomposition process. The lower the interfacial energy the more the process is dominated by the free Helmholtz energy density FðcÞ. As can be seen from Fig. 4 the mixture decomposes much slower for a higher value of k (upper row). 3.2. Long-term investigations Here we simulate the full process of phase decomposition and coarsening in a strain free case. As outlined in Section 2 the diffusion process shows two stages. Initially we observe a (relatively) fast spinodal decomposition, see t ¼ 0; . . . ; 0:05 in Fig. 6. Consecutively follows a (very) slow process of phase coarsening up to full phase separation.
1363
By now both components have the same concentration, cI ¼ cII , and we compute lamellae-like shapes of separated phases therefore. At next we assume component I with a mass fraction of 35% to be embedded in component II. Hence, the initial concentration is set c0 ðxÞ ¼ 0:65 with a small random fluctuation, see Fig. 7. During phase separation we observe now spherical islands of the aphase emerging in the second-component matrix phase. Long-term simulations of phase separation are a challenge for numerical approximation. On one hand a large number of time steps is necessary to cover the full process, on the other hand the size of time step Dt is limited. Moreover, to compute the non-linear problem in every time step a Newton–Raphson iteration is required. One way to improve the numerical performance is, of course, a time-adaptive numerical scheme. However, a significant role plays also the size of the system of equations in the spatial dis-
Fig. 6. Long-term investigation of a specimen with equal phase fractions.
Fig. 7. Long-term investigation of a specimen with unequal phase fractions.
1364
D. Anders, K. Weinberg / Computational Materials Science 50 (2011) 1359–1364
Table 2 Material specific parameters for an eutectic Ag-Cu alloy. g1
GJ m3
5.20
g2
GJ m3
7.27
g3
mol m3
1.11105
v1
GJ m3
2.97
cretization. A comparison of the necessary degrees of freedom nicely illustrates the advantage of our approach. Our 33 33 element mesh has 1225 unknown coefficients for the concentration field. The same mesh with Hermitian polynomials would need more than 18,000 degrees of freedom. A mixed formulation with piecewise linear basis functions still requires 2312 degrees of freedom, not yet guaranteed a fine resolution of phase boundaries. 4. Comparison of simulation and experiment Binary solder alloy is commonly utilized with eutectic concentration of components. The solid mixture then decomposes into phases a and b of specific concentration 0 < ca;b < 1 and not into its components. The application of our model problem to such practical problems is straightforward. However, it requires the material data and Helmholtz free energy function for the specific solder alloy. In the following we analyze a brazing Ag–Cu solder. The aging experiments of the eutectic alloy (71 wt.% Ag, 29 wt.% Cu) have been performed at high temperature, T = 1000 K. The specimen is assumed to be free of mechanical stresses, see [12] for more details. Because of the long time-scale the experiments focus on phase separation and the early stages of coarsening only. In order to simulate the aging process by our numerical approach the thermal component of the free Helmholtz energy density F th , is modelled with a Margules ansatz,
F th ðcÞ ¼ g 1 c þ g 2 ð1 cÞ þ g 3 RT½c lnðcÞ þ ð1 cÞ lnð1 cÞ þ v1 c2 ð1 cÞ þ v2 cð1 cÞ2 ; where R ¼ 8:314 ðJ=mol KÞ is the gas constant. The calculation of the parameters g 1 ; g 2 ; g 3 ; v1 ; v2 from values of a chemical database has been performed in [12] (Table 2). Parameter k is determined
v2
GJ m3
3.01
k ð NÞ
ca
cb
3.431010
0.063
0.945
from the higher gradient coefficients calculated in [12]. Assuming a phase boundary width of about 20 nm the value corresponds well to the known values of surface tension of Ag–Cu alloys. The size of the computational domain is about 1/3 of the experimental one displayed in Fig. 5. The micrographs in Fig. 5 (second row) clearly show the different phases and the observed coarsening process on an area of about 10 lm 10 lm. The light areas denote the Ag-rich a-phase ðca ¼ 0:063Þ and the dark areas represent the Cu-rich b-phase ðcb ¼ 0:945Þ. In the pictures generated from our simulation, see Fig. 5 (first row), the blue areas depict the a-phase and the red areas represent the b-phase. The state shortly after solidification with its fine lamellae-like topology is well reproduced by our simulation results. The phase separation and Ostwald ripening (in early stages) is portrayed in accordance to the correlating time-scale as well. In conclusion, our numerical model shows a good agreement with the experimental results. References [1] J. Cottrell, T.J.R. Hughes, Y. Bazilevs, Wiley, 2009. [2] H. Gomez, V.M. Calo, Y. Bazilevs, T.J.R. Hughes, Comput. Methods Appl. Mech. Eng. 197 (2008) 4333–4352. [3] R.L.J.M. Ubachs, P.J.G. Schreurs, M.G.D. Geers, J. Mech. Phys. Solids 52 (2004) 1763–1792. [4] J.W. Cahn, J.E. Hilliard, J. Chem. Phys. 28 (1958) 258–267. [5] J.W. Cahn, J. Chem. Phys. 30 (1959) 1121–1124. [6] J.W. Cahn, J.E. Hilliard, J. Chem. Phys. 31 (1959) 688–699. [7] J.W. Cahn, Acta Metall. 9 (1961) 795–801. [8] N. Moelans, B. Blanpain, P. Wollants, CALPHAD 32 (2008) 268–294. [9] H. Garcke, Proc. R. Soc. Edinb., Sec. A, Math. 133 (2) (2003) 307–331. [10] A. Fratzl, O. Penrose, J.L. Lebowitz, J. Stat. Phys. 95 (1999) 1429–1503. [11] O.C. Zienkiewicz, R.L. Taylor, McGraw–Hill Book Company, 2003. [12] T. Böhme, Dissertationsschrift, Technische Universität Berlin, 2008. [13] C.M. Elliott, Math. models for phase change problems, in: J.F. Rodrigues, (Ed.), Int. Series Num. Math., vol. 88, 1989, pp. 35–73. [14] C.M. Elliott, S. Larsson, Math. Comput. 58 (1992) 603–630.