Mesoscopic Constitutive Model for Predicting Failure of Bulk Metallic Glass Composites Based on the Free-Volume Model

A meso-mechanical damage model is developed to predict the tensile damage behaviors of bulk metallic glass composites (BMGCs) toughened by ductile particles. In this model, the deformation behaviors of the BMG matrix and particles are described by the free volume model and Ludwik flow equation, respectively. Weng’s dual-phase method is used to establish the relationship between the constituents and the composite system. The strain-based Weibull probability distribution function and percolation theory are adopted in characterizing the evolution of shear bands leading to the progressive failure of BMGCs. Moreover, the present model is performed under strain-controlled loading. Comparing to experiments on various BMGCs, the predictions are in good agreement with the measured results, which confirms that the present model successfully depicts the composite properties, such as yield strength, uniform deformation and strain softening elongation.


Introduction
To improve the poor damage tolerance of pure bulk metallic glass (BMG), many kinds of composite (BMGC) systems have been prepared, and many important conclusions were reached. However, an in-depth understanding of the inherent synergic effect among different constituents in BMGCs is still lacking. In comparison to simulations and experiments, theoretical models are more efficient and convenient in explaining their micro-deformations and composite effect. Moreover, quantitative relations are more efficient in optimizing the ductility/toughness of BMGs via rapidly tuning their microstructures. It is imperative to understand the correlations among processing, microstructures and properties for such composites. According to the thermodynamics and free energy principle, Marandi et al. [1] advanced an elastic-viscoplastic constitutive model for describing the finite deformation behaviors of BMGCs. They [2] further extended their model to better predict the stress-strain relations of in-situ BMGCs. Qiao et al. [3] were the first to consider the work-hardening ability of dendrite phase and softening of metallic glass matrix from the point of micromechanics view, and the predictions are in a fairly good agreement with tensile experiments. It should be noted that the interaction among the constituents is not fully reflected. Yang et al. [4] also established an analytical model to describe the deformation kinematics, free volume evolution, hardening, softening and viscosity of BMGs. Recently, Sun et al. [5] improved Qiao's previous micromechanics model to predict the tensile behaviors of in-situ BMGCs more accurately based on the measured data via nano-indentation. Rao et al. [6] proposed a new meso-mechanical constitutive model to predict the monotonic tensile/compressive deformation of BMGCs with toughening phases, but their analytical model has very complicated formulas. Jiang et al. [7][8][9][10] regarded the shear bands as micro-cracks, Materials 2018, 11, 327 2 of 11 established their equivalence relation, and finally developed two micromechanics models based on the incremental tangent stiffness and secant modulus, respectively.
These analytical models can successfully reflect some main features, such as yield strength, strain hardening and stress softening elongation, of ductile particles filled BMGs. However, they cannot fully take account of the inherent microstructure evolution and deformation features of BMG matrix. It is expected that shear bands will gradually transform into micro-cracks with the deformation increasing, and correspondingly the stress-strain curve of BMGCs will enter the stress-softening stage. To the author's knowledge, the damage effect in the BMGCs was not addressed in these analytical models, and a simple micromechanics model is always required to describe their intriguing mechanical response.
This paper aims to build an analytical damage model for predicting the tensile failure of BMGCs toughened by ductile particles. The deformation behaviors of BMG matrix and particles are described by the free volume model and Ludwik flow equation, respectively, and Weng's homogenization frame is adopted to establish the interaction between the constituents and composites. As compared to the other models, the present model is more convenient to apply, and more readily to be expanded. The developed model is performed under strain-controlled loading, and verified by modeling the monotonic stress-strain relations of particle toughened BMGCs.

Analytical Model of BMGCs
The BMGCs are filled with ductile particles, and the stress-strain relations of the constituents should be described with the proper constitutive equations. For such dual-phase composites, where both phases can undergo plastic flow, Weng [11] developed an analytical model to depict the stress-strain relations of the composites, and later extended by Zhu [12]. Their formula will be used as the basis of a new micromechanics method for BMGC, and a perfect interfacial bonding between two phases is assumed. For a dual-phase composite, the particle phase is referred as phase 1 and the BMG matrix as phase 0, and those of the composite are expressed by symbols without any script. All the tensors and vectors are written in boldface letters. The volume fractions for the particle and matrix phases are denoted by c 1 and c 0 , respectively, and should satisfy the condition c 1 + c 0 = 1.

