- Email: [email protected]
Numerical analysis in nonlinear viscoelasticity
NUMERICAL ANALYSIS IN NONLINEAR VISCOELASTICITY Y. M. HADDAD Westinghouse Canada Limited, Hamilton, Ontario, CanadaUN
3K2
Ah&a&-The well-known Volterra integral equation in linear viscoelasticity is treated to characterize the response of nonlinear viscoelastic materials. In this context, a procedure is presented which includes a differential approximation of the time-dependent kernel and an optimization procedure and leads to an optimum selection of certain parameters such that the deviation of the theoretical results form the ex~~rnent~ data is minimal. The formulations are then applied nume~c~ly to the case of natural cellulose tibres by using available experimental data on the relaxation of such fibres. The results of the investigation prove the close accuracy of the predictive ability of the model within the range of strain levels considered. 1. INTRODUCTION
FROMA phenomenofogical point of view, the concept of the relaxation function in linear viscoefasticity can be expressed in terms of Voftena f I] response equation as follows:
~(t)=+‘(r)-j-’ R”(t-+?(r)dr] -m
(1)
where f(t), e’(t) are the stress and strain tensors respectively at time t, E” is the elastic tensor modulus and R”(c - T) is the relaxation function which is usually a positive decreasing function of time. 2. NDNLINEAR
RHEOLOGICAL
RESPONSE
While eqn (1) represents the linear viscoelastic response, the nonlinear response behaviour of real material is significant [2,3]. Thus, for the formulation of the response in this case, the method first outlined by Distefano and Todeschinif41 is applied such that the Volterra equation read now for the uniaxiaf case as follows:
TO)= de(t), h, kz,. . .) +[
h(e(T), bl, b2,. . .)R(t - 7) dr
(2)
in which t(t), e(t) are the scalar components of stress and strain respectively at time t in the sample subjected to uniaxiaf loading. In this equation, the function g(.) corresponds to the elastic response in a parametric form, whereas the function h(.) accounts for the not&near hereditary effects and is also given in a parameteric form. In these two functions, unknown constants kg, kZ,. . . and b,, b2,. . . are respectively involved. The choice of the functions gf.) and hf.) is guided by a qualitative knowledge of the viscoefastic behaviour of the material. me function R(.) is the relaxation function for the uniaxiaf case. In a strict sense, this function is assumed to satisfy an Nth order differential equation with constant coefficients of the following form:
aa -I- alR"'
+ u2R"' + -
’ * +
u~_~R’~-“+ RcN’= 0
(3)
in which a~, al,. . . , @,_I are unknown coefficients. In order to simplify the analysis, one can introduce for the instan~neous response a single elastic constant, E, only. Let, also, the constant strain being applied to a specimen during a relaxation experiment be denoted by *ei, where i = 1,2,. . . , II (n corresponding to the number of experiments). Then, eqn (2) can be rewritten in a reduced form for a ffxed time t within the range of the whole relaxation experiment in the following manner: *(i(t) = E*ei + h(*ei, bf, 62,. . .,[ where the oversign * is used to indicate the theoretical stress, 325
R(ddT
(4)
Y.MHADDAD
326
In order, however, to find the unknown constants b,, bZ,. . . and to solve for the relaxation kernel R(T), it is convenient to adopt a minimization procedure as outlined below. In this context, the adopted procedure will be presented whilst an actual numerical evaluation for the nonlinear rheological response of natural cellulose fibres is presented in Section 4.
3.APPROXIMATESO~UTlONOFTH~RE~AXATiONKERNEL Whilst eqn (4) has been shown in a general form for the relaxation of a nonlinear viscoelastic material, a solution of the kernel R(T) can only be given in an approximate form. For this purpose eqn (5), below, will be used which is based on experimental observations of the relaxation behaviour of the material. Thus, denoting the experimental stress by *& one can express analytically the experimental relaxation curves as a function of the applied stress and time for a number of ex~~ments i = f,2,. . . , n in the following manner
in which *& = E*ei is the applied stress at t = 0, Gj are constants and fli(t, m/l, M~Z,.. .) are functions of time. The shape of the experimental relaxation curve suggests the form of the function &(.) for the type of material under consideration. In eqn (5) the constants Gi and mil, mi2,. . . are required to be determined for each experiment i by a standard fitting procedure. In order to use the experimentally available data of the relaxation of the material by means of eqn (5), it is evident that the coefficients in (3) have to be found by an optimization procedure such that the kernel &R(r) d7 in eqn (4) is approximated by functions &(t, mil, mi2,. . .) of eqn (5). Thus, using this approximation, it may be convenient to identify d R(7) * z 0i(t, mil, Rli2,.. . where the time T represents the to&I time for each relaxation experiment. By substituti~ the values of I’i co~es~ndi~ to eqn (6) into (3) and using the method of differential approximation[5], we require that the functional:
8 j{CZ(Jri t Q,rpt a2ri2J t . a. + +,rjN-lJ+rfNq2 dt be a minimum with respect to all possible choices of the coefficients flo, aI,. . . , &&l. The limits of the integral in the above equation and the sequel are taken from 0 to T. The minimization of expression (7) with respect to the Uju = 0,1, . . , , N - 1) leads to a system of N-simultaneous linear algebraic equations which can be written in the following form:
2 \ rji){aori t @fl)+
a2rj2)t e. . + aN_,rjN-y
dt = - 2 i-l
I
prfN)dt
wherei, ,..., n,andj=l,2 ,,.., N-l. This system of N-simultaneous algebraic equations is to be solved for the coefficients Uj(j=O,l,+*., N - 1) of eqn (3). Having determined the coefficients of the Iinear differential eqn (3), one solution of this equation may be assumed to be of the form R(T)= exp(E). In order that this solution satisfies the linear differential eqn (3) identically, one must have the characteristic equation:
uota~Ftu~F2t***+~N_~F~-'+F~=O.
(9)
This expression yields N-roots FI(I = 1,2,. , . , N) which correspond to the N-solutions of the original linear differential equation. Hence, the general solution can be written as a linear
327
Numerical analysis in nonlinear viscoelasticity
combination of the N-solutions of the form: (10)
N) are constants to be determined. where4(1=1,2,..., Now, by combining eqns (4) and (10) it follows that *&(t) =
E*ei + h(*ei, bl, b2, . . .)
‘Q=$!$-,I=l,2
,...,
iI ^DI{exp(Frf)- 11
Nandi=1,2
,...,
n).
(11)
The identification problem of the nonlinear viscoelastic response of the material can now be formalized in the following manner: Given A independent experimental functions *&(i = 1,2,. . . , n), eqn (3, corresponding to applied strains *ei, we wish to find the constants bl, bz,. . . and *& (I = 1,2,. . . , N) in the constitutive eqn (11) such that the functional:
Wr, b2,. . . , ^&I= 2 Yi I( “4(t)- *&(t))‘dt i=l
f
(12)
is minimized. For the purpose of minimizing this functional, a quadratic expression of the difference between the theoretics stress ^&t), (i = 1,2,. . . , n) which is given by eqn (If) and the experimental stress *[i(t) as given by eqn (5) is used and where +yi(i = 1,2,. . . , n) represent suitable positive weighting factors. Once the values of the unknown constants of eqn (11) have been determined by the above procedure, this equation may be used as a “model equation” for the rheological behaviour of nonlinear viscoeiastic materia1. 4. NUMERICAL
EVALUATION
The foregoing procedure is appfied for the case of natural cellulose fibres. The relaxation data as derived from the experimental relaxation curves of these fibres are first fitted to a convenient form of the analytical expression (5). Secondly, numerical schemes are developed for the differential approximation of the relaxation kernel as discussed above and, hence, for the minimization of the objective function (12). 4.1 Fitting the experimentaldata of relaxation The ex~~rnen~ curves, Meredith[2], showing the relaxation stress against the loga~thm of time are represented in Fig. 1 for natural cellulose fibres (cotton) in tension.
IO‘ ’
10 103 TIME,s
IOf
Fig. I. Experimental relaxation curves of natural cellulose fibres (cotton), from Meredith@].
328
Y. M. HADDAD
Referring to expression (5) and the shape of the relaxation curves of Fig. 1, one may express the relaxation stress of natural cellulose fibres as a function of time in the following manner *[i(t) = *& - Gitmi
(i = 1,2,. . . ) 5)
(13)
where *& is the instantaneous stress at time t = 0 and Gi and mi are constants. It is thus required to fit for each relaxation experiment the parameters Gi and mi in eqn (13) such that the relaxation stress *g(t) gives a reasonable representation of the outcome of the relaxation experiment. The values of these constants, Haddad[3], are given in Table 1. 4.2 Differential approximation of the relaxation kernel With the foregoing obtained information from the experimental relaxation data of natural cellulose by means of eqn (13), one can proceed using the method of differential approximation indicated in Section 3 to determine the coefficients of the linear differential eqn (3). Comparing eqns (5) and (13) it is evident that for the case of natural cellulose the function ei(.) takes the following form fli(.) = - tmi.
(14)
Consequently, one can write in view of eqns (6) and (14) that d Ii =z(-tmi)
(15)
from which the jth derivative required for eqn (8) can be written as rY) =
-j$ (-
t”i) = - mi(mi - l)(mi
- 2) . . .
trni _ j)tmi-i-1
(j = 0, 1,2, . . . , N - 1). (16)
A straightforward matrix inversion method has been used to solve the system of the simultaneous algebraic eqns (8) for the coefficients aj(j = 0, 1,. . . , N - 1). However, it should be mentioned that with reference to eqn (16) and in view of the values of Mi (i = 1,2,. . . ,5) indicated in Table 1, the singularity of the terms of (8) prevents the use of the lower limit of the integration of this equation to be zero. This problem has been solved here by taking the lower limit of the integral 0.09 hr. The results of the computations are presented in Table 2 for N = 3,4, . . . ,8. Table I. Material parameters characteristic eqn (13)-cotton fibres
i
*ei cm/cm
*E 108dyne/cm*
lo8dyne/(cm* hr mi)
mi
1 3 4 5
0.01 0.02 0.03 0.04 0.05
9.7 4.5 13.6 17.9 22.05
4.28641 1.91883 5.68172 6.90827 6.90621
0.06405 0.06351 0.06612 0.05839 0.06527
Gi
Table 2. Coefficients of the linear differential eqn (3)-cotton fibres N ai a0
01 a2 a3 a4 as a6 al
3
4
5
6
7
0.426757x Iti 0.441519x lo’ 0.4%572x 102
0.629616x 10’ 0.988387x ltr 0.193098x lti 0.891147x 102
0.116839x 106 0.249313x 106 0.7265%x 10’ 0.562153x 10” 0.140015x lo-’
0.261782x 10’ 0.708934x 10’ 0.281851x IO’ 0.318768x 106 0.130205x 10s 0.202358x 103
0.687618x 10s 0.225801x IO9 0.115419xlo9 0.176208x 10s 0.103647x 10’ 0.260282x lo-’ 0.276144x lo’
8 0.207167x 10” 0.799109x 1O’O 0.503567x 10” 0.984242x lo9 0.773721x 108 0.276726x IO’ 0.469378x Iti 0.361374x IO’
329
Numerical analysis in nonlinear viscoelasticity
Having determined the coeficients Uj (j = 0, 1, . _. , N - 1) it is possible to continue to find the roots of the characteristic eqn (9). In this context, the Secant method has been employed and the numerical values of the roots F, (I = 1,2, . . . , N) of eqn (9) are given in Table 3. 4.3 Minimization of the objective function (12) We continue to determine the values of the unknown constants of the model eqn (11)where the function h(.) appearing in this equation is assumed to take, for the present case of natural cellulose, the following form h(.) = exp (be).
(17)
Substituting for h(.) and &(.),as given by (17)and (14),respectively, into (11)and (5), one can write the functional (12) in the following manner Il{b, ‘I& (I = 1,2,. . . , N)} =
g %
I{ GP - @xp (b*ed)$,*DI (exp (F$) - 1))2df.
The computations analysis for minimizingthe objective function (18) has been carried out by using the steepest descent method whereby the lower limit of the integral of this equation has been taken as 0.09hr as previously mentioned. The values of the parameters b and “Dr (I = 1,2,. . . , N) which minimize eqn (18)are shown in Table 4 for N = 3,4, . . . ,8. 4.4 Testing the predictive ability of the model
To test the predictive ability of the model, the predicted values of stress have been computed for each experiment i, at different values of t between 0.09 and 2 hr and for N=3,4,..., 8. This has been carried out by using eqn (I 1)whereby the function h(.) appearing in this equation is given by eqn (17).The predicted values of stress have been then compared with the experimental ones of the same time. The mean square error, [M.S.E.], has been Table 3. Roots of the characteristic eqn (9)-cotton fibres N fi l/hr
4
3
5 6
- I .099525 - 10.~
fi Fl Fs F6 fi Fs
- 38.~35
6
5
-0.740515
- 0.554445 - 4.413529 - 15.347084 - 38.149942 - 81.550000
-6.30&X% - 22.7~725 - 59.359363
8
I
- 0.443661 - 3.325409 - 11.220736 - 27.226782 - 55.507862 - 104.633550
- 0.370971 - 2.640529 - 8.658291 - 20.638786 -41.212553 - 74.250533 - 128.372337
Table 4. Values of the parameters characteristic eqn (ll)-cotton
-0.319748 -2.179361 - 6.953434 - 16.316638 - 32.142377 - 56.795098 -94.051185 - 152.616166
fibres
N Dr 108d~e/cmz
3
4
5
a,
64.848773
- 1.279786
0.681772
:;
61.622924 85.071193
99.129241 756aO531
98.657634 75.206427
34.74769
34.276117 17.350358
:; :;
6 - 0.191855 74.94373s 98.424634 34.043117 17.350358 20.684165
7 - 1.472124 74.828227 98.368221 33.988159 17.295414 20.629207 24.108952
*& b
- 6.49227
-7.01154
- 9.51815
- I3.04031
- 14.58507
8 - 2.125503 74.793207 98.355531 33.976572 17.283834 20.617627 24.097379 27.570524 - 14.92895
Y. M. HADDAD
330
computed for each experiment at each value of N through the relation: [M s E 3 = 2 { ^I(G) - *‘D/)12 . . . (19) X (I = 1,2,. . . , X, tx within the domain of T). Table 5 shows the mean square error for each experiment i(i = 1,2,. . . ,5) as associated with an order of optimization N. The total mean square error reaches a minimum value of N = 4. The predictive ability of the model is also demonstrated in Fig. 2 for strain levels 0.01 and 0.02.
8.0
t
-
EXPERIMENT
-----
PREDICTION
N=3
4.0 N=7
N=B
3.0
O.OL 0.0
I
I
0.5
I
I
1.5
1.0
2.0
TIME, h
Fig. 2. Relaxation curves of naturalcellulose iibres (cotton).
Table 5. Mean square error Mean square error, [ IO8dyne/cm212
i
*et
I 2 3 4 5
0.01 0.02 0.03 0.04 0.05
Total mean square error
N=3
N=4
N=5
N=6
N=l
N=8
0.088020 0.506500 3.302700 5.690000 3.598600
0.026856 0.061094 1.177800 2.24380 0.46755
0.034400 0.073541 1.524300 3.598600 I.816100
0.060531 0.1 II610 2.355300 6.396000 5.457600
0.162840 0.135730 I.890500 6.001600 5.644300
0.42276 0.43198 0.89025 3.7243 3.41%
2.637164
0.79542
I .4093882
2.8762082
2.166994
I .I17116
5. CONCLUDING
REMARKS
The nonlinear viscoelastic response of real material has been considered in a generalized manner by using the hereditary relaxation relation due to Volterra. The analytical approach due to Distefano and Todeschini has been adopted. In this context, a procedure has been presented that leads to an optimum selection of certain parameters that characterize the relaxation response of the material so that the deviation of the theoretical results from the experimental data would be minimal. The theoretical model has been applied to the relaxation of natural cellulose fibres in tension where the experimental data are available due to Meredith. The results of the investigation prove the close accuracy of the predictive ability of the model within the range of strain levels considered, i.e. *ei = 0.01,0.02, . . . ,0.05. In this range of strain, the total mean square error reaches a minimum value at N = 4 for the time range between 0.09 and 2 hr.
Numerical analysis in nonlinear viscoelasticity
331
Acknowledgements-The results in this paper are taken from the author’s research work at McGill University, Montreal, Canada. This work was supported by the Pulp and Paper Research Institute of Canada. The negotiated development grant in building research at the Centre for Building Studies, Concordia University, Montreal is also gratefully acknowledged. REFERENCES 111V. VOLTERRA,~o~c?jonsdeLigaes. Gauthier Villard,Paris (1913).See Y. N. RABOTNOV,Creepjobless in ~~c~u?a~ Members, p. 123.Wifey, New York (1969). 121R. MEREDITH. J. Text. fnns.45,438 (1954). [31 Y. HADDAD, Response behaviour of a two-dimensional fibrous network, Ph.D. Thesis. McGill University, Montreal (1975). (41 N. DISTEFANO and R. TODESCHINI, Inf. 1. Solids Stmctures 9,t305(1973). [5] R. BELLMAN, Methods of nonlinear Analysis. Academic Press, New York (1970).
3K2
Ah&a&-The well-known Volterra integral equation in linear viscoelasticity is treated to characterize the response of nonlinear viscoelastic materials. In this context, a procedure is presented which includes a differential approximation of the time-dependent kernel and an optimization procedure and leads to an optimum selection of certain parameters such that the deviation of the theoretical results form the ex~~rnent~ data is minimal. The formulations are then applied nume~c~ly to the case of natural cellulose tibres by using available experimental data on the relaxation of such fibres. The results of the investigation prove the close accuracy of the predictive ability of the model within the range of strain levels considered. 1. INTRODUCTION
FROMA phenomenofogical point of view, the concept of the relaxation function in linear viscoefasticity can be expressed in terms of Voftena f I] response equation as follows:
~(t)=+‘(r)-j-’ R”(t-+?(r)dr] -m
(1)
where f(t), e’(t) are the stress and strain tensors respectively at time t, E” is the elastic tensor modulus and R”(c - T) is the relaxation function which is usually a positive decreasing function of time. 2. NDNLINEAR
RHEOLOGICAL
RESPONSE
While eqn (1) represents the linear viscoelastic response, the nonlinear response behaviour of real material is significant [2,3]. Thus, for the formulation of the response in this case, the method first outlined by Distefano and Todeschinif41 is applied such that the Volterra equation read now for the uniaxiaf case as follows:
TO)= de(t), h, kz,. . .) +[
h(e(T), bl, b2,. . .)R(t - 7) dr
(2)
in which t(t), e(t) are the scalar components of stress and strain respectively at time t in the sample subjected to uniaxiaf loading. In this equation, the function g(.) corresponds to the elastic response in a parametric form, whereas the function h(.) accounts for the not&near hereditary effects and is also given in a parameteric form. In these two functions, unknown constants kg, kZ,. . . and b,, b2,. . . are respectively involved. The choice of the functions gf.) and hf.) is guided by a qualitative knowledge of the viscoefastic behaviour of the material. me function R(.) is the relaxation function for the uniaxiaf case. In a strict sense, this function is assumed to satisfy an Nth order differential equation with constant coefficients of the following form:
aa -I- alR"'
+ u2R"' + -
’ * +
u~_~R’~-“+ RcN’= 0
(3)
in which a~, al,. . . , @,_I are unknown coefficients. In order to simplify the analysis, one can introduce for the instan~neous response a single elastic constant, E, only. Let, also, the constant strain being applied to a specimen during a relaxation experiment be denoted by *ei, where i = 1,2,. . . , II (n corresponding to the number of experiments). Then, eqn (2) can be rewritten in a reduced form for a ffxed time t within the range of the whole relaxation experiment in the following manner: *(i(t) = E*ei + h(*ei, bf, 62,. . .,[ where the oversign * is used to indicate the theoretical stress, 325
R(ddT
(4)
Y.MHADDAD
326
In order, however, to find the unknown constants b,, bZ,. . . and to solve for the relaxation kernel R(T), it is convenient to adopt a minimization procedure as outlined below. In this context, the adopted procedure will be presented whilst an actual numerical evaluation for the nonlinear rheological response of natural cellulose fibres is presented in Section 4.
3.APPROXIMATESO~UTlONOFTH~RE~AXATiONKERNEL Whilst eqn (4) has been shown in a general form for the relaxation of a nonlinear viscoelastic material, a solution of the kernel R(T) can only be given in an approximate form. For this purpose eqn (5), below, will be used which is based on experimental observations of the relaxation behaviour of the material. Thus, denoting the experimental stress by *& one can express analytically the experimental relaxation curves as a function of the applied stress and time for a number of ex~~ments i = f,2,. . . , n in the following manner
in which *& = E*ei is the applied stress at t = 0, Gj are constants and fli(t, m/l, M~Z,.. .) are functions of time. The shape of the experimental relaxation curve suggests the form of the function &(.) for the type of material under consideration. In eqn (5) the constants Gi and mil, mi2,. . . are required to be determined for each experiment i by a standard fitting procedure. In order to use the experimentally available data of the relaxation of the material by means of eqn (5), it is evident that the coefficients in (3) have to be found by an optimization procedure such that the kernel &R(r) d7 in eqn (4) is approximated by functions &(t, mil, mi2,. . .) of eqn (5). Thus, using this approximation, it may be convenient to identify d R(7) * z 0i(t, mil, Rli2,.. . where the time T represents the to&I time for each relaxation experiment. By substituti~ the values of I’i co~es~ndi~ to eqn (6) into (3) and using the method of differential approximation[5], we require that the functional:
8 j{CZ(Jri t Q,rpt a2ri2J t . a. + +,rjN-lJ+rfNq2 dt be a minimum with respect to all possible choices of the coefficients flo, aI,. . . , &&l. The limits of the integral in the above equation and the sequel are taken from 0 to T. The minimization of expression (7) with respect to the Uju = 0,1, . . , , N - 1) leads to a system of N-simultaneous linear algebraic equations which can be written in the following form:
2 \ rji){aori t @fl)+
a2rj2)t e. . + aN_,rjN-y
dt = - 2 i-l
I
prfN)dt
wherei, ,..., n,andj=l,2 ,,.., N-l. This system of N-simultaneous algebraic equations is to be solved for the coefficients Uj(j=O,l,+*., N - 1) of eqn (3). Having determined the coefficients of the Iinear differential eqn (3), one solution of this equation may be assumed to be of the form R(T)= exp(E). In order that this solution satisfies the linear differential eqn (3) identically, one must have the characteristic equation:
uota~Ftu~F2t***+~N_~F~-'+F~=O.
(9)
This expression yields N-roots FI(I = 1,2,. , . , N) which correspond to the N-solutions of the original linear differential equation. Hence, the general solution can be written as a linear
327
Numerical analysis in nonlinear viscoelasticity
combination of the N-solutions of the form: (10)
N) are constants to be determined. where4(1=1,2,..., Now, by combining eqns (4) and (10) it follows that *&(t) =
E*ei + h(*ei, bl, b2, . . .)
‘Q=$!$-,I=l,2
,...,
iI ^DI{exp(Frf)- 11
Nandi=1,2
,...,
n).
(11)
The identification problem of the nonlinear viscoelastic response of the material can now be formalized in the following manner: Given A independent experimental functions *&(i = 1,2,. . . , n), eqn (3, corresponding to applied strains *ei, we wish to find the constants bl, bz,. . . and *& (I = 1,2,. . . , N) in the constitutive eqn (11) such that the functional:
Wr, b2,. . . , ^&I= 2 Yi I( “4(t)- *&(t))‘dt i=l
f
(12)
is minimized. For the purpose of minimizing this functional, a quadratic expression of the difference between the theoretics stress ^&t), (i = 1,2,. . . , n) which is given by eqn (If) and the experimental stress *[i(t) as given by eqn (5) is used and where +yi(i = 1,2,. . . , n) represent suitable positive weighting factors. Once the values of the unknown constants of eqn (11) have been determined by the above procedure, this equation may be used as a “model equation” for the rheological behaviour of nonlinear viscoeiastic materia1. 4. NUMERICAL
EVALUATION
The foregoing procedure is appfied for the case of natural cellulose fibres. The relaxation data as derived from the experimental relaxation curves of these fibres are first fitted to a convenient form of the analytical expression (5). Secondly, numerical schemes are developed for the differential approximation of the relaxation kernel as discussed above and, hence, for the minimization of the objective function (12). 4.1 Fitting the experimentaldata of relaxation The ex~~rnen~ curves, Meredith[2], showing the relaxation stress against the loga~thm of time are represented in Fig. 1 for natural cellulose fibres (cotton) in tension.
IO‘ ’
10 103 TIME,s
IOf
Fig. I. Experimental relaxation curves of natural cellulose fibres (cotton), from Meredith@].
328
Y. M. HADDAD
Referring to expression (5) and the shape of the relaxation curves of Fig. 1, one may express the relaxation stress of natural cellulose fibres as a function of time in the following manner *[i(t) = *& - Gitmi
(i = 1,2,. . . ) 5)
(13)
where *& is the instantaneous stress at time t = 0 and Gi and mi are constants. It is thus required to fit for each relaxation experiment the parameters Gi and mi in eqn (13) such that the relaxation stress *g(t) gives a reasonable representation of the outcome of the relaxation experiment. The values of these constants, Haddad[3], are given in Table 1. 4.2 Differential approximation of the relaxation kernel With the foregoing obtained information from the experimental relaxation data of natural cellulose by means of eqn (13), one can proceed using the method of differential approximation indicated in Section 3 to determine the coefficients of the linear differential eqn (3). Comparing eqns (5) and (13) it is evident that for the case of natural cellulose the function ei(.) takes the following form fli(.) = - tmi.
(14)
Consequently, one can write in view of eqns (6) and (14) that d Ii =z(-tmi)
(15)
from which the jth derivative required for eqn (8) can be written as rY) =
-j$ (-
t”i) = - mi(mi - l)(mi
- 2) . . .
trni _ j)tmi-i-1
(j = 0, 1,2, . . . , N - 1). (16)
A straightforward matrix inversion method has been used to solve the system of the simultaneous algebraic eqns (8) for the coefficients aj(j = 0, 1,. . . , N - 1). However, it should be mentioned that with reference to eqn (16) and in view of the values of Mi (i = 1,2,. . . ,5) indicated in Table 1, the singularity of the terms of (8) prevents the use of the lower limit of the integration of this equation to be zero. This problem has been solved here by taking the lower limit of the integral 0.09 hr. The results of the computations are presented in Table 2 for N = 3,4, . . . ,8. Table I. Material parameters characteristic eqn (13)-cotton fibres
i
*ei cm/cm
*E 108dyne/cm*
lo8dyne/(cm* hr mi)
mi
1 3 4 5
0.01 0.02 0.03 0.04 0.05
9.7 4.5 13.6 17.9 22.05
4.28641 1.91883 5.68172 6.90827 6.90621
0.06405 0.06351 0.06612 0.05839 0.06527
Gi
Table 2. Coefficients of the linear differential eqn (3)-cotton fibres N ai a0
01 a2 a3 a4 as a6 al
3
4
5
6
7
0.426757x Iti 0.441519x lo’ 0.4%572x 102
0.629616x 10’ 0.988387x ltr 0.193098x lti 0.891147x 102
0.116839x 106 0.249313x 106 0.7265%x 10’ 0.562153x 10” 0.140015x lo-’
0.261782x 10’ 0.708934x 10’ 0.281851x IO’ 0.318768x 106 0.130205x 10s 0.202358x 103
0.687618x 10s 0.225801x IO9 0.115419xlo9 0.176208x 10s 0.103647x 10’ 0.260282x lo-’ 0.276144x lo’
8 0.207167x 10” 0.799109x 1O’O 0.503567x 10” 0.984242x lo9 0.773721x 108 0.276726x IO’ 0.469378x Iti 0.361374x IO’
329
Numerical analysis in nonlinear viscoelasticity
Having determined the coeficients Uj (j = 0, 1, . _. , N - 1) it is possible to continue to find the roots of the characteristic eqn (9). In this context, the Secant method has been employed and the numerical values of the roots F, (I = 1,2, . . . , N) of eqn (9) are given in Table 3. 4.3 Minimization of the objective function (12) We continue to determine the values of the unknown constants of the model eqn (11)where the function h(.) appearing in this equation is assumed to take, for the present case of natural cellulose, the following form h(.) = exp (be).
(17)
Substituting for h(.) and &(.),as given by (17)and (14),respectively, into (11)and (5), one can write the functional (12) in the following manner Il{b, ‘I& (I = 1,2,. . . , N)} =
g %
I{ GP - @xp (b*ed)$,*DI (exp (F$) - 1))2df.
The computations analysis for minimizingthe objective function (18) has been carried out by using the steepest descent method whereby the lower limit of the integral of this equation has been taken as 0.09hr as previously mentioned. The values of the parameters b and “Dr (I = 1,2,. . . , N) which minimize eqn (18)are shown in Table 4 for N = 3,4, . . . ,8. 4.4 Testing the predictive ability of the model
To test the predictive ability of the model, the predicted values of stress have been computed for each experiment i, at different values of t between 0.09 and 2 hr and for N=3,4,..., 8. This has been carried out by using eqn (I 1)whereby the function h(.) appearing in this equation is given by eqn (17).The predicted values of stress have been then compared with the experimental ones of the same time. The mean square error, [M.S.E.], has been Table 3. Roots of the characteristic eqn (9)-cotton fibres N fi l/hr
4
3
5 6
- I .099525 - 10.~
fi Fl Fs F6 fi Fs
- 38.~35
6
5
-0.740515
- 0.554445 - 4.413529 - 15.347084 - 38.149942 - 81.550000
-6.30&X% - 22.7~725 - 59.359363
8
I
- 0.443661 - 3.325409 - 11.220736 - 27.226782 - 55.507862 - 104.633550
- 0.370971 - 2.640529 - 8.658291 - 20.638786 -41.212553 - 74.250533 - 128.372337
Table 4. Values of the parameters characteristic eqn (ll)-cotton
-0.319748 -2.179361 - 6.953434 - 16.316638 - 32.142377 - 56.795098 -94.051185 - 152.616166
fibres
N Dr 108d~e/cmz
3
4
5
a,
64.848773
- 1.279786
0.681772
:;
61.622924 85.071193
99.129241 756aO531
98.657634 75.206427
34.74769
34.276117 17.350358
:; :;
6 - 0.191855 74.94373s 98.424634 34.043117 17.350358 20.684165
7 - 1.472124 74.828227 98.368221 33.988159 17.295414 20.629207 24.108952
*& b
- 6.49227
-7.01154
- 9.51815
- I3.04031
- 14.58507
8 - 2.125503 74.793207 98.355531 33.976572 17.283834 20.617627 24.097379 27.570524 - 14.92895
Y. M. HADDAD
330
computed for each experiment at each value of N through the relation: [M s E 3 = 2 { ^I(G) - *‘D/)12 . . . (19) X (I = 1,2,. . . , X, tx within the domain of T). Table 5 shows the mean square error for each experiment i(i = 1,2,. . . ,5) as associated with an order of optimization N. The total mean square error reaches a minimum value of N = 4. The predictive ability of the model is also demonstrated in Fig. 2 for strain levels 0.01 and 0.02.
8.0
t
-
EXPERIMENT
-----
PREDICTION
N=3
4.0 N=7
N=B
3.0
O.OL 0.0
I
I
0.5
I
I
1.5
1.0
2.0
TIME, h
Fig. 2. Relaxation curves of naturalcellulose iibres (cotton).
Table 5. Mean square error Mean square error, [ IO8dyne/cm212
i
*et
I 2 3 4 5
0.01 0.02 0.03 0.04 0.05
Total mean square error
N=3
N=4
N=5
N=6
N=l
N=8
0.088020 0.506500 3.302700 5.690000 3.598600
0.026856 0.061094 1.177800 2.24380 0.46755
0.034400 0.073541 1.524300 3.598600 I.816100
0.060531 0.1 II610 2.355300 6.396000 5.457600
0.162840 0.135730 I.890500 6.001600 5.644300
0.42276 0.43198 0.89025 3.7243 3.41%
2.637164
0.79542
I .4093882
2.8762082
2.166994
I .I17116
5. CONCLUDING
REMARKS
The nonlinear viscoelastic response of real material has been considered in a generalized manner by using the hereditary relaxation relation due to Volterra. The analytical approach due to Distefano and Todeschini has been adopted. In this context, a procedure has been presented that leads to an optimum selection of certain parameters that characterize the relaxation response of the material so that the deviation of the theoretical results from the experimental data would be minimal. The theoretical model has been applied to the relaxation of natural cellulose fibres in tension where the experimental data are available due to Meredith. The results of the investigation prove the close accuracy of the predictive ability of the model within the range of strain levels considered, i.e. *ei = 0.01,0.02, . . . ,0.05. In this range of strain, the total mean square error reaches a minimum value at N = 4 for the time range between 0.09 and 2 hr.
Numerical analysis in nonlinear viscoelasticity
331
Acknowledgements-The results in this paper are taken from the author’s research work at McGill University, Montreal, Canada. This work was supported by the Pulp and Paper Research Institute of Canada. The negotiated development grant in building research at the Centre for Building Studies, Concordia University, Montreal is also gratefully acknowledged. REFERENCES 111V. VOLTERRA,~o~c?jonsdeLigaes. Gauthier Villard,Paris (1913).See Y. N. RABOTNOV,Creepjobless in ~~c~u?a~ Members, p. 123.Wifey, New York (1969). 121R. MEREDITH. J. Text. fnns.45,438 (1954). [31 Y. HADDAD, Response behaviour of a two-dimensional fibrous network, Ph.D. Thesis. McGill University, Montreal (1975). (41 N. DISTEFANO and R. TODESCHINI, Inf. 1. Solids Stmctures 9,t305(1973). [5] R. BELLMAN, Methods of nonlinear Analysis. Academic Press, New York (1970).