Fractional calculus and viscoelasticity

Material Information

Fractional calculus and viscoelasticity
Welch, Samuel Williamson John
Publication Date:
Physical Description:
vi, 42 leaves : illustrations ; 29 cm


Subjects / Keywords:
Viscoelasticity -- Mathematical models ( lcsh )
Calculus ( lcsh )
Calculus ( fast )
Viscoelasticity -- Mathematical models ( fast )
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )


Includes bibliographical references (leaf 42).
General Note:
Submitted in partial fulfillment of the requirements for the degree of Master of Science, Department of Mechanical Engineering.
Statement of Responsibility:
by Samuel Williamson John Welch.

Record Information

Source Institution:
University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
18003981 ( OCLC )

Full Text
Samuel Williamson John Welch
B.S., University of Colorado at Denver, 1982
A Thesis submitted to the
Faculty of the graduate school of the
University of Colorado at Denver in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Mechanical Engineering

This thesis for the Master of Science degree by
Samuel Williamson John Welch
has been approved for the
Department of
Mechanical Engineering

Welch, Samuel Williamson John
Fractional Calculus and Viscoelasticity
Thesis directed by Professor Ralph Koeller
The shear relaxation moduli for materials obeying
constitutive equations of the fractional calculus type are
found. Two simple solutions to quasistatic boundary value
problems are then presented using one of these moduli. A method
of fitting experimental data to the relaxation modulus is
discussed with some examples of data fits then presented. The
difficulties of the use of the fractional calculus models in
finite element analysis is then discussed. Alternative
representations of the relaxation moduli are used in a finite
element program and comparison of the finite element results
with the analytic solutions obtained earlier are shown

I would like to thank Professor Koeller for providing
the inspiration to attempt this thesis, Professor Clohessy for
his patience, Professor Trapp for his assistance with the finite
element program and Faye Waitman for her assistance with word

I. INTRODUCTION....................................... 1
II. CONSTITUTIVE EQUATIONS............................. 3
III. ANALYTIC SOLUTIONS................................ 10
Viscoelastic Cylinder.......................... 11
Viscoelastic Wall.............................. 15
IV. FITTING DATA.................................. .. 23
V. FINITE ELEMENT.................................... 32
BIBLIOGRAPHY.......................................... 42

3.1 Radial displacement of a viscoelastic cylinder
under internal pressure.......................... 14
3.2 Drawing of the viscoelastic wall................. 15
3.3 X direction displacement of a viscoelastic wall
pushed into a rigid boundary..................... 21
3.4 Y direction displacement of a viscoelastic wall
pushed into a rigid boundary........"............ 22
4.1 Graph of a seven parameter curve fit to
polyisobutylene data.............................: 28
4.2 Graph of a four parameter curve fit to
polyurethane data................................. 29
4.3 Graph of a seven parameter curve fit to
polyurethane data.................................. 30
4.4 Graph comparing ..the seven parameter
Mittag-Leffler fit with an exponential fit.......... 31
5.1 Graph comparing the finite element solution with
the analytic solution for the radial displacement
of a viscoelastic cylinder under internal
pressure.......................................... 39
5.2 Graph comparing the finite element solution with
the analytic solution for the x direction
displacement of a viscoelastic wall pushed
against a rigid' boundary........................... 40
5.3 Graph comparing the finite element solution with
the analytic solution for the y direction
displacement of a viscoelastic wall pushed
against a rigid boundary............................ 41

This work is an attempt to investigate some of the
practical aspects of the use of fractional calculus in the linear
theory of viscoelasticity. The materials considered will be
homogeneous and isotropic while the problems considered will be
two dimensional and quasistatic.
First, some relationships between fractional calculus and
the linear theory of viscoelasticity will be examined. Some
simple solutions will then be found using these relationships.
The problem of fitting experimental data to the functional forms
required by the above theory is considered next and some examples
are provided. The last topic considered will be the use of the
fractional calculus models in a finite element program. Test
problems for the program will be the simple solutions obtained
Viscoelastic behavior is being studied in such diverse
fields as polymer science, biomechanics and geology. One of the
problems facing analysts in these fields is how to effectively
characterize the viscoelastic behavior of materials. With the. use
of the fractional calculus models materials can be characterized

