1 Introduction

Energy dissipative behaviour has been observed in many materials. A large class of these behaviours falls under a subcategory of inelastic material behaviour called viscoelastic response. A viscoelastic response is characterised by a combination of elastic (time independent) and viscous (time dependent) responses. Characterisation of viscoelastic material response in the framework of Classical Continuum Mechanics (CCM) has received and continue to receive extensive research interest within the theoretical and computational mechanics community. However, certain mathematical and physics-based limitations associated with the classical theory motivated the development of numerous continuum models [1,2,3,4] to overcome the limitations of the CCM. One of these models is Peridynamics (PD) which was proposed in [5] as a nonlocal reformulation of the CCM theory. Although itself, a continuum theory, PD circumvent the mathematical limitations of the classical theory by replacing the spatial derivatives present in the equation of motion of the CCM theory with integral operators. As a consequence, PD acquired the advantage of not requiring a continuous differentiable displacement field. Therefore, evolving discontinuity, such as the initiation and propagation of cracks, is treated as an inherent characteristic of the material and is thus modelled using the PD fundamental equation of motion without the need for additional assumptions or modifications. This advantage has motivated a growing interest in using PD to characterise materials with strong or evolving discontinuous response [6,7,8,9,10,11,12,13,14,15,16,17,18].

From a physics standpoint, the integral operator necessarily introduces nonlocality to the PD formulation. This nonlocal feature of PD makes it suitable as a modelling paradigm to study the physics of bodies in which the response at a point is influenced by response of material points located at some finite distance away. This has resulted in extensive research efforts to be spent in pushing the scope of application of PD in many frontiers of numerical characterisation of materials response where nonlocal response is envisaged. In a research [19], it was demonstrated that PD can accurately predict size effect in geometrically equivalent structures across a whole range of sizes. PD has also been employed to study other nonlocal physical processes such as metal machining process [20], diffusion and other transport phenomena [21,22,23]. Peridynamic frameworks have also been proposed for homogenisation of heterogeneous materials [24,25,26,27,28,29,30,31] as well as multiscale [32,33,34,35,36,37,38] and Multiphysics modelling [39,40,41,42].

Notwithstanding the growing research interests in exploiting peridynamics for different applications, relatively few research has been undertaken in the direction of developing models and material characterization in the framework of PD viscoelasticity. This is especially needed in emerging fields such as virtual reality simulation that aims to allow realistic visuals and haptic feedbacks for its users in the field of medical training as well as the use of viscoelastic materials as vibration dampers and acoustic absorbers in aerospace, automobile, marine transport, and engineering of nanodevices. The attainment of this goal will require the development of physics-based models of tissues, biomaterials, and polymers. It is known that these materials often exhibit viscoelastic response [43,44,45] when subjected to stimuli.

In addition, viscoelastic materials are known to exhibit complex behaviour such as variation in material properties with change in sample size [46, 47], a phenomena known as size/scale effect. Scale effect is a very important consideration in the field of miniaturisation engineering where viscoelastic nanomaterials are increasingly being used as candidate materials for development of nanodevices [48]. As the scale of observation becomes smaller, intermolecular interaction becomes too important to ignore [49]. This effectively introduces a notion of nonlocal interaction and hence an internal length scale in the analysis. This effectively invalidates the local hypothesis of the CCM [50], and emphasises the need to develop a nonlocal viscoelastic modelling framework that will incorporate size dependent length scale parameter.

Other complex phenomena encountered in the characterization of viscoelastic materials include initiation and propagation of cracks [51] and its dependence on speed of propagation [52] as well as dynamic crack induced pressure and shear waves [53]. These are a few of the motivating reasons that necessitate the development of PD viscoelasticity framework with capability to model continuous, discontinuous, and nonlocal material response. To this end, characterisation of viscoelastic materials in the framework of PD theory is beginning to gain attention. A Bond-based Peridynamic (BBPD) viscoelastic material model was proposed in [54]. To derive the integral representation of the nonlocal viscoelastic model, Fourier- and Laplace-transform were utilized. It was demonstrated through analytical and numerical solution of the viscoelastic model that the proposed formulation can predict viscoelastic response.

A thermodynamically consistent viscoelasticity model in the framework of Ordinary state-based peridynamics (OSBPD) was proposed in [55] that relies on the additive decomposition of the extension state into a volumetric and deviatoric components. The deviatoric extension state is further decomposed into elastic and back extension components. A key assumption in the formulation is that the volumetric response is assumed to be elastic while the deviatoric response is viscoelastic. The viscoelastic response was assumed to be governed by a combination of the standard linear solid and Maxwell model. The proposed framework was implemented to study the relaxation response of a single bond. In [56], an explicit time integration algorithm for this model was proposed. This OSBPD viscoelastic framework was extended in [57] to incorporate Prony series representation of the viscoelastic material property functions. This allows the possibility of fitting the Prony series with viscoelastic material functions such as the relaxation modulus and creep compliance data obtained from experiment. A viscoelastic characterisation of fracture in membranes was undertaken in [58]. A generalization of the OSBPD viscoelastic framework was proposed in [59] to account for viscoelastic response in both spherical and deviatoric response of materials.

A viscoelastic model based on a NOSBPD correspondence formulation called Bond-Associated Peridynamic [60] framework was proposed in [61]. The framework was utilised to simulate viscoelastic response in the interface region of two dissimilar materials. Other rate dependent formulations proposed in the framework of PD include viscoplasticity [62] and visco-hyperelasticity [63, 64].

The aim in this contribution is to propose a viscoelastic framework for Non-ordinary State-based Peridynamics (NOSBPD). The motivation is to leverage on the ability of the NOSBPD formulation to admit state-of-the-art constitutive models of the CCM theory through correspondence principle, thus allowing for viscoelastic modelling of generalised materials within the peridynamic framework.

This proposed formulation differs from the previously proposed method on several fronts. The first point of departure is that the previous formulations [54, 55, 57] were built on the framework of either BBPD or OSBPD. This poses restriction on the range of materials that can be modelled using these frameworks since BBPD can only admit materials with a Poisson’s ration of \({\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 4}}\right.\kern-0pt} \!\lower0.7ex\hbox{$4$}}\), while OSBPD is based on the assumption that the bond force between two interacting material points is parallel to the bond [65].

An advantage to be drawn from this formulation which is based on NOSBPD is that the NOSBPD is a generalisation of the BBPD and the OSBPD thus allowing for far reaching generalization in the scope of materials that can be modelled. Also, a particular formulation of the NOSBPD called correspondence model uses familiar quantities such as stress and strain based on the principle of constitutive correspondence. This allows the proposed framework to feed off from the state-of-the-art viscoelastic constitutive model from the CCM framework. Another point of departure between this proposed framework and the previously developed frameworks is that this proposed framework allows for both partial and total viscoelastic formulation to be implemented.

