Cohesive zone modeling for crack propagation in polycrystalline NiTi alloys using molecular dynamics

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

Contents lists available at ScienceDirect

Theoretical and Applied Fracture Mechanics journal homepage: www.elsevier.com/locate/tafmec

Cohesive zone modeling for crack propagation in polycrystalline NiTi alloys using molecular dynamics

T

Min Lua, Fang Wanga, , Xiangguo Zengb, Wenjun Chenb, Junqian Zhangc ⁎

a

School of Materials and Energy, Southwest University, Chongqing 400715, China College of Architecture and Environment, Sichuan University, Chengdu 610065, China c Department of Mechanics, Shanghai University, Shanghai 200444, China b

ARTICLE INFO

ABSTRACT

Keywords: NiTi alloys Voronoi tessellation Molecular dynamics simulations Cohesive zone model Intergranular and transgranular cracking

Nickle-titanium (NiTi) alloys are widely accepted to be one of the most favored engineering materials for use in the structural community. In this work, the formulation and application of a multiscale approach combined with cohesive zone model (CZM) and molecular dynamics (MD) were proposed to investigate the crack propagation of NiTi alloys. The Voronoi tessellation was employed to represent microstructures of such polycrystalline materials. MD simulations were performed to achieve the traction-separation law in the CZM with various initial micro-cracks. The resultant characteristic parameters were embedded in cohesive elements along grain boundaries and in grains, and in turn both intergranular and transgranular fracture were simulated by implementing the finite element analysis with failure criteria. Consequently, the stress intensify factor of compact tension specimens was predicted, and further the relationship between the microstructure and fracture toughness was explored. The results show that there is a reasonable agreement between numerical results and published experimental data, highlighting the adequacy of the new analytical method for reproducing the cracking behavior across micro-, meso-, and macro-scales.

1. Introduction It has been widely accepted that nickel-titanium (NiTi) alloys are currently considered amongst the most advanced materials, because they have excellent super-elasticity, long fatigue life, good corrosion resistance, and superior thermo-mechanical performance [1]. Recently, demands on NiTi alloys have been increasing in numerous structural applications, such as automotive, aerospace, and medical industries [2], which is due to their merits characterized by sustaining large elastic strains and recovering the original shape under heating conditions. It was also reported that Nitinol possesses recoverable strains with more than an order of magnitude greater than traditional alloys [3]. In fact, mechanical properties of polycrystalline materials are primarily dominated by the microstructure composed of grains with different sizes and distributions, grain boundaries, components, as well as inherent defects [4,5]. Thus, material damage and failure exhibit a multi-scale mechanical behavior. At the micro-scale, pre-existing defects nucleate and grow, producing micro-cracks. Afterwards, resulting cracks propagate in grains and along grain boundaries. This occurs at the meso-scale. With the ongoing accumulation of damages, the ultimate failure of materials occurs at the macro-scale. It has been proved

⁎

that local damage evolution determines material fracture toughness, stiffness loss, and other mechanical properties [6]. The complexity of microstructures significantly increases the difficulty in studying the effect of structural topographies on mechanical properties. There is no denying that reliability concerns regarding of NiTi alloys in structural applications have motivated an urgent need for providing a comprehensive knowledge of material deformation mechanism and its fracture behavior from the perspective of micro-structure [7]. There are many numerical and analytical models for solving this case, which are mainly divided into two types: one is based on the classical fracture mechanics and the other is established under the framework of damage mechanics. The former is typically characterized by the analytical solid mechanics method utilization for calculating the driving force on a crack, and the method of experimental solid mechanics for investigating material fracture resistance. For example, Robertson and Ritchie [8] studied the NiTi crack-resistance behavior with a dependence on crack-growth angle experimentally. However, most experimental studies focused on the exceptional characteristics of phase transitions in the NiTi system [9]. Also, the complexity of underlying mechanism of microstructural evolutions brings an enormous challenge to experimental investigation.

Corresponding author. E-mail address: [email protected] (F. Wang).

https://doi.org/10.1016/j.tafmec.2019.102402 Received 3 April 2019; Received in revised form 31 October 2019; Accepted 31 October 2019 Available online 04 November 2019 0167-8442/ © 2019 Elsevier Ltd. All rights reserved.

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

Fig. 1. Schematic of Voronoi tessellation with single-line boundary.

