An extended peridynamic approach for deformation and fracture analysis

Accepted Manuscript An extended peridynamic approach for deformation and fracture analysis Dan Huang, Guangda Lu, Chongwen Wang, Pizhong Qiao PII: DOI: Reference:

S0013-7944(15)00203-9 http://dx.doi.org/10.1016/j.engfracmech.2015.04.036 EFM 4603

To appear in:

Engineering Fracture Mechanics

Received Date: Revised Date: Accepted Date:

28 August 2014 16 February 2015 30 April 2015

Please cite this article as: Huang, D., Lu, G., Wang, C., Qiao, P., An extended peridynamic approach for deformation and fracture analysis, Engineering Fracture Mechanics (2015), doi: http://dx.doi.org/10.1016/j.engfracmech. 2015.04.036

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An extended peridynamic approach for deformation and fracture analysis Dan Huanga, Guangda Lua, Chongwen Wanga, Pizhong Qiaob,c,∗ a

Department of Engineering Mechanics, Hohai University, Nanjing, 210098, P. R. China

b

State Key Laboratory of Ocean Engineering and School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai, 200240, P.R. China

c

Department of Civil and Environmental Engineering, Washington State University, Pullman, WA, 99164-2910, USA

A short and self-contained abstract: An extended bond-based peridynamic approach for the quantitative analysis of elastic deformation and simulation of both quasi-static and dynamic crack propagation is proposed by uniquely introducing a continuous function in the commonly used peridynamic constitutive model and considering the internal length effect of long-range forces. A local damping is also incorporated into the peridynamic equations of motion, and the step-by-step loading algorithms and non-equilibrium criterion for particle systems are developed specifically for quasi-static problems.

Analysis of quantitative

elastic deformation and quasi-static fracture behavior of the model is verified and demonstrated through simulating several common tests of beams and plates.

∗

Corresponding author at: Department of Engineering Mechanics, Shanghai Jiao Tong University,

Mulan Building A903, Shanghai, 200240, PR China. E-mail addresses: [email protected]; [email protected] (P. Qiao). 1

An extended peridynamic approach for deformation and fracture analysis Dan Huanga, Guangda Lua, Chongwen Wanga, Pizhong Qiaob,c,∗ a

Department of Engineering Mechanics, Hohai University, Nanjing, 210098, P. R. China

b

State Key Laboratory of Ocean Engineering and School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai, 200240, P.R. China

c

Department of Civil and Environmental Engineering, Washington State University, Pullman, WA, 99164-2910, USA

Abstract: An extended bond-based peridynamic approach is proposed for quantitative analysis of elastic deformation and simulation of both quasi-static and dynamic crack propagation.

A continuous function is uniquely introduced in the commonly used

peridynamic constitutive model by considering the internal length effect of long-range forces. Moreover, a local damping is incorporated into the peridynamic equations of motion, and the step-by-step loading algorithms and non-equilibrium criterion for particle systems are developed specifically for quasi-static problems.

The effective calculation field of material

points in the numerical peridynamic discretization is discussed as well.

Analysis of

quantitative elastic deformation and quasi-static fracture behavior of the model is verified through simulating the uniaxial tension and compression of a slab, bending and failure of a cantilever beam, and quasi-static propagation of mode I and mode II cracks in a plate. The numerical examples show that the present approach is capable of quantitatively calculating the elastic deformation and accurately predicting the critical failure load and crack

∗

Corresponding author at: Department of Engineering Mechanics, Shanghai Jiao Tong University,

Mulan Building A903, Shanghai, 200240, PR China. E-mail addresses: [email protected]; [email protected] (P. Qiao). 2

propagation paths. To further demonstrate the capabilities of the proposed model, dynamic fracture of a brittle plate under high-rate loading is simulated, in which the experimentally-observed phenomena, such as successive branching, are captured naturally as a consequence of the solution. Keywords: quasi-static response; dynamic fracture; crack branching, peridynamic theory.

1. Introduction Peridynamics (PD) is a nonlocal formulation of continuum mechanics originally proposed to handle problems involving discontinuities and long-range forces [1].

In

contrast to the classical theories of continua but also to other nonlocal approaches which require additional treatment or remedy methods across discontinuities such as cracks and interfaces, the peridynamic equation of motion is in integral form and free of any spatial derivative, which makes it a promising way to analyze the failure of materials and structures by spontaneously predicting the crack nucleation and propagation with arbitrary paths without any extra numerical techniques. With the unique capability of naturally describing discontinuities within a single continuum framework [2], the peridynamic approach has been successfully employed in a variety of problems in solid mechanics at both macro- and micro-scales, including deformation of a bar [3,4], stationary and dynamic fracture [5-7], structural stability [8], phase transformation [9], heat transfer [10-12], crack propagation and branching in brittle materials [13-15], and mechanical response of nano-scaled materials and composites [16,17]. Further researches on peridynamic constitutive modeling for nonlinear behavior of materials [18-21] and coupling peridynamic model with molecular dynamics (MD) and finite element method (FEM) [22-24] were conducted as well. Peridynamic formulations may be generally divided into two categories, i.e., the

3

so-called bond-based formulation which is comparatively simple but limited to model materials with the fixed Poisson’s ratio, and the state-based formulation, which makes it possible to incorporate very general material models.

In the recent years, more focuses were

paid on the latter formulation, especially on the so-called non-ordinary state-based peridynamic modeling which appears to allow for far more general constitutive responses [2, 25]. Classical nonlinear constitutive models (e.g., viscoplasticity [26] and ductile damage [27]) have been incorporated to non-ordinary state-based peridynamic formulation recently. Nevertheless, the bond-based PD model is still widely used as it is comparatively much simple to understand and implement, and it is more suitable for coupling with the molecular dynamic models in multi-scale analysis as well.

When coupling the bond-based PD models

with the finite element models, the pair-wise bonds in PD model can be treated as the beam elements directly. On the other hand, there still appears to be many interesting phenomena that can be investigated by using the bond-based PD formulation without the enhancement, especially in the linear elastic deformation and brittle fracture analysis.

Meanwhile, several

