Double horizon peridynamics

Abstract

In this study, a new double horizon peridynamics formulation was introduced by utilising two horizons instead of one as in traditional peridynamics. The new approach allows utilisation of large horizon sizes by controlling the size of the inner horizon size. To demonstrate the capability of the current approach, four different numerical cases were examined by considering static and dynamic conditions, different boundary and initial conditions, and different outer and inner horizon size values. For both static and dynamic cases, it was observed that as the inner horizon decreases, peridynamic solutions converge to classical continuum mechanics solutions even by using a larger horizon size value. Therefore, the proposed approach can serve as an alternative approach to improve computational efficiency of peridynamic simulations by obtaining accurate results with larger horizon sizes.

1. Introduction

Although classical continuum mechanics (CCM) has been widely used for the analysis of materials and structures under external loading conditions, equations of CCM are not suitable to represent situations including discontinuities such as cracks since spatial derivatives in CCM equations are not defined along the crack boundaries. Silling [1] introduced peridynamics with the intention to overcome limitations of CCM formulation. By using integro-differential equations without spatial derivatives, peridynamics has become a powerful approach for predicting crack evolution in materials and structures [2]. Rather than using mesh-based or semi-analytical approaches [3], peridynamics formulation is usually numerically implemented by utilising meshless approach. There has been significant progress on peridynamics especially during recent years. Among numerous studies in the field, Vazic et al. [4] demonstrated the superior capability of peridynamics in fracture analysis by investigating the interaction of macro and micro cracks. De Meo et al. [5] predicted pit-to-crack transition process within peridynamic frameworks. Imachi et al. [6] performed dynamic crack arrest analysis by using peridynamics. Ozdemir et al. [7] used peridynamics to investigate functionally graded materials and their dynamic fracture behaviour. Liu et al. [8] analysed fracture of zigzag graphene sheets by creating an ordinary state-based peridynamic model. Huang et al. [9] extended the capability of peridynamics for visco-hyperelastic materials. Diehl et al. [10] presented a review of benchmark experiments for peridynamic models. Zhou and Yao [11] proposed a new concept of smoothed bond-based peridynamics. Prakash and Stewart [12] demonstrated how to use a multi-threaded method to assemble a sparse stiffness matrix for quasi-static problems of peridynamics. Naumenko et al. [13] compared experimental and peridynamic results for damage patterns in float glass plates. Yan et al. [14] used peridynamics to model soil desiccation. Hidayat et al. [15] provided a review about relationships between meshfree methods and peridynamics. Guski et al. [16] utilised peridynamics to investigate plasma sprayed coatings for Solid Oxide Fuel Cell (SOFC) sealing applications. Kefal et al. [17] demonstrated how to use peridynamics for topology optimisation of cracked structures. There have also been many studies presenting peridynamic formulations for beams and plates to model isotropic materials [18–22], functionally graded materials [23–28], and composite materials [29,30]. Peridynamics has also been extended to model other physical fields. Mikata [31] presented peridynamic formulations for fluid mechanics and acoustics. Diyaroglu et al. [32] introduced a peridynamic wetness approach to be utilised for the analysis of moisture concentration in electronic packages. Martowicz et al. [33] created a thermomechanical model to investigate phase transformation in shape memory alloys.

An important aspect of peridynamics is its length scale parameter, horizon, which does not exist in CCM formulation. Dorduncu and Madenci [34] presented a finite element implementation of ordinary state-based peridynamics having variable horizon. Madenci et al. [35] developed a weak form of bond-associated non-ordinary state-based peridynamics formulation for uniform or non-uniform discretization without experiencing zero-energy mode problem. A comprehensive investigation on how to choose horizon size in bond-based peridynamics and state-based peridynamics is given in Bobaru and Hu [36] and Wang et al. [37], respectively. In this study, a new peridynamic formulation, double horizon peridynamics (DH PD), is proposed, which utilises two horizons for each material point instead of one horizon as in traditional peridynamics formulation. Since simulation-based traditional peridynamics can take significant computational time, for some cases, if a large horizon size and/or refined discretization are utilised, this new formulation can allow improvement of numerical accuracy in peridynamic simulations with less computational time. This formulation is different than “DH PD” formulation, which uses two different horizons for two interacting material points [38–41]. The details of the “DH PD” formulation is given in Section 2. How to treat boundary conditions is explained in Section 3. Analytical solutions for various different boundary conditions under static or dynamic conditions are given in Sections 4 and 5. Some numerical results are presented in Section 6, and conclusions of the study are given in Section 7.

2. DH PD formulation

2.1. Traditional peridynamics

The equation of motion (EoM) for a one-dimensional (1D) rod in CCM can be written as,

$$\rho \frac{{\partial }^{2}u}{\partial t^2}(x,t)=E\frac{{\partial }^{2}u}{\partial x^2}+b(x,t).$$

(1)

peridynamics (PD) EoM for 1D bar can be obtained by converting the local term in equation (1) into nonlocal form by utilising Taylor expansion. Performing Taylor expansion with respect to displacement function $u$ about a particular material point $x$ and absorbing higher order terms within $O(*)$ function yields,

