An Improved Ordinary State-Based Peridynamic Model for Granular Fractures in Cubic Crystals and the Effects of Crystal Orientations on Crack Propagation

Abstract

Material anisotropy caused by crystal orientation is an essential factor affecting the mechanical and fracture properties of crystal materials. This paper proposes an improved ordinary state-based peridynamic (OSB-PD) model to study the effect of arbitrary crystal orientation on the granular fracture in cubic crystals. Based on the periodicity of the equivalent elastic modulus of a cubic crystal, a periodic function regarding the crystal orientation is introduced into peridynamic material parameters, and a complete derivation process and expressions of correction factors are given. In addition, the derived parameters do not require additional coordinate transformation, simplifying the simulation process. Through convergence analysis, a regulating strategy to obtain the converged and accurate results of crack propagation paths is proposed. The effects of crystal orientations on crack initiation and propagation paths of single-crystal materials with different notch shapes (square, equilateral triangle, semi-circle) and sizes were studied. The results show that variations in crystal orientation can change the bifurcation, the number, and the propagation path direction of cracks. Under biaxial tensile loading, single crystals with semi-circular notches have the slowest crack initiation time and average propagation speed in most cases and are more resistant to fracture. Finally, the effects of grain anisotropy on dynamic fractures in polycrystalline materials under different grain boundary coefficients were studied. The decrease in grain anisotropy degree can reduce the microcracks in intergranular fracture and the crack propagation speed in transgranular fracture, respectively.

1. Introduction

Cubic crystalline structures have been widely found in engineering materials. Their macroscopic mechanical and fracture characteristics are significantly influenced by microscopic properties such as grain orientation, size, and grain boundary strength [1]. For example, because of the symmetry of the cubic crystal, its equivalent elastic modulus changes periodically with the crystal orientation [2]. The crack propagation path of single-crystal silicon is significantly affected by the angle between the low-index crystal plane and the direction of the notch front [3,4].

Polycrystalline materials are composed of many grains with different orientations and grain boundaries between adjacent grains. The grain orientation distribution can elucidate the physical mechanisms of microstructure formation [5]. A grain boundary is a two-dimensional planar defect where various atomic reactions and processes are promoted or accelerated [6]. As a result, according to the difference in grain boundary strength, the main fracture modes observed in polycrystalline materials are intergranular fracture and transgranular fracture [7]. In general, as the grain boundaries become weaker and more misorientated, the transition from transgranular to intergranular fracture is likely to occur, i.e., intergranular-type cracks are more likely to occur at the grain boundaries where there is a large difference between the orientations of two grains [8]. Therefore, it is vital to study the effects of crystal orientation and grain boundary characteristics on the mechanical and fracture properties of crystal materials.

The development of experimental techniques has led to a more in-depth understanding of the microstructure and properties of crystalline materials. However, in addition to the high cost caused by a large number of sample manufacturing and complex post-processing, many parameters during dynamic fracture are challenging to measure [9]. In order to provide a more effective method, many scholars have used different numerical methods to solve the fracture problems in crystalline materials. These include the cohesive zone model [10,11,12], extended finite element method [13,14,15], boundary element method [16,17,18], and phase field [1,19,20]. However, these methods are based on continuum mechanics theory and fundamentally have limitations in solving discontinuity problems [21].

Silling et al. [22] proposed a non-local theory, peridynamics (PD), which uses integral equations instead of partial differential equations in classical continuum mechanics (CCM) so it is not limited by material continuity and can spontaneously predict crack initiation and propagation. Recently, PD has been successfully applied to fracture analysis of polycrystalline materials. De Meo et al. [23] proposed a bond-state peridynamic (BB-PD) model for polycrystalline materials based on the cubic crystal. Similar to the definitions of matrix bonds and fibre bonds in the PD model of fibre-reinforced composite lamina [24], two types of bonds are defined in the direction of the relative position vector of a pair of material points in the BB-PD model for cubic crystals. One holds for all pairs of material points (Type-1), and the other is present only for the direction of some particular angles (Type-2), which is when the difference between the relative position vector direction of a pair of material points and the crystal orientation satisfies π/4 ± kπ/2 (k = 0,1,2,3). Subsequently, Li et al. [25] extended the application of the BB-PD model to the temperature field and studied the effects of grain size, grain boundary strength, and composite of materials. Lu et al. [21] investigated an application to the granular fracture behaviours of polycrystalline ice under dynamic loading conditions. However, the BB-PD model has limitations such as invariable Poisson’s ratios and does not support the constitutive model of plastic and incompressible materials [26]. Zhu et al. [27] developed an ordinary state-based peridynamics (OSB-PD) model to address the limitations of the BB-PD model and studied the effect of grain boundary strength, grain orientation, and grain size. Because of the arbitrariness of the orientation of polycrystalline materials, the direction of the relative position vector of all material point pairs easily mismatches the direction of the Type-2 bond in this type of model, making the model degenerate a model for isotropic materials.

