Ductile-brittle transition at edge cracks (001)[110] in BCC iron under different loading rates in mode I: a 3D atomistic study

The comparison of the stress field by LFM and MD was performed at applied stress level of 1 GPa. The applied stress is within the linear (elastic) region of the theoretical stress—strain dependence for the potential used, presented in figure 1 in [14]. The faster linear loading 5 PH (2.92 kN s⁻¹) showed no change in the total number of interactions LINT monitored via MD code. Stress calculations on the atomistic level were performed using the inter-planar concept presented in [22]. We stopped the linear loading at time step 6332 h at the level of σ_A = 1 GPa and using the critical viscous damping in Newtonian equations of motion we obtained the static atomic configuration in the 3D atomistic sample loaded in mode I. Comparison of MD static stress components S_x , S_y and T-stress on the x-axis (angle θ = 0) with the continuum predictions by anisotropic LFM [6] for the (001)[110] edge crack

$$\begin{eqnarray}&&{\sigma }_{x}=-\mathrm{Re}({\mu }_{1}{\mu }_{2}){K}_{I}/\sqrt{2\pi r}+T,\,{\sigma }_{y}={K}_{I}/\sqrt{2\pi r}\,and\,\mathrm{Re}({\mu }_{1}{\mu }_{2})=-1.2868\end{eqnarray} \tag{ 1 }$$

is presented in figures 3(a) and (b). Here, ${K}_{I}={F}_{I}(a/W){\sigma }_{A}\sqrt{\pi a},$ further ${\rm{T}}={f}_{I}(a/W){\sigma }_{A}$ and the function Re(μ₁ μ₂ ) of the complex variables μ_1, μ₂ depends on plane strain or plane stress state in the atomistic sample. Figure 3(c) shows the instantaneous stress component S33 from MD on the x-axis at the applied stress of 1 GPa. Here it is obvious that plane strain state with σ_z > 0 prevails only in the nearest vicinity of the crack tip where the peak occurs, while in the larger part of the sample plane stress state prevails with σ_z ∼ 0. For that reason we consider $\mathrm{Re}({\mu }_{1}{\mu }_{2})=-1.2868$ valid for plane stress [35]. The boundary correction factor is ${F}_{I}=2.2434$ by [37] for ${(a/W)}_{MD}=0.4267,$ while ${f}_{I}=-0.5482$ from [18, 19]. The LFM solution is denoted in figures 3(a) and (b) by lines, MD results via circles (o) and r means the distance from the crack tip on the x-axis, oriented in the $[\bar{1}10]$ direction. Regarding MD, some initial stress arises at the crack front in the atomistic sample after the surface relaxation. The maximum values at the crack tip of the longer crack are $S10=2.266{\rm{GPa}},$ $S20=-0.863{\rm{GPa}},$ $S30=-0.839{\rm{GPa}}$ and they are relatively small in comparison with the peak values at the crack tip under the loading. The initial stresses have to be excluded when comparing MD results with LFM solution. Thus, the values from MD, noted (o) in figures 3(a) and (b), are taken as follows: ${S}_{x}=S11-S10,$ ${S}_{y}=S22-S20,$ and ${T}_{MD}={S}_{x}+\mathrm{Re}({\mu }_{1}{\mu }_{2}){S}_{y}={S}_{x}-1.2868{S}_{y}.$ Relation for T_MD follows from equation (1). The agreement between the static atomic stress S_x , S_y and the LFM components σ_x , σ_y is very good in the interval $r/do=\left\langle 1,10\right\rangle .$ Here $r/do$ is a dimensionless distance where $do=a{\rm{o}}\sqrt{2}.$ As for T-stress, in the nearest vicinity of the crack tip, an anisotropic solution ${\rm{T}}={\sigma }_{A}(\mathrm{Re}({\mu }_{1}{\mu }_{2})+{S}_{x}/{S}_{y})$ for the biaxial loading from [6] rather prevails, while the isotropic static solution ${\rm{T}}={f}_{I}(a/W){\sigma }_{A}=-0.548\,{\rm{GPa}}$ prevails with increasing distance r/do from the crack tip toward the right free surface BL. The agreement mentioned above is not trivial; it confirms the self-similar LFM character (no information on boundary correction factors exists in MD simulations) and justifies the concept used in this study.

Zoom In Zoom Out Reset image size

