Size effects in molecular dynamic simulations of fracture in bcc iron crystals

3.1.1. Edge crack a = 64 d₁₁₀, sample B30

Calculations were performed at an applied stress σ_A = 1 GPa. At this loading level, there was no change in the total number of interactions in MD. Comparison of the MD (o) static stress components Sx, Sy and T-stress at the individual atoms lying on the x-axis (angle θ = 0) with the continuum predictions (─) σ_x, σ_y by anisotropic LFM [15] is presented in figures 2(a) and (b). The commonly used notation x, y, z in LFM is related to our notation (needed e.g. in the appendix) via relations x ≡ x₁, y ≡ x₂, z ≡ x_3. The components ${\sigma }_{{\rm{x}}}=-\mathrm{Re}({\mu }_{1},{\mu }_{2}){K}_{{\rm{I}}}/\sqrt{2\pi r}+T$ and ${\sigma }_{{\rm{y}}}={K}_{{\rm{I}}}/\sqrt{2\pi r}$ for the angle θ = 0 are given by equation (B.2) in the appendix, where r means the distance from the crack tip. The stress intensity factor is ${K}_{{\rm{I}}}={F}_{{\rm{I}}}(a/W){\sigma }_{{\rm{A}}}\sqrt{\pi a}$ and the T-stress is defined as $T={f}_{{\rm{I}}}(a/W){\sigma }_{{\rm{A}}}.$ Here F_I(a/W) and f_I(a/W) stand for the boundary corrections factors: F_I(a/W) = 2.2454 by table 1 and f_I(a/W) = −0.5482 by [12, 13]. The value Re(μ_1, μ₂) = −1.2868 valid for plane stress [16] is considered, since by figure 2(c) plane stress state with S33 ∼ 0 prevails in the larger part of the sample B30. The quantity d₀ = a₀√2 in figure 2 denotes the distance between atoms in the $[\bar{1}10]$ direction. The circles (o) denote MD results and the lines (─) the LFM results. After the surface relaxation, some initial stresses arise in MD at the crack front. They differ somewhat from [5] due to the different sample length L here and in [5]. The initial values at the crack tip are S10 = 2.1850 GPa, S20 = −0.9352 GPa, S30 = −0.8587 GPa and they are relatively small in comparison with the maximum values at the crack tip under a loading. The initial stress is subtracted when comparing MD results with LFM solutions (that do not consider an initial stress). Thus, the values from MD (o) in figures 2(a) and (b) are taken as Sx = S11−S10, Sy = S22−S20, and T_MD = Sx + Re(μ₁, μ₂) Sy = Sx–1.2868 Sy. The relation for T_MD follows from equation (B.2) in the appendix B. The next figure 2(c) shows the instantaneous value of the stress component S33 from MD on the entire x-axis (at B/2) with the help of the continuous lines, at the applied stress of 1 GPa. The dots (S33) are calculated from the stress–strain relation for plane strain when ε₃₃ = 0. It leads to S33 = −s₃₁ S11/s₃₃−s₃₂ S22/s₃₃ = 0.6197 S22−0.0365 S11 with the compliances s_ij for the used potential—see appendix A. Figure 2(c) for the sample B30 also shows how the plane strain state with S33 > 0 (i.e. σ_z > 0) is approached in the nearest vicinity of the crack tip where the peak occurs. In figures 2(a) and (b), one may see very good agreement between the static atomic stress Sx (o), Sy (o) and LFM components σ_x (─), σ_y (─) in the interval r/ d₀ = <1, 10>. As to T-stress, the constant isotropic static solution T = f_I (a/W) σ_A = −0.5482 GPa (the horizontal blue line in figure 2(b)) from LFM complies with MD results at larger distances r/d₀ from the crack tip. The anisotropic MD (o) results for the T-stress differ only near the crack tip, where three-axial stress state exists - see figure 2(c). In spite of the different sample length L here and in [5], the agreement between MD and LFM in figures 2(a), (b) is very good, similar as for the longer crystal in [5].

Zoom In Zoom Out Reset image size

3.1.2. Plane strain conditions in larger samples

Figure 3(a) for edge crack a = 128 d₁₁₀ in sample B60 shows that under linear loading of 4.5 kN s⁻¹ at the applied level of 1 GPa the instantaneous stress component S33 even overcome the value S33 (dots) following from the stress–strain relation ε₃₃ = 0 for plane strain. However, in this case, a small change of interatomic interactions LINT was detected in MD (the relative number of broken bonds in the entire 3D atomistic system is ∼ 0.00004%). The situation under lower loading of 0.5 GPa (with no change of LINT) illustrates figure 3(b) from the sample B70. Figures 3(a) and (b) illustrate that with increasing thickness B, plane strain conditions are better obeyed at the crack front (in comparison with figure 2(c) for B30).

