Mathematical modeling of the crack growth in linear elastic isotropic materials by conventional fracture mechanics approaches and by molecular dynamics method: crack propagation direction angle under mixed mode loading

The crack growth directional angles in the isotropic linear elastic plane with the central crack under mixed-mode loading conditions for the full range of the mixity parameter are found. Two fracture criteria of traditional linear fracture mechanics (maximum tangential stress and minimum strain energy density criteria) are used. Atomistic simulations of the central crack growth process in an infinite plane medium under mixed-mode loading using Large-scale Molecular Massively Parallel Simulator (LAMMPS), a classical molecular dynamics code, are performed. The inter-atomic potential used in this investigation is Embedded Atom Method (EAM) potential. The plane specimens with initial central crack were subjected to Mixed-Mode loadings. The simulation cell contains 400000 atoms. The crack propagation direction angles under different values of the mixity parameter in a wide range of values from pure tensile loading to pure shear loading in a wide diapason of temperatures (from 0.1 К to 800 К) are obtained and analyzed. It is shown that the crack propagation direction angles obtained by molecular dynamics method coincide with the crack propagation direction angles given by the multi-parameter fracture criteria based on the strain energy density and the multi-parameter description of the crack-tip fields.


Introduction
One of the fundamental ideas in the fracture assessment of brittle fracture is the local mode I concept [1]. The proposal of mode I dominance was suggested when dealing with cracked plates under plane loading and transverse shear, where the crack grows in the direction almost perpendicular to the maximum tangential stress (MTS) in radial direction from its tip [1]. This theory is one of the widely used theories for mixed mode crack growth [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. In more detail the criterion states that the crack propagation starts along the direction on which the tangential stress becomes maximum. The fracture occurs when the maximum tangential stress reaches a critical value for the material equal to the fracture stress in uniaxial tension. For the calculation of the tangential stress a critical distance from the crack tip must be introduced. But the question is what one can use as a critical distance from the crack tip? To overcome this difficulty the concept of a core region surrounding the crack tip has been proposed by Sih [2]. The idea is the continuum mechanics solution, as well as experimental measurements, stop at some distance from the crack tip. The distance serves as a scale size of analysis at the continuum level. Together with the MTS criterion the strain energy density (SED) has been used to formulate failure criteria for materials exhibiting both ductile and brittle behavior. Dealing with the strain energy concept it is worth to note that it is necessary to introduce the radius of "core region" surrounding the crack tip [1,2]. The key idea is that the continuum mechanics stops short at a distance from the crack tip. The strain energy density factor ( S ) was defined as a product of the strain energy density by a critical distance from the point of singularity. Therefore it is necessary to introduce the critical distance from the crack tip [3][4][5][6][7][8][9][10][11][12][13][14][15][16].
In [5,6] a discussion on the need to consider several more initial terms of the Williams series expansion during fracture behaviour assessment is presented. The paper [6] shows that the classical (one-parameter) MTS criterion gives significantly different crack propagation angle values in comparison to data calculated numerically by means of the finite element method. It is recommended that several more higher-order terms of the Williams expansion should be used for the description of the stress field as an input for the fracture criterion, especially when the criterion is applied at a large distance from the crack tip. In [10] it is noted that types of structural materials differ in the size of the whole nonlinear zone around the crack tip and also in the proportions of its individual subsets. Due to these effects it is clear that the stress field at a greater distance from the crack tip must be described more precisely. In [12] fracture behavior (crack deflection) of a crack in a silicate-based composite is studied. The three-point bending test is simulated numerically by means of the finite element method and the influence of existence of a stiff aggregate ahead of the crack tip is studied and discussed. The maximum tangential stress criterion is applied in order to calculate the initial crack propagation angle. Various geometrical configurations are investigated (the aggregate eccentricity and depth are varied) and the generalized MTS criterion based on the Williams expansion is tested considering various numbers of initial terms of the series. The authors [12] show that using the generalized MTS criterion is strongly limited by the elastic mismatch caused by the presence of the stiff aggregate near the crack tip. Therefore, as it is concluded in [12] additional extensive investigations are unavoidable in this field in order to find out how the stress field ahead of a crack tip behaves when a material interface exists at a small distance from it.
A new variant of the Element-Free Galerkin (EFG) method that combines the diffraction method, to characterize the crack tip solution, and the Heaviside enrichment function for representing discontinuity due to a crack, has been used in [13] to model crack propagation through nonhomogenous materials. In the case of interface crack propagation, the kink angle is predicted by applying the maximum tangential principal stress (MTPS) criterion in conjunction with consideration of the energy release rate (ERR). The MTPS criterion is applied to the crack tip stress field described by both the stress intensity factor (SIF) and the T-stress, which are extracted using the interaction integral method. The proposed EFG method has been developed and applied for 2D case studies involving a crack in an orthotropic material, crack along an interface and a crack terminating at a bimaterial interface, under mechanical or thermal loading; this is done to demonstrate the advantages and efficiency of the proposed methodology. The computed SIFs, T-stress and the predicted interface crack kink angles are compared with existing results in the literature and are found to be in good agreement. An example of crack growth through a particle-reinforced composite materials, which may involve crack meandering around the particle, is reported.
Nowadays one can estimate the crack propagation direction angle using molecular dynamics simulation [17][18][19][20][21][22]. Molecular dynamics (MD) is a powerful tool in characterizing the inception and evolution of plastic deformation and associated fracture mechanism at the atomic scale [17]. Applied to the problem of crack propagation and growth, the method has been used extensively to analyze the fundamental mechanisms of material's fracture in the past [17][18][19][20][21][22]. Thus, in [18] molecular dynamics (MD) simulations are performed to investigate a crack propagation behavior in ferrite iron with a bodycentered cubic (BCC) crystal structure containing a solute nanocluster under pure and mixedmode loadings. In the MD simulations, coherent Cu and/or Ni nanoclusters are considered to be the solute nanoclusters formed in ferrite iron during neutron irradiation. A MD simulation to assess the effect of crack length on the ultimate tensile strength of infinitely large armchair and zigzag graphene sheets is presented in [19]. The authors of [19] note that the strength of graphene is inversely proportional to the square-root of crack length as in continuum fracture theories. Further comparison of the strength given by MD simulations with Griffith's energy balance criterion demonstrates a reasonable agreement. Armchair and zigzag graphene sheets with 2.5 nm long crack exhibit around 55% of the strength of pristine sheets. Investigation of the influence of temperature on the strength of grapheme performed in [19] indicates that sheets at higher temperatures fail at lower strengths, due to high kinetic energy of atoms. The main result of [19] is that strength obtained from molecular dynamics simulations shows inverse square-root proportionality with crack length as in continuum fracture mechanics theories.
In [20] crack propagation at the Cu/SiC interface under Mode I and mixed mode loading conditions is analyzed by MD simulations. The thickness of the Cu/SiC interface can be defined with monitoring the atomic energy change across the interface and crack growth is observed at strong Cu/SiC interface. For the mode I loading conditions, the crack propagates in an asymmetrical manner. The dislocations generated by lattice mismatch are believed to be responsible for the asymmetric crack propagation. The Rice and Thomson (R-T) model is adopted here as a criterion to identify the crack propagation style by comparing the energy for dislocation nucleation and the decohesion energy of the interface. The style of crack propagation predicted by R-T model is cleave and asymmetrical, which fits well with the results of MD simulations. For the mix mode loading conditions, six different loading angles are simulated by MD method. For the six mixed mode loading conditions, the style of the crack propagation is dependent on the loading angle. The MD simulations proved that cleave and blunt are two typical behaviors of the interfacial cracks, which are consistent with theoretical prediction by R-T model. Such experimental observations and numerical predictions represent a qualitative nanoscale understanding of the fracture processes, providing valuable insights into the crack tip plasticity and associated deformation mechanisms-dislocation nucleation and multiplications [17].
The present study is aimed at the determination of the crack propagation direction angle in a wide range of mixed-mode loading using 1) the multi-parameter maximum tangential stress criterion; 2) the multi-parameter strain energy density criterion; 3) molecular dynamics simulation of the crack propagation behavior under mixed-mode loading in a full range of the mixity parameter.