Anisotropy due to crystal orientation can be introduced into other ways in the PD model. Based on the OSB-PD model for isotropic materials, Li et al. [28] used a non-spherical influence function to characterise grain anisotropy, but this model requires additional calibration of the auxiliary parameters. Zhang et al. [29] added a periodic function with a period of π/2 to the PD strain energy density. Liu et al. [30] introduced a period of π/3 into the force density function to develop a chiral-dependent PD model for graphene sheets. In the BB-PD model, micro-modulus (or bond constant) is one of the critical parameters regarding the elastic properties of materials, the PD material parameters in the OSB-PD model have the same nature [26].

Therefore, this work aims to improve the OSB-PD model based on cubic crystals and introduce a periodic function containing crystal orientation into the PD material parameters. In addition, the material stiffness matrix after the coordinate transformation based on the crystal orientation is used, which allows the derived PD material parameters to contain the orientation angle, so no additional coordinate transformation is required during the simulation process.

The structure of this paper is arranged as follows. Section 2 outlines the fundamentals of OSB-PD theory for isotropic elastic materials. Section 3 establishes the improved OSB-PD theory for cubic crystals with arbitrary orientations and gives the derivations of PD material parameters, surface and volume correction factors, and the damage model. Section 4 uses the finite element method to verify the proposed PD model at first and then carries out the m-convergence and δ-convergence analysis on the static tensile and dynamic fracture problems. Section 5 presents the applications of the proposed PD model to the fracture behaviours of single crystals and polycrystalline materials, involving the effects of notch shape, size, and anisotropy degree. Section 6 summarises the main contributions of this work.

2. Fundamentals of the OSB-PD Theory

Peridynamics is a non-local theory of continuum mechanics, in which the mechanical responses of one material point ${\mathit{x}}_{(k)}$ in an object are related to other material points ${\mathit{x}}_{(j)}$ inside its neighbouring region $H_{\mathit{x}}$, called the horizon. In the two-dimensional plane problem, the horizon domain is a disk, and the horizon size is defined by a radius $\delta$, as shown in Figure 1. The equation of motion at ${\mathit{x}}_{(k)}$ in the PD model is expressed by:

$$\rho \ddot{\mathit{u}}({\mathit{x}}_{(k)},t)={\int }_{H_{{\mathit{x}}_{(k)}}}\left(\mathit{t}\left({\mathit{u}}_{\left(j\right)}-{\mathit{u}}_{\left(k\right)},{\mathit{x}}_{\left(j\right)}-{\mathit{x}}_{\left(k\right)},t\right)-\mathit{t}({\mathit{u}}_{\left(k\right)}-{\mathit{u}}_{\left(j\right)},{\mathit{x}}_{\left(k\right)}-{\mathit{x}}_{(j)},t)\right)d{H}_{{\mathit{x}}_{(k)}}+\mathit{b}({\mathit{x}}_{(k)},t),$$

(1)

where $\mathit{b}$ is the external force density; $\rho$ is the material density; and $\mathit{u}$ is the displacement vector. The connection between ${\mathit{x}}_{(k)}$ and ${\mathit{x}}_{(j)}$ is called the PD bond, so the vector ${\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}$ is the bond vector. $\mathit{t}\left({\mathit{u}}_{\left(j\right)}-{\mathit{u}}_{\left(k\right)},{\mathit{x}}_{\left(j\right)}-{\mathit{x}}_{\left(k\right)},t\right)$, or $\mathit{t}$, is the PD force density vector from ${\mathit{x}}_{(j)}$ to ${\mathit{x}}_{(k)}$ in the horizon domain of ${\mathit{x}}_{(k)}$, and $\mathit{t}({\mathit{u}}_{\left(k\right)}-{\mathit{u}}_{\left(j\right)},{\mathit{x}}_{\left(k\right)}-{\mathit{x}}_{(j)},t)$, or ${\mathit{t}}^{\prime }$, is the one from ${\mathit{x}}_{(k)}$ to ${\mathit{x}}_{(j)}$ in the horizon domain of ${\mathit{x}}_{(j)}$.

In the OSB-PD model, the two pairwise force density vectors $\mathit{t}$ and ${\mathit{t}}^{\prime }$ are parallel and opposite in direction but different in magnitude. Additionally, the force density vector is parallel to the deformed bond vector ${\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}$, where ${\mathit{y}}_{(k)}$ and ${\mathit{y}}_{(j)}$ are the positions of ${\mathit{x}}_{(k)}$ and ${\mathit{x}}_{(j)}$ in the deformed configuration, and $\mathit{y}=\mathit{x}+\mathit{u}$. Thus, the two pairwise force density vectors can be calculated by

$$\mathit{t}({\mathit{u}}_{(j)}-{\mathit{u}}_{(k)},{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)=\frac{1}{2}A\frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|}$$

