Quantification of mesh induced anisotropy effects in the phase-field method

https://doi.org/10.1016/j.commatsci.2005.02.017Get rights and content

Abstract

Phase-field modelling is one of the most powerful techniques currently available for the simulation from first principles the time-dependant evolution of complex solidification microstructures. However, unless care is taken the computational mesh used to solve the set of partial differential equations that result from the phase-field formulation of the solidification problem may introduce a stray, or implicit, anisotropy, which would be highly undesirable in quantitative calculations. In this paper, we quantify this effect as a function of various computational parameters and subsequently suggest techniques for mitigating the effect of this stray anisotropy.

Introduction

The dendrite is a solidification morphology characteristic of many metallic melts and certain other systems with a low entropy of fusion. It is a structure of both great practical and theoretical interest. Dendrites are parabolic needle crystals that develop complex, time-dependent shapes, normally as the result of extensive branching which gives rise to a tree-like structure.

During the processing of metallic components, solidification from the parent melt is almost invariably the first step. Remnants of the dendritic microstructures formed during solidification often survive subsequent processing operations, such as rolling and forging, and the length scales established by the dendrite can influence not only the final grain size but also micro- and hence macro-segregation patterns. This can have a wide-ranging influence on both the properties of finished metallic products, affecting for instance mechanical properties, corrosion resistance and surface finish, and on the formability of metallic feedstock, such as the ability to resist hot tearing during rolling.

Theoretical interest stems from the fact the dendrite is a prime example of a pattern forming system where complex morphologies arise from initially homogeneous conditions due to the highly non-linear response of the controlling system. Moreover, although the governing equations for dendritic growth have been known for many decades, finding solutions to the free-boundary problem, even in the tip region has proved enormously complicated. Indeed, finding analytical steady-state solutions for the radius of curvature of the dendrite at its tip has proved to be beyond all orders of perturbation theory [1].

The problem of predicting the operating point of the dendrite first became apparent in 1947 when Ivantsov [2] showed that an isothermal paraboloid of revolution, growing with radius of curvature R and (tip) velocity V into an isotropic undercooled melt was a shape preserving solution to the diffusion equation, thus giving rise to the idea of the needle dendrite. However, the Ivantsov analytical solution for such a crystal is degenerate in that it relates the Peclet number, and not the growth velocity, to undercooling, where the Peclet number is defined asp=VR2Dlwhere Dl is the diffusivity in the melt. Consequently, at a given undercooling an infinite set of solutions are admissible, subject to the condition VR = constant. Such degeneracy is not observed in nature, where a well-defined growth velocity can always be associated with a given set of initial conditions.

A rigorous approach to the problem of selecting the operating point of a needle crystal is provided by the theory of microscopic solvability [3], [4]. The principal physical insight of solvability theory is that surface tension acts as a singular perturbation which resolves the degeneracy found in the macroscopic problem. Perhaps counter to intuition, it turns out that in the case of an isotropic system the equations have no solution. The principal prediction of solvability theory is that capillary forces break the Ivantsov degeneracy via the relationshipR2V=2Dd0σwhere d0 is the capillary length, which is given byd0=σTmcL2where σ is the interfacial energy between the solid and liquid phases, L is the latent heat per unit volume, c the specific heat per unit volume and Tm is the equilibrium melting temperature of the solid. σ is the anisotropy dependant eigenvalue for the problem, where for a cubic system anisotropy is usually introduced by writingd=d0(1-γdcos4θ)where γd is the anisotropy strength. For small Peclet numbers σ is found to vary as σ(γd)γd7/4 in the limit p, γd  0.

Although it has allowed great advances in our understanding of steady-state dendritic growth, solvability theory is not suited to the time-dependant growth problem. One of the central advances in the ability to predict non-steady dendritic growth in the last 20 years has been the advent of phase-field modelling. First proposed by Langer [5] and subsequently developed by, among others, Caginalp [6] and Penrose and Fife [7], the basis of the phase-field method is the definition of a phase variable (say ϕ) the value of which describes the phase state of the material. At it simplest, for a single-phase solid, this might for instance equate ϕ  1 with the solid and ϕ  −1 with the liquid. The central assumption of the phase-field method is that the interface between phases is diffuse, with a finite width δ, and that ϕ is thus a continuous variable over the whole domain Ω on which the problem is defined. For the simple single-phase system discussed above this would allow ϕ to take values intermediate between the solid and liquid end members, that is, −1  ϕ  1. Unlike the classical Stefan problem, in which the solidification front is treated as a mathematically sharp interface at which boundary conditions are applied, the continuity of ϕ over Ω allows the equations that govern the evolution of ϕ to be written in a differential form. These are usually derived from either a free-energy or entropy functional which incorporates conservation of energy and positive entropy production, with the Gibbs–Thomson condition relating the local interface temperature to the thermodynamic equilibrium temperature, local interface curvature and interface kinetics. This phase-equation, coupled to a transport equation for either heat or solute, can then be solved using conventional techniques for partial differential equations.

