EXPERIMENTAL DETERMINATION OF VIBRATION ENERGY AND SEA COUPLING LOSS FACTORS OF A T-SHAPED BEAM

EXPERIMENTAL DETERMINATION OF VIBRATION ENERGY AND SEA COUPLING LOSS FACTORS OF A T-SHAPED BEAM Edmilson R. O. Santos1 Khaled M. Ahmida2 José Roberto F. Arruda1 José Maria C. Dos Santos1 1Departamento de Mecânica Computacional, FEM, C.P. 6122, Campinas - SP CEP 13083-970, [email protected] 2Centro de Componentes Semicondutores, CCS, C.P. 6061, Campinas - SP CEP 13083-970, [email protected] Abstract. Coupling loss factors are of essential importance in Statistical Energy Analysis (SEA). These factors characterize the dynamic coupling between subsystems and are required to predict vibration and noise levels of structures and acoustical cavities. In this paper, the coupling loss factors of a round-robin T-beam structure are obtained using the Power Injection Method (PIM), which corresponds to the inverse SEA problem. In the estimation of the coupling loss factors by PIM, the subsystem energies were obtained from response measurements. The T-beam was subdivided into six subsystems, three representing the energy of flexural waves and the three others representing energy of longitudinal waves. The subsystem energies were calculated from equally-spaced acceleration measurements in the transverse and longitudinal directions along the beams. Two different techniques were used. The first consists of using the lumped mass assumption and the hypothesis that the total energy is twice the kinetic energy. In the second, the displacement field is interpolated via the Spectral Element Method (SEM) and the kinetic and potential energies are obtained by integration of the analytical expressions. Keywords. Coupling loss factors, Statistical Energy Analysis, Power Injection Method, Subsystem Energies, Spectral Element Method. 1. Introduction The interest in identifying the vibration transmission paths in dynamic systems comes from the needed to predict structural vibration levels at high modal densities. SEA is an approach to vibration prediction at high frequencies developed in the late 1960s by Lyon, Smith, Maidanik and others. Such approach estimates the distribution of vibrational energy in a structure, which is represented in terms of a set of connected subsystems (group of “similar” energy storage modes). Each subsystem has a quantity of stored vibrational energy that is related with the input power through some parameters known as coupling loss factors (CLFs) and internal loss factors. If the coupling loss factors and internal loss factors are not known or are only partly known or are to be verified for the purpose of confirmation of some prediction technique, then it may be desirable to measure them in a test structure (Bies and Hamid, 1980). One way to do it is to measure energy levels in many points of the structure due to different power inputs. Then, for a large number of input power applied to structure, the resulting linear system of equations may be inverted to determine the coupling loss and internal loss factors. Bies and Hamid (1980) seem to have been the first to use approach that is currently called power injection method (PIM). 2. Review of SEA One of main methods for vibro-acoustic analysis of coupled structures in the high frequency range is SEA. This method aims to predict the vibration energy stored in different regions of the structure, that is, the state of vibration is expressed in terms of stored, dissipated and transmitted energies between substructures (or subsystems). Every single subsystem can be modeled as a set of modes. For example, in the case of beams, the longitudinal modes in a frequency band can be considered as a subsystem, and the flexural modes as another subsystems. Longitudinal and flexural behavior should be considered separately (De Langle, 1996). Then, by the principle of conservation of energy, a power balance matrix equation for the connected subsystems can be derived (Cimerman, 1997) in two equivalents forms, the symmetric form:   Proceedings of COBEM 2003 17th International Congress of Mechanical Engineering COBEM2003 - 0483 Copyright © 2003 by ABCM November 10-14, 2003, São Paulo, SP            =                                ∑ +−− −−      ∑ + ≠ ≠ k i k k k ki ikkkk k ki i P P n E n E nnn nnn ............... 1 1 1121 11112111 ηηηη ηηηη ω (1) where the loss factors matrix of Eq. (1) is symmetric because of the reciprocity relationship: jjiiij nn ηη = (2) and the non symmetric form:           =                            ∑ +−− −−      ∑ + ≠ ≠ k i k ki kikkk k ki i P P E E ............... 1 21 12111 ηηηη ηηηη ω (3) Eq. (1) or Eq. (3) are typical models for SEA applications, which result from the relation, {} .,....,1 NiEEP N ij jjii N ij ijii =      ∑ −      ∑ += ≠≠ ηηηω (4) Eq. (4) represents the power balance equations for N coupled subsystems that form the basis for SEA. In the equations above, the parameter ni represents the modal density (number of resonant modes per frequency band), η ii the internal loss factor (damping parameter), η ij the coupling loss factor and the Pi the input power (average power introduced in the driven subsystem) is given by: )( * iii VFP ℜ= (5) where Fi is the point force excitation and * i V represent the complex conjugate of the velocity. Theoretical estimates of the coupling loss factor (CLF) for beams at right angles are available in the literature. These estimates are given as functions of the transmission coefficients ( τ ij) between two subsystems. The theoretical CLF may be determined by, (Cremer, et al., 1988), )2( iijiij Lc ω τ η = (6) where ci is the group velocity of the wave in beam i, Li is the length of beam i, and τ ij is the transmission coefficient across the joint relating the incident waves in subsystem i to the transmitted waves in subsystem j. The coefficients for each wave type may be calculated via the following expressions, 269 ; 269 58 ; 269 12 2 2 2 2 2 2 ++ = ++ + == ++ + = ββ β τ ββ ββ ττ ββ β τ LLLBBLBB (7) where LB cc= β with cB being the speed of the flexural waves, cL the speed of longitudinal waves, τ BB the transmission coefficient between incident flexural waves and transmitted flexural waves, τ BL the transmission coefficient between incident flexural waves and transmitted longitudinal waves, and τ LL the transmission coefficient between incident longitudinal waves and transmitted longitudinal waves. 3. Power Injection Method The power injection method (PIM) is usually used to determine the SEA parameters from experimental or numerical given data. PIM is based on the measurement of the power input into the subsystems and the kinetic energy, as an estimate of total energy, of these subsystems. Eq. (4) can be rearranged to express the CLF as unknown parameters. The time-averaged total energy of, say, subsystem i for an excitation in subsystem j is calculated from: 2 iiij VmE = (8) The derivation of subsystem energies from acceleration measurements of a set of points along the subsystem could present some difficulties. Firstly, Eq. (8) assumes that the total energy is equal to twice the kinetic energy. This statement is valid only if the energy is stored by modes excited at resonance (Cimmerman et al., 1997), and is also valid at higher frequency vibrations. Pavic (2001) showed that the kinetic and potential energies are similar close to resonance. This assumption can be verified through numerical modeling by SEM of a straight free-free beam excited transversally at one end. The material properties and dimensions of the beam are given in Tab. 1. Figure 1 shows the rate of kinetic and potential energy for the beam and the driving-point frequency response function (FRF). It can be observed that, at high frequencies, the kinetic and potential energies tend to be equal, and that they are equal at resonance. Secondly, the measurements are made at discrete locations that may not yield an accurate representation of the subsystem space averaged velocity. Table 1. Properties and parameters geometrical of straight beam free-free. Poisson ( ν )0.30 Area 1.000000e-4 m2 Density ( ρ )7800 kg/m3 Length 2 m Inertial Moment 8.333333e-10 m4 Elastic modulus (E)2.1e11 N/m 2 0100 200 300 400 500 600 700 800 900 1000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fre q uenc y, Hz Energy Ratio Kinetic/Total Potential/Total 0100 200 300 400 500 600 700 800 900 1000 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 Figure 1. Energies and driving-point FRF in a multi-dof system. The PIM can be used with confidence for a small number of subsystems but, for a larger number of subsystems, the energy matrix could be numerically ill conditioned. This is usually related to the different degrees of strength of coupling between the subsystems. Some strategies have been used to minimize the errors associated with the calculation of the CLFs using PIM. The technique used in this paper utilizes a sensitivity analysis of the CLFs (Stimpson and Lalor, 1991). Starting from the PIM matrix of two connected subsystems, and calculating the sensitivity of one CLF with respect to another CLF, one ends up with the relation jiE jjiiijji ≠ ≅∆ , η η η (9) where ji η ∆ is the variation in the factor η ji and Eij is the energy of subsystem i when power is input into subsystem j. Furthermore, the use of this relation results in a reduction of the PIM matrix and, thus, the CLF can be easily calculated using i in ij n ij n jj n ii n ij SL ji P E E EE E ω η =≅ , (10) Thus, the fact that the condition number of the PIM energy matrix (if used) reaches high values is not an issue when Eq. (10) is used. This expression is used in this paper for the estimation of the different CLFs for frame structures where beam members are connected at a right angles. 4. Review of the Spectral Element Method The main advantage of the Spectral Element Method is that the element dynamic stiffness is computed from the exact analytical solution in the frequency domain. Two different types of elements can be used in this method: 2-noded and throw-off. A spectral frame element consists of a combination of a bar (traction) element, a shaft (torsion) element and a beam (flexure) element. For the T-beam example, vibrations only in the plane x-y of the structure were considered. Therefore, torsion was not included, although it is straightforward to assemble a frame element including all six degrees of freedom per node. The spectral elements used here obey the following equations of motion for bar and beam elements, respectively, 2 2 t u A x u EA x∂ ∂ = ∂ ∂ ∂ ∂ ρ (11) 2 2 2 2 2 2 t I x v GA x EI x t v A xx v GA ∂ ∂ =      − ∂ ∂ + ∂ ∂ ∂ ∂ ∂ ∂ =       ∂ ∂ − ∂ ∂ φ ρφκ φ ρ φ κ (12) where EA is the axial stiffness, EI is the bending stiffness, u, φ and v are the axial, torsional, and flexural displacements, respectively, GA κ is the shear stiffness, ρ A and ρ I are the corresponding inertia terms, and κ is a geometrical constant that depends on the shape of the cross-section (5/6 in this case). Spectral analysis presents a solution of the form, ∑ =−− )( ˆ ),( tkxi eutxu ω (13) where k is the wave number. For the bar the equations of motion have a two-coefficient solution, )( ),( ˆxlkikx BeAexu −−− += ω (14) After the application of boundary conditions to a uniform wave-guide of finite length L with loads applied only to both ends, we end up with a system of equations for the bar element [] {} uk u u ee ee e Lik L EA F F L L L kiL L ik L L ikL L ki L L ki Lˆ ˆ ˆ ˆ 12 21 )1( ˆ ˆ 2 1 2 2 2 2 1=             +− −+ − =           −− −− − (15) where L k ˆ is the dynamic stiffness matrix for the bar element, F ˆ is complex amplitude of the applied force, u ˆ is the vector of complex amplitudes of the node displacements and L k is the wave number, defined as 2 1 2         =EA A kL ρω (16) In throw-off elements, waves propagate in one direction only. Thus, its dynamic stiffness matrix can easily be obtained by eliminating the term B in Eq. (14), which represents the reflected waves. Hence, the dynamic stiffness for the single-node infinite element is given by [iEAkL]. For the Timoshenko beam element, the four-coefficient exact solution is given by, )( 2 2 )( 1 1 2 2 1 1 )( 2 )( 121 ),( ˆ ),( ˆ xLikxLikxikxik xLikxLikxikxik DeRCeRBeRAeRxv DeCeBeAex −−−−−− −−−−−− −−+= +++= ω ωϕ (17) where R1 and R2 defined as the amplitude ratios and k1, k2 are the wave numbers defined as, 2 1 2 1 4 2 2 2 3 2 1 2 2 2 2 2 2 1 2,1 1 4 111 2 1 )(                                       −         +         ±                   +         = ω ω ωω c c cccc k (18) with the constants defined as 2 1 3 2 1 2 2 1 1,,         ≡         ≡         ≡A I c A EI c A GA c ρ ρ ρρ κ These solutions can be written in terms of the nodal displacements, and a relation between the applied shear forces and moments and the nodal degrees of freedom can be established as, {} [] {} uk L EI FBˆ ˆ ˆ 3 = (19) where B k ˆ is the dynamic stiffness matrix. It is symmetric and generally complex. The individual elements of this matrix can be found in Doyle (1997). An internal loss factor η can be included in all these wave numbers by using a complex Young modulus E(1+i η ). Only one element is needed between any two discontinuities, independently of its length. This plays the role of making the number of these elements in a 3-D structure relatively small. Thus, the response at different nodal degrees of freedom can be recovered with less computational cost by solving this system of equations at each frequency. These responses are then used to predict the total energy of a given structural element, for a certain wave type (longitudinal or flexural). 4.1. Energy in element bar at traction and compression The kinetic energy at bar element is caused by longitudinal waves propagation. This energy can be calculated through the equation, ∫ =L K RdxuAE 0 2 2 1& ρ (20) and calculated as time-averaged energy by: ∫ =L tt K RdxuAE 0 2 2 1& ρ (21) where {} {}{} )( ˆ )( ˆ 2 1 ) ˆ )( ˆ ( 2 1 ˆˆ 2 1*2**2 xuxuuiuiuuu tℜ=ℜ=ℜ= ωωω &&& (22) and thus, Eq.(21) can be written in the form: {} dxxuxuAE L t K R∫ℜ= 0 *2 )( ˆ )( ˆ 4 1 ρω (23) The solution for displacement )( ˆxu at any arbitrary point along the bar can be calculated using the shape functions (Ahmida, 2001). Thus, the displacement at point x, given the nodal solutions, is given by: 2211 ˆ )( ˆˆ )( ˆ )( ˆuxguxgxu += (24) and the expression to calculate of the kinetic energy is given by: { } 2 2 * 1 21 2 * 21 12 2 2 22 2 1 112 ˆˆˆˆˆˆ 4 1uuEuuEuEuEAE RRRR t K R+++ℜ= ρω (25) using the integral function: dxggEdxggEdxggEdxggE L R L R L R L R * 1 0 2 21* 2 0 1 12* 2 0 2 22* 1 0 1 11 ˆˆ ; ˆˆ ; ˆˆ ; ˆˆ ∫ = ∫ = ∫ = ∫ = (26) where * represent the complex conjugate. Note that 2112 RR EE =. This symmetry is seen in the dynamic stiffness matrix of the element spectral. Equations for the calculation of the potential energy can be developed following a similar procedure. In this paper, the Timoshenko beam element was used for the calculation of kinetic and potential energy. This theory is considered a high order theory, and the expressions for the calculation of kinetic and potential energies are given by Ahmida (2001). 5. Setup of Measurement System A T-beam made of Lexan, and which is continuos at the joint, was suspended by nylon fish lines at three points and excited in the transverse and axial directions at the end of branches A and C, see Fig. 2. The force was applied using an electromechanical shaker, which is driven by a periodic chirp signal in the frequency range of 0-8 kHz with intervals of 0.5 Hz. A signal generator and a data acquisition system (HP 3314A) were used. The shaker was connected to a stinger, which was connected to a force transducer PCB model 208A02. The inertance FRFs at 111 points along the structure were measured using three PCB 353B68 accelerometers connected to PCB 482AC5 conditioner. The data acquisition was conducted by LMS CADA-X software. A subsequent change of format was done in order to process data using Matlab. The analyses were done at 1/3-octave bands. 6. Application of the PIM In the estimation of the coupling loss factors by PIM, the subsystem energies were obtained from response measurements, using the Eq. (8). The T-beam was subdivided into six subsystems, three representing the energy of flexural waves and three representing the energy of longitudinal waves. The two kinds of waves coexist in the beam, the surface displacements on either side are a combination of the effects of them (Szwerc and Hambric, 1996). In this experiment, the measurements of longitudinal components were obtained by averaging the top and bottom measurements collected at a certain point along a beam, Fig. 3. This procedure was used and verified by Linjama and Verheij (1980). The flexural components were obtained only with one transversal measurement at each point, Fig. 3. The distance between the measurement points along the structure is 0.0251 m, see Fig. 4. Figure 2. Round-Robin T-beam. Figure 3. Measurement all consisting of 3 accelerometers. The planar behavior of the beam was investigated (no torsional waves) and four excitations were used, Fig. 4. For each excitation, 6 equations are derived. Using the structure symmetries, the number of the unknown coupling loss factors was reduced from 26 to 22, thus resulting in an over-determined set of linear equations. The subsystem energies were calculated from equally-spaced acceleration measurements in the transverse and longitudinal directions along the three branches (beams). Two different techniques were used. The first consists of using the lumped mass assumption and the hypothesis that the total energy is twice the kinetic energy. In the second, the displacement field is interpolated via the SEM and the kinetic and potential energies are obtained by integration of the analytical expressions. The former is shown to be valid only at frequencies where the modal density is high. Given that the set of linear equations is over-determined, a least squares solution was used via the Singular Value Decomposition (SVD). The experimental coupling loss factors in 1/3-octave bands are compared with the analytical expressions used for the three beams. Table 2.Properties of Lexan and Geometry of the T-beam. Density ( ρ ): 1280 kg/m3 Poisson ( ν ): 0.25 Area: 1.7118000e-3 m2 Jpolar: 3.6588149e-7 m4 Iz: 1.4334755e-7 m4 Loss factor 1e-2 Young’s modulus (E): 2.62 GN/m2 x Pin 3 A B C Pin 1 x y z Pin 2 Pin 4 Length Branch A: 0.7779 m Branch B: 1.0827 m Branch C: 0.9303 m x: 0.0251 m 13132 74 75 111 2 Figure 4. Positions of the excitations in the T-beam. 7. Results and discussion The kinetic energy used in the estimation of the coupling loss and internal loss factors was calculated through the lumped mass assumption and through SEM for a Timoshenko beam. The lumped mass corresponds to the hypothesis that the total energy is twice the kinetic energy. In spite of the assumptions made, the subsystem energies obtained with the lumped model were very close to the energies computed with the SEM model. In Fig. 5, Eij represent the energy of subsystem i due to excitation in the subsystem j. 102103104 -40 -35 -30 -25 Kinetic Energy E11 [dB ref. 1 joule ] Frequency [Hz] Lumped mass SEM-Timoshenko w/o φ 102103104 -60 -55 -50 -45 -40 -35 -30 -25 Kinetic Energy E61 [dB ref. 1 joule ] Frequency [Hz] Lumped mass SEM-Timoshenko w/o φ Figure 5. Kinetic energy by lumped mass assumption and via SEM for a Timoshenko beam without rotational kinetic energy. The coupling loss factors and internal loss factors were determined through the PIM. Fig. 6 and Fig. 7 show the CLFs and internal loss factor calculated theoretically and the experimentally using the lumped mass kinetic energy expression and via SEM. The majority of the experimentally determined CLFs compare well to the theoretically determined ones. However, some CLFs are overestimated, Fig. 7. An alternative approach, firstly mentioned by Stimpson and Lalor (1991), was used to fit these experimentally determined CLFs, Fig 8. The CLFs of Fig. 6 and Fig. 7 are defined as follows: • η 15 represent the CLF between flexural waves incident at branch A and longitudinal waves transmitted to branch C; • η 45 represent the CLF between longitudinal waves incident at branch A and longitunial waves transmitted to branch C; • η 54 represent the CLF between longitudinal waves incident at branch C and longitunial waves transmitted to branch A. 102103104 -40 -30 -20 -10 0 10 20 30 40 Coupling Loss Factor η 11 Frequency [Hz] CLF-Lumped CLF-Timoshenko CLF-Theoretical 102103104 -40 -30 -20 -10 0 10 20 30 40 Coupling Loss Factor η 15 Frequency [Hz] CLF-Lumped CLF-Timoshenko CLF-Theoretical Figure 6. CLFs obtained through PIM with the kinetic energy calculated via lumped mass assumption and via SEM expressions. 102103104 -40 -30 -20 -10 0 10 20 30 40 Coupling Loss Factor η 45 Frequency [Hz] CLF-Lumped CLF-Timoshenko CLF-Theoretical 102103104 -40 -30 -20 -10 0 10 20 30 40 Coupling Loss Factor η 54 Frequency [Hz] CLF-Lumped CLF-Timoshenko CLF-Theoretical Figure 7. CLFs obtained through PIM with the kinetic energy calculated via lumped mass assumption and via SEM expressions. 102103104 -40 -30 -20 -10 0 10 20 30 40 Coupling Loss Factor η 45 Frequency [Hz] CLF-Normalized-Lumped CLF-Normalized-SEM CLF-Theoretical 102103104 -40 -30 -20 -10 0 10 20 30 40 Coupling Loss Factor η 54 Frequency [Hz] CLF-Normalized-Lumped CLF-Normalized-SEM CLF-Theoretical Figure 8. CLFs obtained through the expressions of Stimpson & Lalor. In inverse problems, an important issue is the ill conditioning of the involved matrices. In the case of linear equations, as appear in SEA, the inverse problem is ill conditioned if the condition number is high. It is useful to observe the calculations through matrix condition number, which represent the ratio of the largest singular value to the smallest. Fig. 9 presents the condition number for the two different techniques used to obtain the energy matrix, which is thus used to determine the coupling loss factor. It can be seen that at low frequencies, the condition number is higher, as observed by Heckl and Lewit (1994). 0500 1000 1500 2000 2500 3000 3500 4000 4500 5000 102 103 104 105 106 107 108 Condition Number Frequency [Hz] CN-Lumped CN-SEM Figure. 9 Condition number in third octave bands. 8. Conclusion An experiment was conducted with the objective of estimating the coupling loss factors (CLFs) and the internal loss factors of a T-beam structure. The CLFs characterize the dynamic behavior at joints between beam-type subsystems. Thus, the prediction of vibration levels in such structures can be realized. For this purpose, the power injection method was used, which is very useful to gain information about the energy transmission paths in structures. Some parameters, such as subsystem energies and input power, are necessary for the estimation of the CLFs and internal loss factors. For the determination of the different subsystem energies and power input to the system, two different techniques were used: a lumped mass assumption and the spectral element formulations. The Timoshenko theory for beams was considered. The results of the subsystem energies obtained through the two techniques are compared and seem to be satisfactory. Although the lumped mass technique is an approximated method, as it treats the waveguide as non-continuous system, it still provided cost-effective and reliable solutions. On the other hand, the spectral element formulations are based on exact solutions and are only limited by the adopted waveguide theory. However, some CLFs seem to be overestimated in relation to the results obtained by analytical expressions given by Cremer (1988). Results were improved by using Stimpson and Lalor’s formulation (1991). The inverse problem can be accompanied by calculation of the condition number of the energy matrix, at each frequency component, which could give a good indication of the error of matrix sensitivity of the results. 9 Acknowledgements The authors are thankful to FAPESP and CNPq, for the financial support, to LMS International at its representative in Brazil, Smarttech, for the use of the CADA-X software, and to Dr. Steve Hambric, from Pennsylvania State University, for sending us the INCE-T-beam specimen for test. 10. References Ahmida, K. M., 2001, “Análise Dinâmica de Pórticos em Médias e Altas Freqüências”, Unicamp, Ph.D. thesis. Bies, D. A. and Hamid S., 1980,” In Situ Determination of Loss and Coupling Loss Factors by the Power Injection Method”, Journal of Sound and Vibration, 70 (2), pp. 187-204. Cimerman, B., Bharj, T., Borello, G., 1997,”Overview of the Experimental Approach to Statistical Energy Analysis”, SAE Noise and Vibration Conference 169, pp. 1-6. Cremer, L and Heckl, M., 1988, “Structure-Borne Sound”, Springer: Verlag, New York, 573 p.. De Langle, K., 1996, “ High frequency vibrations: contributions to expermental and computation SEA parameter identification techniques”, ISBN 90-73802-50-4, K.U. Leuven, Department of Mechanical Engineering, Belgium, thesis. Doyle, J. F., 1996, “Wave Propagation in Structures: a spectral analysis approach”, Springer: Verlag, New York, 324 p.. Heckl, M., and Lewit, M., 1994, “Statistical Energy Analysis as a tool for Quantifying Sound and vibration Transmission Paths”, Phil. Trasns. R. Soc. Lond. A., 346, pp. 449-464 . Lyon, R. H., Dejong, G. R., 1975, “Theory and Application of Statistical energy analysis”, MIT Press, USA, 277 p.. Pavic, G., 2001, “Vibration Damping, Energy and Power Flow”, Proceedings of the 17th International Congress on Acoustics, Vibrations and Structural Acoustics, vol I, Rome, Italy, 2 p.. Stimpson, G., Lalor, N., 1991,”Practical Noise Modelling of a Car Body Structures Using Energy Flow Analysis”, Proceedings Internoise, pp. 1233-1236. Szwerc, R. P. and Hambric, S. A., 1996,“The Measurement of Intensity of Longitudinal and Flexural Waves in Intersecting Beams”, Proceeding NOISECON 96, pp. 473-478. Verheij, J. W. , 1980, “Cross Spectral Density Method for Measuring Structure Borne power Flow on Beams and Pipes”, Journal of Sound and Vibration, 70 (1), pp. 133-139.