$$u(x+\xi )-u\left(x\right)=\frac{\partial u}{\partial x}\xi +\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}{\xi }^2+\frac{1}{3!}\frac{{\partial }^{3}u}{\partial x^3}{\xi }^3+O\left({\xi }^4\right).$$

(2)

Considering $x$ as fixed, multiplying each term of equation (2) by the influence function $1/\left|\xi \right|$ and integrating over the PD horizon results in,

$$\begin{array}{c}\underset{-\delta }{\overset{\delta }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi =\frac{\partial u}{\partial x}\underset{-\delta }{\overset{\delta }{\int }}\mathrm{sgn}\left(\xi \right)d\xi +\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}\underset{-\delta }{\overset{\delta }{\int }}\left|\xi \right|d\xi \hfill \\ +\frac{1}{3!}\frac{{\partial }^{3}u}{\partial x^3}\underset{-\delta }{\overset{\delta }{\int }}{\xi }^2\mathrm{sgn}\left(\xi \right)d\xi +\underset{-\delta }{\overset{\delta }{\int }}O\left(\left|{\xi }^3\right|\right)d\xi ,\hfill \end{array}$$

(3)

which gives,

$$\frac{{\partial }^{2}u}{\partial x^2}=\frac{2}{{\delta }^2}\underset{-\delta }{\overset{\delta }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +O\left({\delta }^4\right).$$

(4)

Substituting equation (4) into equation (1) yields,

$$\rho \frac{{\partial }^{2}u}{\partial t^2}(x,t)=\frac{2E}{{\delta }^2}\underset{-\delta }{\overset{\delta }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +b(x,t)+O\left({\delta }^4\right).$$

(5)

It reduces to the classic PD EoM if we eliminate the residual term in the above equation such that,

$$\rho \frac{{\partial }^{2}u}{\partial t^2}(x,t)=\frac{2\mathrm{EA}}{{\delta }^2}\underset{-\delta }{\overset{\delta }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +b(x,t).$$

(6)

2.2. DH PD

One can observe that the classical PD EoM, equation (6) identically converges to that of CCM, equation (1) if and only if when the horizon size, $\delta$ , tends to zero. However, horizon is an important parameter which gives PD a nonlocal character and should take up a value of finite length in general. A numerical approach is usually utilised to solve PD EoM. In order to consistently match PD with CCM, a very small horizon size can be chosen, but this brings the price of losing nonlocality and significant computational time. On the contrary, choosing a large horizon may avoid these issues but cause less accurate results. In order to overcome this contradiction, DH PD formulation can be introduced in the PD EoM.

Consider a small inner horizon inside the original horizon as shown in Figure 1.

Figure 1.

The inner horizon ${\mathrm{\Omega }}_{\mathrm{I}}$ and the outer horizon ${\mathrm{\Omega }}_{\mathrm{O}}$ in the double horizon peridynamic formulation.

First, let us consider the integration over the inner horizon ${\mathrm{\Omega }}_{\mathrm{I}}$ . By making analogy with the derivation of equations from equations (2)–(4) by replacing $\delta$ with $\epsilon$ , one can obtain,

$$\frac{{\partial }^{2}u}{\partial x^2}=\frac{2}{{\epsilon }^2}\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +O\left({\epsilon }^4\right).$$

(7)

Note that the inner horizon size $\epsilon$ can be chosen arbitrarily small, such that $O({\epsilon }^4)<<O({\delta }^4)$ .

Next, let us consider the integration over the outer horizon ${\mathrm{\Omega }}_{\mathrm{O}}$ . Multiplying each term in equation (2) by the influence function 1/ $\left|\xi \right|$ gives,

$$\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}=\frac{\partial u}{\partial x}\mathrm{sgn}\left(\xi \right)+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}\left|\xi \right|+\frac{1}{3!}\frac{{\partial }^{3}u}{\partial x^3}{\xi }^2\mathrm{sgn}\left(\xi \right)+O\left(\left|{\xi }^3\right|\right).$$

(8)

Note that when $\xi$ varies over the outer horizon, the residual term ranges as,

$$\begin{array}{cc}\left|O\left(\left|{\epsilon }^3\right|\right)\right|\le \left|O\left(\left|{\xi }^3\right|\right)\right|\le \left|O\left(\left|{\delta }^3\right|\right)\right|\hfill & \forall \xi \in {\mathrm{\Omega }}_O\hfill \end{array}.$$

(9)

In order to reduce the residual to be consistent with that of the inner horizon, multiplying equation (8) by $|{\epsilon }^3/{\xi }^3|$ implies,

$$\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}=\frac{\partial u}{\partial x}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\mathrm{sgn}\left(\xi \right)+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\left|\xi \right|+\frac{1}{3!}\frac{{\partial }^{3}u}{\partial x^3}{\xi }^2\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\mathrm{sgn}\left(\xi \right)+O\left(\left|{\epsilon }^3\right|\right).$$

(10)

Now considering $x$ as fixed and integrating over the outer horizon results in,