by fewer parameters. The cost of this is increased mathematical
complexity. In many cases it may be necessary to forego the
fractional models for simpler models. The choice of models may be
decided by how well one can fit experimental data with the
various forms for the material moduli. In many cases one need not
consider what model he is using but simply concern himself with
fitting the data to functions.
There are situations where the analyst is aided by the
compactness of the various moduli expressed as Mittag-Leffler
functions. There are other situations where the added
mathematical complexity of the fractional models makes them
impractical for use. The two simple analytic solutions obtained
in chapter three are examples of the former while the development
of a finite element program considered in chapter five is an
example of the latter.

In this chapter we state the form of the integral
constitutive equations to be used. These constitutive equations
will be used explicitly in the finite element program. We also
will define a specific fractional model for these constitutive
equations. The result will be the determination of the form for
the shear relaxation modulus. If need be, a similar result could
be derived for the bulk relaxation modulus. We will also state a
dual form of the constitutive equation and show that the
viscoelastic correspondence principle can still be used for the
fractional models.
The general form of the constitutive equation for a
linear viscoelastic body may be written in terms of Stieltjes
convolutions as
cj i j C t)
Now, since these functions are in the Heaviside class
Oij(t) = Gijkl(t)eki(0+)

Use of the well known theorem for representation of a fourth
order isotropic tensor leads to the form
5ij {[k(t)-2_u(t)]enn(0+) +
0+ 3 dx
Note that in this form p(t) is the relaxation modulus in
pure shear and k(t) is the relaxation modulus in pure
dilatation. Note also that this is a form for the constitutive
equation of a linear isotropic viscoelastic body. In what
follows we will consider a state of pure shear and using a four
parameter fractional model we will find what form the relaxation
modulus must take.
We consider a material that obeys the following
constitutive equation in pure shear.
E0(D^+1/tj)eij(t) = (DB+1/tf)cij(t) i*j (1.2)
Tl = tid+Eo/Ej
and D
indicates fractional
i .e.
Def(t) = DDB_1f(t)

For a detailed discussion of the use fractional calculus in
viscoelasticity along with the above model see [2,31
Use of the Laplace transform on (1.2) yields
- ,_x , En(SB+1/T?) r /.x
aij(s) = o * eij(s)
Where the Laplace transform pair
£[Def(t)] = s6f(s)
has been used. In pure shear the constitutive equation (1.1)
ij )
2y(t)eij(0+) + 2
Jo+ d-r
Taking the Laplace transform of this yields
aij(s) = 2sy(s)Iij(s)
Comparison of (1.3) and (1.4). indicates
y(s) =
Inversion of this yields

u(t) = 1_ E|_ EgC-Ct/tj)2] + EnE,
2 E j+E0 2(E0+Ei)
Where Eg(x) is the Mittag-Leffler function i.e.
Eg(x) = E x^n
i=o rd+en)
To simplify the notation let
G0 = 1_ E^
2 Eo+E,
and Gi = L JL
2 E0+E1
And we have, finally
y(t) = G0 + G^gC-Ct/tJ6] (1.5)
It should be noted that a more direct way to obtain this form
would be to set e-jj(t) = h(t), where h(t) is the Heaviside unit
step function, and solve (1.2) for Oij(t) which by (1.M) is
twice the relaxation modulus.
If we wish to use a model with more parameters the
simplest approach is to assume a relaxation modulus of the form
U(t) = G0 + E GiEg.E-U/ti)^] (1.6)
i='1 1
This can then be used directly in (1.1). If the fractional model
in the form of (1.2) is desired it can be found by a reverse
process to that outlined above. This method is desirable as one
must often seek a form of (1.6) that will fit relaxation data

obtained by experiment. This will be considered in detail in a
later chapter.
The dual form of the constitutive equation will now be
arrived at by using the commutation property Stieltjes
convolutions. The commutation property is stated as
This leads directly to the dual form of the constitutive
Oij(t) = 6 i j {[k(0 + )-2p(0+) ]enn(t) + enn(t--r) [8k( t)-2_ 3y( t ) ]dx }
3 Jo+ 3t 3 3t
+2p(0+)eij(t)+2 e^j(t-t)3u(T)dr
It is customary, once the relaxation modulus is known,
to find the creep compliance. The creep compliance and the
relaxation modulus are related in the Laplace transform domain
G(s)J(s) = 1/s2
Where G(t) is twice the shear relaxation modulus and J(t) is the
creep compliance. Application of this to (1.5) yields the creep