Although the formulation proposed in [61] share some similarities with this proposed framework in that they are both based on the general framework of NOSBPD correspondence model, the two formulation differs fundamentally in their kinematics. Whilst the method proposed in [61] uses a bond-associated method to compute deformation gradient, in this proposed framework, the nonlocal deformation gradient is computed using the method proposed in the original development of the NOSBPD [66]. When computing the nonlocal deformation gradient in the bond-associated method, only a fraction of the material points within the horizon of a material point are utilised. This can lead to inaccurate prediction of strain energy density hence the need for a correction factor to compensate for this inadequacy [60].

Another feature of this formulation that distinguishes it from the prior formulations is that it permits a more generalised viscoelastic implementation. This formulation allows partial viscoelastic formulation, unlike the prior formulations, where the assumption of bulk elastic response and shear viscoelastic response is imposed, as well as a total viscoelastic formulation where both bulk and deviatoric viscoelastic responses are accounted for.

In the reminder of this presentation, a concise review of the NOSBPD theory will be presented in 0 to be followed by a discussion on nonlocal kinematics relevant to the state-based formulation in 0. Constitutive correspondence principle linking constitutive models from CCM with the NOSBPD is presented in 0. Discussion on a state variable viscoelastic constitutive model is presented in 0 and implementation strategy for the proposed framework is discussed in 0. Results of numerical simulation to validate the proposed framework is presented in 1.0 and concluding remarks in 2.0.

2 State-based peridynamic theory

One of the foundational concepts in PD just as in CCM is the continuum hypothesis, which asserts that a body consists of an infinite number of particles, each of which shares the same physical characteristics as the bulk. A key distinction between the two is that while in CCM, a material point interacts only with its immediate neighbours, PD allows interaction over a finite distance. Consider a deformable body \(\mathcal{B}\) as shown in Fig. 1. A point P occupying the position \(\mathbf{x}\) in \(\mathcal{B}\) interacts with infinitely many other points located within a finite domain. If this domain is taken to be a ball of radius \(\delta >0\) centred at \(\mathbf{x}\), such that

$${\mathcal{B}}_{\delta }\left(\mathbf{x}\right)=\left\{{\mathbf{x}}^{\mathrm{^{\prime}}}\in \mathcal{R}:\left|{\mathbf{x}}^{\mathrm{^{\prime}}}-\mathbf{x}\right|<\delta \right\}$$
(1)

then \({\mathcal{B}}_{\delta }\left(\mathbf{x}\right)\) is called the family of \(\mathbf{x}\) and contains all the material points Q occupying the position \(\mathbf{x}\mathbf{^{\prime}}\) that interacts with P (for the rest of this presentation, points P and Q will be identified by their position vectors \(\mathbf{x}\) and \({\mathbf{x}}^{\mathbf{^{\prime}}},\) respectively). The interaction between \(\mathbf{x}\) and \(\mathbf{x}\mathbf{^{\prime}}\) is referred to as a bond and the distance \({\varvec{\xi}}={\mathbf{x}}^{\mathrm{^{\prime}}}-\mathbf{x}\) in the reference or undeformed configuration is called the bond length. The family of bonds \({\mathcal{H}}_{\mathbf{x}}\) for point \(\mathbf{x}\) is thus given as

Fig. 1
figure 1

A continuum body showing a point and its domain of influence

Full size image
$${\mathcal{H}}_{\mathbf{x}}=\left\{{\varvec{\xi}}\in \left({\mathbb{R}}\backslash 0\right)|\left({\varvec{\xi}}+\mathbf{x}\right)\in \left({\mathcal{B}}_{\delta }\left(\mathbf{x}\right)\cap \mathfrak{B}\right)\right\}$$
(2)

Let \({\varvec{u}}(\mathbf{x})\) and \({\varvec{u}}\left({\mathbf{x}}^{\mathbf{^{\prime}}}\right)\) be the respective displacement undergone by \(\mathbf{x}\) and \(\mathbf{x}\mathbf{^{\prime}}\) in the deformed configuration. Then the relative displacement between \(\mathbf{x}\) and \(\mathbf{x}\mathbf{^{\prime}}\) is given by: \({\varvec{\eta}}={\varvec{u}}\left({\mathbf{x}}^{\mathbf{^{\prime}}}\right)-{\varvec{u}}\left(\mathbf{x}\right)\). Consequently, the positions of \(\mathbf{x}\) and \(\mathbf{x}\mathbf{^{\prime}}\) in the deformed configuration are given as \(\mathbf{y}=\mathbf{x}+\boldsymbol{ }{\varvec{u}}\left(\mathbf{x}\right)\) and \(\mathbf{y}={\mathbf{x}}^{\mathbf{^{\prime}}}+\boldsymbol{ }{\varvec{u}}\left({\mathbf{x}}^{\mathbf{^{\prime}}}\right)\), respectively.

In the context of state-based PD [66], functions that are defined on bonds are called states. There are three (3) important states that are deserving of specific mention here. These are the reference position vector state, the deformation vector state, and the force vector state. For the point \(\mathbf{x}\), the reference position vector state \(\underset{\_}{\mathbf{X}}\) is a function that maps a bond vector in the family of \(\mathbf{x}\) in the undeformed configuration unto itself, that is

$$\underset{\_}{\mathbf{X}}\left[\mathbf{x}\right]\langle {\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x} \rangle ={\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x}={\varvec{\xi}}$$
(3)

The deformation vector state is a function that operates on a bond vector in the undeformed configuration, whose value is the image of the bond vector in the deformed configuration. That is,

$$\begin{aligned} {\mathbf{\underline {Y} }}\left[ {{\mathbf{x}},t} \right]{\mathbf{x^{\prime}}} - {\mathbf{x}}{ } & = {\mathbf{y^{\prime}}}\left( {{\mathbf{x^{\prime}}},t} \right) - {\mathbf{y}}\left( {{\mathbf{x}},t} \right) \\ & = {\varvec{u}}\left( {{\mathbf{x^{\prime}}},t} \right) - {\varvec{u}}\left( {{\mathbf{x}},t} \right) + {\varvec{\xi}} \\ \end{aligned}$$
(4)

A force vector state \(\underset{\_}{\mathbf{T}}\) is a function that maps the deformation vector state to a force density vector for each bond vector in \({\mathcal{H}}_{\mathbf{x}}\) such that:

$$\underset{\_}{\mathbf{T}}\left[\mathbf{x},t\right]\langle {\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x}\rangle =\mathbf{t}\left({\mathbf{x}}^{\mathbf{^{\prime}}},\mathbf{x},t\right),$$
(5)