Under the linear loading of 5 PH (2.92 kN s⁻¹) up to a high stress level the crack initiation in MD was monitored already in the vicinity of the lower Griffith stress 1.84 GPa, corresponding to the Griffith stress intensity K_G  = 0.835MPa m^1/2 for plane stress, i.e. at ${\sigma }_{G}={K}_{G}/({F}_{I}\sqrt{\pi a})$ = 1.84 GPa which follows from LFM for cleavage fracture. To investigate the microscopic processes at crack initiation and to avoid discussions of dynamic effects, such as flight time corrections of the loading waves, etc, we use a quasi-static ramp loading presented in figure 13 for temperature of 295 K. In linear region, the loading rate corresponds to 5 PH (2.92 kN s⁻¹), which means a stress rate 1.84 GPa/11650 h for our atomistic 3D specimen. At temperature of 0 K, the loading starts at t = 0 h. For t ≥ 11650 h the applied stress in MD is kept at a constant level SA = 1.84 GPa. Brittle fracture in this run is illustrated in figure 4 by overall shape of the fractured atomistic sample at the end of the simulation. The dependencies RK—time step and RK—COD in figures 5(a) and (b) show that the sample was fractured at a time step ∼ 17650 h (when RK ∼ 0) and when COD corresponds to ∼ 40 Å. The peak value RK_max = 0.336 × 10^–6 N corresponds to a global critical stress in the atomistic sample ${\sigma }_{MD}=R{K}_{\max }/(BW)$ =1.82 GPa, which agrees very well with the independent LFM prediction σ_G = 1.84 GPa for cleavage fracture. This means that the peak value RK_max can be related to fracture toughness of the atomistic sample via the relation ${K}_{MD}={F}_{I}(a/W){\sigma }_{MD}\sqrt{\pi a},$ similar to P_max in fracture experiments. The values K_MD and ${G}_{MD}=C{({K}_{MD})}^{2}$ are K_MD = 0.825 MPa m^1/2 and G_MD = 3.54 J m⁻² using the elastic stiffness coefficient C = 5.198 × 10^–12 m²/N for plane stress. These data are close to Griffith values K_G  = 0.835MPa m^1/2 and 2γ₀₀₁ = 3.612 J m⁻² expected for cleavage fracture. Nevertheless, graphical treatment of MD results shows that this is not a pure cleavage fracture.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The crack was initiated in the middle layer 16 (figure 6(b)) already at time step t_i = 10978 h, while no crack advance is monitored in the first surface layer 1 in the [110] direction due to dislocation emission on $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ oblique slip systems, shown in figure 6(a). The situation in the last surface layer 30 is similar to layer 1. A scheme for this emission is in the right part of figure 6(a). The emission can be thought like a block like shear (BLS) process (see [23, 24]), where the upper plane (011) containing the blue crack tip atom A is rigid and the lower plane (011), containing the red atom B, is the moving plane in the direction of the Burgers vector b oriented in the $[11\bar{1}]$ direction. 3D visualization of the dislocation emission from the edge crack is presented in figure 6(c) using the local numbers of interactions KNT(li) of the individual atoms li from MD with the used short ranged potential [20]. Dislocation cores in figure 6(c) are shown via the red atoms with KNT = 16, in agreement with BLS simulation of dislocation generation in perfect crystals. The atoms lying on free {110} surfaces with KNT = 10 are blue and the atoms lying on the free crack surfaces (001) with KNT = 9 are yellow. Note that the KNT numbers are valid just for the short ranged potentials, where the coordination number at one atom in the perfect bcc lattice is KNT = 14, as in our case. The coordination numbers KNT at the atoms lying on free surfaces are lower due to missing interatomic interactions.

Zoom In Zoom Out Reset image size

As it follows from the geometry of the bcc lattice, the nonzero relative displacements at dislocation emission in figure 6(a) are the components ${\rm{\Delta }}V(A,B)=V(A)-V(B)$ in the y = [001] direction and ${\rm{\Delta }}W(A,B)=W(A)-W(B)$ in z = [110]. The relative displacement of the lower part of the crystal with respect to the upper part is

$$\begin{eqnarray}&&U10={\rm{\Delta }}V(A,B)\cos \,\beta +{\rm{\Delta }}W(A,B)\cos \,\alpha \end{eqnarray} \tag{ 2 }$$

where $\cos \,\beta =1/\sqrt{3}$ and $\cos \,\alpha =\sqrt{2}/\sqrt{3}$ for oblique slip $\left\langle 11\bar{1}\right\rangle \left\{011\right\}.$ In agreement with Rice model [9], we consider the emission to be complete when the relative shear displacement $U10$ in the $[11\bar{1}]$ direction reaches the value $b=ao\sqrt{3}/2,$ i.e. when $U10/b=1,$ which is shown in figure 7 by the horizontal line at the level 1. Figure 7 comes from the continuous monitoring of the displacements of the two atoms A and B in MD at the lower crack face in the first surface layer (110). The time where the horizontal line 1 intersects the $U10$-curve from MD data is the time of the emission t_e = 10836 h. This means that dislocation emission $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ slightly precedes the crack initiation in the middle, i.e. ${t}_{e}\lt {t}_{i}.$ However, figures 6(a) and (c) show that the oblique slip patterns later deviate to horizontal direction, which indicate the change of the slip plane.

Zoom In Zoom Out Reset image size

