Steady-state vibrational response of high-speed flexible mechanisms

MechamsmandMac~neTheoryVol 19,No 4/5,pp 417423, 1984 Pnntedm the US A

0094-114](/84 $300+ 00 PergamonPressLtxl

STEADY-STATE VIBRATIONAL RESPONSE OF HIGH-SPEED FLEXIBLE MECHANISMS W. L. CLEGHORN~" Department of Mechanical Engineering, Umversity of Mamtoba, Winnipeg, Manitoba, Canada R3T 2N2 and

R. G. F E N T O N , and B. T A B A R R O K , Department of Mechanical Engineering, Universtty of Toronto, Toronto, Ontario, Canada M5S 1A4 (Received 28 January 1983) Abstract--A procedure is developed for determtmng the steady-state solution for the govermng finite element equations of high-specd flexible mechanisms, operating with constant input rotational speed. The method revolves representing the periodically varying components of the govermng equations by truncated Fourier series. By matching harmonic terms, a system of linear equataons is obtained and subsequently solved for the harmonic coefficaentsof the response. The procedure is tested by comparing the calculated responses with prewously published expenmental results.

1. INTRODUCTION Unlike linear vibrating systems where small oscillations are considered about a fixed equilibrium positaon, models which take into account flexibihty of a moving mechamsm have the problem of deahng wath time dependent sUffnesses and masses to model the system One method of determining the steady-state wbrataonal response under these condittons would be by imposing arbitrary initial conditions, and solving the governing equations in time until the transient solution dies out and only the steady-state solution would remain. However, other techniques have been developed whereby the steady-state solution is found directly. Sadler and Sandor[1] determined the steadystate vibrational response of a point mass located on a flexible coupler of an otherwise massless, ngld 4-bar mechanism Bahgat and Willmert[2, 3] used a, harmonic series techmque for mechamsms wath all links hawng distributed mass and flexibility. The governing equations were formulated using the fimte element method, and took the followang form [M[{~} + [B]{~} + [K]{4~} = {P}

(1)

where [M], [B] and [K] are the global mass, damping, and stiffness matrices, respectavely, and {~b} and {P} are the global deflection and load vectors. The procedure revolved the representation of periodically varying elastic deflections and the load vector by truncated Founer senes, and thereby expressions were derived for the harmomc components of the global deflection vector in terms of global load vector harmonic coeflicaents and the instantaneous values of the global matrices. The variations of elastic deflections tAsmstant Professor ~Professor 417

throughout a cycle are found by calculatang the global matrices at several locations, and substituting into the derived expressions. A numerical problem arises using tlus technique. The denved expressions contain the term [[K] - j 2072[M]]- 1

(2)

where j is an integer employed as an index of summation in the truncated Fourier series, and 07 ts the constant input rotational speed. As the mechamsm moves, the global matrices change. It is possble that one or more of the elgenvalues of the problem [KI{~} = 2[MI{O} will attain the value ofj207 2 somewhere m the cycle. When this happens, expression (2) becomes mdeterrmnat¢ and there is no solution for the deflections. This problem can occur for any mechamsm. For mechanisms with stiffer hnks, higher eigenvalues m the above equation are obtained In this instance, the equahty j2072 = ~. can stall occur, but at larger./values. Another procedure of determining the steady-state solution for a mechanism in which all links are considered flexible was presented by Midha[4]. A cycle of motion was divided into a large number of equal time intervals. For each ttme interval, a system of differential equations was generated by assembling mass and stiffness matrices and a load vector assocaated with the corresponding instantaneous structure. The matrices were assumed to be constant throughout the time interval To ensure a steady-state solution,

418

W. L. CLEGHORN et al.

additional equaUons were generated by matching deflections and the first time derivative of each deflection at the beginning of each cycle wth those at the end of the previous cycle This tecbmque is statable only m instances where the differential equations can be decoupled. Nath and Ghosh[5] employed a very specmhzed harmonic analysis in which the truncated Fourier series were introduced at the element level early m the analysis The element equations were assembled by using a set of dependent global coordinates which reqmred an additional system of constraints to be satisfied. An elimination of these constraints m the end, has led to a non-banded matrix in the system of equattons. This shortconung is avoided m the present analysts which ymlds a banded matrix Thus, much larger systems can be solved with the same computing effort. Recently, Wiederrich[6] examined the functions of kinetic energy, potenUal energy and energy disstpataon for a mechanism with 1 degree of freedom. From these functions were determined relataonships regarding the operating response of the system The procedure presented below for determining the steady-state solution for flerable mechanisms ~s formulated w~th finite elements. It combines the advantages of representing all equation components by continuously varying functions, while at the same tame employing procedures of element assembly simtlar to those cases of an instantaneous structure formulatton In addition to the properties required for performing a vthrataonal analysis of an equivalent stataonary structure, the only reformation required ts a rigad body dynamm analysts of the mechanism. A s~mflar algorithm for determimng the periodic steadystate solution to that outlined in this paper was covered by Bolotin [7] for the Mathieu-Hfll equataon.

any function f ( ~ ) d

dt (f(cz)) = co

df=ogf,

d -d-ti ( f (0~)) = co2f,,.

(4)

(5)

Applying the above transformation to the time dependent vectors in eqn (1), and implementing eqn (3), gives

o92[M]{(;b"} + eJ[Bl{~b'} + [[Ks] + [KD]I{~b} = {P}. (6) Some of the matrices as well as the global load vector m eqn (6) are dependent on o9 as follows [BI = co[~1

(7)

[Ko] = o9~[£~1

(8)

{P} = ca2{ff}

(9)

Circumflexed quanUUes in eqns (7)-(9) have all c~ dependence removed. Substituting eqns (7)-(9) in eqn (6), gives

[M]{~"} --I-[J~]{~'} + []~{¢ } = {if}

(lO)

[£] = ~ [Ks] + [£.].

(11)

where

In order to determine the steady-state solution, global matrices and the global load vector ofeqn (10) are represented by Fourier series of order N1, and the global deflection vector is represented by a Fourier series of order N2 as follows

2. MATHEMATICAL FORMULATION NI

2.1 Solutton of governing equattons The governing equataons employed in the analysts have been derived in[8] and take the form ofeqn (1). It has been assumed that there are no bearing clearances, and there are no frictaonal losses. Components of matrix [B] arise solely due to gyroscopic accelerations, and no physical dampmg is added to the system The stiffness matrix is considered to be the sum of two components, as follows

[ M ] = ~ ([~]cosmo~ + [ ~ s i n m ~ )

(12)

m--0 NI

[~l = ~' ( [ ~ ] c o s m a +[~=]sinma)

(13)

m~0 NI

[/~]= ~ ([/~]cosm~ + [ ~ ] s m m ¢ )

(14)

mffi0 NI

[KI = [Ks] + [g~,].

(3)

Matrix [Ks], the structural component of the stiffness matrix, is equivalent to the stiffness of a corresponding instantaneous structure. Matnx [KD], the dynamic component of the stiffness matrix, gwes changes of the global stiffness matrix as a result of dynamic effects. For the case of constant input rotational speed, tt ~s expedient to use the input angular displacement, a, as the independent variable, instead of time, t. For

{P} = ~ ({/~} cosja + {~} sin j00

(15)

1=0

N2

{q~} = ~' ({q~.) cos M + {~.} smng).

(16)

~E0

Coefficients in eqns (12)-(15) are found by the method presented in the Appendix. The method requires that the mechanism be positioned at 2NI + I equally spaced locatlons of the input I n k throughout a cycle and the global matnces and the global load vector be calculated.

419

Steady-state vibrational response

and the axial deflection, u~, by a constant. The flexural stress distributaon at the outer fiber of an element is

Substituting eqns (12)-(16) m eqn (10) gives a system of NT linear equations in terms of the unknown deflection vector coefficient components, where

Eh ~2uyE

NT = (2N2 + I)N~

(21)

o = ---2 c~x2

(17)

where E is the modulus of elasticity. Relating the nodal deflections to the transverse deflection distribution and substituting in (21), the stress distribution may be put in the form

and N u is the number of elastic degrees of freedom associated with a mechanism. The equations, derived in detail in Ref. [9] for the case that N2 > 2 N . can be put in the following matrix form

@

[AT]{~br} = {Pr}

(18)

o = a0 + olx + o2x 2 + o3x 3

(22)

Eh O'o= -~- ml

(23)

where

where

{,T} = ({S0} T, {$,}T, {~,V, {$~}T, x {~2} T, .... {$2} T, {~n2}7] T {PT} -~ [{~I$0}T,{/61}T, {P1}T, "" , {JSNt}T,

(19)

× {~,V' {0V, .... {0}T.

(20)

at = Eh

- ~ (v, - v9 -

(301 + 202)

_ 3 (3ml - m,)~

If the equations are arranged such that they are based on harmonic comparisons in the following order

2/.

(24)

_1

1, cos % sin a, cos 2% sin 2 % . . . , cos N2ot, sin N : , then matrix [At] will be banded, unsymmetric, with a band-width of 2(Nt + I)NM. The substitution of eqns (12)-(16) in eqn (10) gives harmonics in the range from N2 + 1 to N I + N2 on the l h.s which have no equivalent terms on the r.h.s These terms are all assocaated only with relatively high harmonic orders. Because the harmonic coefficients for both the equation components and solution converge to zero when sufficmntly large expansions are considered, the effect of unmatched terms is small, and, therefore, can be ignored when calculating the solution. 2.2 Stress and strata analysis Once the global variables, {~}, have been determined, nodal deflections for an element at any partacular time can be found[9]. For the examples presented below, the links were represented by urnform beam elements The nodal variables are shown in Fig 1 The transverse deflection distribution of an element, uye, is represented by a qmntic polynormal,

C

roIolron

U

l

o3 = Eh

-

(v~ - v,.) - ~ (01 + 0,)

L3 (ml -- m2) .

(26)

Having determined the normal stress at the desired locataons, the corresponding normal strain at the outer fiber of an element is given by

= ale

(27)

where the longitudinal strain has been neglected. 3. C O M P A R I S O N

WITH

EXPERIMENTAL

RESULTS

A number of researchers have used Alexander's experimental results [10] of planar 4-bar angnlar function generating mechanisms for comparisons with the

~

~.vm 2

~curvoture I

IP

x-O x=L x Fxg 1. Nodal variables of a typical beam fimte element.

420

W L CL~SHORNet al

results of their analyacal solutions. In the experimental setup, the polar moment of inertia on the shaft which drives the input rink was assumed large enough to prox~de a constant input rotational speed. All hnks were made of aluminum. Experimental results were obtmned by placing strain gauges on the upper and lower surfaces of the coupler and output hnk at various locat, ons, mchidmg rmdway along their lengths. Links of various lengths having rectangular crosssections were consMered. One mechanism m particular has been used by researchers for comparison The properaes of this mechanmm, referred to in this paper as the experimental mechanism, are hsted in Table l_ As recommended by Alexander, the mass densities of the coupler and output links have been increased by 8% to account for the additional mass of the wires and strain gauges attached to these links in the experimental setup In addRlon to the properties listed in this table, beanng masses at pins joining the input link and coupler, as well as the coupler and output link, existed in the experimental set up. The addmonal mass at each end of the two pros was 0 06 Ibm. Unless otherwzse stated in this paper, these masses have been neglected in the mathematmal model Results were presented by plotting steady-state strain measurements throughout a cycle of motion For these plots, the deflected link shapes which give negative strains are illustrated in Fig 2 Angular displacements of the input link, 0~, are measured posmve clockwise from the horizontal line joimng the two stationary bearings Alexander, by companng strain measurements on the upper and lower link surfaces, concluded that longitudinal strains are small as compared to bending strains. In his experiments, longitudinal strains were generally in the order of 5 ~ of bending strains Alexander[10] and Mldha[4], compared their analytical results with these experimental results In both cases, the governing equaaons were based on the instantaneous structure formulation[11] Each mech-

'-4'7

/

Fig 2 Global deflectmn variables of a 4-bar mechamsm. anlsm link was modeled w~th one fimte element which based transverse deflections on a cubic polynomial. Sadler[12], used a number of lumpedmasses along each mechanism link Each lumped mass had its transverse deflection measured from the undeflected state of the mechanism link on which it was located. Results presented in Figs. 3 and 4 employed the flexible mechanism model shown in Fig. 2, assembled from beam finite elements (Fig. 1), using the soluUon

,o

- -

AnalyllC~fl

----

El[perlmenlal

50 --

0

9OJ,

kl

ISO~

~

~27D

~O

o: -5

0

-~o Fig 3. Analyttcal and experimental coupler midpoint strains neglecting lumped bearing masses in analyacal model (ca* = 41.88 rad/sec)

Table 1. Properties of experimental mechamsm

Lxnk

Modulus

of

Elastlclty, (Ib/mn 2 )

Mass E

Lxnk

Denslty,

p

Height,

( I b - s e c Z / i n ~)

Link h

Width,

Llnk b

Length,

(In)

(in)

(in)

Input

i0.

x 106

2.52

x i0 -~

0.157

1.00

4.00

CoupleP

I0.

x 106

2.73

x 10 -4

0.063

1.00

II.0

Output

i0.

x 106

2.73

x I0 -~

0.063

1.00

10.5

length

of base

llnk

= 10.0

in

L

421

Steady-state vibrattonal response

of

Pi

--

Anolyhcal

----

Experirnenlol

5O

I

I

--

i,

0,

,'

I

I

t~

\/

-50

'l

,' ,~ ',,., /"--'\~,o'~, ~do/ ~ ,,-',\ / 0z

91o ,,.,,i

,

/

"" ",,,X.J \

,

-I0

%xJ

Rg. 4 Analytncal and experimental output hnk nudpomt strmns neglecting lumped beanng masses in analytical model (co* = 41 88 rad/sec).

procedure outlined in Section 2.1. Also shown in these figures are the expenmental results obtained by Alexander. One reason for the poor comparison between the analytical and experimental results is that lumped bearing masses at the pins joining the mowng hnks were neglected For the analytical results shown in Figs. 5 and 6, the effect of lumped masses has been included, but only in a bruited sense. Only the energy related to the translations of the lumped masses was considered, and the rotational energy was neglected. In order to include the rotational energy of the bearing masses, it would be necessary to introduce additional global variables. These would be curva-

tures between the lumped masses and the ends of the links. Such a model was not pursued• One portion of the plots comparing analytical and experimental midpoint strains of special interest zs that related to the output link in the vicinity of 07 = 340 °. Large experimental midpoint strains were found to occur dunng tlus portaon of the cycle (see Fig. 4). Th~s has been referred to by Alexander as an inertia pulse. Researchers have experienced difficulty reproduong the amphtude of tl'us pulse using thezr analytical models. One reason for tlus is the number of simphficatzons made when denying the governing equations. For instance, Figs. 7 and 8 show midpoint strains calculated by using the governing equataons based on the instantaneous structure formulation which corresponds to the case in eqns (1) and (3) of [BI = [0l, [K.I = [0].

Also shown in the same figures are the results obtained using the govermng equations derived in[8]. Lumped beanng masses were not used m either case Both cases used the same set of global variables (Fig. 2), the same input runmng speed (co* =41.88rad/sec), the same solution techmque (Section 2 1), and the same orders of truncated Fourier series expanszons (N] = 11, N2 = 23). A companson of the results in Fig. 8 indicates that a better

--

Employing Equahons of Reference 8

---

Inslonfoneous Slruclure Formulohon

I0 ]I0

--

Analyhcal

---

ERperumental

50

SO

,,",

.,

-.

\

o w

'-,?\

w

,y

~

',W V

'

~

,8ox

/~7o

.... ~'~

//

~'-'..?~..

,

~-

~90 ~ - - ' . ~ . -'---180

270 "~

,

360

oz

sto

k../',w

o;

-.5

0

-.50 -I0 -IO

Fig 5 Analytical and expenmental coupler m~dpolnt stratus including lumped bearing masses in analyUcal model (m* = 41.88 rad/sec)

I0.~

.50

;~

,,

--

Anolyhcal

---

~'eperlrnlnlal

Fig 7 Coupler midpoint strains of experimental mechamsm, comparison of two analytical formulations (m* ffi41 88 tad/see)

/~

-.50

-IO

Employ,ng Equohons of Reference 8

-----

Inslontaneous SlrueTure Formulation

50

b " v

--

--

"-'

/

" """"._~',,\

0

02

•

SO

.

180

~_-..

270

1

360

-50

\

i \./

Fig. 6_ Analytical and experimental output link midpoint strmns including lumped bearing masses in analyucal model ( m * = 41.88 rad/sec)

-I0

Fzg 8. Output hnk rmdpomt strains of expenmental mechamsm, comparison of two analyUcal formulations (aJ* = 41.88 rad/sec).

422

W. L CI,mmOgN et al.

representation o f the experimental merUa pulse 1S

achieved using the equations derived in[8]. The reason for this is that there are compressive longitudinal loads m the output hnk dunng the mertaa pulse When using the equations derived in [8], compressive Iongctudinal loads reduce flexural stiffnesses associated with the transverse direction. This leads to larger flexural strain amplitudes. However, when using the instantaneous structure formulation, the effect of longitudinal loads is neglected. Figure 9(a) shows the diagonal stiffness matrix component related to global variable 6 (see Fig. 2), using both the instantaneous structure formulation and the equations derived in[8]. The midpoint longitudinal load of the output link is shown in Fig. 9(b). For all mechamsm models considered in this paper, one finite element was used to represent each flexable link. There is of course nn principle no reason why more than one element could not be used to model such systems It is anticipated that with a multi-element idealization, calculated parameters, including strain, would improve. Some previous investigators[13-15] have compared results using various numbers of fimte elements to model a flexible mechanism. For all calculated results presented, the sets of Lmear equations were solved using subroutine L E Q T I B obtained from the I M S L library of mathematical and statistical subroutine[16] The C P U time required using an IBM 3033 computer to generate the equations and detemune their solution was approximately 0.85 mm. 4. SUMMARY

A procedure has been presented whereby the steady-state solution to the equations of a moving and vibrating mechanism is found directly. The govermng equations developed in[8] and those based on

--

g 2 oo

K66

- - - -

1 0

KS,66 I 180

t 270

I 360

(o)

~° I

90

180

270~"N,,~60

1 J_ P Sadler and G N Sandor, IOneto-clastodynamlc harmomc analys~s of four-bar path generating mechanisms A S M E Paper No. 70-Meeh-61, presented at the llth ,4SME B~enn. Mech Conf., Columbus, Ohio (1970) 2 B M Bahgat, General Finite Element Vibraoonal Analysis of Planar Mechanisms Ph.D. Dissertation, Clarkson College of Technology (1973). 3. B M. Bahgat and K_ D. Willmert, Fimte element vibrational analysis of planar mechanisms. Mechan~m and Machine Theory 11(1), 47-71 (1976). 4. A Midha, A. G Erdman and D. A Frohrib, A closed form numerical algorithm for the periodic response of high-speed elastic linkages Trans. ASME, J. Mech_ Design 101, 154-162 (1979). 5. R. K Nath and A. Ghosh, Steady state response of mechamsms vath elastic hnks by finite element methods Mechamsm and Machine Theory 15(3), 199-211 (1980)

6 J. L. Wiederrich, A theory for the identification of a single degree of freedom machine from ,ts pencdic operating characteristics. A S M E Paper No. 82-DET-i. 7 V V. Bolotin, The Dynamic Stability of Elastic Systems (in Russian), Gostekhizdat, Moscow, 1956 (Translation by V I Weingarten, et al.). Holden-Day, New York (1964). 8 W. L Cleghorn, R. G. Fenton and B Tabarrok, Flmt¢ element analysis of high-speed flexthle mechanisms Mechan~m and Machine Theory 16(4), 407-424 (1981). 9 W. L. Cleghorn, Analysis and Design of Htgh-Speed Flex:ble Mechana'ms Ph D. Dissertation, University of Toronto (1980). 10. R M. Alexander, An Experimental and Analytical lnvesttgatwn of the Dynamic Response of an Elasttc Four-Bar Mechanism. Ph.D Dissertation, University of Texas at Arlington (1975). I 1 I Imam, G. N. Sandor and S N Kramer, Dcflecuon and stress analyms m high-speed planar mechanisms with elastic links Trans ASME, J Engng Ind 95B, 541-548 (1973). 12 J P Sadler, On the analytical lumped-mass model of an elastm four-bar mechanism. Trans. ASME, J Engng Ind_ 97B, 561-565 (1975) 13 A Midha, M. L. Badlam and A G. Erdma,, A note on the effects of multi-clement idealiTation of planarelastic-linkage members Proc 5th OSU Appl. Mech. Conf., Stlllwater, Oklahoma, 27 1-27.11 (1977). 14 D. A Turcic, A General Approach to the Dynam:c Analysts of Elastic Mechamsm Systems. Doctoral Dissertation, The Pennsylvania State Umverslty (1982) 15 W L Cleghorn and C J Konzelman, Comparative analysis of fimte element types used In flexible mechanlsrn models Proc. 8th OSU Appl. Mech. Conf , St Louis, Missouri (1983)

,

02 ~5o S ~

REFERENCES

6th Edn, Vol. II, Houston, Texas (1977).

J

~

the instantaneous structure formulataon were employed with this procedure. A better representation o f experimental results was obtained using the equations developed in Ref. [8]_

(b)

Fig. 9 Global stLffness mamx components K~ and Ks,~ and output link rmdpomt Ionlctedlnal load of experimental mechanism (o~*=41 88 tad/see), (a) Global i g s ~ components. (b) Output link rmdpomt longitudinal load.

APPENDIX

Representatwn of a perwdtc function by a truncated Fourier series The problem of describing penochcaUy varying fuactaons m flus paper suggests the use o f a Fourier series, and the practical inev~tabifity of trunmung the series implies that the denved expression will be approxtmate. However, the ap-

423

Steady-state vibrational response proxtmation can be made as close as desired, recognmng that the senes is complete and convergent The task of determJmng the coeffim~ of a truncated Fourier senes can be simplified considerably by requmng that the truncated senes represents the penod,c fancuon exactly at 2N + 1 equally spaced points As shown in the following, each of the coeffacients of the truncated Fourier senes is obtained directly by a stmple summation, as opposed to techniques which involve numerical integratmn or the solution of sets

where

21tk 2N+I'

'%

(k=0,.

,2N).

(A3)

The method of solwng these equatmns is given in[9], and the result is 1

2M

of equations. 2

~=2N+lk.

Determination o f Fourier series coefficients Consider a periodically varying function F(u) having a penod of 2~t. The function may be represented by a truncated Fourier series as follows N

(A1)

F(a)= E (~cosJ" + ~slnja) J-O

In order to determine the 2 N + 1 coeificmnts ( ~ , j = 0, _ , N, ~ , j = 1. . . . N), the function is evaluated at 2 N + 1 equally spaced positions throughout ,ts period and is then equated to the truncated Fourier series at the same location. That is N

F(u,)= E (~cosj=,+~sinj~,),

REGIME VIBRATOIRE

W.L Cleghron,

R~sum~ bles

-

R.G.

PERMANENT

aux ~ l ~ m e n t s

DES M E C A N I S M E S

,N) (A5)

2/V+I

2,!

2~jk

(J = 1,

, N).

FLEXIBLES

ces m o u v e m e n t s .

varlant p~rlodlquement

qu~es.

L'assocxatxon

n~axres

qul sont r ~ s o l u e s pour

FONCTIONNANT

A TRES G R A N D E V I T E S S E

le r ~ g l m e p e r m a n e n t des m ~ c a n l s m e s

La m ~ t h o d e d ~ v e l o p p ~ e

condult

les c o e f f x c x e n t s h a r m o n x q u e s

les r ~ p o n s e s

flexl-

~ p a r t l r des ~ q u a t l o n s conslste

A representer

dans ces ~ q u a t l o n s par des s~rles de F o u r l e r tron-

des t e r m e s de m@me h a r m o n l q u e

est t e s t ~ e en c o m p a r a n t

(A6)

Using eqns (A3)--(A6), components of F ( - ) are determined exactly by the truncated Founer series at 2N + 1 values of a spe~fied by eqn (A3) Between these ~ values, components of F(-) computed from the truncated Fourier senes are only approxanate. The accuracy with which a truncated Fourier series of order N represents a function between the 2N + 1 values of % depends on the nature of F(-). The example problem presented m thin paper reqmred that an expanmon of 1lth order be used so that the absolute maximum difference throughout a cycle was brought within one percent of its maximum absolute value, for all of the equation components

~ v l t e s s e de r o t a t l o n d ' e n t r ~ e c o n s t a n t e ,

f~nls r ~ g l s s a n t

ant~rleurement.

(A2)

21tjk

o

= 2 N + 1 ,~.o F(o%) sin 2 - - f f - ~ '

xndique comment d~termxner

fonctlonnant

les p a r a m ~ t r e s

,2N)

...

L x't%) cos _..--=----7, (J = I,

F e n t o n et B. T a b a r r o k

Cet a r t x c l e

rapldes

(k =0,

2

~

axns~ c a l c u l ~ e s

~ un s y s t ~ m e d ' ~ q u a t i o n s de la r~ponse.

~ des r ~ s u l t a t s

ll-

La m ~ t h o d e

exp~rlmentaux

publ1~s