J(t) = 1_ [ 1-_Gi_ EgE-Ct/Ti)8]
G0 G1+G0
T? = tf Gp+G,
Models with more parameters, as in (1.6), become
difficult to apply this operation to. We note that the creep
compliance is useful if the constitutive equation is to be used
in the form
^ijkl (t-T)d . 00
Now consider the form of the constitutive equation
(1.1). Taking, the Laplace transform we have
Oij(s) = s[k(s)-2 p(s)]6ijinn(s) + 2sy(s)eij(s)
This form of the constitutive equation is the basis for the
viscoelastic elastic correspondence principal. Once the elastic
solution is known we replace the stresses, strains and
displacements with Laplace transformed variables. We then
replace the elastic moduli with Laplace transformed moduli
multiplied by s to obtain the Laplace transformed solution of
the corresponding viscoelastic problem. This principal along

with restrictions on the boundary conditions is discussed in
many references [1,7,8]. This principal works equally well with
the Fourier transform. The key point here is that use of the
fractional calculus models does not remove this useful tool.
The viscoelastic moduli used in this work are the shear
modulus u and the bulk modulus k. If one wishes to characterize
a material with moduli other than these the procedure is similar
to that in elasticity. The most convenient way to do this is in
the Laplace transform domain. This is also done in many
references. The shear modulus was chosen to work with because
data was available. It is also worth noting that many materials
behave elastically under pure dilatation thus y and k seemed a
good choice.

In this chapter are presented analytic solutions for
two simple problems. The outstanding feature of these solutions
is their relative compactness. If the relaxation modulus is
represented by a series of exponentials these solutions would
become much more cumbersome. The advantage of using a series of
exponentials is that more complicated problems may be handled in
a systematic way with the theory of residue calculus. The method
of inversion is straightforward as the resulting Laplace
transformed solution is expressed as a ratio of polynomials.
Residue calculus is difficult with the fractional models because
of the existence of branch points. There are, however, examples
in the literature where this has been handled successfully.
There is some motivation in the literature for the parameter 0
to be one half [4]. If this is so then some of the difficulties
involved in obtaining analytic solutions are more easily
overcome. None of the data used in this work appeared to have
this desirable characteristic. Another simplification would
arise in models with more then one Mittag-Leffler function if
the different B's were equal. Again, none of the data used had
this characteristic.

Viscoelastic Cylinder
Consider a thick circular cylinder with inside radius a,
outside radius b and an internal pressure PjCt) = Pih(t) where
h(t) is the Heaviside unit step function. Assume plane strain
conditions and negligible body forces. We seek the quasistatic
viscoelastic response for a material which behaves elastically
under dilatation and has a shear relaxation modulus of the form
u(t) = G0 + GjEgC-Ct/tJ3]
The elastic solution is given by
rr = PiaZ [ 1 b2/r2 ] ,
ee = Pia2 [ 1 + b2/r2 ]
a77 = piaZ and
ur = Pi*2 [ 3r' + b2 ]
2(a2-b2) 3k+y yr
We note that the stresses are independent of the
material properties and thus are the same for both the elastic
and the viscoelastic cases. To obtain the viscoelastic solution
for the displacement we make use of the viscoelastic
correspondence principle. We start by substituting the following
into the elastic solution

ur(s) for ur ,
sk(s) for k ,
sp(s) for p and
P^Cs) for Pi
Where the superscript bar indicates the Laplace transformed
variable. For this case we have
= k
= G0 +
Pi-(s) = Pi/S
Where the Laplace transform of the Mittag-Leffler function has
been used. Substituting these into the expression for the
displacement we obtain the Laplace transformed displacement
ur(s) = pi &2
3r(sB+1/tJ )
b2(s6+1/t? )
rCCGo + G^s^+Go/tf-]
The Laplace transform of the Mittag-Leffler function is
£[ EgE-Ct/tJ3] ] =
6+1 /t

Another useful Laplace transform pair is
£[ 1-EgC-Ct/t!)3] ] = 1
Use of these transform pairs allows us to invert the expression
for the Laplace transformed displacement back to the time
domain. The result is
_3r_.[ 1
+ b2 [ 1
GiEg[-(t/Tl)g] ]
G!Eg[-(t/T2)B] ]
1 _ 3k+G0
t B tfOk+Go+Gj
t^Go + Gj
Note that as t--> we recover the elastic solution for the
Up = Pi?2 \ 3r + b2 ]
2(a2-b2) 3k+G0 G0r
This solution is plotted on the following page using material
properties obtained from a curve fit in a later chapter.

