Simulation of Solidification Parameters during Zr Based Bulk Metallic Glass Matrix Composite’s (BMGMCs) Additive Manufacturing

After a silence of three decades, bulk metallic glasses and their composites have re-emerged as a competent engineering material owing to their excellent mechanical properties not observed in any other engineering material known till date. However, they exhibit poor ductility and little or no toughness which make them brittle and they fail catastrophically under tensile loading. Exact explanation of this behaviour is difficult, and a lot of expensive experimentation is needed before conclusive results could be drawn. In present study, a theoretical approach has been presented aimed at solving this problem. A detailed mathematical model has been developed to describe solidification phenomena in zirconium based bulk metallic glass matrix composites during additive manufacturing. It precisely models and predicts solidification parameters related to microscale solute diffusion (mass transfer) and capillary action in these rapidly solidifying sluggish slurries. Programming and simulation of model is performed in MATLAB. Results show that the use of temperature dependent thermophysical properties yields a synergic effect for multitude improvement and refinement simulation results. Simulated values proved out to be in good agreement with prior simulated and experimental results.


Introduction
Bulk Metallic Glasses [1] and their Composites [2] have emerged [3] as competitive structural engineering material [4] during last two decades and have attracted the attention of several major research clusters [5]- [19] around the globe to further investigate into scientific reasons behind their formation [20], microstructural evolution [21], property development and structure-property relationship [22] [23]. Main areas of research activity have been focused around probing into mechanical properties and their improvement as these materials have high hardness, strength and elastic strain limit but suffer from lack of ductility which make them brittle and they fail abruptly [24] [25] under the application of tensile and impact loading. The mechanism behind this is the formation and rapid movement of shear bands [26]- [33] in the volume of material by virtue of which material does not exhibit any yielding. In fact, they exhibit strain softening under tensile loading [34] [35] [36]. This renders them useless in practical structural engineering applications [37]. In addition to this, it is very difficult to scale up their size which could satisfactorily retain 100% monolithic glassy structure in whole volume [38]. It is also difficult to process chemical compositions which have high glass forming ability. Some of them are toxic and others are extremely expensive. Multicomponent nature of alloy itself makes it difficult to cast it into useful shapes as it has very high viscosity even at high temperature and is sluggish during processing. Various methodologies have been proposed to overcome this poor exhibition of mechanical objects which may include increasing quantity and complexity of shear bands by allowing their multiplicity due to a) self-interaction or b) at sites of foreign particles purposefully introduced to create more sites of shear bands multiplication [39]. Other mechanisms proposed are to study liquid to solid phase transformation in them by investigation under synchrotron light [40] [41] using container less levitation solidification techniques. The same effect is studied more elaborately in micro and zero gravity conditions on board International Space Station (ISS) by NASA and ESA (CETSOL) [42]. Various theories such as Liquid-Liquid Transition (LLT) [43], and phase separation prior liquid-to-solid transformation [44] have been proposed to explain their microstructural evolution but still, our knowledge about their exact mechanisms of formation and evolution at microstructural level [45] is scarce and limited. With experimental methodologies, efforts have also been focused to utilize well established solidification theories [46] [47] [48] to investigate and predict microstructural evolution during solidification in additive manufacturing using advanced multi-scale, multi-physics simulation and parallel programming strategies. Numerous algorithms [49] [50] [51] [52] have been proposed which describe and detail microstructure evolution in multicomponent alloys [53]. In this study, an effort has been made to further extend this approach and use one of these methods namely deterministic model to calculate numerical parameters required to understand microstructural evolution of ductile phase in bulk metallic glass matrix composites during solidification of liquid melt pool in additive manufacturing. An indigenous, in house model have been developed using famous KGT [47] and Rappaz [46] models to predict microstructural evolution during solidification of Zr based bulk metallic glass matrix composites. This is first effort of its kind which utilizes strengths of both approaches to formulate a comprehensive model for the prediction of microstruc-