aspects of the mathematical framework of the bond-based PD formulation, including the constitutive model and the numerical algorithms, are still waiting for further researches to obtain higher efficiency and accuracy. For example, more attention is paid to dynamic fracture simulation, but seldom to quantitative quasi-static deformation and fracture analysis. In reality, the external loads, which lead to elastic deformation and crack initiation and growth, are usually given in term of velocities or displacements of particles in most of the earlier PD simulations.

What’s more, the boundary bonds are usually set to be no-fail zones

in order to prevent tearing of the boundary particles from the rest nodes of the model, and the micro-modulus of the pair-wise bond is generally set to be a constant, regardless of the distance between two particles within the horizon. In the present contribution, the aforementioned problems is solved by introducing a

4

continuous function to the expression of bond micro-modulus and implementing the new numerical algorithms including local damping, step-by-step loading, and a non-equilibrium criterion for the bond-based PD simulation. This paper proceeds as follows. Section 2 briefly reviews the bond-based peridynamic theory.

Section 3 presents an improved

bond-based PD constitutive model and numerical approaches including the local damping, step-by-step external force loading, discretization and convergence implementations. Section 4 demonstrates the extended model with typical numerical examples including quantitative elastic deformation calculation, quasi-static propagation of mode I and mode II cracks, and dynamic crack branching.

Finally, Section 5 summarizes the conclusions

resulted from this study. 2. Formulation of peridynamics The peridynamic method is in the category of non-local models, but it is different from the conventional non-local methods by replacing the spatial differential equations with the integral operations, rather than using the average of the derivatives, thus eliminating the assumption on the continuity or differentiability of the displacement field thoroughly. The peridynamic equation of motion at a material point x and time t according to the Newton’s law is given as

ρ ( x )u( x, t ) = ∫ f ( x, x′, u( x, t ) , u( x′, t ) , t )dV ′ + b( x, t ) = Lu ( x, t ) + b( x, t ) Η

(1)

where u denotes the displacement of a material point, ρ is the mass density, b represents the external force per unit reference volume. H is a neighborhood of point x which is usually taken to be a sphere in 3D or a circle in 2D problems centered at x , and it is given as

H = H ( x, δ ) = { x′ ∈ R : x′ − x ≤ δ }

(2)

where δ is the horizon, which refers to the “size” of nonlocal interaction. It indeed plays 5

the role of a length-scale parameter, and one recovers the classical theory of nonlinear elasticity in the limit of vanishing length scale [28]. The peridynamic equation of motion in Eq. (1) reduces to the classical equation of motion in the limit δ → 0 . The pair-wise force function f in Eq. (1) contains all the material-specific information so that the description of f would influence the calculation accuracy to a great extent. A micro-elastic material [1,3] is defined when the pair-wise force derives from a micro-elastic potential w(η , ξ ) =

f (η, ξ ) =

c(ξ , δ )s 2ξ , and then the force function can be expressed as 2

∂w(η, ξ ) ξ + η sc(ξ , δ ) ∀η, ξ = ∂η ξ +η

(3)

where ξ = x′ − x and η = u′ − u denote the relative position and relative displacement vector in the reference configuration, respectively (as shown in Fig. 1), and c (ξ , δ ) = c (0, δ ) g (ξ , δ ) is the micro-elastic modulus function denoting the stiffness of a pair-wise bond. ξ = ξ , and s stands for the stretch of a bond, s=

η +ξ − ξ

(4)

ξ

Damage of a peridynamic material at point x can be defined locally in term of the ratio of amount of broken bonds to total amount of interactions [1,3] as,

D ( x, t ) = 1 −

∫

Hx

μ ( x, ξ , t )dVx′

∫

Hx

(5)

dVx′

where μ ( x, ξ, t ) is a factor mapping the breakage of bond, ⎧1 ⎩0

μ ( x, ξ, t ) = ⎨

s (t ′, ξ ) < s0 , 0 < t ′ < t else

(6)

and s0 denotes the critical stretch or failure of the bond, which can be correlated with the classical fracture parameters of materials, such as fracture energy [3] or fracture strength 6

[29,30].

3. Extended peridynamic approach 3.1. Improved peridynamic constitutive model The mathematical framework of the peridynamic method leads to a key dependence of simulation accuracy on the formulation of the micro-modulus function c (ξ , δ ) which contains all material-specific information. Except that Kilic [8, 31] and Ha [14] have tried to use different micro-modulus functions, in most of the previous studies, the function c (ξ , δ ) is reduced to a constant c(δ ) by simplifying the kernel function g (ξ , δ ) to

⎧1 g (ξ , δ ) = ⎨ ⎩0

ξ ≤δ ξ >δ

(7)

In fact, the kernel function g (ξ , δ ) describes the spatial distribution of the intensity of long-range forces in materials, and it holds the following relations ⎧ g (ξ , δ ) = g ( −ξ , δ ) ⎪lim g (ξ , δ ) = max g ⎪ ξ →0 ⎪ ⎨lim g (ξ , δ ) = 0 ⎪ξ →δ ∞ ⎪ ∞ g (ξ , δ ) = ∫ Δ (ξ ) dx = 1 ⎪⎩ ∫−∞ lim −∞ δ →0

(8)

where Δ (ξ ) is the Dirac Delta function. In this study, Δ (ξ ) , rather than the common form of δ (ξ ) , is used to avoid confusion with the horizon δ . Eq. (8) reflects the non-local nature of the peridynamic theory exactly, and all the conditions given in Eq. (8) are satisfied when using the simplified form of g (ξ , δ ) shown in Eq. (7). However, it is clear that the effect of the distance between particles on the stiffness of the bond, i.e., the internal length effect is ignored so that the intensity of the long-range forces between particles keeps the same within the horizon. Thus, it does not match the physical nature of long-range forces much well.

7

In the present work, another form of the kernel function is uniquely proposed as follows.

ξ 2 2 ⎧ ⎪(1 − ( ) ) g (ξ , δ ) = ⎨ δ ⎪⎩0

ξ ≤δ

(9)

ξ >δ

It meets all the requirements described in Eq. (8) and reflects the weakening of the intensity of the long-range force when the distance between two particles increases. By following the same procedure performed to calculate the micro-modulus function through the equivalent of elastic strain energy in elasticity and PD theory [3,4], the expressions of the micro-modulus function in 2D plane stress conditions are obtained as 9E

