PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of materialsLink to Publisher's site
 
Materials (Basel). 2016 December; 9(12): 977.
Published online 2016 December 2. doi:  10.3390/ma9120977
PMCID: PMC5456998

Modelling of Granular Fracture in Polycrystalline Materials Using Ordinary State-Based Peridynamics

Timon Rabczuk, Academic Editor

Abstract

An ordinary state-based peridynamic formulation is developed to analyse cubic polycrystalline materials for the first time in the literature. This new approach has the advantage that no constraint condition is imposed on material constants as opposed to bond-based peridynamic theory. The formulation is validated by first considering static analyses and comparing the displacement fields obtained from the finite element method and ordinary state-based peridynamics. Then, dynamic analysis is performed to investigate the effect of grain boundary strength, crystal size, and discretization size on fracture behaviour and fracture morphology.

Keywords: polycrystalline materials, ordinary state-based peridynamics, transgranular fracture, intergranular fracture

1. Introduction

Polycrystalline materials are widely used in many different industrial applications. Amongst the various existing polycrystalline materials, metals and ceramics are common examples. Polycrystalline materials are composed of individual crystals that have a particular crystal orientation and are separated from neighbouring crystals via grain boundaries. Microscopic features of polycrystalline materials such as crystal orientation, grain boundary strength, etc. may have a significant effect on the overall macroscopic behaviour of the material, especially on the fracturing behaviour of these materials. Therefore, it is essential to analyse this type of material at the microscopic scale. Although experimental approaches can be useful for this purpose [1,2,3], it is not always possible to perform such experiments and they can also be very expensive. On the other hand, numerical approaches can be a very good alternative. There are various numerical studies available in the literature. For the majority of these studies, cohesive zone elements (CZM) [4,5,6], extended finite element methodology (XFEM) [7,8], and boundary element method (BEM) [9,10] were utilized. Although these techniques are widely used and powerful, they may have limitations for some specific cases and conditions. Instead, a new continuum mechanics formulation, peridynamics (PD) [11], can be considered.

As opposed to partial differential equations that traditional approaches are based on, peridynamics utilizes integro-differential equations without containing any spatial derivatives. Hence, these equations are always applicable regardless of discontinuities such as cracks. Peridynamics has been used for the fracture analysis of many different types of materials and material behaviours [12,13,14,15,16,17,18,19]. It has also been applied for the analysis of polycrystalline materials [20,21,22]. However, these studies used either original bond-based formulation (BB) [11] or non-ordinary state-based (NOSB) [23] formulations. Bond-based formulation has limitations on material constants whereas non-ordinary state-based formulations may encounter the zero-energy mode problem. In order to overcome all these issues, an ordinary state-based (OSB) peridynamic formulation [23,24] can be utilized. The numerical solution of peridynamics is generally obtained by using a meshless scheme. Therefore, the formulation does not suffer from issues such as mesh distortion, especially for problems involving large deformations and fractures. There are also other mesh-free methods available in the literature and used for modelling shear bands in metals [25,26,27,28]: concrete fragmentation [29,30,31], dynamic fracture in thin shells [26,32], and fluid–structure interaction [33]. A good example of a meshless approach used for fracture is the cracking particles method (CPM) [29,34,35]. In CPM, the crack path is represented by a set of cracked particles. CPM is especially useful for complex fracture patterns such as crack branching and coalescence. Other promising techniques used for fracture modelling include smoothed particle hydrodynamics [36,37,38], molecular dynamics [39], the discrete element method [40], and the force potential-based particle method [41]. Although the aforementioned approaches may have certain advantages for particular conditions, in this study peridynamics was chosen for modelling granular fracture in polycrystalline materials. Moreover, the ordinary state-based formulation was utilized for the first time in the literature in order to overcome the limitations of bond-based formulation and zero-energy mode problem of non-ordinary state-based theory. After validating the formulation, several demonstration cases are considered to investigate the effect of grain boundary strength, crystal (grain) orientation, and grain size.

2. Ordinary State-Based Peridynamic Formulation for a Cubic Crystal

The equation of motion of Ordinary State-Based (OSB) peridynamics (PD) can be written as:

ρ(x)u¨(x,t)=H(t(uu,xx,t)t(uu,xx,t))dH+b(x,t),
(1)

where ρ(x) represents the density and u¨(x,t) is the acceleration of material point x at time t. Moreover, t(uu,xx,t) and t(uu,xx,t) denote the force density vectors of the material points x and x, and uu represents the difference of displacements of the material points x and x at time t. In Equation (1), H represents the peridynamic horizon that defines the range of interaction of a particular material point, as shown in Figure 1. In the peridynamic literature, the size of the horizon is usually represented by the symbol δ.

Figure 1
The horizon of the material points located at x and x and the peridynamic forces between them in ordinary state-based peridynamics.

