A sharp interface model for deterministic simulation of dendrite growth

Highlights

• The sharp front and capillary effects are modelled using a level set method.

• A robust ghost-fluid technique is used to apply boundary conditions on the front.

• Speed-curvature relation is quantified in terms of a probability distribution.

• The proposed strategy captures coarsening phenomenon and dendritic structures.

Abstract

A high order level set model is developed for deterministic simulation of dendritic growth in unstable solidifying systems. The model captures motion of the front implicitly on a structured finite difference grid, enables calculation of its geometric properties and also applies boundary conditions on the immersed interface. Interfacial capillary effect is incorporated in the model through the Gibbs-Thomson condition. Canonical problems for evaluating grid convergence of the numerical method and validation tests for stability of a growing nucleus in the presence of isotropic surface tension are presented. The growth morphology of solidifying nuclei in undercooled metallic melts is quantitatively analyzed. Effects of crystal anisotropy and melt undercooling on the front geometry, propagation speed and formation of branched dendritic structures are examined. The complex morphological changes, such as remelting of secondary perturbations under specific conditions of undercooling, are also captured during the quantitative analysis.

Introduction

Phase change phenomenon plays a key role in numerous scientific problems observed in nature. It is of particular interest in materials processing and manufacturing. The first attempt at a mathematical description of phase change was intended at modeling melting of icebergs in huge water bodies [1], popularly known as the Stefan problem. It formulates the motion of a phase front in a continuum domain by solving for the energy balance at the front. Although related, phase transformation in materials involves additional complexities. Microstructures of metals and alloys are known to exhibit complex morphology while undergoing solidification. These microscopic structures are important from an engineering perspective as they strongly influence the mechanical behavior of the manufactured part. Understanding the phase transformation process with particular focus on solidification morphology is extremely important for the development of advanced manufacturing and material processing strategies.

Stefan’s analysis provided a very basic interpretation of the phase change kinetics on a macroscopic scale which is useful in certain systems. However, it is not adequate to model microscopic interfacial effects in metals and alloys undergoing solidification [2]. Stefan’s model addresses stable solidification where the interface is flat without protrusions and the front motion is related to the rate of latent heat evolution at the interface. When the molten metal is maintained at a temperature below the equilibrium (i.e. negative temperature gradient imposed on the interface) it results in unstable solidification. Small perturbations on the solid/liquid interface have a tendency to grow faster, resulting in dendritic and columnar microstructures. Modeling the development of such structures requires understanding the physics coupling the front geometry with bulk phase. Various theories have been proposed to account for such pattern formations in unstable solidifying systems. Ivantsov’s analysis [3] was the first attempt to account for the non-planar structures found experimentally. He assumed specific shapes of the structures and proposed a relation between the growth rate and geometry for an isothermal tip. His criterion was later used by other researchers to develop models for columnar growth [4]. These models were suited to predict tip speeds of growing columnar structures but did not address mechanics of nucleation of the structures. The Mullins and Sekerka [5] instability theory describes a criterion for a perturbed front to transform into columnar structure. They assumed sinusoidal perturbations and obtained the critical wavelength for a metastable state in terms of transport coefficients for diffusion governed growth. Whether the perturbations amplify or decay in time can be predicted by comparing them to this critical value.

Ivantsov’s analysis based on isothermal tip is useful to account for columnar structures but not the branched dendritic morphology observed in thermal dendrites in metals. Capillary effects [6] on solidification of an undercooled melt play a key role in determining the front morphology. For a perturbation to grow, the change in volumetric free energy must overcome the surface energy barrier. The effective driving force for growth of a perturbation is less than that corresponding to an isothermal tip assumed by Ivantsov. The Gibbs-Thomson relation [6] describes the effective local undercooling in terms of interface curvature. The undercooling can be transformed into an interface temperature which acts as a Dirichlet boundary condition relating the front geometry with bulk transport equations. The condition can also be modified to incorporate crystalline anisotropy to govern the orientation of dendrites [7]. Details regarding the functional form of anisotropy will be discussed in later section. In addition to curvature effects, atomic mobility and varying material properties also influence the local undercooling and can be incorporated in the Gibbs-Thomson equation [8].

Jaafar et al. [8] outlines various numerical approaches to model the process based on Gibbs-Thomson theory. The models could be broadly classified as implicit or explicit based on the method of tracking the solid-liquid boundary. Explicit methods typically require interface construction from a collection of markers distributed along the front while implicit approaches track the front using distribution of a scalar function throughout the domain. Front tracking using marker particles [9], [10] are examples of explicit tracking methods where interface can be constructed using high order polynomials from a collection of particles. Phase field [7], [11], [12] and level set methods [13], [14], [15] are examples of implicit tracking methods where the front is represented as an iso-contour of some indicator function. The function itself is governed by a transport equation to account for local front motion using appropriate physics. Application of boundary conditions on the front is usually straightforward with explicit tracking methods while topological changes in the front are easily handled by implicit tracking approaches as these are embedded into a scalar function. Jaafar et al. [8] also outline other implicit tracking methods such as the enthalpy-porosity technique [16], [17] and volume of fluid technique [18], [19], [20]. However, in the context of sharp interface problems, they are less suitable due to numerical challenges.

