Unification of local and nonlocal models within a stable integral formulation for analysis of defects

International Journal of Engineering Science 132 (2018) 45–59

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Unification of local and nonlocal models within a stable integral formulation for analysis of defects Mohsen Nowruzpour, J.N. Reddy∗ Department of Mechanical Engineering Texas A&M University, College Station, TX 77845, United States

a r t i c l e

i n f o

Article history: Received 19 May 2018 Accepted 14 June 2018

Keywords: Nonlocal continuum mechanics Non-ordinary state based peridynamics Fracture in rock Semi-circular bend specimen Mode-I Mode-II Mixed mode Fracture toughness MTS GMTS

a b s t r a c t Employing the discrete Cauchy-Born rule with the principle of virtual work done, a generalized model is formulated with the hope to unify local and nonlocal continuum frameworks. Despite the conventional mapping of the microscopic bond from the undeformed configuration, the consistent derivation requires a transformation on the Average Deviation of Lattice (ADL) vector in the region of influence. The new conversion proffers flexibility to the framework for the analysis of nonuniform distribution of particles in the field. We also found a compact mapping matrix which converts surface-based forces (stresses) to the nonlocal body-based forces. The transformation matrix allows reconstructing continuum models at a lower length scale in a discrete setting. To see the credibility of the model, fracture evolution in SCB specimen made of Polymethyl methacrylate (PMMA) is simulated, and the results are compared with experiments which admit an acceptable agreement. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction The recent advancements in developing new materials with a broader application made scholars revisit conventional constitutive models and frameworks. Besides, the description of empirical data challenges them to have a closer look at materials in a lower scale. Classical Continuum Theories (CCTs) allow a continuous spread of matter in the body and sets the equations of motion holding the local action solely. In additions, CCTs are based on partial differential equations (PDE). These assumptions narrow the validity of the CCT to those macro-level responses where the loading length scale is much larger than the physical length scales. But in the case of loading with a microscopic length scale, the classical prediction differs from experimental observation. Short wavelength excitations, analysis of porous media and state-of-the-art nanomaterial such as carbon nanotube are some examples that CCT failed to describe them accurately. Eringen (1976). Moreover, several experiments showed that mediums with smaller cracks have higher fracture resistance than a body with a larger one, while the CCT does not consider the effect of crack dimensions. To avoid such restrictions of the CCT, Voigt (1887) added a couple- forces to the conventional force-traction to model nonlocal interaction. Eringen verified that a nonlocal description is proficient at predicting a broad range of wavelengths Eringen (1972). The nonlocal theory was improved to be capable of predicting crack growth. Eringen pointed out that unlike the CCT, the stress distribution close to the crack tip is bounded. Eringen (1972) proposed a failure criterion by comparing cohesive stress to atomic bonds strength. Although the suggested nonlocal theory points the bounded stresses at the crack tips, the derivatives of the field variables are preserved ∗

Corresponding author. E-mail addresses: [email protected] (M. Nowruzpour), [email protected] (J.N. Reddy).

https://doi.org/10.1016/j.ijengsci.2018.06.008 0020-7225/© 2018 Elsevier Ltd. All rights reserved.

46

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

