TRANSACTIONS ON ELECTRICAL ENGINEERING

CONTENTS Zablodskiy, N., Lettl, J., Pliugin, V., Buhr, K., Khomitskiy, S.: Induction Motor Design by Use of Genetic Optimization Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65 – 69

Homišin, J.: New Ways of Controlling Dangerous Torsional Vibration in Mechanical Systems . . . . . . . . . . . . . . . . .

70 – 76

Folta, Z., Hrudičková, M.: Experience with Torque Measurement on Rotating Shaft . . . . . . . . . . . . . . . . . . . . . . . . . .

77 – 81

Veleba, J.: Performance of Steady-State Voltage Stability Analysis in MATLAB Environment . . . . . . . . . . . . . . .

82 – 88

Ondrášek, J.: Cross Slide Mathematical Model for Solving Chatter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89 – 96

Vol. 2 (2013)

No.

3

ERGO NOMEN

pp.

65 - 96

TRANSACTIONS ON ELECTRICAL ENGINEERING Publisher:

ERGO NOMEN, o.p.s., K13114 FEE CTU in Prague, Technicka 1902/2, 166 27 Praha 6, Czech Republic E-mail: [email protected]

Editorial Office:

PIVONKA Pavel BAUER Jan HAVLICEK Radek MERICKA Jiri NOVA Ivana VONDRICH Jiri ZDENEK Jiri

Periodicity: Language: Scope: On-line version:

Quarterly English International scientific journal of electrical engineering www.transoneleng.org ISSN 1805-3386

Each paper in the journal is evaluated by two reviewers under the supervision of the International Editorial Board. International Editorial Board Editor in Chief: Prof. LETTL Jiri, Czech Technical University in Prague, Czech Republic Members: Prof. BAUER Palo, Delft University of Technology, Netherlands Prof. BRANDSTETTER Pavel, VSB-Technical University of Ostrava, Czech Republic Prof. DOLEZEL Ivo, The Academy of Sciences of the Czech Republic, Czech Republic Prof. DUDRIK Jaroslav, Technical University of Kosice, Slovakia Prof. NAGY Istvan, Budapest University of Technology, Hungary Prof. NOVAK Jaroslav, University of Pardubice, Czech Republic Prof. ORLOWSKA-KOWALSKA Teresa, Wroclaw University of Technology, Poland Prof. PEROUTKA Zdenek, University of West Bohemia, Czech Republic Prof. PONICK Bernd, Leibniz University of Hannover, Germany Prof. RICHTER Ales, Technical University of Liberec, Czech Republic Prof. RYVKIN Sergey, Russian Academy of Sciences, Russia Prof. SKALICKY Jiri, Brno University of Technology, Czech Republic Prof. VITTEK Jan, University of Zilina, Slovakia Prof. WEISS Helmut, University of Leoben, Austria

Responsibility for the contents of all the published papers and technical notes is upon the authors. Template in MS WORD and basic typographic rules to be followed see www.transoneleng.org.

Copyright:

©2013 ERGO NOMEN, o.p.s. All right reserved.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

65

Induction Motor Design by Use of Genetic Optimization Algorithms prof. N. Zablodskiy 1), prof. J. Lettl 2), doc. V. Pliugin 3), ing. K. Buhr 4), stud. S. Khomitskiy 5) 1)

Donbas State Technical University/Automation of electro-technical systems, Alchevsk, Ukraine, [email protected] Czech Technical University in Prague/Faculty of Electrical Engineering, Prague, Czech republic, [email protected] 3) Donbas State Technical University/Automation of electro-technical systems, Alchevsk, Ukraine, [email protected] 4) Czech Technical University in Prague/Faculty of Electrical Engineering, Prague, Czech republic, [email protected] 5) Donbas State Technical University/Automation of electro-technical systems, Alchevsk, Ukraine, [email protected] 2)

Abstract — The problem of the automated calculation and optimal design of an induction motor is presented . The problem of optimization by use of genetic algorithms is set and solved. The analysis of the obtained results is executed. Keywords — induction motor, optimization, varied variables, genetic algorithm, criteria, limitations, effective variant, EvoJ library

I. INTRODUCTION Neuron networks, being one of perspective trends of researches in the artificial intelligence area, as a result of watching processes going on in the nervous system of man were created. Approximately by the same way genetic algorithms were also «invented», but watched over the man nervous system, by the process of living organisms evolution. Genetic algorithms - one of research trends in the artificial intelligence area, engaging in creation of the simplified evolution models of living organisms for the of optimization task decision [1]. A classic genetic algorithm (GA) consists of following steps: 1) initializing, or choice of initial chromosomes population; 2) an estimation of chromosomes adjustment in a population - calculation of adjusted function for every chromosome; 3) verification of algorithm stop condition; 4) chromosomes selection - choosing of chromosomes, participated in descendants for a next population creation; 5) application of genetic operators are mutations and crossing; 6) forming of new population; 7) choosing of the «best» chromosome. The block-diagram of GA is represented in Fig. 1. A simple GA generates an initial population by a random way. Working of the GA is an iteration process which proceeds until the generations set number or some another stop criterion will not be executed. On every

generation, a proportional selection on adjusted, crossing and mutation is realized. The simplest proportional selection is roulette. The wheel of roulette contains one sector for every member of population. The size of every sector is proportional to the corresponding size of adjusted function. At such selection, members of population with higher adjustment will be chosen with greater probability than individuals with subzero adjustment. The next step is using of crossing and mutation. A previous population, obtained after a mutation, is overwritten and the cycle of one generation is completed. Subsequent generations i.e. selection, crossing and mutation obtained as a GA working result are processed in the same way .

II. THEORY AND PROGRAM REALIZATION In the examined task a GA provides one criterion of optimality only, by virtue of the program realization of a calculation function minimum search [2 - 4]. The order of optimization will be different from considered Carthesian product (CP) in the previous article [5]: 1) 2) 3) 4)

setting the range of the varied variables; setting limitations; choosing the criterion of optimality; calling the CA function for optimization and getting the optimal varied variables set; 5) calling the function of the induction motor (IM) automatic calculation for the found set. We will consider an example of the IM optimization programmatic realization using GA on Java in NetBeans IDE [2]. For a decision the problem we will use free Javalibrary EvoJ (Evolution Java) [3]. A project EvoJ is designed as upgradable framework of Java classes for the decision of various optimization tasks by use of evolutional (genetic) algorithms. For the use of EvoJ a programmer must implement one simple interface, consisting of one method only. All other steps undertake EvoJ algorithm.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

66

START

Initialization – choosing of initial chromosome population

Estimation of chromosomes and population adjustment

YES

NO Algorithm stop condition achieved?

The best chromosome choosing Chromosome selection END Applying of genetic operators

Creating of new population

Fig. 1. Magnetization as a function of applied field.

In the example two varied variables will be considered: the internal diameter of the stator core and the length of the stator core. We shall create Java-interface with the name “Solution” in which we set the range (minimum and maximum values) of the varied variables. The code of the Solution interface is down in the text. EvoJ is able to change the variables without a setting of a range. However if it is needed to implement an own mutation strategy, one have to declare setters. In other case we shall not have a possibility to change the variable range. Pay attention to annotation @of Range - it sets the range of values which a variable can accept. Variables of random values from the set range are initiated.

