**E**ff**ect of Cutting Crystal Directions on Micro-Defect Evolution of Single Crystal** γ**-TiAl Alloy with Molecular Dynamics Simulation**

**Jianhua Li 1,2, Ruicheng Feng 1,2,3,\*, Haiyang Qiao 1,2,\*, Haiyan Li 1,2,3, Maomao Wang 1, Yongnian Qi <sup>1</sup> and Chunli Lei 1,2**


Received: 30 October 2019; Accepted: 25 November 2019; Published: 28 November 2019

**Abstract:** In this work, the distribution and evolution of micro-defect in single crystal γ-TiAl alloy during nanometer cutting is studied by means of molecular dynamics simulation. Nanometer cutting is performed along two typical crystal directions: 100 and 101 . A machined surface, system potential energy, amorphous layer, lattice deformation and the formation mechanism of chip are discussed. The results indicate that the intrinsic stacking fault, dislocation loop and atomic cluster are generated below the machined surface along the cutting crystal directions. In particular, the Stacking Fault Tetrahedron (SFT) is generated inside the workpiece when the cutting crystal direction is along 100 . However, a "V"-shape dislocation loop is formed in the workpiece along 101 . Furthermore, atomic distribution of the machined surface indicates that the surface quality along 100 is better than that along 101 . In a certain range, the thickness of the amorphous layer increases gradually with the rise of cutting force during nanometric cutting process.

**Keywords:** molecular dynamics; nano-cutting; crystal direction; γ-TiAl alloy; stacking fault

### **1. Introduction**

With the development of modern industry, nano-processing technology has been widely used in aerospace, biomedical, national defense and military and other high-tech fields [1]. In the fabrication of nanodevices, the basic understanding of the material removal mechanism and the evolution of defects is becoming more and more important. Nanofabrication involves material deformation of only a few atomic layers of about 2–10 nm, and the removal of atoms is discrete and discontinuous. Obviously, experimental and continuous theory is inaccurate. Since the 1990s, Molecular Dynamics (MD) simulation technology has been successfully applied to study the removal mechanism of various materials [2]. The MD can reveal the phenomena that cannot be observed by the traditional method, and capture the structure of material, position and velocity of atoms.

Over past decades, many scholars have done a lot of work by MD method. These studies were mainly focused on monocrystalline silicon, single crystal copper, silicon carbide and monocrystalline nickel that are described by Embedded Atom Method (EAM), Morse and Tersoff potentials. Wang et al. [3] studied the formation mechanism of stress-induced SFT and results show that the complex SFT is nucleated due to tensile and compressive stress. Chuvashia et al. [4] analyzed the mechanism of chip formation under high temperature conditions, and the key conclusion was that less

heat is released and larger chips are formed when cutting at higher temperatures or on (111) crystal faces. Besides, the cutting force is reduced by 24% when the temperature at 1173 K compared with room temperature. Xie et al. [5] found that, as cutting-edge radius increases, both the dislocation slip zone and chip formation zone are suppressed while the elastic deformation zone tends to continually grow in monocrystalline copper nano-cutting process by MD simulations. Dai et al. [6] discussed the effects of different tool structures on nanometer cutting of single crystal silicon and found that some structural tools can cause a lower temperature. Fung et al. [7] researched the effect of surface defects of tool on tool wear during nano-machining and the statistical results showed that atoms loss in the defective region of tool is one order of magnitude greater than that in the non-defective region. Ren et al. [8] analyzed the grinding process of monocrystalline nickel at various speeds and depths with the molecular dynamics simulation, and revealed that, owing to the effect of the elastic deformation, a part of the defects can be eliminated after grinding and the dislocations and the variation in the number of HCP atoms cause the change of grinding force. Xu et al. [9] investigated the effect of hard particles on surface generation in nanometer cutting. Liu et al. [10] expounded the evolution of SFT and the work hardening effect in monocrystalline copper. Xu et al. [11] explained the effects of recovery on surface generation in nanometric cutting. Zhu et al. [12] simulated monocrystalline Nickel nano-cutting and found that a stacking fault is main cause of work hardening. Wang et al. [13] found that the tangential and normal forces of the cutting are different in different lattice orientations, and the stacking fault induced by cutting propagates in the basal plane and depends on the angle between the cutting direction and the c-axis through the MD simulations of four off-axis 4H-SiC nano-cutting processes.

The above literature focused on revealing the mechanism of material removal, microdefect evolution, the correlation between chip removal and cutting speed or temperature and cutting force change. However, few studies have focused on the effect of cutting crystal orientation on the evolution of micro-defect in the materials during nano-cutting, and, in particular, the evolution of a micro-defect along different cutting crystal orientation of γ-TiAl alloy has not been found. Additionally, γ-TiAl alloy has the advantage of low density, good high temperature strength and oxidation resistance, so it is widely used in aeroengine and automobile industries [14,15]. γ-TiAl alloy has face-centered tetragonal (FCT) with an L10 structure [16], as shown in Figure 1. Ti atoms and Al atoms are alternately arranged along the [001] direction. Furthermore, the crystal axis ratio of γ-TiAl is c/a = 1.02, which reduces the crystal symmetry and reduces the slip system. Therefore, the mechanical properties are different when cutting in different crystal directions. The purpose of this paper is to research the evolution of micro-defect by cutting along different crystal orientations of γ-TiAl alloy using MD methods. Section 2 will introduce the cutting model and methods. Section 3 will discuss the results of simulation and analyze the evolution of the micro-defect. Finally, the results will be summarized in Section 4.

**Figure 1.** L10 structure of γ-TiAl.

### **2. Model and Methods**

### *2.1. The MD Simulation Model*

As shown in Figure 1, the *x*, *y*, and *z* coordinate axes correspond to the orientations [100], [010] and [001], respectively. The lattice constants are a0 = b0 = 0.4001 nm and c0 = 0.4181 nm, which are consistent with the experimental values a0 = b0 = 0.4005 nm and c0 = 0.40707 nm [17]. The nano-cutting model is shown in Figure 2, consisting of a monocrystalline γ-TiAl workpiece and a diamond tool. The workpiece contains 131,175 atoms and has a size of 20 <sup>×</sup> 10 <sup>×</sup> 11 nm3 along the [100], [010] and [001], respectively, and the diamond tool contains 7793 atoms. A periodic boundary condition along the *z* direction is adopted to eliminate the influence of size effect. In addition, five layers of the right and the bottom edges are fixed, and five layers next to them are a thermostat region which surrounds the Newton atoms. The thermostat region serves as a heat reservoir at a constant temperature 293 K for Newton atoms [18]. Newton atoms are used to study all aspects of the cutting process, the motion of Newton atoms obeys Newton's second law. Moreover, the selection of cutting speed is based on references [19,20], which combines with the characteristics of the γ-TiAl alloy.

**Figure 2.** The MD simulation model of nano-cutting.

### *2.2. Interatomic Potential*

The reliability of any MD simulation greatly depends on the quality of the interatomic potential employed. For metallic systems, the EAM is a widely used technique [21], which is fitted from a large database of both experiment and abinitio data. The total potential is expressed by

$$E = \sum\_{i} F\_{i}(\rho\_{i}) + \frac{1}{2} \sum\_{ij(i \neq j)} \rho\_{ij}(r\_{ij}) \tag{1}$$

$$\rho\_{\bar{i}} = \sum\_{j(j\neq \bar{i})} f\_{\bar{j}}(r\_{\bar{i}\bar{j}}) \tag{2}$$

where *Fi* is the embedding energy of atom *i*, whose electron density is ρ*i*, φ*ij* is the relative potential energy of atom *i* and atom *j*, *rij* represents the distance between atom *i* and atom *j*, *Fi*(ρ*i*) is the sum of electron cloud density produced by all other atoms, extranuclear electron to atom *i*, and *fj*(*rij*) is the electron density produced by atom *j*.

The interaction between Al-C and Ti-C is described by Morse potential, and the expression is given by

$$V(r\_{ij}) = D\_c \left[ e^{2\alpha(r\_\ell - r\_{ij})} - 2e^{\alpha(r\_\ell - r\_{ij})} \right] \tag{3}$$

where *De* is the cohesive energy between atoms *i* and *j*, α is the elastic modulus of the material, *rij* represents the instantaneous distance between atoms *i* and *j*, and *re* is the equilibrium distance between atoms *i* and *j*. The parameter of Morse potential can be obtained from Zhu et al. [22]. As *DTi-C* = 0.982 eV, α*Ti-C* = 22.83 nm<sup>−</sup>1, *rTi-C* = 0.1892 nm, *DAl-C* = 0.28 eV, α*Al-C* = 27.8 nm<sup>−</sup>1, *rAl-C* = 0.22 nm.

As the rigidity of diamond is much higher than that of γ-TiAl, the tool is modeled as a rigid body. Table 1 shows the MD simulation parameters.


**Table 1.** The MD simulation parameters.

In this paper, Common Neighbor Analysis (CNA) [23] is utilized to identify the defect type, such as atomic cluster, dislocation, stacking fault etc.

### **3. Results and Discussion**

### *3.1. Material Removal Mechanism in Nanometric Cutting Process*

In the nano-cutting process, substrate atoms adjacent to the tool tip are subjected to the high compressive energy, resulting in highly displaced and disordered atoms. Figure 3 shows the diagram of material deformation; (a) and (b) represent the cutting crystal direction along 101 and 100 , respectively. As shown in Figure 3, in the first place, the lattice orientation is 45◦ to the tool when the cutting crystal direction is along 101 . However, the lattice orientation is 0◦ or 180◦ to the tool along 100 . In the second instance, the shape of the chip is also different: The chip along 100 is higher than that along 101 , and the atomic packing is more compact. During the machining process, the elastic and plastic deformation of γ-TiAl occurs with the increase of cutting distance, which is consistent with the literature [24]. The atoms of workpieces are distributed symmetrically on both sides of the tool to form the side-flow, and the others pile up in front of the tool to form chips. The chips are formed by extrusion force from the tool, and is removed in the form of atomic group. Then, the machined surface is generated with a certain elastic recovery.

