Phase field study of the tip operating state of a freely growing dendrite against convection using a novel parallel multigrid approach

https://doi.org/10.1016/j.jcp.2013.10.004Get rights and content

Abstract

Alloy dendrite growth during solidification with coupled thermal-solute-convection fields has been studied by phase field modeling and simulation. The coupled transport equations were solved using a novel parallel-multigrid numerical approach with high computational efficiency that has enabled the investigation of dendrite growth with realistic alloy values of Lewis number 104 and Prandtl number 102. The detailed dendrite tip shape and character were compared with widely recognized analytical approaches to show validity, and shown to be highly dependent on undercooling, solute concentration and Lewis number. In a relatively low flow velocity regime, variations in the ratio of growth selection parameter with and without convection agreed well with theory.

Introduction

In the shaped casting industry the growth behavior of dendrites in a solidifying alloy controls the as-cast microstructure and has a strong influence on final component mechanical properties. Cast structures in practice, even in simple binary alloys, are complex and rarely conform to easy classification as homogeneously columnar or equiaxed, and frequently present complex cellular/dendritic patterns that vary from place to place. There has been significant effort to better understand the underlying physics controlling the shape, length scale and solute redistribution processes occurring at a growing dendrite tip in an attempt to control the factors that determine final cast microstructure [1], [2], [3], [4], [5]. Both analytical and numerical approaches have been developed, but despite the well-known strong influence of liquid movement and convection on final microstructure in practice, only a small number of recent studies have begun to account for its influence on the prior, more developed thermal-solute approaches.

The operating state of a growing dendrite can be defined by the tip radius Rtip and the tip velocity vtip. By assuming the tip to be a parabola (in 2-D) or a paraboloid of revolution (in 3-D) with parabolic tip radius Rp and the steady dendrite is isothermal with the solid at the melting temperature, Ivantsov [6] proposed the most widely quoted relationship for dendrite operating state for a purely thermally-controlled growing dendrite, comprising the relationship between external imposed undercooling Δ=(TmT)/(L/Cp) and the thermal Peclet number at the tip PeT=Rpvtip/(2α) as Δ=Iv(PeT), where Tm is the melting temperature, T is the temperature of the undercooled melt, L is the latent of fusion, Cp is the specific heat, α is the thermal diffusivity and Iv(x)=πxexp(x)erfc(x) is the Ivantsov function (in 2-D). Ivantsovʼs theory predicted that for a given undercooling, there were infinite pairs of (Rtip, vtip) for the solution of the linked expression since only their product (Peclet number) could be determined. While convenient, this implication is in conflict with experiment where Rtip and vtip are invariant for a fixed undercooling. Two subsequent approaches were then developed by the introduction of a selection constant defined as σ. The marginal stability theory was developed by Langer and Muller-Krumbhaar [7] and involved another relationship between the tip radius and velocity given as Rp2vtip/(d0α)=(1/σ)2 where d0 is the material thermal capillary length. Drawing on a stability analysis [8] based on the allowable shape of a perturbed, non-flat solid–liquid interface, they proposed σ=1/(2π). Ben-Jacob et al. [9], [10] and Kessler et al. [11], [12] developed other approaches that allowed for anisotropic surface energy to give a single, paired solution for Rp and vtip, deduced from the fastest growing mode of perturbed solid–liquid interface, which led to an expression similar to the one given by Langer and Muller-Krumbhaar [7] i.e. Rp2vtip=constant. Kessler and Levine [13] extended this idea and found that the dendrite tip shape computed in this way generally displays a cusp (non-zero slope) at the tip and at a unique (Rp, vtip) pair; the cusp reduces to a smooth shape with zero slope at the tip, which is called the microscopic solvability condition. Further numerical experiments revealed that the selection constant σ was dependent on the strength of the surface anisotropy ε i.e. Rp2vtip=f(ε). Nevertheless, experimental validation of these increasingly complicated analytical/numerical approaches has been difficult since they rely on controlling stable and well-characterized growth conditions, generally far from the more dynamic conditions expected in practice [14].