package MotorClasses; import net.sourceforge.evoj.core.annotation.MutationRange; import net.sourceforge.evoj.core.annotation.Range; public interface Solution { //Diameter of stator core String smin1 = "165"; String smax1 = "205"; //Length of stator core String smin2 = "115"; String smax2 = "145"; @Range(min = smin1, max = smax1) double getX(); //return the optimal diameter @Range(min = smin2, max = smax2) double getY(); //return the optimal length }

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

However as a result of mutation they potentially can go out from the indicated range. It can be prevented using the parameter of strict=«true», that will not allow a variable to take on an impermissible value, even if we make an effort to propose them using setter-method. Another case on which it is needed to pay attention here, is that all parameters of all limitations in EvoJ are strings. It allows both to specify the value of parameter directly and to specify the name of property instead of concrete value, to specify the value of limitation parameters at the compile-time. Now we have an interface with variables and we shall write a fitness-function. Fitness-function in EvoJ is implemented as the following interface: public interface SolutionRating { Comparable calcRating(T solution); }

67

@Override public Comparable doCalcRating(Solution solution){ double fn = calcFunction(solution);//call motor function boolean flag = mot.control();//control of limitations if (Double.isNaN(fn) | flag == false){ return null;//sift-out false variant } else { return - fn;//return an effective variant } } }//end of class

All code lines above are obviously enough. We simply take and count our function, using variables from the Solution interface: double fn = calcFunction(solution);//call motor function

Here parameter is our interface with variables. The greater value return (according the contract Comparable) the more suitable solution is considered. Null can be returned and it is the smallest value of a fitness-function. It is recommended to implement this interface as mediated using helper-classes. They undertake some service functions: elimination of the old decisions (if the maximal life term of decision is set), cashing of function value for decisions which were not sifted from in the previous GA iteration. A Fitness-function for our case will look like the following (we create a new class with the name Rating): package MotorClasses; import net.sourceforge.evoj.strategies.sorting.AbstractSimpleRating; public class Rating extends AbstractSimpleRating { static AMotor mot;//motor object static int krit;//index of optimality criterion static int iter_numb;//number of iterations //Constructor public void set_motor(AMotor mot, int krit){ this.mot = mot; this.krit = krit; this.iter_numb = 0; } //reception of iterations number public int get_iter(){return this.iter_numb;}; public static double calcFunction(Solution solution){ iter_numb++;//increase of iterations count double x = solution.getX();//reception of new diameter double y = solution.getY();//reception of new length mot.stator.set_D(x/1000);//setting of new diameter mot.stator.set_ld(y/1000);//setting of new length double fn = mot.auto(krit);//automatic motor calculation return fn;//return a criterion of optimality }

Because we search minimum and the contract of class supposes that the best decisions must have a greater rating, we return the value of function, multiplied by 1: return - fn;//return an effective variant

The population with the highest value of variable fn will be considered as the most close to optimum result. In addition, we sift-out false decisions (if NaN turned out or motor limitations were not passed), returning null. Override of function Comparable realize the mechanism of genetic populations, using as the achieved result the value returned by the function calcFuction(). In the IM class Motor we create the function of automatic calculation, which accepts as an argument the index of optimality criterion and returns the got criterion after the motor calculation: double auto(int krit){ int res = 0; //Code body of motor calculation //… switch (krit){ case 1://1 is efficiency res = 1/kpdnr; break; case 2://2 – power factor res = 1/cosFinr; break; case 3://3 – starting current res = I1pn; break; case 4://4 – starting torque res = 1/Mpo; break; } return res;//return a criterion depending on its index }

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

68

Further in the motor class Motor we create the function of GA realization (a code structure is explained in comments):

It is obvious from Tab. I, an optimal motor efficiency is higher than a base value and other parameters are satisfying limitation range.

void optimization(int krit){ DefaultPoolFactory pf = new DefaultPoolFactory(); //creation of populations with the amount “populations” GenePool pool = pf.createPool(populations, Solution.class, null); Rating rtg = new Rating();//Constructor of Rating class rtg.set_motor(this, krit);//getting of Motor object and criterion //factory of initial decisions set generation DefaultHandler handler = new DefaultHandler(rtg, null, null, null); //implementation of iterations number “iterations” //over the population “populations” handler.iterate(pool, iterations); //reception of the best found decision Solution solution = pool.getBestSolution(); D_opt = solution.getX()/1000; //optimal diameter L_opt = solution.getY()/1000; //optimal length int iter = rtg.get_iter(); //reception of iterations count number }

TABLE I GENETIC ALGORITHM: TABLE OF PARAMETERS BEFORE AND AFTER OPTIMIZATION Name Base value Optimal value Induction in the air-gap, T 0.748 0.807 Internal diameter of stator core D, mm 185 194 130 115 Length of stator core Lδ, mm 0.895 0.755 Relative size λ = Lδ/τ (τ = πD/2p, where p - number of pole pairs) Height of stator slot, mm 21.9 14.6 Height of rotor slot, mm 32.2 33.2 Width of the upper line of stator slot, 7.7 7.8 mm Width of the down line of stator slot, 10.2 9.3 mm Upper diameter of rotor slot, mm 7.9 7.8 Down diameter of rotor slot, mm 3.7 3.4 Efficiency 0.885 0.891 Power factor 0.893 0.9 Starting current relative value 5.84 6.52 Starting torque relative value 1.4 1.62 Overload torque capability 2.65 2.88 Overall stator winding temperature, C 93.25 95.69

From a NetBeans form, the code of the GA optimization implementation consists of two lines: motor.limits();//setting of limitations motor.optimization(krit);//optimization with the criterion “krit”

Functions limits() consist restrictions on motor geometric sizes, temperature limits, starting currents and etc. In the fitness-function limitations are checked by the control function boolean flag = mot.control(); //control of limitations

This function return false if even one restriction will be broken. In control() function there are 16 motor variables limits have been set. In particular we can avoid high temperatures of the stator that is probably resulted because of increasing of the stator winding cross-section. If a decision will not arranged (do not satisfy to motor restrictions according to limitation function) it is possible to continue the GA iterations (increasing the populations number “populations” and iterations “iterations”), while the desired quality of decision will not be attained. So to solve a GA task using EvoJ it is necessary: 1) to create an interface with variables; 2) to implement the interface of the fitness-function; 3) to create the population of decisions and carry out the necessary amount of the GA iterations above them, using a code, given above. The results of the GA solution at a choice of maximum efficiency as a criterion of optimality are shown in Tab. 1.

III. CONCLUSIONS Algorithm of the previous considered CP [5] in comparison with the GA, allows to execute multi-criterion optimization, that is its undoubted advantage. In addition, the CP always gives only the synonymous best variant among the existing ones. However, in the range of varying of two variables ± 20 % from a base value (3976 combinations) the calculation time is approached up to 48 min. Implementation of CA gives stunning results. At the same varied variables and range of their change ± 100 % (!) from a base value, the calculation time is only 40 sec! However, GA, at least, in the present article task, does not allow to execute optimization for a few criteria. In the GA number of the varied variables and a range of their change is not important from the point of view of the productivity, because a set of the varied variables is created dynamically, but not beforehand, as in the CP method. In addition, all combinations of the variables and values of objective function are realized in a binary form. However, time of the GA work is very critical to the number of the created populations and number of iterations in the populations. The choice of population’s number and iterations realized by an experienced way increases until an acceptable result will not be obtained. A result of the optimization with the use of the GA will always be the best for the chosen criterion, but there is not a guarantee, that a better variant can exist. Actually, herein there is a genetic selection logic - we get a result, approaching the best among the random created populations. A most reliable result depends on population’s number. Thus, the degree of authenticity can be estimated by the variation of

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

the obtained results in repeated calculations at the same population’s number. Therefore, the GA productivity at effective variant populations number is determined but not by the number and range of the varied variables. The result of the optimization is ambiguous and close to the best. When production of approximate calculations in maximum compressed terms is needed and quality of the obtained results is written with a permissible error, then using of the GA optimization will be the irreplaceable instrument for designers.

69

REFERENCES [1] Yemelyanov V.V., Kureychik V.V., Kureychik V.M. Theory and practice of evolution modeling. М.: PHYSMATLIT, 2003. – 432 p. [2] N. Zablodskiy, V. Pliugin, K. Buhr. CAD of electromechanic devices: educational tutorial, part 2, 2013. - 330 p. (will be printed). [3] http://evoj-frmw.appspot.com [ONLINE]. [4] N.K. Vereshchagin, A. Shen. Lectures on mathematical logic and theory of algorithms. Beginning of sets theory. MCNMO, 2008. – 198p. [5] N. Zablodskiy, V. Pliugin, K. Buhr, S. Khomitskiy. Asynchronous motor optimal design with using of Cartesian product (will be printed).

REFERENCES ON RUSSIAN: [1] Емельянов В.В., Курейчик В.В., Курейчик В.М. Теория и практика эволюционного моделирования. М.: ФИЗМАТЛИТ, 2003. – 432 с. [4] Верещагин Н.К., Шень А. Лекции по математической логике и теории алгоритмов. Начала теории множеств. МЦНМО, 2008. – 198с.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

70

New Ways of Controlling Dangerous Torsional Vibration in Mechanical Systems Jaroslav Homišin Department of Design, Transportation and Logistics– Section of Design and Machine Elements Technical University of Košice, Faculty of Mechanical Engineering, Letná 9, 040 01 Košice Slovak Republic e-mail: [email protected] Abstract — In general terms the mechanical systems (MS) mean the systems of driving and driven machines arranged to perform the required work. We divide them into MS operating with constant speed and MS working within a range of speed. In terms of MS dynamics we understand the system of masses connected with flexible links between them, i.e. systems that are able to oscillate. Especially piston machines bring heavy torsional excitation to the system, which causes oscillation, vibration, and hence their noise. Based on the results of our research, the torsional vibration control as a source of given systems excitation, can be achieved by application of pneumatic couplings, or pneumatic tuners of torsional vibration. Existence of pneumatic couplings and pneumatic tuners of torsional vibration developed by us, creates the possibility of implementing new ways of torsional oscillating mechanical system tuning. Based on the above, the aim of the article is to present new ways of dangerous torsional vibration control of mechanical systems by application of pneumatic couplings and pneumatic tuners of torsional vibration developed by us. Keywords — torsion, oscillating mechanical system, pneumatic coupling, pneumatic tuner of torsional vibration, ways of torsional vibration control

I. INTRODUCTION Any MS, in terms of dynamics, we understand the system of masses connected with flexible links between them, i.e. systems that are able to oscillate. Piston machines, either drivers or driven units, bring to the system significant torsional vibration. This means that MS with internal combustion engines, compressors and pumps can be characterized as torsional oscillating mechanical systems (TOMS). In the range of operating speed there can be a very intense resonance between the driver frequencies (reciprocating machines) and the natural frequencies of the system. Consequently, there comes to an excessive vibration and related excessive stress of the whole MS. The excessive dynamic stress often causes malfunction of various parts of the system, such as: fatigue fractures of shafts, gear failures, deformation failures of flexible couplings and the like. Therefore, it is necessary to control their dangerous torsional vibration. Currently, the torsional vibration is reduced to a permissible degree by appropriate adjustment, respectively tuning the system by application of an appropriately selected flexible coupling, based on a dynamic calculation. Thus the principle of tuning is an appropriate adaptation of the basic dynamic properties, particularly the dynamic torsional stiffness of the flexible coupling to the system.

The general characteristics of flexible couplings include their dynamic torsional stiffness and damping coefficient. It should be noted that, they are affected by material (metal, rubber, plastics), shape, number and dimensions of the flexible elements. It follows that they depend on various factors, which are divided into stable and unstable factors [1]. The shape, number, size and various structural modifications of flexible elements can be added to the group of stable factors, while the material of flexible members to the group of unstable factors, as a result of fatigue and aging which are changing their original characteristics. By changing original properties there is a change of coupling characteristics Mk = f(φ) (with respect to initial characteristics), and thus a change of its basic characteristics, which has a largely positive impact on the magnitude of the dangerous torsional vibration of the mechanical system [1], [2], [3]. It should be emphasized that any linear or nonlinear flexible coupling is only one characteristic, tightly coupled to the used flexible element. In the case of a linear coupling it is only one characteristic of a constant dynamic torsional stiffness. Dynamic torsional stiffness of the nonlinear coupling varies in some extent of its characteristics, obviously dependent on the working mode of the system. Changing the characteristics of coupling due to appropriate dynamic tuning of TOMS means the use of an other coupling flexible element or other flexible shaft coupling. Influences such as: temperature of flexible coupling elements and the number of cycles causes that, by effect of external forces any flexible member of the coupling is exposed to fatigue and aging. Consequently, there is a change of coupling characteristics, and thus a change of its basic characteristic properties. This leads to the fact that a suitably tuned TOMS becomes untuned. A flexible coupling in this case does not act as a tuner, but rather as a driver of torsional oscillations. It should be noted that this method of tuning will be suitable only in cases where no previously unforeseen (random) effects occur during the operating mode, particularly in the turbo-machinery and reciprocating machinery [2], [3]. In case of random failure in an operating mode of MS a very intense resonance of lower harmonic excitation occurs, which is usually unexpected. Due to this fact, an intense torsional excitation causes increased torsional vibration, mechanical vibration and hence a noisy mechanism. The torsional vibration control, based on the results of our research, is achieved by use of a pneumatic flexible coupling as well as by application of pneumatic flexible

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

coupling with auto-regulation – pneumatic tuner of torsional systems. The aim of this scientific paper will be presentation of the possibility to control dangerous torsional vibration of mechanical systems by ways suggested by us, by tuning and continuous tuning. Tuning of the given system is provided by application of the pneumatic flexible coupling while continuous tuning (tuning during operation) of the system will be ensured by application of the pneumatic flexible coupling with auto-regulation – pneumatic tuner of torsional vibration.

71

sure p of the gas media in it. Varying pressure will ensure the fact that the coupling always works with different characteristics (Fig. 2), which is defined by the formula (1).

II. PROPOSED WAYS TO TORSIONAL VIBRATION CONTROL OF MECHANICAL SYSTEMS

A change of the pneumatic couplings torsional stiffness can be realized by changing the pressure of gaseous media, out of operation (Fig. 3) or during operation (Fig. 5) of mechanical systems. This leads to two proposed ways to the torsional vibration control of mechanical systems [4]: the torsional vibration control of mechanical systems out of operation, ensuring the so called tuning of the system [5], [6], [7], [8], the torsional vibration control of mechanical systems during operation, ensuring the so called continuous tuning of the system [5], [9], [11]. Under the tuning of torsionally oscillating mechanical systems with pneumatic coupling we understand the inflation space of the coupling compression suitable to pressure value of the gaseous medium out of the operation system. The appropriate pressure value of the gaseous medium, and hence the appropriate value of the dynamic torsional stiffness of the coupling is based on the previously realized dynamic systems in terms of calculation of the torsional dynamics. The mechanical systems run during their entire operation with such inflated pneumatic coupling. The principle of the torsional vibration control of the mechanical system during its operation at steady state by application of torsional oscillations pneumatic tuner [2], [3] shows the adaptation of the basic dynamic properties, particularly the dynamic torsional stiffness of the tuner to the system dynamic . The basic principle of the pneumatic tuner is the ability to auto-regulate the twist angle due to a current change of the load torque on a predetermined constant angular value φK . This will ensure the autoregulation of gas pressure in the compression space of the tuner, thus adapting it to the current value of the load torque. III. CHARACTERISTICS OF PNEUMATIC FLEXIBLE SHAFT COUPLINGS

The differential pneumatic coupling (Fig.1) consists of the driving part (1), driven part (2), between them there is located the compression space filled with gaseous medium (air in our case). The compression space consists of three circumferentially spaced and interconnected differential elements. Each differential element consists of a compressed (3) and expanded pneumatic-flexible element (4). Interconnection of differential elements is provided by the interconnecting hose (5). The filling of compression space of coupling through the valve (6) changes the pres-

Fig. 1. Differential pneumatic flexible shaft coupling.

To another characteristic other characteristic properties always belong, thus still different torsional stiffness and damping coefficient. Therefore, each pneumatic coupling depending on the pressure is always defined by another course of the torsional stiffness in Fig. 3, as described by the formula (2). (1) M = a0 .ϕ + a3 .ϕ 3 , k = a0 +

3 a3 .ϕ 3 , 4

(2)

where: ϕ – twist angle of the coupling, a0, a3 – constants of the coupling characteristics. The orsional stiffness, as the main component in the field of the mechanical system tuning has a decisive influence on the natural frequency of the system

Ω0 = k / I red ,

(3)

where: Ired – reduced mass moment of inertia of the mechanical system. It therefore follows the basic principle of the mechanical system tuning by pneumatic couplings. Its basic principle is to customize the natural frequency of the system Ω0 to the angular excitation frequency ω, so that in the range of the system working mode there is no resonance condition ω = Ω0 and hence dangerous torsional vibration. The pneumatic tuner of the torsional vibrations (Fig. 4), which basic principle results from the patent claims [9], [11] is compared with the differential pneumatic coupling on a common structural base. The main difference is that it does not have the valve, but the controller (6) to ensure a coupling constant twist angle φk. The basic principle of the tuner is the ability to auto-regulate the twist angle due to the torque current load change on a predetermined constant angular value φk. This will ensure auto-regulation of the gaseous media pressure in the compression space of the tuner, thus adapting it to the current value of load torque.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

72

In Fig. 5 there are in general terms shown the traces of the pneumatic tuner of the torsional vibrations and torsional stiffness depending on the load torque. To each constant twist angle φk1, φk2, φk3 and φk4 based on a calculation one course of torsional stiffness labelled a, b, c, d is given.

Fig. 2. Courses of static characteristics of the differential pneumatic coupling a, b, c, d, e, f, g shown in the general version belong to pressure of the gaseous medium at p = 100 ÷ 700 kPa with 100 kPa.

Obr. 5. Courses of the pneumatic tuner of the torsional vibration and torsional stiffness depending on the load torque M.

IV. THE RESULTS OF THE INVESTIGATION OF THE PROPER Fig. 3. Courses of the torsional stiffness k of the differential pneumatic coupling a, b, c, d, e, f, g shown in the general version belong to pressure of the gaseous medium at p = 100 ÷ 700 kPa with 100 kPa.

Fig. 4. Pneumatic tuner of the torsional vibration.

Auto regulation of he pressure of the gaseous media has a direct effect on the characteristics change of the pneumatic tuner (Fig. 2) Of course, for change of the torsional stiffness value (Fig. 5), as a result, we can tune the natural frequency of the system.

TUNING AND CONTINUOUS TUNING OF THE TORSIONALLY OSCILLATING MECHANICAL SYSTEM

When investigating an appropriate tuning, or any continuous tuning of the torsionally oscillating mechanical system we mostly start from the Campbell diagram showing the position of the critical speed nK (or the position of the critical angular frequency ωK) depending on the rotational speed frequency N (or natural angular frequency Ω0). Magnitude of the torsional vibration for the tuning and continuous tuning of the system is mostly presented by: courses of the dynamic torque amplitude excited by the torsional vibration in the mechanical system and hence to the pneumatic coupling, depending on the speed. A. Characteristics of the realized torsionally oscillating mechanical system The examination of an appropriate tuning and continuous tuning was performed on a realized torsionally oscillating mechanical system (Fig. 6). The realized system is composed of the driving part (1), pneumatic flexible shaft coupling (3) and driven part (2). The driving part, formed by a DC electric motor type SM 160 L with a power of 16 kW and an auxiliary thyristor controller of the rotational frequency (4) type IRO with the possibility of continuous control from n = 0 ÷ 2000 min-1, using a pneumatic coupling that drives the exciter of the torsional vibrations represented by the three-cylinder compressor type 3–JSK–S. To increase the impact of the compressor tor-

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

sional excitation to the mechanical system we used a compressor without a flywheel. The load of the torsional oscillating mechanical system by the compressor will be adjusted (regulated) with throttle valve (6) integrated into the output pipe of the compressor. The analyzis the pneumatic coupling load at work of the mechanical system in the steady state will be investigated by a schematic model of the realized torsional oscillating mechanical system (Fig. 7).

7

6

73 I 1 .ϕ&&1 + b.(ϕ& 1 − ϕ& 2 ) + k .(ϕ 1 − ϕ 2 ) = M i . sin .(i.ω .t + γ i ),

I 2 .ϕ&&2 − b.(ϕ& 1 − ϕ& 2 ) − k.(ϕ 1 − ϕ 2 ) = 0. n

M d = ∑Mi. i =1

4

I2 .θ . sin[(i.ω .t + γ i ) + β i + ϑi ] I1 + I 2

(6)

while dynamic coefficients θ , ξ and phase angles βi , υi are described by formulas

θ=

3

(5)

ξ=

 i.ω 1 +   Ω0   i.ω 1 −    Ω 0

  

2

  

2

2

 2χ .  Ω0

  i.ω  +    Ω 0

  

   2

2

 2χ .  Ω0

  

1   i.ω 1 −    Ω 0

  

2

2

  i.ω  2  2 χ  .  +    Ω 0   Ω 0

,

(7)

,

(8)

2

  

2

Where for the damping coefficient 2.χ, natural angular frequency of the system Ω0 and separation margin η it is applied

5

2

1 1 2.χ = b. + ,  I1 I 2 

1

1 1 Ω 0 = k . +  ,  I1 I 2 

η=

i.ω i.n = . (9) Ω0 N

B. Results of a proper tuning of the torsionally oscillating mechanical system The Campbell diagram according Fig. 8 describes the tuning of the realized mechanical system by the tangential pneumatic coupling type name 4–1/70–T–C in the speed range n = 0 ÷ 2000 min-1. Operating mode of the system is defined in the speed range n = 750 ÷ 1500 min-1.

Fig. 6. Realized torsional oscillating mechanical system.

Fig. 7. Schematic model of the realized torsional oscillating mechanical system.

When calculating the loads for the steady-state mechanical system within its operational mode, we suppose that the mechanical system is rotating at an angular speed, which varies in a wide range. On mass (1) with the mass moment of inertia I1 a load torque acts . From the above it is clear that M N + ∑ M i . sin(i.ω .t + γ i ) pneumatic coupling and thus the whole torsional oscillating mechanical system is loaded by both unvariable with time medium torque MN in steady state and excitation of harmonic components Mi. On this basis an additional component of dynamic torque Md is introduced in the pneumatic coupling. Thus pneumatic coupling will be in this case loaded with load torque MS that causes the maximum twist angle φS: MS = MN + Md,

(4)

The magnitude of the additional dynamic torque calculated from the equations of motion (5) can be described by equation (6).

Fig.8. The Campbell diagram of the mechanical system with the applied tangential pneumatic coupling at the constant pressure p = 100 ÷ 700kPa.

The diagram shows the position of the critical speed, depending on the natural speed frequencies. Those given pneumatic couplings are nearly linear, natural speed frequencies are shown by the horizontal straight lines a, b, c, d, e, f, g for the entire range of the gaseous medium pressure p = 100 ÷ 700 kPa. Based on the diagram it is possi-

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

ble to say that the pneumatic coupling is capable to operate at all pressures of the gaseous medium (p = 100 ÷ 700 kPa). On the other hand, in terms of the dynamic tuning we conclude that the pneumatic coupling is suitable for the given system in the pressure range p = 200 to 600 kPa. This is due to the fact that at the beginning of the operation mode at p = 100 kPa there is a resonance with the first harmonic component of the load torque speed n = 820 min-1, while at p = 700 kPa a resonance occurs also with the first harmonic component but at the speed n = 1480 min-1. Based on the above we shall in following focus on the dynamic tuning characteristics of the realized system by the tangential pneumatic coupling at the pressure range p = 200 to 600 kPa. Based on the figure of the Campbell diagram it is possible to say that a given coupling extends from the range of operating speed harmonic series i = 2 ÷ 12 The critical speed due to the main harmonic component (i = 3) at pressures p = 200, 300, 400, 500 and 600 kPa appears at the speed nK = 330, 360, 405, 440 and 460 min-1. These values indicate that the realized mechanical system is welltuned with regard to the operating mode beginning. This is confirmed by the separation margin η = i. n / N, which for the investigated pressures has relatively high values of η = 2,3 to 1,6 for the investigated pressures. At the same time we can see in the figure that the harmonic series i = 1 extends in the range of the gaseous medium pressure p = 200 to 600 kPa into operating speed range (OSR). It follows that when using the pneumatic coupling there is a resonance with a harmonic component at these pressures. Specifically, for the pressure p = 200, 300, 400, 500 and 600 kPa resonances occur at speeds nK = 980, 1090, 1220, 1330 and 1430 min-1. The tuning system realized by the pneumatic coupling for one disabled cylinder with regard to the main (i = 3) and secondary (i = 2, 1) harmonics within the operating speed (n = 750 ÷ 1500min-1) is characterized in Fig. 9.

Fig. 9. The dependence of the dynamic component of the torque Md at speeds in the range n = 0 ÷ 2000min-1 of the mechanical system on the tangential pneumatic coupling application with constant pressure at p = 100 ÷ 700 kPa.

74

The overall analysis shows that differential pneumatic coupling can be applied in the torsional oscillating mechanical system with a range of speed only in a fault-free case of work of the piston device. In case of faults caused mainly by the piston device (unbalanced excitation of engine cylinders, one disabled cylinder) it is to use the linear coupling in mechanical systems with a range of unsuitable speed . The reason of this is that, in this case, particularly lower harmonics cause increased amplitude of the mechanical load across the system (Fig. 9). The results indicate that the linear differential pneumatic coupling would be particularly suitable for the mechanical system operating with constant operating speed. C. The results of the continuous tuning of the torsionally oscillating mechanical system The result of continuous tuning of the system realized by a pneumatic tuner of the torsional oscillation type 4– 1/70–T–C is presented by the Campbell diagram in Fig. 10. In the figure there are represented eight waveforms of the natural speed frequencies marked a, b, c, d, e, f, g, h, corresponding to a constant twist angle of the pneumatic tuner φK = 0,5 °; 1 °; 1,5 °; 2 °; 2,5 °; 3 °; 3,5 °; 4 ° and are characterized by a broken line. Based on the Campbell diagram it is possible to say that critical speeds by φK = 0,5 ° and 1 ° from the main harmonic i = 3 of the load torque are in a sufficient distance with regard to the start of the operating mode (n = 750 min-1) by the lowest pressure p = 100 kPa and also by the highest pressure p = 700 kPa. This fact is confirmed by the separation margin, which for that case has values in the range η100 = 2,3 and η700 = 1,51. At the same time we can see that within the operating speed range of the system, particularly for the speed n = 1480 min-1, a resonance is at the harmonic component series of i = 1 by φK = 0,5° and 1° of the pneumatic tuner with the maximum pressure value of the gaseous medium p = 700 kPa. Based on the above it can be concluded that the constant twist angles of the pneumatic tuner φK = 0,5 ° and 1 ° are not suitable for the realized system. By the minimum pressure value of the gaseous medium of the pneumatic tuner p = 100 kPa with φK = 1,5 ° a resonance occurs at the harmonic component i = 1 at the operating speed n = 850 min-1 . With rising pressure up to the maximum value no resonance is at the harmonic component i = 1. For example at the maximum pressure the separation margin for i = 1 has a value η = 1,34. It indicates that the pneumatic coupling with φK = 1,5 ° is appropriate for the realized system except the beginning of the operating mode. When using the pneumatic tuner with constant angles φK = 2 °; 2,5 °; 3 °; 3,5 ° and 4 ° no resonance is within the operating speed range of the realized system from any harmonic components of the load torque. The results of the torsional vibration magnitude of the realized mechanical systems in the case of a disabled cylinder are shown in Fig.11. They are characterized by courses of the dynamic torque amplitudes Md depending on the operating speed in the range n = 0 ÷ 2000min-1.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

Fig. 10. The Campbell diagram of the realized mechanical system by use of the pneumatic tuner of the torsional vibration type 4–1/70–T–C.

It results from the overall analysis in Fig. 11 that the lowest dynamic loads of the mechanical system in the range of the operating speed n = 750 ÷ 1500min-1 , are obtained at constant twist angles of the pneumatic tuner of torsional vibration φK = 2°; 2,5°; 3°; 3,5° and 4°. It is caused by the fact that the pneumatic tuner in that case acts as a highly flexible pneumatic coupling, thus coupling with a relatively low torsional stiffness. The results indicate that the pneumatic tuner of torsional vibration will be especially suitable for mechanical systems working within a range of operating speed.

75

characteristic. The change of the flexible coupling characteristics, due to appropriate adaptation of its dynamic properties to the system dynamics means to use a different element of the flexible coupling or using a different flexible shaft coupling. In any case, it is not possible to forget the fatigue and aging of flexible materials, which finally have a major impact on the initial dynamic properties. Thus the unsteadiness of flexible coupling dynamic properties caused by aging and fatigue of their flexible elements and as well as the frequent failure rate of some other elements of the system causes the detuning of the tuned torsional oscillating mechanical system. In this case its tuning element, the flexible coupling, has no possibility to remove or reduce the increasing dangerous torsional vibration. Taking into account the given facts we propose to use the pneumatic flexible shaft couplings developed by us in order to reduce dangerous torsional vibration by optimal tuning or rather optimal continuous tuning of torsional oscillating mechanical systems. Based on the presented results it is possible to say that presented differential pneumatic coupling, as well as the pneumatic tuner of the torsional vibration, fulfil all the requirements for their application in torsional oscillating mechanical systems. Based on the detailed analysis of the realized mechanical system we can say that linear pneumatic couplings are especially suitable for mechanical systems operating with constant operating speed. On the other hand, the pneumatic tuners of torsional vibrations will fulfil all the requirements of mechanical systems within a range of operating speeds.

ACKNOWLEDGMENT This paper was written in the framework of Grant Project VEGA: „ 1/0688/12 – Research and application of universal regulation system in order to master the source of mechanical systems excitation”. REFERENCES

Fig.11. Courses of the dynamic torque amplitudes Md depending on the operating speed in the range n = 0 ÷ 2000min-1 for the realized mechanical system with application of the pneumatic tuner of torsional vibration

V. CONCLUSION Based on presented results we can say that negative impact of the dangerous torsional vibration is possible to reduce by application of classical flexible couplings. On this occasion it is necessary to note, that each linear or nonlinear presently used flexible coupling has only one

[1] J. Homišin a kol., “Súčasné trendy optimalizácie strojov a zariadení”, C–Press Košice, 2006, ISBN 80-7099-834-2. [2] J, Böhmer, “Einsatz elastischer Vulkan-Kupplungen mit linearer und progressiver Drehfeder-charakteristik”, MTZ, 44/5, 1983. [3] V. Zoul, V, “Torzní vibrace v pohonech a způsob jejich snižování”, Praha, ČSVTS 1984. [4] J. Homišin, J., “Methods of tuning torsionally oscillating mechanical systems using pneumatic tuners of torsional oscillations”, Transactions of the TU of Košice, 3/4, England, 1993, pp. 415 [5] J. Homišin, J., “Mechanická sústava optimálne vyladená pneumatickou spojkou”, UV SR/5274/2009. [6] J. Homišin, “Plynulo riadená mechanická sústava”. UV SR/ 5275/ 2009. [7] J. Homišin, “Pneumatická pružná hriadeľová spojka”. Patent č. 222411/86. [8] J. Homišin, “Pneumatická pružná hriadeľová spojka s diferenčnými členmi”. UV SR/5278/2009. [9] J. Homišin, “Regulačný systém pre zabezpečenie plynulej zmeny charakteristiky pneumatických spojok”. P ČSSR/259225/87. [10] J. Homišin, “Regulačný systém pre realizáciu plynulého ladenia mechanickej sústavy”. P SR/276927/92. [11] J. Homišin, “Pneumatická pružná hriadeľová spojka so schopnosťou autoregulácie”. P ČSSR 278025/95. [12] J. Homišin, M. Jurčo, “Aplication of differential pneumatic dutches voith and without autoregulation in torsionally oscillating me-

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

chanical systems”. The shock and vibration digest, 29/3, 1997, USA, pp.44, (80%/20%). [13] J. Homisin, M. Jurčo, “Application of differential pneumatic clutch with an additional regulating system”. The Shock and Vibration Digest, USA, 30/6, 1998, pp.490, (80%/20%). [14] J. Homišin, “Dostrajanie ukadów mechanicznych drgajacych skretnie przy pomocy sprzegie pneumatycznych. Kompendium wyników pracy naukowo-badawczych autora”. Bielsko-Biała, ATH, 2008, [106 p]. ISBN 978-83-60714-55-3. [15] R. Grega, “Prezentácia výsledkov dynamickej torznej tuhosti pneumatickej pružnej spojky s autoreguláciou na základe experimentálnych meraní”. Acta Mechanica Slovaca, 2/2002, ročník 6, s. 29 – 34. [16] P. Kaššay, M. Urbanský, “Úvod do problematiky prechodových dejov v torzne mechanických sústavách”. Zborník 51. MVK KČSaM, 2010. s. 130 – 136.

76

[17] P. Kaššay, R. Grega, J. Krajňák, “Determination of objective function for extremal control of torsional oscillating mechanical system“. Transactions of the Universities of Košice. Nr. 3, 2009, pp. 17–20. ISSN 1335-2334. [18] P. Kaššay, “Effect of Pneumatic Flexible Shaft Coupling on the Size of Torsional Vibration“. 2 Międzynarodowa Konferencja Studentów oraz Młodych Naukowców. Bielsko-Biała, Wydawnictwo naukowe Akademii techniczno-humanistycznejj, 2012 pp. 99‒104. ISBN 978-83-63713-23-2. [19] P. Lacko, V. Lacko, “Continuously Drisen Rezonance”. Strojárstvo 42 (3/4), 2000, pp. 127‒135, Zagreb, Croatia. [20] L. Pešík, P. Němeček, “Identification of the dynamic system of a machine with an elastic base”. McNU 97, Chicago, USA, 1997. [21] L. Pešík, “Aplikace převodového mechanizmu v úlohách vibroizolace strojů a zařízení”. Acta Mechanica Slovaca, ročník 6, 2/2002, s. 75–78.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

77

Experience with Torque Measurement on Rotating Shaft Zdeněk Folta1), Milena Hrudičková 2) 1)

