Abstract
The constitutive modeling and numerical implementation of a nonordinary state-based peridynamic (NOSB-PD) model corresponding to the classical elastic model are presented. Besides, the numerical instability problem of the NOSB-PD model is analyzed, and a penalty method involving the hourglass force is proposed to control the instabilities. Further, two benchmark problems, the static elastic deformation of a simple supported beam and the elastic wave propagation in a two-dimensional rod, are discussed with the present method. It proves that the penalty instability control method is effective in suppressing the displacement oscillations and improving the accuracy of calculated stress fields with a proper hourglass force coefficient, and the NOSB-PD approach with instability control can analyze the problems of structure deformation and elastic wave propagation well.
1. Introduction
The most common numerical implementation technique of peridynamics (PD), a novel nonlocal continuum mechanics theory [1–3], is a mesh-free particle method [4]. Due to its intrinsic nonlocality [5], advantages in analyzing discontinuous problems of various physical fields [6, 7], and convenience for multiscale modeling [8], peridynamics has received extensive attentions from many fields. Currently, there are two mainstream peridynamic theories. The state-based peridynamic (SB-PD) theory [2, 9–16] is an extension of the original bond-based peridynamic (BB-PD) theory [1, 3]. In SB-PD theory, the response of an individual particle depends collectively on its interactions with neighboring particles through the state variables. Furthermore, the SB-PD can be divided into the ordinary state-based PD (OSB-PD) [9–11] and the nonordinary state-based PD (NOSB-PD) [12–16]. Most of the traditional constitutive models in classical continuum mechanics (CCM) can be incorporated into the NOSB-PD model with an approximate nonlocal deformation gradient [12–15]. The numerical stability/instability problem of the corresponding formulation of NOSB-PD model is the focus of the present paper.
The numerical stability/instability is significant in developing a robust numerical method. The particle discretization with one-point quadrature and the inaccurate enforcement of boundary conditions are two main sources of the numerical instability. Firstly, the one-point quadrature of many methods (e.g., dynamic finite element method, element-free Galerkin method [17], and mesh-free particle method [18]) will lead to unstable hourglass deformation and zero-energy mode, typically manifesting as the rapid jagged oscillations of nodes. Therefore some instability control algorithms are proposed to mitigate these instabilities, such as the viscous damping hourglass control algorithm and the elastic stiffness hourglass control algorithm [19]. Secondly, the inaccurate enforcement of essential boundary conditions of many mesh-free methods will reduce their accuracy and stability [20], and consequently several methods are developed to accurately exert the essential boundary conditions [21].
Similarly, the corresponding formulation of NOSB-PD may meet instability problems which are probably due to the particle discretization and the inaccurate calculation of variables within boundary regions. Beyond the two sources of instability, the presence of zero-energy mode instabilities in the NOSB-PD theory itself may be the main source of instability. Tupek [15] adopted two numerical implementations, namely, the particle discretization with one-point quadrature and the finite element discretization with Gauss quadrature, and confirmed that the discretization method could only result in slight instabilities. It shows that the particle discretization is not the main source of instabilities in NOSB-PD. Bessa et al. [22] noticed the problems of inaccurate enforcement of essential boundary conditions and the instability in NOSB-PD, yet the solving strategies were not given in their study. Wu and Ren [23, 24] pointed out that the numerical instability is profound near boundary and is apparently caused by the inappropriate boundary condition enforcement. By using the technique of mixed local/nonlocal gradient approximations and a stabilized force state, the instabilities of NOSB-PD were suppressed. Littlewood [25] first noticed that the zero-energy mode instability in NOSB-PD is due to the nonlocal deformation gradient, which allows nonphysical material response. Breitenfeld [14] pointed out that the regions with steep gradients suffer severe oscillations and introduced three modified force states to reduce the oscillations. Tupek [15] also studied the zero-energy mode and the nonphysical deformations due to the inaccurate calculation of nonlocal deformation gradient when the horizon size is either too large or too small. Madenci et al. [26] pointed out that the source of such oscillations is the inadequate approximation in the deformation gradient and the force density vector.
Considering the last two sources of the numerical instabilities in NOSB-PD and the method to suppress these instabilities, this study first presents the constitutive modeling and the numerical implementation of nonordinary state-based peridynamics corresponding to the classical elastic model, and then the penalty instability control method with a decreasing influence function is employed and studied in detail. The remainder of the paper is organized as follows: Section 2 introduces the correspondence formulation of NOSB-PD. In Section 3, the penalty method and the numerical implementation are introduced. Section 4 presents two benchmark problems, namely, the elastic deformation of a simple supported beam and the elastic wave propagation of a 2D rod, to study the sources of instabilities and the effect of present instability control method. Concluding remarks are provided in Section 5.
2. Corresponding Formulation of NOSB-PD
2.1. Basic Theory
Peridynamics states that a configuration is composed of massive material points. As shown in Figure 1, let , be the coordinate vectors of material points in the reference configuration, , the coordinate vectors of material points in the current configuration, and , the displacement vectors. Since then, the initial deformation vector state and the current deformation vector state can be defined aswhere denotes the initial relative position vector between point and point , and the angle bracket, , indicates the vector the state field operates on.
A material point interacts with all other material points within its horizon. Therefore, the integral-differential governing equation of SB-PD is given aswhere is the mass density, is the acceleration, and is a prescribed body force density. is volume of material point . The horizon size of is related to a characteristic length-scale of materials and/or of the considered phenomenon, and is a set of points within the horizon of point . Denote as , which represents the total force density of point . is a force state, which denotes a general constitutive model between the force state and the deformation state in SB-PD.
2.2. Constitutive Modeling
There are three main steps in the corresponding formulation of NOSB-PD modeling: describing structure deformation with the nonlocal deformation gradient tensor and then converting it to the strain tensor or the strain rate tensor, conducting strain analysis and stress analysis to obtain the stress tensor according to the utilized constitutive relation, and converting the stress into force state.
For the first step, the nonlocal shape tensor and the nonlocal deformation gradient tensor of point are defined aswhere the symbol denotes the dyadic product of two vectors. is the total number of material points within horizon. is the influence function. It should be noted that the shape tensor and the deformation gradient tensor should be both positive and reversible [12].
Furthermore, approximating the nonlocal gradient tensor as the traditional gradient tensor in CCM, the strain and stress analysis can be conducted. Besides, the principle of material frame indifference requires that a constitutive law, which is used during peridynamic modeling, must be objective [12]. For elastic infinitesimal deformation problems, the Lagrange/Green elastic strain tensor is given by
Then the second Piola-Kirchhoff stress can be obtained using the following elastic constitutive equation:where is the elastic stiffness matrix. Further the first Piola-Kirchhoff stress is defined bywhere is the Cauchy stress. Finally, the force vector state can be obtained with the following relationship:
Then the explicit dynamic solver is used to solve the motion equations and update the information of material points (e.g., acceleration, velocity, position, and force density). Hereto, we finish the analysis process in a time-step.
3. Numerical Instability Control and Implementation
We adopt the method of adding penalty force into the bond force to mitigate the zero-energy mode (or hourglass mode) instability [25]. A penalty force is applied to associate a nonzero quantity of work with hourglass modes of deformation, resulting in the elimination of spurious solutions of PD motion equations. The penalty force density or hourglass force density is defined as follows:where is the hourglass constant coefficient, taking the value of 10−3~101. is the micromodulus of the prototype microelastic brittle model [27]. , , , are the same variables as previously mentioned. is the projection of hourglass vector onto . Besides the hourglass vector can be defined aswhere is the predicted position of point by direction application of the approximate deformation at current configuration.
In addition to adding the penalty force, we adopt a decreasing influence function in Gaussian form instead of the common used constant one, , as follows:Compared with the constant one, the decreasing influence function in Gaussian form can reduce the computational errors of the deformation gradient tensor and the force vector state of points, especially those within the boundary regions. Therefore, it improves the enforcement accuracy of the boundary conditions and further improves both the computational accuracy and stability.
The penalty force is added to (2) to control instabilities, and an artificial damping term is added to obtain the quasistatic solutions of (2). Hence using explicit time-stepping and one-point quadrature, the discretization form of governing equation can be stated aswhere is the artificial damping coefficient, which can be assigned as to solve dynamic problems. The other variables are the aforementioned variables at th time-step.
The uniform spatial discretization, with grid spacing of and particle volume of in 2D, is usually adopted. Besides, a rough volume approximation of particles near the edge of the horizon circle is used to improve the accuracy [28]. In addition, there must be the initial displacement and initial velocity in the explicit dynamic method. The Verlet-Velocity scheme is used for temporal difference, stated as
4. Numerical Examples
4.1. Elastic Deformation of a Simple Supported Beam
Under the plane stress condition, as shown in Figure 2, the 2D simple supported beam, with its length of 1 m, height of 0.1 m, and unit thickness, suffers a concentrated force of 1 kN at the midpoint of the upper surface. The elastic modulus is GPa, and Poisson’s ratio is . In PD simulations, the beam is discretized into particles with gird spacing of mm and horizon size of . The artificial damping coefficient is assigned as kg/(m3·s). The time-step is s for a duration of s. The same model is simulated by the finite element method (FEM) using ANSYS software as well.
The vertical displacement solutions under different simulation conditions, with or without the instability control, are illustrated in Figures 3(a)–3(f). Comparing Figures 3(a) and 3(b) with 3(d), it indicates that although the PD solution and FEM solution match well in the global distribution, there exist obvious jagged oscillations in NOSB-PD solution, especially the displacement field near the loading region; the instability first occurred at the loading points and then spreads very quickly to the surrounding regions; and the decreasing influence function can reduce the instabilities or oscillations to some degree.
(a) PD solution (, )
(b) PD solution (, )
(c) PD solution (, )
(d) FEM solution
(e) Deflection of neutral axis
(f) Displacement of upper and lower surfaces
Figures 3(d), 3(e), and 3(f) show the PD solutions, in which the decreasing influence function and the penalty force are used to control instabilities. The PD results in Figure 3(c) and the FEM result in Figure 3(d) are in good agreement, indicating that the penalty method can suppress the instability of NOSB-PD effectively. Figure 3(e) shows the PD results with different hourglass force constant coefficients ranging from 10−3 to 101, illustrating that within a wide range of the coefficient, the penalty method can obtain relatively accurate results. However, the influences of penalty force on the accuracy and stability of solutions cannot be neglected. As shown in Figure 3(e), with a small hourglass force constant (), the instabilities still exist, especially near the loading regions; with a little larger hourglass force constant (), although the instabilities are removed, the solutions have lager errors; with a too large hourglass force constant (), the PD simulation will fail quickly. Therefore, when adopting the penalty method to control instabilities, the hourglass force coefficient should be determined properly with a recommended range of 10−3~101.
Figure 3(f) exhibits the vertical displacement of the upper and lower surfaces, indicating that no matter what value the hourglass force coefficient takes, there exists displacement mutation at the junction of the loading region and the free region of the upper surface, which results from the inaccurate calculation of the deformation gradient and the force state within the boundary region.
4.2. Elastic Wave Propagation in a 2D Rod
The analytic solution of a thin elastic rod, which is fixed at the right end and suffers an uniform pressure at the left end, is given as follows [29]:where is the wave speed.
This problem can be analyzed by the present method under the plane stress condition. As shown in Figure 4, the length of the 2D rod is and the width is . The mass density is , the elastic modulus is , and Poisson’s ratio is . All parameters are dimensionless. In PD simulations, the rod is totally discretized into 25653 particles with and . Three layers of particles at the right end are fixed, and an uniform pressure is enforced on three layers of particles at the left end for a duration of .
Figure 5 illustrates the longitudinal displacement obtained from PD simulations with different hourglass force coefficients at the time of , which is compared with the analytical solution. It can be seen that even with small time-step the displacement obtained from NOSB-PD simulation without instability control is severely oscillating, especially that near the loading region. On the contrary, even with large time-step the displacement obtained from NOSB-PD with instability control agrees well with the analytical one. The penalty instability control method is proved to be effective again. Figure 5(b) shows that the hourglass force coefficient will influence the accuracy and stability of NOSB-PD solution. Specifically, with a large hourglass constant (), the oscillations are removed, but the errors become large; with a little smaller hourglass constant ( and ), the slight oscillations still exist, especially near the loading region (seen from the local magnification). Therefore with the accuracy improvement of the deformation gradient and the force state within boundary region, the range of hourglass force coefficient will be wider, and the method will be more robust.
(a)
(b)
Similarly, Figure 6 illustrates similar longitudinal stress results compared with the displacement results obtained from the same simulations at the time of . However, the stress has obvious oscillations, which emerge because the second-order derivatives (strain) are approximated by applying the first-order derivative approximation (deformation gradient) twice (see (4) and (5)). Hence the accuracy of stress is influenced by the type of influence function and the horizon size.
Figures 7(a) and 7(b) show the curves of longitudinal displacement and longitudinal stress at the midpoint along with time. It also indicates that the hourglass force coefficient can affect the amplitude and phase of the elastic wave. Specifically, with a little smaller hourglass force constant the amplitude has slight error and the phase is accurate; with a too large hourglass constant both the amplitude and the phase distort away from the analytical solutions; and with a proper hourglass force constant an accurate and stable result can be obtained from the NOSB-PD.
(a) Longitudinal displacement
(b) Longitudinal stress
5. Concluding Remarks
This paper presents and numerically implements the nonordinary state-based peridynamics corresponding to the classical elastic model and studies the penalty method to control numerical instability in detail. Analysis of two benchmark problems proves that the penalty instability control method is effective in suppressing the displacement oscillations and improving the accuracy of calculated stress fields. Although a range of hourglass force coefficients is recommended, it still needs further study on determination of proper hourglass force coefficient for general problems. Besides, the accuracy of the deformation gradient tensor and the force state vector near the boundary or loading region need to be further improved, not only by choosing a proper influence function and horizon size, but also by calibrating the deformation gradient theoretically.
Competing Interests
The authors of the paper declare that there is no conflict of interests regarding the publication of this paper. The authors do not have a direct financial relation with the commercial identity that might lead to a conflict of interests for any of the authors.
Acknowledgments
This work is financially supported by the National Natural Science Foundation of China (nos. 11372099 and 11672101), Jiangsu Provincial Natural Science Foundation (no. BK20151493), and National Science and Technology Support Program of China (no. 2015BAB07B10).