Constitutive Model of BMGs
The shear band evolution controls the fundamental deformation mechanisms in BMGs. At the microscopic level, shear band formation is accompanying with the evolution of the local structural order. One atomistic mechanism capturing shear band formation and evolution in BMGs is the free volume model proposed by Spaepen [13] and further extended by Steif [14]. From the continuum mechanics point of view, the shear band is regarded as a consequence of strain softening and acts as a strain-localization phenomenon. This model considers free volume as an internal state variable to characterize the structural evolution of BMGs at the atomic level.
Following a J 2 -type, small strain visco-plasticity framework, the free volume model is adapted into the multi-axial stress state. The total strain rate in the BMG matrix is written as which includes the elastic part, σ kk δ ij , and the plastic part, . ε p ij . For the BMG matrix, the plastic strain rate, i.e., the flow equation is expressed as where σ ij = σ ij − σ kk δ ij /3 is the deviatoric stress tensor and σ eq = (σ ij σ ij ) 1/2 is the von Mises stress. f is the flow stress, which is defined by where f 0 is the frequency of atomic vibration; ∆G m is the activation energy; k B is the Boltzmann constant; T is the absolute temperature; Ω is the atomic volume; and ξ is the concentration of free volume. The free volume evolution equation under multi-axial stress state is written as where α 0 is a geometrical factor of order unity; ν * is a critical volume; S is the Eshelby modulus, given by ; v is Poisson's ratio; µ is the shear modulus; and n D is the number of atomic jumps needed to annihilate a free volume equal to ν * and is usually taken to be 3-10.

Constitutive Model of Ductile Phases
The Ludwik equation is adopted for ductile particles in terms of von Mises effective stress and plastic strain as where ε p eq = (2ε p ij ε p ij /3) 1/2 ; σ y , h and n are the initial yield stress, strength coefficient and the work-hardening exponent, respectively; and these material parameters will be determined by fitting with a measured stress-strain curve. Moreover, Hencky's flow rule is adopted,

Homogenization Method for BMGCs
For dual-phase composites, Weng's model is used to establish the relationship among ductile particles, matrix and the resulting composites under monotonic uniaxial tension. The detailed derivations are found in their original work [11]. The relationship between the hydrostatic and deviatoric strains of BMGCs are defined by where α s 0 and β s 0 are the components of the classical Eshelby's tensor for spherical inclusions, and given as and κ and µ denote the bulk and shear moduli, and are written as follows to satisfy the isotropic relations, where E and ν are the Young's modulus and Poisson's ratio, respectively; and E and ν with superscript "s" denote the secant modulus and secant Poisson's ratio, respectively, defined by The relationship between the hydrostatic and deviatoric strains of the constituents and those of BMGC are given as

Failure of the BMG Matrix
During the deformation of BMGCs, shear bands gradually transform into micro-cracks with the applied loading, and therefore are simplified as micro-cracks. Moreover, shear bands lead to the stress softening behavior, which is similar to the effect of micro-cracks on the mechanical behaviors [15]. From this viewpoint, shear bands are equivalent to micro-cracks, and then some analytical models for micro-cracks are also applied to the strain localization effect induced by shear bands [16].
The representative volume element (RVE) is often utilized to account for the micro-crack orientation statistics in the inhomogeneous materials. The micro-cracks generated by shear bands are supposed to be random, and thus the corresponding effective moduli are given by [17] where the subscript "in" denotes the intact materials with no micro-cracks. For the ductile phases, the failure criterion based on statistical probability is associated with strain levels. The strain-based Weibull distribution function is introduced to characterize the shear band induced fracture as where ε p is the plastic strain and ε 0 is the reference strain, and since there is no data available for parameter m, which is need to be determined by fitting from a final stage with damage so that the predicted stress-strain relations can duplicate the experiments. Then, the density of shear-bands in the BMGCs is defined by where ρ 0 denotes the saturate density of shear-bands. After introducing the percolation threshold of shear-band propagation in the BMG matrix, the shear-band density is given by where c cr and χ are constants, which will be assigned by the experimental data. In fact, almost the precipitates are randomly distributed in real BMGCs. Therefore, some necessary statistics experiments should be performed to examine the flaw sensitivity and reliability of BMGs.

