**1. Introduction**

26 Viscoelasticity

[39] Roy, A., Morozov, A., van Saarloos, W., & Larson, R.G. (2006). Mechanism of polymer drag reduction using a low-dimensional model, *Physical Review Letter*, Vol. 97, 234501,

[40] Schmid, P.J. & Henningson, D.S. (2001). *Stability and Transition in Shear Flows*,

[41] Sureshkumar, R. & Beris, A.N. (1995). Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows,

[42] Suzuki, H., Ishihara, K., & Usui, H. (2001). Numerical study on a drag reducing flow with surfactant additives, *Proceedings of 3rd Pacific Rim Conference on Rheology*, Paper No. 019,

[43] Takeuchi, H. (2012). Demonstration test of energy conservation of central air conditioning system at the Sapporo City Office Building., *Synthesiology, English edition*, Vol. 4, 136–143. [44] den Toonder, J.M.J., Hulsen, M.A., Kuiken, G.D.C., & Nieuwstadt, F.T.M. (1997). Drag reduction by polymer additives in a turbulent pipe flow: numerical and laboratory

[45] Tsukahara, T., Ishigami, T., Yu, B., & Kawaguchi, Y. (2011a). DNS study on viscoelastic effect in drag-reduced turbulent channel flow, *Journal of Turbulence*, Vol. 12, No. 13, 13 pp. [46] Tsukahara, T., Kawase, T., & Kawaguchi, Y. (2011b). DNS of viscoelastic turbulent channel flow with rectangular orifice at low Reynolds number, *International Journal of*

[47] Tsukahara, T., Kawase, T., & Kawaguchi, Y. (2011c). Numerical investigation of viscoelastic effects on turbulent flow past rectangular orifice, In: *Proceedings of the 22nd International Symposium on Transport Phenomena*, Paper #129 (USB), 7pp., Delft, The

[48] Tsukahara, T. & Kawaguchi, Y. (2011). Turbulent heat transfer in drag-reducing channel flow of viscoelastic fluid, *Evaporation, Condensation and Heat transfer* (ed., A. Ahsan),

[49] Tsukahara, T. & Kawaguchi, Y. (2012). DNS on turbulent heat transfer of viscoelastic fluid flow in a plane channel with transverse rectangular orifices, *Progress in Computational*

[50] Tsurumi, D., Kawada, S., Kawase, T., Tsukahara, T., & Kawaguchi, Y. (2012). Experimental analysis of turbulent structure of viscoelastic fluid flow in downstream of two-dimensional orifice. In: *Turbulent, Heat and Mass Transfer 7*, Begell House, Inc., in

[51] White, C.M. & Mungal, M.G. (2008). Mechanics and prediction of turbulent drag reduction with polymer additives, *Annual Review Fluid Mechacnics*, Vol. 40, 235–256. [52] Yu, B. & Kawaguchi, Y. (2004). Direct numerical simulation of viscoelastic drag-reducing flow: a faithful finite difference method, *Journal of Non-Newtonian Fluid Mechanics*, Vol.

[53] Yu, B., Li, F., & Kawaguchi, Y. (2004). Numerical and experimental investigation of turbulent characteristics in a drag-reducing flow with surfactant additives, *International*

[54] Zakin, J.L., Myska, J., & Chara, Z. (1996). New limiting drag reduction and velocity profile asymptotes for nonpolymeric additives systems, *AIChE Journal*, Vol. 42,

Springer-Verlag, ISBN 978-0-387-98985-3, New York.

*Journal of Non-Newtonian Fluid Mechanics*, Vol. 60, 53–80.

experiments, *Journal of Fluid Mechanics*, Vol. 337, 193–231.

Vancouver, Canada, 8–13 July, 2001.

*Heat and Fluid Flow*, Vol. 32, 529–538.

Netherlands, 8–11 November 2011.

InTech, Rijeka, Croatia, pp. 375-400.

*Journal of Heat and Fluid Flow*, Vol. 25, 961–974.

*Fluid Dynamics*, in press.

press.

116, 431–466.

3544–3546.

4pp.

Viscoelasiticity refers to the phenomenon in which a material body, when deformed exhibits both elastic (akin to solids)and viscous (akin to liquids) behavior. The body stores mechanical energy (elastic behavior) and dissipates it simultaneously(viscous behavior). Linear theory of viscoelasticity treats the body as a linear system which when subjected to an excitation responds with a response function. If the excitation is a stress, the response is a strain and if the excitation is a strain, the response is a stress. Mechanical models involving a spring-mass connected to a dashpot have been used to explain the elastic and viscous behavior.The mathematical structure of the theory and the spring-dash-pot type of mechanical models used and the so called Standard Linear Solid have all been only too well known [1-4]. In recent years methods of fractional calculus have been applied to develop viscoelastic models especially by Caputo and Mainardi [5,6] , Glockle and Nonenmacher [7], and Gorenflo and Mainardi [8]. A recent monograph by Mainardi[9] gives extensive list of references to the literature connecting fractional calculus, linear viscoelasticity and wave motion.All these works treat the phenomenon of viscoelasticity as a macroscopic phenomenon exhibited by matter treated as an elastic continuum albeit including a viscous aspect as well. It should be recognized that matter has an atomic structure and is fundamentally discrete in nature. A microscopic approach would recognize this aspect and a theoretical model would yield the results for the continuum as a limit.

It is well established that lattice dynamics provides a microscopic basis for understanding a host of phenomena in condensed matter physics, including mechanical, thermal, dielectric and optical phenomena, which are macroscopic, generally described from the continuum point of view. Since the pioneering work of Born and Von Karman [10], lattice dynamics has developed into a veritable branch of condensed matter physics. Lucid treatment of various

© 2012 Achar and Hanneken, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

aspects of lattice dynamical theory of condensed matter can be found in renowned treatises [11-13]. Of particular interest to the present paper is the work of Askar [14] on the lattice dynamical foundations of continuum theories of elasticity, piezoelectricity, viscoelasticity and plasticity.

Microscopic Formulation of Fractional Theory of Viscoelasticity 61

(3)

(5)

(6)

(7)

(8)

*s e* = + (2)

*s e <sup>s</sup>* + = (4)

() () *<sup>d</sup> t mt b*

<sup>1</sup> / () 1 (a)

*e t*

*<sup>t</sup> Jt e m Gt m b t*

= +

where *b m*/ *<sup>e</sup> t* = is referred to as the retardation time.

characterized by the constitutive relation

where *a <sup>s</sup> t* = is referred to as the relaxation time.

() () (b)

A model consisting of a spring in series with a dashpot known as the Maxwell model is

( ) *d d ta b dt dt*

/

( ) (b) *<sup>t</sup>*

Another model introduced by Zener [4], known as the Standard Linear Solid model is a combination of the above two models and is characterized by the constitutive relation

> 1 () 1 () *d d at t dt dt s e*

> > /

*t*

*e*

*t*


/ () 1 (a)

*t*

*s*

*t*

( ) (b)

Here *gJ* and *Ge* are known as the glass compliance and the equilibrium modulus respectively, and *c*+ and *c*- are the limiting values of creep for *t* ¥ and of relaxation

More recently methods of fractional calculus have been used [5-8] to generalize these models leading to the so-called Scott-Blair model characterized by the operator equation

1 ( ) 1 ( ),0 1 *d d at t*

*a a s ea* <sup>é</sup> ùé ù ê úê ú + =+ <£

*dt dt a a*

ë ûë û

<sup>é</sup> ù éù ê ú êú + =+ ë û ëû

*g*

*Gt G e*

= +

*e*

*Jt J e*

*c*

*c*

+ - -

<sup>é</sup> <sup>ù</sup> =+ -ê ú ë û

( ) (a)

*a t J t b b*

= +

*<sup>b</sup> Gt e a <sup>s</sup>* - *t*

=

*d* é ù - = -ê ú ë û

and by the material functions

and by the material functions

and by the material functions

for *t* = 0 respectively.

*dt e*

The objective of the chapter is to develop a lattice dynamical model based on the methods of fractional calculus so as to provide a microscopic basis for the theory of viscoelasticity. The plan of the chapter is as follows. First a brief review of the mechanical models of viscoelasticity is given. This is followed by an account of conventional lattice dynamical theory of viscoelaticity based on a model of linear chain of coupled oscillators with dissipative elements is given [14]. In the next section, lattice dynamical methods are extended to the model of linear chain of coupled fractional oscillators (requiring no additional dissipative elements) developed by the authors [15]. The response to sinusoidal forcing of the linear chain of fractional oscillators starting from a quiescent state is studied in the continuum limit and expressions for phase velocity, absorption and dispersion and specific dissipation function are derived. The results are discussed and numerical applications are presented in the final section

### **2. Mechanical models of viscoelasticity**

According to the linear theory of viscoelasticity, the body may be considered to be a system responding linearly to an excitation. A fundamental aspect of the theory both from a mathematical and a physical point of view is the response of the system to an excitation usually applied as a Heaviside step function. If a unit step of stress σ(t), is applied the material responds by undergoing a strain ε(t). The test is called creep test and the material function characterizing the response is called the creep compliance and is denoted by J(t). If on the other hand, if a unit strain is applied and the system responds with a stress σ(t), the test is called a relaxation test and the material function characterizing the system response, G(t), is called the *relaxation modulus.* The creep compliance and the relaxation modulus are defined through linear hereditary integrals of the Stieltjes type:

$$\begin{aligned} \varepsilon(t) &= \int\_0^t f(t-\tau)d\sigma(\tau) \\ \text{o}^- \\ \sigma(t) &= \int\_0^t \mathbf{G}(t-\tau)d\varepsilon(\tau) \end{aligned} \tag{1}$$

Simple mechanical models consisting of springs and dashpots incorporating elastic and dissipative functions have been advanced to explain the behaviour of viscoelastic solids and a good summary of the models of viscoelasticity can be found in the works of Caputo and Mainardi [5,6], Gorenflo and Mainardi[8], and Mainardi[9]. A simple model constituted by a spring in parallel with a dashpot, known as the Kelvin-Voigt model is characterized by the constitutive relation

Microscopic Formulation of Fractional Theory of Viscoelasticity 61

$$
\sigma(t) = m\varepsilon(t) + b\frac{d\varepsilon}{dt} \tag{2}
$$

and by the material functions

60 Viscoelasticity – From Theory to Biological Applications

applications are presented in the final section

**2. Mechanical models of viscoelasticity** 

defined through linear hereditary integrals of the Stieltjes type:

and plasticity.

constitutive relation

aspects of lattice dynamical theory of condensed matter can be found in renowned treatises [11-13]. Of particular interest to the present paper is the work of Askar [14] on the lattice dynamical foundations of continuum theories of elasticity, piezoelectricity, viscoelasticity

The objective of the chapter is to develop a lattice dynamical model based on the methods of fractional calculus so as to provide a microscopic basis for the theory of viscoelasticity. The plan of the chapter is as follows. First a brief review of the mechanical models of viscoelasticity is given. This is followed by an account of conventional lattice dynamical theory of viscoelaticity based on a model of linear chain of coupled oscillators with dissipative elements is given [14]. In the next section, lattice dynamical methods are extended to the model of linear chain of coupled fractional oscillators (requiring no additional dissipative elements) developed by the authors [15]. The response to sinusoidal forcing of the linear chain of fractional oscillators starting from a quiescent state is studied in the continuum limit and expressions for phase velocity, absorption and dispersion and specific dissipation function are derived. The results are discussed and numerical

According to the linear theory of viscoelasticity, the body may be considered to be a system responding linearly to an excitation. A fundamental aspect of the theory both from a mathematical and a physical point of view is the response of the system to an excitation usually applied as a Heaviside step function. If a unit step of stress σ(t), is applied the material responds by undergoing a strain ε(t). The test is called creep test and the material function characterizing the response is called the creep compliance and is denoted by J(t). If on the other hand, if a unit strain is applied and the system responds with a stress σ(t), the test is called a relaxation test and the material function characterizing the system response, G(t), is called the *relaxation modulus.* The creep compliance and the relaxation modulus are

0

*t*

ò


*t*

ò

() ( ) ( )

*t Jt d*

= -

*e t st*

() ( ) ( )

Simple mechanical models consisting of springs and dashpots incorporating elastic and dissipative functions have been advanced to explain the behaviour of viscoelastic solids and a good summary of the models of viscoelasticity can be found in the works of Caputo and Mainardi [5,6], Gorenflo and Mainardi[8], and Mainardi[9]. A simple model constituted by a spring in parallel with a dashpot, known as the Kelvin-Voigt model is characterized by the

(1)

*t Gt d*

= -

*s t et*

0


$$f(t) = \frac{1}{m} \left[ 1 - e^{-t/\tau\_c} \right] \tag{3}$$

$$G(t) = m + b\delta(t) \tag{b}$$

where *b m*/ *<sup>e</sup> t* = is referred to as the retardation time.

A model consisting of a spring in series with a dashpot known as the Maxwell model is characterized by the constitutive relation

$$
\sigma(t) + a \frac{d\sigma}{dt} = b \frac{d\varepsilon}{dt} \tag{4}
$$

and by the material functions

$$\begin{aligned} f(t) &= \frac{a}{b} + \frac{t}{b} \\ G(t) &= \frac{b}{a} e^{-t/\tau\_{\sigma}} \end{aligned} \tag{5}$$

where *a <sup>s</sup> t* = is referred to as the relaxation time.

Another model introduced by Zener [4], known as the Standard Linear Solid model is a combination of the above two models and is characterized by the constitutive relation

$$
\left[1 + a\frac{d}{dt}\right] \sigma(t) = \left[1 + \frac{d}{dt}\right] \varepsilon(t) \tag{6}
$$

and by the material functions

$$\mathbf{J}(t) = \mathbf{J}\_{\mathcal{g}} + \chi\_{+} \left[ \mathbf{1} - e^{-t/\tau\_{\varepsilon}} \right] \qquad \text{(a)}$$

$$\mathbf{G}(t) = \mathbf{G}\_{\varepsilon} + \chi\_{-} e^{-t/\tau\_{\sigma}} \qquad \text{(b)}$$

Here *gJ* and *Ge* are known as the glass compliance and the equilibrium modulus respectively, and *c*+ and *c*- are the limiting values of creep for *t* ¥ and of relaxation for *t* = 0 respectively.

More recently methods of fractional calculus have been used [5-8] to generalize these models leading to the so-called Scott-Blair model characterized by the operator equation

$$\left[1+a\frac{d^{\alpha}}{dt^{\alpha}}\right]\sigma(t) = \left[1+\frac{d^{\alpha}}{dt^{\alpha}}\right]\varepsilon(t), 0<\alpha \le 1\tag{8}$$

The details of the characteristic material functions and their behavior with respect to time are discussed in detail in reference [8].

## **3. Lattice dynamical theory of viscoelasticity**

#### **3.1. The elastic continuum**

Consider a one dimensional string of mass density ρ per unit length and E the Young's modulus. If the string is homogeneously stressed and released, elastic forces set up vibrations in the string corresponding to elastic waves propagating in the string. The equation governing this propagation is given by 9a and these waves propagate with a velocity given by 9b

$$\begin{aligned} \rho \frac{\partial^2 u}{\partial t^2} &= E \frac{\partial^2 u}{\partial x^2} \\ c^2 &= \frac{E}{\rho} \end{aligned} \quad \text{(a)} \tag{9}$$

Microscopic Formulation of Fractional Theory of Viscoelasticity 63

2 2 *u x* ¶ ¶

