Dynamic Preisach modelling of hysteresis for the piezoceramic actuator system

Dynamic Preisach modelling of hysteresis for the piezoceramic actuator system

Mechanism and Machine Theory 37 (2002) 75±89 www.elsevier.com/locate/mechmt Dynamic Preisach modelling of hysteresis for the piezoceramic actuator s...

268KB Sizes 0 Downloads 52 Views

Mechanism and Machine Theory 37 (2002) 75±89

www.elsevier.com/locate/mechmt

Dynamic Preisach modelling of hysteresis for the piezoceramic actuator system Yunhe Yu a, Zenngchu Xiao b, Nagi G. Naganathan a, Rao V. Dukkipati a

c,*

Mechanical, Industrial and Manufacturing Engineering Department, The University of Toledo, Toledo, OH 43606, USA b Mathematics Department, The University of Toledo, Toledo, OH 43606, USA c Department of Mechanical Engineering, School of Engineering, Fair®eld University, Fair®eld, CT 06430-5195, USA Received 4 April 2000; accepted 3 August 2001

Abstract Rate-dependent hysteresis property is a common phenomenon in various hysteretic systems including the piezoceramic material system. Dynamic Preisach model is needed to describe the rate-dependent hysteresis. This paper proposes a new dynamic Preisach model by introducing the dependence of the Preisach function on the input variation rate. An input variation rate function was introduced to adjust the relationship of hysteresis loop on the input variation rate for di€erent hysteresis systems. A detailed numerical implementation procedure is also presented. Experiments were conducted to study the hysteresis behavior of piezoceramic actuator system and to verify the proposed model. Ó 2002 Elsevier Science Ltd. All rights reserved.

1. Introduction In smart structures, piezoceramic materials can be used as actuators and as sensors. Among nonlinearities present in piezoceramic material systems, hysteresis has been particularly identi®ed to be sensitive to the varying ®eld conditions. Since piezoceramic materials are ferroelectric, they are fundamentally nonlinear in their response to an applied electric ®eld, exhibiting a hysteresis e€ect between the electric ®eld and the displacement or the force. Without modelling and incorporating hysteresis in the controller design, the hysteresis will act as an unmodelled phase lag presence, which has the potential to cause instability in a closed-loop system. Hence, reliable modelling and predictions of hysteresis would be a valuable tool when these piezoelectric actuators are employed as part of closed-loop system for purposes of motion control such as an active *

Corresponding author. Tel.: +1-203-254-4147; fax: +1-203-254-4013. E-mail address: [email protected]®eld.edu (R.V. Dukkipati).

0094-114X/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 4 - 1 1 4 X ( 0 1 ) 0 0 0 6 0 - X

76

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

control and micro-positioning. If hysteresis e€ects of these material systems can be predicted, then actuator controllers can be designed to correct these e€ects and the whole controller system can be made to appear as a device with a single valued output function. There is a vast literature on the hysteresis modelling using Preisach technique, especially in ferromagnetic materials. Recently, the Preisach technique has been applied to the piezoceramic material system. Sreeram and Naganathan [12,13] carried out the application of classic Preisach model (CPM) to a piezoceramic material system. In their work, Mayergozy's identi®cation and numerical simulation procedures were used to model the hysteresis exhibited by a piezoceramic bimorph. The minor loop patterns and enclosed areas from simulation and experimental data matched well for symmetrical minor hysteresis loop. Hughes and Wen [6] discussed and veri®ed the applicability of CPM to piezoceramic and shape memory alloy systems. The parameter identi®cation based on the input/output data and hysteresis compensation via direct inversion of the identi®ed model was also demonstrated. Freeman and Joshi [4] extended the CPM to a twoinput model to handle both the applied electric ®eld and the mechanical stress e€ects on the output strain. Ge [5] used both the CPM and generalized Preisach models to described the hysteresis of a piezoceramic actuator and developed a tracking control approach by incorporating the hysteresis model in the feed-forward loop to increase the tracking control precision. CPMs are basically static and rate-independent models and do not describe the dynamic hysteresis phenomena. Hysteresis with rate-dependent property has been observed in many experiments. Pokines et al. [10] and Ge [5] showed that hysteresis has rate-dependence properties for piezoceramic materials. Holman's experimental results show that the hysteresis curves become tilted if the frequency of the excitation is increased for a piezo material as shown in Fig. 1. It is considered that the CPM can be used to predict the hysteresis loop of a piezoceramic actuator under harmonic input excitation at frequency well below the resonance frequency of the actuator [5]. Experimental results showed that even in a low frequency range, the hysteresis loop changes with the