where \(\mathbf{t}\) is a force density at \(\mathbf{x}\) due to its interaction with \({\mathbf{x}}^{\mathbf{^{\prime}}}\). The equation that tracks the motion of point \(\mathbf{x}\) is given by the following integro-differential equation

$$\rho \frac{{\partial }^{2}{\varvec{u}}\left(\mathbf{x}\right)}{\partial {t}^{2}} ={\int }_{{\mathcal{B}}_{\delta }\left({\varvec{x}}\right)}\left(\underset{\_}{\mathbf{T}}\left[\mathbf{x},{\varvec{t}}\right]\langle {\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x}\rangle -\underset{\_}{\mathbf{T}}\left[{\mathbf{x}}^{\boldsymbol{^{\prime}}},{\varvec{t}}\right]\langle \mathbf{x}-{\mathbf{x}}^{\mathbf{^{\prime}}}\rangle \right)d{\mathbf{x}}^{\mathbf{^{\prime}}}+{\varvec{b}}\left(\mathbf{x},t\right),$$
(6)

where \({\varvec{b}}\left(\mathbf{x},t\right)\) is a body force density and \(t\) denotes time, and the angle brackets, \(\langle \bullet \rangle\) in (3–6) signify the vector on which the state operates.

3 Constitutive correspondence and nonlocal kinematics

As mentioned in Sect. 1 , the NOSBPD correspondence model utilizes the constitutive model of the CCM and thus allowing the familiar concepts of stress and strain tensors to be used in the framework of PD. To be able to leverage this capability, nonlocal kinematic quantities that parallel those in the CCM will be derived in the proceeding passage. Utilising the definition of the nonlocal gradient operator given in [67, 68], a Taylor series expansion of (4) about \(\xi =0\) will yield

$$\underset{\_}{\mathbf{Y}}\left[\mathbf{x},t\right]\langle {\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x} \rangle ={\varvec{\xi}}+{{\mathcal{G}}_{\omega }}_{\mathbf{x}}{\varvec{u}}\left(\mathbf{x}\right)\bullet{\varvec{\xi}}+\frac{1}{2}{\varvec{\xi}}\bullet \left({{\mathcal{G}}_{\omega }}_{\mathbf{x}}^{2}{\varvec{u}}\left(\mathbf{x}\right)\bullet{\varvec{\xi}}\right)+\dots \underset{\_}{\mathbf{Y}}\left[\mathbf{x},t\right]\langle {\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x} \rangle \approx \left({{\mathcal{G}}_{\omega }}_{\mathbf{x}}{\varvec{u}}\left(\mathbf{x}\right)+\mathbf{I}\right)\bullet{\varvec{\xi}}+\mathcal{O}\left({\left|{\varvec{\xi}}\right|}^{2}\right),$$
(7)

where \({{\mathcal{G}}_{\omega }}_{\mathbf{x}}\) is the nonlocal gradient operator evaluated at \(\mathbf{x}\), \(\mathbf{I}\) is the second order isotropic tensor and \(\mathcal{O}\left({\left|{\varvec{\xi}}\right|}^{2}\right)\) represent terms of order higher than one. If small deformation assumption is made, then \(\mathcal{O}\left({\left|{\varvec{\xi}}\right|}^{2}\right)=0\) and (7) reduces to

$$\underset{\_}{\mathbf{Y}}\left[\mathbf{x},t\right]\langle {\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x} \rangle =\left({{\mathcal{G}}_{\omega }}_{\mathbf{x}}{\varvec{u}}\left(\mathbf{x}\right)+\mathbf{I}\right)\bullet{\varvec{\xi}}$$
(8)

The parenthesized term in the right hand side of (8) is the nonlocal deformation gradient [68]. If we denote the nonlocal deformation gradient as \(\mathbf{F}\), then (8) can be written as,

$$\underset{\_}{\mathbf{Y}}\left[\mathbf{x},t\right]\langle {\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x} \rangle =\mathbf{F}{\varvec{\xi}}$$
(9)

The nonlocal deformation gradient is given by [66]

$${\mathcal{G}}_{{\omega {\mathbf{x}}}} \left( {\mathbf{y}} \right) = {\mathbf{F}}\left( {\mathbf{x}} \right) = \left[ {\int\limits_{{{\mathbb{R}} ^{n} }} {\omega \left( {\left| {\varvec{\xi}} \right|} \right)\left( {{\mathbf{y}}\left( {{\mathbf{x}}^{^{\prime}} ,t} \right) - {\mathbf{y}}\left( {{\mathbf{x}},t} \right)} \right) \otimes {\varvec{\xi}}d{\mathbf{x}}^{^{\prime}} } } \right]{\mathbf{K}}^{ - 1} ,$$
(10)

where \(\mathbf{K}\) is a shape tensor defined as

$${\mathbf{K}} = \mathop \smallint \limits_{{{\mathcal{H}}_{{\mathbf{x}}} }} \underline {\omega } {\varvec{\xi}}{\mathbf{\underline {X} }}{\varvec{\xi}} \otimes {\mathbf{\underline {X} }}{\varvec{\xi}}dV_{{\varvec{\xi}}}$$
(11)

A nonlocal strain tensor can be defined by considering the difference between the squares of bond length \(\xi =\left|{\varvec{\xi}}\right|=\left|{\mathbf{x}}^{\mathbf{^{\prime}}}-\mathbf{x}\right|\) in the undeformed configuration and \(\lambda =\left|{\varvec{\lambda}}\right|=\left|\mathbf{y}\left({\mathbf{x}}^{\mathbf{^{\prime}}}\right)-\mathbf{y}\left(\mathbf{x}\right)\right|\) in the deformation configuration as shown in Fig. 2. From (9), we can write

Fig. 2
figure 2

Deformation of a Peridynamic body

Full size image
$${\lambda }^{2}-{\xi }^{2}={\varvec{\lambda}}\bullet{\varvec{\lambda}}-{\varvec{\xi}}\bullet{\varvec{\xi}}={{\varvec{\xi}}}^{T}{\mathbf{F}}^{T}\mathbf{F}{\varvec{\xi}}-{{\varvec{\xi}}}^{T}{\varvec{\xi}}={{\varvec{\xi}}}^{T}\left({\mathbf{F}}^{T}\mathbf{F}-\mathbf{I}\right){\varvec{\xi}}$$
(12)

If a deformation or strain tensor \(\mathbf{E}\) is defined as

$$\mathbf{E}=\frac{1}{2}\left({\mathbf{F}}^{T}\mathbf{F}-\mathbf{I}\right)$$
(13)

then (13) becomes

$${\lambda }^{2}-{\xi }^{2}=2{{\varvec{\xi}}}^{T}\mathbf{E}{\varvec{\xi}}$$
(14)