VŠB-Technical University of Ostrava/Faculty of Mechanical Engineering, Ostrava, Czech Republic, e-mail: [email protected] 2) VŠB-Technical University of Ostrava/Faculty of Mechanical Engineering, Ostrava, Czech Republic, e-mail: [email protected] Abstract — This article deals with problem of torque measurement on the turning shaft and the signal transmission to the measuring equipment. There is described the experience with contact transmission and both non-contact and high frequency transmission (with ESA Messtechnik equipment). Also the examples of the practical realization are presented. Keywords — Torque measurement, telemetry, strain gauge measurement, experiment

I. INTRODUCTION Torque measurement on rotating shaft by the method of strain gauge measurements, especially with the existing equipment, brings usually two problems: - the first problem is the installation of the strain gauge to places giving representative data. These are usually areas without mounting, the groove for tongue area and similar discontinuities – simply such location on the shaft, where the relative deformation is in accordance with the torque value (under known formulas); - the second problem is the transmission of the measured signal from the rotating shaft to the recording apparatus. Our article deals with these issues. II. CONTACT TRANSMISSION The contact transmission by means of rings and brushes used to be the only option for the signal transmission in the past. There were usually the direct signals from strain gauges transmitted, that is why the apparatus had to be reliable and with only very low transition resistance values. Despite the use of high quality materials – especially precious metals – such devices were not only liable to failure but also unfit for measurements on the existing equipment. Contact transmission of signal is still used at present – the Hottinger sensor is one of the examples (see Fig. 1 and Fig. 2). The torque measurement by means of the existing equipment is very often required in practice. For this purpose one very simple, cheap and at the same time reliable signal transmission method was developed by professor Dejl and his colleagues – the signal transmission by means of rings and „wire“ brushes. There are rings with peripheral groove mounted on the insulation lining. The signal is then transmitted by a wire encircling the rings which is tighten by means of springs or coils.