(2)

$$\mathit{t}({\mathit{u}}_{(k)}-{\mathit{u}}_{(j)},{\mathit{x}}_{(k)}-{\mathit{x}}_{(j)},t)=-\frac{1}{2}B\frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|}$$

(3)

where the symbol ‘|…|’ is the magnitude or norm of the vector,

$$A=4ad\frac{\delta }{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}{\Lambda }_{(k)(j)}{\theta }_{(k)}+4b\delta s_{(k)(j)}$$

(4)

$$B=4ad\frac{\delta }{\left|{\mathit{x}}_{(k)}-{\mathit{x}}_{(j)}\right|}{\Lambda }_{(j)(k)}{\theta }_{(j)}+4b\delta s_{(j)(k)}$$

(5)

where ${\theta }_{(k)}$ is the volume dilation. Normally, ${\theta }_{(k)}$ and ${\theta }_{(j)}$ are not equal.

$${\theta }_{(k)}=d{\int }_H\frac{\delta }{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}(\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|-\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|)dH$$

(6)

where ${\Lambda }_{(k)(j)}$ is the cosine of the angle between the relative position vectors of the initial configuration and the deformed configuration. ${\Lambda }_{(k)(j)}$ and ${\Lambda }_{(j)(k)}$ are equal in magnitude.

$${\Lambda }_{(k)(j)}=\frac{{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}}{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}·\frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|}$$

(7)

where $s_{(k)(j)}$ is the stretch ratio, which is similar to the elastic strain in CCM theory.

$$s_{(k)(j)}=\frac{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|-\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}$$

(8)

For isotropic materials, if the temperature is neglected, the strain energy density at point ${\mathit{x}}_{(k)}$ is described by

$$W_{(k)}=a{\theta }_{(k)}^2+b{\int }_H{\omega }_{(k)(j)}(\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|-\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|{)}^{2}dH$$

(9)

where ${\omega }_{(k)(j)}$ is the influence function

$${\omega }_{(k)(j)}=\frac{\delta }{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}$$

(10)

The variables $a,b,\mathrm{and} d$ are the PD material parameters to be solved, which can be obtained by comparing the strain energy densities under the simple loadings in PD to those in CCM.

3. OSB-PD Theory for Cubic Crystals and PD Material

3.1. Coordinate Transformation

Because of the crystal orientation, the material stiffness matrix of the cubic crystals in the global coordinate system changes with the orientation angle. In CCM theory, the Voigt notations of the material stiffness matrix for cubic crystals are described below.

$$\mathit{C}=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{12}& 0& 0& 0\\ C_{12}& C_{11}& C_{12}& 0& 0& 0\\ C_{12}& C_{12}& C_{11}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& C_{44}\end{array}\right]$$

(11)

In the assumption of plane stress, the material stiffness matrix $\mathit{C}$ can be simplified by

$$\mathit{Q}=\left[\begin{array}{ccc}{Q}_{11}& Q_{12}& 0\\ Q_{12}& Q_{11}& 0\\ 0& 0& Q_{44}\end{array}\right]$$

(12)

where $Q_{11}=\frac{C_{11}^2-C_{12}^2}{C_{11}}$, $Q_{12}=\frac{C_{11}{C}_{12}-C_{12}^2}{C_{11}},$ and $Q_{44}=C_{44}$.

When considering a crystal orientation angle $\gamma$, there is a transformation between the crystal coordinate system and the global coordinate system, as shown in Equation (13) and Figure 2,

$${\left[\begin{array}{c}X\\ Y\end{array}\right]}^{global}=\left[\begin{array}{cc}\mathit{cos}\gamma & -\mathit{sin}\gamma \\ \mathit{sin}\gamma & \mathit{cos}\gamma \end{array}\right]{\left[\begin{array}{c}x\\ y\end{array}\right]}^{crystal}$$

(13)

so that in the global coordinate system, the reduced material stiffness matrix with an orientation angle $\gamma$ after the coordinate transformation is calculated as

$$\mathit{R}=\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{14}\\ R_{12}& R_{11}& R_{24}\\ R_{14}& R_{24}& R_{44}\end{array}\right]=\left[\begin{array}{ccc}\frac{C_{11}{D}_{11}-C_{12}^2}{C_{11}}& \frac{C_{11}{D}_{12}-C_{12}^2}{C_{11}}& D_{14}\\ \frac{C_{11}{D}_{12}-C_{12}^2}{C_{11}}& \frac{C_{11}{D}_{11}-C_{12}^2}{C_{11}}& D_{24}\\ D_{14}& D_{24}& D_{44}\end{array}\right]$$

(14)

where

