Traction-Associated Peridynamic Motion Equation and Its Verification in the Plane Stress and Fracture Problems

Abstract

How to prescribe traction on boundary surface is still an open question in peridynamics. This problem is investigated in this paper. Through introducing the induced body force defined by boundary traction, the Silling’s peridynamic motion equation is extended to a new formulation called the traction-associated peridynamic motion equation, which is verified to be compatible with the conservation laws of linear momentum and angular momentum. The energy conservation equation derived from the traction-associated peridynamic motion equation has the same form as that in the original peridynamics advanced by Silling. Therefore, the constitutive models of the original peridynamics can be directly applied to the traction-associated peridynamic motion equation. Some benchmark examples in the plane stress problems are calculated. The numerical solutions agree well with the classical elasticity solutions, and the volume correction and the surface correction are no longer needed in the numerical algorithm. These results show that the traction-associated peridynamic motion equation not only retains all advantages of the original peridynamics, but also can conveniently deal with the complex traction boundary conditions.

1. Introduction

Peridynamics (PD) is a nonlocal continuum theory of mechanics developed in the recent two decades [1,2,3,4,5]. Its core consists in that a weighted integral of relative displacement over spatial domain is used instead of the gradient of displacement (strain) in the governing equations of deformation. Therefore, peridynamics can be used to conveniently and effectively analyze deformation companied with evolution of discontinuities. Peridynamics was firstly proposed by Silling [1] and then further improved by Silling and his collaborators [2]. Since then, it has been applied to investigate various problems associated with wave, damage, fracture, and impact breakage [3,4,5,6,7,8,9]. However, the traction boundary condition is incompatible with peridynamics because the governing equation of original peridynamics appears in the form of an integro-differential equation and does not involve the gradient with respect to spatial position. Therefore, boundary conditions cannot be imposed as naturally in the original peridynamics as in Classical Continuum Mechanics (CCM). A question therefore arises: how to incorporate suitable traction boundary condition?

Recalling the ways of the traction boundary condition is imposed in peridynamics, we can roughly divide them into three types. The first method [10] is to convert the traction on the boundary surface into the body force in an inner boundary layer according to the static force equivalence. The body force is usually supposed to uniformly distribute in the boundary layer, and the thickness of the inner boundary layer is taken as the spacing between material points. The second type of method is still based on the original peridynamics motion equation, but some modifications or operations will be performed to achieve the purpose of directly imposing traction boundary conditions. For example, a method only applicable to the non-ordinary correspondence models in peridynamics is proposed in [11], a weak form of peridynamic governing equations is proposed in [12,13], and a way based on Taylor expansion strategy to impose traction boundary conditions in ordinary state-based peridynamics is proposed in [14,15]. As the third method, Huang [16] suggested modifying the peridynamic motion equation to introduce the surface traction.

The traction boundary condition based on the first method is very popular in theoretical analysis and numerical calculation. Using this traction boundary condition, the well-posedness of the linear peridynamics with a given nonlocal kernel has been proved in mathematics [17,18]. Many numerical computation and analysis for practical problems can be found in [4,5,8]. However, it is unnatural and unhandy to convert the traction on the boundary surface to body forces in the inner boundary layer. Specially, when the material particles are close to the boundary surface or interface, the constitutive parameters need to be corrected [4,10,19,20]. Consequently, the traction boundary condition based on the first method is impractical for the sophisticated loading cases and geometrical surface.

If the non-ordinary correspondence models are adopted, the traction boundary condition can be specified directly. However, the non-ordinary correspondence models are limited due to the zero energy mode [21,22], and there are still some problems to be solved. The weak form of peridynamic governing equations has some changes comparing with the original peridynamic governing equation. Although the method of imposing traction boundary conditions based on Taylor expansion strategy can be successfully used in bond-based peridynamics, it still needs to set fictitious nodes. As for the model established by the third method, it is new and promising. Zhou [23] simplified the boundary transfer functions and proposed a method of imposing traction boundary conditions suitable for two-dimensional problems for bond-based peridynamics. However, since the three scalar-typical boundary transfer functions involved in this idea are difficult to determine, it has not been widely applied yet. Overall, how to characterize the traction boundary condition in peridynamics is still an open question and thus requires further investigation.

