Probabilistic Modelling of Microstructural Evolution in Zr Based Bulk Metallic Glass Matrix Composites during Solidification in Additive Manufacturing

Bulk metallic glass and their composites (BMGMCs) are a new class of materials which possess superior mechanical properties as compared to existing conventional materials. Owing to this, they are potential candidates for tomorrow’s structural applications. However, they suffer from poor ductility and little or no toughness which render them brittle and they manifest catastrophic failure under applied force. Their behavior is dubious, unpredictable and requires extensive experimentation to arrive at conclusive results. In present study, an effort has been made to design bulk metallic glass matrix composites by the use of modeling and simulation. A probabilistic cellular automaton (CA) model is developed and described in present study by author which is used in conjunction with earlier developed deterministic model to predict microstructural evolution in Zr based BMGMCs in additive manufacturing liquid melt pool. It is elaborately described with an aim to arrive at quantitative relations which describe process and steps of operations. Results indicate that effect of incorporating all mass transfer and diffusion coefficients under transient conditions and precise determination of probability number play a vital role in refining the model and bringing it closer to a level that it could be compared to actual values. It is shown that proposed tailoring can account for microstructural evolution in metallic glasses.


Introduction
Bulk metallic glass matrix composites [1] 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 [2]- [8].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 [2], a probabilistic model of nucleation and growth [9] [10] 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 [11], LAMP [12] on 316L SS [10], two dimensional CAFE [13], cladding [14] [15]), selective laser melting (SLM) using CAFE [14] [15] [16] [17], CAPF, CAFVM [18] and modified CAFE [19] 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 [20] [21] 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 [22] 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.

Model
In order to formulate model in present case, following assumptions are drawn.

Assumptions
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).
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 [9] 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).

Procedure of Cellular Automaton Model
Cellular automaton model consists of [24] 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 t ij a 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.

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 [11] [23]; ( ) . Schematic of modified 2D decentred square CA growth algorithm [23].
respectively, where C L is solute concentration in liquid, C S is same in solid and D L is solute diffusion coefficient in liquid and D S is same in solid.Similarly, solute diffusion in a mushy cell is given by below equations ( ) ) where L D and s D are the diffusion coefficients of the liquid phase and solid phase, respectively.s f is the solid fraction, E C is the equivalent concentra- tion, and E D 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 [25] [26].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 o s f = solid fraction at previous time step Time step itself may be calculated by where Δx is cell size, and V max is maximum growth velocity of cell during each time step.
Note: solid fraction (f s ) and solid fraction increment ( s f δ ) 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 s f is the solid fraction where Thomson coefficient, k = mean curvature of S/L interface.

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 [24]; Step 1: During this equilibrium at solid liquid interface is determined by 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 1without replacing any subscript or term.
Step 4: Partitioning of solute 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 C E , which may be defined as; ( ) ( ) when f s = 1, C E = C S and when f s = 0, C E = C L Thus, diffusion equation may be expressed as; ( ) where C E is equivalent concentration defined by above equation and D E is equivalent diffusion coefficient [23].This is schematically represented in Figure 3

Interpolation
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 [27].
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

3) Crystallography
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 (T m ), its state changes to liquid.If it does not pass T m , it stays as solid and if it stays at T m , 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 [10] is determined by equation This equation is tailored to take into account nucleation occurring at where s N = Number of new nucleation sites at the molten pool wall and as N = Total number of liquid cells at the molten pool wall.Similarly ( ) where v N = Number of new nucleation sites inside bulk liquid and av N = 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 ( j i W ) account for this effect.The growth length of cell i with regards to its liquid neighbor j (one of possible 26 cells) at time t c can be calculated by: ( ) W 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 j i L ).Orientation weight coefficient itself may be calculated by [28]; , , where X w , Y w and Z w may be calculated by solving below matrix where ( ) ( ) , , , , , p q r p q r x x x y y y and ( ) , , p q r z z z are direction cosines of [100] and [010] and [001] dendrite arms relating to the coordinate x, y, z axes respectively and ( ) , , j j j i i i p q r are direction cosines of vector j i L relating to the coordinate x, y, z axes respectively.
In detail, j i l 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 cell L l = if j is one of 6 nearest neighbors, if j is one of twelve second-nearest neighbors, and If cell j is captured by i, it will be assigned same orientation as cell i.Interested reader may find further detail in [28].

Conclusion
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.

Figure 1 .
Figure 1.Schematic of crystallographic orientation of one direction. below.

Figure 3 .
Figure 3. Schematic diagram of growing in a cell.

Figure 4 .
Figure 4. Representation of interpolation in one cell.
Bulk liquidNucleation 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