Crack Growth Prediction Using Cohesive Zone Model (CZM)
Although MD has proven to be an effective tool in explaining the atomic-scale failure analysis, but it cannot capture the crack propagation dynamics and the constitutive behavior (traction-displacement relationship) by itself. In this regard, CZM, along with MD, has been extensively employed, attributed to its capability of avoiding stress singularity at the crack tip and the ability to represent the physics of fracture at the atomic scale. The crack tip stress and the crack propagation

FIGURE 7.16 Schematic of FCC SC Ni with edge crack considered for CZM-based crack propagation (where a symbolizes the lattice constant of Ni).
are highly dependent on the crystal orientation, attributed to the difference in the fracture mechanism that transpires during the crack growth process, thereby affecting the crack opening and the local stress concentration. The specimen with three different orientations, orientation 1 (X [100]. Y[010], Z[001]), orientation 2 (X[l 10], Y[IlO], Z[001]), orientation 3 (X[l 11], Y[T 10], Z[Il2]), is stretched along the Y-direction, as shown in the schematic (refer to Figure 7.16) [18]. In case of orientation 1, initially, the propagation of crack is delayed; however, in the later time period, ductile-to-brittle transition of the crack propagation behavior transpires, attributed to the nucleation and growth of void, leading to a change in the site of stress concentration. As a result, a rapid propagation of the new crack is witnessed, leading to a drastic failure of the material. However, in orientation 2, the crack propagation occurs easily without any resistance from the material in a small span of time, suggesting brittle fracture behavior. However, in case of orientation 3, the crack propagation is extremely slow' on account of formation of slip bands at the crack tip, leading to crack tip blunting. Moreover, the crack growth rate and the crack opening displacement vary significantly owing to the difference in the orientation of the SC Ni [18]. However, CZM has been extensively employed to study the interfacial fracture behavior of materials.
The CZM has been extensively employed to study the interfacial fracture behavior of materials. The atomistic approach to CZM (using MD) is largely employed to derive traction separation law characterizing the properties at the interface and also to explain the fracture mechanism associated with the interface, as showm in Figure 7.17 [2]. Pristine materials seem to exhibit excellent mechanical properties. Mutation in the strength behavior is apparent in the presence of a defect; for instance, grain boundary plays a significant role in the initiation and propagation of a crack. The grain boundary acts as a transition point where the direction of crack propagation is altered. The crack grows along the closed-pack plane; as the orientation variation between two grains, the direction of crack changes at the grain boundary when crack crossovers from

FIGURE 7.17 Schematic representation of alteration in propagation of crack at grain boundary: (a) Single crystal, (b) bicrystal, and (c) tricrystal. (From Zhang, Y. et al., Results in Phys., 7, 1722-1733, 2017. With permission.)
one grain to another. In a nutshell, grain boundary plays a vital role in altering the direction of the propagation of crack.
At a triple junction, an orientation difference exists between three grains; therefore, when subjected to tensile stress, the triple junction not only nucleates dislocation but also emits dislocations from every grain boundary. Dislocations with different orientations interact with each other to result in the formation of sessile locks such as the Lomer-Cottrell lock, which impedes any further movement of dislocation, leading to dislocation pile-up, followed by the initiation and propagation of crack. The free movement of dislocation is also impeded at grain boundaries, leading to the concentration of stresses. Such a high concentration of stresses and intersection of the Lomer-Cottrell locks at the triple junction leads to the formation of void at the triple junction, leading to nucleation. The grain boundary plays a vital role in aiding the crack propagation by acting as an emission resource of dislocation [2].
The interfacial strength of metal-ceramic composite is greatly dependent on both metal and ceramic phases; therefore, the crack propagation in the interface exhibits a different behavior from that of either parent phase. The metal-ceramic interfacial fracture can be subdivided into four types, as shown in Figure 7.18 [7]. Type I