The tensor \(\mathbf{E}\) in (14) is the nonlocal equivalent of the Green–Lagrange strain tensor and \(\mathbf{F}\) is the nonlocal deformation gradient tensor defined in (10). Under the assumption of infinitesimal deformation, it has been shown in [68] that the nonlocal deformation gradient tensor (13) can be written as

$$\mathbf{E}\approx {\varvec{\upvarepsilon}}=\frac{1}{2}\left(\mathbf{F}+{\mathbf{F}}^{T}\right)-\mathbf{I},$$
(15)

where \({\varvec{\upvarepsilon}}\) is the infinitesimal strain tensor and \(\mathbf{I}\) is the second order isotropic tensor.

4 Correspondence model

For a given deformation tensor \(\mathbf{F}\), let \(\Phi \left(\bullet \right)\) be a function that maps \(\mathbf{F}\) to the corresponding strain energy density in CCM. Also, let \(W\) be the nonlocal PD strain energy density. It has been shown in [66] that constitutive correspondence between PD and CCM is achieved if for some choice of the deformation tensor \(\mathbf{F}\) given by (9) the PD strain energy density is equal to the strain energy density from CCM, that is,

$$\mathrm{W}\left(\underset{\_}{\mathbf{Y}}\right)=\Phi \left(\mathbf{F}\right)$$
(16)

The particular form of the PD force density vector state \(\underset{\_}{\mathbf{T}}\) is obtained by evaluating the Fréchet derivative of the strain energy function \(W\left(\underset{\_}{\mathbf{Y}}\right)\) with respect to \(\underset{\_}{\mathbf{Y}}\). This operation [66] yields the expression for the PD force density vector state as

$$\underset{\_}{\mathbf{T}}\langle {\varvec{\xi}}\rangle =\underset{\_}{\omega }\left(\left|{\varvec{\xi}}\right|\right)\mathbf{P}{\mathbf{K}}^{-1}{\varvec{\xi}},$$
(17)

where \(\underset{\_}{\omega }\) is a weight function and \(\mathbf{P}\) is the first Piola stress tensor. Under the assumption of small deformation, (17) is rewritten as

$$\underset{\_}{\mathbf{T}}\langle {\varvec{\xi}}\rangle =\underset{\_}{\omega }\left(\left|{\varvec{\xi}}\right|\right){\varvec{\sigma}}{\mathbf{K}}^{-1}{\varvec{\xi}},$$
(18)

where \({\varvec{\sigma}}\) is the Cauchy stress tensor.

5 Viscoelastic constitutive model

The correspondence model (17) or (18) allow PD to admit a wide spectrum of constitutive models from the CCM framework. In this section, a linear viscoelastic model will be derived that is compatible with the PD correspondence model. Two equivalent approaches have generally been utilised to develop viscoelastic constitutive models. The first approach uses ordinary differential equations to establish the constitutive relation while the second establishes constitutive relation between stress and strain using integrals via Boltzmann superposition principle. The three basic mechanical models typically used to model linear viscoelastic response are the Maxwell, Kelvin, and the Standard Linear Solid (SLS) models.

In this work, we limit our focus to linear viscoelasticity based on the integral formulation. The integral formulation in one-dimension is often stated as

$$\sigma \left(t\right)={\int }_{0}^{t}E\left(t-\tau \right)\frac{d}{d\tau }\varepsilon \left(\uptau \right)d\tau ,$$
(19)

where \(E\) is the time dependent relaxation modulus and the variable of integration \(\tau\) is a time variable measured backward from the current time \(t\). Several mathematical representations of the relaxation modulus exist. A very popular representation among these is the Prony-Dirichlet series or Prony series for short. The Prony series consist of a finite sum of decaying exponential often expressed [69] as

$$E\left(t\right)={E}_{0}+{\sum }_{j=1}^{n}{E}_{j}\mathrm{exp}\left(-\frac{t}{{\tau }_{i}}\right)$$
(20)

It is worth noting that in (20), the extreme values of the relaxation modulus can easily be obtained as

$$E\left(0\right)={E}_{0}+{\sum }_{j=1}^{n}{E}_{j}, E\left(\infty \right)={E}_{0},$$
(21)

where \(E(0)\) is the instantaneous Young’s modulus and \(E(\infty )\) is the long term or equilibrium Young’s modulus. Introducing (20) into (19)

$$\begin{aligned} \sigma \left( t \right) & = \mathop \smallint \limits_{0}^{t} E_{0} \frac{d}{d\tau }\varepsilon \left( {\uptau } \right)d\tau + \mathop \sum \limits_{j}^{n} \mathop \smallint \limits_{0}^{t} E_{j} \exp \left( { - \frac{t - \tau }{{\tau_{j} }}} \right)\frac{d}{d\tau }\varepsilon \left( {\uptau } \right)d\tau \\ & = E_{0} \varepsilon \left( t \right) + \mathop \sum \limits_{j}^{n} \mathop \smallint \limits_{ - \infty }^{t} E_{j} \exp \left( { - \frac{t - \tau }{{\tau_{j} }}} \right)\frac{d}{d\tau }\varepsilon \left( {\uptau } \right)d\tau \\ & = \sigma_{0} \left( t \right) + \mathop \sum \limits_{j = 1}^{n} h_{j} \left( t \right) , \\ \end{aligned}$$
(22)

where \({\sigma }_{0}\left(t\right)={E}_{0}\varepsilon \left(t\right)\) is the elastic contribution and

$${h}_{j}\left(t\right)={\int }_{0}^{t}{E}_{j}\mathrm{exp}\left(-\frac{t-\tau }{{\tau }_{j}}\right)\frac{d\varepsilon \left(\uptau \right)}{d\tau }d\tau$$
(23)

is referred to as the internal state variable and account for the time dependent contribution to the viscoelastic response. For this reason, this methodology can be considered to belong to a class of methodologies called internal variable methods in the literature. If the strain in the integral of (23) is substituted by \(\varepsilon \left(\uptau \right)={\sigma }_{0}\left(t\right)/{E}_{0}\), we obtain

$${h}_{j}\left(t\right)={\int }_{0}^{t}{\gamma }_{j}\mathrm{exp}\left(-\frac{t-\tau }{{\tau }_{j}}\right)\frac{d{\sigma }_{0}\left(t\right)}{d\tau }d\tau ,$$
(24)

where \({\gamma }_{j}\) is a normalised relaxation function given by \({\gamma }_{j}={E}_{j}/{E}_{0}\). To implement (22) in a numerical scheme such as in the framework of NOSBPD correspondence model will require a strategy for numerical integration of (23) or (24). This can be achieved using an algorithm [69] that discretises for example (24) in terms of a finite time interval. Let the deformation history be split into two periods \(0\le \tau \le {t}_{n}\) of known deformation and \({t}_{n}\le \tau \le {t}_{n+1}\) of unknown deformation, then the integral in (24) can accordingly be split and written as