Zoom In Zoom Out Reset image size

Here we present MD results under linear loading in time up to the final fracture for the two loading rates of 4.5 kN s⁻¹ and 5.85 kN s⁻¹ corresponding to 5 mm min⁻¹ in experiment and either to Young's modulus E = 1/s₂₂ = 135 GPa for plane stress or to E = 1/A₂₂ = 175.5 GPa for plane strain state. This corresponds to the used potential with the compliance coefficients s₂₂ = 0.7410 × 10⁻¹¹ m²/N and A₂₂ = 0.5698 × 10⁻¹¹ m²/N for our crystal orientation from figure 1. The results are summarized in tables 2 and 3, where the atomistic samples of the different self-similar size (i.e. with the same ratio L/W, B/W and ∼ a/W) are marked according to the different thickness B in table 1 (i.e. as B30, B40, B50, B60 and B70). The time of crack initiation t_i and dislocation emission t_e are defined in section 2 and illustrated in figures 5(a), (b) by the arrows. When final fracture of the atomistic sample occurs, at this moment RK decreases to zero in the dependencies RK-time and RK-COD—see e.g. Figures 11(a), (b) for the sample B60. This moment defines the time t_f and the magnitude COD_f included in table 2.

Similar to [5], our linear loading is of a quasi-static character, described in table 3 by the ratio ∆t_G /(T/2) where ∆t_G is the critical time up to Griffith level and T/2 is the half period of the longest free longitudinal vibrations of the atomic sample, given by the relation λ/2 = L = C_L001 (T/2), where C_L001 is the velocity (5550 m s⁻¹) of the longitudinal waves in the [001] direction for the used potential - see table 2 in [17] . The quasi-static character of loading will be illustrated at the end of this section in the framework of the concept for virial stress. Note that the quantity T/2 also represents the flight time of the loading waves needed to travel the length L of the atomic sample.

In all the cases, the crack initiation below Griffith level in the middle of the crystal (at B/2) was of a cleavage character and this was supported by the surface dislocation emission on oblique slip systems $\left\langle 11\bar{1}\right\rangle \{011\}$ and by subsequent cross slip of the emitted dislocations into oblique slip systems $\left\langle 11\bar{1}\right\rangle \{112\}.$ In the case of the thin sample B30 it leads to a fast brittle fracture, with many cleavage facets on fracture surfaces. For thicker samples B40, B50, B60 and B70, in addition to the surface emission, we newly observed dislocation emission from the crack front also inside the tested crystals (i.e. bulk crystal emission) on the slip systems $\left\langle 11\bar{1}\right\rangle \{011\},$ $\left\langle 11\bar{1}\right\rangle \{123\}$ and twin generation in the slip system $\left\langle 11\bar{1}\right\rangle \{112\}.$ All the slip systems are oblique to the crack front [110]—see figure 4. This leads to larger atomistic fracture toughness shown in table 3 and to more ductile fracture in comparison with the specimen B30.

Zoom In Zoom Out Reset image size

A scheme for the above mentioned oblique slip systems observed in MD is shown in figure 4.

The nominal shear (slip) stress τ from the applied stress σ_A (Schmid factors) for the slip systems is given by the following relations in our case:

$$\begin{eqnarray}&&\tau =0.41{\sigma }_{{\rm{A}}}\,{\rm{for}}\,\left\langle 11\bar{1}\right\rangle \{011\},\,\tau =0.47{\sigma }_{{\rm{A}}}\,\text{for}\,\left\langle 11\bar{1}\right\rangle \{112\},\,\tau =0.46{\sigma }_{{\rm{A}}}\,{\rm{for}}\,\left\langle 11\bar{1}\right\rangle \{123\}\end{eqnarray} \tag{ 2 }$$

The stress barriers for dislocation generation in the individual slip systems with the many-body potential in use are:

$$\begin{eqnarray}&&{\tau }_{{\rm{disl}}}=14.5{\rm{GPa}}\,{\rm{for}}\,\left\langle 11\bar{1}\right\rangle \{011\}\,,{\tau }_{{\rm{disl}}}=16.3{\rm{GPa}}\,{\rm{for}}\,\left\langle 11\bar{1}\right\rangle \{112\},\,{\tau }_{{\rm{disl}}}=16.2{\rm{GPa}}\,{\rm{for}}\,\left\langle 11\bar{1}\right\rangle \{123\}\end{eqnarray} \tag{ 3 }$$

