Transient Effect of Fluid Flow on Dendrite Growth Direction in Binary Fe-C Alloys Using Phase Field in OpenFOAM

An open-source fluid flow and phase-field coupled solidification model has been developed in OpenFOAM to investigate the transient nature of the interface growth direction at different degrees of undercooling for an Fe-0.15 wt.% C binary alloy under isothermal conditions. Though there are works on melt convection effects in binary alloys, none reported the transient nature of the dendrite growth direction since thermodynamic driving force decreases with time at a particular undercooling. Developing a theoretical relation will be helpful in understanding the competition between the crystallographic growth direction and solute transport. Flow decoupled simulation results have a good quantitative agreement with the literature. The bending angle formulations on the effects of flow velocity and growth speed were separated. At the end, improved theoretical formulations for estimation of the bending angle based on the anisotropy in interface energy were put forward compared with only few available empirical correlations.


INTRODUCTION
Steel is one of the most widely used materials in today's world, manufactured through the continuous casting route. The microstructure that results is important for upstream processes and final product quality. The practice of continuously cast near-netshape profiles has gained in popularity because of energy savings, but these limit the subsequent inline thermo-mechanical processing possibilities; therefore, the need to control the cast microstructure becomes vital. The solidified microstructure depends on the alloy chemistry and those conditions in the mould that control the thermal profile and fluid flow. To understand the links and identify the operating conditions that control the evolution of the solidification structure, computational models offer ways to simulate the process and identify sensitivities to distinct conditions. Phase-field 1-3 and cellular automaton 4 are currently used to simulate the dependence of the dendritic solidification microstructure on process conditions. Based on the continuum thermodynamic principles, the phase-field method has received much consideration for simulating the complicated interface structure because of its implicit nature of tracking the solid/liquid front via the diffuse interface approach. The real essence of the method stems its origin from the well-known Cahn-Hilliard equation. 5 Wang et al. 6 deduced phase-field models that had a strong resemblance to that of Kobayashi 7 . Starting from a pure metal, Kim et al. 1,2,8 and Boettinger et al. 9 deduced the phase-field model for alloy solidification in a thermodynamically consistent way. Kim in his model 1,2 of binary alloy solidification assumed the interface to be a mixture of solid and liquid phases with different compositions but with same chemical potential. The model 1 under the vanishing interface kinetics coefficient condition maintains local equilibrium at the interface.
The thermodynamics of excess solute rejection during the alloy solidification governs the nonuniformity of the chemical composition in the cast product commonly known as segregation. Segregation near a growing front will alter the solidification kinetics and have serious consequences on the final mechanical properties. 10 Bulk fluid flow during alloy solidification washes the rejected solute away (thereby changing the local composition) and will thus define the segregation and consequently the dendrite growth direction. 11 Ultimately, this influences the evolution of the solidification structure. The solute layer at the dendrite tip, if driven by bulk flow over long distances, may give rise to segregation at the macro-scale level.
The phenomenon of dendrite deflection is defined as the biased growth of the solid phase towards the upstream flow direction. Investigations of the modelling of the melt convection effect on the solidification microstructure [11][12][13][14][15][16][17] have been focussed on the effect of inlet flow magnitude on increased dendrite bending and branching of dendrites in the upstream direction, which is primarily due to the asymmetric solute profile ahead of the deflected dendrite. Not much information is available on the effect of undercooling on dendrite deflection. Most of the research on melt flow has been focussed on nonferrous systems such as Ni-Cu, 14,15 Al-Si 18 and Al-Cu 12,13 while a few exist for Fe-based alloys such as Fe-C 19 and Fe-Mn, 20 which can improve the understanding of the microstructure formation under industrial conditions. To the best of our knowledge, the quantitative dependence of the growth direction on the magnitude of the flow velocity, solute content and solidification speed is limited to two previously reported studies 11,16 -one showing the empirical dependency of dendrite deflection and in the other study anisotropy in interface energy was not incorporated. Okano et al., 21 also based on limited experiments in a continuously cast steel slab, put forward quantitative relations of the bending angle. For an Fe-0.15 wt.% C alloy at relatively high undercooling conditions, using the phase-field method, Natsume et al. 19 investigated the singledendrite deflection at three different growth speeds only by considering anisotropy in interface energy. During the progress of casting, levels of undercooling at various positions of the solidification front within the melt change with time depending on the process conditions. This may give rise to different deflection behaviours in combination with the incoming fluid flow. Developing theoretical correlations will help in understanding this behaviour. For pure materials, a few reported 22 3D simulations showed that the 3D flow effect is more pronounced than the 2D but did not incorporate the growth direction effect.
In the commercially available CFD software, it is not easy to modify the available numerical models. OpenFOAM 23 is an open source computational software, capable of handling a broad spectrum of partial differential equations along with the advantage of having a set of precompiled libraries that can be customised by the user as per the requirements. Despite the wealth of publications 24,25 on CFD, a limited amount of literature can be found wherein coupling of OpenFOAM to phase-field methods has been carried out. 26,27 Thereby, developing a phase-field method based on an open-source numerical  model will not only be useful in understanding the  solute advancement behaviour due to fluid convection in Fe-C binary alloys but also can be helpful for  the material research community in the development of new