$$\begin{array}{c}\multicolumn{1}{c}{\underset{{\mathrm{\Omega }}_O}{\int }\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi =\frac{\partial u}{\partial x}\underset{{\mathrm{\Omega }}_O}{\int }\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\mathrm{sgn}\left(\xi \right)d\xi +\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}\underset{{\mathrm{\Omega }}_O}{\int }\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\left|\xi \right|d\xi }\\ \multicolumn{1}{c}{+\frac{1}{3!}\frac{{\partial }^{3}u}{\partial x^3}\underset{{\mathrm{\Omega }}_O}{\int }{\xi }^2\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\mathrm{sgn}\left(\xi \right)d\xi +\underset{{\mathrm{\Omega }}_O}{\int }O\left(\left|{\epsilon }^3\right|\right)d\xi .}\end{array}$$

(11)

Ignoring residual gives,

$$\frac{{\partial }^{2}u}{\partial x^2}=\frac{\delta }{{\epsilon }^2(\delta -\epsilon )}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ].$$

(12)

Introducing two weight functions ${\omega }_I$ and ${\omega }_O$ for inner and outer horizon, respectively, such that,

$${\omega }_I+{\omega }_O=1,$$

(13)

and combining equation (7) with equation (12) gives,

$$\begin{array}{c}\frac{{\partial }^{2}u}{\partial x^2}={\omega }_I\frac{2}{{\epsilon }^2}\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi \hfill \\ +{\omega }_O\frac{\delta }{{\epsilon }^2(\delta -\epsilon )}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ],\hfill \end{array}$$

(14)

in which the weight functions can be chosen by considering each horizon size in proportion to the total horizon size as,

$${\omega }_I=\frac{\epsilon }{\delta }{\omega }_O=\frac{\delta -\epsilon }{\delta }.$$

(15)

Coupling equation (15) with equation (14) and substituting into equation (1) yields the refined PD EoM for 1D bar as,

$$\begin{array}{c}\rho \frac{{\partial }^{2}u}{\partial t^2}(x,t)=\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi \hfill \\ +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ]+b(x,t).\hfill \end{array}$$

(16)

Note that when the inner horizon size equals to the outer horizon size, $\epsilon =\delta$ , equation (16) reduces to the traditional PD form.

3. Boundary conditions

From a PD point of view, except for the damage situation, each material point must be completely embedded in its PD influence domain. Moreover, for those material points adjacent to the boundary, whose domain is incomplete, it is necessary to introduce fictitious regions outside the boundary such that the completeness of PD EoM is ensured. The width of fictitious region can be chosen as equal to the horizon size $\delta$ , as shown in Figure 2. Two types of BCs and their implementation in PD framework are explained below.

Figure 2.

Real and fictitious domains.

3.1. Fixed boundary conditions

Recall the EoM in classical elasticity:

$$\rho \frac{{\partial }^{2}u}{\partial t^2}(x,t)=E\frac{{\partial }^{2}u(x,t)}{\partial x^2}+b(x,t)0\le x\le L.$$

(17)

Suppose the body is constrained at $x=0$ such that $u(0,t)\equiv 0$ , the representation of equation (1) at this point is,

$$E{\frac{{\partial }^{2}u(x,t)}{\partial x^2}|}_{x=0}=0.$$

(18)

One can obtain by performing central difference for equation (18) that,

$$E\frac{u(-\mathrm{\Delta }x,t)-2u(0,t)+u(\mathrm{\Delta }x,t)}{{\left(\mathrm{\Delta }x\right)}^2}=0,$$

(19)

where $\mathrm{\Delta }x$ is incremental distance. Simplifying equation (19) and swapping the difference notation $\mathrm{\Delta }x$ by $\xi$ leads to,

$$\begin{array}{cc}u(-\xi ,t)=-u(\xi ,t)\hfill & 0<\xi <\delta \hfill \end{array}.$$

(20)

Here, the material point $x=-\xi$ lies in the fictitious region, and equation (20) ensures that the fixed Boundary Condition (BC) is satisfied for PD EoM.

3.2. Neumann boundary conditions

Suppose the body is subjected to a concentrated load of $p(t)$ at $x=L$ . In classical elasticity theory, the boundary condition can be represented as,

$$\sigma (L,t)=E{\frac{\partial u(x,t)}{\partial x}|}_{x=L}=p\left(t\right).$$

(21)

Again, performing central difference with respect to $u$ gives,

$$E\frac{u(L+\mathrm{\Delta }x,t)-u(L-\mathrm{\Delta }x,t)}{2\mathrm{\Delta }x}=p\left(t\right).$$

(22)

After performing some algebra and rearranging the central difference notation according to PD convention, one can obtain that,

$$u(L+\xi ,t)=2\frac{p\left(t\right)}{E}\xi +u(L-\xi ,t)0<\xi <\delta .$$

(23)

Note that, $p(t)$ can be eliminated from equation (23) and absorbed in equation (16) if it is considered as the body force operated by Dirac delta. According to this point of view, equation (23) reduces to,

$$\begin{array}{cc}u(L+\xi ,t)=u(L-\xi ,t)\hfill & 0<\xi <\delta \hfill \end{array},$$

(24)

and equation (24) is called the free boundary condition relationship in the PD framework.

4. Analytic solutions for static problems

4.1. Fixed–fixed rod

Consider a rod subjected to fixed–fixed boundary condition and arbitrary distributed loading as shown in Figure 3.

