Analytical model for equiaxed globular solidification in multicomponent alloys

Available online at www.sciencedirect.com

ScienceDirect Acta Materialia xxx (2015) xxx–xxx www.elsevier.com/locate/actamat

Analytical model for equiaxed globular solidification in multicomponent alloys ⇑ Gildas Guillemot and Charles-Andre´ Gandin

MINES ParisTech, CEMEF UMR CNRS 7635, CS10207, 06904 Sophia Antipolis, France Received 27 January 2015; revised 11 April 2015; accepted 16 April 2015

Abstract—We present a solidification model for equiaxed globular microstructure based on an analytical solution of the diffusion fields in the liquid with cross-diffusion terms for multicomponent alloys. It is demonstrated that the solute profiles for any element can be defined after computation of the eigenvalues and the eigenvectors of the liquid diffusion matrix. The interface velocity and the solute profiles are compared and validated with a numerical solution developed with a front tracking approach. A comprehensive application to the Al–7 wt.%Si–1 wt.%Mg alloy shows the effect of interdiffusion phenomena on the solidification path. This application is developed with an integral approach for the computation of the liquid composition far from the interface. It appears to be the only relevant choice in order to compute the expected driving force. Analytical solutions of the diffusion profile can also be computed for planar and cylindrical geometries. Coupling of this microsegregation model with thermodynamic equilibrium calculations is finally proposed and discussed. Ó 2015 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Solidification; Equiaxed globular microstructure; Microsegregation; Multicomponent alloys; Cross diffusion

1. Introduction In most casting processes of alloys, solidification starts by the development of a primary solid phase from the liquid phase. Redistribution of solute species takes place at the interface between the solid and liquid phases. When the interface follows thermodynamic equilibrium, its growth is mainly controlled by diffusion in the phases. This evolution is of prime importance for the determination of the amount of the liquid phase remaining at lower temperature when secondary solid phases can form. Indeed, the nature and amount of these secondary phases largely influence the final properties of the cast parts. Approaches have initially been proposed in the literature with complete mixing of the liquid phase. This is relevant for small domain in metallic alloys and when the full solidification time is of interest as the ratio of diffusion in solid to liquid is generally very small. These first models rely on the large diffusion or no diffusion in the solid phase, respectively corresponding to the lever rule (LR) [1] or Gulliver–Scheil (GS) [2,3] approximations. Thereafter, models considering finite diffusion in the solid phase have been proposed [4,5] and eventually compared with an analytical solution [6] for binary alloys. These models were extended to account for equiaxed solidification with finite diffusion in the liquid and solid

⇑ Corresponding

author; e-mail addresses: Gildas.Guillemot@ mines-paristech.fr; [email protected]

phases [7–9], multicomponent alloys [10–15], as well as a succession of secondary phases formed from the melt [16]. Full numerical treatments are also made available [17– 19]. However, they are not favored as the main objective of these solidification models is to reach prediction of structure and segregation in casting processing [20–22], thus requiring coupling with macroscopic heat and mass transfers numerical solutions as well as coupling with thermodynamic equilibrium calculations. Analytical solutions are thus preferred in this context in order to keep tractable computational times. The models generally rely on a combination of assumptions in order to estimate the solid/liquid interface velocity, among which: (i) approximate diffusion lengths for the solute species, as a function of the actual grain size and/or the final grain radius, or predetermined shape of the solute profile, (ii) no interdiffusion or cross diffusion phenomena between elements, (iii) no coupling with thermodynamic equilibrium calculations, linearized phase diagrams and fixed diffusion coefficients being assumed. The purpose of this paper is to propose a general analytical model for the growth of equiaxed globular grains in multicomponent alloys with cross diffusion in the liquid phase. We first propose an analytical solution for steady state growth of a single globular grain in semi-infinite media. This solution is compared and validated with a numerical solution based on a one-dimensional (1D) spherical front tracking method. The analytical solution is then extended and coupled with an energy balance in order to

http://dx.doi.org/10.1016/j.actamat.2015.04.030 1359-6454/Ó 2015 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

2

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

Nomenclature B cp D g K L n N r Rg R Rf S t T TL TE TM U v V w  w wl1

eigenvalue vector (m2 s1) mass heat capacity (J kg1 K1) diffusion matrix (m2 s1) volume fraction (–) vector of eigenvectors norm (wt.%) mass latent heat (J kg1) number of points (–) number of added elements (–) radial coordinate (m) molar gas constant (J mol1 K1) growth radius (m) maximal radius (m) surface (m2) time (s) temperature (K) liquidus temperature (K) eutectic temperature (K) melting temperature (K) unitary eigen matrix (–) growth velocity (m s1) volume (m3) mass composition (wt.%) average mass composition (wt.%) mass composition at large distance from interface (wt.%)

model globular growth while heat extraction is taking place. A methodology is introduced to compute the composition far from the interface; being validated with the numerical solution. The analytical model is then applied for various diffusion matrices in the liquid phase, demonstrating the effect of species interdiffusion. Evolution of the solidification path is shown and discussed in the last part. The ternary Al–7 wt.%Si–1 wt.%Mg alloy is chosen for the simulations. 2. Modeling A globular grain is considered in the present model, as schematized in Fig. 1. It consists of a spherical core of radius R made of the primary solid phase (s) growing into a shell of liquid phase (l) delimited by the spherical coordinate r ¼ Rf ð> RÞ. The volume fraction of the solid and liquid phases, respectively gs and gl , satisfy the relation gs þ gl ¼ 1 with gs ¼ ðR=Rf Þ3 . A set of N solute elements i are present into the domain for a specific alloy (N ¼ 2 in Fig. 1). The model is based on solute redistribution at the solid/liquid interface and solute diffusion in both the solid and liquid phases as the main phenomena controlling the radial growth velocity of the globular grain, v ¼ v  ^r, where v is the velocity vector of the interface and ^r is the unit vector directed along the radial direction. Hence it can be seen as a classical 1D microsegregation problem that aims to describe solute profiles in the two phases. As depicted in Fig. 1, the solute compositions evolve from wsi;r¼0 to wsi;r¼R in the solid phase and from wli;r¼R to wli;r¼Rf

~ l1 w

equivalent mass composition at large distance from interface (wt.%)

Greek symbols k parameter defining radius evolution (m s1/2) uT surface heat flux (W m2) UT mass heat flux (W kg3) Dt time step (s) n normalized Landau coordinate (–) DT undercooling (K) q density (kg m3) Superscripts s solid l liquid s=l solid/liquid interface Subscripts i, j component 0 initial value Mathematical notations erfcðÞ complementary error function E1 ðÞ exponential integral function kk norm of vector

in the liquid phase, where interfacial compositions wsi;r¼R and wli;r¼R are denoted ws=l and wl=s i i , respectively. Density of the phases being assumed equal and constant in all phases, the conservation equation for solute species i in phase a relates the temporal variation of its mass fraction wai at radial coordinate r to its diffusion flux  PN  a a : D grad w ij j j¼1 " #  N   X @wai  a a ¼ div Dij gradwj with a 2 fs;lg and i 2 f1;N g @t r j¼1 ð1Þ

with the same assumption for the density of the phases, the solute balance at the s=l interface writes: N  N    X  X  l=s s=l v wi  wi ¼ Dlij grads=l wlj  Dsij grads=l wsj j¼1

j¼1

with i 2 f1;N g

ð2Þ

where the gradients of solute in liquid and solid phases are computed at the s=l interface (i.e., at position r ¼ R). As can be seen from the above equations, species interdiffusion is accounted for, with Daij the cross diffusion coefficient that represents the effect of the solute gradient of element j on the solute flow of element i in phase a. Due to the spherical geometry of the system, solute flux is null at the center of the domain (r = 0). The average composition of any element over the two-phase system does not

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