Recognition of slip patterns in 3D can be performed in small perfect bcc crystal via independent block like shear (BLS) simulation mentioned above. The BLS simulations are described in detail e.g. in [24] and presented here in figure 8. Figure 8(a) is related to the lower oblique pattern at the crack tip A in figure 6(a), coming from the oblique slip plane (011). The slip plane (011) in figure 8(a) intersects our observation plane (110) in the $[\bar{1}1\bar{1}]$ direction under an angle ∼ 35° ($\tan \,\theta =1/\sqrt{2}$) with respect to the horizontal axis of crack extension. After BLS in the b direction $[11\bar{1}]$ it creates the pattern shown in the middle in figure 8(a). This is the same pattern as the lower oblique slip trace in figure 6(a) after dislocation emission on the oblique plane (011). The upper slip patterns in figure 6(a) come from the mirror oblique plane $(0\bar{1}1).$ The oblique dislocation is initially emitted from the surface crack tip as a short curved dislocation that extends gradually under the shear stress on the (011) plane. When such dislocation reaches the straight screw character, parallel with the Burgers vector b = $[11\bar{1}]$ shown in figure 8 (this vector lies both in the (011) and (112) slip planes), the screw dislocation may change the slip plane to (112) via the cross slip. BLS patterns arising in our observation plane (110) after the slip $[11\bar{1}](112)$ are shown in figure 8(b). This slip creates a horizontal trace given by the horizontal intersection line of planes (112) and (110), parallel with the $[\bar{1}10]$ direction of the crack extension. Figure 8(b) illustrates that the slip leads to separation of planes (001), which can help with cleavage along these planes. This is illustrated in figure 9 for the cracked MD sample at a later time 13000 h after the crack initiation and multiple dislocation emissions at the free WL surfaces (110). The multiple surface emissions to oblique $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ slip systems and the subsequent cross slip of the emitted dislocations to oblique slip planes {112} in figures 9(a) and (c) causes at the interior cleavage planes (001) the mentioned separation, which assist the crack extension in the interior of the crystal—see the detail in the middle layer 16 along the crack front in figure 9(b). Again, the scheme in figure 9(b) is valid for the lower slip patterns. The upper horizontal slip patterns come from the mirror oblique plane $(1\bar{1}2).$ The process described above is a new (unpublished) mechanism for a brittle extension of the edge crack (001)[110].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The process is possible also in anisotropic LFM at Griffith level K_G , where K_G is considered a material parameter that does not depend on sample dimensions.

In the oblique slip system $[11\bar{1}](011)$ slip stress τ_b can be obtained by transformation of the stress tensor from the original coordinate system in figure 1: x₁  = $[\bar{1}10],$ x₂  = [001], x₃  = [110] into a new coordinate system: 'x₁  = [011], 'x₂  = $[2\bar{1}1],$ 'x₃  = $[11\bar{1}]$ where the component 'σ₁₃  = τ_b .

$$\begin{eqnarray}&&{\tau }_{b}=1/\sqrt{6}({\sigma }_{12}/\sqrt{2}+{\sigma }_{22}-{\sigma }_{33})\,for\,slip\,[11\bar{1}](011)\,where\,{\sigma }_{ij}={f}_{ij}(\theta ,{\mu }_{1},{\mu }_{2}){K}_{I}/\sqrt{2\pi r}\end{eqnarray} \tag{ 3 }$$

Here θ ∼ 35° as described above, and the other quantities for plane stress conditions at the free surfaces (110) by [35] are: σ₃₃  = 0, μ₁  = +0.6863 + 0.9032i, μ₂ = −0.6863 + 0.9032i, f₁₂  = 0.234, f₂₂  = 1.209, f₃₃  = 0.715. At the nearest vicinity of the crack tip r = b/2 one obtains with K_G  = 0.835MPa m^1/2 for plane stress the value τ_b  = 16.8 GPa by LFM. This value exceeds the stress barrier τ_disl = 14.5 GPa for the slip $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ with the used potential. Thus, the oblique dislocation emission observed in MD is possible also by LFM. Moreover, the LFM value of 16.8 GPa is greater also by independent Frenkel prediction, where the stress barrier is ${{\tau }}_{F}={\mu }b/(2{\pi }\,{d}_{110})=13.8{\rm{G}}{\rm{P}}{\rm{a}},$ using for shear modulus in the 〈111〉 direction the relation $\mu =({C}_{11}+{C}_{22}+{C}_{33})/3.$

As to oblique slip system $[11\bar{1}](112),$ the plane (112) intersects the observation plane (110) under the angle θ = 0 (see figure 8(b)) and the new coordination system is 'x₁  = $[\bar{1}10],$ 'x₂  = [112], 'x₃  = $[11\bar{1}].$ The slip stress 'σ₂₃  = τ_b is:

