Review of Nonlocal-in-Time Damping Models in the Dynamics of Structures

Abstract

In the present paper, the nonlocal-in-time damping models, called “damping-with-memory” models, are reviewed. Since such models do not involve the distribution along the longitudinal coordinate, they are easily adjustable for the FEM (Finite Element Model) algorithm, which is a big advantage due to the fact that FEM is the most-used method in engineering calculations. Within damping-with-memory models, the internal damping of a structure at the current time, is assumed to be dependent not only on the instant strain-rate magnitude or displacement-velocity magnitude but also on the strain-rate or velocity values along the previous time history. The greater the gap between the two time points, the lower the influence that one of them has on the other. To implement a composite beam vibration simulation involving damping with memory, the equation of motion of a structure written in the matrix form could be solved using the central difference method. The models constructed could be calibrated based on 3D numerical simulation data with the least squares method. It has been shown that the results obtained using the implementation of a calibrated damping-with-memory model within the 1D finite-element beam algorithm are in good correlation with those given by the 3D-FEM numerical simulation data.

1. Introduction

The problem of modeling the vibrations in structural elements made of modern composites and nanomaterials is a significant challenge. Unlike traditional building materials, such as metals, composites are highly heterogeneous and anisotropic [1,2].

One of the features of such structural materials is that the physical properties of matrices and binders that comprise their basis differ significantly. By using various combinations of parameters and orientations of inorganic and organic materials in manufacturing composites, it is possible to control their properties and to achieve the desired characteristics and parameters, allowing one to solve different problems in a wide variety of manufacturing industries [1,2].

Classic hypotheses are often not enough to describe the behavior of such materials. Hence, it is suggested that mathematical models are especially flexible and controllable to simulate the dynamic response of structures considering their orthotropic or other types of anisotropy properties [3,4].

The aim of this review is to overview recent publications and results addressing the application of the models with nonlocal-in-time damping for solving dynamics problems of structural mechanics.

2. Mathematical Modeling of the Internal Friction

The problem of modeling the internal friction in composite or even homogeneous materials is especially challenging. Currently, this problem does not have an unambiguous solution. Internal friction in a material means energy-irreversible processes that accompany the cyclic deformation of bodies [5,6,7]. It is believed that the irreversibility of these processes is characterized by the fact that part of the deformation energy for each cycle is converted into heat or electricity and dissipated. In the literature, there exist several measures of internal friction [5], among them: the tangent of the mechanical loss angle; specific energy dissipation; and the inverse value of the quality factor, etc.

There have been over one hundred and fifty years of research into the issues of internal friction in a deformable material [8,9,10,11,12,13]. During this period, many different and often contradictory theories and models have been developed that describe the process of the dissipation of the energy of vibrations. For example, according to the Kelvin–Voigt hypothesis [14], energy losses due to internal friction in a material are proportional to the rate of change in cyclic deformations. On the contrary, in the theory developed by Maxwell, the internal absorption parameter is inversely proportional to the vibration frequency [15], while Sorokin’s hypothesis of complex stiffness [7] is frequency independent.

None of the currently existing hypotheses is universal or applicable to various materials within a wide frequency range. It is clear that only a stable correspondence to actual dynamic processes could confirm the reliability of one theory or another.

The monograph in Reference [16] was devoted to the problems of vibration and noise in various fields of mechanical and civil engineering, wherein different ways of describing various types of damping were considered. Particular attention was paid to models of viscoelastic damping, which determine the behavior of composites: polymer and vitrified materials. At the same time, it was noted in Reference [16] that during the manufacturing of composites, as a rule, producers do not ascribe to the goal of achieving high damping characteristics. Therefore, in laboratory or ambient experiments, it is difficult to separate the pure damping effects from the effects arising from the nonlinear behavior of the material. Thus, it is necessary to describe the influence of viscoelastic damping on the dynamic behavior of structures.

Adhikari [17] presented an analysis of complex systems with damping. The main object of research was a multiple-degree-of-freedom system with viscous and/or non-viscous internal friction. The theory of a viscoelastic-hereditary medium was used as the theory of inviscid friction, and its particular case was the model of viscous friction in Reference [18].

3. Non-Classical Models of Composite and Viscoelastic Materials

Before we focus on internal friction modeling, the non-classical models of composite materials considering general deformation should be discussed. Various conceptual and mathematical models have been proposed to take into account the peculiarities of composite and viscoelastic materials’ response under loading, among them: gradient models [19]; mathematical models involving the apparatus of fractional calculus [20,21,22]; nonlocal models [23]; and nonlocal models via operators of fractional order [24,25,26,27].

Flügge, in his book on viscoelasticity [28] (see p. 74), noted that “The response q(x₁) at some point x₁ obviously depends, of course, not only upon the local value w(x₁) of the deflection but also upon that at neighbor points x₂. Their influence decreases as the distance |x₁ − x₂| increases.” Thus, Flügge came to the conclusion that in order to solve practical engineering problems, it is necessary to study a class of models that are now called nonlocal.

Lei et al. [4] emphasized that from a physical point of view, the need to take into account nonlocal elastic or damping properties arises when it is expedient to model a two-dimensional or spatial element as a one-dimensional one. Therefore, such an approach has many potential applications in engineering structures with nonlocal energy dissipation mechanisms.

In 1972, Eringen and Edelen [29] proposed a nonlocal model of the elastic properties of a structural material, which was further developed by many researchers, among them References [30,31,32].

In 1975, Ahmadi studied a linear theory of nonlocal viscoelastic material in Reference [33] and obtained the constitutive equations for nonlocal Kelvin–Voigt, Maxwell, and Boltzmann–Volterra viscoelastic materials. It was shown that the simple viscoelastic materials, in fact, are nonlocal-in-time properties.

Baretta et al. [34] suggested a combination of Eringen’s nonlocal elasticity model and the fractional derivative model of viscoelasticity to solve the problem of bending of Euler–Bernoulli beams in porous viscoelastic materials.

Another way to generalize the linear relationship between stress and strain was proposed by Boltzmann [35] and Volterra [36], who developed the theory of a viscoelastic medium with heredity, according to which the strain at any given time depends on all that has happened before, i.e., on the whole stress history. The final deformation at the considered moment of time is represented by the sum of the deformations caused by each of the previously acting forces, considering the decrease in their influence over the current time. The degree of such a decrease is described by a memory function that monotonically decreases with an increase in the argument, which characterizes the hereditary properties of the material. The type of hereditary kernel should be selected based on experimental data.