3

Fig. 1. Schematic illustration of the globulitic grain geometry and its solute profile for component 1 and 2 in solid (s) and liquid (l) phases (i.e. N = 2 in the present case). Finite diffusion in both (s) and (l) phases is considered with a fixed average composition respectively equal to w1;0 and w2;0 . A homogeneous surface heat extraction with flow rate uT is assumed in the present model.

change with time for a close system with respect to mass exchange. This also corresponds to a no flux condition at the periphery of the liquid domain (r = Rf ). The two boundary conditions write: N  X

 Dsij gradr¼0 wsj ¼ 0

with

i 2 f1; N g

ð3aÞ

j¼1 N  X

 Dlij gradr¼Rf wlj ¼ 0 with

i 2 f1; N g

ð3bÞ

j¼1

Thermodynamic equilibrium is also assumed at the s=l interface, meaning that the compositions at the solid/liquid interface are linked by the tie-lines of the multicomponent phase diagram. This condition is expressed using the definition of the liquidus surface, T L , depending from the interfal=s cial composition of all N solute elements fwi g16i6N , and the partition coefficient at the s=l interface, k i : h i l=s ð4aÞ T L ¼ F L fwi g16i6N s=l

l=s

k i ¼ wi =wi

with

i 2 f1; N g

ð4bÞ

At time t = 0 s, an almost fully liquid domain is considered with a uniform composition equal to the nominal alloy composition, fwli;t¼0 s g16i6N ¼ fwi;0 g16i6N , with i 2 f1; N g. The solid is yet assumed to be already present as a very small nucleus center at position r = 0 m with uniform composition for each solute species i, wsi;t¼0 s ¼ k i wi;0 , with i 2 f1; N g. Considering the size of the system used to represent a single globular equiaxed grain, a uniform temperature, T , can be assumed at any time. This is justified by the small magnitude of the Biot number defined by the ratio of the heat flux through the system boundary to the thermal exchange within the system. As for the density, q, the heat capacities per unit mass, cp , are assumed constant and equal in all phases. Furthermore considering a constant latent heat per unit mass for the phase transformation, L, a simple energy balance is expressed as:   @T @gs S uT qcp  qL ¼ ð5Þ @t V @t where uT ¼ uT  ^r is the radial component of the outward directed heat extraction rate vector per unit surface, S is the surface of the system where heat extraction takes place

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

4

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

and V is the volume of the system. Eq. (5) permits to compute a temperature when a heat extraction rate uT is imposed. However, a simpler temperature history is also used in this contribution, simply imposing a fixed temperature, T ¼ T 0 with T0 < TL where T L ¼ F L ½fwi;0 g16i6N . 2.1. Analytical

where K ij and Bj are unknown fixed coefficients still to be determined and the function F is defined by Eq. (A1b) in Appendix A. As detailed in Appendix B, introduction of this mathematical expression for wli ðrÞ in Eq. (1) with a ¼ l and rearrangement gives: ! " # N N X X r r l 0 Bj K ij  Dik K kj F pffiffiffiffiffiffi ¼ 0 3=2 Bj t j¼1 2ðBj tÞ k¼1 i 2 f1; N g

N X

Dlik K kj ¼ Bj K ij

with

ði; jÞ 2 f1; N g

ð8Þ

k¼1

An analytical solution is hereafter given for the problem defined with the previous set of equations. For that purpose, we have to add a series of approximations. A steady state growth regime is sought for the composition profiles, that requires adopting a constant temperature of the system, T 0 . The temperature is chosen below the liquidus of the alloy, T L , to permit the coexistence of the solid and liquid phases as well as a driving force for the growth of the interface. Note that this condition is not the same as imposing an adiabatic condition (i.e. uT ¼ 0 W m2 ). Thus Eq. (5) is not considered in the present situation as temperature is simply kept unchanged and is no longer an unknown of the problem. The steady state growth regime also requires a semi-infinite liquid medium or an infinite value of the domain size Rf . Under such conditions and further neglecting the effect of the interface curvature on thermodynamic equilibrium, fixed compositions are indeed expected at the s=l interface and for the far-field composition corresponding to Rf ! 1, wli;1 . According to the description of the model, the initial value for this composition is the nominal alloy composition, wi;0 . It is worth noting that because the steady regime tends toward constant interfacial compositions, solute fluxes in the solid are null. Finally, constant diffusion matrices also have to be assumed. Kirkaldy et al. [23] proposed a solution to a similar problem written in the 1D Cartesian coordinate system. This solution was later extended by Hunziker [24] for a stability analysis of a planar solid/liquid interface in a multicomponent system. The same approach can be further developed to solve the diffusion equation (Eq. (1)) and compute the solute profile around a spherical geometry for a multicomponent alloy with cross diffusion. Thus, the solute profile for any element i is given by a N -linear combination of the solution corresponding to the 1-component expression. In the 1D spherical coordinate system of interest in the present study, the 1-component solution was provided by Aaron et al. [25] in continuity with the heat diffusion solution given by Carslaw and Jaeger [26]. This solution is reported in Appendix A. The solute profile for any element i of the N -component problem writes: ! N X r l l with i 2 f1; N g ð6Þ wi ðrÞ  wi;1 ¼ K ij F pffiffiffiffiffiffi Bj t j¼1

with

 pffiffiffiffiffiffi The derivative functions F 0 r= Bj t and the first coefficient r=ðBj tÞ3=2 correspond to a set of independent functions non equal to 0. Consequently, the only solution to the previous equation corresponds to:

ð7Þ

This equation is no more than the one defining the eigensystem of the diffusion matrix Dl . The value K ij is the component i of the j eigenvector, K j , for matrix Dl . The associated eigenvalue of vector K j is the value Bj . It should be noticed that eigenvectors are defined with a specific unknown multiplicity factor. Thus the norm of each eigenvector, kK j k, is still to be determined. It is proposed to define later K ij by the product kK j kU ij where U ij is the known i-component of the j-unitary eigenvector, U j , and kK j k an unknown coefficient still to be determined. The set of U ij coefficients define the eigenmatrix U. As for the 1-component problem, the only possible expression for the growth rate in a multipffi component alloy is expressed by Eq. (A2) (i.e. R ¼ k t). Consequently Eq. (6) can be estimated at the s=l interface as a function of the k parameter which corresponds to the following equality: ! N X k l=s l wi  wi;1 ¼ U ij F pffiffiffiffiffi kK j k with i 2 f1; N g ð9Þ Bj j¼1 This is nothing but the i-components of a matrix equation that associates vector Dwl corresponding to the differences l=s wi  wli;1 to the vector K corresponding to the kK i k unknown coefficients. The inverse of the matrix relation gives a direct relation between the two vectors if the k parameter is determined: CðkÞDwl ¼ K

ð10Þ

where the component CðkÞij of matrix CðkÞ can be simply expressed with the inverse, U 1 , of the eigenmatrix U as: U 1 ij  CðkÞij ¼  k F pffiffiffiffiffi Bi

with

ði; jÞ 2 f1; N g

ð11Þ

The factor k is still to be determined with the solute balance at the interface written in Eq. (2). It leads to a relation similar to Eq. (A3a) previously obtained for 1-component:   N X N 2X U jk 0 k s=l l ffiffiffiffiffi p ffiffiffiffiffi p wl=s  w ¼  D kK k F k i i k k¼1 j¼1 ij Bk Bk with i 2 f1; N g

ð12Þ

Due to the relation defining vectors K as eigenvector, Eq. (8), we obtain: s=l wl=s i  wi ¼ 

  N pffiffiffiffiffi 2X k kK k kU ik Bk F 0 pffiffiffiffiffi k k¼1 Bk