The extension of the microscopic solvability theory to binary alloys, where both solute and thermal diffusion are important, was performed by Lipton, Glicksman and Kurz (LGK) [15] and Lipton, Kurz and Trivedi (LKT) [16]. These approaches are also characterized by the use of a selection constant σ (σLGK and σLKT will be used for the LGK and LTK theories respectively):σ=d0RpPeT[ξT+2ξcLe(Mc1(1k)ΔC)] where M=|m|(1k)/(L/Cp) is the scaled dimensionless liquidus slope, m is the actual liquidus slope from the phase diagram, k is the solute partition coefficient, Le=α/D is the Lewis number, D is the solute diffusivity in liquid, c is the solute concentration and c is the equilibrium solute concentration. ΔC=(ctipc)/((1k)ctip) is the dimensionless solutal undercooling and ctip is the solute concentration at dendrite tip. For the LGK theory, both ξT and ξc are unity but for the LKT theory:ξc=1+2k12k1+1σ(Pec)2 andξT=111+1σ(PeT)2

The overall undercooling is then given by:ΔT=LCpΔT+kΔT0Δc1(1k)Δc+ΓRp where the three terms on the right correspond to thermal, solutal and capillary undercooling, respectively. ΔT0=|m|c(1k)/k is the equilibrium freezing range corresponding to c and Γ is the Gibbs–Thomson coefficient. ΔT=(TtipT)/(L/Cp) is the dimensionless thermal undercooling. Eqs. (1) and (4) together uniquely determine the tip radius and tip velocity.

Convection in the melt – almost always significant in practice – has long been realized to have a profound effect on dendritic growth [17]. But it is presently unclear how the preceding theories (LGK and LKT) for binary alloys may remain valid or how they might be modified when convection is present. Ananth and Gill [18] and Saville and Beaghton [19] studied the motion of the freezing front between a needle-shaped crystal and a supercooled liquid for situations where there is forced convection aligned along the crystal growth. Analysis was conducted by modeling the transport problem for a pure material solidifying as a paraboloid of revolution in an infinite undercooled melt. The imposed external undercooling could be characterized by the thermal Peclet number PeT, the flow Peclet number Pef=Rtipv/2α (where v is the imposed external flow velocity) and the Prandtl number Pr=υ/α (ratio between kinematic viscosity and thermal diffusivity) i.e. Δ=Δ(PeT,Pef,Pr). Through a so-called linear solvability analysis, Bouissou and Pelce [20] considered the stability of this solution, and found that the ratio of the selection parameters with convection (σ) and without convection (σ0) could be characterized by a dimensionless parameter χe:σ0σ=1+bχe11/14 where b is a numerical constant,χe=a(Re)d0v/((15ε)0.75Rpvtip) anda(Re)=2Re/πexp(Re/2)/erfc(Re/2) where Re =vRp/υ is the Reynolds number and ε is the anisotropic strength. The theory predicts that parameter b in Eq. (5) is very close to zero i.e. σ0σ when χe1, and therefore the flow only has an effect on the tip-selection parameter if d0v/(Rpvtip)=Pefσ is of the order of unity or greater [20].

The applicability of this approach has not yet been confirmed conclusively by numerical modeling or experiments. Beckermann et al. [21] found a weak dependence of the selection parameter on flow. Tong et al. [22] adopted a 2-D thermal convection phase field model and studied the tip operating state corresponding to a pure material dendrite growing against forced convection. Within the parameter range of 0<χe<0.2, the simulation results agreed well with theoretical predictions and the ratio of the selection parameters with and without convection was very close to unity, as suggested by Bouissou and Pelce. Jeong et al. [17] performed 3-D phase field simulations and compared the results with available theories and experimental results: the selection constant σ decreased slightly in the presence of convection, agreeing with the solutions of both Saville and Beaghton [19] and again with Bouissou and Pelce [20], although the idealized conditions required in the theory were unlikely to be present during the experiments. A similar 3-D simulation approach was also adopted by Lu et al. [23] who showed again reasonable agreement with the Bouissou and Pelceʼs theory that the selection constants under different conditions were similar whether convection was, or was not, present.