$$\begin{eqnarray}&&{\tau }_{b}=\sqrt{2}({\sigma }_{22}-{\sigma }_{33})/3\,for\,slip\,[11\bar{1}](112)\,with\,{\sigma }_{22}={K}_{I}/\sqrt{2\pi r}\end{eqnarray} \tag{ 4 }$$

since f₂₂ (θ, μ_1, μ₂ ) = 1 for θ = 0. At the free surface where the cross slip begins, the normal stress component σ₃₃  = 0 and the slip stress are given by a simple relation

$$\begin{eqnarray}&&{\tau }_{b}=\sqrt{2}{\sigma }_{22}/3\,for\,slip\,[11\bar{1}](112)\,with\,{\sigma }_{22}={K}_{I}/\sqrt{2\pi r}\end{eqnarray} \tag{ 5 }$$

In slip system $[11\bar{1}](112)$ the Schmid factor is ∼ 0.47, thus it is larger than in the original slip system $[11\bar{1}](011)$ where Schmid factor is ∼ 0.41. This is a qualitative explanation why the cross slip is realized in MD. Figure 6(a) shows that the cross slip begins at the distance of about r = 8b from the crack tip under the critical level of loading K_I  = K_G  = 0.835MPa m^1/2, which gives by LFM relations (4) and (5): τ_b  = 3.51 GPa. This value is larger than Peierls stress τ_P  ∼ 2 GPa needed for motion of screw dislocation in a perfect bcc Fe crystal at temperature of ∼ 0 K, as presented in a recent comparative study [25]. On the other hand, the LFM relation (5) shows, that for the distance r = b/2 from the crack tip and the same loading K_G  = 0.835MPa m^1/2, the stress concentration in the $[11\bar{1}](112)$ slip system is lower (τ_b  ∼ 14 GPa) than in the slip system $[11\bar{1}](011)$ treated above. Moreover, the value of ∼ 14 GPa is not sufficient for a direct generation of dislocation in the slip system 〈111〉{112} where the stress barrier is τ_disl  = 16.3 GPa [6] with the used potential. In other words, the anisotropy of the stress field described by LFM relations (3), (4) and (5) explains why the surface dislocations are initially generated in the oblique slip systems $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ and later perform the cross slip to oblique system $\left\langle 11\bar{1}\right\rangle \left\{112\right\}.$ Since the atomistic results agree with LFM predictions that cover crack behavior also on macro-scale, realization of the presented micro-mechanism is probable also in macroscopic specimens.

Slower loading of 1 PH (i.e. dP/dt = 0.584 kN s⁻¹ for 1 mm min⁻¹ in experiment) is given in MD by the relation $dP/dt=0.584\,\mathrm{kN}/s={(BW)}_{MD}{(d\sigma /dt)}_{MD}.$ For the short crack we used a stress rate 4.19 GPa/132580 h that obeys the equation above, since the loaded area BW of the atomistic sample is the same as in the case of the longer crack in section 3. 1.1, i.e. (BW)_MD  = 1.8488 × 10^–16 m². Figure 10 illustrates that the loading is sufficiently slow to achieve an agreement with the static LFM prediction at an applied stress level of 1 GPa without any damping in Newtonian equations of motion. In other words, dynamic effects in MD are negligible in this case and figure 10 confirms that the same boundary correction factors can be used for the atomistic 3D sample as predicted by LFM for macro-specimens. The other stress components have a similar character as in figure 3.

Zoom In Zoom Out Reset image size

With increasing linear loading 1 PH, the crack was initiated in the middle of the crystal (layer 16) at the applied stress σ_i  = 3.013 GPa. Corresponding stress intensity is ${K}_{i}={F}_{I}{\sigma }_{i}\sqrt{\pi a}=0.866\,{{\rm{MPam}}}^{1/2}$ with F_I  = 1.68 by [37] for a/W ∼ 0.3 and G_i  = 3.898 J m⁻². Similar to the fast loading 5 PH in section 3.1.1 the crack initiation is accompanied by the surface dislocation emission on the oblique slip systems $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ which occurs at the lower crack face at time step t_e  = 95394 h under applied stress σ_e  = 3. 0205 GPa. The corresponding stress intensity and energy release rate are K_e  = 0.868 MPa m^1/2 and G_e  = 3.916 J m⁻². Unlike 5 PH, by slow loading of 1 PH the prevailing type of plastic processes is oblique twinning generated by the crack itself or by the emitted oblique dislocations $\left\langle 11\bar{1}\right\rangle \left\{011\right\}.$ This is illustrated in figures 11 and 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The atomistic fracture curves RK—t and RK—COD are shown in figures 11(a) and (b). The atomistic sample was fractured at time t_f  ∼ 103300 h (when RK ∼ 0) after a relative large plastic deformation, characterized by a large final COD_f  ∼ 123 Å. Maximum RK corresponds to 0.56 × 10^–6 N, corresponding critical stress is σ_MD  = RK_max/(BW)_MD  = 3.029 GPa and fracture toughness of the atomistic sample is K_MD  = 0.871 MPa m^1/2 with G_MD  = 3.942 J m⁻². These values are greater than in the case of the longer crack under the fast loading 5 PH. Figure 12(a) shows a plastic zone created in front of the crack in our observation plane (110), layer 16. Figure 12(b) is a detail in a perpendicular plane $(\bar{1}10)$ and it clearly shows that the plastic zone is caused by twinning on the oblique planes (112). No significant cross slip of the emitted dislocations to {112} planes was monitored in MD at temperature of ∼ 0 K and loading rate 1 PH. The oblique twins in front of the crack removes the LFM stress concentration at the crack, which hinders the further crack advance. The final character of fracture under slow loading 1 PH is more ductile, which illustrates the larger value of COD_f  ∼ 123 Å in comparison with the faster loading 5 PH, where COD_f  ∼ 40 Å.