The phase-field method has facilitated significant progress in our understanding of a number of problems associated with time-dependant dendritic solidification including the study of dendritic shapes at high undercoolings [8], a model for spontaneous grain refinement [9] and the inclusion of electric currents through the solidifying material [10]. However, there are a number of problems with the technique, one of which is the implicit anisotropy introduced into the solution by the mesh on which the differential equations are solved. Although most workers in the field recognise this as a potential pitfall when using phase-field methods the problem has received little systematic study. The problem arises for the following reasons:

  • The finite difference (FD) and finite element (FE) routines generally used to solve the differential phase and transport equations normally discretise the computational domain using a regular mesh. The regularity of the mesh introduces a directionality which is equivalent to an implicit or mesh induced anisotropy to the solution.

  • Dendritic solidification is very sensitive to small amounts of anisotropy. Consequently, the small level of anisotropy introduced by using FE or FD solvers with a regular mesh can have a significant effect in phase-field simulations when, in many other computational models, this effect would be insignificant.

This work was largely motivated by our observation that mesh induced anisotropy was a far more significant problem in solutal phase-field models than in thermal models. In this paper, we set out to systematically analyses the origins of mesh induced anisotropy, accounting for why it is more problematic in solute based systems, to quantify the mesh induced anisotropy as a function of the relevant computational parameters and to suggest mitigating strategies.

Section snippets

Computational model

Our investigation was conducted using the single-phase model of Warren and Boettinger [11]. As with most phase-field models the basis of the Warren and Boettinger model is the definition of an entropy functional [12], [13], which in this case isS=Ωs(ϕ,e,c)-12ε2(ϕ)2dΩwhere s is the thermodynamic entropy density, ϕ is the phase variable taking values 0 in the solid and 1 in the liquid, e is the internal energy density and c is the concentration of solute (chemical species B) in the solvent

Determination of the mesh induced (implicit) anisotropy

In the absence of anisotropy no stable solutions exits for a solid nucleus growing out into its undercooled or supersaturated parent melt. Consequently, an initial circular nucleus will grow to a critical radius, determined usually by the surface energy, σ, before breaking up into random fingers. However, until the critical radius is approached, an initially circular nucleus will, to a very good approximation, retain its circular shape. When anisotropy is present the evolution of the

Results

As mentioned above this work was motivated by our observation that mesh induced anisotropy is a far more severe problem in solutal models than it is in thermal models. Indeed, we have applied the methodology described above to measure the mesh induced anisotropy in a thermal-phase field model [16], [17] which uses an essentially identical phase-equation and found γm4 to be <0.002.

Our initial feeling was therefore that the observed difference in the strength of the mesh induced anisotropy was

Discussion

The results presented here indicate that, conceptually, the simplest way to minimise the implicit computationally induced anisotropy is to ensure that there is a high density of mesh cells within the interface region. However, if a regular meshing is utilised this can rapidly lead to an unmanageably large grid with correspondingly long computation times. One potential route to introducing a very fine grid in the interface region while preserving reasonable computational efficiency over the rest

Summary and conclusions

Phase-field modelling is one of the most powerful techniques currently available for the simulation from first principles of solidification microstructures. Yet, there is an inherent problem that the nature, and indeed existence, of steady-state solutions is dependant on a small parameter, namely the crystalline anisotropy. Moreover, the effective value of this parameter can be influenced by the computational mesh used to solve the set of partial differential equations that result from the

References (34)

  • O. Penrose et al.

    Physica D

    (1990)
  • R. Gonzalez-Cinca

    Physica A

    (2002)
  • A.M. Mullis et al.

    Acta Mater.

    (2001)
  • L.N. Brush

    J. Cryst. Growth

    (2003)
  • J.A. Warren et al.

    Acta Metall. Mater.

    (1995)
  • S.-L. Wang et al.

    Physica D

    (1993)
  • R. Kobayashi

    Physica D

    (1993)
  • A.A. Wheeler et al.

    Physica D

    (1993)
  • A.M. Mullis

    Acta Mater.

    (2003)
  • R.J. Braun et al.

    J. Cryst. Growth

    (1997)
  • N. Provatas et al.

    J. Comput. Phys.

    (1999)
  • R. Tonhardt et al.

    J. Cryst. Growth

    (2000)
  • E. Burman et al.

    J. Comput. Phys.

    (2004)
  • M. Plapp et al.

    J. Comput. Phys.

    (2000)
  • A. Karma et al.

    J. Cryst. Growth

    (1997)
  • S. Chen et al.

    J. Comput. Phys.

    (1997)
  • E. Brener et al.

    Phys. Rev. E

    (1995)
  • Cited by (0)

    View full text