NASA 

Tech n i ca I 

Memorandum 

NASA TM-86579 

I 

(NASA- Ttl- E6 5 7 9) 

ON THE ANALYTICAL DETERMINATION OF RELAXATION 

MODULUS OF VISCOELASTIC MATERIALS BY PRONY'S 

INTERPOLATION METHOD 

By Pedro I. Rodriguez 

Structures and Propulsion Laboratory 

Science and Engineering Directorate 

December 1986 

CN 3fiE ANAIkIICAZ 

DETEREFNATIGN CE EELAXBTIC~ ~CEULUS OF 

IN2EREOLASlCU MEIHOD (LASA) ,I F CSCL ILA Unclas 

VlSCOELASTIC BAIIEEIALS BY EIiCYY '5 

G3/b4 43E55 

National Aeronautics and 

Space Administration 

George C. Marshall Space Flight Center 

MSFC - Form 3190 (Rev. May 1983) 
ACKNOWLEDGMENTS 

The author would like to express his appreciation to Mrs. Caroline Wimer for 

her enthusiasm and patience in typing, proofreading, and the retyping the manuscript. 

To her I am truly grateful. 

c 
TABLE OF CONTENTS 

INTRODUCTION .............................................................. 

DESCRIPTION OF PRONY 'S METHOD .......................................... 

COMPUTER IMPLEMENTATION ................................................ 

NUMERICAL EXAMPLES ....................................................... 

TP-H1148 Solid Propellant .............................................. V- 747- 75 Fluorocarbon 0-Ring Material .................................. 

DISCUSSION OF RESULTS .................................................... 

CONCLUSIONS ............................................................... 

APPENDIX A . CURVE FITTING PRONY'S EQUATION BY THE METHOD 

OF LEAST SQUARES ......................................... 

REFERENCES ................................................................. 

Page 

1 

2 

6 

7 

7 

12 

16 

19 

21 

25 

iii 
LIST OF ILLUSTRATIONS 

Figure 

1 . 

2a . 

2b . 

3 . 

4 . 

5 . 

6a . 

6b . 

7 . 

8 . 

Title 

TP-H 1148 Relaxation Modulus ......................................... 

TP-H 1148 Modulus versus actual time ................................. 

TP-H1148 B-Spline Modulus versus actual time ........................ 

TP-H 1148 Modulus (first 0.020 min) ................................... 