c(ξ , δ ) =

(10)

πδ 3

when the kernel function is simplified as shown in Eq. (7). And the micro-modulus function is 315 E ξ (11) (1 − ( ) 2 ) 2 3 δ 8πδ when the proposed kernel function g (ξ , δ ) in the form of quartic polynomial as shown in c (ξ , δ ) =

Eq. (9) is employed. The expression given in Eq. (10) corresponds to the commonly-used constant of bond stiffness in previous studies [3, 5, 29], i.e., c=

6E πδ (1 −ν )

(12)

3

when the Poisson’s ratio is set as ν = 1/ 3 . Note that Eqs. (10)-(12) are exclusive for the 2D plane stress bond-based PD analysis. For 3D or other analysis, the expression of the micro-modulus function c (ξ , δ ) can be derived in a similar fashion. A comparison of the commonly used constant micro-modulus in the previous studies

8

with the form proposed in this study is shown in Fig. 2.

With the commonly used constant

bond stiffness, an identical force intensity for all bonds within the horizon is defined, and it comes down to zero suddenly when the length of a bond is beyond the horizon size. Comparatively, the bond stiffness in the present model is defined to be a continuous function of the distance between the particles in the bond.

The bond will gradually get weaker as the

distance increases, and its stiffnesseventually becomes zero when the length of bond reaches the horizon size.

It appears to fit the physical nature of long-range force well, and Kilic and

Madenci [8] and Kilic [31] discussed the similar internal length effects. In order to introduce damage and failure into the peridynamic model, Silling and Askari [3] and Ha and Bobaru [14] related the critical stretch s0 in Eq. (6) with the fracture energy

G0 which denotes the energy per unit fracture (crack) length in 2D (fracture area in 3D) required for complete separation of a body into two halves, or for breaking of all the bonds that initially connect material points in the opposite halves.

Following the same procedure

performed in [3,14] in the 2D plane stress case, the critical stretch s0 can be derived as

s0 =

960G0 (120π − 133)cδ 4

(13)

when the constant bond stiffness in Eq. (10) is used. And the critical stretch s0 is

s0 =

1024π G0 7(120π − 133) Eδ

(14)

when the improved micro-modulus function given in Eq. (11) is employed. It is notable that when considering c = 9E πδ 3 , the expression of s0 derived in Eq. (13) is exactly equivalent to the approximate value given in [14], i.e., s0 =

9

4π G0 . 9 Eδ

3.2. Local damping and loading rule for quasi-static analysis Engineering materials and structures exhibit some degree of damping especially under dynamic conditions.

The peridynamic equations of motion given in Eq. (1) could not be

employed directly to analyze static or quasi-static problems as material points in the model will keep moving all through.

Accordingly, up to date, quite fewer bond-based peridynamic

studies have focused on the quasi-static problems, such as quantitative elastic response of materials and structures; while much more investigations have been taken into dynamic fracture. In the present work, a local damping or dynamic relaxation is introduced into the peridynamic equations of motion for quantitative elastic deformation calculation and quasi-static fracture analysis.

It could be quite useful without bringing non-physical

phenomena since the dynamic effects play very small role on the behavior of materials in the quasi-static cases. The local damping is introduced into the PD equations of motion in term of the body forces as

ρ ( x)u( x, t ) + Cu = Lu ( x, t ) + b( x, t )

(15)

where C denotes the damping coefficient, and it is a positive real constant. In the previous bond-based peridynamic simulations, a key challenge exists in loading where the boundary particles will disperse from the rest of nodes with high acceleration when the external force is applied to the surface of the model. Therefore, the boundary bonds are usually set to be no-fail zones in order to prevent the tearing of boundary particles from the model.

In the earlier PD simulations, the external loads are usually given in term of

10

velocities or displacements of particles, which thus limits the application of peridynamic approach to some degree, e.g., in the quasi-static analysis of materials and structures. In this study, a step-by-step loading rule to avoid the tearing of PD model in a sudden rush is presented. For quantitative elastic analysis, the total external load F could be divided into m increments by introducing a loading factor λ

0 = λ0 < λ1 < λ2 < ⋅⋅⋅ < λm = 1

(16)

with an increment of Δλm = λm − λm −1 . In order to avoid the breaking of bonds or dispersion of particles resulting from a sudden rush in quasi-static simulation, the load increment should hold the condition 2 ρ s0 Δx . (17) Δt 2 where s0 is the critical stretch of bond, Δx is the grid spacing of the PD model, and Δt Δλm F <

is the time step size. In case of continuously simulating the quasi-static elastic behavior and fracture of materials and structures, the external load could be applied one increment by one increment since the extreme load or the ultimate load that a structure can hold is unknown in advance. In order to predict whether a PD model can sustain the equilibrium or crack would propagate unstably, a non-equilibrium criterion is proposed as follows. ψ ( u) ≤ β F

(18)

where ψ ( u) denotes the non-equilibrium force which exists in the particle system, F is the external load, and β is a given convergence standard which describes the allowable ratio of the non-equilibrium force to external load. Eq. (18) can also be written in detail as

11

N

∑ (ψ

n 2 j

)

j =1

≤β

N

∑ (F

(19)

n 2 j

)

j =1

where N is the total number of degrees of freedom, ψ nj = Lnuj + b nj

describes the

non-equilibrium force along the direction of the jth degree of freedom at time tn . In the present study, the non-equilibrium criterion is presented to judge the equilibrium of the particle system and decide whether another load increment can be applied on the system, so as to predict the critical failure load of a structure.

For example, after an external

load increment is applied on some surface particles at time t, the particles near the loading point will move and accordingly affect their surrounding particles; thus, the whole system will need some time steps to reach a new balance so that one more external load increment can be applied on the system.

However, when the non-equilibrium criterion cannot be

