Mechanics of Time-Dependent Materials 8: 105–118, 2004. © 2004 Kluwer Academic Publishers. Printed in the Netherlands.
105
New Methods for Identifying Rheological Parameter for Fractional Derivative Modeling of Viscoelastic Behavior T. BEDA and Y. CHEVALIER1 E.N.S.P. (National Advanced School of Engineering of University of Yaounde I), B.P. 8390, Yaounde, Cameroon; E-mail:
[email protected]; 1 Institut Supérieur des Matériaux & de la Construction Mécanique (ISMCM Paris), LISMMA – Viboacoustic Group, 3, rue F. Hainaut, F-93407 Saint-Ouen Cedex, France; E-mail:
[email protected] (Received 11 February 2002; accepted in revised form 29 April 2003) Abstract. The modelling of viscoelastic behaviour by fractional derivatives, the graphical method called the TIBI Diagram that is developed in this article, is a new approach, completely different from the classical methods, which are somewhat inspired by the BODE diagram technique. The TIBI technique permits the reduction to a strict minimum of two break frequencies, which is the minimum number necessary, while the other graphical methods of identification require up to ten parameters to determine linear viscoelastic behaviour. The second method developed in this paper is a numerical approach which is very flexible with regard to the experimental complex modulus permitting the determination of rheological characteristics. One needs a group of frequencies limited only to the transition region of the material. These two original methods are combined in certain cases to obtain semi-graphical techniques. Experimental and simulation results are presented for each method with comparative analysis. Key words: fractional derivative, modelling, updating, viscoelasticity
1. Introduction The modelling of linear viscoelastic behaviour using fractional derivatives is a principal theme in modern research into structural calculus. The interest in this approach lies in the fact that this mathematical model optimizes the number of rheological parameters that define the viscoelastic behaviour, unlike the classical integer order derivatives. Recent contributions of Bagley (1983) and Bagley and Torvik (1983, 1985, 1986) , Koeller (1984) and Rogers (1983, 1984, 1989) have established a solid foundation for fractional derivative models, leadint to their increasing use in studying the rheology of viscoelastic materials. The problem that this model poses arises from the determination of relative rheological parameters, some of which are fractional power of frequency (or time). The methods presented in the literature settle the value of the power factors a priori. The method mostly used, known as the BODE Diagram, however, is not
106
T. BEDA AND Y. CHEVALIER
recommended for high damping materials since it requires the evaluation of too many parameters: about 20 (Beda et al., 1991; Soula et al., 1997), which introduce a source of errors. Even methods drawing their inspiration from BODE Diagram, methods which a priori fix the fractional power factors because of the difficulties of their evaluation, certainly reduce the parameters sought, but the number required remains high (Roger, 1983). The methods developed in this article are original graphical and numerical methods that permit the evaluation of all the parameters, including the fractional power factors, without necessarily establishing any of the rhelogical parameters. 2. The Fractional Derivative Viscoelastic Model 2.1. C ONSTITUTIVE E QUATIONS Bagley (1983), Bagley and Torvik (1983), Koeller (1984) and Roger (1983) have shown that for many thermorheologically simple and macroscopically homogenous viscoelastic materials, including a large family of macromolecular, thermohardening and thermoplastic materials and elastomers, the constitutive equation for linear behavior can be defined by the following mathematical formula: σ (t) + aD α σ (t) = E0 ε(t) + bD α ε(t),
(1)
where d 1 D f (t) = (1 − α) dt
t
α
0
f (u) du, (t − u)α
0 < α < 1, t ≡ time, (x): Gamma function . Using the Laplace–Carson transform domain, taking the Fourier transform of Equation (1), we obtain [1 + a(j ω)α ]σ ∗ (ω) = [E0 + b(j ω)α ]ε ∗ (ω),
(2)
where j 2 = −1, and σ ∗ (ω) and ε ∗ (ω) are respectively the Fourier transforms of the stress σ (t) and the strain ε(t).
NEW METHODS FOR IDENTIFYING RHEOLOGICAL PARAMETER
107
2.2. C OMPLEX M ODULUS Equation (2) gives σ ∗ (ω) =
E0 + b(j ω)α ∗ ε (ω), 1 + a(j ω)α
σ ∗ (ω) = E ∗ (ω)ε ∗ (ω),
(3)
where E ∗ (ω) is the complex Young modulus. It is also the Carson transform of the relaxation function of the material. Then E ∗ (ω) =
E0 + b(j ω)α 1 + a(j ω)α
= E0
1+ 1
b (j ω)α E0 + a(j ω)α
can be rewritten with different parameters to give the following relationship: α 1 + j ωω0 α , E ∗ (ω) = E0 1 + j ωω1 E ∗ (ω) = E0
1 + (j τ0 ω)α . 1 + (j τ1 ω)α
(4)
(5)
The Young modulus is therefore a function that depends on a fractional power of frequency. For the multiform function (j ω)α , which is a rational fraction between 0 and 1, we use the part of the principal resolution 0 ≤ θ ≤ 2π because of the property of the relaxation function (Beda, 1990). For such a material, the linear viscoelastic behaviour is completely defined by four rheological parameters, E0 , ω0 , ω1 and α, which are characteristic of the material. Our work concerns the evaluation of these four parameters based on the experimental complex Young modulus obtained over a certain discrete range of frequencies. 3. Mathematical Study of the Complex Modulus Function It is given that Log E(ω) = Log |E ∗ (ω)| = Log E0 + Log |1 + (j τ0 ω)α | − Log |1 + (j τ1 ω)α |,
(6)
where E(ω) = |E ∗ (ω)| is the gain (or magnitude) of the complex Young modulus.
108
T. BEDA AND Y. CHEVALIER
Figure 1. Magnitude of function h(ω) versus reduced frequency (- - - asymptotic value).
Equation (6) shows that the study of Log |E ∗ (ω)| can be deduced from the mathematical study of Log |1 + (j τ ω)α |. 3.1. T HE F UNCTION h(ω) = 1 + (j τ ω)α h(ω) = 1 + (j τ ω)α = 1 + (τ ω)α cos so
with
τ =
1 ωc
απ απ + j (τ ω)α sin , 2 2
(7)
απ + (τ ω)2α . 1 + 2(τ ω)α cos |h(ω)| = 2 Many curves of different values of exponent α ranging between 0.5 and 1 are drawn in Figure 1. The abscissa represents reduced frequency τ ω. Analysis of these curves, recorded in Beda (1990) and Beda et al. (1991), reveals the existence of an empirical value of frequency ωf called the fixed point, which is a frequency at which Log |1 + (j ωf )α | varies insignificantly. We conclude therefore that, for this frequency ωf , the value of the function Log |h(ωf )| remains constant whatever the value of α (0.5 < α < 1) and is equal to 7.1 dB.
NEW METHODS FOR IDENTIFYING RHEOLOGICAL PARAMETER
109
The fixed point ωf is situated at 6.16 dB from the break frequency, ωc , from which all asymptotes of different functions in bilogarithmic drawing of diagrams originate (Figure 1). Remark. Beyond two decades of frequencies, the curves are mingled with the asymptotes. 3.2. T HE F UNCTION H (ω) = (1 + (j τ0 ω)α )/(1 + (j τ1 ω)α ) With H (ω) =
1 + (j τ0 ω)α 1 + (j τ1 ω)α
(8)
if we consider the gain |H (ω)| = G(ω), we obtain Log G(ω) = Log |1 + (j τ0 ω)α | − Log |1 + (j τ1 ω)α | = Log |h0 (ω)| − Log |h1 (ω)|, where τ0 = 1/ω0 , τ1 = 1/ω1 , and ω0 and ω1 are the break frequencies of the function G(ω). 3.2.1. Inflexion Point of the Gain Curve Given Log G(ω) =
απ 1 Log 1 + 2(τ0 ω)α cos + (τ0 ω)2α 2 2 απ 1 + (τ1 ω)2α , − Log 1 + 2(τ1 ω)α cos 2 2
taking the second derivative of (9) gives N1 (ω)D2(ω) − N2 (ω)D1(ω) d2 Log G(ω) , = 2 dLog ω) D1 (ω)D2 (ω) with απ απ N1 (ω) α 2α 3α + 2(τ , = (τ ω) cos ω) + (τ ω) cos 0 0 0 2α 2 2 2 απ απ N2 (ω) + 2(τ1 ω)2α + (τ1 ω)3α cos , = (τ1 ω)α cos 2 2α 2 2 2 απ + (τ0 ω)2α , D1 (ω) = 1 + 2(τ0 ω)α cos 2 2 απ + (τ1 ω)2α . D2 (ω) = 1 + 2(τ1 ω)α cos 2
(9)
110
T. BEDA AND Y. CHEVALIER
We see d2 Log G(ω) = 0 for (dLog ω)2
ωinf = √
1 √ = ω0 ω1 . τ0 τ1
The frequency ωinf is an inflexion point for the function G(ω) in the bilogarithmic drawing of diagrams. 3.2.2. Maximum Loss Factor ηmax The complex Young modulus is expressed as E ∗ (ω) = E (ω) + j E (ω) = E (ω)[1 + j η(ω)],
(10)
where η(ω) = E (ω)/E (ω) is the loss factor. If we consider E ∗ (ω) = E0 H (ω), we obtain 1 + (τ0 ω)α cos απ + j (τ0 ω)α sin απ 2 2 (11) E ∗ (ω) = E0 α sin απ 1 + (τ1 ω)α cos απ + j (τ ω) 1 2 2 and (a) Maximum of the loss factor function η(ω) (τ0α − τ1α )[1 − (τ0 τ1 )α ω2α ] sin απ dLog η(ω) 2 = , α ω 2α ]2 dLog ω [1 + (τ0α + τ1α )ωα cos απ + (τ τ ) 0 1 2 we have dLog η(ω) = 0 for dLog ω
1 √ ωinf = √ = ω0 ω1 , τ0 τ1
so the inflexion frequency ωinf for the curve of the gain G(ω) is also the frequency where the loss factor η(ω) is maximum. (b) Calculus of ηmax ηmax
1 = η(ωinf ) = η √ , τ0 τ1
so ηmax
(1 − ξ α )tg απ 2 1 + ξα −
2ξ α/2 cos απ 2
where
ξ=
τ1 ω0 = . τ0 ω1
(12)
ξ can be neglected for materials which have large transition region (beyond 5 decades of frequencies). In this case: απ with less than 10% of relative error. ηmax ∼ = tg 2
NEW METHODS FOR IDENTIFYING RHEOLOGICAL PARAMETER
111
Figure 2. TIBI diagram.
So α≈
2 Artg ηmax . π
(13)
4. Rheological Parameter Identification 4.1. G RAPHICAL M ETHOD FROM 7.1 AND −6.16 D B: T HE TIBI D IAGRAM This method consists of sketching the two horizontal asymptotes and finding out where points V0 and V1 are situated on an asymptote at a distance of 7.1 dB from the experimental gain curve (Figure 2). The fixed points of the magnitude E(ω) are the projections of these two points on the frequency-axis. The break frequencies ω0 , ω1 are the projections on the frequency axis of the points M0 and M1 obtained by displacing −6.16 dB from V0 and 6.16 dB from V1 on the asymptotes (Figure 2). The exponent parameter is given by the slope of the line D connecting the points M0 and M1 . The parameter E0 which corresponds to the static Young modulus of the material is given by the value of the inferior asymptote. All four characteristics E0 ,
112
T. BEDA AND Y. CHEVALIER
Figure 3. Inflexion point ωinf and magnitude Einf determination.
ω0 , ω1 and α are thus determined. Consequently, we realise that the evaluation of the function E ∗ (ω) which is complex, and has two real functions, necessitates the knowledge of only one of the two real functions, i.e. the gain function E(ω). The experimental loss factor curve could be unknown. If the loss factor is known at certain experimental points, it will serve to validate the parameters calculated. 4.2. N UMERICAL M ETHOD This method is used only if relation (13) is satisfied. In this case, it will be applicable in the hypothesis of the previous method. But it is especially applied when the experimental complex modulus are known only on the transition region (no asymptotic part). The static Young modulus E0 is also known. The exponent α is given by relation (13). The inflexion frequency ωinf is determined from the experimental curve of the loss factor by the maximum ηmax . And the magnitude Einf is deduced from the gain curve at the frequency ωinf : Einf = E(ωinf ) (Figure 3). 4.2.1. Evaluation of the Break Frequencies ω0 and ω1 The equation defining the straight line passing through the points M0 and Minf is given as y(ω) = αLog ω + Log
Einf α . ωinf
NEW METHODS FOR IDENTIFYING RHEOLOGICAL PARAMETER
113
Figure 4. Determination of power factor α and inflexion point ωinf .
So Log E0 = αLog ω0 + Log
Einf α . ωinf
The break frequencies (see Figure 4) are given by ω0 =
E0 Einf
1/α ωinf ,
ω1 =
Einf E0
1/α ωinf .
(14)
4.3. S EMI -G RAPHICAL M ETHOD This method requires a knowledge of the experimental complex modulus E ∗ (ω) on the rubbery plateau (inferior asymptote) and on the transition region of the material. First Technical Method: The coefficient α is given by Equation (13): α≈
2 Artg ηmax . π
The inflexion point is given by the abscissa of the point where the loss factor is maximum: η(ωinf ) = ηmax . The static Young modulus E0 is the value of the inferior asymptote. The break frequency ω0 is the intersection of the asymptote and the right gradient α passing through the inflexion point Minf (Figure 4).
114
T. BEDA AND Y. CHEVALIER
Figure 5. Simulation in order to evaluate parameters ω0 , ω1 and α.
The second break frequency ω1 is deduced from ω1 =
2 ωinf . ω0
Second Technical Method: The break frequency ω0 is determined by the (7.1 dB, −6.16 dB) method. The maximum loss factor permits the determination of the inflexion freqency ωinf which in turn permits the graphical detection of the inflexion point Minf . The exponent coefficient α is the slope of the line connecting the points M0 and Minf . The static Young modulus E0 and the second break frequency ω1 are determined in the same way as in the first method (Figure 4). 5. Results 5.1. S IMULATION : C ASES OF F IGURES 5 AND 6 Real parameters α = 0.72, ω0 = 102 , ω1 = 108 , E0 = 1. 5.1.1. For the TIBI Diagram Technique Calculated parameters α = 0.711, ω0 = 902 , ω1 = 908 , E0 = 1. The relative error on the exponent is 0.4%. Figure 6a shows the magnitude curves and Figure 6b the loss factor curves for the real and calculated values. The loss factor curves are almost identical, which validates the calculated parameters which were obtained independently of the experimental values of loss factor.
NEW METHODS FOR IDENTIFYING RHEOLOGICAL PARAMETER
115
Figure 6. Comparative curves: - - - real curves; estimate curves.
5.1.2. For the Numerical Method Calculated parameters α = 0.717, ω0 = 102 , ω1 = 1.113×108 , E0 = 1. Estimated parameters ηmax = 2.05, Einf = 150, ωinf = 105 . The relative error of the exponent α is 1.25%, and of the break frequencies ω0 : 10% and ω1 : 11.3%. 5.2. E XPERIMENTAL R ESULTS Using the experimental magnitude of complex Young modulus of polyisobutylene (PIB) at 20◦ C obtained by Nashif et al. (1975) (Figure 7), the rheological parameters characterising the material we evaluated are E0 = 147 psi, α = 0.666, ω0 = 40 Hz, ω1 = 2.5 × 106 Hz. So E ∗ (ω) = 147
1 + 0.08571(j ω)0.666 psi, 1 + 5.482 × 10−5 (j ω)0.666
E ∗ (ω) = 1.0135 × 106
1 + 0.08571(j ω)0.666 Pa. 1 + 5.482 × 10−5 (j ω)0.666
(15)
116
T. BEDA AND Y. CHEVALIER
Figure 7. Complex modulus of polyisobutylene (PIB) at 20◦ C: (a) obtained by Nashif et al. (1975); (b) calculated by the TIBI diagram technique.
Figure 8 shows comparative curves of experimental and estimated complex Young modulus. These graphs validate the rheological characteristics determined for polyisobutylene (PIB) at 20◦ C. 6. Conclusions The TIBI Diagram Technique constitutes an efficient method to find rheological characteristics of viscoelastic materials by an appropriate graphical approach, which is a special optisation method. It permits one to reduce to a strict minimum the number of rheological parameters characterising the viscoelastic behavior. However, this grahical method requires a knowledge of experimental values of the magnitude of the complex Young modulus covering a large range of frequencies, including the rubbery plateau, the glassy plateau and the transition region of the material.
NEW METHODS FOR IDENTIFYING RHEOLOGICAL PARAMETER
117
Figure 8. Comparative results of complex Young modulus of polyisobutylene (PIB) at 20◦ C. Experimental results (Nashif et al., 1975): •, . Present results (TIBI diagram technique): , .
The numerical method developed permits one to alleviate the difficulties inherent in covering a large range of experimental frequencies, despite the frequent resort to the temperature shift function in viscoelasticity. The shortcomings of this numerical method lie in the fact that it cannot be used for viscoelastic materials with a weak transition region.
References Bagley, R.L., ‘A theoretical basis of the application of fractional calculus to viscoelasticity’, J. Rheol. 27(3), 1983, 201–210. Bagley, R.L. and Torvik, P.J., ‘Fractional calculus – A different approch to the analysis of viscoelastic damped structures’, AIAA J. 21(5), 1983, 741–748. Bagley, R.L. and Torvik, P.J., ‘Fractional calculus in the transient analysis of viscoelastically damped structures’, AIAA J. 23(6), 1985, 918–925. Bagley, R.L. and Torvik, P.J., ‘On the fractional calculus model of viscoelastic behavior’, J. Rheol. 30(1), 1986, 133–155. Beda, T. and Vinh, T., ‘Technique d’identification des paramètres rhéologiques du modèle à dérivées fractionnaires des matériaux viscoélastique (Identification technique of rheological parameters of viscoelastic materials)’, Mécanique Matériaux Electricité 441, 1991, 1–6 [in French]. Beda, T., ‘Modules complexes des matériaux Viscoélastiques par essais dynamiques sur tiges – Petites et grandes déformations (Dynamic testing of complex moduli of viscoelastic beams – Small and large strains)’, Ph.D. Thesis, CNAM, Paris, 1990 [in French]. Koeller, R.C., ‘Application of fractional calculus to the theory of viscoelasticity’, J. Appl. Mech. 51, 1984, 299–307. Nashif, A.H., Jones, D.I.G. and Henderson, J.P., Vibration Damping, John Wiley, New York, 1975. Rogers, L., ‘Operators and fractional derivates for viscoelastic constitutive equations’, J. Rheol. 27(4), 1983, 351–372.
118
T. BEDA AND Y. CHEVALIER
Rogers, L., ‘A new algorithm for interconversion of the mechanical properties of viscoelastic materials’, in Proceedings of the AIAA Dynamic Specialists Conference, Palm Springs, CA, May 17–18, 1984, Paper No. 84-1038 CP. Rogers, L., ‘An accurate temperature shift function and a new approach to modeling complex modulus’, in Proceedings of the 60th Shock and Vibration Symposium, Cavalier Hotel, Virginia Beach, VA, 14–16 November, 1989, 1–19. Soula, M., Vinh, T. and Chevalier, Y., ‘Transient responses of polymers and elastomers deduced from harmonic responses’, J. Sound Vibration 205(2), 1997, 185–203.
