Phase-field lattice-Boltzmann study on eutectic growth with coupled heat and solute diffusion

International Journal of Heat and Mass Transfer 145 (2019) 118778

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Phase-field lattice-Boltzmann study on eutectic growth with coupled heat and solute diffusion Ang Zhang a, Fengyuan Liu a, Jinglian Du a, Zhipeng Guo a,⇑, Qigui Wang b, Shoumei Xiong a,c,* a

School of Materials Science and Engineering, Tsinghua University, Beijing 100084, China Materials Technology, GM Global Propulsion Systems, Pontiac, MI 48340-2920, USA c Key Laboratory for Advanced Materials Processing Technology, Ministry of Education, Tsinghua University, Beijing 100084, China b

a r t i c l e

i n f o

Article history: Received 3 December 2018 Received in revised form 16 July 2019 Accepted 22 September 2019 Available online 28 September 2019 Keywords: Eutectic growth Thermosolutal diffusion Phase-field simulation

a b s t r a c t The thermosolutal interaction influences eutectic evolution and thus properties of solidified materials. Limited by numerical capability, however, very few studies are performed on the eutectic growth with coupled heat and solute diffusion. Combining the phase-field lattice-Boltzmann approach and a parallel-adaptive mesh refinement algorithm, a novel numerical scheme is developed to efficiently simulate the thermosolutal multiphase eutectic evolution. The contact angles at the triple point agree well with the analytical solution, and the temperature always obtains the local extreme at the solid/liquid interface due to the release of latent heat. The effects of the Lewis number and imposed heat sink on the eutectic evolution are discussed in detail, which includes the change of the lamellar growth velocity and the form of lamellae creation. The dimension is further extended to 3D to investigate the evolution of plate- and rod-like eutectics. How to simulate eutectic evolution with a larger Lewis number (e.g., 106) and how to develop a more general eutectic model are also discussed. As the first attempt to solve the thermosolutal eutectic evolution, our investigation paves a way for further quantitative analysis and direct comparison with experiments. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Prediction and control of the formation of the eutectic morphology have long been a challenging subject in fundamental physics and physical metallurgy for decades [1–3]. When the temperature is below the eutectic point in the phase diagram, coexisting solid phases nucleate from the liquid phase simultaneously, forming a nearly periodic pattern with long-rang order in systems with low melting entropies, e.g., Al-Cu, Cu-Ag, and Pb-Sn alloys [4]. The pattern formation and the size distribution of the frozen-in two-phase structures are dependent on the interfacial thermodynamic equilibrium, i.e., the interaction between capillarity and thermosolutal diffusion during solidification [5]. However, performing experiments to investigate the micro thermosolutal diffusion is infeasible due to a rather large Lewis number in metals (i.e., Le = a/D ~ 104, the ratio of thermal to solutal diffusivity). As another alternative, numerical modelling, such as ab initio first principle calculations [6–9] and phase-field modeling ⇑ Corresponding authors at: School of Materials Science and Engineering, Tsinghua University, Beijing 100084, China (Z. Guo). E-mail addresses: [email protected] (Z. Guo), smxiong@ tsinghua.edu.cn (S. Xiong). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118778 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

(PFM) [10–12], is becoming increasingly popular. Still, most theoretical models on microstructure evolution are only based on pure mass transport. Neglecting thermal effect is generally accompanied by assuming that the thermal propagation is accomplished instantaneously. In fact, the pattern evolution is largely dependent on the coupled thermosolutal diffusion during solidification. Akamatsu et al. [5] observed that the interfacial thermal effect significantly altered the lamellar eutectic pattern under laser perturbations. Perrut et al. [13] found the curve of the eutectic isotherms would induce a dilatation of the growth pattern and further affected the growth stability. Loginova et al. [14] demonstrated that the nonisothermal effect would become significant for larger cooling rate. Thus, removing the isothermal assumption or so-called frozen temperature approximation [15] is essential to reveal the real thermosolutal interaction during eutectic solidification. To our best knowledge, however, there is no report on the coupled thermosolutal eutectic phase transition, and how the released latent heat and the thermal diffusion affect the eutectic evolution remains unknown. Rooted in the diffuse-interface model, the phase-field method (PFM) is becoming a standard computational tool to investigate the phase-change problems, particularly for problems with interface evolution [10–12,16]. By introducing order parameters

2

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

(i.e., the phase field variables) varying smoothly across the diffuse interface, explicit interface tracking is eliminated, and the topological change of the phase interfaces is handled without extra procedures. With its thermodynamic rigor, the PFM can accurately recover the Gibbs-Thomson effect and has been successfully employed to investigate the eutectic phase transition [17–20]. To couple the thermal diffusion, the nonlinear heat-conduction equation (see Eq. (4)) needs to be discretized. Limited by the numerical stability, the time step in the explicit finite difference method (FDM) is inversely proportional to the thermal diffusivity. In other words, the time step needs to be reduced by four orders of magnitude in comparison with that under pure solute diffusion, leading to a very long simulation time to approach the real case. To efficiently solve Eq. (4), Lan et al. [21] employed an adaptive finite volume method, Rosam et al. [22] used an adaptive, fully implicit multigrid algorithm, and the current authors [23] previously adopted an implicit parallel multigrid computing scheme. However, all those algorithm structures are rather complex, and the numerical pointcut is to solve the second-order heat-conduction equation by implicit iteration. As another excellent tool to simulate transport phenomena, the lattice Boltzmann model (LBM) can solve the thermal problem by establishing a dynamic model that incorporates the physics of microscopic transport process. The heat-conduction equation with second-order spatial derivatives is locally accessible in the latticeBoltzmann (LB) framework, and the phase-change term (e.g., the latent-heat source term) can be specified as a force term which is directly added to the evolution equation [24]. With its inherent parallel-processing structure and good stability in handling complex geometries, the LBM is becoming a popular alternative to the traditional continuum models [25]. The current authors adopted the LBM to solve the convective transport during dendritic [26,27] and eutectic [28,29] solidification rather than solving the Navier-Stokes equations. Yin et al. [30] employed the LBM to simulate dendrite growth in a thermo-solute-flow situation. Dorari et al. [31] proposed a multiple-grid-time-step LB technique to achieve a stable and high-efficiency solution on the dendritic evolution with a high Lewis number. Moreover, the simplicity of explicit formulations in the LBM offers enhanced efficiency, enforces the conservation law by improving the round-off precision [32], and avoids the solution of the expensive matrix problem during implicit time marching. In this work, a phase-field lattice-Boltzmann (PFLB) approach is developed to investigate the coupled thermosolutal eutectic phase transition. In detail, the PFM is used to simulate the eutectic growth, and the LBM is employed to determine the thermal evolution. Furthermore, a parallel-adaptive mesh refinement (ParaAMR) algorithm based on our recent work [15,33] is combined with the PFLB approach to improve the computing efficiency. Results show that the thermal distribution varies with the evolving solid/liquid (S/L) interface, and the temperature always obtains the local extreme at the S/L interface. The simulated eutectic evolution under thermosolutal interaction is discussed in detail, especially for the effects of the Lewis number and imposed heat sink. The lamellae creation under large undercooling is predicted, and the dimension is extended to 3D to investigate the evolution of plate- and rod-like eutectics. In the next section, the mathematical method, including the phase-field model, the lattice-Boltzmann model, and the paralleladaptive mesh refinement algorithm, is introduced. In Section 3, the numerical accuracy is evaluated in terms of the thermal evolution, the contact angles at the triple point, and the eutectic pattern under different grid levels, followed by the efficiency test. In Section 4, the typical lamellar eutectic simulation result is presented, and the influence of the Lewis number and imposed heat sink is discussed in detail. Then the dimension is extended to 3D to