with i 2 f1; N g

ð13Þ

The kK k k parameters are expressed in Eq. (10) and can be replaced in the previous equation in order to obtain the following relation:

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx s=l wl=s i  wi

 N X N  X k  l=s 1 ¼ U ik U kj G pffiffiffiffiffi wj  wlj;1 Bk j¼1 k¼1 with

i 2 f1; N g

ð14Þ

where the function G is given by Eq. (A3b). This relation is no more than the definition of a matrix equation between the vector of the composition variation at the interface, l=s s=l Dwl=s , corresponding to the differences wi  wi , and Dwl : KðkÞDwl ¼ Dwl=s

ð15Þ

where the component KðkÞij of the matrix KðkÞ is expressed as:   N X k ffiffiffiffiffi p KðkÞij ¼ U ik U 1 G with ði; jÞ 2 f1; N g ð16Þ kj Bk k¼1 It should be emphasized that Eq. (15) can be directly solved to compute the unknown k parameter. However such resolution requires considering thermodynamic equilibrium at l=s the interface. Consequently the set of values fwi g16i6N s=l and fwi g16i6N are linked and correspond to the same tie-line in the phase diagram through Eq. (4b) valid for each i element. Moreover the current temperature of the domain, T , corresponds to the liquidus one, for the compositions fwl=s i g16i6N , Eq. (4a). Due to its non-linearity, the resolution of Eq. (15) is based on a simplex method [27] in order to minimize the difference between the N computed k values on each line (i.e. solute). The simplex here is the N + 1 vertices polyhedron. Each vertex corresponds to a set of N compositions in each element. The computation of the growing steady state of the multicomponent alloy for temperature T finally corresponds to the following steps: (i) diagonalization of the diffusion matrix, Dl , in order to compute the set of eigenvalues, Bj , and normalized eigenvectors, U j , (ii) computation of the lambda parameter, k, with the resolution of Eq. (15) considering the tie-lines of the alloys and liquidus temperature, Eq. (4), thus providing with an estimation for the interface velocity, (iii) direct computation of the coefficients kK i k with Eq. (10), thus accessing the composition profiles in the liquid phase for all elements. 2.2. Numerical A direct numerical solution of the set of Eqs. (1)–(5) has also been developed. It is based on the Landau transform [28], which has largely been applied to solidification problems [17,18,29,30]. As sketched in the bottom part of Fig. 1, radial coordinates, r, are transformed in normalized coordinates, n, corresponding to the values ns ¼ r=R and nl ¼ ðr  RÞ=ðRf  RÞ in the solid and liquid regions, respectively. Considering a spherical geometry, Eq. (1) can be rewritten with the normalized coordinates for each domain as:     s N @wsi  ns v @wsi  X 1 Dij @ r2 @wsj  ¼ þ R @ns t j¼1 r2 R @ns R @ns t @t ns with

i 2 f1; N g

5

!   l N @wli  ð1  nl Þv @wli  X 1 Dij @ r2 @wlj  ¼ þ  Rf  R @nl t j¼1 r2 Rf  R @nl Rf  R @nl  @t nl t



ð17aÞ

with

i 2 f1; N g

ð17bÞ

The boundary conditions at the limits of the domains, Eq. (3), then take the form:  N X @wsj  Dsij s  ¼ 0 with i 2 f1; N g ð18aÞ @n t;ns ¼0 j¼1 N X j¼1

Dlij

 @wlj   @nl 

¼ 0 with

i 2 f1; N g

ð18bÞ

t;nl ¼1

Eqs. (17) and (18) are solved to determine composition fields everywhere but on the interface. Indeed, the interfacial compositions are computed by solving the solute balance, Eq. (2). In Landau coordinates, this is expressed by the following equation:   N N   @wlj  @wsj  1 X 1X l=s s=l l s  v wi  wi ¼  D þ D  Rf  R j¼1 ij @nl  l R j¼1 ij @ns t;ns ¼1 t;n ¼0 with i 2 f1; N g

ð19Þ

As previously explained for the determination of the value of the analytical model parameter k, the N -equations (19) are solved considering simultaneously the set of N s=l l=s tie-lines linking interfacial compositions fwi g and fwi g as well as values for the interfacial liquid compositions given by the temperature, Eq. (4). The following steps are successively developed in a time stepping procedure: (i) resolution of Eqs. (17) and (18) to compute the composition fields in the solid and liquid domains except at the solid/liquid interface using a regular 1D mesh and the finite difference method in order to estimate first and second spatial derivatives, together with an implicit time-stepping scheme, (ii) resolution of Eq. (19) to compute the interface velocity, v, considering tie-lines and liquidus temperature expressions, Eq. (4), (iii) in case a heat extraction rate is imposed at the boundary of the system (i.e., r ¼ Rf ), the remaining fraction of solid entering Eq. (5) is computed and the temperature is updated; otherwise a temperature history is simply imposed. Steps (i)–(iii) are computed without iteration. This is satisfying providing sufficiently small time steps are used. The model parameters are the time step, Dt, the number of mesh in the solid and liquid domains, respectively ns and nl , as well as the initial solid domain size, R0 , taken equal to an arbitrarily small value, 107 m. 2.3. Comparisons The analytical solution is compared with the numerical solution for the multicomponent alloy Al–7 wt.% Si– 1 wt.% Mg. The phase diagram of this ternary alloy is linearized. Constant values are proposed for the liquidus slopes, mi ¼ @T L =@wl=s i , and partition coefficients, k i (Eq. (4b)), with respect to each solute species i, with i 2 f1; N g. These values, listed in Table 1, are extracted at

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

6

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

Table 1. Properties of the ternary Al–7 wt.%Si–1 wt.%Mg alloy [31]. Melting temperature Nominal composition in Si Nominal composition in Mg Liquidus slope for Si element Liquidus slope for Mg element Partition coefficient for Si element Partition coefficient for Mg element Liquidus temperature Eutectic temperature Heat capacity Latent heat Density

TM wSi;0 wMg;0 mSi mMg k Si k Mg TL TE cp L q

Table 2. Parameters for comparison of the numerical and analytical solutions (Figs. 2 and 3). Radius of the droplet Fixed temperature of the droplet Number of divisions in the solid domain Number of divisions in the liquid domain Time step

Rf T0 ns nl Dt

1  103 874 10 000 10 000 105

(m) (K) (–) (–) (s)

the liquidus temperature of the alloy from a thermodynamic equilibrium database [31]. With the purpose to reach comparison, the system temperature is kept unchanged at T = 874 K and a large value of the system size, Rf ¼ 1 mm, is used. Parameter values for the numerical simulation are given in Table 2. Small values are chosen for the diffusion coefficients in the solid phase, Ds , required for the numerical solution. In fact, these values are sufficiently small to consider that no-diffusion occurs in the solid phase. Table 3 also gives the values of the full diffusion matrix for the liquid, Dl , used hereafter.

936.38 7 1 6.67 4.24 0.107 0.221 885.45 830 1144.2 3.437  105 2552

(K) (wt.%) (wt.%) (K wt.%1) (K wt.%1) (–) (–) (K) (K) (J kg1 K1) (J kg1) (kg m3)

Fig. 2 presents the computed evolution of the interfacial compositions and position of the s=l interface as a function of the square root of time. As expected, steady state is reached after a short transient regime during which the interfacial compositions for silicon and magnesium evolve toward constant values. These values are the same as the one deduced from the analytical solution with the resolution of Eqs. (4) and (15) and given in Fig. 2(a) and (c). Fig. 2(b) shows the linear dependence of the s=l interface position, R, with the square root of time expected from Eq. (A2). The slope of the curve also rapidly approaches the analytical prediction of the proportionality growth parameter given in Fig. 2(b). Simultaneously, the velocity evolves as a function of the inverse of the square root of time as shown in Fig. 2(d). The present comparison is further developed in Fig. 3 where the solute profiles in the liquid phase are plotted at three times: 0.25 s, 5 s and 25 s. The corresponding interface positions are marked with crosses in Fig. 2(b) showing that they correspond to the steady regime. The analytic solutions corresponding to Eq. (6) perfectly superimpose onto the numerical

