The Effect of Adding Boron in Solidification Microstructure of Dilute Iron-Carbon Alloy as Assessed by Phase-Field Modeling

Alloying element like boron, even in small addition, is well known to improve hardenability of steels. Its application can improve mechanical properties of steels and reduce alloying costs. Despite these benefits is not easy to cast boron steels, mainly in dynamical solidification process like continuous casting, due to their crack susceptibility. The strategy of using Phase-Field simulation of the solidification process is based on its proved capacity of predicting realistic microstructure that emerge during solidification under conditions even far from equilibrium. Base on this, some comparative simulations were performed using a three component dilute alloy in a two dimensional domain under unconstrained (isothermal) and constrained (directional) solidification. Simulation results suggested two fragile mechanisms: one related to a deep dendritic primary arms space and other due to the remelting of this region at low temperature. Both resulted mainly from the high boron segregation in interdendritic regions.


Introduction
Phase-Field model is an approach based on the concept introduced by Gibbs and Van der Waals, which states that solid-liquid interface has a physical dimension of nanometer magnitude 3,6,7 .It is widely accepted for prediction of solidification microstructure of pure substance as well as alloys of different solute contents [3][4][5] .For two phases system it uses an order parameter (φ) as a state variable that identify the phases, i.e., it is constant inside the bulk phase (ex: 0 for liquid and 1 for solid) and change continuously between these values along the interface 3 .It does not have the disadvantage of tracking the interface all the time.In fact the model solution determines its position.
Developments of Phase-Field models can be addressed to Van der Waals who perceived that the addition of gradient term in free energy functional has an effect of creating interface with certain thickness 6 .Cahn and Hilliard 8 , after rediscovering this idea, associated it to the interfacial tension and solid-liquid interface thickness, allowing Phase-Field model parameters to be related to them.
Important contribution came from Langer 6 who proposed the first Phase-Field model for solidification of pure substances.No less important was the pioneer work from Kobayachi 9 , whose first numerical simulation introduced the formalism to incorporate anisotropy that is used until present date.
In general, Phase-Field models of solidification of three component alloys without convection in liquid or solid state tension effects are composed by four diffusion equations: one for the order parameter; one for temperature and two for the solutes.These can be phenomenological, based on free energy functional or developed from entropy functional 3,7 .Despite both approaches reproduce similar results 3 , the most used one is the phenomenological, mainly due to its faster convergence and its facility to obtain model parameters based on material properties 3 .Therefore, according to Cahn and Hilliard 8 and Landau Ginsburg 7 the free energy functional (F) can be formulated as follows: , , where f (φ, N i , T) is a free energy density that is function of order parameter (φ), concentrations of solutes (N i ) and temperature (T); ξ is a energy density contribution from the gradient of order parameter and V is system's volume.
The key point of the phenomenological Phase-Field model is to consider that free energy decreases monotonically as solidification advances, i.e: where t is time and M is a Phase-Field mobility term.From the derivative of this functional it follows the general expression for the diffusion of order parameter: In order to describe f (φ, N i , T) it has been used a weighing rule 3 , defined as followed: where h p (φ) is an interpolation function which varies along solidliquid interface from a fixed value in solid (ex.: h p (1) = 1) to a fixed value in liquid (ex.: h p (0) = 0), W is the interface energy density and Φ(φ) is a double well function.In the above expression (Equation 4) the first term of right side is related to the interface energy and the others are the bulk free energy of solid (f S ) and liquid (f L ) phases, respectively.
and isotropic material without noise in an one dimension domain, it is possible to relate these parameters to the interfacial tension (σ) and interface thickness (2λ) as follows 11 : The Phase-Field mobility term (M 0 ) formulation was developed for alloys by Karma 12 using asymptotic analysis and independently by Kim et al. 13 using solidification kinetic consideration.Echebarria et al. 14 , based on work of Karma 12 , developed a model for 2 components alloy without diffusion in solid (one side model).Ode et al. 15 proposed another approach for a three component alloy, however, considering the same simplification and vanishing inverse of kinetic coefficient (β).On the other side, Ramires et al. 16 presented an alternative formulation for 2 components alloy, considering a nonisothermal domain, no diffusion in solid and vanishing inverse of kinetic coefficient.More recently, Kim 10 reviewed this development and proposed a new formulation for multicomponent alloy, however, considering the same simplifications of Ode et al. 15 .In this context, a recent work by Ohno et al. 17 for a 2 component alloys was developed to deal with diffusion in solid.They 17 considered an isothermal case and vanishing inverse of kinetic coefficient, obtaining results 17 similar to the one side model.
Based on Kim 10 , Ramires at al 16 and Ohno et al. 17 the present work employs a custom developed approach to three component alloy, for non-isothermal conditions, non-vanishing kinetic coefficient and diffusion in solid described by the following expression: where D L i (i = 1, 2) are liquid solute diffusivities and N i Le are equilibrium solute concentrations in liquid phase, which can be obtained from phase diagrams, K T is the thermal conductivity and A is a constant.This formulation reduces to the one from Ramires et al. 16 for two components alloy as well as to the one from Karma Rappel 18 for pure materials.
The solute diffusion equations are based on conservative Onsager extrapolation without cross effect on diffusivities terms (D i,l ), i.e.: . 1 Using Kim 10 approach and the 'antitrapping' correction term from Karma 12 ( ) i j  the expression for non-dimensional solute concentrations (U i ) can be written as follows: . . 1 and Hence, it is possible to obtain a general formulation for the order parameter diffusion equation, i.e.: ( ) where Φ´ is the derivative of Φ in relation to φ and source depends on the bulk free energy of solid and liquid phases.Following the multicomponent approach from Kim 10 and considering constants solutes activity coefficients of infinite dilution, it is possible to obtain the following expression (see detailed development in Appendix A): where where T m is the melting temperature of pure solvent, R is the universal gas constant, DH f is the enthalpy of fusion of pure solvent, V m is the alloy molar volume, k i are the solutes equilibrium solute partition coefficients, N i∞ are the initial alloy solute concentrations and h r (φ) is another interpolation function with the same properties of h P (φ).
In order to simulate complex microstructure, anisotropy was included in the model as proposed by Kobayachi 9 , where the parameters ξ and M were made to vary according the angle (θ) between solid-liquid interface normal vector and horizontal axis (x) in the present domain of two dimensions (2D).The final diffusion Equation for the order parameter is as follows: .
In this case the following expressions were used 7 : where δ ξ and δ M are the strength of anisotropy, j is the number of anisotropy, θ 0 is the main growth angle, A is the amplitude of noise and ν is random number obtained from a uniform distribution.Predictions for ξ 0 and W follows from the work of Cahn and Hilliard 8 on surface tension.In that case and considering the analytical solution of Equation 8, simplified to describe the equilibrium of pure avoids joint thermodynamic and Phase-Field calculations inside each time step, reducing the numerical processing time.
Isothermal solidification constitutes the maximum rate of solidification of unconstrained dendrite growth.The alloys were submitted to different undercoolings and their solid-liquid interface evolutions were determined until reaching steady state condition.The interface thickness (2λ) was set up to 8 × 10 -8 m and the inverse of kinetic coefficient was estimated as an arithmetic average between pure iron, determined by collision limited teory 22 , and the value reported by Ogushi e Suzuki 23 for a higher carbon Fe-C alloy, i.e.: 0.54 s.m/K.The properties used in simulation are summarized in Table 1.The results presented in Figure 3 show that boron addition has the effect of narrowing the dendrite primary arm thickness.That phenomenon can be attributed to the tip radius selection, associated to the boron solute gradient around the solid-liquid interface.
In order to determine the dendrite tip radius, a parabolic function (fy) was fitted to data close to the main dendrite tip arm.This is a plausible approach which is based on experimental evidences 24 .
Thus tip radii (r) formed in different solidification velocity were determined using the following formulas: where ( ) Here fy´ and fy´´ are the first and the second order spatial derivative of fy in relation to x direction (horizontal direction) and K is the curvature of the dendrite main tip arm.The results of this analysis are summarized in Figure 4.It shows that the dendrite tip radius decrease as solidification velocity increases until a limit where it starts to increase.This point is related to the high velocity dendritecellular morphological transition 24 .It is a common behavior of all three alloys.In these cases, the increase of boron addition reduced the dendrite tip radius and also it seems to reduce the velocity where this transition should happen.
According to the marginal stability theory, the tip selection is a result of stabilizing forces (e.g.: solute gradients) and nonstabilizing forces (interfacial tension), modulated by the selection parameter 24 , which is commonly used in one dimensional analytical tip model as a easier tool to map the morphological transition of materials under solidification.However is not easy to determine this parameter because it can depend on alloy composition and anisotropy.Nonetheless, Phase-Field simulation can be used to estimate it, and for the three alloys under consideration the results are summarized in Table 3.It seems that the higher the boron content the higher the selection parameter.
where D S i are solid solute diffusivities and h d (φ) is another interpolation function with the same properties of h P (φ) and h r (φ).
In order to describe the temperature field in directional solidification, a quasi-steady state approximation was used 14 .It takes in consideration a sample traveling in one direction (ex.z) with a constant pulling velocity (ex.V P ) under a fixed thermal gradient (G).This approach can be implemented as follows: where T 0 is the initial temperature.Equations 8 and 18 were solved using finite volume method 19 with uniform grid.For the time derivative, the explicit scheme was used.Despite some limitation regarding the definition of the time step, it is easy to implement and has the advantage that the Phase-Field Equation does not need to be solved outside the interface region.