These values follow from the results of block like shear simulations in perfect bcc iron crystals, presented in [18, 19, 20]. Further details are given below.

As mentioned above, the individual quantities in table 2 mean: t_i is the time of crack initiation (see figure 5(a)) detected in the middle of the crystal (i.e. at B/2), t_e is the time of the surface dislocation emission(see figure 5(b)) detected at the crack tip in the first surface layer (WL in figure 1) and t_f is the time of final fracture of the crystal under prescribed loading rate dσ/dt at the lower and upper BW borders in the third column in table 2. The time of fracture t_f is determined from the atomistic tension curves RK-time (see e.g. figure 11(a) for B60) as the time when the resultant force RK decreases to zero at the final fracture. The dependence RK - COD (figure 11(b) for B60) serves to detect the crack opening COD_f at the final fracture, illustrating the brittle or ductile behavior of the individual atomistic specimens in table 2: the small COD_f for B30 indicates the brittle fracture, in larger samples the larger COD_f shows that the character of the final fracture was ductile, unlike the brittle crack initiation in all the cases. The quantities σ_i, σ_e and σ_f correspond to the external stress at the BW borders at the corresponding time t_i , t_e and t_f during the linear time loading, described in the third column in table 2. As mentioned in section 2, the maximum of the resultant force RK in the atomistic tension curves represents the fracture resistance of the atomic sample in MD simulation and serves to determine the stress ${\sigma }_{{\rm{MD}}}=R{K}_{\max }/BW$ and subsequently the fracture toughness of the atomic sample ${K}_{{\rm{MD}}}={F}_{{\rm{I}}}{\sigma }_{{\rm{MD}}}\sqrt{\pi a}$ inserted in table 3. The quantities K_i , K_e and K_f in table 3 are determined in the same way from σ_i, σ_e and σ_f defined above.

Zoom In Zoom Out Reset image size

The surface emission and 3D visualization of the emitted dislocations is shown in figure 6. One may see in table 2 that t_e < t_i , i.e. the surface dislocation emission on to oblique slip systems $\left\langle 11\bar{1}\right\rangle \{011\}$ always precedes and facilitates crack initiation under the applied stress σ_i, smaller than the static Griffith stress σ_G in table 1. When the red dislocation lines reach screw character (parallel with $\left\langle 11\bar{1}\right\rangle$ directions) like in figure 6, they make cross slip to oblique {112} planes, leading to separation of the (001) planes. Later after multiple surface emissions, the cross slip causes a pre-cleavage in the middle layers (at B/2) of the sample B30, similar to figure 9(b) in [5]. This enables the fast brittle fracture, with many cleavage facets on fracture surfaces in figure 7. This failure mechanism is described in detail in [5], including the stress analysis. As mentioned above, the surface dislocation emission and the subsequent crack initiation are similar in all the tested atomistic specimens. This is in agreement with the stress analysis according to anisotropic LFM at Griffith level of loading, presented in [16] and [5].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The next figure 8 shows the cross sections along the thickness B in the 3D atomic samples B30 and B60. The figure 8 enables to explain why the cross slip $\left\langle 11\bar{1}\right\rangle \{112\}$ cannot contribute to cleavage in thicker samples so effectively as in the case of the thin sample B30. The slope of the oblique plane {112} with respect to our observation plane (110) is tg α = d₀₀₁/d₁₁₀ = 1/√2 - see figure 4. The first cross slips are visible as the horizontal slip patterns in figures 8(a), (b) for B30 and B60. The first patterns from the cross slip are always at the distance 4a₀ in the Y = [001] direction from the crack plane after the first emission. Due to the slope, they reach the crack plane (001) near by the layer 8 in the Z = [110] direction parallel with thickness B and from the opposite site near by the layer B-8—see figures 8(c), (d) for B30 and B60. (The number n = 8 follows from the relation tg α = 4a₀ /n d₁₁₀ = 4a₀ /(n a₀√2/2) = 1/√2 ↔ n = 8.) That's obvious that with the increasing applied stress, the failures at the positions 8 and B-8 (caused by the cross slip) extend to the crack plane position B/2 easier in the thin sample B30 in comparison with the thicker specimen B60. Moreover, after the multiple surface emissions like in figure 9(a) in [5], the lowest and upper horizontal lines, coming from the cross slip, are at the distance of about 8a₀ from the crack plane and they reach the crack plane (001) in the layers 15-16 due to the slope of the {112} planes (tg α = 8a₀ /(n a₀√2/2) = 1/√2 ↔ n = 16). This is the same in the samples B30 and B60. In B30, this supports the brittle crack development in the middle (at B/2), but in B60 not, since the layers 15-16 do not lie in the middle at B/2. Thus, the cross slip $\left\langle 11\bar{1}\right\rangle \{112\}$ cannot contribute to cleavage in thicker samples so effectively as in B30. This is similar for B40 or B50 and it is one source of the increased atomistic fracture toughness K_MD in the thicker specimens.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Probably more important reasons for more ductile behavior (and the increased fracture toughness K_MD) with the increasing sample thickness are the new bulk shear processes generated inside thicker samples, presented below.