$$\begin{aligned} h_{j} \left( {t_{n + 1} } \right) & = \mathop \smallint \limits_{0}^{{t_{n} }} \gamma_{j} \exp \left( { - \frac{{t_{n + 1} - \tau }}{{\tau_{j} }}} \right)\frac{d\varepsilon \left( \tau \right)}{{d\tau }}d\tau + \mathop \smallint \limits_{{t_{n} }}^{{t_{n + 1} }} \gamma_{j} \exp \left( { - \frac{{t_{n + 1} - \tau }}{{\tau_{j} }}} \right)\frac{{d\sigma_{0} \left( t \right)}}{d\tau }d\tau \\ & = \mathop \smallint \limits_{0}^{{t_{n} }} \gamma_{j} \exp \left( { - \frac{{t_{n} + \Delta t - \tau }}{{\tau_{j} }}} \right)\frac{{d\varepsilon \left( {\uptau } \right)}}{d\tau }d\tau + \mathop \smallint \limits_{{t_{n} }}^{{t_{n + 1} }} \gamma_{j} \exp \left( { - \frac{{t_{n + 1} - \tau }}{{\tau_{j} }}} \right)\frac{{d\sigma_{0} \left( t \right)}}{d\tau }d\tau \\ & = \mathop \smallint \limits_{0}^{{t_{n} }} \gamma_{j} \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)\exp \left( { - \frac{{t_{n} - \tau }}{{\tau_{j} }}} \right)\frac{{d\sigma_{0} \left( t \right)}}{d\tau }d\tau + \mathop \smallint \limits_{{t_{n} }}^{{t_{n + 1} }} \gamma_{j} \exp \left( { - \frac{{t_{n + 1} - \tau }}{{\tau_{j} }}} \right)\frac{{d\sigma_{0} \left( t \right)}}{d\tau }d\tau \\ & = \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)\mathop \smallint \limits_{0}^{{t_{n} }} \gamma_{j} \exp \left( { - \frac{{t_{n} - \tau }}{{\tau_{j} }}} \right)\frac{{d\sigma_{0} \left( t \right)}}{d\tau }d\tau + \mathop \smallint \limits_{{t_{n} }}^{{t_{n + 1} }} \gamma_{j} \exp \left( { - \frac{{t_{n + 1} - \tau }}{{\tau_{j} }}} \right)\frac{{d\sigma_{0} \left( t \right)}}{d\tau }d\tau \\ & = \exp \left( { - \frac{\Delta t}{{\tau_{i} }}} \right)h_{j} \left( {t_{n} } \right) + \mathop \smallint \limits_{{t_{n} }}^{{t_{n + 1} }} \gamma_{j} \exp \left( { - \frac{{t_{n + 1} - \tau }}{{\tau_{j} }}} \right)\frac{{d\sigma_{0} \left( t \right)}}{d\tau }d\tau \\ \end{aligned}$$
(25)

If we make the simplifying assumption that over the time increment \(\Delta t\), \({\sigma }_{0}\) varies linearly with time, then

$$\frac{{d\sigma_{0} }}{dt} = \frac{{\Delta \sigma_{0} }}{\Delta t} = \frac{{\sigma_{0}^{n + 1} - \sigma_{0}^{n} }}{\Delta t}$$
(26)

Introducing (26) into (25) yields

$$\begin{aligned} h_{j} \left( {t_{n + 1} } \right) & = \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)h_{j} \left( {t_{n} } \right) + \mathop \smallint \limits_{{t_{n} }}^{{t_{n + 1} }} \gamma_{j} \exp \left( { - \frac{{t_{n + 1} - \tau }}{{\tau_{j} }}} \right)\frac{{\sigma_{0}^{n + 1} - \sigma_{0}^{n} }}{\Delta t}d\tau \\ & = \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)h_{j} \left( {t_{n} } \right) + \gamma_{j} \tau_{j} \left( {1 - \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)} \right)\frac{{\sigma_{0}^{n + 1} - \sigma_{0}^{n} }}{\Delta t} \\ & = \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)h_{j} \left( {t_{n} } \right) + \gamma_{j} \frac{{\left( {1 - \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)} \right)}}{{\frac{\Delta t}{{\tau_{j} }}}}\left[ {\sigma_{0}^{n + 1} - \sigma_{0}^{n} } \right] \\ & = \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)h_{j} \left( {t_{n} } \right) + \gamma_{j} \frac{{\left( {1 - \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)} \right)}}{{\frac{\Delta t}{{\tau_{j} }}}}\left[ {\sigma_{0}^{n + 1} - \sigma_{0}^{n} } \right] \\ \end{aligned}$$
(27)

If we write \(A_{j} = \left( {1 - \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)} \right)/\left( {\frac{\Delta t}{{\tau_{j} }}} \right)\), then

$$h_{j} \left( {t_{n + 1} } \right) = \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)h_{j} \left( {t_{n} } \right) + \gamma_{j} A_{j} \left[ {\sigma_{0}^{n + 1} - \sigma_{0}^{n} } \right]$$
(28)

Thus the stress at time \({t}_{n+1}\) is computed as

$$\begin{aligned} \sigma^{n + 1} & = \sigma_{0}^{n + 1} + \mathop \sum \limits_{j = 1}^{n} h_{j}^{n + 1} \\ & = E_{0} \varepsilon^{n + 1} + \mathop \sum \limits_{j = 1}^{n} \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right)h_{j} \left( {t_{n} } \right) + \gamma_{j} A_{j} \left[ {\sigma_{0}^{n + 1} - \sigma_{0}^{n} } \right] \\ \end{aligned}$$
(29)

Observe that (29) is a recursive function that requires the variables \({h}_{j}^{n}\) and \({\varepsilon }^{n}\) from the previous times step to be able to compute the stress state at the current time step. A generalisation of (29) to higher dimension is achieved by rewriting the Prony series (20) in terms of tensor quantities (using Voigt notation) as:

$$C_{ij} \left( t \right) = C_{ij0} + \sum\limits_{m = 1}^{n} {C_{ijm} } \exp \left( { - \frac{t}{{\tau_{m} }}} \right),$$
(30)

where the transient relaxation modulus \({C}_{ij}(t)\) is a fourth order tensor. Consequently, each term in (30) is also a forth order tensor which in the context of viscoelastic materials such as polymers represents a particular mechanism by which the material responds to an applied load. The first term, \({C}_{i{j}_{o}}={C}_{ij}(t=\infty )\) is the equilibrium relaxation modulus which represents the elastic component of the material response. The component of the relaxation tensor \({C}_{i{j}_{m}}\) and the relaxation time parameter \({\tau }_{m}\) in the m\(\mathrm{th}\) term of the series characterises the viscous response. Introducing (30) into (19) would yield the generalised form of (29) as

