Boundary Layer Effects in a Finite Linearly Elastic Peridynamic Bar

Abstract

The peridynamic theory is an extension of the classical continuum mechanics theory. The peridynamic governing equations involve integrals of interaction forces between near particles separated by finite distances. These forces depend upon the relative displacements between material points within a body. On the other hand, the classical governing equations involve the divergence of a tensor field, which depends upon the spatial derivatives of displacements. Thus, the peridynamic governing equations are valid not only in the interior of a body, but also on its boundary, which may include a Griffith crack, and on interfaces between two bodies with different mechanical properties. Near the boundary, the solution of a peridynamic problem may be very different from the classical solution. In this work, we investigate the behavior of the displacement field of a unidimensional linearly elastic bar of length L near its ends in the context of the peridynamic theory. The bar is in equilibrium without body force, is fixed at one end, and is subjected to an imposed displacement at the other end. The bar has micromodulus C, which is related to the Young's modulus E in the classical theory and is given by different expressions found in the literature. We find that, depending on the expression of C, the displacement field may be singular near the ends, which is in contrast to the linear behavior of the displacement field observed in the classical linear elasticity. In spite of the above, we show that the peridynamic displacement field converges to its classical counterpart as a length scale, called peridynamic horizon, tends to zero.

Keywords
Peridynamics; Linear Elasticity; Fredholm Integral Equation; Constitutive Modeling

1 INTRODUCTION

The effects of long-range forces in the peridynamic theory have motivated the study of one-dimensional bars of infinite length made of linearly elastic peridynamic materials. The corresponding problems are formulated in terms of linear Fredholm integral equations, which are solved by means of Fourier transform techniques; thus, providing Green's functions for general loading applications. These solutions exhibit features that are not found in their counterparts in the classical linear elasticity theory and are shown to converge to these classical solutions in the limit of short-range forces.

One feature of interest was observed by Silling et al. (2003) in the analysis of a bar with constant micromodulus in $\left(-\delta ,\delta \right)$, where $\delta$ is the horizon, under a single concentrated load. Even though the displacement field obtained from the solution of the corresponding static problem in the classical linear theory is bounded and continuous at all points of the bar, its counterpart in the linear peridynamic theory is unbounded at the point where the force is applied.

Additional studies on one-dimensional bars of infinite length made of linearly elastic peridynamic materials consist of the works of Weckner and Abeyaratne (2005) on the dynamics of bars considering four different types of micromoduli, Mikata (2012), Weckner et al. (2009), and Wang et al. (2017a, 2017b) on static and dynamic Green's functions, Bobaru et al. (2009) on adaptive refinement strategies for the numerical simulation of the corresponding peridynamic problems, and Seleson and Parks (2011) on the numerical investigation of the role of spherical influence functions in the propagation of waves.

Investigation of peridynamic bar of finite length has an additional complication due to the imposition of boundary conditions, which cannot be imposed locally, as in the case of the classical theory, due to the long-range interactions between material points of the bar. To circumvent this complication, Nishawala and Ostoja-Starzewski (2017) assume a particular form for the solution of the peristatic problem of a linearly elastic bar of length $2L$ and find the body force that is necessary to satisfy the equilibrium equation. Their peridynamic solution has the same general form of the classical elastostatic solution of a linearly elastic bar of length $2L$ being subjected to a constant body force density and axial end loads at the ends of the bar.

On the other hand, Chen et al. (2016) have attacked the above complication directly by considering a class of peridynamic kernels parameterized by an exponent $n$ in the context of a bond-based formulation of the equation of motion with no body force. They have investigated the convergence of the peridynamic model to the classical solution of an elastic wave propagation problem consisting of a Gaussian wave propagating along a one-dimensional bar with fixed ends. They have considered two types of Dirichlet boundary conditions: A local condition, which consists of keeping the displacement at the end nodes equal to zero at all times and a nonlocal condition, which consists of adding extra nodes to the ends of the bar and extending the classical solution and its time derivative beyond the ends of the bar. In general, they obtain much better convergence results with the second type of boundary condition. They note, however, that this second type may not be implemented in two- and three-dimensional problems and in problems where new surfaces are created.

Recently, Queiruga and Moridis (2017) have considered three state-based peridynamic models and five influence functions in numerical experiments designed to investigate the convergence of these models to the classical linear elasticity model. There is no attempt to study the effect of nonlocality near the boundary and the authors acknowledge the need of future work to impose proper Dirichlet boundary conditions.

In this work, we consider finite linear elastic peridynamic bars being pulled at the ends and having different micromoduli that were either introduced or adapted from expressions proposed by the aforementioned authors. Differently from Chen et al. (2016), the boundary conditions are imposed on extended parts of the bar (and not as an extension of the classical solution). Depending on the micromodulus, numerical results indicate that the displacement field is discontinuous, as in the case of the constant micromodulus mentioned above, or, has unbounded derivatives at the ends, which is in contrast to the homogeneous deformation of the bar predicted by the classical linear theory. Nevertheless, in all cases, the numerical results converge to results obtained from the classical linear theory as the horizon $\delta$ tends to zero.

In Section 2 we formulate the one-dimensional problem of motion of a linear elastic peridynamic bar and find a relation between the micromodulus of the peridynamic material and the Young’s modulus of classical linear elasticity. In Section 3 we formulate the problem of an elastic peridynamic bar of finite length in equilibrium without body force being subjected to imposed displacements at its ends. The problem is then recast in terms of an inhomogeneous Fredholm equation. In Section 4 we present the numerical scheme used to obtain approximate solutions for different expressions of the micromodulus. These solutions are highly nonlinear in a boundary layer near the ends and, as observed earlier, are discontinuous, or, singular at these ends. We then concentrate our attention on a particular micromodulus that yields a singular behavior of the solution near the ends. Using this micromodulus, we study both convergence of the proposed numerical scheme and convergence of the nonlocal model to the classical linear elastic model as the horizon $\delta$ tends to zero. In Section 5 we discuss the results of this work and present some concluding remarks.