Figure 9 shows the situation in the sample B40. Figure 9(a) illustrates that beside the surface dislocation emission, the oblique dislocations $\left\langle 11\bar{1}\right\rangle \{011\}$ are generated also inside the crystal in the middle layer 20. Later, formation of twins in oblique slip systems $\left\langle 11\bar{1}\right\rangle \{112\}$ is observed at the crack front and also crack deflections—see figure 9(b). These processes hinder crack growth in the original crack plane (001). When the oblique twins $\left\langle 11\bar{1}\right\rangle \{112\}$ (oTwins) reach the right borders in figure 9(c), new stress concentrators arise on the right free BL surface and new slip processes are generated from this surface, which increases COD. This is the main reason of the larger COD_f in table 2 for the specimens B40, B50, B60 and B70, compared to B30. The last figure 9(d) shows that the atomic shear stress τ₁₁₂ in the lower slip system $[11\bar{1}](112)$ overcome the stress barrier τ_twin = −9.3 GPa for twinning with the used potential [18], even before crack initiation. (Note that after crack initiation all the stresses vanish, including τ₁₁₂ in figure 9(d), since the former crack tip atom already lies on the free crack surface.) Oblique twinning generated by our crack (001)[110] is possible under Griffith loading (and above) as well according to anisotropic LFM as presented in [16]. It is well known that twinning in bcc is also possible through the dissociation of the dislocations according to continuum models [21, 22], because it is energetically favorable. This process was observed as well in our simulations, as figure 9(c) indicates.

In addition to the all above mentioned slip processes, new oblique dislocation emissions of the type 〈111〉{123} were observed at the crack inside the crystals B50 and B60. Figure 10(a) show the slip patterns (the blue dark atoms) arising after dislocation emission $[\bar{1}\bar{1}1](123)$ at the crack front toward the upper part of the crystal. Recognition of the slip patterns from figure 10(a) is possible via the block like shear (BLS) simulations in small perfect 3D cubic crystals - see figure 10(b) from BLS for the case in figure 10(a). As well the emission on other equivalent oblique slip systems 〈111〉{123} (with the same Schmid factor of 0.46) were detected, as shown in figure 10(c). The lower slip patterns (below the crack face) correspond to scheme in figure 4, i.e. apparently to the $[11\bar{1}](123)$ slip system. Occurrence of these slip processes correlates with the second highest maximum in the dependences RK - time, which illustrates figure 11(a) for the atomic sample B60. The dependence RK versus COD in figure 11(b) is an atomistic analogy to the known dependencies P - ∆L from fracture experiments. Figure 11(c) shows the fracture surface of the sample B60 after the final fracture under loading rate 5.85 kN s⁻¹, where the individual observed microscopic processes are marked by the lines. The oblique emissions 〈111〉{123} have been observed before for different crack orientations, e.g. in [20] for a crack orientation (110)[001] (crack plane / crack front). Comparison of the fracture surfaces in figure 7 for B30 with figure 11(c) for B60 illustrates well the difference in microscopic processes participating during fracture of the thin (smaller) and thick (larger) atomistic specimens used for MD simulations. The line {011} in figure 11(c) denotes atomic configuration in an oblique {011} crystalline plane arising after crack deflection along the oblique slip system $\left\langle 11\bar{1}\right\rangle \{011\}.$ During dislocation emission in this slip system also normal separation between the planes {011} often arises (see e.g. lower part in figure 9(b)), causing the formation of normal stress σ_n between the planes {011}. When σ_n reaches the cohesive strength σ_coh in the normal direction [011], then decohesion in the slip system may occur, leading to crack deflection onto a {011} plane as marked in figures 11(c) or 9(c). Such decohesion is also possible according to a continuum model presented in [23].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

According to stress analysis in the appendix B for the slip system $[11\bar{1}](123)$ viewed in figure 4, the slip stress on the (123) plane, acting in the direction of the Burgers vector ${\rm{b}}={{\rm{a}}}_{0}/2\,[11\bar{1}]$ below the crack tip in the lower part of the crystal is given by the relation (derived in the appendix B):