Fig. 1: Hottinger torque sensor

Fig. 2: Hottinger torque sensor - detail of the signal transmission

This system was used successfully in many applications – for example for signal transmission of strain gauge data of the industrial automatic washer drum and shaft (Fig. 3 and Fig 4), or articulated cardan shafts of a rolling mill (Fig. 5 and Fig 6).

Fig. 3. Wired signal transmission on industrial automatic washer drum

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

78

gauge measurement on the automatic washer drum serves as an example (Fig. 7 and Fig. 8).

Fig. 7. Some strain gauges placed on the automatic washer drum Fig. 4. Wired and wireless transmission on shaft of an industrial automatic washer drum

Fig. 8. Detail of strain gauges Fig. 5. Torque measurements on the cardan shaft of rolling mill vertical cylinders

III. TELEMETRY WITH DISC AERIAL Modern electronics has brought a miniaturization and also possibilities of non-contact signal transmission telemetric systems. The telemetric system by ESA Messtechnik GmbH, Munich is an example. Its principle is described in Fig. 9, one of the practical realizations (experimental stand – Faculty of Mechanical Engineering, STU Bratislava) is in Fig. 10. The strain gauges are connected in the module with strain gauge amplifier, bridge charger and transmitting block. By means of a fixed aerial and a rotational aerial system both the power supply of strain gauge module and the signal transmission from the strain gauge bridge to the evaluation unit are secured.

Fig. 6. Torque measurements on the cardan shaft - detail

The application of this system depends on certain experience; it is necessary to set the angle of wire encircling, the tension of springs and also the dimensions and material for the wire and the groove in a appropriate ratio in accordance with the shaft diameter and its peripheral velocity. By means of this equipment it is possible to transmit not only the torque measuring strain gauge signal, but also signals from the strain gauge measuring any other quantities or from quite different sensors. The strain

Fig. 9: Telemetry block diagram

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

Fig. 10: Torque sensoring on the experimental stand

79

Above described system is relatively simple as for installation and usage; it is necessary to have space for placing of the rotational and fixed aerials. Maximum possible diameter of the rotational aerial depends on the power supply unit output. The output 1 W is enough for the diameter of 400 mm, the output 3 W is suitable for the diameter of 750 mm. Despite some problems with sufficient space this equipment was installed into the gearbox of fork-lift truck Jungheinrich Hamburg (Fig. 13) and performed measurements during test operation with the option of switching between strain gauges for the tensile-pressure strength and strain gauges for torque measurements on the conical gear pinion (Fig. 14).

The example of such telemetry usage in practice may be torque measurement of a conveyor belt drive. The space for the sensor locating was very limited and confined (see Fig. 11) and besides this fact it was necessary to place the strain gauges near the shaft mounting, which required additional specification of the scale between the measuring voltage and torque by means of FEM (Fig. 12).

Fig. 13: Location of aerials and strain gauges in the gearbox of fork-lift truck

Fig. 11: Torque measurement of a conveyor belt drive

Fig. 14: Location of the amplifier and connecting bar for strain gauges

Fig. 12: Calculating of conversion scale „measuring voltage – torque“ by means of FEM

This equipment can be used (after applying some suitable insulation measures) even in environment with running water – we have tested it during measurements on a running mill in Budapest (Fig. 15). The conditions of cylinders fixing unfortunately caused such shifts of the cardan shaft both in horizontal and vertical directions that the fixed aerial could not be installed properly – so the data obtained were only informative.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

Fig. 15: Telemetry on the cardan shaft of a rolling mill in Budapest

80

V. TELEMETRY WITH HIGH FREQUENCY TRANSMISSION This telemetry is used in more difficult conditions. The transmission is realized by means of high frequency of 433 MHz. The strain gauge (or other sensor) signal is filtered for values over 1 kHz (antialiasing) and digitized by the 12-bit converter. The digital signal is then transmitted at a distance of minimum 50 m (in free space), then decoded – and the final output is an analogue signal again. The block diagram of this telemetry is in Fig. 18, its realization in Fig. 19. This system is now being tested and will be used for measurements on an experimental car in Josef Božek Competence Centre for Automotive Industry.

IV. TELEMETRY WITH AERIAL RIGHT ON THE SHAFT If there is not enough space in the vicinity of the shaft, it is possible to install the rotational aerial right on the shaft. In our case this solution was used during measurements of torque on Škoda Fabia half-axles. The rotational aerial is installed right on the shaft (see Fig. 16 and Fig. 17)) and the fixed aerial is in this case closer to the shaft.

Fig. 16. Torque sensoring on a car half-axle Fig. 18: Block diagram of the high frequency telemetry

Fig. 17. Torque sensoring on a car half-axle

The used aerial incorporates also the revolution sensor. At each shaft revolution there is a short rectangular pulse on the output connector. Then it is possible to calculate the current revolutions from the pulse frequency – and in this case also the vehicle velocity (speed).

Fig. 19. High frequency telemetric system real appearance

The device can be also used for the signal scanning from moving objects. An example is the measurement of the go-karts frame loading, which realized a student of our university as part of his thesis.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

81

ACKNOWLEDGMENT This research has been realized using the support of Technological Agency, Czech Republic, programme Centres of Competence, project # TE01020020 Josef Božek Competence Centre for Automotive Industry. This support is gratefully acknowledged. REFERENCES

Fig. 20 –Telemetry on the go-car frame

Fig. 21 – Receiving point

CONCLUSION Despite this article looks rather as an advertising leaflet for ESA Messtechnik products, its aim is to inform colleagues in our branch. Especially those of them who need to measure the torque or other tensile strength of a rotating equipment – we would like to inform them of some options we have experience with. While the strain gauge system is based on the evaluation of the voltage on the Wheatstone bridge measurement diagonal, it is possible to connect any other sensor (which output is voltage) to the amplifier - and then transmit this signal to the evaluation unit. For example temperature, force, vibrations and other quantities can be transmitted in this way.

[1] Dejl, Z., Folta, Z., Zieschang, T. Zkrácené životnostní zkoušky kuželočelní převodovky vysokozdvižného vozíku Retrak. (Shortened life tests bevel gear forklift reach truck Retrak) .In Proceedings konference Trwalošč elementow i wezlow konstrukczyjnych maszyn górniczych. Gliwice-Ustoň: TEMAG, 2004. ISBN 83-917265-4-1. [2] Folta, Z., Polák, J., Nečas, J., Vrána, V. Ověření účinnosti pohonu pásového dopravníku. (Verification of the effectiveness of belt conveyor drive) Liberec: TU Liberec, 2005. International conference of Departments of Machine Parts and Mechanisms anthology 2005. [3] Dejl, Z., Folta, Z. Měření zátěžného spektra automobilové převodovky. (Measurement of the automotive gearbox loading spektra.) Bratislava: STU Bratislava, 2007. International conference of Departments of Machine Parts and Mechanisms anthology 2007. ISBN 978-80-227-2708-2. [4] Folta, Z., Hrudičková, M. Zpracování zátěžných spekter převodovky osobního automobilu. (Processing of loading spectra of a car gearbox.) Košice: TU v Košiciach, 2010. International conference of Departments of Machine Parts and Mechanisms 2010. ISBN 978-80-970-294-1-8. [5] Dejl, Z., Němček, M.: Research results in the field of automotive transmissions on dept. of Machine Parts and Mechanisms VŠB – Technical University of Ostrava. In: Proceedings of abstracts of Conference „Mechanical Engineering 2006, Bratislava, Slovak University of Technology, p. 46, ISBN 80-227-2513-7. [6] ESA Messtechnik GmbH Munich, Technical documentation.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

82

Performance of Steady-State Voltage Stability Analysis in MATLAB Environment Jan Veleba 1) 1)

University of West Bohemia in Pilsen, Regional Innovation Centre of Electrical Engineering, Pilsen, Czech Republic, e-mail: [email protected]

Abstract — In recent years, electric power systems have been often operated close to their working limits due to increased power consumptions and installation of renewable power sources. This situation poses a serious threat to stable network operation and control. Therefore, voltage stability is currently one of key topics worldwide for preventing related black-out scenarios. In this paper, modelling and simulations of steady-state stability problems in MATLAB environment are performed using authordeveloped computational tool implementing both conventional and more advanced numerical approaches. Their performance is compared with the Simulink-based library Power System Analysis Toolbox (PSAT) in terms of solution accuracy, CPU time and possible limitations. Keywords — Steady-state voltage stability, continuation load flow analysis, predictor-corrector method, voltage stability margin, voltage-power sensitivity, Power System Analysis Toolbox

I. INTRODUCTION Steady-state voltage stability is defined as the capability of the system to withstand a small disturbance (e.g. fault occurrence, small change in parameters, topology modification, etc.) without abandoning a stable operating point [1]-[5]. Voltage stability problems are generally bound with long "electrical" distances between reactive power sources and loads, low source voltages, severe changes in the system topology and low level of var compensation. However, this does not strictly mean that voltage instability is directly connected only with low voltage scenarios. Voltage collapse can arise even during normal operating conditions (e.g. for voltages above nominal values). Moreover, variety of practical situations can eventually lead to voltage collapse, e.g. tripping of a parallel connected line during the fault, reaching the var limit of a generator or a synchronous condenser, restoring low supply voltage in induction motors after the fault. All these cause the reduction of delivered reactive power for supporting bus voltages followed by increases of branch currents and further voltage drops to even lower reactive power flow or line tripping until the voltage collapse occurs. This entire process may occur in a rather large time frame from seconds to tens of minutes. To prevent voltage collapse scenarios, several types of compensation devices are massively used worldwide both shunt capacitors/inductors, series capacitors, SVCs, synchronous condensers, STATCOMs, etc. To reduce voltage profiles (in case of low demand), var consumptions must be increased by switching in shunt reactors, disconnecting cable lines (if possible), reducing voltage-independent MVAr output from generators and

synchronous condensers, etc. To increase bus voltages, opposite corrective actions are to be taken. These include reconfigurations (connecting parallel lines / transformers / cables), power transfer limitations and activations of new generating units at most critical network areas. Furthermore, the voltage load shedding of low-priority loads (usually by 5, 10 or 20 % in total) is usually realized at subtransmission substations using undervoltage relays. These relays work similarly as on-load tap-changing (OLTC) transformers. They are activated by long-term voltage dips (in region between 0.8 and 0.9 pu) and as the result, they trip the load feeders - typically in steps of 1 to 2 % of the load at any given time (with time delays of 1-2 minutes after the voltage dip). The larger voltage dip, the faster and larger response of the relay [2]. Low voltage profiles are usually averted by actions of OLTC transformers. However, each tap position corresponds to an increase of the load which eventually leads to higher branch losses and further voltage drops [1]-[2],[4]. Therefore, OLTC transformers should be blocked during low voltage stability scenarios. Negative effects of OLTC actions during low voltage conditions are presented in many studies with voltage stability margin calculation from synchrophasor measurements [6]-[7]. This paper is organized as follows. Chapters II and III describe conventional Cycled Newton-Raphson (N-R) and more robust Continuation Load Flow (CLF) methods for the voltage stability analysis, respectively. Independent tool - Power System Analysis Toolbox (PSAT) - is briefly introduced in Chapter IV. In Chapter V, key properties of both of the author-developed codes are discussed. Chapters VI and VII show the results of individual approaches when solving voltage stability of a broad variety of test power systems. Chapter VIII closes the paper with some concluding remarks and the evaluation of each technique applied. II. CONVENTIONAL NUMERICAL CALCULATION OF THE VOLTAGE STABILITY PROBLEM When increasing the loading (or loadability factor λ) of the system, its bus voltages slowly decrease due to the lack of reactive power. At the critical point (called singular or bifurcation), characterized by maximum loadability factor λmax and critical bus voltages, the system starts to be unstable and voltage collapse appears (system black-out). From this point on, only lower loading with low voltage values leads to the solution. The dependence between bus voltage magnitudes and λ is graphically represented by the V-P curve (also referred to as the nose curve). Unfortunately, the current (so-called base-case) position of the system operating point on the V-P curve is

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

