Fatigue degradation strategies to simulate crack propagation using peridynamic based computational methods

Abstract

The aim of this paper is to develop new computational tools to study fatigue crack propagation in structural materials. In particular we compare the performance of different degradation strategies to study fatigue crack propagation phenomena adopting peridynamic based computational methods. Three fatigue degradation laws are proposed. Two of them are original. Initially a cylinder model is used to compare the computational performance of the three fatigue laws and to study their robustness with respect to variations of discretization parameters. Then the fatigue degradation strategies are implemented in a peridynamic framework for fatigue crack propagation analyses. Both cylinder model and peridynamic simulations show that the third proposed degradation law is unique in its combination of high accuracy, high stability and low computational cost.

Keywords
high cycle fatigue; peridynamics; fatigue degradation strategies

1 INTRODUCTION

Fatigue is a major cause of structural failure in many structures. Prediction of fracture propagation under cyclic loading is a challenging problem for two main reasons: from a theoretical point of view the discontinuous nature of fracture phenomena conflicts with the underlying mathematical structure of classical continuum mechanics, and from a more practical side the slow damaging effect of a large number of load cycles (high cycle fatigue) is not easy to model (Weinert and Gergely 2011). Traditional finite element method (FEM), or various modified versions of the FEM, have been used to simulate fatigue crack growth (Allix and Corigliano 1996; Branco et al. 2015; Point and Sacco 1996; Schellekens and Borst 1993). FEM models are based on the classical continuum mechanics theory, which cannot directly be used to simulate problems with discontinuities, such as cracks, requiring extra relationships that control the crack growth speed, direction and remeshing techniques. For example, the FEM with a remeshing technique is widely used to model fatigue crack growth (Branco et al. 2015). However, implementation of automatic remeshing techniques is difficult in 3D cases (Baydoun and Fries 2012; Boroomand et al. 2016; Boroomand et al. 2017), and is computationally costly to remesh geometries when multiple cracks exist. Interface elements equipped with cohesive zone models have been used (Alfano and Crisfield 2001; Alfano and Crisfield 2003; Blackman et al. 2003; Chen et al. 1999; Crisfield and Alfano 2002; Ostergaard et al. 2011; Qiu et al. 2001) in the cases in which the crack path is predictable. The extended finite element method (XFEM) (Belytschko and Black 1999; Bidokhti and Shahani 2015; Dolbow and Belytschko 1999; Shi et al. 2010) allows cracks to pass through the elements leading to better approximations of crack paths without remeshing. Nevertheless, one common problem for XFEM is that additional criteria are needed to study crack coalescence or branching.

Peridynamics, a newly presented continuum theory, has been developed to deal mainly with dynamic fracture problems by Silling in (Silling 2000). It is a non-local continuum theory proposed with the specific goal to tackle the shortcomings of the classical theory. Peridynamics consists of integral equations and does not need any spatial differentiation, for these reasons it is proper to describe problems affected by discontinuities (Bazazzadeh et al. 2018; Bobaru et al. 2016; Madenci et al. 2016; Shojaei et al. 2019). One of the significant features of the peridynamic approach is that it allows for spontaneous formation, interaction, and growth of discontinuities in a consistent framework (Xu et al. 2008). Using this theory, promising results have been attained in elastic bar (Aguiar et al. 2018) dynamic brittle fracture in glass (Ha and Bobaru 2010, 2011; Bobaru and Ha 2011; Hu et al. 2013), fiber-reinforced composites (Ghajari et al. 2014; Hu et al. 2012), functionally graded materials (Cheng et al. 2015) and corrosion damage (Chen and Bobaru 2015; Chen et al. 2016).

A peridynamic model for fatigue cracking has been proposed in the recent past by Silling and Askari (Silling and Askari 2014), which enables the simulation, in a single model, of the three phases of fatigue failure: crack initiation, fatigue growth and final failure controlled by quasi-static crack growth. In (Silling and Askari 2014), the dynamic relaxation method is used to obtain the static solution, and fatigue crack growth. References (Zhang and Bobaru 2016; Zhang et al. 2016) study fatigue crack propagation obtaining the static solution adopting a conjugate gradient energy minimization method, which is matrix-free and much faster than dynamic relaxation method, see (Le 2014). In a recent work (Hu et al. 2015; Hu and Madenci 2017) peridynamics is applied to simulate the fatigue life and residual strength in composite laminates. Moreover, by using peridynamic approach, mixed-mode fatigue crack propagation has been investigated in (Jung and Seok 2017). Other studies to simulate fatigue phenomena based on the peridynamic theory can be found in (Baber and Guven 2017; Jung and Seok 2016; Oterkus et al. 2010).

In the present paper a cylinder model (Galvanetto et al. 2009) is used to efficiently compare the computational performance of three fatigue laws. Three fatigue degradation strategies have been implemented in conjunction with a bilinear constitutive law. The computational performance of the three fatigue degradation strategies is examined by making use of various discretization parameters. Then the same constitutive laws and fatigue degradation strategies are implemented in a peridynamic based computational code. Due to the non-local nature of peridynamic methods, they are computationally more expensive than those based on classical continuum theories. In the conventional peridynamic approach, a constant horizon size and a uniform grid are assumed all over the domain. Accordingly, the efficiency of the peridynamic approach is strongly affected when a fine grid spacing is used. There have been many attempts, so far, to maintain the computational accuracy and reduce the cost by coupling peridynamics with local theories (Bobaru and Ha 2011; Seleson et al. 2013; Shojaei et al. 2016; Galvanetto et al. 2016; Shojaei et al. 2017; Zaccariotto et al. 2017; Zaccariotto et al. 2018), or by using parallel computing techniques (Lee et al. 2017; Mossaiby et al. 2017) and adaptive refinement algorithms (Dipasquale et al. 2014; Gu et al. 2017; Ren et al. 2016; Shojaei et al. 2018). Among the aforementioned techniques the FEM-PD coupling enhanced with adaptivity represents an efficient approach the use of which will be explored in future works.

The paper is organized as follows: in section 2 kinematics and statics of the cylinder model are defined. Then, the bilinear interface constitutive law is recalled in section 3. In section 4, three fatigue degradation strategies are described. In section 5 the computational performances of the three fatigue laws are compared. Finally the implementation of the three fatigue laws in a peridynamic code is illustrated in section 6 where the results they provide are also compared. Section 7 gives the conclusions.

2 THE CYLINDER MODEL

The mathematical model consists of a cylinder of radius R which rolls on a horizontal surface initially composed of two coinciding zero-thickness layers, see Figure 1(a) and ref. (Galvanetto et al. 2009). One of the layers is fixed and the other adheres to the surface of the cylinder at the contact point C. The upper layer remains adhered while the cylinder rotates and the contact point moves to the right (see Figure 1(b)). In the initial state the contact point is at $x=0$ (see Figure 1), the rotation, θ, of the cylinder is zero and the line OA, on the cylinder cross section, is vertical. The value of the rotation of the cylinder is defined as the angle between the line OA and the vertical line connecting O to the current contact point.

Figure 1
(a) Cylinder model in its initial configuration (b) The cylinder model after a rotation $\theta$ has occurred, adhesion between the layers is represented as a discrete system of springs.

If a clockwise moment is applied, the cylinder will rotate in the same direction and it will move to the right of its initial position and one layer will be lifted and separated from the other, as shown in Figure- 1(b). An adhesion stress acting between the two layers opposes the motion of the cylinder. The model represents one of the two symmetric halves of the tip of a crack opening in mode I as depicted in Figure- 2.

Figure 2
Pure moment loading for mode I specimen.

2.1 Continuum equations of the model

The model has only one degree of freedom, the rotation of the cylinder $\theta$ . The coordinate of the contact point $X_{con}$ represents the crack length, which is linked to the rotation $\theta$ as

$X_{con}=R\theta$ (1)

The applied moment,$M_a$, is in equilibrium with the reaction moment, $M_t$, due to the stresses in the interface

$M_t={\displaystyle \underset{0}{\overset{X_{con}}{\int }}\sigma (x)(X_{con}-x)dx}$ (2)

2.2 Discretization of the model

The motion of the cylinder detaches the two layers. The adhesion between them can be represented as a discrete system of vertical springs, with uniform spacing Δl, that oppose the separation, Figure 1.b. The springs are initially unstressed but become stressed as the two layers separate behind the advancing contact point.

A constitutive law for the springs has to be defined to describe the mechanical behavior of the interface which is progressively damaged: the force in the spring eventually reduces to zero (representing complete failure) as the layer separation increases.