$$\begin{eqnarray}&&{{\tau }_{{\rm{b}}}}^{123}=-\displaystyle \frac{1}{2\sqrt{21}}{\sigma }_{12}\left(r,\theta ,{\mu }_{1},{\mu }_{2}\right)-\displaystyle \frac{\sqrt{3}}{\sqrt{14}}{\sigma }_{22}\left(r,\theta ,{\mu }_{1},{\mu }_{2}\right)+\displaystyle \frac{\sqrt{3}}{\sqrt{14}}{\sigma }_{33}\left(r,\theta ,{\mu }_{1},{\mu }_{2}\right)\end{eqnarray} \tag{ 4 }$$

where σ_ij (r, θ, μ₁, μ₂) are the LFM stress components related to our original coordinate system x₁ = $[\bar{1}10],$ x₂ = [001], x₃ = [110], the angle θ (measured from the x₁-axis in the anti-clockwise direction) is +θ = 3/2 π + β = 2π–74.207° for the lower slip system $[11\bar{1}](123)$ in figures 4 and in 10(c), (r) is the distance from the crack tip in this slip system and μ₁, μ₂ are the complex roots following from the compatibility equations in anisotropic continuum owing the cubic symmetry. Using the roots μ₁, μ₂ for plane strain from [16] (i.e. μ₁ = 0.6167+0.8653i, μ₂ = −0.6167+0.8653i) and the angle θ = 2π−α, α = π/2−β = 74.207°, then the stress components in the original system x₁ = $[\bar{1}10]$], x₂ = [001] and x₃ = [110] are:

$$\begin{eqnarray}&&{\sigma }_{22}=1.0957{K}_{{\rm{I}}}/\sqrt{2\pi r},\,{\sigma }_{11}=0.4866{K}_{{\rm{I}}}/\sqrt{2\pi r},\,{\sigma }_{12}=0.2796{K}_{{\rm{I}}}/\sqrt{2\pi r}\end{eqnarray} \tag{ 5 }$$

The component ${\sigma }_{33}=-({s}_{31}/{s}_{33}){\sigma }_{11}-({s}_{32}/{s}_{33}){\sigma }_{22}=-0.0376{\sigma }_{11}+0.6197{\sigma }_{22}$ follows from the relation ε₃₃ = 0 for plane strain. Contribution from the T-stress (acting along the x₁-axis) to the slip stress τ_b ¹²³ is zero according to Schmid law because the direction x₁ = $[\bar{1}10]$ is perpendicular to the direction of the Burgers vector ${\rm{b}}={{\rm{a}}}_{0}/2\,[11\bar{1}].$ (For more detail see appendix B). If we consider a near distance from the crack tip, e.g. r = b/2 (where b = a₀√3/2), and Griffith level of loading K_I = K_G = 0.9115 MPa m^1/2 for plane strain, then magnitude of τ_b ¹²³ ∼ 19.7 GPa which overcomes the stress barrier τ_disl = 16.2 GPa for dislocation emission in the $\left\langle 11\bar{1}\right\rangle \{123\}$ slip system with the potential in use. This means that the emissions observed in MD (and presented in figure 10 for the samples B50 and B60) are possible as well from the point of the anisotropic linear fracture mechanics.

The upper slip system $[\bar{1}\bar{1}1](123)$ shown in figures 10(a), (b) with the angle θ = 1/2 π +β is less favorable for the emission in comparison with the lower slip system $\left\langle 11\bar{1}\right\rangle \{123\}$ presented above (figure 10(c)), in spite the fact that Schmid factor (0.46) for the upper and the lower slip system is the same. This 'non-Schmid behavior' can be explained in our case by the anisotropic LFM stress field at the crack tip, described by equation (B.2) in the appendix. Nevertheless, the magnitude of τ_b ¹²³ (square roots of the components in x₂ and x₃) at r = b/4 corresponds to ∼17. 4 GPa, which exceeds the stress barrier τ_disl = 16.2 GPa in the slip system, i.e. this emission at Griffith level of loading is also possible according to the anisotropic LFM. The details are given in appendix B. Note that the distance b/4 corresponds to the position where the theoretical stress barrier for the slip system $\left\langle 11\bar{1}\right\rangle \{123\}$ (see figure 7 in [20]) reaches the maximum during BLS simulations in the perfect bcc iron crystal.

3D visualization of the emitted dislocations in the specimen B60 is shown in figure 12 via a detail near the crack. Here the dislocation cores are visualized again via the red atoms with coordination number KNT = 16, the yellow atoms show the crack on (001) planes with coordination number KNT = 9 and the blue atoms (KNT = 10) are the surfaces atoms lying on {110} free surfaces, similar to figure 6. The more complex dislocation structure in figure 12 corresponds to the situation near the maximum RK in figure 11(a) at the time step 18800 h, where the slip systems $\left\langle 11\bar{1}\right\rangle \{123\}$ are active.

