NASA / TM-2000-210123

ARL-TR-2206

Determining a Prony Series for a

Viscoelastic Material From Time Varying
Strain Data

Tzikang Chen
U.S. Army Research Laboratory
Vehicle Technology Directorate
Langley Research Center, Hampton, Virginia

May 2000
 The NASA STI Program Office ... in Profile

CONFERENCE PUBLICATION.
Since its founding, NASA has been
dedicated to the advancement of Collected papers from scientific and
aeronautics and space science. The NASA technical conferences, symposia,
Scientific and Technical Information (STI) seminars, or other meetings sponsored
Program Office plays a key part in helping or co-sponsored by NASA.
NASA maintain this important role.
SPECIAL PUBLICATION. Scientific,
The NASA STI Program Office is operated technical, or historical information from
by Langley Research Center, the lead center NASA programs, projects, and missions,
for NASA's scientific and technical often concerned with subjects having
information. The NASA STI Program Office substantial public interest.
provides access to the NASA STI Database,
the largest collection of aeronautical and TECHNICAL TRANSLATION. English-
space science STI in the world. The Program language translations of foreign
Office is also NASA's institutional scientific and technical material
mechanism for disseminating the results of pertinent to NASA's mission.
its research and development activities.
These results are published by NASA in the Specialized services that complement the
NASA STI Report Series, which includes the STI Program Office's diverse offerings
following report types: include creating custom thesauri, building
customized databases, organizing and
TECHNICAL PUBLICATION. Reports publishing research results ... even
of completed research or a major providing videos.
significant phase of research that
For more information about the NASA STI
present the results of NASA programs
and include extensive data or theoretical Program Office, see the following:
analysis. Includes compilations of
significant scientific and technical data • Access the NASA STI Program Home
and information deemed to be of Page at http'//www.sti.nasa.gov
continuing reference value. NASA
counterpart of peer-reviewed formal • E-mail your question via the Internet to
professional papers, but having less [email protected]
stringent limitations on manuscript
length and extent of graphic • Fax your question to the NASA STI
presentations. Help Desk at (301) 621-0134

TECHNICAL MEMORANDUM. • Phone the NASA STI Help Desk at

Scientific and technical findings that are (301) 621-0390
preliminary or of specialized interest,
Write to:
e.g., quick release reports, working
papers, and bibliographies that contain NASA STI Help Desk
minimal annotation. Does not contain NASA Center for AeroSpace Information
7121 Standard Drive
extensive analysis.
Hanover, MD 21076-1320
CONTRACTOR REPORT. Scientific and
technical findings by NASA-sponsored
contractors and grantees.
NASA / TM-2000-210123
ARL-TR-2206

Determining a Prony Series for a

Viscoelastic Material From Time Varying
Strain Data

Tzikang Chen
U.S. Army Research Laboratory
Vehicle Technology Directorate
Langley Research Center, Hampton, Virginia

National Aeronautics and

Space Administration

Langley Research Center

Hampton, Virginia 23681-2199

May 2000
The use of trademarks or names of manufacturers in the report is for accurate reporting and does not constitute an

and°fficialspaceend°rsement'Administrationeitheror
theeXpressedu.s.
Army.°r
implied, of such products or manufacturers by the National Aeronautics

Available from:

NASA Center for AeroSpace Information (CASI) National Technical Information Service (NTIS)
7121 Standard Drive 5285 Port Royal Road
Hanover, MD 21076-1320 Springfield, VA 22161-2171
(301) 621-0390 (703) 605-6000
 Determining a Prony Series for a Viscoelastic Material from Time
Varying Strain Data

ABSTRACT

In this study a method of determining the coefficients in a Prony series representation

of a viscoelastic modulus from rate dependent data is presented. Load versus time test data

for a sequence of different rate loading segments is least-squares fitted to a Prony series

hereditary integral model of the material tested. A nonlinear least squares regression

algorithm is employed. The measured data includes ramp loading, relaxation, and unloading

stress-strain data. The resulting Prony series, which captures strain rate loading and

unloading effects, produces an excellent fit to the complex loading sequence.

KEY WORDS: hereditary integral, viscoelasticity, weighted nonlinear regression, Prony

series, multiple loading segments

INTRODUCTION

In order to determine the time dependent stress - strain state in a linear viscoelastic

material, under an arbitrary loading process, the deformation history must be considered.

The time dependent constitutive equations of a solid viscoelastic material include these