crack propagation comprises the brittle behavior of crack propagation with smooth crack surface (refer to Figure 7.18a). Such kind of crack propagation behavior is observed when the atomic configuration at the interface is regular, and as a result, the effect of lattice mismatch at the interface is trivial. As a result, a butterflyshaped stress distribution along the interface is observed. Type II interfacial crack propagation constitutes nucleation of a twin at crack tip, leading to crack blunting, thereby impeding any further propagation of crack. The Shockley partial emitted at the crack tip traverses to a stable position and leads to the formation of twin boundaries of several atomic widths, thereby blunting the crack and impeding any further motion (refer to Figure 7.18b). Type III interfacial crack occurs when a daughter crack appears in front of the parent crack and coalesces to propagate in a brittle manner. Such a scenario arises due to strong lattice mismatch at the interface, leading to the generation of a daughter crack ahead of the parent crack, eventually leading both the cracks to coalesce and thereby extending in a brittle manner with rough crack surface, as shown in Figure 7.18c. Type IV crack propagation constitutes the formation of stacking fault due to emission of full edge dislocation, thereby blunting the crack tip and impeding any further propagation (refer to Figure 7.18d). The adhesive strength coefficient (ASC) plays a vital role in determining the type of crack propagation that will take place at the ceramic-metal interface. For Al-SiC interface with an ASC parameter with a value greater than 0.9, the crack propagates with either Type I or Type II with clean, smooth crack surfaces. However, Type III and Type IV are dominant with rough crack surface when the ASC parameter is below 0.9 for Al-SiC interfaces.
The MD-based CZMs are extensively employed to parameterize the cohesive model. The relationship between the cohesive traction and separation displacement is mathematically presented as follows:
where cr, is the maximum cohesive strength; / is a dimensionless function that relates to the shape of cohesive traction-separation displacement curve; A, B, and c are constants, which are calculated based on the simulations; T denotes the temperature in K. A is non-dimensional parameter that establishes a relationship between normal (и„) and shear (m5) separations to maximum normal (<5„) and shear (<5() separations, given in Eq. 7.14
Failure occurs in a material when the value of A equals 1. When Al-SiC specimen with crack at its interface is subjected to Mode I loading, it fractures in a brittle manner via large-scale debonding across the interface. However, when subjected to Mode II loading, shear bands start appearing, leading to blunting of the crack tip, thereby leading to crack healing and sliding of the SiC crystal. The traction-response

FIGURE 7.19 Traction-separation relationship for (a) Mode I failure and (b) Mode II failure in Al-SiC. (From Dandekar, C.R., and Shin, Y.C.. Compos. Part A: Appl. Sci. Manuf, 42. 355-363, 2011. With permission.)
behaviors of Al-SiC interface for Mode I and Mode II loading have been presented in Figure 7.19. The traction-separation law for Al-SiC interface based on the MD-based simulation [29] is established as follows:

This established traction separation relationship for Al-SiC interface can be further employed to predict the tensile behavior of the composite.
Ti6A14V metal comprises two phases, that is, a-TiAl and /З-TiV. The crack propagation behavior of Ti6A14V-TiC interface consists of a-TiAl-TiC interface behavior and /З-TiV-TiC interface behavior. Generally, Ti6A14V deforms along the prismatic planes, but the introduction of crack at the Ti6A14V-TiC interface increases the critical resolved shear stress (CRSS) of the prismatic planes, attributed to the generation of stacking faults along these planes, as crack boundaries can act as free surfaces [37]. However, since no stacking fault is present along the basal planes, the crack propagates freely along the basal plane above the interface. In case of the /3-TiV-TiC interface, the crack propagates along the interface. This is due to the reason that the cohesive energy of the interactions between the atoms in TiC and along the interface is higher than that in /З-TiV. Moreover, stacking fault formation transpires along the close-packed direction, thereby increasing the CRSS of /З-TiV and leading to crack propagation along the interface [37]. However, in case of Mode II loading, ductile fracture behavior is found to be dominant. Incoherent interfaces are formed due to lattice mismatch at the interface, which act as major sources of edge dislocation. Under such a scenario, slip phenomenon is observed along the interface under shear loading. The crack-opening displacement for Mode I and Mode II loading is provided in Figure 7.20 [37]. The traction-separation law established for Ti6A14V-TiV interface for Mode I (F„(A)) and Mode II (F,(A)) [37] loading is given as follows:


FIGURE 7.20 Traction-separation relationship for (a) Mode I failure and (b) Mode II failure in Ti6AI4V-TiC. (From Elkhateeb, M.G.. and Shin, Y.C.. Mater. Des.. 155, 161-169. 2018. With permission.)
where F„ and F, correspond to Mode I and Mode II loading, respectively, and
The traction-separation relationship established for A120,-A1 interfaces from MD-based CZM is provided here below [38]
The presence of porosity in metals (such as Al) reduces the ductility of the material, whereas in case of ceramics (such as Al203), it imparts ductility to the material [38]. The failure stress reduces irrespective of where the pores are present, that is, only in metal or only in ceramic or in both the phases. However, the failure strain is highly dependent on the position of porosity. The presence of porosity in the metallic phase reduces the failure strain. On the other hand, the presence of porosity in the ceramic phases increases the failure strain. Nevertheless, the effect of porosity on failure strain is nullified when the porosity is present in both the phases.
In case of bimetallic interfaces such as Al-Si bimetallic specimen, crack propagation transpires in a brittle manner irrespective of the crystal orientation. However, the smoothness of the fracture surface is subjected to the orientations of the bimetallic specimens. Clean and flat crack surfaces are obtained for Al(100)/ Si(100) and Al( 111)/Si(l 11) specimens, and rough surface are obtained for Al(110)/ Si(110) specimens.
Application of Crack Tip Opening Displacement in Predicting Fracture Behavior
The CTOD test of a material measures the resistance of the material to the propagation of the crack [39]. Materials with some plasticity prior to failure can be subjected to CTOD. A schematic representing the definition of CTOD is presented in Figure 7.21. Figure 7.21a corresponds to the opening displacement of the original crack tip, as indicated by 8. On the other hand, 8 in Figure 7.21b represents the displacement at the intersection of the vertex (with an angle of 90°) with either surface of the crack.
An accurate prediction of CTOD involves its division into two parts: elastic and plastic, as provided in Eqs. 7.20 and 7.21 and exemplified in Figure 7.22. In the provided figure, a is length and b is the width of the rest of the specimen for which CTOD is calculated, p is the dimensionless rotational factor, which is used to locate the center of the hinge. Using all these factors, CTOD can be calculated as follows [40]:



FIGURE 7.22 Calculation of CTOD.
where К symbolizes the stress intensity factor, u denotes Poisson’s ratio, E indicates Young’s modulus, cr„ is the yield strength of the specimen, and Vp symbolizes the plastic component of CTOD.
The J-integral has proven to be an effective tool in defining the elastic-plastic fracture mechanics of materials. The J-integral is a vector with three components, with each component belonging to each dimension of the Cartesian coordinate system. The value of J-integral can be calculated by integrating along a trajectory, for instance, Г, as shown in Figure 7.23. The specific elastic energy is computed from the corresponding stress values and strain values along the trajectory. Moreover, the estimation of the J-integral is also dependent on the unit normal vector n, the stress vector F, and the spatial derivatives of the displacement vector. However, at the atomic scale, the J-integral can be expressed in terms of energy, as follows [41]:


where U is the potential energy of the system.
The relationship between CTOD and J-integral is defined as follows [33]:

where in is a constrain factor, oy is effective yield strength, // is the plastic eta factor for J-integral calculation, Ap, stands for plastic area under the load-CTOD curve, and В symbolizes the thickness of the specimen.
The fracture behavior of SC Si is characterized at the atomic scale by employing the CTOD approach, thereby identifying the fracture mechanism responsible for the crack propagation behavior in Si [42]. A transition in the fracture behavior from brittle to ductile is observed as the temperature increases from 500 to 1200 K. The crack propagation eventuates via the crack plane (100) at 500 К by breaking of the covalent bonds, followed by the amorphization of the specimen and ring formation. On the other hand, at 1200 K, the crack propagation constituted of two minor brittle jumps, followed by the emission of dislocation along the slip plane [111], leading to crack tip blunting and thus restricting any further propagation of crack. CTOD90degree gives a better estimation of the fracture toughness of Si, as opposed to CTOD blunting [42]. Although CTOD blunting overestimates CTOD, it provides a detailed description of the hinge rotation mechanism. The overestimation of CTOD in case of CTOD blunting is associated with the bending moment that is introduced during the deformation. The stress intensity factor is calculated for mode I loading from the measured value of CTOD via the following equation [43]:
where ry is the radius of the plastic zone at the crack tip. The stress intensity factor increases with temperature, based on the calculations of CTOD. The prediction of the stress intensity factor based on CTOD at 500 К is found to be somewhat higher than the continuum-mechanics-based prediction of the stress intensity factor. The increment in the stress intensity factor with increasing temperature is attributed to the increasing amorphization of Si, which resists failure. The high stress intensity factor at low temperatures is also influenced by lattice trapping mechanism [42], which is highly dependent on temperature. The dependency of lattice trapping mechanism is attributed to the fact that thermal fluctuations are vital to pull out atoms from the trapped state. However, CTOD90degree estimation for stress intensity factor is not accurate in case of bcc iron. The CTOD90degree method provides lower stress intensity value than the continuum mechanics, based on the prediction of the stress intensity factor. In order to correctly estimate the stress intensity factor via CTOD, the deformation mechanism of the crack propagation behavior is vital [44].
The procedure of calculating CTOD blunting and CTOD90degree is provided in this section. First, a crack profile is to be obtained from the calculation of position of surface atoms in x-axis versus у-axis graph. A python script is developed to obtain surface atoms by iterating over each atom in the system and is denoted as atoms with bond order less than four, as surface atoms and the coordinates of such atoms are stored in the file, which was then used for obtaining the crack profile plot. The CTOD90degree is obtained via the intersection of lines extending from the tip of the crack at an angle of 45° (as shown in Figure 7.22b). The distance between the two intersects of these curves is calculated as follows:
The CTOD blunting method determines the CTOD based on the transition of linear- to-polynomial behavior of the crack profile via the measurement of the angles between the linear approximations associated with the surface atoms of the side fronts in the crack profile graph and the gradient of the crack profile function showing polynomial behavior.
The calculation of the J-integral at the atomic scale in an energy form has been already exemplified in Eq. 7.22. According to the equation, the J-integral is treated as the difference in the potential energy of a specimen loaded identically with crack lengths of a and a + 8a. The term dU indicates the difference in the potential energy of the specimen with crack lengths a and a + da. In other words, the calculation of the J-integral at the atomic scale consists of the calculation of the potential energy difference between the two cracked models at a critical tension state, which can be easily obtained through MD route via potential functions of the systems [42]. Tw'o identical strained models are taken into consideration to compute the potential energy difference of crystals. The two models with different initial crack lengths are initially in relaxed conditions (i.e., state of no load). Thereafter, the models are subjected to load in Mode I loading. The potential energy difference is computed corresponding to a model when the crack just starts to propagate (also known as the critical point of the model) with the other model strained to the same value as the previous model.
The J-integral can be defined as the boundary integral of the Eshelby stress tensor 5, as follows [45]:

where W is internal energy density, / is the identity tensor, H denotes the displacement gradient, and P symbolizes the Piola-Kirchoff stress. In Eq. 7.28, the Eshelby stress acts on N, the outward normal to the surface 8Q enclosing the region Q in the initial configuration. If a system changes its configuration from Q| to Cl2, then the J-integral is given by the difference in the tw'o configurations, as discussed earlier.

In order to define the J-integral at the atomic scale W, H, and P have to be derived from the atomic data. In order to define these parameters, localization function or kernels (y (*)) are used, whose basic properties are being positive and normalized.
W can be aefined as a local weighted average of the potential energy per atom. Mathematically, it is represented as

where 0* corresponds to the potential energy of the reference configuration, and W(X) is the potential energy density corresponding to the zero temperature. H (H = Vxm) can be defined in terms of displacement (и) of the material from its reference or initial position (denoted by X).
The definition of H requires the definition of и (displacement). This is done in mass-weighted fashion, given as follows:
Thus, the displacement gradient can be computed as follows:
The computation of P is influenced by the force between the two atoms and the difference between their positions in the initial or reference configuration.
where XaP = Xa- Xp, BaP indicates the bond function, which is given by
The calculation of the J-integral based on this method provides an almost-accurate estimation of the J-values for both bcc and fee configurations.