not known along with its distance from the voltage collapse (so-called voltage stability margin). Thus, location of the singular point must be found during the voltage stability analysis. Note: Values of λmax and critical voltages are rather theoretical since they do not reflect voltage/flow limits of network buses/branches. When incorporating these practical restrictions, the real maximum loadability λmax* can be found (i.e. the maximum value for keeping all network buses and branch loadings within limits). Traditional approach for finding the maximum system loadability is to apply the standard N-R method [8] for the base-case load flow calculation (i.e. for λ = 1.0). When obtaining current position on the V-P curve, network loading (i.e. loads/generations in selected network buses) is increased in defined manner by a certain step and the load flow is computed repetitively along with a new position on the V-P curve. This process continues in an infinite loop until the singular point is reached. However, total number of iterations in each V-P step is gradually increasing so that when close to the singular point, the NR method fails to converge, i.e. no solution is provided. This relates to the fact, that Jacobian J becomes singular (i.e. det J ≈ 0) and its inverse matrix cannot be computed for successful numerical convergence. For speeding up the calculation, a variable step change is applied. Usually, a single default step value is used. When obtaining the divergence of the N-R method, the step size is simply divided by two and the calculation for the current V-P point is repeated until the convergence is achieved. When the current step size value reaches the pre-set minimum value, the calculation is stopped. Despite of the relatively simple procedure, the Cycled N-R method enables the completion of the stable V-P part only. The unstable part including the singular point cannot be examined. Also, high CPU requirements prevent this method from being employed for larger power systems. In this paper, the Cycled N-R algorithm was developed and further tested on wide range of test power systems.

83

for the selected CP in each network bus k at the current point on the V-P curve. Remaining elements in (1) are the newly computed Jacobian J and step size σ of the CP. The tangent predictor is relatively slow, anyway shows good behaviour especially in steep parts of the V-P curve. Unlike the tangent predictor, secant predictor is simpler, computationally faster and behaves well in flat parts of the V-P curve. In steep parts (i.e. close to the singular point and at sharp corners when a generator exceeds its var limit) it computes new predictions too far from the exact solution. This may eventually lead to serious convergence problems in the next corrector step. Therefore, the tangent predictor is more recommended to be applied. The corrector is a standard N-R algorithm for correcting state variables from the predictor step to satisfy load flow equations. Due to one extra parameter λ, additional condition (2) must be included for keeping the value of the CP constant in the current corrector step. This condition makes the final set of equations non-singular even at the bifurcation point.

λ x k − x kpredicted = 0 , x =  V

if CP is λ (2) if CP is V

Difference between both types of predictors and the entire process of the predictor/corrector algorithm is graphically demonstrated in Fig. 1. Horizontal/vertical corrections are performed with respect to the chosen CP type.

III. CONTINUATION LOAD FLOW ANALYSIS The CLF analysis [1],[9] suitably modifies conventional load flow equations to become stable also in the bifurcation point and therefore being capable of drawing both upper/lower parts of the V-P curve. It uses a two-step predictor/corrector algorithm along with the new unknown state variable called continuation parameter (CP). The predictor (1) is a tangent extrapolation of the current operation point estimating approximate position of the new point on the V-P curve. −1

θ  V     λ 

predicted

M   0   J M K    θ 0   M    (1) = V0  + σ  M   0   λ0  L L L L   1     ek

Vector K contains base-case power generations and loads. Variables θ0, V0, λ0 define the system state from the previous corrector step. The vector ek is filled with zeros and certain modifications (see [1],[9]) are implemented

Fig. 1. Predictor/Corrector Mechanism for the CLF Analysis [10].

As the CP, state variable with the highest rate of change must be chosen (i.e. λ and V in flat and steep parts of the V-P curve, respectively). When the process starts diverging, parameter σ must be halved or parameter CP must be switched from λ to V. The step size should be carefully increased to speed-up the calculation when far from the singular point or decreased to avoid convergence problems when close to the peak. The step size modification based on the current

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

position on the V-P curve (i.e. as a function of the line slope for previous two corrected points on the V-P curve) is recommended in [11]. This approach belongs to socalled rule-based or adaptive step size control algorithms. In [10], several voltage stability margin indices (VSMIi, VSMIik) are presented along with relative var reserve coefficient and voltage-load sensitivity factors (VSFi) for comprehensive voltage stability analysis and location of weak or sensitive system buses/branches/areas. In these regions, preventive or remedial actions should be taken. Procedures for allocating individual compensation devices and possible effects are also discussed. The CLF analysis still remains very popular for highspeed solving of voltage stability studies. Due to its reliable numerical behaviour, it is often included into the N-R method providing stable solutions even for illconditioned load flow cases. Moreover, it is applied in foreign control centres for N-1 on-/off-line contingency studies with frequencies of 5 and 60 minutes [2], respectively. IV. POWER SYSTEM ANALYSIS TOOLBOX (PSAT) PSAT [12] is a Simulink-based open-source library for electric power system analyses and simulations. It is distributed via the General Public License (GPL), its download and use is free of charge. However, there is no warranty that the Toolbox will provide correct and accurate results. All corrections and possible repairs or improvements are to be done on the customer side. It contains the tools for Power Flow (busbars, lines, two-/three-winding transformers, slack bus(es), shunt admittances, etc.), CLF and OPF data (power supply/demand bids and limits, generator power reserves and ramping data), Small Signal Stability Analysis and Time Domain Simulations. Moreover, line faults and breakers, various load types, machines, controls, OLTC transformers, FACTS and other can be also modelled. User defined device models can be added as well. All studies must be formulated for one-line network diagram only - either in input data *.m file in required format or in graphical *.mdl file, where the schema is manually drawn. For the former option, input data conversions from and to various common formats (PSS/E, DIgSILENT, IEEE cdf, NEPLAN, PowerWorld and more others) are available. When compared to another MATLAB-based opensource tool MATPOWER [13], PSAT is more efficient and highly advanced by providing more analyses, problem variations, possible outputs and other useful features in its user-friendly graphical interface. MATPOWER does not support most of advanced network devices, entirely omits CLF analysis and has no graphical user interface or graphical network construction ability. Also, it does not consider var limits in PV buses. Incorrect interpretation of reactive power branch losses can be also observed. V. PROPERTIES OF AUTHOR-DEVELOPED CODES IN MATLAB ENVIRONMENT Both Cycled N-R and CLF procedures were developed in MATLAB environment for providing fundamental examination of medium-sized and larger power systems in terms of steady-state voltage stability. Several key aspects of these codes are discussed below.

84

1] Predictor: Despite of computationally more complex algorithm, the tangent predictor was used for finding reliable estimations of new V-P points especially around the singular point. It is applied in CLF algorithm only. 2] Corrector: First, a corrector step is used at the start of the CLF program to find the base-case point for further calculations. Due to possible weak numerical stability at this point (for badly-scaled power systems), the One-Shot Fast-Decoupled (OSFD) procedure is implemented to the standard N-R method for providing more stable solutions and thus preventing numerical divergence. Moreover, voltage truncation (SUT algorithm) is also included into the state update process at every N-R's iteration. Both of these stability approaches were introduced in [14] and further tuned and tested in [15]. Both were also applied to the Cycled N-R algorithm to increase the loading range, for which the stable load flow solutions can be obtained (i.e. closer proximity to the singular point can be reached). 3] Step size: Largest-load PQ network bus is chosen for computing the angle α between the horizontal and the line interconnecting two adjacent V-P points. Based on this, the step size evaluation function (3) is applied - see Fig. 2.

σL   σ= σU  A / sin 2 α + B 

for

α ≥ π/ 8

for

α ≤ π/ 32 otherwise

(3)

The upper and lower step limit constants σU and σL define the step size for the flat part of the V-P curve and for close vicinity to the singular point, respectively.

Fig. 2. Step size evaluation function [10].

For the Cycled N-R algorithm, this is a rather too complex concept of the step size control. Therefore, only a single step size is chosen at the start and a simple stepcutting technique (dividing by 2) is applied in case of divergence. 4] Ending criterion: Only stable part of the V-P curve (incl. exact singular point calculation) is computed by the CLF code. Thus, if the computed value of λ begins to decrease, the process is stopped. For the Cycled N-R code, the calculation is terminated when the step size falls below a certain small value (e.g. 1×10-8). For each load flow case, maximum number of iterations and permitted tolerance for convergence is set to 20 and 1×10-8, respectively.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

5] Calculation speed and accuracy: For excessively accurate voltage stability solutions, values of σU and σL 2.5×10-2 and 6.25×10-4 are used in the CLF algorithm. Rather compromise values of 5×10-2 and 1×10-2 are also used to obtain fast and fairly accurate solutions for any of tested power systems. For the Cycled N-R algorithm, initial step size of 2.5×10-2 is chosen as sufficient. 6] Code versatility: Both Cycled N-R and CLF procedures are programmed so that the user directly specifies an arbitrary group of network buses for an load/generation increase. From this set of buses, only those non-slack buses with non-zero active power loads/generations are activated for the analysis. In all studies performed, a load/generation increase in the entire network (i.e. all network buses selected) is considered. Two scenarios can be activated by the user. a) L scenario increases both P/Q loads in selected PQ/PV buses with a constant power factor (i.e. with identical increase rate). b) L+G scenario increases both P/Q loads in selected PQ/PV buses and P generations in selected PV buses (with identical increase rate). 7] Var limits: In both approaches, bus-type switching logics are applied to iteratively computed reactive powers QGi in PV buses when exceeding the var limit (4) or to relevant bus voltages when returning the vars back inside the permitted var region (5). Symbol (p) denotes the current iteration number.

if QGi( p ) > QGi max Q QGi( p ) =  Gi max ( p) QGi min if QGi < QGi min

Vi = Vi sp

(4)

QGi( p ) = QGi max AND Vi > Vi sp    if  OR  (5) ( p ) sp Q = Q  Gi min AND Vi < Vi   Gi

The terms QGimax and QGimin are the upper and lower var limits, the term Visp determines the specified value of the voltage magnitude for each PV bus. 8] Code limitations: a) With increased loading, lower/upper var limits in PV buses should not be fixed but variable proportionally to the generated active power. In both codes, constant var limits are used for more pessimistic V/Q control. b) Only identical increase rate is applied. However, implementing user-defined increase rates for each load/generation would not pose any serious problem. 9] Sparse programming: Sparsity techniques along with smart vector/matrix programming are used in both Cycled N-R and CLF codes to significantly decrease the CPU time needed for each load flow case. 10] Outputs: Theoretical value of λmax and V-λ data outputs for V-P curves are computed and stored or graphically projected. Respective values of λ for switching some of PV buses permanently to PQ are also recorded. Voltage and power flow limits were not considered for the evaluation of the real maximum loadability λmax*.

85

VI. TESTING OF CYCLED N-R AND CLF ALGORITHMS FOR SOLVING VOLTAGE STABILITY LOAD FLOW PROBLEMS Total number of 50 test power systems between 3 and 734 buses were analyzed using developed Cycled N-R and CLF algorithms in the MATLAB environment. Identical increase rate was applied to all network buses (before filtering those with non-zero active power loads or generations). For both L and L+G scenarios, only stable part of the V-P curve was calculated with included var limits. Settings of both codes are as introduced in Chapter V, paragraphs 4] and 5]. In Tab. I., voltage stability solutions of several test cases are shown. Presented results contain the maximum loadability, numbers of stable V-P points and CPU times in seconds needed. For each case, the first two rows show the outputs of the CLF code for excessive accuracy and compromise accuracy, respectively. For comparison purposes, the third row provides the results of the Cycled N-R code. TABLE I. VOLTAGE STABILITY SOLUTIONS USING CYCLED N-R AND CLF ALGORITHMS - L AND L+G SCENARIOS Scenario L λmax [-] pts CPU [s] 1.302632 331 0.5616 IEEE9II 1.302632 27 0.1404 1.302632 23 0.4056 1.760331 658 1.2012 IEEE14 1.760331 87 0.2340 1.760331 43 0.5460 1.536905 854 1.9500 IEEE30 1.536905 88 0.2808 1.536905 37 0.6396 1.406778 891 2.9016 IEEE57 1.406778 229 0.6864 1.406778 27 0.8112 1.079959 1640 12.9169 IEEE162 1.079960 464 3.1044 1.079960 13 1.7628 1.024573 8457 103.8655 IEEE300 1.024573 529 7.0044 1.024573 16 2.4180 3.104162 139 4.5864 EPS734II 3.104083 46 1.8720 3.104162 96 8.2369 Case

Scenario L+G λmax [-] pts CPU [s] 1.162053 215 0.3900 1.162052 24 0.1248 1.162053 20 0.4212 1.777995 506 0.9360 1.777995 59 0.2028 1.777995 45 0.6396 1.546751 726 1.6536 1.546752 124 0.4212 1.546751 37 0.6552 1.616845 399 1.3884 1.616845 57 0.2652 1.616845 37 0.8112 1.138996 1185 9.3913 1.138996 65 0.8112 1.138996 16 1.8408 1.058820 311 4.0092 1.058819 94 1.4508 1.058820 17 2.5584 3.104162 139 4.8360 3.104083 46 1.8408 3.104162 96 8.1745

