## Abstract

We present a multiphase binary alloy phase-field-crystal model. By introducing density difference between solid and liquid into a previous alloy model, this new fusion leads to a practical tool that can be used to investigate formation of defects in late-stage alloy solidification. It is shown that this model can qualitatively capture the liquid pressure drop due to solidification shrinkage in confined geometry. With an inherited gas phase from a previous multiphase model, cavitation of liquid from shrinkage-induced pressure is also included in this framework. As a unique model that has both solute concentration and pressure-induced liquid cavitation, it also captures a modified Scheil–Gulliver-type segregation behaviour due to cavitation. Simulation of inter-dendritic channel solidification using this model demonstrates a strong cooling rate dependence of the resulting microstructure.

This article is part of the theme issue ‘From atomistic interfaces to dendritic patterns’.

## 1. Introduction

Phase-field-crystal (PFC) methodology has been gaining popularity due to its unique capability in modelling crystalline materials [1–4]. In the past decade, it has been used to investigate the physical mechanisms of many interesting problems such as solid-state structural transformation [5–7], crystal and gas phase nucleation [8,9], dislocation dynamics [10], interface instability [11,12] and electromigration [13]. The methodology has also been extended to binary and multicomponent alloys [7], multiple crystal structures [5] and multiphase systems [9].

This work further extends the PFC phenomenology by introducing a binary PFC model that supports a low-density vapour phase. This is done by fusing features of two previous PFC models, the binary PFC of Elder *et al.* [2] and the PFC model of pure materials of Kocher *et al.* [9]. The introduction of a vapour phase into solidification modelling is relevant to processes ranging from solidification shrinkage and cavitation to CVD growth. As such, the model formulated here is demonstrated in the context of phenomena relevant to late-stage solidification.

Solidification of alloy has been a topic of both theoretical and practical importance for many years [14,15]. Computational modelling of crystallization microstructure evolution has led to significant progress in our understanding of solidification [16,17]. At low solid volume fraction, classic solidification theories have been developed to predict the properties of dendritic patterns, whose fundamental length scale is tied to the Mullins–Sekerka instability wavelength [14], which is on the order of micrometres. At very high solid volume fraction (99% and above), the size of interdendritic/intergranular liquid films are narrowed down to the order of nanometers, which is comparable to the solid–liquid interface width. This occurs during the late stage of solidification, where liquid transport is gradually shut down by increasing solid confinement. This regime is a transition from a diffusion controlled process at low solid fraction to a surface energy-controlled process at very high solid fraction. Theories based on interface disjoining potentials have been developed recently to explain grain coalescence and bridging of grains in the last stages of solidification [18–20]. Some of these are based on the PFC methodology, which operates at 10–100 nm, the length scales required to investigate the late-stage solidification process.

Another important behaviour of late stage solidification that has long been known is associated with the formation of defects and second-phase particles [15]. Gas-phase defects, such as microporosity and hot tears, in particular, form at this stage and are imprinted into the final solidified microstructure. However, our knowledge on late-stage solidification of low concentration alloys is still very limited despite its strong impact on final material properties [18].

In the classic picture of late-stage solidification, a partially solidified region with a small liquid-volume fraction is called mushy zone and is heuristically considered a ‘semi-solid’ [15]. Many experimental and theoretical studies have been dedicated to understand materials’ behaviour in the mushy zone [21–26]. As there is typically a density difference between the solid and the liquid, solidification in this confined environment has been known to generate local pressure drop due to shrinkage of the solid caused by density change. It has been shown in experiments that the pressure induced from solidification shrinkage may lead to porosity formation, hot tearing failure, deformation of the dendritic network and even recrystallization in solidified structure [27–30]. Niyama’s criterion was developed to explain the formation of gas porosity [31] in narrow liquid channels. Rappaz, Drezet and Gremaud (RDG) developed a two-phase theory that takes into account both elastic and plastic deformation of the solid confining the narrow liquid channels that form hot tears [32]. The nucleation of the gas phase (or cavitation) in these continuum-based models is typically determined by setting up a pressure threshold which can be obtained from atomic-level models [33,34].