$$\left\{\begin{array}{l}{D}_{11}=C_{11}-2\mu {\mathit{sin}}^2\gamma {\mathit{cos}}^2\gamma \\ D_{12}=C_{12}+2\mu {\mathit{sin}}^2\gamma {\mathit{cos}}^2\gamma \\ D_{14}=\mu ({\mathit{sin}}^3\gamma \mathit{cos}\gamma -\mathit{sin}\gamma {\mathit{cos}}^3\gamma )\\ D_{24}=\mu (\mathit{sin}\gamma {\mathit{cos}}^3\gamma -{\mathit{sin}}^3\gamma \mathit{cos}\gamma )\\ D_{44}=C_{44}+2\mu {\mathit{sin}}^2\gamma {\mathit{cos}}^2\gamma \\ \mu =C_{11}-C_{12}-2C_{44}\end{array}\right.$$

(15)

Similarly, in the condition of plane strain, the reduced matrix with a crystal orientation $\gamma$ is

$$\mathit{T}=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{14}\\ T_{12}& T_{11}& T_{24}\\ T_{14}& T_{24}& T_{44}\end{array}\right]=\left[\begin{array}{ccc}{D}_{11}& D_{12}& D_{14}\\ D_{12}& D_{11}& D_{24}\\ D_{14}& D_{24}& D_{44}\end{array}\right]$$

(16)

The matrices $\mathit{R}$ and $\mathit{T}$ are used to derive the PD material parameters under the assumptions of plane stress and plane strain, respectively.

3.2. PD Material Parameters and Force Density Vectors

In BB-PD, the micro-modulus is an intrinsic parameter of the material that represents the stiffness of the bond between a pair of material points and is directly related to elastic constants in CCM theory [31]. In the OSB-PD model for isotropic materials, d is independent of the material type. a and b are proportional to the elastic constants of materials, but a is relevant to the volume dilatation. Thus, b is the critical PD parameter of the associated intrinsic material properties, similar to the micro-modulus.

For cubic crystals, the elastic properties are different depending on the axial angle $\phi$ in the crystal coordinate system, as shown in Figure 2b. In the two-dimensional cases, the equivalent elastic modulus $E_{\phi }$ at any direction of the crystal axis is calculated by

$$1/E_{\phi }=s_{11}+(2s_{12}-2s_{11}+s_{44})({\mathit{sin}}^2\phi {\mathit{cos}}^2\phi )$$

(17)

where $s_{11}=\frac{C_{11}+C_{12}}{C_{11}^2+C_{11}{C}_{12}-2C_{12}^2}$, $s_{12}=\frac{-C_{12}}{C_{11}^2+C_{11}{C}_{12}-2C_{12}^2},$ and $s_{44}=\frac{1}{C_{44}}$ are the components of the compliance matrix of cubic crystals [2].

It is evident in Equation (17) that the equivalent elastic modulus for cubic crystals has a period of π/2. In order to make the PD material parameters match the corresponding anisotropic material parameters at an arbitrary crystal orientation and axial angle of the crystal, the idea is to correlate the PD material parameter b with the axial angle of the crystal $\phi$. Assume that b depends on $\phi$ and has a period of π/2, by comparison with Equation (17), it is constructed by

$$b=b_1+b_2{\mathit{sin}}^2\phi {\mathit{cos}}^2\phi$$

(18)

By substituting it into the equation of the strain energy density at ${\mathit{x}}_{(k)}$ (Equation (9)), the following is obtained

$$W_{(k)}=a{\theta }_{(k)}^2+\left(b_1+b_2{\mathit{sin}}^2\phi {\mathit{cos}}^2\phi \right){\int }_H{\omega }_{(k)(j)}(\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|-\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|{)}^{2}dH$$

(19)

Similar to isotropic materials, the PD material parameters $a,b_1,b_2,\mathrm{a}\mathrm{n}\mathrm{d} d$ for cubic crystals can be solved by comparing the strain energy density under simple loading in PD and CCM theory. The derivation process is shown in Appendix A for the sake of brevity. Then, the solutions of the four PD material parameters are