$${\varvec{\sigma}}^{n + 1} = {\mathbf{\underline {C} }}_{0} {\varvec{\varepsilon}}^{n + 1} + \mathop \sum \limits_{j = 1}^{n} \exp \left( { - \frac{\Delta t}{{\tau_{j} }}} \right){\varvec{h}}_{j}^{n} + \gamma_{j} A_{j} \left[ {{\varvec{\sigma}}_{0}^{n + 1} - {\varvec{\sigma}}_{0}^{n} } \right],$$
(31)

where \({\varvec{\sigma}}\) and \({\varvec{\varepsilon}}\) are the vector representation of the second order stress and strain tensors respectively, and the fourth order tensors \({\underset{\_}{\mathbf{C}}}_{0}\) and \({\underset{\_}{\mathbf{C}}}_{j}\) are as defined in (30).

For practical and numerical consideration, it is sometimes necessary to make the assumption that the material behave elastically in dilatation and viscoelastically in its shear response. To this effect, the stress tensor is split into a spherical and deviatoric components and given by

$${{\varvec{\sigma}}}^{n+1}=\mathbf{I}{\sigma }_{m}^{n+1}+{\widetilde{{\varvec{\sigma}}}}^{n+1},$$
(32)

where \(\mathbf{I}\) is a second order isotropic tensor, \({\sigma }_{m}\) and \(\widetilde{{\varvec{\sigma}}}\) are the mean stress and deviatoric stress, respectively. If the assumption of elastic bulk response is enforced, (32) is shown [69] to take the form

$${{\varvec{\sigma}}}^{n+1}=K{\varepsilon }_{m}^{n+1}\mathbf{I}+{\widetilde{{\varvec{\sigma}}}}^{n+1}$$
(33)

In this case, the procedure for the incremental determination of the state variable is applied only on the deviatoric component \(\widetilde{{\varvec{\sigma}}}\) of the stress tensor such that

$${\widetilde{{\varvec{\sigma}}}}^{n+1}={\widetilde{{\varvec{\sigma}}}}_{0}^{n+1}+\sum_{j=1}^{n}{{\varvec{h}}}_{j}^{n+1},$$
(34)

where

$${{\varvec{h}}}_{j}^{n+1}=\mathrm{exp}\left(-\frac{\Delta t}{{\tau }_{i}}\right){{\varvec{h}}}_{j}^{n}+{\gamma }_{j}{A}_{j}\left[{\widetilde{{\varvec{\sigma}}}}_{0}^{n+1}-{\widetilde{{\varvec{\sigma}}}}_{0}^{n}\right]$$
(35)

6 Implementation strategy

The framework proposed in this contribution is capable of simulating both dynamic and static response of materials undergoing linear viscoelastic deformation. The numerical solution strategy adopted in this formulation consists of a nested spatial and time numerical integration procedures. The spatial integration procedure evaluates the total force acting on a material point at a given time while the time integration procedure tracks the position of material points with passage of time. To numerically implement the spatial integration, the meshfree method proposed in [70] is utilised. In this method, the equation of motion (6) is approximated by

$${\rho }_{i}{\ddot{{\varvec{u}}}}_{{\varvec{i}}}=\sum_{j=1}^{N}\left[\underset{\_}{\mathbf{T}}\left[{\mathbf{x}}_{i},t\right]\langle {\mathbf{x}}_{j}-{\mathbf{x}}_{i}\rangle -\underset{\_}{\mathbf{T}}\left[{\mathbf{x}}_{j},t\right]\langle {\mathbf{x}}_{i}-{\mathbf{x}}_{j}\rangle \right]{V}_{j}+{{\varvec{b}}}_{{\varvec{i}},}$$
(36)

where \(\rho_{i} : = \rho \left( {{\mathbf{x}}_{i} } \right)\), \({\ddot{\user2{u}}}_{{\varvec{i}}} = \frac{{\partial {\varvec{u}}_{i} }}{\partial t}\) with \({\varvec{u}}_{i} : = {\varvec{u}}\left( {{\varvec{x}}_{i} } \right)\) and \(N\) is the number of points in the horizon of node \(i\). The time integration procedure employs forward Euler method to compute the current values of accelerations, velocities, and positions of points. With this method, the explicit time integration of (36) yields the acceleration \({\ddot{{\varvec{u}}}}_{{\varvec{i}}}\), velocity \({\dot{{\varvec{u}}}}_{{\varvec{i}}}\) and displacement \({{\varvec{u}}}_{i}\) of a point \(i\) at current time \(t={t}^{n}\) as follows;

$$\begin{gathered} {\ddot{\user2{u}}}_{i}^{n} = \frac{{{\mathbf{\mathcal{L}}}_{i}^{n} + \user2{b}_{\user2{i}} }}{{\rho _{i} }} \hfill \\ {\ddot{\user2{u}}}_{i}^{{n + 1}} = {\dot{\user2{u}}}_{i}^{n} + {\dot{\user2{u}}}_{i}^{n} \Delta t \hfill \\ \user2{u}_{i} ^{{n + 1}} = \user2{u}_{i}^{n} + {\dot{\user2{u}}}_{i}^{{n + 1}} \Delta t \hfill \\ \end{gathered}$$
(37)

where \({\mathcal{L}}_{i}^{n}=\sum_{j=1}^{N}\left[\underset{\_}{\mathbf{T}}\left[{\mathbf{x}}_{i},t\right]\langle {\mathbf{x}}_{j}-{\mathbf{x}}_{i}\rangle -\underset{\_}{\mathbf{T}}\left[{\mathbf{x}}_{j},t\right]\langle {\mathbf{x}}_{i}-{\mathbf{x}}_{j}\rangle \right]{V}_{j}\). In the case of dynamic problems, the solution strategy is implemented as outlined above. However, for static problems, the numerical procedure is given some modifications. The modified computational scheme consists of integration over two time parameters, viz. a numerical time and a viscoelastic or material time. The numerical time is the time frame in which a steady-state solution is recovered from a dynamic problem by damping out the transient component. As the transient solution is damped out, the solution converges to the steady-state solution. In this approach, the quasi-static transient problem is first converted into a dynamic problem. Then fictitious inertia and damping forces are added to damp out the artificial dynamic motion. The method utilised in this work to attenuate the transient solution is called Adaptive Dynamic Relaxation (ADR) method which is an explicit iterative procedure that relies on mass proportional damping (see [71] for details).

7 Suppression of zero mode instability