In the latter model, certain state variables are usually utilized for representing the damage effect on mechanical properties [10]. It has been proved that the cohesive zone model (CZM) is deemed as one of the effective methods that could be applied to characterize the crack propagation across the perspective of macro and micro scales [11,12]. An initial work was made by Tvergaard and Hutchinson [13], and Xu and Needleman [14]. They presented a CZM model with traction-separation (T-S) law for controlling the crack propagation on the interface. What the nicest feature this model has is that it could predict transgranular and intergranular fracture through continuum interface separation constitutive relations by internal state variable theory [15]. Combined with CZM and finite element model (FEM), Zhu et al. [16] revealed the failure mechanism through the crack growth simulation of Ni-based superalloys. However, it was found that some adhesive parameters were crucial to determining the T-S law [17], like as cohesive strength and toughness, indicating that those parameters had a significant impact on the finite element simulations. As yet, such parameters involved in the T-S law are always determined through nanocrystalline experiments or unreliable assumptions [18]. However, it is rather difficult for achieving the cohesive zone law through a direct experimental quantification [12]. Consequently, the T-S law parameterization has attracted increasing attention by many scholars. There is evidence that atomistic simulations can be considered as an effective tool for investing the T-S behavior by providing micro-structural details and mechanic insights at the atomic scale [15]. Thus, molecular dynamics (MD) exhibits a possibility of simulating the traction-displacement response, leading to the achievement of resultant parameters. In the last few years, many attempts have been devoted to the extraction of relevant CZM parameters for describing the decohesion law from MD simulations. Komanduri et al. [19] carried out MD simulations on certain single-crystal cubic metals to investigate the nature of deformation and fracture under uniaxial tension at the nanolevel. This work demonstrated that the failure of these materials was due to the coalescence of void, resulting in nanocracks. The subsequent fracture or separation observed was similar to the corresponding behavior at the macroscale. Yamakov et al. [15] reported that a tractiondisplacement relationship, extracted from MD simulations, could be embedded into a CZM for intergranular fracture problem solution in aluminum at the micro-scale. Besides, Zhou et al. [12] established a cohesive surface constitutive model consisting of separate normal and shear T-S relations, promoting the development of a qualitative picture of the T-S behavior as well as the functional forms and parameters possible. Recently, Zeng et al. [20] studied the influence of several MD models with initial defects on the extracted CZM parameters, and conducted FEM simulations to predict the fracture toughness of titanium alloys. Zhou et al. [21] established a polycrystalline model through Voronoi tessellation that represents material microstructures. They investigated the influence of grain size, grain boundary strength, and microcracks on cracking patterns, which was beneficial to the understanding of crack growth behavior at the mesoscale. Unfortunately,

the knowledge of fracture behavior and its microscale mechanism is fairly scarce. An atomic-based CZM jointed with structural characters would be more sufficient for material performance research from the perspective of multi-scales. Contrary to previous approaches, the primary objective of this work is to propose a cross-scale approach combined with the CZM and MD to investigate the cracking behavior of NiTi alloys subjected to tensile loading. Compact tension (CT) specimens are prepared to study the crack growth in microstructures represented by the Voronoi tessellation. The T-S law parameters are achieved through MD simulations, and then they are embedded in zero-thickness cohesive elements along grain boundaries and in grains. Subsequently, the finite element calculation is conducted on the ABAQUS for simulating intergranular and transgranular fracture in polycrystalline structures. Besides, the influence of various micro-defects on the mechanical performance is examined as well. The physical mechanism obtained, combined with mechanical properties, serves as a guide to the development of NiTi alloys for advanced engineering applications. 2. Modeling 2.1. Microstructure Voronoi tessellation is a kind of space division mathematically. For example, if random points are given in the plane, the Voronoi diagram of a set of sites partitions the plane into several Voronoi regions according to the nearest-neighbor rule [22]. Fig. 1 represents the schematic diagram of Voronoi tessellation. Numerous points are randomly generated as the nuclei of Voronoi polygons in the plane region. After linking two points closest to each other forms a straight line, its midperpendicular constructs a polygon boundary. Thus, each point is associated with the region of the plane closest to it [23]. It has been justified that such a method is applied to characterize the microstructure of polycrystalline materials due to its structural randomness and computational efficiency [24,25]. Thus, this work employs the closed Voronoi tessellation to characterize polycrystalline microstructures of NiTi alloys. As seen in Fig. 1, the generated Voronoi tessellation boundary is in the form of single line, and then is taken as grain boundary. The left graph represents the partial enlarged drawing after mesh partition. The meshes inside grains are generated with linear triangular volume finite elements. Afterwards, zero-thickness cohesive elements with fracture criteria are embedded along grain boundaries and in grains as potential cracking paths. The closed single-line boundary Voronoi tessellation and generated meshes are imported into ABAQUS to conduct FEM calculations for simulating intergranular and transgranular fracture. One point should be emphasized that the zerothickness cohesive zone element formulation has been proved to be an efficient method for representing the intersection between grains before failure [26,27]. The Voronoi tessellation with double-line boundary will be considered in our future work. 2

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

2.2. Numerical simulations

Table 1 Calculated physical properties of NiTi alloys with B2 structure.

2.2.1. Molecular dynamics Usually, interactions between atoms contained in materials are represented through a potential function, and thus an appropriate functional form combined with physical parameters is needed when conducting molecular simulations for each material type. It is well known that the embedded-atom method (EAM) can be used to describe the interatomic potential applicable to many metals and their alloys [28]. In the EAM, the total energy (E ) of the crystal is given as

E=

Fi ( ¯i ) + i

1 2

ij (Rij ) i, j ( i )

Property Elastic constants (GPa) Lattice constant (Å) Coheisve energies (eV) Bulk modulus (GPa)

Fi ( ¯i ) + i

1 2

Sij

ij (Rij )

j ( i)

=

0

l

( rRe 1)

DFT simulations [37]

C11 C12 C44 a0 Ecoh

144 128 80 2.999 −5.05

162 [35] 129 [35] 35 [35] 3.015 [36] −5.02 [36]

183 146 46 3.019 −5.15

B

134

140 [35]

158