in the formulation. Eringen, Speziale, and Kim (1977) improved their nonlocal theory to model Griffith crack and later, he Ari and Eringen (1980) noted that the results of nonlocal Griffith crack model are in a good agreement with Elliott’s lattice model (1947). Nevertheless, the governing equations were written based on the partial differential equations which are still ill-defined on discontinuity. Generally, the nonlocality comes into the picture by adding strain derivative to the standard constitutive relation or defining strain averaging (Eringen et al., 1977). Despite other nonlocal theories which use derivative of field variable, Rogula (1982) proposed a nonlocal theory based on field variable. But the model was written for one-dimensional problems. Silling (20 0 0) proposed a derivative-free framework capable of analyzing multi-dimensional problems. The PD framework may be considered as an intermediate route between the classical and molecular dynamics (MD) approaches. Since it characterizes spatial interaction via integration, a PD equation can solve problems with discontinuities without resorting to any special treatment Silling (20 0 0). Based on this advantage, the PD framework finds its application not only in mechanics (Silling, 20 0 0) but also in areas like thermo-mechanics (Kilic & Madenci, 2010), electromigration (Gerstle, Silling, Read, Tewary, & Lehoucq, 2008), heat conduction in a body involving discontinuity (Bobaru & Duangpanya, 2012) etc. However, its original bond-based (BBPD) version faces a serious limitation because of its restriction on Poisson’s ratio. Besides, the BBPD does not distinguish between volumetric and distortional deformation. The reason behind such limitations in BBPD is traced back to its assumption of equal and opposite pairwise forces between two particles within a bond. As an important step forward, Silling came up with a modification of the BBPD formalism and proposed statebased peridynamics (SBPD) (more precisely ordinary SBPD) in Silling, Epton, Weckner, Xu, and Askari (2007), which could resolve many of the issues associated with the original BBPD approach. Unlike the BBPD, the forces in a bond are unequal in ordinary SBPD (Warren et al., 2009). However, the interaction forces within a bond are still considered as collinear. The SBPD framework is successfully applied in different areas of mechanics, e.g. plasticity (Silling et al., 2007), visco-elasticity (Kilic, 2008), visco-plasticity (Foster, Silling, & Chen, 2010; Taylor, 2008), dynamic brittle fracture (Ha & Bobaru, 2011), delamination in composite material (Xu, Askari, Weckner, & Silling, 2008), branching phenomena (Agwai, Guven, & Madenci, 2011) etc. But owing to its assumption of collinear forces along a bond, the ordinary SBPD is not applicable to non-linear anisotropic materials (Warren et al., 2009). Such limitation has led to further development, and non-ordinary SBPD has been proposed (Warren et al., 2009). Unfortunately, the non-ordinary SBPD is also scourged with difficulties in implementations. It may suffer from instability arising from the weak coupling of particles in the definition of deformation gradient. Responses via the non-ordinary SBPD may also show zero energy modes (Breitenfeld, Geubelle, Weckner, & Silling, 2014). Despite all the effort to improve laws and models (Chang, 2010; Ghosh & Arroyo, 2013; Maranganti & Sharma, 2007; Zhang, Jiao, Sharma, & Yakobson, 2006) a consistent connection between microscale and continuum level is not well defined. Most of the developed continuum models are not capable of capturing the evolution of microscopic defects such as crack within the framework. The existence of defects, such as imperfection and voids may affect the susceptibility of structure because of inception and growth of defects from an atomic level to the macroscopic scale. Several studies show that evolution of fracturing of brittle materials may not be presented via linear fracture mechanics (Buehler, van Duin, & Goddard III, 2006; Mattoni, Colombo, & Cleri, 2005; Mura, 1993; Remmers, de Borst, & Needleman, 2008; Yukutake, 1989). It started with bondbased peridynamics (Silling, Zimmermann, & Abeyaratne, 2003) which brought up some obstacles to have realistic modeling due to the extreme simplification of particle interactions in the domain. The concept of state in peridynamics framework was introduced to remove some of the restrictions (Silling et al., 2007; Silling & Lehoucq, 2010). However, the new concept required advanced mathematical tools to be implemented in the framework. Moreover, it has been shown that there is a huge instability in the response of the problem (Breitenfeld et al., 2014). Our recent study proposed a discrete Lagrangianbased framework to characterize the response of an elastic body at a length scale of interest (Sarkar, Nowruzpour, Reddy, & Srinivasa, 2017). A great benefit of adopting the Lagrangian description is its flexibility in characterizing coupled multibodyinteraction. In this paper, we present a systematic approach to constructing new nonlocal derivative-free formulation from classical constitutive models. This study could be an effort for the unification of local and nonlocal framework to analyze problems incorporating discontinuities and defects. In Section 2, we present nonlocal deformation in a body based on Discrete Cauchy-Born Rule (DCBR). The new definition compacts the state of deformation at a point of a discrete system into a second-order tensor. In Section 3, we derive a nonlocal derivative-free energetically-conjugate pair using nonlocal rate of work done in the continuum. In Section 4, a transformation matrix is introduced to map surface-based forces (conventional stress) to body-based ones. This transformation makes the continuum particle capable of interacting with the nonlocal region. In Section 5, using peridynamics approach for calculation of strain energy release rate (G), we present a new suitable energy-based criterion for the framework to predicts failure in brittle materials. In Section 6, a numerical investigation is done to show the credibility of the new development in analyzing mixed-mode fracture toughness of the semi-circular bend specimen made of PMMA. Finally, some concluding remarks are drawn in Section 7.

2. Directionality operator Here, we present a projection principle for upscaling microscopic information, which allows non-trivial directional information to develop a variable in a continuous or discrete system. The principal was obtained by employing a stochastic

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

47

projection method, uses certain microscopic details. The directionality operator can be defined as;

G = [ − ] · ϒ

(1)

ϒ = −1

(2)

The gain matrix, G might be described as a nonlocal equivalent of the deformation gradient that dictates the motion of a continuum body under external loading. Numerically it can be shown that such definition converges to conventional deformation gradient for the local limit. A discretized version of each term in Eq. (1) can be given as;



 = nc

m 

κ

n=1

nc

nc



1 = m  n=1

 ϒ = nc

vn

m 



u¯ ni

m 

−

κ

u¯ ci

 nc



xnj

  n − v ei  e j xcj

 m     κ nc xnj − xcj vn ei  e j u¯ ni − u¯ ci vn

n=1

κ

nc



xni

(3)

n=1

−

xci



xnj

 n − v ei  e j

(4)

−1

xcj

(5)

n=1

where, u¯ ni can be given; m 

1

u¯ ci =

m 

vn

n=1

κ nc uni vn

(6)

n=1

Superscripts n and c indicate the neighboring and the central particle. u and x show the displacement and location of particles. vn is the associated volume of the nth particle. The subscripts i and j are showing the direction of Cartesian coordinates. κ is the non-local smoothing function. A stabilized displacement field of the cth particle is defined in Eq. (6) that is dependent on the displacement of the neighbors to avoid weak coupling of particles; Such definition stabilizes the framework and avoids zero-mode energy in the solution. A continulized version of the ,  and ϒ can be given as;

i j = i j = ϒi j =



S

1 S  S

καi β j ds  S

καi ds

κβi β j ds





(7)  S

κβ j ds



(8)