Another important phenomenon at the late stages of solidification is the precipitation of non-equilibrium or second phases [35]. Continuous enrichment of the confined liquid melts due to solute rejection at the solidification front (commonly known as microsegregation) leads to the formation of small secondary solid phases below the eutectic point [36–38]. Materials’ strength can be strongly affected by the size and distribution of those small particles [39,40]. Based on local solute partitioning, the Scheil–Gulliver (SG) equation has been widely used to explain the solute enrichment in eutectic liquid at the late stage.

A comprehensive microscopic continuum model for late-stage solidification is presently lacking because it requires a model that is capable of capturing solute segregation, solidification shrinkage, gas-phase nucleation, and elastic and potential plastic deformation in the solid under confinement [41]. Existing PFC models in the literature have been shown to contain some of these physics separately, but none within a single framework. The focus of this work will be to extend the original PFC phenomenology to consistently capture the aforementioned physics of late-stage solidification in a single theory.

The paper is structured in four sections. Section II presents the model formulation, and some of its alloy and multiphase equilibrium properties are demonstrated. Dynamical evolution for a multiphase alloy system are then discussed. In §3, the model is shown to reproduce typical shrinkage-induced pressure drop in the liquid near the solidification front. A modification of the late-stage SG-type microsegregation behaviour is then predicted due to gas-phase nucleation. Following this, solidification of interdendritic liquid at different cooling rates is also discussed, and a new effect of a gas-induced mechanism for precipitation is presented. The last section concludes the work and discusses future outlook on the use of the model.

## 2. Model formulation

This section begins by introducing a minimal PFC model that can be used to study cavitation in alloy solidification; we begin from previous alloy models [2,7] and a recently developed gas-phase model [9]. The parameters in this new model can be traced back to either of the aforementioned alloy models or the gas-phase model. Using basic thermodynamics, classic binary alloy behaviour is demonstrated in the model. Transitions from high-density phases (*α*, *β* solid and liquid) to a low-density gas phase, and corresponding density change during solidification, are then demonstrated. Dynamics for the microscopic PFC density and concentration are proposed, with the dynamics of the latter field taking into account inhomogeneity in the average density.

### (a) Free energy functional

The dimensionless free energy functional is written in terms of a dimensionless PFC density field , where is the reference fluid density around which the theory is expanded, and a concentration filed *c*(**r**) (as in previous alloy models [2,7]). Its form is the sum of an ideal contribution *F*_{ideal}, the pair correlation contribution *F*_{pair} and the multi-point correlation contribution *F*_{MP}, i.e.
2.1The ideal part of the free energy is expanded as
2.2where *ϵ*, *W*_{1} and *W*_{2} are constants of the model. Here, the spatial coordinate is scaled by a reference solid-phase lattice constant. The pair correlations term has three parts
2.3where *F*_{nn} denotes density correlations, *F*_{cc} correlations in concentration and *F*_{nc} density–concentration correlations. The term *F*_{nn} is written as
2.4For the density pair correlation functions in the above equation, we consider here single peak approximations, i.e. and , where *q*_{1} and *q*_{2} set the first peak positions of the crystal structures corresponding to the limits *c*=0 and and *c*=1, respectively, in reciprocal space. The function *ξ* is smooth and interpolates crystal structure between *c*=0 and *c*=1, as in Greenwood *et al.* [7]. The solid-phase lattice constant *a* is given by . The *F*_{nc} term is a modification of a term in Koher *et al.* [9] given by
2.5where *ζ*(*c*) is a function specified below and the density–concentration pair correlation is approximated by *χ*_{c}(**k**)=*e*^{−k2/2λc} in reciprocal space, where λ_{c} is a constant. The *F*_{cc} term follows the standard form
2.6where *α* is a constant. The contribution from multi-point correlations follows Kocher *et al.* [9],
2.7where *χ*_{3}=(*a*_{1}*ϵ*+*b*_{1})*χ*(**r**_{1}−**r**_{2})*χ*(**r**_{1}−**r**_{3}) and *χ*_{4}=*c*_{1}*χ*(**r**_{1}−**r**_{2})*χ*(**r**_{1}−**r**_{3})*χ*(**r**_{1}−**r**_{4}), with *χ*(**k**)=*e*^{−k2/2λ} in reciprocal space. The parameters λ,*a*_{1},*b*_{1} and *c*_{1} are constants. Assuming the concentration field *c* changes slowly on the scale of the atomically varying density *n*(**r**), equation (2.5) is approximated by replacing *c*(**r**_{2}) with *c*(**r**_{1}) and choosing
2.8where *u*_{1}, *v*_{1} are constants and *c*_{0} is a reference concentration.