cohesive energies of the B2-NiTi compound phase. A comparison between computational results and experimental data and first-principle calculations demonstrates that the MEAM potential could be used in our subsequent simulations. A three-dimensional (3D) approach was employed to study the crack propagation under Mode-I loading condition. Importantly, two types of mechanical models are considered for the modeling. One is for intergranular facture, the characteristic of which is symmetric about the grain boundary of two adjacent crystals. The other model consists of one single crystal for transgranular fracture. Zeng et al. [20] reported that the size of MD model could be determined through the geometrical similarity. Thus, an initial model was set as the B2 structure with the length of 100 cells, the width of 100 cells, and the thickness of 5 cells. In other words, all the MD models have an entire dimension of 425 Å ×425 Å × 21.2 Å, a total of 100,536 atoms. The crystal direction was along the orientation of X [1 −1 0], Y [1 1 0], Z [0 0 1]. An initial crack being located at the left side was inserted in the model edge, as seen in Fig. 3, and the relevant dimensions were presented as well. In addition, four models with different micro-defect patterns are considered for each type when performing MD simulations. To restrain boundary effects, the crystal right side was fixed for motions along the x-direction to be avoided. It is noteworthy that Hu et al. [38] argued that strain rate effect became insignificant when the strain rate was beyond 109 S−1 and thus the corresponding loading could be regarded as quasi-static loading. For the purpose of analytical simplicity, the uniaxial tension load was applied as a pre-determined displacement rate with 0.45 Å/ps along Y axis, which was determined from the simulated strain rate on the order of 109 S−1. Besides, the time step was 0.001 ps with a periodic boundary condition along z direction. Meanwhile, the aperiodic boundary conditions were applied along both Xand Y-directions. All simulations were performed by using a microcanonical (NVE) ensemble and there were 60,000 steps in the whole calculation procedure. The purpose of uniform stretching is to avoid the shock wave during the simulations.

(2)

where Sij is the screening function accordingly. It is noteworthy that ¯i could be determined by integrating several partial electron density terms for different angular contributions with weighting factors. Indeed, each partial electron density is taken as a function of atomic configuration and atomic electron density. The latter is written as [30] a (R )

Experiments

Note: DFT, Density functional theory.

(1)

where ij (Rij ) refers to the pair interaction between atoms i and j separated by a distance Rij . Fi is the embedding energy as a function of background electron density ¯i . In fact, a modified embedded-atom method (MEAM) combined with an angular dependent term was proposed by Baskes [29] and Ko et al. [30]. In particular, what the advantage this potential has are that it contains additional angular terms accounting for the bonding character and allows for characterizing a wide range of phases using a common function in comparison with the EAM potential [31]. Then, the MEAM potential is expressed as follows:

E=

Present potential

(3)

where 0 is the atomic electron density scaling factor. l denotes the decay length. re represents the nearest-neighbor distance in the equilibrium reference structure. A detailed formulation of the MEAM formalism is also available in the literature [32,33]. Herein, the MEAM potential is employed to describe the interatomic interactions in this work. It is accepted that the B2-NiTi is a cubic structure with tetragonal cells, the lattice constants of which are a = a0 and b = c = 2 a0 [34]. Fig. 2 illustrates the atom positions in the unit cell for B2 structure. In order to verify the accuracy of above-mentioned MEAM potential, some physical properties of NiTi alloys are calculated with the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS). Table 1 presents the lattice parameters, elastic constants, and

2.2.2. Cohesive zone model It is well justified that the analysis of both crack onset and its ongoing propagation can be accomplished by finite element modeling combined with the cohesive-zone model, which assumes that crack surfaces carry tractions that resist normal separation and tangential sliding before chemical bonds break [21]. Thus, the traction is dependent upon the normal and tangential displacements. In the CZM, three typical characteristics are crucial for describing the cracking evolution at the front of crack, which are the initial stiffness, damage initiation criterion, and damage evolution law, respectively. It is suggested to be preferably characterized through the relationship between tractions (T) and separations (S), i.e., the T-S law, for damage evolution. There are two kinds of T-S law used to describe the properties of cohesive elements for brittle and ductile fractures. In general, the bilinear traction-separation law is used to analyze the brittle fracture behavior, just like as NiTi alloys [39], seen in Fig. 4. In fact, such a law is usually determined by three following

Fig. 2. B2-NiTi unit cell structure. 3

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

Fig. 3. MD model with (a) blunt crack, (b) blunt crack and void, (c) sharp crack, and (d) sharp crack and void.

According to previous work [40], the stiffness (K ) in FEA code could be determined by

K=

Nmax

(5)

ini

For Mode-I fracture, isotropic materials are assumed to have equal normal and tangential stiffness components [40]. Then, Eq. (6) is expanded to

K = Kn = K s = K t =

parameters: one is the displacement when cohesive element surfaces begin to separate ( ini ), and the other corresponds to the maximum interface traction (Nmax ) , and the last denotes the separation at the element failure ( fai ) . The interface fracture energy ( ) equals to the area enclosed by the T-S curve

1 Nmax × 2

fai

ini

(6)

where Kn , Ks , and Kt denote the stiffness coefficient along the normal direction, the first shear direction, and the second shear direction, respectively. To parameterize the T-S law, the T-S region in the front of pre-crack tip was divided into two vertical layers (up and down), as presented in Fig. 3. Since previous experiences [20,41] revealed that a reasonable TS area was approximately six times larger than the initial crack length, 360 Å seems to be an effective size of cohesive zone. Thus, the finite domain for achieving the T-S curve is chosen as 360 Å × 20 Å for the purpose of computational accuracy and efficiency. For each time step, the atom upper and lower positions were required to be preserved. Thus, the separation could be availably determined by the distance from lower to upper atoms. Moreover, the traction was achieved through the average atomic tensile stress of inclusive atoms. In fact, the stress at atomic scale is defined as the strength measurement of in(i ) teratomic interactions [42]. The atomic stress tensor ( , = x , y, z ) at an atom i is calculated as:

Fig. 4. Bilinear traction-separation law for brittle fracture.

=

Nmax

(4) 4

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

i

=

2

statistical character with average size of 330 um2 and its variance of 43. According to the standard specified in the ASTM E399 [45], the specimen size was 1.25 W = 500 μm in length and 1.2 W = 480 μm in width, where W = 400 μm. In addition, the pre-crack length was a = 0.5 W (200 μm), being the same as the path length embedded with cohesive elements. The uniaxial tension loading was applied at a constant displacement under plane strain conditions. The displacements P along the opposite direction were loaded on the reference points, being set in the center of two circle holes with the diameter of 0.25 W = 100 μm, namely, the PR-1 and PR-2. The material is assumed to be isotropic. The 2D model was meshed with the CPE4R elements and CPE3 cohesive elements. To guarantee meshing efficiency, the mesh size for FEM was 3 μm. In addition, the material properties are E = 67 GPa and v = 0.33 [46]. It needs to be mentioned that this work focuses on conducting a comparative investigation for intergranular and transgranular fracture of polycrystals of B2-NiTi alloys. In brief, the simulated procedure is outlined as follows: first, the Voronoi tessellation is employed to generate massive polycrystalline microstructures with different dimensions, but they are assumed to have uniform properties. Then, cohesive elements are placed along potential crack trajectories. After obtaining the characteristic parameters described by the T-S law, which is achieved by MD simulations, the cracking behavior would be reproduced by finite element analysis. Fig. 6 presents a flowchart for modeling the cracking behavior across the micro-, meso-, and macro-scales.

N

1

(i ) =

f (i , j ) r (i , j ) i j (i )

4 Ri3 3

(7a) (7b)

where i and Ri are the volume and the radius of atom i, respectively. Ri could be achieved through its lattice constant [43], thereby producing i . Besides, N is the number of atoms in a limited region around the atom within the cut-off distance of the MEAM potential. f (i , j ) represents the vector component form of the interaction force applied by atom j to i, whereas r (i, j) is the vector component form of the distance between the two atoms. Taking an average over the volume around the atom i within the potential cut-off distance, the average atomic stress tensor ¯ at the atom is given as [42]:

¯ (i ) =

1 N

N

(j )

(8)

j=1

Most interestingly, the T-S relation of cohesive elements is achieved with the help of MD simulations, and consequently the damage parameters are determined from the T-S law. This method might be called the atomic-based CZM, which is beneficial to useful recommendations for deepening the comprehension of crack propagation behavior. Afterwards, the maximum nominal stress criterion (MAXS) is chosen as the failure rule, as follows:

MAX

n

,

s

,

t

Nmax Smax Tmax

=1

3. Results and discussion

(9)

3.1. T-S parameter calculation

where Nmax , Smax , and Tmax represent the maximum stress along the normal direction, the first shear direction, and the second shear direction, respectively. n is defined as the nominal stress in the pure normal mode. s and t are the nominal stress in the first and second shear directions, respectively. This rule implied that crack initiation and propagation occur when n = Nmax is satisfied for Mode-I cracking. Thus, the CZM has been constructed to perform the crack propagation modeling.

To investigate the influence of initial crack on intergranular and transgranular fracture, four MD models with various micro-defect patterns for each type are constructed with the help of LAMMPS, which are characterized by (a) blunt crack, (b) blunt crack and void, (c) sharp crack, and (d) sharp crack and void. The model visualizations are accomplished through OVITO, which is a professional analysis system for atomistic and particle simulations, as showed in Figs. 7 and 8. All the TS curves derived from MD simulations are fitted by the bilinear T-S law. In addition, it is noteworthy that the tensile traction and separation are taken as the average atomic stress and average opening displacement of the T-S area, respectively. Fig. 7(a) exhibits a progressive cracking behavior for the intergranular fracture with blunt crack (Model A). Point ① indicates an initial state. After that, the stress is found to rapidly increase as the crack tip separation expands (point ②). During this period, this curve is

2.2.3. Finite element method The finite element calculation is conducted on the ABAQUS computing. It has been proved that both two-dimensional (2D) and 3D models produce consistent results [40,44]. For the purpose of analysis, a geometrical model for 2D CT specimen with polycrystalline structures under static loading was employed to simulate the cracking behavior, as shown in Fig. 5. It is observed that the grain diameter exhibited a

Fig. 5. Schematic of CT specimen of polycrystal B2-NiTi alloys for FEM modeling. 5

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

Fig. 6. Flow chart for performing FEM simulations using an atomic-based CZM.

which are 28.61, 25.22, 25.36, and 24.76 J/m2, respectively. Correspondingly, the transgranular-fracture energies are 23.95, 19.56, 24.84, and 23.30 J/m2, respectively. Interestingly, it is discovered that the addition of void could decrease the fracture energy, which is attributed to the fact that there are less atomic bonds when voids appears, leading to an easier probability for crack propagation. Also, the fracture energy with the condition of sharp crack is found to be lower in comparison with blunt crack, which is due to stress concentration.