(9)

where α i and β j are; 

αi = ui − ui

(10)



βj = xj − xj

Here, S is the scope or influence domain of the point of interest. Note that the prime symbol ( ) implies the points in the neighborhood of central one. 3. Nonlocal energetically-conjugate pair To get the new nonlocal energetically-conjugate pair suitable for the new framework, we go through a consistent energy approach. The rate of internal work done in a continuum in the current configuration can be equated to the rate of stored energy in the bonds connected from neighbors to the central particle. The statement might be mathematically expressed as;

W =

1 2

 V

(σ : d )dv =

 

1 .α˙ dS d; 2  S

d=

1 ∇ v + (∇ v )T 2

(11)

where d is the symmetric part of the velocity gradient tensor, σ is the Cauchy stress tensor,  is the force vector developed between the particle of interest and the neighboring particle. The lattice deformation rate is shown by α˙ . Since  creates the micro-potential density of a particle, and micro-potentials are defined with respect to a volume,  may be called BodyBased (BB) forces. On other hands, classical stresses might be renamed to Surface-Based (SB) forces, because they are always described with respect to a surface. Due to the symmetry of the stress tensor;

W =

1 2



V

σ : ld v =

 1 P : F˙ dv 2 

(12)

48

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

where l is velocity gradient and could be determined via l = F˙ F−1 ; Here F is deformation gradient and P is the local first Piola-Kirchhoff stress. To build up a nonlocal description of the deformation, we approximate F with the Gain function, G. Hence, the rate of work done in a continuum can be rewritten in a nonlocal setting as;

W =

 1 ψ : G˙ dV 2 

(13)

Accordingly, ψ might be interpreted as the nonlocal counterpart of First Piola–Kirchhoff Stress (NFPS). Since ϒ is calculated in the initial configuration, the time derivative of the Gain function can be given;





˙ − ˙ ϒ G˙ = 

(14)

For the simplicity of summation over same indices, we use Einstein’s notation;





  ˙ = 1 ˙ −  κ α˙ i γ j ds ei  e j ; S

S

γj = βj −

1



τ

τ

β j dτ 

(15)

where γ is a new lower scale component that naturally rises in the formulation. The parameter quantifies the Average Deviation of Lattice (ADL) vector in the region of influence. Unlike the usual practice of mapping of the undeformed bond vector into the deformed configuration (Breitenfeld et al., 2014; Mitchell, 2011; Warren et al., 2009), here the ADL parameter will be evolved into the new configuration. Therefore, in our derivation, γ may represent the notion of bond in a more generic manner. For a symmetric selection of τ , γ j is equal to β j (lattice vector) which is widely used in previous studies, mainly, Nonordinary State-Based Peridynamic (NSBP) (Amani et al., 2016; Breitenfeld et al., 2014; Wu & Ren, 2015). But in the case of selection of asymmetric influence domain, or an asymmetric distribution of particle in a symmetric scope like boundaries or discontinuities, γ j differs significantly from β j . Note that τ is independent of S and could be equal to or smaller than S. Using Eqs. (9) and (15) in Eq. (14) to get the rate of gain function and substituting the results in Eq. (13) gives;

W =

   



1 1 ψki ϒ jk κ α˙ i γ j dS dV = i α˙ i dS dV 2  2  S S

(16)

From Eq. (16), the structure of nonlocal BB force, , could be found. The transformation equation could be written in an indicial format as. Eq. (17)

i = ψki ϒ jk γ j

(17)

The obtained equation provides a connection from nonlocal surface-based (SB) interaction system to a nonlocal body-based (BB) one. It takes the undeformed ADL vector γ from the initial configuration with the nonlocal stress field ψ and returns the corresponding nonlocal BB force. 4. Body-based transformation matrix To facilitate the use of achieved results in previous section, we intend to prepare a matrix that could turn the nonlocalized derivative-free model from classical theory to the nonlocal BB model. Here, it should be noted that this framework cannot take every nonlocalized model, instead, the models that are constructed based on Gain function are suitable to be implemented in this framework. The vectorial form of the Eq. (17) can be presented as:

 T  T  = ϒ · ψ ·γ = ψ ·T ϒ T ·γ = γ T · ϒ · ψ

(18)

Given the fact that definition of ψ is built on the Gain function that converges to the classical counterpart in a localized view, to measure stresses in large deformation problems, NFPK could be written in terms of; ψ = G.ϕ where ϕ might be interpreted as Nonlocal Second Piola Kirchhoff (NSPK) stress. Rewriting Eq. (18) gives;

T = γ T · ϒ · G · ϕ

(19)

To find a transformation matrix capable of mapping of NSPK stress into BB forces, the following matrices are introduced;



T1

T2



T3 =



γ1

γ2

ϒ11 γ3 ϒ21 ϒ31

ϒ12 ϒ22 ϒ32

ϒ13 ϒ23 ϒ33



(20)

It is noteworthy to mention that T vector solely depends upon the initial configuration of the body. Matrix Tˆ that follows from the union of the T and G can be created as;



Tˆ1

Tˆ2





Tˆ3 = T1

T2

T3

G11

G21 G31

G12 G22 G32

G13 G23 G33



(21)

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

49

Fig. 1. Energetically representation of continuum element in a discrete model.