Zoom In Zoom Out Reset image size

The same slip process and similar crack behavior as for B60 have been observed in the atomistic sample B70. The fracture characteristics are included in tables 2 and 3. The atomistic fracture toughness K_MD for B70 slightly overcomes K_G for plane strain.

We should note that whatever the stress barriers for the slip systems 〈111〉 {123} and 〈111〉{112} are almost the same - see equation (3), one cannot expect for our crack orientation (001)[110] the emission of the blunting dislocations 〈111〉{112} because the inclined systems are oriented in the easy twinning direction where τ_twin (9.3 GPa) is significantly smaller than τ_disl (16.3 GPa) and thus, twinning is preferred. Emission of the blunting dislocations 〈111〉{112} is favored (see e.g. [10]) at the crack $(\bar{1}10)[110]$ where the slip systems are oriented in the hard anti-twinning direction with τ_anti-twin ≫ τ_disl—see [15]. The topic twinning versus anti-twinning in bcc metals was recently studied also in [24]. Generally, twinning play very important role in mechanical properties of the bcc metals - see e.g. [25].

The results presented above illustrate that MD simulations can bring new information on microscopic processes generated by the crack itself. On the other hand, there are also some limitations of the model.

One limitation concerns the wave processes generated by crack initiation, dislocation emission, twin generation, etc Each such a stress relaxation at the crack induces emission of the elastic waves, often called as the acoustic emission (AE) or radiation [26, 27]. Propagation of these waves in anisotropic medium is very complex. When the waves arrive to the borders of the finite MD samples (or experimental specimens), they reflect back to the bulk crystal and may influence the situation at the crack front. After a first bond breakage in the system, MD results are continuously influenced by the back wave reflections along the crack front in the [110] direction parallel with the small thickness B. When comparing MD results with continuum predictions, such a situation have to be avoided because the continuum predictions mostly do not include the back-wave reflections. (That's why people often use periodic boundary conditions, but it also brings also some wrong side effects, influencing the brittle versus ductile behavior in MD, as mentioned e.g. in [5].) When comparing MD with experiment, the topic is not so important, because the back wave reflections act as well in experimental specimens. However, their influence in the experiment is probably less intensive in comparison with the small atomistic 3D samples used in MD simulations. The velocity of the longitudinal and transversal (shear) waves in the important (pure) crystalline directions for our crack orientation is given in [17], our sample dimensions are given in table 1. Simple analysis of the wave motion shows that the data for K_e and K_i are not influenced by the back wave reflections from the BW loaded borders and the right free BL border in the [001] and $[\bar{1}10]$ direction respectively. They are influenced only by continuous back wave reflections along the crack front in the [110] direction parallel with the small thickness B, starting from the first bond breakage in the sample. As to K_MD (derived from the maximum RK - see figure 11(a)), it can be already influenced slightly by the back wave reflections in the [001] and $[\bar{1}10]$ directions, similar to fracture toughness K_IC derived from the maximum P in a real fracture experiment (because the velocity of the elastic waves are the same and the geometry of MD samples and experimental samples from [10] is similar).

As to reproducibility of MD results, it is excellent when sequence processing is used for B30. Parallel and sequence processing of the sample B40 led practically to the same MD results. However, the repeated parallel runs with B60 for loading rate 4.5 kN s⁻¹ showed somewhat different results in the non-linear region of the tension curve RK - time, behind the maximum RK, which also lead to different values COD_f and K_f for final fracture. This can probably be explained by a fact that our interatomic forces are nonlinear. This enables a partial absorption of the local kinetic energies of the individual atoms into random thermal atomic motion, namely after the back wave reflections behind the maximum RK. This could explain the wrong reproducibility of parallel processing behind the maximum RK. We believe that the presented MD results are relevant up to K_MD, but behind it, only up to the moment when the oblique twins reach the free BL borders as in figure 9(c). This can be recognized in the dependence RK - time via the flat part, marked in figure 11(a) by the arrow 'BL border'.

As stated in the Introduction and section 2, this study is limited to one particular relative geometry L/W, B/W and to one ratio a/W ∼ 0.42 because of our next fracture experiments on bcc iron crystals with the same relative ratios. Under the same relative geometry, the size of the finite 3D atomic samples has been varied. Our question (before modeling the experiment) was how the different magnitudes L, W and B (see figure 1) affect MD results generally, and especially behavior of the crack and the atomistic fracture toughness. It should be noted that there are also other MD studies, e.g. [28] where the influence of different ratios a/W on fracture toughness is examined, or in [6] where only the different thickness is tested, etc In the self-similar LFM continuum model, the influence of the relative ratios, including a/W, on the stress intensity factors at a crack is treated via the boundary corrections factors presented in [11–13] and for our particular case in table 1. Our study utilizes this self-similar LFM approach. Figure 2 (comparing the local atomic interplanar stress with LFM solution) confirms that this continuum LFM model, designed originally for macroscopic bodies, is acceptable even for a (sufficiently large) 3D atomic sample of nanoscopic dimensions, which is fascinating. (Let us note that the safety assessment in nuclear power plants is also based on the self-similar LFM concept - see e.g. [1]).

Further, when a fast loading or an increased temperature in MD simulations is applied, then by virial theorem the two components of the virial stress should be considered: coming from the potential energy and from the kinetic energy in the system. It should be noted here that the definition of the stress on the atomic level coming from the potential energy in the system is not unambiguous [29]. The different topics may require the different concepts [30] of the stress calculations on the atomic level. This problem is discussed also in [9]. For the part of the stress originating from the potential energy we used here the concept of the local interplanar stress from [9], because compared to the local average volume stress derived from the strain energy density, the interplanar stress better describes situations when a sharp stress gradient or inhomogeneous strain distribution appears within the range of ​​cohesive forces. Important examples are interfaces such as the near stress field at the crack tip presented in figure 2, further crack initiation and dislocation emission, investigated in this study. (At free surfaces of the sample, the normal component of the local average volume stress oscillates, while the corresponding local interplanar stress decays smoothly to zero, as expected.) Under homogenous strain within the range of cohesive forces in 3D bcc lattice, the two different definitions of the atomic stress lead to the same results [9].

The situation related to dislocation emission and crack initiation is illustrated by the virial stress in figures 13(a), (b) for the surface dislocation emission $[11\bar{1}]\,\left(011\right)$ and via figures 14(a), (b) for the crack initiation in the bulk of the crystal. The components coming from the kinetic energy at an atom l_i are given e.g. in [3b] by Lustko. In this study they are calculated according to [3b] as SIJ_V (l_i ) = m_Fe V_I(l_i )V_J(l_i ) / v₀, where V(l_i ) means the velocity of the atom l_i , v₀ = a₀ ³/2 is the volume per one atom in the bcc lattice and m_Fe is the mass of one iron atom. In the case of the crack initiation, the atom l_i is the crack tip atom in the middle layer at B/2, while in the case of dislocation emission $[11\bar{1}]\,\left(011\right),$ the chosen atom l_i is the moving atom below the crack tip in the surface layer 1 (marked in figure 8(a) via green). The interplanar shear stress τ_b ⁽⁰¹¹⁾ in figure 13(a) is calculated from the local instantaneous stress components SIJ (l_i ) in MD according to equation (3) in [5], i.e. τ_b ⁽⁰¹¹⁾ = (S22 - S33 + S12 /√2) /√6. The velocity dependent part coming from the kinetic energy is calculated analogically as τ_bV ⁽⁰¹¹⁾ = (S22_V - S33_V + S12_V /√2) /√6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 13 shows that up to the time of the surface dislocation emission, where the peak in figure 13(a) occurs, the emission process is driven mainly by the interplanar stress component coming from the potential energy. The peak value in figure 13(a) slightly precedes the time t_e = 6869 h in table 2 and its magnitude corresponds to τ_MD = 15.1 GPa , which overcomes the stress barrier τ_disl = 14.5 GPa in equation (3) for the slip system 〈111〉 {011} following from BLS simulations in [19] for this slip system in a perfect bcc iron crystal with the potential in use. This implies that the dislocation emission in MD is realized in accordance with the Rice model [7b] for 0 K. The level of τ_bV ⁽⁰¹¹⁾ in figure 13(b) coming from the kinetic energy before the time t_e is very low. This illustrates that the loading in MD has a quasistatic character (since ∆t_G > T/2 in table 3) and that the kinetic component τ_bV ⁽⁰¹¹ cannot contribute to dislocation generation significantly. The first smaller peak of magnitude 0.125 GPa is located on the time axis at t ∼ 6820 h. It describes the sharp increase of the relative interplanar shear displacement U10 in figure 5(b) in the emission interval (0.5b - 1b) during ∼ 400 time steps, i.e. ∼ 4 ps on time scale. The larger second peak at time step ∼ 7400 signalizes generation of the second dislocation emission according to U10/b and the graphical treatment of MD results at time step 7400, etc.

As to crack initiation, LFM theory assumes that a significant role in the brittle or ductile behavior of the crack plays the three-axial (hydrostatic) stress—see e.g. [13]. In MD it means S = (S11+S22+S33)/3 for the interplanar stress component and S_V = (S11_V + S22_V + S33_V)/3 for the part coming from the kinetic energy. Time development of the three-axial stress at the crack tip in the middle (at B/2) of the crystal is shown in figures 14(a), (b). This is related to the crack initiation presented in figure 5(a).

Figure 14 confirms that the crack initiation in MD is also driven by the interplanar stress components in the part a) and not by the component S_V coming from the kinetic energy, since S_V is negligible before the time step 6966 in the part b). The largest peak in S_V on the time scale coincides well with the time step at the crack initiation in figure 5(a). It indicates that the peak of ∼ 0.5 GPa arises after the sudden transition of the crack tip atom from the bound state to the free state of an vibrating atom lying on the free crack surface (001). This generate a strong pulse (with duration of about 2 ps by figure 14(b)) of the elastic waves that can be visualized e.g. via mapping of the local kinetic energies of the individual atoms l_i in the 3D crystal. (Note that to recognize which part of the free vibrations behind the peak S_V represents the propagating elastic waves [31] and which part represents already thermal vibrations is not easy and it is out of topic of this study.) The crack initiation in figure 14(a) is easy to recognize since when the crack tip atom li already occurs on the free surface (001) (with S22(li) = 0. and KNT = 9), the all normal stress components are set (artificially) to zero in the MD code. But before it, the true peak value S_MD = 32.48 GPa at the moment of the crack initiation in figure 14(a) approaches the ideal cohesive strength σ_coh = 32.42 GPa in the perfect bcc iron crystal under tension I in the [001] direction and plane strain conditions, where the stress state is also three-axial. This theoretical ideal stress–strain curve for the perfect bcc lattice under plane straining is presented in figure 15. Here, beside the prescribed straining ε₂ (STRAIN) in the x₂ = [001] direction, also the lateral contraction ε₁ = –ν ε₂ is considered in the x₁ = $[\bar{1}10]$ direction, where Poisson ratio is ν = –A₁₂/A₂₂ = 0.4675 and A_IJ are the reduced compliances for plane strain (ε₃ = 0) presented e.g. in [16] and also mentioned in the appendix A. The theoretical cohesive stress in the positive [001] direction (STRESS) is then derived from the strained interatomic forces, related to the area per one atom on (001) plane.