### (b) Thermodynamic properties

The terms comprising equation (2.1) are parametrized by the set {*ϵ*,*W*_{1},*W*_{2},*B*_{x1},*B*_{x2},*q*_{1},*q*_{2},*u*_{1},*v*_{1},*c*_{0},λ_{c},λ,*a*_{1},*b*_{1},*c*_{1}}. These constants can be determined by matching the model’s behaviour to the thermodynamic and material properties of an experimental system. This can be done, for instance, by following an analogous approach to that of Asadi & Zaeem [42], which matches the parameters of a PFC model of a pure material to the elastic constants, solid–liquid coexistence densities, density expansion and compressibility of Fe. For an alloy system, we require the additional properties, including the determination of density expansion and coexistence densities associated with the liquid–vapour and solid–vapour phase changes, all of which vary with *c*. These can be uniquely determined by exploiting the larger parameter set associated with the present model. This endeavour requires considerable effort to implement and will be the topic of a future study and publication. In this work, we will restrict ourselves to selecting parameters in ranges where the desired physical properties of a 2-solid–liquid–vapour alloy system are qualitatively demonstrated.

We begin by choosing the model parameters to reproduce the thermodynamic features of a typical binary alloy system at different values of the temperature parameter *ϵ*. The model parameters used in this work are shown in the caption of figure 1. The free energy plot in figure 1*a* shows two separate solid–liquid coexistence regions as a function of concentration. With the temperature parameter *ϵ* sitting above the eutectic point, there is no overlap between the coexistence regions along the concentration axis. At a lower temperature, the two coexistence regions start to merge, as shown in figure 1*b*. Below the eutectic point, there is a large coexistence region between alpha and beta solids, as shown in figure 1*c*.

A crucial physical mechanism that has long been accepted to be responsible for many late-stage solidification phenomena is the change of density associated with solid–liquid phase transition. As expected, the current model parameters are chosen such that both *α* and *β* solids have a higher density than the liquid phase, for a given concentration. This behaviour of the model is shown in figure 2.

To study cavitation due to solidification shrinkage, one of the fundamental causes of subsequent cracking, the model’s parameter set must be in a range that allows the model to simultaneously capture both solid–vapour and liquid–vapour phase transitions in an alloy system. The vapour phase is obtained by introducing a low-density minimum in the uniform branch of the PFC free energy using the technique in [9]. As the model is intended to be used for solidification, and it is known that nucleation of the gas phase during solidification is related to shrinkage-induced pressure drop in the liquid (demonstrated in a later section), the uniform branch of the free energy must have a liquid–vapour coexistence region with dropping temperature parameter *ϵ*. This property is shown in figure 3.

While the shrinkage-induced pressure in solidification may lead to liquid cavitation as described above, it can also lead to plastic deformation, voids and cracks in the solid state. This is particularly true in rapid solidification. Therefore, it is necessary to correctly model the thermodynamics that capture the coexistence of *α* and *β* solids with a low-density vapour phase below the eutectic point. This physical property of the model is shown in figure 4.

