Non-linear analytical modeling and simulation of microstructural evolution in Zr based Metallic Glass Matrix Composites during solidification using MATLAB® – further results

Metallic glass matrix composites are unique new materials which have superior mechanical properties as compared to existing conventional materials. Owing to this, they are potential candidates for various future structural applications. However, they suffer from disadvantages of brittleness, poor ductility and little or no toughness and they manifest catastrophic failure on application of force. Their behavior is dubious and require extensive experimentation to draw conclusive results. In present study, which is continuation of previous study by author, an effort has been made to overcome this pitfall by simulation. A quantitative mathematical model based on KGT theory has been developed to describe nucleation and growth of crystalline phase dendrites in glassy matrix during solidification. It yields information about numerical parameters to understand the behaviour of each individual element in a phase in multicomponent sluggish slurry and their effect on final microstructural evolution. Model is programmed and simulated in MATLAB®. Its validation is done by comparison with identical curves reported in literature previously for similar alloys. Results indicate that effect of incorporating all heat transfer coefficients at macroscopic level and diffusion coefficients at microscopic level play a vital role in refining the model and bringing it closer to actual experimental observations. Two types of hypoeutectic and eutectic systems namely Zr65Cu15Al10Ni10 and Zr47.5Cu45.5Al5Co2 respectively were studied. Simulation results were found to be in good agreement with prior simulated and experimental values.


INTRODUCTION
Bulk Metallic Glasses [1] have emerged [2] as competitive structural engineering material [3] during last two decades and have attracted the attention of several major research groups [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] around the word to further probe into the science and engineering behind their formation [19], microstructural evolution [20], property development and structureproperty relationship [21,22]. Main areas of focus have been investigation of evolution of mechanical properties in these materials as despite their high hardness and very high elastic strain limit, they do not exhibit any tensile ductility and fail catastrophically [23,24] under tensile and impact loading. This happens due to rapid movement of shear bands [25][26][27][28][29][30][31][32] in the volume of materials by virtue of which material does not exhibit any yielding. Infact, they exhibit strain softening rather than strain hardening upon deformation in tensile loading [33][34][35]. These behaviors render them useless in practical structural engineering applications [36,37]. In addition to that, they are limited by the size in which 100% monolithic glassy structure can be achieved [38]. Chemical compositions in which it is possible are rather limited and difficult to process into useful shapes due to their multicomponent nature which make the material thick / slurry like and sluggish during processing. Various mechanisms have been proposed to overcome these drawbacks. Some of them include; increasing the number 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 sites of multiplication of shear bands [39]. On the other hand there have been efforts to study evolution of phases during liquid to solid transformation of these materials by their investigation in synchrotron light [40,41] using container less levitation solidification techniques. Efforts have also been made to study same phenomena in micro and zero gravity conditions on board International Space Station (ISS). Various theories such as Liquid -Liquid Transition (LLT) [42], Phase separation prior liquidtosolid transformation [43] have been proposed to explain their microstructural evolution but still, there is dearth of knowledge about their exact mechanisms of formation and microstructural evolution [44]. With experimental methodologies, efforts have also been focused to utilize well established solidification theories [45][46][47] to investigate and predict microstructural evolution during solidification using advanced multi scale and parallel modeling and simulation strategies. Various algorithms [48][49][50][51] have been proposed which are associated with development of microstructure in multicomponent alloys [52]. In this study, an effort has been made to further extend this approach and use deterministic methods to calculate numerical parameters necessary to understand microstructural evolution of ductile phase in BMGMCs during solidification. An indigenous, in house model have been developed using famous KGT [46] and Rappaz [45] model to predict microstructural evolution during solidification of Zr based BMGMC. This is first effort of its kind which utilizes strengths of both approaches to formulate a comprehensive model for the prediction of microstructure primarily explaining dendrite tip temperature and dendrite tip radius as a function of growth rate / dendrite tip velocity.