1 i*
figure 3.1
Analytic solution for cylinder with internal pressure.
r = .1 inch, b = .15 inch, a = .05 inch, = 100 psi
and p(t) = G0 + G^gO (t/^) ] where G0, Gp 6 and t1 are
found in figure H.2 and the bulk modulus is 10000 psi.

Viscoelastic Wall
Consider the cross section of a viscoelastic wall shown
The wall is pushed into a rigid boundary by a uniform
pressure P0(t) = P0h(t). The wall is considered very long in the
z direction and conditions of plane strain are assumed.
e13 = e2 3 = e33 = 0
13 = 23= 0
The surface at y=0 is considered lubricated. The
material is such that it.has an elastic bulk modulus and a shear
modulus of the form
y(t) = G0 + G^gH-(t/t!

We neglect body forces and assume quasistatic
conditions. Under these assumptions the equations of motion
3a j j 3a x 2
lx" + ~3y"" 0
3o j 2 3o 2 2
~3x" + TT =
The boundary conditions are
Lateral faces
Oji(a ,y,z) = a12(a,y,z) = 0
Rigid boundary
a12(x,0,z) = 0
u2(x,0,z) = 0
Bottom face
a22(x,-L,z) *3 Pq
a12(x,-L,z) = 0
We assume elastic displacements of the form
Uj = ax + by
u2 = cx. + dy

Application of the boundary conditions and use of the
constitutive equations for linear isotropic elastic solids leads
to the expression
0 2 2 = P 0
0ij = 0 i,j ^ 2
u, = (3k-2y)P0x
Ui = (3k+Hy)P0y
Now, since we have satisfied the equations of motion the
boundary conditions and the equations of compatibility we have
the elastic solution. To obtain the viscoelastic solution we
again use the viscoelastic correspondence principal.
sk(s) = k
Sy(s) = G0 + GisB
P0(s) = P0/s
Note that the elastic stresses are independent of the
material properties hence the stresses in the viscoelastic case
are the same as those in the elastic case. We now write the
Laplace transformed displacements substituting in the above
expressions for the Laplace transformed quantities.

Uj (s) =
2g S 2
s +2as +b
s20"1 (3K-2GO-2GJ
+ s (6k-MGn-2G,)
+ (3k-2G0)
u2(s) = -'l'(y)
28 8 ^2
s +2as +b
s2B '(Bk+HGo+'ilGj)
+ a___(Sk+SGo+HGj)
^ 1
+ (3k+UG0)
i|)(z-) =
b2 =
a = (G0+G1)(3k+G0)+(3k+G0+G1)
Now define

Now, the following Laplace transform pairs are needed
2 0 6,2
s +2as +b
^ EgC-Ct/xJ6] + EB[-(t/T2)3]
2 c 2 c

20 0.2
s +2as +b
J_ EgC-Ct/xJ6] +J_ E6[-(t/x2)e]
2 c 2 c
20 0,2
s +2as +b
_L EgC-Ct/xJ6] +J_Eg[-(t/x2)6]
2c(c-a). 2c(c-a) a2-cz
These transform pairs are obtained.using a technique
similar to partial fraction expansions. The time domain
displacements are
ux(t) =
^ Eg[-(t/Tl)e] + Eg[-(t/x2)B]
2 c
2 c
+ ip(x) (6k-^G0-2G1)
2ct i
+ \|)(x) (3k-2G0)
EgH-tt/xJ0] - Eg[-(t/x2)6]
Eg[ (t/xj^] + Eg[-(t/x2)6] + 1

u2 (t)
2^ EglXt/O6]
2 c
+ 211 Eg[-(t/T2)
2 c
- ^(y)(6k+8G0^G1)
- ^(x)(3k+MG0)
Eg[-(t/Tl)B] + Eg[-(t/T2)B] + J_
_ 2c(c-a) 2c(c+a) b2
Note that as t-> we recover the elastic solution
u, = P0x(3k-2G0)
u, = P0y(3k+MG0)
This solution is plotted on the following two pages using
material properties obtained from a curve fit in a later

figure 3.3
Analytic solution for the viscoelastic wall with x = 1
inch, L = 2 inches, a = 1 inch and P = 300 psi. The relaxation
function is
p(t) = G0 + G^gC-Ct/tJ6]
where G0, Glf 3 and tj are found in figure *1.2 and the bulk
modulus is 10000 psi.