The present alloy PFC model is specifically constructed to augment the physical properties of the original alloy model of Elder *et al.* [2] with the demonstrated vapour-phase thermodynamics above. The properties incorporated into the alloy PFC model allow for a vapour phase that can coexist with both liquid and solids, and thus can serve to study the mechanism of cavitation/voiding that is induced by solidification shrinkage. The parameter range in which the above thermodynamic behaviour is exhibited was observed to be robust enough to be able to also capture other interesting phenomena.

### (c) Evolution equations

Dynamical evolution of PFC density *n* follows conserved variational dynamics given by
2.9while the concentration *c* evolves according to
2.10The corresponding mobilities *M*_{n} and *M*_{c} control the diffusive timescale of the *c* and *n* fields, respectively. The variable is the locally averaged spatial density, and the stochastic noise *η* satisfies the fluctuation–dissipation theorem , where is the reference density about which the PFC model is expanded. As *c* is not, in general, a conserved field as in previous PFC models, equation (2.10) approximates *c* dynamics by weighing changes in solute concentration with the local average density. It is clear that this generalized concentration dynamics for spatially inhomogeneous is reduced to the classic concentration dynamics when is assumed to be constant in the system as in previous PFC literature. It has been shown that nucleation-related critical fluctuations are essentially on the interface length scale, therefore, we only included a noise term in the atomically resolved PFC order parameter field *n* and ignored fluctuations in the slow varying concentration field *c*.

Numerical implementation of the dynamics in equation (2.9) and (2.10) is based on Fourier space semi-implicit algorithm. Details of the algorithm is given in appendix A.

## 3. Demonstrations and predictions of the model

This section begins by benchmarking our PFC model against a well-known phenomenon associated with the late-stage solidification processes. Namely, we show that the model qualitatively captures the pressure drop induced by solidification shrinkage inside an inter-dendritic liquid channel. Then, our model is used to predict how the SG microsegregation relation is modified in the presence of a gas phase. Finally, we examine cooling-rate-dependent microstructures that demonstrate a mechanism for void formation due to the appearance of second-phase precipitates.

### (a) Shrinkage-induced pressure

At the late stage, solidification proceeds mainly within the inter-dendritic liquid and separated eutectic liquid pools. Because of the large solid fraction in this regime, transportation of liquid is very limited. Local mass shortage due to solidification shrinkage cannot be effectively compensated by liquid transport as in traditional solidification processes. A pressure drop in the liquid appears then as the result of local mass shortage.