$$\left\{\begin{array}{l}a=\frac{1}{2}\left(R_{12}-R_{44}\right)\\ b_1=\frac{3}{2\pi h{\delta }^{4}k}\left[\left(10\pi k_1+3\pi k_3\right)(R_{11}-R_{12})-2\left(14\pi k_1+\pi k_3\right)R_{44}\right]\\ b_2=\frac{24\left[2R_{44}-(R_{11}-R_{12})\right]}{h{\delta }^{4}k}\\ d=\frac{2}{\pi h{\delta }^3}\end{array}\right.$$

(20)

where $k=\pi k_3-2\pi k_1$, $k_1={\mathit{sin}}^2\gamma {\mathit{cos}}^2\gamma ,$ and $k_3={\mathit{cos}}^4\gamma +{\mathit{sin}}^4\gamma -4{\mathit{sin}}^2\gamma {\mathit{cos}}^2\gamma$.

For polycrystals, if two material points ${\mathit{x}}_{(k)}$ and ${\mathit{x}}_{(j)}$ belong to different grains, the equivalent bond constant $b_{\alpha }$ are calculated by the weighting method.

$$\frac{2}{b_{\alpha }}=\frac{1}{b_{\alpha (k)}}+\frac{1}{b_{\alpha (j)}}(\alpha =\mathrm{1,2})$$

(21)

Thus, the PD force density vectors of the interaction between a pair of material points ${\mathit{x}}_{(k)}$ and ${\mathit{x}}_{(j)}$ in the equation of motion (Equation (1)) can be summed as follows.

$$\mathit{f}={\mathit{t}}_{\left(k\right)\left(j\right)}-{\mathit{t}}_{\left(j\right)\left(k\right)}=\frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|}\left[2ad\frac{\delta }{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}{\Lambda }_{\left(k\right)\left(j\right)}\left({\theta }_{(k)}+{\theta }_{(j)}\right)+4\delta s_{\left(k\right)\left(j\right)}\left(b_1+b_2{\mathit{sin}}^2\phi {\mathit{cos}}^2\phi \right)\right]$$

(22)

3.3. Surface Correction and Volume Correction

For uniform discretisation, since the horizon domain is a circle for the two-dimensional plane problem, the material points located at the horizon boundary will have a part of its volume beyond the boundary. The volume correction coefficient ${\nu }_{(j)}$ should be introduced to reduce the error [26].

$${\nu }_{(j)}=\left\{\begin{array}{ll}1& \left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|\le (\delta -\varpi )\\ (\delta +\varpi -\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|)/2\varpi & (\delta -\varpi )<\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|\le \delta \\ 0& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$$

(23)

where $\varpi$ is half of the discrete spacing of material points, $\varpi =∆x/2$.

The same issue arises with the material points near the surface/boundary of the computational domain. The integration of the strain energy density in Equation (19) is performed for material points within a complete horizon domain. For a material point ${\mathit{x}}_{(k)}$ close to the boundary, it lacks part of its family points because the horizon domain is incomplete, causing the inaccuracy of the simulation results near the boundary. Therefore, the surface correction factors are introduced to deal with this problem when calculating the volume dilatation and the force density vector. The detailed derivation process is shown in Appendix B for the sake of brevity.

The modified volume dilatation and force density vector at point ${\mathit{x}}_{(k)}$ are expressed below.

$${\theta }_{(k)}=d\sum _{j=1}^N{\mathit{G}}_{(d)(k)(j)}\frac{\delta }{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}\left[\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|-\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|\right]{\Lambda }_{(k)(j)}{\nu }_{(j)}{V}_{(j)}$$

(24)

$$\begin{gathered} {\mathit{t}}_{(k)(j)}=\left[2ad\frac{\delta }{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}{\Lambda }_{(k)(j)}{\theta }_{(k)}+2b^{\prime }\delta s_{(k)(j)}\right]\frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|} \\ \mathrm{with} b^{\prime }={\mathit{G}}_{(b)1(k)(j)}{b}_1+{\mathit{G}}_{(b)2(k)(j)}{b}_2{\mathit{sin}}^2\phi {\mathit{cos}}^2\phi \frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|} \end{gathered}$$

(25)

where ${\mathit{G}}_{(d)(k)(j)}$ is the surface correction factor for volume dilatation; ${\mathit{G}}_{(b)1(k)(j)}$ and ${\mathit{G}}_{(b)2(k)(j)}$ are the surface correction factors for the PD material parameters $b_1$ and $b_2$, respectively.

3.4. Damage Model

In the PD damage model for isotropic materials, the damage of the bond between two material points can be determined by either the critical stretch ratio [27] or the critical energy release rate [32]. The critical stretch ratio $s_0$ is used in this work, and its expression in the OSB-PD model for two-dimensional problems is:

$$s_0=\sqrt{\frac{G_c}{\left(bh{\delta }^5+\frac{8}{9}a{d}^2{h}^2{\delta }^7\right)}}$$

(26)

where the parameters a, b, and d can be obtained by Equations (18) and (20) to import the influence of orientation angle; $G_c$ is the critical energy release rate, which can be estimated by the fracture toughness $K_{Ic}$ and the equivalent elastic modulus $E_{\phi }$.