satisfied, it means that the system has not reached a new equilibrium state. When Eq. (18) or (19) is satisfied, it indicates that the structure can withstand the existing external loads. Otherwise, crack will initiate or propagate to reach a new balance or propagate unstably. In addition, a lower value for β will lead to a higher calculation accuracy but a higher computation cost as well. 3.3. Numerical discretization and iteration As usual, a numerical approximation is implemented by discretizing the peridynamic model into particles and reducing the spatial integration into a finite sum in order to solve the equations of motion. In this study, the peridynamic model is discretized into particles with the uniform grid spacing Δx . For dynamic analysis, the integral in Eq. (1) is replaced with a finite sum as

12

p

ρ uin = ∑ f (unj − uin , x j − xi )Vij + bin

(20)

j =1

where n denotes the number of time step, p is the number of particles within the horizon of point xi , and uin = u( xi , t n ) .

Vij stands for the involved calculation volume of particle

x j within the horizon of point xi , and it is usually given as Vij = VΔx = ( Δx )3 for 3D

simulations and Vij = VΔx = ( Δx )2 in 2D cases when the model is discretized with uniform lattices.

This simplification ignores the effect of current configuration on the effective

calculation volume of neighboring particles, thus correspondingly leading to the quantitative inaccuracy.

While in the present study, the boundary particles within the horizon are taken