THEORY (MATHEMATICAL MODEL)
Durig this study only deterministic apporach is adopted to arrive at numeical quantitative parameters necessary to describe solidification and microstructre evolution during initial stage of alloy cooling. Further visual obervation of real time microstructure may be described by coupling [47,53,54] of present model with well established probablistic models [45,55,56] based on cellular automaton theory [57,58]. For present case, considering a free dendrite tip of parabolic shape, KGT model [46] is taken to be responsible for microstructural evolution. Its basic assumptions are following; 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 / partitition coefficient, k, takes into account solute trapping; i-e, k is (dendrite tip) velocity dependent. 5. Initial partition coefficient (ko) is temperature dependent and binary alloy (Zr -Cu) is assumed to behave as multicomponent alloy. 6. The undercooling of tip (AT) is the sum of solute undercooling and the curvature undercooling. 7. The effect of convection is ignored.
In present study, hoewever, a further practical approach is adopted which takes into account the calculation of supersaturations of individual consitituents / components in alloy which rules out that their diffusion fields superimpose and binary alloy system is assumed to behave as ternary or multicomponent system (BMGMC). This is a contradiction to assumption 5 above in basic KGT model-a unique approach adopted here to explain microstructural evolution in detail.
There are three main velocities of interest here (Figure 1 (ac)) when a beam of high energy (electron or laser) travels on surface of specimen in additive manufacturing setup (Figure 2

Nucleation
This is based on Oldfield theory of heterogeneous nucleation which describes a relationship between undercooling (dT) and grain density at each segment of interest (bulk liquid, mold wall and potent nuclei) in terms of Gaussian distribution. Two most important parameters namely, maximum nucleation density (new) and Grain density (n(AT)) are sought after to be determined. Maximum nucleation density may be obtained by integral of nucleation distribution from zero undercooling to infinite undercooling. dci dbT' Similarly, grain density is given by following equation n(A7') -exp [1] dbT' [2] where AT» and ATR 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 ) as described by Prof. Rappaz in his famous article [45]. p r [3] i-e if at any instant of time t, p exceeds r, nucleation will occur. pbrim.VcA where dn, = grain density increase and VcAone cell volume (measure by noting all dimensions of cell assuming it to have square shape) (Figure 3(a)). A change in state index of a cell represents its growth as depicted by cells in Figure 3 [47] Dendrite growth orientation Second part of the model deals with determining dendrite growth orientation i-e the direction in which some dendrites preferably grow faster and longer as compared to others due to balance between geometrical and kinetic variables. This also highlights and points towards grain competition and selection mechanisms. Two important parameters are considered for assigning and determining grain orientation. a) Growth of first grain as a result of heterogeneous nucleation at mold wall or potent nuclei and b) Location of further subsequent new grain(s) and their crystallographic orientation. For example, for cubic metals, the preferential growth directions of dendrites are given by direction of easy heat flow which is along <100> crystallographic direction / orientation. During early stage of solidification, a nucleus grows at the surface of mold or potent nuclei in the form of hemispherical surface. This surface becomes unstable and then dendritic after a certain incubation time and growth occurs with main trunk and arms coinciding with <100> crystallographic direction. The location of new grains is assumed to be governed by random process. This specific orientation is described to be governed by three Euler angles 8, ‹p and irrespective of grain nucleated at the surface of mold, potent nuclei or bulk liquid. The probability dp(8, ‹p, y) that a newly nucleated grain has its main trunk orientation in the range [d, d + d8] and [p, p + dp] and one of its set of secondary branches within the orientation [p, p + dp] is given by [51] dp(8, ‹p, y) = A . sin 8 . d8 . d‹p . dy [4] 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 determines the time during which grain can survive competitive growth of its neighbors.
Grain growth / Growth kinetics / Dendrite stability theory This section concerns kinetics associated with growth of already formed grains in bulk liquid, mold surface and potent nuclei. A unique feature adopted here constitutes the determination of supersaturation of individual elements in multicomponent alloy (BMGMC) systems. This approach originates from the fact that in contrast to conventional castings in which undercooling related with thermal diffusion, attachment kinetics and curvature is small, in multicomponent systems basic KGT model must be used with certain modifications which not only accounts for superimposition of solute fields around each dendrite tip but also incorporate determination of supersaturation for each individual component (Zr, Cu, Ni, Al and Co) of alloy system. This supersaturation D, is a function of Peclet number, Pe, Putting in [5] Pe; [10] where n,= c;*(1 -k;) [11] c;'concentration of constituent i in liquid at dendrite tip (to be found), c ;initial concentration of constituent i, £t= partition coefficient for this constituent i (velocity dependent) Comparing [10] and [11] e' D( en -u du [12] However, k; [13] [14] where n = length scale related to interatomic distance and is estimated to be between 0.5 -5 rim and D, --D e(" ' R • ) [15] where D -Proportionality constant, Q = Activation energy, Rg' Gas constant, T* = Tip temperature calculated by iterative method (described below) [50] In a linearized phase diagram Using eq. (13) and (18), eq. (16) becomes where R = Dendrite tip radius, V = Dendrite tip velocity This model is iterative model which is based on assigning final values to original value generating loop whose explanation will be given in next section. 100 iterative cycles are used to generate homogeneous and normalized data. In general while writing the program, reading it and executing it, D depends on /p{Pe;) and c;', c;*depends on k, and /p(Pe;), k, depends on D" D, depends on Tip temperature and finally, Tip temperature depends on 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.
The criteria used to specify radius of dendrite tip (R) is assumed to be given by marginal stability wavelength of planar wave front (as given in Mullins G " = solute gradient of constituent i in the liquid near tip which can be written as

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 [50], this program account for changing Peclet number with the change of each state value of system. 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 [46]. A correlation for the use of thermophysical data for major alloying elements of BMGMC systems was developed based on their presence in same reactive group (transition metals) in periodic table and nearest possible commonly studied element (Cr) [50] was chosen for generation of first set of data. Based on this, the parameters used in calculations are listed in Table 1. The thermal gradient (G) which is a function of additive manufacturing (AM) condition is taken as free number and only one value (100 K / mm) is used for calculations. A unique feature of program is it can work for; and be tailored according to the need of any material (alloy) and manufacturing process (additive manufacturing). Data on thermophysical properties of metals and alloys (especially exotic materials) as they change with the change of time and temperature is scarce and more efforts are needed in this front. BMGMC samples were casted in two ways. Firstly, they are made in form of wedge using vacuum suction casting system in lab scale Vacuum Arc Melting (VAM) button furnace at CSIRO -Manufacturing. The process consists of carefully calculating raw material based on weight percentage of each element in the alloy system. These powders / granules / chucks are subsequently mixed using hand spatulas to a homogeneity observable by naked eye. For their positioning, handling and control inside enclosed chamber of VAM, they are wrapped in an aluminum foil which not only protects the powders but also serve to fulfill the requirement of alloying element in sufficient quantity in original mix. This Aluminum foil wrapped toffee is placed in horizontal slot in the Copper hearth of water cooled furnace at appropriate time after which, it is melted to get solid chuck / button for subsequent research. During second approach, casted wedge samples were subjected to laser solid forming (LSF) [62] in Additive Manufacturing Setup. Model theory was developed and tailored to describe microstructure evolution during both processes.

RESULTS AND DISCUSSION
Model works by explaining dendritic growth in cast alloys during solidification by manipulating physical process parameters with the change of heat and mass transfer coefficients. Its unique feature is it explains the behavior of multicomponent alloys in terms of transient state variables. An effort is made to keep constant values to a minimum to get real picture of actual physical processes. Boundary conditions of solidification phenomena are kept open which makes model more rigorous and robust and it is possible to apply this for various alloys systems under various conditions. Following results and graphs have been generated after writing script of solidification code and running it in MATLABO R . Figure 5 is the plot of evolution of solutal Peclet number of different elements of BMGMC system against their dendrite tip growth velocities. It shows the transport behavior of individual elements with the change and evolution of state variables. It is clearly evident that highest growth velocity has been attained by Zirconium with very little heat loss from system to surrounding. One of the reasons this happens is due to higher atomic weight and size of Zirconium atoms by virtue of which they can travel higher and longer in slurry like thick multicomponent system with causing very little or no interaction with surrounding material / elements. That is they; owing to their large size attain a more favorable position and orientation while navigating through thick fluid which consists of elements of various size ranges and undergoing rapid heat loss and due to solidification. Thermal conductivity of Zirconium is also not very high which suggest that very little heat loss will occur as a result of Zirconium, its presence and movement through fluid. Overall large atomic weight remains as the major factor for its increased dendrite tip velocity through thick slurry. A small effect of gravity also contributes towards this overall effect (not considered in detail here). Second element of interest in this multicomponent system is Ni, which is not been able to attain enough velocity in growing dendrite tip due to its lower atomic weight. Thus, it remains ineffective in creating shear zone around it which can allow it to penetrate and gain speed in thick slurry around it. It also does not have very high thermal conductivity thus not a lot of heat gets dissipated because of it to surrounding and again it reflects its inability to create low viscosity zone around it which may facilitate it to gain high speeds. A drastically high speed and very high Peclet number over a range is depicted by Copper. This was not astonishing, as it is the element of highest thermal conductivity with an intermediate position with respect to its atomic weight in periodic table of elements. These both features; specially high thermal conductivity gives it an edge to create a region of very high heat dissipation around it facilitating its motion at a high speed over a range of Peclet number in a thick slurry of multicomponent alloy systems. The atomic size and configuration of copper is also favorable for attaining this high speed. Thus, for this thermal gradient, it remains as the most important element contributing towards overall speed / dendrite tip growth velocity of system. The last element is Aluminum, which despite of its high thermal conductivity does not attain very speed. This happens because; the atomic size and weight of aluminum is not very high which does not help it to gain high velocity in the thick slurry of multicomponent systems at tip of dendrite. Also, as the heat keeps on getting dissipated from its surrounding and it has relatively low atomic weight as compared to others, it does not attain high speed and gets stuck in its own region. This reasoning is based on atomic size, weight, electronic configuration and thermal conductivity of individual elements. However, no data exists about actual dendrite growth velocity of this complex BMGMC system and further experimental research to measure this is needed.

Effect of tip temperature on dendrite tip velocity
Below graphs ( Figure 6 a and b) shows relationship between dendrite tip velocity as a function of dendrite tip temperature for all elements in Zr6sCui 5Ali0Niio system. Figure 6 b is an extrapolated graph as compared to Figure d a which is actual plot. It clearly shows three different regions which are distinctive of velocity evolution of alloying elements in the system as temperature change. It depicts behavior of alloying elements in sluggish slurry at different stages. There is a sharp increase in T r at slow V¿because dendrites are equiaxed in nature and due to their rapid mechanical interaction with each other a lot of heat is accumulated in small area. This is more evident in multicomponent alloys. In addition, for Zirconium only (orange curve), its large atomic size and weight also contribute towards increase. This is early / initial stage of solidification. Another reason for this is there is planar wave front during initial stages which does not allow the development of high surface area i-e surface area / volume ratio remain low and not a lot of heat gets dissipated. This is also shown in previous works by Rappaz et. al. [50] 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«r with increase in velocity is observed. This happens as dendrite of alloying elements (Zr and Cu) gain speed; they start creating regions of very high surface / volume ratio around them in the body of melt and at the tip of advancing dendrite. In other words, faster moving dendrites become source of high heat dissipation. Owing to this, temperature comes down. This effect is absent in Al and Ni. This effect is predominantly not observed in actual plots (Figure 6 a) where dendrite tip temperature is plotted against dendrites growth velocity in actual regime. It can be clearly observed that growth velocity reaches maximum values at a certain temperature and then saturates there. No further drop in temperature is observed beyond that value. This is actual phenomena observed in real casting of Zr6 CuisAli oNiio. As the very nature of alloy does not allow dendrites to become elongated or develop very high surface area for rapid heat dissipation after saturation of equiaxed to columnar transition (ECT). This is actual situation and must be taken into account while further development.

Effect of § (a function of Pei) on Pei
Below graph (Figure 7) is representation of § with respect to Peclet number. It briefly shows that § (which is a function of Pe) 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 gives measure that dissipation was proper and homogenous. The curves are generated by plotting solution of present model by the use of indigenous MATLAB code incorporating different transient thermo-physical data of each individual alloy system from literature [63][64][65][66][67][68]. Simulation results are in nice agreement with previously reported behavior of elements in alloy systems [61]. It clearly indicates 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 8) 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 which you will choose to calculate the values. The smaller the temperature gap, the 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 is kept constant at 100 K / mm.

Effect of dendrite tip growth velocity on supersaturation
This is final and second most important graph (Figure 9 a) of this study. It describes the relationship of supersaturation for Zirconium only with dendrite tip growth velocity. This describes how a dendrite tip evolves and what happens to solute field around it as time passes and temperature decreases. As was expected and seen in experimental studies, it can be observed in this simulation as well that solute field around dendrite tip decreases almost in linear uniform homogenous fashion over a temperature range. This means that, solvent (alloy) keeps on rejecting solute as it solidifies, and its dendritic structure evolves. This is also consistent with experimental observations made for a range of alloy systems previously [46,59]. The drop in supersaturation is very sharp in first region as equiaxed grains are forming during this stage which are large in number and form very quickly followed by achievement of stability (a very small flat region in curve) because of equiaxedcolumnar transition followed by almost linear drop of slope of curve due to formation of columnar dendrites which are large in shape, move fast thus become source of large surface / volume ratio and faster rate of heat and mass transfer. Similar observations and simulation results are expected for other constituents provided temperature gradient is kept constant at 100 K / mm. These, namely for Ni, Cu, and Al are described in Figure 9 (bd).