Bulk metallic glass matrix composites  are unique materials of future having superior mechanical properties which surpasses any other material known till date. Their high hardness, strength and elastic strain limit make them potential candidates for future extreme engineering applications  - . However, they suffer from little or no toughness under tensile loading. In fact, they exhibit tensile softening. This is a detrimental property which renders them unusable for any practical use. This can be overcome by various means. One of which is in-situ introduction of primary phase ductile particles (e.g. spheroidal B2-CuZr intermetallic) during solidification which serves as an obstacle to movement of shear band thus pin them and increase their toughness. This phenomenon is difficult to control, and a lot of rigorous experimentations are required before conclusive results could be drawn. In present study, which is continuation of author’s previous work , a probabilistic model of nucleation and growth   is developed to supplement deterministic model to describe this phenomena in these alloys without recurrence to experimentation. Much research has been carried out previously on modelling and simulation of solidification phenomena in various processes employing cellular automata or modified cellular automata combined with finite element such as modified additive manufacturing (HDMR , LAMP  on 316L SS , two dimensional CAFE , cladding  ), selective laser melting (SLM) using CAFE    , CAPF, CAFVM  and modified CAFE  but none has been conducted on the use of cellular automata method on bulk metallic glasses or their composites. On modelling and simulation of bulk metallic glasses themselves; only two studies have been reported   but again they are not aimed at explaining microstructure evolution. At present scenario, there is no single hybrid/combined model which explains phenomena of heat transfer (liquid melt pool formation as a result of laser—matter interaction and its evolution—solidification) coupled with mass transfer (nucleation and growth (solute diffusion  and capillary action driven)) to predict microstructure and grain size in bulk metallic glass matrix composites as melt cools in additive manufacturing liquid melt pool. Present study is one step forward to answer second (i.e. microscale mass transfer) question. Aim is to model nucleation of primary phase ductile crystalline particles in glassy matrix. Constitutive equations governing microscale mass transfer and diffusion phenomena in individual cell of cellular automata grid are developed which describe evolution of solid fraction and solid fraction increment in each cell as a function of time. Together with neighbourhood transition rules, this will give visualization of microstructure evolved in each cell. Model is described as follows.
In order to formulate model in present case, following assumptions are drawn.
1) Model is developed under steady state conditions.
2) Regular square shape is chosen as shape of one cell out of square, hexagon or triangle.
3) Decentred rectangular shape is chosen to counteract mesh anisotropy.
4) One of crystallographic orientations (e.g. one of <100> directions for fcc metals) of grain stands perpendicular to simulation plane. This direction is diagonal of square and labelled simply as <10> direction (Figure 1).
5) Incubation time (time in which active side branches are influenced by solute field of their neighbours) is neglected.
6) Nucleation radius of grains is neglected.
A summary of steps to be followed by cellular automaton simulation is described below;
1) Choose regular shape for populating cellular automaton (CA) grid (e.g. square).
2) Define decentred square envelope.
3) Grow grain in it under steady state conditions.
4) Grain A nucleates at centre of cell v (it has inherent/intrinsic misorientation of its crystallographic axis to that of CA axis).
NOTE: Misorientation in present case is intrinsic feature (i.e. it is caused by features of grains). To remove it, either choose different lattice (hexagonal, octagonal) “or” refine the algorithm.
In reality, this growth of dendritic front has to be tackled at the level of secondary dendritic arm spacing. This may be done by two ways;
1) Introducing dendrite tip correction or
2) Use of 2D rectangular algorithm  at secondary dendrite arm spacing level.
Both these approaches are difficult to implement and practice as “first” is not a CA model at all while “second” is difficult to extend to 3D as rectangle would correspond to non-regular octahedron shape. That’s why, a new approach known as 2D decentred square growth algorithm is developed and adopted (present study) (Figure 2).
Figure 1. Schematic of crystallographic orientation of one direction.
2.2. Procedure of Cellular Automaton Model
Cellular automaton model consists of  discretizing
1) Simulated area into finite cells while
2) Time into small time steps
At a certain time, a cell has a “state” that can be described by number of variables. Most important of which are;
1) Temperature and
2) Solute concentration
During simulation, the state of each cell is determined by the state of its nearest neighbors through a transition rule known as neighborhood transition rule. Two types of neighborhood transition rules are commonly employed.
1) Moore’s law (8 neighborhood)
2) Von Neuman’s law (4 neighborhood)
In traditional CA models, the state of all cells is simultaneously updated. Cell switch states depending on the outcome of the transition rules. In present case, Von Neumann neighborhood is used (4 nearest neighbors). CA rule of Von Neumann neighborhood can be determined by
where is state of cell, which has mark number (i, j) at the time t; is the evolvement rule and it is vital that this rule is carefully chosen to reflect the objective physics of the simulated phenomena. It may be chosen out of 256 rules.
Now, solute diffusion is explained for spheroidal B2—CuZr intermetallic in glassy matrix. This solute diffusion happens in two steps a) Solute diffusion in single cell and b) Solute diffusion between cells.
2.2.1. Solute Diffusion in Single Cell
At a particular time, each CA cell may have three possible states; liquid (solid fraction is 0), mushy (solid fraction between 0 and 1), or solid (solid fraction is 1). Cell state changes from liquid to mushy when nucleation occurs. Solute diffusion in a cell containing “liquid” and “solid” are described by  ;
Figure 2. Schematic of modified 2D decentred square CA growth algorithm .
respectively, where CL is solute concentration in liquid, CS is same in solid and DL is solute diffusion coefficient in liquid and DS is same in solid. Similarly, solute diffusion in a mushy cell is given by below equations
where and are the diffusion coefficients of the liquid phase and solid phase, respectively. is the solid fraction, is the equivalent concentration, and is equivalent diffusion coefficient.
1) Fraction Solid in Single Cell
To calculate solid fraction increment, first growth velocities of solid liquid interface in x and y directions are calculated  . This is briefly described below;
once calculated, fraction solid evolution is calculated by
where a = mesh size (uniform and constant) for both x and y direction and δt is time step.
In a single CA cell, solid fraction can be calculated using
where = solid fraction at previous time step
Time step itself may be calculated by
where Δx is cell size, and Vmax is maximum growth velocity of cell during each time step.
Note: solid fraction (fs) and solid fraction increment ( ) are two different things and should not be confused. Former means actual numerical increase in amount of solid evolved in a cell during time step δt while later is difference of value between old and new.
2) Mean Curvature of Solid Liquid Interface
By a similar procedure, mean curvature of solid liquid (S/L) interface k may be determined
where N is the neighboring cells and is the solid fraction
where is liquidus temperature at , is local temperature and are constitutional and curvature under cooling temperatures respectively.
Putting in (2)
where = initial concentration, = liquidus slope, = Gibbs – Thomson coefficient, k = mean curvature of S/L interface.
2.2.2. Solute Diffusion between Cells
Procedure adopted to calculate solute diffusion between cells is similar to the one adopted for calculation of solute diffusion in a single cell. Primarily, diffusion occurs owing to concentration gradient. This gradient is created by the release of solute in both liquid and solid. This may be summarized in following four steps ;
Step 1: During this equilibrium at solid liquid interface is determined by
where is solute concentration in solid at solid/liquid interface. is counterpart in liquid and k is partition coefficient. It is utmost important to determine k in a proper manner preferably keeping in view transient conditions such that actual process conditions are simulated.
Step 2: Partitioning of solute in liquid is determined. This may be done by the use of equation 1 by replacing all subscripts by L.
Step 3: Partitioning of solute in solid is determined by again using Equation 1 without replacing any subscript or term.
Step 4: Partitioning of solute in growing cell may be determined as follows;
Neighbors of a cell could be anything between liquid, solid or growing. For a liquid cell, it could be solid or growing. For a solid cell, it could be liquid or growing and for a growing cell, it could be liquid, solid or even growing. When the computing center is a growing cell, regardless of state (i.e. solid, liquid or growing) of the neighbor of growing cell, the solute concentration of growing cell may be expressed as an equivalent concentration CE, which may be defined as;
when fs = 1, CE = CSand when fs = 0, CE = CL
Thus, diffusion equation may be expressed as;
where CE is equivalent concentration defined by above equation and DE is equivalent diffusion coefficient . This is schematically represented in Figure 3 below.
Figure 3. Schematic diagram of growing in a cell.
Since thermal calculations are done at finite element mesh. There is need to relate finite element mesh to cellular automaton mesh. This is done by defining interpolation equation .
Physical domain in three-dimensional space is divided into cells of equal size 0.1 µm × 0.1 µm × 0.1 µm. Each cell is characterized by three variables.
1) State (solid, mushy cell or liquid)
2) Temperature and
Temperature values are calculated on the macro nodes and interpolation equation is needed for use in CA cell. This equation may be described as
where d = distance between 8 nearest nodes (moors model) in finite element model and center of cellular automaton cell (Figure 4).
At the start of simulation, all cells states are considered solid and all crystallographic variables are zero (non-active cells). As simulation starts, temperature of each cell is given by thermal simulation performed on finite element mesh nodes.
If temperature of cell passes the melting temperature (Tm), its state changes to liquid. If it does not pass Tm, it stays as solid and if it stays at Tm, it becomes mushy.
Calculation of Nucleation Density of New Cells and Assigning p Number
Solute diffusion is calculated only for cells which experience melting and solidification (active cells). Liquid cell can again become solid if one of following condition happens. Nucleation or Capturing by another active—solid CA cells.
In each time step density of new nuclei  is determined by equation
This equation is tailored to take into account nucleation occurring at
Figure 4. Representation of interpolation in one cell.
1) Molten pool wall and
2) Bulk liquid
Nucleation densities are then multiplied by the total number of cells for Molten pool wall and bulk liquid in order to calculate number of new nucleation sites
where = Number of new nucleation sites at the molten pool wall and = Total number of liquid cells at the molten pool wall. Similarly
where = Number of new nucleation sites inside bulk liquid and = Total number of liquid cells inside bulk liquid.
After calculating the number of new nucleation sites, a random number (0 < P < 1) is assigned to each liquid cell. The nucleas is only formed in a cell i which is located at molten pool wall if following condition is satisfied.
Similarly, this condition must be satisfied for cells located in bulk liquid. This ensures that state of cell i changes from liquid to solid as nucleation occurs whether it is at molten pool boundary or inside bulk liquid.
Another important parameter necessary to describe probabilistic phenomena is assigning of “preferential growth orientation”. Again, this is determined by a random process. Certain numbers of growth orientations are chosen in three-dimensional spherical space around each nucleated cell as possible growth orientation. State of cell i is chosen randomly from these orientations as nucleation occurs. Different growth velocities are attained by different grains because they have different preferential growth direction. A parameter known as orientation weight coefficient ( ) account for this effect. The growth length of cell i with regards to its liquid neighbor j (one of possible 26 cells) at time tc can be calculated by:
where is the orientation weight coefficient related to angle between preferential growth direction of cell i and vector from cell i linking to cell j (called vector ). Orientation weight coefficient itself may be calculated by ;
where Xw, Yw and Zw may be calculated by solving below matrix
where and are direction cosines of  and  and  dendrite arms relating to the coordinate x, y, z axes respectively and are direction cosines of vector relating to the coordinate x, y, z axes respectively.
In detail, is growth length between two points i and j (i has solid state and j has liquid state) when crystallographic orientation of cell i is considered. Neighbor j is captured by cell i if below condition is satisfied.
where if j is one of 6 nearest neighbors, if j is one of twelve second—nearest neighbors, and if j is one of eight third—nearest neighbors.
If cell j is captured by i, it will be assigned same orientation as cell i. Interested reader may find further detail in .
In present study, a cellular automata method is described for describing nucleation and growth of primary ductile phase particles in glassy matrix in bulk metallic glass matrix composites. Model is built on constitutive equations of microscale mass transfer and diffusion phenomena. Model is built assuming steady state conditions. However, use of transient conditions is recommended and as it brings model and thus simulation results to actual values. Coupling of finite element and cellular automata method is described by the help of interpolation equations. It is found that proper use of transport equations and calculations of random probability number play pivotal role in describing microscale solute diffusion and solid fraction evolution in solidifying alloy. Use of moderate size simulation grid (cartesian) to counter mesh anisotropy along with selection of decentered square algorithm also helps in model refinement and optimization.
Rafique, M.M.A. (2018) Probabilistic Modelling of Microstructural Evolution in Zr Based Bulk Metallic Glass Matrix Composites during Solidification in Additive Manufacturing. Engineering, 10, 130-141. https://doi.org/10.4236/eng.2018.104010
 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.
 Kaufman, L., Perepezko, J. and Hildal, K. (2007) Synthesis and Performance of Fe-Based Amorphous Alloys for Nuclear Waste Repository Applications. Lawrence Livermore National Laboratory (LLNL), Livermore.
 Perim, E., et al. (2016) Spectral Descriptors for Bulk Metallic Glasses Based on the Thermodynamics of Competing Crystalline Phases. Nature Communications, 7, 12315.
 Rappaz, M. and Gandin, C.A. (1993) Probabilistic Modelling of Microstructure Formation in Solidification Processes. Acta Metallurgica et Materialia, 41, 345-360.
 Zhang, J., et al. (2013) Probabilistic Simulation of Solidification Microstructure Evolution during Laser-Based Metal Deposition. Proceedings of 2013 Annual International Solid Freeform Fabrication Symposium—An Additive Manufacturing Conference, Austin, August 2013, 739-748.
 Zhou, X., et al. (2016) Simulation of Microstructure Evolution during Hybrid Deposition and Micro-Rolling Process. Journal of Materials Science, 51, 6735-6749.
 Zinoviev, A., et al. (2016) Evolution of Grain Structure during Laser Additive Manufacturing. Simulation by a Cellular Automata Method. Materials & Design, 106, 321-329.
 Wang, Z.-J., et al. (2014) Simulation of Microstructure during Laser Rapid Forming Solidification Based on Cellular Automaton. Mathematical Problems in Engineering, 2014, 9.
 Tian, F., Li, Z. and Song, J. (2016) Solidification of Laser Deposition Shaping for TC4 Alloy Based on Cellular Automation. Journal of Alloys and Compounds, 676, 542-550.
 Markl, M. and Korner, C. (2016) Multiscale Modeling of Powder Bed-Based Additive Manufacturing. Annual Review of Materials Research, 46, 93-123.
 Fan, Z., et al. (2007) Numerical Simulation of the Evolution of Solidification Microstructure in Laser Deposition. Proceedings of the 18th Annual Solid Freeform Fabrication Symposium, Austin, August 2007, 256-265.
 Zhu, M.F., Lee, S.Y. and Hong, C.P. (2004) Modified Cellular Automaton Model for the Prediction of Dendritic Growth with Melt Convection. Physical Review E, 69, Article ID: 061610.
 Li, H.Q., Yan, J.H. and Wu, H.J. (2009) Modelling and Simulation of Bulk Metallic Glass Production Process with Suction Casting. Materials Science and Technology, 25, 425-431.
 Browne, D.J., Kovacs, Z. and Mirihanage, W.U. (2009) Comparison of Nucleation and Growth Mechanisms in Alloy Solidification to Those in Metallic Glass Crystallisation—Relevance to Modeling. Transactions of the Indian Institute of Metals, 62, 409-412.
 Wang, C.Y. and Beckermann, C. (1993) A Multiphase Solute Diffusion Model for Dendritic Alloy Solidification. Metallurgical and Materials Transactions A, 24, 2787-2802.
 Wang, W., Lee, P.D. and McLean, M. (2003) A Model of Solidification Microstructures in Nickel-Based Superalloys: Predicting Primary Dendrite Spacing Selection. Acta Materialia, 51, 2971-2987.
 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.
 Nastac, L. (1999) Numerical Modeling of Solidification Morphologies and Segregation Patterns in Cast Dendritic Alloys. Acta Materialia, 47, 4253-4262.
 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.
 Dezfoli, A.R.A., et al. (2017) Determination and Controlling of Grain Structure of Metals after Laser Incidence: Theoretical Approach. Scientific Reports, 7, Article ID: 41527.