Meshfree discretization has been widely used to solve PD problems [24,25,26]. However, the standard discretization scheme is still very expensive for calculating the large-scale problems such as 3-dimensional problems. A new meshfree scheme [27] is expected to solve these difficulties.

In PD theory, there is no complete non-local neighborhood of a material point near a material boundary, resulting in the so-called skin or surface effects. Generally, the surface effect needs to be corrected to obtain the correct physical results. Various surface effect correction methods were studied in [10]. In addition, there are also many attempts to eliminate the surface correction [28,29,30,31,32,33].

The outline of the paper is as follows. In Section 2, through introducing the induced body force defined by boundary traction, we propose the traction-associated peridynamic motion equation and show that it is compatible with the conservation laws of linear momentum and angular momentum. In Section 3, the energy conservation law is derived from the traction-associated peridynamic motion equation. Two kinds of the bond-based constitutive models are discussed. The concrete form of the induced body force is determined. The prototype microelastic (PMB) constitutive model with the local damage is briefly introduced. The numerical algorithm is discussed in Section 4. In Section 5, the traction-associated peridynamic motion equation is used to calculate three benchmark examples in the plane stress problems. A fracture problem is simulated to verify the effectiveness of traction-associated peridynamic motion equation in the failure analysis. Finally, we close this paper with summary and comment.

2. The Induced Body Force and Extension of Peridynamic Motion Equation

Let p(x′′′, t) is a traction exerted at the point x′′′ on the boundary surface ∂Ω_p ⊂ ∂Ω of a peridynamic media Ω. Due to the nonlocality of peridynamic media, p(x′′′, t) will permeate the interior of Ω and induce a body force b_i(x, t) acting at the material particles. We call b_i(x, t) the induced body force, which is represented as

$$b_i\left(x,t\right)={\displaystyle {\int }_{\partial {\Omega }_p}\left[G\left(x,x^{‴},t\right)\cdot p\left(x^{‴},t\right)\right]d{A}_{x^{‴}}}.$$

(1)

In Equation (1), G (x, x′′′, t) is a second-order tensor field called the transfer function of boundary traction and with the dimension of 1/m³, which reads

$$G\left(x,x^{‴},t\right)=P\left(x,x^{‴},t\right)\cdot S\left(x^{‴},t\right),$$

(2)

where