Numerical Implementation
The developed model is performed under strain-controlled loading, and the detailed algorithm is explained here. For a time interval from t n to t n+1 (∆t n+1 = t n+1 − t n ), the necessary variables at time t n , such as σ n , ε n , ε n and σ (1) n , are known, and a uniform strain increment ∆ε n+1 is given. The average strain increments ∆ε (r) n+1 (r = 0, 1) in each phase are determined by Equations (12)-(15), and then the secant modulus, secant Poisson s ratio and plastic strain can be computed by the constitutive models for each phase. The overall stress increment corresponding to the current strain increment can be solved.
The key issue in the modeling procedure is to fix ∆ε (r) n+1 (r = 0, 1). At first, the initial tentative value of ∆ε (1) n+1 is given by ∆ε Then, Then, the compatibility of a strain increment ∆ε (1) n+1 in particles is checked by the residual R as follows R represents the difference between particle s tentative average strain increment and that obtained by Weng s model. If R < TOL (R→0), the iteration stops. Otherwise, ∆ε n+1 is updated with another new iteration, ∆ε In the uni-axial strain-controlled loading, only the strain increment in loading direction (∆ε 11 ) is exactly given, and the others should be determined by the overall-stress constraints. An additional iteration procedure should be performed to obtain the values of ∆ε 22 (∆ε 33 ) as For a given strain increment (∆ε 11 , ∆ε 22 = ∆ε 33 ), determining the exact value of ∆ε 22 should be a key procedure in the computation. A simple method is explained here. The iteration of updating ∆ε 22 is also performed from an initial value of zero, and a new ∆ε 22 = ∆ε 22 + ∆ is assigned with a very small ∆. The value of ∆ is adjusted correspondingly based on the computation precision. The iteration stops if Equation (26) is satisfied; otherwise, ∆ε 22 is updated according to the above equation for the next iteration.
Based on the completed stress-strain curves by the above iteration, the shear-band density ρ is determined, and then the overall elasticity is given. Based on Equations (16)- (18), the stress-strain curves for the BMGCs can be re-evaluated by involving the damage effect. The above-mentioned numerical implementation procedure is illustrated by the flowchart in Figure 1, and a Fortran code was programmed to predict the stress-strain relations of BMGCs. Based on the completed stress-strain curves by the above iteration, the shear-band density ρ is determined, and then the overall elasticity is given. Based on Equations (16)-(18), the stress-strain curves for the BMGCs can be re-evaluated by involving the damage effect. The above-mentioned numerical implementation procedure is illustrated by the flowchart in Figure 1, and a Fortran code was programmed to predict the stress-strain relations of BMGCs.

Comparisons with the Experiments
Szuecs et al. [18] prepared a Zr-based BMGCs with dendrite volume fraction of c 1 = 20%, and measured their mechanical properties under tension. The dendrites are described by Equation (5)

