no 10 year 2008 volume 54 JOURNAL OF MECHANICAL ENGINEERING STROJNIŠKI VESTNIK Ol fltf.O ПЗ > (U > 47.5 t) £ 47 IS -419 g CO 33 ■41S Ц "S > 431 > u ё 4НШ Sound Presurre Level (all): 56.9 [dB(A)] Siusd lertKall) И J dt(A) K^iL-nt irv*:(n!i; a 2 ОДА) Sound Presurre Level: 52.2 [dB(A)] Strojniški vestnik - Journal of Mechanical Engineering ISSN 0039-2480 Editorial Office University of Ljubljana Faculty of Mechanical Engineering Journal of Mechanical Engineering Aškerčeva 6 SI-1000 Ljubljana, Slovenia Phone: 386-(0)l-4771 137 Fax: 386-(0)l-2518 567 E-mail: info@sv-jme.eu http://www.sv-jme.eu Founders and Publishers University of Ljubljana - Faculty of Mechanical Engineering University of Maribor - Faculty of Mechanical Engineering Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia - Metal Processing Association Editor Andro Alujevič University of Maribor, Faculty of Mechanical Engineering Smetanova 17, SI-2000 Maribor, Slovenia Phone: +386(0)2-220 7790, E-mail: andro.alujevic@uni-mb.si Publishing Council Jože Duhovnik, chairman Niko Samec, vice chairman Ivan Bajsić Jože Balič Iztok Golobic Mitjan Kalin Aleš Mihelič Janja Petkovšek Zoran Ren Stanko Stepišnik International Advisory Board Imre Felde, Bay Zoltan Inst, for Materials Science and Techn. Bernard Franković, Faculty of Engineering Rijeka Imre Horvath, Delft University of Technology Julius Kaplunov, Brunei University, West London Milan Kljajin, J.J. Strossmayer University of Osijek Thomas Lübben, University of Bremen Miroslav Plančak, University of Novi Sad Bernd Sauer, University of Kaiserlautern George E. Totten, Portland State University Nikos C. Tsourveloudis, Technical University of Crete Toma Udiljak, University of Zagreb Arkady Voloshin, Lehigh University, Bethlehem Deputy Editor Editorial Board Vincenc Butala Anton Bergant University of Ljubljana, Faculty of Mechanical Engineering Franci Cuš Aškerčeva 6, SI-1000 Ljubljana, Slovenia Matija Fajdiga Phone: +386(0)1-4771 421, E-mail: vincenc.butala@fs.uni-lj.si Jože Flašker Janez Grum Technical Editor Janez Kopač Darko Švetak Franc Kosel University of Ljubljana, Faculty of Mechanical Engineering Janez Možina Aškerčeva 6, SI-1000 Ljubljana, Slovenia Brane Širok Phone: +386(0)1-4771 137, E-mail: info@sv-jme.eu Leopold Škerget JOURNAL 0F MECHANICAL ENGINEERING STROJNIŠKI VESTNIK Cover plate: Figures represent sound patterns of a larger industrial plant, recorded with an acoustical camera. Figures show identified main noise source (in total - above) and two noise sources at a particular frequency range (below). Visualisation of noise sources in different frequency ranges is necessary for a successful noise control. Print Littera Picta, Medvode, printed in 600 copies Yearly subscription companies individuals students abroad single issue 100,00 EUR 25,00 EUR 10,00 EUR 100,00 EUR 5,00 EUR © 2008 Journal of Mechanical Engineering. All rights reserved. SV-JME is indexed in SCI-Expanded, Compendex and Inspec. The journal is subsidized by the Slovenian Research Agency. Strojniški vestnik - Journal of Mechanical Engineering is also available on http://www.sv-jme.eu, where you access also to papers' supplements, such as simulations, etc. Contents Strojniški vestnik - Journal of Mechanical Engineering volume 54, (2008), number 10 Ljubljana, October 2008 ISSN 0039-2480 Published monthly Papers Cernetic J., Prezelj J., Cudina M.: Simple Feedback Structure of Active Noise Control in a Duct 649 Vladić J., Malešev P., Šostakov R., Brkljac N.: Dynamic Analysis of the Load Lifting Mechanisms 655 Alic G., Širok B., Schweiger V.: Impact of the Guard Grill on the Integral and Local Characteristics of an Axial Fan 662 Pavlović M., Stanojević M., Ševaljević M., Simić S.: Influence of The Waste Oil Concentration in Water on The Efficiency of The Aeration Process in Refinery Wastewater Treatment 675 Balicević P., Kozak D., Mrcela T.: Strength of Pressure Vessels with Ellipsoidal Heads 685 Braut S., Žigulić R., Butković M.: Numerical and Experimental Analysis of a Shaft Bow Influence on a Rotor to Stator Contact Dynamics 693 Žula T., Kravanja Z., Kravanja S.: MINLP Optimization of a Single-Storey Industrial Steel Building 707 Galeta T., Kljajin M., Karakašić M.: Geometric Accuracy by 2-D Printing Model 725 Instructions for Authors 734 Paper received: 06.09.2007 Paper accepted: 07.07.2008 Simple Feedback Structure of Active Noise Control in a Duct Jan Černetič* - Jurij Prezelj - Mirko Čudina University of Ljubljana, Faculty of Mechanical Engineering, Slovenia An active noise control is usually constructed with the use of electronic filters. For sufficient noise attenuation, electronic filters are not always needed. In this case, the electronic controller and appropriate software is not required so the system can be much easier. This paper deals with the use of a feedback structure of active noise control in an experimental ventilation duct. Simulation was performed to investigate the efficiency of a simple analoguous system for active noise control without incorporating electronic filters. The transfer function of the entire analogous system can be used to predict the maximum attenuation level. Tests were made to verify the simulation and to show what noise attenuation level can be achieved in an experimental duct. It has been shown that in a specific frequency range this kind of a system is efficient enough for use in some ventilation ducts. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: active noise control, feedback control, ANC simulation, phase shift 0 INTRODUCTION The beginning of active noise control (ANC) is in year 1936, when Lueg patented his idea about realization of the active noise control [1]. It uses the principle of interference and absorption. Because electronics in those years was not enough advanced to meet the needs of a controller, no real system was produced. Today, two basic methods are used for active noise control; feed-back and feed-forward. They are using all capabilities of contemporary computers and other electronics. These systems are used mainly in conjunction with adaptive filters, which make them capable to cope with bad system response. But there is a question, if it is possible (in some cases) to use a simpler system, which is efficient enough and at the same time cheaper and more reliable because of less electronic components. This is the purpose of the experimental ventilation duct that was made in this research. The characteristic of ventilation ducts' fans is a constant rotating speed. This means that the emitting noise is constant and tonal in lower frequency band. The intention of this test is to achieve good noise attenuation with feedback structure of an ANC system without incorporating adaptive or other electronic filters. These filters are capable to compensate bad system response, and at the same time, they contribute to considerable phase shift of the original signal. This part of phase shift can be eliminated if the filter is not included in the feedback loop. It should be aware of the fact that, for example, an amplifier itself also contains filters. In this research, no additional filter was used. 1 FEEDBACK STRUCTURE OF ACTIVE NOISE CONTROL Feedback structure was found in 1953, when feedback loop was constructed of a microphone, an amplifier and a loudspeaker [2]. Authors have named this scheme as the »sound absorber«. In this scheme, microphone is used as an error microphone, where error means deviation from the theoretical attenuation. Signal is travelling through the amplifier and the controller, where the amplitude is adjusted and phase is shifted for 180°. Modified signal goes to the loudspeaker, where the attenuation of noise occurs. Schematic representation of a feedback structure is shown in Figure 1. pnmary source errar microphone secondary source W" controller Fig. 1. Schematic representation of a feedback structure of ANC Today, a feedback structure is sometimes used as a way to avoid acoustic feedback in feed- Corr. Author's Address: University of Ljubljana, Faculty of Mechanical Engeenering, Aškerčeva 6, Ljubljana, Slovenia, jan.cernetic@fs.uni-lj.si forward structure. A feedback ANC configuration is also very suitable when there is no available reference signal [3]. To achieve maximum level of attenuation, signal from the microphone should be amplified as much as possible, but should not exceed the stability range. Block diagram, corresponding to a feedback control system in the complex frequency s-plane, is shown in Figure 2, where E(s) is a transfer function of acoustic path between the primary source and the microphone, F(s) is a transfer function of acoustic path between the secondary source and the microphone, M(s) a transfer function of the microphone, C(s) a transfer function of the controller, N(s) a transfer function of the electronics (except the controller) and L(s) a transfer function of the secondary source. Typical characteristics of this method is a feedback loop, which leads signal S(s) from the secondary source back to the microphone through the acoustic path. There it is added to the reference signal R(s) and the error signal D(s) is therefore increased. Signal from the error microphone D(s) can be written as: D(s) = R(s) + F (s)S (s) (1) or D(s) = E (s) P(s) + +F (s)M (s)C (s) N (s) L(s)D(s) (2) It leads to a simplified form of a transfer function between a microphone signal D(s) and a signal of primary source P(s): D(s) E (s) (3) P(s) 1 - M(s)C(s)N(s)L(s)F(s) Because the error signal should be as low as possible (theoretically vanish), it is D(s) = 0 . Further derivation shows that for the best attenuation, the transfer function of the controller should have an infinite amplitude. Because of the stability problem, this is not possible. When a feedback loop is amplified too much, other elements may cause the system to become unstable. Therefore an optimum operational point should be found, which is normally near the stability limit. 2 SYSTEM STABILITY When constructing the ANC system with a feedback loop, the stability problem should always be observed. Position of the microphone and the secondary source has strong influence on stability. In theory, it is considered that both elements are infinitely close to each other. But this is not possible to reach because of the following reasons. First, it should be taken into account that both elements can not be placed close to each other because a near-field influence would be too high to get good accuracy. Second, dimensions of a loudspeaker and a microphone are preventing to define their optimum placement. For this reason, there is always some gap between the microphone and the loudspeaker, which affects the efficiency of the whole system. If the microphone and the loudspeaker would be so close to each other that the sound waves around a microphone would be defined solely by the movement of the loudspeaker, the microphone would get the same signal as the loudspeaker produces (no phase shift). This is not possible and the gap causes some time delay in the signal, which results in phase shift. If the wavelength of the signal is significantly greater than a gap between the microphone and the loudspeaker (intermediary length, l), phase shift is relatively small. Hs) E(s) Rs) + ►v D(s) * M(s) C(s) N(s) L(s) S(s) F(s) Fig. 2. Block diagram of a feedback control system of ANC But if the wavelength is similar to the intermediary length (this occurs at the higher frequencies), phase shift becomes relatively large and efficiency of the system drastically decreases. If the frequency is increasing (at some fixed intermediary length), the phase shift is therefore also increasing, while the efficiency of the system is decreasing. Besides, it should be also taken into account the property of the real acoustical and other components, that each of them causes some phase shift in the signal. This is more significant for loudspeakers and filters. When (for a specific sound frequency) the intermediary length corresponds to the equation (4), the phase shift between the sound, produced by the loudspeaker, and the sound, captured by the microphone, is equal to 180°: l = * (4) 2 This means that the work of the controller, which reverses the phase for 180°, is totally ineffective and instead of attenuation it causes amplification of the primary noise. Real noise, which is wanted to be attenuated, consists of many narrow frequency bands. Therefore, at every single intermediary length, one narrow frequency band exists, which causes the noise at this frequency to be amplified, and not attenuated. When the amplification is increasing, the amplitude is also increasing and at the specific point it exceeds the stability limit. Unstable working conditions occur and the system is not working appropiately anymore. The longer the intermediary length, the lower the limit frequency for stable operation. Many systems of the ANC are constructed with the use of electronic adaptive filters [4] and [5]. They are able to partly compensate non-ideal response of each component and other causes of reduced efficiency. One of the filter problems is, that it causes phase shift in the signal, which results in lower efficiency. Therefore it is interesting to find out how these filters can be avoided. 3 MEASUREMENTS The purpose of this experiment was to construct an ANC system with a feedback loop in an experimental ventilation duct without any additional compensation filter, and to get the noise attenuation, which is good enough. The setup was constructed of plywood with thickness of 2 cm (Figure 3). For simulating a fan, the loudspeaker with generated broadband pink noise was used. Narrowband or tonal noise seems to be a better choice for simulating a fan, but the experiment was intended to examine the efficiency of the system through the wider frequency spectrum. /R4* M2 M1 Fig. 3. Experimental ventilation duct scheme A feedback loop, used for generating anti-noise, was constructed of a microphone Bruel&Kjaer (BK) 4165 (M2 in Fig.3), a measuring amplifier BK 2636 and a loudspeaker Visaton W170S8 with reversed polarity (L2 in Fig. 3). For primary noise generation (pink noise) the following equipment was used: a sine random generator BK 1027, a power amplifier BK 2706 and a loudspeaker JVC CS-HX621 (L1 in Figure 3). 3.1. Frequency Response of a Loudspeaker The loudspeaker for anti-noise generation was built in an open box. According to its casing, theoretical response of the loudspeaker was calculated by the program WinISD 0.44. From the theoretical frequency response in Figure 4 it can be evident that because of resonant frequency at around 30 Hz, the loudspeaker is suitable for the use in a frequency range above 40 Hz. Frequency [Hz] Fig. 4. Theoretical frequency response of the loudspeaker Frequency range, appropriate for noise attenuation, is seen in Figure 5, which shows theoretical group delay of the loudspeaker. At 60 Hz and lower, a group delay is too long for the system to work properly. 3.2. Frequency Response of a System Many components of the ANC system have influence on the effectiveness of noise attenuation and each of them has its own transfer function. The main reason for delayed and modified signal in a simple ANC system (without electronic filters) is the loudspeaker, other components (microphone, amplifier, cables, ect) contribute less. Measurement of the input and output was performed in the open loop (Figures 6 and 7), from which the transfer function can be calculated. It is obvious that the signal is very modified at the end of the loop, after it passes all the components of the ANC system. Such a simple system is very hard to compensate, because some frequencies are filtered out by the transfer function of the system. For appropriate compensation, the amplification of these particular frequencies with zero amplitude should theoretically be infinite, which is not possible. 3.3. Noise Attenuation Measurement Then the noise attenuation with the ANC system was measured. The microphone, which is used for measurement of the attenuation level (M1 in Figure 3), was placed near the exit of the experimental duct, 45 cm from the end. Noise spectrum with and without use of an ANC, measured with the microphone M1, is seen in Figure 8 and the attenuation level in Figure 9. It is obvious that the useful frequency range of the system is approximately between 40 and 140 Hz. 90 85 = 65 ci E < 60 55 50 45 Frequency [Hz] Fig. 8. Attenuated (black line) and non-attenuated (grey line) noise in the duct A V/ J\ \ \ л V w и Frequency [Hz] Fig. 6. Input signal into the system Frequency [Hz] Fig. 9. Attenuation level in the duct 80 75 0) 70 10 83 81 79 77 7b 4 E 73 71 69 67 10 1000 With regard to the fact that some ventilation fans produce noise mainly in lower frequency ranges, this ANC system could be very useful in that cases, despite of its narrow working range. Better attenuation could be measured if the primary noise would be narrowband or tonal noise, which is in fact more significant for ventilation fans. 4 SIMULATION When constructing an ANC system, it is important to know its working limits. Simulation of the ANC was performed, using Matlab 7.1, to predict the effectiveness of the system. It shows the possibilities of a particular ANC system in a sense of how much the noise can be attenuated. For a simulation to perform, the input and output of the system with the open loop were measured (as mentioned before). These two signals were used to calculate the coefficients of the impulse response, with the help of an LMS filter (Figure 10). From the impulse response coefficients, the response of a closed loop can be calculated. For a closed loop, the following equation can be written: x(m) = p(m) + f (m), where f (m) = K [x(m)a1 + x(m - 1)a2 +... ... + x(m - (L - 1))aL ] (5) (6) K is gain of a feedback loop, m goes from L to the end of the signal and L is the number of filter coefficients. Extended equation is: x(m) = p(m) + K [x(m)a1 + + x(m - 1)a2 +... + x(m - (L — 1))aL ] (7) The final equation of the output (attenuated) signal can be written as: x(m) =-[p(m) + Kx(m - 1)a2 + 1 - Kal + Kx(m - 2)a3 +... + Kx(m - (L - 1))aL ] (8) In this way the transfer function of the whole system was simulated by the computer program. To represent the efficiency of the system, the input signal is compared with the output signal of the closed loop, which comes from the computer algorithm. The results of the simulation are represented in Figure 12 (attenuated and original signal) and in Figure 13 (attenuation level). As expected, the simulated noise attenuation level is higher than the real measurement, because this is the best possible (ideal) attenuation for this specific ANC system and the equipment used. Coefficient number Fig. 10. Impulse response coefficients According to Figure 11, p(m) is measured input and x(m) is output of a closed loop. Filter coefficients are marked as ai, where index i means its instantaneous number. P(m) x(m) 4+ f(m) Fig. 11. Block diagram of a feedback loop L 100 Frequency [Hz] Fig. 12. Original input (grey line) and simulated attenuated output signal (black line) 10 a But the useful frequency range is more similar to that in measurement, which confirms that the simulation is correct. A Jrtfrl ------ -J- - -1Ш /Улш4 Frequency [Hz] Fig. 13. Attenuation level of the simulation When evaluating the results, some other influences should also be taken into consideration. The result is in strong correlation with the intermediary length in a feedback loop and also with the position of the microphone with regard to the loudspeaker of the secondary source. This is because sound reflections are present in the duct and the sound is travelling in different directions. The more the microphone is far from the near-field, the more the reflections and other phenomena influence on the signal, consequently lower the efficiency of the system is. Besides, the result also depends on the place, where the attenuation is measured (microphone Ml in Figure 3). Because of the impedance mismatch between the duct and the surroundings, a part of the sound waves is reflected back to the primary source, which leads to the standing waves phenomena [6] and [7]. On the specific points of the duct at the specific frequency, nodes are formed, which mean that at that points the measured noise level would be very low. 5 CONCLUSIONS The experimental ventilation duct with a feedback method of an ANC was constructed. The possibilities of the ANC without additional electronic or other filters were investigated, because an important part of time delay of the signal may be caused by the filter. The LMS filter and the measured input and output signal were used to do the simulation of an ANC system. It shows the maximum possible attenuation for a specific ANC system. Then the measurement of the attenuation level was performed. The experiment showed that the additional filters can be avoided, if only a specific frequency range of noise (about 40 to 140 Hz) must be attenuated. Some ventilation fans can meet this requirement. In this case, the ANC system becomes very simple and robust, because the feedback loop is constructed of fewer components. It contains just a microphone, an amplifier and a loudspeaker. During a design process, an engineer should be aware of the fact that the efficiency of the system depends on the position of the error microphone with regard to the loudspeaker for the anti-noise generation. They should be located close to each other, but not too close to become influenced too much by the near-field. The noise, used for measurements and simulations, was broadband pink noise. Better attenuation is expected in case of using narrowband or tonal noise, which are more similar to the real ventilation fan noise. 6 REFERENCES [1] Tokhi, M. O.Leitch, R. R. (1992) Active Noise Control, Clarendon Press, Oxford, USA. [2] Olson, H. F., May, E. G. (1953) Electronic Sound Absorber, J. Acoust. Soc. Am. 25(6), p. 1130-1136. [3] Prezelj, J. Čudina, M. (2007) Dipole in orthogonal direction as a secondary source for active noise control in ducts, Acta Acustica united with Acustica 93, p. 63-72. [4] Sakai, H., Miyagi, S. (2003) Analysis of the adaptive filter algorithm for feedback-type active noise control, Signal processing 83, p. 1291-1298. [5] Sun, X., Kuo S.M., Meng, G. (2006) Adaptive algorithm for active control of impulsive noise, Journal of Sound and Vibration 291, p. 516-522. [6] Munjal, M. L. (1987) Acoustics of ducts and mufflers, John Wiley and sons, New York. [7] Romeu, J. Saluena, X. Jiménez, S. Capdevila, R. Coll Ll. (2001) Active noise control in ducts in presence of standing waves. Its influence on feedback effect, Applied Acoustics 62, p. 3-14. 25 20 15 10 10 Paper received: 08.04.2008 Paper accepted: 07.07.2008 Dynamic Analysis of the Load Lifting Mechanisms Jovan Vladić* - Petar Malešev - Rastislav Šostakov - Nikola Brkljac University of Novi Sad, Faculty of Technical Sciences, Serbia This paper deals with the problem of dynamic behaviour of load lifting mechanism (such as elevators). In the case of considerable lifting heights, high velocity devices are applied, with the purpose of shortening cycle duration and increasing the capacity. In such case, the standard procedure of dynamic analysis is not applicable. In the paper, the procedure of establishing the appropriate dynamic model and corresponding equations is proposed. It enables the analysis of the relevant influences, such as variation of the rope free length, slipping of the elastic rope over the drum or pulley and damping due to the rope internal friction. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: elevators, pulleys, dynamic models, rope length variation, motion stability analysis 0 INTRODUCTION In the most types of load lifting mechanisms with considerable lifting heights operating under usual conditions, the steel wire-rope and the driving drum (in cranes) or the driving pulley (in elevators) are applied. The paper deals with the influence of these components upon the dynamic behaviour of the lifting mechanism, while other influences were left out of the scope of the research. These influences have already been investigated in details, e.g. dynamic phenomena in the gear box, in [5], the dynamic characteristics of the driving electromotor, in [12], etc. The research is based on the system with many degrees of freedom, also used in [13]. The analysis of the lifting mechanism with relatively low lifting velocity is usually performed by using dynamic models based on the longitudinal vibrations of elastic, homogeneous bar of a constant length, with or without mass, corresponding to certain boundary conditions, e.g. in [10]. The stability of beam acceleration and the effect of various stiffness values have been analysed in [6] and [11]. Boundary condition at the lower rope end is defined by the load mass and at the upper rope end by the moment of inertia of all driving mechanism components and by the excitation force corresponding to the mechanical characteristic of the driving electromotor, both reduced to the shaft of the driving drum or pulley. But, in the case of higher lifting velocity, (e.g. in the case of fast passenger elevator, mine shaft winding system, harbour crane, etc.) the numerous influences upon the dynamic behaviour of the lifting mechanism are also to be taken into account, such as: • rope length variation during load lifting and lowering, • steel wire-rope slipping over the driving pulley or drum, • mechanical properties of the rope (such as elasticity modulus, rope internal friction, rope design, etc.), • influence of dynamic processes on the incoming rope side upon the dynamic processes on the outgoing rope side in the systems with driving pulley, • mechanical characteristic of the driving motor, • influence of friction processes between the elevator car guiding shoes and rails. Basic influence of the rope length variation appears directly through the variation of the basic dynamic parameter - stiffness of the rope free length. There is a considerable possibility that in the case of realistic devices with higher velocities and minor energy losses due to internal friction in elastic elements, the variable stiffness of the rope free length could induce the occurrence of parametric vibrations. 1 EQUATIONS OF LOAD LIFTING SYSTEM A load lifting system with the driving pulley is presented in Fig. 1. It represents a more complex system for dynamic analysis in comparison with the system which contains a drum as a driving element. *Corr. Author's Address: University of Novi Sad, Faculty of Technical Sciences, Trg D. Obradovića 6, 21000 Novi Sad, Serbia, vladic@uns.ns.ac.yu Fig. 1. Elevator driving mechanism The influence of the elasticity of shafts, gears, etc. can be neglected in comparison with the steel wire-rope elasticity. This enables the establishment of a simplified dynamic model of the driving mechanism, as shown in Fig. 2. The model takes into account the influence of the rope free length variations on both the incoming and outgoing pulley sides. The variations in rope free length (shortening or lengthening) directly affect the rope stiffness, and therefore also the dynamic behaviour of the rope, which is of significant importance in the systems with high lifting velocities. The influence of elastic rope slipping over the driving pulley is considered through the adequate boundary conditions. According to [9] and Fig. 2, on the basis of the equilibrium conditions for the elementary rope length and for the driving pulley, the following system of equations can be established: q д uJxJ) d2 —---1- = E- A-- dt2 dx2 +q-(1± a ) g u1(x,t)+b- du^ (x,t ) dt (1) q du2(x,t) g dt2 E- A- dx2 du (x,t) u (x,t)+b—2- dt (2) +q- (1± a) Mm = -E- A-dd i-ц dx u1(l1,t) ~u2(l2,t) "[[[о ~u2(l2,t)\ ] ■ ^ (3) with: u1, u2 E A a Mm rope elastic deformations, m elasticity modulus, Pa rope cross-section, m2 driving mechanism acceleration, ms-2 driving motor torque, Nm i - gear ratio, r\ - driving mechanism efficiency, Jr - moment of inertia of rotating masses, reduced to the pulley shaft, kgm2 R - driving pulley radius, m q - rope weight per meter, Nm-1 g - gravity acceleration, ms-2 In the presented Equations (1) to (3), the wire-rope is considered to be a viscous-elastic body (Kelvin's model), so the internal force can be defined as: Fig.2. Dynamic model S ( x, t ) = E - A — dx u ( x, t ) + b du ( x, t ) dt (4) According to [2], dependence of the damping coefficient upon the rope stress can be adopted as: b = 1 0.5 + 2300 350 + a 10" with: a - rope stress, MPa In the course of solving the Equation system (1) to (3), it is necessary to define the electro-mechanical model of the driving motor, or its mechanical characteristic in the form Mm=f(v). Fig. 3. Diagram of velocity and acceleration for the elevator with double-speed motor For the passenger elevators, the parameters defining the "driving comfort" are specified, limiting the values of maximum acceleration to: a = 1.4 ms-2 and its time derivate to: (da/dt)max = 1.0 ms-3, [7]. This makes it possible to simplify the analysis by excluding the Equation (3) from the previous system. Diagram of velocity and acceleration for the elevator with double-speed motor is presented in Fig. 3. When considering the whole driving, the analysis becomes very complex. Therefore, the paper presents only the analysis of dynamic behaviour concerning the rope incoming side of the driving pulley. In this case, the further analysis is to be carried out only with respect to the Equation (1). 2 BOUNDARY CONDITIONS FOR THE ROPE INCOMING SIDE OF THE DRIVING PULLEY In Fig. 4, the characteristic parameters satisfying the boundary conditions for the rope incoming side on the driving pulley or drum are presented. Boundary conditions at point C, where the rope makes the first contact with the pulley, are dependent on whether the winding is with or without rope slipping. The change of force in the wound rope part can be maintained only if this change is smaller than the adhesive force making it possible. Fig. 4.d shows different cases of force distribution over the wound rope length as a function of elevator velocity, [8]. Fig. 4. Boundary conditions for: a) the pulley without slipping, b) the pulley with slipping, c) the elevator car (or counterweight), d) the distribution of forces and slipping over the wound rope length vi>v2>v3 The limiting value of the rope winding velocity is determined from the condition: dS (l1, t) x for y > x and displacement as u(x) =K (x, y, l)-Fi. In case of several forces acting in points at x = yi, the displacement has the following form: n u ( x) = Z K ( x yi,l ) ■ Fi i=1 For the load case with evenly distributed rope load mass, the expression takes the form: L u(x) = JK(x,y,l) ■ q ■ dy i After applying this procedure to the differential Equation (1), replacing floating argument with (y), multiplying it with the function K (x, y, l1) = K and performing necessary mathematical transformations, e.g. [1] and [4], the rope deformation becomes as follows: 4 u1(xt) = -JK ■qi y) ■ h + Jdu1(lh t) 0 d1u1 d2 — g-a ■dy- dx d± V d (Al1) dt J dt ■ dt - + Al + E A k ■ _d_uL_ ■ds - dx J " 2 (11) ds2 ■ dt - b ■ e ■ A-^u^ ■ K (x, L1, lO dx dt Applying the method of particular integrals, e.g. [10], given as: Ü1 (x, t)= X (x)T1 (t) and considering only the first mode of vibrations, form of which - in the case of a rope with the weight at its lower end - can be adopted as the straight line, i.e. xl(x) = x-1, the simple differential equation in the following form is obtained: T ■T1 =£■ R(l1,v)■ T (12) with: ®2(l1) = g ■ E ■ A Q+ q ■( L1-11) ( L - lo rope fundamental frequency, s- R ((1, v ) = L -1 v - b - E ■ A ■ g Q + q ■(L1-11 ) s = L1 (l1 ) Differential Equation (12) describes nonlinear vibrations with viscous friction when the parameters are "slow-time functions", with the solution, according to [10], in the form: T1 (t ) = a1 (t )■ cos e(t ) whereby, after some mathematical transformations, the obtained values are: (t ) = h, 10 L1 l10 M -1 \ 0.25 1 J with: 91 (t)= Ja\ (()■ d t + 910 0 h0 Q . a h10 = E ■ A g amplitude value of the function, ^0 - length of wound rope at t=0, m m = 2 ■ ъ-м-'Ьгт 2 dt 0 L1 — /1 Jft12(/) ■ dt After substituting previously mentioned values into Equation (9), and after performing mathematical transformations, the equation obtaines the form: v Mj(x, t) = (x —11) ■ a 10 ^ L l10 L1 -11 i.e.: ■ cos \a1(l )dt Q- q ■ (2 ■ L + x +11) a I x — li 1 + — I--1 E ■ A .(I1, t ) (13) On the basis of Equations (4) and (13), the force in the rope is expressed as: 51 (x, t ) = E ■ A ■ a1 (t )■ cos 61 (t ) + + [Q + q ■(L — x)](1 + a1 V g J (14) Analyzing Equation (14), it can be concluded that, during the rope free length decreasing phase (load lifting), under certain conditions, its deformation can be increased, causing a permanent growth of rope dynamic load. Such phenomenon, usually described as unstable lifting, appears in the case that da1 (t) > 0. On the basis of this condition, the dt critical lifting velocity can be determined as: I o К v > vS (15) К + vs v < vs The critical lifting velocity in the case without rope slipping over the pulley at the incoming point C is: o 2 ■ b ■ g ■ E ■ A Q q with: for: is: ■■2 ■ b 3 q ■ E ■ A Q (16) (17) Q >> q ■( L - h ) Limit slipping velocity, according to (5), •(1 ) = k ■ R H v 2 ■ b ■(L —11 ) (18) k ■ R q ■ E ■ A ß IQ '(L — I If the elevator velocity exceeds vc in the lifting phase, the dynamic load will increase (unstable movement). Under these conditions two cases are possible: — If the rope slipping at the incoming point is prevented (pulleys with special clamps for rope gripping), the unprevented increase of the dynamic load appears in the lifting phase, — At standard pulleys, the increased dynamic load also appears in the lifting phase, but only until the moment when the slipping begins in the rope incoming point, v < vs. During the rest of the lifting process, the load decreases, i.e. the stable movement occurs. The value of critical velocity depends on the rope damping characteristics (due to internal friction), slipping conditions in the rope incoming point (C) and the basic characteristics of the elevator. When elevators with high load capacity, high velocity and lifting height (fast passenger elevators and mine elevators) are concerned, the lifting velocity can exceed the critical velocity, making it necessary to check the stability of the lifting process already in the design phase. 4 CONCLUSION For driving mechanisms used for vertical load lifting (elevators and cranes), it is necessary to provide adequate conditions for the correct dynamic analysis of their parts' behavior and especially for the rope as the basic element, already in the design phase. Due to a significant influence of rope free length changes, due to its slipping over the driving pulley and for the reason of its mechanical characteristics, it is impossible to apply a classical dynamic model of longitudinal oscillations for homogenous stick of a constant length, especially by the elevators with large velocities and large lifting heights (express and mine elevators). The introduced dynamic model for elevator driving mechanism provides the analysis of relevant parameters and a base for the assessment of lifting process stability through the critical velocity level, defined with the help of expressions (15) to (19). v v 5 REFERENCES [1] Abramowitz, M., Stegun, I. (1979) Handbook of Special Mathematical Functions, Applied Mathematics Series-55, Mir Publishers, Moscow. [2] Goroshko, A., Artjuchova, E. (1968) Dependence coefficient of damping in rope of its load, Steel wire-rope, 5, p. 57-59, (in Russian). [3] Goroshko, A., Savin, G. (1971) Introduction to mechanics of elastic bodies with variable lenght, Science, Kiev, (in Russian) [4] [4] Mitrinović, D.S. (1975) Introduction to the Special Functions, Civil Engineering Book, Belgrade. [5] Hlebanja, J., Duhovnik, J. (1981) Power Transmission Between Teeth at Small Loads, Proceedings of the International Symposium on Gearing and Power Transmissions, Tokyo, August 30 - September 3, Vol. II. [S. l.: JSPE, 1981], pp. 43-47. [COBISS.SI-ID 7709657]. [6] Pakdemirli, M., Ulsoy, G. (1997) Stability Analysis of an Axially Accelerating String, Journal of Sound and Vibration 203(5), p. 815-832. [7] Thiemann, H. (1979) Elevators, Veb Verlag Technik, Berlin, (in German). [8] Vladić, J. (2002) Characteristics of Power Transmission in Systems with Driving Pulley, Journal of Mechanical Engineering Design 5(1), p. 19-27. [9] Vladić, J., Zlokolica, M. (1998) On the Dynamic Analysis of the Load Lifting Machines, Proceedings of ICVE98, Dalian, China, p. 100-103. [10] Vujanović, B. (1977) Theory of Oscillations, University of Novi Sad, Faculty of Technical Sciences, (in Serbian). [11] Wickert, J.A., Mote, C.D. Jr. (1990) Traveling Load Response of an Axially Moving String, Journal of Sound and Vibration 149, p. 267-284. [12] Slavic, J., Nastran, M., Boltezar M. (2006) Modeling and Analyzing the Dynamics of an Electric-motor Brush, Journal of Mechanical Engineering 52(12), p. 126-142. [13] Kulvietiene, R., Kulvvietis, G., Tumasoniene, I. A (2006) Symbolic/ Numeric Vibrations Analysis of System with Many Degrees of Freedom, Journal of Mechanical Engineering 52(12), p. 309-317. Paper received: 10.01.2008 Paper accepted: 07.07.2008 Impact of the Guard Grill on the Integral and Local Characteristics of an Axial Fan Gregor Alic 1 - Brane Širok2 - Vlado Schweiger1 1 Institute Hidria Klima, Godovic, Slovenia 2 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia This paper presents the impact of the guard grill (GG) on the integral characteristic (IC) and local aerodynamic properties of the outflow of an axial fan (AF). The study of the GG's impact was conducted using measurements and numeric modelling at the integral and local levels. Measurements of IC were performed according to the ISO 5801 standard and the local velocity field was measured by using a five-hole probe. Numeric simulations of the machine's local characteristics were made in the same working conditions of the fan. The comparison of results was used in the analysis of the GG's impact on the fan's efficiency and in the analysis of the feasibility of the CFD modelling of flow conditions of the GG. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: axial flow fans, guard grids, integral characteristic, velocity fields 0 INTRODUCTION In industrial practice the development of an AF is mostly performed on the assumption of the decisive impact of the rotor and the inlet of the fan on operating conditions, while the impact of the GG and the supporting construction is empirically evaluated or ignored. No research has yet been performed on the GG's actual impact and the supporting construction of the AF on the precise determining of the fan's operating characteristics and outflow kinematics. Eck [2] and [3] presented the analytical procedures of designing axial and centrifugal fans, while Wallis [1] discuss the analytical procedures for designing AF together with inlet and outlet channels. The programming and analytical design of a fan based on a similar basis as described in the literature [1] to [3] was given by Zhou Dugao et al. where the results of the fan's optimisation were also presented [17]. A common feature of the presented literature [1] to [3] and [17] is that it fails to take into account the impact of the GG on operating characteristics. Several papers have also been written on experimental and numerical studies of flow conditions in the AF. The paper by Yang Li et al. [4], studying the impact of the direction of the blade curvature of similar fans on IC and by five-hole probe and hot wire anemometry the impact on the local vector velocity field, represents the basis for the experimental procedures used in the presented paper. Studies of local kinematics in the flow field of the AF are also discussed in papers [5] with hot wire anemometry [6], with laser Doppler anemometry [9] and [10], computer visualisation and [7] and [8] numeric modelling which reveals the extensive research in this field. The significance of the GG's impact on the IC of the fan is also evident from the number of patents related to its design [12] to [14]. The question arises of the characteristics of the impact of the GG on the IC of the AF and on the flow kinematics in the flow tract of the fan. Further on, the paper focuses on the GG's impact through a numeric and experimental approach and presents the feasibility of evaluating the impact of the GG on the IC of the fan. 0.1 Description of the Fan with and without the Guard Grill The analysis of the GG's impact on the fan's characteristics was performed on the AF with seven unprofiled blades curved forward. The basic geometrical characteristics of the fan were given by: diameter of the outlet Dmax = 510 mm and the dimensionless hub to tip Rpe/Rob = 0.28. The fan with the GG is presented in Fig. 1 with characteristic ratia adjusted to the maximum radius of the fan's rotor: distance of the GG from the outlet area to the outlet radius s/Rmax = 0.235, with the start of the rectangular area of the GG to the axis of rotation at s/Rmax = 0.898 and the end *Corr. Author's Address: Institute Hidria Klima d.o.o., Godovic 150, 5275 Godovic, Slovenia, gregor.alic@hidria.com of the rectangular area of the GG to the axis of rotation at s/Rmax = 0.592. A fan with no GG is fixed on the outlet via the supporting construction with four legs as presented in Fig. 2. The deviation of four supporting consoles from the level of fixing the rotor is significantly bigger than the distance of the GG which permits an assessment of the hypothesis that the supporting consoles in Figure 2 do not affect the flow conditions at the outlet from the rotor. Fig. 1. Fan with the guard grill Fig. 2. Fixing of the fan to the outlet without guard grill An experimental comparative study was performed on the presented versions of the AF with regard to the GG's impact on the IC and an experimental and numeric study with regard to the GG's impact on the velocity field at the fan's outlet in the same marginal working conditions. 0.2 Measuring the Integral Characteristics Integral measurements of the fan with and without the GG were performed at the fan measurement station at the Hidria Institute Klima according to standard ISO 5801, Fig. 3. The operating characteristics of both fans given as the dependence of pressure difference from the volume airflow on the fan were determined in the entire working range. The volume flow of fans qV was measured on measurement nozzles with a sensor of differential pressure and static pressure in front of the nozzle. The Mensor 6100 sensor was used in both cases. The static pressure difference Ap at the fan was measured with the same sensor. The Vaisala HMT 330 measuring device for temperature and relative humidity was used to determine the ambient conditions and pertaining air density, while the barometric pressure was measured with the Vaisala PTB 220 measuring device. The connecting electrical power Pel and the rotating frequency of the fan's rotor were measured in addition to the aerodynamic parameters. Electrical power was measured with the Zimmer LGM450 digital power analyser and the rotating frequency n was measured with an HIK IJ measuring device. The efficiency of the fan is given by the Equation (1) which, together with the aerodynamic characteristic Ap(qv, Pei), represents the IC of the fan: A p q v n = Pe (1) The measurement uncertainty for the volume flow is ±2%, for differential pressure ±1.4% and for electrical power 1% of the measurement range. As the measurement procedures and measuring equipment used did not change for either of the two tested fans, we can assume that the measurement uncertainty was the same in both cases. 0.3 Description of Measuring the Local Velocity In order to study the GG's impact on the flow conditions, we later conducted measurements of the outflow from the fan. Measurements of local velocity were performed using a five-hole probe which is primarily suitable for determining the intensity and direction of a velocity field; however, it is unsuitable for a reliable assessment of the fluctuation characteristics of the velocity field which was not expected in the analysis. Fig. 3. Test station for fans made occording to standard ISO 5801 and its section view [11] 1 - axuliry fan, 2 - flow straightener, 3- seattling mesh, 4-testedfan Assumptions for the measurement were that the fluid flow is steady for each set working point and that it is symmetrical to the axis of the rotor's rotation. To meet the condition of axial symmetry, the experiment was conducted in a coaxial channel allowing an unhindered entry and exit velocity field with no swirl flow structures. We provided the on-flow length of 3xRmax, where R is the radius of the entry cylindrical channel, and the free exit of the airflow of the fan in the semi-area with a radius of 4xRmax, in line with the requirements of the ISO 5801 standard. The outline of the station is presented in Figure 4. Velocity values and the pertaining velocity components were calculated by using ratia of the measured pressure on the five-hole probe, the measured angle a and calibration algorithms of the probe. For this paper, the algorithms were provided by the manufacturer of the five-hole probe, United Sensor. Figure 5 presents the coordinate system selected for the experiment and enabling the appropriate presentation of the results. The differences in static pressure on the individual pressure connections of the probe were acquired simultaneously. PVC tubes of a 5 mm diameter were used to guide the pressure signals to the pressure sensors of the Endres + Hauser PD 235 type. The sampling time interval was identical for all measured points at 60s. The sampling frequency by channel was 100 Hz with 6000 samples per channel. Measuring of each working point was repeated until a dissipation of measured values has fallen under 1.5% of an average measured pressure. The ambient conditions for the measurements determined by the temperature of the environment, barometric pressure and air humidity were measured before each measurement point. computer positioning table diferential pressure transmiters for five-holes probe pressure transmiter for static pressure measuring area at fan exit settling chamber Fig. 4. Outline of the station for measuring local characteristics with a five-hole probe Vr,d + Vr,d - air flow direction Vi, + Vt_ - Vt,n + Fig. 5. Direction of axial vaks, tangential vtan and radial vrad component of velocity vector The measuring of the ambient conditions was undertaken in an identical manner and with identical equipment as measurements of the IC. The overall measurement uncertainty in measuring the velocity by the five-hole probe [15] and [16] depended on the measurement uncertainty of the differential pressure and on average equalled ±3% of the measured value at velocities above 10 m/s [16]. The relative error was increasing with the reduction of velocity, as expected. Given the fact that areas with a smaller velocity were less important for the analysis and that the trend in changes of the velocity field is maintained, we can ignore the increased measurement uncertainty in that area. The uncertainty of setting the five-hole probe is limited by the uncertainty of the positioning table and setting of the initial point ±1.5 mm. Measurements of local characteristics were made in three 'stable' working points of the fan on the outlet side of the fan (points A, B and C in the diagram in Figure 8). They were selected on the basis of a preliminary analysis of the integral measurements of the properties of both fans. As the working points at the selected number of revolutions of the fan's rotor were given in couples (qV, Ap), the selection of the working point for local measurements was already determined by the measurement of the pressure difference on the fan to simplify the experiment. Measurements of the local characteristics of the fan with and without the GG were made 7 mm behind the fan's orifice, Figure 6. Ten measurements were made within the radii ratios r/Rmax from 0.412 to 0.941 with the step Ar/Rmax 0.059. 0.4 Numeric Model Description The numeric analysis included a definition of the geometric model, setting up the topological model, Figure 7 and a calculation of the flow field by using the SST (Shear Stress Transport) turbulent model of the software CFX 10.0 with identical marginal conditions as in the measurements. SST turbulence model used Reynolds averaged Navier-Stokes equations for a steady conditions: - mass conservation equation, Eq. 2, - momentum conservation equation, Eq. 3 - along with equations of k-ю SST turbulence model, Eq. 4 and 5 are forming closed system of Eqs [8]: d( u ) S{pujui ) dp d ■— +- dx, dx1 dx, j 1 J = 0 duL +duL dx, dxi V L 1 J -put u d( pk uL д dxL dxL rk -d— +Gk - Yk L _ (2) (3) (4) Measuring area 7mm from inlet in flow direction. Fig. 6. Area of outlet velocity measurement with measuring point position d(paui )_ д dx, dx, да dx, + G., - Ym+ Da (5) Gk is the turbulence kinetic energy source, because of velocity gradients, Ga is the source of a. rk and Га are the effective diffusivity of k and ю. Yk and Yk and Ya represent dissipation of k and ю, because of turbulence. Da is the cross-diffusion term. The topological model is based on identical assumptions as in the experiment in local measurements of the flow field: the fluid x flow is steady for each calculation point and the fluid flow is axially symmetrical with regard to the axis rotation. Air density was 1.2 kg/m3 and the dynamic viscosity was 1.79x10-5 Pas. We used those assumptions to establish a topological model with 1/7 segment of the geometric model [8] with a hexagonal grid for the fan with and without the GG, Figure 7, with the specifications presented in Table 1. The assumption in setting the topological grid on the fan with the GG was that the concentric protective wire covering roughly 35% of the available outflow area of the outlet has a much greater impact than the radial supporters of the fan, covering roughly 3% of the available outflow area of the outlet. An identical assumption was made for the fan without the GG, where the radial supporters were further distanced from the outflow outlet by the distance sno/Rmax = 0.76, which is characteristically more than in the case of a GG. Marginal conditions of the numeric simulation were determined on the basis of the preliminary measurement of the fan's local velocities with and without the GG. The convergence criterion was set with regard to the development of pressure and mass flow. The convergence criterion was met with a decrease in the average remains of each quantity in the entire calculation domain by 1x 10-6, which required at least 350 iteration steps. 1 MEASUREMENTS AND NUMERIC SIMULATION RESULTS Below is a presentation of the results of measuring the integral and local characteristics of the fan and results of the numeric model for the fan with and without the GG. A comparison of the results of the experiment and the model is possible on those diagrams where integral and local characteristics were given separately. As regards IC, the analysis of the GG's impact on the fan's characteristics is primarily directed at assessing the GG's impact on the operating characteristic Ap(qV) and efficiency j while at the local level the analysis is directed at the GG's impact on the velocity distribution at the exit of the fan's airflow. 1.1 Integral Characteristic of the Fan The diagrams in Figs. 8 and 9 present the measured IC of the fan. Values were set to constant a rotational speed of n = 1392 rpm and the air density of p = 1.2 kg/m3 so that the fans' characteristics with and without the GG are comparable. In identical set conditions the diagram also gives the modelled values (MV) for both version of the fans. Fig. 7. Overall and detailed 3D calculation model of fan with the guard grill Table 1. Topological data of the calculation area offan with and without the guard gril. Fan with the guard grill Fan without guard grill Number of nodes 832512 527760 240 ■ 180 150 90 30 ar * * . -и.. * * * * 0' 4 \ % Ф Sjj/Д . - - *i i- • % a * +- я % * % % t ß C * » ■ \ % % f? B % » « • • v . д t A V \ t v\. \»T\ 0,9 0,8 ; 0,7 ■ 0,6 0,5 0,4 0,3 0,2 0,1 0,0 10000 2000 2800 3600 4400 5200 6000 6800 7600 volume flowrate qV [m /h] -test without GG + CFD with GG test with GG 8400 9200 CFD without GG Fig. 8. Comparison of the integral characteristic differential pressure-flow offan with and without the guard grill 1200 1130 1060 990 „ 920 3 I 850 О Ci S. 780 С 710 640 570 500 □ * и. * * i » . * S." 1. * •a«. » * . * i C » » * * i Л t 2000 2800 3600 4400 5200 6000 6800 7600 volume flowrate q v [m3/h] -B-test without GG -И-test with GG 8400 9200 10000 Fig. 9. Comparison of the integral characteristic electrical power-flow offan with and without the guard grill 60 0 The abscissa presents the volume airflow qV, and the ordinates present the static pressure difference Ap and the efficiency of the fan which is set for the maximum efficiency of the machine. The functional dependence of the static pressure difference from the flow Ap(qV) (the thin full curve) presented in Figure 8 and the dependence of efficiency from the flow rj(qV) (the thin dotted curve) enable an assessment of the GG's impact. Figure 5 presents the fan's IC -consumed power-flow Pel(qV) (the dotted curve). Both figures present the characteristics of the fan with the GG marked with full squares and characteristics of the fan without the GG marked with empty squares. The modelled values of the IC of the fans Ap(qV) were given by three points for each version on the fan in Figure 8. The mark (-) represents the points for the version with the GG while the mark (+) represents points for the version without the GG. In Fig. 8 we see that the fan with the GG produces a smaller flow at the same static pressure difference because the net represents a resistance to airflow. From the maximum flow with the difference of 4.2%, the difference of flows is decreasing with the flow reduction as the flow pressure loss on the protective grid is falling, except for the medium area where the difference increases. Further, a distinctive difference can be established when comparing the efficiencies where the difference is significantly bigger in higher flows. The increase in deviations in the area of the optimal working point and its movement in the direction of the increased airflow volume in the case of the fan without the GG is the result of the impact of the GG on the kinematics of the flow from the fan. The results of numeric calculations of the IC with and without the GG in the three operating areas marked in the diagram are virtually the same. In both cases, the calculated points are located very near the measured characteristic of the fan without the GG. This leads to the conclusion that the GG's impact is inadequately evaluated in integral results of the numeric model, and that the model does not include all the flow characteristics on the coaxial supporters of the net. Fig. 8 presents the dependence of the connecting electrical power on the airflow volume on the fan. Flow Pel(qV) is related to the integral aerodynamic characteristic Ap(qV) and r(qV), Figure 9. The measurement results point to the division of the entire operating area in four parts. The first part between 2000 and 3600 m3/h, where the efficiency of both fans is very similar with a difference in power consumption which increases with an increased flow and is smaller in the fan with the net as it provides for smaller flows. The second part between 3600 and 6000 m3/h is where the difference between both efficiencies as well as the flow begin to increase. However, the difference between the power consumption levels starts to fall. The third part between 6000 and 7600 m3/h is where the difference between the efficiencies partly falls, yet both fans consume similar power. And the last, fourth part between 7600 m3/h and qV,max is where the difference between efficiency again increases but the difference between power consumption at an identical flow decreases. This division was used as the basis for selecting three working points of the fan pertaining to individual parts. We expect to find the described characteristics to be the result of the flow kinematics and distinctly differing structures of the velocity field at the fan's outlet to be present at the selected points. Figure 8 presents the selected points A, B and C in which measurements of the velocity field were performed as described in Section 2.2. 1.2 Local Characteristics of the Flow Field Diagrams Figs. 10 to 15 present the results of the local characteristics. They are given in the form of a velocity distribution in the meridian plane in dependence on the selected working point (A, B and C; Fig. 8) and the position of the five-hole probe R in the radial direction from the hub to the external housing of the fan. Presented are the measured values (the dotted line) and numerically modelled values of the velocity (the full line) of the fans with and without the GG for individual working points. Absolute velocity curves are marked as full circles, the axial components of the velocity vector with full triangles, the tangential components of the velocity vector with crosses and the radial components of the velocity vector with empty circles, depending on the selected operating point. Marks in the charts Figs. 10 to 12: - Va CFD bm: Numerically calculated axial velocity of the fan without the GG; - Va mer bm: Measured axial velocity of the fan without the GG; - Vt CFD bm: Numerically calculated tangential velocity of the fan without the GG; - Vt mer bm: Measured tangential velocity of the fan without the GG; - Vr CFD bm: Numerically calculated radial velocity of the fan without the GG; - Vr mer bm: Measured axial velocity of the fan without the GG; - Vabs CFD bm: Numerically calculated absolute velocity of the fan without the GG and Vabs mer bm: Measured absolute velocity of the fan without the GG; and for the charts in Figure 13 to Figure 15: Va CFD zm: Numerically calculated axial velocity of the fan with the GG; Va mer zm: Measured axial velocity of the fan with the GG; Vt CFD zm: Numerically calculated tangential velocity of the fan with the GG; Vt mer zm: Measured tangential velocity of the fan with the GG; Vr CFD zm: Numerically calculated radial velocity of the fan with the GG; Vr mer zm: Measured axial velocity of the fan with the GG; Vabs CFD zm: Numerically calculated absolute velocity of the fan with the GG and Vabs mer zm: Measured absolute velocity of the fan with the GG. 20,0 18,0 16,0 14,0 12,0 10,0 (Я E 8,0 > 6,0 4,0 2,0 0,0 -2,0 __. чГ— -A-V ..r^c \ s \ ,0 G... G- —e- -в- —e- —e- Va CFD bm - A— Va mer bm Vt CFD bm Vt mer bm Vr CFD bm Vr mer bm -Vabs CFD bm ■ • - - Vabs mer bm 100 115 130 145 160 175 190 205 220 235 250 R [ mm ] Fig. 10. Local velocities offan without the guard grill at operating point The comparison between the measured and calculated components of the velocity vector at working point A shows an excellent match between the results as regards the fans with and without the GG, Figs. 10 and 13. In both cases, the axial component of the velocity vector increases linearly with the radius - position of the measuring probe, equally in both versions (with and without the GG) its value is reduced at the circumference, which is the result of the impact of the boundary layer and impact of the gap between the blade and outlet. Va mer bm Vt mer bm - O - Vr mer bm -Vabs CFD bm - • - - Vabs mer bm 100 115 130 145 160 175 190 205 220 235 250 R [ mm ] Fig. 11. Local velocities of fan without the guard grill at operating point 16 14 12,0 ] 10,0 W E Va mer bm - X - - Vt mer bm - О - - Vr mer bm Vabs CFD bm -Vabs mer bm 100 115 130 145 160 175 190 205 220 235 250 R [ mm ] Fig.12. Local velocities offan without the guard grill at operating point C Va CFD bm Vt CFD bm Vr CFD bm Va CFD bm Vt CFD bm Vr CFD bm The differences between the measured absolute and modelled velocity in the area near the outlet of the fan are bigger in the case of the free outflow without the GG. These differences can be linked to deviations of measured values of the radial component of velocity for both versions of 20,0 -, 18,0 16,0 14,0 12,0 „ 10,0 w 8,0 > the fans and the pertaining modelled values. In the area around the housing there is a trend of an increase in the radial component of the velocity, which is more prominent in the measured values, leading to the conclusion that the flow is directed diagonally in the peripheral area. 4,0 2,0 0,0 -2,0 ■ * ■ - Va mer zm - X - - Vt mer zm - O - Vr mer zm -Vabs CFD zm - • - - Vabs mer zm 100 115 130 145 160 175 190 R [ mm ] 205 220 235 250 Fig. 13. Local velocities offan with the guard grill at operating point A 20,0 18 16,0 14,0 12 ] 10,0 w E Va mer zm Vt mer zm - - О - - Vr mer zm Vabs CFD zm - • - - Vabs mer zm 100 115 130 145 160 175 190 205 220 235 250 R [ mm ] Fig. 14. Local velocities offan with the guard grill at operating point B Va CFD zm Vt CFD zm Vr CFD zm Vt CFD zm E > 26,0 24,0 22,0 20,0 18,0 16,0 8,0 6,0 4,0 2,0 -A---.. .A..... -L-'f. ^- Ж 'x' A" / -e---- -e---- / —e-- —© * / / - - X - - Vt mer zm -Vabs CFD zm - - • - - Vabs mer zm 160 175 190 R [ mm ] Fig. 15. Local velocities offan with the guard grill at operating point C It is a familiar phenomenon resulting from the swirl of the flow and the related centrifugal forces [1] and [2]. At point A there is a minimal difference between the components of the tangential velocity for both versions of the fans and a relatively good match with the results of the numeric model. There is an increase in the tangential component of the velocity, which is line with the expectations [1] and [2]. The differences between both versions increase at working point B, as evident from Figs. 11 and 14. By comparing the absolute velocities we can find that the trend of increasing velocity accelerates with radius R in the entire measurement area. The distinct reduction in velocity on the fan's circumference does not appear in measured values. However, it is present on the curve presenting the modelled velocity values. At this working point and as opposed to working point A there is bigger difference between the measured and modelled axial component of velocity. The match between the measured and modelled values is also poorer in the tangential component of velocity, while in both versions of the fan there is a distinct increase in the tangential component of the velocity vector, which is in line with the literature [1] and [2] and a decrease in the axial component of the velocity vector resulting from the smaller flow. A comparison with the results for working point A shows a significant increase in the radial component of the velocity vector at the circumference which is the result of an increase in the tangential component of the velocity vector. The trend of growing differences between the modelled and measured values increases with transition to working point C. The measurement uncertainty was greater at this point due to the smaller velocities and unregulated flow which should be taken into account in the analysis of results. For both fans, the numeric calculation showed the occurrence of a backflow in the form of a negative component of the axial velocity on the circumference which was confirmed by the measurements of the fan with the GG, Fig. 15. In comparison with the previous two working points, there is a strong increase in the tangential component of the velocity vector due to the bigger produced pressure, which is greater in the fan with the GG as the airflow has less interaction with the surroundings because of the additional GG. A significant increase in comparison with the previous measurements is also in the radial component of the velocity resulting from the increase in the tangential velocity. The analysis of local velocities together with the integral measurements yielded the Va mer zm Vt CFD zm 14,0 12,0 10,0 Vr mer zm 2,0 100 145 220 235 250 following results for the individual operating points. In the fourth part, the GG primarily affects the size of the velocity field while no change in direction occurs. This is reflected in the IC with an almost parallel efficiency-flow curve and pressure-flow curve. The almost identical direction of the velocity field is presented in Figs. 10 and 13. As regards the fan with the GG, the trend in velocity is the same as in the fan without the GG, only it is significantly smaller. The GG affects both the size and direction of the velocity field at the third operating point, Figs. 11 and 14. The impact on the change in the direction of the velocity field is reflected in the increasing difference between the pressure-flow and efficiency-flow curves, Fig. 8, reaching the biggest difference exactly as the part ends. Around this point, the impact of the GG on the IC starts to decrease while the impact on the direction of the flow field is substantial, Figs. 12 and 15. 2 CONCLUSIONS On the basis of the presented results, we may conclude as follows: (1) GG has a distinct impact on the aerodynamic properties of the AF. In addition to the aerodynamic resistance present due to the flow around the GG structure, impacts of the change in the velocity field should also be taken into account that form because of the effect of the geometry of the structure of the GG on the outflow from the fan's rotor. (2) By further development it is necessary to adjust the GG to the optimal on-flow velocity field from the fan without the GG to the meshed structure for decreasing the aerodynamic resistance and increasing the efficiency in the optimal working range. (3) Given the discrepancy between the model results and the measured values at the integral and local level, we can conclude that the GG has a distinct impact and should be included in the modelling of flow conditions on the fans. More attention should be paid to local phenomena in the flows around the GG. (4) Experimental measurements of the velocity field have a very good match with the numeric calculations in the operating area of the AF. However, as the difference of the static pressure of the AF increases the irregularity of the flow rises and so the uncertainty in measuring the velocity field with a five-hole probe and in the numeric calculation grows. 3 REFERENCES [1] Wallis, R. A. (1983) Axial flow fans and ducts; John Wiley & Sons. [2] Eck, B. (1973) Fans; Design and operation of centrifugal, axial-flow and cross-flow fans; Perganom press Ltd., Headington Hill Hall, Oxford. [3] Eck,B.(1978) Technische Strromungslehre; Springer Verlag; Berlin Heidelberg New York. [4] Li, Y. (2007) Experimental research on aerodynamic performance and exit flow field of low-pressure axial flow fan with circumferential skewed blades; Journal of Hydrodinamics. [5] Kergourlay, G. Koudri, S., Rankin G. W., Rey R. (2006) Eksperimental investigation of the 3D unsteady flow field of axial fans; Flow measurements and instrumentation. [6] Abdel-Rahman, A. A., Chakroun, W., Al-Fahed S. F, (1997) LDA measurements in the turbulent round jet; Mechanics research comunications. [7] Corsini, A., Rispoli F.(2005) Flow analyses in a high-presure axial fan with a non-linear eddy-viscosity closure, Heat andfluidflow. [8] Oh, K-J., Kang, S-H.(1999) A numerical investigation of th,e dual performance characteristics of a small propeller fan using viscosius flow calculations, Computers&fluids. [9] Estevadeordal, J., Gogineni, S., Copenhawer, W., Bloch, G., Brendel, M.(2000) Flow field in a low-velocity axial fan: a DPIV investigation; Experimental thermal and fluid science. [10] Yoon, J-H., Lee S-J.(2004) Stereoscopic PIV measurements of flow behind an isolated low-velocity axial fan; Experimental thermal and fluid science. [11] International organization for standardization; Industiral fans - Performance testing using standardized [15] airways, ISO 5801,(1997). [12] Nobusuke, H., Naoki, S., Shigehisa, F. (2001) Outdoor unit of air condetioner [16] JP2003172528; European patent office. [13] Toshiya, F. (2003) Fan grill and outdoor equipment for air-conditioner [17] JP2003172528; European patent office. [14] Wilhelmus, E. (2004) Outlet grill for use in an air-blowing device having axial fan W0200406839; European patent office. Chondrokostas, C. (2005) Calibration of pneumatic five-hole probes in the free-jet wind tunel; Tehnische Universitaet Wien. Zilliac, G.G. (1989) Calibration of seven-hole pressure probes for use in fluid flows with large angularity; NASA. Dugao, Z., Jing, Z., Juan, S. (1996) Optimization design of an axial flow fan used for mining local ventilation, Computers ind. Engng. Paper received: 23.01.2007 Paper accepted: 07.07.2008 Influence of The Waste Oil Concentration in Water on The Efficiency of The Aeration Process in Refinery Wastewater Treatment Milan Pavlović1, - Miroslav Stanojević2 - Mirjana Ševaljević3 - Stojan Simić4 1 Technical Faculty, University of Novi Sad, Zrenjanin, Serbia 2 Faculty of Mechanical Engineering University of Belgrade, Serbia 3 High Technical School, Zrenjanin, Serbia 4 Oil refinery Ltd., Modrica, Bosnia and Herzegovina Process, aeration system and aeration method for biological treatment of wastewater with activated sludge in bio-aeration tanks are chosen based on the flow parameters, composition of the wastewater and required characteristics of the purified water. Choosing an aeration system is a very complex question, as the capacity of oxygen introduced into the wastewater should correspond to the oxygen consumption in order to achieve the most efficient purification. This paper presents the results of experimental aeration analysis of water contaminated with waste motor oil with variations of airflow introduced into the water. The processes of clean water aeration and water with different oil content (mass concentration from 5 to 10 mg/L). The purpose of these investigations performed on an experimental setup were analysis of the aeration process technical indicators in correlation with the thermodynamic and kinetic parameters, related to the presence of oil in the water in order to achieve more efficient refinery wastewater purification. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: aeration system, refinery wastewater treatment, oxygen transport, relaxation processes 0 INTRODUCTION Almost all production lines and auxiliary installations in refineries create certain amounts of wastewater containing oil during the production process. Refinery wastewater contains large amounts of oil that is removed in tanks using scrapers during pre-treatment processes. Flocculation-sedimentation and flocculation-flotation processes are used for complete removal of free hydrocarbons from water. After that, the water is biologically purified in tanks with activated sludge or in aerated lakes and lagoons. Efficiency of aerobic biological processes in wastewater depends on the amount of dissolved oxygen introduced by aeration system. When aeration systems are designed, it is necessary to use aeration devices with an actual oxygen introduction capacity higher or equal to the actual need for oxygen in a certain phase of the treatment process of wastewater. Introduction of the necessary amount of air into bio-aeration tank and mixing of the whole suspension in the tank requires knowledge of technical characteristics of the aeration device (coefficient of oxygen transport, oxygen introduction capacity, oxygen transport efficiency and energy efficiency of oxygen transport). The coefficient of oxygen transport, kLa is the basic parameter for defining efficiency of the refinery wastewater aeration process. Some of the technical indicators of aeration process installations stated in the literature and in the technical documentation of the producers for conditions of clean water aeration are supplied. Technical characteristics of an aeration process in the wastewater were determined in an experimental installation for clean water and water with 5 and 10 mg/L motor oil concentration and presented in this paper. Also the correlation of oxygen transport coefficient with thermodynamic and kinetic parameters related to the presence of oil in the water and air flow velocity were examined in order to achieve more efficient refinery wastewater purification. *Corr. Author's Address: University of Novi Sad, Technical Faculty "Mihajlo Pupin", Đure Đakovića bb, 23000 Zrenjanin, Serbia, mstanojevic@mas.bg.ac.yu 1 EXPERIMENTAL INSTALLATION AND DESCRIPTION OF THE INVESTIGATED PROCESS The actual coefficient of oxygen transport was determined in an experimental installation by enriching water with oxygen. Dissolved oxygen was first chemically extracted from water and then waste oil with a defined mass concentration, was introduced in the water. A polypropylene column with dimensions 700x700x2200 mm with accompanying connections and construction frame, was used for experimental work in batch conditions. The cross-section surface of the column was determined according to recommendations for an air distributor. Experimental work was performed for batch working conditions and varying air flow of 2 and 10 m3/h. The water level in the column was 2m high and total volume was 980 liters. Water aeration with waste oil content of 5 and 10 mg/L was performed. Dissolved oxygen was previously removed using a chemical method of introducing sodium-sulfate in the presence of a cobalt-chloride hex-hydrate catalyst. Investigation of clean water aeration not containing dissolved oxygen was also performed, in order to compare process parameters for air distribution in standard research conditions with the ones obtained in actual conditions for wastewater with corresponding characteristics. Fig. 1 shows the experimental installation for investigating the influence of oil in water on the efficiency of the aeration process for different air flow values [1]. Table 1 contains a list of measured values and instruments on the described installation. The column was initially filled with previously prepared water from which oxygen was extracted using a chemical method. After that, a determined amount of waste oil was added into the water. Before experiment started, kinematical viscosity, density and surface tension of water in the experimental installation, had to be determined. Figure 1. Diagram of experimental installation 1 - low pressure compressor; 2 -air inflow pipe valve; 3 - relief valve; 4 - air flow regulator; 5 - air flow measuring orifice plate; 6 - column with corresponding connections and framework; 7 - disk-shaped membrane air distributor; 8 - water supply; 9 - sampling connection Complete experimental system for defining process parameters of water aeration with certain characteristics, started by measuring the temperature of the surrounding air and water in the column. The over-pressure value before the distributor (pm)d and the orifice plate i.e. the pressure difference before and after the orifice plate (App), was measured, when the first air bubbles entered in the water. Air flow regulation was performed using a flow regulator and relief valve until a set value for the applied analysis system was attained. When the flow is stabilized, water samples were taken from the column in equal time intervals starts (At = 60 s) and the dissolved oxygen content was measured until the same value appeared three times. After one regime was analyzed, the compressor was switched off and the relief valve was opened completely. Water from the column was released into drains via a draining valve. The column was then filled with a fresh amount of water. Thus, the installation was ready for a new investigation regime, i.e. the described procedure was repeated. 2 DETERMINING TECHNICAL INDICATORS OF THE AERATION PROCESS Oxygen transport coefficient kLa is a parameter used to determine the transport intensity of oxygen in water, i.e. the value at which the equilibrium state is reached. Oxygen transport coefficient is calculated as the product of the oxygen transport coefficient in water kL and the specific surface of the contact between air and water in the aeration process (a). Based on the material balance of oxygen flow through the batch reactor with complete mixing the general equation of material balance is: ■ ■ Ac Qcn + V G Cul - VlR (t ) = Qc + V G Cz +Vl — , where: VG - air flow, m3/s, Q - water flow (for batch process conditions Q = 0), m3/s, Table 1. Description of measuring equipment (list of measured values and instruments) Measured value Denotation Unit Instrument 1. Air temperature tG oC Mercury thermometer 2. Air over-pressure (pm)p mm Hg U-pipe with mercury 3. Pressure loss through the orifice plate App mm Alc Micro manometer with alcohol 4. Air over-pressure in front of the air distributor (pm)d mm Hg U-pipe with mercury 5. Water temperature in the column tL oC Mercury thermometer 6. Kinematical viscosity of water in the column OL m2/s Modified Oswald viscosimeter 7. Water density in the column Pl kg/m3 Laboratory pycnometer 8. Water surface tension OL N/m Microbalance 9. Starting mass concentration of dissolved oxygen in water co mg/L "HANNA INSTRUMENTS 9142" 10. Mass concentration of dissolved oxygen in water c(t) mg/L cin - mass concentration of oxygen in the influent, kg/m3, c - mass concentration of oxygen in the effluent, kg/m3, cul - mass concentration of oxygen in air at the input, kg/m3, ciz - mass concentration of oxygen in air at the output, kg/m3, R()) - specific oxygen consumption during biological treatment, kg/(m3s), Vl - water volume, m3. The gas phase material balance gives: VG (cul - c,z ) = AKL [c* (ciz )- c] , (2) where: c — c (ciz ) - equilibrium mass concentration of oxygen in dependence on the mass concentration of oxygen in air at the output, kg/m3, A = aVL- total contact surface between air and water, m2, a - specific surface of contact between air and water, m2/m3. The basic resistance to oxygen transport through the system occurs in a liquid hence the total transport coefficient of oxygen by water (KL) is approximately equal to the oxygen transport coefficient through water (kL). Mass concentration of oxygen in air at the output, based on modified Henry's constant (Hc), is determined as follows: ciz=Hcc (ciz). (3) Modified Henry's constant (Hc), is determined by equation: Hc = Ha/(RUTGCL), (4) where: Ha- Henry's constant, Pa kmol (02+L)/kmol 02, Ru - the universal gas constant, J/(kmolK), TG - absolute air temperature, K, CL - molar concentration of oxygen in water, kmol/m3. Equations (3) and (4) give the expression for calculating the output mass concentration of oxygen in air: Ha(Vg cul + VLkLac) VG Ha + Vlkba(R^TgCJ where: kLa - coefficient of oxygen transport, 1/s. Time dependence of dissolved oxygen concentration is presented by the following equation: Y=mX, (6) where: Ha Y — In pul M02 -2---co RuTG RuT0KCL pul M0 Ha RuTG Ru T0K CL (7) m - slope coefficient of the equilibrium curve, Ha . — k a i L --R kLa RT°K Cl KG (8) m — - RuT0K CL Ha V X = ). (9) The value of Y is given for each value of c()) determined by experiments for defined time intervals. The previous equations result in formulae for determining the coefficient of oxygen in wastewater as: kL a —- Ha (m + R ()))) Fg Ha -(m + R()))ВД*Cl (Ш) The coefficient of oxygen transport for standard conditions is determined in the following way [2]: n л (kL a)tL (kLa)s , (5) (11) where: (kLa)^ - experimentally obtained oxygen transport coefficient, 1/s, tL - water temperature, oC , в —1.024 - temperature correction factor. Based on the known value for the coefficient of oxygen transport (kLa)s the following technical characteristics of aeration systems for aeration of wastewater are calculated [3]: actual capacity of oxygen introduction OC', actual efficiency of oxygen transport E' and actual energy efficiency of oxygen transport Ee. Actual capacity of oxygen introduction is calculated as the product of the standard capacity of oxygen introduction and corresponding correction factors that convert standard investigation conditions to real ones: OC = aOC ßC* — Co 6(tL—20), (12) where: OC - standard capacity of oxygen introduction into wastewater, kg/h, OC = (kLa)s c*VL, Cs* - equilibrium mass concentration of dissolved oxygen in clean water, for standard conditions, kg/m3, a = 0.8 0.94 - relative transport degree of oxygen in wastewater, ß = 0.90 0.97 - relative saturation degree of wastewater with oxygen, ch - equilibrium mass concentration of dissolved oxygen in clean water in actual conditions, corrected in relation to the height of the liquid column above the air distributor and molar participation of oxygen in air, mg/L, 1 + Pl g Pn (H — h ) Pn - pressure corresponding to standard conditions, Pa, H - total height of the water column, m, h - height from the bottom of the tank to the distributor, m, PL - water density, kg/m3. The actual capacity of oxygen introduction OC' into wastewater should correspond to oxygen consumption during biological treatment. The actual efficiency of the transport system is expressed in percentages and represents the ratio between the real capacity of oxygen introduction and the total oxygen flow brought by the aeration system [4]: E = OC OC Go2 I vG I Pg уо. where: GO2 - mass flow of oxygen introduced into the water by the aeration system, kg/h, Va volume air flow (forpn=101.3 kPa, tL = 20oC), m3/h, , pG - air density for standard conditions (tG = 20oC), kg/m3, Уо = 0.232 - mass participation of oxygen in air, kg/kg. The real energy efficiency of oxygen transport represents the ratio between the actual capacity of oxygen introduction and the applied power needed to drive the aeration device: E. = OC X p (14) where: X P - sum of power off all electrical motors (aerator, pump, blower, etc.), kW. Based on all the experiments performed on the described installation for varying air flow between 2 and 10 m3/h and waste motor oil content from 0 to 10 mg/L using the investigated model (kLa)s values given in Table 2 were determined. The values obtained for (kLa)s were used to determine the values of actual capacity of oxygen introduction OC', real oxygen transport efficiency E' and actual energy efficiency of oxygen transport Ee also given in Table 2. Technical characteristics of the aeration system in dependence on oil presence in wastewater, are given as diagrams in Fig. 2. An increase of polluting components' concentrations (oil) in the water, led to a reduction of oxygen transport coefficient values and affecting other technical characteristics of the investigated aeration system. Strojniški vestnik - Journal of Mechanical Engineering 54(2008)10, 675-684 Table 2. Technical characteristics of the aeration system for corresponding experimental regimes Waste oil concentration in Regime water ( VG)n (kL a)s OC' E' Ee number mg/L m3/h 1/h g/h % g/kWh 1 - 2.361 2.887 11.754 1.465 35.078 2 - 5.843 3.052 12.425 0.617 15.024 3 - 10.768 5.267 21.876 0.550 14.851 4 5 2.361 2.248 9.528 1.118 28.442 5 5 5.862 2.378 9.519 0.463 11.510 6 5 10.735 4.527 18.999 0.625 12.898 7 10 2.379 1.993 8.204 1.002 24.489 8 10 5.899 2.160 8.837 0.426 10.686 9 10 10.724 2.577 10.730 0.271 7.284 c) (Vc)n. m3/h \ S \ л N \ \ s * ч \ .л • s •ч. (V0)„, m3/h b) d) (Vc)n, m3/h \ \ \ \ S \ (Ve)n, m3/h ■Content in water 0 mg/l ■Content in water 10 mg/l ■Content in water 5 mg/l Figure 2. Influence of disc-shaped membrane air distributor's technical characteristics on air flow for different mass concentrations of oil in water (0; 5 and 10 mg/L) a) coefficient of oxygen transport for standard conditions (ka)s; b) actual capacity of oxygen introduction OC'; c) actual oxygen transport efficiency E'; d) actual energy efficiency of oxygen transport Ee' 30 25 20 15 10 5 0 40 1,4 35 1,2 30 1,0 25 0,8 20 0,6 0,4 0,2 3 DEPENDENCY OF TECHNICAL INDICATORS ON THERMODYNAMIC INDICATORS AND PRESENCE OF OIL IN THE WATER DURING AERATION PROCESS The law of energy balance can explain reduction of oxygen transport coefficient in the water during aeration with increasing waste oil concentration: dEtot _ dEkin dE pot dz dz dz (15) Change velocity of the total oxygen molar energy, E0 tot , is determined with oxygen transport coefficient, kLa and influenced by the diffusion rate constant, kd , the kinetic molar energy, E0kin , as well as the adsorption rate constants, ka in potential field through the free liquid surface, a [5]: 170 = kd2E kin + kaEfp kLaE tot = kd2 (16) (17) kin 1 ^a^ pot > or in the form: kLa = kd + (ka /E0tot) AadsG0 The potential molar energy, E pot determined with oxygen bond strength, is equal with the oxygen molar adsorption Gibbs energy, dads Ge, and can be affected by the presence of the oil oxidation products and radical initiators (OH-, H2O2, Fe2+, HCO3-/CO3 2-, PO4 3-, humid acids, aryl-R etc) [6]. Standard oxygen adsorption Gibbs energy in this work is determined with the Gibbs-Helmholtz - van, t Hoff equation [7] , applied for saturated water: MO2)a q = MO2)g . (18) The chemical potential, ,m(O2), is determined with the standard value, ,m8(O2) and with oxygen concentration (c(O2)): H (O2) = / (O2) + RT ln c(Ü2). (19) Determination of adsorption molar energy on the base the Equations (18) and (19) is performed using the obtained linear function (20): ln c(O2)a q = ln c(O2)g - Aads G0 (O2)/RT, (20) and presented in the Table 3 as: y = yo - tg p-x, where: y = ln c (O2)a q , yo = ln c (O2)g , tg p = Aads G0 (O2) /R, x=1/T. The calculated oxygen adsorption Gibb-s energies for the investigated regimes on the basis of the obtained linear functions to equation (21) are presented in the Table 4: Aads G0(O2) = R'g p , (21) The calculated oxygen kinetic parameters (the adsorption rate constants and diffusion rate constants) in aeration process: ka = 3RTtg p, (22) kd = yo , (23) are presented in the Fig. 3 and the Table 5. Table 3. Functions derived from the Equation (20) in conditions defined with c-h-q, c (oil concentration), m3/h Investigated regimes defined with: c-h-q y = yo - tg p-x were is to eqation (20) y= ln c(O2,aq) and x = 1/T Coefficient correlation, R2 0-2-2, 0-2-6, 0-2-10 y = -3663.1x + 4.4554 0.9219 0-2-2, 0-2-6, 5-2-10, 5-1-6 y = 4367x - 23.624 0.9751 0-2-10, 5-2-10, 10 -2-10 y = 165691x - 585.58 0.7455 5-2-2 0-2-10, 5-1-2, 5-1-10 y = 8959.3x - 39.575 0.6137 5-2-6, 10-2-2, 10-2-6 y = 7726.9x - 35.569 0.9939 10-2-10, 10-1-2, 10-1-6, 10-1-10 y = 8492.3x - 38.096 0.9393 Table 4. Thermodynamic and technical characteristics of examined regimes Regime number Waste oil mg/L ( VG)n m3/h kLa 1/min c (O2,l) mg/dm3 T K A ads G (O2) kJ/mol ^Aads G (O2) kJ/mol 1 - 2.361 0.0419 7.5 286 30.4 -36.3 -5.9 2 - 5.843 0.0443 7.6 286 30.4 -36.3 -5.9 3 - 10.768 0.0779 7.9 287 30.4 -1377 -1346 4 5 2.361 0.0339 6.4 288.0 -74.5 -74.5 5 5 5.862 0.0339 6.7 285.1 -64.2 -64.2 6 5 0.735 0.0677 7.0 287.5 -64.2 -1377 -1441.2 7 10 2.379 0.0292 5.9 286.7 -64.2 -64.2 8 10 5.899 0.0315 6.1 286.9 -64.2 -64.2 9 10 10.724 0.0382 6.4 287.1 -70.6 -1377 -1447.6 kLa, 1/min = f (Д ads G0 (O2 ) kJ/mol c (waste oil) = 0 mg/l 0,1 0,c 0,06 kLa, 1/min = f (A ads G0 (O2 ) kJ/mol c (waste oil) = 5 mg/l y = -2E-05x + 0,0322 R2 = 1 0,08 0,06 «004 -0,02 0 -2000 -1500 -1000 -500 0 kLa, 1/min = f (A ads G0 (O2) kJ/mol c (waste oil) = 10 mg/l y = -6E-06x + 0,03 R2 = 0,9395 -2000 -1500 -1000 -500 0,02 0,01 0 0 Figure 3. Diagrams obtained using Equation (17), for the regimes with the constant oil concentrations Table 5. Theoretical values of the relaxation processes rate constants in the aeration system for corresponding investigation regimes Regime number (vg )n m3/h dissG& (O2 ,g/aq) kJ/mol kLa 1/min kd 1/min k a 1/min 1/r * m 1/min Compared measured with the calculated aeration rate constants 1 2.361 -5.9 0.0419 0.043 - 0.214 0.042 1/rm = kLa = kd (1/Tm)/ka = 0.19 2 5.843 -5.9 0.0443 0.045 1/rm = kLa = kd (1/Tm)/k a = 0.21 3 10.768 1346 0.0779 0.048 1/rm /kLa = 0.62 1/rm /kd = 1.11 1/rm /ka = 0.22 4 2.361 -74.5 0.0339 0.032 - 0.143 0.04 (1/im)/kLa =1.17 (1/rm)/kd = 1.25 (1/rm)/ka= 0.28 5 5.862 -64.2 0.0339 0.043 (1/rm)/kLa = 1.27 (1/rm)/kd = 1.34 (1/rm)/ka = 0.3 6 0.735 -1441.2 0.0677 0.053 (1/rm)/kLa = 0.78 (1/rm)/kd = 1.65 (1/rm)/ka = 0.37 7 2.379 -64.2 0.0292 0.03 -0.043 0.042 (1/rm)/kLa = 1.42 (1/rm)/kd = 1.39 (1/rm)/ka = 0.97 8 5.899 -64.2 0.0315 0.032 (1/rm)/kLa = 1.4 (1/rm)/kd = 1.4 (1/rm)/ka = 1.0 9 10.724 -1447.6 0.0382 0.053 (1/rm)/kLa = 1.38 (1/rm)/kd = 1.75 (1/rm)/ka = 1.22 * measured period of aeration time enabling that constant oxygen concentration was achieved The oxygen transport coefficient as technical indicator of aeration kinetic in correlation with the rate constants, obtained on the base of the thermodynamic study of aeration kinetic as in correlation with directly measured saturation time is presented in the Table 5. In clean water, the oxygen transport coefficient is equal with the diffusion rate constant and with parameter calculated on the base saturation time (kLa = kd = 1/rm), but adsorption rate constants were five time greater, ka/(1/rm) = 4.8 -г 5.3. In the presence of added oil, relative adsorption rate constants were from 3.6 г 2.7 for 5 mg/l oil concentration up to 1.34 г 0.82 for 10 mg/L oil concentration. Equal adsorption rate constants with the saturation rate constants reduced aeration efficiency due to activated oxygen in the slow adsorption processes at 10 mg/L oil concentration. 4 CONCLUSION Experimental research shows that 5 mg/L concentration of oil in the water, leads to reduction of the dissolved oxygen concentration from 10 to 15% compared to clean water, while 10 mg/L concentration of oil in the water reduces concentration of dissolved oxygen from 20 to 25% when compared to clean water. Aeration process technical indicators (oxygen introduction capacity, oxygen transport efficiency and energy efficiency of oxygen transport) for the clean water and water with 5 mg/L oil concentrations were significantly better with the maximum air flow, 10 m3/h, than with the less flow of 2 m3/h and 6 m3/h. For an 011 concentration of 10 mg/L, the air flow of 10 m3/h did not have a significant influence on the aeration process technical indicators. Based on the research of the dependence of technical indicators on thermodynamic parameters and presence of oil in the water during aeration process it can be concluded that: 1. The oxygen transport coefficient is equal with the diffusion and saturation rate constants for clean water and contaminated water at the air flow of up to 6 m3/h. 2. The oxygen transport coefficient for 10 m3/h air flow: - is increased compared to diffusion rate constants 1.8 times and for water with 5 mg/L oil concentration 2.1 times for clean water - is only 1.2, except for water with 10 mg/l oil concentration when the relative oxygen transport coefficient was compared to diffusion rate constants 3. The study of the Gibbs energy adsorption dependence from oil concentrations shows: - oxygen adsorption is activated with parallel endothermic and exothermic oxygen processes at the air flows up to 2 m3/h and 6 m3/h in the clean water, - oxygen adsorption is activated with spontaneous, exothermic electrons transitions with the most adsorption affinity at the air flow of 10 m3/h in the contaminated water. 4. Adsorption rate constants with increased oil concentrations become less negative and diffusion rate constants less positive values of up to 10 mg/L oil concentration, when equal adsorption and saturation rate constants prevent the air flow influence on increasing the oxygen transport coefficient. 5 REFERENCES [1] Stanojević, M., Simić, S., Radić, D., Jovović, A. (2006), Wastewater aeration, ETA, pages 116, ISBN 86-85361-07-9, Belgrade, Serbia. [2] Ashley, K., Hall, K., Mavinic, D. (1991), Factors influencing oxygen transfer in fine pore diffused aeration, University of British Columbia, Water Research 25 (1991), pp. 1479-1486., Vancouver, Canada. [3] Stanojević, M., Radić, D., Simić, S. (2004), Determining the technical characteristics of the aeration systems for oil refinery's wastewater treatment, 16th International Congress of Chemical and Process Engineering, CHISA 2004., Praha, Czech Republic, 22-26.08.2004., pp. 3863-3873. [4] Wagner, M., Pöpel, J. (1998), Oxygen transfer and aeration Effciency-Influence of Diffuser submergence, Diffuser Desity and Blower Type, Water Quality International, IAWQ 19th Bienial international Conference, Vancouver, Canada. [5] Ševaljević, M.M. (2000) Ph D Thesis, Development of the method of galvanostatic arsenic hydride and successive lead and cadmium pre-concentration for the AAS determination, Faculty of Technology, Novi Sad. [6] Wall, F.T. (1965) Chemical thermodynamic, Frezenius and co., San Francisco. [7] Von Gunten, U. (2003) Water Research, Ozone reaction mechanisms, 37 (2003) pp. 1469-1487. [8] Dular M., Bachert, R., Širok, B., Stoffel, B. (2005), Transient simulation, visualization and PIV-LIF measurements of the cavitation on different hydrofoil configurations, Strojniski vestnik -Journal of Mechanical Engineering 51 (2005)1, pp.13-27. Paper received: 02.01.2008 Paper accepted: 07.07.2008 Strength of Pressure Vessels with Ellipsoidal Heads Pavo Balicević1, - Dražan Kozak2 - Tomislav Mrcela3 1 Josip Juraj Strossmayer University of Osijek, Faculty of Agriculture, Croatia 2 Josip Juraj Strossmayer University of Osijek, Mechanical Engineering Faculty in Slavonski Brod, Croatia 3 Josip Juraj Strossmayer University of Osijek, Faculty of Electrical Engineering, Croatia A method for stress analysis in cylindrical pressure vessels with ellipsoidal heads, based on the axisymmetric shell theory, was proposed. The starting point were the approximate solutions of the differential equation system that were used to get mathematical expressions for determining internal forces, moments and displacements in the vessel walls. Final expressions that can be applied were acquired by joining the membrane and moment theory and by setting and solving equations of boundary conditions. Diagrams that show distribution of internal forces and moments are in dimensionless form which enables their use regardless of dimensions and load. These expressions were used to develop a method for testing strength of pressure vessels with ellipsoidal heads in the design phase. Application of the method was shown on a selected numerical example, while a special computer programme was created for calculation purposes. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: pressure vessels, ellipsoidal heads, internal forces, analytical solutions 0 INTRODUCTION Because of the increasing use of pressure vessels in all types of industry (e.g. oil-refining industry, gas-processing industry, power plant industry, etc.), in transportation, and other processing facilities, and because of their many side effects on the environment, there is a need for accurate analysis of the shells in the design phase. The designer must know the real distribution of stresses, i.e. he must make calculations and analysis of stress in all critical parts of the vessel in order to make an adequate pressure vessel design. Stress is highest at the points of geometry parameter change, e.g. the transition from the shell of the vessel to the head, where higher wall bending usually occurs. Simple empirical forms determined by norms are used for making calculations of stress in the walls of pressure vessels. This type of calculation does not give the real distribution of stress, the highest values or places where stress occurs. More accurate stress analysis can be done by using finite element method (FEM) or analytical method. Analytical method for calculating strength of the thin-walled pressure vessels is based on the theory of axisymmetric shells. This theory has a system of the fourth order differential equation with internal forces, moments and displacements as unknown functions. The shell theory was established by famous researchers (Lourye 1947, Goldenveiser 1953, Novozhilov 1962 and others) who developed analytical methods for determining general solutions of the system. Accurate solutions of the system can be acquired in general form only for pressure vessels of simple geometry (e.g. cylindrical, conical and spherical). In order to apply those solutions in the process of designing, we must calculate numerous complex functions and infinite orders with complex variables, and we must also set and solve an equation of boundary conditions. These solutions can be used to make calculations on a computer, but the calculation procedure and the algorithm structure are very complex. Far simpler analytical procedure for calculating pressure vessel strength can be developed on the basis of approximate solutions of axisymmetric shell theory equation system [1]. It is of vital importance because the approximate solutions contain an error of the same order of magnitude as do the correct solutions. Relative error of those solutions is of the order h/R (ratio of the wall thickness and the curvature radius), which is less than 1/20 at thin-walled pressure vessels, i.e. 5%, which can be neglected in engineering calculations. Approximate solutions of the shell theory can be applied to pressure vessels with more complex shells. This paper presents solutions for a cylindrical vessel with ellipsoidal heads. Final mathematical expressions for calculating internal *Corr. Author's Address: University of Osijek, Faculty of Agriculture, Trg Svetog Trojstva 3, HR-31000 Osijek, Croatia, pavo.balicevic@pfos.hr forces, moments and displacements in the walls of head and cylindrical part were derived from approximate solutions. By using those expressions, and by connecting membrane and moment theory, a method for determining strength of pressure vessels with ellipsoidal heads, which is suitable for designing, was developed. A special computer programme was created for the application of this method. Computer calculation was done on a selected numerical example, and the analysis results were shown in a diagram. 1 MEMBRANE FORCES IN ELLIPSOIDAL HEAD Stresses in the walls of pressure vessels occur due to different types of loads, depending on the purpose of the vessel and on different influences that a vessel is subjected to during exploitation. Internal pressure has the biggest influence on the amount of stress, so all other types of loads are considered to be less important. Solutions of the shell theory equation show that internal forces which occur in the walls of the vessel can be, under certain conditions, determined by superimposing of two kinds of loads, membrane and moment. These conditions given in [2] are true for curvature radius derivations of vessel meridian R1 and surface load qn. At membrane stress state, it is assumed that only normal forces occur in the wall, while shear forces and bending moments are not considered. In the design phase, pressure vessels should be formed so that the real stress is approximately the same as the membrane state in order to avoid bending of the walls and high stresses due to bending. Cylindrical pressure vessels usually have semi-spherical, toroidal-spherical or ellipsoidal heads. Stress is best distributed in ellipsoidal heads. The middle plane of ellipsoidal head is shaped like a half of the ellipsoid of revolution which is the result of ellipse arc rotation around its minor axis. A major ellipsoid semi-axis is equal to cylindrical part radius (a = R), and a minor semi-axis is equal to the head height (b = H). A parameter determining ellipse form is у = a2 / b2 -1, and R0 = a -J 1 + y is curvature radius of head vertex. Figure 1 shows the basic measures of cylindrical pressure vessels, closed by heads shaped like revolution ellipsoid half. Fig. 1. Pressure vessel with ellipsoidal heads Internal forces in the head walls are determined according to the membrane theory [2] by using the following expressions: \1/2 N (0), = p_r ' 2 1 + y 1 + ^ sin 2 в Nt(o)' = (1 -rsin2 в) Ns(o)' (1) (2) Membrane component of circular force has (in places of shell connection) tensile character on the cylinder, and compressive character on the head, which leads to difference in radial displacements according to [3]: ur(0)'' - ur()' = P*- (1 + r), r r 2 E h У '' (3) where: Ns(0)' - meridian normal force in the head wall, n( ur(0)' - radial displacement of the wall (at the head), ur(0)'' - radial displacement of the wall (at the cylindrical part), E - modulus of elasticity of the pressure vessel wall, h - thickness of the wall, p - pressure in the vessel, в - meridian angle, R - radius of cylindrical part of the vessel. (0)' - circular normal force in the head wall, Superscript(0) denotes values pertaining to a membrane stress condition. Corresponding boundary condition values will be denoted by the superscript(l), and those without the superscript refer to total solutions. As can be seen from the expression (1), meridian normal force of ellipsoidal head changes without interruption in the points of transition from the ellipsoid to the cylinder, which is very important for achieving the membrane state of stress. Both normal forces (meridian and circular) reach their maximum value in the head vertex. Value у = 3 should be chosen for ellipse parameter, because maximum values of internal forces in the head are the same as values of circular force in the cylindrical part of the vessel. There is a discontinuity of the membrane component of circular force in the places where head and the cylinder are connected. Since total displacements of the vessel walls should stay intact, additional forces and moments occur in the vessel walls, which lead to bending of the walls. Expressions for determining their values can be acquired from general solutions of shell theory differential equation system. 2 FORCES AND MOMENTS OF BOUNDARY EFFECT Fig.2. Boundary load of the long axisymmetric shell With those conditions, general solutions of axisymmetric shell theory differential equation system, acquired by approximate solutions of the system by [4] can be reduced to: ( ) F ur' ) ---sin2 00 - e~x cos x r 2Dß3 0 M0 . n _x / .4 --- sin 00 - e ( cos x _ sin x ), 2Dß2 (5) 2Dß0 sin00 - e x (cos x + sin x) M, Dßc — e xcosx (6) Thin-walled axisymmetric shells can be divided according to the moment theory into short and long. The shells are considered long if: X- jß-ds > 3, ß- fì (71 V R2 - h (4) With long shells we can disregard the influence of the load of one end of the shell on the internal forces and displacements on the other end. Each end of the shell can be observed independently (i.e. without considering the conditions at the opposite end). Bending of the shell occurs due to radial forces F0 and bending moments M0 at the boundary of the shell, which is called boundary load. Boundary load of the long axisymmetric shell is shown in Figure 2. The condition (4) is regularly achieved for the heads and cylindrical parts of thin-walled pressure vessels. Fr(1) - F0 - e x (cos x _ sin x) 2M0 ß0 _ siné>n Ms - —- sin 00 - e x sin x s ßc 0 M0 - e~x (cos x + sin x) Mt-v-Ms Nt- — -u( N(1) - Fr(1)- cos 0 s r where: iß - ds , D - - E h3 12 (1 _ (1 _ ) (7) (8) (9) (10) x = In expressions (4) to (10) the following was used: Fr(1) - radial force in the wall at boundary load, Ns(1) - meridian normal force in the wall, Nt(1) - circular normal force in the wall, ur(1) - radial displacement of the wall, 3(1) - wall twist angle, Ms(1) - meridian moment of bending, Mt(1) - circular moment of bending, v - Poisson's ratio, s - meridian arc length from the shell boundary to a certain point, s1 - total shell length, r - radius of the middle plane in a certain point of the wall, R- circular radius of wall curvature. All expressions with a subscript 0 are values for the shell boundary. For example, 60 and ß are values of 6 and ßin points x = 0, at the shell boundary. In order to apply expressions (5) to (8) it is necessary to determine the values of the force F0 and moment M0 of boundary load. They are, by their nature, internal forces of moment theory at the joint of two different shells that superimpose forces of membrane theory, in order to achieve inner mechanical balance. Their values for the example under consideration (cylindrical part of the vessel with ellipsoidal heads) in Figure 3 will be determined by setting equations of boundary conditions. Strain continuity at the joint (i.e. place where cylinder and ellipsoid are connected) depends on equal total displacements and twist angles of joined parts, which is expressed by the following equations: ur(0)' + F0 ■Зц' — M0 ■5,2' = = ur(0)'' — Fn ■S„" — M0 ■512' —F0 ■ V + M0 ■ 522 ' = = —Fn ■S^'—M0 ■S 0 22 (11) (12) Membrane components of the wall twist angle, both for ellipsoidal and cylindrical shell, equal null at the angle в=п/2, which can be shown by using membrane theory [5]. Coefficient values of generalized forces influence (F0 and M0) on generalized strain (ur and 3 ), with long shells, according to [4] are: S11' =S11" = sin 60 2Dß03 S ' S '' sin60 Sn =Sn = - 2Dß0 S22 = S22 " = Dß0 (13) Fig.3. Forces and moments of boundary effect where: ß0 - M (3 ■JTh By solving equations system (11) and (12), and by considering values (3) we get: M о = 0, Fe = P 8 ß0 (1+r). (14) By substituting solution (14) into expressions (5) to (10), we get formulas for internal forces and displacements that occur due to boundary effect in the walls of ellipsoidal (') and cylindrical ('') part of the vessel: 1 (1)' = ERL 4 E h uvw' = --- ( 1 + y) e-x cos x (15) v 4 3 ( 1 ) where: x — ß ds , ß — —— — 0J (23) ,(1) - —- p R r 4 E h (1+r) Ms' — 8 ß 2 ( 1 + y) e x sin x, M " — —i— s „ „ 2 8 ß0' (1 + j) (1) ' — ^ (1 + r) e-x cos x , (1) '' — - pR Fr(1)' — - TT (1 + ^)e^x cos (x + ^4), 8 ß0 Fr(1) '' — pT(1 + r)e-xcos (x + nl4), 8 ß0 (22) (16) (17) (18) (19) - — ( 1 + r) e x cos x, 4 (20) (21) Distribution of internal forces and moments along the meridian length of the pressure vessel wall, according to expressions (15) to (22), is shown in the diagrams (Figures 4, 5, and 6). According to the diagram, circular force Nt (Fig.4) and meridian bending moment Ms (Fig. 5) have the biggest influence on stress. Circular moment Mt , according to (9), has smaller values, and components of radial force Fr are of lower order of magnitude so that they can be disregarded in stress calculation. 4 METH0D F0R TESTING STRENGTH Total values of internal forces in the pressure vessel walls can be acquired by adding components of membrane and moment theory: Nt — Nt (0) + N(1) Ns — N„ (0) + N, (1) (24) If we analyse expressions (17) to (22), the following can be concluded. For designing purposes values of internal forces and moments can be observed as dimensionless magnitudes, e.g. Ns/(p• R) and Ms/(p• R• h). Beside the position of the point on the vessel wall, they depend also on design parameters, i.e. on the head shape j and ratio R/h. Fig.4. Distribution of circular forces along the vessel wall determined by moment theory & лШ/3 Fig.5. Distribution of meridian moments along the vessel wall determined by moment theory Fig.6. Distribution of radial forces along the vessel wall determined by moment theory Therefore, by choosing values of these parameters we can test strength before determining final dimensions of the vessel. In order to perform this procedure successfully, a special computer programme (in Fortran 77) [6], based on mathematical expressions, was created. This programme can be used to calculate internal forces and moments by using expressions (1) and (2), (15) to (24) in a number of points distributed along the meridian of the vessel. It can be also used to determine values of main stresses on the inner and outer surface of the wall and to calculate equivalent stresses. The programme is made in such a way that calculation is performed for the chosen ratio R/h and factor of head shape y It was carried out on the example R/h = 87 and y = 3, and calculated values of internal forces and moments are shown in Table 1 only for a limited number of points. Distribution of total circular force, according to calculated values, is shown in Fig. 7. Equivalent stress according to HMH theory of strength [7] was chosen as a criterion for testing strength of pressure vessel. This programme can be used to calculate equivalent stresses at different points of the wall and to find critical points. Maximum value of equivalent stress was determined for the above mentioned example: K ) = 1. 42. v e/ max h (25) which occurs in the points at the inner wall surface at meridian angle в = 69.61° (which is equal to polar angle p = 5.31°). On the outer surface of the wall, maximum stress occurs at the head vertex, and it is much lower in the mentioned example. Stress determined by the expression (25) is the basis for the choice of pressure vessel material (concerning the necessary strength), after the shape and dimensions of the pressure vessel have been determined. These two steps in the design phase can be done iteratively, i.e. by varying design R/h and у , as input data, we can get values of maximum equivalent stresses for different dimension ratia. 5 CONCLUSION Approximate solutions of axisymmetric shell theory used for expressions (15) to (22) are valid for steep shells, i.e. those with big в angle. This condition is found in the vicinity of the place where the head and the cylinder are joined (i.e. where angle в slightly differs from 90°). It is better to use these expressions instead of exact solutions because they consist of simple mathematical functions, so that there is no need to consider boundary conditions, since they are already included in final expressions. These expressions clearly show the law of distribution of internal forces and moments in the walls by using diagrams, which was done for ellipsoidal heads with R/H=2 ratio, since they are commonly used in praxis. These diagrams can be used for analysis of stress and for choosing dimensions during the design phase of such heads. Computer programme enables quick acquisition of data and representation of such diagrams for other R/H ratia. Table 1. Values of internal forces and moments in the ellipsoidal head of cylindrical pressure vessel Cylindrical vessel with ellipsoidal heads, у = 3, Rh = 87 в x Nt Ns Mt Ms p ■ R p ■ R p ■ R ■ h p ■ R ■ h 0 11.59 1.00 1.00 0 0 л/ 32 9.950 0.9574 0.9859 0.000 0.000 л/16 8.417 0.8391 0.9474 0.000 0.000 3л/ 32 7.063 0.6682 0.8934 0.000 0.000 л/ 8 5.910 0.4699 0.8335 0.000 0.000 5л/ 32 4.947 0.2599 0.7746 0.000 0.0021 3л/16 4.145 0.0448 0.7206 0.0012 0.0040 7^ 32 3.472 -0.1689 0.6731 0.0009 0.0031 Z 4 2.901 -0.3696 0.6325 -0.0012 -0.0040 9^ 32 2.408 -0.5412 0.5984 -0.0055 -0.0182 5^ 16 1.975 -0.6672 0.5704 -0.0116 -0.0386 11л 32 1.588 -0.7338 0.5477 -0.0186 -0.0618 3л 8 1.234 -0.7309 0.530 -0.0249 -0.0831 13л 32 0.9056 -0.6531 0.5166 -0.0289 -0.0963 7 л 16 0.5946 -0.4996 0.5073 -0.0281 -0.0935 15 л 32 0.2946 -0.2764 0.5018 -0.0196 -0.0654 л 2 0 0 0.5 0 0 The described method for calculating strength, performed by a computer programme, gives data on critical stress in the form (25), in a short period of time, for arbitrarily chosen geometry parameters of ellipsoidal head, which makes this method very suitable for designing purposes. 6 REFERENCES [1] Senjanović, I. (1973) Theory of Shells of Revolution, Naval Architect Institute, Zagreb. [2] Biderman, V. L. (1977) Mechanics of Thin-Walled Structures (in Russian), Mashinostroenie, Moscow. [3] Balicević, P., Kozak, D., Kraljević, D. (2007) Analytical and Numerical Solutions of Internal Forces by Cylindrical Pressure Vessel with Semi-Elliptical Heads, Proceedings of 1st International Congress of Serbian Society of Mechanics, Kopaonik, April 10-13, 2007. [4] Alfirević, I. (2003) Linear analysis of structures (in Croatian), Faculty of Mechanical Engineering and Naval Architecture, Zagreb. [5] Balicević, P. (2006) About Reliability and Terotechnology in Design of Pressure Vessels, Doctoral dissertation (in Croatian), Mechanical Engineering Faculty, Slavonski Brod. [6] Dovedan, Z., Smilevski, M., Zalokar, J.D. (1987) Fortran 77 with programming technics (in Slovenian), Association for technical culture of Slovenia, Ljubljana. [7] Fryer, D.M., Harvey, J.F. (1998) High Pressure Vessels, Chapman & Hal, International Thomson Publishing, New York. Paper received: 14.06.2007 Paper accepted: 07.07.2008 Numerical and Experimental Analysis of a Shaft Bow Influence on a Rotor to Stator Contact Dynamics Sanjin Braut 1 - Roberto Žigulić1 - Mirko Butković2 1 University of Rijeka, Faculty of Engineering, Croatia 2 Polytechnic of Karlovac, Croatia The shaft bow problem presents a real situation especially in case of slender rotors. This paper investigates the shaft bow influence on the rotor-stator contact dynamics. For this purpose the rotor is described as a simple Jeffcott model and the stator as an elastically suspended rigid ring. To test the numerical model, except a usual run down analysis, an emergency shut down after the sudden rotor unbalance increase is also analyzed. Numerical integration is carried out by the fourth-order Runge-Kutta method. Two different normal force models for the rotor-stator interaction are analyzed. For analyzed parameters, both linear and nonlinear (Hunt and Crossley), normal force models gave similar rotor and stator responses. To confirm some of the results and to tune the numerical model, the experimental investigation on the test rig was conducted. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: rotor-stator contact, shaft bow, numerical simulations, rotor vibration, vibration measurements 0 INTRODUCTION Contact or rub between rotor and stator is one of the most intensive research subjects in rotor dynamics. Most papers which are considering the rotor-stator contact phenomena can be classified in rigid rotor disc-rigid stator contact [1] to [3], bladed disc-stator contacts [4] and rotor-stator contact in retainer bearings when the rotor is supported by active magnetic bearings [5] and [6]. Choy and Padovan [1] developed a general analytical rub model using the following assumptions: simple Jeffcott rotor model, linear stiffness and damping characteristics, rigid casing supported by springs acting in the radial direction, mass inertia of the casing small enough to be neglected, simple Coulomb friction and onset of rub caused by unbalance. Bartha [2] performed an extensive numerical and experimental research of the backward whirl of rotors considering the rigid and elastically suspended stator. Von Grol and Ewins compared measurements and simulations for a windmilling imbalance in aero-engines which is very similar to the classical rotor to stator contact. They have revealed the rotor response rich in different subharmonics. Influence of torsion on the rotor-stator contact was introduced by Edwards et al. [7]. Attention was paid to the effects of torsion on a steady state response of a system experiencing rotor to stator contact. The analyzed Jeffcott rotor model had only 3 degrees of freedom - d.o.f. because the stator was assumed to be rigid. Žigulić et al. [8] considered nonlinear dynamics of a rotor supported by two dry-friction bearings. The rotor was modeled with FEM-Finite element method and the obtained nonlinear system was integrated with Hilbert-Hughes-Taylor - HHT a method. Karpenko et al. [9] have shown how the preloading of the snubber ring could stabilize the dynamic response of the rotor. Theoretical predictions of the two-degrees-of-freedom Jeffcott rotor model with the preloaded snubber ring subjected to out-of-balance excitation was studied in [10], while an experimental verification with special attention to the analysis of the shaft rotational speed (excitation frequency) as well as different eccentricities and preloading was given in [11]. This paper considers the contact between the rotor, with a bowed shaft, and the stator at the position of rotor disc (seal), after the appearance of sudden unbalance due to blade loss. Linear and nonlinear models of the normal contact force are compared for different shaft bow - mass unbalance combinations. In this study, not only lateral but also torsional d.o.f. are taken into account allowing the additional extension of the *Corr. Author's Address: University of Rijeka, Faculty of Engineering, Vukovarska 58, HR-51000, Rijeka, Croatia, sanjin.braut@riteh.hr model with mechanical model of the induction motor. Results of numerical simulations and experiments are presented and compared using the rotor and stator responses in time domain, spectral maps and rotor orbits 1 MODELING OF THE ROTOR-STATOR SYSTEM A rather simple mechanical model for the rotor-stator interaction description is used. The rotor is described by Jeffcott model with a shaft bow while the stator is modeled as an elastically suspended rigid ring. Except traditional consideration of the rotor lateral d.o.f. (xr, yr), torsional d.o.f. (щ) is taken into account allowing an extension of the model with torsional d.o.f. of the induction motor (щ) as well as appropriate mechanical model. To keep consistency of the described consideration the stator had additional three d.o.f., two lateral (xs, ys) and one torsional (щ). Fig. 1 shows the considered rotor-stator system. Next sections briefly discuss the differential equations of motion, the contact model and mechanical model of the induction motor. Fig.1. Extended Jeffcott rotor model 1.1. Differential equations of motion The equations of motion for the system described above can be derived using the Lagrange equations and have the following form: mpr + crxr + krxr = mrep2 cos^+mre(psiny>+krxr0 -FCx my + cy + k y = m^2 sin^-m^cos^k уй - F^ JM + crt (щ -Щ ) + krt Щг -Щ ) = Mfr1 - MC mXs + CsXs + K Xs = FCx msy + Csy + ks ys = FCy JÄ + CA + = MC Jm9m + Crt [Pm-ф )+ Krt (Щ-Щ ) = Mfr2 + Mm (1) The Equation (1) is taken from literature [12] and presents a further development of the models found in [2], [7] and [13], where mr, ms, c, cs, k, ks are mass parameters, lateral damping coefficients and stiffnesses of the rotor and stator respectively, crt, cst, krt and kst are torsional damping coefficients and stifnesses of the rotor shaft and stator, e is the rotor mass eccentricity, Jr, Js and Jm represent mass moments of inertia of the rotor, stator and induction motor respectively. Mfr1 and Mfr2 represent the torques of total losses (friction in bearings and losses in motor fan). These torques were identified from the measurements conducted on the test rig. To simplify the already complicated model, rotor and stator are assumed as isotropic regarding lateral stiffness and damping coefficient, what is still common practice in literature [14]. Vector rr0 from the Fig. 2 represents the static equilibrium position of the bowed shaft [13], where xr0 and yr0 are its X,Y components and are given by xr0 = r0x cos Щ - r0y sin Щ Уг0 = r0x sin Щ + r0x cos Щ (2) Components r0x and r0y represent the static equilibrium position of the bowed shaft in the rotor fixed x,y,z coordinate system. Fig.2. Jeffcott rotor with a bowed shaft Contact forces FCx and FCy as well as the contact moment MC are defined, according to Fig. 3, by the following expression: FCx = fn cos 7 — ft sin 7 FCy = FN sin 7 + FT cos 7 Mc = ft R (3) where the normal FN and the tangential FT contact forces have positive value when they are acting on the stator. Fig. 3 presents the geometry of the rotor - stator contact model with force definition. 1.2. Contact Models Traditionally, for the rotor-stator impact modeling, two different methods have been applied, namely the Newton's restitution coefficient model and the Contact forceindentation model. The first model is based on the simplifying assumption of perfectly rigid bodies. Actual physical objects are compliant and hence the impact duration is strictly greater than zero. This more realistic view of impact phenomena led many researchers to consider the continuous-dynamics models of collision where bodies deform during impact and the collision dynamics is treated as a continuous-time dynamic phenomenon. In its general form, the forceindentation relationship looks like [15]: Fn = Fk (0 + Fc (S, j) + Fp (S, Š) (4) where, Fk is the elastic (conservative) part of the normal contact force FN, Fc the viscous damping part and Fp the dissipative part due to plastic deformation. In the paper, plastic dissipations have been neglected. Fig. 4 shows the normal contact force model presented as a spring-damper system. In the development of the elastic relation Fk(ö), the Hertz theory has to be mentioned. A very commonly used expression is the forceindentation relation for sphere to sphere contact, according to Hertz: Fk = S (5) where, Fk is the normal force pressing the solids together, ö represents the approach of two spheres, i.e. total deformation of both surfaces while ks is a constant depending on the sphere radii and elastic properties of the sphere materials. Rotor Fig.4. Normal contact force model - a spring-damper system For the case of a cylindrical joint, Rivin [16] based his analysis on a half sine distribution of contact pressure between shaft and bearing. To calculate the pressure distribution, he considered two possible models for the force - deformation relation, linear and quadratic: Fk = kxS Fk = k2S2 (6) (7) Fig.3. Geometry of the rotor - stator contact model Results of his experiments carried out with cast iron and hardened steel showed that the quadratic model offered a better correlation of the force deflection in a journal bearing when a high load per unit length was applied. For low loads, the linear model represented a better fit. The simplest model of viscous dissipation Fc (s, j) which is known from the literature [1] is a linear damper, F = cS (8) where c is the viscous damping coefficient. Non-linear expression for visco-elastic force, originally proposed by Hunt and Crossley [17], is given by F = ÄSnS (9) where Ä and n are model parameters. Fig. 5 shows responses of two analyzed models of the contact force regarding different combination of the elastic force models Fk and damping force models Fc. The area bounded by hysteresis loop represents a loss of work during an impact of rotor and stator and can be expressed by, W = JFcd: (10) The linear contact force model is described by the well known linear equation: mS + cS + kS = 0 (11) where c is the viscous damping coefficient, k is the stiffness, S = Г - Cr is the local deformation of the rotor and stator in the contact point in normal direction, Si is the relative velocity of the rotor and stator along normal direction during the contact, |r| = -J(xr - xs )2 +(yr - ys )2 is the magnitude of vector difference between rotor and stator displacement according to Fig. 3 and Cr is the radial clearance between rotor and stator. Fig.5. Models of normal contact force In order to eliminate the negative forces that appear in the linear model and with the improvement in the sense of avoiding the discontinuity, Hunt and Crossley formulated the equation: miS + (äS11)S + kSn = 0 (12) Their theory gives a relationship between restitution coefficient and dissipation in contacts via assumption of linear dependency of the restitution coefficient s and initial impact velocity v0 i.e. s = 1-av0 (13) Thus, according to [17] Equation (12) can be written in the following form mS + kSn -aS +1 0 (14) where a is the constant dependent mainly on material and geometry of bodies in contact. Based on Equation (14) two separate parts of normal contact force can be identified, Fn = Fc + Fk =- akSn S + kSn (15) According to [16], index n can vary between 1 and 2 for cylindrical joints, so the choice of n = 3/2, originally proposed by Hertz (see Equation (5)), seems fairly reasonable. In this paper both normal force models (linear-Equations (6), (8) and nonlinear-Equation (15)) have been used and compared. Tangential contact force has been represented by the Coulomb dry friction law. The definition of the contact model described above can be expressed as follows FM 0 for r < Cr Fc for r > Cr FT = McFN sgn(vsl) (16) where m is the coefficient of friction and vsl is the sliding velocity between rotor and stator surfaces in contact. The angle у from equation (3) represents the normal force angle or direction of the smallest gap when rotor and stator are not in contact i.e. tan у ; Уг - ys (17) To determine S and vsi further vector calculation should be taken. Fig. 6 shows the velocity definition for the rotor and stator in contact. If the normal direction is known, x n0 =M , r = [rx, ry] = [xr - xs, Уг - ys] (18) as well as the velocities of rotor and stator centers [vxs, vys], one can easily Vr = K^ vyr] and Vs obtain the relative normal rotor-to-stator velocity, Š = vrn - vsn (19) where vrn = vr • n0 and vsn = vs • n0 represent the rotor and stator velocities in normal direction. Fig.6. Velocity definition for rotor and stator in contact The sliding velocity vsl can be further defined as, "sl — "rt - Kst (20) where vrt = vrt0 + pr R and vst = vst0 + ps (R + Cr ) represent the rotor and stator velocities on its periphery in tangential direction. R is the radius of the rotor while Šr, Šs are angular velocities of the rotor and stator. 1.3. Mechanical Model of the Induction Motor The driving moment of the induction electric motor has been modeled according to the Kloss expression [18] and [19] Mm = M„ 2 + 2sn sp ■ + 2sn sp (21) where Mmax and smn are the maximum motor torque or breakdown torque and slip for maximum motor torque, s is the actual slip of motors rotor with respect to the synchronous electromagnetic field and p = ff / f is the ratio of actual frequency given from frequency inverter and nominal frequency given by electrical network (f =50 Hz). If we replace Mmax in Equation (21) with Mmaxb =-Mmax I1 + smn )/(1-smn ), the Kloss expression for regenerative braking is obtained. The regenerative braking starts when actual speed of the motor's rotor is greater than the synchronous speed. The motor is then behaving as an electromagnetic brake. Fig. 7 shows the speed torque characteristics of an induction motor. Fig.7. Speed torque characteristics of an induction motor 2 NUMERICAL SIMULATIONS The numerical integration has been carried out by the fourth-order Runge-Kutta method. Time step of At = 2 • 10 5 s , has been applied in all simulations where the contact between rotor and stator was reasonably expected. Otherwise At = 1 • 10-4 s has been used. In some earlier simulations a greater time step was used but it turned out that it couldn't describe some specific rotor behaviors. Time step during one simulation has been fixed, so the subsequent frequency analysis could be done. The parameters used in all calculations are as follows: rotor mass mr = 4.25 kg, stator mass ms = 3.838 kg, mass moments of inertia Jr = 4.532-10-3 kgm2, Jm = 2.115-10-3 kgm2, Js = 2.11 • 10-2 kgm2, rotor disc radius R = 0.06 m, stiffnesses kr = 131540 N/m, krt = 1080 Nm/rad, ks = 1.237-106 N/m, kst = 7917 Nm/rad, damping coefficients cr = 22.43 Ns/m, crt = 0.0374 Nms/rad, cs = 11.96 Ns/m, cst = 0.0258 Nms/rad, radial clearance Cr = 0.4 mm, motor breakdown torque Mmax = 10.25 Nm and motor speed of max. breakdown torque (at power network frequency f = 50 Hz) nm = 2400 rpm. The identified dependence between angular deceleration and angular velocity of the rotor, obtained from free run down tests (induction motor switched off), had the following form = -4-10-6®3 + 0.0020 -0.53760-10.513 rad /s2 (22) 2.1. Numerical Simulation of Motor Controlled Run-Down (RD) Analysis The numerical RD analysis has been performed to identify the basic dynamic characteristics of the rotor-stator system and to correlate them with the experimental results. Knowing, from the modal testing (rotor response on the impulse force excitation at standstill) that the first rotor natural bending frequency is equal to f = 28.0 Hz and according to the Jeffcott (or Laval) theory it is the only natural frequency, RD analysis has been focused on speed range n = 40 to 0 Hz. The second rotor natural bending frequency obtained by the modal testing has been fr2 = 148.5 Hz hence it has been far away enough from the first natural frequency and its influence has been negligible. The basic goal of this part of numerical analysis was to discover the appropriate relation of the mass unbalance eccentricity vector e and the static equilibrium position of the bowed shaft rr0 according to the measured maximum lateral displacement amplitude of the rotor disc and the qualitative shape of the rotor lateral response in time domain, Fig. 8. Final parameters have been e = 0.046 mm, rr0 = 0.05 mm and phase lag of the e in relation to rr0 i.e. a0 = n rad (see Fig. 2). RD analysis has been controlled by a linear frequency ramp (speed law) of the synchronous electromagnetic field of induction motor, ff = 40 to 0 Hz in the time period of 34 s, while its rotor behaved according to the presented mechanical model. At the beginning of the numerical simulation, 2 s have been taken additionally for disappearance of any transients caused by initial conditions. Fig.8. Simulated lateral rotor response in horizont direction for RD analysis and speed law Fig. 9 shows the spectral map of the lateral rotor response in the horizontal direction where as expected only the first order can be seen. On the same figure at the beginning of the simulation, the rotor critical speed is excited but only because of the numerical transient due to initial conditions. 200 Fig.9. Spectral map of the lateral rotor response in horizontal direction xr 2.2. Numerical Simulations of Sudden Rotor Unbalance Increase (SRUI) Analysis The basic assumption has been that the rotor is running at a nominal speed above the first rotor critical speed and that the rotor is well balanced. To achieve the initial condition which is similar to the steadily running speed and to avoid any transients, the sudden unbalance was introduced 1 second after the start of simulation. After sudden appearance of the additional unbalance madd it has been assumed that the permanent monitoring system of the turbo machine has recorded unallowable vibrations and as a precaution measure it starts an immediate shut down procedure. This was accomplished by switching off the induction motor. Then the whole rotor starts to decelerate due to given torques of losses, expressed by the identified rotor angular deceleration mentioned above. The rotor is then approaching to its critical speed and because of additional unbalance muadd it contacts the stator. When the additional unbalance is big enough, the rotor contacts the stator immediately after appearance of the additional unbalance i.e. before switching off the motor. The time period for response of the shut down procedure after sudden unbalance appearance has been set to 1 s. Two different contact normal force models have been compared i.e. linear and nonlinear. Linear model has contact stiffness kC = 7.7-107 N/m and contact damping cC = 3.2-103 Ns/m and nonlinear model has contact toughness kC = 1.8-108 Pa/m12 and parameter a = 3 s/m. The output inverter frequency before unbalance appearance was set to nsf = 40 Hz. It turns out that the frictional losses in test rig are quite big so the rotor decelerates, with small unbalance and without contact, in some 5 to 6 seconds from 40 Hz to standstill. In Figs. 10 to 12 the rotor response without contact to stator is presented, without any additional unbalance, madd = 0 (or presented via rotor mass eccentricity eadd = 0). If someone compares Fig. 9 and Fig. 12, except differences in the shape of first order (linear in Fig. 9 and nonlinear in Fig. 12), on Fig. 12 exists additional horizontal line at 28 Hz. The reason for this mainly lies in fast rotor speed change while passing though its critical speed. 0.5 У """ X \ / \ / \ t ) l \ \ 35*, \ о I 1 §mЯ шу~\ ) \ 1Я Ж1 if...................1...... ЩМ 1 \ \ \ i J / \ / S У 0.5 -0.5 0 0.5 a . mm Fig. 10. Rotor center orbit during free RD (red dashed line represents a clearance circle) 0.5 ■0.5-i-i-i-i-i-i-i- 012345678 50 j-1-,-,-,-,-1-,- 40-;-k. N 3020.........................................;.................................... 10............;...................................................... o I-i-i-i-;-i- --- 012345678 t, s Fig. 11. Simulated lateral rotor response in horizontal direction for free RD analysis with speed law eadd = 0 Fig. 12. Spectral map of rotor displacements in horizontal direction, free RD, eadd = 0 ■0.4 -0.2 0 0.2 0.4 0.6 :шг Fig.13. Rotor center orbit during RD, rotorstator contact appeared, eadd = 0.124 mm 0.2 -E "I E 0-0.1 --0.2 4 t, s Fig.14. Lateral rotor andstator responses in horizontal direction, contact appeared, eadd = 1.24-10' m, rr0 = 0.05 mm, a0 = n rad .......';. \\ ---Without contact \ \ -Contact \ \ \ \ \ ч \ N \ ч v \ 1: 0.5 0: 3 4 5 f, s Fig.15. Rotor speed laws for simulations without and with contact 200 150 I 100 Fig. 17. Spectral map of rotor relative torsional angles pm - фп contact appeared 200 150 I 100 Fig. 18. Spectral map of stator lateral velocities in horizontal direction, contact appeared 200 150 I 100 Fig. 16. Spectral map of rotor displacement in horizontal direction, contact appeared 200 150 I 100 Fig.19. Spectral map of stator torsional angles, contact appeared Figures 13 to 19 show rotor and stator response after the appearance of the additional mass unbalance (presented in rotor mass eccentricity) eadd = 0.124 mm, so the total rotor mass eccentricity is equal to e = ein + eadd =1.17 mm. The magnitude of the vector rr0 was rr0 = 0.05 mm and the phase lag was a0 = к rad. In this simulation the linear model of normal contact force has been applied and coefficient of friction in contact was /лс = 0.18, according to Bartha [2]. In Fig. 11 we had the supercritical self-balancing effect because the mass eccentricity e was smaller than the magnitude of the shaft bow vector rr0. On the contrary in Fig. 14 the self balancing effect is placed subcritically at the t1 = 3.8 s of simulation because after sudden unbalance increase, the total mass unbalance eccentricity e became greater than rr0. Fig. 15 shows the difference between speed laws for simulations with and without the rotor-stator contact appearance. As expected a rotor deceleration for simulation where the rotorstator contact appeared, was more intensive especially in a time period with established rotorstator contact (see contact indication in Fig. 15). After introduction of SRUI, the rotor made several impacts to the stator and then a permanent contact was established. In Figs. 16 to 19 there is a narrow time period around t1 = 1 s of simulation where multiple harmonics can be seen. This is due to intermittent rotor-stator contact. It is interesting to see after the rotor-stator separation (t1 = 2.75 s) that the rotor and stator continued to vibrate with their own flexural natural frequencies, fr = 28 Hz (Fig. 16), f = 90 Hz (Fig. 18) and torsional frequencies ft = 137.9 Hz (Fig. 17), ft = 97.7 Hz (Fig. 19). In Fig. 20, the influence of the phase lag of the mass unbalance eccentricity vector e in relation to the static equilibrium position of the bowed shaft rr0 i.e. a0 has been analyzed. Presented contact forces are related to a0 = 180°, 135° and 90°. In the same time normal contact forces of the linear and nonlinear models have been compared, so the diagrams in the left column correspond to the linear and in the right column to the nonlinear model of the normal contact force. The main difference between the linear and nonlinear model of the contact normal force, presented in Fig. 20, can be seen in the region immediately after appearance of SRUI. Nonlinear model generally had a greater force response for the first impact, while in permanent contact both linear and nonlinear model had almost the same values. Details about the differences between analyzed normal force models and influence of the phase lag a0 are presented in Table 1. The normal force for the permanent rotorstator contact (with sliding) FCperm has showed expected increasing trend with decreasing phase lag a0, because a0 smaller than 180° means greater Table 1. Results of contact normal force analysis regarding phase lag of the mass eccentricity vector e in relation to shaft bow radius vector rr0 а„, FCmax / FCperm , N Linear Nonlinear 180° 327.7 / 96.8 611.9 / 97.1 135° 279.8 / 99.3 462.2 / 99.5 90° 203.7 / 104.8 277.1 / 105.0 unbalance excitation. The reason for contradictory trend of decreasing of the FCmax (max. value of the normal force for the first impact) with decreasing a0 (increasing unbalance excitation) lies in the fact that a0 of the ein and eadd is the same. This caused an increasing tendency of the radius of initial rotor orbit and therefore the decrease of the maximum normal velocity at the first rotor-stator impact. 3 EXPERIMENTAL ANALYSIS The analysis has been performed on the test rig which can be seen on Fig. 21. The test rig has been specifically designed for the rotor-stator contact measurements. It consists of the following main components: robust foundation, mounting plate, roller bearings with their supports, rotor, stator, flexible coupling, induction motor with speed controller, measuring system based on Bruel and Kjaer non-contacting displacement sensors and National Instruments data acquisition card PCI NI 4472 with adequate software (LabVIEW, Matlab). Measurements showed that the first critical speed of the rotor was 28 Hz, while the first natural frequency of the stator was 90 Hz. Radial clearance between the rotor and stator was 0.4 mm. Before the experiments the rotor had been balanced. a) b) с) Linear model Non-linear model Fig. 20. Normal forces in contact between rotor and stator; left column - linear model, right column nonlinear model; a) a0 = 180°, b) a0 = 135°, с) a0 = 90° Also before every single experiment the contact surfaces of the rotor and stator were lubricated with WD40 spray [2] to decrease the coefficient of friction to the value of approximately ß = 0.18 and to minimize the possibility of destruction of sensors and other structural parts of the test rig. 3.1 Experimental, Motor Controlled, RunDown (RD) Analysis To verify the numerical results presented in chapter 2.1, an experimental RD analysis has been performed. Fig.21. Test rig for rotor-stator contact investigation The speed change has been controlled via frequency inverter Micromaster 440 from Siemens. The speed rate has been the same as in the numerical simulation only the starting speed has been 70 Hz. Figs. 22 and 23 show measured lateral responses of the rotor in horizontal x direction with speed law and spectral map of the rotor lateral response in the horizontal direction. This figures can be compared to Figs. 8 and 9 respectively, for a speed range n = 40 to 0 Hz. In Fig. 23 except for the first harmonic, other higher harmonics can be also seen due to various imperfections of the experimental model like radial and angular misalignment etc. Fig. 24 shows a spectral map of stator displacements in horizontal direction. Although there wasn't direct contact between the rotor and stator, the stator was excited as it can be seen on Fig. 24. This is due to the fact that rotor bearing supports and stator supports are rigidly connected on the mounting plate, so this situation allows an indirect stator excitation through the mounting plate. Measured flexural natural frequencies were, fr = 28 Hz (Fig. 23), f = 90 Hz (Fig. 24) and stator torsional frequency fst = 102.5 Hz (Fig. 24). To excite effectively both stator natural frequencies (with second harmonic), the 70 Hz as a starting speed of this RD analysis was chosen. Fig.22. Measured lateral rotor xr response in horizontal direction and speed law Fig.23. Measured spectral map of rotor displacement in horizontal direction ж ' « ж mm в ■ • Л. - 7 т с,, я- ■„-■ ,, . -....... . ■■ ■■ i'. .-.■ - - ■ .. 'ì,1 j .....i^^vi Fig.24. Measured spectral map of stator displacement in horizontal direction 3.2 Experimental Sudden Rotor Unbalance Increase (SRUI) Analysis Because this test rig couldn't simulate the SRUI scenario during the rotor operation, additional mass has been added instead of being removed but with rotor at standstill. The rotor has been quickly accelerated to rotational speed of 40 Hz to reduce the rubbing while passing through rotor critical speed. The power supply to the motor was then switched off and a whole rotor started a free deceleration. Two measurements have been presented: free rotor RD with good rotor balance condition as in previous analysis (without appearance of the rotor-stator contact), Figs. 25 and 26 and free RD of the rotor with additional mass unbalance (presented in the rotor mass eccentricity) eadd = 0.124 mm, Figs. 27 to 30. The first measurement (Figs. 25 and 26) can be compared directly with the numerical simulation presented in the Figs. 11 and 12. Looking just on Figs. 11 and 25 a great correspondence can be noticed. Comparing spectral maps of simulated and measured lateral rotor displacements, i.e. Figs. 12 and 26, some differences regarding presence of higher harmonics in measured response can be observed. The reason for their presence is explained in the previous chapter and their influence is reflected on extended life of rotor bending natural frequency present in the response till the rotor standstill. The second measurement (Figs. 27 to 30) can be compared with the numerical simulation presented in Figs. 13 to 19 with remark that in simulation the motor is switched off in t, = 2 s (Fig. 14) while in measurement in t, = 1.65 s. Furthermore, rotor and stator were, in measurement, in permanent contact from the very beginning, while in numerical SRUI simulation they got into contact after additional mass unbalance activation i.e. in t, = 1.0 s of simulation. Taking into account the aforementioned reasons for differences between simulation and measurement, a great correspondence between Figs. 14 and 27 can be noticed again. In measured spectral map of the rotor lateral response, Fig. 28, vibrations have strong first harmonic component i.e. 1x rotating frequency, while there are higher harmonics too. They were present also in the rotor lateral response without appearance of the rotor-stator contact (Fig. 26), because of the above mentioned imperfections of the rotor experimental model. Both lateral and torsional natural frequencies of the stator are present in each of its depicted responses, i.e. in the lateral response (Fig. 29) and in the torsional response (Fig. 30), after the final rotor-stator separation. This happened mainly because of the vicinity of stator natural frequencies and due to indirect measurement of torsional stator rotations via 2 parallel non-contacting displacement sensors in horizontal (lateral) direction. Although both stator frequencies can be seen in two separated responses, lateral frequency is more emphasized in lateral response while torsional frequency is more emphasized in its torsional response. Thus in the case of torsionally elastic systems susceptible to the rotor-stator contact it is important to combine the measurements with results of the numerical simulations to make a decision of the character of a particular natural frequency. 4 CONCLUSIONS The presented numerical model of the rotor with a bowed shaft - stator system has proved to be capable for analyzing different contact situations. Linear and nonlinear normal contact force models were compared in the same simulated situations and they showed similar responses for permanent rotor-stator contact, while nonlinear force model generally had a greater response for the first rotor-stator impact. This didn't have a big influence on a subsequent rotor and stator response and their stability, for the analyzed parameters. Both normal force models, for permanent rotor-stator contact (with sliding) FCperm, showed increasing tendency with decreasing phase lag a0, because if a0 is smaller than 180° the unbalance excitation will be greater. Experimental analysis showed more realistic rotor and stator responses with more different harmonics present than in the results of numerical simulations. This was mainly due to tested model imperfections and measuring principle. In the real practice measured signals are even more unclear and mixed with different noise, so it is important to combine the measurements with adequate numerical simulations to recognize the needed information. Fig.25. Measured lateral rotor xr response in horizontal direction with speed law, contact-free situation, free RD, eadd = 0 mm DS 1 IS Fig.28. Spectral map of rotor lateral displacements in horizontal direction, free RD, contact appeared, eadd = 0.124 mm 200180- - . " 160 —Г 1.40 -120 ■ m X 100- Fig.26. Spectral map of rotor displacements in horizontal direction, free RD, eadd = 0 mm Fig.27. Measured lateral rotor xr and stator xs responses in horizontal direction, free RD, contact appeared, eadd = 0.124 mm Fig.29. Measured spectral map of stator lateral displacement in horizontal direction, free RD, eadd = 0.124 mm Fig.30. Measured spectral map of stator torsional rotations, free RD, eadd = 0.124 mm 5 REFERENCES [1] Choy, F. K., Padovan, J. (1987) Non-linear [10] transient analysis of rotor-casing rub events, Journal of Sound and Vibration, 113(3), pp. 529-545. [2] Bartha, A. R. (2000) Dry Friction Backward Whirl of Rotors, PhD. Thesis, Swiss Federal [11] Institute of Technology Zurich. [3] Von Groll, G., Ewins, D. J. (2002) A Mechanism of Low Subharmonic Responce in Rotor/Stator Contact - Measurement and Simulation, Journal of Vibration and [12] Acoustics, vol. 124 (2002), pp.350.-358. [4] Ahrens, J., Jiang, J., Ulbrich, H., Ahaus, G. (2001) Experimentelle Untersuchungen zum Schaufelanstreifen, Schwingungen in [13] rotierenden Maschinen - SIRM V Tagung, Wien, Austria, pp. 97-108, (in German). [5] Fumagalli, M., Schweitzer, G. (1994) Impact [14] dynamics of high speed rotors in retainer bearings and measurement concepts, 4th International Symposium on Magnetic Bearings, ETH Zurich Switzerland. [15] [6] Orth, M., Nordmann, R. ANEAS (2002) A modeling tool for nonlinear analysis of active magnetic bearing systems, 2nd IFAC Conference on Mechatronic Systems, [16] Berkley, USA, pp. 357-362. [7] Edwards, S., Lees, A. W. and Friswell, M. I. (1999) The Influence of Torsion on [17] Rotor/Stator Contact in Rotating Machinery, Journal of Sound and Vibration, 225(4), pp. 767-778. [8] Žigulić, R., Butković, M., Braut S. (2002) Nonlinear dynamics of multi-disc rotor in [18] dry friction bearings, IFToMM Sixth International Conference on Rotor [19] Dynamics, vol. 2, Sydney, Australia, pp. 960-967. [9] Karpenko, E.V., Pavlovskaia, E.E., Wiercigroch, M. (2003) Bifurcation analysis of a preloaded Jeffcott rotor, Chaos, Solitons & Fractals, 15, pp. 407-416. Pavlovskaia, E.E., Karpenko, E.V., Wiercigroch, M. (2004) Nonlinear dynamics of a Jeffcott rotor with a preloaded snubber ring, Journal of Sound and Vibration, 276 (1-2), pp. 361-379. Karpenko, E.V., Wiercigroch, M., Pavlovskaia, E.E., Neilson, R.D. (2006) Experimental verification of Jeffcott rotor model with preloaded snubber ring, Journal of Sound and Vibration, 298, pp. 907-917 Braut, S (2006) Analysis of the rotor-stator contact dynamics, PhD thesis (in Croatian), University of Rijeka, Faculty of Engineering, Rijeka. Childs, D. W. (1993) Turbomachinery Rotordynamics, Phenomena, Modeling, and Analysis, John Wiley & Sons, New York. Chu, F., Lu, W. (2007) Stiffening effect of the rotor during the rotor-to-stator rub in a rotating machine, Journal of Sound and Vibration, 308(2007), 758-766. Faik, S., Witteman, H. (2000) Modeling of Impact Dynamics: A Literature Survey, 2000 International ADAMS User Conference, Orlando, Florida, USA, pp. 1-11. Rivin, E. I. (1999) Stiffness and Damping in Mechanical Design, Marcel Dekker, Inc. New York. Hunt, K. H. and Crossley, F. R. E. (1975) Coefcient of restitution interpreted as damping in vibroimpact, Transactions of the ASME, Journal of Applied Mechanics, pp. 440-445. Jurković, B. (1990) Electric power drivers, Školska knjiga, Zagreb, (in Croatian). Boldea, I., Nasar, S. A. (2002) The Induction Machine Handbook, CRC Press, Boca Raton. Paper received: 07.02.2007 Paper accepted: 07.07.2008 MINLP Optimization of a Single-Storey Industrial Steel Building Tomaž Žula1, - Zdravko Kravanja2 - Stojan Kravanja1 1 University of Maribor, Faculty of civil engineering, Slovenia 2 University of Maribor, Faculty of chemistry and chemical engineering, Slovenia The paper presents the topology and standard sizes optimization of a single-storey industrial steel building, made from standard hot rolled I sections. The structure consists of main portal frames, connected with purlins. The structural optimization is performed by the Mixed-Integer Non-linear programming approach (MINLP). The MINLP performs a discrete topology and standard dimension optimization simultaneously with continuous parameters. Since the discrete/continuous optimization problem of the industrial building is non-convex and highly non-linear, the Modified Üuter-Approximation/Equality-Relaxation (ÜA/ER) algorithm has been used for the optimization. Alongside the optimum structure mass, the optimum topology with the optimum number of portal frames and purlins as well as all standard cross-section sizes have been obtained. The paper includes the theoretical basis and a practical example with the results of the optimization. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: industrial buildings, topology optimization, sizing optimization, non-linear programming, MINLP 0 INTRODUCTION Single-storey frame structures are extensively used for industrial, leisure and commercial buildings. In order to obtain efficient frame designs, researchers have introduced various optimization techniques, suitable either for continuous or discrete optimization. O'Brien and Dixon [1] have proposed a linear programming approach for the optimum design of pitched roof frames. Guerlement et al. [2] have introduced a practical method for single-storey steel structures, based on a discrete minimum weight design and Eurocode 3 [3] design constraints. Recently, Saka [4] has considered an optimum design of pitched roof steel frames with haunched rafters by using a genetic algorithm. One of the latest researches reported in this field is the work of Hernandez et al. [5], where the authors have considered a minimum weight design of the steel portal frames with software developed for the structural optimization. It should be noted that all the mentioned authors deal with the discrete sizes optimization only at fixed structural topologies. This paper discusses the simultaneous topology, standard sizes and continuous parameter optimization of an unbraced single- storey industrial steel building. The optimization of the portal frames and purlins was performed by the Mixed-Integer Non-linear Programming approach (MINLP). The MINLP is a combined discrete and continuous optimization technique. In this way, the MINLP performs the discrete topology (i.e. the number of frames and purlins) and the standard dimension (i.e. the standard cross-section sizes of the columns, beams and purlins) optimization simultaneously with the continuous optimization of the parameters (e.g. the structure mass, internal forces, deflections, etc.). The MINLP discrete/continuous optimization problems of frame structures are in most cases comprehensive, non-convex and highly non-linear. The optimization is proposed to be performed through three steps. The first one includes the generation of a mechanical superstructure of different topology and standard dimension alternatives, the second one involves the development of an MINLP model formulation and the last one consists of a solution for the defined MINLP optimization problem. The objective of the optimization is to minimize the mass of the single-storey industrial building. The mass objective function is subjected to the set of equality and inequality constraints *Corr. Author's Address: University of Maribor, Faculty of Civil Engineering, Smetanova 17, 2000 Maribor, Slovenia, tomaz.zula@uni-mb.si known from the structural analysis and dimensioning. The dimensioning of steel members is performed in accordance with Eurocode 3. The Modified Outer-Approximation /Equality-Relaxation algorithm is used to perform the optimization, see Kravanja and Grossmann [6], Kravanja et al. [7] and [8]. The two-phase MINLP optimization is proposed. It starts with the topology optimization, while the standard dimensions are temporarily relaxed into continuous parameters. When the optimum topology is found, the standard dimensions of the cross-sections are reestablished and the simultaneous discrete topology and standard dimension optimization of the beams, columns and purlins is then continued until the optimum solution is found. 1 SINGLE-STOREY INDUSTRIAL BUILDING The paper presents the topology and standard sizes optimization of unbraced rigid single-storey industrial building steel structures, Fig. 1. The columns, beams and purlins are proposed to be built up of standard hot rolled steel I sections. The considered portal frame structures are optimized under the combined effects of the self- weight of the frame members, a uniformly distributed surface variable load (snow and wind), a concentrated horizontal variable load (wind) and an initial frame imperfection. The purlins are designed to transfer the permanent load (the self-weight of the purlins and the weight of the roof) and the variable load (snow and wind). The internal forces are calculated by the elastic firstorder method. The dimensioning of the steel members is performed in accordance with Eurocode 3 for the conditions of both the ultimate limit state (ULS) and the serviceability limit state (SLS). When the ULS is considered, the elements are checked for the axial, shear and bending moment resistance, for the interaction between the bending moment and the axial force, the interaction between the axial compression/buckling and the buckling resistance moment. The total deflection Smax subject to the overall load and the deflections S2 subjected to the variable imposed load are calculated to be smaller than the limited maximum values: span/200 and span/250, respectively. The frame horizontal deflections are also checked for the recommended limits: the relative horizontal deflection of the portal frame should be smaller then the frame height/150. Fig. 1. Single-storey industrial steel building 2 MINLP MODEL FORMULATION FOR MECHANICAL SUPERSTRUCTURES It is assumed that a non-convex and nonlinear discrete/continuous optimization problem can be formulated as a general MINLP problem (MINLP-G) in the form: min z = cT y + f (x) s.t. h(x) = 0 g(x )< 0 (MINLP-G) By + Cx < b x e X" = {x e R": xLO < x < x^} m У e Y={0,1} where x is a vector of continuous variables specified in the compact set X and y is a vector of discrete, mostly binary 0-1 variables. Functions f(x), h(x) and g(x) are non-linear functions involved in the objective function z, the equality and inequality constraints, respectively. All functions fx), h(x) and g(x) must be continuous and differentiable. All functions fx), h(x) and g(x) must be continuous and differentiable. Finally, By+Cx< b represents a subset of mixed linear equality/inequality constraints. The above general MINLP-G model formulation has been adapted for the optimization of mechanical superstructures. The resulting MINLP formulation for mechanical superstructures (MINLP-MS) that is more specific, particularly in variables and constraints, can be used also for the modelling the steel industrial buildings. It is given in the following form: min z = cT y + f ( x) s.t. h(x ) = 0 g(x )< 0 A(x )< a Ey < e (MNLP-MS) Dye + R ( x )< r Kye + L(dcn)< k Py + S (dst )< s x e X = {x e R": xLO < x < x"?} m y e Y={0,1} The MINLP model formulation for mechanical superstructures is proposed to be described as follows: - Included are continuous variables x={d, p} and discrete binary variables y={ye, yst}. Continuous variables are partitioned into design variables d={d°n, d5*} and into performance (non-design) variables p, where subvectors dcn and dst stand for continuous and standard dimensions, respectively. Subvectors of the binary variables ye and yst denote the potential existence of structural elements inside the superstructure (the topology determination) and the potential selection of standard dimension alternatives, respectively. - The mass or economical objective function z involves fixed mass or cost charges in the linear term cTy and dimension dependant mass or costs in the term fx). - Parameter non-linear and linear constraints h(x)=0, g(x) < 0 and A(x) < a represent a rigorous system of the design, loading, resistance, stress, deflection, etc. constraints known from the structural analysis. - Integer linear constraints Ey < e are proposed to describe the relations between binary variables. - Mixed linear constraints Dye+R(x) < r restore interconnection relations between currently selected or existing structural elements (corresponding ye=1) and cancel relations for currently disappearing or nonexisting elements (corresponding ye=0). - Mixed linear constraints Kye+L(dcn) < k are proposed to define the continuous design variables for each existing structural element. The space is defined only when the corresponding structure element exists (ye=1), otherwise it is empty. - Mixed linear constraints Py+S(dst) < s define standard design variables dst. Each standard dimension dst is determined as a scalar product between its vector of i, iel, discrete standard dimension constants ^={q1, q2, q3,..., qi} and its vector of subjected binary variables yst={ystb у\ у\.., ysti}, see Eq. (1). Only one discrete value can be selected for each standard dimension since the sum of the binary variables must be equal to 1 Eq. (2): ,st V^ st d =E ql Ji i e I E yt =1 (1) (2) 3 OPTIMIZATION MODEL FRAMEOPT The MINLP optimization model FRAMEOPT (FRAME OPTimization) for the optimization of the single storey industrial steel buildings has been developed with relating to the above MINLP model formulation for mechanical structures. The following assumptions and simplifications have been defined in the model FRAMEOPT and considered in the optimization: - Considered was a single load case only, where the partial safety factors and combination of actions were defined according to Eurocodes. The optimization of the structure was performed under the combined effects of: - the self-weight of the structure (the line uniform load of columns, beams and purlins) and the weight of the roof (the vertical surface load) plus - snow and vertical wind (the uniformly distributed vertical surface variable load) plus - horizontal wind (the horizontal force at the top of the columns). - Equal steel portal frames and equal purlins were proposed to compose the structure. - Steel members were proposed to be made from standard hot rolled European wide flange I sections (HEA sections). - The global portal frame geometry including the span, height and the beam inclination was proposed to be fix through the optimization. - Vertical and horizontal bracing systems as well as wall sheeting rails were not included in this calculation/optimization. - The internal forces and deflections were calculated by the elastic first-order method. - The portal frames were classified as non-sway steel portal frames. The ratio between the design value of the total vertical load NSd and the elastic critical value for failure in a sway mode Ncr was constrained: NSd/Ncr< 0.1. - The portal frame was calculated as a laterally supported frame. Hereby, the steel members were checked only for the in-plane instability. Columns were designed for the compression/buckling resistance plus the lateral torsional buckling. Beams were checked for the in-plane bending moment resistance. - Buckling lengths of columns were calculated as the in-plane buckling lengths for the non-sway mode. As an interface for mathematical modelling and data inputs/outputs GAMS (General Algebraic Modeling System), i.e. a high level language, was used [9]. The proposed optimization model includes the structure's mass objective function, parameter structural nonlinear and linear constraints, integer and mixed integer logical constraints, sets, input data (constants) and variables. 3.1. Mass objective function The mass objective function of the industrial building structure is defined by Eq. (3). The mass of the structure MASS comprises the masses of columns, beams and purlins. AC, AB and AP represent the cross-section areas of the column, beam and purlins, respectively. h denotes the height of the column, LB is the length of the frame beam and LL is the length of the industrial building (and purlins). NOFRAME represents the number of portal frames and NOPURLIN denotes the number of purlins. Each portal frame is constructed from two columns and two beams, see Fig. 2. MASS = 2• (AC • h • p) NOFRAME + 2 • (AB • LB • p) • NOFRAME + (Ap • Ll • p) NOPURLIN (3) 3.2. Parameter structural non-linear and linear constraints The first constraints of the model represent the constraints (4) to (7) which determine the relations between the continuous cross-sectional dimensions and the cross-sectional height of the column hC. These equations accelerate the convergence of the optimization when standard dimensions are re-established. They define the section breadth bC, the flange thickness tf,C, the i e I Ac Ab Ib Ic Ab Ib Ac Ic л A C : be A в : Fig. 2. Portal frame and cross-sections of elements web thickness twC and the cross-section area AC constant Iro, C for the frame column are given by (see Fig. 2) for the column. The second moments Eqs. (8) to (11). Similar cross-sectional constraints of the area about the y-y and z-z axis, Iy, C and are defined for the frame beam, Eqs. (12) to (16), Iz,C, the torsional constant It, C and the warping and for the purlins, see Eqs. (17) to (24). bC = -8.7681-10 ■ hC + 3.5913• 10 ■ hC - 5.9883-10-7 ■ hC5 + 5.1897•10-6 ■ hC - 2.457840-3 ■ hC + tf,C = 1.58C1•1C-8■ hC; + 3.4958•1C-6 ■ hC + 2.3488■Ю-4 ■ hC-1.9322■Ю-3 ■ hC + 0.76681 tw, C =-1.0598 ■ 10-5 ■ h C+2.4652 ■Ю-3 ■ hC + 0.23804 AC = 2 ■ be ■ tf e + (he - 2 ■ tf (4) (5) (6) (7) 1 y, C = 2■ be ■ fe , tw, e■(he -2■ tf, e) 12 12 he tf, (8) 2 ■ tf, e ■ bC " 12 (he 2 tf, e )■ tw, e 12 11, e = 3 ■(2 ■ bc ■ tf3 e j -■ (( - 2■ tr L ,= с л \hC 2 ■ tf, с) (9) (10) (11) bB = -8.7681 •10-12 ■ hB + 3.5913■Ю-9 ■ hB -5.9883•1C-7 ■ h£ + 5.1897■Ю-6 ■ hB --2.4578■Ю-3 ■ hB + 6.007■Ю-2 ■ hB -5.8757■ hB + 29.294 tf, B = 1.58C1•1C-8■ hj4 + 3.4958■Ю-6 ■ hB + 2.3488■Ю-4 ■ hB-1.9322■Ю-3 ■ hB + 0.76681 tw, B =-1.0598 •Ю-5 • h B +2.4652 •Ю-3 • hB + 0.23804 AB = 2 • bB • tf, B + (hB - 2 • tf, B )• tw, в (12) (13) (14) (15) z y t L I z, C Iy, в - 2 ■ Ьв ■ f в , tw, в■(( - 2 ■ f 12 12 -2■ bB ■tf>в ■ К bp = 8.7681 ■ 10-12 ■ hp7 + 3.5913■Ю-9 ■ hp6 -5.9883 10-7 ■ hp5 + 5.1897■Ю-6 ■ hp4 --2.4578■Ю-3 ■ hp + 6.007■Ю-2 ■ hpr -5.8757■ h7 + 29.294 tf, 7 = 1.5801 ■Ю-8■ hp + 3.4958■Ю-6 ■ hp + 2.3488■10-4 ■ hp2 -1.9322■Ю-3 ■ h7 + 0.76681 tw„ Ap = 2 ■ bp ■ t, p + (bp - 2 ■ t, Iy, p =- 12 12 -2 ■ bp ■ tf p ■ 2■ tf, p ■ bp3 , (bp - 2t 12 12 1p =—■ (2■ bp ■ fp ) + -■( - 2■ t,p )■ tw I,.,» =■ 'z, p 4 ' , - 2 ■ tr (16) (17) (18) (19) (20) (21) (22) (23) (24) z, p The length of the frame beam Lв is calculated according to Eq. (25) and the angle of the inclination of the beam a is defined by Eq. (26). L represents the frame span and f denotes the overheight of the frame beam: L =V(L/ 2)2 + f1 (25) a = arctan (//(L/ 2)) (26) The uniformly distributed vertical surface variable load qz, the uniformly distributed horizontal surface variable load qy, the self-weight per unit length of the portal frame g, the concentrated design horizontal variable wind load P (for the ULS) and wind load Pw (for the SLS) are defined by Eqs. (27) to (31): qz =(s ■ cos2 (a) + wv )ef (27) qy = s ■ cos(a) ■ sin(a) ■ ef (28) g = Ав ■ p + (Ap ■ ^ ef )/ep (29) P=7ll ■ Wh ■ ef ■ h/2 (30) Pw = Wh ■ ef ■ h/2 (31) Where s, wv and wh represent snow, the vertical and horizontal wind per m2 (the variable imposed load); ef stands for the intermediate distance between the portal frames, p is the density of steel, yq is the partial safety factor for the variable load and h represents the height of the columns. The number of the portal frames NOFRAME, the number of purlins NOPURL1N and the maximal intermediate distance between the purlins ep are determined by Eqs. (32) to (38), where LL represents the length of the industrial building, MINNOtame and MAXNOtrame denote the minimal and maximal number of defined portal frames, and MINNOpurlin and MAXNOpurlin stand for the minimal and maximal number of purlins. NOFRAME=L^l ef +1 (32) NOFRAME > MINNOframe (33) NOFRAME < MAXNOfam° (34) NOPURLIN = 2 ■ ((/ e +1) (35) NOPURLIN > M7NNOpurlin (36) NOPURLIN < MAXNCT^ (37) ep < 250 [cm] (38) Eqs. (39) to (45) represent the constraints which determine the portal frames to be non-sway frames. The column stiffness coefficient KC, the effective beam stiffness coefficient KB, the distribution factor of the column for the sway frame rfì, and the plane buckling length of the column for a sway frame mode ßsway are calculated by Eqs. (39) to (42). The value of the distribution factor r2 is taken to be 1 because of the pinned connection of the columns. K = ^ Kb = ~ r K Kc +1.5-Kb ßswüs 1-0.2 - rr +r) - 0.12-rr-r 1-0.8- rr + r) + 0.6-rS-r (39) (40) (41) (42) Eq. (43) represents the elastic critical load ratio (Nsd/Ncr) which defines the steel portal frame to be a non-sway frame. The distribution factor of the column for the non-sway frame r1NS and the plane buckling length of the column for the non-sway frame mode ß non-sway are given by Eqs. (44) to (45): P - h + + L 2 л2 - E -Iy < 0.1 (43) Г = K KC + 0.5 - KB (44) Äon-Sway =0.5 + 0.14- (rNS + r ) + 0.055- (rNS + r ) (45) The ULS constraints for the frame columns are defined by Eqs. (46)-(52). Eq. (46) represents the condition for the design bending moment resistance of the column (MsdL j3 + 5 A) + P•h• (B + C) -+- • A • fy/ Tm + 16 •(B + A • C) 2 •(B + A • C) *"lT • 2 • Iy, C • fy/ (hC •Tmö ) < 1.0 T • qz + Tg • g )•L •C3 +5-A) + P • h •(B + C ) < 2 • Iy, B • f 16 •(B + A • C ) P• h (Tq • qz +Tg • g)•L cos(a) — 2•(B + A• C) hB • Tmo Tq • qz +Tg • g )• L2 ^(3 + 5 • A) P • h •(B + C) 16 •(B + A^ C) 2^(B + A • C) /h • sin(a) (1.04• h, • t.,b )• fy S- TM0 ( Tq • qz +Tg • g )• L2 j3 + 5 • A) P • h •(B + C ) 16 •(B + A • C) 2 •(B + A • C) (Tq • qz + Tg • g )• L2 •( + 5 • A) P • h •(B + C ) 16 •(B + A • C) + 2 •(B + A^ c) h + Ph + (q1 + r£g))L_ sin(a) < B fy TM0 L p• h (Tq• qz + Tg• g)•L . , ч /h + +-2 —--sin(a) AB • fy/ TM! ( Tq • qz + Tg • g )• L2 •( + 5 • A) P • h •(B + C ) 16 (B + A• C) + + A• c) 2• Iy, в • fy/(hB •Tmo) <1.0 (52) (53) (54) (55) L 1 1 L 2 0.1057 • yg • (AP•p + mr ) • ef2 + 0.1057 • yq • (s• cos(«)• ep + wv • ep)• ef2 < 2 •Iy, p • fy hp • ум0 0.1057•у •(s + mr)• sin(a)• e • ef2 < 2 • Iz, p • fy hP • yM0 0.1057•yg •(APp + mr)• ef2 + 0.1057•у, •(s• cos(o)• ep + wv• ep)• ef 2 •iy, p • fy/hP Ум0 0.1057• yq •(s + mr)• sin(a)• ep • ef 2•iz, p • fy/hp •уме < 1.0 0.567•yg •(A^p)^ + 0.567•у,•(s• cos(a)• ep + wv• ep)• ef < 1.04 • hp • p • f S • yM0 TTZ (Ì-1 -U> hl E • I 1 • L •(V - U )• h -1 • L • 6 3 (q + g )• L2 • h U= (q + g)■ L2 .(3 + 5■ A) 16 ■( + A ■ C ) A< h/150 V = P ■ h ■(B + C ) 2 ■ (B + A ■ C ) 4 =- (q+g) 24 • E • I,,„ L 2V 2 (q + g)• L L 2 •V 2 (q + g)• L 12 •(-U) (q+g) 2 • L + ■ 8 •V q+g )• L L 2 •V L 2 •V 2 (q + g) L + L + 2 (q + g)• L 12 •U - 4 •V )• L (q+g) 4 =L/ 250 s • cos(«) • ep + ( wv • ep + AP • p • ep + mr • ep )• cos(«)) • (269 / 42000) • ef4 / (E • Iy, P ) I • cos(«) + +[(s • ep + AP • p • ep + mr • ep )• sin(«) • (269 / 42000) • ef4 / (E • I, P )) • sin(«) < ^^ 3.3. Integer and mixed integer logical constraints The logical constraint in Eq. (66) defines the number of portal frames, where yn denotes the binary variable which is subjected to each portal frame. Eq. (67) defines only one possible vector of binary variables for each frame topology. Eq. (68) calculates the even number of purlins, where the binary variables ym are subjected to the purlins. Eq. (69) defines only one possible vector of the binary variables for each purlin topology. NOFRAME = ^ yn (66) n Уп < Уп-1 (67) (58) (59) (60) (61) (61 a,b) (62) (63) (64) (65) NOPURLIN = 2 ym (68) m ym < ym-1 (69) Eqs. (70) to (79) calculate the standard cross-sections for the columns with their discrete dimensions and characteristics. The latter are determined as scalar products between their vectors of i, ieI, the discrete standard constants (qAC, qh...) and their vector of subjected binary variables_yi, see Eqs. (70) to (78). Only one discrete value is then selected for each standard section since the sum of the binary variables must be equal to 1, see Eq. (79). Ac=Eq^c *yi ie I (70) he=E qhC • yi bc=E qbC • yi tw, e=E q'w' C • y tf, e=E qf C • y i Iy, e=E q'" C • y i Iz, e=E qiz> C • y i It, e=E qiItC • y I m,e =E qi°,C • y E yi=1 Similarly, Eqs. (80) to (86) determine the discrete values of the cross-sectional characteristics for the frame beams and Eqs. (87) to (96) for purlins. i e I (71) Iz, p =E• yk k e K (93) k i e I (72) It, P =E qktP • yk k e K (94) k ieI (73) Im, P=E q'r • yk k eK (95) k ieI (74) E yk=1 (96) k i e I (75) 3.4. Sets, input data (constants) and variables ieI (76) The following sets, input data (constants: ieI (77) scalars and parameters) as well as continuous and binary variables are involved in the optimization ieI (78) model FRAMEOPT: (79) Sets: AB=E qjAB • yj j h =E j • yj j bB=E qbB • yj j tw, b =E jB • yj j tf, B =E qfB • yj j h. B=1 j • yj j E y=1 j AP =E qkAP • yk k K=E qhP • yk k bP=E qkbP • yk k tw, P=E qkw p • yk k tf, P =E qkf p • yk k Iy, P=E qky P • yk j e J j e J j e J j e J j e J j e J k e K k e K ke K ke K ke K ke K (80) (81) (82) (83) (84) (85) (86) (87) (88) (89) (90) (91) (92) i set for the standard dimension alternatives for columns, ieI j set for the standard dimension alternatives for beams, jeJ k set for the standard dimension alternatives for purlins, keK m set for the number of purlins, meM n set for the number of portal frames (columns and beams), neN Scalars (constants, input data): f denotes the overheight of the frame beam [cm] fy yield the strength of the structural steel [kN/cm2] h height of the column [cm] k effective length factor [-] k?, effective length factor [-] mr mass of the roof plates [kg/cm2] s snow (variable imposed load) [kN/cm2] wv vertical wind (variable imposed load) [kN/cm2] wh horizontal wind (variable imposed load) [kN/cm2] C1,C2 equivalent uniform moment factors [-] E elastic modulus of steel [kN/cm2] G shear modulus of steel [kN/cm2] L frame span [m] Ll length of the industrial building [m] minimum number of defined portal frames [-] maximum number of defined portal frames [-] minimum number of purlins [-] M7NNOframe MAXNÜframe MINNÜpurlin k MAXNOpullla maximum number of purlins [-] ab imperfection factor [-] aLT imperfection factor [-] yq partial safety factor for the variable load [-] 7g partial safety factor for the permanent load [-] yM0 resistance partial safety factor [-] distribution factor [ 1 ] !1 slenderness [-] к Ludolfs number [-] p density of steel [kg/m3] Parameters (constants, input data): q AC vector of i, ieI, discrete standard constants for cross-section area of the column q AB vector of j, jeJ, discrete standard constants for cross-section area of the beam q AP vector of k, keK, discrete standard constants for cross-section area of the purlin q bC vector of i, ieI, discrete standard constants for overall breadth of the column q bB vector ofj, jeJ, discrete standard constants for overall breadth of the beam q bkP vector of k, ke K, discrete standard constants for overall breadth of the purlin q itf, C vector of i, ie I, discrete standard constants for flange thickness of the column q j B vector of j, j e J, discrete standard constants for flange thickness of the beam q tkf, P vector of k, ke K, discrete standard constants for flange thickness of the purlin q itw, C vector of i, ie I, discrete standard constants for web thickness of the column q t"- B vector of j, jeJ, discrete standard constants for web thickness of the beam q tkw, P vector of k, ke K, discrete standard constants for web thickness of the purlin q iIt, C vector of i, ie I, discrete standard constants for torsional constant of the column q Ikt, P vector of k, ke K, discrete standard constants for torsional constant of the purlin q ^ C vector of i, ieI, discrete standard constants for second moment of area about the y - y axis of the column q j'B vector ofj, jeJ, discrete standard constants for second moment of area about the y - y axis of the beam q Iky, P vector of k, ke K, discrete standard constants for second moment of area about the y - y axis of the purlin q j4 C vector of i, ieI, discrete standard constants for moment of area about the z - z axis of the column q Ikz, P vector of k, ke K, discrete standard constants for moment of area about the z - z axis of the purlin q S C vector of i, ieI, discrete standard constants for warping constant of the column q k"P vector of k, keK, discrete standard constants for warping constant of the purlin Continuous variables: bB overall breadth of the beam [cm] bC overall breadth of the column [cm] bP overall breadth of the purlin [cm] ef intermediate distance between the portal frames [cm] ep intermediate distance between the purlins [cm] g self-weight of the portal frame [kN/cm] hB cross-sectional height of the beam [cm] hC cross-sectional height of the column [cm] hP cross-sectional height of the purlin [cm] qz uniformly distributed horizontal surface variable load [kN/cm] qy uniformly distributed vertical surface variable load [kN/cm] tf,B flange thickness of the beam [cm] tfC flange thickness of the column [cm] tf,P flange thickness of the purlin [cm] tw,B web thickness of the beam [cm] tw,C web thickness of the column [cm] twP web thickness of the purlin [cm] AB cross section of the beam [cm2] AC cross section of the column [cm2] AP cross section of the purlin [cm2] Дс torsional constant of the column [cm4] It,P torsional constant of the purlin [cm3] Iy,B second moment of area about the y - y axis of the beam [cm4] IyC second moment of area about the y - y axis of the column [cm4] IyP second moment of area about the y-y axis of the purlin [cm4] IzC second moment of area about the z - z axis of the column [cm4] IzP second moment of area about the z-z axis of the purlin [cm4] I