Analytic solution for viscoelastic wall at y = 1 inch
and all other parameters the same as the previous figure..

In this chapter the method of fitting data obtained by
measurement to the Mittag-Leffler function representation of the
relaxation modulus is described. For ease with notation the four
parameter model will be considered with the generalization to
models with more parameters being obvious. We will consider the
case where we have obtained data from a relaxation test in pure
Let {(yi,tj);i=1...n} be the obtained data. We seek to
fit this data to a four parameter model of the form
U(t) = G0 + GiEgC-Ct/t!)3]
Where Eg is the Mittag-Leffler function with parameter (3. Define
the least squares error as
n , ,2
Eis =i=i iPi"y(tj_)} r
Els Jt {yi-G^EgC-Ui/t,)6]}2
We note that G0 may be fixed from the data as t->. Minimization
of Ejg with respect to the three remaining parameters yields the

three equations.
=il1 Uyi-Go-G^gC-Cti/tJ6])
=i|1 KPi-Go"G1Eg[-(ti/t1 )6])
6 ,
^EgE-Cti/tJ ]} = o
6 ,
EgE-Cti/tJ ]} = 0
3 ,
L.EgE-Cti/t!) ]} = 0
Noting that this is a system of nonlinear equations we attempt
to find a solution using Newton's method. This scheme may be
written as : .
xk xk 1 j(x^ 1) ^F(xk))
M 'Fi(x) '
X ={ Gl ' l = , F2(x)
' W Fs(x) > i
/ \
3F, 3F, 3F,
36 3GX 3tj
J(x) = ^ 3F, 3F? 3F,
36 3GX 3tx
3F, 3F, 3F,
l 93 3GX 3ti J
and the superscript k denotes the kth iteration. The method is
implemented by first solving the linear system of equations
J(xk1)z = F(xk1)

The kth iterate is then
xk xk-1 z
The first guesses for Gj and tt are obtained with the
aid of- a graphics computer. For models with a larger number of
parameters this can be time consuming but in view of the
difficulty involved in programming a general first iteration
algorithm this method seemed the most desirable.
In solving the linear system (3.1) care must be taken
that the Jacobian matrix does not become algorithmically
singular. The matrix is scaled in such a way that the largest
element in each row is of the same order of magnitude. This is
a necessary expedient if one wishes to make use of the condition
number returned by most good linear equations software. This
number gives one an idea of how algorithmically singular a
matrix may be. When the Jacobian matrix did become
algorithmically singular it was found that a slight perturbation
of the current iterate sometimes corrected this and the
iteration could proceed. If this did not work then a new first
guess had to be tried. For a complete discussion of the
condition number estimate see [9].
The derivatives required to fill the Jacobian matrix
are, with the exception of those with respect to G!, calculated
numerically. The analytic expressions for these derivatives are
easily obtained but were not used in the program. The reason for

this is that an investigation into the asymptotic properties of
the various derivatives of the Mittag-Leffler function may be
needed and are beyond the scope of this work. The step sizes
needed for the various derivatives are, in general, different
for different sets of data. These step sizes are found by
evaluating the derivatives for a wide range of step sizes and
selecting step sizes that are approximately in the middle of a
range of stability.
It should be mentioned that the Mittag-Leffler function
is an expensive function to evaluate. The above procedure for
finding the best fit to the data can involve quite a bit of
computer time. One way to help alleviate this is to use the same
Jacobian matrix for a number of iterations. It was also found
that a necessary step for convergence with some data is to use
logarithmic scaling.
On the following three pages are three plots showing the
results of fitting both a four parameter model and a seven
parameter model to relaxation data. The first plot is using data
from Tobolsky and Catsiff. The outstanding feature of this fit
is that fifteen decades of data are fit reasonably, well with two
Mittag-Leffler functions. It should be noted that this material
is essentially a viscoelastic fluid thus G0 is zero. The second
plot shows the result of fitting a four parameter model to the
relaxation function found in Christensen. The outstanding
feature here is that one Mittag-Leffler function fits the data

over a range in which eight exponentials were required. This
relaxation function is used in both the analytic solutions
obtained earlier and in the finite element solutions in a later
chapter. The third plot shows the result of fitting a seven
parameter model to the same data. It is difficult to discuss
accuracy here as what we have is a curve fit to a curve fit. The
fourth plot is a blown up view of the previous plot. The
outstanding feature to note here is the smoothness of the
Mittag-Leffler functions as compared to the roughness of the
In conclusion we note that an advantage to using the
Mittag-Leffler function representation for the relaxation moduli
is that the material properties can be characterized by a fewer
number of parameters. It is believed that the smoothness of
these functions over many decades results in a better fit to the
data. For a discussion of the physical basis of these models see

figure 4.1
Seven parameter fit to data from Tobolsky and Catsiff
y(t) = G0 + GjEgjC-Ct/t!)31] + G2Eg2[-(t/t2)6z]
31 =.6474 B2=-^^96
G :=2.876x10 6 G2=7.2975x10 psi
t !=2.96x10~12 t2=.359 sec.