increasing linearly, which means that the deformation is in elastic period and nearly no metallic bond breaks [39]. When the stress reaches the peak value (16.67 GPa), the initial crack begins to propagate at the crack tip (point ③). Then, the stress decreases gradually until the propagation ends (point ⑥). Evidently, the cracking occurs along the grain boundary. As observed in Fig. 7(b), Model B represents a blunt crack (the main defect) combined with a void (the secondary defect). It is evident that the presence of micro-void is favor of crack growth. Thus, the maximum tensile stress is reduced to 15.11 GPa in comparison with Model A. After reaching the peak value (point ③), the main defect propagation initiates till it is connected with the secondary defect, leading to a perforative crack (point ④). Subsequently, the crack tip move forward. Fig. 7(c) demonstrated the cracking of Model C, the pre-crack type of which is a sharp crack. Owing to stress concentration, the sharp crack is more easily beneficial to propagation compared to the blunt crack, and then the maximum stress decreases from 16.67 to 15.68 GPa. A void as the secondary defect is introduced to Model D. For this case, the maximum stress is 14.85 GPa, being less than that of Model C. It is suggested that a void could facilitate the crack growth more easily, which is attributed by the less requisite energy for potential cracking. Similarity, Fig. 8 exhibits the transgranular cracking under four conditions by MD simulations. The maximum stress is 18.67 GPa in Model A and the minimum value is 17.58 GPa in Model D. Compared with transgranular fracture, the maximum tensile stress for intergranular fracture is found to be lower, suggesting that the crack initiation is more easily to occur in the latter case. However, crack propagation is dominated by not only maximum tensile stress but also fracture energy [40]. As shown in the figures, the bilinear fitting for T-S curves of the simulated results is conducted to determine several critical parameters, like as the maximum traction, the maximum separation, and stiffness coefficients. Table 2 provides the cohesive element parameters achieved from MD simulations. According to Eq. (4), the fracture energies for four intergranular-facture cases are calculated availably,

3.2. Crack propagation simulation To analyze the crack growth under a higher length scale, FEM simulations are performed for CT specimens under isothermal Mode-I loading. The cohesive elements obeying the MAXS criterion are embedded in the pre-crack tip. CT models with polycrystal grains are employed to investigate the intergranular and transgranular fracture. Since cohesive element zones are considered as potential cracking paths, the crack tip begins to move forward along the potential path, leading to a new crack once a cohesive element satisfies the failure criterion. Taking Model A as an example, the intergranular cracking of the polycrystalline specimen is presented in Fig. 9, which also depicts the evolution of the reaction force with regard to time. It is clearly observed that the stress near the crack tip gradually increases as the loading continues at the initial stage. When the force reaches the maximum value, the top and bottom of cohesive elements begin to separate, indicating as a crack initiation. Subsequently, the cracking accumulates continuously till its propagation ends. Thus, a real-life cracking procedure, including crack initiation, growth, and ultimate failure, is reproduced through finite element analysis from a micromechanics perspective. For the comparison purpose, Fig. 10 reproduces the cracking of the specimen by embedding cohesive elements along grain boundaries and in grain. It is noteworthy that those cohesive element parameters, 6

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

Fig. 7. T-S curve derived by MD model with (a) blunt crack, (b) blunt crack and void, (c) sharp crack, and (d) sharp crack and void, fitted by the bilinear law, for the intergranular fracture modeling.

which are derived from grain boundaries and grains, respectively, are completely different. Interestingly, it is found that the transgranulardominated fracture occurs in the polycrystalline specimen, indicating

that the transgranular fracture is prior to intergranular fracture in such materials. Such a behavior was observed in previous study [47]. It is discovered from Figs. 9 and 10 that the relationship between the 7

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

Fig. 8. T-S curve derived by MD model with (a) blunt crack, (b) blunt crack and void, (c) sharp crack, and (d) sharp crack and void, fitted by the bilinear law, for the transgranular fracture modeling.

reaction force and time seems to be a similar trend. The stress profiles detected in the area of pre-crack tip demonstrate that the crack always propagates along the orientation where the stress concentration occurs.

In comparison with intergranular fracture, the crack path is found to be rougher for the transgranular fracture. Meanwhile, for the latter one, the main crack grows accompanied by several tiny cracks, and thus the 8

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

crack path is characterized by a zigzag. This is explained by the presence of microcracks that deflect the growth direction of main crack. Thus, it makes crack propagation more difficult, leading to improved fracture toughness.

Fracture toughness is a key property for describing the ability of a material to resist fracture, and it is also one of the most important properties of any materials for design applications [48]. In fact, the fracture toughness of linear-elastic materials can be determined from the stress intensity factor K at which a thin crack in materials begins to grow [49], which is expressed as follows [45]:

5.36 4.69 5.24 9.99

5.47 5.26 5.64 6.04

3.3. Fracture toughness prediction

K a W

2+

(1

a W

23.95 19.56 24.84 23.30

) 3

3 2

0.886 + 4.64 ×

5.60 ×

a W

a W

13.32 ×

a W

2

4

(7)

25.66 21.68 26.7 26.51