Phase field function (order parameter) [7], [11] is typically similar to a step function smoothed across a small transition thickness near the front. The function takes a value of $-1$ in one domain and 1 in the other. The zero iso-contour represents the front. A governing equation for order parameter can be written by imposing minimization of total energy of the system. Interfacial energy, latent heat and varying material properties can be incorporated in the governing equation. The model asymptotically reduces to a sharp interface assumption as the transition thickness vanishes. A thin transition region also requires fine grid sizes and ultimately results in higher computational cost. On the other hand, level set method [21] maintains the front to be a truly sharp interface. The level set function is governed by a simple advection equation that transports a signed distance function using a modeled continuum velocity. This approach enables straightforward computation of interface geometry and application of boundary conditions on the interface. A disadvantage of the level set method is that it is known to suffer from unphysical mass gain/loss [15], [22], [23], [24]. Through various canonical tests, it has been shown in Ref. [22], [25], [26], [27], [28] that a high order accurate formulation of the function offers superior conservation property as well as accurate computation of front geometry. The level set method will be briefly outlined in later section. A brief comparison of the phase field and level set approaches can be found in Ref. [29].

Our objective is to use the level set model to simulate solidification process and understand the branched dendritic structures formed in undercooled metalic melts. As the level set model enables accurate computation of interface geometric properties, it becomes a more suitable approach to model local curvature effects on propagating fronts. Our aim is also to quantitatively study the relationship between local front propagation and its geometry under various environments. Below is a brief review of previous attempts to model dendritic microstructure with particular focus on the challenges encountered.

A wide range of literature is available on the use of phase field method for simulating branched solidification structures. The works of Kobayashi [7] and Karma and Rappel [30], [31] provide a detailed view of how the method can capture complex dendritic features with relative simplicity. The method requires addition of thermal noise to trigger the nucleation and cause morphological instability on the growing solidification front to transform it into dendritic microstructure. The dendritic features seemed to be affected by the amplitude of thermal noise [31]. A recent study by Mullis [12] is aimed at producing deterministic branched structure using the phase field method. Moreover, the smallest resolvable feature using this method is limited by the thickness of diffuse interface.

Juric and Tryggvason [9] presented a detailed sharp interface model for solving solidification problems. Their approach was based on front tracking using marker points to represent the interface. Udaykumar et al. [10] also used similar technique for computation of solidification fronts in the sharp interface limit. They showed linear and quadratic interpolation techniques for obtaining gradients on a finite difference grid near the interface. The method was able to capture branched dendritic microstructure. One of the major challenges in front tracking methods is to maintain uniform point density along the front as it undergoes expansion. Addition of interpolated markers is frequently required. Other explicit interface tracking schemes use finite element approach [32] where body fitted grid can be generated and interface boundary condition can be applied exactly at the triangulated front. Frequent mesh generation in expanding fronts is a limitation. Moreover, both front tracking and unstructured mesh methods pose challenges in handling topological changes of the front.

Various groups have used the level set method to implicitly capture the sharp interface. Gibou et al. [13], [33] developed a second order accurate discretization scheme for solving the heat equation in presence of embedded boundaries and applied it to simulate Stefan problems. They used the classic level set formulation [21] for capturing the front. The application was limited to simulating canonical problems and did not account for the formation of branched dendritic structures. Chen et al. [14] used the level set method [21] to simulate the unstable solidification problem in different domains. The Mullins-Sekerka instability, four and six fold anisotropic growth were simulated. Comparison of canonical test problems with analytical solutions was also shown. To account for immersed boundary, they developed a complex four-coordinate scheme for near front discretization of the gradients. Later, Tan and Zabaras [34] formulated a hybrid front tracking and fixed-domain approach based on level set function. The front itself was tracked using marker particles while the level set information was used to solve energy and flow equations on a finite difference grid. An artificial mushy region (diffuse interface) was introduced to solve the flow equations while one-sided gradient interpolation was applied to obtain fluxes near the solidification front. The Dirichlet boundary condition for interface temperature was incorporated in the interpolation scheme. They were able to simulate dendritic solidification in the presence of forced convection. However, the method does not maintain a sharp solidification front. A high order discretization method in the presence of an interface has also been outlined by Baeza et al. [35] which seems to be a challenge in all works pertaining to phase change problems.

The literature review indicates that level set technique has immense potential in simulating solidification behavior in materials under the sharp interface limit. The phase change behavior itself has been thoroughly studied but the problem of simulating branched dendritic structures using a sharp interface model still remains a challenge. As noted by theoretical studies, the dendritic growth is highly dependent on the interface geometry and its coupling with bulk phase physics. On a structured finite difference grid, addressing the immersed boundary problem also seems to pose challenges and a robust method is yet to be devised. In the present work, our aim is to formulate a high order level set method which captures the location and geometry of the front accurately without the need for explicit construction. A robust formulation for applying boundary condition at an immersed interface using level set information is also presented. Our objective is also to use the model for deterministic simulation of thermal dendrites in pure metals undergoing solidification and study the effect of front geometry on local propagation speed. The use of a second order level set method ensures that geometric properties of the interface can be calculated accurately [25]. Although the problem of dendritic solidification has been studied extensively, a quantitative strategy to describe the morphology and its dependence on front propagation is fewer in the literature. In the present discussion we use level set approach for modeling interface motion as well as to extract information about the dendrite shape.

The numerical formulation is described in the next section, followed by the results and discussion. Standard validation tests from literature are considered to evaluate accuracy of the numerical model. Validation tests relating stability of a solidifying nucleus to capillary effects and undercooling are also presented. Solidification under various conditions of undercooling and crystal anisotropy is simulated and quantitatively analyzed with particular focus on formation of dendritic microstructure.