Fig. 1. Hysteresis curves at higher operating frequencies.

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

77

increase of input frequency. Bertotti's work [1±3] showed that the hysteresis loop area (the power loss) increases with the increase of magnetizing frequency for soft magnetic materials. Hysteresis curves of ferroelectric liquid crystal displays show that with the increase in frequency, an increase in width in the loop and a decrease in the height of the loop are noticed [11]. Dynamic Preisach models were proposed to deal with such a rate-dependent hysteresis. Two typical models are proposed by Mayergozy [8] and Bertotti [1]. The Preisach function of the dynamic term can be obtained by experiment through measuring the relaxation time. The detailed procedure of the numerical implementation was given by Mayergozy [9]. Bertotti's model is useful in predicting the power loss in a ferromagnetic material. However, it is not so e€ective in the analysis of the internal magnetic ®eld in the ferromagnetic device [14]. In this paper, a new dynamic Preisach model is proposed to deal with the rate-dependent property for piezoceramic materials by introducing the dependence of the Preisach function on the input variation rate. As a given parameter, the input variation rate is easy to measure and control during the experiments. The output variational speed of the system is dependent on the input variation rate and the hysteresis loop is directly related to the input variation rate. Hence, it is reasonable to use the system input variation rate as a parameter of the dynamic model. The details are discussed in the following sections. 2. Proposed dynamic Preisach model 2.1. Classic Preisach model The CPM of hysteresis can be represented in mathematical form as follows: Z Z l…a; b† cab u…t† da db; f …t† ˆ

…1†

aPb

where u…t† is the input and f …t† is the output of a system, cab is the elemental hysteresis operator. The output of each elemental hysteresis operator traces a rectangular loop on the input±output diagram switching from )1 to +1 when the input increases above the threshold a. The output switches from +1 to )1 when the input decreases below the value of b as shown in Fig. 2. l…a; b† is the distribution function, which is also called the Preisach function. The integral takes the past history values of u…t† to determine the current output f …t†. The past input history is accumulated in the a; b domain by storing the maximum and minimum values of the input variation. If the Preisach function l…a; b† is known, then the model (1) can be solved directly by integration. If the function l…a; b† is not known explicitly, experimental data are needed to estimate the Preisach function. However, the CPM represents only static and rate-independent hysteresis properties. Modi®cations to the CPM needed are to handle rate-dependent hysteresis phenomena. 2.2. Proposed dynamic Preisach model The dynamic models can be initiated by introducing the dependence of l-function on input variation rate du=dt as follows:

78

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

Fig. 2. Elemental hysteresis operator.

Z Z f …t† ˆ aPb

  du l a; b; c u…t† da db: dt ab

…2†

In the above integral, the time variable t is independent of the corresponding a and b values and so the integral is still a function of variables a and b. Practically and mathematically, the u-function in the above expression may be ill-behaved on the third variable du=dt, because the variation of du=dt can be great, but the output does not have so much variation. Hence the use of some pre-relationship between the l-function and du=dt can introduce more accurate mathematical description of the model. Therefore, we propose the following dynamic model:   Z Z  du l a; b; v …3† cab u…t† da db; f …t† ˆ dt aPb

where v…du=dt† is a function of the input variation rate, which describes the relationship between the input variation rate and hysteresis loop. The v…du=dt† function can be used to adjust the relationship of the hysteresis loop on the input variation rate for di€erent hysteresis systems. The choice of v-function is based on the experimental data. If the v…du=dt† function is so chosen that the amplitude of v…du=dt† has a small variation, then the power series expansion can be used for lfunction with respect to v…du=dt†. Hence      du du …4† ˆ l0 …a; b† ‡ v l1 …a; b† ‡    l a; b; v dt dt By retaining only the ®rst two terms of these expansions, we arrive at the following dynamic model: Z Z   Z Z du l0 …a; b† cab u…t† da db ‡ v cab u…t† da db: …5† f …t† ˆ l1 …a; b† dt aPb