2 ONE-DIMENSIONAL LINEARLY ELASTIC PERIDYNAMIC MODEL

Let $\mathcal{B}$be the reference configuration of a continuous body in its natural state and let ${\mathcal{N}}_{\delta }(x)$ be a $\delta$-neighborhood of $x\in \mathcal{B}$, where $\delta >0$ is called horizon. Let also $x\text{'}\in {\mathcal{N}}_{\delta }(x)$ and define $\xi \triangleq x\text{'}-x$, which is called a bond of $x\text{'}$ to $x$. The collection of all bonds to $x$ is ${\mathcal{H}}_{\delta }(x)$ and $\mathcal{A}(x,t)\langle ·\rangle :{\mathcal{H}}_{\delta }(x)\to {\mathcal{L}}_m$ is a peridynamic state of order $m$ at $x\in \mathcal{B}$ and time $t\in {ℝ}^{+}$, where ${\mathcal{L}}_m$ is the set of all tensors of order $m=\mathrm{0,1,...}$

The governing equation of motion in the peridynamic state-based theory is given by (Silling et al., 2007)

${\displaystyle {\int }_{{\mathcal{N}}_{\delta }(x)}\left[\underset{\_}{T}(x,t)\langle \xi \rangle -\underset{\_}{T}(x^{\prime },t)\langle -\xi \rangle \right]}d{V}_{x\text{'}}+b(x,t)=\rho (x)\ddot{u}(x,t),\forall x\in \mathcal{B},$ (1)

where $b:\mathcal{B}\times {ℝ}^{+}\to {ℝ}^3$ is a prescribed body force density field, $\rho :\mathcal{B}\to {ℝ}^{+}$ is the mass density in the reference configuration, $u:\mathcal{B}\times {ℝ}^{+}\to {ℝ}^3$ is the displacement vector field, $\ddot{u}={\partial }^{2}u/\partial t^2$ is the acceleration field, and $\underset{\_}{T}(x,t)\langle \xi \rangle$ is a force vector state at $(x,t)\in \mathcal{B}\times {ℝ}^{+}$ evaluated on bond $\xi$. To be clear, if $\chi :\mathcal{B}\times {ℝ}^{+}\to {ℝ}^3$ is the deformation field, such that $u(x,t)=\chi (x,t)-x$, then $\underset{\_}{T}(x,t)\langle \xi \rangle$ is the force acting at the point $\chi (x,t)$ in the distorted configuration $\chi (\mathcal{B},t)$ and has the dimension of force per unit reference volume squared. Below, we omit writing the time variable $t$.

Next, consider a cylindrical bar of cross-sectional area $A$ in its reference configuration. Consider also that ${\mathcal{N}}_{\delta }(x)$ is a penny-shaped region of length $2\delta$ and cross-sectional area $A_{\delta }$, which could be taken as $\pi {\delta }^2$. Both the axis of the penny-shaped region and the generators of the cylindrical bar are parallel to the $x$-axis and $x$ is at the center of the penny-shaped region. Integrate the equation of motion (1) over the cross section at an arbitrary point $x$ having axial coordinate $x$ and assume that the only equation that does not vanish identically is given by

${\displaystyle {\int }_{-\delta }^{\delta }\left[\underset{\_}{T}(x)\langle \xi \rangle -\underset{\_}{T}(x\text{'})\langle -\xi \rangle \right]}d\xi +b(x)=\rho \ddot{u}(x),\forall x\in \Omega ,$ (2)

where $\xi \triangleq x\text{'}-x$, $\Omega \triangleq [a,b]$ for some scalars $a$, $b$, such that $L\triangleq b-a>0$ is the length of the bar in its reference configuration,

$u(x)\triangleq \frac{1}{A}{\displaystyle {\int }_{\mathcal{A}}u}(x)·edA,b(x)\triangleq \frac{1}{A}{\displaystyle {\int }_{\mathcal{A}}b}(x)·edA,$

with $e$ being a unit vector parallel to $x$-axis, and

$\underset{\_}{T}(x)\langle \xi \rangle \triangleq \frac{1}{A}{\displaystyle {\int }_{\mathcal{A}}{\displaystyle {\int }_{{\mathcal{A}}_{\delta }}\underset{\_}{T}}}(x)\langle \xi \rangle ·ed{A}_{\delta }dA.$

The above assumption is quite strong because $\underset{\_}{T}(x)\langle \xi \rangle$ depends on the direction of $\xi$ and there is no reason to believe that this dependence is of a particular form. The situation is different in the classical theory, where it is possible to derive a one-dimensional theory of bars by making plausible assumptions in the three-dimensional theory. Nevertheless, we follow the approach of Silling et al. (2003) and use the equation (2) to formulate the problem of motion of a peridynamic bar. A solution of this one-dimensional problem gives valuable insight on the analysis of higher-dimensional problems.

Let $\widehat{\underset{\_}{A}}\left[\underset{\_}{u}\right]$ denote the dependence of the state $\underset{\_}{A}$on the state $\underset{\_}{u}$ and define the difference displacement quotient scalar state at $x\in \Omega$ through

$\underset{\_}{h}\langle \xi \rangle =\frac{u(x+\xi )-u(x)}{\left|\xi \right|},$

where $|·|$ denotes the absolute value of the argument. Then, a constitutive relation for the force state $\underset{\_}{T}(x)$ in (2) is given by

$\underset{\_}{T}(x)=\widehat{\underset{\_}{T}}(x)\left[\underset{\_}{h}\right],$ (3)

where $\widehat{\underset{\_}{T}}(x)\left[\underset{\_}{h}\right]$ is the response function state of the peridynamic material.

In the case of a linear elastic material, we are motivated by the three-dimensional constitutive theory introduced by Aguiar and Fosdick (2014) and propose the response function

$\widehat{\underset{\_}{T}}(x)\left[\underset{\_}{h}\right]\langle \xi \rangle =c_1{\displaystyle {\int }_{-\delta }^{\delta }\frac{\omega (|\xi |,|\eta |)}{|\xi {|}^2}}\left[u(x+\xi )-u(x)\right]d\eta ,$ (4)

where $c_1$ is a material constant and $\omega (|\xi |,|\eta |)$ is a symmetric influence function. A more general expression than (4) could be proposed. The expression (4) is, however, simple enough to grasp the essential features of the modeling process.

Substituting (3) together with (4) into the equation of motion (2) and considering that the bar is in equilibrium, we obtain

$2c_1{\displaystyle {\int }_{-\delta }^{\delta }{\displaystyle {\int }_{-\delta }^{\delta }\frac{\omega (|\xi |,|\eta |)}{|\xi {|}^2}}}\left[u(x+\xi )-u(x)\right]d\eta d\xi +b(x)=\mathrm{0,}\forall x\in \Omega .$ (5)

Following Aguiar (2016), we introduce the multiplicative decomposition

$\omega (|\xi |,|\eta |)=\tilde{\omega }(|\xi |)|\xi {|}^2\tilde{\omega }(\left|\eta \right|)|\eta {|}^2$ (6)

and define the weighted volume

$m\triangleq {\displaystyle {\int }_{-\delta }^{\delta }\tilde{\omega }}(|\zeta |)|\zeta {|}^{2}d\zeta .$ (7)

Substituting (6) into (5) and using (7), we obtain

${\displaystyle {\int }_{-\delta }^{\delta }C}(\xi )\left[u(x+\xi )-u(x)\right]d\xi +b(x)={\displaystyle {\int }_{{\Omega }^{\prime }}C}(x^{\prime }-x)\left[u(x^{\prime })-u(x)\right]d{x}^{\prime }+b(x)=\mathrm{0,}\forall x\in \Omega ,$ (8)

where

$C(\xi )=\{\begin{array}{c}2c_{1}m\tilde{\omega }(|\xi |)\text{if}\left|\xi \right|\le \delta ,\\ \mathrm{0,}\text{otherwise},\end{array}$ (9)

is the micromodulus function, which represents the spring stiffness of the material, and ${\Omega }^{\prime }\triangleq [a-\delta ,b+\delta ]$ is an extension of the domain $\Omega =[a,b]$ introduced above, as illustrated in Figure 1.

Clearly, if $a\to -\infty$ and $b\to \infty$, then $\Omega \text{'}=\Omega$ and the equilibrium equation (8) has the same form of the one-dimensional equilibrium equation presented by Silling et al. (2003) for an infinite bar made of a linear elastic material with a micromodulus $C_s(\xi )$. There, however, $C_s(\xi )$ may not vanish for $\left|\xi \right|>\delta$. Also, the expression (8) is the governing equation of a bond-based peridynamic material, since it does not take into account the collective deformation of points in a $\delta$-neighborhood of a point $x\in \Omega$.

Figure 1
Domain $\Omega$ and its extension ${\Omega }^{\prime }$.

Even though the concept of stress is not essential in peridynamics, we define below a normal stress function $\sigma (x)$, which is used to relate the micromodulus function $C(\xi )$, given by (9), to the Young’s modulus $E$ in the classical linear elasticity. For this, consider a perpendicular cross section of the cylindrical bar at an arbitrary point $x\in \Omega$ and let $A_1$ and $A_2$be, respectively, the left and the right cross sections, as illustrated in Figure 2. In addition, let $\Omega \text{'}={\Omega }_1\cup {\Omega }_2,$such that

${\Omega }_1=\left\{x^{(1)}\in {\Omega }^{\prime }|x^{(1)}<x\right\},{\Omega }_2=\left\{x^{(2)}\in {\Omega }^{\prime }|x^{(2)}>x\right\},$ (10)

and let ${\mathcal{B}}_i$, $i=\mathrm{1,2}$, be the part of the bar corresponding to ${\Omega }_i$.

Figure 2
The domains ${\Omega }_1$ and ${\Omega }_2$ that arise from the imaginary cut at $x$.

Recall from classical continuum mechanics that the normal stress at the surface $A_1$ is the normal force per unit area exerted across the cross section upon the material on the surface $A_1$ by the material on the surface $A_2$. Since peridynamics is a non-local theory, not only points in $A_2$ but also points in ${\mathcal{B}}_2$ exert forces in ${\mathcal{B}}_1$.

In the one-dimensional theory, we take the normal stress $\sigma (x)$ to be the sum of all resultant forces of the form $F(x^{(2)},x^{(1)})$ that a point $x^{(2)}\in {\Omega }_2$ exerts on a point $x^{(1)}\in {\Omega }_1$. Here,

$F(x^{(2)},x^{(1)})\triangleq \widehat{\underset{\_}{T}}(x^{(1)})\left[\underset{\_}{h}\right]\langle x^{(2)}-x^{(1)}\rangle -\widehat{\underset{\_}{T}}(x^{(2)})\left[\underset{\_}{h}\right]\langle x^{(1)}-x^{(2)}\rangle ,\forall x^{(1)}\in {\Omega }_1,\forall x^{(2)}\in {\Omega }_2.$ (11)

It then follows from (4), (6), (7), (9), and (11) that

$F(x^{(2)},x^{(1)})=C(x^{(2)}-x^{(1)})\left[u(x^{(2)})-u(x^{(1)})\right],\forall x^{(1)}\in {\Omega }_1,\forall x^{(2)}\in {\Omega }_2.$ (12)

The normal stress function of a peridynamic bar is defined by

$\sigma (x)\triangleq {\displaystyle {\int }_{{\Omega }_1}{\displaystyle {\int }_{{\Omega }_2}F}}(x^{(2)},x^{(1)})d{x}^{(2)}d{x}^{(1)}.$ (13)

Taking $x^{(2)}=x+r$ for $r\in [\mathrm{0,}b+\delta -x]$ and $x^{(1)}=x-s$ for $s\in [\mathrm{0,}x-(a-\delta )]$, recalling the definitions in (10), and substituting (12) into (13), we obtain

$\sigma (x)={\displaystyle {\int }_0^{x-a+\delta }{\displaystyle {\int }_0^{b+\delta -x}C}}(r+s)\left[u(x+r)-u(x-s)\right]drds.$ (14)

We now consider a homogeneous deformation, for which $u(x)=\epsilon x$, and rewrite (14) as

$\sigma (x)=\epsilon {\displaystyle {\int }_0^{{\gamma }_1}{\displaystyle {\int }_0^{{\gamma }_2}C}}(r+s)\left(r+s\right)drds,$ (15)

where ${\gamma }_1\triangleq x-(a-\delta )$and ${\gamma }_2\triangleq b+\delta -x$. We present computational results in Section 4 that justify the use of homogeneous deformation in (14). By making the change of variable $v(r)=r+s$ for a fixed $s$ in the inner integral of (15), we rewrite this expression as

$\sigma (x)=\epsilon {\displaystyle {\int }_0^{{\gamma }_1}{\displaystyle {\int }_s^{{\gamma }_2+s}C}}(v)vdvds.$ (16)

Observe from (16) that the domain of integration is given by

$\mathcal{I}=\{(v,s)\in {ℝ}^2|(v,s)\in \left[s,s+{\gamma }_2\right]\times \left[\mathrm{0,}{\gamma }_1\right]\},$

which can be written as

$\mathcal{I}={\mathcal{I}}_1\cup {\mathcal{I}}_2\cup {\mathcal{I}}_3,$

with ${\mathcal{I}}_i$, $i=\mathrm{1,2,3,}$ being a subdomain falling into one of the two cases depicted in Figure 3, which are given by $a){\gamma }_1\le {\gamma }_2$ and $b){\gamma }_1>{\gamma }_2$. We then have from (16) that

$\sigma (x)=\epsilon {\displaystyle \sum _{i=1}^3{\displaystyle {\int }_{{\mathcal{I}}_i}C}(v)vdA}.$ (17)

In case $a)$ of Figure 3, ${\gamma }_1\le {\gamma }_2$, and

$\begin{array}{l}{\mathcal{I}}_1=\left\{(v,s)\in {ℝ}^2|(v,s)\in \left[s,{\gamma }_1\right]\times \left[\mathrm{0,}{\gamma }_1\right]\right\},\\ \\ {\mathcal{I}}_2=\left\{(v,s)\in {ℝ}^2|(v,s)\in \left[{\gamma }_1,{\gamma }_2\right]\times \left[\mathrm{0,}{\gamma }_1\right]\right\},\\ \\ {\mathcal{I}}_3=\left\{(v,s)\in {ℝ}^2|(v,s)\in \left[{\gamma }_2,{\gamma }_2+s\right]\times \left[\mathrm{0,}{\gamma }_1\right]\right\}.\end{array}$

Since $a<x<b$, ${\gamma }_1\ge \delta$ and, because $C(\xi )=0$ for $\left|\xi \right|>\left|\delta \right|$ from (9), the last two integrals in (17) vanish. Taking $dA=vdv$ in the first integral, we arrive at

$\sigma =\epsilon {\displaystyle {\int }_0^{\delta }{\xi }^2}C(\xi )d\xi ,\forall x\in [a,b].$ (18)

In case $b)$ of Figure 3, ${\gamma }_1>{\gamma }_2$, and

$\begin{array}{l}{\mathcal{I}}_1=\left\{(v,s)\in {ℝ}^2|(v,s)\in \left[s,{\gamma }_2\right]\times \left[\mathrm{0,}{\gamma }_2\right]\right\},\\ \\ {\mathcal{I}}_2=\left\{(v,s)\in {ℝ}^2|(v,s)\in \left[{\gamma }_2,{\gamma }_1\right]\times \left[v-{\gamma }_2,v\right]\right\},\\ \\ {\mathcal{I}}_3=\left\{(v,s)\in {ℝ}^2|(v,s)\in \left[{\gamma }_1,{\gamma }_2+s\right]\times \left[{\gamma }_1-{\gamma }_2,{\gamma }_1\right]\right\}.\end{array}$

Again, since $a<x<b$, ${\gamma }_2\ge \delta$ and, because $C(\xi )=0$ for $\left|\xi \right|>\left|\delta \right|$, we arrive at (18).

Figure 3
Domains of integration for the stress function $\sigma (x)$ considering $a){\gamma }_1\le {\gamma }_2$ and $b){\gamma }_1>{\gamma }_2$.

On the other hand, $\sigma (x)=E\epsilon$, $\forall x\in [a,b]$, from classical elasticity, where $E$ is the Young's modulus. It then follows from (18) that

$E={\displaystyle {\int }_0^{\delta }{\xi }^2}C(\xi )d\xi ,$ (19)

which establishes a relation between the Young’s modulus from classical theory and the micromodulus function from peridynamics.

Silling et al. (2003) have obtained an analogous relation between an even micromodulus function, denoted here by $C_S(\xi ),$ and the Young's modulus of an infinite bar. Their relation is given by

$E={\displaystyle {\int }_0^{\infty }{\xi }^2}{C}_S(\xi )d\xi .$ (20)

If $C_S(\xi )$ has the form given by (9), then both expressions (19) and (20) are equivalent in the case of an infinite bar.

Substituting (9) into (19) and using (7), we find that $c_1=E/m^2$, which can be inserted back into (9) to yield

$C(\xi )=\{\begin{array}{ll}\frac{2E\tilde{\omega }(|\xi |)}{m}\hfill & \text{if}\left|\xi \right|\le \delta ,\hfill \\ \mathrm{0,}\hfill & \text{otherwise}.\hfill \end{array}$ (21)

3 FINITE BAR PULLED AT THE ENDS

Now, let us consider a linearly elastic bar of length $L$ in its natural state and take $a=0$, $b=L$, so that $\Omega =(\mathrm{0,}L)$ and ${\Omega }^{\prime }=[-\delta ,L+\delta ]$. The bar is in equilibrium with no body force, so that the equation (8) becomes

${\displaystyle \underset{-\delta }{\overset{L+\delta }{\int }}C}(x^{\prime }-x)[u(x^{\prime })-u(x)]\text{d}{x}^{\prime }=0,\forall x\in (\mathrm{0,}L).$ (22)

On the extended parts of the bar, we impose displacement boundary conditions given by

$u(x^{\prime })=\overline{u}(x^{\prime })for{x}^{\prime }\in [-\delta \mathrm{,0}]\cup \left[\text{L},\text{L}+\delta \right].$ (23)

We then have that the problem of equilibrium without body force of the peridynamic bar consists of finding the axial displacement $u:(\mathrm{0,}L)\to ℝ$ that satisfies the equilibrium equation (22) together with the boundary conditions given by (23).

Now, let

$C_0={\displaystyle \underset{-\delta }{\overset{\delta }{\int }}C}(\xi )\text{d}\xi$ (24)

and observe from (22) that

$\begin{array}{l}0={\displaystyle \underset{-\delta }{\overset{L+\delta }{\int }}C}(x^{\prime }-x)u(x^{\prime })\text{d}{x}^{\prime }-C_{0}u(x)\\ ={\displaystyle \underset{0}{\overset{L}{\int }}C}(x^{\prime }-x)u(x^{\prime })\text{d}{x}^{\prime }-C_{0}u(x)+{\displaystyle \underset{-\delta }{\overset{0}{\int }}C}(x^{\prime }-x)\overline{u}(x^{\prime })\text{d}{x}^{\prime }+{\displaystyle \underset{L}{\overset{L+\delta }{\int }}C}(x^{\prime }-x)\overline{u}(x^{\prime })\text{d}{x}^{\prime }.\end{array}$

This can be rewritten in the form

$u(x)-\frac{1}{C_0}{\displaystyle \underset{0}{\overset{L}{\int }}C}(x^{\prime }-x)u(x^{\prime })\text{d}{x}^{\prime }=\frac{g(x)}{C_0},\forall x\in (\mathrm{0,}L),$ (25)

where

$g(x)={\displaystyle \underset{-\delta }{\overset{0}{\int }}C}(x^{\prime }-x)\overline{u}(x^{\prime })\text{d}{x}^{\prime }+{\displaystyle \underset{L}{\overset{L+\delta }{\int }}C}(x^{\prime }-x)\overline{u}(x^{\prime })\text{d}{x}^{\prime }.$ (26)

The integral equation (25) together with (26) yield an inhomogeneous Fredholm equation of the second kind, where the role of the kernel is played by the micromodulus function and the function $g(x)$ depends only on the boundary data. Conditions for existence of solutions are associated with the form of the kernel (see, for instance, Porter and Stirling (1990)). Granted that a solution exists, the equation can be easily solved numerically. Recall from Section 1 that a Fourier transform technique is used to find a closed form solution for the case of infinite domains. In the next section we introduce different expressions for the micromodulus function and solve numerically the problem defined by the governing equation (22) together with the boundary conditions (23).

4 NUMERICAL RESULTS

We consider four different expressions for the influence function $\tilde{\omega }(|\xi |)$ that appears in (21) together with (7). They are listed in Table 1 together with the corresponding non-zero parts of the micromoduli $C(\xi )$.

Table 1
Expressions for the micromodulus $C(\xi )$, where $\varphi =(-2+e\sqrt{\pi }Erf(1))/3\cong 1.7593$ ¹.

Micromodulus	$\tilde{\omega }(|\xi |)$	$({\delta }^3/3E)C(\xi )$	Ref.
Constant	1	1	Silling et al. (2003)
Exponential	$e^{-{(\xi /\delta )}^2}$	$\varphi e^{-{(\xi /\delta )}^2}$
Triangular	$(1-|\xi |/\delta )$	$4(1-|\xi |/\delta )$	Chen et al. (2016)
Singular	$(\delta /|\xi |-1)$	$2(\delta /|\xi |-1)$	Bobaru et al. (2009)

Observe from Table 1 that we have also included the reference where, to the best of our knowledge, the corresponding expression first appeared. Concerning the exponential micromodulus, Mikata (2012) has introduced an analogous expression for a bar of infinite length, where $\varphi =4/(3\sqrt{\pi })$. Exponential functions are also used in the work of Kilic and Madenci (2010). The constant and the triangular micromoduli correspond to, respectively, the constant and the triangular micromoduli with $n=0$ in Chen et al. (2016) and the singular micromodulus corresponds to the triangular micromodulus with $n=1$ also in Chen et al. (2016).

In all cases, the micromoduli are proportional to $E$ in the interval $(-\delta ,\delta )$ and zero elsewhere. It follows from either (22) or (25) together with both (24) and (26) that the numerical results do not depend on the value of $E$ and, therefore, this value will be ignored in the following computations.

The particularly simple form of the expression (22) allows a straightforward implementation of a numerical scheme, which can be extended to higher dimensions, for the numerical calculation of the displacement field of the peridynamic bar. An alternative numerical scheme², based on the Fredholm equation (25) together with both (24) and (26), has also been implemented, but will not be presented here. Considering the boundary conditions

$\overline{u}(x)=\{\begin{array}{l}0\text{if}x\in (-\delta \mathrm{,0}),\hfill \\ \Delta \text{if}x\in (L,L+\delta ),\hfill \end{array}$ (27)

we first introduce the new variables $x\to xL$ and $u\to u/\Delta$, where, now, $x\in (\mathrm{0,1})$ and $u(x)\in (\mathrm{0,1})$, into (22) and rewrite this expression as

${\displaystyle \underset{-l_0}{\overset{1+l_0}{\int }}C}((x^{\prime }-x)L)[u(x^{\prime })-u(x)]\text{d}{x}^{\prime }=0,\forall x\in [\mathrm{0,1}],$ (28)

where we have introduced the definition $\delta \to \delta L$.

In the simplest approach, one can discretize the interval $(-\delta \mathrm{,1}+\delta )$ with $N+1+2n$ equidistant nodes such that $x_{-n}=-\delta$, $x_0=0$, $x_N=1$, and $x_{N+n}=1+\delta$. We then have that the distance between nodes is constant and given by $h=1/N=\delta /n$ and that the discretized form of the equilibrium equation (28) becomes the system of $N+1$ equations

${\displaystyle \sum _{j=-n}^{N+n}C}((x_i-x_j)L)[u(x_i)-u(x_j)]h=0,i=\mathrm{0,...,}N.$ (29)

In this work we consider geometric sequences of nodes starting with both $N=1000$ and $n=20$ and having a factor of $2$. Substituting $C(\xi )$, given by one of the expressions in Table 1, into the system of equations (29), we can solve for $u(x_i),i=\mathrm{0,}\dots ,N,$ and obtain an approximate expression for the displacement field $u(x)$.

To investigate the influence of the horizon $\delta$ on the solution of the integral equation (28), we consider the constant micromodulus given by the first expression in Table 1 and the four horizons $\delta =1/2,1/5,1/10,1/50$. In Figure 4 we show the displacement $u(x)$ plotted against the position $x\in (\mathrm{0,1})$ for these four horizons and $N=16000$ nodes, which corresponds to the most refined mesh in this work. Because of the relation $1/N=\delta /n$ introduced above, we have that $n=320$ nodes. We see from this figure that (i) the end points of the graphs tend to the points $(\mathrm{0,0})$ on the left end and $(\mathrm{1,1})$on the right end of the bar as $\delta$ decreases; (ii) the solution is highly nonlinear in boundary layers of length $\delta$ near the ends and almost linear outside these boundary layers. In fact, the sequence of solutions tends to the solution of the corresponding problem in the classical linear elasticity theory as $\delta \to 0$, which is linear and satisfies the boundary conditions exactly. This classical solution is represented by the thin solid line in the figure. Thus, at least outside of the boundary layers, we are justified to use homogeneous deformation in (14) to find a relation between the Young’s modulus and the micromodulus function.

We see from these observations that displacement discontinuities appear at the ends of the bar for a fixed $\delta >0$. The jumps are finite and decrease as the horizon decreases. Displacement discontinuities have already been observed for infinite peridynamic bars under concentrated loads (see, for instance, Silling et al. (2003), Weckner and Abeyaratne (2005), Mikata (2012)), but have not been reported yet for finite bars subjected to displacement conditions at the ends.

The reason for the formation of the discontinuities is clear from the analysis of the field equation (28), which expresses, for each particle $x\in (\mathrm{0,1})$ and the corresponding set of particles in its $\delta$-neighborhood, the equilibrium of the forces exerted by the particles on the left- and right-hand sides of $x$. Thus, for a particle located at $x=0$, the particles on the interval $(-\delta \mathrm{,0})$ cannot move and the resultant of the forces on the left hand side must be in equilibrium with the forces on the right-hand side, producing a strain localization. Then, as $x\to 0^{+}$, we have from the left-hand side that ${\displaystyle {\int }_{-\delta }^{0}C}((x^{\prime }-x)L)[u(x^{\prime })-u(x)]\text{d}{x}^{\prime }=-[3E/{(\delta L)}^3]u(x)\delta$ and from the right-hand side that ${\displaystyle {\int }_0^{\delta }C}(x^{\prime }-x)[u(x^{\prime })-u(x)]\text{d}{x}^{\prime }=[3E/{(\delta L)}^3]({\displaystyle {\int }_0^{\delta }u}(x^{\prime })\text{d}{x}^{\prime }-u(x)\delta )$. We then have that $u(0^{+})={\displaystyle {\int }_0^{\delta }u}(x^{\prime })\text{d}{x}^{\prime }/(2\delta )$, that is, the limit value of $u(x)$ at the origin from the right-hand side is half of the average of the displacement field evaluated on the right-hand side of the origin. An analogous argument holds for $x\to 1$.

Figure 4
Displacement $u(x)$ versus position $x\in (\mathrm{0,1})$ obtained from classical linear elasticity and from peridynamics for the constant micromodulus, $N=16000$, and four values of the horizon $\delta$.

We now show that the field equation (28) converges to the second order ordinary differential equation of classical linear elasticity as $\delta \to 0$ for the case of the constant micromodulus. Such demonstration is usually done by using linear displacements as test functions (see, for instance, Silling et al. (2003), Mikata (2012)). Here, however, it is done directly. For $x\in (\mathrm{0,1})$, straightforward applications of the l'Hôpital rule yield

$\begin{array}{l}0=L^3\underset{\delta \to 0}{\lim }{\displaystyle \underset{-\delta }{\overset{1+\delta }{\int }}C}((x^{\prime }-x)L)[u(x^{\prime })-u(x)]\text{d}{x}^{\prime }=\underset{\delta \to 0}{\lim }\frac{3E}{{\delta }_{}^3}({\displaystyle \underset{x-\delta }{\overset{x+\delta }{\int }}u}(x^{\prime })\text{d}{x}^{\prime }-2\delta u(x))\\ =\underset{\delta \to 0}{\lim }\frac{E}{{\delta }^2}([u(x+\delta )-u(x-\delta )]-2u(x))=\underset{\delta \to 0}{\lim }\frac{E}{2\delta }[u^{\prime }(x+\delta )-u^{\prime }(x-\delta )]\\ =\underset{\delta \to 0}{\lim }\frac{E}{2}[u^{″}(x+\delta )+u^{″}(x-\delta )]=E{u}^{″}(x).\end{array}$ (30)

Next, we consider the four micromoduli given in Table 1, the horizon $\delta =1/50$ and, again, $N=16000$ nodes. We then show in Figure 5 the displacement $u(x)$ plotted against the reference position $x\in (\mathrm{0,0.006})$. For comparison purposes, we also show a graph for the displacement $u(x)$ obtained from the exact solution of the bar problem in the context of the classical linear elasticity theory. Observe from this figure that all the curves obtained from peridynamics are above the curve obtained from the classical theory, are almost parallel to each other away from the origin, and tend to different values as we approach the origin from its right-hand side. Following arguments that are similar to those presented above in the analysis of the constant micromodulus case, we have verified that these values are approximations of the expression

$u(0^{+})={\displaystyle \underset{0}{\overset{\delta }{\int }}C}((x^{\prime }-x)L)u(x^{\prime })\text{d}{x}^{\prime }/\left(2{\displaystyle {\int }_0^{\delta }C}((x^{\prime }-x)L)\text{d}{x}^{\prime }\right)$

In the case of the singular micromodulus, we can rewrite this expression in the form

$u(0^{+})={\displaystyle \underset{0}{\overset{\delta }{\int }}C}((x^{\prime }-x)L)[u(x^{\prime })-u(0^{+})]\text{d}{x}^{\prime }/{\displaystyle {\int }_0^{\delta }C}((x^{\prime }-x)L)\text{d}{x}^{\prime },$

and, provided that $u(0^{+})$ is finite, observe that the integrand in the numerator is less singular than the integrand in the denominator, which tends to infinity as $x\to 0^{+}$. Thus, we have that $u(x)\to 0$ as $x\to 0^{+}$. Also, at least outside of the boundary layers, we are again justified to use homogeneous deformation in (14).

From now on, we consider the singular micromodulus only, which is given by the last expression in Table 1. To investigate convergence of the numerical scheme, we show in Figure 6 the displacement $u(x)$ plotted against the reference position $x\in (\mathrm{0,0.003})$ for increasing numbers of nodes and for the horizon $\delta =1/50$. Again, for comparison purposes, we also show the curve obtained from the exact solution of the bar problem in the context of the classical linear elasticity theory. We see from this figure that all the curves obtained from peridynamics are nearly parallel to the curve obtained from the classical linear theory away from the origin, the displacement fields obtained from peridynamics converge to a limit field as the number of nodes increases geometrically, and this limit field seems to have unbounded derivative at $x=0$. Convergence to a limit field that is not the classical solution has been reported in the literature. See, for instance, Bobaru et al. (2009) and Chen et al. (2016).

To verify the unbounded derivative of the limit function, we have considered a node at $x>0$ and its immediate neighbor on its right-hand side and evaluated the relative position $\Delta x$ between the two nodes as well as their relative displacement $\Delta u$ for the fixed value $\delta =1/50$ and for each $N$ shown in Figure 6. We have then obtained Figure 7 showing graphs of the ratio $\Delta u/\Delta x$ versus $x\in (\mathrm{0,0.003})$ for increasing values of $N$. A diamond symbol on the curves yield the position of the first node to the right of the origin corresponding to a discretization with the number of nodes shown next to the symbol. We then see from the figure that all the curves are on the top of each other for $x\ge 0.001$, which corresponds to the position of the first node in a discretization with $1000$ nodes. We also see that the ratio $\Delta u/\Delta x$ becomes unbounded as the mesh becomes refined, with the curves corresponding to less refined meshes being indistinguishable from the curve corresponding to a more refined mesh.

Figure 5
Displacement $u(x)$ versus position $x\in (\mathrm{0,0.006})$ obtained from classical linear elasticity and from peridynamics for $\delta =1/50$, $N=16000$, and four different micromoduli $C$.

Figure 6
Displacement $u(x)$ versus position $x\in (\mathrm{0,0.003})$ obtained from classical linear elasticity and from peridynamics for the singular micromodulus, $\delta =1/50$, and an increasing number of nodes $N$.

In Figure 8 we hold $N=16000$ fixed and show the displacement $u(x)$ plotted against the reference position $x\in (\mathrm{0,1})$ for increasing values of the horizon $\delta$. Again, for comparison purposes, we also show the curve obtained from the classical linear elasticity theory. Observe from this figure that the peridynamic solutions tend to the classical solution as $\delta$ decreases. To investigate the behavior of these curves near the origin, we have plotted Figure 9, which shows $u(x)$ against $x$ in the interval $(\mathrm{0,0.006})$. Observe from this figure that $u(x)$ obtained from peridynamics is highly nonlinear near the origin, tends to $u(x)$ obtained from classical theory as $\delta \to 0$, and, as observed earlier, seems to have unbounded derivative at $x=0$.

Figure 7
Ratio between relative displacement $\Delta u$ and relative position $\Delta x$ versus position $x\in (\mathrm{0,0.003})$ obtained for the case of singular micromodulus, $\delta =1/50$, and an increasing number of nodes $N$.

Figure 8
Displacement $u(x)$ versus position $x\in (\mathrm{0,1})$ obtained from classical linear elasticity and from peridynamics using $N=16000$ and decreasing values of the horizon $\delta$.

5 DISCUSSION AND CONCLUSIONS

We have used the state-based peridynamic theory to formulate the one-dimensional problem of motion of a linear elastic peridynamic bar, which may, or, may not be finite. The governing equation contains a micromodulus function $C(\xi )$, which represents the stiffness of the material at a point $x$ and vanishes outside a $\delta$- neighborhood of $x$. This feature together with homogeneous deformation are used to obtain a relation between $C(\xi )$ and the Young’s modulus $E$ of classical linear elasticity. Analogous relation is available in the literature for an infinite bar only.

Figure 9
Displacement $u(x)$ versus position $x\in (\mathrm{0,0.006})$ obtained from classical linear elasticity and from peridynamics using $N=16000$ and decreasing values of the horizon $\delta$.

We have then investigated the one-dimensional problem of a linear elastic peridynamic bar of finite length $L$ in equilibrium with no body force being pulled at the ends by imposed displacements on extended parts of the bar, which have length $\delta$. We have shown that the problem can be formulated in terms of an inhomogeneous Fredholm equation of the second kind. We have then used a numerical scheme to obtain approximate solutions for the peridynamic problem using different expressions of micromodulus, horizons, and discretizations. Except for the singular micromodulus, all the other micromoduli yield discontinuous displacements at the ends of the bar. Discontinuous displacements have been reported in the study of infinite bars, but not of finite bars. The singular micromodulus gives continuous displacements with unbounded derivatives at these ends. In addition, in all cases we have observed homogeneous deformations outside a boundary layer near the ends of the bar. This observation justifies our approach of using homogeneous deformation in (14) to find a relation between the Young’s modulus and the micromodulus function.

We have then focused our attention in the investigation of the singular case and found that the approximate solutions of the peridynamic problem converge to the classical solution of linear elasticity for vanishing horizon and a fixed discretization. Also, for a fixed horizon, these solutions converge to a limit function as the grid spacing, $h$, tends to zero. The information obtained from the investigation of the one-dimensional peridynamic problem can now be applied to better understand the behavior of solutions in higher dimensions and to design effective numerical schemes for higher-dimensional problems.

6 ACKNOWLEDGEMENTS

The first and fourth authors acknowledge the support of CNPq, under grant 444896/2014-7. The second author acknowledges the support of FAPESP, under grant 2016/12529-4. The third author acknowledges the support of FAPESP, under grant 2016/12217-2, and the partial support of the Italian Ministry, under grant MIUR-PRIN voce COAN 5.50.16.01 code 2015JW9NJT.

References

• Aguiar, A.R. (2016). On the Determination of a Peridynamic Constant in a Linear Constitutive Model. Journal of Elasticity, 122(1):27-39. Erratum to: On the Determination of a Peridynamic Constant in a Linear Constitutive Model. Journal of Elasticity, 122(1):41-42.
• Aguiar, A.R., Fosdick, R. (2014). A constitutive model for a linearly elastic peridynamic body. Mathematics and Mechanics of Solids, 19:502-523.
• Bobaru, F., Yang, M., Alves, L., Silling, S., Askari, E., Xu, J. (2009). Convergence, adaptive refinement, and scaling in 1D peridynamics. International Journal for Numerical Methods in Engineering, 77(6):852-877.
• Chen, Z., Bakenhus, D., Bobaru, F. (2016). A constructive peridynamic kernel for elasticity. Comput. Methods Appl. Mech. Engrg. 311:356-373.
• Kilic, B., Madenci, E. (2010). An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory. Theoretical and Applied Fracture Mechanics 53:194-204.
• Mikata, Y. (2012). Analytical solutions of peristatic and peridynamic problems for a 1D infinite rod. International Journal of Solids and Structures, 49(21):2887-2897.
• Nishawala, V., Ostoja-Starzewski, M. (2017). Peristatic solutions for finite one- and two-dimensional systems. Mathematics and Mechanics of Solids, 22(8):1639-1653.
• Porter, D., Stirling, D. (1990). Integral Equations: A Practical Treatment, from Spectral Theory to Applications. Vol. 5. Cambridge University Press, Cambridge, UK.
• Queiruga, A.F., Moridis, G. (2017). Numerical experiments on the convergence properties of state-based peridynamic laws and influence functions in two-dimensional problems. Comput. Methods Appl. Mech. Engrg. 322:97-122.
• Seleson P., Parks M.L., (2011). On the role of influence function in the peridynamic theory. International Journal for Multiscale Computational Engineering, 9(6):689-706.
• Silling, S., Epton, M., Weckner, O., Xu, J., Askari, E. (2007). Peridynamic states and constitutive modeling. Journal of Elasticity, 88(2):151-184.
• Silling, S., Zimmermann, M., Abeyaratne, R. (2003). Deformation of a peridynamic bar. Journal of Elasticity, 73(1):173--190.
• Wang, L., Xu, J., Wang, J. (2017a). Static and dynamic Green's functions in peridynamics. Journal of Elasticity, 126(1):95--125.
• Wang, L., Xu, J., Wang, J. (2017b). Erratum to: Static and dynamic Green's functions in peridynamics. Journal of Elasticity, 126(1):127--128.
• Weckner, O., Abeyaratne, R. (2005). The effect of long-range forces on the dynamics of a bar. Journal of the Mechanics and Physics of Solids, 53(3):705-728.
• Weckner, O., Brunk, G., Epton, M. (2009). Green's functions in non-local three-dimensional linear elasticity. Proceedings of the Royal Society A, 465(2111):3463-3487.

• Available online: March 19, 2018

• 1
  $Erf(x)=\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-x}^x{e}^{-t^2}}dt$ is the Gauss error function.

• 2
  The scheme is based on a computational code available at http://mathematica.stackexchange.com/questions/104667/solving-fredholm-equation-of-the-2nd-kind, which is written in the program Mathematica®.

Publication Dates

• Publication in this collection
  2018

History

• Received
  31 July 2017
• Reviewed
  26 Nov 2017
• Accepted
  14 Dec 2017