Similar to the bond-based (BB) peridynamic model presented in De Meo et al. (2016) [22], the ordinary state-based model for a cubic crystal can be represented using two types of interactions (bonds), as shown in Figure 2.

Figure 2
Type 1 bonds (green dashed lines) and Type 2 bonds (red solid lines) for the OSB PD cubic crystal model for a grain orientation of φ=π4.

These are:

  1. Type 1 bonds (green dashed lines)—interactions along all directions (ϕ=0~2π),
  2. Type 2 bonds (red solid lines)—interactions along the directions of ϕ=14π,34π,54π,74π,

where ϕ represents the angle between the orientation of the bond and the crystal (grain) orientation. As an example, bonds within the horizon of a particular material point for a grain orientation of φ=π4 are shown in Figure 2.

According to OSB PD theory, the strain energy density of a material point can be written as [24]:

W(k)=aθ(k)2+bT1Hδ|xx|(|yy||xx|)2dH+bT2j=1Jδ|x(j)x(k)|(|y(j)y(k)||x(j)x(k)|)2V(j),
(2)

where J is the total number of material points within the family of material point x(k).

By using the strain energy density expression given in Equation (2), the peridynamic force densities t and t can be obtained as:

t(uu,xx,t)=12Ayy|yy|,
(3)

where

A=4adδ|xx|Λθ+4δ(bT1+μT2bT2)s
(4)

with