Mathematical Model
In order to develop a model consider a dendrite tip evolving out of liquid as it starts cooling in melt pool of additive manufacturing. Its shape could be considered to be resembling a parabola. Using KGT model [47] following assumptions could be made about this type of microstructural evolution.
1) The solute field around the dendrite tip is given by Ivantsov solution.
2) The dendrite tip grows at marginal stability limit.
3) The diffusion coefficient d, is (tip) temperature dependent.
4) The segregation/partition coefficient, k, takes into account solute trapping; i-e, k is (dendrite tip) velocity dependent.

5) Initial partition coefficient (k o ) is temperature dependent and binary alloy
(Zr-Cu) is assumed to behave as multicomponent alloy.
6) The undercooling of tip (ΔT) is the sum of solute undercooling and the curvature undercooling.
7) The effect of convection is ignored.
In present study, however, a further practical approach is adopted which takes into account the calculation of supersaturations of individual constituents/components in an alloy which eliminates the fact that their diffusion fields superimpose and supports the actual conditions in which binary alloy system does not behave as ternary or multicomponent system (BMGMC). This is a con-

Nucleation
This is based on Oldfield theory of heterogeneous nucleation which describes a relationship between undercooling (ΔT) and grain density at each segment of interest (bulk liquid, mold wall and potent nuclei) in terms of Gaussian distribution to explain solidification. Two most important parameters sought after to be determined are, maximum nucleation density (n max ) and grain density (n(ΔT)).
Maximum nucleation density may be obtained by integral of nucleation distribution from zero undercooling to infinite undercooling.   Similarly, grain density is given by following equation where ΔT n and ΔT σ are mean undercooling and standard deviation of grain density distribution respectively.
With this, probability of happening of one event (nucleation) is given by nucleation probability (p v ) as described by Prof. Rappaz in his famous article [46]. cation of new grains is assumed to be governed by random process. This specific orientation is described to be controlled by three Euler angles θ, φ and ψ irrespective of grain nucleated at the surface of mold, potent nuclei or bulk liquid. This is described in Figure 4. In this figure, a reference frame, Oxyz, is attached to laboratory (e.g. the mold) in which direction Oz is perpendicular to mold wall. is given by [52] ( ) where A = constant which takes into account the fourfold symmetry of the dendrite along its trunk axis and the possible permutations of the <100> directions.
In general, dendrite growth direction of grain nucleated at the mold surface Putting in [5] 2 where, * i c = concentration of constituent i in liquid at dendrite tip (to be found), Comparing (10) and (11) ( ) where, a o = length scale related to interatomic distance and is estimated to be between 0.5 -5 nm and where, o D = Proportionality constant, Q = Activation energy, R g = Gas constant, T * = Tip temperature calculated by iterative method (described below) [51].
In a linearized phase diagram where, m i = slope of liquidus surface with respect to constituent i  (17) T L = Liquidus temperature for initial alloy composition, Γ = Gibbs-Thomson Another term can be generated from linearized phase diagram known as fictitious melting point of pure constituent [51].
This model is iterative model which is based on assigning final values to original value thus generating a loop whose explanation will be given in next section. c . Thus, a loop is generated which accounts for "to and fro" motion of information and iterative handling of data. This is the essence of generation of refined outputs and results. A schematic flow chart describing the working of model and interdependence of parameters is presented in Figure 5.
Accordingly, one has for a ternary system ( ) where, G c,i = solute gradient of constituent i in the liquid near tip which can be written as where, Γ = Gibbs-Thomson coefficient, m i = the slope of liquidus, c ξ = a function of the Peclet number and segregation coefficient, G = Thermal gradient.