aPb

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

79

It is clear that in the case of very slow input variations, the above models can be reduced to the corresponding static models. This means that the l0 -function in Eq. (5) should coincide with the l-function of static models in Eq. (1). Hence, we have Z Z   du v cab u…t† da db; …6† l1 …a; b† f …t† ˆ f…t† ‡ dt aPb

where the term f…t† stands for the static components of hysteresis. The second term in the Eq. (6) represents the dynamic term of the hysteresis. The static term in Eq. (6) has been discussed in detail by Mayergozy [9]. The purpose here is to determine the dynamic term using the experimental data. The corresponding numerical implementation method needs to be developed for applications. 2.3. Numerical implementation of the dynamic Preisach model Letting f~…t† ˆ f …t† f…t†, which is the di€erence between the dynamic and static outputs. Then Eq. (5) becomes Z Z   du ~ v …7† l1 …a; b† cab u…t† da db: f …t† ˆ dt aPb

In the limiting triangle, the input history is recorded as positive S ‡ …t† and negative S …t† sets that are separated by the staircase interface with vertices whose a and b coordinates are equal to …ak ; bk † as shown in Fig. 3. The Eq. (7) can be expressed as Z Z Z Z ~ f …t† ˆ l…a; b† da db l…a; b† da db: …8† ‡ Su…t†

Su…t†

Eq. (8) can also be written as Z Z   Z Z   du du 0 f …t† ˆ v v l1 …a; b† da db ‡ 2 l1 …a; b† da db: dt dt T



Fig. 3. Positive sets divided into n trapezoids Qk .

…9†

80

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

Since the integral of the Preisach function over the limiting triangle is equal to the output value when the input reaches the positive saturation value, we have Z Z   du ~ ~ v …10† f …t† ˆ f0 ‡ 2 l1 …a; b† da db; dt S‡

where the f~0 is the output of f~…t† when the input is equal to the positive saturation value. The positive sets S ‡ …t† can be subdivided into n trapezoids Qk with vertices whose a and b coordinates are equal to …ak ; bk † and …ak ; bk 1 † as shown in Fig. 3. Hence, we obtain  Z Z   n…t† Z Z  X du du v v l1 …a; b† da db ˆ l1 …a; b† da db: dt dt kˆ1 S‡

…11†

Qk

Note that the value of v…du=dt† depends on a and b for each Qk . If we divide each trapezoid Qk into nkj sub-trapezoids qkj as shown in Fig. 4, then Eq. (11) can be written as nkj Z Z   Z Z   n…t† X X du du v v l1 …a; b† da db ˆ l1 …a; b† da db: dt dt kˆ1 jˆ1 S‡

…12†

qkj

We can apply the Mean Theorem of Integration on each sub-trapezoid qkj , so that, Eq. (12) can be written as nkj   Z Z Z Z   n…t† X X du du v v l1 …a; b† da db; l1 …a; b† da db ˆ dt dt tˆnj kˆ1 jˆ1 S‡

…13†

qkj

where t ˆ fj corresponds to the point of a; b in qkj , which is determined by the Mean Theorem of Integration.

Fig. 4. Each trapezoid Qk divided into nkj sub-trapezoids qkj .

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

81

If the value of l1 -function is known, then the output of the system can be calculated by using Eq. (13). The l1 -function can be obtained from experiment as follows. In the experimental condition, Eq. (12) is also valid for the same sub-trapezoid qkj . Hence  e  Z Z Z Z  e du du v l1 …a; b† da db: …14† l1 …a; b† da db ˆ v dt dt tˆne j

qkj

qkj