To assess generation of the oblique twins in the oblique slip system $[11\bar{1}](112)$ by the crack itself we may use the LFM prediction (5) and consider a distance r = b = 2.4825 Å and K_I = K_MD = 0.871MPa m^1/2 corresponding to RK_max mentioned above. We obtain larger slip stress τ_b than the stress barrier τ_twin for twin generation presented in [24] for the used potential:

$$\begin{eqnarray}&&{\tau }_{b}=(\sqrt{2}/3){\sigma }_{2}=(\sqrt{2}/3){K}_{I}/\sqrt{2\pi r}=10.4{\rm{GPa}},\,{\rm{while}}\,{\tau }_{twin}=\,9.3GPa\end{eqnarray} \tag{ 6 }$$

Estimate of the stress barrier by Frenkel model ${\tau }_{F}=\mu b/(2\pi d)$ gives a reasonable value τ_F = 8.4 GPa, if one considers the distance $d=3{d}_{112}=3ao/\sqrt{6}$ since 3-layers twin is the minimum for its stability [26, 27]. This short discussion shows that oblique twinning on planes {112} in front of the crack is possible also by continuum predictions, which implies its possibility also in experimental macro-specimens. Note that twins can arise also via dissociation of the core both in the case of the screw [28] and edge dislocations [29]—see e.g. Figures 9(a) and (c).

MD simulations at ∼ 0 K also show that the threshold values for emission of the oblique dislocations $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ and oblique twins at the crack front are somewhat different:

$$\begin{eqnarray}&&{K}_{disl\_o\{011\}\left\langle 111\right\rangle }=0.93{K}_{G}(RN)\,\lt \,{K}_{twins\_o\{112\}\left\langle 111\right\rangle }=1.043{K}_{G}(RN)\,at\,\,0\,\,K.\end{eqnarray} \tag{ 7 }$$

The threshold values can be changed by thermal activation in the system which is presented in the next section.

MD simulations with initial temperature of 0 K in section 3.1 showed that the initiation of edge crack in the interior of the 3D crystal was correlated with the surface oblique emission $\left\langle 11\bar{1}\right\rangle \left\{011\right\}.$ Since dislocation emission is a thermally activated process, we start with searching of threshold values for fracture at elevated temperature of 295 K using the ramp loading from figure 13, where the linear phase of loading corresponds to $dP/dt=2.92\,\mathrm{kN}/s$ via the relation ${(BWd\sigma /dt)}_{MD}=2.92\,\mathrm{kN}/s,$ similar to simulations at 0 K. The results are summarized in table 1.

Zoom In Zoom Out Reset image size

When the constant applied stress of 1 GPa was reached during 6332 h, no dislocation emission and no crack initiation were monitored in the time interval up to 30000 h. In the case 1.1 GPa/6965 h, only one transient attempt to emit the dislocation was detected via monitoring of U10/b mentioned above. Similarly, at the level of 1.15 GPa only a transient emission was monitored. Complete emission (only at the surface layer 1) was monitored in the case 1.17 GPa/7408 h at the time step 12128 h, while in last surface layer 30 only generation of the oblique dislocation was observed. As a consequence, only a transient crack initiation was monitored in the middle layer 16. Symmetrical oblique dislocation emission $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ was observed at both surfaces 1 and 30 under loading 1.18 GPa/7471 h at the time 12400 h and later (at 13500 h) also the cross slip to the oblique system $\left\langle 11\bar{1}\right\rangle \left\{112\right\}.$ The crack growth (monitored in the middle layer 16) was slow and no fracture was monitored up to time step 30000 h. The reason is the low applied stress level and beginning of twin formation in the oblique slip systems $\left\langle 11\bar{1}\right\rangle \left\{112\right\}.$ Similar situation was monitored at the loading 1.20 GPa/7598 h—see table 1.