$$G_c=\left\{\begin{array}{cc}\frac{K_{Ic}^2}{E_{\phi }}& \mathrm{plane} \mathrm{stress}\\ \frac{K_{Ic}^2}{E_{\phi }}\left(1-{\nu }^2\right)& \mathrm{plane} \mathrm{strain}\end{array}\right.$$

(27)

The break of the bond is described by a scalar function $\mu \left(\mathit{x},{\mathit{x}}^{\prime }, t\right)$ correlated with the time t. When the stretch ratio between a pair of material points satisfies $s\le s_0$, the bond breaks, that is, $\mu \left(\mathit{x},{\mathit{x}}^{\prime }, t\right)=0$; otherwise, $\mu \left(\mathit{x},{\mathit{x}}^{\prime }, t\right)=1$. In this way, a material point is determined to fail when the bonds between all the family points in its horizon domain and itself are broken. This accumulative damage process can be calculated by the failure function $\psi \left(\mathit{x},t\right)$ below.

$$\psi \left(\mathit{x},t\right)=1-\frac{{\int }_{H_{\mathit{x}}}\mu \left(\mathit{x},{\mathit{x}}^{\prime },t\right)d{V}^{\prime }}{{\int }_{H_{\mathit{x}}}d{V}^{\prime }}$$

(28)

In addition, the grain boundary strength has an important effect on the fracture behaviours of polycrystalline materials. De Meo et al. [23] defined a grain boundary strength coefficient β to study different fracture modes.

$$\beta =\frac{s_{0GB}}{s_{0GI}}$$

(29)

where $s_{0GI}$ and $s_{0GB}$ are the critical stretch ratios for the transgranular and intergranular fractures, respectively. When the intergranular fracture is dominated, then $\beta <1$, while when the transgranular fracture is dominated, then $\beta >1$.

In the following sections, the displacements in static and dynamic problems are solved by the adaptive dynamic relaxation (ADR) and the explicit time integration in Ref. [26], respectively.

4. Model Validation and Convergence Analysis

In this section, the finite element method (FEM) is used to verify the PD model. Through m-convergence and δ-convergence, the displacements under uniaxial tension and dynamic fracture under biaxial tension are analysed. δ is the horizon size, and m can present the density of material points in one horizon region with a relationship of m = δ/∆x (∆x is the spacing of material points). For all the cases in this section, the material is set as α-Fe, which is a body-centred cubic lattice, and its elastic stiffness constants are C₁₁ = 231.4 GPa, C₁₂ = 134.7 GPa, and C₄₄ = 116.4 GPa [2]. The fracture toughness is 58.4 $\mathrm{MPa}·\sqrt{m}$ [33].

4.1. Static Analysis for the Fe Single Crystal with Different Orientations

Figure 3 shows a thin rectangular plate under uniaxial tensile loading. It has a length L of 20 μm and a width W of 10 μm. The thickness H is 0.1 μm. On both sides of the plate, three groups of material points are added along the x-direction to impose boundary conditions, where the left boundary is completely fixed, and the pressure P applied to the right boundary is 150 MPa.

In order to analyse the effects of the m value and δ value on the static problem, the displacements UX and UY, respectively, along the centre lines at y = 0 and x = 0 are studied when the orientation angle γ = 0° and compared with FEM results, as shown in Figure 4a,b,d,e. The convergence criterion is $1\times 10^{-8}$. By comparing the slopes of the curves, it can be seen in Figure 4c,f that the increase in the m value makes the PD results closer to the FEM results when δ is constant. When m is constant, decreasing δ has little effect on the results.

Comprehensively considering the computational time and accuracy, $\delta =6\mathsf{\Delta }x$ and medium mesh size (200 × 100) are taken to compute the uniaxial tensile results at different orientation angles. The displacements at different orientation angles $\gamma \in \left[\mathrm{0,90}°\right]$ are simulated at an interval of 15°. As shown in Figure 5, the PD results are consistent with the corresponding FEM results at each orientation angle, indicating the effectiveness of the proposed PD model in simulating the static problem of a single cubic crystal with a given orientation angle.

4.2. Static Analysis for Polycrystalline Materials

The polycrystalline microstructure is generated by the Voronoi tessellation method. Voronoi seeds are calculated by the pseudo-random Sobel sequence to obtain the microstructure with relatively uniform grain size and distribution. Figure 6 shows a polycrystalline rectangular thin plate with 20 randomly oriented grains. The range of the orientation angles $\gamma$ is $\left[\mathrm{0,90}°\right]$. The computational domain has a length L of 5 mm and a width W of 2.5 mm. The plate thickness H is 0.01 mm. The numbers of material points in the x- and y-directions are 240 and 120, respectively. The horizon size $\delta$ is $6\mathsf{\Delta }x$. The loading condition is the same as the one for single crystals in Section 4.1. Figure 7 illustrates the displacement contours in the x- and y-directions from PD and FEM simulations, respectively. The convergence criterion is $1\times 10^{-8}$. The results of the two methods are almost consistent, which validates the effectiveness of the proposed PD method in calculating the static behaviours of polycrystalline materials.

4.3. Convergence Analysis of Dynamic Fracture of Polycrystalline Materials

Many studies have shown that mesh size and horizon size can significantly influence the crack propagation path in polycrystalline materials [23,27], but most size-independent convergence analyses are focused on static problems. It is important to obtain the convergence results of a crack propagation path for predicting and controlling the dynamic fracture of polycrystalline materials. Like static analysis, m-convergence and δ-convergence are also good methods to study the influence of meshing and horizon size on the crack propagation path. Three GBC values of 0.5, 2, and 1 are selected to represent intergranular, transgranular, and mixed fracture modes.

Figure 8 shows a 5 mm × 5 mm rectangular plate with preset cracks at the top and the bottom. It contains 100 grains with random orientation angles within [−45°, 45°]. The length of each crack Lc is 0.4 mm. The left and right sides of the plate are subjected to a rapid tensile velocity V of 15 m/s in opposite directions.

Figure 9, Figure 10 and Figure 11 illustrate the m-convergence results of different GBC values (fracture modes) at the time steps of 3000 (3 μs). It can be seen that when m is higher than 6, the changes in the crack propagation path in the three fracture modes decrease with the increase in the m value and tend to be converged.

In comparison, the δ-convergence results show more uncertainties. The change in the δ value affects the crack propagation path in the intergranular fracture more than in the transgranular fracture, as shown in Figure 12 and Figure 13. The decrease in δ indicates a reduction in the integral area when calculating the stretch ratio, resulting in the difference in damage determination. This difference will have a more prominent influence on the boundary of adjacent grains with large differences in the orientation angle, which may lead to different crack propagation paths.

In addition, when the horizon size is higher than the size of a single grain, the results will also be different. For a meshing size of 100 × 100, when m = 6, there are 112 material points in a complete horizon domain, while the average number of material points in one grain is about 100. For an intergranular fracture, the number of local microcracks will be significantly reduced (Figure 12e), while for a transgranular fracture, local damage may increase (Figure 14e).

Overall, the horizon size (δ) has a greater impact on crack propagation paths than the density of material points in one horizon region (m). When conducting dynamic fracture analysis of polycrystalline materials by the PD method, it is suggested to determine a proper δ value first and then increase the m value to achieve a faster convergence result of the crack propagation path. It is also necessary to ensure that the number of material points in a single grain is larger than that in a horizon region.

5. Applications on Dynamic Fractures

This section mainly studies and discusses the applications of the proposed model in two aspects. One is the effect of crystal orientations on the crack initiation and propagation of single-crystal materials with different notch shapes and sizes. The other is the effect of the degree of anisotropy in polycrystalline materials on propagation growth paths under different fracture modes.

5.1. Single-Crystal Materials with Different Notch Shapes and Sizes

Notches of various kinds exist in part/component design, such as chamfer, fillet, and thread. However, several studies have indicated that the stress concentration around the notch tip would significantly reduce the load capacity of parts/components [34]. For a single crystal, it has been found that the orientation difference between the notch front and crystal also makes the crack propagation direction change [3]. In order to study the effects of notch shape and size on the crack propagation path, three notch shapes including square, equilateral triangle, and semi-circle are designed on a rectangular plate, as shown in Figure 15. Dimension b is the side length, height, and radius of the square, triangle, and semi-circle, respectively.

The plate’s length, width, and thickness are 5 mm, 15 mm, and 0.1 mm, respectively. The single cubic crystal is silicon (C₁₁ = 165.6 GPa, C₁₂ = 63.9 GPa, C₄₄ = 79.5 GPa). The fracture toughness is 2 $\mathrm{MPa}\cdot \sqrt{m}$ for all orientation angles because it has little dependence on the tensile direction in the same crystal plane [3]. The velocity imposed on the top and the bottom side of the plate is ±0.5 m/s. The orientation angles are 0°, 30°, 45°, and 60°. The value of dimension b is 0.5, 1, and 1.5 mm.

Figure 16, Figure 17 and Figure 18 show the crack propagation paths of equilateral triangular, square, and semi-circular notches at different orientation angles when b is 1 mm. Under the same working conditions, the equilateral triangular notch is less likely to bifurcate than the others. At the orientation angles of 30° and 60°, the crack propagation paths of the three notch shapes are obviously shifted relative to the horizontal line and are basically symmetrical, which reflects the force shift caused by the orientation angle. For the square notch, when the orientation angles are 0° and 45°, the two right angles induce two crack routes (Figure 17a,c), while at the crystal orientation of 30° and 60°, the crack only occurs from one right angle (Figure 17b,d). As for the semi-circular notch, bifurcation is more likely to occur at a crystal orientation of 0° than at other angles, as shown in Figure 18.

Additionally, the notch dimension has a greater effect on the crack propagation path of the square notch, as shown in Figure 19. At the crystal orientation of 0°, when the notch dimension gradually decreases, the number of crack routes reduces from two to one. However, it has smaller changes in the pattern and number of crack routes in the triangular and semi-circular cases, so it is not shown for simplicity.

The crack initiation time and the average propagation speed are also studied. As shown in Figure 20a, when the crystal orientation is 0°, that is, the loading direction is parallel or vertical to the crystal orientation, the crack initiation is the slowest compared with the other angles for all notch shapes. As the notch size increases, crack initiation occurs faster. Among the three shapes, the semi-circular notch takes the longest time to initiate cracks.

However, the later occurrence of the crack does not mean that the crack propagation speed will be slow, as shown in Figure 20b. When γ = 0°, the propagation speeds of the semi-circular and square notches are the slowest, while the speed of the triangular notch is the fastest at the size of 1.5 mm and 0.5 mm. For different notch shapes, the semi-circular notch with a larger size has the lowest speed at γ = 0°.

In short, the numerical results showed that the crack propagation path of single crystals is sensitive to crystal orientation, notch shape, and size. To be more specific, semi-circular notches are the most resistant to fracture in most cases, which is why rounded corners are often used in structural design to reduce the stress concentration at edges and corners. Rectangular and triangular notches are more likely to crack because of stress concentration caused by sharp corners, which can be improved by modifying the tip to a rounded corner [35].

In addition, the change in crystal orientation results in a change in the crack bifurcation, the number of crack propagation paths, and the direction deviation. It needs to be analysed in combination with the material property and structure. For single-crystal silicon, experimental results showed that cracks occur along several specific cleavage planes. For example, the fracture plane of Si(100) sheets after biaxial tensile loading is mainly along the plane $\left\{110\right\}$ [3]. The results in Figure 18a,c may be closer to this experimental trend. However, further improvements in the theoretical part are required for the fracture with specific cleavage planes.

5.2. Effects of Anisotropy Degree in Polycrystalline Materials

The anisotropy degree caused by grain orientations has a significant influence on the fracture properties of polycrystalline materials. Based on the results in Section 4.3, the plate with a crystal orientation range of [−10°,10°] and a grain number of 200 is further simulated for comparison. The crack distribution results under different fracture modes are shown in Figure 21, Figure 22 and Figure 23.

Figure 21 shows that when intergranular fracture is dominant, the reduction in the anisotropy degree will make the microcrack around the main propagation routes disappear to a large extent. When the two fracture modes are not distinguished, the decrease in the anisotropy degree reduces the bifurcation in the crack propagation path, even approaching a straight line, as shown in Figure 22. When transgranular fracture dominates, the crack propagation with a low degree of anisotropy is significantly slower than that with a high degree of anisotropy (see Figure 23). These simulation results indicate that cracks in polycrystalline materials tend to occur at locations with larger differences in orientation angles, which has also been found in experimental studies [8].

Practically, the grain orientation and size of a polycrystalline material may be random or regular, depending on the manufacturing method. In the application of the model, its fracture toughness and crack propagation path can be calibrated through small amounts of experiments. Then, a more fracture-resistant grain structure and orientation can be obtained through optimisation design in simulations to guide the structural design and manufacturing of polycrystalline materials.

6. Conclusions

In view of the shortcomings of existing PD models in simulating crystal orientation, this work developed an improved OSB-PD model to simulate the granular fractures of cubic crystals involving arbitrary crystal orientations and delivered a complete derivation process and expressions of correction factors. Two numerical experiments are used to study and discuss the effect of crystal orientation on crack propagation. The main findings and innovations of this work are summarised as follows.

(1)

The periodic characteristics of the equivalent elastic modulus of cubic crystals are introduced into the PD parameters to achieve the simulation of arbitrary crystal orientations.

(2)

Convergence analysis is carried out in static and dynamic problems to obtain a proper density of material points in one horizon region (m value) and horizon size (δ value) to ensure computational effectiveness and accuracy. For a static problem, the m value has a greater effect on convergence than the δ value, while for a dynamic fracture problem, the δ value influences the crack propagation path more than the m value, especially in the intergranular fracture mode.

(3)

For convergence analysis on dynamic problems, a regulating strategy to obtain the converged and accurate results of crack propagation paths is proposed as follows: select an appropriate horizon size first and then increase the m value until the accuracy is satisfied.

(4)

In the numerical examples, the influence of crystal orientation on single-crystal materials with different notch shapes and sizes is mainly reflected in bifurcations, numbers, and propagation path directions of cracks. Under biaxial tensile loading, the single crystal with a semi-circular notch is more resistant to fracture than the crystal with square or triangular notches in most cases.

(5)

For polycrystalline materials, the decrease in the degree of grain anisotropy reduces microcracks in intergranular fracture and the crack propagation rate in transgranular fracture.

Therefore, the OSB-PD model proposed in this work can simulate the granular fracture behaviours of cubic crystals with arbitrary orientation angles and evaluate the fracture properties of single crystals and polycrystalline materials. The size parameters obtained by convergence analysis can ensure the accuracy of the crack propagation path and provide a reference for the microstructure design and manufacturing of polycrystalline materials. However, improvements are still needed for materials with specific fracture cleavage planes and three-dimensional problems.