A great contribution to the development of the hereditary theory of elasticity was made by Academician Rabotnov [13,37,38]. The fractional–exponential kernel of the hereditary elasticity proposed by him [37] enabled results that are in good agreement with the results of testing polymeric materials to be obtained.

The monograph in Reference [39] was devoted to the construction of nonlocal models (Kunin chains). It has been noted that nonlocality could be of a physical or geometric nature. In the first case, the scale parameter corresponds to the characteristic size of the material’s microstructure. In the second case, it is due to the approximate consideration of such parameters as the thickness of a structural element (rod or plate), while the nonlocal model serves for an effective approximate description of the behavior of a homogeneous three-dimensional medium. The book [39] presented in detail the construction of the mathematical apparatus of nonlocal theories and also noted a number of inherent ambiguities, such as the ambiguity of the certainty of energy density, stresses, etc.

Nonlocal elasticity models were used for theoretical studies of nanorod and nanoplate systems in References [40,41]. The nonlocal elasticity model was applied to the wave propagation modeling in Reference [42].

Bažant and Jirásek [23] reviewed the progress in nonlocal models of the integral type and discussed their physical justifications, advantages, and numerical applications. It has been noted that “the term nonlocal has in the past been used with two senses, one narrow and one broad. In the narrow sense, it refers strictly to the models with an averaging integral. In the broad sense, it refers to all the constitutive models that involve a characteristic length (material length), which also includes the gradient models. This broad sense stems from the realization that some gradient models are derived as approximations to the nonlocal averaging integrals and that for all the gradient models, the gradient, in fact, includes a dependence on the immediate (infinitely close) neighborhood of the point under consideration”. An interested reader could find a state-of-the-art in the field of application of such nonlocal models in mechanics in Reference [23].

Potapov [43,44] considered the problem of the dynamic stability (in the sense of Lyapunov) of a composite rod with the simultaneous application of the hypotheses of nonlocal elasticity and nonlocal damping. Although this approach, on the one hand, makes the model very flexible, which is very important when working with orthotropic and anisotropic materials, on the other hand, the simultaneous consideration of both nonlocal elasticity and nonlocal damping makes the model quite complex from a mathematical point of view.

Suppose that in order to achieve the required calculation accuracy, it is sufficient to use either only the hypothesis of nonlocal elasticity or only the hypothesis of nonlocal damping. Further, we would focus on the hypothesis of nonlocal damping.

4. Nonlocal Damping Models

For the first time, a damping model, nonlocal-in-spatial coordinates, was proposed by Russell [45] for the dynamic analysis of a composite beam.

Banks and Inman [3] examined four damping models using the example of vibrations of a composite cantilever beam with a concentrated mass at the free end. The beam was made of fiberglass-reinforced plastic. In so doing, the reinforcement was arranged in the longitudinal and transverse directions. It was noted that this material possesses damping properties that are different from those of homogeneous materials. It was assumed that the behavior of the beam corresponded to the Bernoulli hypothesis, and the torsional, shear, and longitudinal displacements were negligible.

Three models presented in Reference [3] describe internal friction: the Voigt model; time hysteresis; and spatial hysteresis. The model called temporary hysteresis in Reference [3] is similar to the theory of the medium with heredity, where the system’s response to external influences that took place in the past affects its dynamic behavior in the present.

Spatial hysteresis is similar to the damping model nonlocal-in-coordinates [45]. This model is based on the assumption that the energy losses during lateral vibrations of the beam are associated with the fact that the internal friction in the material is caused by the rotation of the sections relative to each other. Moreover, not only the considered cross-section of the beam affects the process of internal damping but also the cross-sections adjacent to it. In this case, a kernel function or an influence function should be introduced, which is the law of a decrease in the influence of cross-sections on each other with an increase in the distance between them.

The fourth model in Reference [3] described external viscous damping. In other words, it assumed that the vibration damping process was associated with the friction of the structure against the external medium. The damping forces were considered proportional to the displacement rates. This hypothesis is convenient for practical applications.

Some combinations of these models are also considered in Reference [3]. It is noted that the described approaches are based on physical considerations. Damping models are involved in the equation of bending vibrations of the Euler–Bernoulli beams. In this case, the boundary conditions were chosen so that they could be used with different damping models.

The selection of the parameters for damping models was carried out in Reference [3] according to the experimental data using the least squares method. The parameters obtained in this way were utilized in integrodifferential equations to simulate the dynamic response of the system. The vibration process modeled in such a way was compared with experimental data. It has been shown that the spatial hysteresis used as a model of internal friction in combination with external viscous damping gives the best agreement with the experiment.

Banks and Inman [3] emphasized that the intuitive idea physically confirmed that external damping is significantly manifested at lower vibration frequencies, while internal friction has a stronger effect on the vibration process at higher frequencies. It is also noted that the proposed damping models cannot be constructed using standard damping coefficients obtained by modal analysis since they completely hide the physics and the mechanism of the damping process. The natural frequencies obtained as a result of mathematical modeling for the considered beam were compared with those measured experimentally. It has been noted that the results of numerical calculations of a beam with plastic fiberglass reinforced, according to the Euler–Bernoulli theory, coincide well with the experimental data, and there is no need to use the Timoshenko beam model.

For the dynamic analysis of Euler–Bernoulli beams and Kirchhoff plates, Lei et al. [4] suggested using a nonlocal damping model, which includes spatial and temporal hysteresis. In contrast to the usual local models, the damping forces are calculated as the average of the velocity field in a certain region determined by the kernel function. The resulting integro-differential equation of motion in partial derivatives is solved by the Bubnov–Galerkin method.

It was emphasized by the authors of [4] that a new nonlocal damping model proposed for analyzing the dynamic characteristics of bending beams under various boundary conditions is the generalization of the classical model of viscous friction and could be applied in the calculations of engineering structures when it is important to take into account the nonlocal mechanism of energy dissipation during vibrations. Such structures include those with a viscoelastic damping coating, structures on viscoelastic foundations, structures made of composite and nanomaterials, etc.