This process of taking the continuum limit is illustrated by considering the term on the right

 <sup>2</sup> 0 1 (2 ( ) ( ) ( ) *nn n*

*ut u t ut*

2 2 1 1

*aa a <sup>w</sup>* <sup>+</sup> - <sup>ç</sup> æ ö æ ö - - ÷÷ <sup>ç</sup> <sup>ç</sup> ÷÷ <sup>ç</sup> <sup>ç</sup> - ÷÷ <sup>ç</sup> ÷÷ çç çè ø÷÷ è ø

In the limit of *N* ¥ and *a* 0 , the chain of atoms can be treated as a continuum and

() ( ,) *<sup>n</sup> u t uxt* , a continuous function. The term in the brackets can be written as

<sup>0</sup> and can be expressed as . In the limit and / / *k ka a a ka E m a m ma*

1 ( () () ( () () *n n n n u t ut ut u t*

*w r* = = The term

reduces to the square of the wave velocity. Thus the equation of motion in the limit

In order to develop a theory of viscoelasticity, the model of linear chain of atoms has been generalized to include dissipative effects [14] by incorporating dashpots either in parallel to the springs as shown in figure 2 to yield the lattice dynamical version of the Kelvin-Voigt model, or the dashpots in series with the springs as shown in figure 3 to yield the lattice

It has been shown [14] that the models in the continuum limit lead to the stress-strain

¶ ¶ <sup>=</sup> ¶ ¶ yielding the equation of motion of a continuum, equation

0

*a*

2 2 2 2 2 *u u c t x*

**3.2. Viscoelasticity: Lattice dynamical approach** 

2 2 2

dynamical version of the Maxwell model.

**Figure 2.** Kelvin-Voigt model for viscoelasticity

**Figure 3.** Maxwell model for viscoelasticity

hand side of equation (10):

the factor

2 2 <sup>0</sup> *a* 

(9a)

relations

can be written as

can be written as

The conventional model of lattice dynamics consists of a chain of N identical masses m connected by springs of spring constant k, shown in Figure 1. The equation of motion for the nth mass is given in the standard notation by

$$\begin{aligned} m \frac{d^2 u\_n}{dt^2} &= -k(2u\_n - u\_{n+1} - u\_{n-1}) \\ \text{or } \frac{d^2 u\_n}{dt^2} &= -\omega\_0^2 (2u\_n - u\_{n+1} - u\_{n-1}) \end{aligned} \tag{10}$$

where *<sup>n</sup> u* is the displacement from the equilibrium position of the nth mass. The right hand side represents the elastic restoring force on the nth atom. The atoms n=1 and n= N define the boundaries.

The frequency of oscillation is given by <sup>2</sup> <sup>0</sup> *w* = *k m*/ .

#### **Figure 1.** Linear chain of coupled oscillators

The elastic behavior of a continuous chain in equation (9a) can be obtained from the equation of motion of the linear chain of masses connected to each other by springs considered above in equation (9). In the limit when the number of masses N and the separation between the masses a0, such that the product Na L, a finite length, the linear chain of coupled oscillators reduces to a continuous loaded string and is referred to as the long wavelength limit.

This process of taking the continuum limit is illustrated by considering the term on the right hand side of equation (10):

the factor

62 Viscoelasticity – From Theory to Biological Applications

are discussed in detail in reference [8].

nth mass is given in the standard notation by

The frequency of oscillation is given by <sup>2</sup>

**Figure 1.** Linear chain of coupled oscillators

**3.1. The elastic continuum** 

velocity given by 9b

the boundaries.

long wavelength limit.

**3. Lattice dynamical theory of viscoelasticity** 

The details of the characteristic material functions and their behavior with respect to time

Consider a one dimensional string of mass density ρ per unit length and E the Young's modulus. If the string is homogeneously stressed and released, elastic forces set up vibrations in the string corresponding to elastic waves propagating in the string. The equation governing this propagation is given by 9a and these waves propagate with a

> 2 2 2 2

¶ ¶ <sup>=</sup> ¶ ¶

*r*

*u u <sup>E</sup> t x E*

The conventional model of lattice dynamics consists of a chain of N identical masses m connected by springs of spring constant k, shown in Figure 1. The equation of motion for the

2 1 1

=- - -

or (2 )

where *<sup>n</sup> u* is the displacement from the equilibrium position of the nth mass. The right hand side represents the elastic restoring force on the nth atom. The atoms n=1 and n= N define

<sup>0</sup> *w* = *k m*/ .

The elastic behavior of a continuous chain in equation (9a) can be obtained from the equation of motion of the linear chain of masses connected to each other by springs considered above in equation (9). In the limit when the number of masses N and the separation between the masses a0, such that the product Na L, a finite length, the linear chain of coupled oscillators reduces to a continuous loaded string and is referred to as the

2 0 11

=- - -

(2 )

*nn n*

*nn n*

+ -

+ -

*uu u*

2

*w*

*d u m ku u u*

(a)

(9)

(10)

(b)

2

=

*c*

2

*dt d u*

2

*dt*

*n*

*n*

*r*

$$
\rho\_0^2 \left( (2\mu\_n(t) - \mu\_{n+1}(t) - \mu\_n(t)) \right)
$$

can be written as

$$a^2 \omega\_0^2 \left| \frac{1}{a} \left( \frac{(\mu\_{n+1}(t) - \mu\_n(t))}{a} - \frac{(\mu\_n(t) - \mu\_{n-1}(t))}{a} \right) \right|$$

In the limit of *N* ¥ and *a* 0 , the chain of atoms can be treated as a continuum and () ( ,) *<sup>n</sup> u t uxt* , a continuous function. The term in the brackets can be written as 2 2 *u x* ¶ ¶ 2 2 2 <sup>0</sup> and can be expressed as . In the limit and / / *k ka a a ka E m a m ma w r* = = The term 2 2 <sup>0</sup> *a* reduces to the square of the wave velocity. Thus the equation of motion in the limit can be written as 2 2 2 2 2 *u u c t x* ¶ ¶ <sup>=</sup> ¶ ¶ yielding the equation of motion of a continuum, equation (9a)

#### **3.2. Viscoelasticity: Lattice dynamical approach**

In order to develop a theory of viscoelasticity, the model of linear chain of atoms has been generalized to include dissipative effects [14] by incorporating dashpots either in parallel to the springs as shown in figure 2 to yield the lattice dynamical version of the Kelvin-Voigt model, or the dashpots in series with the springs as shown in figure 3 to yield the lattice dynamical version of the Maxwell model.

**Figure 2.** Kelvin-Voigt model for viscoelasticity

**Figure 3.** Maxwell model for viscoelasticity

It has been shown [14] that the models in the continuum limit lead to the stress-strain relations

$$\begin{aligned} \sigma &= E\varepsilon + \zeta \frac{d\varepsilon}{dt} \\ \frac{d\varepsilon}{dt} &= \frac{\sigma}{E} + \frac{1}{E} \frac{d\sigma}{dt} \end{aligned} \tag{11}$$

Microscopic Formulation of Fractional Theory of Viscoelasticity 65

fractional oscillator behaves like a damped harmonic oscillator [15-18], the oscillations of the linear chain of fractional oscillators can be sustained only if it is driven by an external force. Hence a periodic force, taken to be sinusoidal without loss of generality, acting only on the first member of the chain is included. Thus the equations of motion are given

( ) <sup>0</sup> <sup>1</sup> (0) (0) 2 ( ) ( ) ( ) ( ) , 2 1 ( ) 1 1

<sup>1</sup> (0) (0) 2 ( ) ( ) ( ) ( )( ) ( ) ( )

*uu ut u u t d f t d*

0 (0) (0) 2 ( ) ( ) ( ) ( ) *t NN N n N uu ut u u t d*

The last term on the right hand side of equation (16) represents the effect of the periodic

The set of equations (15) through (17) can be solved formally by using the Laplace transform

<sup>0</sup> { ( )} ( ) () . *st*

The Laplace transform of a fractional integral of order *a a* , (Re 0) > defined by [20] and can be evaluated by considering the fractional integral as a convolution of two functions (19b),

<sup>1</sup> ( ) ( )( ) (a) ( )

*tt t <sup>a</sup>*

( ) and ( ) (b) ( )

<sup>1</sup> ( ) ( )( ) ( ) \* ( ), ( )


<sup>+</sup> = -= <sup>G</sup> ò , (20)

<sup>0</sup> *LI f t L t f t s f s* { ( )} { ( ) \* ( )} ( ) ( ); *<sup>a</sup> f f* <sup>+</sup> = = (21)

*I ft f t d t ft a a t t tf <sup>a</sup>*

0

= - <sup>G</sup>

ò

*t I ft f t d*

*a a*

1


*a*

*<sup>f</sup> <sup>a</sup>*

= G

*<sup>t</sup> <sup>t</sup> f t*

0

and using the formula for the Laplace transform of the convolution. Thus

*t*

technique [20-24]. In the standard notation, the Laplace transform is defined by


0 0

*t t*

( ) <sup>0</sup> 1 1

G G ò ò (16)

*<sup>a</sup> <sup>w</sup> a a t t t t tt t a a*

( ) <sup>0</sup> <sup>1</sup> 1

*<sup>a</sup> <sup>w</sup> <sup>a</sup> t t tt <sup>a</sup>* - =+ - - - - <sup>G</sup> ò (17)

*nn n L u t u s e u t dt* ¥ - = = ò (18)

(19)

1


*aw <sup>a</sup> t t tt t <sup>a</sup>* - =+ - <sup>ò</sup> - - - ££ - <sup>G</sup> + - (15)

*u u u t u u u t d nN*

0

*nn n nn n*

11 1 1 2

0

 <sup>1</sup> 0

+

forcing as already noted.

**4.2. Formal solution** 

i.e.,

*t*

by

respectively for the two models. Harmonic waves are found to propagate with decaying amplitude indicating dissipation. The details can be found in reference [14].

#### **4. Fractional calculus approach**

In this paper, the lattice dynamical approach is generalized by using the methods of fractional calculus, leading to the so-called linear chain of coupled fractional oscillators. However, this generalization does not require the dissipative elements, namely the dashpots to be explicitly considered, for the dissipation is intrinsic to the fractional oscillators [16-19]. As the method of approach is quite different, further reference will be made not to the work of Askar [14], but to the work of the authors on a linear chain of coupled fractional oscillators [19] and the approach is outlined below.

#### **4.1. Linear chain of fractional oscillators**

Instead of starting with equation (10), we consider the integral equation of motion for the nth mass interacting with only its nearest neighbors, which can be written without loss of generality as:

$$u\_n(t) = u\_n(0) + \dot{u}\_n(0)t - \omega\_0^2 \int\_0^t (2u\_n(\tau) - u\_{n+1}(\tau) - u\_{n-1}(\tau))(t-\tau)d\tau \ 2 \le n \le N-1 \tag{12}$$

Here (0) *<sup>n</sup> u* and (0) *<sup>n</sup> u* refer to the displacement from equilibrium and velocity of the nth mass at t=0. For the masses at the ends of the chain, the equations of motion are given by

$$u\_1(t) = u\_1(0) + \dot{u}\_1(0)t - \omega\_0^2 \int\_0^t (2u\_1(\tau) - u\_2)(t-\tau)d\tau + \int\_0^t f(\tau)(t-\tau)d\tau\tag{13}$$

and

$$
\mu\_N(t) = \mu\_N(0) + \dot{\nu}\_N(0)t - \omega\_0^2 \int\_0^t (2u\_N(\tau) - u\_{N-1}(\tau))(t-\tau)d\tau,\tag{14}
$$

respectively. Of course, it is <sup>2</sup> <sup>0</sup> *w* = *k m*/ , in the usual notation. The last term in equation (13) indicates a periodic forcing on the end atom.

The integrals on the right hand side of equations (11) -(13) are generalized to fractional integrals of order *a* (defined by equation (19) below) to yield the equations of motion of a chain of coupled fractional oscillators. However, as has been well known that a fractional oscillator behaves like a damped harmonic oscillator [15-18], the oscillations of the linear chain of fractional oscillators can be sustained only if it is driven by an external force. Hence a periodic force, taken to be sinusoidal without loss of generality, acting only on the first member of the chain is included. Thus the equations of motion are given by

$$u\_n = u\_n(0) + \dot{u}\_n(0)t - \frac{\omega\_0^\alpha}{\Gamma(\alpha)} \int\_0^\alpha (2u\_n(\tau) - u\_{n+1}(\tau) - u\_{n-1}(\tau))(t-\tau)^{\alpha-1} d\tau, 2 \le n \le N-1 \tag{15}$$

$$u\_1 = u\_1(0) + \dot{u}\_1(0)t - \frac{\omega\_0^\alpha}{\Gamma(\alpha)} \int\_0^t (2u\_1(\tau) - u\_2(\tau))(t-\tau)^{\alpha-1} d\tau + \frac{1}{\Gamma(\alpha)} \int\_0^t f(\tau)(t-\tau)^{\alpha-1} d\tau \tag{16}$$

$$u\_N = u\_N(0) + \dot{u}\_N(0)t - \frac{\omega\_0^\alpha}{\Gamma(\alpha)} \int\_0^t (2u\_n(\tau) - u\_{N-1}(\tau))(t-\tau)^{\alpha-1} d\tau \tag{17}$$

The last term on the right hand side of equation (16) represents the effect of the periodic forcing as already noted.

#### **4.2. Formal solution**

64 Viscoelasticity – From Theory to Biological Applications

**4. Fractional calculus approach** 

oscillators [19] and the approach is outlined below.

**4.1. Linear chain of fractional oscillators** 

generality as:

and

respectively. Of course, it is <sup>2</sup>

indicates a periodic forcing on the end atom.

(a)

(11)

<sup>1</sup> (b)

*<sup>d</sup> <sup>E</sup> dt*

*es s*

respectively for the two models. Harmonic waves are found to propagate with decaying

In this paper, the lattice dynamical approach is generalized by using the methods of fractional calculus, leading to the so-called linear chain of coupled fractional oscillators. However, this generalization does not require the dissipative elements, namely the dashpots to be explicitly considered, for the dissipation is intrinsic to the fractional oscillators [16-19]. As the method of approach is quite different, further reference will be made not to the work of Askar [14], but to the work of the authors on a linear chain of coupled fractional

Instead of starting with equation (10), we consider the integral equation of motion for the nth mass interacting with only its nearest neighbors, which can be written without loss of

> ( ) <sup>2</sup> 0 11

*nnn nn n ut u u t u u u t d n N w t t t tt* = + - - - - ££ - ò + - (12)

0 0

( ) <sup>2</sup> 0 1

*NNN NN ut u u t u u t d w t t tt* =+ - - - ò - (14)

<sup>0</sup> *w* = *k m*/ , in the usual notation. The last term in equation (13)

*t t ut u u t u u t d f t d* = + - --+ - *w t tt t tt* ò ò (13)

( ) (0) (0) 2 ( ) ( ) ( ) ( ) 2 1

Here (0) *<sup>n</sup> u* and (0) *<sup>n</sup> u* refer to the displacement from equilibrium and velocity of the nth mass at t=0. For the masses at the ends of the chain, the equations of motion are given by

( ) <sup>2</sup>

0 ( ) (0) (0) 2 ( ) ( ) ( ) , *t*

The integrals on the right hand side of equations (11) -(13) are generalized to fractional integrals of order *a* (defined by equation (19) below) to yield the equations of motion of a chain of coupled fractional oscillators. However, as has been well known that a

( ) (0) (0) 2 ( ) ( ) ( )( )

0

1 1 1 01 2

*t*

*d d dt E E dt*

= +

amplitude indicating dissipation. The details can be found in reference [14].

*<sup>e</sup> s ez*

= +

The set of equations (15) through (17) can be solved formally by using the Laplace transform technique [20-24]. In the standard notation, the Laplace transform is defined by

$$L\{\boldsymbol{\mu}\_n(t)\} = \tilde{\boldsymbol{\mu}}\_n(s) = \int\_0^\infty e^{-st} \boldsymbol{\mu}\_n(t) dt. \tag{18}$$

The Laplace transform of a fractional integral of order *a a* , (Re 0) > defined by [20] and can be evaluated by considering the fractional integral as a convolution of two functions (19b),

$$I\_{0+}^{\alpha}f(t) = \frac{1}{\Gamma\{\alpha\}} \int\_{0}^{t} f(\tau)(t-\tau)^{\alpha-1}d\tau \qquad \text{(a)}$$

$$\phi(t) = \frac{t^{\alpha-1}}{\Gamma\{\alpha\}} \text{ and } f(t) \tag{b}$$

i.e.,

$$I\_{0+}^{\\\alpha}f(t) = \frac{1}{\Gamma(\alpha)} \int\_0^t f(\tau)(t-\tau)^{\alpha-1}d\tau = \phi(t)^\* \, f(t) \, . \tag{20}$$

and using the formula for the Laplace transform of the convolution. Thus

$$L\{I\_{0+}^{\alpha}f(t)\} = L\{\phi(t)^{\*}\}\,\,f(t)\} = \tilde{\phi}(s)\,\tilde{f}(s);\tag{21}$$

But as discussed in [23],

$$\bar{\phi}(\mathbf{s}) = L\{\frac{t^{\alpha - 1}}{\Gamma(\alpha)}\} = \mathbf{s}^{-\alpha} \text{ } \text{Re}\,\alpha > 0\tag{22}$$

Microscopic Formulation of Fractional Theory of Viscoelasticity 67

Then the set of linear equations can be explicitly written down and the inverse Laplace transform applied to these equations on both sides term by term to obtain the displacements as functions of time. This formal procedure appears to be straight enough, however, it is not possible to carry out further simplifications in closed form of the expressions in equation. (28). Nevertheless, this set of equations can be numerically solved for specific values of N, the number of oscillators in the chain. Such calculations have been carried out and the details of numerical applications may be found in [19]. In the next section, the continuum

It would appear to be a straight forward procedure to extend the same considerations to the equations of motion of a chain of coupled fractional oscillators given in equations (15)-(17).

However, it is not so straight forward and due caution has to be exercised in extending such a procedure. The reason is the occurrence of time derivatives of fractional order. Douglas has shown [25] that in the context of fractals, inhomogeneity in space results in the appearance of fractional order time derivatives and inhomogeneity in time implies fractional order derivatives of space. Since fractional order time dependence is being investigated in the present work, the question arises whether inhomogenneity in space can be ignored and whether space derivatives can be taken. The question is one of scale. The long wavelength limit implies that the lengths involved are very much larger than the separation distance between masses. On this scale the inhomogeneity in space can be ignored and space derivative can be interpreted as an average of some sort and a hand waving justification can be made. The limiting form of the equations (12) through (14) can be

0 11 ((2 ( ) ( ) ( ) *uu u nn n* ) *<sup>a</sup> wt t t* - -- + -

2 1 1

*aa a <sup>a</sup> t t t t <sup>w</sup>* <sup>+</sup> - <sup>ç</sup> æ ö æ ö - - ÷÷ <sup>ç</sup> <sup>ç</sup> ÷÷ <sup>ç</sup> <sup>ç</sup> - ÷÷ <sup>ç</sup> ÷÷ <sup>ç</sup> ç è ø ç ÷÷ è ø

2 2

*a a*

2 2 *u x*(,) *x* ¶ *t* ¶

1 ( () () ( () () *n n n n u u u u*

<sup>0</sup> / *k ka*

In the limit a 0, ( ) *<sup>n</sup> u t* becomes a continuous variable *u x*(,) *t* and the expression in

*m ma*

*<sup>a</sup> w* = = (31)

and *ka m a k r* , / represent the tension and

(30)

limit of these equations is studied.

Thus the following factor in equation (15)

0

*a*

**5. The continuum limit** 

considered.

can be written as

parenthesis in (30) reduces to

and

Taking Laplace transforms on both sides of the set of equations (15) through (17) yields

$$
\tilde{u}\_n(\mathbf{s}) = \boldsymbol{u}\_n(\mathbf{0})\mathbf{s}^{-1} + \dot{\boldsymbol{u}}\_n(\mathbf{0})\mathbf{s}^{-2} - \omega\_0^\alpha \mathbf{s}^{-\alpha} \left(2\bar{\boldsymbol{u}}\_n(\mathbf{s}) - \bar{\boldsymbol{u}}\_{n+1}(\mathbf{s}) - \bar{\boldsymbol{u}}\_{n-1}(\mathbf{s})\right), \\
2 \le n \le N - 1 \tag{23}
$$

$$
\tilde{\boldsymbol{\mu}}\_{1}(\mathbf{s}) = \boldsymbol{\mu}\_{1}(\mathbf{0})\mathbf{s}^{-1} + \dot{\boldsymbol{\mu}}\_{1}(\mathbf{0})\mathbf{s}^{-2} - \omega\_{0}^{\alpha}\mathbf{s}^{-\alpha} \left(2\tilde{\boldsymbol{\mu}}\_{1}(\mathbf{s}) - \tilde{\boldsymbol{\mu}}\_{2}(\mathbf{s})\right) + \tilde{f}(\mathbf{s})\mathbf{s}^{-\alpha} \tag{24}
$$

$$
\tilde{\boldsymbol{\mu}}\_{N}(\mathbf{s}) = \boldsymbol{\mu}\_{N}(\mathbf{0})\mathbf{s}^{-1} + \dot{\boldsymbol{\mu}}\_{N}(\mathbf{0})\mathbf{s}^{-2} - \boldsymbol{\omega}\_{0}^{\alpha}\mathbf{s}^{-\alpha} \left(2\tilde{\boldsymbol{\mu}}\_{N}(\mathbf{s}) - \tilde{\boldsymbol{\mu}}\_{N-1}(\mathbf{s})\right) \tag{25}
$$

The set of equations can be rewritten in the following matrix form:

$$
\tilde{A}(\mathbf{s})\tilde{\mathcal{U}}(\mathbf{s}) = \mathcal{U}(\mathbf{0})\mathbf{s}^{-1} + \dot{\mathcal{U}}(\mathbf{0})\mathbf{s}^{-2} + \tilde{F}(\mathbf{s})\tag{26}
$$

where the N x N matrix *A*( )*s* is given by

$$
\tilde{A}(\mathbf{s}) = \begin{pmatrix}
1 + 2\omega\_0^\alpha s^{-\alpha} & -\omega\_0^\alpha s^{-\alpha} & 0 & 0 & \cdots & 0 \\
\cdot & -\omega\_0^\alpha s^{-\alpha} & 1 + 2\omega\_0^\alpha s^{-\alpha} & -\omega\_0^\alpha s^{-\alpha} & 0 & 0 & \cdots \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & \cdot & 0 \\
\cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\
0 & 0 & \cdot & \cdot & 0 & -\omega\_0^\alpha s^{-\alpha} & 1 + 2\omega\_0^\alpha s^{-\alpha} \\
0 & 0 & \cdot & \cdot & 0 & -\omega\_0^\alpha s^{-\alpha} & 1 + 2\omega\_0^\alpha s^{-\alpha}
\end{pmatrix}
(27)
$$

and the N x 1 column vectors are given by

$$
\bar{\boldsymbol{U}}(\boldsymbol{s}) = \begin{pmatrix} \bar{\boldsymbol{u}}\_{1}(\boldsymbol{s}) \\ \bar{\boldsymbol{u}}\_{2}(\boldsymbol{s}) \\ \bar{\boldsymbol{u}}\_{3}(\boldsymbol{s}) \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \bar{\boldsymbol{u}}\_{N-1}(\boldsymbol{s}) \\ \end{pmatrix} \circ \boldsymbol{U}(\boldsymbol{0}) = \begin{pmatrix} \boldsymbol{u}\_{1}(\boldsymbol{0}) \\ \boldsymbol{u}\_{2}(\boldsymbol{0}) \\ \boldsymbol{u}\_{3}(\boldsymbol{0}) \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \end{pmatrix} \circ \boldsymbol{U}(\boldsymbol{0}) = \begin{pmatrix} \boldsymbol{\bar{u}}\_{1}(\boldsymbol{0}) \\ \boldsymbol{\bar{u}}\_{2}(\boldsymbol{0}) \\ \boldsymbol{\bar{u}}\_{3}(\boldsymbol{0}) \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \end{pmatrix} \quad \text{and} \ \bar{\boldsymbol{F}}(\boldsymbol{s}) = \begin{pmatrix} \bar{f}(\boldsymbol{s})\boldsymbol{s}^{-\alpha} \\ \boldsymbol{0} \\ 0 \\ 0 \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \end{pmatrix} \tag{28}
$$

respectively. It may be noted that among the four column vectors, only the first and the last refer to Laplace transformed quantities, but the middle two refer to constants describing the initial conditions. The solution to Eq. (25) can be formally written as

$$\bar{M}(\mathbf{s}) = \bar{A}^{-1}(\mathbf{s}) \left( \mathcal{U}(\mathbf{0}) \mathbf{s}^{-1} + \dot{\mathcal{U}}(\mathbf{0}) \mathbf{s}^{-2} + \bar{F}(\mathbf{s}) \right) \tag{29}$$

Then the set of linear equations can be explicitly written down and the inverse Laplace transform applied to these equations on both sides term by term to obtain the displacements as functions of time. This formal procedure appears to be straight enough, however, it is not possible to carry out further simplifications in closed form of the expressions in equation. (28). Nevertheless, this set of equations can be numerically solved for specific values of N, the number of oscillators in the chain. Such calculations have been carried out and the details of numerical applications may be found in [19]. In the next section, the continuum limit of these equations is studied.

## **5. The continuum limit**

66 Viscoelasticity – From Theory to Biological Applications

where the N x N matrix *A*( )*s* is given by

0 0

*w w*

*A s*

*s s*

*aa aa*

and the N x 1 column vectors are given by

 


0 00

*w ww*

*s ss*

*aa aa aa*


3 3

æö æ ö ç ç ÷ ÷

*U s U U*

*us u u s u* - -

( ) (0) ( ) . . (0) (0) . . ( ) (0) ( ) (0) *N N N N*

= ==

initial conditions. The solution to Eq. (25) can be formally written as

1 1

ç÷ ç ÷ ç ç ÷ ÷ èø è ø

1 ( ) { } , Re 0 ( )

Taking Laplace transforms on both sides of the set of equations (15) through (17) yields

( ) 1 2

1 2 0 0. . 0

. . . ( ) . . ..

1 2 . . ..

<sup>0</sup> 0 . . 1 2

<sup>÷</sup> <sup>ç</sup> -+- <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> - +- <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

<sup>æ</sup> <sup>ö</sup> <sup>ç</sup> + - <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

. . 0 . . ..

<sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> - + - <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> - + ÷÷ <sup>è</sup> <sup>ø</sup>

<sup>÷</sup> <sup>ç</sup> <sup>÷</sup> =ç <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

(27)

0 0 .. 0 1 2

respectively. It may be noted that among the four column vectors, only the first and the last refer to Laplace transformed quantities, but the middle two refer to constants describing the

0 . . 1 2 0

3

1

( ) 1 12 *Us A s U s U s Fs* ( ) ( ) (0) (0) ( ) - -- = ++ (29)


*N N*

*u u*

 

( ) 1 2

<sup>0</sup> 1 1 ( ) (0) (0) 2 ( ) ( ) ( ) ,2 1 *nn n nn n us u s u s s us u s u s n N a a <sup>w</sup>* - -- =+- - - ££ - + - (23)

1 1 1 0 12 *u s u s u s s u s u s f ss* ( ) (0) (0) 2 () () () *a a <sup>a</sup> <sup>w</sup>* - -- - = + - -+ (24)

( ) 1 2 0 1 ( ) (0) (0) 2 () () *NN N N N us u s u s s us u s a a <sup>w</sup>* - -- =+- - - (25)

(22)

1 2 *<sup>A</sup>*( ) ( ) (0) (0) ( ) *sU s U s U s F s* - - =++ (26)

0 0 0

*w w w*

(0) 0 . and ( ) . . . (0) 0 (0) 0

æ ö æ ö <sup>ç</sup> ÷ ç <sup>÷</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> =ç <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> ç ÷ ç ÷ <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> ÷÷ è ø è ø

*F s*

*s s s*

*a a aa aa*


0 0

*w w*

*s s*

*a a a a*

( ) 0 -*a*

(28)

*f s s*


*<sup>t</sup> sL s <sup>a</sup> <sup>a</sup> f a <sup>a</sup>* - - == >

G

The set of equations can be rewritten in the following matrix form:

0 00

*w ww*

*s ss*

1 11 2 22

*u s u u u s u u u s u u*

 

( ) (0) (0) ( ) (0) (0)

*aa aa aa*


But as discussed in [23],

It would appear to be a straight forward procedure to extend the same considerations to the equations of motion of a chain of coupled fractional oscillators given in equations (15)-(17).

However, it is not so straight forward and due caution has to be exercised in extending such a procedure. The reason is the occurrence of time derivatives of fractional order. Douglas has shown [25] that in the context of fractals, inhomogeneity in space results in the appearance of fractional order time derivatives and inhomogeneity in time implies fractional order derivatives of space. Since fractional order time dependence is being investigated in the present work, the question arises whether inhomogenneity in space can be ignored and whether space derivatives can be taken. The question is one of scale. The long wavelength limit implies that the lengths involved are very much larger than the separation distance between masses. On this scale the inhomogeneity in space can be ignored and space derivative can be interpreted as an average of some sort and a hand waving justification can be made. The limiting form of the equations (12) through (14) can be considered.

Thus the following factor in equation (15)

$$-\omega\_0^\alpha \left( (2\mu\_n(\tau) - \mu\_{n+1}(\tau) - \mu\_{n-1}(\tau)) \right)$$

can be written as

$$a^2 \omega\_0^\alpha \left| \frac{1}{a} \left( \frac{\mu\_{n+1}(\tau) - \mu\_n(\tau)}{a} - \frac{(\mu\_n(\tau) - \mu\_{n-1}(\tau))}{a} \right) \right| \tag{30}$$

and

$$a^2 \omega\_0^\alpha = a^2 \frac{k}{m} = \frac{ka}{m/a} \tag{31}$$

In the limit a 0, ( ) *<sup>n</sup> u t* becomes a continuous variable *u x*(,) *t* and the expression in parenthesis in (30) reduces to 2 2 *u x*(,) *x* ¶ *t* ¶ and *ka m a k r* , / represent the tension and

the mass density respectively. The expression in (30) can be written in terms of a quantity 0*c <sup>a</sup>* , where <sup>0</sup>*c* has the dimensions of 'velocity'. Assuming that at t=0 the displacement at the free end is subject to sinusoidal forcing, *u t ft A t* (0, ) ( ) sin( ) = = *w* , the equation of motion reduces in the continuum limit to

$$u(\mathbf{x},t) = \frac{c\_0^\alpha}{\Gamma(\alpha)} \int\_0^t \frac{\partial^2 u(\mathbf{x}, \tau)}{\partial^2 \mathbf{x}} (t-\tau)^{\alpha-1} d\tau \tag{32}$$

Microscopic Formulation of Fractional Theory of Viscoelasticity 69


0

0

*w pa w*

*c*

<sup>÷</sup> <sup>ç</sup> <sup>÷</sup> è ø

<sup>÷</sup> <sup>+</sup> <sup>ç</sup> <sup>÷</sup> è ø (39)

(40)

The transient part, ( ,) *tr u xt* , arises from the Hankel loop consisting of the small circle and

( ,) (, , , ) *rt tr u x t e K r c x dr <sup>a</sup> <sup>a</sup> w*

> /2 0

(, , ,) sin ( ) sin( / 2) ( )

*r c*

The steady part, ( ,) *st u xt* , arises from the residues of the poles of the integrand in eq. (36) at

/2 cos( /4)

( ,) sin ( sin( / 4))

*<sup>w</sup> <sup>a</sup> pa <sup>w</sup>*

This completes the formal analysis of the response to sinusoidal forcing of a chain of fractional oscillators in the continuum limit. Numerical results will be presented in the next

The results of the last section in equations (36)-(40) are new. Applications of fractional calculus to oscillation and dissipation problems and solutions to integral and fractional order differential equations can be found in Gorenflo and Mainardi [8] Impulse response functions are also discussed in their work. The same methodology has been extended to sinusoidal forcing of a linear chain in the present work. The response of a chain of fractional oscillators in the continuum limit, when subject to a sinusoidal forcing starting from a quiescent state consists of a transient part and a steady part. The transient part decays in time and approaches zero as *t* ¥. In the limit *a* 2 the transient part vanishes entirely. Furthermore, it exhibits attenuation as a function of distance from the end as indicated by the spatial dependence of the kernel in eq. (39). No simple closed form expressions can be obtained, the only recourse is through numerical integration. Closed form solutions are not available even for simple mechanical models such as the standard linear solid, which are not based on fractional calculus, for the reason that the solutions are deemed mathematically cumbersome as has been discussed by Caputo and Mainardi [5]. These authors have applied methods of fractional calculus to the linear theory of viscoelasticity based on mechanical models and have discussed the propagation of viscoelastic waves[6]. Their work is macroscopic in its approach and is not based on the continuum limit of lattice dynamics. However, some concepts developed by these authors are quite insightful and the results of

æ ö <sup>ç</sup> <sup>÷</sup> = - <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

*w pa*

0

( ) cos( /2) /2


two lines parallel to the negative x-axis (as shown in figure 1in [17]) and is given by

0

*<sup>r</sup> <sup>x</sup> A r <sup>c</sup> Kr cx <sup>e</sup> <sup>x</sup>*

*<sup>a</sup> pa a a*

0 2 2

*w*

/2

*<sup>x</sup> u x t Ae <sup>t</sup>*

*<sup>x</sup> <sup>c</sup>*

*a*

0

æ ö <sup>ç</sup> <sup>÷</sup> -ç <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>÷</sup> çç <sup>÷</sup> è ø

*p w*

*a*

*st*

**6. Results and discussion** 

*s i* = *w* and is given by

with

section.

¥

with the initial conditions

$$
\mu(\mathfrak{x},0) = 0 \text{ and } \dot{\mathfrak{u}}(\mathfrak{x},0) = 0 \text{ for } \mathfrak{x} > 0 \text{ and } \mathfrak{u}(0,t) = f(t).
$$

Taking the Laplace transform on both sides of equation. (32) yields

$$
\tilde{u}(\mathbf{x}, \mathbf{s}) = \frac{c\_0^\alpha}{s^\alpha} \frac{d^2 \tilde{u}(\mathbf{x}, \mathbf{s})}{d\mathbf{x}^2} \tag{33}
$$

The equation. (33) can be rewritten as

$$\frac{d^2\tilde{u}(\mathbf{x},s)}{d\mathbf{x}^2} - \frac{s^\alpha}{c\_0^\alpha}\tilde{u}(\mathbf{x},s) = 0 \text{ for } \mathbf{x} \neq \mathbf{0} \text{ and } \tilde{u}(0,s) = \tilde{f}(s) \tag{34}$$

This is an ordinary differential equation and can be solved to yield

$$
\bar{u}(x,s) = \bar{f}(s)e^{-\sqrt{\frac{s^\alpha}{c\_0^\alpha}}x} \tag{35}
$$

Substituting for 2 2 ( ) *<sup>A</sup> f s s w <sup>w</sup>* <sup>=</sup> <sup>+</sup> corresponding to a sinusoidal forcing *f*( ) sin( ) *tA t* = *w* ,

The eq. (35) can be inverted as a Bromwich integral

$$u(\mathbf{x},t) = \frac{A\omega}{2\pi i} \int\_{Br} \frac{e^{st}e^{-\int\_{t\_0}^{\alpha} \overline{\mathbf{x}}\_0}}{\left(s^2 + \omega^2\right)} ds\tag{36}$$

The Bromwich integral in eq.(36) can be evaluated by appealing to the theory of complex variables [22-24]. By considering a Hankel-Bromwich path, it can be evaluated as the sum of two contributions:

$$
\mu(\mathbf{x},t) = \mu\_{tr}(\mathbf{x},t) + \mu\_{st}(\mathbf{s},t) \tag{37}
$$

representing a transient part and a steady state part respectively.

The transient part, ( ,) *tr u xt* , arises from the Hankel loop consisting of the small circle and two lines parallel to the negative x-axis (as shown in figure 1in [17]) and is given by

$$\mu\_{tr}(\mathbf{x},t) = \int\_0^\infty e^{-rt} K\_\alpha(r,\omega,\mathbf{c}\_0^\alpha,\mathbf{x}) dr \tag{38}$$

with

68 Viscoelasticity – From Theory to Biological Applications

reduces in the continuum limit to

The equation. (33) can be rewritten as

Substituting for 2 2 ( ) *<sup>A</sup> f s*

two contributions:

2

*s*

The eq. (35) can be inverted as a Bromwich integral

*w <sup>w</sup>* <sup>=</sup> <sup>+</sup>

2

with the initial conditions

0*c*

the mass density respectively. The expression in (30) can be written in terms of a quantity

*<sup>a</sup>* , where <sup>0</sup>*c* has the dimensions of 'velocity'. Assuming that at t=0 the displacement at the free end is subject to sinusoidal forcing, *u t ft A t* (0, ) ( ) sin( ) = = *w* , the equation of motion

> 2 0 1 2

*ux ux x u t f t* ( ,0) 0 and ( ,0) 0 for 0 and (0, ) ( ). = => =

( ,) (,) *<sup>c</sup> duxs*

( ,) ( , ) 0 for 0 and (0, ) ( ) *duxs s uxs x u s f s dx c*

<sup>0</sup> ( ,) () *<sup>s</sup> <sup>x</sup> <sup>c</sup> uxs f se*

2 2 ( ,) <sup>2</sup> ( )

*w*

*Br <sup>A</sup> e e ds uxt i s*

The Bromwich integral in eq.(36) can be evaluated by appealing to the theory of complex variables [22-24]. By considering a Hankel-Bromwich path, it can be evaluated as the sum of

*p w*

*a*

*x <sup>a</sup> <sup>t</sup> <sup>a</sup> t t <sup>a</sup>*

> 2 0

*s dx*

2

*a <sup>a</sup>* -

corresponding to a sinusoidal forcing *f*( ) sin( ) *tA t* = *w* ,


0

*<sup>s</sup> <sup>x</sup> st c*

*a a*

*<sup>a</sup>* - =¹ = (34)

¶ - = - <sup>G</sup> ¶ <sup>ò</sup> (32)

*<sup>a</sup>* <sup>=</sup> (33)

= (35)

<sup>=</sup> <sup>+</sup> ò (36)

( ,) ( ,) (,) *tr st uxt u xt u st* = + (37)

0 (,) ( ,) ( ) ( ) *<sup>t</sup> <sup>c</sup> u x uxt t d*

Taking the Laplace transform on both sides of equation. (32) yields

0

This is an ordinary differential equation and can be solved to yield

representing a transient part and a steady state part respectively.

*a*

*uxs*

$$K\_{\alpha}(r,\omega,c\_{0}^{\alpha},\mathbf{x}) = \frac{A\omega}{\pi(r^{2}+\omega^{2})}e^{-\left(\frac{r}{c\_{0}}\right)^{\alpha/2}\mathbf{x}\cos(\pi\alpha/2)}\sin\left((\frac{r}{c\_{0}})^{\alpha/2}\mathbf{x}\sin(\pi\alpha/2)\right)\tag{39}$$

The steady part, ( ,) *st u xt* , arises from the residues of the poles of the integrand in eq. (36) at *s i* = *w* and is given by

$$\mu\_{\rm sf}(\mathbf{x},t) = A e^{-\left(\frac{\omega}{c\_0}\right)^{\alpha/2} \times \cos(\pi\alpha/4)} \sin\omega(t - \frac{\mathbf{x}}{\omega}) \frac{\omega}{c\_0} \Big|^{\alpha/2} \sin(\pi\alpha/4) \tag{40}$$

This completes the formal analysis of the response to sinusoidal forcing of a chain of fractional oscillators in the continuum limit. Numerical results will be presented in the next section.

#### **6. Results and discussion**

The results of the last section in equations (36)-(40) are new. Applications of fractional calculus to oscillation and dissipation problems and solutions to integral and fractional order differential equations can be found in Gorenflo and Mainardi [8] Impulse response functions are also discussed in their work. The same methodology has been extended to sinusoidal forcing of a linear chain in the present work. The response of a chain of fractional oscillators in the continuum limit, when subject to a sinusoidal forcing starting from a quiescent state consists of a transient part and a steady part. The transient part decays in time and approaches zero as *t* ¥. In the limit *a* 2 the transient part vanishes entirely. Furthermore, it exhibits attenuation as a function of distance from the end as indicated by the spatial dependence of the kernel in eq. (39). No simple closed form expressions can be obtained, the only recourse is through numerical integration. Closed form solutions are not available even for simple mechanical models such as the standard linear solid, which are not based on fractional calculus, for the reason that the solutions are deemed mathematically cumbersome as has been discussed by Caputo and Mainardi [5]. These authors have applied methods of fractional calculus to the linear theory of viscoelasticity based on mechanical models and have discussed the propagation of viscoelastic waves[6]. Their work is macroscopic in its approach and is not based on the continuum limit of lattice dynamics. However, some concepts developed by these authors are quite insightful and the results of

the present study can be related to these concepts as follows. Equation (38) can be rewritten by changing the variable *r* 1 / *t* as

$$\mu\_{tr} = \frac{A}{c\_0^{\alpha}} \int\_0^{\infty} R(\tau) e^{-t/\tau} d\tau \tag{41}$$

Microscopic Formulation of Fractional Theory of Viscoelasticity 71

**0.2 0.4 0.6 0.8 1.0 1.2**

 **= 1.7**

 **= 1.9**

 **= 2.0**

 **= 1.7**

**012345**

 **= 1.9**

 **= 2.0**

**Figure 5.** Phase velocity as a function of driving frequency (in arbitrary units).

 **= 1.5**

**0**

**0.6**

**0.8**

**1.0**

**v**

**1.2**

**1.4**

**1.6**

**Figure 4.** Relaxation Spectrum for different values of *a*

 **= 1.5**

**1**

**2**

**R(**

**)**

**3**

**4**

**5**

to bring out clearly the decay in time. It is obvious that the decay in time is characterized by a distribution of relaxation times. The relaxation spectrum is given by

$$R(\tau) = \frac{\omega}{\pi (1 + \omega^2 \tau^2)} e^{-\frac{\mathbf{x}}{\left(\mathbf{c}\_0 \tau\right)^{\alpha/2}} \cos(\pi \alpha/2)} \times \sin(\frac{\mathbf{x}}{\left(\mathbf{c}\_0 \tau\right)^{\alpha/2}} \sin(\pi \alpha/2)) \tag{42}$$

It may be noted that the relaxation spectrum depends also on the distance *x* , which appears in a similarity combination /2 <sup>0</sup> *x c* / ( )*<sup>a</sup> <sup>t</sup>* and that the factor in the exponent, cos( / 2) *pa* , may become negative. Hence caution should be exercised in a strict interpretation in terms of relaxation times. The relaxation spectrum is displayed in figure 4 for several values of the parameter *a* . The values of the other parameters have been chosen to be 0 *x c* = = 1.0, 1.0 and *w* =1.88 for this display.

The steady state solution given by equation (40) can be written as

$$
\mu\_{\rm st}(\mathbf{x}, t) = \frac{A}{c\_0^{\alpha}} e^{-\delta \mathbf{x}} \sin \omega (t - \mathbf{x} / v) \tag{43}
$$

explicitly showing that it as an attenuated wave, propagating with a velocity (phase velocity)

$$v = \frac{\omega c\_0^{\alpha/2}}{\omega^{\alpha/2} \sin(\pi \alpha/4)}\tag{44}$$

and characterized by an absorption parameter

$$\delta = \left(\frac{\omega}{c\_0}\right)^{\alpha/2} \cos(\pi \alpha / 4) \tag{45}$$

In the limit 0 *a d* 2 , and 0 *v c* as can be seen from equations (44) and (45).

Phase velocity as a function of frequency for several values of *a* is shown in Fig. 5. The variation of the absorption parameter as a function of the driving frequency is shown in Fig.6.

**Figure 4.** Relaxation Spectrum for different values of *a*

by changing the variable *r* 1 / *t* as

in a similarity combination /2

and *w* =1.88 for this display.

velocity)

Fig.6.

the present study can be related to these concepts as follows. Equation (38) can be rewritten

0 0

*c*

/2 0

0

*a w*

/2

Phase velocity as a function of frequency for several values of *a* is shown in Fig. 5. The variation of the absorption parameter as a function of the driving frequency is shown in

*a <sup>w</sup> <sup>d</sup> pa* æ ö <sup>ç</sup> <sup>÷</sup> =ç <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

0

*c*

In the limit 0 *a d* 2 , and 0 *v c* as can be seen from equations (44) and (45).

*c*

*v*

*<sup>a</sup> pa <sup>t</sup>*

*<sup>w</sup> <sup>t</sup> pa p wt t*

( )


*x <sup>c</sup> <sup>x</sup> R e*

*<sup>A</sup> u Re d*

to bring out clearly the decay in time. It is obvious that the decay in time is characterized by

cos( /2)

2 2 /2

( ) sin( sin( / 2)) (1 ) ( )

It may be noted that the relaxation spectrum depends also on the distance *x* , which appears

become negative. Hence caution should be exercised in a strict interpretation in terms of relaxation times. The relaxation spectrum is displayed in figure 4 for several values of the parameter *a* . The values of the other parameters have been chosen to be 0 *x c* = = 1.0, 1.0

( , ) sin ( / ) *<sup>x</sup>*

*<sup>A</sup> u xt e t x v*

*d*

explicitly showing that it as an attenuated wave, propagating with a velocity (phase

/2 0 /2 sin( / 4) *c*

cos( / 4)

*a*

*tr*

a distribution of relaxation times. The relaxation spectrum is given by

The steady state solution given by equation (40) can be written as

*st*

and characterized by an absorption parameter

/

*t*

0

*a*

<sup>+</sup> (42)

<sup>0</sup> *x c* / ( )*<sup>a</sup> <sup>t</sup>* and that the factor in the exponent, cos( / 2) *pa* , may

*<sup>a</sup> <sup>w</sup>* - = - (43)

*w pa* <sup>=</sup> (44)

<sup>÷</sup> <sup>ç</sup> <sup>÷</sup> è ø (45)

*c*

*<sup>a</sup> t t* ¥ - <sup>=</sup> ò (41)

( ) *<sup>t</sup>*

**Figure 5.** Phase velocity as a function of driving frequency (in arbitrary units).

**Figure 6.** Absorption parameter *d* as a function of driving Frequency (in arbitrary units).

According to equation. (40), the amplitude of the response wave decays with distance as /2 0 cos( /4) ~ *<sup>x</sup> <sup>c</sup> Ae a <sup>w</sup> pa* æ ö <sup>ç</sup> <sup>÷</sup> -ç <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>÷</sup> çç <sup>÷</sup> è ø and hence the energy as /2 0 2 cos( /4) <sup>2</sup> *<sup>x</sup> <sup>c</sup> A e a <sup>w</sup> pa* æ ö <sup>ç</sup> <sup>÷</sup> - <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup> <sup>÷</sup> çç <sup>÷</sup> è ø . The specific dissipation function can be defined in terms of the distance over which the energy decays by a factor <sup>1</sup> *e* - , which is given by

$$\infty = \frac{c\_0^{\alpha/2}}{2\omega^{\alpha/2}\cos(\pi\alpha/4)}\tag{46}$$

Microscopic Formulation of Fractional Theory of Viscoelasticity 73

(47)

The specific dissipation function is given by <sup>1</sup> *Q*- where *Q* is defined (in analogy with the quality factor 'Q' of harmonic oscillators) as the phase change in radians experienced as the

<sup>1</sup> sin( / 4) tan( / 4) <sup>2</sup> *Q x*

An alternate definition of ' <sup>1</sup> *Q*- ' has been given [6] by Caputo and Mainardi, who introduced in analogy with the optical case, a complex 'refractive index, *r i n n in* = - . The

/2

*v c a*

/2 1 0

*a w pa*

/2

*a <sup>d</sup> <sup>w</sup> pa w w*

0 cos( / 4) *<sup>i</sup>*

<sup>1</sup> <sup>2</sup> 2cot( / 4) *<sup>i</sup>*

Kolsky's definition [26, equation 5.22] of <sup>1</sup> *Q v* ,(2 / *d w* - = , in our notation), also yields the

As explained earlier, the classical theory of viscoelasticity is based on the so called Standard Linear Solid Model characterized by the constitutive relation given in equation (6). This model gives a good qualitative description of the deformation behavior of viscoelastic solids, but has been found to be inadequate to give quantitative account as has been discussed in detail by Gorenflo and Mainardi [8]. More complex mechanical models lead to difficulties in formulation and solution of the resulting complicated differential equations and no satisfactory resolution has been possible [8]. Quantitative agreement is secured only when recourse is made to the so called fractional viscoelastic models. The stress-strain relations are then expressed in terms of fractional order integrals in both creep- and relaxation- representations mentioned earlier and further details can be obtained in reference [9]. In the space-time domain, the approach based on

*c*

<sup>÷</sup> <sup>ç</sup> <sup>÷</sup> è ø

sin( / 4)

*<sup>w</sup>* - = = (48)

*pa* - = = (50)

(49)

real and imaginary parts of the refractive index for the wave in eq. (40) are given by

0

0 0

æ ö <sup>ç</sup> <sup>÷</sup> = = <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

According to Caputo and Mainardi [6], the specific dissipation function is given by

*r*

*c c*

*<sup>n</sup> <sup>Q</sup> n*

*<sup>w</sup> pa pa* æ ö <sup>ç</sup> <sup>÷</sup> = = <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

wave travels the distance *x* in equation. (45) and is given by

/2

*a*

0

*r c*

*n*

*n*

which is equivalent to the expression in equation (47).

same expression as in equation (50).

**6.1. Further discussion** 

and

*c*

<sup>÷</sup> <sup>ç</sup> <sup>÷</sup> è ø

The specific dissipation function is given by <sup>1</sup> *Q*- where *Q* is defined (in analogy with the quality factor 'Q' of harmonic oscillators) as the phase change in radians experienced as the wave travels the distance *x* in equation. (45) and is given by

$$Q = \text{tr}\left(\frac{\omega}{c\_0}\right)^{\alpha/2} \sin(\pi \alpha / 4) = \frac{1}{2} \tan(\pi \alpha / 4) \tag{47}$$

An alternate definition of ' <sup>1</sup> *Q*- ' has been given [6] by Caputo and Mainardi, who introduced in analogy with the optical case, a complex 'refractive index, *r i n n in* = - . The real and imaginary parts of the refractive index for the wave in eq. (40) are given by

$$m\_r = \frac{c\_0}{v} = \frac{\omega^{\alpha/2} \sin(\pi \alpha/4)}{\omega c\_0^{\alpha/2 - 1}}\tag{48}$$

and

72 Viscoelasticity – From Theory to Biological Applications

**0.00**

**0.25**

**0.50**

/2 0


*<sup>x</sup> <sup>c</sup> Ae a <sup>w</sup> pa* æ ö <sup>ç</sup> <sup>÷</sup> -ç <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

~

<sup>1</sup> *e*

cos( /4)

<sup>÷</sup> çç <sup>÷</sup> è ø and hence the energy as

**0.75**

 **= 1.5**

 **= 1.7**

**1.00**

**1.25**

**012345**

 **= 2.0**

 **= 1.9**

/2

2 cos( /4)

*<sup>a</sup> w pa* <sup>=</sup> (46)

<sup>÷</sup> çç <sup>÷</sup> è ø . The specific dissipation

0

<sup>2</sup> *<sup>x</sup> <sup>c</sup> A e a <sup>w</sup> pa* æ ö

<sup>ç</sup> <sup>÷</sup> - <sup>ç</sup> <sup>÷</sup> <sup>ç</sup> <sup>÷</sup>

**Figure 6.** Absorption parameter *d* as a function of driving Frequency (in arbitrary units).

*x*

According to equation. (40), the amplitude of the response wave decays with distance as

function can be defined in terms of the distance over which the energy decays by a factor

/2 0 /2 2 cos( / 4) *c*

*a*

$$n\_i = \frac{c\_{0s}}{\omega} = \frac{c\_0}{\omega} \left(\frac{\omega}{c\_0}\right)^{\alpha/2} \cos(\pi \alpha / 4) \tag{49}$$

According to Caputo and Mainardi [6], the specific dissipation function is given by

$$Q^{-1} = \frac{2n\_i}{n\_r} = 2\cot(\pi\alpha/4)\tag{50}$$

which is equivalent to the expression in equation (47).

Kolsky's definition [26, equation 5.22] of <sup>1</sup> *Q v* ,(2 / *d w* - = , in our notation), also yields the same expression as in equation (50).

#### **6.1. Further discussion**

As explained earlier, the classical theory of viscoelasticity is based on the so called Standard Linear Solid Model characterized by the constitutive relation given in equation (6). This model gives a good qualitative description of the deformation behavior of viscoelastic solids, but has been found to be inadequate to give quantitative account as has been discussed in detail by Gorenflo and Mainardi [8]. More complex mechanical models lead to difficulties in formulation and solution of the resulting complicated differential equations and no satisfactory resolution has been possible [8]. Quantitative agreement is secured only when recourse is made to the so called fractional viscoelastic models. The stress-strain relations are then expressed in terms of fractional order integrals in both creep- and relaxation- representations mentioned earlier and further details can be obtained in reference [9]. In the space-time domain, the approach based on

fractional calculus leads to the propagation of viscoelastic waves and the details can be found in Mainardi[9]. The approach developed in the present work is based on lattice dynamics and hence is a dynamical model. It is microscopic in approach and hence macroscopic properties are described in the so called long wave limit. Further deformation properties which are 'static' in nature are described in the zero frequency limit. As has been shown above the microscopic approach leads to the macroscopic properties in the appropriate limit. It has been noted [26] that the presence of spatial inhomogeneity in polymeric materials gives rise to fractional differential operators in time in the relevant evolution equations, while temporal in-homogeneity leads to fractional differential operators in space. Since only fractional order operation in time has been considered in the present work, there is an implicit in-homogeneity in space which results in a tendency of the particles to cluster and move in a collective fashion. Such collective motion can be considered as elastic waves of very large wavelength, much larger than the scale of the in-homogeneity. This provides a sort of justification for the limiting long wavelength approximation employed.

Microscopic Formulation of Fractional Theory of Viscoelasticity 75

[4] Zener C M (1948) Elasicity and Anelasticity of Metals Chicago: Chicago University

[5] Caputo M and Mainardi F (1971) Linear Models of dissipation in Anelastic solids Riv.

[6] Caputo M and Mainardi F (1971) A new dissipation model based on memory

[8] Gorenflo R and Mainardi F (1997) Fractional Calculus: integral and differential equations of fractional order. In: Carpinteri A and Mainardi F , editors. Fractals and Fractional Calculus in Continuum Mechanics. Wien: Springer Verlag. pp 223-276. [9] Mainardi F (1997) Fractional Calculus: Some basic problems of continuum and statistical Mechanics. In: Carpinteri A and Mainardi F , editors. Fractals and Fractional

[10] Mainardi, F (2010) Fractional calculus and Waves in Linear Viscoelasticity. London:

[12] Born M and Huang K (1954) Dynamical Theory of Crystal Lattices. Oxford,: Clarendon

[13] Maradudin A A, Montroll E M, Weiss G H and Ipotova I P (1971) Theory of Lattice Dynamics in the Harmonic Approximation (Solid State Physics Suppl. 3 2nd edition

[14] Bottger H (1983) Principles of the Theory of Lattice Dynamics. Weiheim: Physik Verlag [15] Askar A (1985) Lattice Dynamical Foundations of Continuum Theories. Singapore:

[16] Narahari Achar B N, Hanneken J W, Enck T and Clarke T (2001) Dynamics of the

[17] Narahari Achar B N, Hanneken J W, and Clarke T (2002) Response Characteristics of the

[18] Narahari Achar B N, Hanneken J W and Clarle T 2004 Damping Characteristics of the Fractional Oscillator *Physica A* 339 311-319. A typographical error in equation. (13) of this reference has been corrected and the corrected equation appears as equation (37) in

[19] Narahari Achar B N, Hanneken J W 2005 Dynamic Response of the Fractional Relaxor-Oscillator to a Harmonic Driving Force *Proc. IDETC/CIE 2005 ASME 2005 Int. Design Engineering Technical Conf. & Computers and Information in Engineering Conf.* (Long Beach,

[20] Narahari Achar B N, Prozny T and Hanneken J W 2007 Linear chain of coupled oscillators: Response dynamics and its continuum limit *Proc. IDET/CIE 2007 ASME 2007 Int. Design Engineering Technical Conf. & Computers and Information in Engineering Conf.* 

this reference. The corrected form of the kernel is given in equation (21a) in reference

missing in equation.(27) of

[7] Glockle W G and Nonenmacher T F (1991) Macromolecules 24: 6426-6434

Calculus in Continuum Mechanics. Wien: Springer Verlag. pp 291-348.

Press.

Press

Nuovo Cimento (Ser. II) 1: 161-198.

Imperial College Press

.New York: Academic Press.

World Scientific.

[16] cited below

Ref. [18]

mechanism Pure and Appl. Geophys. 91: 134-147.

[11] Born M and Von Karman T (1912) Phyzik, Z. 13: 297.

Fractional Oscillator. Physica A 297: 361-367.

Fractional Oscillator *Physica A* 309 275-288 There is a factor *r*

CA , 22-24 September 2005) DETC2005- 84345 pp 1-7

(Las Vegas, Nevada, 4-7 September 2007) DETC2007-35403 pp 1-7

[21] Podlubny L 1999 *Fractional Differential Equations,*(San Diego, CA:Academic Press)