Physical Fundamentals
The dendritic solidification structure is the result of competition between the interface energy and thermodynamic driving force due to undercooling. The phase-field variable / gradually changes its value from 0 in the fully liquid phase to 1 in the fully solid phase across a distance of finite value called the interface thickness (2k). It is to be noted that this interface is not a physical one present in the actual solidification microstructure but a mathematical construction for solving the governing equations. 1 For the case of solidification, the freeenergy functional of the system can be represented as a sum of the free energy of the bulk phases and that of the excess energy associated with the creation of the interface. The free-energy functional is assumed to be a continuous function of the phasefield variable across the diffuse interface. The interface energy 9 constitutes the energy barrier at the interface (where w is the barrier height) and gradient energy (where gradient energy coefficient is denoted by e) associated with the gradient of the phase-field variable across the interface. These are two counteracting effects-the gradient energy tries to widen the interface region to reduce the energy while the barrier height ''w'' tries to make the bulk phases stable by sharpening the interface region. This interface energy vanishes in the bulk phases and is non-zero in the interface region only. It is to be noted that both w and e are computational parameters of the phase-field model. For isotropic interface energy, the equilibrium shape is a sphere. Formation of solid crystals during solidification takes place in a faceted manner. Solidification of cubic symmetry crystals in dendritic mode has been Transient Effect of Fluid Flow on Dendrite Growth Direction in Binary Fe-C Alloys Using Phase Field in OpenFOAM assumed to take place in a weakly anisotropic 28 manner at high undercooling. The anisotropic nature is taken into account by assuming the interface energy ðrÞ to be dependent on the angle h between the interface normal and a reference crystalline axis. This anisotropy dictates the preferred growth direction for the dendrites to grow. The governing equation for evolution of the solid phase with time under dilute solution approximation 1 used in the present study is given as where M is the phase-field mobility, h is the orientation of the normal to the solid/liquid interface with respect to the positive x-axis, T is the temperature, R is the universal gas constant, k e is the equilibrium partition coefficient of carbon in diron, T m is the melting point of pure iron, v m is the molar volume, m e is the equilibrium liquidus slope, C L is the solute concentration in the liquid phase in mole fraction, hð/Þ is a monotonically increasing interpolating function, and wg 0 ð/Þ is the imposed double-well potential. 2 The term eðhÞ is associated with the gradient of the phase-field variable, i.e., r/, as mentioned before. The above governing equation shows that the interface movement is the interplay between smoothing by the diffusion term (which tends to make the interface diffuse to reduce gradient energy) against the thermodynamic driving force (which tends to make the interface sharp to minimise the excess material within the interface and stabilise the bulk phases). Thus, the most widely used method in phase-field models to include the anisotropy is to assume e (gradient energy coefficient) as a function of h represented as eðhÞ. For materials with cubic symmetry undergoing dendritic solidification, the weakly anisotropic nature of the interface energy in 2D is taken into account by assuming the following dependency of e as 19 where d M is the anisotropy constant and j is the mode of anisotropy taken as 0.03 and 4 (fourfold symmetry), respectively. d M determines the strength of anisotropy and e is the gradient energy coefficient. The above equation for anisotropy is not valid for highly faceted solidification, i.e., for materials having strong anisotropy. Under the condition of local thermodynamic equilibrium at the interface and a thin interface limit, the model parameters w, e are related to the material properties, interface energy (r 0 ) and interface width 2k 0 (where / changes from 0.05 to 0.95) as 1,8 It is to be noted that r 0 is the orientationally averaged material interface energy. Thus, since e is related to the material interface energy (Eq. 3), the assumed dependency of e on h also relates the material interface energy to h. The material properties have been assumed to be independent of composition.
The polynomial functions for hð/Þ and gð/Þ are given as 1 Since the simulations were performed isothermally, the governing equation for solute transport by both convection and diffusion for binary alloy solidification is given as 1,29 where Dð/Þ is the diffusivity of the solute defined as Dð/Þ ¼ D S / þ ð1 À /ÞD L with D S and D L (Table I) being the diffusivity coefficients in solid and liquid respectively and V is the liquid velocity. The mixture solute composition C is defined as the weighted average of the C L and C S (the solute concentration in the solid phase) as Each point within the interface is assumed to be a mixture of solid and liquid phases with different compositions of C S and C L . Local equilibrium has been assumed at the interface.
The phase-field (interface) mobility M in the thininterface limit with vanishing kinetic coefficient for binary alloy is defined as 2, For incompressible fluid flow, Eq. 10 represents the conservation of mass assuming solid velocity to be zero r Á ð1 À /ÞV ¼ 0 ð10Þ and that for the conservation of momentum 29 is given as where p is the pressure, m is the kinematic viscosity of liquid steel, and d is the interface thickness.

