Abstract
Peridynamics is based on integro-differential equations and has a length scale parameter called horizon which gives peridynamics a non-local character. Currently, there are three main peridynamic formulations available in the literature including bond-based peridynamics, ordinary state-based peridynamics and non-ordinary state-based peridynamics. In this study, the optimum horizon size is determined for ordinary state-based peridynamics and non-ordinary state-based peridynamics formulations by using uniform and non-uniform discretisation under dynamic and static conditions. It is shown that the horizon sizes selected as optimum sizes for uniform discretisation can also be used for non-uniform discretisation without introducing significant error to the system. Moreover, a smaller horizon size can be selected for non-ordinary state-based formulation which can yield significant computational advantage. It is also shown that same horizon size can be used for both static and dynamic problems.
Similar content being viewed by others
1 Introduction
Solid mechanics is an important area of engineering dealing with deformations of materials and structures under external loading conditions. Continuum mechanics has been widely used for this purpose for the last two hundred years, and there are currently different continuum mechanics formulations available in the literature with different advantages and limitations. The most common continuum mechanics formulation was developed by Cauchy, and the equation of motion of the fundamental object of continuum mechanics, i.e. “material point”, is expressed in the form of a partial differential equation. Since the analytical solution of this equation is limited to particular geometries, boundary conditions and material systems, different numerical techniques including finite element method have been developed to obtain solution for numerous problems of interest. However, due to the spatial derivatives in the partial differential equations, the standard solution procedures become invalid if discontinuities such as cracks exist in the solution domain. In such cases, additional steps should be taken to get around the discontinuity problem.
As an alternative approach, Silling [1] developed a new continuum mechanics formulation and named it as “peridynamics”. There are several fundamental differences between peridynamics and Cauchy’s continuum mechanics (CCM) formulation. First of all, the equation of motion of a material point in peridynamics is in integral form and does not contain any spatial derivatives. Therefore, it does not suffer from the discontinuity problem mentioned earlier. Moreover, peridynamics [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] is a non-local continuum mechanics formulation where a material point can interact with other material points which are located at a finite distance with respect to each other rather than only material points in its nearest neighbourhood as assumed in Cauchy’s formulation. An extensive review on peridynamics is given in [21]. The range of interactions between material points is denoted as “horizon” which is a length scale parameter in peridynamics. Such a parameter does not exist in Cauchy’s formulation. “Horizon” is a fundamental concept in peridynamics, and the term “peri” in the name basically corresponds to “horizon” in Greek language.
Although horizon is a very important parameter, research on how to choose this parameter has been rather limited and mainly depends on suggestions made in the influential paper written by Silling and Askari [2]. They suggested to use a horizon size equivalent to three times of the grid spacing between material points based on the experiences of these researchers for their simulations. However, their conclusion was obtained based on original peridynamic formulation, named as bond-based peridynamics. Although bond-based peridynamics is an effective approach, it has certain limitations. To overcome these limitations, advanced peridynamic formulations, such as ordinary state-based peridynamics and non-ordinary state-based peridynamics, were developed. Moreover, current peridynamic implementations are no longer limited to uniform discretisation. Therefore, it is critical to determine optimum horizon sizes for ordinary state-based peridynamics and non-ordinary state-based peridynamics formulations, so that sufficiently accurate results can be obtained in reasonable amount of time. Hence, this study focuses on determination of optimum horizon size for state-based peridynamic formulations using uniform and non-uniform discretisation under static and dynamic conditions.
2 Peridynamic theory
2.1 Bond-based peridynamics
By taking into these issues regarding CCM, Silling (2000) [1] proposed a new continuum mechanics formulation and named as peridynamics (PD). Silling relaxed the condition of interaction of material points which are directly in contact with each other in Cauchy’s formulation. Instead, all material points inside the structure can interact with each other. In peridynamics, an influence domain is defined to limit the range of interactions which is called horizon, \(H_{{\mathbf {x}}}\) (Fig. 1). Moreover, the equation of motion does not contain any spatial derivatives. Therefore, these equations are always valid regardless of discontinuities. The equation of motion of peridynamics can be expressed as
where \({\mathbf{f}}\left( {{\mathbf{{u}'}}-{\mathbf{u}},{\mathbf{{x}'}}-{\mathbf{x}}} \right) \) represents the interaction (bond) force between material points \({\mathbf{x}}\) and \({\mathbf{{x}'}}\) and \({\mathbf{u}}\) is the displacement of material point \({\mathbf{x}}\).
The definition of the interaction force depends on the material behaviour. For linear elastic isotropic materials, it can be considered as a spring force. However, it is represented in a slightly different form as
where \({\mathbf{y}}={\mathbf{x}}+{\mathbf{u}}\) is the position of the material point \({\mathbf{x}}\) in the deformed configuration. In Eq. (2), c denotes bond constant and s is the stretch of the bond which can be defined as
The interaction force is calculated in the deformed configuration. In the original peridynamic formulation (bond-based peridynamics), it is assumed that the force between two material points is equal and opposite to each other (Fig. 2). Bond constant is the material parameter of peridynamics and can be expressed in terms of material constants of CCM. These relationships can be established by considering a common parameter used in both approaches. For instance, strain energy density of a material point can be used for this purpose. An imaginary structure can be considered, and this structure can be subjected to a simple loading condition. A material point inside this structure should be identified, and its strain energy density can be calculated by using both peridynamics and CCM. Since the value of strain energy density calculated from both approaches should be the same, a relationship between peridynamics parameters and material constants of CCM can be established. For a linear elastic isotropic material, this relationship can be expressed as
where \(\delta \) is the size of the horizon. Note that there is only one peridynamic parameter used in the bond-based peridynamics formulation as opposed to two independent material parameters of CCM which can be chosen as elastic modulus, E and Poisson’s ratio, \(\nu \). Because of this mismatch, one of the parameters of CCM is constant. In other words, it is not possible to freely define Poisson’s ratio value since peridynamic formulation automatically captures a constant Poisson’s ratio value which is 1/4 for 3-dimensional geometries. For materials which have Poisson’s ratio value different than this value, advanced version of peridynamic formulations should be used. Currently, the most common advanced peridynamic approaches are ordinary state-based peridynamics and non-ordinary state-based peridynamics.
2.2 Ordinary state-based peridynamics
As mentioned earlier, the original bond-based peridynamic formulation encounters a limitation on material constants since it was assumed that the peridynamic force between two material points is equal and opposite to each other. Moreover, the peridynamic force between two material points only depends on the motion of associated material points. In order to eliminate the limitation of bond-based peridynamic formulation, its assumptions should be relaxed. In other words, it can be assumed that the magnitude of the force between two materials does not have to be equal to each other (Fig. 3). In addition to this, this force can depend on not only the motion of associated material points, but also the motion of their family members. Based on these new assumptions, the equation of motion of a material point can be written as
or
where \({\underline{\mathbf{T}}} \left\langle \bullet \right\rangle \) represents the force state which is a new terminology in state-based peridynamics [4] and t is the peridynamic force that material point \({\mathbf{{x}'}}\) is exerting on \({\mathbf{x}}\).
State is basically an infinite dimensional array or matrix and stores information about a particular parameter of peridynamic bonds associated with a particular material point. Therefore, force state stores the peridynamic forces belonging to peridynamic bonds. When a state operates on a particular bond, it only returns the stored value for that particular bond. For a linear elastic isotropic material, the force state can be expressed as
where a, b and d are peridynamic parameters and \(\theta \left( {{\mathbf{x}},t} \right) \) is peridynamic dilatation term which can be defined as
2.3 Non-ordinary state-based peridynamics
As mentioned earlier, in ordinary state-based peridynamics formulation, although the peridynamic forces between two material points can have different magnitudes, their directions are assumed to be along the bond direction. This assumption can be disregarded by allowing the direction of peridynamic forces in arbitrary directions (Fig. 4). By doing this, it is essential to explicitly impose a condition on conservation of angular momentum since it is not automatically conserved as in bond-based and ordinary state based formulations. Therefore, the following relationship must hold:
Moreover, the peridynamic force between two material points can be expressed in terms of appropriate stress definitions of CCM. Such an approach can allow direct integration of material models in CCM into peridynamics. For instance, a peridynamic force can be related to the first Piola–Kirchhoff (Lagrangian) stress tensor, \({\mathbf{P}}\) as
where \({\mathbf{K}}\) is the shape tensor which is defined as
where \({\underline{\mathbf{X}}} \left\langle \bullet \right\rangle \) is the position state which contains the relative position of materials points associated with a particular bond and \({\underline{\textit{w}}} \left\langle \bullet \right\rangle \) state contains the influence function information which defines the strength of interactions.
In order to incorporate material models of CCM in peridynamics, it is essential to relate the stress components to associated strain components. The definition of deformation gradient will be necessary to calculate the strain tensor. The deformation gradient, F, can be approximated in peridynamic framework as
where \({\underline{\mathbf{Y}}} \left\langle \bullet \right\rangle \) is the deformation state containing the relative position of bonds in the deformed configuration. Non-ordinary state-based peridynamics is a useful and practical approach. However, it encounters zero-energy mode problem and several techniques were proposed to overcome this problem [5]. In this study, the approach proposed by Silling [6] is used to remove the zero-energy mode problem by adding a stabilisation term to non-ordinary state-based peridynamic force term given in Eq. (10) as
where G is a positive constant, c is the bond-based peridynamic bond constant given in Eq. (4), \(\delta \) is the horizon size, and \(w_{0}\) is defined as
\({\underline{\mathbf{z}}}\) in Eq. (13) is the non-uniform part of the deformation state which can be expressed as
3 Numerical implementation
3.1 Spatial integration
Solution of peridynamic equation of motion by using analytical techniques is usually not possible. Therefore, numerical techniques are utilised and meshless approach is widely used for this purpose. The solution domain is discretised into finite number of volumes and each volume can be represented by a point located at its centre (Fig. 5). The peridynamic equation of motion in integral form can then be expressed in finite summation form as
where \(N_{k} \) is the number of points inside the horizon of the material point k.
For the material points close to the horizon boundary, only part of the material point is inside the horizon as shown in Fig. 6. Therefore, a correction parameter, \(\upsilon _{kj} \), is introduced and defined as
where \(r=\Delta /2\) with \(\Delta \) being the discretisation size and \(\xi _{jk} =\left| {{\mathbf{x}}_{j} -{\mathbf{x}}_{k} } \right| \).
Then, the peridynamic equation of motion given in Eq. (16) can be modified as
In addition, material points close to the surfaces do not have complete horizon as opposed to the material points located inside the solution domain. This causes these material points to have less stiffness with respect to internal material points and introduces error to the solution. To reduce this numerical error, surface correction factors are introduced to compensate the deficiencies due to lack of interactions. Hence, the peridynamic equation motion given in Eq. (18) can be further modified as
where \({\bar{\mathbf{{t}}}}_{kj} \) and \({\bar{\mathbf{{t}}}}_{jk} \) are corrected peridynamic forces if one of the interacting material points is influenced by the surfaces.
3.2 Time integration
3.2.1 Dynamic problems
For dynamic problems, either explicit or implicit integration schemes can be utilised. Although explicit integration scheme is easy to implement and does not have large memory requirement, it is only conditionally stable. The time step size should be chosen smaller than a critical time step size value which is given for ordinary state-based peridynamics as [3]
The value of the critical time step size is usually very small. Therefore, it is more suitable for short duration analysis. For long duration analysis, implicit time integration is more preferable. Please note that implicit time integration may have significant memory requirement for large number of material points in the system. In this study, explicit time integration is utilised for time integration of dynamic problems. After solving the acceleration of a material point at time step n, the velocity, \({\dot{\mathbf{{u}}}}\) and displacement, \({\mathbf{u}}\), of the same material point at the next time step \(n+1\) can be calculated as
where \(\Delta t\) is the time step size.
3.2.2 Static problems
For static problems, adaptive dynamic relaxation (ADR) is a common technique which artificially enforces the system to reach the steady-state condition [3]. Another approach that can be used is to assign the inertia term given in Eq. (19) to 0 and perform implicit solution by solving a matrix system. Again, such system may require large memory and special solution techniques to obtain the solution. In this study, ADR is used based on explicit type integration technique which does not cause any memory problems although reaching the steady-state condition can take certain number of time steps. According to ADR approach, new fictitious inertia and damping terms are introduced, so that the peridynamic equation of motion can be written as
where \({\mathbf{D}}\) is the fictitious diagonal density matrix, c is the damping coefficient,
and
for M material points.
The ith component of the force vector, \({\mathbf{F}}\) can be written as
By utilising the central-difference explicit integration scheme, displacement and velocities at the next time step can be obtained as
The time step size, \(\Delta t\), can be chosen as 1. Since explicit time integration is utilised, there is a stability condition to be satisfied to obtain a stable solution. In ADR, this can be satisfied by using sufficiently large values for the components of the fictitious diagonal density matrix, \({\mathbf{D}}\).
3.3 Application of boundary conditions
Application of boundary conditions in peridynamic theory is different than classical continuum mechanics. Rather than applying loading as point forces or distributed load, loads are exerted to a volume since the peridynamic equation motion is in integral form. In this study, boundary volumes are chosen as fictitious regions, \(\mathfrak {R}_{f} \) outside of the main solution domain, Re, as shown in Fig. 7, which is an effective procedure suggested in [7].
3.3.1 Displacement constraints
In this study, the size of the fictitious region is specified as twice the horizon size, i.e. \(2\delta \). The prescribed boundary value of the displacements \(U^{*}\), \(V^{*}\) and \(W^{*}\) in the \(x-\), \(y-\) and \(z-\) directions can be applied by specifying the displacements of the material points in the fictitious region, \(u_{f} \), \(v_{f} \), \(w_{f} \) in terms of displacements of the material points in the actual domain, u, v, w as (see Fig. 7),
3.3.2 Traction boundary conditions
Traction boundary conditions can also be applied similar to displacement constraints by using a fictitious region as shown in Fig. 8. The displacements of the material points in the fictitious region depend on the dimensionality of the problem (2-D or 3-D) and the unit normal of the boundary on which the traction boundary condition, \(\sigma _{ij}^{*}\), with \(i,j=x,y,z\), is exerted.
For 2-dimensional problems and the traction boundary with a unit normal in the \(x-\)direction, they can be written as [7]
Similarly, for the traction boundary with a unit normal in the \(y-\)direction, they can be written as
4 Determination of horizon size
In this study, we considered a special condition that there is no existence of damage in the structure and non-local effects are insignificant. For such a condition, peridynamic solution should converge to the classical continuum mechanics solution as horizon size converges to 0 [3]. Therefore, in this case, classical continuum mechanics solution can serve as a reference solution for peridynamics. Comparing peridynamic results against analytical and finite element method solutions will be sufficient to make decisions on the suitable value of horizon size. Several important aspects are explored to determine the horizon size for both uniform and non-uniform discretisation.
4.1 Horizon size for uniform discretisation
In this task, we explored how to determine the optimum size of the horizon with respect to the grid spacing and the size of the solution domain for uniform discretisation. As mentioned earlier, peridynamic equations are in the form of integro-differential equations and analytical solutions for such equations are limited. Therefore, numerical solution based on meshless discretisation is a common practice. For simplicity, uniform discretisation is utilised by dividing the solution domain into equal volumes with finite size and each volume is represented with a point located at its centre. Therefore, in practice, it is not possible to obtain a condition where horizon size becomes infinitely small. In such cases, Silling and Askari (2005) [2] suggested to use a horizon size equivalent to three times of the discretisation size, i.e. smallest distance between two neighbouring points based on their experience and observations of the results that they obtained using bond-based peridynamics (Fig. 9).
To make decisions on the horizon size, there are several important aspects to be considered. The first one is to take into account sufficient number of interactions between material points to capture all possible deformation modes. For instance, if a material point is only interacting with its nearest neighbours except the ones located at its diagonals, it is not possible to capture shear deformation especially in bond-based peridynamic formulation. The second important aspect is the dependence of horizon size on the discretisation size. Using smaller discretisation size will increase the accuracy of the numerical calculations. However, this will also increase the computational time. Therefore, it is important to determine an optimum discretisation size providing sufficient accuracy and leading to reasonable computational time. Since horizon size should be bigger than discretisation size, achieving efficiency in computational time introduces additional restriction on the horizon size approaching to its ideal size, i.e. becoming infinitely small. Currently, a common approach is to perform m- and \(\delta \)- convergence analysis by changing the discretisation size and horizon size, respectively. However, since such convergence studies are time-consuming especially for large-scale problems, it is essential to determine a horizon size which can be safely used in all applications and leading to minimum computational time. Another important aspect is the nature of the problem being static or dynamic. It is necessary to determine if the same size of horizon is suitable for both conditions. Moreover, dynamic problems may experience unphysical wave reflection which can lead to inaccurate results as the waves are travelling inside the solution domain. Finally, it is important to investigate if the same horizon size is applicable for all 1-dimensional (1-D), 2-dimensional (2-D) and 3-dimensional (3-D) geometries. Moreover, both ordinary state-based and non-ordinary state-based formulations are utilised and it is important to investigate if same horizon size can be used in both formulations. This investigation can also show if state-based formulations will have a different tendency with respect to bond-based formulations since in state-based peridynamic formulations the effect of family members are taken into account which is ignored in bond-based peridynamic formulation.
To determine the horizon size, first, simple geometries, boundary and loading conditions are considered under both static and dynamic conditions. To demonstrate the general applicability of the determined horizon sizes, more complicated problem cases are considered. Finite element solutions for the same cases are also evaluated as a reference solution.
4.2 Horizon size for non-uniform discretisation
In this task, we explored how to determine the optimum size of the horizon that we can use for non-uniform discretisation. Uniform discretisation is simple and widely used in peridynamic simulations. However, for certain problems, uniform discretisation becomes unfeasible such as a sub-region with a very small thickness. Moreover, as in the finite element simulations, it may be computationally advantageous to use different grid sizes at different parts of the domain (Fig. 10). Note that such a case is more prone to numerical errors.
5 Results and discussion
To determine the optimum horizon size, several benchmark problems are considered. To investigate the dynamic behaviour, a square plate is subjected to uniaxial strain condition as initial condition. For static analysis, same geometry is subjected to uniaxial tension condition. Results are obtained by using ordinary state-based peridynamics and non-ordinary state-based peridynamics with either uniform or non-uniform discretisation. Finally, more complex problems are considered to check if the determined optimum horizon size values can produce accurate results for these more complex problems.
5.1 Vibration of a plate with uniform discretisation
In the first problem case, a square plate with dimensions of \(L = W =\) 1m is considered. The plate has a thickness of 0.01 m and subjected to an initial condition in the form of uniaxial strain as shown in Fig. 11. Elastic modulus and density of the plate are specified as 200 GPa and \(7850\hbox { kg/m}^{3}\), respectively. The left edge of the plate is fully fixed by using a fictitious region as demonstrated in Fig. 12. The solution is obtained using explicit time integration with a time step size of \(1\times 10^{{-7}}\) s. The plate is uniformly discretised with a grid size of \(\Delta x = 0.01\hbox { m}\). Five different horizon sizes are considered as \(\delta =1\Delta x\), \(2\Delta x\), \(3\Delta x\), \(4\Delta x \) and \(5\Delta \textit{x}\).
5.1.1 Ordinary state-based peridynamics
Ordinary state-based peridynamics is a more general case of bond-based peridynamics which allows specification of any Poisson’s ratio. In this case, Poisson’s ratio is specified as 0.25. Comparison of the horizontal and vertical displacement values of the point located at (0.255m, 0.255 m) is shown in Figs. 13 and 14. Although horizon displacements of all horizon size cases agree well with finite element method (FEM) results obtained by using ANSYS, a commercially available finite element software, best agreement is obtained with the horizon size values of \(\delta =3\Delta x\) and \(\delta =4\Delta x \) for the vertical displacements. Horizon size values of \(\delta =1\Delta x\) and \(\delta =2\Delta x \) could not capture accurate displacements for the vertical displacements.
5.1.2 Non-ordinary state-based peridynamics
Vibration of a plate analysis is also considered by using non-ordinary state-based peridynamics and Poisson’s ratio of 0.25. As shown in Figs. 15 and 16, different than ordinary state-based cases, better agreement is obtained for the horizon size values of \(\delta =1\Delta x\) and \(\delta =2\Delta \textit{x}\). Other horizon size values \(\delta =3\Delta x\), \(4\Delta x \) and \(5\Delta x\) also yield accurate horizontal and vertical displacement values.
5.2 Plate under tension with uniform discretisation
In the second problem case, the square plate considered in the previous section is subjected to a uniaxial tension loading of \(\sigma \)* \(= 200\hbox { MPa}\) at the right edge as shown in Fig. 17. The applied loading is specified by creating a fictitious region at the right edge as demonstrated in Fig. 18. The steady-state solution is obtained by using the adaptive dynamic relaxation technique.
5.2.1 Ordinary state-based peridynamics
First, ordinary state-based formulation is utilised and horizontal and vertical displacements along the central axes are obtained as shown in Figs. 19 and 20. According to these results, it can be concluded that horizon size values of \(\delta =3\Delta x\) and \(4\Delta x \) provide better agreement with FEM results. Similar to the dynamic case, horizon values of \(\delta =1\Delta x\) and \(2\Delta x\) could not provide accurate displacement values.
5.2.2 Non-ordinary state-based peridynamics
The plate under tension problem is also investigated using non-ordinary state-based peridynamics formulation. As in the vibration case, all horizon size values yield accurate results as shown in Figs. 21 and 22. The best match with FEM results is achieved with the horizon size values of \(\delta =1\Delta x\) and \(2\Delta \textit{x}\).
5.3 Vibration of a plate with non-uniform discretisation
In Sects. 5.1–5.2, two different cases are considered to determine the suitable horizon size for 2-dimensional structures using ordinary state-based peridynamics and non-ordinary state-based peridynamics under static or dynamic conditions. Although uniform discretisation is widely used in the peridynamic studies available in the scientific literature, it can be computationally advantageous to have flexibility to use different grid sizes at different parts of the solution domain which is a common procedure used in numerical calculations such as in FEM. Hence, in this part of the study, we present the capability of peridynamics using variable grid sizes. As the first demonstration case, vibration of a plate problem considered is investigated. As shown in Fig. 23, the solution domain is split into two regions with a possibility that each of these regions can have different grid sizes.
By keeping the grid size in Region 2 constant, i.e. \(\Delta x_{2} = 0.01\hbox { m}\), three different grid sizes for Region 1 are considered, i.e. \(\Delta x_{1} = 0.01\hbox { m}\) (\(k=1\)), \(\Delta x_{1} = 0.005\hbox { m}\) (\(k=2\)) and \(\Delta x_{\mathrm {1}} = 0.025\hbox { m}\) (\(k=4\)) for the horizon size value of \(\delta _{2}=3\Delta x_{2}\).
5.3.1 Ordinary state-based peridynamics
The solution is first obtained by using ordinary state-based peridynamics. As shown in Figs. 24 and 25, although all three cases yield close results for the horizontal displacements, \(k=2\) and \(k=4\) cases yield better results since the number of material points inside the horizon is higher for these cases which reduces the error for numerical integration.
5.3.2 Non-ordinary state-based peridynamics
Next, non-ordinary state-based peridynamic formulation is utilised to obtain the solution. As shown in Fig. 23, by keeping the grid size in Region 2 constant, i.e. \(\Delta x_{2} = 0.01\), three different grid sizes for Region 1 are considered, i.e. \(\Delta x_{1} = 0.01 (k=1)\), \(\Delta x_{1} = 0.005 (k=2)\) and \(\Delta x_{1} = 0.025 (k=4)\) for the horizon size value of \(\delta _{2}=2\Delta x_{2}\). As shown in Figs. 26 and 27, all three cases yield close results regardless of discretisation size which shows that non-ordinary state-based peridynamics is not sensitive to discretisation being uniform or non-uniform.
5.4 Plate under tension with non-uniform discretisation
As the second demonstration case, plate under tension problem is considered to investigate the effect of non-uniform discretisation under static condition. As shown in Fig. 23, the solution domain is split into two regions with a possibility that each of these regions can have different grid sizes.
5.4.1 Ordinary state-based peridynamics
First, ordinary state-based peridynamics is utilised to obtain the solution. By keeping the grid size in Region 2 constant, i.e. \(\Delta x_{2} = 0.01\), three different grid sizes for Region 1 are considered, i.e. \(\Delta x_{1} = 0.01 (k=1)\), \(\Delta x_{1} = 0.005 (k=2)\) and \(\Delta x_{1} = 0.025 (k=4)\) for the horizon size value of \(\delta _{2}=3\Delta x_{2}\). As shown in Figs. 28 and 29, although all three cases yield close results for the horizontal displacements, \(k=2\) and \(k=4\) cases yield better results since the number of material points inside the horizon is higher for these cases which reduces the error for numerical integration.
5.4.2 Non-ordinary state-based peridynamics
Next, non-ordinary state-based peridynamics is used to obtain the solution. As shown in Fig. 23, by keeping the grid size in Region 2 constant, i.e. \(\Delta x_{2} = 0.01\), three different grid sizes for Region 1 are considered, i.e. \(\Delta x_{1} = 0.01\) (\(k=1\)), \(\Delta x_{1} = 0.005 (k=2)\) and \(\Delta x_{1} = 0.025 (k=4)\) for the horizon size value of \(\delta _{2}=2\Delta x_{2}\). As shown in Figs. 30 and 31, all three cases agree well with FEM results although there exist some oscillations at the interface region between Regions 1 and 2.
5.5 Complex cases
In addition to the several benchmark problems considered for the determination of the optimum horizon size, two more complex cases are selected to check the performance of the selected horizon sizes.
5.5.1 Plate with a square cut-out under tension
In the first problem case, plate with a square cut-out problem is considered as shown in Fig. 32. The plate is subjected to displacement loading at the right edge. The left and upper boundaries of the plate are constrained in x- and y- directions, respectively.
The plate with a square cut-out problem is first investigated by using ordinary state-based peridynamics with a Poisson’s ratio of 0.25 and horizon size of \(\delta = 3\Delta x\). Both FEM and peridynamic results are given in Figs. 33 and 34. According to these figures, it can be concluded that a very good agreement is observed between peridynamic and FEM analysis results.
Next, plate with a square cut-out problem is also analysed by using non-ordinary state-based peridynamics with a horizon size of \(\delta = 2\Delta x\). By comparing the reference FEM results given in Fig. 33, peridynamic results given in Fig. 35 have a very good agreement with FEM results.
5.5.2 Plate with a square cut-out vibration
In the second problem, plate with a square cut-out problem considered in the previous section is analysed under dynamic loading conditions by applying uniaxial strain condition to the entire domain as initial condition as shown in Fig. 36.
The problem is analysed both by using ordinary and non-ordinary state-based peridynamics and peridynamics results are compared with FEM results given in Fig. 37. According to the comparison, ordinary and non-ordinary peridynamic results shown in Figs. 38 and 39 match very well with FEM results under dynamic conditions for the selected horizon sizes of \(\delta = 3\Delta x\) and \(\delta = 2\Delta x\), respectively.
5.6 Discussions of results
Since numerous results are generated to determine the optimum horizon size for peridynamics and presented in the previous sections, several important observations are summarised and discussed in this section.
First, uniform discretisation is utilised to discretise the solution domain and two benchmark problems are considered including vibration of a plate and plate under tension. Hence, both dynamic and static conditions are explored. Solutions are obtained by using ordinary state-based and non-ordinary state-based peridynamics. Dynamic analysis and static analysis results yield similar conclusions about the horizon size selection. Based on numerical results, it was confirmed that horizon size of \(\delta = 3\Delta x \) is an optimum horizon size for ordinary state-based peridynamic analysis as for bond-based peridynamic analysis that Silling and Askari [2] suggested earlier. Although bond-based peridynamics is a special case of ordinary state-based peridynamics, the outcome is not intuitive since the peridynamic force in ordinary state-based peridynamics is dependent on motions of family members rather than only motions of interacting material points as in bond-based peridynamics formulation. For non-ordinary state-based peridynamics, it was concluded that \(\delta = 2\Delta x\) is the most optimum choice for horizon size. This will result in significant reduction in computational time especially for 3-dimensional analysis since there are much less material points inside a horizon of \(\delta = 2\Delta x \) than a horizon with a size of \(\delta = 3\Delta x\). This conclusion should also be confirmed with the existence of the damage in the solution domain. Based on the numerical results, it should also be noted that for non-ordinary state-based formulation the horizon should be small since the non-ordinary state-based model used in this study, i.e. correspondence model, has a direct connection with classical (local) continuum mechanics and requires determination of deformation gradient in peridynamic framework as the first step.
After determining optimum horizon sizes for all ordinary and non-ordinary state-based peridynamic formulations, the performance of these horizon sizes is investigated for non-uniform discretisation by considering the same benchmark problems for the uniform discretisation analysis. By considering different scenarios having refined discretisation on one side of the domain with same horizon size, it was concluded that if there are more material points within a horizon, it reduces the numerical error due to numerical integration and leads to a more accurate solution
Finally, the suggested horizon sizes are checked by considering more complex problems and comparison against reference FEM solution shows that suggested horizon sizes can provide accurate solutions with the advantage of computational efficiency.
6 Conclusions
In this study, optimum horizon sizes and shapes for two widely used peridynamic formulations, i.e. ordinary state-based peridynamics and non-ordinary state-based peridynamics, are investigated. After numerous numerical tests, it is concluded that horizon size of three grid size is suitable for ordinary state-based peridynamics. On the other hand, it is sufficient to use a horizon size of two grid size for non-ordinary state-based peridynamics. Moreover, it is shown that peridynamic implementation can be done using non-uniform discretisation. Although non-uniform discretisation can introduce numerical error especially for dynamic analysis, accurate results can be obtained if the suggested horizon sizes are used for such cases.
References
Silling, S.A.: Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids 48(1), 175–209 (2000)
Silling, S.A., Askari, E.: A meshfree method based on the peridynamic model of solid mechanics. Comput. Struct. 83(17–18), 1526–1535 (2005)
Madenci, E., Oterkus, E.: Peridynamic Theory and Its Applications. Springer, Berlin (2014)
Silling, S.A., Epton, M., Weckner, O., Xu, J., Askari, E.: Peridynamic states and constitutive modeling. J. Elasticity 88(2), 151–184 (2007)
Gu, X., Madenci, E., Zhang, Q.: Revisit of non-ordinary state-based peridynamics. Eng. Fract. Mech. 190, 31–52 (2018)
Silling, S.A.: Stability of peridynamic correspondence material models and their particle discretizations. Comput. Methods Appl. Mech. Eng. 322, 42–57 (2017)
Madenci, E., Oterkus, S.: Ordinary state-based peridynamics for plastic deformation according to von Mises yield criteria with isotropic hardening. J. Mech. Phys. Solids 86, 192–219 (2016)
De Meo, D., Zhu, N., Oterkus, E.: Peridynamic modeling of granular fracture in polycrystalline materials. J. Eng. Mater. Technol. 138(4), 041008 (2016)
Gao, Y., Oterkus, S.: Ordinary state-based peridynamic modelling for fully coupled thermoelastic problems. Continuum Mech. Thermodyn. 31(4), 907–937 (2019)
Wang, H., Oterkus, E., Oterkus, S.: Predicting fracture evolution during lithiation process using peridynamics. Eng. Fract. Mech. 192, 176–191 (2018)
Alpay, S., Madenci, E.: Crack growth prediction in fully-coupled thermal and deformation fields using peridynamic theory. In: 54th AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference (p. 1477) (2013)
Yang, Z., Oterkus, E., Nguyen, C.T., Oterkus, S.: Implementation of peridynamic beam and plate formulations in finite element framework. Continuum Mech. Thermodyn. 31(1), 301–315 (2019)
Diyaroglu, C., Oterkus, E., Oterkus, S.: An Euler–Bernoulli beam formulation in an ordinary state-based peridynamic framework. Math. Mech. Solids 24(2), 361–376 (2019)
Diyaroglu, C., Oterkus, S., Oterkus, E., Madenci, E.: Peridynamic modeling of diffusion by using finite-element analysis. IEEE Trans. Components Packag. Manuf. Technol. 7(11), 1823–1831 (2017)
Imachi, M., Tanaka, S., Bui, T.Q., Oterkus, S., Oterkus, E.: A computational approach based on ordinary state-based peridynamics with new transition bond for dynamic fracture analysis. Eng. Fract. Mech. 206, 359–374 (2019)
Imachi, M., Tanaka, S., Ozdemir, M., Bui, T.Q., Oterkus, S., Oterkus, E.: Dynamic crack arrest analysis by ordinary state-based peridynamics. Int. J. Fract. 221, 155–169 (2020)
Basoglu, M.F., Zerin, Z., Kefal, A., Oterkus, E.: A computational model of peridynamic theory for deflecting behavior of crack propagation with micro-cracks. Comput. Mater. Sci. 162, 33–46 (2019)
Dell’Isola, F., Andreaus, U., Placidi, L.: At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: an underestimated and still topical contribution of Gabrio Piola. Math. Mech. Solids 20(8), 887–928 (2015)
Ozdemir, M., Kefal, A., Imachi, M., Tanaka, S., Oterkus, E.: Dynamic fracture analysis of functionally graded materials using ordinary state-based peridynamics. Compos. Struct. 244, 112296 (2020)
Kefal, A., Sohouli, A., Oterkus, E., Yildiz, M., Suleman, A.: Topology optimization of cracked structures using peridynamics. Continuum Mech. Thermodyn. 31(6), 1645–1672 (2019)
Javili, A., Morasata, R., Oterkus, E., Oterkus, S.: Peridynamics review. Math. Mech. Solids 24(11), 3714–3739 (2019)
Acknowledgements
This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-18-1-7004.
Author information
Authors and Affiliations
Corresponding author
Additional information
Communicated by Luca Placidi.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Wang, B., Oterkus, S. & Oterkus, E. Determination of horizon size in state-based peridynamics. Continuum Mech. Thermodyn. 35, 705–728 (2023). https://doi.org/10.1007/s00161-020-00896-y
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00161-020-00896-y