Figure 3.

Rod subjected to fixed–fixed boundary condition and arbitrary distributed loading.

The length of the bar is denoted as $L$ and Young’s modulus is $E$ . As explained earlier, the PD EoM and BCs can be expressed as,

$$\begin{array}{c}\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi \hfill \\ +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ]+b\left(x\right)=0.\hfill \end{array}$$

(25)

$$\begin{array}{cc}\{\begin{array}{c}u(-\xi )=-u\left(\xi \right)\hfill \\ \multicolumn{1}{c}{u(L+\xi )=-u(L-\xi )}\end{array}\hfill & 0\le \xi \le \delta \hfill \end{array}.$$

(26a, b)

According to the completeness characteristic of a trigonometric system, it is reasonable to assume that the displacement function accommodates in the vector space spanned by trigonometric functions, i.e.,

$$u\left(x\right)\in \mathrm{Span}\left\{\begin{array}{ccccc}\frac{1}{\sqrt{2}},\hfill & \cos \left(s_{1}x\right),\hfill & \cos \left(s_{2}x\right),\hfill & \cos \left(s_{3}x\right),\hfill & ...\hfill \\ \multicolumn{1}{c}{}& \sin \left(s_{1}x\right),\hfill & \sin \left(s_{2}x\right),\hfill & \sin \left(s_{3}x\right),\hfill & ...\hfill \end{array}\right\},$$

(27)

which can be expressed as a linear combination in terms of the bases as follows:

$$u\left(x\right)=\frac{a_0}{2}+\sum _{n=1}^{\infty }{a}_n\cos \left(s_{n}x\right)+b_n\sin \left(s_{n}x\right).$$

(28)

Substituting equation (28) into equation (26a) yields,

$$\frac{a_0}{2}+\sum _{n=1}^{\infty }{a}_n\cos (-s_n\xi )+b_n\sin (-s_n\xi )=-[\frac{a_0}{2}+\sum _{n=1}^{\infty }{a}_n\cos \left(s_n\xi \right)+b_n\sin \left(s_n\xi \right)].$$

(29)

Without considering rigid body motion, one can obtain that,

$$u\left(x\right)=\sum _{n=1}^{\infty }{b}_n\sin \left(s_{n}x\right).$$

(30)

Substituting equation (30) into equation (26b) yields,

$$\sum _{n=1}^{\infty }{b}_n\sin \left[s_n(L+\xi )\right]=-\sum _{n=1}^{\infty }{b}_n\sin \left[s_n(L-\xi )\right],$$

(31)

which gives,

$$\sum _{n=1}^{\infty }{b}_n\sin \left(s_{n}L\right)\cos \left(s_n\xi \right)=0\Rightarrow s_n=\frac{n\pi }{L}.$$

(32)

Plugging equation (32) back into equation (30) implies,

$$u\left(x\right)=\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi }{L}x\right).$$

(33)

One can obtain that,

$$\begin{array}{c}\multicolumn{1}{c}{\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi =\sum _{n=1}^{\infty }{b}_n\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[\sin \left(\frac{n\pi (x+\xi )}{L}\right)-\sin \left(\frac{n\pi x}{L}\right)]d\xi }\\ \multicolumn{1}{c}{=\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi x}{L}\right)\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[\cos \left(\frac{n\pi \xi }{L}\right)-1]d\xi .}\end{array}$$

(34)

Similarly,

$$\begin{array}{c}\multicolumn{1}{c}{\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi }\\ \multicolumn{1}{c}{=\underset{-\delta }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi -\underset{-\epsilon }{\overset{\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi }\\ \multicolumn{1}{c}{=\sum _{n=1}^{\infty }{b}_n\underset{-\delta }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|[\sin \left(\frac{n\pi }{L}(x+\xi )\right)-\sin \left(\frac{n\pi }{L}x\right)]d\xi -\sum _{n=1}^{\infty }{b}_n\underset{-\epsilon }{\overset{\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|[\sin \left(\frac{n\pi }{L}(x+\xi )\right)-\sin \left(\frac{n\pi }{L}x\right)]d\xi }\\ \multicolumn{1}{c}{=\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi x}{L}\right)\underset{-\delta }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|[\cos \left(\frac{n\pi \xi }{L}\right)-1]d\xi -\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi x}{L}\right)\underset{-\epsilon }{\overset{\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|[\cos \left(\frac{n\pi \xi }{L}\right)-1]d\xi }\\ \multicolumn{1}{c}{=\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi x}{L}\right)\{\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{n\pi \xi }{L}-1)d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{n\pi \xi }{L}-1)d\xi \}.}\end{array}$$

(35)

Substituting equations (34) and (35) into equation (25) yields,

$$\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi x}{L}\right)\{\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[\cos \left(\frac{n\pi \xi }{L}\right)-1]d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{n\pi \xi }{L}-1)d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{n\pi \xi }{L}-1)d\xi ]\}=-b\left(x\right).$$

(36)

According to the orthogonal property of trigonometric functions, the coefficients can be determined as,