As can be seen, exact solutions of maximum loadability were obtained for both of tested methods and each of the three accuracy settings. The first setting was definitely too much focused on producing exact results. Therefore, numbers of V-P points and CPU times were pushed often above 200 and 1 second, respectively. When using fair compromise setting, the maximum error for λmax from all 50 test power systems was only 0.0185 percent, while numbers of points and CPU times were decreased on average by 75.27 percent and 64.11 percent, respectively. The Cycled N-R code obtains highly accurate results in terms of solution accuracy. In majority of cases, it provides even better solutions than CLF algorithm with compromise accuracy. Surprisingly, it always computes slightly higher maximum loadability values than by the high-accurate CLF code. This seems to be one visible drawback of the Cycled N-R method. Only low numbers of V-P points are needed for reaching close proximity to the singular point. These numbers are well comparable to those needed for the compromise CLF code.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

Unfortunately, each divergence case (between 22 and 28) significantly delays the entire computation process of the Cycled N-R method. Therefore, the Cycled N-R code suffers from being extremely time-dependent on computing of each V-P point. When compared with the compromise CLF code, the CPU time needed by the Cycled N-R method is on average about 167 % higher. Therefore, the compromise CLF code seems to be the best method for providing fast and highly accurate voltage stability solutions. Stable V-P curves of the IEEE 30-bus power system (L+G scenario) are computed by both Cycled N-R and CLF methods and shown in Figs. 3 and 4, respectively. For the CLF method, the V-P curves are extended to demonstrate numerical stability of the CLF algorithm around the singular point. Extension of V-P curves in the unstable region is provided for 0.97×λmax < λ < λmax.

Voltage Magnitude [pu]

1

0.9

0.8

0.7

0.6

0.5

0.4

1

1.1

1.2

1.3

1.4

Loadability Factor λ [-]

1.5

Fig. 3. V-P curves for the IEEE 30-bus system (Cycled N-R method).

86

VII. TESTING OF PSAT FOR SOLVING VOLTAGE STABILITY LOAD FLOW PROBLEMS Before solving voltage stability in PSAT, the load flow analysis of a system must be performed. Therefore, large number of load flow studies is solved using PSAT to detect any of its possible weaknesses. Results were compared with the author-developed N-R code in MATLAB. Despite of unconstrained network size to be solved, several limitations of PSAT were found during the testing stage. 1] Inefficient PV-PQ bus type switching logic is applied. Probably, reverse switching logic (5) is not used and the need for convergence is requested to activate forward switching logic (4). As a result, unnecessarily more PV buses are being switched permanently to PQ. Furthermore, the switching logic completely fails to switch PV buses to PQ for larger systems with high numbers of PV buses. 2] Nominal voltages must be defined in the input data file or the error message 'Divergence - Singular Jacobian' is obtained during the simulation. This seems to be entirely illogical since nominal voltages should not be necessary for the 'in per units defined' problem. 3] It seems that no advanced stability techniques are applied for the N-R method in PSAT because of severe numerical oscillations appearing in several studies. 4] PSAT intentionally neglects transformer susceptances and thus causes errors in final load flow results. A column for shunt susceptances is available for power lines only. For transformers, this column is filled with zeros by default. Under these limitations, load flow results show very good congruity between the author-developed N-R method and PSAT. Higher total numbers of iterations are needed by PSAT due to missing stability technique(s). Also, CPU times are higher in PSAT due to combining the codes with other analyses and related tool features. As an example, the load flow and voltage stability analysis of the IEEE 14-bus system is done by PSAT (Figs. 5-9).

Voltage Magnitude [pu]

1

0.9

0.8

0.7

0.6

0.5

0.4

1

1.1

1.2

1.3

1.4

Loadability Factor λ [-]

1.5

Fig. 4. Extended V-P curves for the IEEE 30-bus system (CLF method).

As Tab. I. indicates, applied version of the CLF method is still not applicable for real-time voltage stability monitoring, but it can be useful for off-line reliability, evaluation or planning studies of even larger networks.

Fig. 5. GUI in PSAT for the load flow analysis of IEEE 14-bus system.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

87

Fig. 8. Settings for the CLF analysis when solving the IEEE 14-bus power system.

Fig. 6. Final voltage magnitudes of IEEE 14-bus power system in PSAT.

The CLF algorithm in PSAT is defined so, that power increases are realized by adding a power increment (loadability factor multiplied by the increase rate) to the base-case loading, i.e. initial λ is zero. In the authordeveloped Cycled N-R and CLF codes, power increases are performed by multiplying the base-case loading with λ. Therefore, the maximum loadability in PSAT must be increased by unity when comparing both codes.

Fig. 7. Final voltage angles of IEEE 14-bus power system in PSAT.

For voltage stability studies, PSAT contains the advanced CLF algorithm, which is combined with contingency and OPF analyses. Load flow data are extended by two matrices specifying the sets of PQ and PV buses, where the loads and generations are to be increased (different increase rates are possible). Thus, various loading scenarios of the system can be modelled. The CLF code is then started via a specialized window (Fig. 8). Calculation can be adjusted by the user for better computational performance - e.g. by setting a more suitable step size, maximum number of V-P points or by checking the option for controlling voltage, flow or var limits. PSAT offers two CLF methods - perpendicular intersection (PI, as in Fig. 1) and local parametrization (LP). Three stopping criteria are available: The complete Nose Curve (computing both stable/unstable parts of the V-P curve), Stop at Bifurcation (when singular point exceeded) and Stop at Limit (when voltage/flow/point limit hit).

Fig. 9. Nose curves for all network buses of the IEEE 14-bus test system in PSAT.

In Tab. II., voltage stability results for medium-sized IEEE test systems are provided by the author-developed Cycled N-R and compromise CLF codes. These outputs are compared to those obtained by PSAT - see Tab. III. In PSAT, both of the CLF modes were tested (i.e. PI with step 0.025 and LP with default step 0.5). TABLE II. VOLTAGE STABILITY ANALYSIS OF MEDIUM-SIZED IEEE TEST SYSTEMS (CYCLED N-R VS. COMPROMISE CLF) L+G

Cycled N-R code

Compromise CLF code

Case

λmax [-]

pts

CPU [s]

λmax [-]

pts

CPU [s]

IEEE9

2.485393

74

0.5460

2.485382

84

0.2964

IEEE13

4.400579

148

0.6708

4.400577 112

0.3120

IEEE14

4.060253

137

0.8268

4.060252

92

0.3276

IEEE24

2.279398

61

0.6396

2.279398

58

0.2496

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

IEEE30

2.958815

88

0.8112

2.958814

57

0.2964 0.3432

IEEE35

2.888962

91

0.6864

2.888950 107

IEEE39

1.999203

57

0.7644

1.999202

30

0.2184

IEEE57

1.892091

47

0.8892

1.892089

92

0.4836

IEEE118

3.187128

100

1.6536

3.187128

66

0.5772

TABLE III. VOLTAGE STABILITY ANALYSIS OF MEDIUM-SIZED IEEE TEST SYSTEMS (PSAT - PI VS LP MODE) L+G

PSAT - PI mode

88

ACKNOWLEDGMENT This paper has been supported by the European Regional Development Fund and the Ministry of Education, Youth and Sports of the Czech Republic under the Regional Innovation Centre for Electrical Engineering (RICE), project No. CZ.1.05/2.1.00/03.0094. This work has been also sponsored by Technology Agency of the Czech Republic (TACR), project No. TA01020865 and by student science project SGS-2012-047.

PSAT - LP mode

Case

λmax [-]

pts

CPU [s]

λmax [-]

IEEE9

2.481220

7

0.2093

2.482000

13

0.3241

IEEE13

4.390420

13

0.3292

4.399570

20

0.4832

IEEE14

4.060100

18

0.4098

4.059420

19

0.4939

IEEE24

2.277550

10

0.2600

2.278670

16

0.4313

IEEE30

2.958550

16

0.8761

2.958250

20

1.5023

IEEE35

2.872940

16

1.1242

2.878420

10

0.2940

IEEE39

1.999110

11

0.2932

1.997840

12

0.3692

pts CPU [s]

IEEE57

1.891920

12

0.9089

1.892090

26

3.9090

IEEE118

3.187100

613 19.1693

3.187120

82

19.7464

Theoretical values of λmax, numbers of stable V-P points and CPU times in seconds are provided for comparison. For all voltage stability studies in PSAT, identical power increase rates were considered. Only the L+G scenario was examined, logics for var limits were deactivated. Both of PSAT modes showed only average accuracy with satisfiable numbers of V-P points and rather lower computational speed. The LP mode was computationally more time-consuming, but needed lower numbers of V-P points and usually provided more accurate results. The compromise CLF code provided the best combination of solution accuracy and CPU time requirements in each of the cases. Although higher numbers of V-P points were needed, CPU times were still significantly smaller than those in PSAT due to optimized sparse programming applied. Identical conclusions can be made when mutually comparing CLF and Cycled N-R codes. VIII. CONCLUSION For solving voltage stability problems, both the Cycled N-R and CLF codes were programmed and comprehensively tested on a broad range of test power systems in MATLAB environment. Various stability techniques, step size approaches and numerical settings were applied and used to upgrade their performance in order to find the algorithm with fair compromise between calculation speed and solution accuracy. The results were compared with outputs obtained from PSAT. The studies imply that the best technique (i.e. best combination of precision level and CPU requirements) is the CLF algorithm with compromise step size settings programmed by Author in MATLAB. However, final technique can be applicable in practice only for off-line planning and development studies of electric power systems. For realtime evaluations of system's voltage stability, a more robust algorithm with minimized numbers of stable V-P points must be developed. Therefore, follow-up research activities will be concentrated especially on this area of interest.

REFERENCES [1] P. Kundur, Power System Stability and Control, McGraw-Hill, 1994. [2] C. Canizares, A.J. Conejo and A.G. Exposito, Electric Energy Systems: Analysis and Operation, CRC Press, 2008. [3] V. Ajjarapu, Computational Techniques for Voltage Stability Assessment and Control, Springer, 2006. [4] I. Dobson, T.V. Cutsem, C. Vournas, C.L. DeMarco, M. Venkatasubramanian, T. Overbye and C.A. Canizares, “Voltage Stability Assessment: Concepts, Practices and Tools - Chapter 2,” Power System Stability Subcomittee Special Publication, IEEE Power Engineering Society, 2002. [5] J.H. Chow, F.F. Wu and J.A. Momoh, Applied Mathematics for Restructured Electric Power Systems - Optimization, Control and Computational Intelligence, Springer, 2005. [6] I. Šmon, M. Pantoš and F. Gubina, “An improved voltage-collapse protection algorithm based on local phasors,” Electric Power Systems Research 78 (2008), pp. 434-440, 2008. [7] Y. Gong and N. Schulz, “Synchrophasor-Based Real-Time Voltage Stability Index,” Proceedings of PSCE conference, pp. 1029-1036, 2006. [8] J.J. Grainger and W.D. Stevenson, Power System Analysis, McGraw-Hill, 1994. [9] M. Crow, Computational Methods for Electric Power Systems, CRC Press, 2002. [10] J. Veleba, “Application of Continuation Load Flow Analysis for Voltage Collapse Prevention,” Journal Acta Technica, vol. 57, pp. 143-163, 2012. [11] P. Zhu, “Performance Investigation of Voltage Stability Analysis Methods,” PhD. thesis, Brunel University of West London, 2008. [12] Homepage of PSAT. [Online]. Available: http://www3.uclm.es/profesorado/federico.milano/psat.htm. [Accessed: 3 Mar. 2013]. [13] Homepage of MATPOWER. [Online]. Available: www.pserc.cornell.edu/ matpower/. [Accessed: 23 Dec. 2012]. [14] Z. Tate, “Initialization Schemes for Newton-Raphson Power Flow Solvers,” [Online]. Available: http://grb.physics.unlv.edu/~zbb/files/upload/29UV3GPCVQWO6 ERVE9RABFQ5M.pdf. [Accessed: 15 Apr. 2010]. [15] J. Veleba, “Acceleration and Stability Techniques for Conventional Numerical Methods in Load Flow Analysis,” Proceedings of ELEN conference, pp. 1-10, 2010.

THE AUTHORS Jan Veleba received his Master degree in Power Engineering at the University of West Bohemia, Pilsen, Czech Republic in 2008. Currently, he is the PhD. student at the Department of Electrical Power Engineering and Environmental Engineering at the University of West Bohemia in Pilsen, Faculty of Electrical Engineering. His main research activities concern load flow, voltage-power control and voltage stability analyses of particularly larger and more complex electric power systems.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

89

Cross slide mathematical model for solving chatter Jiří Ondrášek VÚTS, a.s./Mechatronics Department, Svárovská 619, 460 01 Liberec XI, Czech Republic, e-mail: [email protected]

Abstract—The paper deals with the issues of creating a mathematical model of self-excited chatter of the cross slide. This type of oscillation occurs in systems where there is still an internal source from which the system consumes energy to maintain or even to increase the amplitude of oscillation. This consumption is controlled by an oscillatory movement of the system itself. Such an energy source greatly affects the dynamic and stability properties of the system.

objects, in the following part of the text there are again provided some basic relationships of continuous cutting machining operation that are based on those assumptions: 1) Cutting force F is proportional to width b and depth of cut yN in the direction of the normal to the machined surface by the relation:

Keywords—chatter, cross slide, chip thickness, lobe diagram