Continuum linear fracture mechanics approach
In the development of fracture mechanics Williams made a major breakthrough in the analysis of asymptotic stress field at the vicinity of the crack tip in isotropic linear elastic plane media. With the well-known eigenfunction expansion method it is possible to establish the separable variable nature of the solution and to obtain asymptotic expressions for the stress field in a plane medium with a tractionfree crack submitted to mode I, mode II and mixed-mode (mode I and mode II) loading conditions: The multi-parameter fracture mechanics concept consists in the idea that the crack-tip stress field is described by means of the Williams expansion (1). In this work the central crack in an infinite plane medium is considered. Analytical determination of coefficients in crack-tip expansion for a finite crack in an infinite plane medium is given in [3]: for a Mode I crack; 21 12 Since all the coefficients in the Williams asymptotic series expansion for the geometry considered are known it is possible to estimate the crack propagation direction angle by means of the multi-parameter fracture mechanics concept. In this paper two fracture criteria were chosen for the estimation of the initial crack growth direction: maximum tangential stress (MTS) criterion and strain energy density (SED) criterion [2]. The most well-known criterion for the estimation of the crack propagation direction is MTS (maximum tangential stress). This criterion is strictly stress-based and its analysis does not depend on the conditions of plane stress and plane strain. It assumes that a crack extends in the direction of the maximum tangential stress. Mathematically the criterion of the maximum tangential stress is written as:  (8) is tested while considering various numbers of the initial terms of the Williams expansion (1). Therefore, the tangential stress has to be expressed via the power series and then its maximum is sought numerically. The results obtained by the MTS criterion are shown in Table 1 where N is the number of terms keeping in the Williams series expansion, / c r r a  is the dimensionless distance from the crack tip. The first column gives the crack propagation direction angles obtained by the MTS criterion when the leading term of the Williams asymptotic expansion is taken into account. While the following columns show the crack propagation direction angles given by the multi-parametric asymptotic expansion of the stress field in the vicinity of the crack tip. Here the multi-parameter crack tip expansion of the stress field contains 100 terms. One can see that the values of the crack propagation direction angle obtained by the oneparameter fracture criterion are close to the values of the second column especially for Mode II loading. Obviously it can be seen the influence of T-stresses at the crack tip. As the distance from the crack tip increases the difference between results given by the one-parameter fracture criterion and the multi-parameter fracture criterion enhances.
The criterion based on the minimum of the strain energy density deals with strain energy and therefore the conditions for predicting crack growth direction depend on whether or not the specimen is subjected to plane strain or plane stress conditions. Specifically, the minimum strain energy density 5 1234567890 ''"" The Williams series expansion (1)     ). It is obvious that the fracture criteria in its classical oneparameter form are independent on the length parameter, i.e. on the radial distance. The following conclusions can be stated based on the results obtained from the MTS criterion: at larger distances from the crack tip, it is very useful to take into account more terms of the Williams series expansion; at small distances from the crack tip, the crack propagation direction is nearly independent of the number of Williams series expansion terms considered in the fracture criterion, which is caused by the strong singular stress field in the region. At larger distances from the crack tip, the initial crack propagation angle converges more slowly and it is necessary to keep high number of terms in the asymptotic expansion. In the present study 100 terms of the Williams series expansion are kept. Although the dependences in Tables 1-5 are not strictly uniform the qualitative trend is clear; with increasing distance from the crack tip a greater number of higher-order terms of the Williams series expansion is needed. As it is noted in [6] the very important question is how the choice of a right distance from the crack tip should be made. When non-singular terms of the stress expansion are taken into account it is impossible to avoid the dependence of the criterion on length parameter [6]. It is common that this value should be constant for each material and should be connected to materials properties. In order to consider this phenomenon in this study atomistic modeling of mixed-mode loading of the plane medium with the central crack is performed.