Table 3. Diffusion matrices in the solid and liquid phases [31,32]. Eigenvalues and eigenvectors liquid diffusion matrices are also provided for the first (Si) and the second (Mg) component. Solid phase Diffusion matrix Ds ð1012 m2 s1 Þ   1 0 0 1 Liquid phase Diffusion matrix Dl

Dlinf

9

ð10 m2 s1 Þ   2:777 3 2 5:914

Dlsup

9

Dldia

9

ð10 m2 s1 Þ   2:777 0 2 5:914

ð10 m2 s1 Þ   2:777 3 0 5:914

ð109 m2 s1 Þ   2:777 0 0 5:914

B linf

B lsup

B ldia

Eigenvalues Bl 9

ð10 m s Þ

ð10 m s Þ

ð10 m s Þ

ð109 m2 s1 Þ

ð 1:045 7:146 Þ

ð 2:277 5:914 Þ

ð 2:277 5:914 Þ

ð 2:277 5:914 Þ

Ul  inf  0:876 0 0:482 1

U lsup   1 0:636 0 0:771

U ldia   1 0 0 1

2

9

1

2

1

9

2

1

Eigenvectors Si Mg

Ul 

0:925 0:525 0:380 0:851 Si Mg



Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

7

Fig. 2. Time evolutions of the interfacial composition in the liquid phase of an equiaxed globular grain growing in an undercooled liquid with fixed composition. Evolutions are shown for both silicon (a) and magnesium (c) elements. The interfacial position and velocity are respectively shown in (b) and (d). Details on solidification conditions are listed in Tables 1–3. Exact values from the analytical steady state regime (dashed curves and labels) are compared with the numerical solution (plain curves).

solutions deduced from the front tracking model. Also note that at time t = 25 s, the solute layer in the liquid has just reached the outer boundary of the liquid domain for the silicon element. The diffusion profile for magnesium also shows a non-monotonic behavior typical of the effect of interdiffusion between species frequently discussed in the literature [23,24]. These comparisons validate the present analytical solution for equiaxed globular growth in a multicomponent alloy with cross diffusion.

3. Results

where D~ ~ li;1 for curwl is the vector of the differences wl=s i w l ~ i;1 is the composition far from the s=l rent time t and w interface. We shall see later various methodologies to estimate this quantity. Once k is determined, all other variables and fields can be computed. The time derivative of Eq. (A2) also permits determining the interface velocity: v¼

k2 2R

ð21Þ

The radius of the globular grain is then increased by (v  Dt) and the solid fraction is also updated from the new interface position. Further neglecting the diffusion in the solid  si , is computed phase, the average composition in the solid, w step by step by a numerical integration of the interfacial composition: Z t s  si ðtÞ ¼ gs w ws=l with i 2 f1; N g ð22Þ i dg

On the basis of the above validated analytical solution, we now propose to develop a model for equiaxed globular solidification in multicomponent alloys for a fixed geometry as previously presented in Fig. 1. This requires extending the solution of the previous equations by assuming a quasi-steady state regime during which cooling takes place as well as solutal interactions in the liquid. In other words, at the periphery of the system corresponding to the radial coordinates r ¼ Rf , not only a heat extraction rate is now present (Eq. (5) holds) but composition also locally varies due to the solute profiles of the chemical species that extend up to this outmost boundary. For that purpose Eq. (15) is solved in a time-stepping scheme in order to re-evaluate the growth rate parameter k when the temperature and the far-field composition vary. It is rewritten as:

The present hypothesis of no solid diffusion is consistent with the value of the solid diffusion coefficient proposed in Table 3. However, as it is demonstrated in the discussion part, few differences are observed between the computations developed with and without diffusion in the solid phase for the present ternary alloy. A global solute balance then permits to access the aver li : age liquid composition w

KðkÞD~ wl ¼ Dwl=s

 si þ gl w  li ¼ wi;0 gs w

ð20Þ

0

with i 2 f1; N g

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

ð23Þ

8

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

Finally the energy balance, Eq. (5), is solved to update the temperature of the system. All variables previously pre~ li;1 that needs to sented are known at any time t, except w be determined. 3.1. Far field approximations ~ li;1 has to be linked to the The far-field composition w other variables in order to estimate the driving force for the growth of the grain. Various approaches have been proposed in the literature [7,8,33,34]. The three choices studied hereafter are schematized in Fig. 4. They are all based on the  li , which is deduced known averaged liquid composition, w from the global solute balance equation, Eq. (23), at any time of the solidification process. According to Eq. (23),  li only depends on the actual fraction computation of w and average composition of the solid phase and is thus the  li is reflected same for the three choices. The uniqueness of w ~ li;1 are schein Fig. 4 while different values of compositions w matized for the far-field composition. Choice 1: Average liquid composition (Fig. 4(a)). Estimation is simply given by the average value of the liquid phase: w ~ li;1 ¼ w  li

ð24Þ

Choice 2: Composition at radial coordinates Rf (Fig. 4(b)). The composition computed at position r ¼ Rf from the solute profile, Eqs. (6) and (10), is set equal to the average  li . This is composition value of the liquid phase: wli;r¼Rf ¼ w illustrated in Fig. 4(b) by the black cross showing the position where solute profile crosses the average liquid composition at Rf . The following expression can thus be obtained ~ li;1 , where the values Dwlk can be estimated at the prefor w vious time step solution: ! Rf k pffiffiffiffiffi F N X N X R Bj l l 1 ! Dwlk ~ i;1 ¼ w i  U ij U jk ð25Þ w k j¼1 k¼1 F pffiffiffiffiffi Bj Choice 3: Integration of the composition profile (Fig. 4(c)). The average composition, w  li , can also be linked to the integration of the approximated solute profile (Eq. (26)) between positions R and Rf , hwli ðrÞi½R;Rf  . Graphically, this means that the blue area under the approximated profile shown in Fig. 4(c) should correspond to the gray area  li . This under the profile defined by the average value w choice is mathematically expressed as: Z Rf 1  li ¼ hwli ðrÞi½R;Rf  ¼ 4 w wli ðrÞ4pr2 dr ð26Þ pðR3f  R3 Þ R 3 where the liquid composition, wli ðrÞ, is expressed as: ! rk pffiffiffiffiffi F N X N X R Bj l l 1 ! Dwlk ~ i;1 þ wi ðrÞ ¼ w U ij U jk ð27Þ k j¼1 k¼1 F pffiffiffiffiffi Bj

Manipulation of the two above equations provides with an ~ li;1 : estimation of w N X N X 1 U ij U 1 jk 3 pðR3f  R3 Þk j¼1 k¼1 3   pffiffiffiffiffi pffiffiffiffiffi I Rf k; R Bj  I Rk; R Bj    Dwlk k F pffiffiffiffi

~ li;1 ¼ w  li  4 w

ð28Þ

Bj

where the function Ið; Þ is defined by: 2 3 x2    p ffiffiffi 2p 6 x 2 7 Iðx; yÞ ¼ 42yð2y 2  x2 Þe 4y þ px3 erfc 5 3 2y

ð29Þ