Numerical Methodology
The code developed to solve the governing equations using the finite volume method in OpenFOAM was written in C++ programming language. The phase-field and solute transport equations were solved using an explicit Euler scheme and the governing equations for fluid flow were solved implicitly with the standard PISO algorithm. 23 A random noise of the form 16ar/ 2 ð1 À /Þ 2 of 1% amplitude was added to the liquid concentration at the interface to simulate the formation of secondary arms where r is the random number between + 1 and À 1 and a is the noise amplitude. 30 The width of the domain was taken to about 11 times the length of the reported secondary dendrite arm spacing and the domain height was taken to be around 4.5 times the secondary arm spacing. A mesh size of 0.5 lm was used to calculate the dendrite arm spacing. Noflux boundary conditions for both the phase-field and concentration were used. As an initial condition, a uniform thin layer of solid was placed at the bottom of the domain. The effect of cooling rate was implemented in form of the equation TðtÞ ¼ T 0 À tDT 3 where T 0 is the initial temperature a few degrees below the liquidus temperature of the alloy and DT is the cooling rate. The initial temperature T 0 was constant for all the cases. Growth of secondary dendrite arms was simulated at three different cooling rates, namely-83.33 K/s, 33.33 K/ s and 8.33 K/s. Two-dimensional simulations were performed for Fe-C binary alloys under isothermal conditions with the thermodynamic data taken from Thermo-Calc software. 31 Simulations for solidification coupled with fluid flow were performed under different levels of undercooling. The mesh size was composed of 1000 9 500 cells with an equal distance of 10 À8 m (= 10 nm). The seed crystal (initial condition) was placed at the centre of the bottom wall of the domain. Liquid melt at a constant inlet velocity along the x-axis enters the left wall (inlet) of the domain and exits at the right wall of the domain with zero-gradient velocity boundary condition. The top wall of the domain was assumed to be a slip wall and the bottom wall was assumed to be in contact with stationary mould wall. The uniform value of zero pressure was taken at the right wall and the rest of the walls were taken as zero-gradient pressure boundary condition. The boundary conditions for the phase-field and concentration were taken as zero gradient except for the left wall of the domain. For isothermal flow-coupled simulations, the k e value was constant up to three significant digits.