where a and W are consistent with Fig. 5. P is the reaction force loaded on the reference points (PR-1 and PR-2) along the vertical direction. With consideration that the thickness effect is neglected in the 2D plane strain condition case, the specimen thickness (B) is chosen as 1 μm. It is evident that the value of P is vital for calculating the stress intensity factor. In this work, the critical stress intensity factor (KQ) that reflects the actual stress intensity at which fracture takes place could be derived from the maximum reaction force [50]. The relationships between the applied displacement and the reaction force of the CT specimens are presented in Figs. 11 and 12. The vertical reaction force and displacement at the center of the circular hole (point PR-1) are extracted from the ABAQUS results. It is discovered that the time when the reaction force reaches the maximum denotes that the crack begins to propagate. Afterwards, the required external load gradually decreases as the cracking evolution continues. Compared to other schemes, the blunt crack model has the maximal reaction force for both intergranular and transgranular fracture, and it is also found that the maximum tensile stress determined by MD simulations is the largest for this case. Accordingly, the crack growth is more difficult to occur in Model A, compared as other models. Interestingly, the reaction force after reaching the peak exhibits a decreasing behavior with a slower manner in the polycrystalline model, which is explained by the existence of microcracks that probably hinder the main crack propagation [21]. The predicted critical stress intensity factors for intergranular and transgranular fracture are presented in Fig. 13(a) and (b). It is found that KQ ranges from 31.6 to 35.3 MPa· m for the former model, while the corresponding value of the latter model is in the range from 32.8 to 37.1 MPa· m . The calculated results indicates that the polycrystalline model seems to be more appropriate for simulating the cracking behavior of NiTi alloys, because this modeling approach can yield a closer fit to the experimental value of 38 ± 3 MPa· m [46]. It further indicates that the Voronoi tessellation is valid for describing polycrystalline structures.

3.48 3.85 3.55 1.76 18.67 18.05 18.61 17.58 Blunt crack Blunt crack and void Sharp crack Sharp crack and void

P B W

a + 14.72 × W

4. Conclusions

A B C D

16.67 15.11 15.68 14.85 Blunt crack Blunt crack and void Sharp crack Sharp crack and void

In this study, a novel numerical approach has been presented to investigate the crack propagation of NiTi alloys across micro-, meso-, and macro-scales. At the atomic scale, MD simulations combined with an effective interatomic potential are conducted to derive the T-S law, and the resulting characteristic parameters are incorporated into cohesive elements along grain boundaries and in grains. Importantly, the Voronoi tessellation is employed to represent the grain distribution in

Transgranular fracture

Intergranular fracture

A B C D

Maximum tensile stress, Nmax (GPa)

3.05 2.87 2.78 2.46

34.33 32.17 33.57 33.35

28.61 25.22 25.36 24.76

=

Micro-defect type Example

Table 2 CZM parameters derived from MD simulations.

Separation at Nmax, δini (Å)

Maximum separation, δfai (Å)

Fracture energy, Г (J/m2)

Stiffness coefficients, Kn, Kt, Ks (×107 uN/um3)

M. Lu, et al.

9

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

Fig. 9. The intergranular-cracking evolution of polycrystal specimen with micro-defect of Model A.

Fig. 10. The transgranular-cracking evolution of polycrystal specimen with micro-defect of Model A.

Fig. 11. The response of reaction force and displacement for polycrystal specimen with various micro-defect types in the intergranular mode.

Fig. 12. The response of reaction force and displacement for polycrystal specimen with various micro-defect types in the transgranular mode.

10

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

