Vol.10 No.03(2018), Article ID:83155,24 pages

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

Muhammad Musaddique Ali Rafique

Eastern Engineering Solutions LLC, Detroit, MI, USA

Copyright © 2018 by author and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

Received: February 5, 2018; Accepted: March 17, 2018; Published: March 20, 2018


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.


Simulation, Metallic Glass, Solidification, Toughness

1. 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 microstructure primarily explaining dendrite tip temperature and dendrite tip radius as a function of growth rate/dendrite tip velocity based on dendrite tip stability theory.

2. 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 (ko) 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 contradiction to assumption 5 above in basic KGT model and provides the basis of development of present model to explain microstructural evolution in detail.

There are three main velocities which are of interest here (Figures 1(a)-(c)) when a beam of high energy (electron or laser) travels on surface of specimen in additive manufacturing setup (Figure 2).

1) Moving heat source velocity (Vb), 2) solidification front velocity (Vs) and 3) dendrite tip velocity (Vhkl)

Model comprises of three main parts and separate mathematical expression is developed for each segment.


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 (nmax) and grain density (n(ΔT)). Maximum nucleation density may be obtained by integral of nucleation distribution from zero undercooling to infinite undercooling.

n max i = 0 d n d Δ T Δ T (1)

Figure 1. Schematic showing relationship between different velocity components (reprinted with permission from Springer) [50].

Figure 2. Schematic of melt pool formation in additive manufacturing (EBM) and velocity of heat source (reprinted with permission from Springer) [50].

Similarly, grain density is given by following equation

n ( Δ T ) = 0 Δ T n max Δ T σ exp [ 1 2 ( Δ T Δ T N Δ T σ ) ] d Δ T (2)

where ΔTn 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 (pv) as described by Prof. Rappaz in his famous article [46].

p v r (3)