Since qkj are narrow enough, the values of v…du=dt† on qkj are almost constant. Hence, we can replace v…due =dt†jtˆnej by the value at the mid-point or the end-point. In Eq. (14), v…due =dt† is a function of input variation speed. Each sub-trapezoid qkj can also be represented by the di€erence of the two triangles T …aj ; bj † and T …aj ; bj 1 † Z Z  e Z Z  e du du v v l1 …a; b† da db ˆ l1 …a; b† da db dt dt qkj

T …aj ;bj †

Z Z

‡ T …aj ;bj



 due v l1 …a; b† da db: dt

…15†



By substituting Eq. (14) into (15), we have 0 1 Z Z Z Z  e Z Z  e 1 du du B C l1 …a; b† da db ˆ due @ v v l1 …a; b† da db ‡ l1 …a; b† da dbA: dt dt v… dt † tˆne qkj

Tj

j

Tj

1

…16† In Eq. (16), the integration of the Preisach function over triangle T …aj ; bj † can be obtained from the experimental data and expressed as Z Z  e du v …17† l1 …a; b† da db ˆ fa0j fa0j bj ; dt T …aj ;bj †

where fa0j is the output value at the ®rst-order curve when the increase in input value is equal to aj , and fa0j bj is the output value when the input increases to aj and then decreases to bj as shown in Figs. 5 and 6. Eq. (16) can be simpli®ed as Z Z   1 l1 …a; b† da db ˆ due fa0j bj 1 fa0j bj : …18† v… dt † tˆn qkj

j

The expression for Eq. (13) can be written as

nkj Z Z   n…t† X  v…du † X du dt tˆnj v fa0j bj l1 …a; b† da db ˆ e du dt v… † e kˆ1 jˆ1 dt tˆn S‡

j

1

 fa0j bj :

…19†

82

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

Fig. 5. Triangular integration.

Fig. 6. Output value corresponding to triangular integration.

So the l-function is obtained by the experimental data. The output of the system can be calculated by Eq. (13).

3. Experimental study on the hysteresis behavior of a piezoceramic actuator system 3.1. Experimental set-up The purpose of this study is to ®nd the hysteresis behavior of a bimorph piezoceramic actuator at di€erent input conditions, so that it is possible to decide on how to model the hysteresis by the Preisach technique. To further understand the rate-dependent properties of hysteresis in the piezoceramic material system and to validate the proposed dynamic Preisach model, experiments were conducted on a piezoceramic actuator. The experimental set-up consists of various components such as actuators, input signal generating system, the output signal measuring system and the data acquisition system. The actuator is a commercial PZT bimorph with 0.019 in. thickness, 0.5 in. width, rigidly clamped using a ®xture at one end like a cantilever beam. The voltage input initiated by a signal generator was sent to the actuator through a high voltage ampli®er. The output displacement is measured by using an MTI Fotonic Sensor. The sensor probe is mounted

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

83

on a ®xture over the free end of the bimorph actuator. The measured displacement signal was converted into an electric signal and was sent to the PC. The graphical program LabVIEW made by National Instruments was used as the data acquisition system. Fig. 7 shows the schematic of the experimental set-up. 3.2. Experimental results Fig. 8 shows the response of the PZT bimorph actuator with sinusoidal input voltage for di€erent input frequencies. When the input frequency changes from 0.01 to 10 Hz, the input vs. output loops have similar patterns of change. Both the output maximum value and loop width decrease with the increase in the input frequency. The loops are normal hysteresis loops with input maximum corresponding to the output maximum value. In between the two extreme values there exist some phase lags. The details of change tendency of hysteresis loop at this input frequency range can be seen in Fig. 9. When the input frequency exceeds 10 Hz, both the displacement amplitude and the loop width increase. At a speci®c frequency, the displacement amplitude and width of the loop reach a maximum value. After that frequency, both the displacement amplitude and the loop width decrease with the increase in input frequency. The response behavior of the system in a relatively high input frequency range can be explained by the system inertia: the change of maximum value and width result from the amplitude and phase lag change associated with the system natural frequency. In this case, a simple one degree of freedom system dynamic equation can be used to explain this output phenomenon A X ˆ q ; ‰1 …f =fn †Š2 ‡ …2nf =fn †2 U ˆ tan

2n…f =fn †

1

1

…f =fn †2

;