Within such an approach, it is assumed that damping at some point of the bar with coordinate $x_1$, measured along its axis, is considered dependent not only on the local value of the rate of deformation change $\dot{\epsilon }(x_1,t)$ at the same point $x_1$ but also on the value of the rate of change in deformations at neighboring points at some area adjacent to this point. The less the degree of influence of the damping properties of the points under consideration on each other is, the greater the distance between them.

The approach proposed in Reference [4] was further utilized for solving different engineering problems, among them: a modal analysis [46] and a finite element analysis of the Euler–Bernoulli beams with internal nonlocal damping [47,48]; a finite element method for elastic beams resting on nonlocal foundations [49]; and Timoshenko beams with internal nonlocal damping [50].

Reference [48] concluded that “the external damping parameters have simple effects on the natural frequencies and the dependence with nonlocal parameter is not strong”. Hence, the nonlocal-internal and -external damping parts in the equation of motion can be modeled independently. Therefore, the modal analysis technique for the nano-scale Euler–Bernoulli beams was proposed with the consideration of damping that is nonlocal both in time and space in Reference [51].

A novel dynamic finite-element approach for bending vibrations of damped nonlocal-in-space beams on an elastic foundation was proposed in Reference [51]. Both internal and external damping mechanisms were employed. The stiffness and mass matrices of the nonlocal beam were obtained using the conventional finite element method. The scale factor of the nonlocal model was not determined by any experimental data, and for numerical examples, it was considered from 0 to 2 nm [51].

Emphasizing the physical significance of nonlocal models, Gonzales et al. [52] utilized the approach by Lei et al. [4,46,47,48,49,50] for treating the Euler–Bernoulli beams with separate areas of nonlocal properties via the Galerkin method. In Reference [53], the models of spatial and spatiotemporal nonlocal-internal damping were applied to the nonlinear seismic analysis of a four-story frame. It was shown that nonlocal damping models do not exhibit spurious damping properties as exhibited by the global damping models on the onset of inelasticity.

We should also mention the series of theoretical papers [54,55,56,57,58,59,60] which are devoted to the nonlocal energy damping model with linear and nonlinear terms, wherein the energy dissipation is considered in Euler–Bernoulli beams and Kirchhoff plates with different boundary conditions. The authors of [54,55,56,57,58,59,60] declared that the model is suitable for flight structures, pipes, and plate dynamic simulation, including supersonic flatter.

The significant contribution to the development and application of the nonlocal damping theory to practical engineering problems was made by Potapov [61,62]. In 2012, he implemented the nonlocal-in-space damping model to the problem of the stability of vertical columns subjected to longitudinal force [61]. The solution was found using the Bubnov–Galerkin and the Runge–Kutta fourth-order methods. Further, this approach was utilized for modeling nonlinear systems, and the problem of shallow arch stability was solved [62].

To this purpose, instead of the traditional Voigt hypothesis

$$\sigma =E\epsilon +E{t}_{\kappa } \dot{\epsilon },$$

(1)

where σ and $\epsilon$ are the normal stress and longitudinal strain, respectively, $\dot{\epsilon }$ is the strain rate, E is Young’s modulus, and $t_{\kappa }$ is the retardation time, the following relationship was used:

$$\sigma \left(x,t\right)=E\epsilon \left(x,t\right)+\chi \underset{0}{\overset{l}{\int }}C\left(\left|x-\theta \right|\right)\dot{\epsilon }\left(\theta ,t\right)d\theta ,$$

(2)

where $C\left(\left|x-\theta \right|\right)$ is a function that describes the decrease in the influence of points with coordinates $x$ and $\theta$ on each other, and $\chi$ is the dynamic viscosity.

The research started by Potapov and his collaborators [61,62,63,64] on modeling internal damping nonlocal-in-spatial coordinates was further developed in References [65,66], resulting in a technique allowing one to calibrate the model parameters using the data of numerical simulations of the vibrations of a composite beam in a three-dimensional configuration, taking its orthotropic properties into account. The results presented in References [65,66] indicate that a one-dimensional nonlocal model could be used instead of a detailed classical spatial computational model.

Although the nonlocal-in-space damping model shows good agreement with the results of numerical experiments [65,66], its embedding into the algorithm of a numerical calculation method, for example, the finite element method, is rather difficult and, to some extent, farfetched. This is due to the fact that the kernel function describing the character of the decrease in the influence of neighboring points of an element on each other can step over the boundaries of finite elements.

But since the finite element method is the predominant method of analysis in engineering practice, it becomes necessary to construct a model that is flexible enough to describe the dynamic behavior of a composite material and, at the same time, is easily integrated into the finite element method algorithm. Such a model similar to that for damping nonlocal-in-time, proposed in Reference [3], was developed further by Sidorov et al. [67], with the only difference that the authors of Reference [67] assumed that damping at the current moment of time $t$ depends not only on the instantaneous value of the rate of change in deformations at this moment of time $\dot{\epsilon }\left(t\right)$ but also on the value of the rate of change in deformations at the previous moments of time $\tau$. Moreover, the less the effect of the value of the rate of change in deformation at a certain moment of time $\tau$ is, the longer the time interval between $\tau$ and the current moment $t \mathrm{i}\mathrm{s}$. That is why we will refer to this model as damping-with-memory.

See

Table 1. Basic models of internal damping.

Basic Classical Models of Damping due to Internal Friction
Kelvin–Voigt model	$\sigma =E\epsilon +t_{\kappa }\dot{\epsilon }$	[3,8,50,63,66]
Sorokin’s complex stiffness model	$\sigma =\left(1+i\gamma \right)E\epsilon$	[3]
Nonlocal models of damping due to internal friction
Internal friction model nonlocal-in-spatial coordinate x	$\sigma \left(x,t\right)=E\epsilon \left(x,t\right)+\chi \underset{0}{\overset{l}{\int }}C\left(\left|x-\theta \right|\right)\dot{\epsilon }\left(\theta ,t\right)d\theta$	[3,50,63,64,65,66]
Internal friction model nonlocal-in-time t	$\sigma \left(x,t\right)=E\epsilon \left(x,t\right)+\chi {\int }_{-\infty }^{t}G\left(t-\tau \right)\dot{\epsilon }\left(x,\tau \right)d\tau$	[3,67,68,69,70]