2) Cutting force constant C0 is independent of speed. 3) Cutting edge geometry does not affect the direct correlation between force F and chip depth yN in the normal direction. 4) Angle β (an angle between cutting force F and the axis perpendicular to the normal of the machined surface) does not change with chip depth yN. 5) Friction forces between tool, workpiece and leaving chip are neglected. In Fig. 1, there is shown a block diagram of the continuous machining process in which the dynamic system can be described by the so-called oriented dynamic compliance Gy(s). This is a transfer function between the Laplace transforms of cutting force F(s) and the movement of the tool group y(s):

I. INTRODUCTION One of the main causes of generating self-excited vibration in mechanical systems is dry friction between two mutually moving parts that are directly related to damping in the system. Chatter is undesirable and there are efforts to avoid it either by increasing positive damping or eliminating the causes of negative damping. The frequency of steady chatter is close to the natural frequency of a mechanical system. This phenomenon often occurs during machining operation when a part of the energy of cutting process during cutting operation can change in the energy oscillating the machine as a whole. Vibration then manifests in a significant waviness of the cut surface and is usually accompanied by noise. Generally, it is theoretically possible to establish a certain range of cutting conditions under which, when applying them, no chatter arises. One of the means of such a designation is speed stability diagram – lobe diagram, which expresses the dependence of chip thickness on the workpiece speed. The methodology of creating lobe diagrams is based on the Laplace transform images of the cutting force and movement of a tool group in the direction of material removal. The very issue of dealing with chatter in cutting machining operation is creating the mathematical models of the following physical objects: Model of the machining process, Model of mechanical system inclusive flexible links and real constraints, Drive model that presents a model of electromotor itself and its control. II. CHATTER A mathematical process of machining process is described in detail in publication [2], see Chapter 5, a brief description of continuous machining is given in article [1], as the case may be. In the case of creating a mathematical model of the cross slide to investigate chatter it is assumed the latter mentioned method of machining. For this reason and for the reason of the logical linking of mathematical models of particular

F = −C0byN .

G y (s) =

y(s) , F(s)

(1)

(2)

where s is the complex variable. If depth of chip y0(t) < 0 is specified, cutting force F is generated which will cause the movement of tool group y(t) that will be superimposed to the specified depth, so instantaneous chip depth yN(t) is given by an expression:

yN (t ) = y0 (t ) + y (t ).

(3)

The chip thickness is negative because chip cutting occurs in the opposite direction of the y-axis orientation. If there is a case of yN(t) > 0, the cutting edge got out from engagement. When working with a fixed time link between the cuts, the freshly machined surface will get again into contact with the tool after a defined time with the socalled transport time delay Td. This fact can be expressed with a relation:

y N (t ) = y0 (t ) + y (t ) − y (t − Td ),

(4)

where y0(t) is the feed per one revolution, term y(t) is valid for immediate waviness and term y(t-Td) applies for waviness from the previous cut which will come under the cutting edge of the tool over a time period Td.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

90

y y(t-Td)

-Tds

F

+y(t1-Td)



y0 yN(t1)

β

y0(t1) R

y(t)

–y(t1)

x

e

+

+

yN

−C 0 b

F

G y (s)

y

+

vobv = 2πRnZ Fig. 1. Continuous machining and its block diagram.

3) From pairs f ↔ nZ and f ↔ bmez, there are assembled pairs nZ ↔ bmez which are shown as a set of curves – lobs for various values N again. 4) The individual „lobs“ may intersect each other and the area of the stable widths is under their lower envelope.

Waviness y(t-Td) from the previous cut together with an immediate tool position y(t) always affects the variation of chip depth yN(t), see Fig. 1, where the instantaneous chip depth yN at time t = t1 is shown. The Laplace image of equation (4) is:

(

y N (s) = y0 (s) + y(s) 1 - e

-Td s

).

(5)

III. EFFECT OF FEED DRIVES ON THE CHATTER

Using equations (1), (2) and (5), the closed loop transfer as to Fig. 1 can be determined in the following form: G Celk (s) =

−C 0 bG y (s) y(s) = y 0 (s) 1 + C b G (s) 1 - e- Td s 0 y

(

)

.

The shape of speed diagram and defining the areas of stability in the diagram depend, in addition to the properties of the coupled mechanical system, also on the characteristics of feed drives. In the same way as the behavior and properties of a mechanical system, also the behavior and properties of a control drive can be considered. They are commonly expressed by the dynamic flexibility of the control when it is again a transfer function. The total dynamic compliance of the associated dynamic system is then determined on the basis of the causal interconnection of particular dynamic systems, i.e., feed drives and a mechanical system. This interconnection is simple since the outputs of the one system are inputs of the second system, as shown in Fig. 2 and Fig. 3. In computer simulations it is assumed that the drive of the feeds of machine tool will be implemented by a 3phase synchronous electromotor with permanent magnets with which an exciting magnetic rotor flow is produced. For the purposes of simulations, the model of this servomotor was replaced by a simplified model which is based on the description of a DC motor, see Fig. 2, in which the block Mechanical system is a motor rotor with mass moment of inertia J. The parameters of a single coil, ie., inductance Ls, electrical resistance Rs and motor voltage constant KE are substituted in this model. Only at torque constant KT it is necessary to substitute the value:

(6)

To ensure the stability of the machining process it is necessary so that chip width b does not exceed the limit chip width bmez at which the amplitude characteristic of the closed loop transfer (6) shall not exceed the unit level at any angular frequency ω. This condition is expressed by the relationship:

b < bmez = −

1 = funkce1(ω ) , (7) 2C0 Reneg Gy ( jω ) min

(

)

in which RenegGy represents negative values of the real component of dynamic compliance Gy, see [2]. The stability lobe diagram expresses the dependence of limit chip width bmez on the speed of workpiece nZ. The speed equation, see [1], is given by the expression:  Re G y ( jω )  1  = funkce 2(ω ) , (8) f = nz  N + 1 − arctg  π Im G y ( jω )  

where f = ω/2π is the frequency of self-excited vibrations on the limit of stability for each value of N individually. Number N indicates the number of whole waves of the workpiece surface ripple to be incurred over period Td = 1/nZ. The full version of the stability lobe diagram is generated in practice as follows: 1) To the theoretically or experimentally determined set of the values of dynamic compliance, there are completed arranged pairs f ↔ bmez as to (7). 2) By means of (8), there are created pairs f ↔ nZ repeatedly for values N = 0,1,2,… which may be further displayed by a set of curves f = f(nZ). U

+ –

1 Lss + R s

UE

I

K TCelk =

3 KT , 2

(9)

whereby the common interaction of all three coils is taken into account, see [4]. Furthermore, in the scheme, there stands U for voltage, I electrical current, M electromagnetic torque, ω rotor angular velocity and φ rotor angular displacement. The above given parameters of the electromotor were substituted in its simulation model from the data sheet stated by the manufacturer.

KTCelk

M

Mechanical system

KE

Fig. 2. A simplified block diagram of a synchronous motor.

φ ω

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

+ y*

Ky



y

v*

+

91

Kv



Motion function

Tiv s + 1 I* Tiv s

+ –

v Mechanical ω system

KI

TiIs + 1 TiIs

I

U Electromotor

M

Fig. 3. A cascade control circuit with a current speed and position feedback.

In the vector control of this type of electric motor, it is almost exclusively used the cascade control circuit with three hierarchically arranged feedbacks: current, speed and position, see Fig. 3. Maintaining the required values of position, revolutions and current is ensured by PID linear controllers. Constant K is a proportional component of the controller, constant Ti expresses the integral time constant of the controller and Td is a derivative time constant, υ = y, v, I. Index υ reflects competence of the various constants to the position, velocity and current feedback. Those constants in the control structure according to Fig. 3 were debugged for the needs of the simulation model by the Ziegler-Nichols method, see [5]. In the diagram in Fig. 3, y expresses the position of the tool group, v its velocity, and * denotes the desired value of a particular variable. IV. MATHEMATICAL MODEL OF A MECHANICAL SYSTEM A. Mathematical description of a coupled mech. system To build equations of motion of a coupled mechanical system with flexible links, it is usually based on the Lagrange equations of mixed type which in matrix form are as follows:

d  ∂Ek  dt  ∂q&

  ∂Ek  −    ∂q

  ∂E p   ∂Rd  +   +    ∂q   ∂q&

 ∂f V   = Q +    ∂q

  λ , (10) 

where Ek and Ep represent the kinetic and potential energy of a mechanical system, Rd the so-called Rayleigh dissipative function and vector Q represents the vector of action generalized forces whose components correspond to appropriate coordinates qj. To describe the coupled mechanical system, there are used generally dependent physical coordinates q of dimension r, which are coupled by a system s of scalar constraint conditions:

Knowing kinetic Ek, potential Ep, and dissipative Rd energies of the coupled mechanical system, it is possible to establish equations of motion. Together with the second time derivative of constraint conditions fV, a summary record of these equations in matrix form can be included: M  J

T &&  p1 − Dq& − Kq  − J  q   =  . p2 0  λ   

(12)

in which:

J=

∂f V T ∂q

(13)

expresses the Jacobi matrix of the system of coupling equations. Symbols M, K and D are gradually mass, stiffness and damping matrices of the mechanical system as a whole. Vectors p1 and p2 are: T

& q& + 1  ∂M q&  p1 = Q − f g − M 2  ∂q  ∂ ∂  ∂f V p 2 = − T (Jq& )q& − 2 T  ∂q ∂q  ∂t

q& ,  ∂ 2f V q& − 2 ,  ∂t 

(14)

where vector fg is the vector of gravitational forces. Equations (12) compose a system (r+s) of algebraicdifferential equations for r unknown generally dependent physical equations q and s unknown Lagrange multipliers λ. These equations are currently assembled and solved by computational mechanics. Comprehensive information about creating mathematical models of multibody systems with flexible bodies is cited e.g. in publication [3].

For i-th flexible (pliable) body, coordinate vector qi can be written as:

B. Frequency response of the coupled mechanical system Transfer functions between the acting forces and the mass displacements of the controlled mechanical system with discretely distributed mass and stiffness parameters are generally determined by the relevant matrix elements:

qi = ri , pi , qei ,

G ( s ) = ( M s 2 + D s + K ) −1 .

in which ri is the vector of the coordinates that define the location of the given body in a fixed coordinate system Oxyz. Orientation of the body in the basic space is determined with Euler parameters pi. The elastic deformations of the body are expressed by the vector of elastic (or standardized modal) coordinates qei. Furthermore, in equations (10), vector λ of Lagrange multipliers is found in the number of s that have a direct connection with the reactionary forces in kinematic constraints. The coupled mechanical system is characterized by i = r – s degrees of freedom.

The method of determining the transfer function of the mechanical system according to relation (15) is applicable only for linear non-conservative mechanical systems. In the case of mechanical systems in which there occur nonlinearities due to passive resistances in kinematic constraints for example, it is advisable to use a method in which the system is excited by the swept signal. It is one of the possible methods that allows the study of nonlinearities in the time domain only. Getting a frequency response of a nonlinear dynamic system is possible only round about the operating point.

f (q, t ) = 0 . V

(11)

(15)

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

H(jω)

u(t)

If any of the eigenvalues should have a positive real part, it would be an unstable mechanical system. The relationship between the imaginary part of the eigenvalue βi and the natural frequency of the system Ω0i is determined by the following conversion:

y(t)

Fig. 4. Time-invariant system block.

Consider a dynamic system described by the frequency response: H ( jω ) =

Y ( jω ) , U ( jω )

(16)

where Y(jω) and U(jω) are the images of output and input signal y(t) and y(t), see Fig. 4. Because it is a complex function, it can be decomposed into the real and imaginary parts:

H ( jω ) = V (ω ) + jW (ω ) ,

(17)

with which the amplitude part A(ω) and the phase part φ(ω) of frequency response can be determined: A(ω ) = H ( jω ) = V 2 (ω ) + W 2 (ω ) ,  W (ω )   .  V (ω ) 

ϕ (ω ) = arg (H ( jω ) ) = arctan 

(18)

Alternatively: Re{H ( jω )} = H ( jω ) cos ϕ (ω ),

Im{H ( jω )} = H ( jω ) sin ϕ (ω ) .

(19)

If a swept sine wave signal is fed to the system input: u (t ) = A1 sin(πf 0 t 2 ),

(20)

then, a signal with variable frequency and amplitude will be stabilized at the system output:

y (t ) = A2 (t ) sin(ω (t )t + ϕ 0 ).

(21)

Subsequently, both signals are converted with Fourier transformation into the frequency domain and by their dividing according to (16), it is set the appropriate frequency response around the given operating point. In other words, there is a linear approximation of the appropriate transfer function. C. The natural frequency of the coupled mech. system If the coupled mechanical system is generally expressed in a system of equations of motion in matrix form:

&&(t ) + Dq& (t ) + Kq (t ) = f (t ), Mq

(22)

it can be solved at f(t) = 0 the problem of eigenvalues, which is given in n-dimensional space by the expression:

(

)

det λ2 M + λ D + K = 0 .

(23)

Roots λi of the characteristic equation (23) are the eigenvalues of the system. Those can be complex coupled in pairs:

λi = −α i + jβ i ,

λ∗i = −α i − jβ i ,

i = 1,2,..., m

or real:

λi = −α i ,

i = 2m + 1,2m + 2,..., n .

92

Ω 0i =

βi . 2π

(24)