Fig. 2. Conversion of local SB forces to nonlocal BB forces.

using the Tˆ , Eq. (18) might be shown in matrix form;

⎡

 Tˆ1 1 2 = ⎣ 0 3 0 

0 Tˆ2 0

0 0 Tˆ3



0 Tˆ3 Tˆ2

Tˆ3 0 Tˆ1

Tˆ

⎡ ⎤ ϕ ⎤ ϕ 11 22 ⎥ Tˆ2 ⎢ ⎢ϕ33 ⎥ ⎥ Tˆ1 ⎦ ⎢ ⎢ϕ23 ⎥ 0 ⎣ϕ ⎦  13 ϕ12

(22)

Tˆ matrix might be named as a transformation matrix that converts nonlocal classical constitutive model to a generalized nonlocal derivative-free one. An evolution of SB force to BB forces are shown in Fig. 2. Note that 1 , 2 , 3 should not be confused with the resultant of stress tensor in the classical model. In fact, they are new quantities which are associated with the volume of the body and resultant of microscopic interaction with the local and nonlocal medium. This is in contrast to the definition of stress as a quantity that can only communicate through the surface (local definition of interaction). The given formulation can be directly applied to the orthotropic linear elastic material;

⎡

⎤ ⎡ ϕ11 C11 ⎢ϕ22 ⎥ ⎢C21 ⎢ϕ33 ⎥ ⎢C31 ⎢ ⎥=⎢ ⎢ϕ23 ⎥ ⎢ 0 ⎣ ⎦ ⎣ ϕ31 0 ϕ12 0

C12 C22 C32 0 0 0

C13 C23 C33 0 0 0

0 0 0 C44 0 0

0 0 0 0 C55 0

⎤⎡

⎤

0 ξ11 0 ⎥⎢ξ22 ⎥ ⎢ ⎥ 0 ⎥ ⎥⎢ξ33 ⎥ 0 ⎥⎢ξ23 ⎥ ⎦⎣ ⎦ 0 ξ31 C66 ξ12

(23)

50

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

Fig. 3. Calculation of nonlocal strain energy dissipated during formation of unite fracture surface area.

Fig. 4. Areas representing particles with broken bonds connected to central particles lying on green axis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

where ξ is nonlocal derivative-free Green strain and it could be written based upon G as follows;

ξ=

 1 T G ·G−I 2

(24)

Ultimately, using the equation of motion (EOM) for a two-body interaction system, developed in the previous study, Sarkar et al. (2017), EOM could be rewritten in the following form:

φ + f = ρ u¨ ; φ =



 S

  −  dS

(25)

where ρ , f and u are density, body force vector and displacement vector respectively.  is generalized nonlocal body-based force vector. One of the major distinctions of such a framework with classical continuum theory is that the new nonlocal forces are not any longer interacting through the surface, since they are volume-dependent quantity and they can fix the restriction of reaching nonlocal region. 5. Failure criteria In this section, we try to develop an energy-based failure criterion for the framework. Zhou and Wang (2016) provided a stress-based failure criterion in a discrete system. They proposed an extra criterion in addition to the tensile failure to incorporate shear failure of the bonds. Silling and Askari (2005) calculated the summation of the required work per unit length to eliminate every interaction between particles on the symmetric line below the crack surface (Pi ) and particles above the crack surface (Pj ) that are in the influence domain of particle (Pi ) (see Fig. 3). Using this calculation, an energybased failure approach which relaxes the need for an extra criterion is proposed. The strain energy release rate can be calculated once all the critical failure energy of bonds that cross the crack surface are summed up (see colored area in Fig. 4). Note that as we move downward, the area above the crack gets smaller. For a thin linear isotropic elastic material, the strain energy can be calculated via:



U=

S

U ∗ dv =

1 2



S

ϕ i j ξi j d v

(26)

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

51

Fig. 5. Scheme of SCB specimen with inclination crack angle, α .

where U∗ is strain energy density. Note that ξ ij is the Green strain. For a two dimensional analysis, strain energy release rate, G can be calculated as follows;

G = U nc =

t 2



y=rs

 θ =α 

y=0

r=rs

θ =−α r= cosy(θ )



 ϕincj ξincj rdrdθ dy

(27)

 where t is the thickness. Assuming a symmetric distribution of particles i.e. τ1 τ β j dτ  = 0 ⇒ γ j = β j , the strain energy release rate for particles under the crack surface could be;

⎛ G

nc

=

y

t nc nc 2 ϕ ξ ⎝rs y cos−1 2 ij ij rs

−

rs y2



1−

3

y2 rs2

−

2 rs3



1−

3

y2 rs2

⎞y=rs ⎠

=

y=0

=

t 3 nc nc r ϕ ξ 3 s ij ij

(28)

t 3 nc nc nc nc nc nc c nc c nc nc nc r [ϕ ξ + ϕ22 ξ22 + ϕ33 ξ33 + 2(ϕ23 ξ23 + ϕ31 ξ31 + ϕ12 ξ12 )] 3 s 11 11

Here, rs is a nonlocal parameter that emerges in the definition of the nonlocal strain energy release rate. Since the calculation of energy can be repeated for particles on the line above the crack surface (i.e., Gcn ), the total strain energy release rate can be given;

G = f1 (γ nc )G

