This work considers a cylindrical rod, acting as a nonlinear viscoelastic spring, with an attached mass that can displace along and rotate about the rod’s axis, thereby inducing an extension and twist in the rod. The coupling between them produces the Poynting effect in the body. The material is modeled by the nonlinear single integral Pipkin–Rogers constitutive equation. The mass has been assumed to be at rest in a long-time equilibrium state and then given a small axial or rotational disturbance and released (plucked). The subsequent motion is governed by a linear Volterra integro-differential equation. For a particular choice of material parameters, analytical expressions for the time-dependent decay of the disturbance are obtained in terms of the extension/torsion coupling, material stiffness, the amount and rate of stress relaxation, and the inertia of the mass. For other choices of material parameters, the results provide insight into how these quantities affect the time-dependent decay. A number of plucking scenarios are treated such as extensional or rotational disturbances from the undeformed state, rotational plucking after a finite axial stretch, and extensional plucking after a finite rotation. Simultaneous rotational and extensional plucking leads to a system of Volterra integro-differential equations whose treatment is deferred to later work.
When a circular cylinder undergoes simultaneous stretching and twisting, each material element undergoes a shear superposed on an elongation. The elongation and the shear both contribute to the normal stress on the cross-section, resulting in a coupling between the angle of twist, axial extension ratio, twisting moment, and axial force on the cylinder. This effect was observed by Poynting who showed that if no compressive axial force was applied to a twisted steel cylinder, it would elongate (the Poynting effect [1]). Freudenthal and Ronay [2] discussed this effect in experiments involving torsion of rods of various materials, both time independent (steel, copper, aluminum, rubber) and time dependent (polyurethane).
Within the context of nonlinear elasticity, Rivlin and Saunders [3] presented results of an experimental study of the Poynting effect in vulcanized rubber cylinders bonded to steel end plates. They showed that when a cylinder was twisted, a compressive axial force was required to keep it at zero elongation. A simple qualitative picture for the origin of the compressive normal stress in a shear deformation was given by Treloar [4] (see also Gent [5]). Treloar’s explanation suggests that the Poynting effect would be observed in experiments on torsion of cylinders for different materials and kinds of mechanical responses. Min et al. [6] derived results based on the Poynting effect for incompressible isotropic nonlinear viscoelastic solids while Goldberg and Lianis [7] demonstrated the effect in an experimental extension/twist program on an SBR cylinder. When extension/torsion experiments on rubber bars are conducted in a thermal environment at temperatures above, say , the rubber will undergo significant time-dependent changes in its macromolecular structure due to continuous scission and re-cross-linking of macromolecular network junctions. The Poynting effect for rubber undergoing these microstructural changes was discussed in Wineman [8–10].
In each of the aforementioned studies, whether experimental or theoretical, the axial force and twisting moment were applied to the cylinder by a hard loading device. Dynamical loading conditions were treated in Beatty and Bhattacharyya [11] (and also the references therein). One end of the cylinder was attached to a fixed support and the other end was attached to a rigid cylindrical mass that could move axially along and rotate about the cylindrical axis. This system, coupling the axial and rotational motions of the mass due to the Poynting effect in the cylinder, was referred to in Beatty and Bhattacharyya [11] as a Poynting oscillator. For an undamped, elastic (neo-Hookean) cylinder, Beatty and Bhattacharyya [11] presented analytical results that show how the oscillatory axial and/or rotational motion of the mass depends on system parameters as well as pre-stretch and/or pre-twist.
The motion of the mass when it is attached to a viscoelastic cylinder is also of interest. In Fosdick [12], attention was restricted to axial motions when the cylinder is modeled by the constitutive equation for finite linear viscoelasticity [13]. A special case was considered which allowed the governing equation to be reduced to a system of linear ordinary differential equations that were amenable to analysis. The stress relaxation function was a one-term exponential. This introduced a relaxation time and a coefficient indicating the amount of stress relaxation. A detailed discussion examined the influence of these parameters on two aspects of the motion: resonance and whether a small disturbance died out as underdamped or overdamped motion. In Fosdick et al. [14], nonlinear oscillatory aspects of the motion were considered. Some aspects of these results can also be found in Wineman and Rajagopal [15] for a mass attached to a standard linear viscoelastic rod.
The interest of the present work is in examining cases when viscoelastic response and the Poynting effect in the cylinder combine with the equations of motion for both axial and rotational motions of an attached mass. To this end, a nonlinear single integral constitutive equation is used that differs from that in Fosdick et al. [12]. The instantaneous initial and long-time equilibrium responses are like that of a Mooney–Rivlin material, but with different moduli at and in the limit as . The rotational/axial motion of the mass is governed by nonlinear Volterra integro-differential equations that do not appear to be convertible to a system of differential equations.
The particular case that is considered is as follows. The mass is first given an axial and/or rotational displacement and then held fixed until the applied loads on the cylinder have relaxed to their long-time equilibrium values. The mass is then given an axial and/or rotational disturbance and released, i.e., it is plucked like a violin string. The purpose is to determine how the time-dependent decay of the disturbance depends on the viscoelastic properties and the underlying deformation.
Some motivation for this work comes from the use of elastomeric elements in soft robotic devices. As these elements are combined with movable stiffer components, their nonlinear viscoelastic properties can be important factors in the design and control of such devices. The Poynting oscillator has a number of features that make it useful in providing insight into these design and control issues. First, the deformation of the cylinder is known a priori and has been shown to be the same for every isotropic incompressible solid material [16–18]. Thus, the study can use any such material model. There are distinct modes of deformation that influence each other, i.e., axial and torsional. Finally, its response depends on the nonlinear viscoelastic properties, modal interaction, and inertia of the attached mass.
The constitutive equation is introduced in section 2 along with a particular choice for the material parameters. The governing equations for the axial and torsional motion are formulated in section 3. In section 4, the mass is given a small axial or rotational disturbance from the undeformed state and released. For particular material parameters, each disturbance is shown to satisfy a linear Volterra integro-differential equation. Solutions are found for the axial and rotational disturbances that show they decay differently with time, both in their time course and whether the motion is overdamped or underdamped. A small rotational plucking superposed on a finite axial displacement is treated in section 5. The disturbance satisfies a linear Volterra integro-differential equation that depends on the base finite extension. Its solution shows how the underlying axial displacement affects the time decay of the rotational disturbance. Section 6 contains an analogous treatment for a small axial plucking superposed on a finite rotation. Concluding comments are in section 7.
2. Constitutive equation
Let denote the current time and denote an earlier time, i.e., . The configuration of a body at time is denoted by . The body is assumed to be in the same configuration, denoted by , for times . This is taken as its reference configuration and is assumed to be a natural or stress free configuration. A material particle at position in is connected to its position in by a relation . Since the body is initially in , for . The deformation gradient is , with , .
The material considered here is assumed to be incompressible which requires that , . The right Cauchy–Green tensor at time is , with invariants and . For notational convenience, this set of invariants is denoted by .
The constitutive equation for an incompressible nonlinear viscoelastic solid has the form:
where is the Cauchy stress, is a Lagrange multiplier arising from the incompressibility constraint (interpreted as a pressure), and the notation indicates dependence on the history of . Nonlinear single integral models generally have the form:
The particular form used here is that introduced by Pipkin and Rogers [19]:
is a tensor-valued deformation-dependent relaxation property. At fixed , monotonically decreases from to a nonzero limit denoted by . For an isotropic material:
and the scalar coefficients are the deformation-dependent relaxation properties.
The scalar functions in equation (4) are selected to have the form:
where and are the stress relaxation functions.
The motivation for these specific forms is seen when they are combined with equations (1), (3), and (4). Evaluation at gives the stress-deformation relation for instantaneous response:
where use has been made of the definition and the Cayley–Hamilton theorem. The term combines in equation (1) with all terms of the form that arise in the calculation. Assuming the material approaches an equilibrium state in the limit as , a similar approach gives the stress-deformation relation for long-time equilibrium response:
where . Thus, the stress-deformation relations at and in the limit as have the form of that for a nonlinear elastic Mooney–Rivlin material, but with different sets of material parameters.
The relaxation functions and can have any form. For the purpose of the numerical examples in later sections, the relaxation functions are:
where and are the characteristic relaxation times. The form of assures that the constitutive equation reduces to that for incompressible isotropic linearized viscoelasticity in the limit of small deformation histories.
3. Formulation
Consider a solid cylindrical rod with radius and length in its natural or stress free configuration. One end is attached to a fixed support and the other is attached to a rigid cylinder with mass and moment of inertia about its central axis. The rigid body undergoes a displacement along and rotation about the cylindrical axis of the rod while the rod undergoes a combined extension/torsion deformation history. An axial force and twisting moment act on each end of the rod while its cylindrical surface is free of traction. It is assumed that the inertia of the rod is small compared to that of the rigid body.
3.1. Kinematics
The motion is described with respect to a cylindrical coordinate system whose -axis coincides with the cylindrical axis of the rod with its origin at the fixed end. The motion connecting the coordinates of a particle in its reference configuration , , to its coordinates in its configuration at time , , has the form:
where is the angle of twist per current length and is the axial stretch ratio. If denotes the twist angle of the rod or the rotation angle of the mass, then .
Since , dependent variables can be regarded as functions of or , whichever is convenient for the development. The components of the deformation gradient with respect to an orthonormal basis oriented along cylindrical coordinate lines are:
Note that as required by the constraint of incompressibility.
The components of the are:
is readily calculated but is omitted for the purpose of brevity. The invariants are:
3.2. General expressions for stresses, axial force, and twisting moment
Note that , and hence by equation (3), . Then, denoting :
where all quantities are evaluated at .
Since depends only on and , or equivalently , a straightforward analysis of the linear momentum equations at time , neglecting inertia, shows that the same holds for , and:
The resultant twisting moment on the rod at is:
With the use of the first of equation (9), this is transformed into an integral over the cross-section in the reference configuration:
The torque–twist–stretch relation for the rod is obtained by combining the shear stress in equation (20), the integral for the twisting moment (16), and the specific choice of material properties in equation (5). On recalling that and introducing :
3.5. Axial force
The axial force–twist–stretch relation for the rod is similarly obtained by substituting the expressions (21) and (22) into equation (19):
where:
and
3.6. Equations of motion for the mass
The equation for the axial motion of the mass, as in Beatty and Bhattacharyya [11], is:
where is the current length of the rod, and is an axial external force applied to the mass. The equation for the rotational motion of the mass is:
where is an external moment about the axis of the cylinder applied to the mass.
These Equations form a system of coupled Volterra integro-differential equations for and , the axial and rotational motions of the mass.
4. Small disturbances from the reference configuration
Let the mass be subjected to small axial and rotational disturbances from the reference configuration, i.e., if , then and . Beatty and Bhattacharyya [11] discussed this motion in detail in the absence of damping, and Fosdick et al. [12] discussed the case with damping for only axial motion. The interest here is the influence of viscoelasticity on the decay with time of both axial and rotational disturbances from the reference configuration when they are imposed and then released, i.e., motions for which initial values are specified and in equation (31) and in equation (31).
Let in the torsion equation (24) and in the axial force equation (25) with equations (26)–(28) and then linearize the equations using and . The equations of motion for the mass uncouple and each has the form:
where is given in terms of the relaxation functions in equation (5) by . For the axial motion:
and for the rotational motion:
As the physical dimension of is , characteristic axial and rotational inertial times can be defined by and . Since the moment of inertia of a cylindrical disc of radius is :
If , then or . This suggests that the damping of rotational disturbances differs from the damping of axial disturbances.
in which are the roots of a quintic polynomial equation and depend on . As it is impractical to explore the dependence of on this parameter space, let , so that now in equations (32) and (33). The solution of equation (31) then has the form:
in which are determined from the cubic polynomial equation:
The initial conditions provide two other equations for .
The values depend on only two parameters, and . Parameter is the characteristic stress relaxation time relative to the characteristic inertial time and is the amount of stress relaxation. They account for the influence of the fundamental features of stress relaxation on the time-dependent decay of a disturbance.
The formulae for the roots of a cubic polynomial equation are well documented (e.g., Abramowitz and Stegun [20]). The roots depend on the parameters in equation (37) through the quantities:
The roots are given by:
If , are real, , in equation (36), and the motion is overdamped. If , is real, , and , are complex conjugates with . The motion (36) can be written as:
and the disturbance is underdamped.
Figure 1 shows the points where . in the small enclosed region where and . For the purpose of a numerical example, let , for which . denotes the point where and denotes the point directly below at where . Hence, both the axial and rotational disturbances decay as underdamped motions. denotes the point where and denotes the point where . Now, the axial disturbance decays as an underdamped motion while the rotational disturbance decays as an overdamped motion. For other choices of the parameters , it may be the case that both disturbances are underdamped, both are overdamped, or one is overdamped and the other is underdamped. In general, even though both the axial and rotational motions may be underdamped or overdamped, their time courses differ, as shown in the following numerical example. For the axial disturbance corresponding to , in equation (41), while for the rotational disturbance corresponding to , . For completeness in presentation, for the axial disturbance , and for the rotational disturbance , .
Regions showing where and . Points correspond to axial motion and points correspond to rotational motion.
The next several figures show the complicated dependence of in equation (41) on stress relaxation, as represented by and . For this purpose, attention is limited to the region of Figure 1 where . Figure 2 shows that decreases rapidly with and then levels off to a limiting value that increases with . In contrast, Figure 3 shows that increases rapidly with and then levels off to a limiting value that decreases with . , shown in Figure 4, increases monotonically with and then becomes independent of .
versus for .
versus for .
versus for .
5. Decay of rotational disturbance due to combined axial pre-stretch and viscoelasticity
Let the mass be axially displaced a specified amount without being rotated. While at a fixed axial displacement, a small amplitude rotational disturbance is applied and released, i.e., the mass is rotationally plucked. This imposes a combined extension/twisting deformation on the cylinder. The influence of the axial displacement on the time-dependent decay of the rotational disturbance arises from the terms in equation (24) that couple the twist and the axial stretch.
Thus, as the mass displaces axially during a time interval , the axial stretch ratio of the cylinder varies from to . The mass is then held at a fixed axial displacement, i.e., , . Meanwhile, for , , and then varies so that , . Time is large enough that the stress required to maintain the fixed stretch has reached a long-time equilibrium value.
The stress required to produce the axial stretch history is found using the Pipkin–Rogers constitutive equation in equations (1), (3), and (4) and the material properties in equations (5) and (8). The details of the calculation are omitted in the interests of brevity of presentation. The long-time equilibrium stress is:
the same relation as for a Mooney–Rivlin material. In developing this expression, the transient stress differs from its long-time equilibrium value by terms of the form:
where , so that for . For each fixed , as increases regardless of the form of the relaxation function. For example, if is one of the functions in equation (8), the integral in equation (43) has the form , , with being a constant and are the relaxation times in equation (8). As another example, consider materials with a finite memory of duration , i.e., for . The integral in equation (43) becomes zero when . Thus, can be chosen so that terms of the form equation (43) can be neglected compared to the long-time equilibrium term.
Next, on letting in equation (20), using for notational convenience, and linearizing using :
where denotes and denotes .
On replacing the general properties by the expressions in equation (5) and simplifying the group of terms in equation (44) that are in the braces and multiplied by become:
with:
The expression in equation (44) that is in the braces and multiplied by contains terms of the form:
By a discussion similar to that associated with equation (43), for times , such terms become very small and the expression in the braces has the limit:
Let in equation (45) and introduce the change of variable . Finally, replace by so that time is now measured from . Equation (44) becomes:
where:
is a stress relaxation function that depends on the stress relaxation properties for linear viscoelastic response and the axial pre-stretch.
Figure 5 shows the sets of points where for . in the small enclosed region where . The shape of this region varies with the axial pre-stretch . The axial pre-stretch also affects the time dependence of a disturbance. This is illustrated in Figure 6 for an underdamped motion for and . increases while and decrease with .
Regions showing where and for .
versus .
6. Small axial and rotational disturbances from a pre-stretched and pre-twisted state
Let the mass be axially displaced and rotated to some state and then be held in that state while the applied axial force and torque relax to long- time equilibrium values. The nonlinear viscoelastic cylinder has a corresponding combined elongation/twisting deformation history from its reference configuration to a deformed configuration . The influence of the initial axial and rotational displacements of the mass on the time-dependent decay of a small amplitude disturbance is to be characterized.
The rotation history of the mass is , where is a known base history and is a superposed disturbance. During the time interval , varies from to and for later times, , . The disturbance for where . For , , where and . Time is large enough that the applied loads required to maintain the fixed deformation have reached their long-time equilibrium values.
The axial position history of the mass is , where is a known base history and is a superposed disturbance. During the time interval , , varies from to and, for later times , . The disturbance for where . For , where and . The following expressions are useful: for the stretch ratios—, and for ; for the twist—, and for .
The twisting moment and axial force histories on the cylinder that are associated with the specified twisting and elongation histories during are calculated by substituting and into equations (24) and (25)–(28). By a discussion similar to that related to equation (43), integrals over the transient portions of and , i.e., of the form:
decrease to zero as increases. The integrals over the time interval when the twist and stretch are held constant have the form:
Time at the onset of a disturbance is assumed to be sufficiently large that integrals of the form in equation (58) can be neglected and that is close to . Thus, the twisting moment and axial force have reached their long-time equilibrium expressions at and are given by:
Now, consider times . Let the disturbed axial stretch and twist be substituted into equations (24) and (25)–(28) and the resulting expressions be linearized. Then, introduce the notation , , . With the resulting expressions for the twisting moment and axial force in terms of the disturbances, the equations of motion for the rotational and axial motion of the mass, (29) and (30), become:
where is a linear Volterra integral operator on function . and are known from equations (60) and (61), and it is assumed that . Specific expressions for are:
Equations (62) and (63) are a system of linear Volterra integro-differential equations for , . As in previous sections, let so that the systems (62) and (63) simplify to:
and
Starting with equation (69), can be expressed directly in terms of and then substituted into equation (70) to give a higher-order linear Volterra integro-differential equation for . In a similar manner, starting with equation (70) leads to a corresponding equation for .
The solution of such equations is beyond the scope of this work. Instead, in the spirit of previous sections, consider the special case when the rod undergoes the time-dependent twisting described above while the axial force . As the twist history appears in the expressions for in equation (27) and in equation (28), these will not be zero. Accordingly, as seen from equation (25), there must be an axial stretch history so that . A numerical example illustrating this time-dependent version of the Poynting effect was given in Wineman [9] for a rod whose material is undergoing time-dependent microstructural change. In the present context, because of the Poynting effect, as the twist increases from zero to , the axial stretch ratio will increase from one to .
For the twisting and axial stretch histories described above that have specific limits and , relations (60) and (61) connect these limits with the long-time equilibrium twisting moment and axial force. When , equation (61) implies that . Consider the special case when the rotation of the mass is undisturbed, so that , but there is an axial disturbance. This case when affects the time decay of an axial disturbance is the complementary situation to that treated in section 5.
In the equation of axial motion for the mass (29), and the axial displacement is now . The axial force is given by equation (63) with and simplified by setting . The resulting equation is:
By applying the Laplace transform, or by direct substitution of the assumed representation (36) into equation (71), it can be shown that has the form in equation (54) if the constants satisfy the cubic equation:
where:
In the interests of brevity of presentation, it is sufficient to point out that the essential features of the discussion in section 5 can be applied here. In particular, since and , is present in the above expressions and it affects how fast the axial disturbance decays and whether it is overdamped or underdamped.
7. Concluding comments
This work has been concerned with how a nonlinear viscoelastic cylindrical rod affects the motion of an attached mass. In the particular situation considered, the mass is axially or rotationally displaced and then held at rest in a long-time equilibrium state. This causes an extension and torsion in the rod and their nonlinear coupling produces the Poynting effect. The mass is then “plucked” extensionally or rotationally, i.e., given a small axial or rotational disturbance and then released. Expressions for the subsequent time-dependent decay of the disturbance due to viscoelasticity have been developed that account for the Poynting effect, material stiffness, the amount and rate of stress relaxation, and the inertia of the mass. A number of different scenarios have been treated such as extensional or rotational plucking from the undeformed state, rotational plucking after a finite axial stretch, and extensional plucking after a finite rotation.
The analysis has been carried out using the Pipkin–Rogers nonlinear single integral constitutive framework for nonlinear viscoelasticity. It has been chosen because it contains the important features of nonlinear viscoelasticity, yet is analytically convenient. Specific material properties have been assumed so that the instantaneous and long-time equilibrium elastic responses have the same form as the Mooney–Rivlin constitutive for nonlinear elasticity. As a result, the time dependence of the disturbance is determined by a linear Volterra integro-differential equation. Under the further assumption that the initial and long-time equilibrium responses act as neo-Hookean materials, each disturbed motion is shown to be a sum of three exponentials with explicit expressions for their exponents as functions of the underlying stretch or twist, material stiffness, the amount and rate of stress relaxation, and the inertia of the mass.
The results presented here are limited by a number of simplifying assumptions, However, they may be of use in giving insight into the behavior of a disturbance when there is added complexity. When , inverting the Laplace transform will involve finding the roots of higher-degree polynomials for which analytical expressions are not available. The roots will need to be determined numerically. The disturbance will be a sum of exponential terms, but their exponents will not be explicit functions of the parameters as before. A method presented by Lazaro [21] may be useful for determining the boundary between underdamped and overdamped behavior. Another complexity occurs when there is both extensional and rotational plucking. The disturbances are governed by a system of linear Volterra integro-differential equations. The Laplace transforms for the individual disturbances will involve high-degree polynomials, and the above comments will again apply. It was assumed that the disturbances were applied at a time large enough that time-dependent quantities had reached their long-time equilibrium values. In some applications, this assumption may not be practical. The presence of such time-dependent terms may require numerical methods of solution. However, the results presented in the simpler cases here may provide useful guidance for considering the influence of various parameters.
Footnotes
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.
ORCID iD
Alan Stuart Wineman
References
1.
PoyntingJH.On changes in the dimensions of a steel wire when twisted, and on the pressure of distortional waves in steel. Proc R Soc London Ser A1912; 86: 534–561.
2.
FreudenthalAMRonayM.Second order effects in dissipative media. Proc R Soc London Ser A1966; 292: 14–50.
3.
RivlinRSSaundersDW. Large elastic deformations of isotropic materials. VII. Experiments on the deformation of rubber. Philos Trans Royal Soc Lond Ser A1951; A243: 251–280.
4.
TreloarLRG.The physics of rubber elasticity. 3rd ed.Oxford: Clarendon Press, 1975.
5.
GentAN.Article 3: elasticity. In: GentAN (ed.) Engineering with rubber: how to design rubber components. 2nd ed.Munich: Hanser Publications, 2001, pp. 37–88.
6.
MinBKKolskyHPipkinAC.Viscoelastic response to small deformations superposed on a large stretch. Int J Solids Struct1977; 13: 771–781.
7.
GoldbergWLianisG.Stress relaxation in combined torsion-tension. J Appl Mech1970; 37: 53–60.
8.
WinemanA.A universal relation for a chemomechanically changing rubber cylinder. Mech Res Commun2018; 94: 28–32.
9.
WinemanA.The Poynting effect in elastomeric bars undergoing chemo-mechanical evolution. Math Mech Solids2020; 26: 683–693.
10.
WinemanA.Combined extension/torsion experiments for assessing a constitutive theory for elastomers undergoing thermally induced microstructural changes. J Elasticity. Epub ahead of print 14October2022. DOI: 10.1007/s10659-022-09944-x.
11.
BeattyMBhattacharyyaR.Poynting oscillations of a rigid disk supported by a neo-Hookean shaft. J Elasticity1990; 24: 135–186.
12.
FosdickRKetemaYYuJ-H.Vibration damping through the use of materials with memory. Int J Solids Struct1998; 35: 403–420.
13.
ColemanBNollW. Foundations of linear viscoelasticity. Rev Modern Phys 196; 33: 239–249.
14.
FosdickRKetemaYYuJ-H.A non-linear oscillator with history dependent force. Int J Nonlin Mech1998; 33: 447–459.
15.
WinemanASRajagopalKR.Mechanical response of polymers: an introduction. Cambridge: Cambridge University Press, 2000.
16.
CarrollMM.Controllable deformations of incompressible simple materials. Int J Eng Sci1967; 5: 515–525.
17.
CarrollMM. Finite deformations of incompressible simple solids. I. Isotropic solids. Quart J Mech Appl Math1968; 21: 147–170.
18.
FosdickRLDynamically possible motions of incompressible, isotropic simple materials. Arch Ration Mech Anal1968; 29: 272–288.
19.
PipkinACRogersTG.A non-linear integral representation for viscoelastic behaviour. J Mech Phys Solids1968; 16: 59–72.
20.
AbramowitzMStegunIA.Handbook of mathematical functions (National Bureau of Standards applied mathematics series). Washington, DC: US Government Printing Office, 1965.
21.
LazaroM.Exact determination of critical damping in multiple-exponent kernel-based viscoelastic single degree-of-freedom systems. Math Mech Solids2019; 24: 3843–3861.