Under a higher loading 1.23 GPa/7788 h a relatively fast crack growth was monitored, accompanied by the oblique emission $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ and subsequent cross slip to slip systems $\left\langle 11\bar{1}\right\rangle \left\{112\right\},$ i.e. the same mechanism of crack growth was monitored both at 0 K and at 295 K. This is illustrated by figure 14 for time 13500 h. When the horizontal slip traces from screw dislocations $\left\langle 11\bar{1}\right\rangle \left\{112\right\}$ cross the crack plane (001) on inner layers (110) (e.g. the layers 6, 16, etc) they prepare a pre-processing zone for easy crack extension via separation of the planes (001)—see figure 14(b) where the position of the original crack tip point is denoted via the arrow. Cleavage crack extension in middle layer 16 at later time is illustrated in figure 15(a). However, generation of the oblique twinning on planes {112} was observed in front of the crack before the final fracture (figure 15(b)), unlike corresponding simulations at ∼ 0 K in section 3.1.1. When the twins reach the right free surface BL, new stress concentrators arise at free surface BL and new oblique dislocations $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ are emitted from the free surface. This increases crack opening and the final phase of fracture is realized via a small plastic necking of the atomistic sample. The final phase of fracture might be influenced by the finite sample dimensions. However, it did not occur at 0 K in figure 4 under the loading 5 PH and it may not occur even at 295 K, as mentioned below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Oblique twins occurred in all runs 5 PH we tested at elevated temperature (see table 1), but randomly at different times, which affected the velocity and character of crack extension. The results are summarized in table 1. The largest portion of cleavage fracture was observed at the loading 1.59 GPa/10067 h and at 1.23 GPa/7788 h, which illustrates a 3D visualization of the fracture surfaces in figure 16.

Zoom In Zoom Out Reset image size

Table 1 is also used to assess thermal activation in the system. As to surface dislocation emission $\left\langle 11\bar{1}\right\rangle \left\{011\right\},$ if we consider the transient emission under σ_e = 1.1 GPa with G_e =1.286 J m⁻² = G(T) as a threshold, then we may estimate the degree m of thermal active tion from the relation:

$G\left(T\right)=G\left(0K\right)-m{k}_{B}T/{A}_{110}$ for $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ slip system (8) by analogy with [8]. For T = 295 k and with the area per one atom ${A}_{110}={\left(a{\rm{o}}\right)}^{2}\sqrt{2}/2$ and $G(0K)=3.135\,{J/m}^{2}$ for the emission at 0 K (see table 1) we obtain the degree of thermal activation m = 26.43. This value does not differ significantly from the prediction of $30{k}_{B}T$ by Rice and Beltz [11] for thermal activation of dislocation emission in bulk crystal of a metal (see the fluctuation of m-values in table 2 in [8]). If we consider the emission from both surfaces 1 and 30 under σ_e  = 1.18 GPa with ${G}_{e}=1.4656\,{J/m}^{2}=G(T),$ then a lower value of m = 16.65 is obtained. The differences between the surface and bulk emissions are treated e.g. in [30, 31].

Important factor for brittle or ductile behavior of the crack (001)[110] is the threshold for generation of the oblique twins. This generation can start as an emission of partial dislocations $\left\langle 11\bar{1}\right\rangle \left\{112\right\}$ with the Burgers vector b/3 and thus, we may consider thermal activation of ∼ $30{k}_{B}T$ for bulk crystal from [11]. Using K_twin from relation (7), then ${G}_{twin}(0K)=C{({K}_{twin})}^{2}=G0.$ Using the relation $G(T)=G0-30{k}_{B}T/{A}_{112},$ we may derive the value of G(T), K_twin (T) and subsequently σ_twin  = 1.32 GPa for T = 295 K. The area per 1 atom on plane {112} is ${A}_{112}=a{o}^{2}\sqrt{3}/\sqrt{2}.$ Close to crack initiation we observed the oblique twins at the crack tip on the inner layer 16 under the ramp loading 1.30 GPa/8231 h at the time step 13500 h (see table 1). This caused that no fracture was monitored up to time step 50000 h. Note that a small change in loading rate to dP/dt = 3 kN s⁻¹ improved the situation and the overall shape of the atomistic sample at the final fracture was similar to the brittle fracture in figure 4 at 0 K. Nevertheless, the inspection of fracture surfaces showed that more cleavage facets occurred in MD simulations using the standard loading rate 2.92kN s⁻¹ (5 PH), namely the ramp loading 1.23 GPa/7788 h–see figure 16.