investigate the evolution of plate- and rod-like eutectics, and how to simulate a larger Lewis number (e.g., 106) and how to develop a more general eutectic model are prospected. Concluding remarks are summarized in Section 5. 2. Mathematical method 2.1. Phase-field model In this work, the eutectic phase-field model developed by Kim et al. [18] is extended by coupling the thermal and solutal diffusion. Based on the interface-field method [34] and the postulate of equal chemical potentials among coexisting phases, the PFM is deduced from variational principle and has been used to investigate the eutectic evolution including the CBr4-C2Cl6 [15] and AlCu eutectic alloys [28,35]. The two solid phases (a and b) and the liquid phase (L) are denoted by three phase-field variables (i.e., /1, /2 and /3) varying smoothly between 0 and 1. The free energy over the domain V is expressed as an integral of the density functional which is split into the interfacial energy density gP and the chemical free energy density gT.

!# Z " X G¼ g P þ g T þ kL /i  1 dV i

V

¼

Z "X X j>i

V

þ

X



i

2 ij

e

2

!

r/i  r/j þ xij /i /j

/i g i ðC i ; T Þ þ kL

i

X

!# /i  1

dV

ð1Þ

i

with the governing equations: n @/i 2X dG dG si sj M ij  ¼ n j–i d/i d/j @t

! ð2Þ

3 X @C /i rC i ¼rD @t i¼1

ð3Þ

@T L @f s ¼ ar2 T þ  q_ @t cp @t

ð4Þ

where kL is a Lagrange multiplier to conserve the mass balance P between the phases i /i ¼ 1, si is a step function, i.e., si(x,t) = 0 if P /i = 0 and si(x,t) = 1 otherwise, and n ¼ i si ðx; tÞ. eij is the gradient energy coefficient, and xij is the height of the double well potential between the i and j phases (i, j = a, b and L, and i – j), both values of which are determined by considering an equilibrium state with a pffiffiffiffiffiffiffiffi flat phase boundary and they satisfy eij ¼ 4 frij =p, and xij ¼ 2rij =f, where rij and f are the interface energy and the half interface width respectively. Mij is the phase-field mobility coefficient which is determined by vanishing the kinetic coefficient P [18]. C ðx; t Þ ¼ i /i C i is the solute concentration of the coexisting phases, where Ci is the concentration of the i phase, fs = 1  /3 is the solid fraction, and D ¼ ð1  /3 ÞDs þ /3 Dl is the diffusivity which is dependent on the phase field, where Ds and Dl are the solute diffusivity in the solid and liquid phases respectively. a is the thermal diffusivity, L is the latent heat, and cp is the specific heat. q_ > 0 is an imposed heat sink to simulate the heat flux out of the computa_ the released latent heat would tional domain [23,36]. Without q, raise the domain temperature and arrest solidification. The functional derivative in Eq. (2) is given by:

dG X ¼ d/i j–i

e2ij 2

!

r /j þ xij /j þ g i  C i g C 2

ð5Þ

3

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

  RT E TE  T RT E with g i  C i g C ¼  CE  ð1  ki Þ  C mi tm tm i

ð6aÞ

ð10Þ

for i = 1 and 2 (i.e., a and b), and

g  C 3 gC ¼  3

RT E

tm

ð6bÞ

C3

where R is the molar gas constant, tm is the molar volume, CE and TE are the eutectic concentration and eutectic temperature respectively. mi = (TE  T)/(CE  C3i) is the liquidus slope, and ki = Ci/C3 is the partition coefficient of the i (i = a and b) phase, where C3i is the liquid concentration in equilibrium with the i phase at the temperature T. It is noted that Eqs. (6a) and (6b) indicate the amount of the supersaturation from the equilibrium concentration. The phase-field equations are discretized using the explicit finite difference method on a Cartesian uniform grid in space and a first-order forward Euler stepping method in time [15]. A ninepoint stencil scheme is employed for the Laplace operator in Eq. (5), while a net flux control volume scheme is used for the divergence operator in Eq. (3). To relax the restriction on the time step and avoid implicit iteration, the heat-conduction equation Eq. (4) is solved using the following lattice Boltzmann equations. 2.2. Lattice Boltzmann model As a mesoscopic kinetic method, the LBM describes the transport process by repeated streaming and collision of a collection of factitious particles, and it keeps consistent with the heatconduction equation (i.e., Eq. (4)) by the Chapman-Enskog expansion [25]. The starting point is the kinetic description of the particle distribution function.

  ! @f i r ; t LBM @tLBM

  ! ! þ c i  rf i r ; tLBM ¼ Xi

ð7Þ

  ! where the particle distribution function f i r ; t LBM represents the ! probability of finding a given particle with velocity c i at a given ! ! position r and time t LBM , i.e., the density of particles with c i at ! r and t LBM [25]. For the employed two-dimensional nine-velocity ! ! (D2Q9) model, c 0 ¼ ð0; 0Þ, c i ¼ cðcoshi ; sinhi Þ for i = 1–4, and pffiffiffi   ! c j ¼ 2c coshj ; sinhj for j = 5–8, where hi = p(i  1)/2, hj = p (2j  1)/4, and c = dx/dt is the lattice velocity, where dx and dt are the lattice spacing and the time step in the LBM respectively, and both are rescaled to 1. Oi represents the collision operator and can be simplified by the Bhatnagar-Gross-Krook (BGK) approximation [37].

h 





eq ! Xi ¼  f i ! r ; tLBM  f i r ; tLBM

i

=sLBM

ð8Þ

where sLBM is the single relaxation time toward the equilibrium distribution and is related to the macroscopic thermal diffusivity a by   ! a ¼ c2 dtð2sLBM  1Þ=6. f eq r ; t LBM ¼ wi T is the equilibrium distrii bution function, where T is the macroscopic temperature, and wi is the weight coefficient which is dependent on the discrete velocity model, i.e., w0 = 4/9, w1-4 = 1/9 and w5-8 = 1/36 for the D2Q9 model. The kinetic equation Eq. (7) with the BGK approximation Eq. (8) is discretized as [24]:

 fi

! r ;tLBM þdt

   ! f i r ;t LBM dt

¼

    ! ! ! ! f i r þ c i dt;tLBM þdt f i r ;tLBM þdt þ ci ! c dt     i ! ! f i r ;t LBM f eq r ;t LBM i sLBM

  ! After canceling out f i r ; tLBM þ dt , Eq. (9) is simplified as

         eq ! ! ! ! ! f i r þ c i dt; t LBM þ dt ¼ f i r ; t LBM  f i r ; t LBM  f i r ; t LBM =sLBM

ð9Þ

To incorporate the phase change and the imposed heat sink in Eq. (4), a source term is added to the right side of Eq. (10) [24,38], i.e.,

  ! ! f i r þ c i dt; tLBM þ dt        eq ! ! ! ¼ f i r ; tLBM  f i r ; t LBM  f i r ; t LBM =sLBM þ F i dt ð11aÞ  with F i ¼ 1 

   L @/3 wi   q_ LBM 2sLBM cp @t LBM 1

ð11bÞ

Accordingly, the temperature is expressed as follows [38].

T¼

X