nc

+ f2 (γ cn )G

cn

(29)

In general, we are differentiating between the connections from n to c particle and c to n particle. To show this difference, we use f1 and f2 function as a coefficient of G. In this paper, for f1 = f2 = 1, we are still able to get acceptable results. Here, the failure criterion states that the connection between two bonds fails if the strain energy calculated based upon two ends of a bond exceeds the critical amount. 6. Results and discussion 6.1. Calculation of mode I/II fracture toughness using SCB specimens Kuruppu and Chong (2012) suggested semi-circular bend (SCB) specimen for measuring mode I of brittle materials such as rock or other geological material. Many factors such as convenient of running test, minimal machining labor for preparation of specimens and stability of results made it a popular test specimen. Fig. 5 depicts a scheme of the SCB of radius R which is resting on two supports. The length of the distance from one support to another is 2S. Each specimen has a crack of length a, that creates an angle with the line of symmetry. With different crack inclination (α ), different crack modes or a mixture of them can be obtained. For instance, when α = 0, the specimen is subjected to pure mode I. As we increase α ,

52

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

Fig. 6. Geometry factor values of mode I and mode II for different inclination angles.

Fig. 7. Algorithm of Verlet-Velocity method.

mode II contributes more to the failure of specimen. To calculate the stress state close to the tip of the crack, stress intensity factor in mode I and II are given;



YI

YII

a S R R

α, ,



a S R R

α, ,



KI = √

2Rt

(30)

π a Pcr

KII 2Rt = √ π a Pcr

(31)

where t and Pcr are thickness and critical applied force on the top of specimen respectively. YI , YII are geometry factors associated with mode I and mode II. The geometry factor parameters are a function of crack angle (α ), crack length ratio, a/R, and support length ratio S/R. Ayatollahi and Aliha (2007) has provided a wide range of experimental data of YI and YII for various geometrical ratios. Fig. 6 shows values of YI and YII for a SCB with dimensional ratio of a/R = 0.3 and S/R = 0.43 for different inclination angles (Ayatollahi, Aliha, & Hassani, 2006) 6.2. Numerical method; Verlet scheme To solve equation of the motion numerically, the body is discretized into number of nodes. a discretized form of the governing equations could be given as;

φc + f c = ρ u¨ c ; φc =  T c = ϒ c · ψ c ·γ c  T n = ϒ n · ψ n ·γ n

m  

 c − n vn

(32)

n=1

(33) (34)

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

53

Fig. 8. Flowchart of the developed code.

Fig. 9. Semi-Circular Bend specimen under mode I (α = 0) or mixed mode loading (α = 0).

γ nc = βnc −

m 1  nc n β v S

(35)

n=1

A flowchart of the developed code is shown in Fig. 8. Since we have a discrete structure, we try to pick one of the powerful methods for numerical integration of the equation of motion. Verlet scheme is a fascinating method particularly for researchers in molecular dynamics area. The method is offered in three version; Velocity, position (Verlet, 1967) and Leapfrog (Swope, Andersen, Berens, & Wilson, 1982) type. This algorithm provides acceptable stability along with simplicity. To update position at time t and velocity at time (t + 12 δt ) and t, Taylor expansion series is employed:

X(t + δt ) = X(t ) + X˙ (t )δt +

1 ¨ (t )δt 2 + · · · X 2

(36)

54

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

Fig. 10. Variation of critical load with crack inclination angles for different rs . Table 1 Mechanical properties of PMMA.



X˙ t +

Parameter

Value

Units

E (Young’s Modulus) ν (Poisson’s ratio) ρ (Density) KIc (Fracture Toughness I)

3.75 0.35 1.41 2.13

GPa – Mg/m3 √ MPa m



1 1 δt = X˙ (t ) + X¨ (t )δt + · · · 2 2

φ (t ) ρ 1 1 ¨ (t + δt )δt + · · · X˙ (t + δt ) = X˙ t + δt + X ¨ (t + δt ) = X

2

2

(37) (38) (39)

The new configuration can be updated using position, velocity, and acceleration at the previous step (see Fig. 7, a). The nonlocal forces can be calculated by updating new position vectors (see Fig. 7, b). Finally, to get velocity at time (t + δt ), the velocity at the middle of time increment is determined (Fig. 7, c), and then the velocity can be determined from the updated position and velocity at previous half-step. 6.3. Numerical simulation of dynamic crack growth and fracture analysis of SCB rock specimen Failure of engineering parts due to crack propagation is a topic that attracts much attention. Maiti and Prasad (1980) studied the prediction of the unstable crack path using the stress distribution before the crack starts propagating in the body. Oliver, Dias, and Huespe (2014) developed a numerical technique by injecting a discontinuous displacement to simulate the evolution of crack in materials. In this model, the evolution of crack is investigated without applying any external criteria to dictate crack path. As a benchmark problem, we implemented a numerical simulation of a cracked SCB made of Polymethylmethacrylate (PMMA) with the different inclination angles. The mechanical properties of PMMA are given in Table 1. The domain is discretized into 16,366 particle uniformly (x = y) with the particle spacing of x = 0.0 0 05. Fig. 9 depicts a discretized SCB with constraints and loading condition. The SCB is constrained vertically on both supports, and a constant displacement rate is applied to the top of the specimen. For the sake of stability and avoiding fracture at supports, all constraints and loading (on the top) are applied to an area of 3x × 3x. The following dimensions and loading condition are considered for the simulation;

