Modeling of coupled motion and growth interaction of equiaxed dendritic crystals in a binary alloy during solidification

Motion of growing dendrites is a common phenomenon during solidification but often neglected in numerical simulations because of the complicate underlying multiphysics. Here a phase-field model incorporating dendrite-melt two-phase flow is proposed for simulating the dynamically interacted process. The proposed model circumvents complexity to resolve dendritic growth, natural convection and solid motion simultaneously. Simulations are performed for single and multiple dendritic growth of an Al-based alloy in a gravity environment. Computing results of an isolated dendrite settling down in the convective supersaturated melt shows that solid motion is able to overwhelm solutal convection and causes a rather different growth morphology from the stationary dendrite that considers natural convection alone. The simulated tip growth dynamics are correlated with a modified boundary layer model in the presence of melt flow, which well accounts for the variation of tip velocity with flow direction. Polycrystalline simulations reveal that the motion of dendrites accelerates the occurrence of growth impingement which causes the behaviors of multiple dendrites are distinct from that of single dendrite, including growth dynamics, morphology evolution and movement path. These polycrystalline simulations provide a primary understanding of the sedimentation of crystals and resulting chemical homogeneity in industrial ingots.


Rising of a circular particle in bulk liquid
Before applying the proposed phase-field model to the case of moving dendrites, it is necessary to validate the model in handling solid-liquid flow in some simple situations that can be benchmarked using analytical solutions or reported numerical simulations. The convergence of the model is studied by simulating a freely floating solid particle in two dimensions. The particle is assumed to be circular with radius a equal to 400d0 (here, d0 = 1.96×10 -6 cm). The size of computational domain is 250a × 250a, which is to guarantee the particle moving freely in bulk liquid and not affected by the wall. The densities of two phases are s = 2.45 g/cm 3 and l = 4.9 g/cm 3 , respectively. Given that / 2 l s    and density variation is not sufficiently small, the averaged density (l + s)/2, instead of pure liquid density, is selected as the background density in the Boussinesq approximation 1 . The kinematic viscosities of liquid and solid are chosen to be 5.0 × 10 -3 cm 2 /s and 5.0 × 10 3 cm 2 /s, respectively.
The value of W0 ranges from 5d0 to 40d0. The simulated rising particle with W0 = 10d0 is presented in Fig. S1. From the distribution of phase-field variable evolving with time in Fig. S1(a), it can be seen that though the particle is not assumed to be rigid, it always keeps its initial shape and size during rising, as well no stretch and deformation are observed. The closer flow field depicted in Fig. S1(b) shows each point in the particle region moves upwards at a same speed. Two vortexes generate symmetrically on both sides of the particle as a result of flotation in Fig. S1(c, d).
These features of the moving particle imply that the proposed model is able to correctly simulate the solid particle motion in the liquid.  Fig. S2. Initially, the particle rises slowly but with the highest acceleration. Due to the viscous friction that is proportional to the square of the rising velocity, the acceleration decreases with the increment of velocity. Eventually, the particle reaches steady state where the driving 0.12s 0.06s force is balanced by frictional force. For all the choices of W0, the rising velocities evolve with the same trend but just differ at the steady state. The drag force acting upon a circular particle by the liquid has been proposed by Van Dyke 2 where CD is drag coefficient, = 0.5772 is Euler constant, V is rising velocity of the particle,  is kinematic viscosity of the fluid. Terminal velocity is achieved when the sum of gravity and drag force is equal in magnitude to buoyance, and yields where g = 980 cm/s 2 is gravitational acceleration. After substituting the parameters of two phases into Eq. (2), Vinf = 0.072 cm/s is calculated through an iterative algorithm.
Comparison of phase-field simulated floating velocity of a particle with the analytical solution is illustrated in Fig. S2. As the interface width parameter used in the simulation decreases, the prediction by the phase-field model shows convergence and approaches the analytical solution. In following simulations, the interface width parameter W0 is chosen to be 10d0 in consideration of both numerical accuracy and computational efficiency.

Rotation of a growing dendrite at a constant angular velocity
The rotation of a fourfold symmetry dendrite growing from undercooled melt is simulated, aiming to test the advection and curl terms in the crystallographic orientation equation. All of the computational parameters for dendrite growth are identical to that proposed by Karma 3 . An angular velocity which is set to be π/(100000) is imposed over the entire computational domain. Initially, the orientation of the dendrite is 0°. After a 45° rotation under the given angular velocity, the simulated dendritic morphology and the imposed rotational velocity field are presented in Fig. S3(a). Generally, the dendrite keeps the fourfold symmetry, as it should be. Nevertheless, when the crystallographic orientation change resulting from rotation is canceled in the phase-field equation, a distorted dendrite is obtained in Fig.   S3(b), which is similar to the swirling dendrite obtained by Yamaguchi 4 0.072 targeted angle calculated from the prescribed angular velocity are plotted for comparison in Fig. S3(c). A wonderful agreement is achieved, which claims the accuracy of the proposed model in handling the crystallographic orientation evolution.
Further quantitative comparison of the tip growth velocity with pure diffusion-controlled case and that by Karma 3 is plotted in Fig. S3(d). The time dependence of the grain orientation is shown in Fig. S5. It exhibits an approximate linear reduction, which means that the dendrite rotates clockwise in a constant angular velocity.