…20†

…21†

where X is the displacement amplitude, U is the phase lag, A is the static displacement amplitude, f is the input frequency, fn is the natural frequency of the system and f is the damping factor. When

Fig. 7. Experimental set-up.

84

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

Fig. 8. Hysteresis loop at di€erent input frequencies.

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

85

Fig. 9. Comparison of hysteresis loop at relative lower input frequency range.

the input frequency is around the system natural frequency, the displacement will reach its maximum value and after that, the displacement amplitude will decrease with the increase in input frequency. Hence, the dynamic equations can be used to describe the input and output relationships in this relatively high input frequency range. However, the hysteresis behavior of the piezoceramic actuator at low input frequency range cannot be explained by the structure dynamic equations in this case. Holman et al. [7] explained this phenomenon as a creep e€ect. If a step voltage is applied to the piezo, then it will in the ®rst instance react almost immediately to the step, followed by a delayed reaction, which moves the extension of the piezo to the desired end value. If the frequency of the excitation voltage is low enough for the piezo to settle down, then there is no problem. However, if the frequency of the input voltage is increased, then it will go to a new position while the previous location has not yet been reached because of the delayed expansion. This will, therefore, result in lower amplitude compared to the lower frequency and will show up as a tilted hysteresis curve. In our case, both the amplitude and loop width decrease with the increase in the input frequency. Another interesting hysteresis behavior is shown in Fig. 10. The hysteresis loop shows the obvious di€erence when the input waveforms are sinusoidal and triangular. When the input waveform is sinusoidal, the hysteresis loop width and amplitude are bigger than those in the triangular input waveform case. This can also be explained by Holman's explanation. This may be better stated in terms of the Fourier analysis of the two waveforms, i.e. the triangular waveform has harmonics of higher frequency, which are e€ectively attenuated by the system. It gives the piezo actuator more expansion time to settle down. Hence the amplitude is higher.

4. Hysteresis modelling and prediction The proposed dynamic Preisach model is used to model and predict the output of the system in the relatively low frequency range. As discussed previously, the v…du=dt† function describes the

86

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

Fig. 10. Comparison of hysteresis loop at di€erent input waveforms.

Fig. 11. Comparison of measured and predicted hysteresis loop at input frequency equal 0.5 Hz: (a) static Preisach prediction; (b) dynamic Preisach prediction.

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

87

relationship between the input variation speed and the change tendency of the hysteresis loop. Normally the form of v…du=dt† function varies in the Preisach domain (a±b half plane) and can be determined by the experimental results. Our experimental results and Holman's explanation show that both the output maximum value and loop width decrease with the increase of the input frequency. This means that v…du=dt† function can be constant in Preisach domain in this case. According to the experimental results shown in Fig. 8, the v…du=dt† function can be expressed in simple form to describe the amplitude di€erence of static and dynamic hysteresis loop as    n du du c2 ; …22† ˆ c1 v dt dt where c1 , c2 and n are the constant coecients and equal to 1.37, 0.286 and 0.39, respectively. Now the proposed dynamic Preisach model can be used to model and predict the hysteresis of PZT bimorph actuator system. In implementing the proposed dynamic Preisach model, a static and a dynamic ®rst-order curve need to be created. Input frequency of 0.01 Hz was considered as the static state since experiments

Fig. 12. Comparison of measured and predicted hysteresis loop at input frequency equal to 5 Hz: (a) static Preisach prediction; (b) dynamic Preisach prediction.

88

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

showed that there is almost no di€erence in the hysteresis loop when the input frequency is lower than 0.01 Hz. The hysteresis loop at an input frequency of 1 Hz was used for the dynamic hysteresis loop. The proposed dynamic Preisach model predicts the hysteresis behavior when the input frequency changes from 0.01 to 10 Hz in this case. The veri®cation of the proposed dynamic Preisach model can be inferred by its ability to predict the higher-order transition curves. To do so, two high-order transition curves were measured at input frequencies equal to 0.5 and 5 Hz, respectively. Then, by using the developed computer codes that implement the proposed dynamic Preisach model, the higher-order transition curves predicted by the model were computed for the same sequence of ®eld reversals as in the experiments. The comparison of the predicted higher-order transition curves with the experimentally measured higher-order transition curves is demonstrated in Figs. 11 and 12. Figs. 11(a) and 12(a) show the classic Preisach prediction results. Since the CPM cannot account for the rate-dependent e€ect, the higher the input frequency, the more the prediction error. The comparison of results also shows that the dynamic predicted results agree well with the experimental results. This veri®es that the proposed dynamic Preisach model is e€ective for the modelling of the hysteresis of the piezoceramic actuator system.