Simulation Using Object Oriented Programming (OOP)
A computer program was written in MATLAB ® . Instead of fixing the Peclet number as was done in previous approaches [51], this program account for changing Peclet number with the change of each state value of system [63]. Furthermore, the program is based on transient transport processes and values are incorporated into original values using "for" loop and are assigned to their initial values using iterative approach. This helps in improving the efficiency of program and generation of fine mesh less results. Program asks for initial values and upon assigning initial values to dynamic variables, it generates first set of data which is repeated and assigned back to original variable to generate a loop. This process is repeated 100 times based on the number of iterations assigned in loop.
It also takes into account temperature dependent diffusion coefficient and velocity dependent partition/segregation coefficient in accordance with KGT model [47]. Exact thermophysical data of BMGMCs derived out of very recent studies [64] [65] was normalized and used in simulations. A correlation for the use of thermophysical data for major alloying elements of BMGMC system was developed. This was done based on their presence in same reactive group (transition metals) in periodic table. Nearest possible commonly studied element (Cr) [51] was chosen for generation of first set of data. Based on this, the parameters used in calculations are listed in Table 1

Results and Discussion
Model works by explaining dendritic growth in cast alloys during solidification dissipated. This is also shown in previous works by Rappaz et al. [51] but as equiaxed to columnar transition (ECT) stage reaches, a stability region evolves.
It happens due to very nice delicate balance between heat loss at dendrite tip and growth velocity of propagating solidification front. After that, again drop in T tip lidifying alloys [67]. These results are consistent with earlier observations [47] [61] and very recently observed phenomena of nanoscale solute partitioning in bulk metallic glass matrix composites [68]. This effect is not as sharp as observed in previous studies [47] because bulk metallic glass and their composites are sluggish alloys inherently and it is difficult to observe sharp interfaces and growth front in these materials specially at low thermal gradients (G).
Effect of ξ (a function of P ei ) on P ei Below graph (Figure 9) is representation of ξ with respect to Peclet number. It briefly shows that ξ (which is a function of P e ) shows almost decreasing linear relationship with the increase of rate of heat transfer from system to surrounding for all major alloying elements of hypoeutectic system at a fixed thermal gradient. It shows the effectivity of heat transfer process for BMGMC and describes that dissipation was proper and homogenous. The curves are generated by plot-  [72]. Simulation results are in nice agreement with previously reported behavior of elements in alloy systems [62].
It clearly shows that ξ value finally turns to zero for all elements generating a nice synergy between simulation and experimental observations in high speed additive manufacturing processes.

Evolution of segregation/partition coefficient with temperature
This is very interesting graph ( Figure 10) which shows the relationship between temperature dependent partition coefficients as a function of increasing temperature itself. It shows that partition coefficient is not uniform in its behavior when studied over a temperature range. It evolves with the change/evolution of temperature. Although assumed to be, and observed to be almost linear, its evolution is highly dependent on the gap chosen to calculate the values. The smaller the temperature gap, better will be the representation of actual behavior or evolution of partition coefficient over that period. For simplicity reasons, this effect is not studied in detail and general assumption that it shows linearity of evolution over temperature and time in the range of interest is made and adopted.
Again, thermal gradient (G) is kept constant at 100 K/mm.

Effect of dendrite tip growth velocity on supersaturation
This is final and second most important graph (Figure 11(a) and Figure   11(b)) of this study.   Apart from these, its salient features are; 1) Supersaturation of individual elements was measured to account for overall behavior of multicomponent system-an approach missing in previous studies 2) Due to scarcity, dispersion and unavailability of data, a correlation from nearest possible element in same group in periodic table was used.

Salient Features of Model
3) An effort was made to remove/reduce error by use of iteration based approach to refine model. 6) A unique approach based on segregation coefficient (k) as a function of temperature was adopted (Previously [47], it was only velocity dependent). 7) Slope of liquids (m) is taken to be concentration (C * ) dependent.

Conclusions
Following conclusions may be drawn out of this study. In essence, model comprises of extension of KGT theory for multicomponent systems beyond BLL model [62] employing real time temperature dependent conditions experienced in Additive Manufacturing. All these considerations must be taken into account while designing an alloy system (present case BMGMC) for a practical application processed by additive manufacturing. Any fault or carelessness will result in erroneous reading and in worst case scenario, catastrophic failure of components, parts and assemblies which must be avoided.