$$b_n=-\frac{2}{L}\frac{\underset{0}{\overset{L}{\int }}b\left(x\right)\sin \frac{n\pi x}{L}\mathrm{dx}}{\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}(\cos \frac{n\pi \xi }{L}-1)d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{n\pi \xi }{L}-1)d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{n\pi \xi }{L}-1)d\xi ]}.$$

(37)

Hence, the analytical PD solution is,

$$u\left(x\right)=\sum _{n=1}^{\infty }\frac{2}{L}\frac{\underset{0}{\overset{L}{\int }}b\left(x\right)\sin \frac{n\pi x}{L}\mathrm{dx}\sin \left(\frac{n\pi x}{L}\right)}{\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}(1-\cos \frac{n\pi \xi }{L})d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{n\pi \xi }{L})d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{n\pi \xi }{L})d\xi ]}.$$

(38)

Note that when $\epsilon =\delta$ , equation (38) reduces to analytical solution for classical PD theory as,

$$u\left(x\right)=\sum _{n=1}^{\infty }\frac{{\delta }^2}{\mathrm{EL}}\frac{\underset{0}{\overset{L}{\int }}b\left(x\right)\sin \frac{n\pi x}{L}\mathrm{dx}}{\underset{-\delta }{\overset{\delta }{\int }}\frac{1}{\left|\xi \right|}(1-\cos \frac{n\pi \xi }{L})d\xi }\sin \left(\frac{n\pi x}{L}\right).$$

(39)

4.2. Fixed–free rod

Consider a rod subjected to fixed–free boundary condition and arbitrary distributed loading as shown in Figure 4.

Figure 4.

Rod subjected to fixed–free boundary condition and arbitrary distributed loading.

The PD EoM and BCs for this case can be expressed as,

$$\begin{array}{c}\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ]+b\left(x\right)=0.\hfill \end{array}$$

(40)

$$\begin{array}{cc}\{\begin{array}{c}u(-\xi )=-u\left(\xi \right)\hfill \\ \multicolumn{1}{c}{u(L+\xi )=u(L-\xi )}\end{array}\hfill & 0\le \xi \le \delta \hfill \end{array}.$$

(41a, b)

Again, suppose the displacement function belongs to the vector space spanned by bases of trigonometric functions and hence can be expressed as,

$$u\left(x\right)=\frac{a_0}{2}+\sum _{n=1}^{\infty }{a}_n\cos \left(s_{n}x\right)+b_n\sin \left(s_{n}x\right).$$

(42)

Using equations (41a, b), one can obtain,

$$s_n=\frac{(2n-1)\pi }{2L},$$

(43)

and,

$$u\left(x\right)=\sum _{n=1}^{\infty }{b}_n\sin \frac{(2n-1)\pi x}{2L}.$$

(44)

Plugging equation (44) back into equation (40) and performing algebraic simplifications results in,

$$\begin{array}{c}\sum _{n=1}^{\infty }{b}_n\sin \frac{(2n-1)\pi x}{2L}\left\{\begin{array}{c}\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi ]\hfill \end{array}\right\}\hfill \\ =-b\left(x\right),\hfill \end{array}$$

(45)

where the coefficients can be determined as,

$$b_n=-\frac{2}{L}\frac{\underset{0}{\overset{L}{\int }}b\left(x\right)\sin \frac{(2n-1)\pi x}{2L}\mathrm{dx}}{\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi ]}.$$

(46)

Substituting equation (46) into equation (44) yields the PD analytical solution as,

$$\begin{array}{c}u\left(x\right)=-\sum _{n=1}^{\infty }\frac{2}{L}\frac{\underset{0}{\overset{L}{\int }}b\left(x\right)\sin \frac{(2n-1)\pi x}{2L}\mathrm{dx}\sin \frac{(2n-1)\pi x}{2L}}{\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}+\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi ]}.\hfill \end{array}$$

(47)

In particular, when $\epsilon =\delta$ , equation (47) reduces to the analytical solution for traditional PD model as,

$$u\left(x\right)=-\sum _{n=1}^{\infty }\frac{{\delta }^2}{\mathrm{LE}}\frac{\underset{0}{\overset{L}{\int }}b\left(x\right)\sin \frac{(2n-1)\pi x}{2L}\mathrm{dx}}{\underset{-\delta }{\overset{\delta }{\int }}\frac{1}{\left|\xi \right|}(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi }\sin \frac{(2n-1)\pi x}{2L}.$$

(48)

5. Analytical solutions for free vibration condition

5.1. Fixed–fixed rod

Let the PD EoM, BCs, and Initial Condition (IC)s to be expressed for the fixed–fixed rod as,

$$\begin{array}{c}\rho \frac{{\partial }^{2}u}{\partial t^2}(x,t)=\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi \hfill \\ +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ].\hfill \end{array}$$

(49)