R = 50 mm, a = 0.3R, S = 0.43R, t = 5 mm,

t = 5 × 10−8 s, Vy (0, R ) = −0.005 mm/min, ux (±S, 0 ) = 0.0

(40)

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

55

Fig. 11. Crack trajectory in SCB specimen made of PMMA with initial crack and inclination angle of α = 0◦ , 10◦ , 20◦ , 30◦ .

Since rs is a new nonlocal parameter, it can be determined by comparison of the fracture results for different rs with the experimental values. Fig. 10 displays the critical load associated with different crack angles. It can be seen that for rs = 2x, the simulation’s results have a closer prediction to experiment. Therefore, the remaining results are presented based on rs = 2x = 0.001 . It is noteworthy to mention that rs is not a mesh dependent parameter, but it is material property, and it has a unique value for each material. Figs. 11 and 12 illustrate the final status of specimens after the crack is fully grown. The color bar in these figures shows the extent of damage to a particle. The damage index of one implies that a particle has missed all its links, whereas the index of zero indicates a particle that preserves all of its connections. This parameter is given as:

μ=1−

R T

(41)

where R is the number of remaining bonds and T is the total number of connected bonds before the loading starts. When α = 0.0◦ (pure mode I), The crack starts propagating on a straight line, along the direction of initial crack till it reaches the top of the SCB. As the inclination angle increases, cracks trajectory deviates from the direction of initial flaw due to the contribution of mode II loading. It can be clearly seen for angles greater than 10°, the propagation starts from the tip, and it eventually aligns with the symmetry line of semi-circle. This is in agreement with the Ayatollahi et al.’s experiments

56

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

Fig. 12. Crack trajectory in SCB specimen made of PMMA with initial crack and inclination angle of α = 40◦ , 43◦ , 47◦ , 50◦ .

(2006) and Xie, Cao, Jin, and Wang’s reports (2017) for SCB with length ratio a/r = 0.3. The mean values of experimental and theoretical fracture toughness along with their corresponding critical load are listed in Table 2. The magnitude of the error for KIC , KIIC and P are relatively low, except at α = 20, 43. According to the trend of Pc , it is expected to get a value more than 2.53 kN at α = 20, but the experimental data is showing lower magnitude. Such inconsistency might be a result of manufacturing flaws or measurement errors. Note that Ayatollahi et al. (2006) has dropped these values from their graphical results to get a decent trend in their curve fitting. Figs. 13 and 14 are showing the theoretical and experimental fracture toughness of mode I and II for each specimen. The results are compared with the maximum tangential stress (MTS) method (Erdogan & Sih, 1963) and generalized MTS (GMTS) (Smith, Ayatollahi, & Pavier, 2001) that are commonly used in the analysis of mixed-mode fracture. Note that MTS only takes the first term (singular term) of the series expansion of tangential stress around the crack tip. However, GMTS includes singular and nonsingular terms of the series to predict the start of crack propagation. In both developments, the crack propagates along a direction perpendicular to the maximum of the tangential stress, σ θ . The simulation from current study predicts mode I fracture toughness at α = 10◦ , 20◦ , 40◦ with better accuracy (see Fig. 13). It also performs better in determining the magnitude of mode II fracture toughness at 20°, 30°, 40° in comparison to both conventional and generalized MTS (Ayatollahi et al., 2006) (see Fig. 14). It is common practice to present the fracture toughness results from mixed mode loading in terms of normalized fracture toughness i.e. KI /KIc and KII /KIc . The mean value of experimental data (Ayatollahi et al., 2006) for pure mode I fracture √ toughness of PMMA is 2.1258MPa m. In Fig. 15, the simulation results of mode I and mode II fracture toughness are plotted along with the experimental data, MTS and GMTS simulation. Each of the green squares is showing an initial crack incli-

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

57

Fig. 13. Comparison of the mode I normalized fracture toughness of different crack angle obtained from experiment and conventional criteria with the predicted results from the present study.

Table 2 Experimental and numerical fracture results for SCB with a length ratio of S/R = 0.43 and a/R = 0.3 subjected to a three-point bending. Experimental Ayatollahi et al. (2006)

Theorical

α (°)

KI

KII

Pc (kN)

KI

KII

Pc (kN)

KI

Error % KII

Pc (kN))

0 10 20 30 40 43 47 50

2.13 2.10 1.63 1.30 0.75 0.53 0.23 0.00

0.00 0.47 0.75 1.07 1.31 1.22 1.24 1.12

2.38 2.53 2.45 3.03 3.73 3.63 4.13 4.24

2.27 2.22 1.91 1.32 0.75 0.53 0.22 0.00

0.00 0.47 0.87 1.13 1.33 1.33 1.32 1.29

2.54 2.68 2.85 3.09 3.79 4.00 4.39 4.73

6.74 5.96 16.44 1.73 1.45 10.13 6.22 11.65

6.66 5.95 17.12 1.24 0.04 0.97 1.61 0.00

0.00 0.05 16.30 4.73 1.33 9.02 6.47 15.74