i-e if at any instant of time t, pv exceeds r, nucleation will occur. p v = δ n v V C A where δnv = grain density increase and VCA = one cell volume (measure by noting all dimensions of a 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(b). Typically, cell dimensions are in nanometers while typical number of cells reported in previous calculations are 30,000 [52].

Dendrite growth orientation

Second part of the model deals with determination of 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 held responsible 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

(a) (b)

Figure 3. (a) Schematic of growth in one cell; (b) Schematic diagram of movement of dendrite tip in a grid of cells represented by change of state index of each cell (reprinted with permission from Taylor and Francis Group) [48].

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 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. First two angles describes the growth and orientation of main [001] trunk while third angle ψ describes the orientation of [100] and [010] secondary branches. A schematic showing their arrangement is shown in Figure 4.

The probability d p ( θ , φ , ψ ) that a newly nucleated grain has its main trunk orientation in the range [ θ , θ + d θ ] and [ φ , φ + d φ ] and one of its set of secondary branches within the orientation [ ψ , ψ + d ψ ] is given by [52]

d p ( θ , φ , ψ ) = A sin θ d θ d φ d ψ (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 concerns the determination of supersaturation of individual elements in multicomponent alloy (BMGMC) systems. This approach arises from the notion that in contrast to conventional castings in which undercooling related with thermal diffusion, attachment kinetics and curvature is small, in multicomponent systems (undergoing additive manufacturing treatment) 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 Ωi is a function of Peclet number, Pei

Figure 4. Schematic of crystallographic orientation of new nucleated grain as defined by three Euler Angles, θ, φ and ψ (reprinted with permission from Springer) [52].

Ω = I v ( P e ) (5)

Ω i = I v ( P e i ) (6)

P e i = R V 2 D i (7)

Putting in [5]

Ω i = I v ( R V 2 D i ) (8)

Ω i = P e i e P e i E 1 ( P e i ) = P e i e P e i P e i e u u d u (9)

Ω i = R V 2 D i e R V 2 D i E 1 ( P e i ) = R V 2 D i e R V 2 D i R V 2 D i e u u d u (10)


Ω i = c i * c o , i c i * ( 1 k i ) (11)

where, c i * = concentration of constituent i in liquid at dendrite tip (to be found), c o , i = initial concentration of constituent i, k i = partition coefficient for this constituent i (velocity dependent).

Comparing (10) and (11)

c i * c o , i c i * ( 1 k i ) = R V 2 D i e R V 2 D i R V 2 D i e u u d u (12)


c i * = c o , i 1 ( 1 k i ) ( I v ( P e i ) ) (13)

k i = k o + ( a o V D i ) 1 + ( a o V D i ) (14)

where, ao = length scale related to interatomic distance and is estimated to be between 0.5 - 5 nm and

D i = D o e ( Q R g T * ) (15)

where, D o = Proportionality constant, Q = Activation energy, Rg = Gas constant, T* = Tip temperature calculated by iterative method (described below) [51].

In a linearized phase diagram

T * = T L + i = 1 2 m i ( c i * c o , i ) 2 Γ R (16)

where, mi = slope of liquidus surface with respect to constituent i

m i = T L c i (17)

TL = Liquidus temperature for initial alloy composition, Γ = Gibbs-Thomson coefficient, 2 Γ R 0 (negligible) (under normal solidification conditions), 2 Γ R = 1 (under rapid solidification conditions (additive manufacturing)).

Another term can be generated from linearized phase diagram known as fictitious melting point of pure constituent [51].

T m = T L i = 1 2 m i c o , i (18)

Using Equations (13) and (18), Equation (16) becomes

T * = T m + i = 1 2 m i c o , i 1 ( 1 k i ) ( I v ( P e i ) ) 2 Γ R (19)

where, R = Dendrite tip radius, V = Dendrite tip velocity.

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. 100 iterative cycles are used to generate homogeneous and normalized data based upon best software engineering practice. In general while writing the program, reading it and executing it, Ω depends on I v ( P e i ) and c i * , c i * depends on ki and I v ( P e i ) , ki depends on D i , D i depends on Tip temperature and finally, Tip temperature depends on c i * . 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.

Finally, total undercooling (ΔT) is related to supersaturation (Ω) by [52]

Δ T i = m i c o , i [ 1 1 1 Ω i ( 1 k i ) ] (20)

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 and Sekerka [60], Langer and Mueller-Krumbhaar [61], KGT [47] and BLL [62] models).

Accordingly, one has for a ternary system

R = 2 π ( Γ i = 1 2 m i G c , i ε c ( P e i ) G ) 1 / 2 (21)

where, Gc,i = solute gradient of constituent i in the liquid near tip which can be written as

G c , i = υ D i C i * ( 1 k i ) (22)

Figure 5. Schematic flow chart describing working of model.

G = Average thermal gradient near tip and ε c ( P e i ) = f ( P e i ) = 1 (low speed/ low Pei).

Putting values of P e i = R V 2 D i and Equation (22) in Equation (21) [51]

4 π 2 Γ R 2 + 2 R i = 1 2 P e i m i C o , i ( 1 k i ) ξ c ( P e i ) [ 1 ( 1 k i ) I v ( P e i ) ] + G = 0 (23)

where, Γ = Gibbs-Thomson coefficient, mi = the slope of liquidus, ξ c = a function of the Peclet number and segregation coefficient, G = Thermal gradient.

3. 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. 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 (e.g. 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

Table 1. Parameters used in the calculations of dendritic growth for Zr based BMGMCs.

front. Performing simulations employing data from other elements in a BMGMC system e.g. Cu, Al, Co, Ni is left for reader as an exercise and home work task.

4. Experimental

BMGMC samples were produced 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 furnace, they are wrapped in an aluminum foil which not only protects the powders but also serve as alloying element in sufficient quantity in original mix. This Aluminum foil wrapped toffee is placed in horizontal slot in the water cooled copper hearth of 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) [66] in Additive Manufacturing setup. Model theory was developed and tailored to describe microstructure evolution during both processes.

5. 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 to a variety of alloys systems under various conditions. Following results and graphs have been generated after writing script of solidification code and running it in MATLAB®.

Effect of heat dissipation on dendrite tip velocity

Figure 6 is the plot of evolution of solutal Peclet number of different elements of bulk metallic glass matrix composite 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

Figure 6. Plot of solutal Peclet number Vs dendrite tip velocity for individual elements of BMGMC system.

heat loss 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 presence and movement of Zirconium through fluid. Overall large atomic weight remains as the most favorable 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. It also does not have very high thermal conductivity thus not a lot of heat gets dissipated to surrounding because of it and again it reflects its inability to create low viscosity zone around it which may facilitate it to gain high speeds. 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 or spheroid. 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 is needed to measure this.

Effect of tip temperature on dendrite tip velocity

Below graph (Figure 7) shows relationship between dendrite tip velocity as a function of dendrite tip temperature for Zirconium only in the systems considering it to be major alloying element. Graph shows three different regions which are distinctive of velocity evolution of Zr in the system as temperature change. There is a sharp increase in Ttip at slow Vf 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. Here, since only Zirconium is under consideration, so the large atomic size and weight of Zr also contribute towards increase of this value. 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. [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 Ttip

Figure 7. Plot of dendrite tip temperature Vs dendrite tip velocity for Zr in BMGMC system.

with increase in velocity is observed. This happens as Zr dendrite gain speed; it starts creating regions of very high surface/volume ratio around it in the body of melt and at the tip of advancing dendrite. In other words fast moving dendrites become source of high heat dissipation. Owing to this, temperature drops. Discrepancy between experimental values and simulation results is due to lesser number of iterative cycles which create low convergence. Hence, a temperature field with a deviation from actual values is observed. This anomaly could be removed by increasing iterations which enhances simulation efficiency and helps in model refinement. A decreasing trend in dendrite tip temperature is also expected at higher thermal gradients (G) as this will provide a large sink for the heat to get dissipated away from liquid melt pool which eventually will lower the tip temperature even at low growth velocities.

Effect of tip temperature on dendrite tip radius

The effect of dendrite tip temperature on the evolved radius is shown in Figure 8 below for Zirconium. This is a typical graph describing transient nature transport phenomena in metallurgical unit operations. At very high melting points liquid and solid fronts are immiscible with each other. Hence, no sharp interface/boundary is observed. As the alloy cools in mold / liquid melt pool of additive manufacturing, a distinct interface describing evolution of dendrite starts appearing. This effectively can be described in terms of curvature known as dendrite tip radius. The graph shows the behavior of this radius at different values of temperature and clearly shows that this radius decreases with the rise of temperature leading towards its complete disappearance near melting point. At lower tip temperatures, more time is available to solute to diffuse out of solvent to accumulate in the form of a mass which evolves in the form of radius of solidifying front. This radius can change from planar to dendrite depending on

Figure 8. Plot of dendrite tip temperature Vs dendrite tip radius for Zr in BMGMC system.

thermal gradient which is final description of microstructural phenomena in solidifying 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 Pei) on Pei

Below graph (Figure 9) 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 describes 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 [64] [65] [69] [70] [71] [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

Figure 9. Plot of ξ Vs Peclet number for various elements in BMGMC system.

Figure 10. Plot of evolution of partition coefficient with temperature for Zr in BMGMC system at constant thermal gradient.

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. 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 observed in experimental studies, it can be witnessed in this simulation as well that solute field around dendrite tip decreases almost homogeneously and linearly 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 [47] [60]. 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 equiaxed-columnar transition followed by almost linear decrease 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 quicker 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.


Figure 11. (a) Plot of evolution of dendrite tip growth velocity as a function of supersaturation for all alloying elements in Zr65Cu15Al10Ni10 BMGMC system at constant thermal gradient; (b) Plot of evolution of dendrite tip growth velocity as a function of supersaturation for Ni in Zr65Cu15Al10Ni10 BMGMC system at constant thermal gradient.

6. Salient Features of Model

Primarily, present model is formulated as a result of combination of following [47] [50] [51] [52] studies

1) Gandin, C.-A., M. Rappaz, and R. Tintillier, Three-dimensional probabilistic simulation of solidification grain structures: Application to superalloy precision castings. Metallurgical Transactions A, 1993. 24(2): p. 467-479.

2) Rappaz, M., et al., Analysis of solidification microstructures in Fe-Ni-Cr single-crystal welds. Metallurgical Transactions A, 1990. 21(6): p. 1767-1782.