$$\begin{array}{cc}\mathrm{BCs}:\{\begin{array}{c}u(-\xi ,t)=-u(\xi ,t)\hfill \\ \multicolumn{1}{c}{u(L+\xi ,t)=-u(L-\xi ,t)}\end{array}\hfill & 0\le \xi \le \delta \hfill \end{array}.$$

(50a, b)

$$\mathrm{ICs}:\{\begin{array}{c}u(x,0)=u_0\left(x\right)\hfill \\ \multicolumn{1}{c}{\frac{\partial u}{\partial t}(x,0)=v_0\left(x\right)}\end{array}.$$

(51a, b)

By using separation of variables approach, the displacement function can be decomposed as,

$$u(x,t)=X\left(x\right)T\left(t\right).$$

(52)

Hence, equation (49) can be rewritten as,

$$\rho X\frac{{\partial }^{2}T}{\partial t^2}=T\{\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{X(x+\xi )-X\left(x\right)}{\left|\xi \right|}d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{X(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{X(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ]\}.$$

(53)

Isolating variables yields,

$$\frac{\rho }{T}\frac{{\partial }^{2}T}{\partial t^2}=\frac{1}{X}\{\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{X(x+\xi )-X\left(x\right)}{\left|\xi \right|}d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{X(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{X(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ]\}=-\lambda ,$$

(54)

which gives,

$$\begin{array}{c}\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{X(x+\xi )-X\left(x\right)}{\left|\xi \right|}d\xi \hfill \\ +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{X(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{X(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ]=-\lambda X,\hfill \end{array}$$

(55a)

and,

$$\rho \frac{{\partial }^{2}T}{\partial t^2}+\lambda T=0.$$

(55b)

By comparing equation (55a) with equation (25), we can consider $X(x)$ as an analogue to $u(x)$ and $\lambda X(x)$ to $b(x)$ . Thus, associating with equations (33) and (36), one can obtain,

$$X\left(x\right)=\sum _{n=1}^{\infty }{b}_n\sin \frac{n\pi x}{L},$$

(56)

and,

$$\begin{array}{c}\sum _{n=1}^{\infty }{b}_n\sin \left(\frac{n\pi x}{L}\right)\{\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[\cos \left(\frac{n\pi \xi }{L}\right)-1]d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{n\pi \xi }{L}-1)d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{n\pi \xi }{L}-1)d\xi ]\}\hfill \\ =-\sum _{n=1}^{\infty }{\lambda }_n{b}_n\sin \frac{n\pi x}{L},\hfill \end{array}$$

(57)

which results in the PD “eigenvalues” as,

$$\begin{array}{c}{\lambda }_n=\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[1-\cos \left(\frac{n\pi \xi }{L}\right)]d\xi \hfill \\ +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{n\pi \xi }{L})d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{n\pi \xi }{L})d\xi ].\hfill \end{array}$$

(58)

The general solution to equation (55b) can be expressed as,

$$T_n\left(t\right)=A_n^{*}\cos \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)+B_n^{*}\sin \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right).$$

(59)

Substituting equations (59) and (56) into equation (52) and rearranging the coefficient notations yields the general PD solution as,

$$u(x,t)=\sum _{n=1}^{\infty }[A_n\cos \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)+B_n\sin \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)]\sin \frac{n\pi x}{L}.$$

(60)

Differentiating equation (60) with respect to time yields,

$$\frac{\partial u}{\partial t}(x,t)=\sum _{n=1}^{\infty }\sqrt{\frac{{\lambda }_n}{\rho }}[-A_n\sin \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)+B_n\cos \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)]\sin \frac{n\pi x}{L}.$$

(61)

Applying the ICs gives,

$$\begin{array}{ccccc}u(x,0)=u_0\left(x\right)\hfill & \Rightarrow \hfill & \sum _{n=1}^{\infty }{A}_n\sin \frac{n\pi x}{L}=u_0\left(x\right)\hfill & \Rightarrow \hfill & A_n=\frac{2}{L}\underset{0}{\overset{L}{\int }}{u}_0\left(x\right)\sin \frac{n\pi x}{L}\mathrm{dx}\hfill \end{array}.$$

(62a)

$$\begin{array}{ccccc}\frac{\partial u}{\partial t}(x,0)=v_0\left(x\right)\hfill & \Rightarrow \hfill & \sum _{n=1}^{\infty }\sqrt{\frac{{\lambda }_n}{\rho }}{B}_n\sin \frac{n\pi x}{L}=v_0\hfill & \Rightarrow \hfill & B_n=\frac{2}{L}\sqrt{\frac{\rho }{{\lambda }_n}}\underset{0}{\overset{L}{\int }}{v}_0\left(x\right)\sin \frac{n\pi x}{L}\mathrm{dx}\hfill \end{array}.$$

(62b)

As a summary, the PD analytical solution for free vibration of fixed–fixed rods can be summarised as,

$$u(x,t)=\sum _{n=1}^{\infty }[A_n\cos \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)+B_n\sin \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)]\sin \frac{n\pi x}{L}.$$

(63a)

$$\begin{array}{c}{\lambda }_n=\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[1-\cos \left(\frac{n\pi \xi }{L}\right)]d\xi \hfill \\ +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{n\pi \xi }{L})d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{n\pi \xi }{L})d\xi ].\hfill \end{array}$$

(63b)

$$A_n=\frac{2}{L}\underset{0}{\overset{L}{\int }}{u}_0\left(x\right)\sin \frac{n\pi x}{L}\mathrm{dx}.$$

(63c)

$$B_n=\frac{2}{L}\sqrt{\frac{\rho }{{\lambda }_n}}\underset{0}{\overset{L}{\int }}{v}_0\left(x\right)\sin \frac{n\pi x}{L}\mathrm{dx}.$$