Simulations are performed with the goal to compare the effect of the above choices for the estimation of the far-field liquid composition. The ternary Al–7 wt.%Si– 1 wt.%Mg alloy is considered, properties in Table 1 being used again. The final radius of the grain, the surface heat flux and the initial temperature are listed in Table 4, together with other model parameters. Diffusion matrices for the solid phase, Ds , and in the liquid phase, Dl , are given in Table 3. Simulation results are shown in Fig. 5. The blue curves labeled [A] correspond to the analytical solution using the three choices to compute the far-field liquid composition. They are also compared with the black curves labeled [N] deduced from the numerical solution, parameters for the simulation being also listed in Table 4. It is clearly demonstrated that a large deviation exists between the three choices. The cooling rate (Fig. 5(a)) changes as soon as a noticeable latent heat is released, i.e. around 2.5 s. However, the recalescence predicted by the numerical solution is only reproduced with choice 3 (plain line). For choices 1 and 2, only a change of the slope of the cooling history takes place without the expected temperature increases. Since the temperature evolution is the result of an energy balance between the latent heat released due to the growth of the globular grain (term @gs =@t proportional to v in Eq. (5)) and the heat extracted through the outer boundary of the system, uT , this difference is directly linked to the s=l interface velocity and hence to the difference l=s ~ li;1 . Fig. 5(b) and (c) shows the time evolution of wi  w the interfacial compositions for Si and Mg, respectively. Comparison with the numerical solutions retrieves the previous observation made for the temperature evolution: introduction of choice 3 in the analytical solution gives l=s an evolution of the interfacial composition, wi , close to the numerical solution. This is the main reason for the prediction of the recalescence with the analytical solution. Finally, Fig. 5(d) and (e) compares the far-field liquid compositions. Deviations from the numerical solution appear between curves, even with choice 3. In fact, the numerical solution is closer to choice 1 at longer times. Indeed, values ~ li;1 for the numerical solution have been extracted from of w the computed composition profile at position, r ¼ Rf , wli;Rf . Since the composition gradient in the liquid phase vanishes upon propagation of the interface, composition wli;Rf neces li and sarily tends toward the average liquid composition w hence the better comparison with choice 1 when

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

9

Fig. 3. Composition profiles in the liquid for (a) Si and (b) Mg at times 0.25 s, 5 s and 25 s comparing (plain lines, [N]) the numerical solutions with (dashed lines, [A]) the analytical solutions. The vertical lines correspond to the radial position of the solid/liquid interface, also shown by crosses in Fig. 2(b).

(a)

(b)

(c)

~ li;1 that defines the driving force for the growth of the s=l interface Fig. 4. Schematic illustration of the possible choices to compute the composition w ~ li;1 ¼ w  li , (b) choice 2: the composition at position r ¼ Rf (black cross), wli;r¼Rf ¼ w  li , with (a) choice 1: the average composition in the liquid phase, w  li . The gray curves represent the actual solute profile while the blue and (c) choice 3: the average composition in the liquid phase, hwli ðrÞi½R;Rf  ¼ w ~ li . The blue dot is the actual curves are the approximated profiles. The horizontal dotted gray lines are the average liquid compositions, w ~ li;1 for the various choices. Please note that this latter composition is obviously badly schematized as it corresponds to the farapproximation for w field value, i.e. at infinite distance from the solid/liquid interface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

~ li;1 . However, these differences reveal that, if considering w the goal is to reach a good estimation of the interfacial composition and hence the growth velocity, the composi~ li;1 is a bad choice. tion at the boundary of the domain for w The comparisons definitely demonstrate that largest driving

force for the interface and thus largest growth velocity is obtained for choice 3. This latter choice is thus the only ~ li;1 in one kept for the estimation of the key parameter w the rest of this paper for the use of the quasi-steady state approximation of the analytical solution.

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

10

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

Table 4. Geometry and cooling conditions used for the globulitic solidification of the domain and simulation parameters proposed both for the analytical model and numerical model (Figs. 5 and 6).

Initial temperature

Rf UT uT uT uT T0

2  104 5000 851 1276 2552 885.45

(m) (W kg1) (W m2) (sphere) (W m2) (cylinder) (W m2) (plan) (K)

Analytical model Time step

Dt

2  104

(s)

Numerical model Number of divisions in the solid domain Number of divisions in the liquid domain Time step

ns nl Dt

2000 2000 2  106

(–) (–) (s)

Radius of the droplet Mass heat flow Surface heat flow (equivalent to UT )

Fig. 5. Evolution of temperature, T , (a), interfacial compositions, wl=s ~ li;1 (d, e). Composition evolutions are i , (b, c) and far-field compositions, w shown both for silicon (b, d) and magnesium (c, e) elements. Results are detailed for the three choices proposed in the present model (blue lines – [A]). The front tracking model (black line – [N]) has also been applied to the cooling simulation and the evolutions are superimposed to those corresponding to the analytical approach. Composition at infinite for the front tracking simulation in (d) and (e) are considered as the ones estimated to the position r ¼ Rf . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

11

~ li;1 , (d, e). Compositions are Fig. 6. Evolution of temperature, T , (a), interfacial compositions, wl=s (b, c), and compositions at a large distance, w i l l l shown both for silicon (b, d) and magnesium (c, e) elements. Curves are represented for D (black), Dinf (red), Dsup (green) and Dldia (blue) diffusion matrices for both the front tracking numerical simulations (plain curves – [N]) and the present analytical approach (plain curve – [A]). Composition at infinite for the front tracking simulation in (d) and (e) is considered as the ones estimated to the position r ¼ Rf . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3.2. Comparisons with the numerical solution The analytical and numerical simulations have been performed for the three other liquid diffusion matrices presented in Table 3. The same small diffusion matrix, Ds , is used for the solid phase. Simulation parameters are also the same as the ones described in Tables 1 and 4. Results are presented in Fig. 6 where the same variables time-evolution as in Fig. 5 are drawn. Reference results shown in Fig. 5 for the simulations corresponding to the full diffusion matrix Dl are indeed reproduced for a better discussion on the effect of the diffusion matrix. Simulations labeled Dldia correspond to a simple diagonal matrix with no cross diffusion terms. It is the usual approximation found in the literature. For a distinct and progressive consideration of the cross diffusion terms, matrices Dlinf and Dlsup are also considered, in which the off-diagonal terms DlSi Mg and DlMg Si are set to zero, respectively

(Table 3). This methodology is extracted from recent analysis on solid-state precipitation in multicomponent alloys [35]. The first observation is that the analytical and numerical simulations are very close to each other when considering the temperature and interfacial composition for a given choice of the diffusion matrix. This demonstrates that the analytical model enables to reproduce the equiaxed globular grain growth whatever the diffusion matrix when considering choice 3 for the estimation of the far-field liquid ~ li;1 lead composition. Simulations with choices 1 and 2 for w to similar conclusions as in Fig. 5 for the others matrices and are consequently not shown here. It is worth noticing that, for Si, the interfacial composil=s tion wSi (Fig. 6(b)) is always higher than the far-field com~ lSi;1 (Fig. 6(d)). This is not the case for Mg as position w shown by comparing curve in Fig. 6(c) with that in Fig. 6(e): interfacial composition wl=s Mg is higher than

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

12

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

~ lMg;1 for cases Dldia and Dlsup and lower for cases Dl and w Dlinf . The latter situation with Dl is in fact shown in Fig. 3 where the solute gradients at the s=l interface in the liquid are negative for Si and positive for Mg. According to Eq. (2), the interface velocity and thus the event of the recalescence are directly influenced by the overall diffusion fluxes. If we consider Si for the case Dlsup , the diffusion flux of Si induced by Mg segregation in the liquid,   DlSi Mg grads=l wlMg , is positive. It increases the positive