Model Validation
The developed numerical code was validated by comparing the results with those of Kim et al. 3 for Fe-0.1 wt.% C alloy. Kim et al. 3 validated his modelled data against experimental reported results. Figure 1 shows the evolution of the secondary dendrite arms for the binary alloy at different time instants and comparison of the secondary dendrite arm spacing data as a function of cooling rate with literature. 3 Figure 1a and b shows the phase-field profile of the growth of secondary dendrite arms and their coarsening behaviour at two different time instants (t 1 and t 2 ), which seems to be in good agreement as shown in Ref. 3. The solid phase is represented by when the phase-field variable takes the value of 1 (color code ''red'') while the liquid phase is represented as the value 0 (colour code ''blue''). The interface region surrounding the growing solid is quite wide because of the coarser mesh size used though there is evidence of entrapped inter-dendritic liquid in between successive secondary arms. It should also be noted that the solute-enriched regions near the bottom still remain liquid because of their reduction in melting point, according to the iron-carbon phase diagram. Some secondary arms grow preferably while some arms stop moving forward thereby giving rise to Table I Figure 1c shows the comparison of the secondary dendrite arm spacing data between the literature and present work as a function of cooling rates. The predicted behaviour of secondary dendrite arm spacing with cooling rate follows the similar quantitative trend as reported. 3 The modelled secondary dendrite arm spacing depends on the thermodynamic parameters used in the system. 3 For validation, both the liquidus slope and the equilibrium compositions along with the phase-field mobility have been put as a function of temperature as opposed to constant values used in Ref. 3. Also, the initial temperature taken in the work is a few degrees below the liquidus temperature. That is why maybe there is a difference in the secondary dendrite arm spacing values shown in Fig. 1c  Fluid Flow Effect Figure 2 shows the effect of fluid flow on the interface growth direction for an isothermal solidification of Fe-0.15 wt.% C binary alloy at an undercooling of 24.3 K. Figure 2a shows the growth of a solid in the form of a dendrite in the absence of fluid flow. The primary arm grows (shown by the black line) in the direction of maximum anisotropy, i.e., perpendicular to the base of the domain. The image is symmetric with respect to the black line. Figure 2d shows the concentration profile along the white dotted line with and without flow. The rejected solute build-up layer ahead of the moving interface is symmetric on both sides of the primary arm. For flow-coupled solidification simulations, the melt inflow velocity is taken as 0.15 m/s. 19 The fluid flow profile (similar to that obtained by Natsume et al. 19 ) at a time instant during the process of solidification is shown in Fig. 2c. The variable 'ULIQ' in Fig. 2c is the product of the liquid phase fraction multiplied by the liquid velocity 'V'. The flow vectors enter from the left wall of the domain. Then, facing the growing solid obstacle, they rise up and after reaching the tip of the primary arm they follow the path towards the exit on the right wall. The magnitude of the flow vectors close to the tip of the primary arm reaches a maximum as the growing tip approaches the upper wall. Figure 2b represents the growth of the solid in presence of fluid flow. The growth direction of the primary arm is shown by the yellow line, which is tilted with respect to the original growth direction in absence of fluid flow. Thus, the growing solid bends towards the upstream direction of fluid flow. The angle between the original growth direction and deviated growth direction is represented as h in subsequent sections (as shown later in the Eq. 12). This is because the incoming fluid on its way towards the exit takes away the rejected solute from the upstream side of the primary arm to the downstream side thereby creating a washing effect as pointed out by previous researchers. 12,32 This washing effect is responsible for asymmetrical solute build-up on either side of the primary dendrite as evident from Fig. 2d and thus defines the growth direction. The increased peak on the downstream side shows the enrichment of the solute, whereas there is depletion of solute on the upstream side. Hence, the interface growth is favored towards the upstream side. The higher the speed of the flow, the higher will be the asymmetry in the concentration profile and hence dendrite deflection angle. Also, the change in the growth direction being an instantaneous phenomenon, it is more relevant to speak of the flow velocity close to the tip rather than the incoming melt velocity. In case of multi-dendrite growth, even if the input melt flow speed is the same, the flow velocity near the dendrite tips will depend on the growth pattern of the multiple dendrites and thus will contribute to the solute profile. During a stable continuous casting operation, the input liquid jet stream will have a stable velocity but the flow velocity in front of the moving solidification front at various depths from the meniscus will be different and thus will contribute differently to the solute profile at various stages of the casting operation. It is therefore imperative to study the influence of the flow velocity near the dendrite tip on the dendrite bending angle rather than the initial melt velocity at the entry point.