Fig. 14. Comparison of the mode II normalized fracture toughness of different crack angle obtained from experiment and conventional criteria with the predicted results from the present study.

58

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

Fig. 15. Comparison of mode I/II normalized fracture toughness from the experiment and conventional criteria with the predicted results from the present study.

nation angle (given in the table. 2) that starts from zero as we move on the curve from right to left. The red dash line is showing the average of the experimental value (red dots) for mode I and mode II normalized fracture toughness. It can be seen from the figure that the current study has better estimation than other two methods (MTS, GMTS), especially, once the specimen is undergoing pure mode II (i.e., KI = 0). 7. Conclusion In this study, a new nonlocal derivative-free framework based on the discrete Cauchy-Born rule was formulated. The new approach granted a systematic procedure to generate all nonlocal version of classical continuum definitions such as deformation gradient, strain, and stresses from the local counterparts. For such a purpose, a compact matrix was constructed to transform surface-based forces to body-based ones. The new integral form of equation enables us to investigate mechanical problems having a discontinuity in the field variables. Instead of using conventional definition of bond, we introduced a more generic term to represent micro-interaction of particles. As an application of this study, a mixed-mode fracture analysis of PMMA semi-circular bend specimens was done. The results showed a remarkable improvement in prediction of fracture toughness in comparison to MTS method. Although GMTS performed much better than MTS, the simulation results showed that the new development still has a better determination of mode II fracture toughness compared to GMTS. Finally, a comparison of the experimental observation with the simulated crack trajectory confirms an acceptable agreement. Acknowledgement The authors gratefully acknowledge the support of this research by the Oscar S. Wyatt Endowed Chair. Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ijengsci.2018.06. 008. References Amani, J., Oterkus, E., Areias, P., Zi, G., Nguyen-Thoi, T., & Rabczuk, T. (2016). A non-ordinary state-based peridynamics formulation for thermoplastic fracture. International Journal of Impact Engineering, 87, 83–94. Agwai, A., Guven, I., & Madenci, E. (2011). Predicting crack propagation with peridynamics: A comparative study. International Journal of Fracture, 171(1), 65–78. Ari, N., & Eringen, A. C. (1980). Nonlocal Stress Field at Griffith Crack.. Technical Report. Princeton Univ NJ Dept of Civil and Geological Engineering. Ayatollahi, M., Aliha, M., & Hassani, M. (2006). Mixed mode brittle fracture in pmmaan experimental study using scb specimens. Materials Science and Engineering: A, 417(1–2), 348–356. Ayatollahi, M., & Aliha, M. (2007). Wide range data for crack tip parameters in two disc-type specimens under mixed mode loading. Computational Materials Science, 38(4), 660–670. Bobaru, F., & Duangpanya, M. (2012). A peridynamic formulation for transient heat conduction in bodies with evolving discontinuities. Journal of Computational Physics, 231(7), 2764–2785.

M. Nowruzpour, J.N. Reddy / International Journal of Engineering Science 132 (2018) 45–59

59