$$\left\{\begin{array}{l}P\left(x,x^{‴},t\right)=\alpha \left(x^{‴},x\right)\frac{\left(y^{‴}-y\right)\otimes \left(y^{‴}-y\right)}{\left(y^{‴}-y\right)\cdot \left(y^{‴}-y\right)}\\ S\left(x^{‴},t\right)={\left[{\displaystyle {\int }_{\Omega }P\left(x,x^{‴},t\right)d{V}_x}\right]}^{-1}\end{array}\right..$$

(3)

In Equation (3), y′′′ = x′′′ + u′′′(x′′′, t) is the position vector of x′′′ in the deformed configuration and u′′′ = u′′′(x′′′, t) is a displacement field on ∂Ω_p. α(x′′′, x) is a weight function determining the characteristic of G(x, x′′′, t) and will be discussed in the following. By Equations (2) and (3), it is easy to verify

$${\displaystyle {\int }_{\Omega }G\left(x,x^{‴},t\right)d{V}_x}=I,$$

(4)

where I is the second order unit tensor. Equation (4) is a sufficient and necessary condition to ensure the compatibility of peridynamic motion equation with total equilibrium of linear momentum and angular momentum. Next, let us discuss this argument.

By Equation (4), the integrals of Equation (1) and y(x, t) × Equation (1) over Ω lead to

$${\displaystyle {\int }_{\Omega }{b}_i\left(x,t\right)d{V}_x}={\displaystyle {\int }_{\mathsf{\partial }{\Omega }_p}p\left(x^{‴},t\right)d{A}_{x^{‴}}}$$

(5)

$${\displaystyle {\int }_{\Omega }y\left(x\right)\times b_i\left(x,t\right)d{V}_x}={\displaystyle {\int }_{\mathsf{\partial }{\Omega }_p}{y}^{‴}\left(x^{‴}\right)\times p\left(x^{‴},t\right)d{A}_{x^{‴}}}.$$

(6)

After the induced body force b_i(x, t) is introduced, the motion equation of peridynamic media subjected simultaneously to boundary traction and external body force can be written as

$$\rho \left(x\right)\ddot{u}\left(x,t\right)=b_i\left(x,t\right)+{\displaystyle {\int }_{H_x}\left\{\mathsf{Ṯ}\left[x,t\right]〈\xi 〉-\mathsf{Ṯ}\left[x^{\prime },t\right]〈-\xi 〉\right\}d{V}_{x^{\prime }}}+b\left(x,t\right),$$

(7)

where ξ = x′ − x, H_x is a spherical neighborhood of x with radius δ, T the force vector state field [2], b(x, t) the external body force, and ρ(x) the mass density. We call Equation (7) the traction-associated peridynamic motion equation. Clearly, Equation (7) is an extension of the Silling’s peridynamic motion equation. Let x′′′∈∂Ω_p and x∈Ω. If we take G(x, x′′′, t) = H(|x′′′−x|)I/V_b where H(|x′′′ − x|) is the dimensionless square wave function and V_b the volume of the boundary layer of ∂Ω_p with the thickness of Δ, then Equation (1) will degenerate into a common formula when ones deal with the boundary traction in peridynamics.

By Equations (5) and (6), the integrals of Equation (7) and y(x, t) × Equation (7) over Ω yield

$${\displaystyle {\int }_{\Omega }\rho \left(x\right)\ddot{u}\left(x,t\right)d{V}_x}={\displaystyle {\int }_{\mathsf{\partial }{\Omega }_p}p\left(x,t\right)d{A}_x}+{\displaystyle {\int }_{\Omega }b\left(x,t\right)d{V}_x}$$

(8)

$${\displaystyle {\int }_{\Omega }\rho \left(x\right)\stackrel{·}{\overline{y\left(x\right)\times v\left(x,t\right)}}d{V}_x}={\displaystyle {\int }_{\mathsf{\partial }{\Omega }_p}y\left(x\right)\times p\left(x,t\right)d{A}_x}+{\displaystyle {\int }_{\Omega }y\left(x\right)\times b\left(x,t\right)d{V}_x},$$

(9)

which describe total equilibrium of linear momentum and angular momentum, respectively. Therefore, Equation (7) is consistent with the conservation laws of linear momentum and angular momentum. In addition, it is easy to see that Equation (7) is form-invariant under the Galileo transformation.

3. Peridynamic Constitutive Model

3.1. Balance Equation of Energy

Let v = v(x) is the velocity field within material, and ε is the internal energy density. Only elastic deformation is concerned, in peridynamics, total energy conservation can be represented as

$$\frac{D}{Dt}{\displaystyle {\int }_{\Omega }\left(\frac{1}{2}\rho v^2+\rho \epsilon \right)dV}={\displaystyle {\int }_{\Omega }{b}_i\cdot vd{V}_x}+{\displaystyle {\int }_{\Omega }b\cdot vdV}.$$

(10)

Equation (10) can be further written as

$${\displaystyle {\int }_{\Omega }\rho a\cdot vdV}+{\displaystyle {\int }_{\Omega }\rho \dot{\epsilon }dV}={\displaystyle {\int }_{\Omega }{b}_i\cdot vd{V}_x}+{\displaystyle {\int }_{\Omega }b\cdot vdV},$$

(11)

where a is acceleration field. In terms of Equation (7), Equation (11) reduces to

$${\displaystyle {\int }_{\Omega }\rho \dot{\epsilon }dV}={\displaystyle {\int }_{\Omega }\left\{{\displaystyle {\int }_{H_x}\left\{\mathsf{Ṯ}\left[x^{\prime },t\right]〈-\xi 〉-\mathsf{Ṯ}\left[x,t\right]〈\xi 〉\right\}d{V}_{x^{\prime }}}\right\}\cdot vd{V}_x}.$$

(12)

Since H_x ⊂ Ω is a compact supported set of T[x′, t]<−ξ> and T[x, t]<ξ>, Equation (12) can be written as

$${\displaystyle {\int }_{\Omega }\rho \dot{\epsilon }dV}={\displaystyle {\int }_{\Omega }\left\{{\displaystyle {\int }_{\Omega }\left\{\mathsf{Ṯ}\left[x^{\prime },t\right]〈-\xi 〉-\mathsf{Ṯ}\left[x,t\right]〈\xi 〉\right\}d{V}_{x^{\prime }}}\right\}\cdot vd{V}_x}.$$

(13)

Exchanging x′ and x, and then using definition of the compact supported set, we have

$${\displaystyle {\int }_{\Omega }\rho \dot{\epsilon }dV}={\displaystyle {\int }_{\Omega }\left\{{\displaystyle {\int }_{H_x}\mathsf{Ṯ}\left[x,t\right]〈\xi 〉\left[v\left(x^{\prime }\right)-v\left(x\right)\right]d{V}_{x^{\prime }}}\right\}d{V}_x}.$$

(14)

By the localized hypothesis [34,35], the balance equation of energy is given as follows

$$\rho \dot{\epsilon }={\displaystyle {\int }_{H_x}\mathsf{Ṯ}\left[x,t\right]〈\xi 〉\left[v\left(x^{\prime }\right)-v\left(x\right)\right]d{V}_{x^{\prime }}}.$$

(15)

Clearly, Equation (15) has the same form as that in the original peridynamics advanced by Silling and his collaborators [1,2]. Equation (15) is a basis to determine the peridynamic constitutive models of hyperelastic material. Therefore, the hyperelastic constitutive models in the original peridynamics can be inherited without modification by the traction-associated peridynamics.

3.2. Bond-Based Constitutive Models

Bond-based (BB) constitutive models are simplified versions of state-based (SB) constitutive models. For brevity, only the BB constitutive models are considered, not concerned the SB constitutive models. The BB constitutive models have been systematically established by Silling [1], and their forms are not unique. The microelastic models [5] are used to describe the elastic deformation of isotropic materials. Two commonly used microelastic models are listed below.

• General microelastic models:

$$f\left(u^{\prime }-u,\xi \right)=\left\{\begin{array}{ll}c\left(\left|y^{\prime }-y\right|-\left|\xi \right|\right)\frac{y^{\prime }-y}{\left|y^{\prime }-y\right|}\hspace{1em}& \left|\xi \right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right..$$

(16)

The linear form of Equation (16) is

$$f\left(u^{\prime }-u,\xi \right)=\left\{\begin{array}{ll}c\frac{\xi \otimes \xi }{{\left|\xi \right|}^2}\cdot \left(u^{\prime }-u\right)\hspace{1em}& \left|\xi \right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right.,$$

(17)

where f(u′−u, ξ) is the force density vector [1] with the dimension of N/m⁶. The relation between f(u′−u, ξ) and the force vector state T in Equation (7) is represented as T[x, t] <ξ> = −T [x′, t] <−ξ> == f(u′−u, ξ)/2. The parameter c is called the spring constant or the bond-constant, which reads [5]

$$c=\left\{\begin{array}{ll}\frac{15E}{\pi {\delta }^5}\hspace{1em}& 3-\mathrm{dimension}\\ \frac{12E}{\pi h{\delta }^4}\hspace{1em}& 2-\mathrm{dimension} \mathrm{plane} \mathrm{stress}\\ \frac{64E}{5\pi h{\delta }^4}\hspace{1em}& 2-\mathrm{dimension} \mathrm{plane} \mathrm{strain}\\ \frac{3E}{h_1{\delta }^3}\hspace{1em}& 1-\mathrm{dimension}\end{array}\right.,$$

(18)

where E is the Young’s modulus, h the thickness of plate and h₁ the cross-sectional area of rod.

2.

The prototype microelastic (PM) model:

PM model is another special form of the microelastic models

$$f\left(u^{\prime }-u,\xi \right)=\left\{\begin{array}{ll}c\frac{\left|y^{\prime }-y\right|-\left|\xi \right|}{\left|\xi \right|}\frac{y^{\prime }-y}{\left|y^{\prime }-y\right|}\hspace{1em}& \left|\xi \right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right..$$

(19)

The linear form of Equation (19) is

$$f\left(u^{\prime }-u,\xi \right)=\left\{\begin{array}{ll}c\frac{\xi \otimes \xi }{{\left|\xi \right|}^3}\cdot \left(u^{\prime }-u\right)\hspace{1em}& \left|\xi \right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right..$$

(20)

The parameter c in Equations (19) and (20) is still the bond-constant and it takes the value below [4,5]

$$c=\left\{\begin{array}{ll}\frac{12E}{\pi {\delta }^4}\hspace{1em}& 3-\mathrm{dimension}\\ \frac{9E}{\pi h{\delta }^3}\hspace{1em}& 2-\mathrm{dimension} \mathrm{plane} \mathrm{stress}\\ \frac{48E}{5\pi h{\delta }^3}\hspace{1em}& 2-\mathrm{dimension} \mathrm{plane} \mathrm{strain}\\ \frac{2E}{h_1{\delta }^2}\hspace{1em}& 1-\mathrm{dimension}\end{array}\right..$$

(21)

It should be noted that in the BB models, the Poisson’s ratio of 3D and 2D plane strain problems are fixed at 1/4, while that of 2D plane stress problem are fixed at 1/3.

3.3. Transfer Function of Boundary Traction

According to Equations (2) and (3), the transfer function G(x, x′′′, t) of the boundary traction characterizes the change of intensity when the traction is transferred from the boundary surface to the interior of the body, which is determined by the weight function α(x′′′, x). In physics, α(x′′′, x) should attenuates to zero as it moves away from the boundary surface. Therefore, α(x′′′, x) can be defined a function with compact support, i.e.,

$$\alpha \left(x^{‴},x\right)=\left\{\begin{array}{ll}q\left(\left|x^{‴}-x\right|\right)\ne 0\hspace{1em}& \left|x^{‴}-x\right|\le \iota \\ 0& \mathrm{otherwise}\end{array}\right.,$$

(22)

where x′′′ ∈ ∂Ω_p while x ∈ Ω. ι is a scale parameter. For simplicity, we take ι = δ. As thus, Equation (22) means that the traction on the boundary surface is dispersed in the boundary layer with the thickness of δ.

Consider the quasi-static uniaxial tension of a rod subjected to tensile force p at two ends. Through the inverse method [36] and the undetermined coefficient method, we find that when the general microelastic constitutive Equations (16) and (17) are adopted, if α(x′′′, x) takes the form below

$$\alpha \left(x^{‴},x\right)=\left\{\begin{array}{ll}k\left({\delta }^2-{\left|x^{‴}-x\right|}^2\right)\hspace{1em}& \left|x^{‴}-x\right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right.$$

(23)

then the classical elasticity solution of the uniaxial tension can be acquired.

Similarly, when the PM constitutive Equations (19) and (20) are used, if α(x′′′, x) is written as

$$\alpha \left(x^{‴},x\right)=\left\{\begin{array}{ll}k\left(\delta -\left|x^{‴}-x\right|\right)\hspace{1em}& \left|x^{‴}-x\right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right.$$

(24)

the same result is also given. It should be noted that k in Equations (23) and (24) is any constant and is not zero. For simplicity, we can take k = 1.

Equations (23) and (24) can be extended to 2-dimensional and 3-dimensional form below

$$\alpha \left(x^{‴},x\right)=\left\{\begin{array}{ll}{\delta }^2{\left|x^{‴}-x\right|}^2\hspace{1em}& \left|x^{‴}-x\right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right.$$

(25)

$$\alpha \left(x^{‴},x\right)=\left\{\begin{array}{ll}\delta -\left|x^{‴}-x\right|\hspace{1em}& \left|x^{‴}-x\right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right..$$

(26)

Thus, the concrete form of the induced body force is determined due to Equation (25) or (26), In the following, they will be directly used to analyze benchmark examples of 2D plane stress problems.

3.4. Prototype Microelastic Brittle Damage Model

In PD theory, local damage at a point is defined as the weighted ratio of the number of eliminated interactions to the total number of initial interactions of the material point with its family members, that is [4,19]

$$\phi \left(x,t\right)=1-\frac{{\displaystyle {\int }_{H_x}\mu \left(\xi ,t\right)d{V}_{x^{\prime }}}}{{\displaystyle {\int }_{H_x}d{V}_{x^{\prime }}}}.$$

(27)

It should be noted that the local damage φ ranges from 0 to 1. When φ = 1, all the interactions initially associated with the point have been eliminated, while φ = 0 means that all interactions exist. The measurement of the local damage value is an indicator of the possible formation of cracks within a body. In Equation (27), μ is a history-dependent scalar-valued function, which reads

$$\mu \left(\xi ,t\right)=\left\{\begin{array}{ll}1\hspace{1em}& s\left(\xi ,t^{\prime }\right)<S_c for all 0\le t^{\prime }\le t\\ 0& \mathrm{otherwise}\end{array}\right.,$$

(28)

where S_c is the critical stretch of bond failure, while s is bond stretch defined by

$$s=\frac{\left|y^{\prime }-y\right|-\left|\xi \right|}{\left|\xi \right|}.$$

(29)

Therefore, a simple way to introduce failure into the constitutive model to allow bonds (springs) to break when they are stretched beyond a predefined limit. After bond failure, there is no tensile force sustainable in the bond, and once a bond fails, it is failed forever (there is no provision for “healing” of a failed bond). As thus, the PMB model to characterize brittle damage [4,5,19] can be written as

$$f\left(u^{\prime }-u,\xi ,t\right)=\left\{\begin{array}{ll}cs\mu \left(\xi ,t\right)\frac{y^{\prime }-y}{\left|y^{\prime }-y\right|}\hspace{1em}& \left|\xi \right|\le \delta \\ 0& \mathrm{otherwise}\end{array}\right..$$

(30)

4. Numerical Algorithm

4.1. Spatial Discretization

Meshfree spatial discrete method [19] is used to discretize continuum into a range of arbitrary shaped subdomains, in which, the collocation points (nodes) are placed. With one-point Gauss quadrature strategy [37], the spatial discrete form of Equation (7) can be written as

$$\begin{array}{l}\rho \left(x_i\right)\ddot{u}\left(x_i,t\right)\\ ={\displaystyle \sum _{\mathsf{\partial }{\Omega }_p}G\left(x_i,x_k,t\right)\cdot \overline{p}\left(x_k,t\right)A_{x_k}}+{\displaystyle \sum _{H_{x_i}}\left\{\mathsf{Ṯ}\left[x_i,t\right]〈x_j-x_i〉-\mathsf{Ṯ}\left[x_j,t\right]〈x_i-x_j〉\right\}{V}_{x_j}}+b\left(x_i,t\right).\end{array}$$

(31)

The same spatial discrete strategy can be applied to acquire the integral value of G(x, x′′′). It is worth noting that the collocation points must also be set on the traction boundary surface due to the introduction of the boundary integration.

4.2. Time Integration

The adaptive dynamic relaxation (ADR) method [38] is used to solve Equation (31). After introducing new fictitious inertia and damping terms, Equation (31) is rewritten as

$${\lambda }_i{\ddot{u}}_i^n+c^n{\lambda }_i{\dot{u}}_i^n=F_i^n,$$

(32)

where n is n-th time (iterative) step, λ_i the fictitious density of the i-th node, and cⁿ the damping coefficient. The vector F_iⁿ in Equation (32) is the summation of internal and external forces, which reads

$$\begin{array}{l}{F}_i^n\\ ={\displaystyle \sum _{\mathsf{\partial }{\Omega }_p}G\left(x_i,x_k,t\right)\cdot \overline{p}\left(x_k,t\right)A_{x_k}}+{\displaystyle \sum _{H_{x_i}}\left\{\mathsf{Ṯ}\left[x_i,t\right]〈x_j-x_i〉-\mathsf{Ṯ}\left[x_j,t\right]〈x_i-x_j〉\right\}{V}_{x_j}}+b\left(x_i,t\right).\end{array}$$

(33)

By the explicit central-difference integration, the iterative scheme of displacements and velocities can be written as

$$\left\{\begin{array}{l}{\dot{u}}_i^{n+\frac{1}{2}}=\frac{\left(2-c^n\Delta t\right){\dot{u}}_i^{n-\frac{1}{2}}+\frac{2\Delta t{F}_i^n}{{\lambda }_i}}{2+c^n\Delta t}\\ u_i^{n+1}=u_i^n+\Delta t{\dot{u}}_i^{n+\frac{1}{2}}\end{array}\right..$$

(34)

As Equation (34) has an unknown velocity field at t^−1/2, the iterative process can be not started. However, if we assume u_i⁰ ≠ 0 and v_i⁰ = 0, the iteration can be completed by using

$${\dot{u}}_i^{\frac{1}{2}}=\frac{\Delta t{F}_i^0}{2{\mathsf{\lambda }}_i}.$$

(35)

In this algorithm, since the fictitious density λ_i, the damping coefficient cⁿ and the time step Δt are not actual physical quantities, their values can be chosen so as to make the convergence of numerical solution as fast as possible. Therefore, we take Δt = 1, while the calculation of λ_i and cⁿ can refer to [4,38,39,40].

5. Some Plane Stress Benchmark Problems

The nonlinear PM constitutive model is adopted to analyze some benchmark examples of the plane stress problems. In calculation, we take the Young’s modulus E = 200 GPa and the Poisson’s ratio ν = 1/3.

5.1. Example 1: A Rectangular Plate with Two Opposite Edges Subjected to Tension

As shown in Figure 1, the upper and lower edges of an isotropic rectangular plate are subjected to uniform tension q = 200 MPa. The length of the plate L = 1 m, the width W = 0.5 m, and the thickness h = Δ.

The rectangular plate (excepting near the boundary) is uniformly discretized into a particle set with the equal spacing Δ in the plane, as shown in Figure 2.

In Figure 2, the black dots represent the collocation points (nodes). Due to the symmetry of the structure and load, the constrains to rigid-body displacement need be imposed at the two symmetrical axes of the plate. It is necessary to collocate nodes on the boundary. The volume of subdomain associated with the boundary node is Δ²h/2 or Δ²h/4.

When the horizon size is specified as δ = 3.015Δ, three different grid sizes Δ = L/50, L/100, and L/200 are used to show the influence of the grid density on the convergence and computational accuracy, as illustrated in Figure 3 and Figure 4.

It can be seen from Figure 3 that the numerical calculation converges at the 500th time step regardless of the grid size. Moreover, Figure 4 shows that the computational results are close to each other for different grid density, and the errors between them and classical solutions are all within 5%. Therefore, the numerical algorithm can be considered to satisfy the δ-convergence requirement [41,42,43,44].

When the horizon size is fixed as δ = mΔ = 0.06 m, three different combinations (Δ = L/50, m = 3; Δ = L/100, m = 6; Δ = L/200, m = 12) are used to study m-convergence. The results are illustrated in Figure 5 and Figure 6.

Figure 5 that the calculation converges very quickly with time step, no matter which combination mode is adopted. From Figure 6, we see that the numerical results of the three different combinations agree with each other, and the errors between them and classical solutions are all within 5%. Therefore, the numerical algorithm is of m-convergence [41,42,43,44].

Comparison between Equation (7) and the original PD is carried out through numerical calculation. Figure 3 and Figure 7 show that the displacements calculated by Equation (7) and the original PD converges very quickly with time step, and the convergence modes are similar. As can be seen from Figure 8, Figure 9 and Figure 10, there exist good matches between predictions of traction-associated peridynamic motion equation and analytical solutions as well as original PD predictions.

In order to balance the computational accuracy and efficiency, in the following, the horizon size is specified as δ = 3.015Δ and grid size Δ = L/100. The total time step is set to 3000. The material parameters of plates in all examples are the same as those of the plate in Figure 1, and the same discretization as Figure 2 is adopted.

5.2. Example 2: A Rectangular Plate Subjected to Bending

As shown in Figure 11, a rectangular plate with the same size as the plate in Figure 1 is subjected to an anti-symmetrically distributed loads with a maximal value q = 200 Mpa. This is a pure bending problem with the stress boundary condition.

The convergence analysis is shown in Figure 12. It can be seen that the displacements calculated by the original PD and Equation (7) converge with time step in the same way, and completely converge at the 1000th time step.

The numerical results are illustrated in Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17, From which, it can be seen that the displacement distribution predicted by Equation (7) agrees with that by original PD, and the two have better matching with analytical solutions.

5.3. Example 3: A Square Plate with A Circular Hole Subjected to Tension by Two Opposite Edges

As shown in Figure 18, a squared plate with central circular hole is subjected to uniform tension q = 200 MPa. The side length of the plate L = 0.5 m and the radius of the circular hole r = 0.05 m.

The convergence analysis is shown in Figure 19. It can be seen that the displacements calculated by the original PD and Equation (7) converge with time step in the same way, and completely converge at the 500th time step.

The results in Figure 20 show that the displacements given by FEA, PD, and Equation (7) are close to each other in distribution. The relative error between them is less than 6.7%.

5.4. Example 4: Failure of A Square Plate with A Circular Hole under Quasi-Static Loading

We continue to investigate the fracture of the plate with central circular hole under tension. As shown in Figure 21, in order to avoid nodes occurring at the propagation path of cracks, all nodes (red and black dots) have been anew collocated. The grid size is taken as Δ = L/100. The nodes on the boundary surface (red dots) are only involved in the integration of the boundary traction and correspond to a volume of 0. In calculation, the critical stretch S_c of bond failure takes 0.0058.

The propagation of crack is characterized by the value of the damage φ. When the load q arrives at 380 MPa, the plate breaks due to cracking. Figure 22 shows the damage φ calculated by the original PD. At the 500th time step, the crack initiates from two sides of the hole. As the time step increases to 700, the crack propagates and the damage φ reaches 0.479. As the time step continues to increase, the plate breaks at the 900th time step.

The damage φ calculated by Equation (7) is illustrated in Figure 23. The results show that the crack initiates at the 700th time step, and at the 900th time step, the damage φ of the crack propagation reaches 0.42. When the time step arrives at 1125, the plate fails due to cracking.

6. Conclusions

Through introducing the induced body force defined by boundary traction, the Silling’s peridynamic motion equation is extended to the traction-associated peridynamic motion equation. From investigation on this equation, the conclusions are summarized as follows.

• The traction-associated peridynamic motion equation is consistent with the conservation laws of linear and angular momentum, and it is form-invariant under the Galileo transformation.
• The constitutive models in the original peridynamics can be inherited without modification by the traction-associated peridynamics. The concrete form of the induced body force is determined by matching with the constitutive models.
• Numerical calculations for the typical plane stress problems are in good agreement with the classical elasticity solutions, and the volume correction and the surface correction are no longer needed in the numerical algorithm.