Effect of Undercooling
During solidification at a particular undercooling, it is known 15 that the growth speed of the dendrite tip initially remains high and then it gradually decreases with time because of a decrease in thermodynamic driving force, thereby achieving an almost steady-state value. It is therefore quite important to study the time-dependent behaviour of the dendrite growth direction in the presence of fluid flow for a particular undercooling. For a constant inlet melt flow speed of 0.15 m/s, simulations were performed at five different levels of undercooling of 24.3 K, 19.3 K, 14.3 K, 9.3 K and 6.3 K respectively. For each level of undercooling, the dendrite bending angle, the maximum flow velocity at the dendrite tip and the average growth speed of the tip of the primary arm were calculated at different time instants. Figure 3 shows the 2D contour plot of the variation in dendrite bending angle with maximum fluid flow velocity at the tip and the growth speed as the axes for the whole set of simulation data points (shown as symbols). The bending angle is represented by the sine value of the angle (h) instead of the value of the angle itself. This is because the sine value gives the displacement between the tip and seed crystal position in the flow direction. Similar growth speeds at high undercooling have also been reported by Natsume 19 . Two distinct features can be seen from Fig. 3-one being the dependency on the maximum flow velocity and the other being the dependency on the growth  speed. At a constant tip growth speed it can be seen that with an increase in the maximum flow velocity at the tip, the bending angle increases. The initial steep increase in the bending angle decreases afterward and has a tendency to saturate towards the end. This situation may be correlated with the initial period of the continuous casting of steel when the liquid steel undercooling is relatively high. The liquid steel jet tends to wash away the rejected solute from the advancing solidification front as a result of which the bending angle may increase. Also reducing the casting speed, i.e., inlet flow speed will reduce the bending angle.
At a constant maximum flow velocity, it can be seen that with a decrease in the tip growth speed, the bending angle increases. 19 The increase in the bending angle is quite steep at lower growth velocities. This is because the fluid passing the dendrite tip has much more time to sweep away the rejected solute from the upstream side of the dendrite tip to the downstream side and thereby contributes to a higher bending angle. However, in this growth speed-dominant regime, the solute is not washed away to further distances and instead creates a trail of an asymmetrical solute boundary layer. Thus, this might be the point where diffusion tries to gain importance over bulk convection and may contribute to macro-segregation in the cast product. The authors thus have made an attempt to separate out the two effects-the flow velocity effect and growth speed effect-on the bending angle. This information on the bending angle might be useful for casting operators in finding out the linkage between the dendrite growth direction and casting parameters at various stages of the casting and thereby taking corrective action.