3) Rappaz, M., et al., Development of microstructures in Fe-15Ni-15Cr single crystal electron beam welds. Metallurgical Transactions A, 1989. 20(6): p. 1125-1138.

4) Kurz, W., B. Giovanola, and R. Trivedi, Theory of microstructural development during rapid solidification. Acta Metallurgica, 1986. 34(5): p. 823-830.

5) Zhang, J., et al. Probabilistic simulation of solidification microstructure evolution during laser-based metal deposition. In Proceedings of 2013 Annual International Solid Freeform Fabrication Symposium—An Additive Manufacturing Conference. 2013.

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.

3) An effort was made to remove/reduce error by use of iteration based approach to refine model.

4) Programming of model was done in MATLAB®—not done elsewhere previously.

5) Temperature dependent properties (transient heat transfer conditions) were used.

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.

8) Peclet number (Pe) & ξ are not taken as constant like previous studies [Bobadilla, M., J. Lacaze, and G. Lesoult, Journal of Crystal Growth, 1988. 89(4): p. 531-544] in which it is assumed

a) ξ = 1 (low growth rate) (low Pe)

b) ξ = 0 (very fast cooling rate—typical Additive Manufacturing conditions)

9) 2Γ/R = 1 (high velocity AM conditions).

10) New relation for dendrite tip temperature was developed.