5 -
4 -
3 '
2 -
1 -
- 1 2
10 10 10
log10(t) (seconds)
figure 4.2
Four parameter fit to Christensens [1] relaxation function.
p(t) = G0 + G1EB[-(t/t1)e]
G0 = 500 psi
Gj = 2500 psi
S = .323
tx = 1.56x10 3 sec

y (t)/G0
5 -
-H -1 2
10 10 10
logio^t) (seconds)
figure 4.3
Seven parameter fit to Christensens [1] relaxation function.
y(t) = G0 + G^g^-Ct/tJ31] + G2Eg2[- (t/t2 )3z]
G0 = 500 psi
G! = 1982 G2 = 1666 psi
0! = .582 02 -375
tj = 1.56x10 t2 = .0109 sec

Christensen's fit
7 parameter Mittag-Leffler fit
Blown up portion of Christensen's fit along with the
seven parameter Mittag-Leffler function fit.

In this chapter the development of a linear viscoelastic
plane strain finite element program will be considered. At an
intermediate point in the development the disadvantage of using
the Mittag-Leffler function representation of the relaxation
modulus will be apparent. The relaxation modulus will therefore
be represented by a series of exponentials. The materials
considered will be isotropic and behave elastically in their
bulk behavior.
We seek approximate solutions to the equations of motion
9ij + p0bi
Subject to the following assumptions
a) Plane strain conditions
b) Quasistatic motion
c) Negligible body forces
The equations of motion then simplify to

The finite element approximations will be formed using
triangular elements with three nodes. This corresponds to
assuming both a linear displacement profile and constant stress
across each element. Galerkin's method will be used to obtain
the approximation equations to be solved. We have
Z 3ij 4 =0 m = 1......M M = 2N
. R 3xj 1
Where the are the linear basis functions, N is the number of
nodes, M is the number of basis functions and R is the region of
interest. Use of the divergence theorem in two dimensions yields

. S
Where S is the bounding curve of the region R. This step is
necessary as the basis functions chosen do not have non-trivial
second derivatives as would be required in the original form of
the Galerkin equations. The stress boundary conditions now
appear naturally in terms of the line integral. Displacement
boundary conditions will be implemented by manipulation of the
final system of equations. A result of using linear basis
functions is that the integrals over the region are simple to
evaluate hence no integral approximation method is needed.
The form of the constitutive equation to be used is
= 0