into account, and the effective calculation field Vij is defined as ⎧ δ − ξ +η 1 + )VΔx (δ − rj ) ≤ ξ + η ≤ δ ⎪( 2rj 2 ⎪ ⎪ ξ + η ≤ (δ − rj ) Vij = ⎨VΔx ⎪ else ⎪0 ⎪ ⎩

where rj =

Δx

2

.

(21)

This modification is illustrated in Fig.3.

With the given initial and boundary conditions, the displacement of particle xi at time

tn +1 can be obtained by approximating the acceleration in Eq. (20) with an explicit central difference formula

in = u

uin+1 − 2uin + uin-1 Δt 2

uin+1 =

(22)

Δt 2 ⎡ n n n⎤ n n-1 ⎢ ∑ f ( u p − ui , x p − xi )V p + Fi ⎥ + 2ui − ui ρ ⎣ p ⎦

(23)

where Δt is the time step size. The requirement for a stable time step was derived by Silling and Askari [3], and this rule is adopted in the present work.

13

For quasi-static analysis with local damping and step-by-step loading, the displacement and velocity of particle xi at time tn +1 can be obtained similarly with the explicit central difference u in +1 = (( ρ / Δt − C / 2) u in + (λm F − Lnui ))( ρ / Δt + C / 2) uin +1 =

λm F + Lnui 2ρ ρ / Δt 2 − C / 2Δt n −1 n − + u u i Δt 2 ( ρ / Δt 2 + C / 2 Δt ) ρ / Δt 2 + C / 2 Δt i ρ / Δt 2 + C / 2 Δt

(24) (25)

The constitutive modeling and numerical approaches is implemented in Fortran using an in-house 2D serial peridynamic code.

4. Numerical examples In order to verify the unique capabilities of the present extended peridynamic approach for solving both quasi-static and dynamic fracture problems as well as performing quantitative elastic analysis accurately, several typical numerical examples including static elastic deformation, quasi-static mode I and mode II crack propagation, and dynamic failure are conducted in this section, and the results are compared with the numerical finite element analysis, available analytical prediction, or experimental observations. To simulate the mechanical behavior of real linearly elastic brittle materials like concrete, in the following examples, the material properties considered are: the Young’s modulus of material E = 22 GPa , density ρ = 2,400 kg/m3, fracture energy G0 = 110 J/m2, and uni-axial tensile strength ft = 1.6 MPa. For the sake of numerical calculation, the models are discretized into nodes with the uniform grid spacing Δx = 1 mm, and the time step size −8

Δt = 10 s according to the stable requirement.

When using the numerical approach

presented in this work, the non-equilibrium factor in Eq. (19) is set to be β = 5 ×10-7 , and the local damping coefficient for quasi-static analysis is given as C = 106 kg / m 3 ⋅ s based on

14

the convergence, efficiency and accuracy of tests. 4.1. Quantitative deformation and fracture analysis In the previous studies, the bond-based peridynamic method has seldom been utilized to calculate the elastic response of materials and structures quantitatively for multiple reasons. In this section, two kinds of typical static elastic calculation, i.e., a brittle slab under uni-axial tensionor compression and a cantilever beam subjected to shearing force, are considered. A 2D peridynamic model for a quasi-brittle (e.g., concrete) slab with the length of 200 mm and the width of 100 mm is shown in Fig. 4. The left side of the slab is fixed, and the right side is free and subjected to uni-axial tension and uni-axial compression, respectively as shown in Fig. 4(a) and (b).

Similar to the conventional bond-based PD simulations

previously performed, in this example, the external load appears in term of a given fixed displacement of d = 0.1 mm. The comparison of displacement contours in the slab between the PD and FEM (executed by ANSYS software) results are shown in Figs. 5 and 6. It is evident that in both the tensile and compressive cases, the distribution of displacement contours in the slab resulted from PD simulation agrees quite well with that of the FEM analysis.

What’s more, there is nearly no

difference on the horizontal displacement between results from two methods (as shown in Figs. 5(a) and 6(a)), and the maximum error on the vertical displacement is no more than 1% (as shown in Figs. 5(b) and 6 (b)).

The close displacement comparisons demonstrate the

high accuracy of the present approach for quantitative elastic calculation when using the conventional loading method. In order to further verify the accuracy of the proposed approach on quantitative calculation of static elastic responses of materials and structures, the PD model for a typical 2D cantilever beam with the size of L = 1.0 m and H = 0.2 m is established as shown in 15

Fig. 7. The presented numerical algorithms, such as local damping, gradual loading and non-equilibrium criterion, are employed in this example. First, a given concentrated load F = 1 kN is applied on the free end of the beam (on the far right layer of particles), and the elastic displacements of the beam are calculated by using the FEM and proposed PD approach, respectively. The comparison of results by the PD and FEM analysis is shown in Fig. 8. As shown in Fig. 8, the displacement contours in the cantilever beam using the proposed PD model are consistent with those by the FEM analysis.

The maximum displacements

along the vertical and horizontal directions given by the PD simulation are 2.376 × 10 −5 m and 3.357 × 10 −6 m, respectively, which closely match the corresponding FEM results of 2.375 × 10 −5 m and 3.384 × 10 −6 m.

Considering the FEM results as an exact solution for the

realistic response of the structure in such a simple static elastic case, the proposed PD approach illustrates its capability of accurately modeling the static elastic problems both qualitatively and quantitatively. The available analytical solutions for 1D peridynamic elastic bar/rod problems [28,32] showed some nonlocal features that would not be present in the classical solution, such as a pattern of wave-like decaying oscillations that spread out to infinity from the loading region, and a repeating pattern of discontinuities from the loading region [28]. For the presented 2D elastic cantilever beam problem, the analytical solution for the vertical displacement of the mid-point at the free end of the beam, when considering the shearing deformation and using the continuum elasticity, can be given as u y ( L )= − F (

L3 L + ) 3 EI κ GA

(26)

16

where F is the applied external load, EI and GA are the bending and shear stiffness of the beam, respectively, and κ = 2 / 3 is the shear correction factor. Fig. 9 shows a comparison of the load-displacement curves by using the present approach with respect to the analytical solution by using the continuum elasticity, finite element solution, and the results using the commonly-used peridynamic kernel function. The horizontal axis denotes the external concentrated load applied on the free end of the beam, and the vertical axis stands for the vertical displacement of the mid-point at the free end. It appears that the results from the present model is consistent with the analytical solution and finite element results.

The present extended peridynamic approach, when

considering the internal length effect of particles, can quantitatively calculate the displacement of the beam with higher accuracy. To investigate the independence of the quantitative elastic analysis to the grid spacing, the other four different grid spacing sizes, i.e., Δx = 2 mm, 4 mm, 8 mm, and 10 mm are chosen to calculate the displacement of the cantilever beam when the external concentrated load keeps the same as F = 1 kN . Fig. 10 shows the solutions for the vertical displacements of points on the beam axis when subjected to the external load of F = 1 kN at the free end. It is clear that the quantitative displacement calculation results keep unchanged with various grid spacing sizes. Further, the identical PD model is implemented to simulate their reversible crack (damage) initiation and crack propagation until ultimate failure, i.e., the whole quasi-static response process of the cantilever beam (see Fig. 7) under the increasing external load. This problem seems to be a long-standing challenge in computational mechanics since although the classical methods, such as FEM, are capable of analyzing the elastic response of 17

structures effectively and satisfactorily, they need extra fracture criteria and numerical remedies to handle discontinuity-concerned problems, such as crack initiation and propagation. For this kind of simulation, the following numerical rules are employed: the concentrated force acting on the particles at the far right end of the cantilever beam increases gradually with a uniform load increment fulfilling the conditions described in Eq. (17). As a consequence of equations of motion, the surface particles will move with the updated acceleration and naturally affect their neighboring particles. Hence, the response of the cantilever beam to the external load increment transmits from the loading location (the free end) to the far-away fixed end, just like wave propagation in solids.

After a period of

iteration time similar to the relaxation in molecular dynamics simulation, the structure will reach an equilibrium state if the condition given in Eq. (18) is met, which means that a new load increment may be applied on the structure for failure simulation.

In this study, the load

increment is set to be 1 kN initially for the elastic deformation simulation, and once damage occurs in the structure, the load increment will reduces to 0.1 kN for the failure simulation. For the problem shown in Fig. 7, the extreme load which will lead to the cracking of the cantilever beam can be approximately estimated through the simplified elastic calculation as Fex =

ft H 2 6L

= 10.67 KN

(27)

and it is predictable that for a homogeneous isotropic model, the crack will initially take place on the upper surface at the fixed end under the downward tip load, where the maximum tensile stress is located. The damage map (crack paths) in the cantilever beam using the present PD approach is shown in Fig. 11. It appears that the beam holds a downward tip load of 10.8 kN and reach

18

the balance of elastic deformation without any damage or bond breaking (Fig. 11(a)). However, when another load increment of 0.1 kN is applied to the structure, damage appears on the upper surface at the clamped end, just as analytically predicted. Furthermore, once damage occurs locally with several bonds being broken at the fixed end, the structure cannot reach a new equilibrium state again.

In contrast, crack will emerge and propagate rapidly,

which leads to the failure of the structure, as shown in Fig. 11(b)-(d). It is notable that no damage occurs in the other region in the beam except the field near the crack. As a whole, the PD results including the location where the crack will initiate and the path along which the crack will propagate show remarkable agreement with the analytical prediction.

Also, the extreme load resulted from the PD simulation is about 10.8-10.9 kN,

which agrees with the analytical results as well. 4.2. Quasi-static crack propagation Crack propagation simulation is an active and persistent challenge in the classical continuum mechanics theory. In this section, the PD simulation on typical mode I and mode II cracks in a finite-sized plate is undertaken. The size of the plate is 400 mm x 400 mm, and a pre-existing notch with the length of 200 mm starts from the middle point of the left boundary horizontally into the plate, as shown in Fig. 12.

The quasi-static step-by-step loading style as mentioned in Section 3 is employed up to a point which does not set the crack propagation, and smaller load increments are then applied one by one, as done in [33], until the last load increment which will lead to unstable crack propagation.

The external loads for mode I and mode II cracksare given in term of

uniformly distributed tensile force and shearing force acting on the top and bottom layers of 19

particles (see Fig. 12), respectively. The damage contours which reflect the crack propagation paths in the plate are shown in Figs. 13 and 14.

It is validated that the present approach inherits the capabilities of

simulating crack growth using the identical model from the usual (conventional) PD model, but provides an enhanced capacity for the quasi-static fracture analysis. The damage initiation and crack propagation under far away mode I tensile load is shown in Fig. 13. With a quasi-static loading rule and the load increment of 0.1 MPa, several material bonds near the crack tip fails and the local damage occurs when the external tensile stress increases to 0.6 MPa (much lower than the tensile strength of the material, i.e., 1.6 MPa), as shown in Fig. 13(a).

However, the particle system can still reach equilibrium

despite of the broken bonds near the crack tip.

At this moment, the load increment is

reduced to 0.01 MPa. When the external load increases to 0.7 MPa (i.e., after ten more load increments of 0.01 MPa), it appears that the plate cannot keep equilibrium longer, and the crack will propagate along the horizontal direction from the crack tip to the right side of the plate rapidly, as shown in Fig. 13(b)-(d). According to the linear elastic fracture mechanics theory and the handbook of stress intensity factor [34], the stress intensity factor for the problem described in Fig. 12(a) is

KI = Wσ π a

(28)

where a denotes the length of the pre-existing crack (i.e., 200 mm in this example), and W is a geometrical factor.

For this example, W can be obtained according to the handbook of

stress intensity factor [34] as W = 1.12 − 0.23(

a a a a ) + 10.6( ) 2 − 21.7( ) 3 + 30.4( ) 4 = 2.84 2a 2a 2a 2a

20

(29)

Considering the relationship between the critical stress intensity factor and the fracture 2 energy, KIC = GIC E , the critical remote stress which will result in crack propagation in this

example can be derived analytically as σ cr =

EGIC

π aW

2

= 0.691 MPa

(30)

Comparing the analytical value of critical stress with the PD results shown in Fig. 11, the present approach illustrates its capability of accurately predicting the quasi-static crack propagation quantitatively. Similarly, the mode II crack propagation in the brittle plate (see Fig. 12(b)) is shown in Fig. 14. When the remote shearing stress reaches 1.8 MPa, the damage occurs at the crack tip with the broken bonds; but the model can still keep in balance. The crack will propagate along some specific orientation rapidly when the shearing stress increases to 1.88 MPa, as shown in Fig. 14(b)-(d).

Again, according to the linear elastic fracture mechanics theory, the pre-existing crack in Fig. 12(b) will propagate along the maximum circumferential stress orientation, which can be derived by solving the following equation KI sin ϕ + K II (3cos ϕ − 1) = 0

(31)

For the mode I crack in Fig. 12(a), K II = 0 , and the solution for the crack propagation 0

orientation is ϕ = 0 , which again validates the crack paths given in Fig. 13. While for the mode II crack in Fig. 12(b), K I = 0 , so the crack will propagates along the orientation which fulfills 3 cos ϕ − 1 = 0 , i.e., ϕ = ±70.5 . In Fig. 14(b), it shows that as a consequence of the 0

PD solution, the mode II crack propagates along the orientation of ϕ = −70.2 , which 0

21

matches surprisingly well with the analytical prediction based on Eq. (31). 4.3. Dynamic fracture and crack branching Crack branching phenomena are observed in various experiments on dynamic brittle fracture [35,36].

However, up to date there is still no analytical solution for crack branching

problem, and it is an open but challenging task to simulate crack branching with the classical numerical methods.

Atomistic models (e.g., the molecular dynamics approach) have been

employed to simulate the crack branching, as done in [37]; but the atomistic simulation size is much limited, and the crack branching angles are much larger than the experimental observations. Kilic [13,31] and Ha and Bobaru [14,15] verified that the PD approach is able to correctly model and simulate the dynamic crack branching in brittle or quasi-brittle materials.

In this study, as a demonstration, the extended PD approach can simulate not

only the elastic response and quasi-static fracture quantitatively and qualitatively as shown in the above examples, but also it has the capability of simulating the dynamic fracture and crack branching problems. A brittle plate with the same size of 400 mm x 400 mm as shown in Fig. 12(a) is modeled, except that the length of the pre-existing notch is reduced to 20 mm in order to track the crack path more clearly.

In this example, the proposed local damping and loading rule for

quasi-static problems are removed from the calculation for dynamic analysis. The remote tensile stress (the mode-I loading) is applied on the top and bottom layers of the particles with a loading rate of 1 GPa/s, since the time step size here is in 10-8 s scale and it is unlikely to conduct an explicit computation over the time length of minutes as in the experiments [35, 36]. The crack propagation and branching process is shown in Fig. 15.

Since the

pre-existing notch in this example is much shorter than that in the above quasi-static case, it is analytically predictable that the higher stress level is required for crack propagation. It is 22

shown that when the remote tensile stress increases to 1.08 MPa, damage occurs at the crack tip (see Fig. 15(a)) and crack will then propagate horizontally as the usual mode I crack behaves (see Fig. 15(b)). The primary crack will branch successively and symmetrically as shown in Fig. 15(c) and (d).

Some branches will arrest but the others will propagate until

they reach the boundary, just as observed in the experimental studies [35, 36]. In this study, the quantitative elastic response analysis and quasi-static fracture problems are the focus of the proposed approach, and the dynamic fracture as in this example is just to demonstrate that the present extended PD model does not lose the capabilities of simulating dynamic fracture inherent from the PD formulation.

On the other hand, the loading style,

material properties and geometry of the plate are different from the experiments. Therefore, the dynamic fracture and crack branching simulation and the comparisons with experimental results are only qualitative in nature. The quantitative and detailed investigations into crack branching and other dynamic fracture observations, such as secondary cracks, crack-path instability, wave transmission, and reflection, by the present extended PD approach will be taken in the follow-up study.

5. Conclusions An extended bond-based peridynamic (PD) approach is developed, and an in-depth evaluation on the capabilities of the proposed numerical model to both qualitatively and quantitatively capture the elastic response and quasi-static fracture problems of materials and structures is performed.

The internal length effects of materials and the distribution

intensity of long-range forces are incorporated into the bond micro-modulus definition by introducing a continuous function in term of quartic polynomials, which improves the accuracy of quantitative calculation remarkably.

In addition, a local damping is introduced

into the PD equations of motion, and a step-by-step loading rule and non-equilibrium criteria 23

are proposed for quasi-static analysis.

The corresponding numerical discretization and

iteration algorithms are provided, and the effective calculation field of particles considering the current configuration is discussed as well. Several typical examples are analyzed with the present extended PD approach, and the numerical results match remarkably well with the numerical finite element predictions, analytical solutions, and experimental observations.

For the static elastic deformation

calculation, the extended PD approach can output the displacement results consistent with the finite element analysis qualitatively and quantitatively under various loading cases. Moreover, the whole response of structures under increasing external load, from elastic deformation to damage occurrence and crack initiation until ultimate failure, can be simulated all through with the same model, and the extreme load that a structure can hold is obtained accurately. For the quasi-static fracture analysis, the simulations on the crack propagation under the typical mode I and mode II loading show that the present model is able to not only analyze various crack growth processes spontaneously, but also predict the critical failure load and the crack paths without deviation.

What’s more, the extended PD approach naturally inherits

the capability of handling dynamic fracture problems and captures the experimental observations (e.g., crack branching) quite well.

Further efforts on the quantitative and

detailed analysis of dynamic fracture and other experiment-observed phenomena will be conducted in the follow-up study.

24

Acknowledgements The support of the National Natural Science Foundation of China (Nos. 11132003, 51209079, 51179064, 11372099, and 51478265) is gratefully acknowledged.

References: [1]

Silling, S.A. 2000, “Reformulation of elasticity theory for discontinuities and long-range forces,” Journal of the Mechanics and Physics of Solids, 48, pp. 175-209.

[2]

Tupek, M.R. and Radovitzky, R., 2014. “An extended constitutive correspondence formulation of peridynamics based on nonlinear bond-strain measures,” Journal of the Mechanics and Physics of Solids, 64, pp. 82-92.

[3]

Silling, S.A. and Askari,E., 2005, “A meshfree method based on the peridynamic model of solid mechanics,” Computers and Structures, 83, pp. 1526-1535.

[4]

Bobaru, F., Yang, M., Alves, L.F., et al., 2009, “Convergence, adaptive refinement, and scaling in 1D peridynamics,” International Journal for Numerical Methods in Engineering, 77, pp. 852-877.

[5]

Breitenfeld, M.S., Geubelle, P.H., Weckner, O., et al., 2014, “Non-ordinary state-based peridynamic analysis of stationary crack problems,” Computer Methods in Applied Mechanics and Engineering, 272, pp. 233-350.

[6]

Silling, S.A., 2003, “Dynamic fracture modeling with a meshfree peridynamic code,” The Second MIT Conference on Computational Fluid and Solid Mechanics, MIT, MA.

[7]

Weckner, O. and Abeyaratne, R., 2005, “The effect of long-range forces on the dynamics of a bar,” Journal of the Mechanics and Physics of Solids, 53, pp. 705-728.

[8]

Kilic, B. and Madenci, E., 2009, “Structure stability and failure analysis using peridynamic theory,” International Journal of Non-linear Mechanics, 44(8), pp.845-854.

[9]

Dayal, K. and Bhattacharya, K., 2006, “Kinetics of phase transformations in the

25

peridynamic formulation of continuum mechanics,” Journal of the Mechanics and Physics of Solids, 54, pp. 1811-1842. [10] Oterkus, S., Madenci, E., and Agwai, A., 2014, “Fully coupled peridynamic thermomechanics,” Journal of the Mechanics and Physics of Solids, 64, pp. 1-23. [11] Oterkus, S., Madenci, E., and Agwai, A., 2014, “Peridynamic thermal diffusion,” Journal of Computational Physics, 265, pp. 71-96. [12] Bobaru, F. and Duangpanya, M., 2010, “The peridynamic formulation for transient heat conduction,” International Journal of Heat and Mass Transfer, 53, pp. 4047-4059. [13] Kilic, B. and Madenci, E., 2009, “Prediction of crack paths in a quenched glass plate by using peridynamic theory,” International Journal of Fracture, 156, pp. 165-177. [14] Ha, Y.D. and Bobaru, F., 2010, “Studies of dynamic crack propagation and crack branching with peridynamics,” International Journal of Fracture, 162, pp. 229-244. [15] Ha, Y.D. and Bobaru, F., 2011, “Characteristics of dynamic brittle fracture captured with peridynamics,” Engineering Fracture Mechanics, 78, pp. 1156-1168. [16] Bobaru, F. and Silling, S.A., 2004, “Peridynamic 3D problems of nanofiber networks and carbon nanotube-reinforced composites,” Materials and Design: Proceedings of Numiform, American Institute of Physics, pp. 1565-1570. [17] Lehoucq, R.B. and Silling, S.A., 2007, Statistical Coarse-Graining of Atomistics into Peridynamics, Sandia National Laboratory Report, SAND2007-6410. [18] Silling, S.A., Epton, M., Weckner, O., et al, 2007, “Peridynamic states and constitutive modeling,” Journal of Elasticity, 88, pp. 151-184. [19] Silling, S.A., 2009. Linearized Theory of Peridynamic States, Sandia National Laboratory Report, SAND2009-2458. [20] O’Grady, J. and Foster, J., 2014, “Peridynamic beams: a non-ordinary, state-based model,” International Journal of Solids and Structures, 51, pp. 3177-3183.

26

[21] Mitchell, J.A., 2011, ANonlocal, Ordinary, State-Based Plasticity Model for Peridynamics, Sandia National Laboratory Report, SAND2011-3. [22] Askari, E., Bobaru, F., Lehoucq, R.B., et al., 2008, “Peridynamics for multiscale materials modeling,” Journal of Physics: Conference Series 125, 2008-012078. [23] Macek, R.W. and Silling, S.A., 2007, “Peridynamics via finite element analysis,” Finite Elements in Analysis and Design, 43, pp. 1169-1178. [24] Oterkus, E., Madenci, E., Weckner, O., et al., 2012,“Combined finite element and peridynamic analyses for predicting failurein a stiffened composite curved panel with a central slot,” Composite Structures, 94, pp. 839-850. [25] Warren, T.L., Silling, S.A., Askari, A., et al., 2009, “A non-ordinary state-based peridynamic method to model solid materialdeformation and fracture,” International Journal of Solids and Structures, 46, pp. 1186-1195. [26] Foster, J.T., Silling, S.A., and Chen,W.W., 2010, “Viscoplasticity using peridynamics,” International Journal for Numerical Methods in Engineering, 81, pp. 1242-1258. [27] Tupek, M.R., Rimoli, J.J., and Radovitzky, R., 2013, “An approach for incorporating classical continuum damage models in state-based peridynamics,” Computer Methods in Applied Mechanics and Engineering, 263, pp. 20-26. [28] Silling, S.A., Zimmermann, M. and Abeyaratne, R., 2003, “Deformation of a peridynamic bar”, Journal of Elasticity, 73, pp.173-190. [29] Gerstle, W., Sau, N., and Aguilera, E., 2007, “Micropolar peridynamic constitutive model for concrete,” 19th International Conference on Structural Mechanics in Reactor Technology (SMiRT19), Toronto, Canada. [30] Gerstle, W., Sau, N., and Sakhavand, N., 2009, “On peridynamic computational simulation of concrete structures,” ACI Special Publication, 65, pp. 245-264. [31] Kilic, B., 2008, PeridynamicTheory for Progressive Failure Prediction in

27

Homogeneous and Heterogeneous Materials, Ph.D. thesis, Tucson: The University of Arizona, Tucson, AZ. [32] Mikata, Y., 2012, “Analytical solutions of peristatic and peridynamic problems for a 1D infinite rod”, International Journal of Solids and Structures, 49, pp.2887-2897. [33] Song, J.H., Wang, H., and Belyschko, T., 2008, “A comparative study on finite element methods for dynamic fracture,” Computational Mechanics, 42, pp. 239-250. [34] Sih, G.C., 1973, Handbook of Stress Intensity Factors, Lehigh University, PA. [35] Bowden, F.P., Brunton, J.H., Field, J.E. et al., 1967, “Controlled fracture of brittle solids and interruption of electrical current,” Nature, 216, pp. 38-42. [36] Ramulu, M. and Kobayashi, A.S., 1985, “Mechanics if crack curving and branching -- a dynamic fracture analysis,” International Journal of Fracture, 27, pp. 187-201. [37] Zhou, S.J., Lomdahl, P.S., Thomson, R., et al., 1996, “Dynamic crack processes via molecular dynamics,” Physical Review Letters, 76, pp. 2318-2321.

28

Figures:

Fig. 1. Illustration of variables in a peridynamic model

Fig. 2. Comparison of micro-modulus function (stiffness of the bond)

Fig. 3. Effective calculation field of material points

29

d

L

L-d

d

Fig. 4. The quasi-brittle slab under (a) uni-axial tension; and (b) uni-axial compression

(a)

(b)

Fig. 5. Displacement contours in the slab under uni-axial tension (unit: m): (a) Horizontal displacement, and (b) Vertical displacement (upper: PD; lower: FEM)

30

(a)

(b)

Fig. 6. Displacement contours in the slab under uni-axial compression (unit: m): (a) Horizontal displacement, and (b) Vertical displacement (upper: PD; lower: FEM)

Fig. 7. A cantilever beam under concentrated load

31

(a)

(b)

Fig. 8. Displacement contours in the cantilever beam (unit: m): (a) Horizontal displacement, and (b) Vertical displacement (upper: PD; lower: FEM)

Fig. 9. Load-displacement curves by various methods

32

Fig.10. Quantitative displacement calculation results with respect to the grid spacing sizes

(b)

(a)

(b)

(d)

Fig. 11. Damage maps (crack paths) in the cantilever beam: (a) F = 10.8 kN ; (b) 1 μs after F = 10.9 kN ; (c) 30 μs after F = 10.9 kN ; and (d) 50 μs after F = 10.9 kN

33

(a)

(b)

Fig. 12. PD models for typical crack propagation: (a) Mode I crack, and (b) Mode II crack

34

(aa)

(b))

