Numerical experiments of phase separation in ternary mixtures

Mathematics and Computers in Simulation 52 (2000) 41–51

Numerical experiments of phase separation in ternary mixtures M.I.M. Copetti LANA, Departamento de Matemática, Universidade Federal de Santa Maria, 97119-900 Santa Maria, RS, Brazil Received 1 November 1999; accepted 24 November 1999

Abstract We consider an explicit finite element approximation of a model for phase separation in a ternary mixture and describe some numerical experiments in two space dimensions. © 2000 Published by Elsevier Science B.V. All rights reserved. Keywords: Phase separation; Ternary mixture; Critical temperature

1. Introduction In this paper, the numerical approximation of a model for phase separation in ternary alloys is considered. When a one-phase homogeneous system composed of N species, at high temperature and thermal equilibrium, is rapidly cooled to a uniform temperature θ below a critical temperature θ c , where it is unstable with respect to concentration fluctuations, spinodal decomposition takes place: the system separates into spatial regions rich in one specie and poor in the other species and evolves into an equilibrium state with lower overall free energy. The evolution is towards an equilibrium state with two, three or more coexisting phases. The Cahn–Hilliard equation [2] is a continuum model for phase separation in binary systems. An extension to this model for ideal mixtures with more than two components was proposed by Elliott and Luckhaus [6]. A system of nonlinear diffusion equations modelling isothermal phase separation was derived and studied by them. This system was approximated by Blowey, Copetti and Elliott [1] using an implicit finite element scheme and some numerical experiments were performed in the ternary case, in one dimension. A discrete model was considered by De Fontaine [4,5] who described the stability regions in the phase diagram of multicomponent solutions. Spinodal decomposition in ternary alloys was investigated by Hoyt [8] using the structure functions (the Fourier transform of the pair correlation functions). The equations of motion for the three structure functions were obtained by Hoyt using an extension of the nonlinear theory of spinodal decomposition E-mail address: [email protected] (M.I.M. Copetti) 0378-4754/00/$20.00 © 2000 Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 9 9 ) 0 0 1 5 3 - 6

42

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51

developed by Langer [10]. Subsequently, these equations were linearised and solved using various choices of average compositions and atomic mobilities [9]. The equilibrium and dynamic behaviour of multicomponent systems were studied by Eyre [7]. The purpose of this work is to present some numerical simulations of phase separation in ternary systems in two space dimensions. The outline of the paper is as follows: in the next section we introduce some notation and present the mathematical model. In Section 3, we consider and analyse an explicit finite element approximation and in Section 4, we describe some numerical experiments. Throughout the paper, we denote the norms of L∞ () and L2 () by || · ||∞ and || · ||0 , respectively and the semi-norm ||∇v||0 is indicated by |v|1 . Here we take || · || to be the matrix norm induced by the 2-norm. For symmetric matrices, it is the spectral radius. 2. The ternary model We shall consider the numerical approximation of the following system of nonlinear diffusion equations w, u t = L 1w

x ∈ ,