diffusion flux created by Si at the interface,  l DSi Si grads=l wlSi . The segregation coefficient being smaller than one, the composition difference at the interface is positive. It is also constant as the T is fixed. Moreover, the composition gradients have the same magnitude as composition evolutions in Fig. 6 are almost superimposed for Dlsup and Dldia . Thus, the velocity is increased in case Dl compared to case Dldia (with  sup  s=l l l DMg Si grad wSi ¼ 0), and recalescence takes place sooner Fig. 6(a). Now considering diffusion of Mg in case Dlinf , the composition gradient at the interface for Si is negative, giving rise to a  positive contribution of the flux

for Mg DlMg Si grads=l wlSi . It is decreased by the negative   diffusion flux created by Mg DlMg Mg grads=l wlMg . Thus,

the interface velocity is now controlled by the diffusion of Mg in the presence of a gradient of Si. This situation is more difficult to compare with case Dldia as the contribution   DlMg Mg grads=l wlMg

is now positive. When the full diffu-

l

sion matrix D is considered, the overall diffusion flux for each solute species are by the summation of a  controlled s=l l l positive flux, D for Si and grad w Si Si Si  

DlMg Si grads=l wlSi for Mg, and a negative flux,     s=l l l DSi Mg grad wMg for Si and DlMg Mg grads=l wlMg for Mg. Consequently, the presence of Si increases the total flux of Mg, while the presence of Mg decreases the total flux of Si. Compared to case Dlinf (same sign for each term with   DlSi Mg grads=l wlMg ¼ 0), case Dl leads to a lower velocity. Note that the effect of cross diffusion terms becomes very clear in Fig. 6. While it can be comprehensively explained for the present ternary system, it becomes more complicated for an alloy with more chemical species. Indeed the cumulative effect of all the solute gradients has to be considered in the analysis. This is possible with the present model.

binary alloy. So the same methodologies previously detailed can be used. The solute profile is expressed with Eq. (6) for any geometry. The kK i k values are computed with Eq. (10) and the key parameter k is still the solution of Eq. (15). The position of the planar front or cylinder radius is still expressed with Eq. (A2) as a linear function of the square root of time. The only differences are the corresponding F(.) and G(.) functions which are expressed as: x2  x  2e 4 x and GðxÞ ¼ pffiffiffi F ðxÞ ¼ erfc 2 px erfc 2 for planar growth ð30Þ  2 x F ðxÞ ¼ E1 4

x2 4e 4  2 and GðxÞ ¼ x x2 E 1 4 

for cylindrical growth

ð31Þ

It should be mentioned that these two expressions of the G(.) function, Eqs. (30) and (31), are similar to the ones encountered for the estimation of the supersaturation in growing condition for a 2D (i.e. parabolic) and 3D (i.e. paraboloid) dendrite tip [39,40]. Thus the G(x) function is equal to the inverse of the supersaturation expressed at a dendrite tip where x2 =4 is considered as the Peclet number pffiffiffiffiffi 2 [40], p.180. Indeed x2 =4, is equal to ðk= Bi Þ =4 (Eq. (16)) considering the computation associated to the eigenvalue i. However R, v and k are linked by the general relation, pffiffiffiffiffi 2 Eq. (21), and consequently ðk= Bi Þ =4 is equal to R v=ð2Bi Þ which is no more than the chemical Peclet number associated to index i. So the expression of the supersaturations at the surface of the two geometries, Eqs. (30) and (31), exhibits the same relation than the one encountered at the surface of a dendrite tip. Horvay and Cahn have previously obtained similar expressions of the supersaturation for a planar front (resp. cylindrical front) and parabolic dendrite (resp. paraboloid dendrite) in their mathematical development for the computation of dendritic and spheroidal growth [39]. Similarly, the methodology to compute the far-field composition using choice 3, i.e., Eqs. (26)–(29), can be repeated here. It leads to: N X N X U ij U 1 1 jk ! ðRf  RÞk j¼1 k¼1 k F pffiffiffiffiffi Bj pffiffiffiffiffi pffiffiffiffiffi  ½IðRf k; R Bj Þ  IðRk; R Bj ÞDwlk

ð32Þ

x2   2y x 2 Iðx; yÞ ¼ pffiffiffi e 4y  x erfc 2y p for planar growth

ð33Þ

~ li;1 ¼ w  li  w

4. Discussions



4.1. Planar and cylindrical evolutions The quasi-steady state analytical model can easily be extended to other geometries. Indeed Eqs. (1)–(5) are valid for the 1D planar (Cartesian) coordinates system as well as for the 1D cylindrical coordinates system, in which case the interface adopts a planar and a cylindrical morphology, respectively. As for the spherical coordinates system, the steady state solute profile for the multicomponent alloy is a linear combination of the solution corresponding to a

with

and

N X N X U ij U 1 1 jk ! ðR2f  R2 Þk2 j¼1 k¼1 k F pffiffiffiffiffi Bj pffiffiffiffiffi pffiffiffiffiffi  ½IðRf k; R Bj Þ  IðRk; R Bj ÞDwlk

~ li;1 ¼ w  li  w

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

ð34Þ

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

13

Fig. 7. Evolution of temperature, T , (a) solid fraction, gs , (b) and interfacial compositions, wl=s i , both for silicon (c) and magnesium (d) elements. Computations are developed for the Dl diffusion matrix for planar (dotted line), cylindrical (dashed line) and spherical (plain line) geometry with the present analytical [A] approach. Same cooling rate per unit mass has been applied for the three geometries. A magnification of the temperature evolution at the beginning of the cooling is proposed in (a).

with

2 3 x2  2 1 6 2  4y 2 x 7 Iðx; yÞ ¼ 48y e  2x2 E1 5 2 4y 2 for cylindrical growth

4.2. Large solute diffusion in the solid phase

ð35Þ

In order to visualize the effect of the coordinate system, simulations are proposed for the same full diffusion matrix Dl (Table 3) and the same alloy composition (Tables 1 and 4). The cooling condition per unit mass is given in Table 4. This value is the same for all geometries and corresponds to interfacial heat flux values uT = 2552 W m2 and uT = 1276 W m2 for planar and cylindrical geometry, respectively. This latter values are equal to UT qV =S where S and V are the area and volume of the respective geometries. Consequently the heats extracted per unit volume are the same in all cases. Fig. 7 shows the evolution of the temperature, fraction of solid and interfacial composition for the three geometries, so the results for the spherical geometry are the same as in Figs. 5 and 6. Large differences are observed between the curves. At the beginning of cooling, solid fraction is the highest for planar coordinates and the lowest for spherical coordinates. This is due to faster development of a planar solidification front. Temperatures are also higher in this later case due to a larger heat release and the corresponding cooling rate is lower. Moreover, no recalescence is observed for both planar and cylindrical coordinates. Similar types of differences are observed for the interfacial composition. As explained earlier, this is due to the fact that the cooling rate is directly linked to the growth velocity, itself dependent on the interfacial composition. If the dimension of the geometry decreases, a larger solid fraction is developed and the cooling rate is lowered.

The present analytical model has also been extended with large solute diffusion in the solid phase. With such hypothesis, the analytical model is modified with the assumption that the average composition in the solid phase  si , is equal to the interfacial composition, ws=l in Eq. (23), w i , for any chemical species. A simulation with large solute diffusion is proposed and compared to previous computations in Fig. 8. The alloy properties, liquid diffusion matrix Dl , and cooling conditions are unchanged (Tables 1 and 4). Only the spherical coordinates system is considered. Computations start at the liquidus temperature. It is stopped when the eutectic temperature, T E , is reached (Table 1). Note that this temperature is chosen arbitrarily here since the eutectic several successive secondary phases can form in multicomponent alloys, at distinct temperatures. Also a unique temperature does not characterize a given secondary phase as this is also dependent on diffusion in both the solid and liquid phases. So T E should be seen here as a selected temperature at which we compared the amount of remaining liquid phase left. The analytical model succeeds to model temperature and solid fraction evolutions when large diffusion in the solid phase is assumed. However, the difference with the simulation using no diffusion in the solid only becomes visible when reaching about 30 s. Before this threshold time, the two curves (Dl ; Ds ¼ 0 and Dl ; Ds ! 1) cannot be distinguished. This small difference is mainly explained by the small values of partition ratio, k Si and k Mg , (Table 1) leading to small deviations between the solidification paths. Higher solubility of alloying elements in the solid phase, higher compositions, lower eutectic temperature or larger nucleation

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