f þ i i

  dt L @/3   q_ LBM 2 cp @t LBM

ð12Þ

It is noted that there is a dimensional conversion between the PFM and the LBM, i.e., t ! tLBM and x ! xLBM for the time and length respectively. In detail, the time and length are scaled by the time step dt and the maximum grid size dxmax (i.e., the grid size at the bottom level) during simulation, where dt is limited by the constraint stability of the explicit algorithm, i.e.,

dt ¼

rdx2min 2Nd D

ð13Þ

where r is an introduced stability coefficient which is less than 1, dxmin is the minimum grid size (i.e., the grid size at the top level), and Nd is the domain dimension (i.e., 2 and 3 for 2D and 3D cases respectively). It is noted that the length scaling is dxmax rather than dxmin because the quantities in the multilevel LBM must be advanced from the coarsest to the finest grid level, i.e., from ‘‘bottom to top” [39]. To match the macroscopic thermal diffusivity, the relaxation time in the LBM should satisfy





sLBM ¼ 3aLBM = c2 dt þ 1=2

ð14Þ

where aLBM ¼ adt=dx2max . Note that dxmax = 2n1 dxmin, and by combining Eqs. (13) and (14), the following is obtained





sLBM ¼ 3rLe= 22n1 Nd þ 1=2

ð15Þ

where n is the number of the grid level in the hierarchical data structure with a refining spacing ratio of 2 (i.e., the length ratio of mesh size between neighboring levels). It is noted that if Eq. (4) is solved using the explicit difference method, the time step will satisfy dta ¼ rdx2min =ð2N d aÞ ¼ ðr=LeÞdx2min =ð2N d DÞ ¼ r a dt, where ra = r/Le. Taking Le = 104 for instance, r ~ 1, ra ~ 104, and thus dta = 104 dt in the one-level structure, i.e., the time step is reduced by four orders of magnitude compared with dt in Eq. (13). On the other hand, based on Eqs. (13)–(15), the time step in the multilevel data structure is still inversely proportional to the solutal diffusivity rather than thermal diffusivity. But this does not mean that the time step is enhanced by four orders of magnitude compared with dta. In the multilevel LBM structure, the stability coefficient r must be properly predetermined, because r is related to the relaxation time sLBM which should be close to unity if possible [25]. A larger sLBM will cause larger numerical error, while employing a sLBM close to 0.5 might lead to numerical instability in the BGK scheme. Accordingly, in 2D, one can choose Nd = 2, n = 5, Le = 104, sLBM = 1.08, and thus r = 0.02, i.e., the time step is 200 times larger than dta. Moreover, the complicate programming on solving the nonlinear equations during

4

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

implicit difference method is avoided, which enforces the conservation laws by improving the round-off precision [32]. 2.3. Coupling AMR with PFLB approach Besides less restriction on the time step, the multilevel data structure is an extension of the concept invoking nonuniform grids to improve the computing efficiency. Grids with higher resolution is preferred at the phase interfaces (including the a/b, a/L, and b/L interfaces) where a more accurate solution is desired. In the following, we use the ‘phase interfaces’ to denote all three interfaces including the a/b, a/L, and b/L interfaces. The prerequisite to build the multilevel data structure is to establish a refinement criterion, i.e., how to find and tag the cells that need refined/coarsened. In this work, a block-structure adaptive-mesh-refinement (AMR) technology based on our recent work [15,33,40] is extended to efficiently solve the thermosolutal multiphase evolution. The refinement criterion is performed according to a gradient condition, i.e.,

max ðjr/i jÞ þ EC jrC j þ ET jrT j P n 16i63

ð16Þ

where EC and ET are the weight coefficients for the solute and temperature respectively, and n is a threshold value determined via numerical tests. After tagging, the grids on each level are subdivided with a predefined refining ratio (i.e., 2), and the refinement process continues until the top (i.e., the finest) grid level is reached. The tagged grids are then grouped into the so-called patchboxes by a cluster algorithm developed by Berger and Rigoutsos [41]. Fig. 1 (a) and 1(b) show the typical hierarchical mesh structure superimposed on the eutectic solute field. The patchboxes are depicted using the coarser solid lines, in which the meshes with same size are filled. The whole domain is first covered by the coarsest meshes, and once Eq. (16) is reached, the cells are tagged, subdivided and clustered sequentially. Accordingly, a multilevel data structure is built, and the nearer to the phase interfaces, the smaller the grid size. The structure of the patchboxes on each grid level must be correctly nested, i.e., the finer grids must be located inside the coarser ones. As shown in Fig. 1(b), all grids on the level n are always inside those on the level n  1, and the phase interfaces can be perfectly characterized by the array of the patchboxes on the top level (see level IV in Fig. 1(b)). To accomplish the data communication among the patchboxes, a layer of ghost cells is added at the external boundaries of each patchbox to collect the data from its closest neighbors, which happens either within one processor or between neighboring processors. For the data on neighboring grid levels, the restriction and interpolation operations are performed to update all variables [33,42].

In the multilevel LBM, a multistep advancing scheme needs to be executed to maintain the lattice velocity to be 1 on any grid level. In detail, after finishing one ‘‘streaming-collision” step on the coarser grid level, 2lvd ‘‘streaming-collision” steps is proceeded on the subsequent finer grid level [43], where lvd is the level difference between the chosen coarser and finer grid levels. Besides, the relaxation time needs to change with the grid level to keep the thermal diffusivity constant [44], i.e.,

sLBMf ¼ 2  ðsLBMc  1=2Þ þ 1=2

ð17Þ

where sLBMf and sLBMc are the relaxation time at the fine and coarse grid levels, respectively. To ensure the continuity of the physical quantities over the coarse/fine interface, the distribution function needs to be modified by summing the equilibrium and nonequilibrium parts [44]. f

     _c ! eq;f ! ! r ; t LBM ¼ a f i r ; t LBM þ ð1  aÞf i r ; tLBM

ð18aÞ

c

     f  eq;c ! ! ! r ; t LBM ¼ 1=a  f i r ; t LBM þ ð1  1=aÞf i r ; t LBM

ð18bÞ

fi fi

where a = 0.5 sLBMf /sLBMc is a weight coefficient, and the superscripts f and c denote the fine and coarse grid levels, respectively. _  f i and f i denotes the spatially and temporally interpolated and restricted distribution function respectively. 2.4. Parallel computation After the adaptive mesh refinement, the computing data is broadcast to all processors, based on which the computing work is constructed, partitioned, and dispatched on a distributed memory level controlled by a message passing library named Message Passing Interface. Consequently, each process has its own but different data which is dependent on the layout of patchboxes, achieving so-called single program multiple data model. Based on a space-filling curve approach, the load assignment of eight processors is depicted by different colors in Fig. 1(c). Corresponding to Fig. 1(a), each patchbox is dispatched to an individual processor, and the configuration of processors is consistent with the layout of patchboxes. 3. Numerical tests The Al-Cu eutectic alloy is simulated using thermo-physical parameters shown in Table 1, which are assumed constant during simulation. Four eutectic couples with the height of 3 lm are initialized at the bottom of the square domain (51.2  51.2 lm2) filled with the undercooled melt (1 K). Unless stated otherwise,

Fig. 1. (a) Typical hierarchical mesh structure superimposed on the eutectic solute field. (b) Mesh structure on each level. (c) Parallelization of the computing work to eight processors which are denoted by different colors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778 Table 1 Thermo-physical parameters of the Al-Cu eutectic alloy. Parameters