Atomistic modeling for mixed-mode loading of the plane medium with the central crack
Face-centered cubic (FCC) materials such as Cu are broadly used in various industries as structural materials [18][19][20]. Thus, in order to improve safety and reliability it is crucial to understand how such materials fail [20]. As noted in [20] even with the same FCC structure different materials exhibit diverse failure patterns due to their dissimilar material properties. The fracture of materials is intrinsically a multiscale phenomenon, which originates at the atomic scale as a result of the breaking of the bonds between atoms. The fracture behaviour of materials is, therefore, strongly dependent on the local atomic environment, such as the atomic structure, lattice orientation, and the discrete nature of matter distribution [18][19][20]. Thus, it is unsuitable to analyze fracture at the atomic scale by a conventional continuum mechanics-based approach. Atomistic scale modeling and simulation are required. Molecular dynamics is one of the most widely used numerical simulation techniques for investigating the fracture behavior of materials at the atomic scale [17][18][19][20][21]. In this paper to model copper plate under mixed loading we used large-scale atomic/molecular massively parallel simulator (LAMMPS) [21,22] in combination with embedded atom method (EAM) [17,24]. EAM potentials are widely used in variety of different simulations, focusing mainly on mechanical properties [17,24]. Potentials for copper were proposed by Foiles et al. [24] and shown decent results in elastic limit, therefore, all the strain considered in this work are elastic. Molecular dynamics simulations were performed with the plate with the central crack in Cu single crystal to investigate the governing deformation mechanisms at the crack tip under mixed-mode (Mode I and Mode II) loadings in a full range of the mixity parameter. Periodic boundary conditions were implemented in all three directions of the cell. To neglect the effect of neighbouring cells we choose the size of the central crack to be relatively small (1:10 ratio) to the size of the simulation cell. Furthermore, we added small noninteracting boundaries to the edge of the plane. Total number of atoms in the cell is 300000. Before mixed loading is applied, plate is optimized to the minimal energy with conjugated gradient method. When minimum energy state is achieved, we apply mixed strain. During all 50000 steps of simulation, we collect data of the state of all atoms in the cell. Results are shown in Figures 1-6, colour coding is obtained by OVITO [26,27] tool. Brighter colours correspond to higher stress.      The stress-strain curve for Cu single crystal is shown in Figure 7. For this system material remains elastic till 5% strain. The subsequent plasticity caused by gliding of dislocations resulting in crack tip blunting (Figures 2, 3, 6). This phenomenon is in consistence with the previous numerical experiments [17]. Figures 1-3 show the deformation process and the crack growth for the mixity parameter . One can clearly see that in the vicinity of the macrocrack microcracks and vacancies appear, evolve and merge. The vacancies coalescence can be seen in Figure 8. The material shows elastic behavior till the first plastic effects begin to develop (Figure 8). and coalescence. Figure 9 and 10 show the effects of temperatures on the crack propagation process. The crack propagation process at 0K is shown in Figure 9. The crack propagation process at 77K is shown in Figure 10. The mixity parameter for calculation illustrated by Figures 9 and 10 is equal to 0.5. It can be seen that the temperature influences on the crack propagation process and increasing the temperature leads to the blunting of the crack.