Although the non-ordinary state-based peridynamic correspondence model offers a lot of advantages some of which have already been highlighted in Sect. 1 of this study, it suffers from a problem commonly referred to as the ‘zero-energy modes’ instability. This is known to cause oscillation in the deformation and stress field. While several techniques have been suggested to address this problem, including [20, 72,73,74,75,76,77], this study will utilise the method proposed in [73] for reason of simplicity to implement. This method supresses the zero-energy mode instability by introducing a supplementary artificial force density vector \({\underset{\_}{\mathbf{T}}}_{a}\left[{\varvec{x}},t\right]\langle {{\varvec{x}}}^{^{\prime}}-{\varvec{x}}\rangle\) on the bond \({{\varvec{x}}}^{^{\prime}}-{\varvec{x}}\) such that (18) becomes:

$$\underset{\_}{\mathbf{T}}\langle {\varvec{\xi}}\rangle =\underset{\_}{\omega }\langle \left|{\varvec{\xi}}\right|\rangle {\varvec{\sigma}}{\mathbf{K}}^{-1}{\varvec{\xi}}+{\underset{\_}{\mathbf{T}}}_{a}\left[{\varvec{x}},t\right]\langle {{\varvec{x}}}^{^{\prime}}-{\varvec{x}}\rangle$$
(38)

The artificial force density vector \({\underset{\_}{\mathbf{T}}}_{a}\left[{\varvec{x}},t\right]\langle {{\varvec{x}}}^{^{\prime}}-{\varvec{x}}\rangle\) in (38) is given by:

$${\underset{\_}{\mathbf{T}}}_{a}\left[{\varvec{x}},t\right]\langle {{\varvec{x}}}^{^{\prime}}-{\varvec{x}}\rangle =\frac{1}{2}\underset{\_}{\omega }\langle \left|{\varvec{\xi}}\right|\rangle {\varvec{C}}\underset{\_}{{\varvec{z}}}\langle {\varvec{\xi}}\rangle ,$$
(39)

where \(\underset{\_}{\omega }\langle \left|{\varvec{\xi}}\right|\rangle\) is the weight function, \({\varvec{C}}\left({\varvec{\xi}}\right)=c\left({\varvec{\xi}}\otimes{\varvec{\xi}}\right)/{\left|{\varvec{\xi}}\right|}^{3}\) is the tensor-valued symmetric micromodulus function, \(c=18k/\pi {\delta }^{4}\) is the bond force constant in the BBPD, and \(\underset{\_}{{\varvec{z}}}\langle {\varvec{\xi}}\rangle =\underset{\_}{\mathbf{Y}}\langle {\varvec{\xi}}\rangle -\mathbf{F}{\varvec{\xi}}\).

8 Model verification

In this section, the NOSBPD viscoelastic framework proposed will be verified against results obtained from analytical solution and Finite Element Method (FEM), either based on results obtained using a commercial FEM software (ANSYS), or results from published works. To model the viscoelastic response using the proposed NOSBPD viscoelastic formulation, computer codes were developed in MATLAB to solve three benchmark problems including a bar subjected to tensile load representing 1-dimensional and a plate subjected to tensile load representing a 2-dimensional quasi-static transient problems as well as axial vibration in a bar representing transient dynamic problem.

The material properties utilized in the following numerical examples is adopted from Abaqus Benchmark Manual (version 6.11). The relaxation function of the material is represented by the Prony series

$$E\left(t\right)={k}_{1}+{k}_{2}\mathrm{exp}\left(-\frac{t}{\tau }\right),$$
(40)

where \({k}_{1}=6.89\mathrm{MPa}\), \({k}_{2}=62.01\mathrm{MPa}\) and \(\tau =1.0\mathrm{sec}\). In this example, the bulk response is assumed to be elastic and the time independent bulk modulus \(K\) is \(689MPa\). From (40), the long-term Young’s modulus is \(E\left(\infty \right)={E}_{0}={k}_{1}\).

8.1 Viscoelastic bar subjected to constant axial load

Consider a bar having a length of \(1000\mathrm{ mm}\) and a square cross-section with dimensions \(10\mathrm{ mm}\times 10\mathrm{ mm}\). Let the bar be fixed at one end and a constant load of \(0.689\mathrm{ MPa}\) is suddenly applied at the other end. The solution strategy adopted is the complete viscoelastic formulation. The load is applied and held constant over a time period of \(50 s\). The time increment is taken to be \(\Delta t=0.0001 s\) which translates to a total time step of \(\mathrm{50,000}\). The displacement of points in the bar at initial and final times are presented in Fig. 3, while the displacement as well as the creep strain of the tip of the bar with respect to time are presented in Fig. 4.

Fig. 3
figure 3

Displacement of all point in the bar at a initial time, \(t=0\), b final time, \(t=50\mathrm{ s}\)

Full size image
Fig. 4
figure 4

a Displacement of bar tip with time b creep strain of bar tip with time

Full size image

Results presented in Fig. 3a, b respectively, shows that the instantaneous and long-term displacements predicted using the proposed framework agree well with analytical results. Results of creep strain presented in Fig. 4 also show good agreement with the analytical results.

8.2 Plate extension

This problem is a plane stress quasi-static transient analysis of a viscoelastic rectangular plate with dimensions \(1.0\mathrm{m}\times 1.0\mathrm{m}\) as shown in Fig. 5. The thickness, \(h\) of the plate is \(0.1\mathrm{ m}\). The left edge of the plate is subjected to a roller-type support. A load of \(1 \mathrm{MPa}\) is suddenly applied to the right edge of the plate and held constant for \(50\mathrm{ s}\) and then removed. The plate is then allowed to recover in a further 50 s thus bringing the total simulation time to \(100 \mathrm{s}\).

Fig. 5
figure 5

Plate geometry and loading

Full size image

To implement the numerical simulation, the plate was discretised into \(10\times 10\) grid giving a total of \(100\) material points. To utilise the recursive function (31) to update the stress at any time \({t}_{n+1}\), a numerical algorithm based on ADR methodology is utilised. This numerical algorithm consists of a two-tier time integration solution strategy. The first-tier time integration computes the viscoelastic deformation of the plate over a finite time increment \(\Delta t\). This time increment would be denoted as the ‘material’ time. The material time increment used in this simulation is \(\Delta t=0.1\mathrm{ s}\). The simulation is done over a total material time of \(t=100\mathrm{s}\). As a quasi-static transient problem, the plate is assumed to attain equilibrium over every material time step. This allows the treatment of the plate at every material time step as a static problem. Hence the ADR time integration methodology is utilised in this tier of the solution strategy to find the steady-state solution of the problem over each material time increment. The time increment in the ADR time integration is specified as 1. The time frame in the ADR algorithm is denoted as ‘ADR time’. Results from this simulation are compared with results from FE implementation in ANSYS. The viscoelastic displacement of the loaded edge of the plate with time using this proposed formulation is shown in Fig. 6, alongside prediction from ANSYS. Also, the displacement field in both \(x\) and \(y\) directions corresponding to loading time \(t=50\mathrm{ s}\) (when the load is removed) and \(t=100\mathrm{ s}\) (at the end of simulation) are, respectively, shown in Figs. 7, 8.