Al-Cu 2

Dl (solute diffusivity in liquid, m /s) Ds (solute diffusivity in solid, m2/s) TE (eutectic temperature, K) CE (eutectic concentration, mol frac) ma (liquidus slope of a phase at TE, K/mol frac) mb (liquidus slope of b phase at TE, K/mol frac) ka (partition coefficient in a phase to L phase) kb (partition coefficient in b phase to L phase) raL(a/L interface energy, J/m2) rbL (b/L interface energy, J/m2) rab (a/b interface energy, J/m2) L (latent heat, J/kg) cp (specific heat, J/(kg K) a

3  109 [45] 3  1013 [45] 821.4 [46,47] 0.173 [46,47] 1050 [46,47] 488 [46,47] 0.14 [45] 1.85 [45] 160.01  103 [46] 88.363  103 [46] 219.484  103 [46] 3.814  105a 4.827  104a

slope of the curve becomes flatter over time, i.e., the domain temperature tends to become uniform. The simulated values using the LBM agree well with the analytical solutions, indicating that the LBM can accurately reproduce the thermal evolution. Fig. 2(b) shows the temperature distribution along the domain diagonal under different grid levels, in which the inserted cloud picture shows the typical four-level grid architecture based on the gradient refinement criterion. The curves overlap, indicating that the multilevel grid configuration does not compromise the accuracy. The stair-step curves can be attributed to the mesh size of the discretized cells, which becomes more significant when a larger-size mesh is spanned. 3.2. Contact angle

TM

The parameter comes from AnyCasting software.

the following parameters are adopted: the heat sink is 2.5 K/s, the initial lamellar spacing is 12.8 lm, the minimum mesh size is 0.2 lm, the number of the grid level is 4, and the maximum mesh size is 1.6 lm. For all variables including the phase field, solute and temperature, a periodic boundary condition is applied at both left and right sides, while a zero-Neumann boundary condition is set at both top and bottom sides. 3.1. Thermal evolution To test the accuracy of the temperature field calculated by the LBM, the conduction-induced heat transfer in a semi-infinite slab is simulated as a benchmark case, in which the solution has been analytically determined [48] as follows.

 pffiffiffiffiffiffiffiffi T ¼ A þ B  erf x= 2at

5

ð19Þ

where A and B are constants dependent on the boundary conditions, and erf (y) is the error function. To make comparison, Eq. (4) is numerically solved using the LBM by setting a = 1, L = 0 and q_ ¼ 0. As shown in the inset in Fig. 2(a), a periodic boundary condition is set at both top and bottom sides, and an adiabatic boundary condition is applied at the right side to satisfy the semi-infinity assumption. A constant temperature TM = 1 is imposed at the left side, and T0 = 0 is initialized elsewhere. Accordingly, the analytical solution of Eq. (19) can be determined by setting x = 0, T = 1 and x = 1, T = 0, and Eq. (19)  pffiffiffiffiffiffiffiffi becomes T ¼ 1  erf x= 2at . Fig. 2(a) shows the comparison of temperature distribution along the x axis at different times. The

The eutectic phase-field model has been exhaustedly studied by direct comparisons with experiments [15,18]. To validate the simulated eutectic pattern, the contact angles are measured at the triple point where the three binary interfaces (the a/b, a/L, and b/L interfaces) are connected by imposing local equilibrium. According to Young’s law, the surface tension at the steady state satisfies (see Fig. 3(b)).

raL cosha ¼ rbL coshb

ð20aÞ

raL sinha þ rbL sinhb ¼ rab

ð20bÞ

where ha (hb) is the contact angle of the a/L (b/L) interface with the x direction and is assumed isotropic during eutectic evolution. The analytical solution based on the values in Table 1 is ha = 69.92° and hb = 51.55°. Fig. 3(a) shows the evolution of interface profiles according to /1 = 0.5 and /2 = 0.5. The interface reaches a steady state after a transition period. The triple point then moves vertically, i.e., the trajectory of the a/b interface keeps vertical, and the contact angles at the steady state remain unchanged. The measured ha and hb in Fig. 3(b) are 68.36° and 52.40°, and the relative deviation is 2.2% and 1.6% comparing with the analytical solution respectively, indicating the accuracy of the numerical results. It is noted that the anisotropy of surface energy is neglected and thus the contact angles can be determined by the force balance, i.e., Eqs. (20a) and (20b). 3.3. Eutectic evolution The thermosolutal interaction determines the eutectic evolution. Considering that the thermal diffusivity is normally four orders of

Fig. 2. Temperature distribution in the conduction-induced heat transfer of a semi-infinite slab. (a) Comparison between the LBM and analytical solution for the temperature along the x axis at different times. The inset shows the initial configuration including the boundary conditions and domain size. (b) Temperature distribution along the domain diagonal under different grid levels. The insert shows a typical four-level grid structure superimposed on the cloud picture of the temperature field.

6

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

Fig. 3. (a) Evolution of interface profiles according to /1 = 0.5 and /2 = 0.5. (b) Local enlarged image of the region circled by the blue line in (a). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

magnitude larger than the solutal diffusivity for metals, the domain temperature will become uniform instantly. For that, we fix the Lewis number at 10 in the present test cases, which weakens the thermal transport capacity and reduces the computing overhead. Fig. 4 shows the comparison of the solute field and temperature field at 1.6 s (i.e., after 600,000 steps) under different grid levels, i.e., 1, 2, 3, and 4 from left to right in column order respectively. In the cloud pictures of the solute field (see the first row), the a, b, and liquid phases are colored using blue, red, and cyan respectively. As the number of the grid level increases, the grid structure adjusts its layout to characterize the phase interfaces. The average eutectic growth velocities of the four cases are all 15.6 lm/s. In the cloud pictures of the temperature field (see the second row), a local high-temperature region exhibits at the S/L interface, and the temperature ranges under different grid levels remain same.

Fig. 5 shows the comparison of the interface profiles according to /1 = 0.5, /2 = 0.5, and /3 = 0.5 between the results with one grid level (depicted by the dashed red line) and those with four gird levels (depicted by the solid blue line). The interface profiles in the two cases overlap, and the largest difference is 0.1% in the middle of the a/L interface. The inset shows the multilevel structure with four grid levels which is superimposed on the interface profiles. The interface region including three binary interfaces and the triple region is covered by the finest grids. The coarser grids are used to cover the region far from the phase interfaces, which ensures the numerical accuracy. Fig. 6 shows the comparison of the distribution of the solute concentration in both solid phase (a) and liquid phase (b) along the direction parallel to the S/L interface, as designated by the white arrows in Fig. 4(a2). The distance away from the S/L interface

Fig. 4. Comparison of the solute field (the first row, a1–a4) and temperature field (the second row, b1–b4) under different grid levels, i.e., 1, 2, 3, and 4 from left to right in column order respectively. The a, b, and liquid phases are colored using blue, red, and cyan respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

7

3.4. Efficiency In Fig. 4, eight processors are called to compute the thermosolutal eutectic evolution, and the corresponding elapsed times are 3882, 2773, 1941, and 1705s from left to right, respectively. The mesh adaptivity shortens the elapsed time by reducing the number of the involving grids. However, it is noted that the mesh generation process also occupies the computing resources, which becomes more significant if the number of the grid level is too large for the limited computing data [33]. Thus the elapsed time of the four-level architecture does not exhibit remarkable improvement comparing with that of the three-level architecture. The parallel computation also significantly enhances the computing efficiency. When the number of the parallel processors changes from 8 to 48, the total elapsed time is shortened to 383 s from 1705 s for the four-level data structure. Comparing with 23436 s by using the common method (i.e., both the number of the parallel processors and the number of the grid levels are 1), the combination of parallel computing and the adaptive mesh refinement can achieve the most attracting efficiency improvement. In detail, the efficiency is accelerated by a factor of 61 for 48 processors and four-level structure in comparison with the common method. Fig. 5. Comparison of the interface profiles according to /1 = 0.5, /2 = 0.5, and /3 = 0.5 between the results with one grid level (depicted by the dashed red line) and those with four gird levels (depicted by the solid blue line). The inset shows the four-level mesh structure superimposed on the interface profiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

is 0.5 k (k is the lamellar spacing, see Fig. 4(a3)). The concentrationdistance curves under different grid levels overlap with one another except very small deviation at the concentration extremes in the liquid (see Fig. 6(b)), further indicating that the adaptive nonunifom meshes do not compromise numerical accuracy. The solute concentration remains constant in the solid phase, while it exhibits sinusoidal-like distribution in the liquid and obtains the local extremes on the top of the midpoint of the a/L and b/L interfaces [2]. The difference of the solute concentration corresponds to different constitutional undercooling. The eutectic lamellae grow into the undercooled melt horizontally, indicating the nearly equal total undercooling between the a- and b-lamellae. Accordingly, the curvature undercooling needs to be adjusted, which causes the nonplanar S/L interface. Besides, the temperature along the two white lines in Fig. 4(a2) is similar to one another among the four kinds of grid levels, i.e., remaining constant as 819.9937 K and 819.9908 K in the liquid and solid phases respectively.

4. Results and discussion 4.1. Typical eutectic growth Fig. 7 shows the evolution of the solute field and temperature field with four-level grid structure during eutectic growth. Each column corresponds to the simulation time 0, 0.53, 1.33, and 2.13 s respectively. The coexisting solid phases grow into the undercooled melt with dynamically adjusted grid structure. Because the field variables change the fastest at the phase interfaces, the finest meshes are generated to cover the phase interfaces. The eutectic couples grow in a nearly planar S/L interface, leading to a significant direction-dependent temperature distribution. In detail, the temperature field exhibits little difference along the x axis, but forms a nonuniform distribution along the growth direction (i.e., y axis). Besides, the combined effect of heat sink and released latent heat during eutectic growth causes the fluctuation of the temperature scale. Fig. 8 shows the instantaneous growth velocity and domain temperatures versus time, where T0.25, T0.5, Tmax, Tmin are the temperature at the points (0.5 L, 0.25 L), (0.5 L, 0.5 L), the minimum temperature, and the maximum temperature, respectively, and L is the domain height and width. As the eutectic lamellae grow, the domain temperature first decreases due to the dominated

Fig. 6. Comparison of the distribution of the solute concentration in the solid phase (a) and liquid phase (b) along the direction parallel to the S/L interface under different grid levels. The direction is along the white arrows in Fig. 4(a2), and the distance away from the S/L interface is 0.5 k (k is the lamellar spacing).

8

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

Fig. 7. Evolution of the solute field (the first row, a1–a4) and temperature field (the second row, b1–b4) with four-level grid structure during eutectic growth. Each column corresponds to the simulation time 0, 0.53, 1.33, and 2.13 s respectively.

Fig. 8. Instantaneous growth velocity and domain temperatures versus time, where T0.25, T0.5, Tmax, Tmin are the temperature at the points (0.5 L, 0.25 L), (0.5 L, 0.5 L), the minimum temperature, and the maximum temperature, respectively, and L is the domain height and width. A and B are two intersections, and the inset shows the local zoom near the intersection A.

effect of heat sink in the initial stage, which corresponds to the increase of the growth velocity. After experiencing 1.368 s (see the vertical dashed line), the accumulated heat due to the release of latent heat becomes comparable with the heat sink and raises the domain temperature slightly. Since then, both domain temperature and growth velocity have changed in a narrow range. When the S/L interface reaches (0.5 L, 0.25 L), T0.25 is equal to Tmax, as shown by the intersection A in the inset, indicating that the maximum temperature is obtained at the S/L interface, which can be attributed to that the S/L interface is the position of the heat source due to the release of latent heat. Similar behavior can also be deduced from the intersection B for T0.5. Fig. 9 shows the distribution of the phase field /3, solute concentration, and temperature along the centerlines of both the a and b phases at 1.6 s. Both the phase field variable and solute concentration remain unchanged in the bulk phases but exhibit rapid change at the S/L interface, where /3 changes smoothly from 0 in

the solid to 1 in the liquid. Because of the release of latent heat, the temperature first increases and then decreases and obtains the extreme at the S/L interface. From the local enlarged view in Fig. 9(b), the changes of both /3 and temperature along the centerline of the a phase lag behind those along the centerline of the b phase, which can be attributed to the different interface profiles between the a- and b-lamellae (see Fig. 5). Because a symmetric model alloy with the same physical parameters has the same interface structure [17], the asymmetric interface profiles between the a/L and b/L interfaces in the Al-Cu alloy is considered to be induced by the difference of both surface energy and atoms mobility of eutectic elements. Furthermore, the transition layer of the temperature field is remarkably larger than that of the solute field. It is expected that if the Lewis number becomes larger, e.g., 104, the thickness of the thermal boundary layer will increase accordingly, and the eutectic growth will be mainly controlled by the solute transport.

4.2. Influence of growth parameters The thermal distribution during thermosolutal eutectic growth can be changed by adjusting the heat sink and the Lewis number. Fig. 10 shows the extreme temperatures and the average growth velocity of the a phase versus the heat sink and the Lewis number. The default values of the heat sink and the Lewis number are 2.5 K/ s and 10 respectively. To make better comparison, the simulation time is fixed at 1.6 s in all cases. The average growth velocity is computed in a time range from the initialization to 1.6 s. As shown in Fig. 10(a), the increase of the heat sink lowers the domain temperature but increases the driving force for solidification, leading to a faster growth velocity. When the value of heat sink reaches 3.0 K/s, the lamellar spacing is halved by creating new lamellae at the center of the a phase (see the insets in which the values of the heat sink are 2.0 and 4.0 K/s, respectively). Similar phenomenon has been observed in eutectic experiments reported by Dantzig and Rappaz [2], in which the lamellae creation occurs when the growth velocity is increased. To accommodate the extra

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

9

Fig. 9. Field variables including the phase field /3, solute, and temperature versus the y axis at 1.6 s. The change paths, starting from domain bottom to domain top, are along the centerlines of both the a and b phases, and (b) is enlarged view of the shaded area in (a).

Fig. 10. Extreme temperatures and average growth velocity versus (a) the heat sink and (b) the Lewis number, respectively.

solute at the center of the wider lamellae (i.e., the a phase designated by blue in the insets), a negative curvature is developed [2], which increases the S/L interface area and leads to a temperature rise due to more released latent heat. After that, the domain temperature continues to decrease linearly, but the extreme difference of temperature (i.e., Tmax  Tmin) is enlarged, indicating the domain temperature becomes increasingly nonuniform. It is noted that the presence of heat sink mimics the cooling of the whole system and has little influence on the thermal diffusivity. Thus the change of the domain temperature exhibits linear change in Fig. 10(a). When the heat sink is fixed, both the domain temperatures and the average growth velocity change with log10Le in a nonlinear way, and further plotting indicates that such variation is also nonlinear in both temperature-Le and growth velocity-Le relations. A larger Le corresponds to faster heat propagation, i.e., tending to uniform the domain temperature, which decreases the thermal accumulation near the S/L interface and extreme difference of temperature. After the lamellae creation occurs, the eutectic evolution reaches a new thermodynamic state. The interaction among the thermal diffusion, the movement of the heat source (i.e., the advancement of the S/L interface), and the heat sink determines the thermal distribution and eutectic evolution. The increase of either heat sink or Lewis number can induce creation or generation of new lamellae which makes the lamellar spacing halved. The eutectic growth is controlled by the interaction between capillarity and diffusion. The increase of either heat sink or Lewis number will increase the growth velocity, leading to a situation where the steady-state lamellar spacing becomes larger than the extreme spacing based on the Jackson-Hunt eutectic theory [3]. Accordingly, the eutectic system is located to a large degree

on the diffusion branch in the undercooling-lamellar spacing profile proposed by Dantzig and Rappaz [2], i.e., dominated by the diffusion effect. More eutectic elements are rejected into the liquid phase, causing the build-up of solute ahead of the a/L inter-

Fig. 11. Evolution of the a/L interface profile according to /1 = 0.5 during the form of new b-lamellae by nucleation in the middle. The lamellar spacing knew = kold/2, and the S/L interface at the last snapshot is depicted by the red line according to /3 = 0.5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

10

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

Fig. 12. 3D plate-like eutectic growth including the solute field and temperature field. (a–c) are the solute field at different simulation time 0.36, 1.42, and 2.13 s respectively. (d) is the temperature field at 2.13 s.

able due to more rejected eutectic elements, and the lamellae creation becomes earlier triggered. The heat sink is introduced to avoid the continuous temperature rise induced by the release of latent heat. To simulate steady-state eutectic evolution, the combination of the heat sink and Lewis number should be kept in an appropriate range, otherwise the unexpected undercooling (too large or too small) will deteriorate the present simulation. The employed phase-field model is established on a small-undercooling assumption, i.e., approaching the equilibrium state. The near-equilibrium assumption means that the solute concentration in the liquid is close to the eutectic point in the phase diagram, and the eutectic growth can fast converge to the equilibrium state [18]. Too large undercooling will deviate this prerequisite, while too small undercooling will cause tedious computing consumption. Fig. 13. Distribution of the solute concentration along the x direction corresponding to Fig. 12(b). The distance away from the S/L interface is 0.5 k, and the position is at z = 0.125 Z, 0.25 Z, and 0.375 Z respectively (see the white arrow in Fig. 12(b)).

face and solvent on the b/L interface. The interface of the wider lamellae (i.e., a) depresses to accommodate the extra elements, i.e., forming a pocket with negative curvature. Fig. 11 shows the evolution of the a/L interface profile according to /1 = 0.5 during the form of new b-lamellae by nucleation in the middle, in which the lamellar spacing knew = kold/2, and the S/L interface at the last snapshot is depicted by the red line according to /3 = 0.5. Le = 50, and the heat sink is 2.5 K/s. As the eutectic grows, the interface depression is aggravated until the other eutectic solid phase (i.e., b) forms. A larger heat sink means a greater undercooling, and a larger Le corresponds to a faster thermal diffusion (i.e., less heat accumulation near the interface), both of which increase the undercooling and the driving force for solidification. As the heat sink or Le increases, the interface depression becomes increasingly remark-

4.3. 3D eutectic growth To generalize the thermosolutal eutectic simulation, the present study is extended to 3D. A so-called three-dimensional nineteenvelocity (D3Q19) discrete model is employed to handle the distribution function in the LBM [25], and the corresponding weight coefficients along different directions are w0 = 1/3, w1-6 = 1/18 and w7-18 = 1/36. Fig. 12 shows the simulated 3D plate-like eutectic including the solute field and temperature field. The thickness along the z axis Z = 12.8 lm = 0.25 L, q_ = 2.5 K/s, Le = 10, and the simulation time is 0.36, 1.42, 2.13, and 2.13 s for (a–d), respectively. A periodic boundary condition is applied at new added sides in 3D, and the boundary conditions of the other sides keep the same as those in 2D. Four eutectic couples are initialized at the domain bottom. The eutectic couples grow stably with a nearly planar S/L interface, and the local high-temperature region exhibits near the S/L interface due to the release of latent heat. As a more appealing progress, the 3D thermosolutal simulation makes it possible to compare with real solidification of bulk eutectics.

Fig. 14. 3D rod-like eutectic growth including solute and temperature. (a) is the initialization, (b–c) are the solute field at 1.42 s and 2.13 s, and (d) is the temperature field at 2.13 s, respectively.

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

Fig. 15. Distribution of the solute concentration along the x direction corresponding to Fig. 14(b). The distance away from the S/L interface is 0.5 k, and the position is at z = 0.125 Z, 0.25 Z, and 0.375 Z respectively (see the white arrows in Fig. 14(b)).

Fig. 13 shows the distribution of the solute concentration along the x direction corresponding to Fig. 12(b). The distance away from the S/L interface is 0.5 k, and the position is at z = 0.125 Z, 0.25 Z, and 0.375 Z respectively (see the white arrow in Fig. 12(b)). Sinusoidal wave distribution is observed, and the three concentrationdistance curves overlap, i.e., the concentration distribution along the x direction does not vary with the z coordinate in the platelike eutectic. This indicates that the plate-like eutectic pattern can be considered as an extension of the 2D lamellar eutectic along one coordinate direction (i.e., the z axis). Fig. 14 shows the 3D rod-like eutectic growth, in which Fig. 14 (a) shows the initiation, and Fig. 14(b–d) show the solute field at 1.42 s and 2.13 s and temperature field at 2.13 s, respectively. The thickness along the z axis is Z = 25.6 lm = 0.5 L, and other parameters and boundary conditions keep the same as those in Fig. 12. The initial radius of the b phase is 3.2 lm. Because the nearly planar interface remains unchanged in the rod-like eutectic, the temperature does not change along the direction parallel to the S/L interface and always obtains the extreme at the S/L interface. Fig. 15 shows the distribution of the solute concentration along the x direction corresponding to Fig. 14(b). The distance away from the S/L interface is 0.5 k, and the position is at z = 0.125 Z, 0.25 Z, and 0.375 Z respectively (see the white arrows in Fig. 14(b)). The solute concentration of the rod-like eutectic behaves differently among the three given directions, indicating the difference of diffusion effect from the plate-like eutectic. Nevertheless, the solute concentration always obtains a local maximum (minimum) on the top of the a (b) phase, and it is the fluctuation of solute concentration ahead of the S/L interface that maintains the cooperative growth of eutectic solid phases. 4.4. Approaching realistic Lewis number The employed Lewis number in this work is relatively small compared with the real value (e.g., ~104 for metals). How to simulate eutectic evolution with a realistic Le and how to generalize the present phase-field model into large-undercooling condition are discussed here. The simulation with a larger Lewis number is computationally intensive. A combination of the PFM, the LBM and the multilevel data structure in this work seems a preferable alternative to solve the problem. The time step is inversely proportional to the solute diffusivity in the present discretization scheme, and the stability coefficient in Eq. (13) needs to be determined. Although the largest Le is 500 in the present work (see Fig. 10(b)), it can be theoretically evaluated that Le being larger than 104 can also be readily realized

11

in the multilevel data structure. For example, if Le = 106, one can establish eight–level hierarchical data structure and set key simulation parameters as sLBM = 0.96, and r = 0.01 according to Eq. (15). It is noted that the number of grid levels can reach 11 in the constraint of both stability restriction and convergence issue [49], and the largest grid levels can be up to 50 theoretically if more sophisticated method is employed to ensure the numerical accuracy [43]. In terms of the algorithm, one of the largest breakthrough is that the time step in an explicit discretization scheme can be enhanced by orders of magnitude by constructing a multilevel LBM structure, i.e., the dissimilar time scale is dealt with in a multilevel structure. Similar concept has been validated by Dorari et al. [31], in which they developed a multi-grid-time-step to accelerate the dendritic simulation with a high Lewis number but in a onelevel grid structure. Comparing with the work performed by Dorari et al. [31], the novelty is that a unified multilevel structure is constructed based on a gradient criterion, and the number of the grid level can be artificially determined in a reasonable range. The finest mesh is only preferred near the phase interfaces, which especially benefits the diffuse-interface method like the PFM. However, a larger Le requires a bigger computational domain to cover the thermal diffusion because a larger Le corresponds to faster thermal diffusivity. In this work, we focus on the thermosolutal interaction during eutectic growth, and a smaller Le is employed in a relatively small computational domain instead. In terms of the physical mechanism, the present work, as the first attempt to reveal the thermosolutal interaction during eutectic growth, can pave the way for further quantitative study. Although there is no direct experimental comparison on small Le or related analytical prediction, the accuracy of simulation results can be ensured from the validation of the eutectic model [15,18], thermal evolution (see Sec. 3.1) and contact angle (see Sec. 3.2). Simulating a larger domain with a larger Lewis number, especially in 3D, is in progress, which can not only show the robustness of the numerical scheme but also make better quantitative comparison with bulk eutectics. The remaining problem is how to extend the present phasefield model to a large undercooling, which is challenging and is an ongoing work. Two speculative ideas can be pursed together. The first is to introduce the anti-trapping current like that developed by Folch and Plapp [17]. The second is to employ the multibinary extrapolation in the phase diagram formulation which is used to determine thermodynamic data from free energy databases [50].

5. Conclusions and remark A phase-field lattice-Boltzmann approach is developed to investigate the coupled thermosolutal eutectic phase transition in both 2D and 3D cases. In detail, the PFM is used to simulate the eutectic growth, and the LBM is employed to solve the thermal evolution. To improve the computing efficiency, a parallel-adaptive mesh refinement algorithm based on our recent work [15] is combined with the PFLB approach, which makes coupled thermosolutal eutectic simulation become tractable. The thermosolutal interaction is controlled by the heat sink and Lewis number, and the most remarkable effect of such interaction manifests in the change of lamellar growth velocity. The increase of the heat sink lowers the domain temperature but increases the driving force for solidification, leading to the increase of growth velocity in a linear manner. A larger Le corresponds to faster heat propagation, i.e., making the domain temperature more uniform by decreasing thermal accumulation and temperature difference, which accelerates the eutectic growth. The lamellae creation under large undercooling is expected, which is attributed to the dominated diffusion effect and is considered to maintain a new thermo-

12

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778

dynamic equilibrium by forming a pocket with negative curvature in the middle of the wider phase. Because of the release of latent heat during solidification, the local high-temperature region exhibits at the S/L interface and evolves with the advancement of the S/L interface. However, the largest Le is 500 in the present work and a larger Le will require a bigger computational domain to cover the thermal diffusion. Restricted by the near-equilibrium assumption (or small undercooling) in the present eutectic model, the temperature range cannot considerably deviate from the eutectic temperature. How to simulate eutectic evolution with a larger Le and how to generalize the present phase-field model into large-undercooling condition are discussed. To be closer to the practical situation, a more robust numerical scheme is expected to handle the thermosolutal eutectic growth in a more general behavior, which is in our future work. Declaration of Competing Interest We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled. Acknowledgments This work is financially supported by the Joint Funds of the National Natural Science Foundation of China (grant number U1537202), the Tsinghua-General Motors International Collaboration Project (grant number 20153000354), the Tsinghua University Initiative Scientific Research Program (grant number 20151080370), and the Tsinghua Qingfeng Scholarship (THQF2018-15). The authors would also like to thank the National Laboratory for Information Science and Technology in Tsinghua University for access to supercomputing facilities. References [1] S. Akamatsu, M. Plapp, Eutectic and peritectic solidification patterns, Curr. Opin. Solid State Mater. Sci. 20 (2016) 46–54. [2] J.A. Dantzig, M. Rappaz, Solidification, EPFL Press, Lausanne, 2009. [3] K.A. Jackson, J.D. Hunt, Lamellar and rod eutectic growth, Trans. Metall. Soc. AIME 236 (1966) 1129–1142. [4] J.D. Hunt, K.A. Jackson, Binary eutectic solidification, Trans. Metall. Soc. AIME 236 (1966) 843–852. [5] S. Akamatsu, K. Lee, W. Losert, Control of eutectic solidification microstructures through laser spot perturbations, J. Cryst. Growth 289 (2006) 331–338. [6] J. Du, A. Zhang, Z. Guo, M. Yang, M. Li, S. Xiong, Atomic cluster structures, phase stability and physicochemical properties of binary Mg-X (X= Ag, Al, Ba, Ca, Gd, Sn, Y and Zn) alloys from ab-initio calculations, Intermetallics 95 (2018) 119– 129. [7] J. Du, Z. Guo, A. Zhang, M. Yang, M. Li, S. Xiong, Correlation between crystallographic anisotropy and dendritic orientation selection of binary magnesium alloys, Sci. Rep-Uk 7 (2017) 13600. [8] J. Du, A. Zhang, Z. Guo, M. Yang, M. Li, F. Liu, S. Xiong, Effect of additional solute elements (X= Al, Ca, Y, Ba, Sn, Gd and Zn) on crystallographic anisotropy during the dendritic growth of magnesium alloys, J. Alloy. Compd. 775 (2019) 322– 329. [9] J. Du, A. Zhang, L. Zhang, S. Xiong, F. Liu, Quantitative and qualitative correlations by atomistic determination for the precipitated phases in Al–Li– Cu system, Intermetallics 112 (2019) 106551. [10] A. Zhang, J. Du, Z. Guo, Q. Wang, S. Xiong, Phase-field lattice-Boltzmann investigation of dendritic evolution under different flow modes, Philosoph. Mag. (2019), https://doi.org/10.1080/14786435.2019.1646437. [11] J. Du, A. Zhang, Z. Guo, M. Yang, M. Li, S. Xiong, Atomistic determination of anisotropic surface energy-associated growth patterns of magnesium alloy dendrites, ACS Omega 2 (2017) 8803–8809. [12] J. Du, A. Zhang, Z. Guo, M. Yang, M. Li, F. Liu, S. Xiong, Atomistic underpinnings for growth direction and pattern formation of hcp magnesium alloy dendrite, Acta Mater. 161 (2018) 35–46.

[13] M. Perrut, S. Akamatsu, S. Bottin-Rousseau, G. Faivre, Long-time dynamics of the directional solidification of rodlike eutectics, Phys. Rev. E: Stat. Nonlinear Soft Matter Phys. 79 (2009) 32602. [14] I. Loginova, G. Amberg, J. Ågren, Phase-field simulations of non-isothermal binary alloy solidification, Acta Mater. 49 (2001) 573–581. [15] A. Zhang, Z. Guo, S.M. Xiong, Eutectic pattern transition under different temperature gradients: a phase field study coupled with the parallel adaptivemesh-refinement algorithm, J. Appl. Phys. 121 (2017) 125101. [16] X. Fan, A. Zhang, Z. Guo, X. Wang, J. Yang, J. Zou, Growth behavior of c0 phase in a powder metallurgy nickel-based superalloy under interrupted cooling process, J. Mater. Sci. 54 (2019) 2680–2689. [17] R. Folch, M. Plapp, Quantitative phase-field modeling of two-phase growth, Phys. Rev. E: Stat. Nonlinear Soft Matter Phys. 72 (2005) 11602. [18] S. Gyoon Kim, W. Tae Kim, T. Suzuki, M. Ode, Phase-field modeling of eutectic solidification, J. Cryst. Growth 261 (2004) 135–158. [19] A. Zhang, J. Du, Z. Guo, Q. Wang, S. Xiong, Dependence of Lamellar Eutectic growth with convection on boundary conditions and geometric confinement: a phase-field lattice-Boltzmann study, Metall. Mater. Trans. B 50 (2019) 517– 530. [20] A. Zhang, J. Du, Z. Guo, Q. Wang, S. Xiong, Abnormal solute distribution near the eutectic triple point, Scr. Mater. 165 (2019) 64–67. [21] C.W. Lan, Y.C. Chang, C.J. Shih, Adaptive phase field simulation of nonisothermal free dendritic growth of a binary alloy, Acta Mater. 51 (2003) 1857–1869. [22] J. Rosam, P.K. Jimack, A.M. Mullis, An adaptive, fully implicit multigrid phasefield model for the quantitative simulation of non-isothermal binary alloy solidification, Acta Mater. 56 (2008) 4559–4569. [23] Z. Guo, J. Mi, P.S. Grant, An implicit parallel multigrid computing scheme to solve coupled thermal-solute phase-field equations for dendrite evolution, J. Comput. Phys. 231 (2012) 1781–1796. [24] W.S. Jiaung, J.R. Ho, C.P. Kuo, Lattice Boltzmann method for the heat conduction problem with phase change, Numer. Heat Transf. B-Fund. 39 (2001) 167–187. [25] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen, The Lattice Boltzmann Method Principles and Practice, Springer, Cham, Switzerland, 2017. [26] A. Zhang, S. Meng, Z. Guo, J. Du, Q. Wang, S. Xiong, Dendritic growth under natural and forced convection in Al-Cu alloys: from equiaxed to columnar dendrites and from 2D to 3D phase- field simulations, Metall. Mater. Trans. B 50 (2019) 1514–1526. [27] A. Zhang, J. Du, Z. Guo, Q. Wang, S. Xiong, A phase-field lattice-Boltzmann study on dendritic growth of Al-Cu alloy under convection, Metall. Mater. Trans. B 49 (2018) 3603–3615. [28] A. Zhang, J. Du, Z. Guo, S. Xiong, Lamellar eutectic growth under forced convection: a phase-field lattice-Boltzmann study based on a modified Jackson-Hunt theory, Phys. Rev. E 98 (2018) 43301. [29] A. Zhang, Z. Guo, S.M. Xiong, Quantitative phase-field lattice-Boltzmann study of lamellar eutectic growth under natural convection, Phys. Rev. E 97 (2018) 53302. [30] H. Yin, S.D. Felicelli, L. Wang, Simulation of a dendritic microstructure with the lattice Boltzmann and cellular automaton methods, Acta Mater. 59 (2011) 3124–3136. [31] E. Dorari, M. Eshraghi, S.D. Felicelli, A multiple-grid-time-step lattice Boltzmann method for transport phenomena with dissimilar time scales: Application in dendritic solidification, Appl. Math. Model. 62 (2018) 580–594. [32] I. Rasin, W. Miller, S. Succi, Phase-field lattice kinetic scheme for the numerical simulation of dendritic growth, Phys. Rev. E 72 (2005) 66705. [33] Z. Guo, S.M. Xiong, On solving the 3-D phase field equations by employing a parallel-adaptive mesh refinement (Para-AMR) algorithm, Comput. Phys. Commun. 190 (2015) 89–97. [34] I. Steinbach, F. Pezzolla, A generalized field method for multiphase transformations using interface fields, Physica D 134 (1999) 385–393. [35] A. Zhang, Z. Guo, S. Xiong, Phase-field-lattice Boltzmann study for lamellar eutectic growth in a natural convection melt, China Foundry 14 (2017) 373– 378. [36] J.A. Warren, R. Kobayashi, A.E. Lobkovsky, W. Craig Carter, Extending phase field models of solidification to polycrystalline materials, Acta Mater. 51 (2003) 6035–6058. [37] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [38] T. Seta, Implicit temperature-correction-based immersed-boundary thermal lattice Boltzmann method for the simulation of natural convection, Phys. Rev. E 87 (2013). [39] X. Zhang, J. Kang, Z. Guo, S. Xiong, Q. Han, Development of a Para-AMR algorithm for simulating dendrite growth under convection using a phasefield–lattice Boltzmann method, Comput. Phys. Commun. 223 (2018) 18–27. [40] J. Du, A. Zhang, Z. Guo, M. Yang, M. Li, S. Xiong, Mechanism of the growth pattern formation and three-dimensional morphological transition of hcp magnesium alloy dendrite, Phys. Rev. Mater. 2 (2018) 83402. [41] M. Berger, I. Rigoutsos, An algorithm for point clustering and grid generation, IEEE Trans. Syst. Man Cybern. 21 (1991) 1278–1286. [42] A. Zhang, Z. Guo, Q. Wang, S. Xiong, Three-dimensional numerical simulation of bubble rising in viscous liquids: a conservative phase-field latticeBoltzmann study, Phys. Fluids 31 (2019) 63106.

A. Zhang et al. / International Journal of Heat and Mass Transfer 145 (2019) 118778 [43] O. Filippova, D. Hänel, Grid refinement for lattice-BGK models, J. Comput. Phys. 147 (1998) 219–228. [44] A. Dupuis, B. Chopard, Theory and applications of an alternative lattice Boltzmann grid refinement algorithm, Phys. Rev. E 67 (2003) 66707. [45] W. Kurz, D.J. Fisher, Fundamentals of Solidification, Trans Tech Publications Ltd., Switzerland, 1992. [46] N. Marasßli, J.D. Hunt, Solid-liquid surface energies in the Al-CuAl2, Al-NiAl3 and Al-Ti systems, Acta Mater. 44 (1996) 1085–1096. [47] M. Gündüz, J.D. Hunt, The measurement of solid-liquid surface energies in the Al-Cu, Al-Si and Pb-Sn systems, Acta Mater. 33 (1985) 1651–1672.

13

[48] A. Solomon, some remarks on the Stefan problem, Math. Comput. 20 (1966) 347–360. [49] J. Töelke, S. Freudiger, M. Krafczyk, An adaptive scheme using hierarchical grids for lattice Boltzmann multi-phase flow simulations, Comput. Fluids 35 (2006) 820–830. [50] J. Eiken, B. Boettger, I. Steinbach, Multiphase-field approach for multicomponent alloys with extrapolation scheme for numerical application, Phys. Rev. E 73 (2006).