TP-H 1148 curve fit (Prony's method) ................................. 

V- 747- 75 Relaxation Modulus .......................................... 

V- 747- 75 Modulus versus actual time .................................. 

V- 747- 75 B-Spline Modulus versus actual time ......................... 

V-747-75 Modulus (first 0.080 min) ................................... 

V- 747- 75 curve fit (Prony's method) .................................. 

Page 

7 

8 

8 

9 

11 

12 

13 

13 

14 

17 
LIST OF TABLES 

Table 

1. 

2. 

3. 

4. 

5. 

A-1. 

Title Page 

FORTRAN Program PRONY ........................................... 6 

10 

V- 747- 75 Fluorocarbon Relaxation Modulus ............................ 15 

TP-H1148 Result Comparison Chart ................................... 18 

V- 747- 75 Result Comparison Chart .................................... 18 

FORTRAN Program CFIT ............................................. 24 

TP-H 1148 Propellant Relaxation Modulus.. ............................. 
TECHNICAL MEMORANDUM 

ON THE ANALYTICAL DETERMINATION OF RELAXATION MODULUS OF 

VISCOELASTIC MATERIALS BY PRONY'S INTERPOLATION METHOD 

INTRODUCTION 

When analyzing materials with time dependent mechanical properties, it is Of 

vital importance that the characterization of such properties be accurately defined. 

For viscoelastic materials, properties such as Relaxation Modulus, Creep Compliance, 

and time-dependent Poisson's Ratio are of utmost importance for numerical and closedform solutions to various problems. 

A popular method of obtaining analytical expressions for these properties consists of obtaining discrete values as a function of the logarithm of time and curve 

fitting the data to an appropriate expression. Due to the decaying nature of such 

properties as Relaxation Modulus and Poisson's Ratio, they are conveniently represented as a series of exponential functions. A widely used form of these functions is 

the so-called Prony series (due to Gaspard Francois Clair Marie Riche de Prony, 

1755- 1839), which can be expressed as, 

i= 1 

When using expressions such as this for curve fitting it is noted that there are 

too many unknowns for the amount of equations available for simultaneous solving. 

This generally leads to a trial-and-error approach where values of yi are assumed, 

and an expression for the time-dependent variable is obtained. 

repeated until a satisfactory curve fit has been obtained. 

This process is 

The intent of this paper is to demonstrate that with the use of the Intergraph 

Interactive Graphics Design System (1.G .D. S .) and its three-dimensional sculptured 

surfaces capability the true Prony method can be efficiently used, and the trial and 

error approach totally eliminated. 

to a third-degree, B-Spline interpolation which can be done quickly and with double 

precision accuracy on the I.G.D.S. 

desired number of equally spaced increments of time for the application of Prony's 

method. 

order to apply the method properly. 

results is not compromised by the introduction of the intermediate step of the B-spline 

interpolation. 

This is accomplished by curve fitting the test data 

The resultant curve can then be divided into the 

It must be stated at this time that equal time increments are necessary in 

It will also be shown that the accuracy of the 

The method is applied in this paper to a solid rocket propellant (TP-H1148) and 

a fluorocarbon elastomer commonly used as O-ring material (V- 747-75). 

The calculations are performed on a Univac 1108 computer and the graphics on 

a VAX 780 computer, both with double precisis: accuracy. 
DESCRIPTION OF PRONY'S METHOD 

This method is well documented [7,1,6] and will be explained here as a means 

of describing the Fortran program PRONY which will solve, without trial and error, 

for the exponents yi. 

Let Prony's equation be of the form 

m 

where k = 0,1,2,. . .n; i = 1,2,3,. . .m; and n L 2m - 1. 

value of Yi 

the successive values of time (abscissa) must form an arithmetic progression, that is 

tk = to + kw, where w is the difference between any two successive values of time. 

Equation (1) can now be expressed as, 

We first intend to obtain the 

and then the values of A and Bi. We must remember at this time that, 

n i= 1 

It is convenient to use the following a,,revia 

YiW 

V. = e 

1 

Yito C. = B. e 

1 1 9 

therefore, equation (1) becomes the set of equations 

m 

T (t)k =A+x C.VF 1 . i= 1 

ions to simp ify equation (2), 

( 3) 

(4) 

(5) 

Willers E71 explains how the process of obtaining the yi can be simplified by 

forming the difference between two successive ordinates in order to conveniently 

eliminate A from the approximating equations. As an example, take 4 equidistant 

measurements and 3 exponential terms or k = 0,1,2,3 and i = 1,2,3; this will lead to 

the following set of equations (5), 

2 
= A + c1v13 + c~v~~ + c3v3 3

f;t), 9 

the differences can be expressed as, 

- 

(tl0 A 1-0 = C1(V1-1) + C2(V2-1) + C3(V3-1) = P(t)l -T 

- 

= C,(V,-1)V12 + c2(v2-l)v22 + C3(V3-1)V,2 = F(t)3 - A 3-2 

- 

A 4-3 = cl(vl-1)v13 + c2(v2-1)v23 + c3(v3-1)v33 = F- (04 - %I3 

If we let the values of Vi be the roots of the equation 

m-2 + vm + s1 v m-l + s2 v ... + smLl v+s m =o 

then we can solve for the values of yi from equations (3) or 

1 = - log, vi . Yi w 

When trying to interpolate with three or more exponential 

tion (I)] there is the possibility of encountering complex roots Vi. 

precludes successful curve fitting. It does, however, change the final form of equation (1) thus forcing the analyst to utilize an expression with a combination of exponential and harmonic terms [ 1,7] . On the other hand, if the value of Vi is a negztive 

This by no means 

(9) 

terms [m 2 3 in equa3 
real value then equation (9) becomes unsolvable due to the fact that the natural 

logarithm of a negative number is undefined. 

impressive accuracy can be obtained, for the materials in question, by using 2 exponential terms in equation (1) , precluding the necessity for value of m greater than 2. 

Our problem now becomes that of obtaining the values of the S. in equation (8). 

Equations (7) can exist simultaneously if we write 

Fortunately, it has been found that 

1 

If we have a total of n+l data points and m exponential terms in the Prony series, we 

can obtain a total of n-m+l equations of the form (10) or 

- 

S + '2-1 S m-1 + ... + a =o '1-0 m (m+l)-m 

- - '3-2 s m +a 4-3 S m-1 + .'* + '(m+3)-(m+2) =o 

Equations (11) can readily be expressed in matrix form as 

r 

- - - ' (m+l)-m ... ' 1-0 '2-1 

- ' 2-1 

or simply 

161 {SI = 0 

'rn 

'm- 1 

1 

=o 

where the 6 matrix contains the ak-(k-l) coefficients and, in this case, k = 1,2 , 3.. .n. 

We can express (13a) in tensor notation as 

- 

6 sr=o 

qr 

4 
where 

q = 1,2,3 ... n-m+l 

r = 1,2,3 ... m+l 3 

- - - etc. From equations (13) we can - 

'3-2' in this manner we have Tl1 = A l-o, 6 22 - 

form a set of normal equations [ 1,7] by the method of least squares which written in 

matrix form are 

In the tensor notation we have 

ajR s = E. a1 

where 

j,R = 1,2,3, ... m 

and 

n-m+l 

q= 1 

n-m+l 

q= 1 

where p = m+l. 

Once the and E. are obtained, equation (14a) can be solved for {SI by 'j 1 

premultiplying both sides of the equation by the inverse of [a] or 

Equation (8) can now be solved for the Vi and from equation (9) the exponents yi can 

be obtained. 

5 
COMPUTER IMPLEMENTATION 

Table 1 is the listing of the FORTRAN program PRONY 

method. The program reads the xi and yi, in this case time 

and forms the a matrix. It is then inverted by means of the 

method and the values of S are obtained as in equation (17). 

Once the Si are obtained, then the yi can be obtained. 

implementing the Prony 

and Relaxation Modulus, 

Gauss-Jordan reduction 

It is then a simple 

more to solve for the constants A and Bi by means of the least-squares method. 

TABLE 1. FORTRAN PROGRAM PRONY 

PEOROBIN197*ME799(1).PRONV 

1 

2 

3 

4 

5 

6 

7 

B 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

PARAMETER N.32, M=2, NM1 =N- 1, MP 1 =M+ 1. M2=2+M 

PARAMETER NM1M=NMl-M.M2PZ12*M+2,M4P4=4*M+4 

DIMENSION JC(12) 

DOUBLE PRECISION ROOT(M2P2).AA(M2P2),6B(MZPZ) 

DOUBLE PRECISION CC(M4P4).W.GAM.GAMMA 

DOUBLE PRECISION V(P).S(MPl) 

DOUBLE PRECISION E(N).T(N).OIFF(NMl).OELTA(NY1M.MP1) 

DOUBLE PRECISION ALPHA(M.YPl).EPS(M) 

00 10 K-l.N 

REA0(5,100)T(K).E(K) 

WRITE( 10.100)T(K),E(K) 

100 FORMATO 

10 CONTINUE 

00 11 K-l.NM1 

KP1=K+1 

OIFF(K)-E(KPl)-E(K) 

00 12 1-1.NM1M 

00 12 J-l.MP1 

JJ=I+J- 1 . 

11 CONTINUE 

OELTA(I.J)=OIFF(JJ) 

12 CONTINUE 

00 14 I-l.M 

00 14 J'1.M 

ALPHA( I, J)*O. 0 

DO 14 K=l.NMIM 

ALPHA(I.J)-ALPHA(I,J)+OELTA(K.I)+DELTA(K.J) 

00 15 I'1.M 

EPS(1)-0.0 

00 15 J=l.NHlM 

EPS(I)=EPS(I)+OELTA(~,I~*OELTA~d,MPl) 

00 16 I.1.M ALPHA(I.MPl)=-l.O*EPS(I) 16 CONTINUE 

V(0'4.0 

CALL OGJR(ALPHA.MP1.M.M.MP1.$300.JC.V) 

14 CONTINUE 

15 CONTINUE 

300 WRITE(6,301)JC(O 

301 FORMAT(' SYSTEM SOLVE0 UP TO EOUATION NO. '.13//) 

S( 1)=1.0 

00 17 I-l.M 

J=M-I+Z 

S(J)=ALPHA(I.MPl) 

00 18 I-1.MP1 

WRITE(6.400) I.S(I) 

17 CONTINUE 

400 FORMAT(' S('.I2.')*'.E20.6) 

18 CONTINUE 

TOL-1.0-16 

CALL ROOTZ(S.M.ROOT.AA.BB.CC.TOL) 

WRITE~6.6l~~ROOT~I~.I-l.M2~ 61 FORMAT(20X.'ROOT'.5X.2E30.8) 

W.T(3)-1(2) 

00 19 I.1.M 

NN- 2 * I - 1 

GAM=ROOT(NN) 

GAMMA.(OLOG(GAM) )/W 

WRITE(6.700)GAMMA 

700 FORMAT(2OX.'EXPONENT'.E30.8) 

19 CONTINUE 

STOP 

EN0 

eBRKPT PRINTS 

6 

1 

step 
NUMERICAL EXAMPLES 

TP-H 1148 SRM Propellant 

Figure 1 shows a plot of the Relaxation Modulus versus log of time as obtained 

from a stress relaxation test [21 . - 1.659583-3 to 1.659583+1 min. In order to properly apply this method, this data must 

be plotted in real time and then divided into equal time increments, w. This would be 

close to an impossible task unless performed on a graphics computer with the capabilities of I.G.D.S. 

As can be seen, the time range extends from 

-1.78 -1.28 -0.78 -0.28 0.22 0.72 122 -2.78 -2.28 

LOO OF TIME IN MINUTES 

Figure 1. TP-H 1148 Relaxation Modulus. 

As can be seen in Figure 2a, we have plotted seven of the nine data points 

presented in Figure 1. 

curve, even though all nine points are used when performing the B-Spline interpolation shown in Figure 2b. As can be seen, the greatest relaxation occurs during the 

first 0.020 min of the test, and therefore we shall concentrate on this portion of the 

curve. 

This is done in order to show with clarity the shape of the 

By llzoorning-inl' to the portion of interest, we can take increments of time of 

0.001 min. 

The I.G.D.S. can accurately locate the intersection between the B-Spline curve and 

These values can be read directly off the screen and tabulated as 1.n Table 2. This 

is the basic information required in order to run PRONY. 

Figure 3 shows the portion between 0.002 and 0.020 min of Figure 2b. 

- the time increments to yield values of Relaxation Modulus, E(t) , at each intersection. 
-lI 

I, 

v) 

W 

I3 

$ 

Q, 

E 

-r--- 

8. 

5 0 

N 

8 
20 

15 

E (t) 

KSI 

10 

5 

0 .005 .010 .015 .020 

TIME IN MINUTES 

Figure 3. TP-H1148 modulus (first 0.020 min). 

Equation (1) shows the general form of the Prony series used in our discussion. 

For this example let us take m = 2 since at m 2 3 we have encountered complex values 

of Vi. By inputing the data from Table 2 into our program we obtain the following 

form of equation (8) : 

V2 - 0.822591 V + 0.095196 = 0 

where, 

S1 = -0.822591 

S2 = 0.095196 

The solution of equation (18) yields the following roots 

9 
TABLE 2. TP-H 1148 PROPELLANT RELAXATION MODULUS 

TIME 

(MINUTES) 

.002 

.003 

.ow 

.005 

.006 

.007 

.008 

.009 

.010 

.011 

.012 

.013 

.014 

.015 

.016 

.017 

.018 

.019 

.020 

V1 = 0.6832668 

V2 = 0.1392342 

20.3704 

18.1835 

16.0257 

14.4590 

13.3989 

12.6508 

12.1197 

11.7508 

1 1 5090 

11.3706 

11.2952 

11.2244 

11.1563 

11.0909 

1 1.0282 

10.9683 

10.9112 

10.8569 

10.8056 

Substitution of equations (19) into equations (9) yield the values of the required 

exponential term constants. 

- loge 0.6832668 = -380.8698 y1- o.001 

loge 0.1393242 = -1971.8239 - 1 

y2 - m 

We can now rewrite our Prony equation as 

-380.8698 t + B2 e -1971.8239 t E(t) = A + B1 e 

10 
where t = time in minutes, and A, B1, and B2 are constants to be determined by the 

method of least-squares [ 51. These constants are, 

A = 10973.3 

B1 = 23397.3 

B2 = -79056.7 

It should be noted that equation (21) is valid only for the time range 0.002 to 0.020 

min and that extrapolation beyond this range can lead to erroneous results. 

shows a plot of equation (21) as compared to the B-splined experimental data. 

can be seen extremely good accuracy is obtained with only 2 terms. 

Figure 4 

As 

E (t) 

KSI 

20 

15 

10 

5 

0 

--- 

I 

t 

I 

I 

I 

I 

I 

I 

I 

I 

---- I 1-- 

RELAXATION MODULUS FROM B-SPLINED TEST DATA : 

I 

I , --- LEAST SQUARES CURVE FIT *..I 

.006 .010 .015 .020 

TIME IN MINUTES 

Figure 4. TP-H 1148 curve fit (Prony's method). 

11 
V- 747- 75 Fluorocarbon 0-Ring Material 

Figure 5 shows a plot of the Relaxation Modulus versus log of time as obtained 

from a stress relaxation test [ 41. Once again the time range extends from 1.659583-3 

to 1.659583+1 min. 

followed in order to obtain values of E(t) as a function of equal time increments w. 

Figures 6a and 6b are plots of actual data points and B-Spline interpolation, respectively. 

The same procedure as for the TP-H1148 solid propellant is 

5000 

4000 

3000 

2000 

1000 

0 

IIIII 

V-747-75 FLUOROCARBON O-RING 

MATERIAL AS MEASURED TEST DATA 

AT 3OoF AND 2% STRAIN 

-2.78 -2.18 -1.78 -1.38 -1.08 -0.90 -0.78 0.22 1.22 

LOG OF TIME IN MINUTES 

Figure 5. V- 747- 75 Relaxation Modulus. 

For this example we shall take equal increments of time of 0.0025 min and 

analyze the portion of the curve that extends to 0.0800 min. 

equal time increments, and Table 3 shows their respective values of Relaxation Modulus 

as obtained from the I .G .D .S. 

Figure 7 shows the 

With this information we are again ready to run PRONY. 

Again we have found that for m 1 3 we encounter complex values of Vi; therefore, for m = 2 we obtain the following form of equation (8). 

2 V - 1.14458 V + 0.255224 = 0 

with 

S1 = -1.14458 

S2 = 0.255224 . 

12 
I 2 4 I 

0 0 I----- 

, I- L d --;-----. I----- 

i 

0 

0 

a 

I 0 

, 

1 

II 

--1------ 8 r---- 

I-___ 

, 

8 a' 8 c.i cu' 88 c 8 

c 

88 

ei cu' 

8 

d 

A 

c- 

-v) 

WY 

- 

c- 

-v) 

WY 

13 
I t; r----------- 

w 

, 

I 

# 

r---------- 

c 

z 

Q 

2 

4 

X 

w 

K 

w 

w 

w 

a 

5 

B 

K 

lA 

c" 

2 

z 

a 

n 

I, 

8 

u 

8 

Pi 

8 

cu' 

n 

c 

E 

.#-I 

0 

00 

0 

0 

+ 

m 

k 

c 

W 

; 

x 

4 

3 

z 

lr 

b 

b 

b 

I 

I 

3 

CI 

al 

k 

'7 

ho 

tz 

14 
TABLE 3. V-747-75 FLUOROCARBON O-RING MATERIAL RELAXATION 

MODULUS AS OBTAINED FROM FIGURE 7, FOR TIME BETWEEN 

0.0025 AND 0.0800 MINUTES 

TIME 

(MINUTES) 

.0025 

.0050 

.0075 

.0100 

.0125 

.0150 

,0175 

.0200 

.0225 

.0250 

.0275 

.0300 

.0325 

.0350 

.0375 

.04W 

.0425 

.0450 

.0475 

.woo 

.0525 

.%50 

.m75 

.moo 

.M25 

.0650 

.0675 

.0700 

.0725 

.07 50 

.0775 

.0800 

3977.9 

3317.1 

3032.1 

2873.1 

2772.4 

2687.4 

2613.8 

2551.3 

2499.9 

2459.7 

2430.5 

241 1.6 

2395.4 

2379.8 

2364.8 

2350.4 

2336.6 

2323.6 

2310.7 

2298.7 

2287.2 

2276.4 

2266.1 

2256.4 

2247.3 

2238.7 

2230.6 

2222.8 

2215.5 

2208.1 

2202.1 

2196.1 

15 
The roots of equation (23) are 

V1 = 0.84113 

V2 = 0.30345 

The required exponential constants are obtained from equation (9) as 

1 1 

loge 0.84113 = -69.204 (25a) I I 

0.0025 y1 = w loge v1 = 

- log, 0.30345 = -477.017 - 1 

y2 - w loge '2 - 0.0025 

Our Prony equation for the V-747-75 material during the range 0.0025 5 t L 0.0800 

is now 

- 477.017t 

-69-204t + B2 e E(t) = A + B1 e 

The constants A, B 1, and B2 can be determined by the method of least squares [ 51, 

A = 2239.88 (27a) 

B1 = 1247.08 (27b) 

B2 = 2253.44 (27c) 

Figure 8 shows a plot of equation (26) as compared to the B-Splined experimental 

data. Again a good correlation between the data and the curve fit is seen. 

DISCUSSION OF RESULTS 

Tables 4 and 5 show the comparison between the test data and the results from 

PRONY. We can see in both cases that there is very good accuracy between test 

points and curve fit points as long as we maintain ourselves between the time range 

under study or as determined by Tables 2 and 3. In both cases the value of E(t) 

departs from the actual values as soon as we leave the time range under study. At 

the beginning of the time range this problem can be overcome by extrapolating the 

the test data on to t = 0 [or t = 0.00001 min to avoid numerical problems when plotting 

Log(t)]. At the end of the time range the accuracy can be maintained by including 

as many time points as required for a specific analysis. 

16 
, 

l a 

a 

n 

r----------- II - - - - - - - - - - 

.__________ 

, 

I---------- 

I---------- 

c, 

w 

In 

w 

z 

J 

Q 

v) 

I 

m 

L 

LL 

v) 

2 

3 

0 

L 

z 

0 

i= 

4 

X 

n 

5 

w 

K 

I 

v) 

W 

c 

3 

z 

f 

E 

- 

w 

I8 

d 

17 
TABLE 4. 

Time (min) 

0.00166 

0.00661 

0.01666 

0.04168 

TABLE 5. 

Time (min) 

0.00166 

0.00611 

0.01666 

0.04168 

0.08318 

TP-H1148 RESULT COMPARISON CHART 

Relaxation Modulus (ksi) (0.002 1. t 5 0.020) 

Test Data B -Spline PRONY 

21.114 21.114 20.412 

12.913 12.913 12.861 

10.992 10.992 11.014 

10.144 10.144 10.973 

V-747-75 RESULT COMPARISON CHART 

Relaxation Modulus (psi) 

(0.0025 I t I 0.0800) 

Test Data B - Spline PRONY 

4247.00 4247.00 2636.05 

3113.00 3113.00 3179.15 

2639.00 2639.00 2634.39 

2341.00 2341.00 2309.58 

2189.00 2189.00 2243.82 

Another point of interest is that if any of the roots of equation (8) are negative 

we will not be able to solve equation (9) for that particular yi since the natural 

logarithm of a negative number is undefined. 

not been addressed by any available source on the implementation of Prony's method 

[ 1,6,7] . require that the values of Vi meet the following requirements : 

A method to overcome this problem has 

It is the author's belief that due to the nature of equation (9) we must 

1. 

2. 

All Vi shall be positive real values or complex conjugates. 

If the data to be curve fitted is of a decaying nature as a function of time, 

the Vi must be less than 1.00 except for complex conjugate values. 

the Vi must be greater than zero, with at least one greater than 1, except for complex conjugate values. 

3. If the data has increasing values of the or'dinate with increasing time, then 

18 
CONCLUSIONS 

A very good tool for the analytical determination of functions which can be 

represented by a series of exponential terms is available. The method, when utilized 

properly, can provide the analyst accurate expressions of the ordinate in question. 

Although the method leaves the unanswered question of what to do when Vi is a real 

negative number, the author feels that enough accuracy for materials such as elastomers and solid propellants can be obtained to perform an accurate stress and strain 

analysis of these materials under load. 

. 

19 
APPENDIX A 

CURVE FITTING PRONY'S EQUATION BY THE METHOD OF LEAST SQUARES 

After obtaining the values of the yi exponents, Prony's equation (1) can be 

solved for the A and Bi constants from the data points by using the method of least 

squares [ 5,7] . 

The method, due to Legendre, is applied to Prony's equation as follows. Let 

Prony's equation be 

n Yit 

F(t) = A + B~ e 

i= 1 

If we want to approximate a curve with f(t) as its discrete data points, we can 

express a deviation for the argument value of t as, 

where m = 1,2,3,. . . , k and k is the total number of data points. 

squares of all the deviations is expressed as, 

The sum of the 

k k 

D = (f(t)m - f(t)m)2 = 6m 2 

m=l m=l 

If we substitute equation (A-1) into equation (A-3), we obtain, 

The 

this 

k 

D=C (A+ 

m= 1 

method of least 

to be true it is 

i= 1 

(A- 3) 

(A-4) 

squares requires that equation (A-4) be a minimum. 

necessary that 

In order for 

From equations (A-5) and (A-4) we can write, 

(A-5) 

21 
m= 1 i= 1 

k n 

aF(t) 

a B, - f(t),) - = 0 Yitm laD - (A+C ~~e2 aB, J ’ m=l i= 1 

where j = 1, 2, 3...n. 

(A-7) 

Equation (A-7) constitutes a set of n equations or one equation for each exponential term. These latter two equations can be rewritten as 

(A-8) 

m=l i= 1 m= 1 m=l 

For the purpose of writing a computer program that will solve for the constants A and 

Bj , equations (A-8) and (A-9) may be expressed in matrix form as follows, 

[Cl {A) = {FI (A-10) 

where the vector {A) holds the sought unknowns A and Bj and the coefficients of IC] 

and {F) are defined as follows, 

aT(tlm 

=k cll = c aA 

m=l 

k Yi-1 t m af(tIm 

aA cli = c e 

m=l 

(A- 11) 

(A- 12) 

(A- 13) 

1 m=l 
JL m=l 

For equation (A-14) only we have 

i = j = 2, 3, ... k+l 

whereas for equations (A-12) and (A-13), 

i = 1, 2, 3, ... k 

The coefficients of the vector {F} are defined as, 

m=l 

(A- 14) 

(A- 15a) 

(A- 15b) 

(A- 16) 

(A- 17) 

and in this case the i is defined as in equation (A-15a). 

Once the C and F coefficients are properly identified, both sides of equation 

(A-10) can be premultiplied by the inverse of [C] to obtain the solution vector {A} 

or , 

CAI = LC1-l {F) (A- 18) 

Table A-1 shows a printout of the program CFIT applying the method described above. 

23 
TABLE A-1. FORTRAN PROGRAM CFIT 

PEOROBIN197*ME799(1).CFIT 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

71 

78 

19 

80 

28 

PARAMETER K=2.KPl=K+l.KP2=K+Z 

PARAMETER M=33,VAL=lOO..MM=100 

DIMENSION JC(lP).E(M).T(M) 

DOUBLE PRECISION B(K).F(KPl).C(KPl.KP2).PRONV(MM) 

DOUBLE PRECISION V(Z),MOD(MM) 

DOUBLE PRECISION TT(MM+l).R(K) 

C INITIALIZE MATRIX AN0 VECTOR 

DO 10 N=l,KPl 

F(N)=O.O 

DO 10 L=l.KPZ 

C(N.L)=O.O 

R(1)=118.159 

R(2 )=661.228 

10 CONTINUE 

C R( 3 1.4 15.608 

C R ( 4 1 =500. 

00 11 I=l,M 

READ(S.lO)T( I 1, E( I 

100 FORMATO 

11 CONTINUE 

C( 1.1 )=l .OtM 

00 20 N=Z.KPl 

DO 20 I=l,M 

C(l.N)=C(l.N)+EXP(-R(N-l)*T(I)) 

C(N,l)=C(N.l)+EXP(-R(N-l)*T(I)) 

20 CONTINUE 

DO 30 J=Z.KPl 

DO 30 N52.KPl 

DO 30 I=l,M 

C(J.N)=C(J.N)+EXP(-R(N00 40 IZ1.M 

F( l)=F( l)+E(I) 

DO 50 N.2.KPl 

DO 50 I-1.M 

F(N)=F(N)+E(I)*EXP(-R(N 

DO 60 N=l.KPl 

C(N.KPZ)=F(N) 

60 CONTINUE 

VI 1 )=4.0 

30 CONTINUE 

40 CONTINUE 

50 CONTINUE 

CALL OG&(C .KP2 . K P 1 . K P 1 .KP2,5300. JC .V) 

300 WRITE(6.301)JC(l) 

301 FORMAT(' SYSTEM SOLVE0 UP TO EOUATION NO.'.I3//) 

A=C(l.KPZ) 

00 70 N=Z.KPl 

B(N-l)=C(N.KPZ) 

WRITE(6.102)A.B(l),B(2),B(3) 

TMAX=T (M) 

TMIN=T( 1) 

DO 80 I=l.MM 

PRONV ( I ) -0.0 

70 CONTINUE 

102 FORMAT(lOX.4E15.6) 

80 CONTINUE 

XINC=(TMAX-TMIN)/VAL 

TT( 1 )=TMIN 

00 90 I=l.MM 

00 101 N=Z.KPl 

PRONV(I)=PRONV(I)+B(N-~)*EXP(-R(N-~)~TT(I)) 

MOO(I)=A+PRONV(I) 

101 CONTINUE 

TT( I+l)=TT( I )+XINC 

90 CONTINUE 

ARG=TMIN 

v1=0.0 

v2=0.0 

ARG 1-0.0 

00 110 I=!.MM 

WRITE(9)ARGl.Vl.VZ 

ARG~=ALoG~O(ARG) 

CALL XNTRP5(ARG.T.E.M.2.Vl.N) 

V2=M00( I) 

ARG=ARG+XINC 

C WRITE(6.105)ARG.Vl.VP 

C 105 FORMAT(lOX.3E20.6) 

110 CONTINUE 

STOP 

EN0 

, 

24 
REFERENCES 

1. 

2. 

3. 

4. 

5. 

6. 

7. 

Buckingham, R. A:: Numerical Methods. 

London, 1957. 

Sir Isaac Pitman and Sons, LTD., 

Harvey, A. R . : 

TP-H 1148 SRM Propellant. TVJR- 14153, Morton Thiokol, Inc. /Wasatch Division, 

1984. 

Strain Evaluation Study and Propellant Characterization for 

Ketter, R. L. and Prawel, S. P. : 

McGraw-Hill Book Co., New York, 1969. 

Modern Methods of Engineering Computation. 

Russell, P. G. and Harvey, A. R.: Space Shuttle Characterization of Fluorocarbon 0-Ring Material. TWR- 300165, Morton Thiokol, Inc. /Wasatch Division, 

1986. 

Sokolnikoff, I. S. and Sokolnikoff, E. S. : 

and Scientists. McGraw-Hill Book Co., New York, 1941. 

Higher Mathematics for Engineers 

Whittaker, Sir Edmund and Robinson, G. : The Calculus of Observations. 

Blackie and SM, LTD., London, 1958. 

Willers, F. A., translated by Beyer, R. T. : Practical Analysis. Dover Publications, Inc., New York, 1947. 

25 
APPROVAL 

ON THE ANALYTICAL DETERMINATION OF RELAXATION MODULUS OF 

VISCOELASTIC MATERIALS BY PRONY'S INTERPOLATION METHOD 

By Pedro I. Rodriguez 

The information in this report has been reviewed for technical content. 

of any information concerning Department of Defense or nuclear energy activities or 

programs has been made by the MSFC Security Classification Officer. This report , in its entirety, has been determined to be unclassified. 

Review 

G. F. McDONOUGH 

Director, Structures and Dynamics 

Laboratory 

26 *U.S. GOVERNMENT PRINTING OFFICE 1987--730-067/40059 

~- ~~ 

~~ 
. 

. 

TECHNICAL REPORT STANDARD TITLE PAGE 

I REPORT NO. 2. GOVERNMNT ACCESSION NO. 3. RECIPIENT'S CATALOG NO. 

NASA TM - 86579 

On the Analvtical Determination of Relaxation Modulus of 

I. TITLE AND SUBTITLE 5. REPORT DATE 

December 1986 

Viscoelastic Materials by Prony's Interpolation Method 6. PERFORMING ORGANIZATION CODE 

1. AUTHOR(S) . 8. PERFORMING ORGANIZATION REPOR r # 

Pedro I. Rodriguez 

b. PERFORMING ORGANIZATION NAME AND ADDRESS 10. WORK UNIT NO. 

George C. Marshall Space Flight Center 11. CONTRACT OR GRANT NO. 

13. TYPE OF REPOR;' 8 PERIOD COVERED 

Marshall Space Flight Center, Alabama 35812 

- 

2. SPONSORING AGENCY NAME AND ADDRESS 

National Aeronautics and Space Administration 

Washington D .C . 20546 

Technical Memorandum 

I 

5. SUPPLEMENTARY NOTES 

Prepared by Structures and Propulsion Laboratory Science and Engineering Directorate. 

6. ABSTRACT 

A computer implementation to Prony's curve fitting by exponential functions is 

presented. The method, although more than one hundred years old, has not been 

utilized to its fullest capabilities due to the restriction that the time range must be 

given in equal increments in order to obtain the best curve fit for a given set of 

data. 

The procedure used in this paper utilizes the 3-dimensional capabilities of the 

Interactive Graphics Design System (I.G.D.S.) in order to obtain the equal time 

increments. The resultant information is then input into a computer program that 

solves directly for the exponential constants that yield the best curve fit. 

exponential constants are known, a simple least squares solution can be applied to 

obtain the final form of the equation. 

Once the 

7. KEY WORDS 18. DISTRIBUTION STATEMENT 

Viscoelastic Prony Interactive Graphics 

Design Unclassified - Unlimited 

I 

19. SECURITY CLASSIF. (d thh roeat) 120. SECURITY CLASSIF. (of thh pago) 121. NO. OF PAGES 122. PRICE 

NTIS I Unclassified Unclassified 32