Conclusions
The paper is focused on the application of the different approaches for the determination of the initial crack propagation angle. Copper plate with the central crack under complex mechanical stresses (Mode I and Mode II loading) is studied by extensive molecular dynamics simulations based on the EAM potential. On the other hand, the complete Williams expansion for the crack tip fields containing the higher-order terms is used. The crack propagation angle is obtained by 1) the multi-parameter fracture mechanics approach based on two fracture mechanics criteria, MTS and SED; 2) atomistic modeling for the mixed-mode loading of the plane medium with the central crack. From our simulations we can get crack propagation directions and crack angles. Calculations of MD method were run for three different values of e M : 0.4, 0.5 and 0.6. Calculated values of crack angles were -51.5°, -46.6° and -42.2° accordingly. All the fracture criteria tested give similar values of the crack growth angle for different values of the mixity parameter. It is shown that the initial crack propagation angles given by the both approaches are very close especially for the case when the higher order terms of the Williams series expansion for the stress/displacement field description are taken into account (Tables 4, 5). Thus one can conclude that the criteria of classical continuum mechanics MTS and SED can give satisfactory predictions for crack initiation direction. The crack propagation direction angles given by the conventional fracture mechanics reasonably agree with the angles obtained from molecular dynamics simulations.
The paper is aimed at the application of the different approaches for the determination of the initial crack propagation angle. The crack propagation angle is obtained by 1) the multi-parameter fracture mechanics approach based on two fracture mechanics criteria, MTS and SED; 2) atomistic modeling for the mixed-mode loading of the plane medium with the central crack. It is shown that the initial crack propagation angles given by the both approaches are very close especially for the case when the higher order terms in the Williams series expansion are kept. One can conclude either that considerable efforts still remain while modeling the mechanical behavior of crystalline materials with empirical potential in the framework of MD simulations. One can conclude that it is crucial and natural to develop multiscale methods [28 -34]. For instance, in [28] the authors try to unify both the continuum mechanics and atomistic approaches and divide the simulation area into atomistic and continuum subdomains such that nanoscale defects are described at atomistic resolution while a continuum mechanics is employed elsewhere. Thus the combination of the continuum mechanics approach and the atomistic simulation allows us to get deeper understanding and realistic description of the crack growth behavior and deformation process under different complex loadings.