5. Conclusion Experiments showed that the dynamic behavior of the piezoceramic actuator system is di€erent at di€erent input frequency ranges. In a relatively low frequency range, the system behavior cannot be described by using dynamic equations. The proposed new dynamic Preisach model is used to model and predict the output of the piezoceramic actuator system in this range. The new dynamic introduces the dependence of the Preisach function on the input variation rate. When a function of input variation rate is also introduced to adjust the relationship of the hysteresis loop on the input variation rate, the proposed model can be used for di€erent hysteresis systems. The results show that the predicted results on the higher-order transition curves agree well with the experimental ones. This veri®es that the proposed dynamic Preisach model is e€ective for the hysteresis modelling of the piezoceramic actuator system.

References [1] G. Bertotti, Dynamic generalization of the scalar Preisach model of hysteresis, IEEE Transactions on Magnetics 28 (5) (1992) 2599±2601. [2] G. Bertotti, F. Fiorillo, M. Pasquale, Measurement and perdition of dynamic loop shapes and power losses in magnetic materials, IEEE Transactions on Magnetics 29 (6) (1993). [3] G. Bertotti, F. Fiorillo, M. Pasquale, Dynamic Preisach model interpretation of power losses in rapidly quenched 6.5% SiFe, IEEE Transactions on Magnetics 30 (6) (1994). [4] A.R. Freeman, S.P. Joshi, Numerical modelling of PZT nonlinear electromechanical behavior, SPIE Conference Proceeding 2715 (1995) 602±613. [5] P. Ge, Modelling and control of hysteresis in piezoceramic actuator, Ph.D. Dissertation, University of Rhode Island, 1996.

Y. Yu et al. / Mechanism and Machine Theory 37 (2002) 75±89

89

[6] D. Hughes, J.T. Wen, Preisach modelling of piezoceramic and shape memory alloy hysteresis, SPIE Proceeding 2715 (1995) 507±528. [7] A.E. Holman, P.M.L.O. Scholte, W.Chr. Heerens, F. Tuinstra, Analysis of piezo actuators in translation constructions, Review of Scienti®c Instruments 66 (5) (1995) 3208±3215. [8] I.D. Mayergozy, Dynamic Preisach models of hysteresis, IEEE Transactions on Magnetics 24 (6) (1988) 2925± 2927. [9] I.D. Mayergozy, in: Mathematical Modelling of Hysteresis, Springer, New York, 1991, p. 1991. [10] B. Pokines, W.K. Belvin, D.J. Inman, Static and dynamic characteristics of a piezoceramic strut, in: Proceedings of Fifth NASA/DoD Controls±Structures Interaction Technology Conference, Lake Tahoe, Nevada, March 3±5, 1992. [11] C. Reynerts, D.V. Alexis, Polarization against voltage hysteresis loop of surface-stabilized ferroelectric liquid crystals, Journal of Physics D 22 (1989) 1504±1509. [12] G. Salvady, P.N. Sreeram, G. Naganathan, Hysteresis model for piezoceramic actuator systems using Preisach models, in: Proceedings of the North American Conference on Smart Structures and Materials, Orlando, Florida, February 13±18, 1994. [13] P.N. Sreeram, N.G. Naganathan, Hysteresis prediction for a piezoceramic material system, in: Proceeding of 1993 ASME Winter Anneal Meeting, New Orleans, LA, Vol. AD-Vol33, 1993, pp. 35±42. [14] F. Vajda, E.D. Torre, Relationship between the moving and the product Preisach models, IEEE Transactions on Magnetics 27 (5) (1991) 3823±3826.