The value of R is to be chosen as very large compared to the displacement at which the stretched springs fail so that the springs can be considered to remain vertical while they are stressed.

2.2.1 Kinematics

Once $X_{con}$ is known and given the spring spacing Δl, it is easy to identify the number of springs that become active (n) and their coordinates $x_i(i=\mathrm{1,}\mathrm{2,}\mathrm{3,}\dots n)$.

Figure 3 shows that simple geometry is sufficient to define the elongation δi of the i-th spring and its distance $d_i$ from the contact point.

Figure 3
Angle, opening displacement and distance from contact point of springs.

Since the springs have a zero initial length, the elongation will be defined as the vertical displacement (${\delta }_i$) of the top of the spring. Moreover it was previously mentioned that the spring eventually degrades as the displacement is increased and the stress finally reduces to zero. If ${\delta }_c$defines the critical displacement at which failure occurs, then the angle $\alpha$ at which such a displacement is achieved is given by (Figure 3):

$\alpha ={\cos }^{-1}\left(\frac{R-{\delta }_c}{R}\right)$ (3)

2.2.2 Equilibrium

As the cylinder advances, a number of springs become active. The interface stress (${\sigma }_i$), which is a function of the opening displacement (${\sigma }_i=f({\delta }_i)$) can be obtained by using an interface constitutive law.

The interface stress of each spring can then be used to find the spring force ($F_i$) generated by multiplying this value ${\sigma }_i$ by the area which the spring ‘represents’. For dimensional reasons the cylinder is assumed to have a width of 1. The interface area associated with a typical spring is therefore 1·Δl. An exception is for the area of the first spring which is taken as 1·Δl/2. The value of the force for a typical spring is given by

$F_i=1\cdot \Delta l{\sigma }_i$ (4)

The moment ($M_i$) generated by the spring force with respect to the contact point ($C$) is obtained by multiplying the force by the distance to the current contact point

$M_i=d_i{F}_i$ (5)

and the total moment ($M_t$) generated by the springs is found by adding the moments $M_i$ of all the active springs

$M_t={\displaystyle {\sum }_{i=1}^n{M}_i},$ (6)

This is the moment which has to be applied to the cylinder to maintain it in equilibrium for a given applied rotation. The critical moment $M_c$, which is the minimum value of the moment to break the interface is (Galvanetto et al. 2009):

$M_c=G_{c}R$ (7)

where G _c is the critical energy release rate. The moment M _c of the cylinder model can be found by imposing static equilibrium for a rotation of magnitude α.

3 INTERFACE CONSTITUTIVE LAWS

As it was mentioned previously, a constitutive law is required to describe the behavior of $\sigma$ as a function of $\delta$. The constitutive law should be defined in such a way that for every value of $\delta$ there is a unique value of $\sigma$. In this research the bilinear interface constitutive law has been considered (Alfano and Crisfield 2001; Alfano and Crisfield 2003; Crisfield and Alfano 2002; Galvanetto et al. 2009; Peerlings et al. 2000; Qiu et al. 2001; Robinson et al. 2005). There are two parameters, additional to the critical displacement value, ${\delta }_c$, that are needed to define an interface constitutive law: the critical energy release rate ($G_c$) and the maximum stress value (${\sigma }_0$). In the present case of the bilinear law the elastic limit relative displacement (${\delta }_0$) is also required (Galvanetto et al. 2009; Peerlings et al. 2000; Robinson et al. 2005).

Figure 4
Mode I bilinear interface constitutive law.

The critical energy release rate is a physical quantity used in Linear Elastic Fracture Mechanics (LEFM) for defining a criterion of propagation of the delamination in the interface. Its value is measurable in lab tests. It is indirectly used in the definition of the constitutive laws by equating the areas under the stress-displacement curves to the critical energy release rates. Therefore, the total energy dissipated during the delamination is correctly computed at each point although it is not released instantaneously as is assumed in LEFM (Alfano and Crisfield 2001).

δ0 is the maximum displacement for which the spring behavior is still linear elastic, whereas the ratio (${\sigma }_0/{\delta }_0$) is clearly linked to the stiffness of the interface.

The constitutive behavior for monotonic increasing relative displacement δ will be described below.

The bilinear interface law, as seen in Figure 4, is divided into three parts: linear elastic, stiffness degradation and failure. In the first part a linear-elastic behavior is assumed for a relative displacement that is increased from zero, where there is no stress ($\sigma =0$), until it reaches the elastic limit ${\delta }_0$. At the elastic limit the stress assumes its maximum value, ${\sigma }_0$. The linear-elastic part of the law is described by the following expression:

$\sigma =K\delta 0\le \delta \le {\delta }_0$ (8)

where $K$ is the initial stiffness of the interface and is given by

$K={\sigma }_0/{\delta }_0$ (9)

The second part of the interface constitutive law defines the behavior between the elastic limit and the critical displacement ${\delta }_c$. Once the elastic limit is exceeded the interface starts to be degraded. This degradation is described by the following relationship:

$\sigma ={\sigma }_0\left(\frac{{\delta }_c-\delta }{{\delta }_c-{\delta }_0}\right){\delta }_0<\delta <{\delta }_c$ (10)

The third and last part of the interface constitutive law is for relative displacement values that are equal or larger than the critical displacement value. When this critical value is reached or exceeded the element fails irreversibly so that no stress can be carried by the interface:

$\sigma =0\delta \ge {\delta }_c$ (11)

As it was stated previously the critical energy release rate $G_c$ is equal to the area under the $\delta -\sigma$curve, so in this case

$G_c=\frac{{\delta }_c{\sigma }_0}{2}$ (12)

Therefore of the four parameters $G_c$, ${\delta }_c$, ${\delta }_0$, ${\sigma }_0$ only three are independent.

The constitutive law of the interface can be written as well in terms of a damage variable D as:

$\sigma =K(1-D)\text{}\delta$ (13)

where the damage variable is:

$\{\begin{array}{cc}D=0& 0\le \delta \le {\delta }_0\\ D=\frac{{\delta }_c}{\delta }\left(\frac{\delta -{\delta }_0}{{\delta }_c-{\delta }_0}\right)& {\delta }_0<\delta <{\delta }_c\\ D=1& \delta \ge {\delta }_c\end{array}$ (14)

The increment of the opening displacement from ${\delta }_1$ to ${\delta }_2$ (${\delta }_2>{\delta }_1$) causes an increment of the damage variable:

$\Delta D=D_2-D_1=\frac{{\delta }_c{\delta }_0}{{\delta }_c-{\delta }_0}\left(\frac{1}{{\delta }_1}-\frac{1}{{\delta }_2}\right),\text{where}{\delta }_0<{\delta }_1<{\delta }_c,{\delta }_0<{\delta }_2<{\delta }_c,{\delta }_2>{\delta }_1$ (15)

4 FATIGUE DAMAGE

If the fatigue phenomenon is involved a time-like independent variable has to be considered, the number of cycles N of the external load. The applied load history consists of two components: quasi-static ramping phase and fatigue phase as shown in Figure 5.

Figure 5
(a) Envelope cyclic load approach and (b) envelope of the displacement for a large number of cycles.

In order to simplify the calculation process some assumptions are made (Peerlings et al. 2000). The first one is related to the cyclic load which is assumed to be sinusoidal and oscillating between zero and a maximum value. Since it would be very computationally expensive to simulate the relative displacement history for each spring, the load numerically applied to the structure, or to the cylinder model in this paper, is taken as constant and equal to the maximum value of the actual cyclic load as shown in Figure 5. Consequently, the relative displacement calculated at the crack tip will be the envelope of its true cyclic variation with time, as shown in Figure 5(b) (Galvanetto et al. 2009; Peerlings et al. 2000; Robinson et al. 2005).

We assume that the total damage rate in a spring can be obtained as the sum of the static and the fatigue damage rates as follows

$\dot{D}={\dot{D}}_s+{\dot{D}}_f$ (16)

${\dot{D}}_s$ is the rate of static damage due only to the increase in magnitude of $\delta$; ${\dot{D}}_f$ represent the fatigue damage increment due to the repetition of the load cycles and is greater than zero even if δ is constant. As the novelty of this work is in the fatigue part we use the same static damage rate based on the bilinear constitutive law for all different fatigue approaches.

Three strategies are considered to estimate the fatigue part of the damage.

4.1. Static damage rate

The static damage rate with respect to time can be obtained from Equation (14) as

${\dot{D}}_s=\frac{\partial D_s}{\partial t}=\frac{{\delta }_0{\delta }_c}{{\delta }_c-{\delta }_0}\frac{\dot{\delta }}{{\delta }^2}{\delta }_0<\delta <{\delta }_c$ (17)

One may integrate Equation (17) over an increment of the number of load cycles ($\Delta N$) as

$D_{sN+\Delta N}-D_{sN}=\frac{{\delta }_0{\delta }_c}{{\delta }_c-{\delta }_0}\left(\frac{1}{{\delta }_N}-\frac{1}{{\delta }_{N+\Delta N}}\right)$ (18)

in which, $D_{sN}$ and $D_{sN+\Delta N}$ denote the static damage component corresponding to the opening displacements ${\delta }_N$ and ${\delta }_{N+\Delta N}$ where ${\delta }_{N+\Delta N}\ge {\delta }_N$. Equation (18) provides the same expression as Equation (15) and therefore clarifies that the static damage component is due only to the increment of the opening displacement.

4.2 First fatigue degradation strategy

In the first strategy, originally presented in (Galvanetto et al. 2009; Robinson et al. 2005) the fatigue rate is assumed to have the following form:

${\dot{D}}_f=\frac{\partial D_f}{\partial N}=C{e}^{\lambda D}{\left(\frac{\delta }{{\delta }_c}\right)}^{\beta }$ (19)

$C$, $\lambda$ and $\beta$ are parameters to be determined by fitting of experimental data. The physical meaning of Equation (19) is that the fatigue damage rate is bigger for bigger relative displacement δ and for higher values of the current total damage.

Integrating Equation (19) over a number of loading cycles ($\Delta N$), one may write the increment of fatigue damage as

$D_{fN+\Delta N}-D_{fN}=\Delta N\frac{C}{1+\beta }{e}^{\lambda D_{\mu }}{\left(\frac{{\delta }_{\mu }}{{\delta }_c}\right)}^{1+\beta }$ (20)

Moreover, based on the reference (Galvanetto et al. 2009), $D_{\mu }$ and ${\delta }_{\mu }$ are taken as

$\begin{array}{l}{D}_{\mu }=(1-\mu )D_N+\mu D_{N+\Delta N}\\ {\delta }_{\mu }=(1-\mu ){\delta }_N+\mu {\delta }_{N+\Delta N}\end{array}$ (21)

Where $\mu$ value is 0.7. Additionally, by putting (18) and (20) into Equation (16), the total damage increment, based on the bilinear constitutive law and stiffness degradation damage (see Figure 6), can be obtained as

$D_{N+\Delta N}-D_N=\underset{static}{\underbrace{\frac{{\delta }_0{\delta }_c}{{\delta }_c-{\delta }_0}\left(\frac{1}{{\delta }_N}-\frac{1}{{\delta }_{N+\Delta N}}\right)}}+\underset{fatigue}{\underbrace{\Delta N\frac{C}{1+\beta }{e}^{\lambda D_{\mu }}{\left(\frac{{\delta }_{\mu }}{{\delta }_c}\right)}^{1+\beta }}}$ (22)

Figure 6
First (and third) fatigue law based on stiffness degradation strategy.

Note that in usual numerical procedures the only unknown variable in Equation (22) is $D_{N+\Delta N}$, it appears on both sides of Equation (22) which makes it implicit and nonlinear.

4.3 Second fatigue degradation strategy

The second fatigue degradation strategy, originally proposed by the authors in (Zaccariotto et al. 2015b), is based on the reduction of the value of the failure displacement, due to the fatigue phenomenon (Figure 7). It can be assumed that the failure elongation of a spring in each increment of load cycle, $\Delta N$, is a function of displacement as follows.

${\delta }_{cN+\Delta N}={\delta }_{cN}-B\text{}{({\delta }_{N+\Delta N}-{\delta }_0)}^n\Delta N$ (23)

Where $B$ and $n$ are parameters related to the material behavior when it is subjected to fatigue load cycles and ${\delta }_{cN+\Delta N}$ is the failure displacement of a spring after $N+\Delta N$ loading cycles.

The total damage at the end of each load cycle increment can be obtained as

$D_{N+\Delta N}=\frac{{\delta }_0{\delta }_{cN+\Delta N}}{{\delta }_{cN+\Delta N}-\delta {}_0}\left(\frac{1}{{\delta }_0}-\frac{1}{{\delta }_{N+\Delta N}}\right)$ (24)

Accordingly, the total damage increment in this approach is given by

$D_{N+\Delta N}-D_N=\frac{{\delta }_0{\delta }_{cN+\Delta N}}{{\delta }_{cN+\Delta N}-\delta {}_0}\left(\frac{1}{{\delta }_0}-\frac{1}{{\delta }_{N+\Delta N}}\right)-\frac{{\delta }_0{\delta }_{cN}}{{\delta }_{cN}-\delta {}_0}\left(\frac{1}{{\delta }_0}-\frac{1}{{\delta }_N}\right)$ (25)

In the above formula if $\Delta N\to 0$, only the static damage takes place due to the increment of external load so that Equation (25) is simplified to

$D_{N+\Delta N}-D_N=\frac{{\delta }_0{\delta }_{cN}}{{\delta }_{cN}-\delta {}_0}\left(\frac{1}{{\delta }_0}-\frac{1}{{\delta }_{N+\Delta N}}\right)-\frac{{\delta }_0{\delta }_{cN}}{{\delta }_{cN}-\delta {}_0}\left(\frac{1}{{\delta }_0}-\frac{1}{{\delta }_N}\right)=\text{}\underset{\text{Static}\text{Damage}}{\underbrace{\frac{{\delta }_0{\delta }_{cN}}{{\delta }_{cN}-\delta {}_0}\left(\frac{1}{{\delta }_N}-\frac{1}{{\delta }_{N+\Delta N}}\right)}}$ (26)

which is consistent with Equation (18). In Equation (24), the dependence on the number of cycles is implicit in ${\delta }_{cN+\Delta N}$ (see Equation (23)) and in the value of the current displacement. Since $D_{N+\Delta N}$ doesn’t appear on both sides of Equation (25), this numerical approach, differently from the first strategy, is explicit.

Figure 7
Second fatigue law, failure stretch shortening.

4.4. Third fatigue degradation strategy

In the third law, which has not presented before, it is assumed that the derivative of fatigue damage with regard to $N$ is given in Equation (27)

${\dot{D}}_f=\frac{\partial D_f}{\partial N}=A^{\prime }{\left(\frac{\delta }{{\delta }_c}\right)}^{m^{\prime }}$ (27)

In the above formula, $A^{\prime }$ and $m^{\prime }$ are constants related to the material behavior while $\delta$ represents the elongation of the spring. Equation (27) differs from equation (19) because in (27) the damage rate does not depend on the current value of total damage D.

Updated fatigue damage for crack growth can be determined as (see references (Bak et al. 2016; Peerlings et al. 2000))

$D_{fN+\Delta N}=D_{fN}+{\displaystyle \underset{N}{\overset{N+\Delta N}{\int }}\frac{\partial D_f}{\partial N}dN}$ (28)

Using a mid-point rule, one may approximately calculate the integral of the fatigue damage rate as

${\displaystyle \underset{N}{\overset{N+\Delta N}{\int }}\text{}\frac{\partial D_f}{\partial N}dN\approx }{A}^{\prime }\Delta N{\left(\frac{1}{2}(\frac{{\delta }_N}{{\delta }_c}+\frac{{\delta }_{N+\Delta N}}{{\delta }_c})\right)}^{m^{\prime }}$ (29)

Accordingly, by combining Equations (29) and (18) in Equation (16), the increment of total damage due to an increment $\Delta N$ of the number of cycles is given by the following expression, explicit in the unknown $D_{N+\Delta N}$:

$D_{N+\Delta N}-D_N=\underset{static}{\underbrace{\frac{{\delta }_0{\delta }_c}{\delta {}_c-{\delta }_0}\left(\frac{1}{{\delta }_N}-\frac{1}{{\delta }_{N+\Delta N}}\right)}}+\underset{fatigue}{\underbrace{A^{\prime }\Delta N{\left(\frac{1}{2}(\frac{{\delta }_N}{{\delta }_c}+\frac{{\delta }_{N+\Delta N}}{{\delta }_c})\right)}^{m^{\prime }}}}$ (30)

Note that in the third law, as in the first one, damage directly affects the stiffness of the springs as shown in Figure 6.

Once the sum of static and fatigue damage increments has been computed, Equation (13) is used to compute the stress in the spring.

5 COMPARISON OF THE PERFORMANCE OF THE THREE FATIGUE DEGRADATION STRATEGIES APPLIED TO THE CYLINDER MODEL

In this section we compare the computational performances of the three fatigue degradation strategies proposed in section 4. The fatigue phenomena considered are described by the following Paris law, based on the experimental data gathered from (Asp et al. 2001), (the specimen is under pure mode I loading):

$\frac{\partial X_{con}}{\partial N}=C_I{\left(\frac{G}{G_c}\right)}^{m_I}{G}_0<G<G_c$ (31)

where $C_I$ and $m_I$ are material parameters, $X_{con}$ and $N$ represent the crack length and number of load cycle respectively. $G$,$G_c$ and $G_0$ express the maximum, critical and threshold energy release rates respectively. For the cylinder model, using Equation (7), one may rewrite Equation (31) as

$\frac{d{X}_{con}}{dN}=C_I{\left(\frac{M_a}{M_c}\right)}^{m_I}$ (32)

$C_I$ and $m_I$ assume the values given in Table 1, provided in reference (Bak et al. 2017). In a log-log graph representing $d{X}_{con}/dN$vs $G/G_c$Paris law is a straight line, as shown by the red line in Figure 8.

Figure 8
Paris law and results produced by the three fatigue degradation strategies.

Table 1
Parameter values of Paris law and cylinder model.

	Parameter		Value
	$C_I$		0.00308 mm/cycle
	$m_I$		5.4
	$G_c$		0.260 Nmm/mm2
	$R$		100 mm
	$\Delta l$		0.005 mm
	$\Delta N$		100
	${\sigma }_0$		30 N/mm2
	${\delta }_0$		10-6 mm
	${\delta }_c$		0.01733 mm

The parameters of the three fatigue degradation strategies are chosen so that they reproduce the Paris law and the relevant results are shown in Figure 8. The parameter values are shown in Table 2.

Figure 9(a) and 9(b) show the excellent agreement among the results obtained with the three fatigue laws for two cases $M_a=0.2M_c$and $M_a=0.3M_c$. Given the initial static ramp of the load phase, the crack fatigue propagation starts from a value different from zero.

Table 2
Parameters used in the three proposed fatigue degradation strategies.

	Fatigue laws		Model specific parameters
	First		$C=2\times {10}^6,\lambda =\mathrm{0.5,}\beta =2.0$
	Second		$B=\mathrm{2540,}n=3.8$
	Third		$A^{\prime }=\text{1}\text{.099}\times {\text{10}}^{-6},m^{\prime }=3.0$

Figure 9
Crack length $X_{con}$ for (a) $M_a=0.2M_c$ and (b) $M_a=0.3M_c$ using the three different strategies.

As depicted in Figure 9, the fatigue crack propagation is divided into two phases. At the beginning of crack propagation, one may observe a slow increase of the crack length, and then the crack grows more quickly with a constant slope. It is also apparent that the slope $d{X}_{con}/dN$is bigger for the case $M_a=0.3M_c$. All three proposed approaches are able to produce the same slope, $d{X}_{con}/dN$, with an excellent agreement.

All the codes are developed using standard Matlab 2016 language.The run times reported later are measured on an Intel Core i7-6700HQ 2.60 GHz CPU, on a 64 bit Windows 10 Enterprise system.

It takes 59.84 CPUs and 57.32 CPUs to solve this example using the second and third fatigue laws respectively when$M_a=0.2M_c$, while this computational time in the first law is 68.01 CPUs which is approximately 15-20% slower than the other two.

Normalized stress vs. displacement graphs are shown in Figure 10 for an energy release rate of$G=0.5G_c$, i.e. the applied moment is $M_a=0.5M_c$. The three graphs show the stress-elongation curve (σ−δ curve) for the same spring in the constant slope propagation phase. The three approaches generate a similar curve due to the simultaneous action of increasing displacement and fatigue cycles. For all cases the area under the σ−δ fatigue curve approximates well the green area which represents half of the area under the bi-linear law.

Figure 10
Single spring behavior under static and fatigue load using (a) the First, (b) the second and (c) the third approach. The green area represents half of the triangular area under the bi-linear law, which is a measure of the energy dissipated during static crack propagation.

In many practical problems, the applied $G$ is not constant. So, following what was done also in (Bak et al. 2016), we compare the performance of the three fatigue degradation strategies when used with a variable G (i.e a $M_a$ in the case of the cylinder model). We alter the applied moment instantly as shown in Figure 11. The response to this energy release change can be observed in the crack length vs. the number of load cycles as depicted Figure 11. It can be noticed that there is a good agreement between all proposed approaches and the Paris law. The results of the Paris law are obtained by integrating numerically Equation (31) to evaluate the crack length as a function of energy release rate as (see Reference (Bak et al. 2016))

$X_{con}={{\displaystyle \underset{}{\overset{}{\int }}C\left(\frac{G(X_{con})}{G_c}\right)}}^{m}dN$ (33)

Figure 11
Blockwise constant energy release rate $G$.

Moreover, in many engineering problems gradual energy release changes occur. Assuming the energy release rate as a function of crack length, one may change the loading condition in a way to obtain an increasing/decreasing energy release rate as depicted in Figure 12.

Figure 12
Simulation linear dependency of the energy release rate ratio $G(X_{con})/G_c$ on the crack length.

For increasing $M_a$ the first fatigue degradation strategy provides results slightly better than the third law, whereas the results of the second law are worse. For decreasing M _a the three laws have a similar performance (see Figure 13).

Figure 13
Crack growth for an linearly variable energy release rate as a function of the crack length.

5.1 Stability of the results with respect to variations of discretization parameters

The results presented in the previous section seem to suggest that the performance of the three degradation strategies are rather similar; the second law has a worse performance in the case of increasing M _a; and the first law requires a longer computational time than the other two, but the overall behavior of the three laws does not exhibit apparent differences. In this section we want to compare the stability of the results provided by the three laws with respect to variations of two important discretization parameters:$\Delta l$, the spring spacing, and $\Delta N$, the increment of load cycles in one step. The best approach is the one for which the slope $d{X}_{con}/dN$ is affected the least by changes in $\Delta l$ and $\Delta N$. All parameters are given in Table 2. The applied moment is set to be $M_a=0.2M_c$ and the “exact slope” is generated by adopting small values of $\Delta N$(=100) and $\Delta l$ (=0.005mm). The ‘exact slope’ is the same for all fatigue degradation strategies. Note that smaller values of $\Delta l$ and $\Delta N$ increase the cost of the procedure but reduce the error in the slope $d{X}_{con}/dN$. We consider values of $\Delta l$ in the range$0.005\text{mm}<\Delta l<0.21\text{mm}$, which is divided into 206 uniform subintervals of amplitude equal to 0.001. $\Delta N$ is considered in the range $100<\Delta N<60000$, which is divided into 600 uniform subintervals. Consequently, 123600 different parameter pairs are taken into account to produce three contour plots of the error, shown in Figure 14. For every pair of values (Δl, ΔN) the relevant slope is computed and compared to the value of the ‘exact slope’. Different colors in Figure 14 represent different percent errors in the computed slope with respect to ‘exact slope’. Analyzing these error-plots, we observe that as $\Delta l$ and $\Delta N$ increase the error in general rises. The region of small error is much wider for the third law and much smaller for the second. Therefore, the third law, not only is computationally cheaper than the other methods for the same values of Δl and ΔN, but also can be used with larger values of $\Delta l$ and $\Delta N$.

Figure 14
Error plots using bilinear constitutive law for (a) first (b) second and (c) third fatigue laws.

Therefore, the cylinder model suggests that the third fatigue law is preferable. In the next section, we verify such a suggestion by implementing the three laws in a peridynamic-based code. In particular, the bilinear constitutive law equipped with the three fatigue degradation strategies will be adopted as constitutive law that defines the force in the bond as a function of the bond stretch.

6 PERIDYNAMICS FATIGUE MODEL

Bond based peridynamics (BBPD) is a non-local theory of continuum proposed in Reference (Silling 2000). Since the formulation of the theory avoids differentiations and is only based on spatial integration BBPD is particularly suited for describing crack propagation problems. For a point, $x$, and time, $t$, the peridynamic static equilibrium equation is written as (see Figure 15)

$0={\displaystyle \underset{H_x}{\int }f(u(\widehat{x},t)}-u(x,t),\widehat{x}-x){\text{dV}}_{\widehat{x}}+b(x,t)forx\in \Omega andt\in [t_0,t_f]$ (34)

Figure 15
Neighborhood of a generic material point $x$.

In which, $\Omega$ is the body domain, $f$ is the pairwise force function in each peridynamic bond between node $x$ and $\widehat{x}$, $u$, is the displacement vector field, $H_x$ is a local spherical region of radius ${\delta }_H$, in which the integral is defined. $\xi$ and $\eta$ are the initial relative position and current relative displacement respectively:

$\begin{array}{l}\xi =\widehat{x}-x\\ \eta =u(\widehat{x},t)-u(x,t)\end{array}$ (35)

The stretch of the bond can be expressed as

$s=\frac{\left|\xi +\eta \right|-\left|\xi \right|}{\left|\xi \right|}$ (36)

6.1. Spatial Discretization of the model

We discretize the solution domain with a grid of points called nodes (see Figure 16), to approximate the peridynamic equation.

Figure 16
A spatial discretization for a peridynamic model.

By adapting the one-point Gauss quadrature rule to proceed with the spatial integration, one can discretize Equation (34) as

$0={\displaystyle \sum _{j=1}^{N}f(u(x_j,t)-u(x_i,t)},x_j-x_i){\beta }_j(\xi )V_j+b(x,t)$ (37)

where, $x_i$ is the point of interest and $x_j$ denotes the family nodes located within the horizon of the node $x_i$. $V_j$ represents the volume of the collocation point $x_j$ . Moreover,$\beta$ is the volume correction factor which determines the portion of $V_j$ within the neighborhood of source node $x_i$ . In this study, $\beta$ value is determined as proposed in (Yu et al. 2011):

$\beta =\{\begin{array}{l}1\text{for}\Vert \xi \Vert \le {\delta }_H-0.5\Delta x\\ \frac{{\delta }_H+0.5\Delta x-\Vert \xi \Vert }{\Delta x}\text{for}{\delta }_H-0.5\Delta x<\Vert \xi \Vert \le \\ 0\text{otherwise}\end{array}{\delta }_H+0.5\Delta x$ (38)

Equation (38) expresses that the volume of nodes close to the boundary of the neighborhood falls only partially within the horizon of the source node. Additionally, the distance between two adjacet nodes is referred to as the grid spacing.

As proposed in (Silling 2000), in the prototype microelastic brittle material (PMB), the pairwise force function $f$ and the stretch $s$ are related according to the following equation

$f=cs\mu (\xi )\frac{\eta +\xi }{\left|\eta +\xi \right|}\cong cs\mu (\xi )\frac{\xi }{\left|\xi \right|}$ (39)

in which $c$ is the micro-modulus that can be interpreted as the bond stiffness, $\mu (\xi )$ is a history-dependent function which takes the value of 1 for active bonds and 0 for broken bonds. Furthermore, for “small deformations” the magnitude of the bond force, $f$, is approximated as the second part of the Equation (39): $f=cs\mu$. Considering the elastic deformation energy, one can define the micro-modulus $E$ in terms of the Young’s modulus $E$ and of the horizon radius ${\delta }_H$ as represented in Equation (40) and Equation (41) for 3D and and planar cases respectively.

$c=\frac{12E}{\pi {\delta }_H{}^4}\text{3D}\text{case}$ (40)

$\{\begin{array}{c}c=\frac{9E}{\pi t_h{\delta }_H{}^3}\text{Plane}\text{stress}\\ c=\frac{48E}{5\pi t_h{\delta }_H{}^3}\text{Plane}\text{strain}\end{array}$ (41)

in which, $t_h$ represents the thickness. The value for Poisson’s ratio for the 3D and plane strain cases is fixed to $\nu =0.25$; however, for the plane stress cases this value is restricted to $\nu =1/3$.

The limitation of Poisson’s ratio is resolved in (Silling et al. 2007) and (Warren et al. 2009) by using the state-based formulation of the peridynamic theory.

6.2 Pairwise force function

Based on Reference (Silling and Askari 2005), for a PMB material the stretch in a bond is related to the pairwise function $f$ by Equation (39). When the stretch reaches its limit value,$s_0$, the bond breaks so that the global behavior of the system is nonlinear. In this section the constitutive law of the bond $f=f(s)$ is assumed to be bilinear of the type presented in Section 3. The stretch of the bond plays the same role as the elongation of the spring. $s_0$ is the elastic limit, and $s_c$ is the failure limit of the stretch, as shown in Figure 17. According to reference (Ha and Bobaru 2010), following expressions are assumed:

$s_0=\sqrt{\frac{5G_0}{6E{\delta }_H}}\text{3D}\text{case}$ (42)

$\{\begin{array}{c}{s}_0=\sqrt{\frac{4\pi G_0}{9E{\delta }_H}}\text{Plane}\text{stress}\\ s_0=\sqrt{\frac{5\pi G_0}{12E{\delta }_H}}\text{Plane}\text{strain}\end{array}$ (43)

where, $G_0$, is the energy required to start the damaging process. The ratio between $G_c$ and $G_0$, called $k_r$, is defined as (see Reference (Zaccariotto et al. 2015a)):

$k_r=\frac{G_c}{G_0}=\frac{s_c}{s_0}$ (44)

Equation (44) provides the value of $s_c$ for the bilinear constitutive law.

Figure 17
Magnitude of the bond force f vs bond stretch $s$ law: (a) elastic brittle material (b) elastic-progressively damaging material.

Using a bilinear constitutive law, one may define the pairwise force function magnitude, $f$, as

$\{\begin{array}{cc}f(s)=cs& 0\le s\le s_0\\ f(s)=c(1-D)\text{}s& s_0<s<s_c\\ f(s)=0& s\ge s_c\end{array}$ (45)

in which, $D$, represents the total damage (sum of the static damage and the fatigue damage) in the bond.

The damage rates are defined according to the equations shown in section 4.The increment of static damage in a bond can be obtained with respect to Equation (18) as follows

$D_{sN+\Delta N}-D_{sN}=\frac{s_0{s}_c}{s_c-s_0}\left(\frac{1}{s_N}-\frac{1}{s_{N+\Delta N}}\right)$ (46)

Moreover, similar to the cylinder model, the increment of total damage in each bond for the first fatigue strategy is given by (see Equation (22))

$D_{N+\Delta N}-D_N=\underset{static}{\underbrace{\frac{s_0{s}_c}{s_c-s_0}\left(\frac{1}{s_N}-\frac{1}{s_{N+\Delta N}}\right)}}+\underset{fatigue}{\underbrace{\Delta N\frac{C}{1+\beta }{e}^{\lambda D_{\mu }}{\left(\frac{s_{\mu }}{s_c}\right)}^{1+\beta }}}$ (47)

where,

$s_{\mu }=(1-\mu )s_N+\mu s_{N+\Delta N}$ (48)

Similarly, for the second fatigue law, the total damage of each bond at the end of a load cycle can be obtained as (see Equation (24))

$D_{N+\Delta N}=\frac{s_0{s}_{cN+\Delta N}}{s_{cN+\Delta N}-s{}_0}\left(\frac{1}{s_0}-\frac{1}{s_{N+\Delta N}}\right)$ (49)

Finally, for the third fatigue law, the increment of total damage in each bond, can be written as (see Equation (30))

$D_{N+\Delta N}-D_N=\underset{static}{\underbrace{\frac{s_0{s}_c}{s{}_c-s_0}\left(\frac{1}{s_N}-\frac{1}{s_{N+\Delta N}}\right)}}+\underset{fatigue}{\underbrace{A^{\prime }\Delta N{\left(\frac{1}{2}(\frac{s_N}{s_c}+\frac{s_{N+\Delta N}}{s_c})\right)}^{m^{\prime }}}}$ (50)

The damage level in a point $x$, at cycle number $N$, can be expressed by

$\varphi (x,N)=1-\frac{{\displaystyle \underset{H(x)}{\int }\mu (\xi ,N)d{V}_{x^{\prime }}}}{{\displaystyle \underset{H(x)}{\int }d{V}_{x^{\prime }}}}$ (51)

In Equation (51), $\varphi$, indicates the ratio of the number of broken bonds to the total number of bonds originally connected to the point $x$. The damage level, $\varphi$, takes a value between 0 and 1, $0\le \varphi \le 1$ . $\varphi =1$ means there is no interaction between the point and all surrounding points within its horizon, while $\varphi =0$ denotes an undamaged state.

6.3. Nonlinearity of the system

The problem becomes nonlinear after the bonds get damaged, therefore, the stiffness matrix of the structure is not constant. Various techniques have been used to satisfy the equilibrium in the whole structure. In this paper, the Newton-Raphson (N-R) method is adapted to solve Equation (37) and update the stiffness matrix in each iteration.

Using the conventional finite element concept, one is able to write the relationship between external and internal forces as follows:

$Ku=F_{ext}$ (52)

in which, $K$ represents the global stiffness matrix of the structure, $u$and $F_{ext}$ stand for the nodal displacements and external force vector respectively.

An out of balance force in the whole system is generated due to the nonlinearity of the problem. Therefore, out of balance force vector is represented by

$g=F_{int}-F_{ext}=Ku-F_{ext}$ (53)

The equilibrium is satisfied when the norm of $g$ is lower than a given tolerance.

According to the Newton-Raphson solution, the derivatives of$g$ with respect to the nodal displacement vector,$u$, have to be calculated as follows

$\frac{\partial g}{\partial u}=\frac{\partial }{\partial u}\left(Ku-F_{ext}\right)=K+\frac{\partial K}{\partial u}u$ (54)

At a given point in which $u=u_0$, the term $\frac{\partial g}{\partial u}{|}_{u=u_0}$ stands for the tangent stiffness matrix, $K_T$. One may calculate each component of this matrix by

$K_{T_{ij}}=K_{ij}+{\displaystyle \sum _{m=1}^N\frac{\partial K_{im}}{\partial u_j}{u}_m}$ (55)

in which, N, represents the total number of degrees of freedom. i and j denote the rows and the columns of the $K$ and $K_T$ respectively.

In the case that the stretch of a bond exceeds the failure stretch, $s_c$, the bond is removed from the whole structure.

6.4. Loading and boundary conditions

The application of boundary conditions in peridynamics is a significant issue: since integral operators, rather than partial differential operators, are used to calculate the peridynamic equilibrium equations, the traction boundary conditions should be applied over a volume (Madenci and Oterkus 2014).

In order to converge to the equilibrium solution, static forces are applied gradually. After static equilibrium is satisfied, the fatigue phase initiates in the model, as shown in figure 5(a).

Numerical models with static damage have been comprehensively investigated in (Zaccariotto et al. 2015a), where deformations of the structures have been compared to the exact and FEM numerical solutions. Since the main contribution of this work is on the fatigue damage, the authors just present the examples with static and fatigue damage to further validate the fatigue degradation strategies in a peridynamic framework.

6.5. Double cantilever beam peridynamic example

In order to verify the indications of the cylinder model, a double cantilever beam is simulated with a BBPD code as depicted in Figure 18. The beam has the length of $L=0.0161$ m and width of $W=0.0029$ m. Its Young’s modulus and pre-crack length are specified as $E=70$ GPa,$a=0.005$ m. The applied moment is $M_a=0.048\text{Nm}$ which is sufficient to introduce static damage in several bonds during the initial static ramp. Furthermore, based on Equation (44): $k_r=20$.

Figure 18
Double cantilever beam model discretized with a BBPD software; (a) Initial structure with a pre crack of $a=0.005$ m (b) Deformed structure after $4\times {10}^6$ load cycles (a scale factor of 20 has been applied to displacement values).

Assuming a constant horizon size, ${\delta }_H=3.75\times {10}^{-4}$ m, and the grid spacing,$\Delta l$, namely $m={\delta }_H/\Delta l$, we are able to examine the efficiency of each fatigue low by using $m$-convergence approach, three different values of m and $\Delta l$are used as: $m=\mathrm{2,3,4}$and $\Delta l=\text{1}\text{.8971}\times {\text{10}}^{-4},\text{1}\text{.25}\times {\text{10}}^{-4},\text{9}\text{.4298}\times {\text{10}}^{-5}$ m respectively. In Figure 19, the slope of the lines ‘crack-length vs number of cycles’ in the third law converge for m values equal to 3 and 4. Due to the fact that $m=3$ is computationally less expensive than $m=4$, $m=3$ is employed for the rest of analysis.

Figure 19
Slope of the lines ‘crack-length vs number of cycles’ by using fixed values of ${\delta }_H$ and $\Delta N$(${\delta }_H=3.75\times {10}^{-4}$ and $\Delta N=1000$) and different m values for the (a) first (b) second and (c) third fatigue laws.

This example has been solved using different grid spacing as $\Delta l=6.27\times {10}^{-5}\text{m}$, $\Delta l=1.06\times {10}^{-4}\text{m}$ and $\Delta l=1.25\times {10}^{-4}\text{m}$; and different increments of load cycles, ΔN=1000, 10000 and 50000 for the first and the third proposed fatigue laws, in contrast, due to convergence problems we used smaller values of 500, 1000 and 10000 for the second fatigue law. (See Table 3).

Table 3
Ranges of $\Delta N$ and $\Delta l$ investigated using the first (purple triangle), the second (red square) and the third (green circle) fatigue laws

$\Delta l(\text{m})$	$\Delta N$
	500	1000	10000	50000
$6.27\times {10}^{-5}$	------	------	■	▲●
$1.06\times {10}^{-4}$	■	▲■●	▲■●	▲●
$1.25\times {10}^{-4}$	■	▲■●	▲■●	▲●

Assuming a constant parameter, $m=3$, one can examine the efficiency of the fatigue laws via ${\delta }_H$ -convergence approach (see Figure 20). The horizon size, ${\delta }_H$, is reduced and the grid spacing is equal to $\Delta l={\delta }_H/m$ .

In this example, the crack length with respect to the number of load cycles, $N$, is computed by using the three fatigue strategies (See Figure 21).

Figure 20
${\delta }_H$ - Convergence with a fixed m ratio.

Figure 21 :
Number of cycles vs. crack length for different values of ${\delta }_H$ and $\Delta N$for (a) the first (b) the second and (c) the third laws.

Figure 21 shows the second law to be much more sensitive to parameter variations than the other two: it does not work for ΔN= 50000 and for smaller values of ΔN provides results that are clearly more scattered than those produced by third law. Moreover the comparison between Figure 21(a) and 21(b) shows that if the lines in Figure 21(a) obtained for ΔN= 50000 are neglected, since there is no equivalent counterpart in Figure 21(b), the results of the second law are worse even than those of the first law. Finally comparing Figure 21(a) and 21(c) confirms the findings of the cylinder model: the third law is much more stable than the other two and provides very similar values of the slope of the line ‘crack-length vs number of cycles’ in a reasonably wide range of values of ${\delta }_H$ and ΔN.

Additionally, Figure 21 shows that the third fatigue law is less sensitive than the other two to the horizon size, ${\delta }_H$.

Both the cylinder model and the DCB discretized with the peridynamic based code confirm that the third law is the best of the three considered in the present work because it is more stable with respect to the variations of the discretization parameters and it is cheaper from a computational point of view.

6.6. Experimental validation of the model

To further validate the presented numerical methods, we simulate a DCB, which is taken from an experimental study (Juntti et al. 1999), using BBPD code. Due to the good convergence and high computational speed of the third fatigue law which represented in previous sections, hereinafter we only use the third law for the rest of the examples. In this example, we investigate the effectiveness of the selected fatigue degradation strategy by comparing with the experimental data (Juntti et al. 1999). The specimen used for this simulation is shown in Figure 22. The beam has the length of $L=0.15$ m, thickness of $W=0.0031$ m and width of $B=0.02$ m. Its Young’s modulus, Poisson’s ratio and pre-crack length are specified as $E=120$ GPa, $\nu =0.3$ and $a=0.035$ m. Although the specimen is 3D we simplify it into a 2D plain strain peridynamic model. The specimen arms are loaded under pure mode I with opposing moments as depicted in Figure 2. (This type of loading with constant applied moment, $M_a$, should cause a constant crack growth length so that the energy release rate is independent of crack length). Bilinear constitutive law parameters are set as:$G_c=260$ J/m² and $k_r=4.38$. Moreover, fatigue constant parameters of the fatigue approach are shown in Table 4.

Figure 22
Geometry of the specimen.

Table 4
Parameters used in the proposed fatigue degradation strategies for the DCB model.

	Fatigue law		Model specific parameters
	Third		$A^{\prime }={\text{10}}^7,m^{\prime }=4.0$

The energy release rate under pure mode I for a DCB, is related to the applied moment, $M_a$, as (Mi and Crisfield 1996)

$G=\frac{M_a^2}{BEI}$ (56)

Where, $I$, is the second moment of area of the specimen arm.

In a log-log graph representing $d{X}_{con}/dN$vs $G/G_c$, the third fatigue law using $m=3$ and $m=4$, is able to produce the linear part of Paris plot (See Figure 23).

Figure 23
Paris plot for mode I with experimental and BBPD numerical results.

6.7. Center cracked beam with a hole

In this section, fatigue crack propagation is simulated in a single notched beam with a hole by using the third fatigue law (see Figure 24). Experimental data on this test can be found in (Miranda et al. 2003). As shown in Figure 24, the beam is loaded in two points at the top and supported in two other points at the bottom. This problem has been solved under plain strain conditions with the physical parameters of: $E=205$GPa, $\nu =0.3$, $a=0.025$m and $P=100$ N. Bilinear constitutive law parameters are set as:$G_c=1.07\times {10}^5$ J/m², ${\delta }_H=1.143\times {10}^{-3}$, $m=4$ and$k_r=4.00$. Moreover, fatigue constant parameters are taken as: ${{\rm A}}^{\prime }=2.1\times {10}^8$ and $m^{\prime }=4.0$.

Figure 24
Beam with a center crack, loaded in two points and supported in two points.

The solution domain is discretized with a 106×438 uniform grid of nodes which results in 45,390 nodes. This example has been solved up to $6.7\times {10}^5$ load cycles using $\Delta N=50000$ cycles. As shown in Figure 25, the crack path grows toward the center hole. Moreover, the crack path is in a very good agreement with experimental data in (Miranda et al. 2003).

Figure 25
Comparison of the crack paths between numerical results and experimental results.

7. CONCLUSIONS

The paper compares three fatigue degradation strategies to be used in the simulation of fatigue crack propagation. The adopted constitutive law of the material is bilinear. The increment of damage due to fatigue is added to that due to the static increment of the opening displacement and is defined in three different ways. A mathematical cylinder model is employed to carry out a reasonable comparison of the computational efficiency of the three fatigue degradation strategies.

Fatigue degradation strategy 1 assumes that the fatigue damage rate is a function of the opening displacement δ and of the current value of the total damage D. In fatigue degradation strategy 2 the damage directly affects the critical value of the opening displacement δ _c. Finally, fatigue degradation strategy 3 assumes that the fatigue damage rate is a function only of the opening displacement δ.

Then the three laws are implemented in a code using bond based peridynamics to simulate fatigue crack propagation, and the sensitivity of the results to variations of the main discretization parameters is assessed. Both the cylinder model and the BBPD code provide the same assessment of the three fatigue degradation strategies. Two main conclusions can be obtained by the current work:

1. The third fatigue degradation strategy is the best, among those investigated, for being used in BBPD codes.

2. The cylinder model appears to provide reliable indications about the computational performance of the fatigue laws.

8. ACKNOWLEDGMENT

The financial support of the University of Padova BIRD2017 NR.175705/17 is gratefully acknowledged.

REFERENCES

• Aguiar, A.R., Patriota, T.V.B., Royer-Carfagni, G., Seitenfuss, A.B., (2018). Boundary layer effects in a finite linearly elastic peridynamic bar, Latin American Journal of Solids and Structures: In press.
• Alfano, G., Crisfield, M.A., (2001). Finite element interface models for the delamination analysis of laminated composites: Mechanical and computational issues, International Journal for Numerical Methods in Engineering 50: 1701-1736.
• Alfano, G., Crisfield, M.A., (2003). Solution strategies for the delamination analysis based on a combination of local-control arc-length and line searches, International Journal for Numerical Methods in Engineering 58: 999-1048.
• Allix, O., Corigliano, A., (1996). Modeling and simulation of crack propagation in mixed modes interlaminar fracture specimens, International journal of Fracture 77: 111-140.
• Asp, L.E., Sjögren, A., Greenhalgh, E.S., (2001). Delamination growth and thresholds in a carbon/epoxy composite under fatigue loading, Journal of Composites Technology & Research 23: 55-68.
• Baber, F., Guven, I., (2017). Solder joint fatigue life prediction using peridynamic approach, Microelectronics Reliability 79: 20-31.
• Bak, B.L.V., Turon, A., Lindgaard, E., Lund, E., (2016). A Simulation Method for High-Cycle Fatigue-Driven Delamination using a Cohesive Zone Model, International Journal for Numerical Methods in Engineering 106: 163-191.
• Bak, B.L.V., Turon, A., Lindgaard, E., Lund, E., (2017). A benchmark study of simulation methods for high-cycle fatigue-driven delamination based on cohesive zone models, Composite Structures 164: 198-206.
• Baydoun, M., Fries, T., (2012). Crack propagation criteria in three dimensions using the XFEM and an explicit-implicit crack description, International Journal of Fracture 178: 51-70.
• Bazazzadeh, S., Shojaei, A., Zaccariotto, M., Galvanetto, U., (2018). Application of the Peridynamic differential operator to the solution of sloshing problems in tanks, Engineering Computations (In press), DOI:10.1108/EC-12-2017-0520.
  » https://doi.org/10.1108/EC-12-2017-0520
• Belytschko, T., Black, T., (1999). Elastic crack growth in finite elements with minimal remeshing, International Journal for Numerical Methods in Engineering 45: 601-620.
• Bidokhti, A.A., Shahani, A.R., (2015). Interaction analysis of non-aligned cracks using extended finite element method, Latin American Journal of Solids and Structures 12: 2439-2459.
• Blackman, B.R.K., Hadavinia, H., Kinloch, A.J., Williams, J.G., (2003). The use of a cohesive zone model to study the fracture of fibre composites and adhesively-bonded joints, International journal of fracture 119: 25-46.
• Bobaru, F., Foster, J.T., Geubelle, P., Silling, S.A. (2016). Handbook of Peridynamic Modeling. Apple Academic Press Inc, Oakville, Canada.
• Bobaru, F., Ha, Y.D., (2011). Adaptive refinement and multiscale modeling in 2D peridynamics, International Journal for Multiscale Computational Engineering 9: 635-659.
• Boroomand, B., Bazazzadeh, S., Zandi, S.M., (2016). On the use of Laplace's equation for pressure and a mesh-free method for 3D simulation of nonlinear sloshing in tanks, Ocean Engineering 122: 54-67.
• Boroomand, B., Bazazzadeh, S., Zandi, S.M., (2017). On the use of Laplace's equation for pressure and a mesh-free method for 3D simulation of nonlinear sloshing in tanks; Reply to the discussion, Ocean Engineering 134: 176-177.
• Branco, R., Antunes, F., Costa, J.A., (2015). Review on 3D-FE adaptive remeshing techniques for crack growth modelling, Engineering Fracture Mechanics 141: 170-195.
• Chen, J., Crisfield, M.A., Kinloch, A.J., Busso, E.P., Matthews, F.L., Qiu, Y., (1999). Predicting progressive delamination of composite material specimens via interface elements, Mechanics of Composite Materials and Structures 6: 301-317.
• Chen, Z., Bobaru, F., (2015). Peridynamic modeling of pitting corrosion damage, Journal of the Mechanics and Physics of Solids 78: 352-381.
• Chen, Z., Zhang, G., Bobaru, F., (2016). The influence of passive film damage on pitting corrosion, Journal of The Electrochemical Society 163: C19-C24.
• Cheng, Z., Zhang, G., Wang, Y., Bobaru, F., (2015). A peridynamic model for dynamic fracture in functionally graded materials, Composite Structures 133: 529-546.
• Crisfield, M.A., Alfano, G., (2002). Adaptive hierarchical enrichment for delamination fracture using a decohesive zone model, International Journal for Numerical Methods in Engineering 54: 1369-1390.
• Dipasquale, D., Zaccariotto, M., Galvanetto, U., (2014). Crack propagation with adaptive grid refinement in 2D peridynamics, International Journal of Fracture 190: 1-22.
• Dolbow, J., Belytschko, T., (1999). A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering 46: 131-150.
• Galvanetto, U., Mudric, T., Shojaei, A., Zaccariotto, M., (2016). An effective way to couple FEM meshes and Peridynamics grids for the solution of static equilibrium problems, Mechanics Research Communications 76: 41-47.
• Galvanetto, U., Robinson, P., Cerioni, A., Armas, C.L., (2009). A Simple Model for the Evaluation of Constitutive Laws for the Computer Simulation of Fatigue-Driven Delamination in Composite Materials, SDHM: Structural Durability & Health Monitoring 5: 161-190.
• Ghajari, M., Iannucci, L., Curtis, P., (2014). A peridynamic material model for the analysis of dynamic crack propagation in orthotropic media, Computer Methods in Applied Mechanics and Engineering 276: 431-452.
• Gu, X., Zhang, Q., Xia, X., (2017). Voronoi‐based peridynamics and cracking analysis with adaptive refinement, International Journal for Numerical Methods in Engineering 112: 2087-2109.
• Ha, Y.D., Bobaru, F., (2010). Studies of dynamic crack propagation and crack branching with peridynamics, International Journal of Fracture 162: 229-244.
• Ha, Y.D., Bobaru, F., (2011). Characteristics of dynamic brittle fracture captured with peridynamics, Engineering Fracture Mechanics 78: 1156-1168.
• Hu, W., Ha, Y.D., Bobaru, F., (2012). Peridynamic model for dynamic fracture in unidirectional fiber-reinforced composites, Computer Methods in Applied Mechanics and Engineering 217: 247-261.
• Hu, W., Wang, Y., Yu, J., Yen, C.-F., Bobaru, F., (2013). Impact damage on a thin glass plate with a thin polycarbonate backing, International Journal of Impact Engineering 62: 152-165.
• Hu, Y.L., Carvalho, N.V.D., Madenci, E., (2015). Peridynamic modeling of delamination growth in composite laminates, Composite Structures 132: 610-620.
• Hu, Y.L., Madenci, E., (2017). Peridynamics for fatigue life and residual strength prediction of composite laminates, Composite Structures 160: 169-184.
• Jung, J., Seok, J., (2016). Fatigue crack growth analysis in layered heterogeneous material systems using peridynamic approach, Composite Structures 152: 403-407.
• Jung, J., Seok, J., (2017). Mixed-mode fatigue crack growth analysis using peridynamic approach, International Journal of Fatigue 103: 591-603.
• Juntti, M., Asp, L.E., Olsson, A., (1999). Assessment of evaluation methods for the mixed-mode bending test, Journal of Composites Technology & Research 21: 37-48.
• Le, Q.V. (2014) Relationship between microstructure and mechanical properties in Bi2Sr2CaCu2Ox round wires using peridynamic simulation. Ph.D thesis, Raleigh (NC): North Carolina State University.
• Lee, J., Oh, S.E., Hong, J.-W., (2017). Parallel programming of a peridynamics code coupled with finite element method, International Journal of Fracture 203: 99-114.
• Madenci E., Oterkus E., (2014). Peridynamic Theory and Its Applications, Springer-Verlag New York.
• Madenci, E., Barut, A., Futch, M., (2016). Peridynamic differential operator and its applications, Computer Methods in Applied Mechanics and Engineering 304: 408-451.
• Mi, Y., Crisfield, M.A. (1996). Analytical derivation of load/displacement relationship for the DCB and MMB and proof of the FEA formulation. IC-AERO Report 97-02, Department of Aeronautics, Imperial College, London, UK.
• Miranda, A.C.O., Meggiolaro, M.A., Castro, J.T.P., Martha, L.F., Bittencourt, T.N., (2003). Fatigue life and crack path predictions in generic 2D structural components, Engineering Fracture Mechanics 70: 1259-1279.
• Mossaiby, F., Shojaei, A., Zaccariotto, M., Galvanetto, U., (2017). OpenCL implementation of a high performance 3D Peridynamic model on graphics accelerators, Computers & Mathematics with Applications 74: 1856-1870.
• Ostergaard, M.G., Ibbotson, A.R., Roux, O.L., Prior, A.M., (2011). Virtual testing of aircraft structures, CEAS Aeronautical Journal 1: 83.
• Oterkus, E., Guven, I., Madenci, E. (2010) Fatigue failure model with peridynamic theory. 12th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems, Las Vegas, NV, USA.
• Peerlings, R.H.J., Brekelmans, W.A.M., Borst, R., Geers, M.G.D., (2000). Gradientenhanced damage modelling of high-cycle fatigue, International Journal for Numerical Methods in Engineering 49: 1547 -1569.
• Point, N., Sacco, E., (1996). A delamination model for laminated composites, International Journal of Solids and Structures 33: 483-509.
• Qiu, Y., Crisfield, M.A., Alfano, G., (2001). An interface element formulation for the simulation of delamination with buckling, Engineering Fracture Mechanics 68: 1755-1776.
• Ren, H., Zhuang, X., Cai, Y., Rabczuk, T., (2016). Dual‐horizon peridynamics, International Journal for Numerical Methods in Engineering 108: 1451-1476.
• Robinson, P., Galvanetto, U., Tumino, D., Bellucci, G., Violeau, D., (2005). Numerical simulation of fatigue-driven delamination using interface elements, International journal for numerical methods in engineering 63: 1824-1848.
• Schellekens, Borst, R., (1993). A non-linear finite element approach for the analysis of mode I free edge delamination in composites, International Journal of Solids and Structures 30: 1239-1253.
• Seleson, P., Beneddine, S., Prudhomme, S., (2013). A force-based coupling scheme for peridynamics and classical elasticity, Computational Materials Science 66: 34-49.
• Shi, J., Chopp, D., Lua, J., Sukumar, N., Belytschko, T., (2010). Abaqus implementation of extended finite element method using a level set representation for three-dimensional fatigue crack growth and life predictions, Engineering Fracture Mechanics 77: 2840-2863.
• Shojaei, A., Mudric, T., Zaccariotto, M., Galvanetto, U., (2016). A coupled meshless finite point/Peridynamic method for 2D dynamic fracture analysis, International Journal of Mechanical Sciences 119: 419-431.
• Shojaei, A., Zaccariotto, M., Galvanetto, U., (2017). Coupling of 2D discretized Peridynamics with a meshless method based on classical elasticity using switching of nodal behaviour, Engineering Computations 34: 1334-1366.
• Shojaei, A., Mossaiby, F., Zaccariotto, M., Galvanetto, U., (2018). An adaptive multi-grid peridynamic method for dynamic fracture analysis, International Journal of Mechanical Sciences 144: 600-617.
• Shojaei, A., Galvanetto, U., Rabczuk, T., Jenabi, A., Zaccariotto, M., (2019). A generalized finite difference method based on the Peridynamic differential operator for the solution of problems in bounded and unbounded domains, Computer Methods in Applied Mechanics and Engineering 343: 100-126.
• Silling, S.A., (2000). Reformulation of elasticity theory for discontinuities and long-range forces, Journal of the Mechanics and Physics of Solids 48: 175-209.
• Silling, S.A., Askari, A. (2014) Peridynamic model for fatigue cracking. Report, SAND2014-18590, Albuquerque (NM, United States): Sandia National Laboratories (SNL-NM).
• Silling, S.A., Askari, E., (2005). A meshfree method based on the peridynamic model of solid mechanics, Computers & Structures 83: 1526-1535.
• Silling, S.A., Epton, M., Weckner, O., Xu, J., Askari, E., (2007). Peridynamic states and constitutive modeling, Journal of Elasticity 88: 151-184.
• Warren, T.L., Silling, S.A., Askari, A., Weckner, O., (2009). A non-ordinary state-based peridynamic method to model solid material deformation and fracture, International Journal of Solids and Structures 464: 1186-1195.
• Weinert, A., Gergely, P., (2011). Fatigue-Strength Surface: basis for structural analysis under dynamic loads, CEAS Aeronautical Journal 2: 243-252.
• Xu, J., Askari, A., Weckner, O., Silling, S., (2008). Peridynamic analysis of impact damage in composite laminates, Journal of Aerospace Engineering 21: 187-194.
• Yu, K., Xin, X.J., Lease, K.B., (2011). A new adaptive integration method for the peridynamic theory, Modelling and Simulation in Materials Science and Engineering 19: 045003.
• Zaccariotto, M., Luongo, F., Sarego, G., Galvanetto, U., (2015a). Examples of applications of the peridynamic theory to the solution of static equilibrium problems, The Aeronautical Journal 119: 677-700.
• Zaccariotto, M., Mudric, T., Tomasi, D., Shojaei, A., Galvanetto, U., (2018). Coupling of FEM meshes with Peridynamic grids, Computer Methods in Applied Mechanics and Engineering 330: 471-497.
• Zaccariotto, M., Sarego, G., Dipasquale, D., Galvanetto, U. (2015b) Strategies for fatigue damage modeling with peridynamics. 8th International Congress of Croatian Society of Mechanics, Opatija, Croatia.
• Zaccariotto, M., Sarego, G., Dipasquale, D., Shojaei, A., Bazazzadeh, S., Mu-dric, T., Galvanetto, U. (2017). Discontinuous mechanical problems studied with a peridynamics-based approach. Aerotecnica Missili & Spazio, 96(1), 44-55.
• Zhang, G., Bobaru, F., (2016). Modeling the evolution of fatigue failure with peridynamics, Romanian Journal of Technical Sciences - Applied Mechanics, 61: 22-40.
• Zhang, G., Le, Q., Loghin, A., Subramaniyan, A., Bobaru, F., (2016). Validation of a peridynamic model for fatigue cracking, Engineering Fracture Mechanics 162: 76-94.

• Available online: January 04, 2019.

Publication Dates

• Publication in this collection
  2019

History

• Received
  09 Apr 2018
• Reviewed
  21 Nov 2018
• Accepted
  03 Jan 2019