Zoom In Zoom Out Reset image size

MD results above indicate that the criteria by three-axial (hydrostatic) stress S_MD ≈ σ_coh for crack initiation and τ_MD ≈ τ_disl for dislocation emission can be transferred to more macroscopic models, like crystal plasticity or finite element method.

The reliable values K_i , K_e and K_MD in dependence on sample thickness B (i.e. on sample size) from our 3D atomistic simulations in bcc iron crystals with edge crack (001)[110] by MD technique are shown in figure 16. The lower line denotes Griffith stress intensity factor K_G for plane stress, and the upper line represents K_G for plane strain. It is obvious that the curves have a tendency to be saturated with the increasing B, similar to what reported for fcc nickel in [6]. The values K_i for crack initiation in thicker samples lie below K_G because the surface dislocation emission facilitates the crack initiation, as mentioned above. Linear fracture mechanics and as well fracture experiments (see e.g. ref. [1b]) consider two different fracture toughness for plane stress and plane strain state This is in a qualitative agreement with our results presented in figure 16: with increasing thickness, the plane strain conditions in the atomic 3D samples are better met (see figure 3) and the atomistic fracture toughness K_MD in figure 16 approaches the higher K_G for plane strain.

Zoom In Zoom Out Reset image size

In relation to experimental observations, the all named slip processes from MD are often observed in bcc iron and ferritic steels - see e.g. the references [1], or recently [32]. As for fracture experiments with edge crack (001)[110] presented in [5] and [10], the MD model with the minimum thickness B30 is closest to reality, because it describes well the brittle fracture under fast loading observed in [10] with the shorter SEN specimens (treated in this study) and as well with the longer SEN crystals (treated in [5]). This may seem paradoxical, but on the other hand, the atomic sample B30 obey LFM predictions concerning the static stress state at the crack tip in elastic region of loading - see figure 2. Under higher dynamic loading, the MD results in figure 13 and figure 14 comply with continuum ideas about crack initiation and dislocation emission. Likewise, the atomic sample B30 already fulfills the condition B/b > 20, which is necessary for plausible modeling the dislocation emission from the crack tip at elevated temperature [7].