A basic theory in the literature that describes the above phenomenon in directional solidification is based on an effective theory that considers the mushy zone as an effective one-dimensional inter-dendritic liquid channel that terminates in a liquid reservoir. It assumes that liquid flow within the channel is driven only by a pressure gradient (Darcy’s Law). The pressure drop in the liquid at distance *x*=*d* from the reservoir is given by [15] where *g*(*x*)=*μf*_{l}(*x*)/*K*(*x*) and where *μ* is the liquid viscosity, *K* is the mushy zone permeability and *f*_{l} is the local liquid volume fraction. The velocity of the liquid normal to the solidification front *v*_{ln} is given by , where *v*_{n} is the velocity of the solid–liquid front and the shrinkage factor where *ρ*_{S} is the solid density and *ρ*_{L} is the liquid density. The boundary condition is , with *f*^{D}_{l} the upper limit of the liquid fraction where the Darcy-type flow becomes dominant. This model has been extended to incorporate elastic and plastic deformations of the solid confinement by Rappaz *et al.* [32]. The salient physics of this model involves the balance of mass flux due to liquid transport and solidification consumption. As liquid flows through microstructure with increasing solid volume fraction, a larger liquid surface area is exposed to the solidification front and more liquid mass is drawn into the solid. Meanwhile, the liquid transport becomes less effective within narrowing liquid channels, which leads to a pressure drop within the channels that cannot be fed by liquid.

The model in the previous paragraph is for pure materials and considers fluid flow. Moreover, it considers mass conservation across sharp interfaces. Nevertheless, its mechanism for pressure drop due to density change should be contained with our PFC model. As the length scale of the PFC model is limited by its atomic density resolution, we test shrinkage-induced pressure drop in a single groove geometry (figure 5) where the solid volume fraction also goes from 0 to 1, instead of a large typical mushy zone microstructure. The essential physics of mass flux balance is still retained in this simple set-up. Namely, as the PFC model uses a conserved ‘density’ field to describe the solidification, we expect that it qualitatively captures pressure drop due to solidification shrinkage. We examine the pressure drop in the test groove geometry shown in figure 5 using the PFC model formulated in equations (2.1)–(2.7) with the dynamics of equations (2.9)–(2.10). A typical pressure drop curve is shown in the bottom frame of figure 5. This resembles the same qualitative form seen in other theories. Simulations are performed in a 4096 by 512 grid box with grid spacing and Δ*y*=*π*/(4*q*_{1}). Pressure in the liquid is measured along the centre line shown in the top frame of figure 5. As shown in the previous thermodynamics section, this model allows a low-density vapour phase to form in eutectic liquid. Thus, if the pressure drop in the liquid channel becomes large enough, it may lead to liquid cavitation. We shall examine this topic below.

### (b) Solute segregation