(63d)

In particular, when $\epsilon =\delta$ , the PD “eigenvalues” reduce to,

$${\lambda }_n=\frac{2E}{{\delta }^2}\underset{-\delta }{\overset{\delta }{\int }}\frac{1}{\left|\xi \right|}[1-\cos \left(\frac{n\pi \xi }{L}\right)]d\xi ,$$

(63e)

which is valid for traditional PD models.

5.2. Fixed–free rod

Let the PD EoM, BCs, and ICs for the fixed–free rod be expressed as,

$$\begin{array}{c}\rho \frac{{\partial }^{2}u}{\partial t^2}(x,t)=\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi \hfill \\ +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{u(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ].\hfill \end{array}$$

(64)

$$\begin{array}{cc}\mathrm{BCs}:\{\begin{array}{c}u(-\xi ,t)=-u(\xi ,t)\hfill \\ \multicolumn{1}{c}{u(L+\xi ,t)=u(L-\xi ,t)}\end{array}\hfill & 0\le \xi \le \delta \hfill \end{array}.$$

(65a, b)

$$\mathrm{ICs}:\{\begin{array}{c}u(x,0)=u_0\left(x\right)\hfill \\ \multicolumn{1}{c}{\frac{\partial u}{\partial t}(x,0)=v_0\left(x\right)}\end{array}.$$

(66a, b)

By using separation of variables approach, i.e.,

$$u(x,t)=X\left(x\right)T\left(t\right),$$

(67)

two characteristic functions can be obtained by substituting equation (67) into equation (64) as,

$$\begin{array}{c}\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{X(x+\xi )-X\left(x\right)}{\left|\xi \right|}d\xi \hfill \\ +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{X(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^3}\right|\frac{X(x+\xi )-u\left(x\right)}{\left|\xi \right|}d\xi ]=-\lambda X,\hfill \end{array}$$

(68a)

and,

$$\rho \frac{{\partial }^{2}T}{\partial t^2}+\lambda T=0.$$

(68b)

By comparing equation (68a) with equation (40), we can consider $X(x)$ as an analogue to $u(x)$ and $\lambda X(x)$ to $b(x)$ . Thus, associating with equations (44) and (45), one can obtain,

$$X\left(x\right)=\sum _{n=1}^{\infty }{b}_n\sin \frac{(2n-1)\pi x}{2L},$$

(69a)

and,

$$\begin{array}{c}\sum _{n=1}^{\infty }{b}_n\sin \frac{(2n-1)\pi x}{2L}\left\{\begin{array}{c}\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[\cos \frac{(2n-1)\pi \xi }{2L}-1]d\xi +\hfill \\ \multicolumn{1}{c}{\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(\cos \frac{(2n-1)\pi \xi }{2L}-1)d\xi ]}\end{array}\right\}\hfill \\ =-\sum _{n=1}^{\infty }{\lambda }_n{b}_n\sin \frac{(2n-1)\pi x}{2L},\hfill \end{array}$$

(69b)

which results in the PD “eigenvalues” as,

$${\lambda }_n=\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[1-\cos \frac{(2n-1)\pi \xi }{2L}]d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{(2n-1)\pi \xi }{2L})d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{(2n-1)\pi \xi }{L})d\xi ].$$

(70)

Similar to the derivations of the previous case, the refined PD analytical solution can be obtained as,

$$u(x,t)=\sum _{n=1}^{\infty }[A_n\cos \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)+B_n\sin \left(\sqrt{\frac{{\lambda }_n}{\rho }}t\right)]\sin \frac{(2n-1)\pi x}{2L}.$$

(71a)

$${\lambda }_n=\frac{2E}{\delta \epsilon }\underset{-\epsilon }{\overset{\epsilon }{\int }}\frac{1}{\left|\xi \right|}[1-\cos \frac{(2n-1)\pi \xi }{2L}]d\xi +\frac{E}{{\epsilon }^2}[\underset{-\delta }{\overset{-\epsilon }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{(2n-1)\pi \xi }{2L})d\xi +\underset{\epsilon }{\overset{\delta }{\int }}\left|\frac{{\epsilon }^3}{{\xi }^4}\right|(1-\cos \frac{(2n-1)\pi \xi }{L})d\xi ].$$

(71b)

$$A_n=\frac{2}{L}\underset{0}{\overset{L}{\int }}{u}_0\left(x\right)\sin \frac{(2n-1)\pi x}{2L}\mathrm{dx}.$$

(71c)

$$B_n=\frac{2}{L}\sqrt{\frac{\rho }{{\lambda }_n}}\underset{0}{\overset{L}{\int }}{v}_0\left(x\right)\sin \frac{(2n-1)\pi x}{2L}\mathrm{dx}.$$

(71d)

In particular, when $\epsilon =\delta$ , the PD “eigenvalues” reduce to,

$${\lambda }_n=\frac{2E}{{\delta }^2}\underset{-\delta }{\overset{\delta }{\int }}\frac{1}{\left|\xi \right|}[1-\cos \frac{(2n-1)\pi \xi }{2L}]d\xi ,$$

(71e)

which is valid for the traditional PD model.