μT2={1Type-2 bonds0otherwise
(5)

and

t(uu,xx,t)=12Byy|yy|,
(6)

where

B=4adδ|xx|Λθ+4δ(bT1+μT2bT2)s.
(7)

In Equations (3) and (6), y and y represent the location of material points x and x after deformation, i.e., y=x+u and y=x+u (see Figure 1). The PD dilatation, θ, for a crystal can be expressed as:

θ(k)=dHδ|xx|(|yy||xx|)ΛdH,
(8)

and the parameter, Λ, is defined as:

Λ=(yy|yy|)(xx|xx|).
(9)

The stretch parameter s can be expressed as:

s=|yy||xx||xx|.
(10)

The PD material parameter a is associated with the deformation involving dilatation, θ(k). The remaining material parameters, bT1 and bT2, are associated with deformation of the bonds along the Type 1 and Type 2 bond directions, respectively, as shown in Figure 2. All PD material constants can be expressed in terms of material constants of a cubic crystal, Qij, from classical theory. The procedure for obtaining these relationships is presented in Section 3.

When the stretch, s(k)(j) between material points k and j exceeds a critical stretch value, sc, the interaction breaks and damage occurs. Hence, there will no longer be any interaction between these two particles. The critical stretch parameter (2D) can be expressed as [24]:

sc=Gc(6πμ+169π2(κ2μ))δ,
(11)

where μ represents the shear modulus and κ is the bulk modulus of the material. According to [42], the critical energy release rate Gc can be obtained from fracture toughness KIc, as:

Gc=KIc2Eplane stress,
(12)

where E is the Young’s modulus.

An “interface strength coefficient” is introduced by [20] to investigate various fracture modes of polycrystalline materials and is defined as:

β=scGBscGI,
(13)

where scGB and scGI denote the critical stretch of interactions that cross the grain boundary and the critical stretch of interactions that are located within the grain, respectively, i.e., GB represents the grain boundary and GI represents the grain interior.

3. Derivation of PD Parameters

The PD material parameters, a,d,bT1,and bT2, that appear in the force density vector-stretch relations for in-plane deformation of a cubic crystal can be related to the engineering constants by considering three different simple loading conditions:

  1. Simple shear: γ12=ζ;
  2. Uniaxial stretch in crystal orientation direction: ε11=ζ,ε22=0;
  3. Biaxial stretch: ε11=ζ,ε22=ζ.

3.1. First Loading Condition (Simple Shear γ12=ζ)

In the first loading condition, a constant simple shear strain is applied as shown in Figure 3 and the corresponding dilatation and strain energy density from classical continuum mechanics (CCM) can be expressed as:

θ(k)CCM=0
(14)

and

W(k)CCM=12Q44ζ2.
(15)

Figure 3
Simple shear loading condition (x: crystal orientation direction).

For this loading condition, the length of the relative position of material points y(j) and y(k) in the deformed state becomes:

|y(j)y(k)|=[1+(sin ϕ(j)(k)cos ϕ(j)(k))ζ]|x(j)x(k)|.
(16)

The PD dilatation term can be evaluated as:

θ(k)PD=dHδξ{[1+(sin ϕcos ϕ)ζ]ξξ}dH=0,
(17)

in which ξ=|x(j)x(k)|. As expected, there is no dilatation for this loading condition. Therefore, the strain energy density can be calculated as:

W(k)PD=a(0)+bT1Hδξ{[1+(sin ϕcos ϕ)ζ]ξξ}2dH+bT2j=1Jδ|x(j)x(k)|((sin ϕ(j)(k)cos ϕ(j)(k))ζ|x(j)x(k)|)2V(j)
(18)

or

W(k)PD=(πhδ4ζ212)bT1+(δζ24j=1Jξ(j)(k)V(j))bT2.
(19)

Equating expressions of strain energy density from CCM and OSB PD formulations, i.e., Equations (15) and (19), results in:

(πhδ412)bT1+(δ4j=1Jξ(j)(k)V(j))bT2=Q442.
(20)

3.2. Second Loading Condition (Uniaxial Stretch in Crystal Orientation Direction: ε11=ζ,ε22=0)

In the second loading condition, a constant strain is applied along the direction of crystal orientation (Figure 4), and the corresponding dilatation and strain energy density from CCM can be expressed as:

θ(k)CCM=ζ
(21)

and

W(k)CCM=12Q11ζ2.
(22)

Figure 4
Uniaxial loading condition in crystal orientation direction, x.

The length of the relative position of material points y(j) and y(k) in the deformed state becomes:

|y(j)y(k)|=[1+(cos2 ϕ(j)(k))ζ]|x(j)x(k)|.
(23)

Due to this deformation, the dilatation of PD can be evaluated as:

θ(k)PD=dHδξ{[1+(cos2 ϕ)ζ]ξξ}dH
(24)

or

θ(k)PD=πdhδ3ζ2.
(25)

By equating expressions of the dilatation term from CCM and PD formulations, i.e., Equations (21) and (25), results in:

d=2πhδ3.
(26)

The PD strain energy density for this loading condition can be calculated as:

W(k)PD=aζ2+bT1Hδξ{[1+(cos2 ϕ)ζ]ξξ}2dH+bT2j=1Jδ|x(j)x(k)|((cos2 ϕ(j)(k))ζ|x(j)x(k)|)2V(j)
(27)

or

W(k)PD=aζ2+(πhδ4ζ24)bT1+(δζ24j=1Jξ(j)(k)V(j))bT2.
(28)

Hence, by equating expressions of strain energy density from Equations (22) and (28), the following relationship can be obtained:

a+(πhδ44)bT1+(δ4j=1Jξ(j)(k)V(j))bT2=12Q11.
(29)

3.3. Third Loading Condition (Biaxial Stretch: ε11=ζ,ε22=ζ)

In the third loading condition, a constant strain is applied in all directions (Figure 5), and the corresponding dilatation and strain energy density from CCM can be expressed as:

θ(k)CCM=2ζ
(30)

and

W(k)CCM=(Q11+Q12)ζ2.
(31)

Figure 5
Biaxial stretch loading condition.

The length of the relative position under this loading condition can be evaluated as:

|y(j)y(k)|=[1+(cos2 ϕ(j)(k)+sin2 ϕ(j)(k))ζ]|x(j)x(k)|.
(32)

Hence, the dilatation term in PD formulation can be expressed as:

θ(k)PD=dHδξ{[1+ζ]ξξ}dH
(33)

or

θ(k)PD=πdhδ3ζ.
(34)

By equating the expressions of dilatation from both CCM and PD, i.e., Equations (30) and (34), the same expression given in Equation (26) can be obtained. Moreover, the PD strain energy density under this loading condition can be evaluated as:

W(k)PD=aζ2+bT1Hδξ{[1+(sin2 ϕ)ζ]ξξ}2dH+bT2δζ2(j=1J|x(j)x(k)|V(j)).
(35)

By equating Equations (31) and (35), a new relationship can be obtained, as:

Q11+Q12=4a+(2πhδ43)bT1+(δj=1Jξ(j)(k)V(j))bT2.
(36)

Hence, the OSB PD material parameters can be expressed in terms of engineering constants of CCM by utilizing the three relationships given in Equation (20), (29), and (36) as:

{a=12(Q12Q44)bT1=3(Q11Q12)πhδ4bT2=2Q44Q11+Q12δ(j=1J|x(j)nx(k)n|V(j))d=2πhδ3.
(37)

For BB PD, the parameter a associated with dilatation should vanish, leading to the constraint equation Q12=Q44, which is a limitation of BB PD in cubic polycrystal analysis.

4. Numerical Results and Discussion

In this section, the results generated from static and dynamic PD analyses are presented, and comparisons with finite element method (FEM) results are also provided. For static analysis (Section 4.1), a single cubic Niobium (Nb) crystal model is considered first (Section 4.2.1) and displacement fields of PD and FEM are compared. Then, a cubic Molybdenum (Mo) polycrystal model with 18 Voronoi grains is analysed (Section 4.2.2), and the PD and FEM displacement fields are compared. For the dynamic analysis, the influence of the discretization size and the interface strength coefficient (β) on the results is considered first (Section 4.3.1). Then, the influence of crystal size on fracture behaviour is investigated (Section 4.3.2).

4.1. Material Data

Two different materials are considered in this study: niobium (Nb) for single crystal static analysis, and molybdenum (Mo) for polycrystal static and dynamic analysis. According to [43], the local stiffness matrix of each individual crystal can be written as:

[C]=[C11C12C12000C12C11C12000C12C12C11000000C44000000C44000000C44].
(38)

However, for plane stress configuration, the material matrix given in Equation (38) can be written by using reduced stiffness matrix following the procedure given in [44] as

[Q]=[Q11Q120Q12Q11000Q44],
(39)

where Qij are the reduced stiffness coefficients and can be expressed in terms of the elements of the stiffness matrix, Cij as

Q11=C112C122C11Q12=C11C12C122C11Q44=C44.
(40)

Therefore, the material properties of Nb and Mo are given in Table 1 as shown below [44]:

Table 1
Material properties of niobium and molybdenum.

The fracture toughness of Mo is specified as KIC=24.2 MPam [45], which corresponds to a critical stretch value of 0.008127 for plane stress configuration.

4.2. Static Analysis

A cubic crystal model with a length of 152.4 μm and a width of 76.2 μm is considered and the number of particles along the horizontal and vertical directions is 240 and 120, respectively. The values of grid spacing and horizon are Δx=0.635 μm and δ=1.915 μm, respectively. A uniform discretization scheme is used throughout this study. However, non-uniform discretization can also be possible by using small grid sizes at critical regions such as interfaces and utilizing the “Dual-horizon peridynamics” concept, as introduced by Ren et al. (2016) [46,47]. A horizontal tension loading of, P, is applied on the right edge of the model, and the left edge is fully fixed as shown in Figure 6. The tension loading is specified as a body load and applied to a single layer of material points at the right edge of the model. The displacement constraint condition at the left edge is also imposed to a single layer of material points. To reach the steady-state condition and perform static analysis, an adaptive dynamic relaxation technique was utilized as described in Madenci and Oterkus (2014) [24].

Figure 6
Crystal model for static analysis.

4.2.1. Static Analysis of Nb Single Crystal

The tension loading applied on the right edge of the model is P=174.4 MPa. Figure 7 and Figure 8 show a comparison of the results obtained from FEM and PD analysis under plane stress conditions for crystal orientations of 0° and 45°, respectively. Particles located along the central x-axis and y-axis are selected and horizontal and vertical displacements are compared, respectively.

Figure 7
Comparison of displacements between FEM and PD analyses for Nb crystal for 0° orientation: (a) horizontal displacements for particles along the central x-axis; (b) vertical displacements for particles along the central y-axis.
Figure 8Figure 8
Comparison of displacements between FEM and PD analyses for Nb Crystal for 45° orientation: (a) horizontal displacements for particles along the central x-axis; (b) vertical displacements for particles along the central y-axis.

Based on the results presented in Figure 7 and Figure 8, good agreement is obtained between PD and FEM analyses. Therefore, it can be concluded that the OSB PD crystal model presented in this study can produce accurate results for different crystal orientations for a single crystal.

4.2.2. Static Analysis of Mo Polycrystal

In the second case, a polycrystal model with 18 randomly orientated grains is generated by using Voronoi tessellation. A uniform discretization is utilized. Depending on the location of the material point, corresponding grain orientation is determined. Hence, Type 2 bonds will exist in many different directions according to the random orientation of the crystals. The average crystal size is 645.16 μm2, and the amount of tension loading applied on the right edge is P=374.5 MPa. The layout of the polycrystal model is shown in Figure 9.

Figure 9
Model for the static analysis of Mo polycrystal, composed of 18 randomly orientated grains with respect to the x–y coordinate system located at the centre of the model.

Displacement distributions for both x and y directions are compared between FEM and PD solutions as shown in Figure 10 and Figure 11 and a good agreement is obtained between two solutions, which confirms that the present model can also accurately represent polycrystal material behaviour.

Figure 10
Displacement field comparison between FEM and PD analyses for Mo polycrystal (x-direction).
Figure 11Figure 11
Displacement field comparison between FEM and PD analyses for Mo polycrystal (y-direction).

4.3. Dynamic Analysis of Mo Polycrystals

For the dynamic analysis, a 5 mm by 5 mm square plate with randomly oriented grains is considered as shown in Figure 12. A horizontal velocity boundary condition of V=5 m/s is applied on both the left and right edges of the model. Three layers of virtual particles are placed along the left and right edges to impose this condition, as suggested in [24]. A no-fail zone is also imposed on virtual particles and their neighbouring particles in order to allow the load to be accurately transferred inside the plate. Two pre-existing cracks with a length of 0.4 mm are applied at the centre of the bottom and top edges, as shown in Figure 13. The time step size is specified as dt=0.05 ns and the total number of time steps is 100,000, i.e., the total simulation time of 5.0 μs. The study considers three different interface strength coefficient, β, values (0.5, 1.0, and 2.0), three different mesh sizes (74 × 74, 150 × 150, and 300 × 300) and three different total numbers of grains (25, 100, and 400; i.e., different crystal sizes).

Figure 12
Polycrystal model for dynamic analysis (100 grains).
Figure 13
Location of the cracks in the model for dynamic analysis.

4.3.1. Effect of PD Discretization Size and Interface Strength Coefficient (β)

The aim of this analysis is to investigate the effect of the peridynamic discretization size on the crack pattern predicted by PD model and the morphology of intergranular and transgranular fracture modes when changing the value of the interface strength coefficient, β. The horizon is specified as δ=3.015Δx, which means that it is controlled by changing the PD discretization (74 × 74 particles, 150 × 150 particles, and 300 × 300 particles). Moreover, three different interface strength coefficient, β, values are considered to investigate the intergranular and transgranular fracture modes of the polycrystal.

Figure 14, Figure 15 and Figure 16 show the fracture pattern of the polycrystal under plane stress configuration at five different times (1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs) for β=0.5 with 74 × 74 particles, 150 × 150 particles, and 300 × 300 particles, respectively.

Figure 14
Fracture pattern of polycrystal when β=0.5 with 74 × 74 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 15
Fracture pattern of polycrystal when β=0.5 with 150 × 150 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 16
Fracture pattern of polycrystal when β=0.5 with 300 × 300 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.

The results show that with an increasing total number of particles, the intergranular crack pattern can be predicted more accurately and in more detail. However, the simulation time will increase rapidly as well. Therefore, it is important to find a good balance between accuracy and time. In this study, 150 × 150 particles can provide appropriate results, which is the reason why most of the simulations in this paper are chosen by using this number of particles.

Figure 17, Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22 show the fracture patterns of the polycrystal at five different times (1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs) for β=1.0 and β=2.0 with 74 × 74 particles, 150 × 150 particles, and 300 × 300 particles, respectively.

Figure 17
Fracture pattern of polycrystal when β=1.0 with 74 × 74 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 18
Fracture pattern of polycrystal when β=1.0 with 150 × 150 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 19
Fracture pattern of polycrystal when β=1.0 with 300 × 300 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 20
Fracture pattern of polycrystal when β=2.0 with 74 × 74 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 21
Fracture pattern of polycrystal when β=2.0 with 150 × 150 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 22
Fracture pattern of polycrystal when β=2.0 with 300 × 300 particles. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.

As described above, similar conclusions can be found in these simulations. For instance, branching of cracks can be obtained more clearly by increasing the total number of particles, but the simulations become more time-consuming. Moreover, the transgranular fracture mode becomes more dominant as the interface strength coefficient increases.

4.3.2. Effect of the Crystal Size

The aim of this section is to investigate the effect of the crystal size on fracture pattern. The plate is discretized by 150 × 150 particles, containing three different numbers of randomly orientated grains (25 grains, 100 grains, and 400 grains). Three different grain boundary strength coefficients, β=0.5, β=1.0 and β=2.0, are considered to investigate the effect of crystal size for different fracture modes.

Figure 23, Figure 24 and Figure 25 show the fracture pattern of the polycrystal at five different times (1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs) for β=0.5 with 25 grains, 100 grains, and 400 grains, respectively.

Figure 23
Fracture pattern of polycrystal when β=0.5 with 25 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 24
Fracture pattern of polycrystal when β=0.5 with 100 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 25
Fracture pattern of polycrystal when β=0.5 with 400 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.

According to the damage plots shown in Figure 23 with 25 grains at 2.0 μs, the propagation does not always occur from pre-existing cracks. Only the top pre-existing crack propagates in the 100-grain model and both pre-existing cracks propagate in the 400-grain model. This is because with an increase in the total number of grains, the probability of the pre-existing cracks being located on a grain boundary increases. In other words, since the grain boundary strength β=0.5 promotes intergranular fracture mode, the crack can more easily propagate if it is located on the grain boundary. However, for the grain boundary strength values of β=1.0 and β=2.0, there is no such difference observed in fracture behaviour and both pre-existing cracks propagate as shown in Figure 26, Figure 27, Figure 28, Figure 29, Figure 30 and Figure 31.

Figure 26
Fracture pattern of polycrystal when β=1.0 with 25 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 27
Fracture pattern of polycrystal when β=1.0 with 100 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 28
Fracture pattern of polycrystal when β=1.0 with 400 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 29
Fracture pattern of polycrystal when β=2.0 with 25 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 30
Fracture pattern of polycrystal when β=2.0 with 100 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.
Figure 31
Fracture pattern of polycrystal when β=2.0 with 400 grains in total. From left to right: time = 1.5 μs, 2.0 μs, 2.5 μs, 3.0 μs, and 3.5 μs, respectively.

5. Conclusions

In this paper, a new ordinary state-based peridynamic formulation is presented and related derivations are provided. The current model does not have any limitations on material constants as in the bond-based peridynamic theory. Static analyses were carried out for validation purposes and a comparison of results between PD and FEM shows that the proposed PD model can accurately capture the deformation behaviour of cubic polycrystalline materials. Then, dynamic analyses were carried out by considering different configurations to investigate the effect of interface strength coefficient, discretization size, and crystal size. The observations based on the evaluated results can be summarized as:

  1. Intergranular and transgranular fracture modes can be captured by changing the interface strength coefficient. As a future study, by comparing the experimental and PD results of crack morphology, actual interface strength coefficients can be estimated.
  2. The accuracy of simulation can be improved by increasing the total number of particles for intergranular fracture. However, the difference is not significant for transgranular fracture. In order to prevent the simulation from being time-consuming, a good balance should be considered between accuracy and simulation time.
  3. Pre-existing cracks can propagate more easily with decreasing crystal size for inter-granular fracture mode, since there is a higher probability of a pre-existing crack interacting with a grain boundary.

As a future study, experimental studies can be used to validate and refine the damage predictions of the proposed PD model. Moreover, as the current study is mainly focused on a 2D model, the formulation can be extended to a 3D model.

Acknowledgments

Results were obtained using an EPSRC-funded ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk). EPSRC grant no. EP/K000586/1.

Nomenclature

CZMcohesive zone model
FEMfinite element method
CCMclassical continuum mechanics
PDperidynamics
BBbond-based
OSBordinary state-based
NOSBnon-ordinary state-based
Hxhorizon of a generic particle x
δradius of the horizon [m]
fmechanical response function [N/m6]
cbond constant [N/m6]
sbond stretch
s0critical stretch
xvector defining the position of a generic particle x
xvector defining the position of a generic neighbour of particle x
yvector defining the position of particle x in the deformed configuration
yvector defining the position of particle x in the deformed configuration
b(x,t)body force density field [N/m3]
KICfracture toughness [MPam]
hplate’s thickness [m]
κbulk modulus [N/m2]
μshear modulus [N/m2]
Gccritical energy release rate [N/m]
[C]local stiffness matrix
[Q]reduced stiffness matrix
Δxgrid spacing [m]
bT1bond constant type-1 [N/m6]
bT2bond constant type-2 [N/m6]
ξijundeformed bond length between particles i and j [m]
ϕbond angle with respect to the crystal orientation angle [rad]
Vjvolume of a generic neighbouring particle j [m3]
u(x,t)displacement field at x [m]
u(x,t)displacement field at x [m]
ρ(x)mass density at x [Kg/m3]
u¨(x,t)acceleration vector field [m/s2]
βinterface strength coefficient

Author Contributions

Author Contributions

Ning Zhu, Dennj De Meo and Erkan Oterkus developed the formulation. Ning Zhu performed numerical simulations. Ning Zhu and Erkan Oterkus wrote the paper.

Conflicts of Interest

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Herbig M., King A., Reischig P., Proudhon H., Lauridsen E.M., Marrow J., Buffiere J.-Y., Ludwig W. 3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast X-ray tomography. Acta Mater. 2011;59:590–601. doi: 10.1016/j.actamat.2010.09.063. [Cross Ref]
2. Yan Y., Noufi R., Al-Jassim M. Grain-boundary physics in polycrystalline CuInSe2 revisited: Experiment and theory. Phys. Rev. Lett. 2006 doi: 10.1103/PhysRevLett.96.205501. [PubMed] [Cross Ref]
3. Gay P., Hirsch P., Kelly A. X-ray studies of polycrystalline metals deformed by rolling. III. The physical interpretation of the experimental results. Acta Cryst. 1954;7:41–50. doi: 10.1107/S0365110X54000060. [Cross Ref]
4. Sfantos G.K., Aliabadi M.H. A boundary cohesive grain element formulation for modelling intergranular microfracture in polycrystalline brittle materials. Int. J. Numer. Methods Eng. 2007;69:1590–1626. doi: 10.1002/nme.1831. [Cross Ref]
5. Kraft R., Molinari J., Ramesh K., Warner D. Computional micromechanics of dynamic compressive loading of a brittle polycrystalline material using a distribution of grain boundary properties. J. Mech. Phys. Solids. 2008;56:2618–2641. doi: 10.1016/j.jmps.2008.03.009. [Cross Ref]
6. Barut A., Guven I., Madenci E. A meshless grain element for micromechanical analysis with crystal plasticity. Int. J. Numer. Methods Eng. 2006;67:17–65. doi: 10.1002/nme.1623. [Cross Ref]
7. NSukumar N., Srolovitz D.J., Baker T.J., Prévost J.-H. Brittle fracture in polycrystalline microstructures with the extended finite element method. Int. J. Numer. Methods Eng. 2003;56:2015–2037. doi: 10.1002/nme.653. [Cross Ref]
8. Sukumar N., Srolovitz D. Finite element-based model for crack propagation in polycrystalline materials. Comput. Appl. Math. 2004;23:363–380.
9. Benedetti I., Aliabadi M. A three-dimensional cohesive-frictional grain-boundary micromechanical model for intergranular degradation and failure in polycrystalline materials. Comput. Methods Appl. Mech. Eng. 2013;265:36–62. doi: 10.1016/j.cma.2013.05.023. [Cross Ref]
10. Benedetti I., Aliabadi M. A three-dimensional grain boundary formulation for microstructural modeling of polycrystalline materials. Comput. Mater. Sci. 2012;67:249–260. doi: 10.1016/j.commatsci.2012.08.006. [Cross Ref]
11. Silling S. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids. 2000;48:175–209. doi: 10.1016/S0022-5096(99)00029-0. [Cross Ref]
12. Panchadhara R., Gordon P.A. Application of peridynamic stress intensity factors to dynamic fracture initiation and propagation. Int. J. Fract. 2016;201:81–96. doi: 10.1007/s10704-016-0124-8. [Cross Ref]
13. Madenci E., Oterkus S. Ordinary state-based peridynamics for plastic deformation according to von Mises yield criteria with isotropic hardening. J. Mech. Phys. Solids. 2016;86:192–219. doi: 10.1016/j.jmps.2015.09.016. [Cross Ref]
14. Oterkus S., Madenci E., Oterkus E., Hwang Y., Bae J., Han S. Hygro-thermo-mechanical analysis and failure prediction in electronic packages by using peridynamics; Proceedings of the 2014 IEEE 64th Electronic Components and Technology Conference (ECTC); Orlando, FL, USA. 27–30 May 2014; pp. 973–982.
15. Oterkus S., Madenci E. Peridynamics for antiplane shear and torsional deformations. J. Mech. Mater. Struct. 2015;10:167–193. doi: 10.2140/jomms.2015.10.167. [Cross Ref]
16. Diyaroglu C., Oterkus E., Oterkus S., Madenci E. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. 2015;69:152–168. doi: 10.1016/j.ijsolstr.2015.04.040. [Cross Ref]
17. Amani J., Oterkus E., Areias P., Zi G., Nguyen-Thoi T., Rabczuk T. A non-ordinary state-based peridynamics formulation for thermoplastic fracture. Int. J. Impact Eng. 2016;87:83–94. doi: 10.1016/j.ijimpeng.2015.06.019. [Cross Ref]
18. Oterkus E., Madenci E. Peridynamics for failure prediction in composites; Proceedings of the 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference; Honolulu, HI, USA. 23–26 April 2012; p. 1692.
19. Oterkus E., Barut A., Madenci E. Damage growth prediction from loaded composite fastener holes by using peridynamic theory; Proceedings of the 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference; Orlando, FL, USA. 12–15 April 2010.
20. Askari E., Bobaru F., Lehoucq R., Parks M., Silling S., Weckner O. Peridynamics for multiscale materials modeling. J. Phys. Conf. Ser. 2008;125:012078. doi: 10.1088/1742-6596/125/1/012078. [Cross Ref]
21. Ghajari M., Iannucci L., Curtis P. A peridynamic material model for the analysis of dynamic crack propagation in orthotropic media. Comput. Methods Appl. Mech. Eng. 2014;276:431–452. doi: 10.1016/j.cma.2014.04.002. [Cross Ref]
22. De Meo D., Zhu N., Oterkus E. Peridynamic modeling of granular fracture in polycrystalline materials. J. Eng. Mater. Technol. 2016;138:041008. doi: 10.1115/1.4033634. [Cross Ref]
23. Silling S.A., Epton M., Weckner O., Xu J., Askari E. Peridynamic states and constitutive modeling. J. Elast. 2007;88:151–184. doi: 10.1007/s10659-007-9125-1. [Cross Ref]
24. Madenci E., Oterkus E. Peridynamic Theory and Its Application. Springer; New York, NY, USA: 2014.
25. Li S., Hao W., Liu W.K. Mesh-free simulations of shear banding in large deformation. Int. J. Solids Struct. 2000;37:7185–7206. doi: 10.1016/S0020-7683(00)00195-5. [Cross Ref]
26. Li S., Hao W., Liu W.K. Numerical simulations of large deformation of thin shell structures using meshfree methods. Comput. Mech. 2000;25:102–116. doi: 10.1007/s004660050463. [Cross Ref]
27. Rabczuk T., Areias P.M.A., Belytschko T. A simplified meshfree method for shear bands with cohesive surfaces. Int. J. Numer. Methods Eng. 2007;69:993–1021. doi: 10.1002/nme.1797. [Cross Ref]
28. Rabczuk T., Samaniego E. Discontinuous modelling of shear bands using adaptive meshfree methods. Comput. Methods Appl. Mech. Eng. 2008;197:641–658. doi: 10.1016/j.cma.2007.08.027. [Cross Ref]
29. Rabczuk T., Belytschko T. A three dimensional large deformation meshfree method for arbitrary evolving cracks. Comput. Methods Appl. Mech. Eng. 2007;196:2777–2799. doi: 10.1016/j.cma.2006.06.020. [Cross Ref]
30. Wu Y., Wang D.D., Wu C.T. Three dimensional fragmentation simulation of concrete structures with a nodally regularized meshfree method. Theor. Appl. Fract. Mech. 2014;72:89–99. doi: 10.1016/j.tafmec.2014.04.006. [Cross Ref]
31. Chuzel-Marmot Y., Ortiz R., Combescure A. Three-dimensional SPH-FEM gluing for simulation of fast impacts on concrete slabs. Comput. Struct. 2011;89:2484–2494. doi: 10.1016/j.compstruc.2011.06.002. [Cross Ref]
32. Rabczuk T., Areias P.M.A., Belytschko T. A meshfree thin shell method for nonlinear dynamic fracture. Int. J. Numer. Methods Eng. 2007;72:524–548. doi: 10.1002/nme.2013. [Cross Ref]
33. Rabczuk T., Gracie R., Song J.H., Belytschko T. Immersed particle method for fluid-structure interaction. Int. J. Numer. Methods Eng. 2010;81:48–71. doi: 10.1002/nme.2670. [Cross Ref]
34. Rabczuk T., Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. Int. J. Numer. Methods Eng. 2004;61:2316–2343. doi: 10.1002/nme.1151. [Cross Ref]
35. Rabczuk T., Zi G., Bordas S., Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Comput. Methods Appl. Mech. Eng. 2010;199:2437–2455. doi: 10.1016/j.cma.2010.03.031. [Cross Ref]
36. Maurel B., Combescure A. An SPH shell formulation for plasticity and fracture analysis in explicit dynamics. Int. J. Numer. Methods Eng. 2008;76:949–971. doi: 10.1002/nme.2316. [Cross Ref]
37. Caleyron F., Combescure A., Faucher V., Potapov S. SPH modeling of fluid-solid interaction of dynamic failure analysis of fluid-filled thin shells. J. Fluids Struct. 2013;39:126–153. doi: 10.1016/j.jfluidstructs.2013.02.023. [Cross Ref]
38. Monaghan J.J. Smoothed particle hydrodynamics. Annu. Rev. Astron. Astrophys. 1992;30:543–574. doi: 10.1146/annurev.aa.30.090192.002551. [Cross Ref]
39. Krivtsov A. Molecular dynamics simulation of impact fracture in polycrystalline materials. Meccanica. 2003;38:61–70. doi: 10.1023/A:1022019401291. [Cross Ref]
40. Tavarez F.A., Plesha M.E. Discrete element method for modelling solid and particulate materials. Int. J. Numer. Meth. Eng. 2007;70:379–404. doi: 10.1002/nme.1881. [Cross Ref]
41. Brighenti R., Corbari N. Dynamic failure in brittle solids and granular matters: A force potential-based particle method. J. Numer. Meth. Eng. 2016;105:936–959. doi: 10.1002/nme.4998. [Cross Ref]
42. Anderson T. Fracture Mechanics—Fundamentals and Applications. 2nd ed. CRC Press; Boca Raton, FL, USA: 1995.
43. Hosford W.F. The Mechanics of Crystals and Textured Polycrystals. Oxford University Press; New York, NY, USA: 1993.
44. Kaw A.K. Mechanics of Composite Materials. 2nd ed. Taylor & Francis Group; London, UK: Boca Raton, FL, USA: 2006.
45. Sturm D., Heilmaier M., Schneibel J., Jehanno P., Skrotzki B., Saage H. The influence of silicon on the strength and fracture toughness of molybdenum. Mater. Sci. Eng. A. 2007;463:107–114. doi: 10.1016/j.msea.2006.07.153. [Cross Ref]
46. Ren H., Zhuang X., Cai Y., Rabczuk T. Dual-Horizon Peridynamics. Int. J. Numer. Methods Eng. 2016;108:1451–1476. doi: 10.1002/nme.5257. [Cross Ref]
47. Ren H., Zhuang X., Rabczuk T. A New Peridynamic Formulation with Shear Deformation for Elastic Solid. J. Micromech. Mol. Phys. 2016;1:1650009. doi: 10.1142/S2424913016500090. [Cross Ref]

Articles from Materials are provided here courtesy of Multidisciplinary Digital Publishing Institute (MDPI)