**Figure 3.** The diagram of cutting distances at 12 nm and the cutting crystal directions along 101 , 100 respectively, (green represents Ti and red represents Al). (**a**) Cutting crystal direction along 101 . (**b**) Cutting crystal direction along 100 .

### *3.2. System Potential Energy and Lattice Deformation in Nanometric Cutting Process*

The system potential energy and lattice deformation can be used as indicators of material deformation during the cutting process. Figure 4 shows the effect of different cutting crystal directions on the potential energy of the system, and the two potential energy curves rise with the increase of the cutting distance and the trend of the grows are the same. However, the potential energy of the system is higher than that along 100 when the cutting direction is along 101 . On one hand, the different lattice arrangements result in different surface energy of the workpiece. On the other hand, the lattice bond breaks more and the internal energy releases more when the cutting direction is along 101 . However, the lattice bond breaks less and the internal energy releases less during nanometric cutting along 100 .

**Figure 4.** Change curve of system potential energy in the nanometric cutting process.

The CNA method is used to identify the crystal structure transformation during nano-cutting process. Figure 5 shows the effect of cutting directions on atomic structure transition. Figure 5a,b represents the number of structural transitions of HCP (Hexagonal Close-Pack) and BCC (Body-Centered Cubic), respectively. Two bar graphs are compared between the number of HCP transition crystal structure and the number of BCC transition along 101 . In the initial stage of cutting, it increases with the increase of cutting distance. The value fluctuates greatly and shows the opposite trend when the cutting distance is 8–16 nm. This is due to the complex stacking faults nucleation, expansion, annihilation. There is a large area of stacking fault at 8–16 nm, which increases the number of HCP transformations and inhibits the number of BCC. The HCP structure plays a dominant role in the whole cutting process. Comparing the number of HCP and BCC transitions along 100 , we can see that the BCC crystal structure predominates when the cutting distance is less than 6 nm and more than 14 nm, while the HCP crystal structure dominates and the value fluctuates sharply between 8–14 nm. As the cutting force increases with the increase of cutting distance, which leads to the formation, expansion and annihilation of SFT, thus increasing the number of HCP and decreasing the number of BCC. This indicates that there is a competitive relationship between the transformation of HCP and BCC crystal structures in nano-cutting. It can be seen from Figure 5a that the number of transitions along 101 is larger than that along 100 , while Figure 5b is just the opposite. As the crystal orientation of the workpiece is 45◦ with the tool along 101 , it is easier to form intrinsic stacking faults; therefore, the number of HCP is greater along 101 .

**Figure 5.** Graph of atomic structure variations with the cutting distance on the surface.

### *3.3. Evolution of Dislocation in Di*ff*erent Cutting Crystal Direction*

In order to clearly understand the evolution process of surface defects during nano-cutting process, Figure 6 shows the effect of different cutting crystal directions on the cutting surface. Figure 6a,b represents dislocation nucleation and dislocation slip along 100 and the corresponding time is 124.5 ps and 125.5 ps, respectively. Figure 6c,d represents dislocation nucleation and dislocation slip along 101 and the corresponding time is 98 ps and 100 ps, respectively. Among them, white, red and green represent surface atoms (dislocation atoms), HCP atoms and FCC atoms, respectively. The nucleation and slip of dislocations are marked in Figure 6. In the early stage of nano-machining, it can be seen clearly the dislocation nucleated, slide and emitted at 45◦ to the cutting tool, which is explained by the Schmid rule; the resolved shear stress reaches the maximum value at ±45◦ in front of the tool and slips in these directions. The dislocation nucleation along 100 is late; however, its dislocation slip region is larger than that along 101 , as a result of the formation of SFT along 100 , which is a reaction between dislocations. Finally, the SFT annihilates to the boundary of the workpiece. Due to the certain image force to dislocation from free surface, the slip area of dislocation increases gradually with the increase of cutting distance, as shown in Figure 6. This is owing to the rise of cutting distance, cutting force and temperature, which contributes to the slip of dislocations during the nanometric cutting process.

**Figure 6.** Variation of cutting surface in nanometric cutting process. (**a**) Dislocation nucleation at 124.5 ps along 100 . (**b**) Dislocation slip at 125.5 ps along 100 . (**c**) Dislocation nucleation at 98 ps along 101 . (**d**) Dislocation slip at 100 ps along 101 .

*Metals* **2019**, *9*, 1278

The evolution of internal defects in the nanometric cutting process is a dynamic process. Figure 7 shows the internal defects evolution in the workpiece along 101 ; atoms are colored according to the CNA. For a better observation of the evolution of defects, ordinary FCC atoms are removed. The white, red and blue represents Other, HCP and BCC atoms, respectively. Two layers of red atom represent an intrinsic stacking fault [25]. In the early stage of the nano-machining, the stacking fault is generated along - 011 crystal planes for the instability of the cutting, as shown in Figure 7a. The stacking fault nucleates and extends at shear-slip zone below the cutting tool, and finally annihilates at free surface of workpiece, as depicted in Figure 7b. Two Shockley partial dislocations extend along (011) and - 011 planes and form a Lomer-Cottrell lock [26], which finally forms a "V"-shaped dislocation loop, as illustrated in Figure 7b. In the middle stage of nanomachining, it can be seen clearly that the complex dislocation loop and the stacking fault are generated beneath the tool attributing to the rise of cutting force and temperature, as shown in Figure 7c. With the increase of cutting distance, the stacking fault and the complex dislocation loop disappear gradually and are removed in the form of chips, as depicted in Figure 7d. However, the defects continue to evolve in the workpiece, and a small number of the atomic clusters are generated inside the workpiece, as illustrated in Figure 7d. The atomic clusters are formed as a result of the combined effects of temperature and plastic deformation during nano-cutting.

**Figure 7.** Internal defect evolution of workpiece along 101 . (**a**) Stacking fault along - 011 . (**b**) "V"-shaped dislocation loop. (**c**) Dislocation loop and stacking fault. (**d**) Atomic clusters.

Figure 8 shows the internal defect evolution inside the workpiece along 100 . It can be seen that the stacking fault, the SFT, the dislocation loop and the atomic cluster are formed in the subsurface. The stacking fault in the shear-slip zone extends along the - 1 11 and - 111 planes, as shown in Figure 8a,b. After processing, there are two stacking faults slipping along each slip direction, which finally annihilate at free surface of the workpiece. As can be seen from Figure 8c,d, the complex SFT consists of three edges, three stacking fault planes and one workpiece surface, which are formed in the subsurface of the workpiece. The order of the nucleation of the intrinsic stacking fault is marked "S1", "S2", "S3", where 1/3<001> Hirth dislocation is formed between "S1" and "S2", 1/6<112> Shockley partial dislocation is produced between "S2" and "S3", and 1/6<110> Stair-rod dislocation is generated between "S1" and "S3" along the respective slip surfaces, in which Stair-rod and Hirth are immobile dislocations to preserve the shape of the SFT. As a result of the instability of external forces during the cutting process, the SFT is eventually removed in the form of chips, and a small number of atomic clusters and dislocation loops are formed, as illustrated in Figure 8e.

**Figure 8.** Internal defect evolution of the workpiece along 100 . (**a**) Stacking fault along - 1 11 . (**b**) Stacking fault along - 111 . (**c**) SFT. (**d**) Nucleation of the intrinsic stacking fault. (**e**) Atomic clusters.

### *3.4. E*ff*ect of Cutting Crystal Direction on Cutting Force, Temperature and Cutting Surface Quality*

The cutting force is an indispensable factor in the analysis of the nanometer cutting process. Figure 9 shows tangential force, *Fx* and *Fy*, in the cases with different cutting crystal directions. It is noticeable that when the tool and the workpiece are close to a certain distance, the workpiece atoms give a certain repulsive force to the tool, where the cutting force has a negative value. Subsequently, the *Fx* increases rapidly with cutting distance and becomes gradually stable. This conforms with reference [27]. The tangential force is basically stable and fluctuates within a certain range when the cutting distance exceeds 6 nm, and the fluctuation of the curve is caused by crystal structure transformation, stacking fault and lattice reconstruction. In addition, there is a consistent trend in the curve of *Fx* for the cutting crystal directions. In Figure 10, the temperature of the cutting zone grows along with increase of cutting distance before the distance reaches 15 nm. Comparing derivation of the temperature and the cutting force, the temperature of the workpiece grows along with the increase of cutting force at the early stage of cutting. During the final of the cutting process, the temperature goes down, which happens as the cutting force drops. The temperature of workpiece with cutting direction of 100 is higher when the cutting distance is between 12 nm and 16 nm than that of the one along the 101 direction. This phenomenon results in the resistance of materials varying along different directions and more heat being generated during the cutting of a workpiece along the 100 direction.

**Figure 10.** Temperature of cutting zone.

In order to clearly distinguish the size of the cutting force, this article details statistics of the average cutting force. Besides, the cutting distance is selected between 6–18 nm, which aims to reduce error of the cutting force. The average cutting force is shown in Table 2. It can be clearly seen that the average cutting force along 100 is greater than that along 101 . This results from the interaction of dislocations along 100 , which activates more dislocation slip systems and intersects in the cutting direction, and, in consequence, the greater cutting force is required for material removal than that along 101 . This is in line with reference [28]. When the cutting crystal direction changes from 101 to 100 , the tangential force and normal force is increased by 13.32%, 10.25%, respectively.