6. Numerical results

6.1. Fixed–fixed rod under static conditions

In the first case, a 1D rod with a length of L = 1 m is subjected to a body load of b(x) = 1000(1 –x)¹⁰ N/m³ under static conditions. Both edges of the rod are fixed (fixed–fixed). Elastic modulus of the rod is specified as E = 200 GPa.

Variation of the displacement field along the rod by changing inner horizon size, ε = δ, δ/2, δ/5, δ/10, for a horizon size value of δ = 0.1 m are given in Figure 5. As shown in this figure, as the size of the inner horizon decreases, DH PD solution converges to the CCM solution.

Figure 5.

(a) Variation of the displacement field along the rod by changing inner horizon size, ε, for a horizon value of δ = 0.1 m and(b) Zoomed view.

In Figure 6, two different scenarios are compared. First, as the horizon size decreases from δ = 0.1 m to δ = 0.05 m, peridynamic solution converges to the CCM solution as expected. On the contrary, even if a larger horizon size is used, i.e., δ = 0.1 m, the new DH PD approach provides closer results to CCM solution for an internal horizon size of ε = 0.01 m.

Figure 6.

(a) Variation of the displacement field along the rod by changing the horizon size and using a traditional or double horizon approach and (b) Zoomed view.

6.2. Fixed–free rod under static conditions

In the second case, the same properties are considered as in the first case except the right edge is left as free condition whereas the left edge is fixed under static conditions. First, the effect of the internal horizon size is investigated by changing the inner horizon size as ε = δ, δ/2, δ/5, δ/10, for a horizon size value of δ = 0.1 m. As shown in Figure 7, as the inner horizon size decreases, DH PD results converge to the CCM solution based on the variation of displacements along the rod.

Figure 7.

(a) Variation of the displacement field along the rod by changing inner horizon size, ε, for a horizon value of δ = 0.1 m and (b) Zoomed view.

Next, two different horizon size values are considered, δ = 0.1 m, 0.05 m, by using traditional peridynamics and as expected, for the smaller horizon size value, peridynamic results are closer to CCM results (see Figure 8). On the contrary, if DH PD formulation is utilised, better agreement can be obtained even for a larger horizon size δ = 0.1 m by considering the inner horizon size as ε = 0.01 m.

Figure 8.

(a) Variation of the displacement field along the rod by changing the horizon size and using a traditional or double horizon approach and (b) Zoomed view.

6.3. Fixed–fixed rod under free vibration conditions

In the third case, the capability of the DH PD formulation under the dynamic free vibration conditions is investigated. The length of the rod is the same as the two previous cases. The horizon size is specified as δ = 0.1 m. Elastic modulus and density are given as E = 200 GPa and ρ = 7850 kg/m³. Initial conditions are imposed as,

$$u(x,0)=0.01x^5{(x-1)}^5\mathrm{m}.$$

(72a)

$$\frac{\partial u}{\partial t}(x,0)=10x^5{(1-x)}^5\mathrm{m}/\mathrm{s}.$$

(72b)

As shown in Figure 9, the variation of the displacement at the centre of the rod (x = 0.5 m) as the time progresses is obtained by changing the internal horizon value as ε = δ, δ/5, δ/10, δ/100. Similar to the static cases, as the internal horizon size decreases, DH PD results capture the CCM results very well for the same outer horizon size value which demonstrates the capability of the double horizon approach.

Figure 9.

(a) Variation of the displacement field at the centre of the rod (x = 0.5 m) by changing inner horizon size, ε, for a horizon value of δ = 0.1 m and (b) Zoomed view.

6.4. Fixed–free rod under free vibration conditions

In the final numerical case, the same properties are considered as in the previous case except by assigning free condition for the right edge and fixed condition for the left edge under free vibration conditions and imposing the initial conditions as,

$$u(x,0)=0.01x^5\mathrm{m}.$$

(73a)

$$\frac{\partial u}{\partial t}(x,0)=0\mathrm{m}/\mathrm{s}.$$

(73b)

As in the previous case, four different internal horizon values are considered as ε = δ, δ/5, δ/10, δ/100, and the variation of the displacement at the centre of the rod (x = 0.5 m) is obtained as the time progresses as shown in Figure 10. Similar to the previous case, DH PD results match very well with the CCM results as the internal horizon size decreases.

Figure 10.

(a) Variation of the displacement field at the centre of the rod (x = 0.5 m) by changing inner horizon size, ε, for a horizon value of δ = 0.1 m and (b) Zoomed view.

7. Conclusion

In this study, a new DH PD formulation was introduced by utilising two horizons instead of one as in traditional peridynamics. To demonstrate the capability of the current approach, four different numerical cases were examined by different boundary (fixed–fixed and fixed–free) and initial conditions under static and dynamic conditions for different outer and inner horizon size values. For both static and dynamic cases, it was observed that as the inner horizon decreases, peridynamic solutions converge to CCM solutions even by using a larger horizon size value. Therefore, the proposed approach can serve as an alternative approach to improve computational efficiency of peridynamic simulations by obtaining accurate results with larger horizon sizes. Potential future directions can be extending the current formulation to two-dimensional and three-dimensional models, and problems including fracture.