history effects. The load (stress) and displacement (strain) history, the loading rate

(displacement rate), and time of load application on the specimen are all needed to determine

the constants in the constitutive equations. A common form for these constitutive equations

N
employs a Prony series (i.e., a series of the form _, _i "e-t _:i ).
i=1

Creep and relaxation tests are most commonly used to determine the viscoelastic

material properties, see Figure 1. In ideal relaxation and creep tests, the displacements or

loads are applied to the specimen instantly. In the real test, especially for a large structural

component, limitations of the testing equipment result in a relatively low strain rate and long
periodof loading. The responseduringthe periodof loadingis typically ignored,andonly
the dataobtainedduringthe periodof constantdisplacementor constantload areusedto
determinethematerialproperties.1'2Ignoringthis long loadingperiodandthe strainrate
effectsin the datareductionintroducesadditionalerrorsin the determinationof thematerial
properties.
Therearenumerousmethodsfor determiningthe Pronyseriesfrom relaxationand/or
creepdata. An earlymethod3involvedconstructinglog-logplots of relaxationdatain
which straightline approximationsfor the dataonthe log-log graphyield the time constants
(i.e. r i 's) from the slopes, and the exponential coefficients (i.e. o_i 's) are obtained from the

intercepts. Other methods have also been employed. For example, Johnson and Quigley 4

determined a relaxation time constant for a nonlinear model which is similar to a one-term

Prony series model. They minimized a least-squares error measure, of the difference

between the nonlinear model and measured data, by iteratively integrating (numerically) an

internal variable equation. When attempting to determine relaxation time constants for

higher fidelity nonlinear models, Johnson, et al. 5 employed trial and error procedures,

similar to early linear methods, 3 due to the complexity of the resulting nonlinear least-

squares problem. More recently, a few authors 6 8 employed nonlinear optimization methods

to obtain a high quality Prony series representation of relaxation data with a minimum

number of terms in the series. The viscoelastic model can also be formulated in differential

form. This is becoming popular recently 9 11 since the differential models can be effectively

incorporated into finite element algorithms. When using these internal variable methods,

each Prony series term is associated with a material internal state variable. In the discrete

finite element model, each term in the Prony series adds a substantial number of global

variables. Thus, a short Prony series, which can accurately represent the material, is

desirable. Nonlinear regression methods can help with determining a short and accurate

Prony series.

The purpose of this paper is to present a method for including the loading and

unloading data, along with the relaxation data, in a nonlinear regression analysis to obtain the

Prony series. The resulting viscoelastic material model is then capable of simulating the

loading segments as well as the relaxation segments. This is an improvement when modeling

hysteretic effects is important. The analytical solution for loading and/or unloading is

determined herein and employed in a nonlinear regression analysis to determine the Prony

series. In addition, data weighting functions are investigated and are shown to improve the

2
fit in thebeginningof relaxationperiod. Again,this methodallowsall the measureddatato
beincludedandresultsin animprovedconstitutivemodel.

Hereditary Integrals for linear viscoelasticity

Detailed descriptions of linear viscoelasticity can be found in the literature. 1'2 An

overview of linear viscoelasticity is provided here in order to introduce the hereditary

integral method which is used below to determine the analytical solution for the loading and

unloading segments. Linear viscoelastic constitutive models are represented by simple

physical models composed of springs and dashpots. The spring is the linear-elastic

component, and its constitutive equation is

o- = E. (1)
The dashpot is the viscous component, and its constitutive equation is

DE
O" =T] .-- (2)
Ot

where r/is the viscosity constant. Linear viscoelastic constitutive models are constructed by

superimposing components with constitutive equations given by equations (1) and (2).

Since the mechanical response of the dashpot is time dependent, the behavior of a

viscoelastic material that is modeled by parallel and/or series combinations of springs and

dashpots is also time dependent.

The creep test consists of a constant stress, C_o, applied to a specimen for a period of

time while its strain is recorded (Figure la). In a relaxation test, the specimen's strain, eo, is

held constant for a period of time while the stress is recorded (Figure lb). In Figure 1, eo and

C_oare the initial strain and stress, respectively. For the relaxation test, a constitutive relation

for the period of constant strain can be written as follows:

(7(t) = Y(t). eo (3)

where Y(t) is a relaxation function. When the material is assumed to be a general Maxwell

solid, the relaxation function is typically modeled with a Prony series as follows,
 tl

Y(t)=Eo(1- ZPi "(1-e-t/ri))

i=1

(4)

where:

Pi is the i'th Prony constant ( i =1,2 ..... )

is the i'th Prony retardation time constant ( i = 1,2 .... )

Eo is the instantaneous modulus of the material

For time t = O, Y(O) = Eo and for t = oo, y(oo) = Eo (1-Xpi).

In the case of a creep test, a creep compliance function, J(t), is defined as follows.

e(t) =J(t). cr o (5)

The compliance function is then determined by procedures analogous to those described

above.

To determine the stress state in a viscoelastic material at a given time, the

deformation history must be considered. For linear viscoelastic materials, a superposition of

hereditary integrals describes the time dependent response a. If a specimen is load free prior

to the time t = O, at which a stress, c_0+ c_(t), is applied the strain for time t > 0 can be

represented as follows.

E(t) = G O • J (t) + I;J (t-_) d°(_).d_ (6)

where J(t) is the compliance function of the material and dcr(_)/d_is the stress rate. A
similar equation can be used for the relaxation model to obtain the stress function introduced
by an arbitrary strain function 1

a(t) = eo.r(t) + Ior(t - de( )d (7)

where Y(t) is the relaxation function (Equation 4) and de (_)/d_is the strain rate. An

example of applying hereditary integrals for a multiple loading segment process is shown in

next section.

4
Hereditary integrals for a multiple loading process

Hereditary integrals with Prony series kernels can be applied to model a loading

process such as the one shown in Figure 2. The process in Figure 2 is divided into four

segments for which strain and strain rate functions are defined. The functions are:

to<t< fi to<t< fi
[ t/(t - to)
fi<t<t 2 de=
[ e /(q 0 - to) t 1 < t _< t 2
(8)
-81. t 3 - t2) t2 <t<t 3 ' dt ]- el/(t 3 - t2) t 2 < t _< t 3
e(t) = l elft
t3 <t <<-t4 L 0 t3 <t<t 4

where go = g (0) = 0 and to = 0.

For a material with a relaxation function in the form of a Prony series (Equation 4),

the stress functions of the loading process can be derived as follows:

Step 1. ( to < t_ h)

Substitute Equations 4 and first strain and strain rate functions of Equation 8 into Equation 7

and obtain:

_r(t) = e o. Y (t) + _oY (t - 4) d_) d4

t n -(t-_r _ E 1
=O+_oEo.(1- Y-J,Pi .(1-e )).--.d 4
i =1 tl

= )i " _-
E _-_pi_ 4-_-_pi'Ti'8 " 0

)i • t-- ___ [git 4- ___[9i " Ti -- _-_ [9i " Ti " e -/_ri

(9)

where n is the number of terms in the Prony series. To simplify the expression, n will not be

shown in following equations. Pi and "ciare the constants in the i-th term of the Prony series.

Step 2. ( tl < t __ t2)

Using the second strain rate function obtain:

 G([) : _0" Y ([) "[- I; 1-Y ([--_) ded_)d_ + it +y (t-_) de(_) d£d_
_

E o .e I ,
(10)
=0+ )i _--_l)i_ +_Pi_i.c +0
0

=E°el(tl-_pitl+_Pivi.e -_PiVi )
tl

Step 3. ( tz < t_ t3)

The third strain rate function yields:

0-([) : EO" Y ([)-1-I;1-Y ([- 4)d_)-d 4`'*f_ -[-If2-Y([-4)d-_)-d 4`'*f_

1
aq aq
de(_)
+f' Y(t-4) -d4
tl

=0+ Eo'el __Zpi_+Zpi.ci. e +0

/1 E -(t-_/_/
0

E--°--_-1 4_Zpi4+Zpi.ci. e
(t 3 -t 2)
t2

-t

- E° "el (tl-_.,pit 1 +_.,Pivi .e -_.,Pivi.e )

[1

-( t-t2_r i
(t-Y_Pit+Y_PiV i - t 2 + _., Pit2 - _., pivY i . e )
(t3-t 2)

(11)

Note, the first portion of Equation 11 is equal to Equation 10. Thus,

-( t-t2_v i
a(t) =%(t) (t-_pit+_pi'ci-t2+_pit2-_Pi'cie )
(t 3 -t 2)

(12)

where (Y2(t) is Equation 10 and it is function of time.

Step 4. ( t3 < t_ t4)

Similarly, the equation for the fourth step can be written as follows:

G(t) = O-2(t ) E---O--'_-J-] -(t-t3_r_

(t3 - t2) (t3-Y_Pit 3 + Y_PiTi .e -t 2 (13)

-( t-t2_r i
+Y_Pitz-Y_Pi'cie )

A numerical example of a multiple loading segment process using MATHCAD 12

software is shown in Appendix A. In the example, the stress function was calculated based

on the strain and strain rate functions shown above, and it employed a two-term Prony series.

The results of the viscoelastic analysis are shown in the stress-time and stress-strain plots.

This worksheet can be used to generate data in a parametric study involving viscoelastic

materials. The worksheet is also used as part of the weighted nonlinear regression algorithm

as is shown in the following sections.

Weighted Nonlinear Regression

The Prony series coefficients and retardation times appearing in Equation 4 need to

be determined in a regression analysis. Here, a standard nonlinear regression method (the

Marquardt-Levenberg Method 13'14) is used to perform the data fitting. In the nonlinear

regression, an error function (Z 2) with respect to the unknown constants is defined as,

x2(a) = .Yi - Y(xi;a). (14)

i=lL °-i

where xi and yi are the experimental data, function y(xi;a) is the model to be fitted, and (_i is

the standard deviation of measurement error of i-th data point. A set of unknown constants

(a) will be determined that minimize the error function 2. The error function (Z 2) is

approximated by its Taylor series with the quadratic form:

2"2 (a) _c-d.a+l/_2.a.D.a (15)

where c is a constant and d is the gradient of Z 2 with respect to the parameters a, which will

be zero where Z 2 is minimum. Matrix D is the second partial derivative of Z 2 with respect to

the parameters a. Initial trial values of a are specified and improved values are determined
7
 2
by the nonlinear regression algorithm. Iteration is continued until the error function, Z ,

effectively stops decreasing.

Since each Prony term includes two variables (pi and "ci)and since the instantaneous

modulus (Eo) must be determined, the total number of variables in the regression is 2n+l.

Based on thermodynamic principles, several constraint conditions must be applied:

P/>0, ZP/<1, r i >0, Eo>0 (16)

In addition, the distribution of the standard deviation of measurement error (c_) is not easily

determined based on the error of data acquisition equipment and the error of test machine, the

error is usually assumed to be uniform for all data points (c_i = 1). As is well known, the

viscoelastic effects are most significant at the beginning of the relaxation period, the fitting

error in this region is significant. Since the percentage of the number of data points at the

beginning of the relaxation period is less, the error function Z 2 is dominated by a long

uniform tail region of the relaxation period. To reduce the error and improve the fit at

beginning of the relaxation period, a weight function (w = 1/c_, 0 < c_ <1) is used. The larger

the weight factor at a data point, the better the curve fit the data point. There is no analytical

method to determine the weight function, thus a trial-and-error method is used. The

acceptance of the weight function is based on a graph of the data and regression model

results. Weighted nonlinear regression requires more iterations than unweighted nonlinear

regression, but it can provide a better fit to the experimental data in the region of most

interest.

Weighted Nonlinear Regression for Relaxation Test

A three-point bending relaxation test of a composite material was performed 15. The

specimen (12 in. x 2 in. x 0.768 in.) was loaded in 22 seconds to a maximum deflection of

0.103 in. at the middle of the span. Then the deflection was held for 11,711 seconds (Figure

3). The Marquardt-Levenberg nonlinear regression method was applied by the commercial

software SigmaPlot 13 to a hereditary integral model using two segments (loading and

holding) to obtain the Prony series coefficients for the data. Since no analytical method

exists to form the weight function, a trial-and-error method based on material properties was

used to obtain a fit curve. The viscoelastic stress decays exponentially in the relaxation test,

8
thereforetheviscoelasticeffectis moresignificantduringthe loadingperiodandatthe
beginning(< 30 seconds)of the holdingperiod.The numberof datapointsin theseperiods
(<100)is muchlessthanthe numberof datapointsin tail regionof the relaxationperiod(>
5000). Theerror function(Z2) will be dominatedby a long uniform tail regionof the

relaxationperiodif a uniform weightfunctionis applied. Therefore,a piecewiseweight

functionwasusedto obtainbetterfits for theseperiodsandimprovethe accuracyof the
regression. Figure4 showsthe loadrelaxationatthe beginningof theprocess. The dots
representthetestdata. Threeregressionresultsareshown. Thedash-dotcurveis the result
of a regressionanalysiswithoutthe weight function (w/o WF) for a two-termPronyseries.
Thelong-dashcurveis the resultfor a two-termPronyserieswith weight functionnumber1
(WF1) shownbelow:

150<t<22
(17)
w = J 10022 < t < 50
[150< t < 11711

The short-dash curve is the result for a three-term Prony series with weight function 2 (WF2)
as follows:

w = 106/(t2)22 _<t _<1000 (18)

10.t0< t < 22
11000 < t < 11711

Since an initial value is required for each of the variables, a trial data set based on the

test data was assumed. The sum of the P constants should be about 0.09 since the load at the

end of relaxation period is 9% lower than the load at the beginning. The retardation time

constants can be set to powers of ten. As long as the initial trial values are reasonable,

convergence will be achieved.

As is well known, the regression data fits better in a particular region if the relative

weight factor (weight factor / sum of weight factors) in that region is greater than the average

value. In Figure 4, the curves with weighted functions were closer to the data points near the

beginning of the relaxation period. Since the weight factor at t = 22 seconds of function 2

(=2066) is greater than the value of the weight function 1 (=100), the curve of function 2 is a

better match to the data than curve of function 1 at the beginning of relaxation.
 However,asshownin the Figure3,weight function2 doesnot fit the dataaswell as
the othercurvesafter1000seconds.It appearsthatfunction2 wasoverweightedatthe
beginningof the relaxationandthe relativeweight factoratthe otherregionwastoo small.
ThePronyconstantsof the regressionareshownin Table1.

Table 1. ThePronyconstantsfor theregression

Modulus P1 _1 P2 _2 P3 _3
(E0) (sec.) (sec.) (sec.)
w/oWF 7259.4 0.0259 134.69 0.0507 4898.64
WF1 7323.8 0.0319 57.23 0.0502 4040.36
WF2 7450.9 0.0262 9.29 0.0180 80.65 0.0460 2028.03

The results showed that a weight function should be selected based on the material

properties and test data distribution. Properly selecting the weight function in the most

significant region can improve accuracy of the regression.

Weighted Nonlinear Regression for a Multiple Loading Process

In a recent study 15, a thick composite panel responded viscoelasticly when it was

tested. In order to characterize the properties of the panel, a three-point bending test with

multiple loading segment processes was performed on it. Since the stiffness of the panel was

quite high, the test was conducted with a large hydraulic testing machine. High deformation

rates were not available with this loading machine. The loading head displacement schedule

of the machine is shown in Figure 5. The schedule is unlike a standard relaxation test.

Though load relaxation periods exist, the time required to apply a full load was clearly long

when compared to the relaxation periods. As mentioned in the previous sections, in order to

include the data of loading and unloading segments, a nonlinear regression combined with

the hereditary integral was used to simulate the data.

Since the deformation was very small compared to the specimen dimensions a linear

model was used for the calculations. The resulting displacements and the loads are linearly

related to the strains and stresses. Thus, the displacement and load data were used as strain

and stress data in the equations.

The displacement (strain) schedule in Figure 5 was divided into 11 steps. The

piecewise continuous function for the first 10 steps was assumed to be approximately linear,

similar to the previous example (Equation 8). The data for step 11 was fit to a cubic
10
polynomialfor a morepreciseregressionmodel. The strainandstrainratefunctionswere
definedasfollows,

8(t) = 81 tl < t _< t 2 (19)

t e r t/(q - to) to < t _ q

E i + a. t + b. t2 +... tlO < t _< t 11

J el/(q -to) to < t <_tl

de = 0 t1 < t <_t 2

dt [ ......
a+2.b.t+.., tlo <t<tll

The parameters of the polynomial used in step 11 were ei = -6.47997, a = 0.022976, b

= -2.516844e-5 and c = 8.772892e-9.

The piecewise continuous stress function for first 10 steps, based on the hereditary

integral, was derived by the same technique as shown from Equation 9 to Equation 13. The

load function for step 11 was derived as follows:

de(_)d£
= + f;,oY(t d_
//

=O-lo(0+ fttmo
Eo(1- Y_2pi(1-e-(t-_)/_'))(a + 2.b._ + 3. c._2).d_
i=1

=o'lo(t)+Eo.[(a. ¢ +b.¢ 2 +c.¢3)-_,pi .(a.¢ +b.¢ 2 +c.¢ 3)

+ _.Pi "e-(t-_)/ri "ri .(a + 2.b. (_-ri) + 3. c. (_ 2 -2._.r + 2 •ri2))] t

rio

= O'10(t ) + E 0 • [(1- ZPi)'( a" (t- tl0 ) + b. (t 2 - t?o ) + c. (t 3 - t3o))

+_.pi.ri .(a+2.b.(t-z'i)+3.c.(t 2 -2.t.r i +2.z'/z))

- _, Pi "e-(t-tm°)/ri "Vi "(a + 2. b. (tlo - v i) + 3. c. (tlo 2 - 2. tlo • r + 2./'i 2 ))

(20)

where _o(t) is the stress function from step 10.

A piecewise weight function was generated in order to increase the accuracy of

regression. Again, the weight function was determined by a trial-and-error method. After

11
severaliterations,the constantweight factorsusedfor eachstepweredetermineas:W(t) =
[1.0, 5.0, 1.0,1.5,1.0,7.0, 1.0,1.5,1.0,5.0, 1.0]andthe analysisresultswereshownin
Figure6. The weight factorswereequalto 1.0for theloading andunloadingsteps,andwere
greaterthan 1.0for theholding steps.
Thetotal numberof variablesin theProny seriesmaterialmodelis 2n+l, wheren is
thenumberof Pronyterms. In this study,a two-termPronyserieswassufficientto fit the
data. Theregressionanalysisresultsandtestdataareshownin Figures6 and7. The only
differencebetweenthe resultsof the two analyses,shownin the figures,wasthe
displacementfunctionof step11. Oneof the analysesassumedthatthe displacement
functionis linearandotherassumedit wascubic.Both modelresultsagreedcloselywith the
experimentaldata. Sincethe cubicpolynomialdescribeddisplacementin step11more
precisely,the onewith cubicdisplacementfunctionfit betterin step11thanthe onewith
linearfunction. The errorfunction(Z2 ) of the onewith cubicfunctionis 10%lower thanthe
onewith linear function.

CONCLUDING REMARKS:

A method of determining the coefficients in a Prony series representation of a

viscoelastic modulus from rate dependent data has been presented. The hereditary integral

method was employed to obtain an analytical representation of material response when it is

subjected to rate dependent loading. The analytical representation was used in a nonlinear

regression analysis, with measured data, to evaluate the Prony series constants. Several

regression analyses were performed using different weight functions. For the data analyzed

in this study, improved simulations of the hysteresis effects were obtained when the data at

the beginning of each relaxation period was appropriately weighted. Note, the data analyzed

here had loading and relaxation regions of similar length in time. Other weighting functions

may be needed for different loading schedules.

The method presented here provided a highly accurate representation of the material

behavior in the rate dependent loading region. It can also represent the response of a

viscoelastic material for other unique loading schedules. For example, it can be used for

schedules in which the material is not allowed to relax between subsequent loading changes.

12
References

1. F10gge, W. "Viscoelasticity", Blaisdell Publishing Co. , Massachusetts, 1975.

2. Cristensen, R. M. "Theory of Viscoelasticity" 2 na Edition, Academic Press, New
York, 1982.
3. Schapery, R. A., "Mechanics of Composite Materials", Vol. 2, Ed. Sendeckyj, G. P.,
Academic Press, New York, pp. 85-169 (1974)
4. Johnson, A. R., and Quigley, C. J., "A Viscohyperelastic Maxwell Model for Rubber
Viscoelasticity", Rubber Chemistry and Technology, Vol. 65, No. 1, pp 137-153
(1992).
5. Johnson, A. R., Quigley, C.J., Young D.G., and Danik,J.A. "Viscohyperelastic
Modeling of Rubber Vulcanizates", Tire Science and Technology, TSTCA,Vol. 21,
No. 3, July-Sept. 1993, pp. 179-199
6. Hill, S.A., "The Analytical Representation of Viscoelastic Material Properties Using
Optimization Techniques", NASA TM-108394, February, 1993
7. Bower, M.V. and Gant, D.F., "Stress Relaxation Functions: Method of
Approximation", NASA-CR-195830, April, 1994
8. "ABAQUS/Standard User's Manual", Hibbitt, Karlsson and Sorensen, Inc 1998
9. Johnson, A.R. "Modeling Viscoelastic Materials Using Internal Variables", The
Shock and Vibration Digest, Vol. 31,No. 2 1999, pp. 91-100
10. Lesieutre, G.A. and Govindswamy, K., "finite elements modeling of frequency-
dependent and temperature dependent dynamic behavior of viscoelastic materials in
simple shear," Int. J. Solids Structures 33(3), 1996,pp. 419-432
11. Johnson, A. R., Tessler, A., and Dambach, M. "Dynamics of thick viscoelastic
beams," ASME J. of Eng. Mat. And Tech., 119, 1997, pp. 273-278
12. "Mathcad User's Guide, Mathcad Plus 6.0", Mathsoft Inc., 1995.
13. "Transforms & Regressions, Reference Manual of Sigmaplot 4.0 for Windows",
SPSS Inc., 1997.
14. Press, W. , Plannery, B., Teukolsky, S. and Vetterling, W., "Numerical Recipes",
Cambridge University Press, 1989.
15. Chen, T.K., Dfivila, C. and Baker, D. "Analysis of Tile-Reinforced Composite Armor.
Part 2: Viscoelastic Response Modeling", Proc. 21st Army Science Conference,
Norfolk, VA. 1998

13
Appendix A Numerical Solution for Multiple Loading Segment Process

The purpose of this appendix is to present an input file for MATHCAD which can be used to
numerically compute results from Equations (9) to (13) for the multiple loading segment process.
Two Prony terms are used in this simulation. The loading process is defined by Figure 2. In the
data file below, text with a bold font represents a comment and with a normal font represents a
command. The result of the simulation is saved to a file (OUTPUT.PRN) which is shown at end
of the appendix.

Input file (Hereditary.MCD)

The relaxation function of the material is as follows:

Y(t) = E ( 1 - P1(1 - e -t/_l). P2 (1-e -t/_2))

Viscoelastic Material Constants:

E := 10 9 P1 := 0.2 "el := 10 P2 := 0.1 "c2 := 100

where E is the modulus, P and • are Prony constants

Loading time: (Seconds)

delt := 5 holdt := 50 tO := 0 tl := tO + delt t2 := tl + holdt
t3 := t2 + delt t4 := t3 + holdt

Piecewise Strain Function: (Multiple Loading Process) t := tO.. t4

el := 0.01 e2 := 0.0
strain function

t 0.01
E(t) := el' if O<t<tl
(tl tO)

(el) if tl<t<t2

t t2
el+(e2 el). if t2<t_t3
(t3 t2)
o i_
(e2) if t3<t<t4 0 55 110
0 otherwise

Hereditary integral for the stress function

Step 1 ( tO < t < tl) Loading to e = el
t := tO..tl

crl(t):=E'Cltl (t P1 t+P1 rl-PI.rl exp(-@l ) P2 t+P2 vZ-PZ.v2 exp(-_2))

Step 2 (tl < t < t2 ) Holding the load for 50 seconds

t := tl..t2

cr2(t) :=E. c1. (tl- P1. tl+ P1. vl. exp(--- -P1. vl. ) -P2. tl
vlt-tl) exp(- _1
tl

t - tl) _ . t))
+P2.v2.exp(--_ P2.v2 exp(- 75

14
 Step 3 (t2<t<t3) Unload to _-=e2
t := t2..t3

o-2(t) := E.(c2-c1).(t_Pl.t+Pl.rl_P2.t+P2.r2_t2+Pl.t2+P2.t2
t3 - t2

t - t2 t - t2))
- el. rl. exp(- --7i-- ) - e2. r2. exp(-

o-3(t) := o-2(t) + o-3_ l(t)

Step 4 ( t3 < t < t4) Holding for 50s

t := t3, t3 + 5..t4
t- t3 t
exp(- t3) _ t2
o-4_ l(t) := E. (t3(c2- -t2)c1). (t3 - P1. t3 + P1. rl. exp(- --7i-- ) - P2. t3 + P2. r2. r2

t
+Pl.t2+P2.t2-Pl.rl.exp( t2) _ P2. r2. exp(- t - t2))
rl r2
o-4(t) := o-2(t) + o-4 _ l(t)
Piecewise stress function stress function
t: t0.. t4 1.10 7

or(t) : (cYl(t)) if 0<t<tl

5.10 6
(cY2(t)) if tl<t<t2

(cY3(t)) if t2<t<t3 o(t)

(cY4(t)) if t3<t<t4 0

0 otherwise

5.10 6
Set up a matrix to print out of data o 50 100

xt,0: t Xt,l: _(t) xt,2: _(t)

Stress vs. Strain

Write the result file (OUTPUT.PRN)
1-10 7
WRITEPRN(output) := x

(End of Hereditary.MCD)

5-106
(Y(t)
Output file (OUTPUT.PRN)
Three columns of data (time, strain and stress) are
included in the output file as shown below:

0.005 0.01
(time) (strain) (stress)
£(t)

0 0 0
1 0.002 1979653.653
2 0.004 3921103.522
3 0.006 5827816.446
4 0.008 7702931.033
5 0.01 9549288.871
6 0.01 9389809.131
7 0.01 9244678.865

15
(a)

Creep Test
£

I • I

to Time to Time

(b)

Relaxation Test
E
G

Eo Go

U"..............
2;...........
I
• I

to Time to Time

Figure 1. Viscoelastic material characterization

tests: a) creep test, b) relaxation test.

16
 0.01

0.008

c
• 0.006
0.004

0.002

o
0 20 40 60 80 100 120

Time (seconds)

Figure 2. Multiple-segment loading process

760 -

75O

740

..d

730
720

710 -

700 - _,____
v"_

690

680 i i i I i i i I i i i I i i i I i i i I i i i I

2000 4000 6000 8000 10000 12000

Time, sec.

Figure 3. The regression results of relaxation test

17
 760 -

755

..d 75O

O
...d
745

740

735

730 i' i i i I i i i i I i i i i I i i i i I i i i i I i i i i I

20 25 30 35 40 45 50

Time, sec.

Figure 4. The regression results at beginning

of relaxation test.

18
 0.30
®
®
0.25

®
¢-
0.20

r-

E 0.15
o

o..
or)
n 0.10

0.05 Test Data

-- Cubic function for step 11
-- -- Linear function for step 11
0.00
I I I I I I

0 200 400 600 800 1000

Time (Sec.)

Figure 5. Test and simulation of loading schedule

19
 0
.d

E
0
z

Test Data
Nonlinear Regression with cubic polynomial strain for step 11
-- - Nonlinear Regression with linear strain function for step 11

I I I I

200 400 600 800 1000

Time (Sec.)

Figure 6. Normalized load vs. time for the test data and
nonlinear regression results.

20
 1.0

0.9

0.8

0.7
"o
o
.d 0.6
"o
U 0.5

E 0.4
o
z 0.3

0.2
Test Data

0.1 -- Cubic function for step 11

-- -- Linear function for step 11
0.0
I I I I I

0.00 0.05 0.10 0.15 0.20 0.25 0.30

Displacement,in.

Figure 7. Test and regression simulation of Normalized

load-displacement for CAV panel

21
 REPORT DOCUMENTATION PAGE Form Approved
OMB NO.0704-0188
Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data
sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other
aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and
Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188),
Washington, DC 20503.

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

May 2000 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

Determining a Prony Series for a Viscoelastic Material From Time Varying

Strain Data WU 706-13-31-01

6. AUTHOR(S)

Tzikang Chen

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

REPORT NUMBER
NASA Langley Research Center U.S. Army Research Laboratory
Hampton, VA 23681-2199 Vehicle Technology Directorate
NASA Langley Research Center L-17978
Hampton, VA 23681-2199

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING

AGENCYREPORTNUMBER
National Aeronautics and Space Administration
Washington, DC 20546-0001
and NASA/TM-2000-210123
ARL-TR-2206
U.S. Army Research Laboratory
Adelphi, MD 20783-1145
11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

Unclassified-Unlimited
Subject Category 64 Distribution: Nonstandard
Availability: NASA CASI (301) 621-0390
13. ABSTRACT (Maximum 200 words)

In this study a method of determining the coefficients in a Prony series representation of a viscoelastic modulus
from rate dependent data is presented. Load versus time test data for a sequence of different rate loading
segments is least-squares fitted to a Prony series hereditary integral model of the material tested. A nonlinear
least squares regression algorithm is employed. The measured data includes ramp loading, relaxation, and
unloading stress-strain data. The resulting Prony series which captures strain rate loading and unloading effects,
produces an excellent fit to the complex loading sequence.

14. SUBJECT TERMS 15. NUMBER OF PAGES

hereditary integral, viscoelasticity, weighted nonlinear regression, Prony 26

series, multiple loading segments 16. PRICE CODE

A03
17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION
OF REPORT OF THIS PAGE OF ABSTRACT OF ABSTRACT

Unclassified Unclassified Unclassified UL

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z-39-18
298-102