NOTE: No physical microstructure was either reported previously or tried in present approach. Physical simulation is “not possible” at deterministic stage as it is numerical/analytical model which defines parameters for next stage probabilistic studies only.

7. Conclusions

Following conclusions may be drawn out of this study.

a) There is significant effect of initial metal temperature, composition, type of alloying elements, temperature gradient and thermo-physical properties on final microstructure developed as a result of heat and mass transfer phenomena.

b) Determination of supersaturation of individual elements yield best possible strategy for its correlation with superimposition of solute field around each dendrite tip.

c) Determination of dimensionless solutal Peclet number is the main factor responsible for accurate quantitative prediction of microstructure in solidifying alloys.

d) Dependence and evolution of ξ on, and with respect to solutal Peclet number is decisive in explaining transient nature transport phenomena in additive manufacturing processes.

e) Employment of iterative process helps in refining the model and generates accurate results.

f) Final microstructure evolution is expressed in the form of dendrite tip temperature and dendrite tip radius as a function of growth rate/dendrite tip velocity and must be carefully measured.

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.


This work is partially supported by scholarship provided by RMIT University for tuition and living of author. Further, author would like to acknowledge moral support and encouragement provided by his primary supervisor (Dr. Dong Qiu) and Prof. Milan Brandt throughout the work and helpful discussion and constructive criticism by Prof. Mark Easton. The authorization of gracious use of experimental facilities during short term stay at CSIRO and above average moral support by Dr. Daniel East is also gratefully acknowledged.

Competing Financial Interests

The author(s) declare no competing financial interests.

Cite this paper