The quality of the cutting surface is crucial during nanometer cutting process. Figure 11 shows the atomic distribution of machined surface under different cutting crystals. Figure 11a,b represents the cutting direction along 101 and along 100 , respectively. The blank area in the figure means the atom lost in cutting, and there is a large atom loss area. The number of atoms remaining on the machined surfaces in Figure 11a,b is 706 and 697, respectively, and these two numbers are close. In order to analyze the surface quality more accurately, the *y*-direction distribution of atoms on the top layer of the machined surface is calculated. As shown in Figure 11, the black and red curves represent the *y*-direction distribution of atoms along 101 and 100 , respectively. It can be seen that the atom distribution span along 101 is 0.01–0.1 nm, with the largest number of atoms at 0.05 nm, and its peak value is 145. However, the atom distribution span along 100 is 0.01–0.1 nm, and the peak value is 211 at 0.06 nm. From the analysis of data and curve distribution, it can be seen that the atomic spans of the two cutting directions are the same. Furthermore, the *y*-direction atom distribution along 100 is more concentrated and compact, which shows that the machined surface quality along 100 is better than that along 101 .

**Figure 11.** Atomic distribution of machined surface under different cutting crystals.

In order to analyze the effect of different cutting directions on the quality of the cutting surface, the atomic distribution in the First layer along the *y*-direction was quantified, that is, the unevenness of the First layer atomic plane. Figure 12 shows the effect of different cutting directions on atomic distribution. It can be seen that both curves have multiple peaks; the maximum number of atoms is 135, 87 at 0.15 nm, 0.145 nm. Secondly, the distribution of atoms along 100 is concentrated, while the distribution of atoms along 101 is dispersed. This indicates that the distance between atoms is smaller and the planeness is better than that along 101 . When the cutting direction is along 101 , the orientation of the lattice is 45◦ to the tool, resulting in smaller fluctuations of atoms in the nano-cutting process.

**Figure 12.** *Y*-direction atom distribution in the top layer of machined surface under different cutting crystals.

As shown in Figure 13, the effect of different cutting directions on the amorphous layer is displayed. The amorphous layer is composed of dislocation and disordered atoms, which is formed below the machined surface. Figure 13a,b represents the cutting direction along 101 and 100 , respectively. It is clear that the thickness of amorphous layer along 100 is slightly larger than that along 101 . The cutting force along 100 is slightly higher than that along 101 . This, in a certain range, the deeper the thickness of the amorphous layer, the greater the cutting force during the nanometric cutting process.

**Figure 13.** Effect of cutting direction on thickness of amorphous layer. (**a**) Thickness of amorphous layer along 101 . (**b**) Thickness of amorphous layer along 100 .

### **4. Conclusions**

In this paper, the nanometric cutting process of single crystal γ-TiAl alloy has been studied by molecular dynamics simulation. The effect of cutting directions on the material deformation, system potential energy, lattice deformation, machined surface, cutting force and the evolution of dislocation are thoroughly discussed. The conclusions are as follows:

(1) The orientation of the lattice is 45◦ to the tool when the cutting crystal direction is along 101 . However, the arrangement of the lattice is 0◦ or 180◦ to the tool when the cutting direction is along 100 , and the chips are removed as atomic groups. The potential energy of cutting orientation along 101 is larger than that along 100 .

(2) The evolution of defects is affected by cutting directions in the workpiece. Dislocation slip starts at the cutting surface. Inside the workpiece, the intrinsic stacking faults, dislocation loops and atomic clusters are generated. The SFT is produced inside the workpiece when the cutting direction is along 100 . However, a "V"-shape dislocation loop is generated when the cutting direction is along 101 .

(3) As the cutting direction changes from 101 to 100 , the tangential force and normal force is increased by 13.32% and 10.25%, respectively. Quantitative analysis of atomic distribution along the *y*-direction proves that the machined surface quality along 100 is better than that along 101 . The thickness of the amorphous layer along 100 is slightly deeper than that along 101 .

**Author Contributions:** Conceptualization, R.F.; methodology, R.F. and J.L.; software, J.L.; formal analysis, H.Q.; investigation, H.Q. and Y.Q.; resources, H.L.; data curation, H.L. and M.W.; writing—original draft preparation, H.Q.; writing—review and editing, R.F. and J.L.; visualization, H.Q. and J.L.; supervision, R.F. and C.L.

**Funding:** This research was funded by a grant from the China Scholarship Council, grant number 201808625035, the National Natural Science Foundation of China, grant number 51865027; the Program for Changjiang Scholars and Innovative Research Team in University of Ministry of Education of China, grant number IRT\_15R30; and the Hongliu First-class Disciplines Development Program of Lanzhou University of Technology.

**Acknowledgments:** The technical support of Gansu Province Supercomputer Centre is gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Fracture Toughness Calculation Method Amendment of the Dissimilar Steel Welded Joint Based on 3D XFEM**

### **Yuwen Qian 1,2 and Jianping Zhao 1,2,\***


Received: 21 March 2019; Accepted: 27 April 2019; Published: 30 April 2019

**Abstract:** The dissimilar steel welded joint is divided into three pieces, parent material–weld metal–parent material, by the integrity identification of BS7910-2013. In reality, the undermatched welded joint geometry is different: parent material–heat affected zone (HAZ)–fusion line–weld metal. A combination of the CF62 (parent material) and E316L (welding rod) was the example undermatched welded joint, whose geometry was divided into four pieces to investigate the fracture toughness of the joint by experiments and the extended finite element method (XFEM) calculation. The experimental results were used to change the fracture toughness of the undermatched welded joint, and the XFEM results were used to amend the fracture toughness calculation method with a new definition of the crack length. The research results show that the amendment of the undermatched welded joint geometry expresses more accuracy of the fracture toughness of the joint. The XFEM models were verified as valid by the experiment. The amendment of the fracture toughness calculation method expresses a better fit by the new definition of the crack length, in accordance with the crack route simulated by the XFEM. The results after the amendment coincide with the reality in engineering.

**Keywords:** undermatched; integrity identification; XFEM; fracture toughness calculation method

### **1. Introduction**

### *1.1. The Undermatched Welded Joint*

High-strength low-alloy (HSLA) steel is commonly used for petrochemical equipment. The development of special electrodes for the new HSLA is expensive and time-consuming. Ready-made lower strength electrodes for HSLA are a good choice to fabricate the undermatched welded joint. The undermatched welded joint has a low crack tendency and is economical within a wide range of matching ratios [1].

### *1.2. The Integrity Identification of the Undermatched Welded Joint*

The integrity identification of BS7910-2013 [2] shows the three pieces of the undermatched welded joint, parent material–weld metal–parent material, which is like a sandwich structure. In reality, the undermatched welded joint geometry is different [3], parent material–heat affected zone (HAZ)–fusion line–weld metal, as Figure 1a shows. Figure 1b shows the different grain sizes of the different areas of the undermatched welded joint. The blue line in the figure is the line between the HAZ and the parent material. The combination of CF62 (parent material) and E316L (welding rod) is the example undermatched welded joint in this study. The fracture toughness of the different areas was researched by experiments and extended finite element method (XFEM) simulations in this paper. The XFEM models of the different areas of the undermatched welded joint were built and verified as valid by the experiments.

**Figure 1.** Different areas of the investigated welded joint. (**a**) Different local areas; (**b**) Grain.

### *1.3. The Research Status of XFEM*

The extended finite element method (XFEM) was put forward by Belytschko [4], an American professor, and was used to research the crack growth route. Compared to the traditional finite element, the XFEM simulated the crack initiation and growth without the need to refresh the meshing, reducing the workload. Sukumar et al. [5] promoted the 2D XFEM crack model to a 3D XFEM crack model. Legay [6] and Ventura et al. [7] researched the blending element (shown in Figure 2) to improve the convergence and precision of the simulation. Menouillard [8] improved integral stability. Zhuang and Cheng [9] realized the fusion line crack extension between the dissimilar steels by the XFEM.

**Figure 2.** Planar crack, cut element, and blending element.

Research of the XFEM began late in China. Xiujun Fang of Tsinghua University researched the cohesive crack model based on the XFEM in 2007 [10]. Hai Xie of Shanghai Jiaotong University developed the XFEM with the ABAQUS User Subroutine in 2009 [11]. Yehai Li of Nanjing University of Aeronautics and Astronautics researched the interfacial crack growth of the dissimilar material in 2012 [12]. Shaoyun Zhang of Zhejiang University researched the helical crack growth by the XFEM in 2013 [13]. Zhifeng Yang of Nanjing Tech University developed the 2D elastic-plastic crack XFEM model and verified its validity by the *J* integral in 2014 [14]. Yixiu Shu of the Northwestern Polytechnical University analyzed the multiple crack interaction problems by the XFEM in 2015 [15]. Yang Zhang of Harbin Engineering University applied a 3D crack to the pressure vessel by the XFEM in 2015 [16]. Zhen Wang and Tiantang Yu researched the adaptive multiscale XFEM model for 3D crack problems in 2016 [17].

### *1.4. The Welding of the Undermatched Joint*

The size of the weld plate is 600 × 400 × 40 mm (length, width, and thickness, respectively) with an X-shaped groove, shown in Figure 3. The welding rod (E316L) is the special electrode for 316L stainless steel. The welding parameters are listed in Table 1, in accordance with the criteria of ASME-VIII [18]. Tables 2 and 3 list the chemical composition tested results of CF62 and E316L. Tables 4 and 5 list the mechanical properties of the CF62 steel and the E316L steel, respectively.

**Figure 3.** The investigated welded joint.

**Table 1.** Welding parameters of the undermatched welded joint adapted from [18], with permission from American Society of Mechanical Engineers, 2019.


**Table 2.** Chemical composition of CF62 (wt.%).