Breitenfeld, M., Geubelle, P., Weckner, O., & Silling, S. (2014). Non-ordinary state-based peridynamic analysis of stationary crack problems. Computer Methods in Applied Mechanics and Engineering, 272, 233–250. Buehler, M. J., van Duin, A. C., & Goddard III, W. A. (2006). Multiparadigm modeling of dynamical crack propagation in silicon using a reactive force field. Physical Review Letters, 96(9), 095505. Chang, T. (2010). A molecular based anisotropic shell model for single-walled carbon nanotubes. Journal of the Mechanics and Physics of Solids, 58(9), 1422–1433. Elliott, H. (1947). An analysis of the conditions for rupture due to griffith cracks. Proceedings of the Physical Society, 59(2), 208. Erdogan, F., & Sih, G. (1963). On the crack extension in plates under plane loading and transverse shear. Journal of basic engineering, 85(4), 519–525. Eringen, A. C. (1972). Linear theory of nonlocal elasticity and dispersion of plane waves. International Journal of Engineering Science, 10(5), 425–435. Eringen, A. C. (1976). Continuum physics. Vol. iii: Mixtures and em field theories. Academic Press. Eringen, A. C., Speziale, C., & Kim, B. (1977). Crack-tip problem in non-local elasticity. Journal of the Mechanics and Physics of Solids, 25(5), 339–355. Foster, J. T., Silling, S. A., & Chen, W. W. (2010). Viscoplasticity using peridynamics. International Journal for Numerical Methods in Engineering, 81(10), 1242–1258. Gerstle, W., Silling, S., Read, D., Tewary, V., & Lehoucq, R. (2008). Peridynamic simulation of electromigration. Comput Mater Continua, 8(2), 75–92. Ghosh, S., & Arroyo, M. (2013). An atomistic-based foliation model for multilayer graphene materials and nanotubes. Journal of the Mechanics and Physics of Solids, 61(1), 235–253. Ha, Y. D., & Bobaru, F. (2011). Characteristics of dynamic brittle fracture captured with peridynamics. Engineering Fracture Mechanics, 78(6), 1156–1168. Kilic, B. (2008). Peridynamic theory for progressive failure prediction in homogeneous and heterogeneous materials. ProQuest. Kilic, B., & Madenci, E. (2010). Peridynamic theory for thermomechanical analysis. Advanced Packaging, IEEE Transactions on, 33(1), 97–105. Kuruppu, M. D., & Chong, K. P. (2012). Fracture toughness testing of brittle materials using semi-circular bend (scb) specimen. Engineering Fracture Mechanics, 91, 133–150. Maiti, S., & Prasad, K. (1980). A study on the theories of unstable crack extension for the prediction of crack trajectories. International Journal of Solids and Structures, 16(6), 563–574. Maranganti, R., & Sharma, P. (2007). A novel atomistic approach to determine strain-gradient elasticity constants: Tabulation and comparison for various metals, semiconductors, silica, polymers and the (ir) relevance for nanotechnologies. Journal of the Mechanics and Physics of Solids, 55(9), 1823–1852. Mattoni, A., Colombo, L., & Cleri, F. (2005). Atomic scale origin of crack resistance in brittle fracture. Physical Review Letters, 95(11), 115501. Mitchell, J. A. (2011). A Nonlocal Ordinary State-Based Plasticity Model for Peridynamics.. Technical Report. Sandia National Laboratories (SNL-NM), Albuquerque, NM (United States). Mura, T. (1993). Micromechanics: overall properties of heterogeneous solids. Oliver, J., Dias, I., & Huespe, A. E. (2014). Crack-path field and strain-injection techniques in computational modeling of propagating material failure. Computer Methods in Applied Mechanics and Engineering, 274, 289–348. Remmers, J. J., de Borst, R., & Needleman, A. (2008). The simulation of dynamic crack propagation using the cohesive segments method. Journal of the Mechanics and Physics of Solids, 56(1), 70–92. Rogula, D. (1982). Introduction to nonlocal theory of material media. In Nonlocal theory of material media (pp. 123–222). Springer. Sarkar, S., Nowruzpour, M., Reddy, J., & Srinivasa, A. (2017). A discrete lagrangian based direct approach to macroscopic modelling. Journal of the Mechanics and Physics of Solids, 98, 172–180. Smith, D., Ayatollahi, M., & Pavier, M. (2001). The role of t-stress in brittle fracture for linear elastic materials under mixed-mode loading. Fatigue & Fracture of Engineering Materials & Structures, 24(2), 137–150. Silling, S. A. (20 0 0). Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 48(1), 175–209. Silling, S. A., Zimmermann, M., & Abeyaratne, R. (2003). Deformation of a peridynamic bar. Journal of Elasticity, 73(1–3), 173–190. Silling, S. A., & Askari, E. (2005). A meshfree method based on the peridynamic model of solid mechanics. Computers & structures, 83(17–18), 1526–1535. Silling, S. A., Epton, M., Weckner, O., Xu, J., & Askari, E. (2007). Peridynamic states and constitutive modeling. Journal of Elasticity, 88(2), 151–184. Silling, S., & Lehoucq, R. (2010). Peridynamic theory of solid mechanics. Advances in Applied Mechanics, 44, 73–168. Swope, W. C., Andersen, H. C., Berens, P. H., & Wilson, K. R. (1982). A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: application to small water clusters. The Journal of Chemical Physics, 76(1), 637–649. Taylor, M. J. (2008). Numerical simulation of thermo-elasticity, inelasticity and rupture in membrane theory. University of California, Berkeley. Verlet, L. (1967). Computer” experiments” on classical fluids. i. thermodynamical properties of lennard-jones molecules. Physical Review, 159(1), 98. Voigt, W. (1887). Theoretische studien fiber die elastizit ∼ itsverh ∼ iltnisse der kristalle. Abh. Geschichte Wissenschaften, 34. Warren, T. L., Silling, S. A., Askari, A., Weckner, O., Epton, M. A., & Xu, J. (2009). A non-ordinary state-based peridynamic method to model solid material deformation and fracture. International Journal of Solids and Structures, 46(5), 1186–1195. Wu, C., & Ren, B. (2015). A stabilized non-ordinary state-based peridynamics for the nonlocal ductile material failure analysis in metal machining process. Computer Methods in Applied Mechanics and Engineering, 291, 197–215. Xie, Y., Cao, P., Jin, J., & Wang, M. (2017). Mixed mode fracture analysis of semi-circular bend (scb) specimen: A numerical study based on extended finite element method. Computers and Geotechnics, 82, 157–172. Xu, J., Askari, A., Weckner, O., & Silling, S. (2008). Peridynamic analysis of impact damage in composite laminates. Journal of Aerospace Engineering, 21(3), 187–194. Yukutake, H. (1989). Fracturing process of granite inferred from measurements of spatial and temporal variations in velocity during triaxial deformations. Journal of Geophysical Research: Solid Earth, 94(B11), 15639–15651. Zhang, X., Jiao, K., Sharma, P., & Yakobson, B. (2006). An atomistic and non-classical continuum field theoretic perspective of elastic interactions between defects (force dipoles) of various symmetries and application to graphene. Journal of the Mechanics and Physics of Solids, 54(11), 2304–2329. Zhou, X.-P., & Wang, Y.-T. (2016). Numerical simulation of crack propagation and coalescence in pre-cracked rock-like brazilian disks using the non-ordinary state-based peridynamics. International Journal of Rock Mechanics and Mining Sciences, 89, 235–249.