355

((c)

( (d)

M ; Figg. 13. 1 Daamaage maaps (modee I crac c ck pro p opaggatiion patths)) in n thee brrittlle pplatte: ((a) σ = 0..6 MPa

0 MP Pa ; an ( 23 (b) 2 μs afte a er σ = 0.77 MPPa ; (cc) 30 μs affter σ = 0.7 nd (dd) 42 4 μs after a r σ = 0.7 MPPa

(a))

(bb)

(c))

(dd)

366

Fig. 14. Damage maps (mode II crack propagation paths) in the brittle plate: (a) τ = 1.8 MPa ; (b) 2 μs after τ = 1.88 MPa ; (c) 12 μs after τ = 1.88 MPa ; (d) 20 μs after τ = 1.88 MPa

(a)

(b)

(c)

(d)

Fig. 15. Damage maps and crack branching in dynamic fracture: (a) t = 1.08 μs; (b) t = 1.12 μs; (c) t = 1.24 μs; (d) t = 1.38 μs

37

Nomenclature:

b : the external force per unit reference volume

β : the convergence standard which gives the allowable non-equilibrium force. c (ξ , δ ) : the micromodlulus function C : the damping coefficient

δ : the horizon Δ (ξ ) : the Dirac Delta function Δ (ξ ) : the common form of δ (ξ ) which is used to avoid confusion with the horizon δ f : the pair-wise force function

g (ξ , δ ) : the kernel function

G0 : the fracture energy H : the neighborhood of point x

λ : a loading factor

μ ( x, ξ, t ) : a factor mapping the breakage of bond n : the number of time step

N : the total numbers of degrees of freedom p : the number of particles within the horizon of point xi

ρ : the mass density

s0 : the critical stretch or failure of the bond Δt : the time step size u : the displacement of a material point 38

v: Poisson’s ratio Vij : the involved calculation volume of particle

x j within the horizon of point xi

Δx : the grid spacing of the PD model

ψ nj = Lnuj + b nj : the non-equilibrium force along the direction of the jth degree of freedom at time tn

39

An extended peridynamic approach for deformation and fracture analysis Dan Huanga, Guangda Lua, Chongwen Wanga, Pizhong Qiaob,c, a

b

Department of Engineering Mechanics, Hohai University, Nanjing, 210098, P. R. China

State Key Laboratory of Ocean Engineering and School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai, 200240, P.R. China

c

Department of Civil and Environmental Engineering, Washington State University, Pullman, WA, 99164-2910, USA

Highlights: -

An extended bond-based peridynamics for quantitative elastic and fracture analysis

-

A continuous function proposed in peridynamic constitutive model

-

A local damping incorporated into the peridynamic equations of motion

-

Quantitative deformation and quasi-static fracture analysis for 2D structures

-

Dynamic fracture and branching of a brittle plate under high-rate loading



Corresponding author at: Department of Engineering Mechanics, Shanghai Jiao Tong University, Mulan Building A903, Shanghai, 200240, PR China. E-mail addresses: [email protected]; [email protected] (P. Qiao).