Apart from these pure material studies, binary alloys have also been considered in which the influence of external flow on dendrite growth was investigated. By employing an adaptive mesh methodology, Lan et al. [24] studied the influence of thermal-solutal convection on morphology changes of a freely growing dendrite. They showed the importance of employing the so-called “anti-trapping” current developed by Karma [25] in the phase field model in order to produce qualitatively acceptable predictions. Recently, Siquieri et al. [26] calculated the shape, using the phase field method, of a freely growing alloy dendrite under isothermal conditions where convection was present.

The phase field method for simulating dendritic solidification and other phase transitions has become popular in recent years due to its thermodynamic rigor and the perceived benefits of avoiding explicit tracking of the complex shaped and evolving solid–liquid interface. The phase field method adopts one or several order parameters to characterize the transition of the phases during solidification. In particular, to model the dendritic growth of alloys usually one order parameter, i.e. the phase filed ϕ, is used to distinguish the solid and liquid phases in which ϕ varies smoothly but steeply from −1 in bulk liquid to 1 in bulk solid over the diffuse interface region of width W0.

As discussed in our earlier work [27], despite apparent benefits, a major problem for phase field modeling remains the enormous computing overhead required to resolve the detail of the solid–liquid interface in complex dendritic geometries. Several methods have been developed and proposed to help address the practical difficulties, of which the most important was the so-called thin interface limit analysis, developed by Karma and co-workers [25], [28], [29]. This methodology enabled a diffuse interface width larger than the capillary length, which has the practical effect of greatly reducing the computing overhead. However, for alloy solidification processes in the casting industry, where multiple dendrite growth proceeds according to coupled thermal, solutal and fluid flow fields, use of the phase field approach has been restricted by the need thus far to make limiting simplifications to the real process physics.

The primary difficulties arise from three aspects: (1) the different length- and time-scales in the coupled energy and mass transport processes (high Lewis number Le=α/D104 and low Prandtl number Pr=υ/α102 for metallic alloys); (2) enormous computing demands arising from the requirement to achieve practically useful and realistic spatial resolution of the solid–liquid interface; and (3) a lack of robust numerical schemes that can solve such a complex multi-scale problem in a reasonable time span.

In this paper, an approach combining a full approximation storage (FAS) multigrid [30] algorithm and a parallel computing architecture is developed for the first numerical solutions of coupled thermal-solute-convection fields in the dendritic growth of a real alloy. The phase field model including convection and the parallel multigrid approach are firstly described and then numerical experiments are performed to investigate convergence and computational efficiency. Predictions of dendrite tip radius, tip velocity and particularly the selection constant are compared with widely accepted theory as a way to validate various behaviors predicted by the model under different liquid flow conditions. In particular, we focus on the effect of convection on the selection constant σ according to Bouissou and Pelce [20] as this requires a study of the linked effects of temperature, solute concentration and flow.

Section snippets

Mathematical model

The following assumptions were made:

  • (1)

    Physical properties of materials including thermal diffusivity α, liquid solute diffusivity D, latent heat of fusion L, alloy specific heat Cp, and kinematic viscosity υ are constant unless stated otherwise.

  • (2)

    Solute diffusion in the solid phase is neglected.

  • (3)

    Only forced convection is considered, i.e. the effect of gravity is neglected.

  • (4)

    The flow is incompressible i.e. the alloy density ρ is constant.

Formulation of the phase field model starts from the

Discretization of equations

It is easy to notice that Eqs. (17), (18), (19), (20), (21) can be re-written in a simple form as:ΦtEt=Φl2E+ΦsE+ΦxEx+ΦyEy+ΦC+ΦE where E is the target variable such as ϕ, U and θ. Φt, Φl, Φs, Φx, Φy and ΦC are constants which are always related to the values from the previous time step, ΦE is a non-linear term related to variable E which is non-zero only for phase field. Eqs. (17), (18), (19), (20), (21) were then all discretized onto a 2-D rectangular domain (as shown in Fig. 1) with