Table 1 shows also other differences between the simulations at 0 K and 295 K. As mentioned in section 3.1.1, the maximum of the resulting force RK in one half of the crystal physicaly represents a global resistence of the atomistic system against fracture and it occurs on the time scale in each case before the dislocation emission and crack initiation. Consequently, σ_MD represents an internal critical stress available in MD at elevated temperature, K_MD is global fracture toughness in MD at elevated temperature and G_MD is the corresponding energy release rate in the heated atomistic sample. All the MD quantities are larger than σ_i (or σ_e ), K_i and G_i for crack initiation, derived from external loading, due to thermal activation. The average value of G_MD at 295 K is 5.35 J m⁻² and the difference ${G}_{MD}\left(295\,\,K\right)-{G}_{MD}\left(0\,K\right)=1.81\,{J/m}^{2}$ is close to an average from continuum predictions of $30\,{k}_{B}\left(295{\rm{K}}\right)/\,{A}_{110}=2.1\,{J/m}^{2}$ and $30{k}_{B}\left(295{\rm{K}}\right)/\,{A}_{112}=1.2\,{J/m}^{2}$ for thermal activation of the oblique slip system $\lt 11\bar{1}\gt \{011\}$ and $\lt 11\bar{1}\gt \{112\}.$ This is because the brittle fracture in MD at 295 K and 0 K is associated with the mentioned slip processes.

During linear loading in time we found out that the shorter edge crack with relative ratio a/W ∼ 0.3 was initiated at the external stress of 2.35 GPa and thus for the ramp loading 1 PH we used the stress loading 2.35 GPa/74400 h by the scheme in figure 13, where the linear phase corresponds to $dP/dt=BW(2.35{\rm{GPa}}/74400\,h)=0.584\,\mathrm{kN}/s.$ Similar to the fast loading 5 PH, the crack initiation was correlated with the emission of dislocations to oblique slip systems $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$—see table 2. Unlike 5 PH, only one-sided cross slip to slip system $\left\langle 11\bar{1}\right\rangle \left\{112\right\}$ was observed. Soon after crack initiation the oblique twins were generated—see figure 17. As a consequence, ductile fracture was monitored in MD, which illustrates figure 18 with the deformed fracture surface after the oblique twining. Such a sloping fracture surface is presented in figure 7 in [4], coming from SEM analysis of the fractured specimen No10 with the edge crack a/W ∼ 0.3 loaded slowly in mode I.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that the oblique twinning and oblique dislocation emission $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ at the edge crack (001)[110] were already briefly reported in [32], where 3D simulation of the fracture experiment [4] on the shorter SEN specimens is presented. Oblique emission $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ was also reported in [23] under cyclic loading in mode I and mode II, where it resulted in crack deflections. No cross slip of the screw dislocations $\left\langle 11\bar{1}\right\rangle \left\{011\right\}$ and their contribution to crack extension is reported in [23, 32].

Table 2 also includes comparison of the slow ramp loading 1 PH and the fast ramp loading 5 PH. The results presented in table 2 and the comparison of fracture surfaces in figure 18 for 1 PH and the fracture surfaces in figure 16 for 5 PH show that the slow loading 1 PH supports more oblique twinning, which lead to a ductile fracture. To verify that, we also performed MD simulations under linear loading in time, briefly presented in the next section.

As mentioned above, one disadvantage of the linear loading in time is that MD results for dislocation emission and crack initiation occurring at x-axis need to be corrected at least by flight time ${\rm{\Delta }}{t}_{x}=(L/2)/{C}_{L[001]}$ of the fastest longitudinal loading waves (C_{L [001]} = 5550 m s⁻¹ [14]) from the loaded borders to the x-axis (where they are generated). In our case Δt_x = 1162 h, which is not negligible namely under faster loading 5 PH. The second disadvantage is that the high stress level before fracture makes the small 3D atomistic samples more ductile than expected in relation to experiments. Nevertheless, the initial crack behavior in the vicinity of the initiation can be described qualitatively well in relation to the experiments.