14

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

Fig. 8. Solidification paths with the analytical model [A] considering (Ds ¼ 0) no solute diffusion and (Ds ! 1) large solute diffusion in the solid phase showing (a) the temperature evolution as a function of time and (b) the solid fraction evolution as a function of the temperature. Analytical evolution considering Lever Rule (LR) and Gulliver–Scheil (GS) approaches is superimposed. Simulation with the analytical solution [A] is developed with the hypothesis of cross diffusion matrix in the liquid phase, Dl (Table 3).

undercooling lead to larger differences when comparing the remaining liquid at the temperature assumed for the formation of a second solid phase. The evolutions expected from the LR [1] and GS [2,3] solidification paths are also reported in Fig. 8. They are deduced considering the phase equilibrium computed with thermodynamical properties detailed in Table 1 and the balance equation of solute for each element with the hypothesis of no-diffusion in the solid phase: for the LR solidification path : wSi;0 wMg;0 þ mMg T ¼ T M þ mSi 1 þ ðk Si  1Þgs 1 þ ðk Mg  1Þgs

ð36Þ

for the GS solidification path: T ¼ T M þ mSi wSi;0 ð1  gs ÞkSi 1 þ mMg wMg;0 ð1  gs ÞkMg 1 ð37Þ where T M is the melting temperature of the alloys, ðmSi ; mMg Þ and k Si ; k Mg are the liquidus slope and partition coefficient for solute element i, respectively. Values are

given in Table 1. As expected, the curves computed with the analytical models retrieve the usual LR and GS approximations at low temperature, i.e. after about 40 s when solute gradients in the liquid phase have vanished. 4.3. Coupling with thermodynamic database The analytical model has been applied so far considering a linear multicomponent phase diagram. However, it is of prime importance to show that full coupling with thermodynamic databases is possible as multicomponent phase diagram can rarely by linearized. Such coupling has been previously demonstrated [14–16]. Hereafter, the Thermo-Calc software and the TCAL2 database are used to compute the interfacial compositions assuming thermodynamic equilibrium (Eq. (15)). Diffusion coefficients are extracted from the MOBAL2 database. However, as off-diagonal coefficients are not well defined in this database, they are still fixed to the same previous constant values listed in Table 3 for Dl . The on-diagonal coefficients are those proposed in this database. They correspond to the

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

15

Fig. 9. Solidification process with the analytical model considering a linear phase diagram and fixed diffusion coefficient in the liquid phase (dashed line) and coupling with thermodynamic database (plain line). Solidification path is shown both for temperature (a) and solid fraction (b) evolutions.

ones published by Du et al. [36]. Thus, the expression of the diffusion matrix becomes: 0

30000 B B 134e Rg T Dl ¼ B @ 

2

1 C C  109 m2 s1 71600 C A  99000e Rg T 3

ð38Þ The evolutions of the solidification process for the two approaches are presented in Fig. 9. For comparison with previous simulations, the simulation from Fig. 8 using no diffusion in the solid is reported. Differences between the curves are due to the decrease of the on-diagonal diffusion coefficients as well as non-constant values of the partition coefficients and liquidus slopes. This new simulation shows reduced solid fraction when reaching T E . Deviation between the curves in Fig. 9 is of the same magnitude as in Fig. 8. This means that the variations of both the diffusion coefficients and the properties of the phase diagram have effect of the same order of magnitude as diffusion in the solid phase and cannot be neglected. It should also be noted that the diffusion coefficients are also function of

the local composition. This means that they not only vary with time due to the temperature evolution but also in space. This effect is not accounted for in the present approach since the analytical solution is only valid with a constant diffusion matrix in space. The effect is yet expected to be smaller compared with the present analyses. The above comparisons also clearly identify the need for more systematic quantitative evaluation of the diffusion coefficients. 5. Conclusion A steady state analytical solution was proposed for isothermal 1D growth of a s=l interface in an undercooled melt controlled by multicomponent diffusion in the liquid phase. The interface velocity and the composition profile in the liquid phase were computed considering diffusion interactions between the chemical solute species. This analytical solution was validated by comparison with 1D numerical simulations developed with a finite difference implicit time stepping front tracking approach. The analytical solution being extended to cooling conditions, the integration of the composition profile was demonstrated as the

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

16

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

only relevant choice to estimate the composition at infinite position which defines the driving force for the growth of the interface. The extended analytical model was also compared and validated with the same numerical model for a set of diffusion matrices in the ternary Al–7 wt.%Si– 1 wt.%Mg ternary alloy. It was also demonstrated that the model could be extended and applied to planar and cylindrical geometries. This analytical model was also applied considering large diffusion in the solid phase. The analytical simulations retrieve the lever rule and Gulliver–Scheil approximations when the solute gradients in the liquid phase have vanished. Finally, full coupling with thermodynamic equilibrium calculations was also demonstrated. The present analytical model is currently limited to globular grain growth with/without diffusion in the solid phase. It has to be further developed for equiaxed dendritic growth. This should be conducted following previous efforts to model equiaxed dendritic solidification [7,8,11– 13,20–22] by using dendrite tip growth kinetics for multicomponent alloys [24]. Similarly, cross diffusion in the solid phase needs to be accounted with consideration of the full diffusion matrix as well as the possibility to account for successive secondary solid phases formed upon peritectic and/or eutectic reactions [16]. Relevance of such advanced microsegregation model is found in modeling of solidification structures and segregation in casting. But the present model could also find some application for modeling of solid-state precipitation [35,37,38]. It should find applications to various diffusion limited phase transformations in the field of physical metallurgy and materials science.

Appendix A. Globular solidification for a binary alloy For a binary alloy, Aaron et al. [25] have derived the solution for the solute diffusion field in a matrix phase ahead of an advancing spherical interface when its interfacial composition is fixed. However it should be pointed out that a similar analysis was previously proposed by Horvay and Cahn [39] for a large set of geometries. In Fig. 1, the situation corresponds to the existence of the first solute profile alone, with N = 1. The following additional considerations are defined to derive an analytical expression (see main text): no diffusion flux in the solid, semi-infinite medium with Rf ! 1, fixed far-field composition wl1 , fixed temperature T , constant value for the diffusion coefficient of the solute species in the liquid Dl . The composition in the liquid phase at time t and radial position r, wl ðrÞ, that satisfies Eqs. (1)–(4) with the above listed conditions is given by:   r ffi F pffiffiffiffiffiffi l wl ðrÞ  wl1 ¼ ðwl=s  wl1 Þ  D t ðA1aÞ k ffiffiffiffiffi p F Dl x2 pffiffiffi x e 4 p with F ðxÞ ¼  erfc ðA1bÞ x 2 2 and where the unknown growth parameter k is provided by the condition pffi ðA2Þ R¼k t 

Eq. (A2) was shown to be the only solution that permits to satisfy the set of Eqs. (1)–(4), and in particular the solute balance at the interface, Eq. (2), since the time derivative of Eq. (A2) gives an expression for the interface velocity. The solution can also be rewritten:   k wl=s  ws=l ¼ G pffiffiffiffiffil ðwl=s  wl1 Þ ðA3aÞ D with GðxÞ ¼ 