t > 0, X u) − Au u − γ 1u u − 1 – (θ φ(u u) − Au u), w = θ φ(u

(1a) x ∈ ,

t > 0,

(1b)

subject to the initial condition u (xx , 0) = u 0 (xx ),

x ∈ ,

(1c)

R 3 , γ and θ are positive with periodic boundary conditions on the square =[0, l]×[0, l]. Here u , w ∈R constants, A is a 3×3 symmetric matrix and L is a 3×3 positive semi-definite matrix satisfying L1=0. The vector u=(u1 , u2 , u3 )t represents the concentration of three species B, C, D in the mixture with 3 X

ui = 1,

0 ≤ ui ≤ 1,

(2)

i=1

and w =(w1 , w2 , w3 )t is a generalised chemical potential satisfying 3 X

wi = 0.

i=1

The surface represented by (2) is called Gibbs triangle. Each of the three corners of the Gibbs triangle represents 100% of one specie and 0% of the other two species. For example, when u =(0, 0, 1) we can say we have pure phase B. R 3 and we set, for i=1, 2, 3, We indicate by vi the i-th component v ∈R {1}i = 1,

3 X 1X – v= vi , 3 i=1

{φ(vv )}i = φ(vi ) = ln vi .

Introducing u) = θ 9 (u

3 X

1 u, ui ln ui − u t Au 2 i=1

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51

43

the free energy of the system, it follows that u9 (u u)}i + {Au u}i − θ. θ φ(ui ) = {∇u We observe that we have mass conservation for this model: Z Z u (xx , t) dxx = u 0 (xx ) dxx = u m || ∀t, 



where u m is the average composition of the system, 0<{u m }i <1, i=1, 2, 3. The existence of a unique solution to the above problem with Newmann boundary conditions was proved by Elliott and Luckhaus [6]. In this work we take     2 −1 −1 0 1 1 1    2 −1  A = −θc  1 0 1  , L =  −1 3 −1 −1 2 1 1 0 with ||L||=1. It follows that u) = θ (u1 ln u1 + u2 ln u2 + u3 ln u3 ) + θc (u1 u2 + u1 u3 + u2 u3 ) 9 (u and problem (1) may be written as u − θ 1(Lφ(u u)) = 0. ut + γ 12u + θc 1u


Note that this free energy has three equal minima and, for θ < (θc /3), a local maximum at u = 13 , 13 , 13 . In this case, the equilibrium theory presented by Eyre [7] predicts that the ternary system will evolve towards an equilibrium state with phase whose concentration values are given by the points of minimum of 9. When θ is decreased the points of minimum approximate the corners of the Gibbs triangle. If the mean concentration has the form u m =(c1 , c1 , 1–2c1 ), 0≤c1 ≤0.5, we find that the two curves θ = c1 θc ,

θ = 3c1 θc (1 − 2c1 )

define four regions in which H9, the Hessian of the free energy 9 with respect to the two independent concentrations, evaluated at um , is positive definite, negative definite or indefinite (see Fig. 1). Thus, depending on θ , θ c and c1 , the system is unconditionally stable, unconditionally unstable or conditionally stable (fluctuations in composition, in some directions, are unstable). Fig. 2 represents the isothermal section θ c =1, θ =0.3 through the phase diagram of the ternary system. The arrows, at selected points, indicate the directions in concentration which are unstable (see Blowey, Copetti and Elliott [1] and De Fontaine [4] for details).

3. Numerical approximation Let us introduce the notation Hp1 () = {v ∈ H 1 () : v is periodic}

44

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51

Fig. 1. Regions of instability. The selected points indicate some values used in the simulations.

and consider a uniform triangulation T h of  with nodes x i , obtained by dividing the square  into N2 squares with side h = l/N and connecting the north-west vertex to the south-east vertex of each square. We define S h ⊂ Hp1 , by ¯ : χ |τ is linear for τ ∈ T h and χ is periodic}, S h = {χ ∈ C() Sh. and, if vi ∈Sh , 1≤i≤3, we write v ∈S 2 N −1 Taking {χi }i=0 , defined by χ i (xx j )=δ ij , i, j=0, 1,. . . , N2 −1, as a basis for Sh , we consider the following numerical integration scheme based on nodal values:

Fig. 2. Instability directions at selected points when θ c =1 and θ =0.3. At the top corner u=(0, 0, 1), at the left corner u=(1, 0, 0) and at the right corner u =(0, 1, 0). When we have 50% of species C and D and 0% of specie B then u =(0.5, 0.5, 0).

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51 h

(χ , η) =

XZ τ ∈T h

τ

Iτh (χ η) dxx

=

2 NX −1

mi χ(xx i )η(xx i ) ∀χ,

45

¯ η ∈ C(),

i=0

where Iτh is the linear interpolation operator on τ and mi =(χ i , χ i )h =h2 . We use the (·, ·)h inner product to approximate de L2 -inner product and indicate | · |h =[(·, ·)h ]1/2 on Sh . S h , an approximation to u0 , satisfying discrete problem we wish to solve is the following: given U 0 ∈S PThe 3 0 x 0 x n +1 n +1 h x ∈ , find U S , n=0, 1,. . . , such that ∀χ χ ∈S Sh ,W ∈S i=1 Ui (x ) = 1, 0 < Ui (x ) < 1 ∀x U n+1 − U n , χ )h + 1t (L ∇W W n+1 , ∇χ χ ) = 0, (U

 X h W n+1 , χ )h = (θ φ(U U n ), χ )h − (AU U n ,χ χ )h + γ (∇U U n , ∇χ χ ) − 1 – (θφ(U U n ) − AU U n ), χ (W

(3a) (3b)

and 3 X

Uin+1 (xx ) = 1,

0 < Uin+1 (xx ) < 1

∀xx ∈ ,

(3c)

i=1

where h

χ , η) = (χ

3 X

h

χ , ∇ηη ) = (χi , ηi ) and (∇χ

i=1

3 X

(∇χi , ∇ηi ).

i=1

We remark that with this explicit scheme we cannot guarantee the a priori bound in (3c). However, provided 1t was sufficiently small, a sequence {U n } satisfying (3c) was always generated. Note that the sequence generated by (3) satisfies 3 X

Win = 0,

(Uin , 1)h = (Ui0 , 1)h .

i=1

Representing the solution of (3) in terms of the basis functions Sh as n

U }i = {U

Uin

=

2 NX −1

cijn χj ,

n

W }i = {W

j =0

we find

that c ni

=

{cijn }

Win

=

2 NX −1

dijn χj ,

1 ≤ i ≤ 3,

j =0

solves, for i=1, 2,

c n+1 = (I + 1t θc M −1 K − 1t γ M −1 KM−1 K)cc ni − 1t θM −1 K{Lφ(cc n )}i i ×cc n+1 = 1 − c n+1 − c n+1 3 1 2 , P where {Lφ(cc n )}i = 3j =1 φ(cc nj )Lij , Mij = (χi , χj )h , Kij = (∇χi , ∇χj ). Theorem 1. The sequences {U n , W n } generated by (3) satisfy   2 2 U n+1 ) − F h (U U n) 8 θc n+1 32γ F h (U − U n h ≤ 0 + 1 − 4 1t − 2 θ Cn+1 1t L1/2W n+1 1 + U 1t h h 21t

χ ) = (9(χ χ ), 1)h + (γ /2) |χ χ |21 , and Cn+1 =maxi =1,2,3 φ 0 (ξin+1 ) ∞ and ξin+1 depends on Uin+1 where F h (χ n and Ui .

46

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51

Proof. Taking χ = W n+1 in (3a) and χ = U n+1 − U n in (3b) we find that 2 γ 2 1 γ W n+1 ) + (A(U W n+1 , ∇W U n+1 − U n ), U n+1 − U n )h + U n+1 1 − U n 1 1t (L ∇W 2 2 2 γ n+1 1 1 2 U n , U n )h + (AU U n+1 , U n+1 )h = − U n 1 − (AU U 2 2 2 2 1 γ U n , U n )h U n ), U n+1 − U n )h = U n+1 − U n 1 − (AU −(θ φ(U 2 2 1 U n+1 , U n+1 )h − (θφ(U U n+1 ), U n+1 − U n )h + θ (φ(U U n+1 ) + (AU 2 U n ), U n+1 − U n )h −φ(U Noting that Uin ln Uin − Uin+1 ln Uin+1 ≥ (1 + ln Uin+1 )(Uin − Uin+1 ) we find U n ln U n − U n+1 ln U n+1 , 1)h U n+1 ), U n − U n+1 )h ≤ θ (U (θ φ(U 9 (U U n ) − 9 (U U n+1 ), 1)h + 21 (AU U n , U n )h − 21 (AU U n+1 , U n+1 )h . = (9 Thus, 2 γ 2 1 γ W n+1 ) + (A(U U n+1 − U n ), U n+1 − U n )h + U n+1 1 − U n 1 W n+1 , ∇W 1t (L ∇W 2 2 2 2 γ n+1 2 9 (U U n ) − 9 (U U n+1 ), 1)h + θ Cn+1 U n+1 − U n h . ≤ − U n 1 + (9 U 2 Recalling that the triangulation is uniform and observing that 2 n+1 U − Uin 1 = (cc n+1 − c ni )t K(cc n+1 − c ni ) i i i we obtain 2 2 n+1 8 U − U n 1 ≤ 2 U n+1 − U n h . h By Eq. (3a)

2

n+1 U W n+1 0 U n+1 − U n 1 − U n h ≤ 1t L ∇W which yields

2 32γ

2 γ n+1 W n+1 0 − U n 1 ≤ 4 (1t)2 kLk L1/2 ∇W U 2 h and

2 8

2 n+1 U W n+1 0 . − U n h 2 (1t)2 kLk L1/2 ∇W h

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51

47

Therefore,   2 γ 2 γ 2 θc 2 8 32γ 1t 1 − 4 1t − 2 θ Cn+1 1t L1/2W n+1 1 + U n+1 1 − U n 1 + U n+1 − U n h h h 2 2 2 n n+1 h 9 (U U ) − 9 (U U ), 1) . ≤ (9 Recalling the definition of Fh we obtain the result.



Remark 2. The result given by Theorem 1 shows that if the stability condition 1t <

h4 32γ + 8h2 θ Cn+1

holds for n=1, 2,. . . , then Fh satisfies U n+1 ) − F h (U U n) F h (U ≤0 1t which is a decreasing energy property. We can use the fact that we expect to observe phases with concentration values close to the points of minimum of the free energy to estimate the time step. Corollary 3. Suppose that the stability condition given above hold. Then, the sequences {U n , W n } generated by (3) satisfy the following stability estimates n n+1 2 i+1 1/2 i 2 X i 2 U ≤ Const U L W + − U + h 1 1 i=0

where Const is a constant that depends on U 0 . Proof. Summing over n the inequality given in Theorem 1 we find that 9 (U U (9

n+1

h

), 1) + 1t

n+1  X i=1

 2 γ 2 32γ 8 1 − 4 1t − 2 θ Ci 1t L1/2W i 1 + U n+1 1 h h 2

n+1 2 θc X i+1 + − U i h ≤ F h (U 0 ). U 2 i=0

U n +1 ), 1)h >−Const, the result follows. Nothing that (9(U



4. Numerical experiments In this section, we display the results of our numerical simulations of phase separation in ternary systems using the numerical approximation (3). In the simulations, the initial conditions were random perturbations of amplitude 0.05 of the uniform um with the average composition of the first component equal to the average composition of the state u =u second component. A 56×56 mesh was used on the square =[0, 56]×[0, 56] and we took θ c =1, θ =0.3,

48

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51

Fig. 3. Long time evolution of a ternary system with u m = respectively, at time t=50, 250, 900 and 1800.

1 1 1 , , 3 3 3




. Columns 1, 2 and 3 show the evolution of u1 , u2 and u3 ,

γ =0.5 and 1t=0.01. Under this conditions, it follows that the points of minimum of 9 are approximately u B =(0.0551, 0.0551, 0.8898), u C =(0.8898, 0.0551, 0.0551) and u D =(0.0551, 0.8898, 0.0551). The experiments are presented in Figs. 3–5 where the long time evolution of the three components in the system is shown. In all figures, time increases down the columns (t=50, 250, 900, 1800) and in

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51

Fig. 4. Long time evolution of a ternary system with um = respectively, at time t=50, 250, 900 and 1800.

1 1 1 , , 4 4 2




49

. Columns 1, 2 and 3 show the evolution of u1 , u2 and u3 ,

50

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51

Fig. 5. Long time evolution of a ternary system with u m = respectively, at time t=50, 250, 900 and 1800.

2 2 1 , , 5 5 5




. Columns 1, 2 and 3 show the evolution of u1 , u2 and u3 ,

columns 1, 2 and 3, the concentration of the first, second and third component, respectively, is plotted against (x, y).
Fig. 3 shows the time evolution of the ternary system with u m = 13 , 13 , 13 inside the region of global instability (point (a) of Fig. 2). As expected, we observe three phases in the early stages of spinodal decomposition.
In the second experiment performed, u m = 41 , 41 , 21 (point (b) in Fig. 2). Initially, we see two phases, one of them dominated by u3 , and, for some time, the evolution of the system is in the direction u1 =u2 . At t=900, we observe three phases with either u1 , u2 or u3 dominating (see Fig. 4).

M.I.M. Copetti / Mathematics and Computers in Simulation 52 (2000) 41–51


1

51

In the third experiment, u m = 25 , 25 , 5 (point (c) in Fig. 2) and we observe, at short times after the quench, two phases dominated by u1 and u2 with u3 being approximately constant. At t=900, the third phase has appeared. The results are presented in Fig. 5. Similarly to the binary case [3], in all experiments we observed a rapid evolution towards a structure with distinct phases followed by a slow process that involves the coarsening of the regions occupied by the three components. At t=1800, in all experiments performed, the system has a three-phase structure with the concentration of the peaks, for some of the phases, very close to 0.89 as predicted by the equilibrium theory given in [7]. Moreover, the directions of decomposition observed in the early stages of phase separation are in agreement with the results obtained by Blowey, Copetti and Elliott [1] on instability directions. Acknowledgements Research funded by CNPq and FAPERGS. References [1] J.F. Blowey, M.I.M. Copetti, C.M. Elliott, Numerical analysis of a model for phase separation of a multi-component alloy, IMA J. Numer. Anal. 16 (1996) 111–139. [2] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy, J. Chem. Phys. 28 (1958) 258–267. [3] M.I.M. Copetti, C.M. Elliott, Kinetics of phase decomposition processes: numerical solutions to Cahn–Hilliard equation, Materials Sci. Technol. 6 (1990) 273–283. [4] D. De Fontaine, An analysis of clustering and ordering in multicomponent solid solutions — I. Stability criteria, J. Phys. Chem. Solids 33 (1972) 297–310. [5] D. De Fontaine, An analysis of clustering and ordering in multicomponent solid solutions — II. Fluctuations and kinetics, J. Phys. Chem. Solids 34 (1973) 1285–1304. [6] C.M. Elliott, S. Luckhaus, A generalised diffusion equation for phase separation of a multicomponent mixture with interfacial free energy, Preprint series 887, IMA, University of Minnesota, 1991. [7] D.J. Eyre, Systems for Cahn–Hilliard equations, SIAM J. Appl. Math. 53 (1993) 1686–1712. [8] J.J. Hoyt, Spinodal decomposition in ternary alloys, Acta Metall. 37 (1989) 2489–2497. [9] J.J. Hoyt, Linear spinodal decomposition in a regular ternary alloy, Acta Metall. Mater. 38 (1990) 227–231. [10] J.S. Langer, Theory of spinodal decomposition in alloys, Ann. Phys. 65 (1971) 53–86.