Results and Discussion
This model was already validated against analytical solution and experimental data for pure material and 2 components alloys 20,21 .However, some new calculations were first performed in order to compare their results with equilibrium thermodynamic data for dilute ternary alloys.As far as the rate of solidification is small, solute concentration should approach the values given by phase diagrams.That constitutes an additional validation procedure.
The first analysis was performed for Fe (BCC)-C-P dilute alloy.Input data are summarized in Table 1.Phase-Field model was set up to use a flat interface thickness (2λ) of 8,0 × 10 -8 m and vanishing inverse kinetic coefficient.The domain was defined with 3 nodes in horizontal direction.Anisotropy effect and noise were set up to zero.A flat solid-liquid interface was allowed to advance in vertical direction until reaching steady state.The resulting solute concentration profiles of carbon and phosphorus are presented in Figure 1 together with three isothermal slices of dilute Fe(BCC)-C-P phase diagram reported by Ode et al. 15 .The agreement is remarkable.
A second Phase-Field calculation was focused on Fe(BCC)-C-B dilute alloy.The procedure was similar to that applied to the Fe(BCC)-C-P dilute alloy.Again, input data are from Table 1.As it can be seen from Figure 2, a good agreement was obtained between the present Phase-Field calculations and phase diagram data obtained from Blasek et al. 1 .
In order to analyze the effect of boron on the solidification microstructure, the present work simulated three alloys with different boron contents (Table 2).Their equilibrium solute concentrations necessary to determine the mobility parameter M (Equation 16) were determined previously as function of temperature.This procedure   The addition of boron in the dilute Fe(BCC)-C alloy should increase the non-stabilizing forces during solidification because it creates an extra solute diffusion field.Therefore one should expect an increase of solid-liquid interface disturbances at a certain solidification velocity.In order to verify this effect, the domain and the simulation time were expanded.The way to perform this is to increase the solid-liquid interface thickness.It has been shown that is possible to obtain convergence of Phase-Field simulation to the analytical solution of Mullins Sekerka instability theory, and therefore to obtain consistent results, using interface thickness smaller than dendrite tip radius 14 .In order to make this simulation feasible with the present model, the interface thickness (2λ) was increased to 2 × 10 -7 m.The grid size and the time step were set up to 5.0 × 10 -8 m and 5.0 × 10 -9 s, respectively.Morphological results on the three alloys under analysis can be compared in Figure 5. Actually, it seems that, for similar solidification time, the increase of boron addition lead to an intensification of secondary arms growth.Also, it can be observed some indication of secondary arms dynamic coalescence phenomena, where larger arms grow faster and smaller ones are melted back.This suggests that boron addition can lead to an increase in the secondary arms spaces due to increasing dynamic coalescence, which is in agreement with the theory of secondary arms growth 22 .
Solute segregation is an important phenomenon that can strongly influence macroscopic properties of materials.Regarding this, Figure 6 and 7 show the carbon and boron segregation resulting from morphologies presented before (Figure 5).It is clear the increase of carbon and boron concentrations in the interdendritic regions.These higher concentrations could lead, according to Blasek et al. 1 , to a low temperature remelting (around 1200 °C) of already solidified interdendritic region.This phenomenon has been associated to the ductility loss and crack formation during continuous casting of boron steel 1 .
In directional solidification (constrained growth) the solidification front is inside the range defined by solidus and liquidus temperature.At steady state, pulling velocity is equal to the solidification velocity and morphology is mainly columnar, which is the most common microstructure in continuous casting process.In order to simulate this process, Equation 21was included in the present Phase-Field model.However, typical continuous casting cooling rate is not feasible to be simulated with the present implementation due to the large domain and time required to resolve the microstructure, which would result in a huge processing time (ex.: a system with 10 8 grid points and 10 7 time steps).Therefore a faster cooling rate was set up based on work of Miettinen 25 .In this case, fixed values of temperature gradient (G), pulling velocity (Vp) and initial temperature (T 0 ) were set up.Two seeds were implemented in both upper and lower left corner of the domain.The distance between them was defined by primary arms space as determined by Miettinen 25 for these conditions, i.e.: 1.5 × 10 -5 m.The dimension of the domain in the other direction (now z direction) was set up in order to have a well developed microstructure and, at same time, ensure a certain superheat at the high side of the domain, i.e.: 7.5 × 10 -5 m.Interface thickness, grid size and time step were defined as 3 × 10 -7 m, 7.5 × 10 -8 m, 2 × 10 -8 s, respectively.Field simulations.Therefore, in order to reduce the susceptibility for crack formation in this steel one should focus on minimizing boron segregation.As product hardenability improvement is due to the presence of solid soluble boron 1 and, as reported before, its solubility in this phase is very low (figures 8 and 9), the first approach should be the reduction of its addition to a minimum defined by its solubility limit in solid phase.Another possible strategy would be the use of another element capable of forming a high temperature third phase with boron.
It is important to point out that present simulation were performed using properties independent of temperature and concentrations.Nonetheless, the excellent agreement of Phase Field results with equilibrium phase diagrams indicates that this simplification can work properly concerning dilute solution.
Another point is on some larger values used in present work to define interface thickness necessary to extend the simulation domains.In this regard, the present Phase-Field model were implemented using the 'antitraping' flux (Equation 20), therefore one can not expect any artificial solute trapping, which could reduce solute segregation close do solid-liquid interface.Besides, these values were always kept lower enough to resolve dendritic arms.
In actual continuous casting process the solidification front is expected to be sensitive to the liquid convection that can affect solute segregation.This phenomenon can be driven by fluid flow of liquid steel entering in the mould in the first stage of solidification and also by shrinkage of solid.Also, change in liquid density and surface tension gradient (Maragoni's convection) can provoke convection, both due to concentration gradient formed at solidification front.As these phenomena were not include in the present Phase-Field model one should expect some differences between present results and actual process.Nonetheless, based on sample presented in Figure 10, The solid-liquid interface growth in directional solidification is driven by constitutional undercooling, which in turn is defined by solute segregation at the interface.Simulation results on solute segregation using the present Phase-Field model can be seen in Figure 8, where the darkest areas correspond to the solid phase and the solute gradients between the two arms comes from segregation phenomena, where solid phase rejects solute to the liquid phase.From these results it is possible to determine the constitutional undercoolings, which are presented in Figure 9 (negative values).As it can be seen, the interdendritic region is not totally undercooled in the present simulation.Additionally it is noticeable that boron addition increases the interdendritic region making the main dendrite arms thinner.On the other side, close to dendrite tip, the constitutional undercoolings become higher as boron addition increases.
As a consequence of the above results, liquid fraction between the thinner arms increases, resulting in more fragile microstructure, which would be under transverse tensile stress common in continuous casting.Indeed, it has been reported severe crack formation parallel to the dendrite main axis during casting of high boron content steel, as shown in Figure 10, which is consistent with present Phase-Field simulations.
Based on simulation of present work, it can be concluded that two mechanisms making high boron alloyed steel to be more sensitive to crack formation during continuous casting could be present: one occurring at solidification front and other at lower temperature (around 1200 °C).The fist one could be related to a deep dendritic primary arms space that would concentrate the tensile stress common to the continuous casting process.The second one could be due to the remelting of this region at low temperature resulting in a large localized ductility loss.Both mechanisms resulted mainly from the high boron segregation in interdendritic regions, predicted by Phase-    one can speculate that convection could be acting to accelerate the primary arm growth in its main axis direction, where convection is more effective, due to an increase of constitutional undercooling (increase in latent heat removal and concentration gradient).On the other side, in interdendritic region, where space limits convection one could expect results similar to the present work, i.e., without convection.Therefore the final results could be a thinner dendritic primary arm than that predicted by the present work and, therefore, a deepest interdendritic region.i.e., the presence of convection could induce an increase in the fragility of solid-liquid front.However, in order to confirm the above comments it is important to include in the Phase-Field model the convection phenomenon.

Conclusions
The present Phase-Field model reproduced remarkably well the equilibrium solute concentration of rich iron size of Fe(BCC)-C-P and Fe(BCC)-C-B ternary phase diagrams.
Isothermal simulations reported a decrease of arm thickness, and dendrite tip radius with incremental addition of boron.
The relation between dendrite tip radius and solidification velocity were reproduced for the three alloys simulated.In this regard incremental addition of boron seems to reduce the velocity where dendrite-cellular transition can occur.
According to the isothermal simulation, boron addition increases solid-liquid interface instability, inducing faster secondary arms formation and its dynamic coalescence.That phenomenon could increase secondary arms spaces.
Interdendritic carbon and boron segregations and the equilibrium phase diagram from Brasek et al. 1 suggest that remelting of already solidified interdendrite region could occur in lower temperature range (around 1200 °C), decreasing hot ductility of the material.
Directional solidification process simulation showed the effect of incremental boron addition in narrowing the dendrite main arms and the increasing of interdendritic liquid fraction.These results can be attributed to the carbon and boron segregation that influence constitutional undercoolings.
Present Phase-Field model predicted two mechanisms that make high boron alloyed steel to be more sensitive to crack formation during continuous casting: one related to a deep dendritic primary arms space that would concentrate the tensile stress common to the continuous casting process and other due to the remelting of this region at low temperature resulting in a large localized ductility loss.Both mechanisms resulted mainly from the high boron segregation in interdendritic regions.
Even though the present model did not consider convection effect present in actual process, the sample of high boron alloyed steel took from it showed cracks consistent with the fragility mechanism present above.
In order to reduce the sensitivity for crack formation during continuous casting of boron steel, segregation of this element should be minimized.The ideal value would be the maximum solubility of boron in solid phase.Another approach could be the use of another element that reacts with boron to yield a third phase.

Figure 1 .
Figure 1.Solid-liquid equilibrium solute concentration in dilute Fe(BCC)-C-P alloy: points determined by the present Phase-Field model; lines are from equilibrium phase-diagram.

Figure 2 .
Figure 2. Solid-liquid equilibrium solute concentration in dilute Fe(BCC)-C-B alloy: points determined by the present Phase-Field model; lines are from equilibrium phase-diagram.

Figure 4 .
Figure 4. Effect of boron addition on tip radius-solidification velocity relation.

Figure 3 .
Figure 3. Phase-Field calculation of the effect of boron addition on dilute Fe-C alloy: a) FeC alloy; b) FeCB1 alloy; and c) FeCB2 alloy.

Figure 5 .
Figure 5.The effect of boron addition on dendrite secondary arms growth: a) alloy FeC; b) alloy FeCB1; and c) alloy FeCB2.

Figure 10 .
Figure 10.Sulfur print of a transversal slice of continuously casted slab of high boron steel 2 .

Table 1 .
Data used in the present Phase-Field calculation.

Table 2 .
Iron-carbon alloys simulated in present work.

Table 3 .
Selection parameter for the 3 alloys simulated by the present Phase-Field model.