Rafique, M.M.A. (2018) Simulation of Solidification Parameters during Zr Based Bulk Metallic Glass Matrix Composite’s (BMGMCs) Additive Manufacturing. Engineering, 10, 85-108.


  1. 1. Klement, W., Willens, R.H. and Duwez, P.O.L. (1960) Non-Crystalline Structure in Solidified Gold-Silicon Alloys. Nature, 187, 869-870.

  2. 2. Hays, C.C., Kim, C.P. and Johnson, W.L. (2000) Microstructure Controlled Shear Band Pattern Formation and Enhanced Plasticity of Bulk Metallic Glasses Containing in situ Formed Ductile Phase Dendrite Dispersions. Physical Review Letters, 84, 2901-2904.

  3. 3. Johnson, W.L. (1999) Bulk Glass-Forming Metallic Alloys: Science and Technology. MRS Bulletin, 24, 42-56.

  4. 4. Ashby, M.F. and Greer, A.L. (2006) Metallic Glasses as Structural Materials. Scripta Materialia, 54, 321-326.

  5. 5. Flores, K.M. and Dauskardt, R.H. (1999) Local Heating Associated with Crack Tip Plasticity in Zr-Ti-Ni-Cu-Be Bulk Amorphous Metals. Journal of Materials Research, 14, 638-643.

  6. 6. Eckert, J., et al. (2007) Mechanical Properties of Bulk Metallic Glasses and Composites. Journal of Materials Research, 22, 285-301.

  7. 7. Das, J., et al. (2009) Designing Bulk Metallic Glass and Glass Matrix Composites in Martensitic Alloys. Journal of Alloys and Compounds, 483, 97-101.

  8. 8. Das, J., et al. (2005) “Work-Hardenable” Ductile Bulk Metallic Glass. Physical Review Letters, 94, 205501.

  9. 9. Choi-Yim, H. and Johnson, W.L. (1997) Bulk Metallic Glass Matrix Composites. Applied Physics Letters, 71, 3808-3810.

  10. 10. Cheng, J.L. and Chen, G. (2013) \ Glass Formation of Zr-Cu-Ni-Al Bulk Metallic Glasses Correlated with L → Zr2Cu + ZrCu Pseudo Binary Eutectic Reaction. Journal of Alloys and Compounds, 577, 451-455.

  11. 11. Chen, M. (2011) A Brief Overview of Bulk Metallic Glasses. NPG Asia Materials, 3, 82-90.

  12. 12. Chen, M. (2008) Mechanical Behavior of Metallic Glasses: Microscopic Understanding of Strength and Ductility. Annual Review of Materials Research, 38, 445-469.

  13. 13. Chen, H.S. (1974) Thermodynamic Considerations on the Formation and Stability of Metallic Glasses. Acta Metallurgica, 22, 1505-1511.

  14. 14. Akihisa, I., et al. (1988) Glass Transition Behavior of Al-Y-Ni and Al-Ce-Ni Amorphous Alloys. Japanese Journal of Applied Physics, 27, L1579.

  15. 15. Johnson, W.L., et al. (2011) Beating Crystallization in Glass-Forming Metals by Millisecond Heating and Processing. Science, 332, 828-833.

  16. 16. Jiang, M.Q., et al. (2010) Fractal in Fracture of Bulk Metallic Glass. Intermetallics, 18, 2468-2471.

  17. 17. Qiao, J., Jia, H. and Liaw, P.K. (2016) Metallic Glass Matrix Composites. Materials Science and Engineering: R: Reports, 100, 1-69.

  18. 18. Schroers, J. (2010) Processing of Bulk Metallic Glass. Advanced Materials, 22, 1566-1597.

  19. 19. Greer, A.L. (2010) Materials Science: A Cloak of Liquidity. Nature, 464, 1137-1138.

  20. 20. Greer, A.L. (1995) Metallic Glasses. Science, 267, 1947-1953.

  21. 21. Yi, J., et al. (2016) Glass-Forming Ability and Crystallization Behavior of Al86Ni9La5 Metallic Glass with Si Addition. Advanced Engineering Materials, 18, 972-977.

  22. 22. Cheng, Y.Q., Sheng, H.W. and Ma, E. (2008) Relationship between Structure, Dynamics, and Mechanical Properties in Metallic Glass-Forming Alloys. Physical Review B, 78, 014207.

  23. 23. Sarac, B. (2015) Microstructure-Property Optimization in Metallic Glasses. Springer, Berlin.

  24. 24. Greer, A.L. (2011) Metallic Glasses: Damage Tolerance at a Price. Nature Materials, 10, 88-89.

  25. 25. Gu, X.W., et al. (2014) Mechanisms of Failure in Nanoscale Metallic Glass. Nano Letters, 14, 5858-5864.

  26. 26. Schroers, J. and Johnson, W.L. (2004) Ductile Bulk Metallic Glass. Physical Review Letters, 93, 255506.

  27. 27. Schuh, C.A., Hufnagel, T.C. and Ramamurty, U. (2007) Mechanical Behavior of Amorphous Alloys. Acta Materialia, 55, 4067-4109.

  28. 28. Donovan, P.E. and Stobbs, W.M. (1981) The Structure of Shear Bands in Metallic Glasses. Acta Metallurgica, 29, 1419-1436.

  29. 29. Dodd, B. and Bai, Y. (2012) Adiabatic Shear Localization: Frontiers and Advances. Elsevier, Amsterdam.

  30. 30. Gao, Y.F., et al. (2011) On the Shear-Band Direction in Metallic Glasses. Acta Materialia, 59, 4159-4167.

  31. 31. Greer, A.L., Cheng, Y.Q. and Ma, E. (2013) Shear Bands in Metallic Glasses. Materials Science and Engineering: R: Reports, 74, 71-132.

  32. 32. Jiang, M.Q., Wang, W.H. and Dai, L.H. (2009) Prediction of Shear-Band Thickness in Metallic Glasses. Scripta Materialia, 60, 1004-1007.

  33. 33. Leng, Y. and Courtney, T.H. (1991) Multiple Shear Band Formation in Metallic Glasses in Composites. Journal of Materials Science, 26, 588-592.

  34. 34. Hajlaoui, K., et al. (2007) Unusual Room Temperature Ductility of Glassy Copper-Zirconium Caused by Nanoparticle Dispersions That Grow during Shear. Materials Science and Engineering: A, 449-451, 105-110.

  35. 35. Zhang, Y. and Greer, A.L. (2007) Correlations for Predicting Plasticity or Brittleness of Metallic Glasses. Journal of Alloys and Compounds, 434-435, 2-5.

  36. 36. Lewandowski, J., Wang, W.-H. and Greer, A. (2005) Intrinsic Plasticity or Brittleness of Metallic Glasses. Philosophical Magazine Letters, 85, 77-87.

  37. 37. Kruzic, J.J. (2016) Bulk Metallic Glasses as Structural Materials: A Review. Advanced Engineering Materials, 18, 1308-1331.

  38. 38. Nishiyama, N., et al. (2012) The World’s Biggest Glassy Alloy Ever Made. Intermetallics, 30, 19-24.

  39. 39. Schroers, J. (2005) The Superplastic Forming of Bulk Metallic Glasses. JOM, 57, 35-39.

  40. 40. Guo, G.-Q., et al. (2015) Detecting Structural Features in Metallic Glass via Synchrotron Radiation Experiments Combined with Simulations. Metals, 5, 2093-2108.

  41. 41. Guo, G.-Q., et al. (2015) How Can Synchrotron Radiation Techniques Be Applied for Detecting Microstructures in Amorphous Alloys? Metals, 5, 2048-2057.

  42. 42. Zimmermann, G., et al. (2011) Investigation of Columnar-to-Equiaxed Transition in Solidification Processing of AlSi Alloys in Microgravity—The CETSOL Project. Journal of Physics: Conference Series, 327, 012003.

  43. 43. Zu, F.-Q. (2015) Temperature-Induced Liquid-Liquid Transition in Metallic Melts: A Brief Review on the New Physical Phenomenon. Metals, 5, 395-417.

  44. 44. Kim, D.H., et al. (2013) Phase Separation in Metallic Glasses. Progress in Materials Science, 58, 1103-1172.

  45. 45. Ott, R.T., et al. (2005) Micromechanics of Deformation of Metallic-Glass-Matrix Composites from in situ Synchrotron Strain Measurements and Finite Element Modeling. Acta Materialia, 53, 1883-1893.

  46. 46. Rappaz, M. and Gandin, C.A. (1993) Probabilistic Modelling of Microstructure Formation in Solidification Processes. Acta Metallurgica et Materialia, 41, 345-360.

  47. 47. Kurz, W., Giovanola, B. and Trivedi, R. (1986) Theory of Microstructural Development during Rapid Solidification. Acta Metallurgica, 34, 823-830.

  48. 48. Wei, Y.H., et al. (2007) Numerical Simulation of Columnar Dendritic Grain Growth during Weld Solidification Process. Science and Technology of Welding and Joining, 12, 138-146.

  49. 49. Rappaz, M. and Blank, E. (1986) Simulation of Oriented Dendritic Microstructures Using the Concept of Dendritic Lattice. Journal of Crystal Growth, 74, 67-76.

  50. 50. Rappaz, M., et al. (1989) Development of Microstructures in Fe-15Ni-15Cr Single Crystal Electron Beam Welds. Metallurgical Transactions A, 20, 1125-1138.

  51. 51. Rappaz, M., et al. (1990) Analysis of Solidification Microstructures in Fe-Ni-Cr Single-Crystal Welds. Metallurgical Transactions A, 21, 1767-1782.

  52. 52. Gandin, C.-A., Rappaz, M. and Tintillier, R. (1993) Three-Dimensional Probabilistic Simulation of Solidification Grain Structures: Application to Superalloy Precision Castings. Metallurgical Transactions A, 24, 467-479.

  53. 53. Zhang, J., et al. (2013) Probabilistic Simulation of Solidification Microstructure Evolution during Laser-Based Metal Deposition. 24th Annual International Solid Freeform Fabrication Symposium—An Additive Manufacturing Conference, Austin, 2013, TX, 739-748.

  54. 54. Zhou, X., et al. (2016) Simulation of Microstructure Evolution during Hybrid Deposition and Micro-Rolling Process. Journal of Materials Science, 51, 6735-6749.

  55. 55. Gu, C., et al. (2017) A Three-Dimensional Cellular Automaton Model of Dendrite Growth with Stochastic Orientation during the Solidification in the Molten Pool of Binary Alloy. Science and Technology of Welding and Joining, 22, 47-58.

  56. 56. Nastac, L. (1999) Numerical Modeling of Solidification Morphologies and Segregation Patterns in Cast Dendritic Alloys. Acta Materialia, 47, 4253-4262.

  57. 57. Laurentiu, N. and Doru, M.S. (1997) Stochastic Modelling of Microstructure Formation in Solidification Processes. Modelling and Simulation in Materials Science and Engineering, 5, 391.

  58. 58. Von Neumann, J. and Burks, A.W. (1996) Theory of Self-Reproducing Automata. University of Illinois Press, Urbana, IL.

  59. 59. Reuther, K. and Rettenmayr, M. (2014) Perspectives for Cellular Automata for the Simulation of Dendritic Solidification—A Review. Computational Materials Science, 95, 213-220.

  60. 60. Mullins, W.W. and Sekerka, R.F. (1964) Stability of a Planar Interface during Solidification of a Dilute Binary Alloy. Journal of Applied Physics, 35, 444-451.

  61. 61. Langer, J.S. and Müller-Krumbhaar, J. (1977) Stability Effects in Dendritic Crystal Growth. Journal of Crystal Growth, 42, 11-14.

  62. 62. Bobadilla, M., Lacaze, J. and Lesoult, G. (1988) Influence des conditions de solidification sur le déroulement de la solidification des aciers inoxydables austénitiques. Journal of Crystal Growth, 89, 531-544.

  63. 63. Rafique, M.M.A., Qiu, D. and Easton, M. (2017) Modeling and Simulation of Microstructural Evolution in Zr Based Bulk Metallic Glass Matrix Composites during Solidification. MRS Advances, 2, 3591-3606.

  64. 64. Wu, K., Li, R. and Zhang, T. (2013) Crystallization and Thermophysical Properties of Cu46Zr47Al6Co1 Bulk Metallic Glass. AIP Advances, 3, 112115.

  65. 65. Yamasaki, M., Kagao, S. and Kawamura, Y. (2005) Thermal Diffusivity and Conductivity of Zr55Al10Ni5Cu30 Bulk Metallic Glass. Scripta Materialia, 53, 63-67.

  66. 66. Yang, G., et al. (2012) Laser solid forming Zr-based bulk metallic glass. Intermetallics, 22, 110-115.

  67. 67. Flemings, M.C. (1974) Solidification Processing. McGraw-Hill, New York.

  68. 68. Yang, L., et al. (2009) Nanoscale Solute Partitioning in Bulk Metallic Glasses. Advanced Materials, 21, 305-308.

  69. 69. Mills, K.C. (2002) Front Matter A2—Recommended Values of Thermophysical Properties for Selected Commercial Alloys. Woodhead Publishing, Cambridge, iii.

  70. 70. Grimvall, G. (1999) Front Matter A2—Thermophysical Properties of Materials. Elsevier Science B.V., Amsterdam, iii.

  71. 71. Valencia, J.J. and Quested, P. (2001) Thermophysical Properties. Modeling for Casting and Solidification Processing, 189.

  72. 72. Choy, C.L., et al. (1991) Thermal Conductivity of Amorphous Alloys above Room Temperature. Journal of Applied Physics, 70, 4919-4925.