It is well known that microsegragation of solute atoms has a strong impact on defect formation and final solidification microstructure. The SG relation has been widely used to describe microsegragation and the formation of metastable solids at high solid volume fraction [15]. To increase the solid volume fraction by df_{s}, the solidification front has to reject solute (*c*_{l}−*c*_{s})df_{s} with *c*_{l}(*c*_{s}) the local liquid(solid) concentrations. The remaining liquid is then enriched due to this rejected solute by an amount dc_{l}=(*c*_{l}−*c*_{s})*df*_{s}/*f*_{l}, assuming the rejected solute atoms are evenly redistributed within the liquid. With *f*_{s}+*f*_{l}=1, partition coefficient *k*=*c*_{s}/*c*_{l} and boundary condition *c*_{l}=*c*_{0} at *f*_{s}=0, the average liquid concentration *c*_{L} can be expressed as a function of the final solid fraction *f*_{S} as
3.1However, equation (3.1) does not take into account the effect of shrinkage-induced pressure drop that occurs in the liquid at the late stage of solidification. Specifically, since a large pressure drop may lead to gas phase nucleation in the interdendritic liquid. A crude modification of equation (3.1) that considers a gas phase fraction *f*_{g} can be written as
3.2where *c**_{L} and *f**_{S} are corresponding liquid concentration and solid fraction before gas phase nucleation. Because the rejected solute has to be distributed within a smaller liquid volume following gas nucleation, equation (3.2) predicts a potentially higher liquid concentration for the same solid volume fraction at the late stage of solidification. Numerical results from our PFC model qualitatively demonstrate such a behaviour in figure 6. The simulation is carried out in a quasi-one-dimensional system of 16 384 by 16 grid (one unit cell along the small dimension) where a nano-size liquid channel is enclosed by two solid slabs. To our knowledge, this is the first model that qualitatively captures this enhanced solute segregation due to gas-phase nucleation. This can lead to interesting solidification structures, as shown in the next section and Wang *et al.* [43].

### (c) Solidification of liquid channel

Following the demonstration of pressure change and void formation on segregation, the PFC model is pushed a little further to examine second-phase formation in late-stage solidification of a confined inter-dendritic liquid. Instead of simulating a long inter-dendritic liquid channel, we choose to go with a solid–liquid–solid sandwich configuration which can be seen as a slab of a long channel. Starting from the top configuration in figure 7, the resulting microstructure can be vastly different depending on the cooling rate, as shown in the two middle configurations in figure 7. For the slow cooling rate, the microstructure at the end is still sandwich-like, with a stretched liquid (SL) region with enriched solute concentration similar to SG segregation. For the fast cooling rate, both gas phase (labelled as G) and *β* solid are observed at the end of solidification. Detailed investigation of the solidification process in this case reveals that the liquid first quickly solidifies into metastable solid precipitates, which produces a large shrinkage-induced pressure that triggers gas-phase nucleation between these precipitates. The metastable solid then decomposes to form *α* and *β* solids. Although the model is only qualitative at this point, the interplay of concentration and mass transport is shown to be an important factor in late-stage solidification.

## 4. Conclusion

A novel alloy PFC model that incorporates solid–liquid and gas-phase interactions was developed. It is the first step towards a self-consistent phase field formalism for investigating elasto-plastic effects in the late stages of rapid solidification. Here, the model is shown to qualitatively reproduce shrinkage-induced pressure drop and cavitation in inter-dendritic liquid. To our knowledge, this is the first phase field/PFC type of model that captures this important physics in solidification.

As a multiphase alloy model, it also captures an interesting microsegregation behaviour in eutectic liquid when cavitation is triggered by shrinkage-induced pressure. To our knowledge, this coupled evolution of solute concentration with shrinkage-induced pressure/cavitation in late-stage solidification has not been examined in previous experiments or theories.

Cooling rate-dependent solidification behaviours in interdendritic liquid were examined. In the case of fast cooling, the model predicted that fast solidification of eutectic liquid into metastable solids leads to void formation (which promotes solidification cracking), while the decomposition of the metastable solids leads to *β* and *α* solids. Although the parametrization of the model is only qualitative at present, the interesting phenomena it demonstrates in this work suggested that further investigation along the lines shown here may lead to new understanding of late-stage solidification processes.

The multiphase binary alloy PFC model developed in this work is also expected to be easily generalized to study the combined effect of segregation and defect nucleation in other multiphase systems. Furthermore, its analysis using coarse-graining schemes can be used to guide the design of a complex amplitude multi-phase phase-field model for use with higher-scale simulations.

## Data accessibility

No experiments were performed in this work. Numerical schemes used to obtain results are described in the main text.

## Authors' contributions

N.W. and N.P. conceived of and designed the study. N.W. performed all calculations and drafted the manuscript. G.K. contributed to the idea that leads to gas phases in the PFC model. N.P. revised the manuscript. All the authors read and approved the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

The article is funded by National Science and Engineering Research Council of Canada, the Canada Research Chairs Program and Compute Canada.

## Acknowledgements

We acknowledge Nathan Smith for helpful discussions early in this work.

## Appendix A. Numerics

Implementation of the evolution equations is based on the Fourier space semi-implicit algorithm. The equation of motion for the PFC density field is relatively straightforward because the mobility parameter is assumed to be a constant. The equation for the concentration field is slightly different from the density field equation because there is a density-dependent prefactor , where is a local average of *n* (defined in the text). Implementing of a semi-implicit algorithm for the concentration proceeds as follows: by dividing on both sides of equation (2.10), the concentration equation can be rewritten as
A 1Retaining only the second term on the right-hand side of equation (4.1) reduces it to the standard PFC concentration dynamics which can be solved directly with the standard semi-implicit algorithm [3]. To evaluate the full equation, the 1st term on the right-hand side is first calculated in real space and then treated as a nonlinear term that is added to the second term in the standard semi-implicit algorithm for the PFC model. This approach also makes it possible to model two-sided solute diffusion in a straightforward way, although that is not considered here.

## Footnotes

One contribution of 16 to a theme issue ‘From atomistic interfaces to dendritic patterns’.

- Accepted October 4, 2017.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.