Comparison of the corrected MD results at temperature of 295 K with linear loading 1 PH and 5 PH is presented in table 3. Under linear loading the values COD_f and W_f /BW are larger than those in table 2 (ramp loading) because of the larger external forces acting during linear loading. This supports more oblique twining in the atomistic samples, which leads to more ductile behavior in comparison with ramp loading. Nevertheless, the values K_i , COD_f and W_f /BW presented in table 3 show that the slow loading 1 PH leads to more ductile behavior in comparison with 5 PH, which is in a qualitative agreement with fracture experiments [3, 4]. In this study it is illustrated in figures 19(a) and (b) for the linear loading 1 PH and 5 PH shortly after the crack initiation. Figure 20 shows the atomistic dependencies RK—t and RK—COD, in analogy with the experimental dependencies P—t and P—ΔL presented in [3, 4] and in other experimental studies, reviewed e.g. in [2].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As already mentioned in Introduction, one may expect only a qualitative agreement between MD and the experimental results, not quantitative. And really, in agreement with the experimental results on bcc metals reviewed in [2], our MD curves RK—t and RK—COD with the increasing quasi-static loading rate are displaced to the left and towards the lower peak values, similar to experiments, i.e. with the increasing loading rate the fracture is faster and less ductile, or more brittle respectively. The same is valid for the ramp loading 1 PH and 5 PH that is more informative to answer our question on the contribution of the crack itself to ductile versus brittle behavior under different loading rates. The faster loading rate 5 PH leads to more brittle fracture, while the slower loading rate 1 PH to more ductile fracture, which was confirmed via inspection of the fracture surfaces in MD under the loading 1 PH (figures 18) and 5 PH (figure 16). The presented MD results comply with our experiments [3, 4] that motivated this study. The grape shape plastic zone in figure 4 from the experiment [3] and in figure 2 from [4] at the short crack a/W ∼ 0.3 under the loading rate 1 mm min⁻¹ is similar to MD figures 12(a) and 19(a) showing the oblique twinning in front the short crack a/W ∼ 0.3 under the slow loading 1 PH. Under the fast loading 5 PH (for 5 mm min⁻¹) such a plastic zone does arise neither in MD nor in the experiment. The average velocity of the crack growth under the ramp loading 5 PH is Vcr = 417 m s⁻¹, as follows from table 2. If we suppose such a crack speed in the experiment [3] with the longer crack a/W = 0.422 (where W ∼ 10 mm and loading rate was 5 mm min⁻¹), then the time needed for fracture will corresponds to about 1.4 × 10^–5 s, which is too short time for monitoring the crack extension via classical optical microscope and the process seems like a sudden fracture, as observed in [3].

As discussed already in [33], qualitative comparison of MD results in pure bcc iron with the experimental results on Fe-Si single crystals is justified in the region of low Si contents (up to about 6wt%Si), where the elastic constants are practically the same as in bcc iron and thus, where the same stress and energy barriers for a given crack tip process are expected by continuum predictions, valid both on macro and atomistic level.

Slip patterns from planes {112}, created by gliding dislocations and twinning nearby and on the fracture surfaces of Fe-3wt%Si crystal No2 (with K_Ic = 24MPa m^1/2) were reported in [15], where the post-fracture SEM and TEM analysis of the experiment [3] was done. The specimen No2 with the short edge crack a/W ∼ 0.3 was loaded with the cross head speed 1 mm min⁻¹, which in our MD study correspond to the case 1 PH. While the fracture surface from specimen No 2 was of ductile character, under higher loading rate (5 PH in this study) prevailing cleavage character of fracture is mentioned in [15]. It complies with presented MD results for 1 PH and 5 PH at temperature of 295 K that show that at least a part of the slip processes revealed in [15] can be generated by the crack itself.

In agreement with this MD study, the activity of the slip systems on the oblique planes {011} was recently reported in [16], where tension fracture experiments with a central crack (001)[110] in thin single crystals Fe-3wt%Si were analyzed in detail via electron channeling, electron back scattering and scanning electron methods (SEM). Main contribution to the observed striations in [16] is attributed to slip patterns from the inclined slip systems $\left\langle 11\bar{1}\right\rangle \left\{112\right\}$ due to a larger Schmid-factor (∼0.47) for this system in comparison with the oblique systems 〈111〉{011}, where it is ∼ 0.41. Generation of blunting dislocations $\left\langle 11\bar{1}\right\rangle \left\{112\right\}$ by the edge crack (001)[110] was neither observed in our MD simulations nor preferred by LFM analysis. This so-called non-Schmid behavior [34] can be explained in our case by the anisotropy of the near stress field at the crack (001)[110] presented in detail in [35]. (As mentioned in section 3.1.1, at Griffith level of loading and at the distance r = b/2 from the crack tip we obtain from equation (3) the slip stress τ_b  = 16.8 GPa for the slip system 〈111〉{011}, while for the inclined slip systems $\left\langle 11\bar{1}\right\rangle \left\{112\right\}$ according to [35]: τ_b  = 7.59 GPa for plane strain and τ_b  = 5.87 GPa for plane stress conditions.) It is not in contradiction with the conclusions from [16], since the slip patterns on the inclined systems $\left\langle 11\bar{1}\right\rangle \left\{112\right\}$ from [16] may concern pre-existing dislocations or easy twinning in the slip systems (see fig.4.14 in [1] or ref. [24, 26]). Note that the emission of blunting dislocations $\left\langle 11\bar{1}\right\rangle \left\{112\right\}$ is preferred by MD and by anisotropic LFM for the crack orientation $(\bar{1}10)[110],$ where the same slip systems are oriented in hard, anti-twinning direction—see [6–8] and [35]. This was confirmed experimentally e.g. in [4] and [33].

Recent MD simulations in fcc Nickel [36] show how fracture toughness may depend on the thickness of the 3D atomistic sample. As to bcc Fe and the edge crack (001)[110], we believe that the mentioned oblique slip processes from MD, generated by the crack itself at the free (110) surfaces (WL in figure 1), can also occur in thicker specimens of our orientation. This is possible according to the independent LFM analysis presented above and it complies with the experimental observations in [16].