2 F 0 ðxÞ 2 1 ¼ 0 x F ðxÞ pffiffiffi x2   x C p B xe 4 erfc x2 @1  A 2 2

ðA3bÞ

The non-linear Eq. (A3) can be solved to extract the value of the k coefficient considering, for chosen values of temperature T and far-field composition wl1 , the interfacial composition wl=s and k s=l wl=s provided by Eq. (4). The interface position and velocity are then accessible by Eq. (A2) and its time derivative, respectively. Appendix B. Computation of the profile composition coefficients In spherical coordinates, for a composition evolving only as a function of the radial component, Eq. (1) writes: # "  N @wlj  @wli  1 @ 2X l with i 2 f1; N g ðB1Þ ¼ Dij r  @t r r2 @r @r  j¼1 t

Assuming a steady far-field composition wli;1 , introduction of the mathematical expression for wli ðrÞ proposed in Eq. (6) in Eq. (B1) leads to: " !# N X @ r K ij F pffiffiffiffiffiffi @t Bj t j¼1 " "  # # N N X 1 @ 2X r l @ r ¼ 2 Dij K jk F pffiffiffiffiffiffiffi r @r @r k¼1 Bk t j¼1

ðB2Þ

Further considering constant values for all coefficients of D and K, we can develop the expression as: ! N X r r 0 K ij pffiffiffiffiffiffiffiffi F pffiffiffiffiffiffi Bj t 2 Bj t 3 j¼1 " "  # # N N X 1 @ 2X 1 r l 0 ðB3Þ ¼ 2 Dij K jk pffiffiffiffiffiffiffi F pffiffiffiffiffiffiffi r r @r Bk t Bk t j¼1 k¼1 If we compute the derivatives, we can write: ! N X r r 0 K ij pffiffiffiffiffiffiffiffi F pffiffiffiffiffiffi Bj t 2 Bj t3 j¼1

    N N XX 1 4r r 2r2 00 r F pffiffiffiffiffiffiffi ¼ Dlij K jk 2 pffiffiffiffiffiffiffi F 0 pffiffiffiffiffiffiffi þ 2r Bk t Bk t Bk t Bk t j¼1 k¼1 ðB4Þ However the function F(x) defined in Eq. (A1b) is the solution of the equation 4F 0 ðxÞ þ 2xF 00 ðxÞ ¼ x2 F 0 ðxÞ. Consequently the previous equality can be simplified for

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030

G. Guillemot, C.-A. Gandin / Acta Materialia xxx (2015) xxx–xxx

pffiffiffiffiffiffiffi x equal to r= Bk t. Finally, if we exchange indices j and k, we obtain Eq. (7): "

N X

r

j¼1

2ðBj tÞ3=2

with

Bj K ij 

i 2 f1; N g

N X k¼1

# Dlik K kj

F

0

! r pffiffiffiffiffiffi ¼ 0 Bj t

[17] [18] [19] [20] [21]

ð7Þ

[22] [23]

References [1] W. Kurz, D. Fisher, Fundamentals of Solidification, fourth ed., Trans Tech Publications, 1998. [2] G. Gulliver, J. Inst. Met. 9 (1909) 120. [3] E. Scheil, Z. Metall. 34 (1942) 70. [4] H.D. Brody, M.C. Flemings, Trans. AIME 236 (1966) 615. [5] T.W. Clyne, W. Kurz, Metall. Trans. A 12 (1981) 965. [6] S. Kobayashi, Trans. Iron Steel Inst. Jpn. 28 (1988) 728. [7] M. Rappaz, Ph. The´voz, Acta Metall. 35 (1987) 2929. [8] C.Y. Wang, C. Beckermann, Metall. Mater. Trans. A 24 (1993) 2787. [9] H. Combeau, J.-M. Drezet, A. Mo, M. Rappaz, Metall. Mater. Trans. A 27 (1996) 2314. [10] U.R. Kattner, W.J. Boettinger, S.R. Coriell, Z. Metall./ Mater. Res. Adv. Technol. 87 (1996) 522. [11] M. Rappaz, W.J. Boettinger, Acta Mater. 47 (1999) 3205. [12] B. Appolaire, H. Combeau, G. Lesoult, Mater. Sci. Eng. A 487 (2008) 33. [13] M. Wu, J. Li, A. Ludwig, A. Kharicha, Comput. Mater. Sci. 79 (2013) 830. [14] H. Zhang, Ch.-A. Gandin, H. Ben Hamouda, D. Tourret, K. Nakajima, J. He, ISIJ Int. 50 (2010) 1859. [15] H. Zhang, Ch.-A. Gandin, K. Nakajima, J. He, IOP Conf. Ser.: Mater. Sci. Eng. 33 (2012) 012063. [16] D. Tourret, G. Reinhart, Ch.-A. Gandin, G.N. Iles, U. Dahlborg, M. Calvo-Dahlborg, C.M. Bao, Acta Mater. 59 (2011) 6658.

[24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40]

17

X. Dore´, H. Combeau, M. Rappaz, Acta Mater. 48 (2000) 3951. L. Thuinet, H. Combeau, Comput. Mater. Sci. 45 (2009) 294. Dictra, Thermo-Calc Software, Stockholm, Sweden. C.Y. Wang, C. Beckermann, Metall. Mater. Trans. A 27 (1996) 2754. H. Combeau, M. Zaloznik, S. Hans, P.E. Richy, Metall. Mater. Trans. B 40 (2009) 289. M. Wu, J. Li, A. Ludwig, A. Kharicha, Comput. Mater. Sci. 92 (2014) 267. J.S. Kirkaldy, D.J. Young, Diffusion in the Condensed State, The Institute of Metals, London, 1987. O. Hunziker, Acta Mater. 49 (2001) 4191. H.B. Aaron, D. Fainstein, G.R. Kotler, J. Appl. Phys. 11 (1970) 4404. H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Oxf. Sci. Publications, 1986. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes – The Art of Scientific Computing, third ed., Cambridge University Press, 2007. H.G. Landau, Q. Appl. Math. 8 (1950) 81. Ch.-A. Gandin, Acta Mater. 48 (2000) 2483. Ch.-A. Gandin, Y. Bre´chet, M. Rappaz, G. Canova, M. Ashby, H. Shercliff, Acta Mater. 50 (2002) 901. TCAL2, TCS Al-Alloys Database v2.1.1, Thermo-Calc Software AB, Stockholm, SE, 2013. MOBAL2, TCS Al-Alloys Mobility Database v2.0, ThermoCalc Software AB, Stockholm, SE, 2011. M.A. Martorano, C. Beckermann, Ch.-A. Gandin, Metall. Mater. Trans. A 34 (2003) 1657. M.A. Martorano, C. Beckermann, Ch.-A. Gandin, Metall. Mater. Trans. A 2004 (1915) 35. L. Rougier, A. Jacot, Ch.-A. Gandin, P. Di Napoli, P.-Y. The´ry, D. Ponsen, V. Jaquet, Acta Mater. 61 (2013) 6396. Y. Du, Y.A. Chang, B. Huang, W. Gong, Z. Jin, et al., Mater. Sci. Eng. A363 (2003) 140. Ch.-A. Gandin, A. Jacot, Acta Mater. 55 (2007) 2539. Q. Chen, B. Sundman, Mater. Trans. 43 (2002) 551. G. Horvay, J.M. Cahn, Acta Metall. 9 (1961) 695. J.A. Dantzig, M. Rappaz, Solidification, EPFL Press, 2009.

Please cite this article in press as: G. Guillemot, C.-A. Gandin, Acta Mater (2015), http://dx.doi.org/10.1016/j.actamat.2015.04.030