Fit Function Analysis
The empirical dependence of deflection angle h for steel ingots put forward by Takahashi et al. 16 is given as where h is in degrees, V is the flow velocity of the bulk liquid in cm/s, and v solidification is the solidification rate in cm/s. The flow velocity close to the moving front can at times be higher than the reported flow velocities of the incoming bulk fluid. The relation shows a logarithmic dependence on solidification rates and liquid flow velocity. Okano's relation 21 has a resemblance to that of Takahashi. For both the relations, the dendrite bending angle is not mathematically defined at extremely low flow velocities and solidification rates because of the logarithmic dependence. Since the relations consider only the flow magnitudes, it is expected that the bending angle will always be positive. However, in some cases, the predicted bending angle may turn out to be negative for the relations. Also, it is often challenging to carry out high-temperature experiments to measure the dendrite bending angle and relate to the industrial fluid flow conditions, which involves cost and effort. Thus, it can be quite useful to have simple bending angle relations with sound theoretical links that can be used to predict the bending angle from a wide range of flow velocities and growth speeds. In future, these relations can also be fine-tuned based on industrial conditions. A theoretical relation developed for one particular alloy system can be extended to other systems-binary (e.g., Fe-Mn, Fe-Ni, Fe-Al, Al-Cu, etc.) or even ternary systems (Fe-C-Mn, Fe-C-Al, etc.). In this section, by revisiting these 16,21 bending angle formulae, the authors have made an attempt to improve the mathematical basis with the obtained phase-field method-based modelling results by incorporating the anisotropy in interface energy and separating the effects of flow velocity and growth speed. Figure 4 below shows the surface fitted (using MATLAB R2018a software) with the data points in the flow velocity-dominated regime and the corresponding fit function (with 95% confidence limit) is given as where V is the maximum fluid velocity at the tip (m/ s) and v tip is the tip growth speed (m/s). The R 2 value for the relation is 0.99. The above relation (Eq. 13) is similar to that of the Takahashi 16 relation (Eq. 12) in the sense that the bending angle is directly proportional to the flow velocity but inversely proportional to the growth speed of the front. The higher exponent value of the flow velocity signifies the dominant effect of flow. The fitted surface shows that the bending angle increases with the flow velocity in an exponential manner and most of the data points lie close to the surface. The fitted surface moves towards a saturation level at very high flow velocities. The fitted equation shows that the bending angle is proportional to 1 À e À0:7551V 0:6 ð Þ , which has a resemblance to that of Takahashi 16 where the bending angle is proportional to logV 2:08 . The exponent of the flow velocity may differ for different alloy systems.
For the second set of data points in the growth speed-dominated regime, using the corresponding flow velocities and growth speeds, the bending angle values were first estimated from the previous fit function and then it was subtracted from the original simulated bending angle values for those points. This was done to remove the base effect of the flow velocity on the bending angle for this set of data points so that the remaining effect would be the contribution from the growth speed. The corresponding fit function for the growth speed-dominated regime after removing the effect of flow velocity (with 95% confidence limit) is given as sinh ¼ 1 À e À0:3505V 0:01 e À8:753v 0:5 tip ð14Þ where V is the maximum fluid velocity at the tip (m/s) and v tip is the tip growth velocity (m/s). The R 2 value for the relation is 0.80. A higher exponent value of tip growth speed signifies the dominating effect of growth speed compared with flow velocity in this regime. The corresponding fitted surface will show the additional effect of growth speed on the bending angle on top of the fluid flow effect. A decrease in undercooling decreases the thermodynamic driving force and hence decreases the growth speed of the interface. Since the growth speed decreases, the fluid flow will have more time to wash away the rejected solute and hence the bending will increase with a decrease in undercooling. This particular fit function can be useful in predicting the bending angle in low growth speed regimes. In this case, also most of the data points lie close to the surface. At very low growth speed, the bulk solute transport will be dominated by the diffusion, which gives rise to high solute build-up ahead of the interface. Growth speed tending to very high values will give a very low bending angle, i.e., the interface will tend to grow as if it is a case of solidification without flow. Thus, the combined fit function for the whole set of data points can be written as sinh ¼ 1 À e À0:7551V 0:6 e À1:851v 0:01 tip þ 1 À e À0:3505V 0:01 e À8:753v 0:5 tip ð15Þ Thus, the above equation shows that the total bending angle is the contribution from the flow velocity along with an additional contribution from growth speed. The above fit equation is mathematically defined for limiting cases, i.e., cases of extremely high velocity or growth speed as well as extremely low velocity and growth speed. Obviously at zero velocity, the equation gives zero bending angle. Figure 5 shows the fitted surface (fit function  Transient Effect of Fluid Flow on Dendrite Growth Direction in Binary Fe-C Alloys Using Phase Field in OpenFOAM being Eq. 15) for the entire data set with an R 2 value of 0.95.

CONCLUSION
Fluid flow-coupled numerical simulations based on the open source model developed have been carried out at different levels of undercooling for an Fe-C binary alloy. The phase-field method-based solidification model seems to be in quantitative agreement with the literature. The following conclusions can be obtained: At a particular level on undercooling, the tip growth speed and flow velocity at the tip changed with time and hence the bending angle. In this way, an array of deflection angle data along with flow velocity at the tip and tip growth speed was obtained for different levels of undercooling. An attempt to separate out the fluid flow velocity and growth speed effects has been made that can be useful in predicting the interface growth direction under industrial casting conditions by studying at what point of casting and up to what extent both the diffusion and bulk fluid flow interact with the solute layer ahead of the solidification interface. Anisotropy in interface energy-based separate fit functions dependent on flow velocity at the tip and growth speed have been postulated, which are an extension of the empirical correlations proposed by Takahashi. 16 In future, they can also be extended to full 3D, which at present is still a formidable challenge.