Acknowledgements The authors gratefully acknowledge the financial support of the Natural Science Foundation Project of CQ CSTC under grant number cstc2018jcyjAX0581, the Fundamental Research Funds for the Central Universities under grant number XDJK2018B002, the Venture & Innovation Support Program for Chongqing Overseas Returnees under grant number cx2018078, and the Science and Technology on Reactor System Design Technology Laboratory under grant number HT-KFKT10-2017001. References [1] X.M. Wang, Y.F. Wang, Z.F. Yue, Finite element simulation of the influence of TiC inclusions on the fatigue behavior of NiTi shape-memory alloys, Metall. Mater. Trans. A 36A (2005) 2615–2620. [2] M. Elahinia, N.S. Moghaddam, M.T. Andani, A. Amerinatanzi, B.A. Bimber, R.F. Hamilton, Fabrication of NiTi through additive manufacturing: a review, Prog. Mater. Sci. 83 (2016) 630–663. [3] S.W. Robertson, A.R. Pelton, R.O. Ritchie, Mechanical fatigue and fracture of Nitinol, Int. Mater. Rev. 57 (2012) 1–36. [4] Z.D. Sha, S.S. Quek, Q.X. Pei, Z.S. Liu, T.J. Wang, V.B. Shenoy, Y.W. Zhang, Inverse pseudo Hall-Petch relation in polycrystalline graphene, Sci. Rep. 4 (2014) 5991. [5] G. Rajasekaran, P. Narayanan, A. Parashar, Effect of point and line defects on mechanical and thermal properties of graphene: a review, Crit. Rev. Solid State 41 (2016) 47–71. [6] H. Zhang, D.H. Fu, H.P. Song, Y.L. Kang, G.Y. Huang, G. Qi, J.Y. Li, Damage and fracture investigation of three-point bending notched sandstone beams by DIC and AE techniques, Rock Mech. Rock Eng. 48 (2015) 1297–1303. [7] S. Daly, A. Miller, G. Ravichandran, K. Bhattacharya, Experimental investigation of crack initiation in thin sheets of nitinol, Acta Mater. 55 (2007) 6322–6330. [8] S.W. Robertson, R.O. Ritchie, In vitro fatigue-crack growth and fracture toughness behavior of thin-walled superelastic Nitinol tube for endovascular stents: a basis for defining the effect of crack-like defects, Biomaterials 28 (2007) 700–709. [9] L.L. Meisner, M.G. Ostapenko, A.I. Lotkov, A.A. Neiman, Surface microstructure and B2 phase structural state induced in NiTi alloy by a high-current pulsed electron beam, Appl. Surf. Sci. 324 (2015) 44–52. [10] J.Q. Zhang, F. Wang, Modeling of progressive failure in ductile matrix composites including local matrix yielding, Mech. Adv. Mater. Struc. 16 (2009) 522–535. [11] X.W. Zhou, J.A. Zimmerman, E.D. Reedy, N.R. Moody, Molecular dynamics simulation based cohesive surface representation of mixed mode fracture, Mech. Mater. 40 (2008) 832–845. [12] X.W. Zhou, N.R. Moody, R.E. Jones, J.A. Zimmerman, E.D. Reedy, Molecular-dynamics-based cohesive zone law for brittle interfacial fracture under mixed loading conditions: effects of elastic constant mismatch, Acta Mater. 57 (2009) 4671–4686. [13] V. Tvergaard, J.W. Hutchinson, The influence of plasticity on mixed mode interface toughness, J. Mech. Phys. Solids 41 (1993) 1119–1135. [14] X.P. Xu, A. Needleman, Numerical simulations of fast crack growth in brittle solids, J. Mech. Phys. Solids 42 (1994) 1397–1434. [15] V. Yamakov, E. Saether, D.R. Phillips, E.H. Glaessgen, Molecular-dynamics simulation-based cohesive zone representation of intergranular fracture processes in aluminum, J. Mech. Phys. Solids 54 (2006) 1899–1928. [16] W. Zhu, L. Yang, J.W. Guo, Y.C. Zhou, C. Lu, Determination of interfacial adhesion energies of thermal barrier coatings by compression test combined with a cohesive zone finite element model, Int. J. Plasticity 64 (2015) 76–87. [17] P.A. Gustafson, A.M. Waas, The influence of adhesive constitutive parameters in cohesive zone finite element models of adhesively bonded joints, Int. J. Solids Struct. 46 (2009) 2201–2215. [18] T.T. Zhou, C.Z. Huang, Simulation of crack propagation in single phase ceramic tool materials, Comp. Mater. Sci. 104 (2015) 177–184. [19] R. Komanduri, N. Chandrasekaran, L.M. Raff, Molecular dynamics (MD) simulation of uniaxial tension of some single-crystal cubic metals at nanolevel, Int. J. Mech. Sci. 43 (2001) 2237–2260. [20] X.G. Zeng, T.X. Han, Y. Guo, F. Wang, Molecular dynamics modeling of crack propagation in titanium alloys by using an experiment-based Monte Carlo model, Eng. Fract. Mech. 190 (2018) 120–133. [21] T.T. Zhou, C.Z. Huang, H.L. Liu, J. Wang, B. Zou, H.T. Zhu, Crack propagation simulation in microstructure of ceramic tool materials, Comp. Mater. Sci. 54 (2012) 150–156. [22] F. Aurenhammer, Voronoi diagrams – a survey of a fundamental geometric data structure, Comput. Surv. 23 (1991) 345–405. [23] S. Fortune, A sweepline algorithm for Voronoi diagrams, Algorithmica 2 (1987) 153–174. [24] J. Fleig, The grain boundary impedance of random microstructures: numerical simulations and implications for the analysis of experimental data, Solid State Ionics 150 (2002) 181–193. [25] H.D. Espinosa, P.D. Zavattieri, A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. Part I: Theory and numerical implementation, Mech. Mater. 35 (2003) 333–364. [26] T. Diehl, On using a penalty-based cohesive-zone finite element approach. Part I: Elastic solution benchmarks, Int. J. Adhes. Adhes. 28 (2008) 237–255. [27] Z.J. Wu, X.Y. Xu, Q.S. Liu, Y.T. Yang, A zero-thickness cohesive element-based

Fig. 13. The calculated critical stress intensity factor of polycrystalline structures with various micro-defect types for (a) intergranular and (b) transgranular fracture.

materials from the mesoscopic scale, and in turn its statistical character is described quantitatively. Afterwards, FEM combined with parameterized CZM are performed to reproduce the macroscopic intergranular and transgranular fracture of CT specimens subjected to uniaxial tensile loading under plane strain conditions, by taking into account various micro-defect types. In addition, the fracture toughness of NiTi alloys is predicted availably. Thus, the relationship between the microstructure and mechanical property is established explosively. A comparison between the predicted values and experimental data demonstrates that the proposed model is more appropriate for studying the cracking evolution of NiTi alloys. The results are expected to contribute to a comprehensive understanding of the crack propagation from the multi-scale perspective. However, the fracture toughness and crack-growth resistance are affected by microstructures and loading conditions, such as grain size, grain boundary strength, microcracks, and strain rate. Most importantly, materials exhibit a fracture statistics, which is associated with grain size dependence [51]. Thus, a more precise calibration of the model, capable of explaining such effects on fracture behavior is the subject of our further work. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