Fig. 6
figure 6

Prediction of displacement of the loaded edge of the plate with time along a x-direction, and b y-direction

Full size image
Fig. 7
figure 7

Displacement field, \(u\), corresponding to \(t=50 \mathrm{s}\) a \({u}_{x}\) computed using proposed framework b \({u}_{y}\) computed using proposed framework c \({u}_{x}\) computed using FEM—ANSYS d \({u}_{y}\) computed using FEM—ANSYS

Full size image
Fig. 8
figure 8

Displacement field, \(u\), corresponding to \(t=100 \mathrm{s}\) a \({u}_{x}\) computed using proposed framework b \({u}_{y}\) computed using proposed framework (c) \({u}_{x}\) computed using FEM—ANSYS d \({u}_{y}\) computed using FEM—ANSYS

Full size image

As can be observed from Figs. 6, 7, 8, the prediction of the transient displacement of the loaded edge of the plate as well as the displacement field in the plate at times \(t=50 s\) and \(t=100 s\) using the proposed method agrees well with FE prediction using ANSYS, thus further reinforcing the validity of the proposed method to model time dependent behaviour of materials.

8.3 Axial vibration of a viscoelastic bar

In this section, the proposed NOSBPD viscoelastic formulation is applied to simulate the axial vibration of a viscoelastic bar. This example will be used as a probe to examine the suitability of the proposed framework in the viscoelastodynamic regime. The bar is clamped at one end and subjected to an initial sudden stretch at the other end for a short period of time after which the load is removed. The viscoelastic material property of the bar is defined by the Prony series coefficients given in Sect. 8. The density of the material is arbitrarily taken to be \(10000 \mathrm{Kg}/{\mathrm{m}}^{3}\). The bar has a length of \(1.0 m\) and a cross sectional area of \(1.0\times {10}^{-6} {\mathrm{m}}^{2}\). To implement this problem in the numerical scheme of the proposed formulation, the bar is discretised into \(1000\) material points, so that each material point has a size of \(\Delta x=1.0\times {10}^{-3} m\). A horizon size, \(\delta =3\Delta x\) is adopted. Numerical time integration solution technique [71] is utilised to track the position of material points with passage of time due to applied load. The numerical integration was implemented over a total simulation time of \(2.8065 s\) with a time interval, \(\Delta t=3.742\times {10}^{-5}\) s.

The transient displacement of two points located at \(x=0.4995\mathrm{ m}\) and \(x=0.9995\mathrm{ m}\) are shown in Fig. 9. For comparison, the response of the bar subjected to the same loading condition but with equilibrium material properties corresponding to time \(t=\infty\) was also computed using analytical solution method. The equilibrium material property yields elastic response of the material. It is observed from the figures that compared with the elastic response; the viscoelastic property progressively damps out the oscillation thereby creating a response that is asymptotic to the steady-state solution. This attenuating effect that viscoelastic properties have on vibration of a medium is called viscoelastic damping [78, 79]. Viscoelastic materials added to structural elements and devices to enhance their damping characteristics are called viscoelastic damping materials (VDM). Characterisation of static and dynamic response of VDM is an active field of research in many disciplines including automotive and aerospace industries.

Fig. 9
figure 9

Displacement at a centre of bar, \(x=0.4995\mathrm{ m}\), and b at the tip of bar, \(x=0.9995\mathrm{ m}\)

Full size image

Introduction of nonlocality into continuum models is known to introduce dispersion in wave propagation [80]. The horizon is the parameter that defines the degree of nonlocality in the NOSBPD theory, and just like other nonlocal continuum models, the degree of nonlocality also plays important role in the behaviour of peridynamic solution. To investigate the influence of horizon size on the dynamic response of viscoelastic materials, the current viscoelastodynamic problem is simulated using varying horizon sizes of \(\delta =3\Delta x\), \(10\Delta x\), \(15\Delta x\), and \(20\Delta x\). The result for this parametric study showing the trend of dispersion behaviour is presented in Fig. 10.

Fig. 10
figure 10

Effect of horizon size on vibration of viscoelastic bar

Full size image

Numerical results presented in Fig. 10 indicate that the horizon influences the period of oscillation while the amplitude of vibration remains almost the same. This dispersive behaviour due to the introduction of nonlocal interaction has also been observed in a one-dimensional, viscoelastic bond-based peridynamic medium [81], as well as in the dynamic analysis of nanobeams [79, 82]. The ability to control the dispersive behaviour of viscoelastic materials through a nonlocal parameter is a very important strength of this proposed framework that promises to find application in engineering of nano materials and devices where size effects due to dimensional constraints for example, affects mechanisms of deformation such as dislocation motion. In such situation where size dependent response becomes an important consideration, the use of local models based on CCM becomes questionable as they do not admit intrinsic size dependence in their solution, and so nonlocal models are increasingly being utilised to characterise such materials and have been shown to better predict response of materials than local models based on CCM theory [83].

9 Conclusion

This paper presents a NOSBPD viscoelastic framework for simulation of materials exhibiting rate dependent behaviour. The proposed formulation leveraged on a constitutive correspondence principle to allow highly established viscoelastic models from CCM to be implemented in NOSBPD framework. Being a nonlocal theory, and endowed with the capability of modelling continuous and discontinuous material response using the same governing equation of motion without the need for modifications, a NOSBPD viscoelastic formulation promises to provide a computational framework with capacity to model viscoelastic materials in the presence of discontinuities such as initiation and propagation of cracks or nonlocal response due to scale effects as observed in small scale material characterisation.

To validate the proposed framework, three numerical problems representing 1-D and 2-D quasi-static transient and 1-D transient dynamic problems were solved. In the quasi-static regime, the proposed framework was utilised to characterise the response of a bar and a plate subjected to a constant tensile load. The viscoelastic transient response predicted by the proposed framework shows good agreement with results from analytical and FEM methods.

In the viscoelastodynamic regime, results of simulating the axial vibration of a bar also shows good agreement with results obtained using FEM. Parametric study indicated that the influence of nonlocality on the dispersion behaviour of the bar is limited to the period of oscillation and not so much influence is observed on the amplitude of vibration. Characterization of nonlocal behaviour is made possible because of the notion of horizon introduced by the NOSBPD framework as a material parameter.

Although one of the benefits to be derived in developing this framework is the effectiveness of PD to model materials with strong discontinuous behaviour such as the existence of cracks, discussion on the application of this framework to characterise materials with cracks is deferred to a future contribution.