for the main classical and nonlocal models of internal friction for solids and structures. The main types of kernel functions $C\left(\left|x-\theta \right|\right)$ and $G\left(t-\tau \right)$ utilized in the spatial and temporal nonlocal models, respectively, are presented in

Table 2. Kernel functions used in nonlocal damping models.

Kernel Functions for Nonlocal-in-Space Damping Models
Exponential kernel function	$C\left(\left|x-\theta \right|\right)=\frac{\mu }{2} e^{-\mu \left|x-\theta \right|}$	[50,63,64,65,66]
Error kernel function	$C\left(\left|x-\theta \right|\right)=\frac{\mu }{\sqrt{2\pi }} e^{\frac{-{\mu }^2{(x-\theta )}^2}{2}}$	[3,50,53]
Hat kernel function	$C\left(\left|x-\theta \right|\right)=\left\{\begin{array}{c}\frac{1}{l_0}, \mathrm{for} \left|x-\theta \right|\le \raisebox{1ex}{$l_0$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ 0, \mathrm{otherwise} \end{array}\right.$	[50]
Triangular kernel function	$C\left(\left|x-\theta \right|\right)=\left\{\begin{array}{c}\frac{1}{l_0}\left(1-\frac{\left|x-\theta \right|}{l_0}\right), \mathrm{for} \left|x-\theta \right|\le l_0\\ 0, \mathrm{otherwise}\end{array}\right.$	[50]
Kernel functions for nonlocal-in-time damping models
Exponential kernel function	$G_v\left(t-\tau \right)=\mu e^{-\mu \left(t-\tau \right)}$	[67]
Error kernel function	$G_{\nu }\left(t-\tau \right)=\frac{2\mu }{\sqrt{\pi }}{e}^{-{\mu }^2{\left(t-\tau \right)}^2}$	[67,69,70,71]

, where $\mu$ and l₀ are the characteristic parameters of the damping material [50], and γ is the coefficient of inelastic reaction due to Sorokin [3].

5. Damping-with-Memory Model

Within the algorithm of the finite element analysis, the equilibrium equation of a structure deformed during its motion could be represented in the matrix form as follows [71]:

$$M\ddot{V}\left(t\right)+D\dot{V}\left(t\right)+KV\left(t\right)=F\left(t\right),$$

(3)

where $V\left(t\right)$ is the displacement vector, $M$ is the mass matrix, $K$ is the stiffness matrix of the finite element model, $F\left(t\right)$ is the load vector, while $D$ is the damping matrix.

The stiffness and mass matrices of the i-th beam element have the form [71]:

$$K_i=\left[\begin{array}{cccc}\frac{12EI}{l^3}& \frac{6EI}{l^2}& -\frac{12EI}{l^3}& \frac{6EI}{l^2}\\ \frac{6EI}{l^2}& \frac{4EI}{l}& -\frac{6EI}{l^2}& \frac{2EI}{l}\\ -\frac{12EI}{l^3}& -\frac{6EI}{l^2}& \frac{12EI}{l^3}& -\frac{6EI}{l^2}\\ \frac{6EI}{l^2}& \frac{2EI}{l}& -\frac{6EI}{l^2}& \frac{4EI}{l}\end{array}\right],$$

(4)

$$M_i=\frac{\rho Al}{420}\left[\begin{array}{cccc}156& 22l& 54& -13l\\ 22l& 4l^2& 13l& -3l^2\\ 54& 13l& 156& -22l\\ -13l& 3l^2& -22l& 4l^2\end{array}\right],$$

(5)

where l is the element length, $I$ is the moment of inertia, and $\rho$ is the material density.

Assuming that the memory model of material damping depends on the values of the strain rate, its damping matrix could be deduced according to the requirement of stationary state of the full energy of the vibrating system. Thus, the dissipation of energy during deformation of a finite element under dynamic loading will be represented by a dissipative function ${\mathsf{\Phi }}_D$ [72,73]:

$${\mathsf{\Phi }}_D=\frac{1}{2}\chi {\dot{\epsilon }}^2,$$

(6)

where the coefficient defining material viscosity $\chi$ could be obtained as follows:

$$\chi =E{t}_{\kappa },$$

(7)

and $t_{\kappa }$ is the retardation time (sometimes it is called the delay time [72]) which could be determined in terms of the damping coefficient (critical fraction) ${\xi }_{cr}$ and the first natural frequency of the system $\omega$ [72,73]

$$t_{ \kappa }=\frac{2{\xi }_{cr}}{\omega }.$$

(8)

Energy dissipation through the material of the whole system could be presented as the summation of energy dissipation ${\mathsf{\Phi }}_{Di}$ through each finite element:

$${\mathsf{\Phi }}_D=\sum _{i-1}^N{\mathsf{\Phi }}_{Di},$$

(9)

where $i$ is the finite element number $(i=1, 2, \dots , N)$, $N$ is the number of elements in the whole FE computational model, which, in the case of a bent beam, has the form [68]

$${\mathsf{\Phi }}_{Di}=\frac{1}{2}\underset{l}{\overset{}{\int }}\underset{A}{\overset{}{\int }}\chi {\dot{\epsilon }}_b^{2}dAdz,$$

(10)

where A is the element’s cross-section area, z is the longitudinal coordinate, l is the element length, and ${\dot{\epsilon }}_b$ is the bending-induced axial strain rate.

Bending-induced axial strain for the Euler–Bernoulli beam is as follows [73]:

$${\epsilon }_b=\frac{1}{R}y=\frac{d^{2}v}{d{z}^2}y,$$

(11)

where R is the radius of curvature of the beam’s neutral layer, y is the distance to the considered beam fiber from its neutral layer, and $v$ is the transverse displacement.

Then Equation (10) is reduced to

$${\Phi }_{Di}=\frac{1}{2}I\chi \underset{l}{\overset{}{\int }}{\left(\frac{d^2\dot{v}}{d{z}^2}\right)}^{2}dz.$$

(12)

Within the FEA, the transverse displacements could be approximated along the beam element with cubic shape function:

$${\left[N_v\right]}^T=\left[\begin{array}{c}1-3{\xi }^2+2{\xi }^3\\ l(\xi -2{\xi }^2+{\xi }^3)\\ 3{\xi }^2-2{\xi }^3\\ l(-{\xi }^2+{\xi }^3)\end{array}\right],$$

(13)

where $\xi =\frac{z}{l}$ is the reduced local longitudinal coordinate.

Considering (13), the Expression (12) could be reduced to

$${\mathsf{\Phi }}_{Di}=\frac{1}{2}I\chi \underset{0}{\overset{1}{\int }}{\left(A_v\left[N_v\right]\right)}^T{A}_v\left[N_v\right]{\dot{v}}_i^{2}ld\xi .$$

(14)

where ${\dot{v}}_i=\left(\begin{array}{c}{\dot{v}}_0\\ {\dot{\phi }}_0\\ {\dot{v}}_l\\ {\dot{\phi }}_l\end{array}\right)$ are nodal velocity vectors during bending. In so doing, ${\dot{v}}_0$ and ${\dot{v}}_l$ are the velocities of the transverse displacements at the initial and terminal points of each i-th element, and ${\dot{\phi }}_0$ and ${\dot{\phi }}_l$ are velocities of rotation displacements, and $A_v=\frac{d^2}{l^2·d{\xi }^2}$.

The requirement for the stationary state of the full energy of the vibrating system

$$\frac{\partial {\mathrm{\Phi }}_D}{\partial \dot{v}}=0$$

(15)

gives us the $D\dot{\overline{V}}$ term in the equilibrium equation in Motion (3).

The beam’s damping matrix $D$ is obtained by topological summation of the element’s damping matrices $D_i$ as follows:

$$D_i=I\chi \underset{0}{\overset{1}{\int }}{\left(A_v\left[N_v\right]\right)}^T{A}_v\left[N_v\right]ld\xi ,$$

(16)

or

$$D_i=\chi ·\left[\begin{array}{cccc}\frac{12I}{l^3}& \frac{6I}{l^2}& -\frac{12I}{l^3}& \frac{6I}{l^2}\\ \frac{6I}{l^2}& \frac{4I}{l}& -\frac{6I}{l^2}& \frac{2I}{l}\\ -\frac{12I}{l^3}& -\frac{6I}{l^2}& \frac{12I}{l^3}& -\frac{6I}{l^2}\\ \frac{6I}{l^2}& \frac{2I}{l}& -\frac{6I}{l^2}& \frac{4I}{l}\end{array}\right].$$

(17)

To construct the damping-with-memory model, Equation (3) should be represented as follows: [67]

$$M\ddot{V}\left(t\right)+D{\int }_0^{t}G\left(t-\tau \right)\dot{V}\left(\tau \right)d\tau +KV\left(t\right)=F\left(t\right),$$

(18)

where $G(t-\tau )$ is the kernel function which describes the decrease in the influence of the motion velocity at the moment τ on damping at the current moment t. The kernel function should be normalized to satisfy the following condition:

$${\int }_{-\infty }^{t}G\left(t-\tau \right)d\tau =1.$$

(19)

Considering (19), the error function could be used as the kernel function

$$G\left(t-\tau \right)=\frac{2\mu }{\sqrt{\pi }}{e}^{-{\mu }^2{\left(t-\tau \right)}^2},$$

(20)

where $\mu$ is a parameter of the model that determines the rate of decay of the influence of the nonlocal damping over distance, i.e., the parameter characterizing the level of damping temporal nonlocality [68]. The higher is $\mu$, the closer the model is to the classic local one (Figure 1).

6. Equilibrium Equation Solving with the Method of Central Differences

To solve the dynamic equilibrium equation, the method of the central differences could be adopted [69], what allows one to reduce Equation (18) to the following form [67]:

$$\frac{1}{\Delta t^2}M\left(V_{i+1}-2V_i+V_{i-1}\right)+\frac{D}{2}Z+\frac{1}{2\Delta t}D\left(V_{i+1}-V_i\right)+K{V}_i\left(t\right)=F_i,$$

(21)

where $i=1, 2, 3\dots$ is a number of the considered moments in time $t$, and $\Delta t$is the time step.

In Equation (21), the vector has the form

$$Z=\sum _{j=1}^i\frac{2\mu }{\sqrt{\pi }}{e}^{-{\mu }^2{\left[t-\left(\tau -\Delta t/2\right)\right]}^2}\left(V_j-V_{j-1}\right),$$

(22)

where $\tau =j \Delta t$, and $t=i\Delta t$.

From Equation (21), it follows that the displacement vector $V_{i+1}$ could be calculated in terms of $V_i$ and $V_{i-1}$ as follows:

$$V_{i+1}=Q{F}_i-Q_1{V}_i-Q_2{V}_{i-1}-Q_{3}Z,$$

(23)

where

$$\begin{array}{l}Q={\left(\frac{1}{\Delta t^2}M+\frac{1}{2\Delta t}\mathrm{D}\right)}^{-1},\\ Q_1=Q\left(-\frac{2}{\Delta t^2}M-\frac{1}{2\Delta t}D+K\right),\\ Q_2=Q\left(\frac{1}{\Delta t^2}M\right),{ Q}_3=\frac{1}{2}QD.\end{array}$$

(24)

At the first step, i.e., for $i=1$, $V_0=0$ and $V_1=0$ are taken as the initial conditions.

The solution obtained by the explicit scheme turned out to be quite unstable and demanding in a sense of calculation time. For the explicit scheme, the maximum allowed time increment was 0.001 s. For the larger time steps, the finite element model fell apart [67]. Thus, it was further decided to use the Newmark method (the implicit scheme) instead, since it has been found in Reference [69] that the implicit scheme allows the maximum time increment of 0.01 s, i.e., 10 times larger in comparison to the central difference method.

7. Modified FEA Model Considering Nonlocal Damping Solved by the Implicit Scheme

In order to make the modal more flexible and controllable, the additional terms have been added to the FEA equation of Motion (3) [69,70]:

$$M\ddot{\overline{V}}\left(t\right)+\alpha D\dot{\overline{V}}\left(t\right)+\left(1-\alpha \right)D\left[\underset{t_0}{\overset{t-∆t}{\int }}G\left(\mu ,t-\tau \right)\dot{\overline{V}}\left(\tau \right)d\tau \right]+K\overline{V}\left(t\right)=\overline{F}\left(t\right).$$

(25)

Similar to Reference [4], the material had separate parts of locality and nonlocality, the damping properties were also divided into nonlocal and local parts, and $\alpha$ was the coefficient of nonlocality, i.e., the parameter assigning the share of local properties of the model ($0<\alpha <1$).

The equation of motion (25) could be solved using the modified Newmark method [71]. The computational scheme for the vertical displacements of the beam is as follows:

$$\overline{V}\left(t_{i+1}\right)=K_{ef}^{-1}{R}_{ef}\left(t_{i+1}\right),$$

(26)

where the effective stiffness matrix $K_{ef}$ and the effective load vector $R_{ef}\left(t_{i+1}\right)$ have the following forms:

$$K_{ef}=\frac{1}{a·∆t}M+\frac{\alpha }{∆t}D+K,$$

$$R_{ef}\left(t_{i+1}\right)=F\left(t_{i+1}\right)+Q_1\left(\overline{V}\left(t_i\right)+∆t\dot{\overline{V}}\left(t_i\right)+∆tb\ddot{\overline{V}}\left(t_i\right)\right)+$$

$$+Q_2\overline{V}\left(t_i\right)-\left(1-\alpha \right)D{\sum }_{j=1}^i\left(G\left(t_i,{\tau }_j\right)\dot{\overline{V}}\left({\tau }_j\right)∆\tau \right),$$

$Q_1=\frac{1}{a·∆t}M,Q_2=\frac{\alpha }{∆t}D,$ $i=1,2,3\dots n$ is the number of the time step, $j=1÷i$, $∆t=\left|t_{i+1}-t_i\right|$ is the time gap between the moments $t_i \mathrm{a}\mathrm{n}\mathrm{d} t_{i+1}$, and $∆\tau$ is considered equal to $∆t$.

Iterative schemes for calculating accelerations and velocities based on the assumption of the change in accelerations over a finite interval $[t_{i-1},t_{i+1}]$ according to a linear law utilizing the Newmark method are the following:

$$\ddot{\overline{V}}\left(t_{i+1}\right)=\frac{\overline{V}\left(t_{i+1}\right)-\overline{V}\left(t_i\right)}{2{\left(∆t\right)}^2}-\frac{\dot{\overline{V}}\left(t_{i-1}\right)}{2∆t}-\frac{a+b}{b+c}\left[\ddot{\overline{V}}\left(t_{i+1}\right)-\ddot{\overline{V}}\left(t_i\right)\right],$$

(27)

$$\dot{\overline{V}}\left(t_{i+1}\right)=\dot{\overline{V}}\left(t_{i-1}\right)+2∆t\ddot{\overline{V}}\left(t_{i-1}\right)+2∆t\frac{a+b}{b+c}\left[\ddot{\overline{V}}\left(t_{i+1}\right)-\ddot{\overline{V}}\left(t_i\right)\right],$$

(28)

where a is the interval between $t_{i-1}$ and $t_i$, $b$ is the interval between $t_i$ and $(t_i+\widehat{t}$), $c$ is the interval between $\left(t_i+\widehat{t}\right)$ and $t_{i+1}$, and, thus, $a+b+c=2∆t$ (Figure 2).

In Formulas (27) and (28), the characteristic time $\widehat{t }$ defines the influence coefficient ($\frac{a+b}{b+c}$). It could be chosen according to the minimal divergence between the calculated and experimental magnitudes of displacements. In the model discussed, it has been taken as a constant value with due account for the minimal error of calculations [69].

Since the influence kernel function is fading with time and after some time becomes insignificantly small, there is no need to remember the whole vibration process to determine the damping forces. For a rather long duration of the vibration processes, the memory length is limited using the mnemonic parameter according to the rule (Figure 3)

$$G\left(t_i-{\tau }_j\right)=\frac{2\mu }{\sqrt{\pi }}{e}^{-{\mu }^2{\left(t_i-{\tau }_j\right)}^2},$$

(29)

where $i=1,2,3,\dots , {\mathrm{T}}_n/\Delta t+1$is the number of the current discrete time step, $j=k,k+1,\dots ,i$ are numbers of all computational steps previous to $i$, ${\mathrm{T}}_n=\left|t_n-t_0\right|$ is the length of the time interval of forced mechanical vibrations until their entire damping, and $\Delta t={\mathrm{T}}_n/n$ is the step of discretization of the time interval $[t_0,t_n]$ divided into $n$segments. In so doing, $k=1,2,\dots ,{\mathrm{T}}_n/(\Delta t{M}_{\eta })$ is the number of the time step that limits the interval of memory history with the length equal to the magnitude of the mnemonic parameter $M_{\eta }$of nonlocal-in-time damping, which characterizes the endurance of material memory about its deformations.

8. Numerical Example and Results for the Beam Model

The mathematical model of beam vibrations considering the nonlocal damping-with-memory was implemented in the MATLAB software package. In order to make this model to be applicable for practical use, it should be calibrated according to the results of experimental data.

The method for determining the parameter μ, describing the level of damping nonlocality of the beam material, based on the data of the numerical experiment via the least squares method is presented in Reference [65]. This technique was used to determine the influence distance for the damping-with-memory model in References [69,70]. For this purpose, a three-dimensional finite element model of the considered beam element was realized within the verified calculation software SIMULIA Abaqus CAE. The obtained data were imported into the MATLAB software package. Determination of the three parameters $\mu$, $\alpha$, $M_{\eta }$, which control the nonlocal damping model, was implemented according to the rule of minimum error based on the least squares method (Figure 4):

$$\frac{100\mathrm{\%}}{N}·\sum _{i=1}^N\left|\frac{{\overline{V}}^{exp}\left(t_i\right)-{\overline{V}}^{model}(\alpha ,\mu ,M\eta ,t_i)}{{\overline{V}}^{exp}\left(t_i\right)}\right|=f^{error},$$

(30)

where $N_t$ = 251 is the number of time increments, ${\overline{V}}^{exp}\left(t\right)$ and ${\overline{V}}^{model}(\alpha ,\mu ,M_{\eta },t)$ are, respectively, the displacement vector components obtained via numerical experiments and due to the nonlocal-in-time model.

As an example, let us consider the dynamic behavior of a composite beam shown in Figure 5, which was made of orthotropic thermoset vinyl ester class 1 GFRP, the properties of which were determined experimentally in References [74,75,76] and are given in

Table 3. Properties of thermosetting vinyl ester fiberglass 1FRP [74,75,76].

Young’s modulus in the longitudinal direction Elong	17.2 GPa
Young’s modulus in the transverse direction Etrans	12.2 GPa
Poisson’s ratio in the longitudinal direction μlong	0.32
Poisson’s ratio in the transverse direction μtrans	0.15
density of the material ρ	1900 kg/m3
Damping coefficient (critical fraction)	0.042

. The beam is clamped at the ends and is subjected to instantly applied uniformly distributed load q = 10 kN/m.

The displacements of a middle span of the beam in time are shown in Figure 6. The solid line corresponds to the beam displacements obtained using the calibrated nonlocal model, and the dashed curve to those via a standard 3D model simulated in SIMULIA Abaqus (Figure 7).

In order to estimate the value of divergence between the results, the root-mean-square error ${\sigma }^{error}$ was calculated.

$${\sigma }^{error}=\sqrt{\frac{\sum _{i=1}^N{\left(y_i-f_i\left(\mu \right)\right)}^2}{N}},$$

(31)

where $y_i$ are the ordinates of the numerical experiment data, and $f_i\left(\mu \right)$ are the ordinates of the results via the nonlocal-in-time damping model.

The root-mean-square error was 5.9% for the deflection at 20 s.

A comparison of the displacements of the middle node obtained using the classical (local in time) one-dimensional damping model based on the Voigt hypothesis with those obtained in SIMULIA Abaqus is shown in Figure 8. In this case, the relative root-mean-square error was 16.5%.

Instead of central differences (explicit scheme), the Newmark method (implicit scheme) could be utilized [69,70]. In this case, the nonlocal model becomes more stable and larger time increments could be implemented. The accuracy of the obtained results is almost the same in comparison to the central differences, but the time of one simulation of the vibration process is shortened more than twenty times [69].

9. Modeling of Composite Frame Vibrations Considering Nonlocal-in-Time Damping Model

Let us show how the nonlocal-in-time damping model could be integrated into the FEA algorithm to study dynamic response of frame structures [77,78]. The equation of motion for the frame structure is similar to Equation (3), while the stiffness matrix for the 2D frame element has the following form [71]:

$$K_i=\left[\begin{array}{cccccc}\frac{EA}{l}& 0& 0& -\frac{EA}{l}& 0& 0\\ 0& \frac{12EI}{l^3}& \frac{6EI}{l^2}& 0& -\frac{12EI}{l^3}& \frac{6EI}{l^2}\\ 0& \frac{6EI}{l^2}& \frac{4EI}{l}& 0& -\frac{6EI}{l^2}& \frac{2EI}{l}\\ -\frac{EA}{l}& 0& 0& \frac{EA}{l}& 0& 0\\ 0& -\frac{12EI}{l^3}& -\frac{6EI}{l^2}& 0& \frac{12EI}{l^3}& -\frac{6EI}{l^2}\\ 0& \frac{6EI}{l^2}& \frac{2EI}{l}& 0& -\frac{6EI}{l^2}& \frac{4EI}{l}\end{array}\right].$$

(32)

The damping matrix could also be obtained according to the requirement of stationary state of the total energy of the vibrating system as was demonstrated earlier. The procedure for developing the damping matrix has been presented in detail in Reference [67]. Unlike the beam element, for the frame element, the axial displacement u should be taken into account. Then, energy dissipation could be written as follows:

$${\mathsf{\Phi }}_{Di}=\frac{1}{2}\underset{l}{\overset{}{{\displaystyle \int }}}\chi {\dot{\epsilon }}_a^{2}Adz+\frac{1}{2}\underset{l}{\overset{}{{\displaystyle \int }}}\underset{A}{\overset{}{{\displaystyle \int }}}\chi {\dot{\epsilon }}_b^{2}dAdz,$$

(33)

where A is the frame element cross-section area, and ${\dot{\epsilon }}_a$ is the tension-induced axial strain rate.

Tension-induced axial strain is as follows [73]:

$${\epsilon }_a=\frac{du}{dz}.$$

(34)

Then Equation (33) will take the form

$${\mathsf{\Phi }}_{Di}=\frac{1}{2}A\chi \underset{l}{\overset{}{{\displaystyle \int }}}{\left(\frac{d\dot{u}}{dz}\right)}^{2}dz+\frac{1}{2}I\chi \underset{l}{\overset{}{{\displaystyle \int }}}{\left(\frac{d^2\dot{v}}{d{z}^2}\right)}^2{}^{}dz.$$

(35)

Considering the linear shape function for the axial displacements $\left[N_u\right]=\left[\begin{array}{cc}1-\xi & \xi \end{array}\right]$ relationship (35) could be transformed to:

$${\mathsf{\Phi }}_{Di}=\frac{1}{2}A\chi \underset{0}{\overset{1}{{\displaystyle \int }}}{\left(A_u\left[N_u\right]\right)}^T{A}_u\left[N_u\right]{\dot{u}}_i^{2}ld\xi +\frac{1}{2}I\chi \underset{0}{\overset{1}{{\displaystyle \int }}}{\left(A_v\left[N_v\right]\right)}^T{A}_v\left[N_v\right]{\dot{v}}_i^{2}ld\xi .$$

(36)

where ${\dot{u}}_i=\left(\begin{array}{c}{\dot{u}}_0\\ {\dot{u}}_l\end{array}\right) \mathrm{is}$ the vector of nodal velocities of axial displacements, $A_u=\frac{d}{l·d\xi }$, and ${\dot{u}}_0$ and ${\dot{u}}_l$ are velocities of the initial and terminal points of the i-th element.

Therefore, the damping matrix of the frame element $D_i$ could be written as follows:

$$D_i=A\chi \underset{0}{\overset{1}{{\displaystyle \int }}}{\left(A_u\left[N_u\right]\right)}^T{A}_u\left[N_u\right]ld\xi +I\chi \underset{0}{\overset{1}{{\displaystyle \int }}}{\left(A_v\left[N_v\right]\right)}^T{A}_v\left[N_v\right]ld\xi ,$$

(37)

or

$$D_i=\chi \left[\begin{array}{cccccc}\frac{A}{l}& 0& 0& -\frac{A}{l}& 0& 0\\ 0& \frac{12I}{l^3}& \frac{6I}{l^2}& 0& -\frac{12I}{l^3}& \frac{6I}{l^2}\\ 0& \frac{6I}{l^2}& \frac{4I}{l}& 0& -\frac{6I}{l^2}& \frac{2I}{l}\\ -\frac{A}{l}& 0& 0& \frac{A}{l}& 0& 0\\ 0& -\frac{12I}{l^3}& -\frac{6I}{l^2}& 0& \frac{12I}{l^3}& -\frac{6I}{l^2}\\ 0& \frac{6I}{l^2}& \frac{2I}{l}& 0& -\frac{6I}{l^2}& \frac{4I}{l}\end{array}\right].$$

(38)

In order to solve the equation of frame motion, the nonmodified Newmark method (an implicit scheme) could be used [71]. In this case, the first- and second-order time derivatives of the displacement vector $V\left(t\right)$ (velocity and acceleration) entered into Equations (3) and (18) could be presented as follows:

$${\dot{V}}_{i+1}=\frac{V_{i+1}-V_i}{\Delta t}, {\ddot{V}}_{i+1}=\frac{2}{\Delta t^2}\left(V_{i+1}-V_i-{\dot{V}}_i\Delta t\right)-{\ddot{V}}_i.$$

(39)

Considering Relationships (39), Equation (18) is reduced to the following form:

$$M\left[\frac{2}{\Delta t^2}\left(V_{i+1}-V_i-{\dot{V}}_i\Delta t\right)-{\ddot{V}}_i\right]+D\frac{1}{\Delta t}{\displaystyle \sum }_{j=1}^{i+1}G\left(t_{i+1},{\tau }_j\right){\dot{V}}_j+K{V}_i=F_i.$$

(40)

Here $G\left(t_{i+1},{\tau }_j\right)$ is the discrete analog of the $G\left(t-\tau \right)$ kernel, which, for Error Function (20), is calculated as follows:

$$G\left(t_{i+1},{\tau }_j\right)=\Delta t\left[\frac{2\mu }{\sqrt{\pi }} e^{-{\mu }^2{\left(t-\tau +\frac{\Delta t}{2}\right)}^2}\right],$$

(41)

where $t=\left(i+1\right)\Delta t$, $\tau =j\Delta t$, and $j=1,2,\dots ,i+1$ is the number of the time step, at which the displacement vector is calculated.

The Newmark scheme is implicit, so one needs to know ${\dot{V}}_{i+1}$ to calculate the displacement vector. To solve this problem, the damping-with-memory term in Equation (40) should be divided into two parts, resulting in

$$M\left[\frac{2}{\Delta t^2}\left(V_{i+1}-V_i-{\dot{V}}_i\Delta t\right)-{\ddot{V}}_i\right]+\beta D+\alpha D+K{V}_i=F_i,$$

(42)

where one of two memory parts, $\alpha$, is related to the i + 1 step, when the other, $\beta$, corresponds to all previous steps:

$$\begin{gathered} \beta ={\displaystyle \sum }_{j=1}^i\frac{2\mu }{\sqrt{\pi }}{e}^{-{\mu }^2{\left(t_j-t_{j-1}\right)}^2}\left(V_j-V_{j-1}\right), \\ \alpha =\frac{2\mu }{\sqrt{\pi }}{e}^{-{\mu }^2{\left(t_{i+1}-t_i\right)}^2}\left(V_{i+1}-V_i\right). \end{gathered}$$

(43)

Equation (42) could be transformed into the computational scheme for the step-by-step calculation of $V_{i+1}$using the vectors $V_i$ and $V_{i-1}$, which are calculated on the previous increments $i$ and $i-1$

$$Z{V}_{i+1}=F_{i+1}+M{\ddot{V}}_i+Q_1{\dot{V}}_i+Q_2{V}_i+\beta D,$$

(44)

where

$$\begin{gathered} Z=\left(\frac{2}{\Delta t^2}M+\frac{1}{\Delta t}\alpha D+K\right), Q_1=\frac{2}{\Delta t}M, \\ {Q}_2=\frac{1}{\Delta t}\left(\frac{2}{\Delta t}M+\alpha D\right). \end{gathered}$$

(45)

10. Numerical Results for the Frame Model

As a numerical example, the finite element model of a vibrating frame was implemented in the MATLAB software package. The frame geometry, boundary conditions, and loading are shown in Figure 9. The load is applied instantly.

The frame is made of the same material as the beam structure considered before.

For the first step of the simulated vibration process $i=1$, we assume $V_1={\dot{V}}_1={\ddot{V}}_1=0$ as the initial conditions.

To determine the scale parameter for the considered structure, once again the computer simulation results obtained in SIMULIA Abaqus were used. The frame was modeled utilizing 3D-solid finite elements considering the orthotropic properties of the material. Using the least squares method, the magnitude of the parameter was calculated as $\mu =0.22$, which points to a significant level of the nonlocal properties of the material. The deflection time history simulated by the calibrated nonlocal model in comparison with the results of 3D modeling via Abaqus is shown in Figure 10.

Even though the orthotropic properties of the material were not taken into account in the 2D-nonlocal-damping model of the frame, the results match the numerical experiment with satisfying accuracy. This is due to the fact that the nonlocal model could be flexibly controlled by its scale parameter $\mu$. The difference between the results provided by the classical local-damping model and Abaqus data is shown in Figure 11.

11. Conclusions

From the review of damping models presented in this paper, it is evident that in comparison to local damping models, the nonlocal models allow one to obtain the main characteristics of the simulated vibration process in a more reliable and flexible way. The proposed models enable one to consider both internal and external damping at a time.

The enlarged flexibility makes it possible to use one-dimensional models of beam elements in the dynamic analysis of structures made of modern composite materials with orthotropic properties. The damping model with memory, which could be calibrated using the least-squares method, allows one to approximate the data of a numerical experiment with satisfactory accuracy.

The so-called nonlocal-in-time damping models are more easily integrated into the finite element method algorithm as compared with the nonlocal-by-coordinate damping models, what makes its application quite effective for solving different engineering problems.