Benchmark simulations

Fig. 3 shows the configuration of the benchmark dendrite growth scenario. The 2-D computation domain is rectangular with M×N cells. Only half of the growing dendrite was simulated because of the symmetry of the case. The dendrite grew from the midpoint of the left side of the domain with an initial solid “seed” radius of R0=30d0. The initial temperature of the seed was set to be zero while for the rest of the domain the temperature was set at an initial undercooling θ=θ0. For phase field,

Operating state and dendrite morphology transition

A parametric study of dendrite growth was performed using the same domain as shown in Fig. 3 and full details are given in Table 2. There were 126 simulations and five main categories of parameters were investigated, including Mc, Δθ, Le, Pr and v, i.e. the main factors that could influence the growth of the dendrite tip. The spatial step was fixed at Δx=Δy=0.8 and the time step at dt=0.8Δx2/(4D). The size of the computing domain was chosen according to the parameters used in the simulation,

Conclusions

Free dendritic growth against forced convection has been investigated using phase field modeling in which fully coupled and non-linear equations for heat, solute and liquid flow were solved using a novel parallel-multigrid numerical approach. Convergence and efficiency tests confirmed that this approach was very robust and enabled the solution of coupled dendrite growth during solidification of alloys of practical interest with Le104 and Pr102. Dendrite growth dynamics were studied by

Acknowledgements

The authors would like to thank the Natural Science Foundation of China (Project No. 51205229), UK Royal Academy of Engineering/Royal Society through the Newton International Fellowship Scheme, and the EPSRC Centre for Innovative Manufacture: Liquid Metal Engineering (EP/H026177/1) for financial support, and the Oxford Supercomputer Centre and the National Laboratory for Information Science and Technology in Tsinghua University for access to supercomputing facilities and support for parallel

References (37)

  • E. Ben-Jacob et al.

    The formation of patterns in non-equilibrium growth

    Nature

    (1990)
  • L. Granasy et al.

    A general mechanism of polycrystalline growth

    Nat. Mater.

    (2004)
  • T. Haxhimali et al.

    Orientation selection in dendritic evolution

    Nat. Mater.

    (2006)
  • W. Kurz et al.

    Fundamentals of Solidification

    (1998)
  • G.P. Ivantsov

    Temperature field around spherical, cylindrical and needle-shaped crystals which grow in supercooled melt

    Dokl. Akad. Nauk SSR

    (1947)
  • W.W. Mullins et al.

    Stability of a planar interface during solidification of a dilute binary alloy

    J. Appl. Phys.

    (1964)
  • E. Ben-Jacob et al.

    Dynamics of interfacial pattern formation

    Phys. Rev. Lett.

    (1983)
  • E. Ben-Jacob et al.

    Boundary-layer model of pattern formation in solidification

    Phys. Rev. A

    (1984)
  • Cited by (36)

    • A review on solidification of alloys under hypergravity

      2023, Progress in Natural Science: Materials International
    • Phase-field lattice-Boltzmann study on dendritic growth of hcp metals under gravity-driven natural convection

      2023, Transactions of Nonferrous Metals Society of China (English Edition)
    • Phase-field-lattice Boltzmann method for dendritic growth with melt flow and thermosolutal convection–diffusion

      2021, Computer Methods in Applied Mechanics and Engineering
      Citation Excerpt :

      Most of those PFM–LBM models can be considered as hybrid models in which finite-difference- or finite-volume-based PFM was applied to simulate the phase field evolution, while the LBM was implemented to model the melt flow and heat and solute transfer. In addition, fully coupled PFM models considering all the effects of melt flow and thermosolutal convection–diffusion in the literature are very rare (e.g., [30]) due to the lack of general, convenient, and efficient numerical schemes. Recognizing the capabilities and advantages of the LBM, there has been growing interest in constructing LB schemes to solve the governing equation for the phase field [31–33]; as a result, the generic LB algorithm, and thus a single grid system, can be applied to model all the transport phenomena as well as the phase field evolution.

    View all citing articles on Scopus
    View full text