Determination of Mechanical Spectra from Experimental Responses Določanje mehanskega spektra na osnovi izmerjenih relaksacijskih krivulj Emri L, Fakulteta za strojništvo, Ljubljana N. W. Tschoegl, California Institute of Technology, Pasadena, California 91125, USA A recursive computer algorithm was developed which generates line spectra from experimental response functions. The method allows storing information on the mechanical properties of polymeric materials in a convenient way. The algorithm also interconverts betvveen relaxation and retardation spectra. From the spectra, any desired response funetion can then be recovered. The algorithm essentially utilizes the fact that the kernel functions resemble step functions. Slightly different codes are used for each kernel funetion. The appearance of negative relaxation or retardation lines is obviated. Mathematically such lines vvould be acceptable, and they do not seriously affect reconstruction of responses vvithin relaxation or retardation behavior. Hovvever, they vvould seriously interfere vvith interconversion betvveen the tvvo types of behavior, and they vvould also pose problems in the interpretation of the spectra. Key vvords: viscoelasticity, relaxation, retardation, spectra, kernel funetion, polymers Razvit je rekurziven računalniški algoritem, ki izračuna spekter iz eksperimentalno dobljenih odzivnih funkcij. Metoda omogoča shranjevanje informacij o mehanskih lastnostih polimernih materialov v prikladni obliki. Algoritem omogoča tudi interkonverzijo med relaksacijskim in retardacijskim spektrom. Iz spektra lahko izračunamo katerokoli materialno funkcijo. Bistvo algoritma je uporaba dejstva, da so jedra podobna koračni funkciji. Malenkostno drugačne metode so uporabljene za vsako od jeder. Pri določevanju relaksacijskega in retardacijskega linijskega spektra, se s tem algoritmom izognemo negativnim spektralnim linijam. Matematično so take linije sprejemljive in nimajo velikega vpliva na rekonstrukcijo materialnih funkcij. Problem zaradi negativnih spektralnih linij se pojavi pri interkonverziji med dvema tipoma obnašanja, poleg tega pa povzročijo nepremostljive ovire pri fizikalni interpretaciji spektra. Ključne besede: viskoelastičnost, relaksacija, retardacija, spekter, jedro, polimeri 1. Introduction In the response of a linearly viscoelastic material to a strain excitation, eomplete information on the time-dependent part of the response is contained in the relaxation spectrum. In the response to a stress excitation the same role is played by the retardation spectrum. The response then consists of the appropriate viscoelastic constants (such as the equilibrium modulus, or the glass compliance and the steady-flow fluidity), in addition to the integral over the spectrum multiplied by a kernel funetion char-acteristic of the type of excitation chosen to elicit the response. If this is a strain or a stress as a step funetion of time, the result is the relaxation modulus in the first. and the creep compliance in the second čase. Thus, once the spectrum is known in addition to the viscoelastic constants, it is possible to generate from it the response to any desired type of excitation. The spectrum itself is not accessible by direct experiment. From a theoretical response curve it can often be calculated1". From experimental data the spectrum is necessarily obtained as an approximation to the true spectrum. A variety of methods have been proposed1"2 to extract approximations to the spectra from experimental data by mathematical manipulation. The well knovvn methods based on numerical or graphical differentiation typify this approach. In another approach an attempt is made to determine a distribution of diserete spectral lines from vvhich the original response curve can be more or less faithfully repro-duced1'. Among the better knovvn older methods are Procedure X of Tobolsky and Murakami3, the collocation method of Schapery4, and an extension of it by Cost and Becker1 vvhich they call the multidata method. The tvvo last named require matrix in-version. Both are likely to generate a negative spectrum lines, vvhich makes virtually impossible to use this spectra for inter-conversion betvveen the relaxation and retardation behavior. In the past few years several new papers have been devoted to this subject6"14. An extensive comparison of these methods along vvith the method introduced here is presented elsevvhere". We present here the method based on an iterative computer algorithm for calculating a distribution of spectral lines vvhich, for a given set of data, is unique vvithin the desired degree of ac-curacy. It must be emphasized, that the discrete distribution of moduli or compliances on respondance (relaxation or retardation) times obtained in this way is stili an approximation. Hovvever, as vvill be shovvn belovv, the calculated distributions appear to yield better results than any of the other methods. In this paper we shall sketch only the application of the algorithm to the relaxation modulus. The power of the algorithm vvill be than demonstrated on the experimental data obtained by Catsiff and Tobolsky16. These data vvere originally reported in terms of the tensile modulus, E(t). For the comparison reasons, presented in' the data vvere first converted to the shear modulus, as reported in1'. For the same reason the units of the modulus vvere converted from dynes/cm2 and hours to N/m2 and seconds. A complete presentation of the algorithm has been published elsewhere-in the series of four papers17 20. The first paper in this series17 describes the algorithm for obtaining a discrete distribution of relation times from simulated relaxation modulus data, or of retardation times from simulated creep compliance data. The second paper"1 deals vvith the determination of the spectra from theoretical storage and loss functions. The third paper1'' takes up the problem of converting betvveen relaxation modulus and creep compliance. In these first three papers the algorithm have been thus applied to data sets obtained by sampling continuous theoretical curves. This has simplified presentation of the details and the povver of the algorithm. The fourth paper20 deals finally vvith the application of our method to experimental data, i.e., data that are not free of experimentai error. 2. Theoretical Essentially the method consists in predetermining a set of respondance times vvhich are equally spaced on the logarithmic time-axis, and then calculating the strength of the spectral line associated vvith each respondance time. Our method shares this feature vvith the collocation method and the multidata method. Hovvever. in contrast to these methods, ours does not require matrix inversion and thus avoids mathematical difficulties associated vvith sueh an inversion (sueh as. e.g., occasionally gen-erating troublesome negative respondance times). We calculate the intensity (i.e. the strength) of the k"' spectral line eorresponding to the k,h respondance time. Tk. from ali source data ly-ing vvithin a fixed time interval ('Window 1) around Tk, Fig. 1. The contributions to the strength of this line arising from the presence of neighboring spectral lines is taken into account by making appropriate assumptions concerning the spread of the neighboring lines vvhich vvill make non-negligible contributions (Window 2). Fig. 2. The spectrum is calculated by proceeding from the lovv end tovvards the high end of the response, again making appropriate assumptions. Finally, the erude distribution of spectral lines obtained in the first pass is improved by itera-tion. As mentioned previously, we shall demonstrate our method on hand of the shear relaxation modulus. The method can be eas-ily adapted to dealing vvith the tensile relaxation modulus or the tensile creep compliance or, for that matter, vvith any other response to the imposition of a strain or a stress as a step funetion oftime1''. 3. The Algorithm The theory of linear viscoelastic behavior describes the shear relaxation modulus by the relation G(t) = {Ge} + (Gg-{Gc})jh(r) exp-t/ r d ln r (1) u vvhere h(T) is the (normalized) continuous relaxation spectrum, and Gg and Ge are the instantaneous and the equilibrium modulus. respectively. The braces signify that (Gf| = Ge vvhen the modulus describes an arrheodictic material, and that {GJ = 0 vvhen the material is rheodictic (The term rheodictic refers to a material shovving steady-state flovv). Dividing by Gf - (GJ yields g(t) = {ge} +Jh(x) exp-t/ t d In t, (2) n vvhere g(t) and ge are the normalized relaxation modulus, and equilibrium modulus, respectively. The source data are assumed to be available as a set of M discrete data points GJ e {tj, G(tj); j = 1, 2..... m}. (3) Each of these data points can be normalized by the difference betvveen the largest. G,, and the smallest point, GM, to yield the set g, e {t,. g(tj); j = l, 2..... M}. (4) Novv. the modulus can be expressed alternatively by a discrete set of (normalized) spectral lines, h,. In terms of these we have G(t) = {G,} + (G(-{Ge})£hi exp-t/r, (5) or, in normalized form, g(.) = K} + Sh, exp-t/x,. (6) i i We intend to determine, from the set of source data, I g,; j=l, 2, 3, ... , M}, a set of spectral lines, j h,; i=l, 2, 3.....n), vvhich vvill faithfully reproduce the modulus, G(t). In proceeding to ex-plain hovv this is done, we initially use the continuous represen-tation (2), instead of its discrete equivalent (6), because this sim-plifies the presentation. We begin by splitting the integral in equation (2) to obtain r, g(t) = {c.} + Jh(r) exp-1/ r d In r + ,n r, +Jh(t) exp-t/r d In T+h(rk) exp-t/rk + r, +Jh(r) exp-t/ t d In z. (7) Let the kernel in the integrals in equation (7) be represented by the funetion K(t)=exp-t/i\ Figure 1 shovvs a plot of Kk = exp-t/Tk as a funetion of ln t/Tk. The broken line represents the tangent to Kk(t) at ln t/Tk = 0. From the behavior of the kernel vve can dravv tvvo immediate conclusions. First, Kk(l00xk) = 3.72x10 44 = 0. (g) and t, =- and In In k+1 r, (9) (10) If P is the number of spectral lines per decade of log t, then we have log t Ul - log Tk and, therefore. t, =- 2.303 P(l0l/p-l) (11) 1.0 - .6 (WIND0W1)K (WIND0W 2)k \ / > \ \ m L / X / \ \ v ' ■ \ ■■l«|II \ ■ i * V- , >— h-i log (t/r) Figure 2. Definition of Window 2 (the Modeling Window) Slika 2. Definicija Okna 2 (Okno modeliranja) 2.303x10'p t„ =-7-r- T. p(lO'" -l) Therefore, the spread of Window 2 is given by [t„ tj. The width of the window decreases as the number of spectral lines inereases. Accordingly. lim Window 2 = 0 p—♦<» (13) -.5 0 1.5 log (II*) Figure 1. Definition of Window I (the Boundary Window) Slika 1. Definicija Okna I (Mejno okno) Second, the logarithmic time interval -0.60.4 defines the region over which Kk(t) shows its largest time dependence. We call this interval Window I. The interval from vvhich data points will be drawn from the set of source data to calculate the k"' spectrum line depends on the number of lines per logarithmic decade of time. We call this interval Window 2. Figure 2 shows dtran(0 = K'(t), the derivative of K(t), as a funetion of In t/-rk for t = tk-i. tk- and tk+1. The interseetions of the central derivative with its neighbors to the left and to the right, designated by t, and t„ (for lower and upper) are given by marks the transition from the diserete to the continuous form of the representation of G(t). The third conclusion we can draw from Fig. 1 states that Window2 must not be largerthan Window 1. Outside of Window I the k"1 spectral line cannot handle data points because its con-tribution on the right is virtually zero, and on the left it ap-proaches a constant value. On the other hand, we must have at least one diserete data point lying vvithin Window 2. This real-ization allovvs us to determine the optimum number of spectrum lines per logarithmic decade to be chosen. This number must be sueh that the width of Windo\v 2 approaches (but does not ex-ceed) the vvidth Window 1. The number can be found by solving the transcendental equations (11) and (12). In equation (11) we let log t,/tk=-0.6, in equation (12) we let it be log t/rk=0.4, and solve for the nearest integer P. The smaller of the tvvo solutions gives the desired optimum number of spectral lines per decade. We can now return to equation (7). Each datum point within Window 2 pertaining to the k"1 spectrum line can be modeled as T gj = {gc} +Jh(x) exp-t|/x d In T + o \ + Jh(t) exp-tj/T d /«x + h(tk) exp-t(/tk + T + |h(x) exp-tj/x d In i. (14) By our first conclusion from Fig. 1. if t, < tk/100, then the first integral in (14) vanishes. Proceeding now from the continuous to the diserete representation. we vvrite gj={gM}+Zhi eKp-tj/T.+h, exp-t|/xk + i=m + ^ h, exp-t|/il +Ar (15) i=k+i In the equation above m is the diserete counterpart of Ta and is given by m=k-2n-l. The term A, has been added to take into account the absolute error introduced by svvitching from the continuous to the diserete representation. Using the abbreviation i=k i ^(tJ) = {gM}+ Xhlexp-tj/r + i=m + hk exp-ti / rk + ^ h, exp-t; / r,. (16) where Z(tj) denotes the theoretical (error-free) value of the nor-malized experimental datum point gr the term A! can be ex-pressed as (17) To avoid the instability problems caused by the large differ-ence in magnitude of modulus on both sides of Window 2 we in-troduce the relative error of approximation (for details see17) (18) The sum. Qk, of the squares of d, within Window 2 is j=\u = 2>2 (19) j=\ where su and skli are the first and the last discrete points in Window 2 belonging to the k"1 spectrum line. Minimizing the er-ror according to d h, (20) where hk is the k"' spectrum line. leads to the expression from which the strength of the k"1 spectral line is to be obtained, g I 2- T—F8' exp-t,/rL =0. (21) " No] Using the abbreviation -^r(tJ) = {gM}+Shi exP-lA + + exp-tj/Ti, (22) i=k+l and equation (16) in addition, we obtain the equation X | + hk exp-tJ/xk] [■nj + hk exp-ti / TkJ r g, exp-tj / Tk = 0. (23) This equation must be satisfied to minimize the sum of squares of the relative quadratic errors. Qk, in the k'h window. It must be solved numerically by iteration for it cannot be made ex-plicit for hk as is possible vvhen minimizing the absolute error17". We therefore now recast (23) in a form in vvhieh it be-comes suitable for recursion. We let hm' = "hf where hk'" is the value of hk obtained in the i"' iteration, ^ [^--^(0] g, exp-t,/xk and J3 h' = X -^(tj + hi'1 exP-t,/Tk [exp-t, / xkj ^(tjj + hf exp-tj/xk (24) (25) (26) The starting set of spectrum lines for the iteration is (V = h("(Tk); k=l,2.....n). This starting set is obtained in a first svveep through the data using the method detailed in the paper17. For subsequent svveeps we then use the method just deseribed. The iteration is broken off vvhen an appropriately chosen limit of the absolute differenee betvveen the new and the preceding square relative error, Qk. has been reached. 4. Results We demonstrate the povver of the algorithm by obtaining the distribution of spectral lines from the experimental data obtained by Catsiff and Tobolsky"' - CT data. The data has been first con-verted as deseribed in the introduetion. In this form they are tab-ulated in'7, and are presented in Fig. 3. The solid line represents a spline funetion21 through the data. Figure 4 shovvs the relax-ation spectrum Hct(t) calculated from these data, using the presented method. The reconstruction of the shear relaxation modulus, G(t), from the calculated spectrum Hct(t), using the relation (5) is compared vvith the spline funetion through the original experi-mental data in Fig. 5. Both curves can not be distinguished with-in the resolving povver of the plot. Data of CT raw data - smoothed data 6 — O O O) o -12 -10 _L I i I _L -6 -4-2 0 2 4 6 Log t - second Figure 3. Relaxation modulus, log G,tU), as funetion of log t Slika 3. Relaksacijski modul, log G(l(t). kot funkcija log t -S -4 -2 0 2 4 S 8 Log r- second Figure 4. Relaxation spectrum, log Hct(t), as funetion of log t Slika 4. Relaksacijski spekter, log H^t), kot funkcija log t Knovving the spectrum H(t) in addition to the viscoelastic constants one can generate from it the response to any desired type of excitation. In order to demonstrate this we first calculate G' O) o —t -10 , | i | , | , | i j , | i | i 1 ' ! 1 1 1 - . Data of FGF - • raw data — smoothed data -t 10 t 2 3 4 5 G 7 8 Log co - 1/second Figure 6. Storage compliance, log J'HilKo), and loss eomplianee. log J"KG[ (w), as functions of log 0) Slika 6. Shranitveni dinamični modul log J(co) in dinamični modul izgub log J"F,0F(w) kot funkciji log co J"( 0)) = [G'H];+[G>f J'ct(C0) and J"CT(co) are compared with the experimental data obtained by Fitzgerald, Grandine, and Ferry22 -FGF data. The comparison is presented in Fig. 7. J'CT(co) and J"CT((0), repre-sented as broken lines, are compared with the spline functions through the experimental data, J'KUF(co) and J"R1F(C0). shown as solid line. The original FGF data are shovvn in Fig. 6. The agree-ment betvveen the prediction and the spline function through the experimental data is excellent. cd O -5 E -6 ll O cn o ll O U_ -10 1 1 1 1 1 1 - \ - 7 Data of FGF - — smoothed data V -— auto - prediction from CT spectra ",l,l,l,i,i,i,i, ,l,l," 2 3 4 5 8 Log co — 1/second 10 Figure 7. Smoothed J<1F(co) and J"FGF(co) data compared vvith J'( ,(co) and J",T((o) derived from H,t(t) Slika 7. Zglajeni eksperimentalni krivulji log J'fgf(w) 'n J"h»(o» primerjani z log J'n(a>) in log J"ct(co) izračunanimi iz Hcx(ti 5. Conclusion In this paper we have presented the algorithm for evaluation of relaxation line spectra from experimental data. The algorithm can be easily modified for the assessment of retardation spectra. The algorithm essentiafly utilizes the fact that the kernel functions resemble step functions. Slightly different codes are used for each kernel function. as presented elsevvhere17. We feel that vve have demonstrated that the proposed algorithm is indeed capable of generating the underlying line spectra from the experimental data vvithout producing a negative lines that are physically unacceptable. Acknowledgements - The authors gratefullv acknovvledge support of this vvork by the Slovene Ministry of Science under Grant P2-1131-782, and partial support by the California Institute of Technology. G V, (CD) and G"c-T(co) can be than interconverted to storage compliance, J'rT(K>). and loss compliance, J", ,(co). using the relations1' r(o>) = cr(co) j'H]j+[g»]j (29) 6. References 1 N. W. Tschoegl, The Thcory of Linear Viseoelastic Behavior. Springer-Verlag, Heidelberg, i 989; a. Chapter6; b. Chapter4; c. Chapter 3, Section 6; d. Chapter 11: e. Chapter 6. Section 3; f. p.40: g. p. 329. : J. D. Ferry. Viseoelastic Properties of Polvmeric Materials, Wiley & Sons. NY, 1980. I A.V. Tobolsky and K. Murakami, J. Pohmer Sei. 40, 443.1955. 4 R. A. Schapery, Proč. Fourth U.S. Nat. Congr. Appl. Mech., 2. 1075. 1962.' 5 T. L. Cost and E.B. Becker, lntern. J. Num. Methods in Eng.. 2. 207, 1970. 6 H. M. Laun, J. ofRheol.. 30, 459. 1986. 7 M. Baumgaertel and H.H. Winter, Rheol. Aeta, 28, 501, 1989. s J. Honercamp, Rheol. Aeta. 28, 363, 1989. 9 J. Honercamp and J. Weese, Macromolecules, 22,4372, 1989. 10 V. M. Kamath and M.R. Mackley, J. Non-Newt. FluidMech.. 32. 1199, 1989. II C. Elster and J. Honercamp, Macromolecules, 24, 310, 1991. 12 C. Elster and J. Honercamp, J. Rheol., 36, 911, 1992. 13 N. Orbey and J.M. Dealy../. Rheol. 35, 1035, 1991. 14 J. Janzen and J. Scanlan, Submitted to Rheol. Aeta.. 1992. 15 NAV. Tschoegl and I. Emri, Submitted to Rheol. Aeta., 1994. 16 E. Catsiff and A.V. Tobolsky, Colloid Sei., 10, 375, 1955. 17 I. Emri and N.W. Tschoegl, Rheol. Aeta. 32, 311, 1993. 1K N.W. Tschoegl and I. Emri, Rheol. Aeta. 32, 322, 1993. 19 N.W. Tschoegl and I. Emri, lntern. J. Polymeric Mater., 18,95, 1992. 20 I. Emri and NAV. Tschoegl, Rheol. Aeta. 33, 60, 1994. 21 C. de Boer, A Practical Guide to Splines. Springer-Verlag, NY, 1978. 22 R. E. Fitzgerald. L.D. Grandine Jr. and J.D. Ferry, J. Appl. Phys„ 24, 650, 1953.