**Table 3.** Chemical composition of E316L(wt.%).


**Table 4.** Mechanical properties of CF62 reproduced from [19], with permission from Taiyuan University of Science and Technology, 2019.



**Table 5.** Mechanical properties of E316L reproduced from [20], with permission from welding technology, 2019.

### *1.5. The Fracture Toughness Calculation Method*

ASTM 1820-18 [21] is the newest standard test method for the measurement of fracture toughness. Fracture toughness is calculated by the relationship between the *J* integral and the crack length, *ai* [21]:

$$J = \frac{K^2 \left(1 - \nu^2\right)}{E} + J\_{pl} \tag{1}$$

$$J\_{pl} = \frac{\eta\_{pl} A\_{pl}}{B\_N b\_o} \tag{2}$$

where:

*Apl* is the area under force versus displacement, as shown in Figure 4;

η*pl* is either 1.9 if the load-line displacement is used for *Apl* or 3.667 − 2.199(*a*0/*w*) + 0.437(*a*0/*w*) <sup>2</sup> if the recorded crack mouth opening displacement is used for *Apl*;

*BN* is the net specimen thickness;

*bo* is W − *a*0.

**Figure 4.** Definition of the area for the *J* calculation using the basic method.

The crack length, *ai*, is calculated by the displaced force point with Equation (3):

$$\frac{a\_i}{W} = 1.000196 - 4.06319\mu + 11.242\mu^2 - 106.043\mu^3 + 464.335\mu^4 - 650.677\mu^5 \tag{3}$$

where:

$$\mu = \frac{1}{\left(B\_{\epsilon}E \ C\_{\mathfrak{c}(i)}\right)^{1/2} + 1} \tag{4}$$

$$B\_{\varepsilon} = B - \left(B - B\_{N}\right)^{2}/B \tag{5}$$

$$C\_{c(i)} = \frac{C\_i}{\left(\frac{H^\*}{R\_i}\sin\theta\_i - \cos\theta\_i\right)\left(\frac{D}{R\_i}\sin\theta\_i - \cos\theta\_i\right)}\tag{6}$$

where:

*Ci* is the measured specimen elastic compliance (Δvm/ΔPm) (at the load-line);

*H*∗ is the initial half-span of the load points (center of the pin holes);

*D* is one-half of the initial distance between the displacement measurement points;

*Ri* is the radius of the rotation of the crack center line, (W+a)/2, where a is the updated crack size.

Figure 5a shows the relationship between the crack rotation angle, θ*i*, and the respective crack length, *ai* [21]. The crack rotation angle, θ*i*, is calculated by Equation (7). Figure 5a shows that the crack length, *ai*, was equal to the length parallel to the center line of the compact tension (CT) samples. The crack deflection has not been considered in the crack length calculation in the criteria.

$$\theta\_i = \sin^{-1}\left\{ \left( D + \frac{V\_{m\{i\}}}{2} \right) \left( D^2 + R\_i^2 \right)^{1/2} \right\} - \tan^{-1}\left( \frac{D}{R\_i} \right) \tag{7}$$

where *Vm*(*i*) is the total measured load-line displacement at the beginning of the *i*-th unloading/reloading cycle.

**Figure 5.** Crack rotation angle, θ, of the compact tension (CT) specimens. (**a**) Initial crack; (**b**) Crack after deflection.

Figure 5b shows the crack route after the crack deflection. It shows the crack length after the deflection, *a* = *ax* <sup>2</sup> + *ay* 2. The crack length in the criteria of ASTM 1820-18 is *ai* = *ay*, which means that *ai* < *a*. Therefore, the crack length, *ai*, calculated by the criteria of ASTM 1820-18, was smaller than the true crack length after the deflection, *a*. The fracture toughness calculation method was then amended by the crack length after the deflection was simulated by the XFEM.

The division of the geometry (parent material–HAZ–fusion line–weld metal) of the undermatched welded joint and the experimental results changed the fracture toughness of the dissimilar steel welded joint. The crack length amendment after the deflection defined a new fracture toughness calculation method of the dissimilar welded joint.

### *1.6. J Integral and the Fracture Toughness Value,Kmat*

The *J* integral was put forward by Rice [22]. Figure 6a shows the integral loop/area around the crack tip. Equation (8) is calculated with the equivalent integral region to deal with the stress and the strain of Figure 6b [22]. Figure 6c shows the integral area, which contains a heterogeneous material interface. In [23], it was verified that the interface integral area *I th <sup>A</sup>*<sup>+</sup> = *<sup>I</sup> th <sup>A</sup>*<sup>−</sup> = *Iinter f ace* = 0 when the crack tip is under a constant temperature with a pure mechanical load. Equation (8) is applicable to calculate the *J* integral for the investigated welded joint in this paper.

$$J = \int\_{A} \left( \sigma\_{i\bar{j}} \frac{\partial \mathbf{u}\_{\bar{j}}}{\partial \mathbf{x}\_{1}} - \alpha \delta\_{1i} \right) \frac{\partial q}{\partial \mathbf{x}\_{i}} dA \tag{8}$$

where *ui* is the component of the displacement vector; *dA* is the tiny increment area on the integral path, Γ; ω is the strain energy density factor; and *q*(*x*,*y*) is a mathematical function, which is *q* = 0 in the outer boundary and *q* = 1 in the other places of the regional; and ω, *Ti* are listed in Equation (9).

$$\begin{aligned} \rho &= \int\_{0}^{\varepsilon\_{ij}} \sigma\_{ij} d\varepsilon\_{ij} \\ T\_{i} &= \sigma\_{ij} n\_{j} \end{aligned} \tag{9}$$

where σ*ij* is the stress tensor and ε*ij* is the strain tensor.

**Figure 6.** Counterclockwise integral loop around the crack tip. (**a**) Integral loop; (**b**) Integral area; (**c**) Interface integral area.

BS7910-2013 [2] gives Equation (10), which shows the relationship between the *Jmat* value and the fracture toughness value, *Kmat*, during the elastic stage:

$$J\_{mat} = \frac{K\_{mat}^2 \left(1 - \nu^2\right)}{E} \tag{10}$$

where *E* is the elastic modulus and ν is the Poisson's ratio.

### **2. Experiments**

### *2.1. The Fracture Toughness Experiments*

The fracture toughness of the different areas of the joint geometry (parent material–HAZ–fusion line–weld metal) was tested with the individual CT specimens that were taken from the undermatched welded joint. The crack initiated and grew along the center line of the CT specimens during the experiments. The center lines in the different areas of the joint present the locations of the CT specimens in the dissimilar joint, as seen in Figure 7. The CT size is shown in Figure 8.

**Figure 7.** Center lines of the CT specimens in the investigated welded joint.

**Figure 8.** Size of the CT specimens.

The fracture toughness was tested with the MTS-809 (MTS Systems Cor., Eden Prairie, MN, USA). The local fracture surfaces taken by SEM of the different areas of the joint geometry are shown in Figure 9. Figure 9b shows that the crack initiated from the high-strength region (parent material) and deflected to the lower strength region (weld metal) in the HAZ. Figure 10 depicts the HAZ fracture surface macro-image, which indicates the parent material (PM) region and the weld metal

region. It means that the crack deflection existed in the dissimilar material during the experiment. The experimental process was in accordance with the criteria of ASTM 1820-18 [21], and the results of the *J* integral and the crack length, Δ*a*, were managed under the regulations.

**Figure 9.** Fracture surfaces after the fracture toughness experiment. (**a**) Parent material (PM); (**b**) HAZ; (**c**) Fusion line; (**d**) Weld metal.

**Figure 10.** HAZ fracture surface macro-image.

### *2.2. The Tensile Experiments*

The maximum stress damage of the parent material and the weld metal was needed to build the XFEM models. It was tested with a flat tensile specimen, which was 0.8 mm in thickness. Figure 11a shows the location of the tensile specimen, which was taken from the parent material to the weld metal and is parallel to the fusion line. Figure 11b shows the shape and dimensions of the flat tensile specimen.

**Figure 11.** The location and size of the flat tensile specimen. (**a**) Location; (**b**) Shape and dimensions.

### **3. XFEM Models**

A pre-crack was needed when the 1/2 CT specimen was tested and the pre-crack length was 1.5 mm, in accordance with the criteria of the ASTM 1820-18 regulation. The XFEM models were built with a 1.5 mm pre-crack.

The four XFEM models were the parent material, weld metal, fusion line (a combination of PM and weld metal), and HAZ. The mechanical properties were tested by the flat tensile specimen. There were two important parameters—the maximum stress damage, −σ*max*, and the equivalent energy release rate, *G*θ*c*. In [14,24], the relationship between, σ*max*, and the tensile strength, σ*b*, is given as σ*max* = σ*b*. When the maximum energy release rate (MERR), *G*θ*max*, is equal to the threshold, *G*θ*c*, the crack begins to initiate. The MERR, *G*θ*max*, is calculated in Equation (11), and the threshold is calculated in Equation (12) [13]:

$$\mathbf{G\_{\mathcal{O\max}}} = \frac{1-\nu^2}{4E}\mathbf{1} + \cos\theta\_0 \mathbf{[K\_I^2[1+\cos\theta\_0]-4K\_\mathrm{I}K\_\mathrm{II}\sin\theta\_0 + K\_\mathrm{II}^2[5-3\cos\theta\_0]]} \mathbf{G\_{\mathcal{O\max}}} \tag{11}$$

where *E* is the elastic modulus, ν is the Poisson's ratio, *K*<sup>I</sup> is the applied tensile stress intensity factor of Mode I, and *K*II is the applied tensile stress intensity factor of Mode II.

There are three kinds of crack mode—the opening crack mode (Mode I), the sliding crack mode (Mode II), and the tearing crack mode (Mode III). The experimental and XFEM simulation load is perpendicular to the crack surface. The research crack in this paper is the opening crack mode(Mode I). The MERR threshold, *G*θ*c*, is related to the critical stress intensity factor of Mode I, which is *K*II = 0:

$$\mathbf{G}\_{\partial \mathbf{c}} = \mathbf{G}\_{\mathbf{l} \mathbf{c}} = \frac{1 - \nu^2}{E} \mathbf{K}\_{\mathbf{l} \mathbf{c}}^2 \tag{12}$$

Since the crack mode in this paper is Mode I, *G*I*<sup>c</sup>* = *G*II*<sup>c</sup>* = *G*III*<sup>c</sup>* = *G*θ*<sup>c</sup>* [13]. The XFEM simulation properties of CF62 and E316L are listed in Table 6.


**Table 6.** Mechanical properties of the crack propagation.

Figure 12 shows the meshing of the XFEM of the 1/2CT. The sparse meshing elements in the XFEM are applicable [25–27]. Figure 13 shows the dissimilar material XFEM models of the joint geometry with the material printed on their parts. The parent material and weld metal XFEM models are not shown since they are the uniform material model, i.e., they are the material printed on the whole model.

**Figure 12.** The meshing model of the extended finite element method(XFEM).

**Figure 13.** The dissimilar material XFEM models. (**a**) Fusion line; (**b**) HAZ.

### **4. The Fracture Toughness of the Undermatched Welded Joint**

### *4.1. Fracture Toughness of the Undermatched Welded Joint*

The test results of the *J* integral and crack length, Δ*a*, are fitted by Equation (13) in the criteria of ASTM 1820-18 [21]:

$$J = \alpha + \beta (\Delta a)^r \tag{13}$$

The experimental data and the resistance equations of the different areas of the joint geometry are listed in Tables 7 and 8, respectively. The resistance curves are shown in Figure 14.


**Table 7.** Experimental data.


**Table 8.** The resistance equations and the *Jmat* and *Kmat* value of the undermatched welded joint.

**Figure 14.** The ensemble of the *J-*Δ*a* resistance curves of the undermatched welded joint.

### *4.2. Amendment of HAZ Fracture Toughness*

Figure 15a shows the fracture surface of the HAZ and shows that there was a 1.04 mm length of the smooth fracture surface and a 1.07 mm length of the parent material fracture surface. The fusion line between the parent material and the weld metal is marked with a red line. The weld metal fracture surface is below the red line.

**Figure 15.** Fracture toughness results of the HAZ. (**a**) Fracture surface; (**b**) Resistance curve regions of the HAZ.

Figure 15b shows the resistance curve regions of the HAZ. It shows that in the initial crack (Δ*a* < 1.04 mm) the *R*-curve fits the resulting data points well. The parent material region is next to the fitting region: 1.04 mm < Δ*a* < 2.11 mm. The *J* integral is larger than the curve point which shows the parent material's fracture toughness. After the crack grew across to the weld metal region, the *J* integral was smaller than the *R*-curve points.

Haitao Wang tested and simulated the crack growth paths in a dissimilar metal welded joint and found that the crack deflected to the soft area during the experiment [28,29]. Longchao Cao [30] and Pavel Podany [31] researched the microstructure and the local mechanical properties of the dissimilar butt joints. Jianbo Li [32] studied the fracture toughness near the fusion line (HAZ and fusion line) of the dissimilar material and, in addition, found the crack deflection. The HAZ experimental results demonstrated the crack deflection when the CT specimens were tested. The resistance curve of the HAZ (Figure 15b) shows that the fracture toughness of the undermatched welded joint needs to be amended. The crack deflection shows that the crack length test method needs to be amended as well.

As there was a misfit of the data points in the two regions and the lower strength region was the weld metal region, which was below the high-strength region, the HAZ was considered to be the high-strength region. The crack initiated in the HAZ high-strength region and then deflected to the weld metal region. During the crack growth process, the crack crossed the fusion line between the high-strength region and the weld metal region. Figure 9b shows the local fracture surface figure taken by SEM, which shows the location of the high-strength region and the weld metal region (lower strength region).

Table 9 lists the resistance equations after the amendment of the different areas (parent material–HAZ–fusion line–weld metal) of the undermatched welded joint geometry. The fusion line crack deflection initiated in the pre-crack stage and the crack growth of the fusion line (the heterogeneous material XFEM model) were always in the weld metal region.

Figure 16 shows the resistance curves of the high-strength regions of the HAZ after amending. The figure shows that the HAZ is more applicable between the resistance curve and the data points after amending.

**Table 9.** The resistance equations and the *Jmat* and *Kmat* values of the undermatched welded joint after amending.


Note: the crack growth of the fusion line was always in the weld metal region with the deflection of the pre-crack.

**Figure 16.** Resistance curves of the HAZ after amending.

### *4.3. J Integral Verification of XFEM*

The XFEM gives the *J* integral of the crack growth, which is used to compare the experimental *J-*Δ*a* resistance curves. The XFEM models of the parent material and the weld metal were verified, respectively, by the experimental results. Figures 17 and 18 are the XFEM results of the parent material and the weld metal. Figures 17a and 18a show the CT stress distribution, and Figures 17b and 18b show the crack route.

**Figure 17.** XFEM result of the parent material. (**a**) CT stress distribution; (**b**) Crack route.

**Figure 18.** XFEM result of the weld metal. (**a**) CT stress distribution; (**b**) Crack route.

There are five *J* integral routes when the XFEM model outputs the *J* integral and the results are the same. The *J* integral result is independent of the route, which has been demonstrated previously [27,33,34].

Figure 19 shows the comparison of the experimental *J-*Δ*a* resistance curves and the *J* integral data point output of the XFEM. The figure shows that the data point output of the XFEM coincides well with the experimental *J-*Δ*a* resistance curves. This means that the XFEM models in this paper are valid for the undermatched welded joint.

**Figure 19.** *J* integral verification of the XFEM. (**a**) Parent material; (**b**) Weld metal.

### **5. Amendment of the Fracture Toughness Calculation Method**

Table 10 lists the resistance equations of the investigated welded joint after the crack length amendment. Figure 20 shows the homogeneous material XFEM models' (PM and weld metal) resistance curves after amendment, while Figure 21 shows the heterogeneous material XFEM models' (fusion line and HAZ) resistance curves after amendment. Figure 21a shows the resistance curve of the fusion line, and Figure 21b shows that of the HAZ high-strength region.

**Table 10.** The resistance equations and the *Jmat* and *Kmat* values of the undermatched welded joint after the crack amendment based on the XFEM.


Note: the crack growth of the fusion line was always in the weld metal region with the deflection of the pre-crack.

**Figure 20.** Resistance curves after the crack length amendment of the homogeneous material. (**a**) Parent material; (**b**) Weld metal.

**Figure 21.** Resistance curves after the crack length amendment of the heterogeneous material. (**a**) Fusion line; (**b**) HAZ.

The crack length amendment results show that the fracture toughness calculation method of the undermatched welded joint needed to be amended, due to the crack growth deflection of the heterogeneous material.

The fracture toughness of the investigated welded joint was given by experimenting with the different areas of the joint (parent material–fusion line–HAZ–weld metal). In addition, the fracture toughness calculation method was amended based on the XFEM models of the heterogeneous material in this paper.

### **6. Conclusions**

The division of the dissimilar steel welded joint geometry of the current integrity identification needed to be amended to four pieces: parent material–HAZ–fusion line–weld metal. In accordance with this new division of the undermatched welded joint, XFEM models of the undermatched welded joint were built to calculate their fracture toughness. The conclusions are as follows:


**Author Contributions:** Y.Q. is mainly on the Methodology, Software, Validation, Formal Analysis, Investigation, Data Curation, Writing-Original Draft Preparation, Writing-Review & Editing. J.Z. is mainly on Writing-Review & Editing, Visualization, Supervision, Funding Acquisition.

**Funding:** This research was funded by the National Key Research and Development Program of China [2016YFC0801905].

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Using DEFORM Software for Determination of Parameters for Two Fracture Criteria on DIN 34CrNiMo6**

### **Ivana Poláková \*, Michal Zemko, Martin Rund and Ján Džugan**

COMTES FHT, Pr ˚umyslová 995, Dobˇrany 334 41, Czech Republic; mzemko@comtesfht.cz (M.Z.); mrund@comtesfht.cz (M.R.); jdzugan@comtesfht.cz (J.D.) **\*** Correspondence: ivana.polakova@comtesfht.cz; Tel.: +420-377-197-321

Received: 27 February 2020; Accepted: 26 March 2020; Published: 28 March 2020

**Abstract:** The aim of this study was to calibrate a material model with two fracture criteria that is available in the DEFORM software on DIN 34CrNiMo6. The purpose is to propose a type of simple test that will be sufficient for the determination of damage parameters. The influence of the quantity of mechanical tests on the accuracy of the fracture criterion was explored. This approach was validated using several tests and simulations of damage in a tube and a round bar. These tests are used in engineering applications for their ease of manufacturing and their strong ability to fracture. The prediction of the time and location of the failure was based on the parameters of the relevant damage model. Normalized Cockroft-Latham and Oyane criteria were explored. The validation involved comparing the results of numerical simulation against the test data. The accuracy of prediction of fracture for various stress states using the criteria was evaluated. Both fracture criteria showed good agreement in terms of the fracture locus, but the Oyane criterion proved more suitable for cases covering larger triaxiality ranges.

**Keywords:** flow stress; stress triaxiality; material damage; FEM simulation

### **1. Introduction**

The aim of this paper is to develop a procedure for failure prediction through simulations when a limited amount of tests are available. The purpose of the paper is to demonstrate quick and sufficiently accurate predictions based on very few—and if possible simple—experiments. From engineering and practical points of view, the number of tests on specimens should be no higher than absolutely necessary. When it comes to calibrating fracture criteria, the challenge lies in choosing the right tests [1]. DEFORM software (developed by Scientific Forming Technologies Corporation, Columbus, OH, USA) was used for all computational modelling tasks presented in this paper. Two fracture criteria available in DEFORM were chosen for validating the procedure of failure prediction: the Normalized Cockcroft-Latham criterion (NCL) and the Oyane criterion.

An accurate description of material behaviour is a prerequisite for successful modelling of cold metal forming processes. Strain hardening (flow stress) and fracture behaviour depend on basic material properties, which are governed by processing history, typically the temperature and strain rate profiles. Such properties can be determined through experiments and testing. In forming processes, where fracture must be avoided, identification of forming limits is a matter of key importance. However, the ability to determine the conditions that lead to fracture is important for processes in which it is undesirable, as well as in processes where material separation is the desired outcome (such as cutting).

Numerous researchers have attempted to predict damage and fracture [2–4] and a multitude of models have been proposed over the years. These include Bridgman [5], McClintock [6] and Rice and Tracey [7], who demonstrated the effect of hydrostatic pressure on fracture. Wilkins et al. postulated that both hydrostatic pressure and deviatoric stress contribute to damage. Johnson and Cook [8] introduced a fracture model that accounted for strain rate and temperature.

Bao and Wierzbicki [9] performed numerous experiments involving a wide range of triaxialities and suggested that the dependence of fracture strain on stress triaxiality was not a monotonic curve. They claimed that there were different fracture mechanisms operating in different stress triaxiality ranges. Later, these authors proposed a cut-off value for triaxiality at negative 1/3, where fracture strain becomes infinity. Kim et al. [10], Gao et al. [11], and Brünig et al. [12] studied triaxiality and Lode angle. This has been followed up by many authors, including our laboratory, to develop testing procedures and the shape of the samples [13–15].

The failure prediction procedure described in this paper consists of defining material behaviour for numerical simulation on the basis of several different tests and finding correlation between them and the accuracy and usability of the prediction. For this purpose, mechanical tests classified into three groups were performed. The first one involved tensile and compression tests through which flow stress data were acquired. The second group comprised notched tensile test, shear test and plane strain test conducted to determine fracture parameters. The last group of tests were performed to validate the fracture criteria. These included radial compression of a ring specimen and axial compression of a notched cylinder.

### **2. Material**

Fracture modelling for tests was carried out on DIN 34CrNiMo6 steel at room temperature and at quasi-static strain rates. This steel was chosen due to its versatile engineering applications [16]. It is used for heavy-duty components, such as axles, shafts, crankshafts, connecting rods, valves, propeller hubs, gears, couplings, torsion bars, aircraft components, and heavy parts of rock drills.

The advantages of DIN 34CrNiMo6 martensitic steel include high ductility, hardenability and strength [16]. In fact, it provides the highest yield strength and ultimate tensile strength from the entire class of high strength steels. Generally, martensitic steels are characterized by a martensitic matrix produced by quenching with small amounts of ferrite or bainite. This steel is strengthened by the precipitation of a fine dispersion of carbides during tempering. Its nominal chemical composition, in weight percent, is presented in Table 1. The material was heat-treated to 34 HRC. Basic measured mechanical data are shown in Table 2.


**Table 1.** Nominal chemical composition of DIN 34CrNiMo6.



### **3. Uncoupled Fracture Models**

Fracture models are divided into coupled and uncoupled categories. Coupled models are complex, require a great many material constants, and are computationally demanding. In uncoupled models, the plasticity model is separate from the failure model. This means a significant reduction in the computation effort. The DEFORM software only offers uncoupled fracture criteria.

Two fracture criteria were selected for this paper. Their applicability to industrial problems was evaluated. The NCL criterion has only one parameter which governs the critical damage. Therefore, the nature of the problem at hand needs to be accounted for. Identifying the applicability range of a criterion is an essential task. Nevertheless, the NCL criterion is based on maximal principal stress, which is dependent on both the triaxiality and the Lode angle. The NCL criterion only averages the

effect of both parameters using maximal principal stress. Generally, fracture criteria with multiple parameters are applicable to a wider range of stress states. One of them is the Oyane criterion, the other criterion explored in this work. The Oyane criterion only includes stress triaxiality dependence. Oyane prefers triaxiality, which seems to be more significant for crack initiation and evolution.

The stress state in a material is affected by external loads. Tensile stress leads to damage more rapidly than compression. Physical damage is activated by microvoid nucleation and propagates through their growth and coalescence. The development of damage is studied in correlation to stress triaxiality. Stress triaxiality is defined as a ratio of the mean hydrostatic stress σ*<sup>m</sup>* to the equivalent stress σ*eq*:

$$
\eta = \frac{\sigma\_m}{\sigma\_{eq}} = \frac{\sigma\_1 + \sigma\_2 + \sigma\_3}{3\sigma\_{eq}},
\tag{1}
$$

where σ*<sup>i</sup>* (*i* = 1, 2, 3) denotes the principal stress. The equivalent stress σ*eq* is defined in terms of three principal stresses. Assuming that the material is isotropic, the equivalent stress is given by the von Mises equation:

$$
\sigma\_{\varepsilon q} = \frac{1}{\sqrt{2}} \sqrt{\left[ \left( \sigma\_1 - \sigma\_2 \right)^2 + \left( \sigma\_2 - \sigma\_3 \right)^2 + \left( \sigma\_3 - \sigma\_1 \right)^2 \right]} \tag{2}
$$

The equivalent strain ε*eq* is expressed as:

$$\varepsilon\_{eq} = \frac{\sqrt{2}}{3} \sqrt{\left[ \left( \varepsilon\_1 - \varepsilon\_2 \right)^2 + \left( \varepsilon\_2 - \varepsilon\_3 \right)^2 + \left( \varepsilon\_3 - \varepsilon\_1 \right)^2 \right]} \tag{3}$$

Many published experimental reports [17] have shown that fracture strain decreases with increasing triaxiality, although the dependence is not monotonic as in Figure 1. A peak is normally found near the 1/3 value, which corresponds to a smooth tensile test. Another factor in the occurrence of fracture is the Lode parameter. This is related to the third deviatoric invariant and reflects shear deformation. The value of the Lode parameter is 0 for plane strain conditions, +1 for axisymmetric tension, and −1 for axisymmetric compression. Data from various experimental tests plotted in a fracture locus graph provide the basis for fracture prediction, see example in Figure 1, where the dependence of the fracture strain on stress triaxiality and the Lode parameter is plotted. The black line in the graph correspond to plane stress. The critical level of accumulated damage is reached when equivalent strain exceeds a limit value.

**Figure 1.** Representation of fracture locus in space.

In computer simulation, the critical value of a damage criterion indicates the initiation of fracture. The actual fracture strain can be determined from a simulation of the real test. In DEFORM code,

the growth and coalescence of voids that lead to the evolution of fracture can be modelled by either the softening method or by element deletion. The critical damage value depends on the stress state in the test and on the fracture model used.

In the NCL criterion, the largest principal stress is the most relevant factor for the initiation of fracture:

$$D\_{\mathbb{C}} = \int\_{0}^{\varepsilon\_{f}} \frac{\sigma\_{1}}{\sigma\_{eq}} d\varepsilon \,, \tag{4}$$

where *Dc* denotes the critical damage value, ε*<sup>f</sup>* is the fracture strain, σ<sup>1</sup> stands for the maximal principal stress and σ*eq* is the equivalent stress. The maximal principal stress can be expressed by stress triaxiality and Lode angle [18,19]:

$$
\sigma\_1 = \sigma\_{\rm eq} \Big( \eta + \frac{2}{3} \cos \frac{\pi}{6} (1 - \overline{\theta}) \Big). \tag{5}
$$

where η is triaxiality and θ is a normalized Lode angle.

The Oyane fracture criterion is based on stress triaxiality as the main reason for fracture:

$$\mathbf{D}\_{\mathbb{C}} = \int\_{0}^{\varepsilon\_{\mathrm{f}}} \left[ 1 + \frac{1}{\mathbf{a}\_{0}} \frac{\sigma\_{\mathrm{m}}}{\sigma\_{\mathrm{eq}}} \right] d\varepsilon,\tag{6}$$

where *a*<sup>0</sup> is a material coefficient to be determined experimentally and the ratio σ*m*/σ*eq* denotes the stress triaxiality.

The dimensionless normalized Lode angle [17] is defined in terms of invariants of the stress tensor:

$$I\_1 = p = -\sigma\_m = -\frac{1}{3}(\sigma\_1 + \sigma\_2 + \sigma\_3),\tag{7}$$

$$I\_2 = q = \sigma\_{eq} = \sqrt{\frac{1}{2} \left[ \left( \sigma\_1 - \sigma\_2 \right)^2 + \left( \sigma\_2 - \sigma\_3 \right)^2 + \left( \sigma\_3 - \sigma\_1 \right)^2 \right]} \tag{8}$$

$$I\_3 = r = \frac{27}{2} [ (\sigma\_1 - \sigma\_m)(\sigma\_2 - \sigma\_m)(\sigma\_3 - \sigma\_m) ]^{\frac{1}{3}} \,\,\, \,\tag{9}$$

$$\overline{\theta} = 1 - \frac{2}{\pi} \arccos\left(\frac{r}{q}\right)^3,\tag{10}$$

where *Ii* (*i* = 1, 2, 3) are invariants of the stress tensor and θ is the normalized Lode angle parameter. Only the NCL criterion has a built-in dependence on the Lode angle. The effect of the Lode angle on Oyane criterion can be done by separating the fracture coefficients for tests with different Lode angle.

### **4. Experiments**

In the following sections, three groups of tests are described. The basic tests include tensile and compression tests. Tests for fracture coefficients determination comprised a tensile test of a notched specimen, plane strain, and shear tests. Validation tests involved notched and unnotched rings and notched cylinder specimens.

### *4.1. Basic Tests*

Tensile testing was carried out at room temperature. The test data were used for flow stress calibration. A drawing of the smooth tensile specimen is shown in Figure 2; the dimensions are in millimetres. The figure also shows the locations and the distance of the landmarks (red crosses) monitored by an optical (laser) extensometer. The rate of loading for smooth tensile specimens was 2 mm/min. The initial strain rate was 0.001 s<sup>−</sup>1. Elongation was measured by optical and mechanical extensometers but only the optical extensometer data were used for calculations.

The other tests described below were recorded using a DIC (Digital Image Correlation) system, except for the tensile tests of smooth and notched specimens, which were measured by an optical

extensometer. A stochastic pattern was applied by airbrush on samples recorded with the DIC system. The actual plastic strain on the specimen surface was calculated using a correlation algorithm implemented in GOM Aramis software [20]. The resolution of cameras was 4248 × 2832 pixels. The frame rate was 5 Hz.

Cylindrical specimens 10 mm in diameter and 15 mm in height were used for compression testing. The locations of and the distance between the landmarks for load-displacement measurement are indicated with red crosses as indicated in Figure 3. The displacement was measured using DIC. The camera was pointed at the specimen and dies, onto which a pattern was applied by spraying to provide contrast for strain measurement and for calculation by means of GOM Aramis software.

**Figure 2.** Drawing of tensile specimen.

**Figure 3.** Drawing of compression specimen. The distance shown in red refers to the dies.

### *4.2. Tests for Fracture Coe*ffi*cient Determination*

The critical value of damage is estimated from plastic strain accumulated prior to fracture, which depends on the stress state. The choice of these tests for calibrating the fracture criteria was intended to cover a wide range of triaxialities and was based on literature survey and previous experience of the authors. They included the notched tensile test (Figure 4), shear test (Figure 5), and plane strain test (Figure 6).

**Figure 4.** Drawing of notched tensile test specimen.

**Figure 5.** Drawing of shear test specimen.

**Figure 6.** Drawing of plane strain test specimen.

Tensile testing of notched specimens at a loading velocity of 0.25 mm/min produced the same strain rate as in the previous tests. Optical and mechanical extensometers were used for measuring the gauge length by the red crosses as indicated in Figure 4. Fractured tensile specimens after testing are shown in Figure 7.

**Figure 7.** Notched tensile test specimens.

The purpose of the shear test and plane strain test was to cover a wider range of stress states for this study. The loading velocity was set at 0.5 mm/min in order to obtain the same strain rate in the relevant region as in notched tensile test. In both tests, extension was measured by DIC. Drawings of shear test and plane strain test specimens are presented in Figures 5 and 6, respectively. The gauge lengths are indicated by red crosses. Some fractured specimens with their surface patterns are shown in Figures 8 and 9.

**Figure 8.** Shear test specimens.

**Figure 9.** Plane strain test specimens.

### *4.3. Validation Tests*

The proposed fracture criteria parameters were validated using two simple tests. One involved radial compression of notched and unnotched rings, see Figures 10 and 11. The other was axial compression of a notched cylinder, as shown in Figure 12.

**Figure 10.** Drawing of notched ring specimen.

**Figure 11.** Drawing of unnotched ring specimen.

**Figure 12.** Drawing of the notched cylinder.

Based on an extensive literature survey, it was found that the presented geometry of test specimens is easy to manufacture and has a strong capacity for fracture initiation. These tests are not standardized, and are used for engineering purposes, such as assessing the damage mechanisms in pipes with and without defect, or the load response of a round bar with a defect.

A loading velocity of 2 mm/min was applied at room temperature. A DIC system was used for measuring the specimen displacement. Photographs of notched and unnotched specimens after testing are shown in Figure 13. Notched cylinders after axial compression are displayed in Figure 14. The cracked specimens clearly demonstrate that the failure was a shear-dominated process in this case.

**Figure 13.** (**a**) Notched and (**b**) unnotched rings after testing.

**Figure 14.** Notched cylinder specimens after testing with fractures in the notch region.

### **5. Numerical Simulation Using DEFORM Software**

### *5.1. Construction of Plasticity Models*

Stress–strain (flow stress) curves for the material were calibrated against data from the smooth tensile and compression tests. Obtaining flow stress curves from uniaxial tensile test data up to the onset of necking is relatively straightforward. This technique has been reported by several authors [21]. However, a disadvantage of tensile testing is that its flow stress data applies at very small strains only.

Much higher strain values are attained in compression than in tension. In compression tests, friction can be a significant issue. It should be reduced as much as possible to prevent specimen barrelling. Using compression test data to expand the flow stress range obtained from tensile tests is a common technique. The resulting flow stress curve is then more accurate at larger strains. The complete set of flow stress data used for simulations in this study is summarized in Table 3 and plotted in Figure 15.


**Table 3.** Flow stress data in a tabular form.

**Figure 15.** Flow stress plot.

Using the flow stress data calculated from tensile and compressive test data, these tests were simulated, replicating the real test conditions. Force-extension diagrams with the actually measured curves and the curves obtained by simulation using fitted flow stress data are shown in Figures 16 and 17 for the smooth tensile test and the compression test, respectively. Both tests were modelled as axisymmetric problems. A fine mesh with elements of 0.25 mm in size was used in the region where fracture was expected to occur. The reason was that strain gradient becomes steep in this region. Remeshing was applied when the specimen shape changed.

**Figure 16.** Comparison of tensile test curves.

None of the smooth tensile specimens fractured at mid-length. Specimens 10\_1 and 10\_2 exhibited almost identical trends, see Figure 16. Specimen 10\_3 fractured earlier than the other two. All the measured compression test curves were almost identical. The simulated load response is in very good agreement with the measured one. The aim of these simulations was to verify the plasticity data to be used for modelling. No fracture criteria were considered at this stage.

**Figure 17.** Comparison of compression test curves.

### *5.2. Identification of Fracture Parameters*

To determine the critical damage threshold for simulation, experimental fracture data are needed [22]. Fracture is assumed to occur when the force-displacement curve begins to drop. This indicates fracture initiation and the amount of fracture strain. However, since the specimen deforms rapidly during necking, the exact fracture strain value it is very difficult to identify from measured test plots. An estimate of fracture strain can also be obtained by measuring the cross-sectional area of a specimen before and after testing [23]. Using FE simulation of tests, the onset of fracture and the fracture strain can be determined more accurately.

The measured and simulated forces are compared in Figure 18 for the notched tensile test, in Figure 19 for the shear test, and in Figure 20. for the plane strain test. These simulations were based on flow stress data from the previous tests. The measured and simulated forces are in very good agreement for all the tests, which means that the flow stress data from smooth tensile and compression tests had been determined accurately. The conditions and constraints used in the simulations of notched tensile test, shear and plane strain tests replicated the actual test conditions. A fine mesh with elements 0.15–0.3 mm in size was used at the mid-length region of the specimen.

The amount of fracture strain in the simulated notched tensile test can be found at the notch tip where failure was expected to occur. The point is indicated by an arrow in Figure 21. In the shear test simulation, a point on the specimen surface was chosen for this purpose where strain reaches the highest values, as in Figure 22. The location of the highest strain in the plane strain specimen is not on the surface but inside the specimen, see Figure 23. The moment when the crack is detected by DIC on the surface of the plane strain specimen may be different. The average strain on the cross-section of the specimen can be smaller than the fracture strain.

**Figure 18.** Comparison of forces for notched tensile test.

**Figure 19.** Comparison of forces for shear test.

**Figure 20.** Comparison of forces for plane strain test.

**Figure 21.** The fracture location in 2D simulation of notched tensile test.

**Figure 22.** The fracture point in simulation of the shear test.

**Figure 23.** The fracture point in simulation of the plane strain test (longitudinal section).

Fracture strain at crack initiation, as determined from simulations calibrated against the actual test data, are listed in Table 4. Cracking was considered to initiate when the measured force curve either dropped off or began to deviate from the simulated force curve. This time instant and location indicate the critical damage value for the NCL criterion which is characterized by normalized Lode angle and triaxiality. The parameters can be obtained directly in the DEFORM simulation data. The critical damage value is identical, accidentally, for the notched tensile specimen and the plane strain specimen.


**Table 4.** Critical damage values identified by simulation.

The calculated triaxiality paths at selected points on the specimens are plotted in Figures 24–26. The triaxiality is changing significantly during the stress evolution, so an average value is usually used and recommended. The average value of triaxiality η*av* is calculated based on the following equation:

$$
\eta\_{\text{av}} = \int\_0^{\varepsilon\_f} \frac{\eta}{\varepsilon\_f} d\varepsilon \,. \tag{11}
$$

The normalized Lode angle paths at identical points in test specimens are plotted in Figures 27–29. The normalized Lode angle for the notched tensile test is close to 1. In the shear test, the normalized Lode angle is close to 0, similar to in the plane strain test. These and the average stress triaxiality values indicate that the choice of these tests was appropriate to reflect the different stress states in the material. All average triaxiality and Lode angle values are given in Table 2.

**Figure 24.** Triaxiality path for notched tensile test.

**Figure 25.** Triaxiality path for shear test.

**Figure 26.** Triaxiality paths for the plane strain test at two points.

**Figure 27.** Lode angle for notched tensile test.

**Figure 28.** Lode angle for shear test.

**Figure 29.** Lode angle for plane strain test.

The Oyane fracture criterion is based on stress triaxiality. The parameters of the Oyane criterion can be evaluated based on literature data [24]. The equation for the criterion can be rewritten as:

$$
\varepsilon\_f = D\_\mathbb{C} - \frac{1}{a\_0} \int\_0^{\varepsilon\_f} \frac{\sigma\_m}{\sigma\_{eq}} d\varepsilon\_\prime \tag{12}
$$

This expression presents a linear relationship between the fracture strain and the integral:

$$\int\_{0}^{\varepsilon\_{f}} \frac{\sigma\_{\text{m}}}{\sigma\_{\text{eq}}} d\varepsilon,\tag{13}$$

which can be obtained from numerical simulation. The graphic representation of Equation (12) is shown in Figure 30. The material constants for the Oyane criterion can be determined from a linear fit to points for the individual tests. Material constant a0 can be derived from the slope of the line as its negative inverse value. *Dc* is in fact the value at the intersection of the line with the ordinate. The values determined in this manner were *Dc* = 0.63 and *a*<sup>0</sup> = 0.26.

When Lode angle is considered for the Oyane criterion, see Figures 27–29**,** the notched tensile test becomes the least relevant of the three tests. Lode angle for the notched tensile test is far from the other two tests. It seems appropriate to have the values of Lode angle close to each other, see the dependence illustrated in Figure 1. In an alternative course of calculation, the notched tensile test was omitted from the evaluation of material constants. The resulting values *Dc* = 0.65 and *a*<sup>0</sup> = 0.27 do not differ substantially from the previous pair of values.

**Figure 30.** Relationship between <sup>ε</sup>*<sup>f</sup>* 0 σ*<sup>m</sup>* <sup>σ</sup>*eq d*ε and ε*<sup>f</sup>* .

Figure 31 shows a graphic representation of the Oyane criterion along with the values from Table 4. in a plot of fracture strain versus triaxiality.

**Figure 31.** Graphic representation of the Oyane criterion.

The damage coefficients for the Oyane criterion are summarized in Table 5.

**Table 5.** Critical damage values for the Oyane criterion.


### **6. Validation of Criteria Coe**ffi**cients and Discussion**

The main purpose of numerical analysis of damage is to predict the failure location and time in a part under load. The validity of the above-calculated critical damage values and material constants can be assessed by comparing numerical simulations with experimental data. This validation was based on the third group of experimental tests: notched and unnotched ring specimens and notched cylinder specimens. All these tests were simulated using NCL and Oyane fracture criteria with accurate replication of the test conditions.

### *6.1. Radial Compression of Ring Specimens*

Two types of ring specimens were tested: notched and unnotched ones. The primary purpose of the notches was to ensure the ring fractures. However, once a simulation of compression of an unnotched ring proved that it would suffer fracture as well, an additional actual test of an unnotched ring was performed to validate the simulation. The simulation models are shown in Figures 32 and 33.

As expected, the numerical analysis confirmed that fracture in the notched ring starts at the notch tip. The fracture propagates inside the ring with the increasing load, opposite the line of contact with the dies. In the unnotched ring, fracture initiates first inside the ring, opposite the outer line of contact with the dies. Then, as loading continues, another fracture occurs on the outer surface where the ring bulges out. Fracture propagation was represented by the softening method implemented in the DEFORM code.

**Figure 32.** Simulation model of compression of a notched ring.

**Figure 33.** Simulation model of axial compression of an unnotched ring.

The simulation model with the NCL criterion had two versions, each for one of the two critical damage values listed in Table 4. The forces from these simulations and from the experiments are compared in Figure 34. The calculated material constants for the Oyane criterion, as presented in Table 5, were substituted into simulation and the resulting simulated forces were compared with the measurement data, as shown in Figure 35.

**Figure 34.** Measured forces in radial compression of rings and the forces simulated with the use of the NCL fracture criterion.

**Figure 35.** Measured forces for notched and unnotched rings and the forces simulated with the use of the Oyane fracture criterion.

The location of the fracture initiation was predicted correctly by simulation with both the NCL and Oyane fracture criteria. For the notched ring test and the NCL criterion, the simulation model with the value of 0.23 predicted a premature fracture, whereas the one with the value of 0.38 predicted failure later than in the experiment. The actual unnotched ring test and its simulation yielded similar results. However, no drop occurred on the simulated force curve for the critical damage value of 0.38. In the simulations, the time of fracture initiation was understood to be the moment the stress decreased, as indicated by the softening method. The point at which the stress begins to fall is indicated by

the red arrow in Figure 34. It illustrates an important finding that although a drop in the simulated force demonstrates that fracture is present, the actual crack initiation might have occurred earlier. During these experiments, the specimen surfaces had to be carefully observed in order to detect fracture. However, the DIC camera was aimed at the specimen face to measure the displacement. The other surfaces were only observed with the naked eye. Hence, the exact time of crack initiation was not determined.

The NCL critical damage value 0.38 in Table 4 applies to the shear test. The value for the notched ring test is the same as the value for the tensile test. The more appropriate value appears to be 0.23, which predicts a slightly earlier failure. The critical damage for the NCL criterion depends significantly on the stress state. This criterion should therefore be used with care.

The simulation of notched ring compression with the Oyane criterion showed very good agreement with the experiment on the notched ring. However, the failure prediction for the unnotched ring was earlier than in the actual test. The calculations with different coefficients gave almost identical results, indicating that the dependence on Lode angle has of little importance for this material with respect to the used triaxiality range from around 0 to 0.7, see Table 4. This result was predictable due to the small difference in the coefficients as shown in graphical representation of Oyane criterion in Figure 31.

### *6.2. Axial Compression of a Notched Cylinder*

Axial compression of a notched cylinder was also used for validating both criteria (Figure 36). The measured and simulated forces for the NCL and Oyane criteria are plotted and compared in Figure 37. In this test, failure is not indicated by a drop in the test force. In the actual tests, the moment of fracture initiation was found by means of DIC measurement. In simulations, the fracture propagation is represented by the softening method. Hence, the crack initiation is identified as a decrease in stress in the relevant region.

**Figure 36.** Simulation model of notched cylinder compression.

In the notched cylinder specimen, the NCL criterion predicted cracking rather early for *Dc* = 0.23. The simulation with *Dc* = 0.38 shows better agreement with the actual test data. With reference to the critical values in Table 4, the NCL value from the shear test (*Dc* = 0.38) appears to better predict fracture in this shear dominated process.

Fracture predictions for notched cylinder obtained with the use of Oyane criterion were very accurate, with the fracture being located in the centre of the range of measured fracture displacement values. The simulation with Oyane criterion was performed for two different coefficients, see Figure 31, that were chosen to assess the influence of the Lode angle. As expected, the difference is small for this material and the used triaxiality range, see Table 4.

**Figure 37.** The comparison between measured and simulated compression forces and fracture initiation points for notched cylinder specimens.

Both validation tests showed quite a good applicability of the Oyane criterion, based on fracture parameters from the simple tests performed in this study. The Oyane criterion covers a wider range of stress states and industrial problems using two inserted parameters, as reflects the dependence on stress triaxiality. The NCL model can also predict failure rather accurately but the choice of the critical damage value is strongly dependent on the stress state in the process. The failure prediction in radial compression of rings was more accurate with *Dc* = 0.23, whereas the simulations of axial compression of notched cylinders gave better results for *Dc* = 0.38. In the rings, tensile stress dominated in fracture initiation locations. By contrast, failure in the notched cylinders was shear dominated. Therefore, a different critical damage value was needed for a reliable failure prediction.

### **7. Conclusions**

Tensile and compression tests were conducted to determine plasticity data for DIN 34CrNiMo6 steel. Notched tensile tests, shear tests and plane strain tests were used for finding critical damage values for the Oyane and NCL criteria. Fracture prediction simulations were carried out using DEFORM code. The choice of tests for sufficiently precise fracture prediction was dictated from an engineering point of view. Plasticity and fracture models were validated using simple tests: radial compression of notched and unnotched rings and axial compression of notched cylinders.

Reducing the inaccuracy of the material flow stress is an important issue for increasing accuracy in subsequent predicting the fracture behaviour. Both fracture criteria predicted the failure locations accurately in the validation tests. The time of fracture was predicted more accurately by the Oyane criterion for both tests. The Oyane criterion can be used for failure prediction in industrial problems without any extra experience of and knowledge with measuring the data for material models.

The NCL criterion lacks general applicability, unlike the Oyane criterion. The effect of both triaxiality and Lode angle in the NCL criterion is hidden in maximal principal stress. The NCL criterion does not allow a damage prediction for arbitrary processes but only for similar processes in which the stress state of material does not differ significantly for the critical damage value. The results indicated that it is essential to identify an applicability range for each critical damage value selected for the NCL criterion.

This paper can serve as a guide for fracture evaluation based on a limited number of simple tests. The effectiveness of this procedure was demonstrated above. Using the fracture criteria available in the DEFORM software, one can obtain the desired results without any programming tasks.

Follow-on investigations will involve validating both fracture criteria for other types of specimens. In addition, other fracture criteria will be examined in order to improve failure prediction in industrial conditions.

**Author Contributions:** Data curation, M.R.; funding acquisition, M.Z.; investigation, I.P.; methodology, M.Z. and J.D.; project administration, J.D.; writing—original draft, I.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** European Regional Development Fund: CZ.02.1.01/0.0/0.0/17\_048/0007350. European Regional Development Fund: CZ.02.1.01/0.0/0.0/16\_019/0000836.

**Acknowledgments:** The paper was developed with a support from the ERDF project Pre-Application Research of Functionally Graduated Materials by Additive Technologies, No. CZ.02.1.01/0.0/0.0/17\_048/0007350 and from ERDF Research of advanced steels with unique properties, No. CZ.02.1.01/0.0/0.0/16\_019/0000836.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