cj^j(t) 61j {kenn(t)- 2p(t)enn(0+) } + (0 + )]j(fc)
26ij y(t-T)denn(T)dT
3 Jo+
u(t-T)^eij (tMt
eij = + 9uJ)
2 9xj 3xi
Now, introduce the approximation
r ux(xr / / s (p (x)
** 2N
< ^ = E An(t) ^
n-1 n
k u2(x)i k 2(x) j
Where = 0 n even
5 = 0 n odd
The basis functions are normalized such that, for
example, at the coordinates for node 1 J = 1, 2 = 0 and
ux = Aj. Substituting (5.3) into (5.2) and the result into (5.1)
we obtain 2N equations for the 2N unknown displacements. In what
follows we will consider what is necessary to advance the scheme

in time. This is the point at which it is convenient to compare
the effect that the two forms for the relaxation modulus have on
the program development.
Consider the integral
. 0+
When t = tn we have
f^n r^n
y(tn-x)3eijdt = u(tn-t)3eijdi+
0+ 9^ t '
0+ -
)3fi jd-r
- 9t
We seek to write the integral
In terms of the integral
. 0+ 3t
This is a critical step because it allows us to use the
integrals up to the previous time step in calculating the
integrals for the current time step. One alternative is to store
the displacements and use them to calculate the integrals over
the full time range at each time step. This alternative is

undesirable as a large amount of memory would be required and a
large number of computations would have to be performed for each
time step.
We now need consider what must be done if the Mittag-
Leffler function representation of the relaxation modulus is to
be used. Using a Taylor series we write
p(tn-T) *= p(tn_i-T) + At9yI, + o(At2)
§u n-1"T
We note that this requires accumulating two integrals and the
integral arising from the second term on the right will also
have to be approximated from the previous time step. This was
tried with limited success. A higher order scheme could be tried
but three difficulties immediately come to mind. First, we would
be approximating higher order derivatives from functions
obtained from a curve fit. Second, the second derivative of the
relaxation modulus is fairly large and at least two more orders
of accuracy seem to be needed in the Taylor series
approximation. Third any higher order scheme would involve
accumulating many integrals increasing greatly the memory
requirements of the program. If there were no other way these
disadvantages could be worked around. Another alternative would
be to try the dual form of the constitutive equation. This was
not tried and possible difficulties would be the need to take
derivatives of the strains in the same Taylor series

approximations as those noted above. This alternative needs to
be tried. A third alternative, chosen here, is to approximate
the relaxation modulus by a series of exponentials, i.e.
y(t) = G0 +.Z1Giexp(-t/ti)
The integral
. 0+
Is now easily written in terms of the integral
. 0+ 3t
Note also that we only need accumulate one integral for this
Having described the method used to accumulate the
required integrals the integral from tn to tn_-| will now be
discussed. We have
There are a number of simple ways to approximate this integral.
We will assume the strain is linear over the interval leading to

Jtn-i 3t
ei j (tn)~ei j (^n-l )
'n n
{go+Z G1exp[-(tn-T)/ti]}dx
So at t = tn we have calculated
aij (tn) = Sij (1 ^n
~26ij y(tn-x)9enn(x)dx
3 o+ 9t
+2 u(tn-x)9eij(x)dx
Jo+ 9t
This is now substituted into (5.1) and we have completed the
development of the finite element approximations. On the next
three pages are plots of the finite element solutions vs. the
analytic solutions obtained earlier.

t (seconds)
figure 5.1
Comparison of the finite element solution and the analytic
solution obtained earlier for the viscoelastic cylinder. The '+'
indicates the finite element solution.

figure 5.2
Comparison of the
solution obtained earlier
finite element solution and the analytic
for the viscoelastic wall.

figure 5.3
Comparison of the finite element solution and the analytic
solution obtained earlier for the viscoelastic wall.

1. Christensen, R. M.: Theory of Viscoelasticity, An
Introduction, Academic Press, New York, 1971.
2. Koeller, R. C.: Application of Fractional Calculus to the
Theory of Viscoelasticity, ASME Journal of Applied
Mechanics, 51 299-307 (1984)
3. Koeller, R. C.: Polynomial Operators, Stieltjes Convolution,
and Fractional Calculus in Hereditary Mechanics, Acta
Mechanica, 58 251-264 (1986)
4. Bagley, R. L., and Torvik, P. J.: A Theoretical basis for
the Applicatiion of Fractional Calculus to Viscoelasticity,
Journal of Rheology, Vol. 27, No. 3, 201-210 (1983)
5. Bagley, R. L., and Torvik, P. J.: A Generalized Derivative.
Model for an Elastomer Damper, Shock Vibr. Bull., No. 49,
Part 2, Sept. 1979, 135-143
6. Tobolsky, A. V., and Catsiff, E.: Elastoviscous Properties
of Polyisobutylene (and Other Amorphous Polymers) from
Stress Relaxation Studies. IX. A Summary of Results, Journal
of Ploymer Science, Vol XIX, 111-121 (1956)
7. Gurtin, M. E., and Sternberg, E.: On the Linear Theory of
Viscoelastitity, Arch. Ration. Mech. and Anal., Vol. 11,
291-356 (1962)
8. Rabotnov, Yu N.: Elements of Hereditary Solid Mechanics, Mir
Publishers, Moscow, 1980
9. Dongarra, J. J., Mole'r.C. B., Bunch, J. R., and Stewart,
G. W.: Linpack Users' Guide, SIAM Publishers, Philadelphia