Comparisons with the Experiments
Szuecs et al. [18] prepared a Zr-based BMGCs with dendrite volume fraction of c1 = 20%, and measured their mechanical properties under tension. The dendrites are described by Equation (5)       Comparisons of the macroscopic stress-strain relations between the prediction and experiments for BMGCs with particle volume fraction of 43% [19]. ε 0 = 0.1 and m = 16 are adopted in the computation. BMGC's deformation increases step by step; a slight work hardening was observed, followed by remarkable improvement in plastic strain level. By including the damage effect, the collapse stage during deformation could be clearly presented. Finally, the predicted results are in good agreement with the measured data. On the other hand, the increase in dendrite loading level will greatly impair the yielding stress. The present comparison confirmed that this model can reflect the dependence of composite strength and ductility on phase volume fractions.   BMGC's deformation increases step by step; a slight work hardening was observed, followed by remarkable improvement in plastic strain level. By including the damage effect, the collapse stage during deformation could be clearly presented. Finally, the predicted results are in good agreement with the measured data. On the other hand, the increase in dendrite loading level will greatly impair the yielding stress. The present comparison confirmed that this model can reflect the dependence of composite strength and ductility on phase volume fractions.
The impacts of model parameters in the present method should be carefully explained. Figures 5  and 6 demonstrate the effect of Weibull modulus m and reference strain ε 0 on the overall stress-strain curves, respectively. In Figure 5, the Weibull modulus m ranges from 2 to 16, and the reference strain ε 0 = 0.06. At a given particle volume fraction, the plastic elongation increases with increasing Weibull modulus, which is controlled by the physics meaning of Weibull modulus. Additionally, the dependence of the overall stress-strain relations on reference strain ε 0 is illustrated in Figure 6. As expected, the uniform stretches decrease with the increase in micro-crack density. For the reference strain less than 0.1, the difference in plastic elongation becomes evident. Strain ε xx (%) Figure 4. Comparisons of the macroscopic stress-strain relation between the prediction and experiments for BMGCs with different particle concentrations [20]. ε0 = 0.11 and m = 8 are adopted in the computation. The impacts of model parameters in the present method should be carefully explained. Figures 5 and 6 demonstrate the effect of Weibull modulus m and reference strain ε0 on the overall stress-strain curves, respectively. In Figure 5, the Weibull modulus m ranges from 2 to 16, and the reference strain ε0 = 0.06. At a given particle volume fraction, the plastic elongation increases with increasing Weibull modulus, which is controlled by the physics meaning of Weibull modulus. Additionally, the dependence of the overall stress-strain relations on reference strain ε0 is illustrated in Figure 6. As expected, the uniform stretches decrease with the increase in micro-crack density. For the reference strain less than 0.1, the difference in plastic elongation becomes evident. 100 Pa) Based on the work hardening coefficient, the ability to endure damage can be judged for engineering materials. The reliance of work-hardening rate dσ/dε on particle concentration is shown in Figure 7. Here, the definition of dσ/dε, is used to describe the resistance to flow localization, i.e., necking. It is noted that the value of dσ/dε is higher for BMGCs with higher filler volume fraction. Moreover, the value reduces rapidly with deformation for dilute composites. These curves also show that uniform elongation becomes larger with the higher initial hardening rate. Since the material with a high dσ/dε can lead to a much uniform plastic flow, which is against the early appearance of deformation localization, and an increased stretch is finally reached.
Based on the work hardening coefficient, the ability to endure damage can be judged for engineering materials. The reliance of work-hardening rate dσ/dε on particle concentration is shown in Figure 7. Here, the definition of dσ/dε, is used to describe the resistance to flow localization, i.e., necking. It is noted that the value of dσ/dε is higher for BMGCs with higher filler volume fraction. Moreover, the value reduces rapidly with deformation for dilute composites. These curves also show that uniform elongation becomes larger with the higher initial hardening rate. Since the material with a high dσ/dε can lead to a much uniform plastic flow, which is against the early appearance of deformation localization, and an increased stretch is finally reached.

Discussion
The comparisons between the predictions and experiments indicate that the present model can describe the damage behaviors of BMGCs accurately. Under the multi-axial stress state, the BMG matrix in the composites exhibits a certain degree of plasticity, which is absent for pure BMG under uniaxial loading. Based on the free volume model and iso-hardening flow law for the constituents, the equivalent plastic behaviors of BMG matrix can be determined by the present model.
Based on the Eshelby's tensor and Mori-Tanaka mean field frame, many micromechanic models have been developed, and are commonly applied under stress-controlled loading. Such models have difficulty describing strain-softening deformation. In the present work, the secant modulus and strain-controlled formula were adopted, and the aforementioned difficulty can be easily avoided. In Rao's mesoscale model, two important terms need to be performed: The first one is the algorithmic tangent operator, obtained by consistent linearization of the time discretized constitutive equations. The second is a new one and called an affine strain increment. Therefore, their model has very complicated in mathematical formulas. On the other hand, the present model does not need to calculate the tangent stiffness instead of the secant modulus, and is very simple in mathematic form, which is readily realized in programming. It is expected that particle volume fraction is usually higher than 30%, even exceeding 50% for many BMGCs. The interaction between particles and matrix cannot be well considered by the Mori-Tanaka method. Additionally, shear bands will finally transform into micro-cracks by increasing the applied deformation, and moreover the microstructure evolution is a complicated process, and their effect on the macroscopic performance is not clear yet. These problems should be deeply studied and tackled by improving the present model in the future.

Conclusions
Based on Weng's method for dual-phase composites, a meso-mechanical damage model was proposed to depict the tensile failure of BMGCs. Free volume theory was adopted to describe the strain softening of BMG matrix. The strain-based Weibull probability distribution function and percolation theory are used to reflect the damage effect induced by the transformation from shear bands to micro-cracks. The displacement controlled loading was applied to predict the collapse stage during the deformation. The final comparisons with the experiments confirm that the present model can replicate the monotonic tensile stress-strain relations of BMGCs until final failure.