11

Theoretical and Applied Fracture Mechanics 105 (2020) 102402

M. Lu, et al.

[28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]

numerical manifold method for rock mechanical behavior with micro-Voronoi grains, Eng. Anal. Bound Elem. 96 (2018) 94–108. M.S. Daw, M.I. Baskes, Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals, Phys. Rev. B 29 (1984) 6443–6453. M.I. Baskes, Modified embedded-atom potentials for cubic materials and impurities, Phys. Rev. B 46 (1992) 2727–2742. W.S. Ko, B. Grabowski, J. Neugebauer, Development and application of a Ni-Ti interatomic potential with high predictive accuracy of the martensitic phase transition, Phys. Rev. B 92 (2015) 134107. Y. Guo, X.G. Zeng, H.Y. Chen, T.X. Han, H.Y. Tian, F. Wang, Molecular dynamics modeling of the effect of nanotwins on the superelasticity of single-crystalline NiTi alloys, Adv. Mater. Sci. Eng. (2017) 7427039. B.J. Lee, M.I. Baskes, H. Kim, Y.K. Cho, Second nearest-neighbor modified embedded-atom-method potential for bcc transition metals, Phys. Rev. B 64 (2001) 184102. B.J. Lee, W.S. Ko, H.K. Kim, E.H. Kim, The modified embedded-atom method interatomic potentials and recent progress in atomistic simulations, Calphad 34 (2010) 510–522. Y. Zhong, K. Gall, T. Zhu, Atomistic characterization of pseudoelasticity and shape memory in NiTi nanopillars, Acta Mater. 60 (2012) 6301–6311. O. Mercier, K.N. Melton, G. Gremaud, J. Hagi, Single-crystal elastic constants of the equiatomic NiTi alloy near the martensitic transformation, J. Appl. Phys. 51 (1980) 1833–1834. Y. Zhong, K. Gall, T. Zhu, Atomistic study of nanotwins in NiTi shape memory alloys, J. Appl. Phys. 110 (2011) 033532. N. Hatcher, O.Y. Kontsevoi, A.J. Freeman, Role of elastic and shear stabilities in the martensitic transformation path of NiTi, Phys. Rev. B 80 (2009) 144203. T. Fu, X.H. Peng, C. Huang, S.Y. Weng, Y.B. Zhao, Z.C. Wang, N. Hu, Strain rate dependence of tension and compression behavior in nano-polycrystalline vanadium nitride, Ceram. Int. 43 (2017) 11635–11641. X.M. Wang, B.X. Xu, Z.F. Yue, X.Y. Tong, Fracture behavior of the compact tension specimens in NiTi shape memory alloys, Mat. Sci. Eng. A – Struct. 485 (2008)

14–19. [40] T.L. Xu, R. Stewart, J.H. Fan, X.G. Zeng, A.L. Yao, Bridging crack propagation at the atomistic and mesoscopic scale for BCC-Fe with hybrid multiscale methods, Eng. Fract. Mech. 155 (2016) 166–182. [41] Y. Sheng, X.G. Zeng, The deformation mechanisms in process of crack propagation for alpha titanium with compounding microdefects, Adv. Mater. Sci. Eng. (2016) 2156936. [42] M.F. Horstemeyer, M.I. Baskes, Atomistic finite deformation simulations: a discussion on length scale effects in relation to mechanical stresses, J. Eng. Mater. – T ASME 121 (1999) 114–119. [43] Y. He, Y.J. Li, C.F. Chen, H.B. Yu, Diffusion coefficient of hydrogen interstitial atom in α-Fe, γ-Fe and ε-Fe crystals by first-principle calculations, Int. J. Hydrogen Energ. 42 (2017) 27438–27445. [44] ABAQUS, ABAQUS/Standard User’s Manual, Hibbitt, Karlsson & Sorensen, Inc, 2001. [45] ASTM E399-90. Standard Test Method for Plane-Strain Fracture Toughness Testing of Metallic Materials. American Society for Testing and Materials, Philadelphia, 1990. [46] B. Haghgouyan, C. Hayrettin, T. Baxevanis, I. Karaman, D.C. Lagoudas, Fracture toughness of NiTi towards establishing standard test methods for phase transforming materials, Acta Mater. 162 (2019) 226–238. [47] D.A. Miller, W.R. Thissell, D.A.S. Macdougall, Dynamic tensile plasticity and damage evolution in shape-memory Ni-Ti, J. Phys. IV 10 (2000) 341–346. [48] A. Ahadi, Q.P. Sun, Grain size dependence of fracture toughness and crack-growth resistance of superelastic NiTi, Scr. Mater. 113 (2016) 171–175. [49] M.D. Wei, F. Dai, N.W. Xu, T. Zhao, Stress intensity factors and fracture process zones of ISRM-suggested chevron notched specimens for mode I fracture toughness testing of rocks, Eng. Fract. Mech. 168 (2016) 174–189. [50] R.L. Holtz, K. Sadananda, M.A. Imam, Fatigue thresholds of Ni-Ti alloy near the shape memory transition temperature, Int. J. Fatigue 21 (1999) S137–S145. [51] S. Ozaki, Y. Aoki, T. Osada, K. Takeo, W. Nakao, Finite element analysis of fracture statistics of ceramics: effect of grain size and pore size distributions, J. Am. Ceram. Soc. 101 (2018) 3191–3204.

12