V. MATHEMATICAL MODEL OF A CROSS SLIDE The aim of the calculations in the mathematical model of the chatter of the cross-slide was to calculate the limit chip width bmez and to create the stability lobe diagram in the machining processes of grooving (machining in udirection) and longitudinal turning (machining in wdirection) while both operations do not occur simultaneously. To establish them, it is necessary to determine the appropriate dynamic compliance Gy(s), y = u, w, in the given direction depending on the method of machining. As stated above, the problems of vibration during machining can be divided into the description of three basic objects: 1) The process of cutting itself. 2) The description of the mechanical system – machine tool or machine group. 3) Description of the drive, i.e., the engine itself and its control. In the first case it was used the knowledge contained primarily in [2], and which are briefly mentioned in the previous Chapter II. Here are two parameters that, on the basis of assumptions, were considered for constants that had to be put into the resulting simulation model. These are: - β = 22.5° cutting angle, - C0 = 2·109 Nm-2 cutting force constant. A cross slide mathematical model is generally expressed in equation (12), see Fig. 5. These equations were compiled based on the following assumptions. This is a spatial system of 24 perfectly rigid bodies (without considering the frame 1) and 4 flexible bodies – a console of U and W axes and a slide of U and W axes. Between the bodies there were defined kinematic pairs in such a way so that the analyzed mechanical system is free of redundant constraints. W axis console seating to the frame, ball screw support, ball screw nut and the ball screw itself of both axes were considered as flexible. Their compliance was defined by stiffness as to TABLE I. Stiffness of both ball screws is given by the ball screw nut distance from the ball screw support. These are the following: lW = 140 mm. lU = 125 mm, Both ball screws are overhung. TABLE I. STIFFNESS OF FLEXIBLE ELEMENTS Torsion stiffness of the coupling Screw support stiffness Nut axial stiffness Screw longitudinal stiffness Screw torsion stiffness

Unit [Nmrad-1] 106 [Nm-1] 106 [Nm-1] 106 [Nm-1] [Nmrad-1]

U-axis 825 450 410 688 17620

W-axis 825 450 410 614 15730

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

U-axis motor

93

Fu

U-axis slide

W-axis slide

W-axis motor

Fw

v U-axis console

w

W-axis console

u Fig. 5. Cross slide model.

The damping of these elements was selected proportional to their stiffness in the simulation model, i.e.: dU,W = 2.5·10-6kU,W in the case of longitudinal stiffness and dtU,W = 2.5·10-3ktU,W in the case of torsion stiffness. The deformation field of the flexible bodies was approximated with deformation modes, see [3]. It is a linear combination of static wave shapes of the body that respect the boundary conditions and the natural vibrations of the given flexible body. The modes of individual elements were established on the basis of the model of the particular body, created by using the FEM. Through those modes, there were determined further generalized mass matrix M and stiffness matrix K of the given flexible body whose dimension was acceptable for the purposes of subsequent calculations already. Damping matrix D of the flexible body was defined on the basis of modal damping coefficients ci. That one was introduced by means of damping ratio bri that expresses the ratio of ci modal damping coefficient in relation to cicr critical damping of the given mode of the appropriate pliable body. The coefficients of proportional damping and damping ratio were verified with respect to the experimentally determined modal damping of the natural wave shapes with the real mechanical system. TABLE II. PROPORTIONAL DAMPING OF THE OWN SHAPES OF PLIABLE BODIES W axis console W axis slide + U axis console U axis slide

bri [%] 10 25 20

The dovetail slides of the cross slide were modeled through the interaction of two bodies. It is ensured by the action of coupling forces as linear dampers with bυ damping coefficients. The origin of forces is spread over

6 points in each functional area of each from the dovetail slides of U and W machining axes. This gives a total of 2x24 force relations between two pairs of bodies. Points were distributed evenly in such a way so that the appropriate pairs of points can lie always opposite each other in the normal to the functional surfaces of the groove. Oil viscosity was estimated at 250 mPa·s and the normal clearance on one side of the groove 10 µm, while these values correspond to the specific damping coefficient δ in a size of 5·1010 Nsm-3. Then, damping coefficients bυ of one damper for each of the functional surfaces of the dovetail slide of U and W machining axes can be determined as follows: bν =

1 δ lν dν , 12

ν = Uh ,Us ,Wh ,Ws ,

(25)

in which U and W indices stand for belonging to U and W axes and h and s indices express the horizontal and inclined functional surface of the dovetail slide. Knowing length lυ and width dυ of the contact surfaces of the dovetail slide groove, damping coefficients bυ of linear dampers were determined and which are given in TABLE III. TABLE III. DAMPING COEFFICIENTS OF LINEAR DAMPERS

U W

h s h s

lυ [mm] 260 260 304 304

dυ [mm] 35.29 19.03 40.30 23.54

bυ·107 [Nsm-1] 3.8 2.1 5.1 3.0

The resulting simulation model was characterized by 108 DOF. In the model of feed drive, there were substituted the numerical values of the relevant parameters and which are given in TABLE IV.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

TABLE IV. DRIVE PARAMETERS pp – Number of pole pairs Rs – Electrical resistance of one phase Ls – Inductance of one phase KE – Motor voltage constant KT – Motor torque constant rI – Current controller proportional component TI – Current controller time integration constant rω – Speed controller proportional component Tω – Speed controller time integration constant rφ – Position regulator proportional component

Unit [-] [Ω] [mH] [Vrad] [NmA-1] [VA-1] [s] [Asrad-1] [s] [s-1]

Value 3 1.5 13.3 0.4444 0.4054 0.8 0.004 0.75 0.05 5

VI. SIMULATION RESULTS At first, the mathematical model of the cross slide was verified on the basis of measured natural frequencies on a real machine installation. By comparing the values of calculated and measured natural frequencies, a quite good compliance between those variables is apparent, see TABLE V. Since according to the measurement the third mode dominates, this model was not further verified in terms of damping due to modal damping 1 of natural oscillation shape. TABLE V. NATURAL FREQUENCIES

Measured modal Calculated Calculated modal damping [%] f0 [Hz] damping [%] 2.8 118.6 1.4 2.6 160.0 2.6 2.6 252.0 2.6

(

)

Fu , v , w = AF sin πf 0t , AF = 50 N , f 0 = 400Hz. (26) 2

To determine inertance, it is necessary to know the progress of acceleration at the site. Then, both signals are converted through Fourier transformation into the frequency domain and with their division, the appropriate course is determined. In Fig. 6 up to Fig. 8, there is shown the course of real and imaginary components of the inertances of the cross slide provided both by measurements and calculations on a mathematical model of this mechanical system further then. By comparing the courses of inertances determined by calculating in different directions with the measured ones, a very similar character of the dynamic properties of the mathematical model of the cross slide with a real object is obvious. For its further refinement it is necessary to better specify the submodel of damping caused mainly with passive resistances in the dovetail slide, which has a significant influence on the courses of inertances.

Measured real inertance component Calculated real inertance component

ReInU [m/Ns2]

1 2 3

Measured f0 [Hz] 118.0 157.0 260.0

In the next step, the mathematical model was verified with experimentally determined transfer functions of a real mechanical system. In this case, it is the course of inertances in the given direction. As an inertation it is called the transfer function between the Laplace images of cutting force and the acceleration of a tool group. In the point of acting of cutting forces, the frequency variable course of harmonic force acts in the appropriate direction according to the following equation:

f [Hz] Measured imaginary inertance component Calculated imaginary inertance component

ImInU [m/Ns2]

Shape

94

f [Hz] Fig. 6. Real and imaginary inertance component in u – direction.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

95

ReInV [m/Ns2]

Measured real inertance component Calculated real inertance component

f [Hz]

ImInV [m/Ns2]

Measured imaginary inertance component Calculated imaginary inertance component

f [Hz]

ReInW [m/Ns2]

Fig. 7. Real and imaginary inertance component in v – direction.

Measured real inertance component Calculated real inertance component

f [Hz]

ImInW [m/Ns2]

Measured imaginary inertance component Calculated imaginary inertance component

f [Hz] Fig. 8. Real and imaginary inertance component in w – direction.

96

Chip width b [mm]

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

bmezmin = 18.9 mm

Spindle speed nZ [rev/min] Fig. 9. Lobe diagram in the machining process of grooving.

By comparing the course of inertances determined by calculations in different directions with the measured ones, a very similar character of the dynamic features of the mathematical model of the cross slide with the real object is obvious. For its further refinement it is necessary to better specify the submodel of damping caused mainly by passive resistances in the dovetail slide that has a significant influence on the course of inertances. A sample of the speed stability diagram of grooving is shown in Fig. 9. The diagram was designed in accordance with the procedure referred to in Paragraph II. For its formation, transfer functions determined on the basis of calculations on a mathematical model of the cross slide were used. The theoretical value of the chip minimum limit width achieves the size bmezmin = 18.9 mm. The stable area of cutting conditions is highlighted with hatching. VII. CONCLUSION Creating a mathematical model of the cross slide as a coupled mechanical system with flexible links and real constraints are not an entirely trivial matter. Especially in the case of the dynamic properties of the dovetail slide it is important to realize that this is a non-linear dynamic system. In order to take into account the dynamic properties of the real system it is necessary to create a relatively detailed computational model which leads to a large number of degrees of freedom and often shows a nonlinear character of dynamic behavior as well. It was selected such a procedure of the creation of the linked system of the cross slide when it is a composition of abstract dynamic subsystems with causal orientation input – output. This assembly is very simple because the outputs of one model are the inputs of another model. Such a model of the related mechanical system can be solved in both time domain and frequency domain.

The course of transfer functions of the analyzed mechanical system is highly dependent on the submodel of damping - above all, on the damping in the dovetail slides of the cross slide, i.e., on passive resistances. The course of transfer functions is not significantly sensitive to the size of the specific damping coefficient δ. The shape of Lobe diagrams depends on: Cutting force constant C0, The direction of cutting force β, The course of the oriented dynamic compliance of the mechanical system. To create Lobe diagrams, it is therefore necessary to identify the internal damping of the system and the damping due to the effect of passive resistances in the given system. ACKNOWLEDGMENT This paper was created within the work on the 2A2TP1/038 – Project supported by the Ministry of Industry and Trade - Ministerstvo průmyslu a obchodu. REFERENCES [1] J. Ondrášek, “Creating a mathematical model for solving chatter and dealing the problems concerning the maximum allowable size of a machining chip”, Advances in Mechanisms Design, Proceedings of TMM 2012, pp. 237-243, Springer Sciens+Business Media Dordrecht, 2012, ISBN 978-94-007-5124-8. [2] P. Souček, A. Bubák, Vybrané statě z kmitání v pohonech výrobních strojů, (in Czech), Česká technika – nakladatelství ČVUT, 2008, ISBN 978-80-01-04048-5 [3] J. Slavík, V. Stejskal, V. Zeman, Základy dynamiky strojů, (in Czech), Vydavatelství ČVUT, 1997, ISBN 80-01-01622-6. [4] P. Souček, Servomechanismy ve výrobních strojích, (in Czech), Vydavatelství ČVUT, 2004, ISBN 80-01-02902-6. [5] M. Valášek, kolektiv, Mechatronika, (in Czech), Vydavatelství ČVUT, 1995, ISBN 80-01-01276-X.

Transactions on Electrical Engineering, Vol. 2 (2013), No. 3

Zablodskiy, N., Lettl, J., Pliugin, V., Buhr, K., Khomitskiy, S.: Induction Motor Design by Use of Genetic Optimization Algorithms The problem of the automated calculation and optimal design of an induction motor is presented. The problem of optimization by use of genetic algorithms is set and solved. The analysis of the obtained results is executed.

Homišin, J.: New Ways of Controlling Dangerous Torsional Vibration in Mechanical Systems In general terms the mechanical systems (MS) mean the systems of driving and driven machines arranged to perform the required work. We divide them into MS operating with constant speed and MS working within a range of speed. In terms of MS dynamics we understand the system of masses connected with flexible links between them, i.e. systems that are able to oscillate. Especially piston machines bring heavy torsional excitation to the system, which causes oscillation, vibration, and hence their noise. Based on the results of our research, the torsional vibration control as a source of given systems excitation, can be achieved by application of pneumatic couplings, or pneumatic tuners of torsional vibration. Existence of pneumatic couplings and pneumatic tuners of torsional vibration developed by us, creates the possibility of implementing new ways of torsional oscillating mechanical system tuning. Based on the above, the aim of the article is to present new ways of dangerous torsional vibration control of mechanical systems by application of pneumatic couplings and pneumatic tuners of torsional vibration developed by us.

Folta, Z., Hrudičková, M.: Experience with Torque Measurement on Rotating Shaft This article deals with problem of torque measurement on the turning shaft and the signal transmission to the measuring equipment. There is described the experience with contact transmission and both non-contact and high frequency transmission (with ESA Messtechnik equipment). Also the examples of the practical realization are presented.

Veleba, J..: Performance of Steady-State Voltage Stability Analysis in MATLAB Environment In recent years, electric power systems have been often operated close to their working limits due to increased power consumptions and installation of renewable power sources. This situation poses a serious threat to stable network operation and control. Therefore, voltage stability is currently one of key topics worldwide for preventing related black-out scenarios. In this paper, modelling and simulations of steady-state stability problems in MATLAB environment are performed using author-developed computational tool implementing both conventional and more advanced numerical approaches. Their performance is compared with the Simulink-based library Power System Analysis Toolbox (PSAT) in terms of solution accuracy, CPU time and possible limitations.

Ondrášek, J.: Cross Slide Mathematical Model for Solving Chatter The paper deals with the issues of creating a mathematical model of self-excited chatter of the cross slide. This type of oscillation occurs in systems where there is still an internal source from which the system consumes energy to maintain or even to increase the amplitude of oscillation. This consumption is controlled by an oscillatory movement of the system itself. Such an energy source greatly affects the dynamic and stability properties of the system.

_____________________________________________________________________________________________ TRANSACTIONS ON ELECTRICAL ENGINEERING VOL. 2, NO. 3 HAS BEEN PUBLISHED ON 30TH OF SEPTEMBER 2013