Strojniški vestnik Journal of Mechanical Engineering Strojniški vestnik - Journal of Mechanical Engineering (SV-JME) Aim and Scope The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis. The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue. The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s). Editor in Chief Vincenc Butala University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Technical Editor Pika Škraba University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Founding Editor Bojan Kraut University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerčeva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www.sv-jme.eu Print: Grafex, d.o.o., printed in 310 copies Founders and Publishers University of Ljubljana, Faculty of Mechanical Engineering, Slovenia University of Maribor, Faculty of Mechanical Engineering, Slovenia Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association President of Publishing Council Branko Širok University of Ljubljana, Faculty of Mechanical Engineering, Slovenia International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, University of Engineering and Technology, Pakistan Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, UL, Faculty of Mechanical Engineering, Slovenia Franci Čuš, UM, Faculty of Mechanical Engineering, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Igor Emri, UL, Faculty of Mechanical Engineering, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Janez Grum, UL, Faculty of Mechanical Engineering, Slovenia Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, UM, Faculty of Mechanical Engineering, Slovenia Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan Julius Kaplunov, Brunel University, West London, UK Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany Jernej Klemenc, UL, Faculty of Mechanical Engineering, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, UL, Faculty of Mechanical Engineering, Slovenia Gorazd Lojen, UM, Faculty of Mechanical Engineering, Slovenia Thomas Lübben, University of Bremen, Germany Janez Možina, UL, Faculty of Mechanical Engineering, Slovenia George K. Nikas, KADMOS Engineering, UK José L. Ocana, Technical University of Madrid, Spain Miroslav Plančak, University of Novi Sad, Serbia Vladimir Popović, University of Belgrade, Faculty of Mech. Eng., Serbia Franci Pušavec, UL, Faculty of Mechanical Engineering, Slovenia Bernd Sauer, University of Kaiserlautern, Germany Rudolph J. Scavuzzo, University of Akron, USA Arkady Voloshin, Lehigh University, Bethlehem, USA Vice-President of Publishing Council Jože Balič University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: Front cover shows a close-up of an optical system, developed for inspection of microlens arrays. The system is capable of measuring the shape of the generated wavefronts and analysing the presence of optical aberration in terms of Zernike polynomial coefficients. Image Courtesy: Laboratory for Digital Systems and Electrical Engineering, Faculty of Mechanical Engineering, University of Ljubljana, Slovenia Photo: Željko Stevanić/ifp ISSN 0039-2480 General information Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue). Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor. Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http://en.sv-jme.eu/. You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content. We would like to thank the reviewers who have taken part in the peerreview process. © 2016 Strojniški vestnik - Journal of Mechanical Engineering. All rights reserved. SV-JME is indexed / abstracted in: SCI-Expanded, Compendex, Inspec, ProQuest-CSA, SCOPUS, TEMA. The list of the remaining bases, in which SV-JME is indexed, is available on the website. The journal is subsidized by Slovenian Research Agency. Strojniški vestnik - Journal of Mechanical Engineering is available on http://www.sv-jme.eu, where you access also to papers' supplements, such as simulations, etc. Contents Strojniški vestnik - Journal of Mechanical Engineering volume 62, (2016), number 10 Ljubljana, October 2016 ISSN 0039-2480 Published monthly Papers Afrasiab Raisi: Natural Convection of Non-Newtonian Fluids in a Square Cavity with a Localized Heat Source 553 Zhenxing Kong, Dawei Pi, Xianhui Wang, Hongliang Wang, Shan Chen: Design and Evaluation of a Hierarchical Control Algorithm for an Electric Active Stabilizer Bar System 565 Yongxiang Li, Xifan Yao, Jifeng Zhou: Multi-objective Optimization of Cloud Manufacturing Service Composition with Cloud-Entropy Enhanced Genetic Algorithm 577 Rade Vasiljević, Milomir Gašić, Mile Savković: Parameters Influencing the Dynamic Behaviour of the Carrying Structure of a Type H Portal Crane 591 Cuneyt Uysal, Kamil Arslan, Huseyin Kurt: A Numerical Analysis of Fluid Flow and Heat Transfer Characteristics of ZnO-Ethylene Glycol Nanofluid in Rectangular Microchannels 603 Xiaozhou Li, Manlin Zhu, Jiancang Xie: Numerical Simulation of Transient Pressure Control in a Pumped Water Supply System Using an Improved Bypass Pipe 614 Strojniški vestnik - Journal of Mechanical Engineering 62(2016)10, 553-564 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.3218 Original Scientific Paper Natural Convection of Non-Newtonian Fluids in a Square Cavity with a Localized Heat Source Afrasiab Raisi Shahrekord University, Faculty of Engineering, Iran This study investigates the natural convection cooling of a heat source placed on the bottom of a cavity filled with non-Newtonian power-law fluids. The heat source produces a uniform heat flux. The remaining parts of the bottom wall of the cavity are insulated, and the vertical, and top walls are kept at a low temperature. The governing equations for the power-law fluid flow are solved with the numerical finite difference method based on the control volume formulation and SIMPLE algorithm. The study investigates the effects of relevant parameters such as the Rayleigh number, the power-law index and the heat source length and heat source location on the thermal performance of the cavity. The results show that as the Rayleigh number increases, particularly when n<1, the thermal performance of the cavity is improved. Furthermore, using non-Newtonian shear thinning fluids, especially at higher Rayleigh numbers, improves the cooling performance of the heat source. The results also show that, depending on the Rayleigh number and the power-law index, the length and position of the heat source have significant effects on the thermal performance of the cavity. Keywords: natural convection, cavity, heat source, non-Newtonian, power-law Highlights • The numerical simulation of natural convection of non-Newtonian power-law fluids. • The thermal performance of an enclosure filled with a non-Newtonian fluid. • The thermal performance of the enclosure is enhanced as the Rayleigh number increases. • The average Nusselt number increases for shear-thinning and reduces for shear-thickening fluids compared to Newtonian fluids. 0 INTRODUCTION Natural convection is often considered to be the main mechanism of heat transfer in cavities. It has a variety of applications in many engineering systems. The fluid inside the cavity can occasionally be regarded as Newtonian; however, in many natural or artificial systems, this fluid has a non-Newtonian behaviour. With their important applications in engineering, fluids such as nanofluids, molten polymers, paints, food stuff, inks, organic matter, and glues may all exhibit strong non-Newtonian behaviours. Examining natural convection in non-Newtonian fluids in cavities is, therefore, highly crucial to engineering and has many applications in it. Food processing, oil drilling, polymer engineering, geophysical systems, electronic cooling systems, and nuclear reactors are examples of natural convection in non-Newtonian fluids. Natural convection in Newtonian fluids has been extensively studied in cavities under various boundary conditions. Catton [1] examined natural convection in vertical, horizontal and tilted rectangular enclosures. Davis [2] studied natural convection inside an enclosure with isothermal sidewalls and adiabatic horizontal walls. A comprehensive study of natural convection literature was done by Ostrach [3]. The applications of natural convection were presented by Bejan [4] and Emery and Lee [5]. Aydin [6] assessed natural convection inside a rectangular enclosure, which was heated from the left vertical wall and was cooled from the top wall; the other two walls were adiabatic. The thermal performance of rectangular enclosures is dependent on the aspect ratio of the enclosures [7]. The conducted studies on the natural convection field have been reviewed by Ganguli et al. [8]. Arpino et al. [9] analysed transient natural convection in tall cavities. They used a dual time-stepping to improve the transient solution. In another study, Arpino et al. [10] numerically studied transient natural convection in porous and partially porous cavities with an emphasis on the dependence of flow behaviour on Ra, porous layer permeability, and cavity aspect ratio [1] to [10]. Ozoe and Churchill [11] were perhaps the first to study forced convection in non-Newtonian fluids in cavities. They investigated natural convection in two non-Newtonian fluids, i.e. Ostwald-de Waele (power-law) and Ellis, in a shallow horizontal cavity heated from the bottom and cooled from above. They found that the critical value of Rayleigh number for the onset of natural convection increases with the fluid index. Cavities with natural convection are divided into two general groups. The first group includes cavities that are heated and cooled from the side walls, and the second group includes cavities that are heated from the bottom. In the second group, natural convection begins when the Rayleigh number exceeds a critical value. Kaddiry et al. [12] conducted a numerical investigation of a Rayleigh- Benard convection of a non-Newtonian power-law fluid in a square cavity. They showed that the critical value of the Rayleigh number for onset natural convection was an increasing function of the power-law index. For example, the results indicated that Rac~600 for n=0.6 and Rac ~ 5000 for n = 1.4. Kim et al. [13] studied transient natural convection in non-Newtonian power-law fluids in a vertical cavity by assuming the horizontal walls to be insulated and taking the cause of the buoyancy force and natural convection inside the cavity to be the simultaneous temperature change at the vertical walls. They found that, at a certain value of the Rayleigh number, both the convection strength and the heat transfer rate change for non-Newtonian fluids in comparison to Newtonian fluids, so that they are both increased at n < 1 and reduced at n > 1. Lamsaadi et al. [14] studied transient natural convection in non-Newtonian power-law fluids in a shallow cavity using numerical and analytical methods. Their cavity consisted of long, insulated, horizontal walls and short vertical walls, with the vertical walls heated and cooled with a constant heat flux. The results of their study showed that heat transfer and flow characteristics were not sensitive to an increase in the Prandtl number or the aspect ratio of the cavity, given that the parameters were sufficiently large. They argued that n (the power-law index) and the Rayleigh number were parameters that affected the rate of heat transfer and flow field in shallow cavities with non-Newtonian fluids and large Prandtl numbers. Lamsaadi et al. [15] also investigated steady natural convection in non-Newtonian power-law fluids in a tilted rectangular slit using the numerical method. A constant heat flux was applied to the side walls, and the remaining walls were insulated. This study was conducted within the Rayleigh number range of 10 < Ra < 105, the power-law index of 0.6 < n < 1.4 and the rotation angle of -180° < Ф < 180°. They found that given a certain value of Rayleigh number, the rotation of the cavity had a dramatic effect on the rate of heat transfer. The highest rate of heat transfer occurred when the cavity was heated from the bottom, and the effect of the cavity rotation increased as the value of n decreased. Several studies have recently investigated natural convection in non-Newtonian fluids in square and rectangular cavities with horizontal insulated walls and vertical walls with constant temperatures. Turan et al. [16] used the FLUENT ANS YS commercial package to simulate power-law fluids. Ternik and Rudolf [17] investigated the natural convection of a non-Newtonian nanofluid in an enclosure with isothermal sidewalls and adiabatic horizontal walls. The natural convection of power-law fluids has been assessed in a rectangular enclosure [18] and in a square enclosure [19] with isothermal vertical walls. The common findings of these studies showed that the heat transfer and flow field were affected by the Rayleigh number and n (the power-law index). However, for any given n and Ra, the average Nusselt number was not affected by the Prandtl number. Moreover, for rectangular cavities with a boundary condition of a constant temperature on the vertical walls, the Nusselt number did not vary uniformly as the aspect ratio increased; however, when the boundary condition was a constant heat flux applied to the vertical walls, this variation was uniform. Alloui [20] studied natural convection in a non-Newtonian fluid in a vertical cavity using the Carreau-Yasuda model and showed that compared to Newtonian fluids, the rate of heat transfer and convection strength increased when n < 1 and decreased when n > 1. The changes in the n value did not greatly affect the rate of heat transfer when the Rayleigh number was small. Natural convection in non-Newtonian fluids in cavities heated from the bottom has received less attention in literature than the studies considering cavities heated from the sides. Lamsaadi et al. [21] studied natural convection in a horizontal cavity heated from the bottom and cooled from above using numerical and analytical methods. They showed that fluid flow, temperature distribution and the rate of heat transfer were sensitive to n (the power-law index); however, for large Prandtl numbers, they were no longer sensitive to the Prandtl number. Compared to Newtonian fluids, non-Newtonian fluids were associated with the onset of natural convection by a single cell flow when n < 1 and behaved inversely when n > 1. Ouertatani et al. [22] conducted a numerical assessment of Rayleigh-Bénard convection in a rectangular cavity; they considered a cavity with the bottom wall at a constant temperature of Th and the top wall at a constant temperature of Tc. They assumed the vertical walls to be insulated. Their study yielded streamlines and isotherms for the Rayleigh numbers in the range of 104 < Ra < 106. They argued that the rate of heat transfer increased with the Rayleigh number. Mahrood Khadangi et al. [23] conducted an experimental study of natural convection of a non-Newtonian nanofluid in a cavity heated from the bottom with a constant heat flux and cooled from above. They used the Al2O3 and TiO2 nanoparticles and a 0.5 weight percent solution of carboxymethyl cellulose as the base fluid for the nanofluid. The results of their study showed improvements in the rate of heat transfer of natural convection in non-Newtonian nanofluids at low concentrations of nanoparticles. However, the rate of heat transfer began to decline as the volume percent of nanoparticles exceeded a certain amount. Khazzar et al. [24] conducted a numerical study of natural convection in a non-Newtonian power-law fluid in a tilted rectangular cavity heated from the bottom and cooled from above. They conducted this study for different tilt angles of the cavity (0 < ф < 90) and calculated the rate of heat transfer and flow field for different aspect ratios, Rayleigh numbers and power-law indices. Their results showed an increasing and a decreasing average Nusselt number for non-Newtonian shear-thinning (n < 1) and non-Newtonian shear-thickening (n > 1) fluids, in respective order. The increase and decrease in the average Nusselt number depended on the Rayleigh number, the Prandtl number, the aspect ratio, and the power-law index. Alloui et al. [25] conducted a numerical and analytical investigation on the onset of natural convection in a shallow rectangular cavity filled with power-law fluids. The cavity was heated from the bottom through a constant heat flux. They calculated the critical value of the Rayleigh number for non-Newtonian fluids given n > 1 and n < 1 using the linear stability theory. Turan et al. [26] conducted a numerical study of a laminar natural convection of power-law fluids in a square cavity heated from the bottom through a constant heat flux. They argued that the average Nusselt number was significantly affected by the Rayleigh number and n (the power-law index); however, it was not affected significantly by the Prandtl number. They also compared bottom-heated and side-heated cavities and concluded that for shear-thinning fluids and at high Rayleigh numbers, the average Nusselt number was higher for side-heated cavities in comparison to bottom-heated ones. However, for shear-thickening and Newtonian fluids, the Nusselt number remained approximately the same for both side-heated and bottom-heated cavities. Yigit et al. [27] studied the effects of aspect ratio on natural convection in non-Newtonian Bingham fluids in bottom-heated and top-cooled rectangular cavities. They varied the aspect ratio in the range of (0.25 1) becomes more complicated due to the formation of the Banard cells. Despite the valuable studies conducted on the natural convection of non-Newtonian fluids in bottom-heated cavities as shown in the review of the literature, no studies were found on cooling a local heat source placed on the bottom of a cavity by power-law fluids. The heat source in mind can be an electronic device producing heat while operating in an electronic circuit, which requires cooling so as not to exceed a certain permitted temperature. As regards the length of the heat source, as well as the position of its placement on the floor of the cavity, there can be significant effects on the heat transfer and flow of fields; in this study, the effects of the length of the heat source and its position are also reviewed. The present study is, therefore, conducted to examine the effects of parameters such as the Rayleigh number, the power-law index and the positioning and length of the heat source on the rate of heat transfer and temperature and flow fields. 1 METHODS 1.1 Governing Equations It is assumed that the enclosure is full of a non-Newtonian power-law fluid. The flow is considered to be steady, two-dimensional and laminar, and the radiation effects are negligible. Also, density changes have been modeled in the momentum equation with the Boussinesq approximation. Based on the stated assumptions, the dimensional governing equations of continuity, momentum, and energy are as follows: du dv — + — = 0, dx dy (1) .du du. dp дтхх дт p(u— + v—) = —— +—^ + —(2) дх ду дх dx dy , dv dv dp dTxy , dTyy , агт т\ dx dy dy dx dy dT dT u--+ v— = a dx dy fd 2T d 2ТЛ ~dx2 +~dy2 (4) For a non-Newtonian fluid which follows the power-law model the viscous stress tensor is given by [29]: Tj = 2HaDj = Ha (du. du — +—L dxj dxi (5) where Dj is the rate of strain tensor for the two-dimensional Cartesian coordinate and fia is the apparent viscosity that is derived for the two-dimensional Cartesian coordinate as: Va = K \ 2 du Y + f dv dx J ^dy rdv + duf Ì 2 dx dy J (6) where K and n are power-law model constants. K is the consistency coefficient and n is the power-law index. Where (n < 1), is for shear-thinning fluids and (n > 1), is for shear-thickening fluids. When (n = 1), a Newtonian fluid is obtained. By using the following non-dimensional parameters: X = f Y = U = UL, V = L L a P = pL2 в = T - T pa' AT gßATL2n+1 AT = a q"L Ra = anK / p Pr = k KL2"-2 pa (7) Eqs. (1) to (4) are converted to dimensionless forms, as follows: dU dV n -+-= 0, dX dY TTdU TrdU dP U-+ V-=--- dX dY dX (8) + Pr 2 d , du 2-1 U„- dX i" dX —Irt,— l+T-l И* ^7 + ^7 I I, (9) _8_( , (dU_ dV_ dY V" i dY + dX jjdV dV dP U-+ V-=--+ Pr dX dY dY d ( * (dU_ d¥_ + dX \Иа i dY + dX 2—( * —1+ dY dY J + + RaPrQ, (10) TTde „dO dO d2в U— + V— = —- + —-dX dY dX2 dY2 where /л* is the apparent dimensionless viscosity and is defined as: / = i 2 dU dX +(dV I dY ( dV dU +1-+- [dX dY n-1 2 i T — + — I } . (12) 1.2 Problem Description A schematic diagram of the problem under examination is shown in Fig. 1. The figure illustrates a two-dimensional cavity containing a non-Newtonian power-law fluid. A heat source is placed on the bottom wall of the cavity, producing a uniform heat flux. The remaining parts of the bottom walls are insulated, and the vertical side walls and the top horizontal walls are kept at a relatively low temperature of Tc. The cavity contains a shear-thinning fluid (n < 1), a Newtonian fluid (n = 1), and a shear-thickening fluid (n > 1). The Prandtl number is assumed constant and equal to 100. Fig. 1. Schematic diagram of the physical model According to the schematic diagram of the physical model (Fig. 1), the dimensionless boundary conditions are as follows: on the left and right walls: U = V= в=0, on upper horizontal wall: U = V= в=0, on the edges of the lower wall: U = V= дв/dY = 0, on the heat source: U= V= 0 and дв/öY=-1. (13) After solving the governing equations numerically, as a measure of the heat transfer rate of the enclosure, the local Nusselt number on the heat source surface can be defined as follows: n-1 + hL q"L Nu = — = —-1——. k (T -T )k (14) In Eq. (14), h is the convection heat transfer coefficient. Using the dimensionless parameters, the following relationship is obtained for the local Nusselt number: Nu = (15) ds (X) where 9s is the dimensionless temperature of the heat source. The average Nusselt number can be obtained by integrating the local Nusselt number along the heat source. S+-W 1 2 1 Nu„ = — J , ^ dX. W S - -W 2 0. (X ) (16) between the predictions obtained with mesh M2 and extrapolated values is extremely good for both the average Nasselt number and the maximum temperature of the heat source. The grid size 100*100 is, therefore, used in the numerical solution. Table 1. Grid independence study (Ra = 106, n = 0.6, W = 0.4, S= 0.5). Nx*Ny NUm Error [%] n (7s,max Error [%] Mesh M1 50*50 26.0491 1.6 0.06431 4 Mesh M2 100*100 25.7902 0.6 0.0669 0.194 Mesh M3 200*200 25.6755 0.15 0.067 0.045 fext 25.6372 0.06703 Fig. 2 shows a quarter of the computational grid which is made up of uniform meshes, and Fig. 3 shows the staggered grid in the x direction. 1.3 Numerical Method, Grid Study and Code Validation For obtaining a numerical solution, the governing differential equations have to be converted into algebraic equations. The governing dimensionless equations (Eqs. (8) to (11)) are integrated over each of the finite control volumes in the flow domain; obtained integrated transport equations along with the relevant boundary conditions are discretized by using finite difference method. The convection-diffusion terms are estimated using the power-law scheme. For obtaining a numerical solution to the discretized equations, SIMPLE algorithm proposed by Patankar [30] is used and a computer program is written in FORTRAN. The convergence criterion is to reduce the maximum mass residual of the grid control volume below 10-8. A uniform staggered grid is used for obtaining the numerical solution. The grid independence is examined for n = 0.6, Ra = 106, S= 0.5 and W = 0.4. The results of grid study for three different uniform meshes are presented in Table 1. Using the Richardson extrapolation, the grid-converged value of a general variable such as ф is obtained as follows. Vex, = Vm 3 V M 2 -V.M 3 rp -1 (17) where фм3 is obtained on the basis of the finest grid and фм2 is thus based on the next level of coarse grid, r=2 is the ratio between the coarse to fine grid spacing and p=2 is the theoretical accuracy. The results of the grid study show that as the grid points increase from 50*50 to 200*200 the error of the numerical solution decreases, and the agreement Fig. 2. A quarter of the computational grid Fig. 3. Staggered grid in the x direction For validation of the numerical code, the results of the present code are compared with the results obtained by Turan et al. [16]. For this purpose, a square cavity filled with a non-Newtonian power-law fluid is considered, as are insulated horizontal walls, and vertical walls that are at constant and different temperatures of Th and Tc. The average Nusselt number is calculated for different Rayleigh numbers and power-law indices. Fig. 4 compares the results and shows a good consistency between the results of the present study and of the one by Turan et al. [16]. 10" 10s Ra Fig. 4. Validation of the present code against Turan et al. [16] 2 RESULTS AND DISCUSSION The results of the present study are shown in the form of the effects of the Rayleigh number (103 1, the average Nusselt number is not affected by changes in n. At high Rayleigh numbers (Ra = 105, 106), the average Nusselt number increases due to an enhanced convection; in these states, an increase in n reduces the average Nusselt number significantly. Fig. 9. Variation of average Nusselt number versus power-law index at various Rayleigh numbers (W = 0.4, S = 0.5, Pr = 100) For a better understanding of the variations of Num with Ra and the power-law index (n), Table 2 presents the percentage changes in the average Nusselt number of non-Newtonian fluids compared to Newtonian fluids at various Rayleigh numbers. Table 2. The percentages of average Nusselt number variation for non-Newtonian fluids (W = 0.4, S = 0.5, Pr = 100) Г Num - Nu Л m,n=l xl00 v Num,n= 1 J Ra 103 104 105 106 n = 0.6 0.02 57.2 114.3 92.44 n = 0.8 0 10.6 45.67 40.03 n = 1.2 -0.01 -0.96 -24.68 -27.26 n = 1.4 -0.01 -1.31 -35.63 -45.41 n = 1.8 -0.02 -1.54 -40.67 -63.28 In Ra= 103, the percentage of changes in the average Nusselt number is negligible in non-Newtonian fluids compared to Newtonian fluids. In general, the percentage of changes in the average Nusselt number is positive for shear-thinning non-Newtonian fluids (n< 1) and negative for shear-thickening non-Newtonian fluids (n > 1) compared to Newtonian fluids. At increased Rayleigh numbers and an enhanced convection inside the cavity, the percentage of changes in the average Nusselt number increases in non-Newtonian fluids compared to in Newtonian fluids. The percentage of changes in the average Nusselt number is highest for shear-thinning fluids compared to Newtonian fluids at Ra = 105, mainly because the dominant mechanism of heat transfer changes from conduction to convection at this Ra value and the reduction in viscosity has the most significant enhancing effect on convection. 3.2 The Effects of the Heat Source Length In this section, the power-law index and the position of the heat source were taken as constant and equal to n = 0.6 and S = 0.5, and the effects of the heat source length (0.2 < W< 0.8) on the flow and temperature fields and thermal performance of the cavity were then investigated. Fig. 10 presents the streamlines (left) and the isotherms (right) at Ra = 105 for various lengths of the heat source. The strength of the rotating cells is increased due to the greater heat generated by the heat source's increased length. The isotherms indicate that convection is the dominant mechanism of heat transfer. Although convection is boosted with the increase in the length of the heat source, the maximum heat source temperature has not been reduced. In contrast, since a greater heat is generated in the heat source by the increase in its length, the heat source temperature has also increased. Fig. 11 shows the effects of the Rayleigh number and the length of the heat source on the average Nusselt number (Num) and the maximum temperature of the heat source. As expected, the average Nusselt number increases with an increase in the Rayleigh number and the subsequently enhanced convection. The average Nusselt number decreases slightly with the increase in the length of the heat source. When the length of the heat source increases, the fluid is in contact with the heat source for a prolonged period, thereby leading to an increase in fluid temperature and a reduction in the temperature gradient next to the heat source. As a result, Num decreases slightly with the increase in the length of the heat source and despite Fig. 10. Streamlines (left) and isotherms (right) for different heat source lengths (Ra = 105, n = 0.6, S = 0.5, Pr = 100) a) 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 --Ra = 10 ----Ra = 105 --Ra = 10 0.5 w b) Fig. 11. Effect of heat source length and Rayleigh number on; a) average Nusselt number; and b) heat source maximum temperature (n = 0.6, S = 0.5, Pr = 100) the enhanced convection. The reduction in the average Nusselt number reduces the cooling power of the fluid; as a result, the maximum temperature of the heat source increases with the length of the heat source. The increase in the maximum temperature of the heat source occurs with greater intensity at Ra= 103, as at this Rayleigh number, heat is transferred through conduction, and an increase in the length of the heat source leads to the generation of more heat, which then has to be let out by conduction. 3.3 Effects of the Positioning of the Heat Source In this section, the Rayleigh number is taken as Ra = 105 and the length of the heat source as W = 0.4, and through changing the position of the heat source (0.2 < S < 0.5) for three values of the power-law index (n=0.6, 1, 1.4), the flow and temperature fields and the rate of heat transfer have been examined. Fig. 12. Streamlines for different locations of the heat source at three different power-law indexes (Ra = 105, W = 0.4, Pr = 100) Figs. 12 and 13 show the streamlines and the isotherms for four different positions of the heat source (S = 0.2, 0.3, 0.4, 0.5) and three values of n. When the heat source is placed next to the cold wall on the left, two asymmetric counter-rotating vortices are formed in the cavity. The strength of the left vortex is lower than the right one due to the limitations of the space in which the fluid rotates and because it is less in contact with the heat source. Through distancing, the heat source from the left wall, the strength of the left vortex increases while the strength of the right vortex decreases following a brief phase of increase. Finally, the strength of both vortices becomes equal when the heat source is in the mid-section of the bottom wall of the cavity. Fig. 13. Isotherms for different locations of the heat source at three different power-law indexes (Ra = 105, W = 0.4, Pr = 100) Isotherms show different behaviours with power-law index changes (Fig. 13). At n = 0.6, isotherms indicate that a considerable convection strength exists inside the cavity for different positions of the heat source and also that convection is the main mechanism for heat transfer. Distancing the heat source from the left wall enhances convection and reduces the maximum temperature of the heat source. At n = 1, both conduction and convection mechanisms play a role in heat transfer. By distancing the heat source from the cold left wall, less heat is conducted out and the maximum temperature of the heat source, therefore, increases, until the heat source reaches the mid-section of the bottom wall of the cavity, where cold fluid is in contact with the heat source from both sides, thereby reducing the maximum temperature of the heat source. Finally, at n = 1.4, where conduction is the main mechanism of heat transfer, distancing the heat source from the left wall increases the maximum temperature of the heat source. Table 3 presents the average Nusselt number for four different positions of the heat source and three values of the power-law index. At n = 0.6, at which a strong convection is generated inside the cavity, the average Nusselt number does not vary significantly with the position of the heat source. At this state, distancing the heat source from the left wall reduces the average Nusselt number slightly and then, when two vortices with equal strength are formed inside the cavity (5 = 0.5), the average Nusselt number increases slightly. At n = 1, where conduction and convection are equally powerful mechanisms for heat transfer, distancing the heat source from the left wall slightly reduces the average Nusselt number at first due to the reduced rate of conduction heat transfer, and when two vortices with equal strength are formed, the average Nusselt number increases. At n = 1.4, where conduction is the main mechanism of heat transfer, the average Nusselt number decreases with the heat source distancing from the left wall and approaching the mid-section of the bottom wall. Table 3. Average Nusselt number for different locations of the heat source at three different power-law indexes (Ra = 105, W = 0.4, Pr = 100) Num n = 0.6 n = 1 n = 1.4 5 = 0.2 14.5487 9.3163 7.9428 5= 0.3 14.5450 7.0014 5.0477 S= 0.4 14.5148 6.9105 4.6306 5= 0.5 14.9846 7.0234 4.5285 4 CONCLUSIONS This study examines the natural convection in a cavity filled with a non-Newtonian power-law fluid and partially heated from the bottom through a source of uniform heat flux. The effects of the Rayleigh number, the power-law index for non-Newtonian fluids and the length and position of the heat source on the thermal performance of the cavity were investigated. The following conclusions are drawn based on the results obtained. Increasing the Rayleigh number enhances the natural convection inside the cavity, resulting in a higher rate of heat transfer and a reduced temperature for the heat source. At n < 1, the natural convection is further enhanced when the Rayleigh number is increased. Reducing the power-law index (n) reduces the apparent viscosity of the fluid and enhances natural convection inside the cavity, which then leads to an increased rate of heat transfer and a reduced temperature for the heat source. Enhancing natural convection as a result of reducing n is more pronounced at higher Rayleigh numbers. The average Nusselt number increases for shear-thinning and decreases for shear-thickening fluids compared to Newtonian fluids. These changes are more pronounced for shear-thinning fluids (n < 1). Increasing the length of the heat source enhances natural convection inside the cavity, but simultaneously increases the temperature of the heat source and also reduces its average Nusselt number. When the heat source approaches the mid-section of the bottom wall from the left wall of the cavity, the thermal behaviour of the cavity varies depending on whether the value of the power-law index is greater than, less than or equal to unity. 5 NOMENCLATURE Dj rate of strain tensor [s1] g gravitational acceleration [ms-2] h convection heat transfer coefficient [Wm-2K-1] k thermal conductivity [Wm-1K-1] K consistency coefficient [Nsnm-2] L Cavity length [m] n power-law index [-] Nu Nusselt number [-] Num average Nusselt number [-] p fluid pressure [Nm-2] P modified pressure [Nm-2] P dimensionless pressure[-] Pr Prandtl number[-] q" heat flux [Wm-2] Ra Rayleigh number[-] 5 heat source distance from the left wall [m] S dimensionless distance of heat source from the left wall [-] T temperature [K] u,v velocity components in x, y directions [ms-1] U, V dimensionless velocity components [-] w length of the heat source [m] W dimensionless length of the heat source [-] X, y Cartesian coordinates [m] X Ydimensionless coordinates^] Greek symbols a thermal diffusivity [m2s-1] ß thermal expansion coefficient [K-1] AT reference temperature difference [K] Tj stress tensor[Nm-2] в dimensionless temperature[-] H dynamic viscosity [Nsm-2] Ha dimensionless apparent viscosity [-] p density [kgm-3] Subscripts a apparent c cold m average 5 surface of the heat source 6 REFERENCES [1] Catton, I. (1978). Natural convection in enclosures. Proceedings of the 6th International Heat Transfer Conference, vol. 6, p. 13-31. [2] De Vahl Davis, G. (1983). Natural convection of air in a square cavity: a bench mark numerical solution. International Journal for Numerical. Methods in Fluids, vol. 3, no. 3, p. 249-264, D0l:10.1002/fld.1650030305. [3] Ostrach, S. (1988). Natural convection in enclosures. Journal of Heat Transfer, vol. 110, no. 4b, p. 1175-1190, D0I:10.1115/1.3250619. [4] Bejan, A. (2004). Convection Heat Transfer, 3rd ed. John Wiley, New York. [5] Emery, A.F., Lee, J.W. (1999). The effects of property variations on natural convection in a square cavity. Journal of Heat Transfer, vol. 121, no. 1, p. 57-62, D0I:10.1115/1.2825966. [6] Aydin, O., Ünal, A., Ayhan, T. (1999). Natural convection in rectangular enclosures heated from one side and cooled from above. International Journal of Heat and Mass Transfer, vol. 42, no. 13, p. 2345-2355, D0I:10.1016/S0017-9310(98)00319-6. [7] Turan, O., Poole, R.J., Chakraborty, N. (2012). Influences of boundary conditions on laminar natural convection in rectangular enclosures with differentially heated side walls. International Journal of Heat and Fluid Flow, vol. 33, no. 1, p. 131-146, DOI:10.1016/j.ijheatfluidflow.2011.10.009. [8] Ganguli, A.A., Pandit, A.B., Joshi, J.B. (2009). CFD simulation of heat transfer in a two dimensional vertical enclosure. Chemical Engineering Research and Design, vol. 87, no. 5, p. 711-727, DOI:10.1016/j.cherd.2008.11.005. [9] Arpino, F., Cortellessa, G., Dell'Isola, M., Massarotti, N., Mauro, A. (2014). High order explicit solutions for the transient natural convection of incompressible fluids in tall cavities. Numerical Heat Transfer, Part A-Application: An International Journal of Computation and Methodology, vol. 66, no. 8, p. 839-862, D0I:10.1080/10407782.2014.892389. [10] Arpino, F., Cortellessa, G., Mauro, A. (2015). Transient thermal analysis of natural convection in porous and partially porous cavities. Numerical Heat Transfer, Part A-Application: An International Journal of Computation and Methodology, vol. 67, no. 6, p. 605-631, DOI:10.1080/10407782.2014.949133. [11] Ozoe, H., Churchill, S.W. (1972). Hydrodynamic stability and natural convection in Ostwald-de Waele and Ellis fluids: The development of a numerical solution. AIChE Journal, vol. 18, no. 6, p. 1196-1207, D0I:10.1002/aic.690180617. [12] Kaddiri, M., Naimi, M., Raji, A., Hasnaoui, M. (2012). Rayleigh-Benard convection of non-Newtonian power-law fluids with temperature-dependent viscosity. ISRN Thermodynamics, vol. 2012, Article ID 614712, DOI:10.5402/2012/614712. [13] Kim, G.B., Hyun, J.M., Kwak, H.S. (2003). Transient buoyant convection of a power-law non-Newtonian fluid in an enclosure. International Journal of Heat and Mass Transfer, vol. 46, no. 19, p. 3605-3617, DOI:10.1016/S0017-9310(03)00149-2. [14] Lamsaadi, M., Naimi, M., Hasnaoui, M. (2006). Natural convection heat transfer in shallow horizontal rectangular enclosures uniformly heated from the side and filled with non-Newtonian power law fluids. Energy Conversion and Management, vol. 47, no. 15-16, p. 2535-2551, DOI:10.1016/j.enconman.2005.10.028. [15] Lamsaadi, M., Naimi, M., Hasnaoui, M., Mamou, M. (2006). Natural convection in a tilted rectangularslot containing non-Newtonian power-law fluids and subject to a longitudinal thermal gradient. Numerical Heat Transfer, Part A: Application: An International Journal of Computation and Methodology, vol. 50, no. 6, p. 561-583, DOI:10.1080/10407780600599513. [16] Turan, O., Sachdeva, A., Chakraborty, N., Poole, R.J. (2011). Laminar natural convection of power-law fluids in a square enclosure with differentially heated side walls subjected to constant temperatures. Journal of Non-Newtonian Fluid Mechanics, vol. 166, no. 17-18, p. 1049-1063, DOI:10.1016/j. jnnfm.2011.06.003. [17] Ternik, P., Rudolf, R. (2013). Laminar natural convection of non-Newtonian nanofluid in a square enclosure whit differentially heated side walls. International Journal of Simulation Modeling, vol. 12, no. 1, p. 5-16, DOI:10.2507/ IJSIMM12(1)1.215. [18] Turan, O., Sachdeva, A., Poole, R.J., Chakraborty, N. (2013). Aspect ratio and boundary conditions effects on laminar natural convection of power-law fluids in a rectangular enclosure with differentially heated side walls. International Journal of Heat and Mass Transfer, vol. 60, p. 722-738, DOI:10.1016/j.ijheatmasstransfer.2013.01.017. [19] Ternik, P., Buchmeister, J. (2015). Buoyacy-induced flow and heat transfer of power law fluids in a side heated square cavity. International Journal of Simulation Modelling, vol. 14, no. 2, p. 238-249, DOI:10.2507/IJSIMM14(2)5.293. [20] Alloui, Z., Vasseur, P. (2015). Natural convection of Carreau-Yasuda non-Newtonian fluids in a vertical cavity heated from the sides. International Journal of Heat and Mass Transfer, vol. 84, p. 912-924, DOI:10.1016/j. ijheatmasstransfer.2015.01.092. [21] Lamsaadi, M., Naimi, M., Hasnaoui, M. (2005). Natural convection of non-Newtonian power law fluids in a shallow horizontal rectangular cavity uniformly heated from below, Heat and Mass Transfer, vol. 41, no. 3, p. 239-249, DOI:10.1007/s00231-004-0530-8. [22] Ouertatani, N., Ben Cheikh, N., Ben Beya, B., Lili, T. (2008). Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure. Comptes Rendus Mécanique, vol. 336, no. 5, p. 464-470, D0l:10.1016/j.crme.2008.02.004. [23] Khadangi Mahrood, M.R., Etemad, S. G., Bagheri, R. (2011). Free convection heat transfer of non Newtonian nanofluids under constant heat flux condition. International Communications in Heat and Mass Transfer, vol. 38, no. 10, p. 1449-1454, D0I:10.1016/j.icheatmasstransfer.2011.08.012. [24] Khezzar, L., Siginer, D., Vinogradov, I. (2012). Natural convection of power law fluids in inclined cavities. International Journal of Thermal Sciences, vol. 53, p. 8-17, DOI:10.1016/j. ijthermalsci.2011.10.020. [25] Alloui, Z., Ben Khelifa, N., Beji, H., Vasseur, P., Guizani, A. (2013).The onset of convection of power-law fluids in a shallow cavity heated from below by a constant heat flux. Journal of Non-Newtonian Fluid Mechanics, vol. 196, p. 70-82, D0I:10.1016/j.jnnfm.2013.01.008. [26] Turan, O., Lai, J., Poole, R.J., Chakraborty, N. (2013). Laminar natural convection of power-law fluids in a square enclosure submitted from below to a uniform heat flux density. Journal of Non-Newtonian Fluid Mechanics, vol. 199, p. 80-95, D0I:10.1016/j.jnnfm.2013.06.002. [27] Yigit, S., Poole, R.J., Chakraborty, N. (2015). Effects of aspect ratio on natural convection of Bingham fluids in rectangular enclosures with differentially heated horizontal walls heated from below. International Journal of Heat and Mass Transfer, vol. 80, p. 727-736, DOI:10.1016/j. ijheatmasstransfer.2014.09.046. [28] Vinogradov, I., Khezzar, L., Siginer, D. (2011). Heat transfer of non-Newtonian dilatant power law fluids in square and rectangular cavities. Journal of Applied Fluid Mechanics, vol. 4, no. 2, p. 37-42. [29] Anderson, H.I., Irgens, F. (1990). Film Flow of Power-Law Fluids, Cheremisinoff, N.P. (ed.) Encyclopedia of Fluid Mecchanics, Polymer Flow Engineering, vol. 9. Gulf Publishing Company, Houston. [30] Patankar, S. (1990). Numerical Heat Transfer and Fluid Flow. Mc Graw Hill Book Company, New York. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)10, 565-576 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3381 Original Scientific Paper Design and Evaluation of a Hierarchical Control Algorithm for an Electric Active Stabilizer Bar System Zhenxing Kong - Dawei Pi* - Xianhui Wang - Hongliang Wang - Shan Chen Nanjing University of Science & Technology, Department of Mechanical Engineering, China This paper describes a novel hierarchical control algorithm for the electric active stabilizer bar (ASB) system, which is applied to a four-wheeled road vehicle. The proposed control algorithm is designed to improve vehicle roll and yaw dynamics. The upper-level controller calculates the target active anti-roll torque via sliding mode control, which is aimed at improving the roll stability. The middle-level controller distributes active anti-roll torque between front and rear axle via fuzzy control, which can improve the handling stability through changing lateral load transfer of front and rear axles in real time. The lower level controller is employed to control the output torque of ASB actuators via PI control and improve the output characteristic of actuators with excellent response and stability. The numerical simulation and hardware-in-the-loop (HIL) experiment are carried out to evaluate the performance of the proposed control algorithm. It is demonstrated that the ASB system based on proposed control algorithm makes a significant improvement in the vehicle roll stability ride comfort and handling stability. Keywords: active stabilizer bar, hierarchical control, electric, sliding mode control, fuzzy control, numerical simulation, hardware-in-the-loop experiment Highlights • Electric actuator modelling of an active stabilizer bar system. • Enhanced roll dynamics and yaw dynamics of vehicles. • Designed sliding model controller for roll dynamics control. • Development of fuzzy controller for yaw dynamics control. 0 INTRODUCTION Vehicle roll motion can reduce the roll stability and handling stability. Improving the vehicle roll dynamics response is important. Previously, the passive stabilizer bar was installed to stabilize the vehicle body, but the performance is limited. Thus, some active roll control (ARC) systems have been developed. One of the effective solutions is active stabilizer bar (ASB) system. ASB system can generate active anti-roll torque in real time, reduce the body roll angle and roll rate, and it is suitable for applying to the mass production of cars because of the low cost and simple structure. Therefore, many researchers have been focusing on ASB systems. The HIL benches for hydraulic ASB are designed and tested in [1] to [3]. An electric active stabilizer suspension system is developed to control vehicle roll motion. The actual vehicle tests proved superior roll stability and ride comfort in [4] to [6]. A control logic for the ASB system with rotary-type hydraulic stabilizer actuators is proposed in [7]. The control logic consists of a feedforward controller and a feedback controller. Through the test, the ASB system demonstrated the successful reduction of the roll motion under all conditions. Moreover, [8] puts forward a linear quadratic (LQ) controller based on two-level architecture. There has also been some research on coordination control of ASB and other chassis control systems. Yim et al. [9] present a method for designing a controller that uses ASB and ESP for rollover prevention. Moreover, an integrated control strategy of the differential braking, the semi-active suspension damper and the active roll moment is analysed in [10]. In this paper, a novel hierarchical control scheme of ASB is put forward. The upper-level controller calculates active anti-roll torque. The middle-level controller distributes active anti-roll torque between front and rear active stabilizer bars. The lower level controller is employed to control the output torque of actuators. The structure of this paper is as follows. The system dynamic models including vehicle and tire model, road input model and ASB actuator model are presented in Section 1. The hierarchical controller of the ASB system is designed in Section 2. In Section 3, numerical simulations and HIL experiments are implemented, together with results. Finally, the conclusions of the paper are summarized in Section 4. 1 SYSTEM DYNAMIC MODELS 1.1 Vehicle and Tire Model A 14 degree-of-freedom (DOF) vehicle dynamics model is adopted in this paper [11], which includes the longitudinal, lateral, roll, yaw, pitch, vertical and Corresponding Author's Address: Nanjing University of Science & Technology, Department of Mechanical Engineering, Nanjing, Jiangsu 210094, China, pidawei@mail.njust.edu.cn, 13770669850@139.com rotational motion of four wheels, vertical motion of the vehicle body. The top, front and left view of a vehicle system model is presented in Fig. 1. Moreover, the Dugoff tire model [12] and [13] is adopted, which can reflect the longitudinal tire force and lateral force variation with the slip ratio, the side slip angle, and the tire vertical load. "m J df f F rwyfl fr v F x 1 wxfl F. wyfr F •yrld. ф в M Hw zqii Fig. 1. Vehicle model 1.2 Road Input Model The road roughness plays a major role in the motion of the vehicle. Considering the road roughness of four wheels, a C-level road model is established based on the filtered white noise method in MATLAB/Simulink [14]. 1.3 ASB Actuator Model The electric ASB actuator consists of left-half and right-half stabilizer bars, a brushless direct current (BLDC) motor, a reduction gear, a housing and two bushings [6]. The force diagram of an active stabilizer bar is presented in Fig. 2. F, M. MActuator b ^—D-EH3—. ^ 1=1 ( Y. <-M—+i Jr--^. p t/2 1 * Fig. 2. Force diagram of an electric active stabilizer bar The relationship between MActuator and MARC is denoted by Eqs. (1) and (2). F _ MARC _ MActuator MARC , MActuator. b (1) (2) Therefore, the output torque of ASB actuators can be calculated by Eq. (3) and (4). Marc , f _ ~M Actuator, f > MARC,r j MActuator br (3) (4) where MActuatori is the output torque of ith ASB actuator, MARC i is the active anti-roll torque of ith ASB system; t is the length of the stabilizer bar, and b is length of the lever. To calculate and debug easily, the BLDC motor model is replaced by the direction current (DC) motor model [15] by Eqs. (5) to (8). U - Ue =(L - M) + IR, U = K co, e e Te - T = Jdw+ Bw, e ' dt t _ KI. (5) (6) (7) (8) Considering the reduction gear, the output torque of front and rear ASB actuators are given by Eqs. (9) and (10). MActuator,f f , MActuator ,r (9) (10) where Ke is back-emf constant, Kt torque constant, i transmission ratio of the reduction gear, I motor current, B motor damping coefficient, J motor inertia, M electromotive damping, R armature resistance, Te electromagnetic torque, Tl load torque, U supply voltage of the armature, Ue back electromotive force (emf), and m motor angular velocity. 2 HIERARCHICAL CONTROLLER DESIGN Based on the hierarchical structure, a design solution for the ASB system is proposed. The block diagram of the control system is presented in Fig. 3. The b t S F r y F wxrl wxri hierarchical control system includes three level controllers: 1. Upper-level controller: to improve the roll stability of vehicles, the upper-level controller is designed to generate the target active anti-roll torque. 2. Middle-level controller: the middle-level controller is designed to distribute the active anti-roll torque between front and rear active stabilizer bars, which can improve yaw stability. 3. Lower level controller: to improve the output performance of the ASB actuator, the lower level controllers are designed to control the output torque of the ASB actuator. Control Algorithm M. Lower Level 1 Front Controller 1 1 ASB Upper Level Controller MARC -► Middl Cont e Level roller M ARC, * Lower Level 1 Rear ASB Controller 2 1 5f, v,, ay r 14 DOF Vehicle Model Steering Input q Road Disturbances Fig. 3. Block diagram of the control system 0), (14) then = ce + c2e - f (ф,ф, t )-Д/ (ф,ф, t )- - b (t )MÄRC - d (t )+ф. (15) The uncertainty of the system is given by Eq. (16). 1-Л/(ф,ф,t)|< F, \-d(t) + ф| < D. (16) + Ö Choose the reaching condition which guarantees the asymptotic rapidity and stability as follows Eq. (17): s = -s sgn(s) - Ks, (s,K > О). The gains are determined as (Eq. (18)): s = F + D + n,K > 0. The control law is designed as follows: Ce + c2è - J (ф,ф, t) +s sgn(s) + Ks MARC = b (t ) (17) (18) (19) Consider a Lyapunov function candidate in Eq. (20): V 1 2 V =— s . 2 (20) The sliding condition is satisfied according to Eq. (21): V = ss = s[cle + c2è - (7(ф,ф, t) + Л/(ф,ф, t)) - -b(t)marc - d7) + ф] = = s [-Л/ (ф,ф, t ) - d (t ) + ф -e sgn(s) - Ks~j < < s [F + D - e sgn(s) - Ks] = = (F + D) s-e s - Ks2 < < (F + D) s-e s - Ks2 =-v\s\-Ks2 < 0. (21) To eliminate the high-frequency chatter further, the sign function is replaced by a saturation function in Eq. (22): sat I ^ I = s Ф sgn(s) < 1 , > 1 (22) where Ф is boundary-layer thickness. The output of sliding mode controller is presented by Eq. (23): MARC = 1 b (t ) c1e + c2e - / (ф,ф, t ) + ) + Ks , (23) that is MARC = (msgK - Кф)Ф- СфФ + mshray - 2 S -(Ix + msh2 ) cxe + c2e + £sat(—) + Ks . (24) Thus the control law obtained above can meet the requirement of system stability and track the target value accurately. All control parameters are listed in Table 1. Table 1. Controller parameters c1 = 140 K = -- 10 Ф = 0.1 С2 = 10 e = 0.1 2.2 Middle Level Controller The middle-level controller is designed to distribute the anti-roll torque between front and rear active stabilizer bars. Moreover, the good coordination of roll stability and yaw stability can be guaranteed [18]. The diagram is shown in Fig. 6. Fig. 6. Diagram of middle-level controller Based on two DOF linear vehicle model [9], the steady-state value of target yaw rate rt ss is achieved as follows (Eq. (25)): ' ' * (25) 1 + - ■■2 ( I 12 j. C S5/ where vx is longitudinal velocity, Sf front wheel angle, Cf / Cr cornering stiffness of front/ rear axle, l wheel base, and f / lr front/ rear share of wheel base. Moreover, the yaw rate transient response is characterized in Eq. (26). 1 0.1s +1 (26) Considering limitation of the road friction, the limited value of yaw rate [17] is given as: 0.85^ (27) where ц is road friction coefficient. Thus, the target yaw rate is computed as follows (Eq. (28)): r = sgn (rti )•min (Ы> Ir 21). (28) The inputs of the fuzzy controller are yaw rate r and yaw rate error dr. The output of the fuzzy r V x controller is the distribution coefficient Л ( X= MARC, f IMARC e [0,1] ). In the fuzzification step, linguistic variables are used to make the input variables r and dr and the output variable Л compatible with the condition of the knowledge-based rules. The linguistic terms are shown in Table 2. The fuzzy sets used for inputs and outputs are defined as: {r}={N, ZE, P} {dr}={NB, NS, ZE, PS, PB} U}={ZE, S, M, B, L}. Table 2. Linguistic terms for middle-level controller N Negative ZE Zero P Positive NB Negative big NS Negative small PS Positive small PB Positive big S Small M Middle B Big L Large During the fuzzy decision process, a list of fuzzy rules is defined on the basis of expert knowledge and extent simulation. The Mamdani fuzzy inference system is adopted, and the weights of the rules are considered to be 1. The fuzzy rules are shown in Table 3. These rules are introduced on the basis of the following criteria. Table 3. Fuzzy rules for middle level controller r dr Л r dr Л N NB ZE ZE PS M N NS S ZE PB M N ZE M P NB L N PS B P NS B N PB L P ZE M ZE NB M P PS S ZE NS M P PB ZE ZE ZE M Criterion 1 : if r is NB and dr is N, then Л is ZE. In this case, the vehicle is in a serious understeer state. The lateral load transfer of front axle is set to increase, and the lateral load transfer of rear axle is set to decrease. Thus, the distribution coefficient tends to zero. Criterion 2: if r is PB and dr is N, then Л is L. In this case, the vehicle is in a serious oversteer state. The lateral load transfer of front axle is set to decrease, and the lateral load transfer of rear axle is set to increase. Thus, the distribution coefficient tends to be large. Criterion 3: if r is ZE or dr is ZE, then Л is M. In this case, the vehicle is in a stable state. The lateral load transfer of the front axle is set to be a little smaller than the lateral load transfer of the rear axle, which can make the vehicle be in a slight understeer state. Thus, the distribution coefficient is set about 0.55. The centre-of-area defuzzification method is used to scale and map the fuzzy output to produce an output value for the control system. Through many simulations and analyses, r is set within the range of [-1 rad, 1 rad], dr is set within the range of [-0.25 rad/s, 0.25 rad/s], and Л is set within the range of [0, 1]. The membership functions of input and output variables are shown in Fig. 7. & j3 2 0.8 и л §0.6 S ° 0.4 2 ^0.2 Q 0 N ZE P -1 a) -0.75 -0.5 -0.25 0 0.25 0.5 0.75 r [rad/s] 0.8 ■o S 0 Л1 о 0 (D ад 0 и v Q b) NB NSZ EPS PB i -0.2 -0.1 0 0.1 0.2 dr [rad/s] ■о s О Л1 v О 0 (D &0 Q Fig. 0 0.2 0.4 0.6 0.8 1 I 7. Membership function of variables; a) r, b) dr and c) A 1 1 1 0 c) Therefore, active anti-roll torques of the front and rear active stabilizer bars are given by Eq. (29) and (30). MARC, f = ШЛКС , M. = (1 -b)MA (29) (30) 2.3 Lower Level Controller The lower level controller is designed with PI control theory to control the output torques of the front and rear active stabilizer bars, which can improve the actuator output performance with excellent response and stability. In practical engineering application, the control of actuator output torque can be realized by controlling the electric current of motors [16]. The control diagram of the system is shown in Fig. 8. If PI а Drive DC Motor Controller Circuit Model T Fig. 8. Diagram of lower level controller Considering the ASB actuator model in section 1.3, the target electric current of motors can be presented in Eq. (31) and (32). I, / = Kt/ -M ARC, /' I.. =- iK,t -M, (31) (32) The input of PI controller is the current error e, (e, = It -I), and the output of PI controller is PWM duty ratio a of the motor current. The control algorithm is given by Eq. (33). The controller parameters are listed in Table 4. а = kpet + ki J etdt. (33) Table 4. Controller parameters kp = 0.1 h = 40 manoeuvres. The parameters of vehicle and ASB actuator adopted are listed in Tables 5 and 6. Table 5. Parameters of the vehicle m 1704.7 kg hr 0.445 m ms 1526.9 kg 53015 Nm/rad Ix 744 kg-m2 3534 Nm/(rad/s) Iz 3048.1 kg-m2 Rw 0.313 m If 1.035 m Iw 0.99 kg-m2 Ir 1.655 m Cf 66000 N/rad df 1.535 m Cr 70000 N/rad dr 1.535 m Table 6. The output characteristic of the actuator Actuator Maximum torque Maximum torque rate ASB actuator 700 Nm 1600 Nm/s The vehicle is set to travel on the C-level road with an initial speed of 80 km/h, no braking, and accelerating. The road friction coefficient ц is 0.8. The C-level road height is presented in Fig. 9. Simulations are conducted with different types of suspensions, including 'without ASB', 'ASB without dynamic distribution' (distribution coefficient X is 0.55) and 'ASB with dynamic distribution' respectively, under J-turn, sine wave and zero input manoeuvres. The front wheel angle under different manoeuvres above is shown in Fig. 10. 0.02 0.01 .c M I ■S 0lfl •c -0.01 -0.02 -Front left wheel Front right wheel Rear left wheel -----Rear right wheel 0 12 3 4 10 Fig. 9. The C-level road height of four wheels 3 SIMULATION AND EXPERIMENT 3.1 Numerical Simulation The numerical simulation is carried out to evaluate the performance of the ASB system under typical 3.1.1 Roll Stability The simulation results of roll angle and roll rate under J-turn and sine wave manoeuvres are shown in Figs. 11 and 12. Compared with 'without ASB' case, roll angle and roll rate of 'ASB with dynamic distribution' b r 6 0.2 2 0.15 0.1 e 0.05 a) 0123456789 10 Time [s] 0.2 T3 2 0.1 £ -0.1 Ич -0.2 b) 0123456789 10 Time [s] Fig. 10. Front wheel angle under different manoeuvres; a) J-turn and b) sine wave 4 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 —— Target value ■ ■ ■ - Without ASB ASB with dynamic distribution 0 1 a) 23456 Time [s] 7 8 9 10 10 -Without ASB "ASB with dynamic distribution b) 0123456789 10 Time [s] Fig. 11. Roll motion under J-turn manoeuvre; a) roll angle and b) roll rate case is reduced by 46.1 % and 45.2 % in average. The roll angle is very close to target roll angle, which guarantees a rapid response and small static error. Under sine wave manoeuvre, the roll angle follows the target value rapidly, even though the steering wheel input frequency reaches 0.7 Hz as shown in Fig. 12a. It is confirmed that the roll angle control of ASB system shows excellent frequency response characteristics. — Target value ■ ■ - Without ASB ■ _ ASB with dynamic distribution a) 0123456789 10 Time [s] Without ASB "ASB with dynamic distribution M ü U 0 о Pi ■ -10 -15 b) 0123456789 10 Time [s] Fig. 12. Roll motion under sine wave manoeuvre; a) roll angle and b) roll rate Figs. 13a and b show vehicle roll motion under zero input manoeuvre. The roll angle and roll rate caused by the road input excitation both decreases. Considering the human body sensitivity to roll vibration (around 0.8 Hz) and roll resonant frequency (around 2 Hz), the PSD of roll rate in the 0.3 Hz to 3 Hz domain is used to evaluate the vehicle ride comfort [5]. From Fig. 13c, the PSD of roll rate in the frequency range of 0.3 Hz to 3 Hz is reduced by 5 dB to 20 dB. Obviously, the vehicle with the ASB system has great ride comfort performance. 0 5 0 ^ 0.1 м и •o 0 "ад § -0.1 -0.2 -0.3 -0.4 a) -Witout ASB ASB with dynamic distribution 0 1 2 3 4 5 6 7 8 9 10 Time [s] b) -Without ASB ASB with dynamic distribution "0 123456789 10 Time [s] 20 m T3 _ -20 IS 2 -40 О Q £ -60 Without ASB ASB with dynamic distribution c) 10 10' Frenquency [Hz] Fig. 13. Roll motion under zero Input manoeuvre; a) roll angle, b) roll rate and c) PSD of roll rate 3.1.2 Handling Stability The simulation results of yaw rate are shown in Fig. 14. Compared with the 'without ASB' case, the yaw rate is closer to the target value in the 'ASB without dynamic distribution' and 'ASB with dynamic distribution' cases. Moreover, the yaw rate peak values of two cases are both reduced. Under the J-turn manoeuvre, the yaw rate of 'ASB with dynamic distribution' case has better transient and steady-state response characteristics in comparison to 'ASB without dynamic distribution' case, which indicates smaller overshot and oscillation of yaw rate. Moreover, yaw rate of 'ASB with dynamic distribution' case is closer to the target value as shown in Fig. 14a. 30 25 20 И 1 5 и и IL £ 10 л iS 5 >- a) 0 -5 -10 0 60 г —Target value -Without ASB ■ - ASB without dynamic distribution ASB with dynamic distribution 1 2 3 4 5 6 Time [s] 7 8 9 10 — Target value Without ASB ASB without dynamic distribution ASB with dynamic distribution _ b) 4 5 6 Time [s] Fig. 14. Yaw rate under different manoeuvres; a) J-turn and b) sine wave However, in Fig. 14b, the 'ASB with dynamic distribution' case has no significant effect on the yaw rate control in comparison to the 'ASB without dynamic distribution' case. Two important reasons are supposed as follows: 1. Under extreme steering manoeuvres, the vehicle is in a serious unstable state. Considering the limitation of the coupling relationship between roll and yaw dynamics, it is very difficult to implement the vehicle yaw rate control. 2. The higher response speed of the actuator is put forward for the high steering input frequency. The performance of the yaw rate control decreased in 0 the high-frequency region due to the limitation of the ASB actuator's output characteristic, such as the torque rate. In general, the ASB system with dynamic distribution provides a more rapid response and smaller oscillation of the yaw rate, which considerably improves vehicle handling stability. 3.2 Hardware-in-the-loop Experiment The experiments are performed based on an HIL system. The HIL system includes the real-time simulation system (AutoBox), ASB electric control unit (ASB ECU), ASB actuators and force sensors. Fig. 15 shows the block diagram of the HIL system. The AutoBox runs a vehicle dynamics model, a Dugoff tire model, and a road input model, and exports vehicle state signals. According to those signals, ASB control algorithm in ASB ECU calculates the target anti-roll torque of front and rear ASB actuator. ASB actuators generate actual anti-roll torque, which is measured by force sensors. The actual anti-roll torque of front and rear ASB actuator are imported to the vehicle dynamics model. 4r 3.53 2.521.5 1 0.50 0.5 a) -Target value ■ ■ ■ -Without ASB --ASB with dynamic distribution 0123456789 10 Time [s] —Without ASB ■_ ASB with dynamic distribution b) 0123456789 10 Time [s] Fig. 16. Roll motion under J-turn manoeuvre; a) roll angle and b) roll rate ASB ECU (ASB Control Algorithm) 7\ Vehicle State AutoBox Target Anti-roll Torque_front Vehicle Dynamics Model /L-N W Dugoff Tire Model Road Input Model Target Anti-roll Torque_rear Actual Anti-roll Torque_front Actual Anti-roll Torque_rear Fig. 15. Block diagram of HIL system The initial setup of the experimental vehicle is the same as that in numerical simulation (Section 3.1). The HIL experiments are conducted under J-turn, sine wave, and zero input manoeuvres. 3.2.1 Roll Stability Figs. 16 and 17 show the vehicle roll motion under J-turn and sine wave manoeuvre. The 'ASB with dynamic distribution' case reduces the vehicle roll angle and roll rate by 45.5% and 43.8%, respectively, W r\ * * Ìi A Mr л 0 1 2 3 4 5 6 7 8 9 10 Time [s] -Without ASB ASB with dynamic distribution b) 0123456789 10 Time [s] Fig. 17. Roll motion under sine wave manoeuvre; a) roll angle and b) roll rate 4 0 in average. Moreover, the roll angle follows the target value quickly with small chatter caused by the physical limitation of ASB actuators output characteristic. The steady-state error of roll angle is less than 0.2 deg. From Fig. 17a, the amplitude of roll angle is just about 0.15 deg bigger than the target value during the entire sine wave maneuver, and the phase of roll angle is the same as a target value, which means good frequency response characteristic of ASB system. It is verified that the proposed ASB control algorithm can improve the vehicle roll stability. a) -Without ASB ASB with dynamic distribution 1 3 4 5 6 7 8 9 10 Time [s] M (D ТЗ b) 1 0 -1 -24 20 0 -20 -40 -60 -80 —Without ASB ASB with dynamic distribution 123456789 10 Time [s] -Without ASB --ASB with dynamic distribution c) 10 10 Frequency [Hz] Fig. 18 shows vehicle roll response under zero input manoeuvre. From Fig. 18a and 18b, compared with 'Without ASB' case, 'ASB with dynamic distribution' reduces the roll angle and roll rate, which improves vehicle roll stability. From Fig. 18c, the PSD of roll rate decreases by 0 dB to -18 dB in the frequency range of 0.3 Hz to 3 Hz, which means better vehicle ride comfort. 3.2.2 Handling Stability Fig. 19 shows the vehicle yaw rate under J-turn and sine wave manoeuvre. 'Without ASB' and 'ASB without dynamic distribution' both have poor performance in vehicle yaw stability. From Fig. 19a, during 1 s to 2.5 s, the vehicle is in an oversteer state, and the anti-roll torque distribution reduces the overshot of yaw rate by 21.1%. During 2.5 s to 10 s, the anti-roll torque distribution reduces the oscillation of yaw rate, makes yaw rate follow the target value, and the steady-state error of yaw rate is less than 3 deg/s. 30 252015 10 50 -5 M "Й I >H -10 L — Target value Without ASB ----ASB without dynamic distribution ■ ■ ASB with dynamic distribution 0 1 2 3 a) 4 5 6 Time [s] 7 8 9 10 60 40 f20 TD rt SI 0 >н -20 -40 -Target value -Without ASB ASB without dynamic distribution ■ ■ ASB with dynamic distribution b) 456 Time [s] Fig. 18. Roll motion under zero Input manoeuvre; a) roll angle, b) roll rate and c) PSD of roll rate Fig. 19. Yaw rate under different manoeuvres; a) J-turn and b) sine wave 0 2 Thus, the vehicle is in a stable state. From Fig. 19b, since the limitation coupling relationship between the vehicle roll and yaw dynamics and physical limitation of ASB actuators output performance, the anti-roll torque distribution does not achieve desired control effect, which just slightly reduces the peak value of yaw rate by about 1deg/s in average. Above all, ASB system with dynamic distribution makes the vehicle yaw rate within a stable region and improves the vehicle handling stability under some extreme conditions. 4 CONCLUSIONS A novel hierarchical control algorithm for an electric ASB system is proposed. Numerical simulations and HIL experiments are implemented to study the roll and handling stability of the vehicle with the proposed ASB control system under typical manoeuvres. Conclusions can be drawn as follows. 1. Three level controllers implemented for ASB system are devised based on the hierarchical control structure. In the upper level, a sliding mode controller is designed to generate active anti-roll torque, which improves vehicle roll stability. In the middle level, a fuzzy controller is proposed to distribute active torque between front and rear active stabilizer bars, which improves vehicle handling stability. In the lower level, a PI controller is applied to generate the actuator output torque. 2. Simulation and experimental results confirmed the reliability and precision of the ASB hierarchical control algorithm. The proposed control algorithm can reduce vehicle roll angle and roll rate effectively and improve vehicle yaw dynamics response. In summary, the good coordination of the vehicle roll stability, ride comfort and handling stability can be achieved by the proposed ASB control algorithm. 5 ACKNOWLEDGEMENTS This work was supported by the National Science Foundation (Grant No. 51205204, No. 51205209), the Fundamental Research Funds for the Central Universities (Grant No. 30915118832), and the Six Talent Program of Jiangsu Province (Grant No. 2014003). 6 NOMENCLATURES ay lateral acceleration of the vehicle [m/s2] Cf cornering stiffness of front axle [N/rad] Cr cornering stiffness of rear axle [N/rad] Сф total suspension roll damping [N/(rad/s)] df front tread [m] dr rear tread [m] FwXa/FMy^ii/FM^zii longitudinal, lateral and vertical force of iith wheel [N] g acceleration of gravity [m/s2] hp distance between pitch axis and centre of gravity [m] hr distance between roll axis and centre of gravity [m] i transmission ratio of reduction gear [-] I motor current [A] Iw moment of inertia of wheel [kgm2] IX roll moment of inertia [kgm2] Iz yaw moment of inertia [kgm2] J motor inertia [kgm2] Ke back-emf constant [V] Kt torque constant [Nm/A] Кф total suspension roll stiffness [Nm/rad] lf front share of wheel base [m] lr rear share of wheel base [m] m total mass [kg] ms sprung mass [kg] M electromotive damping [H] MActuator f output torque of front ASB actuator [Nm] MActuator r output torque of rear ASB actuator [Nm] MARC,f active anti-roll torque of front ASB system [Nm] MARC, r active anti-roll torque of rear ASB system [Nm] r yaw rate [rad/s] rt target yaw rate[rad/s] R armature resistance [ß] Rw rolling radius of wheel [m] Te electromagnetic torque [Nm] Tl load torque [Nm] U supply voltage of the armature [V] Ue back electromotive force (emf) [V] vX longitudinal velocity of the vehicle [m/s] vy lateral velocity of the vehicle [m/s] Zb displacement of vehicle body [m] Zqii road height of iith wheel [m] Zwii displacement of iith wheel [m] df steering angle of the front wheel [rad] H road friction coefficient[-] ф roll angle [rad] ф( target roll angle [rad] m motor angular velocity [rad/s] 7 REFERENCES [1] Sorniottia, A., Morgandoa, A., Velardocchiaa, M. (2006). Active roll control: system design and hardware-in-the-loop test bench. Vehicle System Dynamics: International Journal of Vehicle Mechanics and Mobility, vol. 44, suppl. 1, p. 489-505, D0l:10.1080/00423110600874677. [2] Cimba, D., Wagner, J., Baviskar, A. (2006). Investigation of active torsion bar actuator configurations to reduce vehicle body roll. Vehicle System Dynamics: International Journal of Vehicle Mechanics and Mobility, vol. 44, no. 9, p. 719-736, D0l:10.1080/00423110600618280. [3] Sorniotti, A., D'Alfio, N. (2007). Vehicle dynamics simulation to develop an active roll control system. SAE Technical Paper, no. 2007-01-0828, D0I:10.4271/2007-01-0828. [4] Buma, S., Ookuma, Y., Taneda, A., Suzuki, K., Cho, J.-S., Kobayashi, M. (2010). Design and development of electric active stabilizer suspension system. Journal of System Design & Dynamics, vol. 4, no. 1, p. 61-76, D0I:10.1299/jsdd.4.61. [5] Mizuta, Y., Suzumura, M., Matsumoto, S. (2010). Ride comfort enhancement and energy efficiency using electric active stabiliser system. Vehicle System Dynamics: International Journal of Vehicle Mechanics & Mobility, vol. 48, no. 11, p. 1305-1323, D0I:10.1080/00423110903456909. [6] Jeon, K., Hwang, H., Choi, S., Kim, J., Jang, K., Yi, K. (2012). Development of an electric active roll control (ARC) algorithm for a SUV. International Journal of Automotive Technology, vol. 13, no. 2, p. 247-253, D0I:10.1007/s12239-012-0021-8. [7] Kim, S., Park, K., Song, H.J., Hwang, Y.K., Moon, S.J., Ahn, H.S., Tomizuka, M. (2012). Development of control logic for hydraulic active roll control system. International Journal of Automotive Technology, vol. 13, no. 1, p. 87-95, DOI:10.1007/ s12239-012-0008-5. [8] Varga, B., Németh, B., Gaspar, P. (2015). Design of anti-roll bar systems based on hierarchical control. Strojniški vestnik -Journal of Mechanical Engineering, vol. 61, no. 6, p. 374-382, D0I:10.5545/sv-jme.2014.2224. [9] Yim, S., Jeon, K., Yi, K. (2012). An investigation into vehicle rollover prevention by coordinated control of active anti-roll bar and electronic stability program. International Journal of Control, Automation, Systems, vol. 10, no. 2, p. 275-287, D0I:10.1007/s12555-012-0208-9. [10] Her, H., Suh, J., Yi, K. (2015). Integrated control of the differential braking, the suspension damping force and the active roll moment for improvement in the agility and the stability. Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, vol. 229, no. 9, p. 1154-1154, D0I:10.1177/0954407014550502. [11] Ghike, C., Shim, T. (2006). 14 degree-of-freedom vehicle model for roll dynamics study. SAE Technical Paper, no. 200601-1277, D0I:10.4271/2006-01-1277. [12] Dugoff, H., Fancher, P.S., Segel, L. (1970). An analysis of tire traction properties and their influence on vehicle dynamic performance. SAE Technical Paper, no. 700377, D0I:10.4271/700377. [13] Guntur, R., Sankar, S. (1980). A friction circle for Dugoff's tyre friction model. International Journal of Vehicle Design, vol. 1, no. 4, p. 373-377, D0I:10.1504/IJVD.1980.061234. [14] Zhang, L.J., Zhang, T.X. (2008). Study on general model of random inputs of the vehicle with four wheels correlated in time domain. Transactions of the Chinese Society of Agricultural Machinery, vol. 27, no. 7, p. 75-78, D0I:10.3969/j. issn.1000-1298.2005.12.008. [15] Zhang, H., Liu, J., Ren, J., Zhang, Y., Gao, Y. (2009). Research on electric power steering with BLDC drive system. IEEE 6th International IEEE Power Electronics and Motion Control Conference, p. 1065-1069, D0I:10.1109/ IPEMC.2009.5157543. [16] Kim, J.-H., Song, J.B. (2002). Control logic for an electric power steering system using assist motor. Mechatronics, vol. 12, no. 3, p. 447-459, D0I:10.1016/S0957-4158(01)00004-6. [17] Pi, D.-W., Chen, N., Zhang, B.-J. (2011). Experimental demonstration of a vehicle stability control system in a split-M manoeuvre. Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, vol. 225, no. 3, p. 305-317, D0I:10.1177/09544070JAUT01541. [18] Rajamani, R. (2005). Vehicle Dynamics and Control. Springer, New York. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)10, 577-590 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3545 Original Scientific Paper Multi-objective Optimization of Cloud Manufacturing Service Composition with Cloud-Entropy Enhanced Genetic Algorithm Yongxiang Li - Xifan Yao* - Jifeng Zhou School of Mechanical & Automotive Engineering, South China University of Technology, China To consider the service-matching degree, the composition harmony degree, and the service composition complexity in cloud manufacturing service composition optimization problems, a new composition optimization approach, called cloud-entropy enhanced genetic algorithm (CEGA), is put forward to solve such problems with multi-objectives. The definitions of service-matching degree, composition harmony degree, and cloud-entropy and the corresponding calculation methods are given. A multi-objective optimization mathematical model of cloud manufacturing service composition is built. The manufacturing task of AGV (automated guided vehicle) is taken as an example to verify the proposed CEGA algorithm on the established composition model. The studied result shows that CEGA converges faster than a standard genetic algorithm with shorter time. Keywords: cloud manufacturing, service composition optimization, cloud-entropy, service-matching degree, composition harmony degree, genetic algorithm Highlights • Cloud-entropy enhanced genetic algorithm (CEGA) is proposed. • Service-matching degree, composition harmony degree, and cloud-entropy are defined, and the corresponding calculation methods are given. • The multi-objective optimization mathematical model of cloud manufacturing service composition is built. • The manufacturing task of AGV is taken as an example to verify the proposed CEGA algorithm on its feasibility and effectiveness. • The proposed CEGA converges faster than a standard genetic algorithm with shorter time. 0 INTRODUCTION With the rise of cloud computing [1] and the development of Internet of things [2] and other related technologies, the cloud manufacturing based on them provides a new way to realize centralized management, unified scheduling, on-demand use of all kinds of distributed manufacturing resources [3]. Cloud manufacturing, in the form of cloud services, for the whole manufacturing life cycle, provides all kinds of manufacturing resources which can obtain and pay for use at any time [4]. A single cloud service is often difficult to meet customer's manufacturing requirements. Therefore, with service composition, a cloud manufacturing service platform assembles a number of fine-grained simple cloud services into coarse-grained complex cloud services to meet customer's complex manufacturing requirements [5]. Implementation of cloud manufacturing is a collaborative process that requires the participation of multiple distributed manufacturing resources. In the cloud manufacturing service platform, it is shown as a series of activities, such as manufacturing demand decomposition, service matching with tasks, service composition optimization, service execution monitoring, etc. [6]. According to the opening degree of services to users, cloud manufacturing services can be divided into public cloud services and private cloud services [7]. The manufacturing resources of private cloud services are limited to internal use, not for providing services directly to users outside the enterprise. The manufacturing resources of public cloud services are not only available for internal users but also directly provide services to users outside the enterprise. In manufacturing processes, there are a variety of information interaction and material transportation among cloud manufacturing services. It is a kind of complex relationship that is full of mutual dependence, constraints and competition [8]. Therefore, assembling appropriate fine-grained cloud services according to the customer's manufacturing needs is an important part of the implementation process of cloud manufacturing. Some scholars have done research on the optimization of service composition. For example, Lartigau et al. [9] proposed an improved artificial bee colony optimization algorithm for solving the cloud manufacturing service composition model based on QoS with geo-perspective transportation; Jeong and Lee [10] proposed a web service composition formal model with business logic process buffer based on CSP (Communication Sequential Process); Castejon et al. [11] proposed multi-objective optimal design procedure applied to a robotic arm for service tasks; Jovanovic et al. [12] described production cycle scheduling algorithm on the grounds of investigations of manufacturing capacity utilization levels and causes of loss for special-purpose products in complex manufacturing environments; Gaaloul et al. [13] proposed a dynamic mining algorithm based on a statistical technique to discover composite web service patterns from execution logs; Stegaru and Stanescu [14] presented a methodological modeling framework for quality driven web service composition that tackles quality from the user, final product, and process views; Omid et al. [15] proposed a framework for context-aware web service composition using planning techniques; Iordache and Moldoveanu [16] introduced a QoS-aware end-to-end web service composition approach that handles all the stages from the web service discovery step to the actual binding of the services, using a genetic algorithm to compute and compare the aggregated QoS of the composite services; Berlec et al. [17] proposed basic and extended models that take into account the tied-up capital in a production to calculate the optimal batch quantity; Avitabile et al. [18] used component modes from unconnected components as projection matrices to identify the system level full field response; Florjanic et al. [19] proposed an artificial neural network model to address the problem of estimating the volume of manufacturing hours in model manufacturing; Suzic et al. [20] studied mass customized production by using group technology and production flow analysis in a panel furniture manufacturing company; Chifu et al. [21] presented an ant-inspired method for selecting the optimal solution in semantic web service composition; Xiao et al. [22] proposed an improved shuffled frog leaping algorithm for solving the optimization problems of multi-objective production transportation scheduling in a cloud manufacturing environment; Dong and Guo [23] proposed a cloud manufacturing service evaluation and selection optimization model based on a global trust degree and composite template; Jing et al. [24] presented a cloud manufacturing service composition optimization algorithm based on discrete particle swarm optimization and execution reliability; Zhao et al. [25] proposed the service capability evaluation model of small and medium enterprise design knowledge resources and the mathematical model of design knowledge resource serialization combination under a cloud manufacturing environment, and a quantum harmony search algorithm for solving the model was given; Bao et al. [26] divided the service composition into three stages: the matching of adjacent nodes, the cleaning of broken branch nodes and the combination of atomic services, and put forward three steps service composition algorithm of the combination of atomic services based on knowledge. The above documents built service composition models and the corresponding solution algorithms mainly based on the analysis of service composition targets, most of which idealized a cloud manufacturing service as a kind of "rigid body", which was separated from the social relationship, seldom considering the relationships between the cloud manufacturing services and manufacturing tasks, and between cloud manufacturing services on the service composition optimization. However, to meet the personalized needs in the cloud manufacturing environment, it requires collaborative participation among customers and cloud manufacturing service providers and others [27]. The matching degree of cloud manufacturing services and manufacturing tasks, the harmony degree of cloud manufacturing services in service composition, and the complexity of service composition have great influence on the completion of customized/personalized product manufacturing tasks. Therefore, it is imperative to consider those factors for studying the optimization of cloud manufacturing service composition. This study constructs a cloud manufacturing service composition optimization mathematical model with the consideration of the matching degree of service and task in the process of cloud manufacturing service composition, the composition harmony degree, and the complexity of service composition, and then proposes a cloud-entropy enhanced genetic algorithm (CEGA) to solve the problem. Finally, the model and the algorithm are verified by an example of an automatic guided vehicle manufacturing task. The remainder of this paper is organized as follows: Section 1 analyzes the objectives and the strategies of cloud manufacturing service composition; Section 2 gives the definitions of service-matching degree, composition harmony degree, and cloud-entropy, as well as the corresponding calculation methods; Section 3 proposes CEGA algorithm; Section 4 analyzes the application case; and Section 5 concludes the whole paper and gives future work. 1 CLOUD MANUFACTURING SERVICE COMPOSITION ANALYSIS In a cloud manufacturing environment, customers can organize multiple distributed manufacturing resources to complete manufacturing tasks using cloud manufacturing service platform service composition modules. The goal of cloud manufacturing service composition is that various product manufacturing tasks are allocated to private cloud services and the most appropriate public cloud services, making the service composition meet the technical requirements of product, process time constraints and other conditions, and to obtain the maximum economic benefit [28]. The matching degree of cloud manufacturing services and manufacturing tasks, the harmony degree of cloud manufacturing services in service composition, and the complexity of service composition have significant impacts on the service composition, so cloud service composition should not only consider the product delivery period, manufacturing cost and other traditional factors, but also address service-matching degree, composition harmony degree and cloud-entropy. In the cloud manufacturing service composition, matching manufacturing tasks with public cloud services and/or enterprise private cloud services is necessary. If the cloud manufacturing service provider cannot understand its customers' manufacturing needs completely, or in the cloud manufacturing service platform it is hard to accurately express customers' requirement, customers and cloud manufacturing service providers need to communicate, coordinate, and collaboratively participate in product manufacturing processes. If the cloud service customer belongs to a branch of a group company, due to the restriction of its equipment manufacturing capacity and technical level, it is difficult for the customer to complete a manufacturing task alone, so it needs to be completed with the enterprise private cloud services or public cloud services. The main effect of customers' participation in manufacturing processes is the verification and feedback of product manufacturing requirements. Such a manufacturing model in which customers are directly involved in the product manufacturing processes is different than traditional closed-door manufacturing and networked manufacturing based on the knowledge of the customers' description of manufacturing requirements. Its outstanding advantage is that customers can directly inspect and provide feedback on a manufacturing process or even some parts and components in the process of product development or manufacturing, and without any need to wait for all the manufacturing tasks to be completed before the product information feedback. To a certain extent, it can reduce a product's defect rate and ensure that the product meets customer's customized manufacturing needs with high standards. Moreover, it also shortens the time of trial production and the rework time of unqualified products, reduces production costs, and ensures the product's quality and delivery time. In the composition process, first of all, determining the service composition strategy based on the decomposition of customer's manufacturing requirements should be considered, specifically, the basic strategies for selecting cloud manufacturing services according to the manufacturing tasks; secondly, service-matching degree, composition harmony degree, cloud-entropy, and other problems should be addressed. Based on the above analysis, the following cloud manufacturing service composition strategies are given: 1. All manufacturing tasks must be allocated to one or more corresponding public cloud services or private cloud services to complete. 2. Manufacturing tasks that include core product technologies and need to be kept strictly confidential must be allocated to enterprise private cloud services to be completed. 3. The comprehensive optimization objective of the cloud manufacturing service composition scheme should be as close to the ideal value as possible, which is composed of service-matching degree, composition harmony degree, cloud-entropy, and so on. Service composition according to the above strategies could provide support for the construction of a scientific and rational organization structure of the cloud manufacturing and the optimization of the use of distributed manufacturing resources. 2 MODELING CLOUD MANUFACTURING SERVICE COMPOSITION 2.1 Problem Description Suppose that a customized product manufacturing project can be broken down into l sub-tasks, namely W= {w1, ..., wl}, it is necessary to select the appropriate public cloud services in n public cloud service sets and enterprise private cloud services in m private cloud service sets from the cloud resource pool to complete. The m private cloud service sets are represented as Q1, ..., Qm. The number of private cloud services contained in each set is represented as k1, ..., km. The ith private cloud service set is represented as Q, = {Q^, ..., Qi:ki}. The n public cloud service sets are represented as Qm+i, ..., Qm+n. The number of public cloud services contained in each set is represented as km+1, ..., km+n. The ith public cloud service set is represented as Q,={Q^, ..., Qi>ki}. The l manufacturing tasks are allocated to the most m+n j appropriate public cloud services in ^ k and the most suitable enterprise private cloud services in m J / k to compete for the task collaboratively. The allocation relationship is shown in Fig. 1. 02 = { Ö2k2 } ßm = {, ßm2,", ßmk, } ßm+i = {ßm+1Д ' ßm+i,2 ' " ' ' ßm+i,km+1 } ßm+2 = {ßm+2,1, ßm+2,2,", ßm+2,k,+2 } ßm+n = {ßm +„Д, ßm+„,2,-, ßm+„,k,+, [ Fig. 1. Cloud manufacturing service composition mapping relationship 2.2 Service-Matching Degree Definition and Calculation In cloud manufacturing, whether manufacturing services are suitable for manufacturing tasks has a significant influence on the cloud manufacturing service composition. The higher the suitability of cloud manufacturing services for manufacturing tasks, the greater the ability to control the manufacturing tasks, the better the completion quality of the manufacturing tasks. Therefore, in the process of service composition, considering the suitability of cloud manufacturing services for manufacturing tasks is of great significance, which is defined as the service-matching degree. • Definition 1: Service-matching degree: the quantitative level of suitability of cloud manufacturing services (including private and public cloud services) for the manufacturing tasks. Many factors may affect the suitability of public and private cloud services for manufacturing tasks, such as the geographic location of manufacturing resources mapped by cloud manufacturing services, the frequency and experience of implementation of similar manufacturing tasks, the software and hardware strength of manufacturing resources mapped by cloud manufacturing services, the interest in undertaking the manufacturing tasks, and so on. The specific description is as follows: 1. Capacity factor. Capacity matrix C = (ci])(m+nyxi is introduced to describe the manufacturing capacity of m private cloud services and n public ones. ci]-(1 f )) f < f , (13) where Enn = Rand(n) x He + En, Ex = f , En = fnax- f )/C3, where C3 is the control coefficient; He = En / C4, where C4 is the control coefficient; In Eqs. (12) and (13), fmax, fmm and f are the maximum, minimum and average fitness values of the population respectively; f is the larger of the two crossover individual fitness values; f is the fitness value of the mutation individual; pcmax and pcmin respectively represent the maximum mutation and minimum crossover probabilities of the population; pmmax and pmmin respectively represent the maximum and minimum mutation probabilities of the population; and къ ..., k4 are constants of [0,1] with k = k2 = 0.5, and k3=k4 = 1. According to Eqs. (12) and (13), the initial values of pc and pm are larger, but gradually reduced with the continuous evolution of the population. In the same generation, different individuals have different crossover and mutation probabilities. The individuals with high fitness value should be protected, and their crossover and mutation probabilities are reduced. The individuals with low fitness value should increase their crossover and mutation probabilities. Therefore, each individual in each generation has different crossover and mutation probabilities, so as to realize the adaptive crossover and mutation. Entropy En, expected value Ex, hyper entropy He and standard deviation Enn are important control parameters of the normal cloud model. Among them, the value of standard deviation Enn is mainly affected by the two parameters of entropy En and hyper entropy He. Entropy En reflects the steepness of the cloud model; expected value Ex reflects the horizontal position of the cloud model; hyper entropy He reflects the cloud particle discrete degree [34]. Cloud particles fluctuate near the expected value, and the fluctuation degree is controlled by He. The stability of the normal cloud model is decreased when the hyper entropy He is too large, and the randomness of the normal cloud model is declined when the hyper entropy He is too small. The algorithm ensures that the best individuals around the maximum fitness value maintain the stable tendency of a normal cloud model, improving the search ability of the individual with a lower fitness value, thus, producing new individuals in larger space, improving the randomness of the algorithm effectively, and restraining premature convergence. The specific steps of the cloud-entropy genetic algorithm are discussed as follows. 3.1 Encoding Strategy The algorithm uses a binary encoding method, in which each chromosome is divided into l segments, and each gene segment represents a manufacturing task. Assuming that the number of the candidate public cloud services to complete the ith manufacturing task is ni, and that of the candidate private ones to complete the ith manufacturing task is mi, then the total number of the candidate cloud manufacturing services to complete the ith manufacturing task is mi+ni. In encoding, the genetic value of 0 indicates that the manufacturing task corresponding to the gene segment is not allocated to any public or private cloud c m services; whereas the genetic value of 1 indicates that the manufacturing task corresponding to the gene segment is allocated to a public/ private cloud service. Such encoding can also be used to decode the chromosome. An example of such encoding is given in Fig. 4, which indicates that the manufacturing task w1 is allocated to the first service of the private cloud service set Q1 with a genetic value of 1. In the same segment, the other genetic values corresponding to the other candidate cloud services are 0. Similarly, the manufacturing task wl is allocated to the third service of the public cloud service set Qm+n. 3.2 Fitness Function Cloud manufacturing service composition is a nonlinear and multi-objective optimization problem, and there is a certain link among multiple objectives. It is difficult to directly give the optimal value of such a multi-objective optimization problem, but the ideal positive point (i.e. the most desired value) and the ideal negative point (i.e. the least expected value) of the single objective can be given easily. Therefore, it is more appropriate to choose the ideal point method to construct the fitness function in the cloud-entropy genetic algorithm. The ideal point method is based on the distance between the objective function value and the ideal point to evaluate the advantages and disadvantages of the composition service scheme. The smaller the distance, the better the scheme. The ideal point is composed of the ideal value of each objective function. The ideal value can be specified by the decision maker, and also can be obtained from the single objective optimal value. According to the above analysis, the evaluation function of the cloud manufacturing service composition problem can be constructed: where ( I1,12,I3 ) is the ideal point, which is composed of the optimal values of three single objective functions. (Ij, I2, I3) is the objective function value of the cloud manufacturing service composition. I is the distance between objective function value of the cloud manufacturing service composition and the ideal point, namely deviation. The quantity level and the important degree of the objective functions I1, I2 and I3 are not the same, so they should be dealt with dimensionless correction and be distinguished with weight coefficients. Thus the evaluation function can be modified as follows: I -1* min I ' = Je1(I--£- )2 +f2(IV^)2 +£3(IV3)2> (15) I -1* I* I I where ej, s2 and s3 are the weight coefficients of the objective functions, and sj + e2+s3 = 1. Their values can be determined by the expert evaluation. According to the above analysis, the fitness function of CEGA is constructed as: f (i) = U -lSl ( ^ )2 + ^ (12ф )2 + £з ( ^ )2, (16) where f(i) is the fitness function of the ith chromosome; Ii, /2 and 13 are the objective function values corresponding to the ith chromosome; and U is a positive and sufficiently large number. 3.3 Selection Operation The roulette method is adopted to perform the selection operation. According to its fitness, each generation has the corresponding probability to be copied to the next generation. Given that the population size is popsize, and the fitness of the individual is f(i), then the selection probability to be copied to the next generation can be calculated as follows. f (i) mm in I (Ii - I*)2 + (12 - I*)2 + (I3 - I*)2 Pi =- (14) Z f (i) (17) mi+ni m+ni Private cloud service sets Public cloud service sets Private cloud service sets Public cloud service sets Qi Q2 Qm+1 Qi Q2 Qm+1 1 0 0 0 0 0 0 Manufacturing task w i Qm+n 0 0 0 0 0 0 0 0 0 1 Manufacturing task wi Fig. 4. Example of the encoding method 2 =1 3.4 Crossover Operation The crossover probability pc of the population is calculated by Eq. (12), and the expected value Ex, entropy En, and hyper entropy He of the crossover operator are calculated. Produce two parent individuals by the normal cloud model. The operation method called "double crossover points" is adopted; this is randomly selecting two gene segments as crossover points, exchanging gene encodes of the crossover points between the first parent individual and the second parent individual, the unchosen gene encodings remaining the same; thus, two offspring individuals are produced. The crossover operation is shown in Fig. 5. Gene segment I Gene segment 2 Gene segment 3 Parent I 10001 1000000001 000100000100110 Randomly select two gene crossover operation. Offspring 1 0 001010000001 10 Offspring 2100010000100001 Fig. 5. Crossover operation 3.5 Mutation Operation Calculate the mutation probability pm of the population according to Eq. (13), the expected value Ex, entropy En, and hyper entropy He of the mutation operator. Produce a new individual by the normal cloud model. When pm is bigger than the given threshold, randomly select a gene segment in the chromosome, one of encoded gene in the selected gene segment is set to 1, the others in the same gene segment are set to 0, then a new chromosome is produced. The mutation operation is shown in Fig. 6. CEGA randomly generates an initial population, repeats steps (Sections) 3.2 to 3.5 until the maximum generation equals the setting value 150 (i.e. maxgen=150) or other termination conditions of the algorithm occur, and outputs the optimal results. Gene segment 1 Gene segment 2 Gene segment 3 Parent 000100000101010 The first gene encoding in the gene segment 2 is set to 1 for mutation Offspring 000101000001010 Fig. 6. Mutation operation 4 APPLICATION EXAMPLE A company in Foshan city has built two manufacturing workshops. A private cloud service platform is built inside the enterprise for sharing the equipment, raw materials and other manufacturing resources for the internal manufacturing departments. The manufacturing tasks of the enterprise are synergistically implemented by the internal manufacturing resources (i.e. the enterprise private cloud services) and the external manufacturing resources (i.e. the public cloud services). Take the production of five sets of Automatic Guided Vehicles (AGVs) as an example to illustrate the specific application of the proposed cloud manufacturing service composition optimization algorithm. The manufacturing task of the AGV can be divided into six sub-tasks: car body task w1, driving device task w2, auxiliary equipment task w3, power supply system task w4, auxiliary control system task w5 and main control system task w6. All of the sub-tasks are allocated according to the proposed CEGA. The available public cloud service sets for the allocation of the manufacturing tasks w1, w2, w3, w4, and w5 are Qb Q2, Q3, Q4, and Q5, respectively. The numbers of public cloud services contained in the five sets are 3, 3, 4, 2 and 2, respectively. The manufacturing task w6 involves the core technology of the enterprise with a high level of confidentiality, and must be implemented by the enterprise private cloud service set Q6, which contains two private cloud services. In summary, the six manufacturing tasks wb w2, w3, w4, w5, and w6 are allocated to 14 public cloud services and two private cloud services. According to the manufacturing task number, the candidate cloud manufacturing services are arranged in a sequence. The influence factor vectors are obtained according to the definition of service-matching degree. Capacity factor vector: C = [0.6, 0.2, 0.6, 0.2, 0.4, 0.6, 0.8, 0.6, 0.4, 0.2, 0.2, 0.6, 0.8, 0.2, 0.6, 0.8] Desire factor vector: D = [0.7, 0.6, 0.6, 0.8, 0.7, 0.6, 0.5, 0.8, 0.8, 0.9, 0.7, 0.8, 0.7, 0.8, 0.8, 0.7] Equipment factor vector: E= [0.6, 0.6, 0.4, 0.4, 0.6, 0.6, 0.6, 0.4, 0.6, 0.6, 0.4, 0.6, 0.4, 0.6, 0.8, 0.6] Position factor vector: P = [1, 0.6, 1, 0.6, 1, 1, 1, 0.6, 1, 1, 1, 0.6, 0.6, 1, 1, 1]. Set a = 0.3, ß=0.2, у = 0.2, and 8 = 0.3, and then the service-matching degree vector can be obtained: V= [0.74, 0.44, 0.68, 0.48, 0.68, 0.72, 0.76, 0.60, 0.70, 0.66, 0.58, 0.64, 0.64, 0.64, 0.80, 0.80]. The execution time [hour] of each task is represented by a vector as follows: texe=[73, 70, 68, 75, 70, 60, 70, 85, 90, 88, 52, 54, 48, 60, 63, 58]. The maximum continuous working time [hour] and the maintenance time [hour] are expressed as follows: tcon=[16, 20, 16, 20, 24, 30, 20, 18, 16, 24, 20, 30, 40, 24, 22, 24]; trep= [2, 3, 2, 2, 3, 1.5, 3, 2, 1, 2, 2, 1, 3, 1.5, 2, 1]. The unit time cost [€ per hour] is expressed as follows: y = [37, 61, 47, 28, 43, 57, 47, 67, 31, 59, 45, 41, 46, 43, 49, 51]. The composition harmony degree matrix of the manufacturing task can be obtained according to the definition of the composition harmony degree. W1 W2 W3 W4 W5 W6 w1 К h12 h13 h14 h15 h16 W2 h21 h22 h23 h24 h25 h26 II К h32 h33 h34 h35 h36 W4 h41 h42 h43 h44 h45 h46 W5 h51 h52 h53 h54 h55 h56 W6 Ai h62 h63 h64 h65 h66 The matrix is a symmetric one, whose diagonal elements hn, h22, •••, and h66 are 1. h12 represents the composition harmony degree of three candidate cloud manufacturing services for manufacturing task w1 and three for manufacturing task w2. Other elements of the matrix are similar to h12. Their values can be calculated as follows: "0.769 0.667 0.589" 0.588 0.667 0.770 0.667 0.626 0.910 h12 = h13 = 1.000 0.589 0.668 0.770 0.625 0.834 0.557 0.770 0.909 0.528 0.627 0.716 h14 = 0.835 0.628 0.669 0.771 0.627 0.668 h15 = 0.594 0.770 0.671 0.834 0.670 0.770 hi6 = 0.770 0.834 0.909 0.770 0.909 0.909 h23 = 0.667 0.625 0.589 0.527 0.769 0.668 0.590 0.627 0.770 0.592 0.671 0.671 h24 = h26 = h35 = h45 = 0.772 0.591 0.669 0.669 0.589 0.770 0.834 0.910 0.910 0.770 0.910 0.910 0.593 0.770" 0.676 0.670 0.568 0.836 0.777 0.593 0.625 0.589" 0.715 0.556 h25 = h34 = h36 = h46 = 0.774 0.716 0.593 0.626 0.834 0.667 0.717 0.590" 0.673 0.672 0.776 0.633 0.535 0.774 0.910 0.770" 0.771 0.911 0.836 0.837 0.910 0.911 0.770 0.834" 0.910 0.833 h56 = 0.910 0.834 0.769 0.910 The cloud-entropy vector of the manufacturing task can be obtained according to the definition of cloud-entropy and Eq. (4): EnC=[1.784, 1.539, 1.658, 1.602, 1.341, 1.622, 1.539, 1.658, 1.946, 1.482, 1.274, 0.764, 0.561, 1.202, 1.307, 1.144]. The delivery time constraint is 480, and the cost constraint is 100000. Specifically, max(tb t2, ..., t6) <480 and С = YÜ'-f • Уг ^ Ю0000. MATLAB R2015a is used to realize the CEGA programming. Set population size Popsize = 50, maximum generation Maxgen = 150, and U = 100. The weight coefficients of the objective functions are given as s1 = 0.4, e2 = 0.3, and s3 = 0.3, respectively. The ideal point obtained by the single objective optimal function is (4.30, 12.768, 9.135). After 48 iterations of the CEGA, the optimal fitness value of the population is 99.934, the service-matching degree is 4, the composition harmony degree is 11.799, the cloud-entropy is 8.708, the delivery time is 375, the cost is 80490, the chromosome encoding of the optimal solution is 1001001000100110, and the distance between the optimal objective function value and the ideal point is 1.101. The chromosome encoding of the optimal solution is shown in Fig. 7, where manufacturing task w1 is allocated to the first cloud service of the public cloud service set task w2 to the first cloud service of the public cloud service set Q2; task w3 to the first cloud service of the public cloud service set Q3; task w4 to the first cloud service of the public cloud service set Q4; task w5 to the second cloud service of the public cloud service set Q5; and task w6 to the first cloud service of the private cloud service set Q6. The average running time of the algorithm is 11.63 s, and the evolution curves of CEGA are shown in Fig. 8. Fig. 8a) shows the optimal individual fitness value evolution curve, Fig. 8b the service-matching degree evolution curve, Fig. 8c the composition harmony degree evolution curve, Fig. 8d the cloud-entropy evolution curve, Fig. 8e the execution time evolution curve, Fig. 8f the execution cost evolution curve, Fig. 8g the deviation evolution curve and Fig. 8h the three-dimensional scatter diagram of CEGA solution process. The population tends to be stable when it evolves to the forty-eighth generation. Q\ Qi Qi Qa QS 06 1 0 0 1 0 0 1 0 0 0 1 0 0 1 1 0 W\ w2 w3 w4 w5 w6 Fig. 7. The optimal solution of the AGV manufacturing task Fig. 9 shows the evolution curves of the proposed CEGA and Standard Genetic Algorithm (SGA) under the same conditions. As can be seen from the figure, CEGA converges to the optimal solution of a) 9 8.5 Q. 0 8 "Ó =5 7.5 0 0 7 6.5 50 100 Generation d) 50 100 Generation cu 4.4 CD DJ cu -о 4.2 D) С 'si 4 га E 8 3.8 '£ cu w 3.6 150 0 b) 350 150 0 50 100 Generation e) 50 100 Generation er 12.5 CT T3 > 12 n o E 11.5 ar -С n 11 oi (Л o 10.5 CP E 10 o 0150 0 10 S9.5 "co 8 9 с о '■g 8.5 о 50 100 Generation 50 100 Generation 3.5 3 с ■В 2 5 co > CD 2 Q 15 1 g) 50 100 Generation h) matching degree Fig. 8. CEGA evolution curves; a) the optimal individual fitness value evolution curve; b) the service-matching degree evolution curve; c) the composition harmony degree evolution curve; d) the cloud-entropy evolution curve; e) the execution time evolution curve; f) the execution cost evolution curve; g) deviation evolution curve; and h) three-dimensional scatter diagram of CEGA solution process 99.95 99.9 99.85 99.8 0 c) 0 f) the problem after 48 generations, and SGA converges after 106 generations. On the same computer with Intel Core i3-3110M, 2.4 GHz CPU, and 4G memory, CEGA takes 11.63 s, and SGA 19.34 s. Thus, in the solution of cloud manufacturing service composition, CEGA converges faster than SGA does. Generation Fig. 9. Comparison of CEGA and SGA evolution curves 5 CONCLUSION In order to optimize the allocation of cloud manufacturing resources and improve the execution efficiency of cloud manufacturing service composition, research on the cloud manufacturing service composition optimization problem has been done with the following contributions: 1. Service-matching degree, composition harmony degree and cloud-entropy are defined, and the corresponding calculation methods are given, which broadens the scope of the influence factor set of cloud manufacturing service composition optimization. 2. The mathematical model of cloud manufacturing service composition optimization is built, which takes the service-matching degree, composition harmony degree and cloud-entropy as the multi-objective functions with the execution time, execution cost, and so on as the constraints. 3. On the introduction of the theory of a normal cloud model to improve the traditional genetic algorithm, a so-called CEGA is proposed to solve the mathematical model of cloud manufacturing service composition with multi-objectives and further verified by taking the task of AGV manufacturing as an example of its feasibility and effectiveness. Cloud manufacturing service composition optimization is a complex problem, which has many influence factors such as the service-matching degree, the composition harmony degree and the cloud-entropy as studied in this paper. Further research can include other influence factors and algorithms. 6 ACKNOWLEDGEMENTS The project was supported by the National Natural Science Foundation of China (51175187, and 51375168), the Science & Technology Foundation of Guangdong Province (2015A020220004, 2016B090918035 and 2016A020228005), and the Science & Technology Foundation of Zhanjiang City (2015A01001). The authors would like to thank the editors and anonymous reviewers for their constructive and helpful comments which helped to improve the presentation of the paper. 7 REFERENCES [1] Jo, M., Maksymyuk, T., Strykhalyuk, B., Cho, C.H. (2015). Device-to-device-based heterogeneous radio access network architecture for mobile cloud computing. IEEE Wireless Communications, vol. 22, no. 3, p. 50-58, D0l:10.1109/ MWC.2015.7143326. [2] Chun, S.-M., Kim, H.-S., Park, J.-T. (2015). CoAP-based mobility management for the internet of things. Sensors, vol. 15, no. 7, p. 16060-16082, D0I:10.3390/s150716060. [3] Li, Y.X., Yao, X.F., Xu, X.M., Jin, H. (2014). Formal verification of cloud manufacturing service composition and BPEL codes generation based on extended process calculus. Information Technology Journal, vol. 13, no. 11, p. 1779-1785, D0I:10.3923/itj.2014.1779.1785. [4] Li, B.-H., Zhang, L., Wang, S.-L., Tao, F., Cao, J.-W., Jiang, X.-D., Song, X., Chai, X.-D. (2010). Cloud manufacturing: a new service-oriented networked manufacturing model. Computer Integrated Manufacturing Systems, vol. 16, no. 1, p. 1-7. [5] Tao, F., LaiLi, Y.J., Xu, L.D., Zhang, L. (2013). FC-PACO-RM: A parallel method for service composition optimal-selection in cloud manufacturing system. IEEE Transactions on Industrial Informatics, vol. 9, no. 4, p. 2023-2033, D0I:10.1109/ TII.2012.2232936. [6] Omid, F.V., Mahmoud, H. (2014). A platform for optimization in distributed manufacturing enterprises based on cloud manufacturing paradigm. International Journal of Computer Integrated Manufacturing, vol. 27, no. 11, p. 1031-1054, D0I:10.1080/0951192X.2013.874582. [7] Dantas, J., Matos, R., Araujo, J., Maciel, P. (2015). Eucalyptus-based private clouds: availability modeling and comparison to the cost of a public cloud. Computing, vol. 97, no. 11, p. 11211140, DOI:10.1007/S00607-015-0447-8. [8] Hu, J.L., Den, J.B., Wu, J.B. (2013). A green private cloud architecture with global collaboration. Telecommunication Systems, vol. 52, no. 2, p. 1269-1279, D0I:10.1007/s11235-011-9639-5. [9] Lartigau, J., Xu, X.F., Nie, L.S., Zhan, D.C. (2015). Cloud manufacturing service composition based on QoS with geo-perspective transportation using an improved artificial bee colony optimization algorithm. International Journal of Production Research, vol. 53, no. 14, p. 4380-4404, D0I:10. 1080/00207543.2015.1005765. [10] Jeong, H.Y., Lee, Y.S. (2012). CS P based web service composition model with buffer at the business logic process level. Journal of Internet Technology, vol. 13, no. 3, p. 501508, D0l:10.6138/JIT.2012.501.508. [11] Castejón, C., Carbone, G., Garcfa-Prada, J.C., Ceccarelli, M. (2010). A multi-objective optimization of a robotic arm for service tasks. Strojniški vestnik - Journal of Mechanical Engineering, vol. 56, no. 5, p. 316-329. [12] Jovanovic, J.R., Milanovic, D.D., Djukic, R.D. (2014). Manufacturing cycle time analysis and scheduling to optimize its duration. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 7-8, p. 512-524, D0I:10.5545/sv-jme.2013.1523. [13] Gaaloul, W., Baina, K., Godart, C. (2008). Log-based mining techniques applied to web service composition reengineering. Service Oriented Computing and Application, vol. 2, no. 2, p. 93-110, DOI:10.1007/S11761-008-0023-6. [14] Stegaru, G., Stanescu, A.M. (2015). Towards quality driven web service composition. UPB Scientific Bullintin, Series C: Electrical Engineering, vol. 77, no. 2, p. 43-54. [15] Omid, M., Safi-Esfahani, F., Nadimi-Shahraki, M.-H. (2014). A framework for context-aware web service composition using planning techniques. Multiagent and Grid Systems, vol. 10, no. 4, p. 185-197, D0I:10.3233/MGS-140222. [16] lordache, R., Moldoveanu, F. (2015). An end to end web service composition based on QoS preferences. UPB Scientific Bullintin Series C: Electrical Engineering, vol. 77, no. 3, p. 3-16. [17] Berlec, T., Kusar, J., Zerovnik, J., Starbek, M. (2014). Optimization of a product batch quantity. Strojniški vestnik -Journal of Mechanical Engineering, vol. 60, no. 1, p. 35-42, D0I:10.5545/sv-jme.2013.1009. [18] Avitabile, P., Nonis, C., Obando, S.E. (2014). System model modes developed from expansion of uncoupled component dynamic data. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 5, p. 287-297, D0I:10.5545/sv-jme.2014.1829. [19] Florjanič, B., Govekar, E., Kuzman, K. (2013). Neural network-based model for supporting the expert driven project estimation process in mold manufacturing. Strojniški vestnik - Journal of Mechanical Engineering, vol. 59, no. 1, p. 3-13, D0I:10.5545/sv-jme.2012.747. [20] Suzić, N., Stevanov, B., Ćosić, I., Anišić, Z., Sremčev, N. (2012). Customizing products through application of group technology: a case study of furniture manufacturing. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 12, p. 724731, D0I:10.5545/sv-jme.2012.708. [21] Chifu, V.R., Salomie, I., Pop, C.B., Niculici, A.N., Suia, D.S. (2014). Exploring the selection of the optimal web service composition through ant colony optimization. Computing and Informatics, vol. 33, no. 5, p. 1047-1064. [22] Xiao, Y.Y., Li, B.H., Zhuang, C.H., Lin, T.Y. (2015). Distributed supply scheduling oriented to multi-variety customization. Computer Integrated Manufacturing Systems, vol. 21, no. 10, p. 800-812, D0I:10.13196/j.cims.2015.03.025. [23] Dong, Y.F., Guo, G. (2014). Evaluation and selection approach for cloud manufacturing service based on template and global trust degree. Computer Integrated Manufacturing Systems, vol. 20, no. 1, p. 207-214, D0I:10.13196/j.cims.2014.01. dongyuanfa.0207.8.20140126. [24] Jing, S.K., Jiang, H., Xu, W.T, Zhou, J.T. (2014). Cloud manufacturing service composition considering execution reliability. Journal of Computer-Aided Design & Computer Graphics, vol. 26, no. 3, p. 392-400. [25] Zhao, N., Niu, Z.W., Guo, W. (2012). Quantum harmony search method for design knowledge resource serialization combination in cloud manufacturing environment. Computer Integrated Manufacturing Systems, vol. 18, no. 7, p. 14351443. [26] Bao, J.P, Song, N., Tao, B. (2013). A knowledge-based service composition algorithm. Journal of Xi'an Jiaotong University, vol. 47, no. 8, p. 1-6+67. [27] Yao, X.F., Lian, Z.T., Li, Y.X., Jin, H., Xu, C., Tan, W., Zhang, J., Lin, Y. (2012). Service-oriented architecture and integrated development environment for cloud manufacturing. Computer Integrated Manufacturing Systems, vol. 18, no. 10, p. 23122321. [28] Huang, B.Q., Li, C.H., Tao, F. (2013). A chaos control optimal algorithm for QoS-based service composition selection in cloud manufacturing system. Enterprise Information Systems, vol. 8, no. 4, p. 445-463, D0I:10.1080/17517575.2013.7923 96. [29] Bao, B.F., Yang, Y., Li, L.T., Li, F., Liu, A.J., Liu, N. (2014). Multi-objective optimization for task allocation of product customization collaborative development. Computer Integrated Manufacturing Systems, vol. 20, no. 4, p. 739-746, D0I:10.13196/j.cims.2014.04.baobeifang.0739.8.2014042. [30] Chen, Y.J., Yao, X.F., Xu, D.L. (2010). Job shop scheduling with profit and entropy as performance measures. Journal of Beijing University of Technology, vol. 36, no. 10, p. 13051311. [31] Zhang, F., Cao, J.W., Hwang, K., Li, K., Khan, S.U. (2015). Adaptive workflow scheduling on cloud computing platforms with iterativeOrdinal optimization. IEEE Transactions on Cloud Computing, vol. 3, no. 2, p. 156-168, D0I:10.1109/ TCC.2014.2350490. [32] Ueda, Y., Horio, K., Kubota, R. (2014). A modified real-coded genetic algorithm considering with fitness-based variability. International Journal of Innovative Computing, Information and Control, vol. 10, no. 4, p. 156-168. [33] Lin, C.C. (2012). Hierarchical path planning for manipulators based on genetic algorithm with non-random initial population. International Journal of Innovative Computing, Information and Control, vol. 8, no. 3, p. 1773-1786. [34] Dai, C.H., Zhu, Y.F., Chen, W.R., Lin, J.H. (2007). Cloud model based genetic algorithm and its application. Acta Electronic Sinica, vol. 35, no. 7, p. 1419-1424. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)10, 591-602 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3553 Original Scientific Paper Parameters Influencing the Dynamic Behaviour of the Carrying Structure of a Type H Portal Crane Rade Vasiljević1'1' - Milomir Gašić2 - Mile Savković2 1 High Technical School of Vocational Studies, Serbia 2 University of Kragujevac, Faculty of Mechanical and Civil Engineering, Serbia This paper deals with the problems investigating of parameters that influence the dynamic behaviour of the carrying structure of type H portal cranes with high carrying capacities excited by the motion of the crane and the load swing. The introduction of the paper presents the importance of the problem of the dynamics of structures. In accordance with the given problem, dynamic spatial models of the boom and portal of the crane were formed and solved. The main part of the paper presents a modal analysis and determines dynamic responses of the space structure with the application of the finite element method and direct integration. The influence of changes in speed and acceleration on the dynamic response of the structure was analysed. The influence of the change of boom rotation angle on the dynamic response of the structure was also analysed. The investigation results provide significant bases and conclusions for designing portal cranes and continuation of research in this field. Keywords: type H portal crane, rotating boom, excitation, analysis of dynamic behaviour, finite element method, dynamic response Highlights • Definition of excitation in a portal crane. • Dynamic model of portal crane boom. • Finite element model of the carrying structure of a type H portal crane. • Dynamic response of the carrying structure to excitation. 0 INTRODUCTION Research into the dynamic behaviour of structures is present in mechanical and civil engineering. Structures of various purposes are exposed to variable influences. These influences can be common [1], which are variable only in time, and influences of moving loads, which are variable in time and space [2] and [3]. The action of these influences on structures is reflected in the occurrence of dynamic displacements, velocities, and accelerations of the structure. In mechanical engineering, this problem is present, above all, in the field of crane engineering. One of the representatives of crane machines, which are the subject of this paper, is the portal-rotating crane. These are crane machines with a boom that is connected to the portal and that has the possibility of a 360° rotation. These cranes are interesting for dynamic analysis due to their high carrying capacities and the fact that they have a small base in proportion to their height. In portal-rotating cranes, there are influences that are variable in time. Load swing contributes significantly to the dynamic loading of cranes. The importance of the problem of load swing is considered and confirmed in papers [1] and [4]. Furthermore, the importance of this problem is defined in standards [5]. The first papers considered planar dynamic models, and later papers published specifically in the field of dynamics of portal cranes mostly treated planar dynamic models [6] and [7]. The authors of this paper found a small number of papers with dynamic spatial models which dealt with the problems of the dynamic behaviour of portal-rotating cranes [8]. Therefore, the search was directed, in a wider sense, to the group of rotating cranes. First, a group of dynamic models that deal with the problems of the dynamic behaviour of rotating cranes was considered [9] to [11]. Next, a group of dynamic models that deal with the problems of the control of rotating cranes was reviewed [12] and [13]. In almost all the above-mentioned papers [6] to [13], the authors modelled portal-rotating cranes as direct or discrete-continuous models. Vasiljević et al. show [14] that more accurate results are obtained by the application of consistent masses. The papers mostly provide an investigation of the dynamic models of only booms of portal cranes, but few of them investigate dynamic models which take into account the carrying structure. In accordance with this fact, this paper should contribute to the knowledge of the dynamic behaviour of the carrying structure of the portal crane excited by the motion of the crane and the load swing. The solution to this problem requires appropriate modelling of the boom and the carrying structure. The approach to modelling of the dynamic model of the portal crane is such that the whole portal crane is divided into two subsystems, i.e. the carrying structure and the boom. The relationship between the carrying structure and the boom is simplified by reducing the influence of load and the dead weight of the boom to the corresponding points of the portal. The carrying structure of the portal crane is modelled as a spatial linear beam. In reality, spatial linear beams are systems with continuously distributed masses, whereas the spatial model of the carrying structure is more complex. For these reasons, the continuously distributed mass is replaced with the consistent mass. Hence, the consistent mass matrix of the system is formulated on the basis of the same interpolation functions that are used for the derivation of stiffness matrices. The spatial linear beam is modelled by means of a space girder element, with two node points and six unknowns in each node. For the adopted dynamic model of the portal-rotating crane, modal analysis of oscillation was conducted in the first step, and it forms the basis for a good understanding of the dynamic characteristics of the carrying structure. In the second step, the dynamic response of the portal-rotating crane to excitation was determined. 1 PROBLEM SETTING Oscillation results from the elasticity of mechanisms and the carrying structure due to the action of inertial forces. The shorter the time of action of those forces, the more unfavourable the oscillations. The first main problem in the investigation of the dynamic behaviour of portal cranes is the definition of excitation. The structure, i.e. the base of the crane considered, can perform only rectilinear motion from point to point. The kinematic values of the motion of the portal crane are the speed of motion, acceleration and the distance passed. The trapeze velocity profile is most frequently used for the calculation of mechanisms for the motion of cranes. Accordingly, in this paper, it is assumed that the crane moves along its track with the trapeze velocity profile over the course of one cycle (Fig. 1). From Fig. 1, it is seen that the motion of the structure consists of three phases. In phase 0 to 1 (time ta), the structure accelerates with the constant acceleration aa so that the velocity linearly increases, and at the end of the first phase, it reaches the velocity vu. The second phase, 1 to 2, is the phase of uniform motion of the structure during the time tu, which is characterized by the constant velocity vu. In phases 2 to 3 (time td), the structure constantly decelerates with the deceleration ad until it stops. The distance passed by the structure is determined based on the profile of velocities (Fig. 1). The presented profile of velocities will be used for obtaining the dynamic response of the carrying structure excited by the motion of the crane and the load swing. v, a, S ad Fig. 1. Kinematic diagrams of motion of the crane The second main problem in the investigation of the dynamic behaviour of portal-rotating cranes is the adoption of the appropriate model. Formation of the portal crane model is approached in such a way that the influence of loads and the boom is reduced to the corresponding points of the carrying structure of the crane. An analogous approach was applied in the investigation of the dynamic behaviour of the reloading bridge. It is shown in [15]. 2 DYNAMIC MODEL OF THE BOOM 2.1 Model Description The boom is a constituent part of the portal-rotating crane. The subsystem of the boom consists of the carrying structure and the rope system. The carrying structure of the boom is connected to the platform by joints. In the system of the portal crane, the boom is observed, by idealization, as an independent subsystem. Accordingly, Fig. 2 presents a dynamic equivalent model of the portal crane boom. The model was formed for the needs of determination of dynamic loads. The boom is observed as an independent subsystem which oscillates. A dynamic equivalent model is formed in such a way that it both keeps the main dynamic characteristics of the boom and that the defined problem could be mathematically solved. From Fig. 2 it is seen that the subsystem of the boom with a load is represented with two lumped masses, two lightweight bars, and a circular disc. Discretization of the carrying structure of the boom was performed on a lightweight bar with a reduced mass on the tip. In other words, the carrying structure of the boom is represented by a lightweight bar with the length Lb and the mass mj, which is reduced to the tip of the boom. The rope system of the boom is represented by a non-elastic lightweight bar with the length Lr and the mass m2, which allows load swing. It means that the tip of the boom, i.e. the mass mj and the mass m2, are joined by the non-elastic lightweight bar with the length Lr. The rotating column is represented by a circular disc with the axial moment of inertia J and the moment of rotation T. The distance between the boom joint and the column axis is r. Fig. 2. Dynamic model of the boom Based on the recommendations in [15], the reduced mass m1 can be determined according to the following relation: m = |4... m (1) where mb is the boom mass. 2.2 Mathematical Formulation The Lagrange equations of the second kind will be used for the setting of the mathematical formulation of the formed dynamic model of the boom. The equations of motion of the boom elements are set based on the dynamic equivalent model presented in Fig. 2. The dynamic equations of motion of the system read: ((t ) + ®;2sm0(t ) = —- x (t )cosp, (2a) y(t ) + ®2siny(t ) = --1 x (t )sin^, (2b) where в is the angle of oscillation of load in the longitudinal direction, щ is the angle of oscillation of load in the side (lateral) direction, ф is the angle of rotation of the column, i.e. the boom, x is the rectilinear motion of the boom, and m is the circular frequency of load oscillation. 2.2.1 Load Oscillation The laws of non-attenuated oscillations of the load as a function of time along the generalized coordinates в and щ due to the crane acceleration given by the diagram according to Fig. 1 will be determined from Eqs. (2a) and (2b). In the first step, the motion of load, i.e. mass m2 in the longitudinal direction along the generalized coordinate в is observed in such a way that it is taken out of the system, Fig. 3. Fig. 3. Load oscillation in the longitudinal direction For determination of the law of load swing, the differential equation of relative motion of the particle, i.e. mass m2, is defined: mA^F™ + F™ + G2 + 5. (3) Jy(t )=T, (2c) Eq. (3) is projected onto the axes t and n of the natural trihedron, which is connected to the mass m2. The details of obtaining the law of load oscillation in the longitudinal direction are given. As there is an analogy between load oscillation in these two directions, only the final form will be given for the side direction. The law of load oscillation in the longitudinal direction obtains the form: e(t ) + — sine(t ) = --1 x (t )cose(t )cos9. (4) In Eq. (4), the relation g/Lr represents the quadrant circular frequency of the load m2. For the case of small oscillations, the approximation that the angle в is small can be introduced, so that sine - в and cose~1. Furthermore, the replacement a = x(t ) is introduced. Eq. (4) is linearized so that it obtains the following form: в (t) + w2e(t)=--1 a(t)cos^. (5) Eq. (5) is most suitably solved by the method of Laplace transformation. In the first step, it is determined that: L{{s )) = -l-L ) L s2 +®2 \cosq. (6) From the kinematic diagrams from Fig. 1, the acceleration diagram is particularly suitable for solving Eq. (6). This diagram is presented in Fig. 4. It is adopted that aa = ad = h and ta = td = т. For previous calculations, it can be adopted that 0.5(ta + td)=3 to 5 s. Fig. 4. Acceleration diagram Now, according to the diagram of change of acceleration in time, Fig. 4, the expression for acceleration in the Laplace domain can be represented in the following form: (s ) = — (l - ^ - ^ + (7) where h is the amplitude of input acceleration and t, is the corresponding time step in the acceleration diagram. Eq. (6) obtains the form: L {(s )} = - h-L I i Us t-,s , tAs I 1 - e2 - e3 + e [ s (s2 +ffl2 ) L fcos^. (8) Finally, the law of load oscillation in the longitudinal direction is obtained by transforming Eq. (8) for 0(s) in the time domain by using the inverse Laplace transform. 0(t) = - ——-cos^{(l - cos mt)- [l - cos m(t -12 )] HeavisideTheta(t -12 ) -[l - cos m(t -13 )] HeavisideTheta(t -13 ) + [l - cos m(t -14 ) HeavisideTheta(t -14 )}, (9) ¥(t ) = - L„ a -sin^ {(1 - cos at )- [l - cos a(t -12 ) HeavisideTheta(t -12 ) -[l - cos aa(t -13 ) HeavisideTheta (t -13 ) + [l - cosa (t -14 ) HeavisideTheta(t -14 )}. (10) The obtained laws of load oscillation (Eqs. (9) and (10)) represent Heaviside (step) functions with a linear polynomial. The load swing contributes the most to the instability and overturning of the crane when the angle of load oscillation reaches its highest point. This case occurs when the impulses of time calculation are mutually in phase. In that case, load oscillation caused by every step of acceleration is constructively added and produces the highest amplitude of load oscillation. Eqs. (9) and (10) show that the maximum angle of oscillation appears when the following conditions are fulfilled: • each of the cosine members within Eqs. (9) and (10) is in phase, • multiplying the cosine members within Eqs. (9) and (10) by the corresponding step of the function HeavisideTheta results in obtaining the value equal to 1 (it is necessary to have a sufficient time of operation, t>t4, in order to achieve acceleration according to the diagram). The cosine members within Eqs. (9) and (10) are in phase when four steps of the acceleration diagram are performed by their constructive adding in perfect time. In the worst case, the angle of oscillation increases four times in relation to the oscillation caused by one step of input acceleration. The maximum angle of load oscillation in the longitudinal and side directions is equal to: \0 \ = \ \ = 4—. max r max (11) To produce maximum load swing, it is necessary to have acceleration steps completely in phase, i.e. it is necessary to fulfil two conditions. The first condition requires that the time interval between the first and the second steps, as well as between the third and the fourth ones should be equal to one half of the oscillation period т: t2 - tx = t2 =(0.5 + n), and t4 -13 =(0.5 + n) т, (12) where n = 1,2,3, ... The time delay is necessary because the steps in sets have different designations. Therefore, the time delay is equivalent to the excitation of the phase delay n according to the opposite signs of function, which excludes the possibility of the change of phases (and brings two functions in phase). The second condition requires that the time interval between the second and the third impulses must contain the time delay of the multiple period т: t3 -t2 = m, (13) where n = 1, 2, 3, ... In other words, this holds because the second and the third steps have the same sign so that the change of phase by 2n requires that the second and the third steps should be in phase. 2.2.2 Dynamic Loads of the Portal In accordance with the adopted generalized coordinates of oscillation of the dynamic model of the boom, the dynamic moment of bending occurs in two directions: • in the longitudinal direction, and • in the side (lateral) direction. The dynamic moment of bending in the longitudinal direction reads, Fig. 2: Mdyn,l =mig (Lb cosa + Г ) + m2 (g + Lß2 cosß) ( sin0 + Lb cosa + r) +m2Lß2 sm0(Lb sina -Lr cosß + m2Lry/2cosy(Lr cosa + r). (14) The dynamic moment of bending in the side direction reads, Fig. 2: Mdyn,s =m2 ( + LA> 2 C°w)Lr SW + m2 Lry/ 2sin^(L4 sina- Lr cos^). (15) 3 FINITE ELEMENT MODEL OF THE STRUCTURE 3.1 Model Description The whole portal-rotating crane is divided into two subsystems: the moving structure, and the boom. The relationship between the structure and the boom is simplified in such a way that the influence of load and the dead weight of the boom is reduced to the points of the upper and lower supports of the boom. The type H carrying structure of the considered portal crane is shown in Fig. 5. Fig. 5. Carrying structure (portal) of type "H" m Fig. 6. Finite element model of the type H structure The carrying structure is a rigid spatial frame. Its base has dimensions L*B. The main structural parts of the carrying structure are the legs, the slanted columns, and the lower and the upper beams. The legs are identical (height H) and stand at the same level. The slanted columns are identical (length C) and they are connected to the legs as well as to the upper and lower beams. The distance between the upper and the lower beams is equal to H0. The upper beam has the function of the upper support of the boom, and it is made of a circular ring with the diameter D. The lower beam has the function of the lower support of the boom, and it is made of a type H frame. The finite element model of the structure of the considered portal crane is shown in Fig. 6. Taking into account that the boom rotates around its axis, it can be concluded that the assumed planar model is not sufficient for describing the dynamic behaviour of the considered portal-rotating crane. The spatial model of the portal crane, to which dynamic load is reduced, is formed in accordance with the spatial dynamic model of the boom. Discretization of the carrying structure of the crane was performed on 18 finite beam elements connected by nodes. The whole system has 82 degrees of freedom of motion. The legs were modelled as one finite element with the characteristics A1, Ix1, Iy\, Iz1. The slanted columns were modelled as one finite element with the characteristics A5, Ix5, Iy5, Iz5. The upper beam (the circular ring with the diameter D) was divided into four finite elements with the characteristics A 9, Ix9, Iy9, Iz9. The lower beam (the type H frame) was divided into six finite elements with the characteristics A13, Ix13, Iy13, Iz13. The supporting structure rests on four elastic supports with the stiffness k. The main mechanical characteristics of the carrying structure, i.e. of all finite elements of the formed spatial model are the module of elasticity E, the slipping module G, and the density p. The dynamic moments, Eqs. (14) and (15), are reduced to the combination of horizontal dynamic forces P(t) and Ps(t) in the nodes of finite elements of the upper beam and the central node of the lower beam. The position of these forces relative to the longitudinal direction (the axis X) is defined by the angle ф, which represents the angle of boom rotation. The dynamic forces P(t) and 4*Pl(t)/4 make a combination of forces that opposes the dynamic moment of overturning in the longitudinal direction Mdyn,l of the platform, i.e. the rotating part of the crane: Mdyni = P (t ) • H0. (16) The dynamic forces Ps(t) and 4*Ps(t)/4 make a combination of forces which opposes the dynamic moment of overturning in the side direction Mdyn,s of the platform: Mdyns = P (t) • H0. (17) The vertical static force Pm1 of the load mass mx is reduced to the central node of the lower beam (the type H frame). The vertical static force Pm2 of the reduced boom mass m2 is also reduced to the central node of the lower beam. The total vertical static force Pm that acts in the central node of the lower support of the boom is equal to the sum of the force of the load mass and the force of the reduced boom mass: P = P. (18) The characteristic nodes of the structure model for obtaining dynamic response of the structure to excitation are the nodes at the point of elastic supports (points of wheel-rail contact). 3.2 Mathematical Formulation In accordance with the formed mathematical model of the portal crane, the differential equation of dynamic equilibrium, i.e. forced oscillations of the carrying structure reads: [M ]{U } + [K ]{U} = {P(t)}, (19) where [M] is the mass matrix of the system, [K] is the stiffness matrix of the system, {Ü} the vector of acceleration of generalized coordinates of the carrying structure, {U} the vector of displacement of generalized coordinates of the structure, and {P(t)} the vector of external load of the structure nodes. Eq. (19) is used for studying the dynamic response of the structure due to the motion of the crane. It can be decomposed to the equilibrium conditions of active and inertial forces in the direction of unknown and known displacements so that it gains the form: Muu Muk " /U ul Kuu Kuk _Mku Mkk _ lü kj _Kku Kkk _ Pu (t )| Pk (t)J .(20) As the displacements and accelerations of the supports in horizontal directions are equal to zero (Uk=0 and Ük=0), Eq. (20) obtains the form: [Мцц ]{U }+[Кцц ]{иц } = { (t) . (21) Natural frequencies of the structure are obtained by solving the algebraic equation: det(Kuu -w2stMuu ) = 0. (22) u In the direction of unknown displacements, the mass submatrix of the whole system is equal to the matrix of the carrying structure of the portal crane: Muu= [Mu (23) In the direction of unknown displacements, the stiffness submatrix of the whole system includes the matrix of the structure and the stiffness matrix of elastic supports: K [K- k82+[k ] (24) In the direction of unknown displacements, the subvector of external dynamic forces of the whole system is equal to the subvector of external dynamic forces of the carrying structure of the portal crane: p (t)=(t)}. (25) The subvector of external dynamic forces Pu(t) can be expressed in the following form: {0 0 0 0}, i = 1...4 {0 0 0 0 0 0}, i = 5...8 f Px(t) Py (0 Pu(t)T = 0 0 0 0},i = 9...12, (26) where: (27) 4 4 -Px (t ) -Py (t ) Pm 0 0 0}, i = 13 {0 0 0 0 0 0}, i = 14,15 Px (t) = p (t)cos^ + Ps (t)sin^, and Py (t) = -p (t)sinp + p (t)cosp. 4 NUMERICAL RESULTS AND DISCUSSION For determination of numerical values of parameters influencing the dynamic behaviour of the carrying structure of the type H portal crane and the formation of appropriate conclusions, the initial data were defined based on the solutions derived and manufacturers' catalogues: vc = 0.6 m/s; mt = 10000 kg; ms = 9200 kg; r = 1 m; ij = 30 m; H = 6 m; H0 = 3.675 m; An = 31500X10-6 m2, ln = 6 m, Iyn = 3.117x10-3 m4, Izn = 1.391X10-3 m4, Ixn = 2.776X10-3 m4 (n = 1 to 4); An = 31500X10-6 m2, ln = 4.243 m, Iyn = 3.117x10-3 m4, Izn = 1.391x10-3 m4, Ixn = 2.776x10-3 m4 (n = 5 to 8); An = 26500x10-6 m2, ln = 3 m, Iyn = 1.303x10-3 m4, Izn = 1.066x10-3 m4, Ixn = 1.630x10-3 m4 (n = 9 to 12); An = 25200x10-6 m2, ln = 3 m, Iyn = 1.115x10 3 m4, Izn = 1.119x10-3 m4, Ixn = 1.611x10-3 m4 (n = 13 to 18); E = 2.1x10ii Nfe2; G = 0.8x10" Nm2; p к = 107 N/m; g = 9.81 m/s2. 4.1 Excitation 7850 kg/m3; The motion of the crane along the tracks is described by the diagram of velocity vc(t), Fig. 7. -- vc (0 CO / / / / / / \ \ \ 4 \ \ \ f у N \ 4 / / / . / / // л' / / у / / \ \ 4 N 4 \ V \ \ 4 \\ Time [sj Fig. 7. Time diagram of the crane velocity vc(t) The change of angle of load oscillation in the longitudinal direction в according to Eq. (9) is presented in Fig. 8. The maximum angle of load oscillation is 0.0815 rad. Й 0.05 0.00 -0.05 0 2 4 6 8 10 12 Time [sj Fig. 8. Angle of load oscillation in the longitudinal direction - VC,(0 -- vc2(/) / / \ / л / » 1 ! 1 If 1 1/ У \ / 1 / nI 1/ V '/ \ \ yi \ t \ / 4 / W V \ A \\jl \ 1 V \ 1 \J \J 6 Time [sj Fig. 9. Angle of load oscillation in the longitudinal direction for different speeds of motion of the crane Fig. 9, based on Eq. (9), presents the angle of load oscillation in the longitudinal direction for two different speeds of the crane. The speed vc1 is the nominal speed of the crane, whereas the speed vc2 represents the extreme performance of the portal crane. Fig. 9 shows that the angle of load oscillation increases with the increase in the speed of the motion of the portal. The change of the angle of load oscillation in the lateral direction у is completely the same. The change of the dynamic bending moment in the longitudinal and side directions (8, у) according to Eqs. (14) and (15) is presented in Figs. 10 and 11. Z 3.49X106 -о £ 3.45Ч06 I-i- ti 2 4 6 8 10 12 Time [sj Fig. 10. Dynamic bending moment in the longitudinal direction 0 2 4 6 8 10 12 Time [sj Fig. 11. Dynamic bending moment in the side direction Detailed analysis of excitation in the portal crane is shown in [16]. 4.2 Response of the Carrying Structure to Excitation 4.2.1 Modal Analysis Table 1 presents the values of the first two frequencies of oscillation of the carrying structure of the portal-rotating crane, Eq. (22). The values of natural frequencies represent the first and most important element of estimation of dynamic stability of cranes in the first phase of the design of new solutions [14]. The crane has good dynamic behaviour if its first frequency of oscillation is high. According to this aspect, the carrying structure of the considered portal crane has favourable dynamic stability (/1 = 2.65 Hz). Table 1. Frequencies Mode no. Period Tst [s] Circ. freq. rnst [rad/s] Frequency/st [Hz] 1 0.378 16.63 2.65 2 0.344 18.78 2.99 The verification of the mathematical model of the carrying structure (Fig. 6) was done by creating an FE model in the programme package SAP2000® [17]. Using the modal analysis in SAP2000, the assumed first 12 frequencies of the carrying structure were determined. The first two frequencies are / = (2.66; 3.01) Hz. By comparing the values of frequencies obtained through the mathematical model with the finite element approach and in a purely numerical way (FEM software SAP2000), coinciding well between the results of the first two frequencies and the relative deviations Л = (0.38; 0.66) % is observed. Other frequencies also coincide well. Thus, for example, for the following four frequencies the relative deviation is up to 2.5 %, which is very good for the spatial model. 4.2.2 Dynamic Displacements The evaluation of the quality of new solutions of cranes is given based on the maximum values of dynamic displacements. Eq. (21) for investigating the dynamic behaviour of the portal-rotating crane was solved by means of the direct, step-by-step integration method. In the Mathematica® program, the original module is written based on the Newmark integration method [18]. The time interval of integration was chosen to be Л?1 = 0.01. Dynamic displacements for all degrees of freedom of the model were obtained, but only characteristic displacements are presented here. The dynamic response of the portal to excitation is, before all, contained in the dynamic displacement of elastic supports of the model for different positions of the boom relative to the direction of motion of the crane. Figs. 12 through 16 present the comparative change of dynamic displacements of nodes 1 to 4 for the boom positions relative to the rectilinear motion of the portal ф = 0°, 45°, 90°, 135°, and 180°. Fig. 12. Dynamic displacements of nodes 1 to 4; ф = 0 Fig. 16. Dynamic displacements of nodes 1 to 4; ф = 180° -0.10 11111 Uz i uz2 uz3 Uz4 A max ------ 4.-V/V, л ..... /чЛ Л,/ ч / \ \ _ ' \А ч v **' .4* / /' \ "ЧХ '"■N /*" \ v' 5.0 5.2 5.8 5.4 5.6 Time [s] Fig. 13. Dynamic displacements of nodes 1 to 4; ф = 45° 0.15 -0.10 5.0 U zi, Uz- UZ2, U 24 •л -ч \ -У \ / 4 'Ч \ / \ / / 1 - ' / V ч/ Ч/ Time [s] Fig. 14. Dynamic displacements of nodes 1 to 4; ф = 90° -0.10 5.0 ,1111 Uz, Uz2 Uz3 uz4 /\ \ / \ Л f '*-' Ч.' ^/Л ; ч v /а >-: ' \ X / / ч— \ / \ / Time [sj Fig. 15. Dynamic displacements of nodes 1 to 4; ф = 135° As the carrying structure is symmetric on one side, and as the dynamic moment is much higher in the longitudinal direction than in the lateral one, on the other side, dynamic displacements of nodes 1 and 2, i.e. 3 and 4 for the angle of the boom 0° and 180° almost coincide. For the angle of the boom 0°, displacements of nodes 1 and 2 are positive, while for the angle of the boom 180° displacements of nodes 3 and 4 are positive. Analogously, displacements of nodes 1 and 3, i.e. 2 and 4 for the angle of the boom 90° almost coincide. Displacements of nodes 1 and 3 are positive. For the angle of the boom 45°, displacement of node 1 stands out because it is positive. Analogously, for the angle 135° the positive displacement of node 3 stands out. Based on the analysis of displacements of elastic supports, Figs. 12 through 16, it was shown that the elastic support - node 1 has the largest positive displacement for the value of the angle of boom rotation of ф = 45°. The maximum displacement of node 1 is 78.93 mm. In this case, lifting of that leg of the portal crane most often occurs in practice. Therefore, the carrying structure of the considered portal crane is most sensitive to the boom position of 45°. Accordingly, this position is most relevant for describing the dynamic state of the carrying structure. In accordance with the parameters included in the mathematical model of the considered carrying structure of the portal crane, the expected value of maximum dynamic displacement was obtained. 5 CONCLUSIONS Analysis of the dynamic behaviour of portal cranes of large carrying capacities is necessary because it provides quality bases for their optimal design. The research conducted allows thorough consideration of the parameters influencing the dynamic behaviour of portal cranes. A combined approach, which combines equations of analytical mechanics and the finite element method, was developed. In accordance with that, dynamic spatial models of the boom and the structure of the portal crane were formed and solved. It was established that the influence of load swing is critical in portal cranes. It was shown that speed and acceleration/ deceleration of the crane influence the dynamic response of the structure. For the nominal value of the crane speed of 0.6 m/s, the obtained maximum angle of oscillation of load is 0.0815 rad/s. The stability of the portal crane excited by the motion of the crane and the load swing was investigated. The dynamic response of the structure to excitation is presented through dynamic displacements of elastic supports. The boom position at which maximum displacements of elastic supports occur was established. The critical elastic support, i.e. the leg that may be lifted, was defined for this position of the boom. The maximum dynamic displacement of the critical elastic support of 78.93 mm was obtained. This paper leaves some space for the continuation of the investigation of parameters influencing the dynamic behaviour of the carrying structure in portal-rotating cranes and the optimization of such types of structures. Further work on the problems considered should lead to the improvement of the created models, which will include a larger number of parameters. Furthermore, this paper presents a quality basis for research into the dynamic behaviour of other types of portal rotating cranes (e.g. type X). 6 NOMENCLATURE [K] structural stiffness matrix [N/m] {P(t)} external force vector [N] {U} displacement vector [m] {Ü} acceleration vector [m/s2] [Muu] mass submatrix of the structure in the direction of unknown displacements [kg] [Kuu] stiffness submatrix of the structure in the direction of unknown displacements [kg] [k] stiffness matrix of elastic supports [N/m] к stiffness of elastic supports [N/m] {Pu(t)} external subvector of forces in the direction of unknown displacements [N] {Uu} displacement vector in the direction of unknown displacements [m] {Üu} acceleration vector in the direction of unknown displacements [m/s2] L length of the portal base [m] B width of the portal base [m] H height of the leg [m] H0 distance between the upper and the lower supports of the boom [m] ln length of the finite element [m] An cross-sectional area[m2] Ixn moment of inertia of cross section the x axis [m4] Iyn moment of inertia of cross-section the y axis [m4] Izn moment of inertia of cross-section the z axis [m4] p mass density of material [kg/m3] E Young's modulus [N/m2] G slipping module [N/m2] wst circular frequency of structure [rad/s] fst frequency of structure [Hz] At time interval of integration [s] p total number of time steps [-] vc speeds of the moving crane [m/s] x(t ) acceleration of the moving crane [m/s2] t time of acceleration (deceleration) of the crane [s] Ш] reduced mass of the boom [kg] m2 load mass [kg] mb boom mass [kg] g gravitational constant [m/s2] Lb boom length [m] Lr rope length [m] r distance between the boom joint and the column axis [m] в angle of load oscillation in the longitudinal direction [rad] щ angle of load oscillation in the lateral direction [rad] Ф angle of column (boom) rotation [°] x rectilinear motion of the crane [m] ю circular frequency of the load [rad/s] T moment of rotation [Nm] [M] structural mass matrix [kg] 7 REFERENCES [1] Jerman, B., Hribar, A. (2013). Dynamics of the mathematical pendulum suspended from a moving mass. Tehnički vjesnik -Technical Gazette, vol. 20, no. 1, p. 59-64. [2] Fryba, L. (1999). Vibration of Solids and Structures under Moving Loads. 3'd ed., Thomas Telford, London, D0l:10.1680/ vosasuml.35393. [3] Gašić, V., Zrnić, N., Rakin, M. (2012). Consideration of a Moving Mass Effect on Dynamic Behaviour of a Jib Crane Structure. Tehnički vjesnik - Technical Gazette, vol. 19, no. 1, p. 115-121. [4] Raubar, E., Vrančić, D. (2012). Anti-Sway System for Ship-to-Shore Cranes. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 5, p. 338-344, D0l:10.5545/sv-jme.2010.127. [5] F.E.M 1.001. (1998). Rules for the Design of Hoisting Appliances. European Materials Handling Federation, Brussels. [6] Osinski, M., Wojciech, S. (1998). Application of Nonlinear Optimisation Methods to effect Input Shaping of the Hoist Drive of an Off-Shore Crane. Nonlinear Dynamics, vol. 17, no. 4, p. 369-388, D0l:10.1023/A:1008333417693. [7] Jovanović, M., Radoičić, G., Milić, P. (2009). Selection of Finite Elements Considering Loadcases and Geometry. Proceedings of the XIXth International Conference MHCL, p. 61-66. [8] Sun, Y., Li, D. (2012). Dynamic analysis and design method study on the combined-boom system of portal crane. Applied Mechanics and Materials, vol. 152-154, p. 1645-1649, D0I:10.4028/www.scientific.net/AMM.152-154.1645. [9] Marinović, I., Sprečić, D., Jerman, B. (2012). A slewing crane payload dynamics. Tehnički vjesnik - Technical Gazette, vol. 19, no. 4, p. 907-916. [10] Maczynski, A., Szczotka, M. (2002). Comparison of Models for Dynamic Analysis of a Mobile Telescopic Crane. Journal of Theoretical and Applied Mechanics, vol. 40, no. 4, p. 10511074. [11] Fujioka, D., Rauch, A., Singhose, W. (2009). Selection of finite elements considering loadcases and geometry double-pendulum payloads. American Control Conference, p. 31363014. [12] Gustafsson, T. (1996). On the design and implementation of a rotary crane controller. European Journal of Control, vol. 2, no. 3, p. 166-175, DOI:10.1016/S0947-3580(96)70042-0. [13] Huimin, O., Naoki, U., Shigenori, S. (2010). Anti-Sway control of rotary crane only by horizontal boom motion. IEEE International Conference on Control Applications, p. 591-595. [14] Vasiljević, R.; Savković, M.; Bulatović, R. (2013). The Approaches to the Mathematical-mechanical Modeling Supporting Construction. IMK-14 - Research & Development in Heavy Machinery, vol. 19, no. 1, p. EN29-38. [15] Gašić, V. (2004). Dynamic behaviour identification of bridge type stacker - reclaimer with bucket chain booms in power plants. M.Sc. thesis, Faculty of Mechanical Engineering, [16] Vasiljević, R., Gašić, M. (2015). The Dynamic Model of the Boom Portal Cranes. IMK-14 - Research & Development in Heavy Machinery, vol. 21, no. 4, p. EN125-130. [17] Computers and Structures Inc. (2009). SAP2000® Analysis reference manual, Berkeley. [18] Bathe, K.J. (1996). Finite Element Procedures. Prentice-Hall, New Jersey. 8 APPENDIX 8.1 Equations (2a) and (2b) 8.1.1 Equations (2a) The differential equation of relative motion of the mass m2 in the longitudinal direction reads: - - - - (A1) m2ar = F'; + f;+G2 + S. Intensity of inerti al forces: Fn=m2Z ((). (A2) (A3) Eq. (A1) is projected onto the axes t and n of the natural trihedron: m2Lß (t) = -G2 sin0 (t) - m2x(t)cos^cos0(t) + +m2z (t ) sin ß(t ), (A4) m2Lß (t) = -G2 cos 0(t)- m2x(t)cos^ sin 0(t) + +m2z (t)cos0(t). (A5) For ž = 0 Eqs. (A4) and (A5) obtain the form: m2Lß (t) = -G2 sin e(t)- m2x (t)cos^cos0 (t), (A6) m2Lß (t) = -G2cosö(t)- m2x (t)cos^sin0 (t). (A7) From Eq. (A6) is determined by the law of oscillation of cargo in the longitudinal direction: 0(t ) + g sin 0 (t ) = --1 i (t )cospcos0(t ). (A8) The approximation sine ~ в and cose ~ 1 is introduced: e(t) + ge(t) = -1x(t)cosp, (A9) where g/Lr represents the quadrant circular frequency of the load ю2. The replacement a = X(t ) is introduced: e(t) + a2e(t) = -—a (t)cosq>. (A10) Lr Laplace transformation: i{0(.s )}(2 +a2 ) = -—L {a (s )}cosp, (A11) L {(s )} = -LL a2 (S ^ f cosp. (A12) ( +ю2 )) The expression for acceleration in the Laplace domain can be represented in the following form: a (s) = — (l - e'25 - ^ + e'*s). (A13) Eq. (A12) obtains the form: ^ЖЬ -H1-^^ U* (A14) Lr I s (s +a> ) I Finally, the law of load oscillation in the longitudinal direction is obtained by transforming the Eq. (A14) in the time domain by using the inverse Laplace transform (Mathematica® program used): e(t ) = cos^{(l - cos at) - -|l - cosa(t -12) )HeavisideTheta(t -12 ) --|l - cosa(t -t3 ) JHeavisideTheta(t -t3 ) + + |l - cosa(t -14 ) )HeavisideTheta(t -14)}. (A15) 8.1.2 Equations (2b) The differential equation of relative motion of the mass m2 in the side (lateral) direction reads: m2ar=F;ny + Fpn + g2 + S. Intensity of inertial forces: F™ = m2 x (t )sinp, F; = m2z (t ). (Al 6) (Al 7) (A 18) Further, the Eq. (A16) solved by analogue Eq. (A1). The law of load oscillation in the side direction reads: h у/(t) = -——2sin^{(l - cosct)- - [l - cos co(t -12 ) HeavisideTheta (t -12 ) - - [l - cos co(t -13 ) HeavisideTheta (t -13 ) + + [l - cos C (t -14 )] HeavisideTheta (t -14 )}. (A 19) 8.2 Equations (24) Eq. (24) reads: Kuu =[KUU ] 82x82+[k ], (A20) In order to execute Eq. (24), the stiffness matrix of elastic supports [k]16x16 is, by adding zero rows and zero columns, extended to the matrix with dimensions 82x82, except for 16 degrees of freedom of four nodes: ' K1,1 + k1,1 K 1,16 + k1,16 K 1,17 - K 1,82 K16,1 + k16,1 " K 16,16 + k16,16 K16,17 - ' ' K 16,82 K 17,1 ' ' K 17,16 K 17,17 - ' ' K17,82 K 82,1 K 82,16 K82,17 - ' ' K82,82 _ (A21) 8.3 FE model in SAP2000® software The FE model of the structure of the portal crane is shown in Fig. A1. Model is created in FE package SAP2000®. The first 2 mode shapes of the structure are shown in Figs. Al and A2. Fig. A1. FE model of the structure of type "H" Fig. A2. 1st mode, f1=2.66 Hz Fig. A3. 2nd mode, f2= 3.01 Hz к = Strojniški vestnik - Journal of Mechanical Engineering 62(2016)10, 603-613 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.3170 Original Scientific Paper A Numerical Analysis of Fluid Flow and Heat Transfer Characteristics of ZnO-Ethylene Glycol Nanofluid in Rectangular Microchannels Cuneyt Uysal* - Kamil Arslan - Huseyin Kurt Karabuk University, Mechanical Engineering Department, Turkey The fluid flow and heat transfer characteristics of ZnO-Ethylene glycol (EG) nanofluid through rectangular microchannels having different aspect ratios (a is 1 to 2) are numerically investigated. The different nanoparticle volume percentages (q> is 1 % to 4 %) of ZnO-EG nanofluid are used. The flow is considered under single-phase, three-dimensional, steady-state, incompressible, thermally developing, laminar flow conditions. As a result, the microchannel with an aspect ratio value of 1 has the highest convection heat transfer coefficient and the lowest pressure drop. It is also observed that the convection heat transfer coefficient and pressure drop increases with an increase in nanoparticle volume fraction value of nanofluid. However, the Nusselt number decreases with increasing nanoparticle volume fraction, while the Darcy friction factor is not affected. Keywords: convection heat transfer, pressure drop, rectangular microchannel, ZnO-EG nanofluid Highlights • ZnO-EG nanofluid flow and convective heat transfer characteristics through rectangular microchannels are numerically investigated. • Thermal and hydrodynamic performances of rectangular microchannels are performed for different aspect ratios. • Results showed that the thermal and hydrodynamic performances decrease with an increase in the aspect ratio of the rectangular microchannel. • The Nusselt number decreases with increase in nanoparticle volume fraction of ZnO-EG nanofluid, while the Darcy friction factor is not affected. 0 INTRODUCTION Nanofluids are prepared with the addition of nano-sized metallic or non-metallic particles to conventional working fluids to enhance their heat transfer ability. The term "nanofluid" was introduced by Choi [1] in 1995. In recent decades, the use of nanofluids as a heat transfer fluid in microchannel applications which are used in electronic device cooling has been investigated. Azizi et al. [2] investigated the thermal performance of Cu-water nanofluid through a rectangular microchannel assembled into a cylindrical geometry. The Nusselt number is enhanced up to 23 % for 0.3 wt % nanoparticle addition. Rimbault et al. [3] reported that the CuO-water nanofluid flow through rectangular microchannels has the higher pressure drop by 70 % with a 4.5 % nanoparticle addition. However, the ratio of the amount of heat transferred to the pumping power required decreases with an increase in nanoparticle volume fraction. Kalteh [4] has compared different nanofluids which are a combination of nine different nanoparticle and three different base fluids. The highest and the lowest heat transfer coefficients were obtained for diamond-water and SiO2-water nanofluid, respectively. However, the obtained pressure drop values for all nanofluid type in same base fluid were almost equal. Halelfadl et al. [5] have investigated fluid flow and heat transfer characteristics of CNT-water nanofluid through the rectangular microchannel. It was found that the nanofluid use reduces the total thermal resistance. Zhang et al. [6] have studied the fluid flow and heat transfer characteristics of Al2O3-water nanofluid through circular microchannels. The Nusselt number increased with increase in nanoparticle volume fraction and the Reynolds number. The maximum heat transfer enhancement was reported to be 10.6 %, while the Darcy friction factor increased by 7.9 %. Nebbati and Kadja [7] found that yAl2O3-water nanofluid use in microchannels increased the local Nusselt number and decreased the bottom surface temperature and shear stress on the wall. Mohammed et al. [8] numerically studied the fluid flow and heat transfer characteristics of Al2O3-water nanofluid having different nanoparticle volume fractions through a rectangular microchannel. It was found that the heat transfer coefficient and wall shear stresses increase with increasing nanoparticle volume fraction, while thermal resistance decreases. Mohammed et al. [9] numerically compared the nanofluids prepared with the addition of Al2O3, Ag, CuO, diamond, SiO2 and TiO2 to water. The highest and the lowest heat transfer coefficients were obtained for the diamond-water and the Al2O3-water nanofluids, respectively. Moreover, the highest and the lowest pressure drops were obtained for SiO2-water and Ag-water nanofluids, respectively. Salman et al. [10] compared the thermal performances of Al2O3-water and SiO2-water nanofluid flows through microtubes, numerically and experimentally. They found that SiO2-water nanofluid has the highest Nusselt number value and followed by Al2O3-water nanofluid and pure water, respectively. Mohammed et al. [11] numerically compared the thermal performances of Al2O3-water, Ag-water, SiO2-water and TiO2-water nanofluids through a square microchannel heat exchanger. The highest heat transfer coefficient was obtained for Al2O3-water nanofluid. Moreover, the lowest pressure drop was obtained for Ag-water nanofluid. The thermophysical properties of working fluid are critical parameters of the thermal performance of working fluid, because the thermal performance parameters such as heat transfer rate, thermal diffusion, the Prandtl number, and the Reynolds number are related to thermophysical properties of working fluid. Therefore, the definition of thermophysical properties of nanofluids is a critical topic in the investigation of the thermal performance of nanofluids. Vajjha and Das [12] proposed a correlation to calculate the specific heat of nanofluids. In their study, Al2O3-EG/ water (60:40), ZnO-EG/water (60:40), SiO2-EG/water (60:40) nanofluids were used. They reported that the nanoparticle addition to base fluid decreased the specific heat of base fluid. Satti et al. [13] measured the specific heat of five different nanofluids containing Al2O3, ZnO, CuO, TiO2, and SiO2 nanoparticles dispersed in a base fluid of 60 % propylene glycol and 40% water by mass. It was found that the specific heat decreases with increase in nanoparticle volume fraction. At 243 K, specific heats of ZnO-PG/water (60:40) nanofluids for 0.5 % and 4.0 % nanoparticle addition reduce by about 28 % and 32 %, respectively. Esfe and Saedodin [14] and [15] proposed correlations to determine the thermal conductivity and viscosity of ZnO-EG nanofluids. Yu et al. [16] found that in low volume concentrations of ZnO-EG nanofluid showed Newtonian behavior, while at higher volume concentrations it showed shear-thinning behavior. Prajapati et al. [17] reported that the heat transfer coefficient increases with an increase in nanoparticle volume fraction of the ZnO-water nanofluid. Salman et al. [18] found that the highest Nusselt number was obtained for SiO2-EG nanofluid, followed by ZnO-EG, CuO-EG, Al2O3-EG and pure EG. EG and EG/water mixtures are used in car engines, heating systems, and solar heating installations as heat transfer fluids in cold regions. Nanofluid-based EG and EG/water mixtures can be used as alternative heat transfer fluids in these fields. Therefore, the fluid flow and heat transfer characteristics of nanofluid-based EG should clearly be investigated. From the above literature survey, it was clearly seen that the studies performed on the fluid flow and heat transfer characteristics of ZnO-EG nanofluid in microchannels are limited. In this study, the fluid flow and heat transfer characteristics of ZnO-EG nanofluid for different nanoparticle volume fractions through rectangular microchannels having different aspect ratios were numerically investigated under laminar flow conditions to investigate the use of ZnO-EG nanofluid as a heat transfer fluid in electronic device cooling and the importance of geometric configuration. The dimensionless temperatures, convective heat transfer coefficient, the Nusselt number, pressure drop, the Darcy friction factor values for ZnO-EG nanofluid were determined. 1 NUMERICAL PROCEDURE 1.1 Model Description The system is considered to be a microchannel attached over a hot surface to dissipate heat generated by any electronic device (e.g. computer chips). The bottom wall of the microchannel is exposed to heat generated by the electronic device. The schematic diagram of the rectangular microchannel is illustrated in Fig. 1. Six different aspect ratios between a is 1 and 2 are defined. The hydraulic diameter and length of the microchannel for each aspect ratio values of the microchannel are assumed 150 ^m and 5 cm, respectively. The other geometric parameters are defined by using the following equation, D _ 4A _ 2HW h _ P ~(H + W) (1) where Dh, H and W are the hydraulic diameter, height and weight of microchannel, respectively. Fig. 1. The schematic diagram of microchannel The obtained results for other geometric parameters belonging to the rectangular microchannel are shown in Table 1. Table 1. The geometric dimensions of microchannels а (H /W) H [um] W [/um] 1.0 150 150 1.2 165 137.5 1.4 180 128.57 1.6 195 121.875 1.8 210 116.67 2.0 225 112.5 Z-momentum equation: U dW+V dW+W dX dY rdW dZ = -dP _L dZ Re d W d W d2 +--^ + - dX2 dY2 dZ2 Energy equation: .30 и^+v О^ —\ = dX dY dZ, 1 ( d20 d20 d20 л Re ■ Pr dX2 dY2 dZ2 (6) 1.2 Governing Equations Fluid flow containing infinitesimal solid particles with diameters less than 100 nm can be modeled as single-phase flow [19]. Therefore, the ZnO-EG nanofluid flow is modelled as single-phase flow in this study. The following assumptions are also adopted for this study: (i) both heat transfer and fluid flow in microchannel are in three dimensional and steady-state; (ii) nanofluid flow is incompressible and laminar; (iii) the physical properties of nanofluid, such as density, specific heat, thermal conductivity are temperature independent; and (iv) the buoyancy effect, viscous dissipation and radiation heat transfer are negligible. The governing equations for the singlephase model can be written as follows under the abovementioned assumptions: Continuity equation: dU dV dW n -+ — +-= 0, (2) dX dY dZ X-momentum equation: TTdU dU , U-+ V-+ W dX dY ,dU_ dZ = -dP _L = dX Re Y-momentum equation: d 2U d 2U d 2U dX2 dY2 dZ2 (3) dV dV dV ' U — + V — + W — I = dX dY dZ dp 1 dY Re f d 2V d V d2V 2 +—t + —г dX2 dY2 dZ2 + + x=r I> (4) where X, Y, and Z are dimensionless distances; U, V, and W are the dimensionless velocity components; Re is the Reynolds number; Pr is the Prandtl number; P is the dimensionless pressure and в is the dimensionless temperature. The dimensionless parameters in the above equations can be expressed as follows: (7) (8) (9) x -, Y = У, Z z ~ D Dh ~ D, u , V v w =-, W ' U in '' Uin = Un Re = PU D и ' Pr CPß k P = в = AP pK ' T - Tin T - T (10) (11) (12) where u, v, w are the velocity components in x, y and z coordinates; p is the density; ц is the viscosity; CP is the specific heat; k is the thermal conductivity; AP is the pressure drop; Uin is the inlet velocity; Tw and Tin are the wall and inlet temperature, respectively. To solve the governing equations, boundary conditions are required. The inlet velocity is obtained by considering the Reynolds numbers. The no-slip condition is imposed at the walls of the microchannel. The fixed heat flux is applied at the bottom surface of the microchannel, and the heat is transferred from the bottom wall to the nanofluid. The side and top surfaces of the microchannel are assumed to be insulated. The pressure outlet boundary condition is used at the outlet of the microchannel. The boundary conditions used in this study are given below in detail: at the inlet: U = 1, в = 1, at the outlet: P = Pout = 0, at the solid-fluid interfaces: U = 0, 0 = -k™=-k. Ъ dn dn дв„ at the bottom surface: qw = -ks —-, dn at the side and top surfaces: qw = 0. 1.3 Thermophysical Properties The governing equations include some thermophysical properties belonging to working fluid. To solve the governing equations, these thermophysical properties should be defined. The density and specific heat of the ZnO-EG nanofluid can be determined using the following equations, which are based on conventional mixture theory, respectively: Pnf =(1 -v)pbf + VPnp , (13) (PC ) =(1 -v)(PC ) MPC ) ' (14) where ф is the nanoparticle volume fraction of the nanofluid and the subscripts nf, np and bf denote nanofluid, nanoparticle, and the base fluid, respectively. To calculate the thermal conductivity of ZnO-EG nanofluid, the following correlation proposed by Esfe and Saedodin [14], which is based on experimental data obtained for ZnO-EG nanofluid can be used; knf = (0.24859Г2'504^ - 0.7492 Ik, 'bf (15) where T is the temperature and in units of °C. The validity of this correlation is in ф = 0.0625 % to 5 % and T = 25 °C to 50 °C intervals. The viscosity of ZnO-EG nanofluid can be calculated by using following correlation proposed by Esfe and Saedodin [15] which is based on experimental data obtained for ZnO-EG nanofluid; Pnf = 0.9118Exp (5.49ф - 0.00001359T +0.0303Zn (T ) Pf, (16) where T is the temperature and in units of °C. This correlation is valid for ф = 0.25 % to 5 % and T = 25 °C to 50 °C. 1.4 Numerical Procedure To carry out the numerical computation, the finite volume method is used. The governing equations are solved along with boundary conditions. The convective terms existing in the governing equations are discretized by using the second order upwind scheme. The standard scheme is employed in the discretization of pressure. A SIMPLE algorithm is used to resolve the velocity and pressure coupling [20]. The SIMPLE-code family is based on the spatial integration of the conservation equations over finite control volumes. The Green-Gauss cell-based method is used to discretize the momentum and energy equations. The governing equations are iterated until the convergence value of 1*10-6. No convergence problems are observed during the calculations. 1.5 Grid Testing and Validation In the modelling of the rectangular microchannel, the hexahedral mesh distribution is used. Finer meshes are used in the region near the microchannel walls to enhance the accuracy of the solution. The mesh distribution of rectangular microchannel is shown in Fig. 2. 75,00 225,00 Fig. 2. The mesh distribution of rectangular microchannel The grid independence test is realized for each microchannel with different aspect ratios to invalidate the effect of grid number on a solution. Six different grid numbers ranging from 40,000 to 1,400,000 are used to model each microchannel. Optimum grid numbers are defined considering the change of result obtained by the solution. The numerical computations are realized for defined optimum models. To validate the accuracy of code, the results obtained by the present study are compared with the results of Lee et al. [21], which are based on experimental data for rectangular microchannels in thermally developing laminar flow under constant heat flux conditions. The comparison of the results obtained by the present study with that of Lee et al. [21] is illustrated in Fig. 3. As can be seen, the local Nusselt number values obtained by the present study showed good agreement with the results of Lee et al. [21]. The maximum deviation from the results of Lee at al. [21] is obtained to be 9.45 %. 0.000 0.005 0.010 0.015 0.020 0.025 0.030 X Fig. 3. Comparison of results obtained by present study with the results of Lee et al. [21] 2 RESULTS AND DISCUSSION The fluid flow and heat transfer characteristics of ZnO-EG nanofluid through rectangular microchannels are numerically investigated. In the calculations, different nanoparticle volume fractions (ф = 1 % to 4 %) of ZnO-EG nanofluid are used. The Reynolds number is in the range of 10 to 100. The calculations are realized for different aspect ratio values (a = 1, 1.2, 1.4, 1.6, 1.8, 2) of the rectangular microchannel. The fixed heat flux of 1000 W/m2 is applied at the bottom surface of the microchannels. 2.1 Dimensionless Temperature The dimensionless temperature is calculated by using Eq. (12). The variation of dimensionless temperature along with axial dimensionless distance for each aspect's ratio values of the rectangular microchannel in the use of pure EG at Re = 10 and Re = 100 are illustrated in Figs. 4a and b, respectively. The variation of the dimensionless temperature along with the axial dimensionless distance for different nanoparticle volume fractions in a = 1 at Re = 10 and Re = 100 are illustrated in Figs. 5a and b, respectively. The ZnO nanoparticle addition to pure EG decreased the bulk and wall temperature of the fluid. However, the decrease in wall temperature is higher than that of the bulk temperature. Therefore, the Fig. 4. The variation of dimensionless temperature along with axial dimensionless distance for different aspect ratios (working fluid is pure EG) at a) Re = 10 and b) Re = 100 dimensionless temperature increases with an increase in the nanoparticle volume fraction of ZnO nanofluid. 2.2 Convective Heat Transfer Coefficient The convective heat transfer coefficient can be written as follows: h _ m CP (outlet Tinlet ) " A (( - Tb ) ' (17) where m is the mass flow rate, and As is the heated surface area of the microchannel. The variation of convective heat transfer coefficient with the Reynolds number for each aspect ratio of the rectangular microchannel in the use of pure EG is illustrated in Fig. 6. As can be seen from Fig. 6, the convective heat transfer coefficient decreases with an increase in aspect ratio of the rectangular microchannel. The convective heat transfer coefficients obtained for a = 1 Fig. 5. The variation of dimensionless temperature along with axial dimensionless distance for different nanoparticle volume fractions (for а = 1) at a) Re = 10 and b) Re = 100 40 60 Re Fig. 6. The variation of convective heat transfer coefficient with the Reynolds number for different aspect ratios (working fluid is pure EG) at Re = 10 and Re = 100 is higher my 19.69 % and 7.18 % compared to that of а = 2.0, respectively. The microchannel having an aspect ratio of а = 1 provides better heat dissipation by convection heat transfer. The microchannel width is decreased with increasing aspect ratio to be able to provide constant hydraulic diameter. This tendency causes decreases in the heated surface area of the microchannel and consequently increases the convective heat transfer coefficient. However, a decrease in the outlet temperature of pure EG with increasing aspect ratio is higher than the decrease in heated surface area. Therefore, the convective heat transfer coefficient decreases with the increasing aspect ratio of the microchannel. The variation of convective heat transfer coefficient with the Reynolds number for different nanoparticle volume fractions in а = 1 is illustrated in Fig. 7. Fig. 7. The variation of convective heat transfer coefficient with the Reynolds number for different nanoparticle volume fractions (for а = 1) The convective heat transfer coefficient increases with the increase in nanoparticle volume fraction of ZnO-EG nanofluid. The increments in the convective heat transfer coefficient for 4 % ZnO-EG nanofluid compared to pure EG is 19.33 % and 16.89 % at Re = 10 and Re = 100, respectively. The ZnO nanoparticle addition to pure EG increases the density, outlet temperature of the nanofluid, and the velocity of the nanofluid at fixed Reynolds numbers. Moreover, the ZnO nanoparticle addition to pure EG decreases the temperature difference between wall and bulk temperatures of nanofluid. This fact is also reported by Tsai and Chein [22]. These facts increase the convective heat transfer coefficient of flow. The specific heat of pure EG decreases with increasing nanoparticle volume fraction. This case has a reducing effect on the convective heat transfer coefficient but is eliminated by other additive effects. The variation of local convective heat transfer coefficient along with axial dimensionless distance for each aspect ratio of the rectangular microchannel Fig. 8. The variation of local convective heat transfer coefficient along with axial dimensionless distance for different aspect ratios (working fluid is pure EG) at a) Re = 10 and b) Re = 100 Fig. 9. The variation of local convective heat transfer coefficient along with axial dimensionless distance for different nanoparticle volume fractions (a = 1) at a) Re = 10 and b) Re = 100 in the use of pure EG at Re = 10 and Re = 100 are illustrated in Figs. 8a and b, respectively. As can be seen from Figs 8a and b, the local convective heat transfer values are much higher in the inlet region of the microchannel. It is also observed that the reduction in local convective heat transfer coefficient after the inlet region of the microchannel at Re = 100 is higher compared to that of at Re = 10. The variation of local convective heat transfer coefficient along with the axial dimensionless distance for different nanoparticle volume fractions in а = 1 at Re = 10 and Re = 100 are illustrated in Figs. 9a and b, respectively. As expected, the local convective heat transfer coefficients for 4 % ZnO-EG nanofluid take higher values in comparison to other nanoparticle volume fractions of ZnO-EG nanofluid and pure EG. It is also observed that the ZnO-EG nanofluid flow is thermally developing a flow for this study. 2.3 The Nusselt Number The Nusselt number is expressed as follows: Nu = k (18) The variation of the Nusselt number with the Reynolds number for each aspect ratio of the rectangular microchannel in the use of pure EG is illustrated in Fig. 10. The Nusselt number decreases with increases in the aspect ratio of the rectangular microchannel. It is found that the Nusselt number obtained for а = 1 at Re = 10 and Re = 100 is higher by 19.69 % and 7.18 % compared to that of а=2, respectively. The variation of the Nusselt number with the Reynolds number for different nanoparticle volume fractions in а = 1 is illustrated in Fig. 11. As can be seen from Fig. 11, the Nusselt number surprisingly decreases with increases in nanoparticle volume fraction of ZnO-EG nanofluid. The Nusselt number values for 4.0 % ZnO-EG nanofluid at Re = 10 and Re = 100 are lower by 2.02 % and 4.02 % compared to that of pure EG, respectively. Fig. 10. The variation of the Nusseit number with the Reynolds number for different aspect ratios (working fluid is pure EG) Fig. 11. The variation of the Nusselt number with the Reynolds number for different nanoparticle volume fractions (for а = 1) The reason for the decrease in the Nusselt number is a higher increment in conductive heat transfer coefficient when ZnO nanoparticle is added to pure EG. The increment in convective heat transfer coefficient is obtained to be 19.33 % for 4 % ZnO-EG nanofluid at Re = 10, while the increment in conductive heat transfer coefficient is 21.78 % for 4 % ZnO-EG nanofluid. This higher increment in conductive heat transfer coefficient compared to convective heat transfer coefficient reduces the Nusselt number of the flow. 2.4 Pressure Drop The variation of normalized pressure drop with the Reynolds number for each aspect ratio of the rectangular microchannel in the use of pure EG is illustrated in Fig. 12. The normalized pressure drop increases with increases in aspect ratio of the rectangular microchannel. It is found that the normalized pressure drop for а=2 at Re = 10 and Re = 100 is 9.05 % and 8.91 % higher compared to that of а = 1, respectively. The variation of normalized pressure drop with the Reynolds number for different nanoparticle volume fractions in а = 1 is illustrated in Fig. 13, respectively. It is found that the normalized pressure drop increases with increases in nanoparticle volume fraction of ZnO-EG nanofluid in that ZnO nanoparticle addition to pure EG increases the viscosity of the nanofluid. The normalized pressure drop values for 4 % ZnO-EG nanofluid at Re = 10 and Re = 100 are 28.96 % and 29.23 % higher compared to that of pure EG, respectively. Fig. 12. The variation of normalized pressure drop with the Reynolds number for different aspect ratios (working fluid is pure EG) Fig. 13. The variation of normalized pressure drop with the Reynolds number for different nanoparticle volume fractions (for а = 1) 2.5 The Darcy Friction Factor The Darcy friction factor is determined as follows: AP f = 2 D L pU„ (19) where L is the length of the microchannel, AP is the pressure drop. The variation of the Darcy friction factor with the Reynolds number for each aspect ratio of the rectangular microchannel in the use of pure EG is illustrated in Fig. 14. Fig. 14. The variation of the Darcy friction factor with the Reynolds number for different aspect ratios (working fluid is pure EG) The Darcy friction factor increases with increases in the aspect ratio of the rectangular microchannel. The pressure drop for а=2 at Re = 10 and Re = 100 is 9.05 % and 8.91 % higher compared to that of а = 1, respectively. The variation of the Darcy friction factor with the Reynolds number for different nanoparticle volume fractions in а = 1 is illustrated in Fig. 15. Fig. 15. The variation of the Darcy friction factor with the Reynolds number for different nanoparticle volume fractions (for а = 1 ) It is observed that the ZnO nanoparticle addition to pure EG does not affect the Darcy friction factor. The same finding is reported for Al2O3-water and Al2O3-EG/water nanofluids by Jung et al. [23]. The local Darcy friction factor can be calculated as follows fx pU (20) where тх is local shear stress. The variation of the local Darcy friction factor along with axial dimensionless distance for each aspect ratio of the rectangular microchannel in the use of pure EG at Re = 10 and Re = 100 are illustrated in Figs. 16a and b, respectively. Fig. 16. The variation of the local Darcy friction factor along with axial dimensionless distance for different aspect ratios (working fluid is pure EG) at a) Re=10 and b) Re=100 As can be seen from Figs 16a and b, the flow is hydrodynamically developed. The velocity of flow and the local shear stress values increase with an increase in the Reynolds number. However, an increase in velocity is higher than that of local shear stress. Therefore, the local Darcy friction values for Re = 10 are higher than that of Re = 100. The highest local Darcy friction factor values along with the axial dimensionless distance are obtained for a=2, while the lowest values are obtained for a = 1. The surface area of a microchannel increases with increases in aspect ratio of the microchannel. A larger surface area leads to increases in local shear stress. 2.6 Thermal Resistance and Pumping Power Thermal resistance can be defined as follows: T - T Rh =- Qin (21) where Qin is the amount of heat transferred to nanofluid. Another important parameter is the pumping power required to drive the nanofluid through the microchannel. The pumping power required can be written as: W = UA AP, (22) where Ac is the cross-section area of the microchannel. The increase in nanoparticle volume fraction of ZnO-EG nanofluid reduces the wall temperature. Therefore, the increase in the nanoparticle volume fraction of ZnO-EG nanofluid caused a reduction in thermal resistance. The nanoparticle addition to working fluid also increases the pressure drop and, therefore, the pumping power increases. The nanofluids should be evaluated by their thermal resistance in a given pumping power. The variation of the thermal resistance with pumping power required for each aspect ratios of the rectangular microchannel in the use of pure EG is illustrated in Fig. 17. In a given pumping power required, the highest thermal resistance values are obtained for a=2, while the lowest values are obtained for a = 1. The variation of thermal resistance with pumping power required for different nanoparticle volume fractions of ZnO-EG nanofluid in a = 1 is illustrated in Fig. 18. 0.5 1.0 1.5 2.0 Pumping power [W] Fig. 18. The variation of thermal resistance with pumping power required for different aspect ratios (working fluid is pure EG) In a given pumping power, the highest thermal resistance values are obtained for pure EG. The thermal resistance decreases with increases in the nanoparticle volume fraction of ZnO-EG nanofluid. The lowest thermal resistance values for a required pumping power are obtained for 4 % ZnO-EG nanofluid. 0.5 1.0 1.5 2.0 Pumping power [W] Fig. 17. The variation of thermal resistance with pumping power required for different aspect ratios (working fluid is pure EG) 4 CONCLUSIONS The fluid flow and heat transfer characteristics of ZnO-EG nanofluid through rectangular microchannels with different aspect ratios are numerically investigated. The results showed that the geometrical configuration of microchannels is a very effective parameter on fluid flow and heat transfer characteristics. Therefore, the hydrodynamic and thermal performances of microchannels can be enhanced with only geometric configuration in the system design. The microchannel with a = 1 gave better results in the sense of both fluid flow and heat transfer. In the analyses, the different nanoparticle volume fractions of ZnO-EG nanofluid are used as a working fluid. The use of ZnO-EG nanofluid enhanced the convection heat transfer; however, it simultaneously increased the pressure drop. In this case, the thermal resistance value at a given pumping power provides insight about the advantage or disadvantage of nanofluid use. ZnO-EG nanofluid has an advantage over pure EG due to their lower thermal resistance values at a given pumping power compared to pure EG. 5 REFERENCES [1] Choi, S.U.S. (1995). Enhancing thermal conductivity of fluids with nanoparticles. ASME FED, vol. 231, p. 99-103. [2] Azizi, Z., Alamdari A., Malayeri, M.R. (2015). Convective heat transfer of Cu-water nanofluid in a cylindrical microchannel heat sink. Energy Conversion and Management, vol. 101, p. 515-524, D0I:10.1016/j.enconman.2015.05.073. [3] Rimbault, B., Nguyen, C.T., Galanis, N. (2014). Experimental investigation of CuO-water nanofluid flow and heat transfer inside a microchannel heat sink. International Journal of Thermal Sciences, vol. 84, p. 275-292, DOI:10.1016/j. ijthermalsci.2014.05.025. [4] Kalteh, M. (2013). Investigating the effect of various nanoparticle and base liquid types on the nanofluids heat and fluid flow in a microchannel. Applied Mathematical Modeling, vol. 37, no. 18-19, p. 8600-8609, DOI:10.1016/j. apm.2013.03.067. [5] Halelfadl, S., Adham, A.M., Mohd-Ghazali, N., Maré, T., Estellé, P., Ahmad, R. (2014). Optimization of thermal performances and pressure drop of rectangular microchannel heat sink using aqueous carbon nanotubes based nanofluid. Applied Thermal Engineering, vol. 62, no. 2, p. 492-499, DOI:10.1016/j. applthermaleng.2013.08.005. [6] Zhang, H., Shao, S., Xu, H., Tian, C. (2013). Heat transfer and flow features of Al2O3-water nanofluids flowing through a circular microchannel -Experimental results and correlations. Applied Thermal Engineering, vol. 61, no. 2, p. 86-92, DOI:10.1016/j.applthermaleng.2013.07.026. [7] Nebbati, R., Kadja, M. (2015). Study of forced convection of a nanofluid used as a heat carrier in a microchannel heat sink. Energy Procedia, vol. 74, p. 633-642, DOI:10.1016/j. egypro.2015.07.799. [8] Mohammed, H.A., Gunnasegaran, P., Shuaib, N.H. (2010). Heat transfer in rectangular microchannels heat sink using nanofluids. International Communications in Heat and Mass Transfer, vol. 37, no. 10, p. 1496-1503, DOI:10.1016/j. icheatmasstransfer.2010.08.020. [9] Mohammed, H.A., Gunnasegaran, P., Shuaib, N.H. (2011). The impact of various nanofluid types on triangular microchannels heat sink cooling performance. International Communications in Heat and Mass Transfer, vol. 38, no. 6, p. 767-773, DOI:10.1016/j.ijheatmasstransfer.2011.03.024. [10] Salman, B.H., Mohammed, H.A., Kherbeet, A.Sh. (2014). Numerical and experimental investigation of heat transfer enhancement in a microtube using nanofluids. International Communications in Heat and Mass Transfer, vol. 59, p. 88100, DOI:10.1016/j.icheatmasstransfer.2014.10.017. [11] Mohammed, H.A., Bhaskaran G., Shuaib, N.H., Abu-Mulaweh, H.I. (2011). Influence of nanofluids on parallel flow square microchannel heat exchanger performance. International Communications in Heat and Mass Transfer, vol. 38, no. 1, p. 1-9, DOI:10.1016/j.icheatmasstransfer.2010.09.007. [12] Vajjha, R.S., Das, D.K. (2009). Specific heat measurement of three nanofluids and development of new correlations. ASME Journal of Heat Transfer, vol. 131, no. 7, 071601, p. 1-7, DOI:10.1115/1.3090813. [13] Satti, J.R., Das, D.K., Ray, D. (2016). Specific heat measurements of five different propylene glycol based nanofluids and development of a new correlation. International Journal of Heat and Mass Transfer, vol. 94, p. 343-353, DOI:10.1016/j.ijheatmasstransfer.2015.11.065. [14] Esfe, M.H., Saedodin, S. (2014). Experimental investigation and proposed correlations for temperature-dependent thermal conductivity enhancement of ethylene glycol based nanofluid containing ZnO nanoparticles. Journal of Heat and Mass Transfer Research, vol. 1, no. 1, p. 47-54. [15] Esfe, M.H., Saedodin, S. (2014). An experimental investigation and new correlation of viscosity of ZnO-EG nanofluid at various temperatures and different solid volume fractions. Experimental Thermal and Fluid Sciences, vol. 55, p.1-5, DOI:10.1016/j.expthermflusci.2014.02.011. [16] Yu, W., Xie, H., Chen, L., Li, Y. (2009). Investigation of thermal conductivity and viscosity of ethylene glycol based ZnO nanofluid. Thermochimica Acta, vol. 491, no. 1-2, p. 92-96, DOI:10.1016/j.tca.2009.03.007. [17] Prajapati, O.S., Rohatgi, N., Rajvanshi, A.K. (2013). Heat transfer behavior of nano-fluid at high pressure. Journal of Materials Science and Surface Engineering, vol. 1, no. 1, p. 1-3. [18] Salman, B.H., Mohammed, H.A., Munisamy, K.M., Kherbeet, A.Sh. (2014). Three-dimensional numerical investigation of nanofluids flow in microtube with different values of heat flux. Heat Transfer-Asian Research, vol. 44, no. 7, p. 599-619, DOI:10.1002/htj.21139. [19] Das, S.K., Choi, S.U.S., Yu, W., Pradeep, T. (2008). Nanofluids: Science and Technology. John Wiley & Sons, Jersey. [20] Patankar, S.V. (1980). Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, Washington. [21] Lee, P.S., Garimella, S.V., Liu, D. (2005). Investigation of heat transfer in rectangular microchannels. International Journal of Heat and Mass Transfer, vol. 48, no. 9, p. 1688-1704, DOI:10.1016/j.ijheatmasstransfer.2004.11.019. [22] Tsai, T.H., Chein, R. (2007). Performance analysis of nanofluid-cooled microchannel heat sinks. International Journal of Heat and Fluid Flow, vol. 28, no. 5, p. 1013-1026, DOI:10.1016/j. ijheatfluidflow.2007.01.007. [23] Jung, J.Y., Oh, H.S., Kwak, H.Y. (2009). Forced convective heat transfer of nanofluids in microchannels. International Journal of Heat and Mass Transfer, vol. 52, no. 1-2, p. 466-472, DOI:10.1016/j.ijheatmasstransfer.2008.03.033. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)10, 614-622 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3535 Original Scientific Paper Numerical Simulation of Transient Pressure Control in a Pumped Water Supply System Using an Improved Bypass Pipe Xiaozhou Li- Manlin Zhu* - Jiancang Xie Xi'an University of Technology, State Key Laboratory Base of Eco-hydraulic Engineering in Arid Area, China The use of a bypass pipe in a pump system is an economical and simple method of preventing water hammer surges to control the extreme transient pressures induced by pump failure. Currently, conventional bypass pipes (CBPs) have certain limitations in engineering applications. To overcome these limitations, in this study, a CBP is modified to create an improved bypass pipe (IBP). Unlike that of a CBP, the control valve in an IBP system is a hydraulic control valve that uses a controller on an uninterruptible power source (UPS) to turn the valve on or off. This modification enables the hydraulic control valve to be precisely controlled and facilitates the prevention of water hammer surges. To study the transient control effect of the IBP using the method of characteristics (MOC), a mathematical model of a complex system that includes an IBP, a check valve, and a water pump is established. Using MATLAB software, a water supply system in an industrial zone is used as an example for numerical simulations of the transient process experienced by this pump after it fails. The extreme pressure variations and the envelope of the pipe pressure in a typical section of the pumped water supply system are analysed and compared. The results show that an IBP can effectively reduce the maximum extreme pressure in a pumped water supply system. Keywords: improved bypass pipe, hydraulic transient, pumped water supply system, UPS controller, numerical simulation Highlights • Based on a conventional bypass pipe, a new, economical, and effective water hammer prevention measure, i.e., an improved bypass pipe, is proposed. • A hydraulic control valve with an on/off status determined by an uninterruptible power source controller is installed in the bypass pipe to control the transient pressure. • The on/off time for the newly developed improved bypass pipe is accurately controlled to facilitate transient pressure control. • Based on the method of characteristics, a mathematical model for a complex system containing an improved bypass pipe is established. 0 BACKGROUND Pump failures induce rapid changes in the flow rate, which causes extremely high-pressure water hammer surges in the water supply system, including extremely high- and low-pressure water hammer surges. Extreme water hammer pressure damages pipes and hydraulic equipment and results in water supply service interruptions. To eliminate or reduce the impact of extreme water hammer pressure on water delivery systems, water hammer prevention measures are normally implemented; these include bypass pipes, air valves, check valves, air tanks, pressure regulating chambers, unidirectional pressure regulating chambers, pressure discharge valves, water hammer prevention valves, etc. [1] to [3]. An air valve is a common, economical, and effective piece of water hammer prevention equipment that is normally installed at the local peak of a pipe. When negative pressure occurs inside the pipe, the air valve opens to allow air in; when the pressure inside the pipe increases, the air valve releases it by discharging air [4] and [5]. A check valve is an economical and reasonable water hammer prevention measure that is widely used in pumping stations [6]. Check valves can prevent water pumps from rotating in reverse due to liquid backflow and thus prevent damage to the motor. However, instant check valve closure results in a catastrophically high pressure [7] and [8]. Therefore, when a check valve is installed in a pumping station to control transients, an air tank should also be installed near the pumping station to prevent the increase in pressure that results from closing the check valve [9]. However, air tanks have drawbacks that include their large volume, high maintenance cost and complex operation (the gas in the tank should be replenished). A pressure discharge valve with the proper parameters can prevent a water delivery system from generating excessively high or low pressures; otherwise, the negative impact of water hammer surges on the water delivery system is exacerbated [10]. To ensure proper protection of a water delivery system by a pressure discharge valve, the pressure discharge valve should have so little inertia that it can be opened promptly in response to a rapid change in the pressure to prevent delays in its opening [11]. A water hammer 614 *Corr. Author's Address: Xi'an University of Technology, NO.5 South Jinhua Road, Xi'an 710048, Shaanxi ,China, zhuml@xaut.edu.cn prevention valve may have inaccurate hydraulic control that causes its opening to be delayed, or that causes it to close excessively fast, which results in water hammer damage. Wylie [12] proposed two types of conventional bypass pipes (CBPs) for preventing water hammer surges in pumping stations. One was to install a check valve in the water pump's outlet and a control valve in the bypass pipe to prevent the pump from rotating in reverse due to water backflow. After pump failure, the bypass pipe control valve opens instantly and then closes gradually to release the high pressure inside the system and eliminate water hammer surges. The other type was designed to prevent low pressure and water column separation in the water pump's outlet from being caused by pump failure. A control valve was installed at the water pump's outlet, and a check valve was installed in the bypass pipe; therefore, the water flowing into the reservoir could be directed to the water pump's outlet through the bypass pipe. Both types of CBPs have certain drawbacks. The first type of bypass pipe cannot accurately control the valve's on/off timing; therefore, it may not be able to prevent water hammer surges as well as is desired. The second type of bypass pipe cannot release the high pressure in the system and is only applicable to pumping station water supply systems that ignore the effect of backflow and have water pump inlets with positive water head pressures. Compared with experimental studies, numerical simulations take less time and are more economical. With advances in computer technology, numerical simulations have become widely employed in water hammer analyses [13] to [15]. To predict pressure variations in a water supply system and to choose appropriate water hammer preventive measures, numerical simulations of water hammer surges are indispensable, simple, and effective [16]. Methods of numerically analysing water hammer surges include arithmetic, graphical, algebraic and linear analysis methods as well as the method of characteristics (MOC). At present, the MOC is the simplest and most popular method for numerically solving transient flow problems. Izquierdo proposed mathematical models for simulating transient boundary conditions in simple and complex hydraulic systems [17] and [18]. Because the results of hydraulic transient numerical simulations are in good agreement with those of test measurements [19] and [20], the results of numerical simulations of water hammer surges are widely accepted. A composite bypass pipe and check valve is installed in a pumping station to reduce the extremely high water pressure associated with water hammer surges in a pumped water supply system. However, this type of CBP has limitations. Therefore, this study has three purposes: 1) To overcome the limitations of CBPs, an improved bypass pipe (IBP)-based water hammer prevention measure is proposed, and the operational principle of an IBP is explained. 2) Based on previous studies and the method of characteristics (MOC), a mathematical model of the boundary conditions of an IBP is established, and a solution for this mathematical model is described. 3) A water supply system in an industrial zone is used as an example in a numerical simulation of the water hammer prevention effect of a pump with an IBP-based protective device on a water supply system and to provide evidence for the use of an IBP-based water hammer prevention measure in a system. 1 THE STRUCTURE AND OPERATIONAL PRINCIPLE OF THE IMPROVED BYPASS PIPE Fig. 1 shows a diagram of the IBP. It is based on a CBP, but the control valve in the CBP is replaced by a hydraulic control valve. The hydraulic control valve consists of a valve body, a hydraulic pressure system, an oil cylinder, an energy storage tank, and an electromagnetic valve. The IBP's hydraulic control valve leverages an uninterruptible power source (UPS), which provides power to the electromagnetic valve, which controls when and for how long the main valve is open. This improvement can prevent the valve from turning on or off in response to man-made or hydraulic interference and provide more accurate control over the valve's status than a conventional manually or hydraulically controlled valve can. Hydraulic control valve Fig. 1. Schematic of the IBP After the pump fails, the UPS-powered controller can predefine a delayed opening time, t1 (the interval between when the pump fails and when the electromagnetic valve powers on), for the electromagnetic valve according to the requirements. After the electromagnetic valve powers on and opens, oil from the energy storage tank enters the rod-less chamber of the oil cylinder and rotates the valve board to open the main valve. The main valve's opening time, t2 (the time it takes the main valve to make the transition from being completely closed to being completely open), is regulated via the pressure at the valve that controls the oil's flow rate according to actual requirements. After the time allotted for the electromagnetic valve to be open expires, the UPS-powered controller closes the electromagnetic valve, and the oil in the energy storage tank enters the rod chamber of the oil cylinder and pushes the piston back to close the main valve. The main valve's closing time, t4 (the time it takes for the main valve to make the transition from completely open to completely closed), is regulated by a valve closure regulation valve according to actual requirements. The difference between the time for which the electromagnetic valve is open and the time required for the main valve to open is defined as the main valve's opening duration, t3 (the time the main valve is completely open). Fig. 2 shows the IBP in operation. Fig. 2. Flowchart of the IBP's operation 2 CHARACTERISTIC EQUATIONS FOR CALCULATING WATER HAMMER SURGES AND THEIR SOLUTION 2.1 Basic Differential Equations for Water Hammer Surges equations comprise a pair of non-linear hyperbolic partial differential equations, which is very difficult to solve using conventional methods. Therefore, the MOC is used to convert the two partial differential equations that take pipeline hydraulic friction into account into ordinary differential equations of a specific form. These are called characteristic equations and are given in Eqs. (3) and (4) for C+ and in Eqs. (5) and (6) for C-: g dh a dt dv fv|v + — + = о ' 2D dx — = a, dt _ gdh + dv + fvJH = 0 a dt dt 2D dx — = -a, dt (3) (4) (5) (6) where, because vsinO is typically much smaller compared with any other item in the equation for the transient process and because wave speed a is far larger than the flow velocity v, v sinö is ignored in Eq. (3) and Eq. (5), and the flow velocity, v, is ignored in Eq. (4) and Eq. (6). Next, the characteristic equations are integrated to generate a set of simplified finite-difference equations, Eqs. (7) to (10). C + : hiJAl = C p - BQiJAl, (7) U, jòt' (8) CP = hi-l,(j-1)At + Щ-Нj-1)At — Щ-Uj-1)At |Qi-l,(j-1)At I' (9) C - : hiJAt = CM + BQU 1,( j-1)At ■'M ~ "i+1,(j-1)At ' C„ ^nBQi+i,(j-i)M RQi+i,(j-y)Al\Qi+i,(j-1)At ■ (10) The basic differential equations for water hammer surges consist of an equation of motion Eq. (1), and a continuity equation, Eq. (2). dh 1 dv v dv f vivi — +--+--+ -—— = 0, (1) dx g dt g dx D 2 g dh .dh „ a2 dv _ — + v(--sm0) +--= 0. (2) dt dx g dx 2.2 Simplified Finite-Difference Equations The numerical calculation of water hammer surges is based on comprehensive and complete basic differential equations, Eqs. (1) and (2). These two In Eqs. (7) to (10) B = a/(gA), R = fAx/(2gDA2), ^г-1,(/-1)ЛЬ ^г+1,(/-1)ДЬ Qi-1,(j-l)At and Qi+lXj-X)^ and are known from the previous step; only and Q^ are unknown. 2.3 Procedure for Solving the Water Hammer Equations Fig. 3 shows the grid diagram used in the MOC. The solution process normally starts with a stable flow at t = 0. Therefore, the initial values of h and Q in each calculation section are known. Based on their values at t=0 in all the sections of a layer, the values of h and Q at t=Лt in all the sections of the layer are calculated. The boundary section is determined by the boundary condition. Next, the values at t = 2Лt in all the sections of the layer are calculated. The same procedure is followed until values at the required time, t =jAt, are calculated. Fig. 3. A characteristic line in space (x) - time (t) 3 BOUNDARY CONDITIONS AND A MATHEMATICAL MODEL FOR THE IBP Fig. 4 shows that the check valve and the controllable bypass pipe composite are installed at the pump's outlet to prevent water hammer surges. The check valve is installed at the water pump's outlet to prevent backflow in the water delivery pipe due to accidental power failure and to facilitate maintenance and repairs. However, a check valve by itself may close suddenly and cause damage if there is a water hammer surge resulting from a pump failure or if there is a severe engineering accident. Therefore, installing a controllable bypass pipe can prevent an increase in pressure resulting from a sudden closure of the check valve. To study the water hammer prevention effect of a controllable bypass pipe, the MOC is employed to perform a numerical analysis and calculations. The MOC is a widely used numerical method that is verified in reference [12]. Before establishing a mathematical model, the following assumptions are made: 1. From when the pump fails to when the flow rate at the water pump's outlet drops to zero, the check valve local coefficient of resistant is a constant, кеъ and when the flow rate at the water pump's outlet is Q6 < 0, the check valve closes instantly. 2. The local water head losses at T-type nodes upstream and downstream of the water pump are ignored. 3. The T-type nodes connecting the pipe and the bypass pipe upstream and downstream of the water pump (the dashed box in Fig. 4) are very short, and loss during the flow is ignored; they are analysed as an integral unit. When a power failure occurs, the water pump in the water supply system suddenly loses power, the pump's outlet pressure drops, and water flows back. When the flow rate at the water pump's outlet is less than or equal to zero (Q^q-^ = 0), the check valve closes instantly (Q6, jAt = 0). When the pump fails at time T < t j, the IBP closes (Q4JAt = Q 5, jAt = 0). When the pump fails at time tj < T < ti+t2+t3+t4 , the IBP opens quickly and stays open. At this point, the backflow water column is discharged from the IBP to the reservoir. When the pump fails at time T > /j+t2+t3+t4 , the IBP closes slowly (Q4, ]At = Q 5, ]At = 0). Fig. 4. A boundary model of the IBP Continuity equation: Q2, jAt = Q3, jAt + Q4, jAt, Q7, jAt = Q5, j At + Q6, j At, Q3jAt = Q6,jAt =UQR > Q4 , jAt jAt ' water balance equation: (11) h3,jAt + HR ( +U2) a0 (n + tan 1 —) + b0 a Q3, jAt Q3, jAt I , -kci ~ л 2- = h 2 gA1 6, jAt ' Q4,jAt |Q4,jAt j , h4, jAt — kC2 ^ 71 = h5, jAt 2gA2 h2, jAt = K, jAt h4, jAt, h5, jAt = h6, jAt =^7, jAt, characteristic equation: ^2, jAt = CP B1-2Q; _ CM + B7-8Q7, 2, j An "7, jAt ^M ' 7-8b^7, jAt (12) (13) water pump's unit inertia equation: (a2 +u2) al (n + tan 1 —) + b a ßo - WRA Nr n . R (a0-a) = 0. g Tr 15 A? Eqs. (11) to (14) are consolidated into (14) Fl = CP - CM - (( + B7-g )(u + Qj ) + +HR (a2 +u2 ) a0 (n + tan 1 —) + b0 a , Qru\u\ -kc = 0, 1 2gA F2 = Cp - CM -(( + B-8){Qru + Q4JAt)- (15) Q4,jAt Q4,jAt\ „ -kc2-!-7—1 = 0, 2 2£Л2 (16) F 3 = (а2 +u2) + ß _ WR2 N n +Po a1 (n + tan 1 —) + b1 а g Tr 15 At (а0 _а) = 0. (17) Eq. (15), Eq. (16) and (17) can be solved using the Newton-Raphson method [12]. = -(Bì_2 + B7_g)QR du +Hr \2u a01 n + tan 1 — I + b0 + a0a !> - kcì QrU (18) F1a~ = h\2a da a0 j n + tan — j + b0 - aov\, (19) FL dF1 = -( B1-2 + B7-s)> 4, jAt F 2v=dF2 = -( B-2 + B7-S)Qr , du F 2 = F = о, F 2 Qaj dF 2 dQT da = -(B1-2 + B7_8) - kc2 (20) (21) (22) (23) 4, jAt ГГЛ 3F3 F 3Ц =-= 2u dv al j n + tan 1 — j + èj gA + a1a, (24) гтл dF3 F 3a=^ = 2a da al j n + tan 1 — j + bl F 3 a,; ERI ^ g Tr 15 At ' dF 3 = 0. (25) dQ (26) The Newton-Raphson equations are represented as matrices: ( Flv Fla Flo JAt F 2u F 2a F 2o °4, JA F 3u F 3a F 3o °4, JAt (Au ^ (-F1 Л Aa = -F 2 l AQ4J A, J -F 3 j (27) The solutions of the above equations are represented as matrices, ^ Au Ì Aa \AQ4,jAt , (Fl F2 F3n + Fl F2n F3 V u a Q4,jAt a QiJAt u +F1Qj a,F 2u F 3a- F1&J a,F 2a F 3u - F1uF 2в^ A,F 3a- F1a F 2uF 3&J A, ) (F2a F- F2 j F3a F1& jlt, F3a - F1a F3Qja, Fl F2„ - FL F 2 ^ a Qa F^A,F3u- F2UF^A, F 2, F 3 - F 2 F 3,, F1„F 3n - F1n F 3,, F1n F 2- F1„F 2n ^JAt u u n4,jAt u a a u n4,jAt n4, jAt F1 F3 - Fl F3„ a u u a F1., F2„ - F1a F2u , u a -F1 -F2 y-F 3 У (28) Initial values are assigned to u and a as follows: = 2U( j-1)At — U( j -2)At, a = 2a( j—1)At — a( j—2)At, Q4, jAt = 2Q4,( j—1)At — Q4,( j—2)At ' (29) Eqs. (15) to (26) are substituted into Eq. (28) and solved for До, Да and AQ4,jM . Then, Eq. (29) is rewritten as: u = u + Au а=а + Aа, ß,,^ = QAJht + AQ4JÄ,. (30) The above procedure is repeated until Au, Да and AQ4,jAt are within the allowable deviation range: |Au| + | Да| + | AQ4,]At\ < e. The allowable deviation e is set to approximately 0.0002 [12]. 4 USE OF AN IBP IN A PUMPED WATER SUPPLY SYSTEM 4.1 System Overview coefficientCd , and its local resistance coefficient, kc, is shown in Eq. (31). kc = 1/Cd2 . (31) As shown in Fig. 5, a pumping station water supply system consists of three identical water pumps operating in parallel, two parallel pipes, and upstream and downstream reservoirs. The pumping station's static head is 142 m, the water pump's rated head is HR = 149.5 m, the rated flow rate is QR = 0.37 m3/s and the water pump's rated rotation speed is Nr = 1450 rpm. The pump's characteristic curve is based on data for the specific rotation speed ns=25 (SI), according to reference [12]. Each pump unit's extreme moment of inertia is WR2 = 450Nm2. The pipe length is L = 1390 m, and the pipe's inner diameter is D = 781 mm. The wave speed is a = 1213 m/s, and the pipe's coefficient of friction is f= 0.02. When the pump stops and backflow starts, the check valve behind the pump closes immediately, and the bypass pipe control valve opens or closes according to the predefined procedure. According to Appendix B of reference [12], when the check valve is completely open (i.e., the valve's opening level is т = 1), its discharge coefficient is Cd = 1.4, and when the check valve is completely closed (i.e., the valve's opening level is т=0), its discharge coefficient is Cd=0. Fig. 6 shows the relationship between the hydraulic control valve's discharge coefficient and opening level. The relationship between the valve's discharge Pump Check valve Upstream -U-w-i O reservoir t' iBP Pipeline Pipeline Downstream reservoir Fig. 5. A diagram of a pumped water supply system 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 Opening ratio Fig. 6. The hydraulic valve's discharge coefficient 4.2 Numerical Simulations Numerical simulations and calculations for two postpump failure scenarios, without an IBP and with an IBP, are performed for the pumped water supply system shown in Fig. 5. 4.2.1 Case 1: Without an IBP When all three parallel water pumps fail due to an accident and backflow occurs, the check valves behind the pumps immediately close, and the IBP is not involved. The transient process that occurs in the pumped water supply system in this situation is numerically simulated. Fig. 7 shows the flow rate at each water pump's outlet (Q6) and the typical section pressure (h6). Fig. 7. Variations in the flow rate and pressure without an IBP Fig. 7 shows that when the three parallel water pumps fail simultaneously after an accident, the flow rate at each water pump's outlet (the pipe between the water pump and the main pipe) drops rapidly until T = 1.99 s, when it becomes negative (i.e., reaches its minimum value), and the typical section pressure at the water head also drops to its minimum value. At this point, the check valve closes instantly to prevent a high volume of backflow, and the flow rate at each water pump's outlet stabilizes at zero. However, the closure of the check valve also stops the transient high-pressure wave from returning and being discharged, which may lead to an extremely high-pressure water hammer surge in the system. 4.2.2 Case 2: With an IBP Fig. 8 shows that the hydraulic control valve's on/ off program opens the IBP 2.5 s after pump failure (t1 = 2.5 s). It takes 3 s (t2 = 3.0 s) to make the transition from being completely closed to being completely open; it remains completely open for 5 s (t3 = 5 s), and then it takes 15 s (t4 = 15 s) to make the transition from being completely open to being completely closed. 25 Time [s] Fig. 8. Operation of the hydraulic valve Fig. 9 shows the flow rate at the outlet of a water pump with an IBP and typical variations in the section pressure. Fig. 9 shows that when the three parallel water pumps fail simultaneously, the flow rate at each water pump's outlet (Q6) also drops rapidly until T = 1.99 s, when it becomes negative (i.e., reaches its minimum value), and the typical section pressure at the water head drops to its minimum value as well. At this point, the check valve closes instantly to prevent a high volume of backflow, and the flow rate at each water pump's outlet stabilizes at zero. Because the IBP opens at a predefined time, the returned transient high-pressure wave is not blocked by the closed check valve. Instead, it is discharged by the IBP, eliminating the extremely high-pressure water hammer surge from the system. During this period, the maximum flow rate in the IBP reaches Q5 = -1.062 m3/s. In addition, the IBP gradually closes at a predefined time to prevent valve closure-induced water hammer surges. 5 10 15 20 25 30 35 40 45 50 Time [s] Fig. 9. The flow rate and pressure variation at the IBP 4.3 Comparison and Analysis of Simulation Results The mathematical model of the IBP presented in this study is used in numerical simulations and calculations of a transient process in a pumped water supply system. The typical variations in the section pressure determined by the simulations are shown in Fig. 10. In the graph, the red dashed line represents typical variations in the section pressure when there is no IBP, and the solid blue line represents typical variations in the section pressure when there is an IBP. The simulation results show that without an IBP, the water head of the maximum transient pressure in a typical section reaches 223.2 m, which is 1.51 times the initial pressure water head (147.5 m). With an IBP, the water head of the maximum transient pressure in a typical section drops to 171.8 m, which is 1.16 times the initial pressure water head (147.5 m). These results show that the presence of an IBP can significantly reduce the water head of a typical section's maximum transient pressure. 0 5 10 15 20 25 30 35 40 45 50 Time [s] Fig. 10. Pressure variations in a typical section Fig. 11. Envelopes of extreme pressure Fig. 11 shows the envelope of the extreme pipe pressure during the transient process in a pumped water supply system. In the graph, the blue dashed-dotted line represents the water supply pipe's central line, the blue solid line represents the water supply system's initial hydraulic grade line, the red dotted line represents the envelope of the extreme pressure during the transient process in a water supply system without an IBP, and the green dashed line represents the envelope of the extreme pressure during the transient process in a water supply system with an IBP. This graph shows that an IBP significantly reduces the maximum extreme pressure in a pumped water supply system, which ensures that the water supply system operates safely and reliably. 5 CONCLUSIONS AND FUTURE WORK In a pumped water supply system, pump failure reduces the pressure in the water supply pipe and results in subsequent liquid backflow in the pipe. An enormous backflow not only results in water column separation but also impacts the water pump's impeller and makes it rotate in reverse. In severe cases, it threatens the pumped water supply system's safe operation. Check valves are widely used in pumped water supply systems to prevent pump failure-induced water backflow and reverse rotation of the water pump. However, rapid check valve closure also results in extremely high-pressure water hammer surges. Previous researchers combined a CBP with a check valve to release the high pressure in the system and eliminate water hammer surges. However, because the timing of the CBP valve's on/off setting cannot be accurately controlled, it may not be possible to achieve the desired water hammer prevention effect. To overcome this limitation of CBPs, this study improved the CBP and proposed a new water hammer prevention measure, the IBP. An IBP is a bypass pipe with a hydraulic control valve controller whose on/ off procedure is predefined by a UPS. Therefore, the valve's setting can be accurately controlled. Using the MOC, a mathematical model of the IBP is established, and calculations are performed. This mathematical model simulates and calculates typical variations in the transient section pressure and the distribution of the pressure in the water supply pipe for a pumped water supply system in an industrial zone. A comparison of the simulation results shows that the IBP (whose hydraulic control valve's setting is changed at t1 = 2.5 s, t2 = 3.0 s, t3 = 5 s and t4 = 15 s) water hammer prevention measure in a pumped water supply system can significantly reduce the maximum extreme pressure of water hammer surges in the pumped water supply system. However, some of the problems encountered in this study require further research: 1) The mathematical model of the IBP does not take into account the bypass pipe's length and diameter; it only considers the bypass pipe valve's local coefficient of resistance, so the next step is to incorporate the bypass pipe's length, l, and diameter, d, into the mathematical model to study the effect of the bypass pipe's length and diameter on water hammer prevention. 2) The numerical simulation and calculation of the pumped water supply system's transient process shows that an IBP can protect the water supply system from damage due to water hammer surges; however, the sensitivity of key parameters affecting the IBP's water hammer prevention effect (the bypass pipe's length, l, and diameter, d, and the bypass pipe hydraulic control valve's delayed opening time, t1, opening time, t2, opening duration, t3 and closing time, t4) will analyzed to provide reference parameters for engineering applications. 6 ACKNOWLEDGMENTS This research program is sponsored by the National Natural Science Foundation of China (Grant Nos. 41471451 and 51479160). 7 NOMENCLATURE t time [s] T time [s] h pressure head [m] v flow velocity [m/s] g acceleration of gravity [m/s2] a wave speed of water hammer surge [m/s] X distance along pipe [m] в pipe slope f Darcy-Weisbach friction factor D pipe inner diameter [m] A area of pipe section [m2] C+ name for forward characteristic line C- name for reverse characteristic line CP known constant in MOC CM known constant in MOC B known constant in MOC R known constant in MOC At spatial step size [m] At time step size [s] i subscript denoting the position of the section j number of time segments N number of pipe segments Aj area of check valve section [m2] A2 area of hydraulic valve section [m2] kc1 minor loss coefficient of check valve kc2 minor loss coefficient of hydraulic valve Qr rated discharge of pump [m3/s] HR rated head of pump [m] v dimensionless pump discharge a dimensionless pump speed ratio a0, b0 parameters for pump head characteristics a1, b1 parameters for pump torque characteristics ß dimensionless torque ratio WR2 combined polar moment of inertia of pump Nr rated rotation speed of pump [rpm] TR rated torque of pump [Nm] n constant F function symbol e assigned precision control t valve opening ratio L pipe length [m] ns specific speed Cd discharge coefficient of valve l bypass pipe length [m] d bypass pipe inner diameter [m] 8 REFERENCES [1] Boulos, P.F., Karney, B.W., Wood, D.J., Lingireddy, S. (2005). Hydraulic transient guidelines for protecting water distribution systems. Journal American Water Works Association, vol. 97, no. 5, p. 111-124. [2] Bergant, A., Kruisbrink, A., Arregui, F. (2012). Dynamic behaviour of air valves in a large-scale pipeline apparatus. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 4, p. 225-237, D0l:10.5545/sv-jme.2011.032. [3] Ramezani, L. (2015). An exploration of transient protection of pressurized pipelines using air valves, PhD thesis, University of Toronto, Toronto. [4] Karadzic, U., Bulatović, V., Bergant, A. (2014). Valve-induced water hammer and column separation in a pipeline apparatus. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 11, p. 742-754, D0I:10.5545/sv-jme.2014.1882. [5] Carlos, M., Arregui, F.J., Cabrera, E., Palau, C.V. (2011). Understanding air release through air valves. Journal of Hydraulic Engineering, vol. 137, no. 4, p. 461-469, DOI:10.1061/(ASCE)HY.1943-7900.0000324. [6] Thorley, A.R.D. (2004). Fluid Transients in Pipeline Systems: A Guide to the Control and Suppression of Fluid Transients in Liquids in Closed Conduits. ASME Press, New York. [7] Sugiyama, H., Mizobata, K., Ohtani, K., Ohishi, T., Sasaki, Y., Musha, H., Miwa, T. (2003). Dynamic characteristics of a check valve for drain pumps. ASME/JSME 4 Joint Fluids Summer Engineering Conference, p. 2851-2856, D0I:10.1115/ FEDSM2003-45256. [8] Thorley, A.R.D. (1989). Check valve behavior under transient flow conditions: A state-of-the-art review. Journal of Fluids Engineering, vol. 111, no. 2, p. 178-183, D0I:10.1115/1.3243620. [9] Purcell, P. (1997). Case study of check-valve slam in rising main protected by air vessel. Journal of Hydraulic Engineering, vol. 123, no. 12, p. 1166-1168, D0I:10.1061/(ASCE)0733-9429(1997)123:12(1166). [10] Zhang, K.Q., Karney, B.W., McPherson, D.L. (2008). Pressure-relief valve selection and transient pressure control. Journal American Water Works Association, vol. 100, no. 8, p. 62-69. [11] Chaudhry, M.H. (1987). Applied Hydraulic Transients, 2nd ed. Van Nostrand Reinhold, New York. [12] Wylie, E.B., Streeter, V.L., Suo, L. (1993). Fluid Transients in Systems. Prentice Hall, New Jersey. [13] Kaliatka, A., Vaišnoras, M., Valinčius, M. (2014). Modelling of valve induced water hammer phenomeno in a district heating system. Computers&Fluids, vol. 94, p. 30-36, D0I:10.1016/j. compfluid.2014.01.035. [14] Skulovich, O., Perelman, L., Ostfeld, A. (2014). Bi-level optimization of closed surge tanks placement and sizing in water distribution system subjected to transient events. Procedia Engineering, vol. 89, p. 1329-1335, D0I:10.1016/j. proeng.2014.11.449. [15] Wan, W., Huang, W., Li, C. (2014). Sensitivity analysis for the resistance on the performance of a pressure vessel for water hammer protection. Journal of Pressure Vessel Technology, vol. 136, no. 1, p. 011303, D0I:10.1115/1.4025829. [16] Wood, D.J. (2005). Waterhammer analysis-essential and easy (and efficient). Journal of Environmental Engineering, vol. 131, no. 8, p. 1123-1131, D0I:10.1061/(ASCE)0733-9372(2005)131:8(1123). [17] Izquierdo, J., Iglesias, P.L. (2004). Mathematical modelling of hydraulic transients in complex systems. Mathematical and Computer Modelling, vol. 39, no. 4-5, p. 529-540, D0I:10.1016/S0895-7177(04)90524-9. [18] Izquierdo, J., Iglesias, P.L. (2002). Mathematical modelling of hydraulic transients in simple systems. Mathematical and Computer Modelling, vol. 35, no. 7-8, p. 801-812, D0I:10.1016/S0895-7177(02)00051-1. [19] Kendir, T.E., Ozdamar, A. (2013). Numerical and experimental investigation of optimum surge tank forms in hydroelectric power plants. Renewable Energy, vol. 60, p. 323-331, D0I:10.1016/j.renene.2013.05.016. [20] Kwon, H., Lee, J. (2008). Computer and experimental models of transient flow in a pipe involving backflow preventers. Journal of Hydraulic Engineering, vol. 134, no. 4, p. 426-434, D0I:10.1061/(ASCE)0733-9429(2008)134:4(426). Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 62, (2016), številka 10 Ljubljana, oktober 2016 ISSN 0039-2480 Izhaja mesečno Razširjeni povzetki (extended abstracts) Afrasiab Raisi: Naravna konvekcija pri nenewtonskih fluidih v pravokotni votlini z lokaliziranim virom toplote SI 97 Zhenxing Kong, Dawei Pi, Xianhui Wang, Hongliang Wang, Shan Chen: Zasnova in vrednotenje hierarhičnega algoritma za krmiljenje sistema električnega aktivnega stabilizatorja SI 98 Yongxiang Li, Xifan Yao, Jifeng Zhou: Večciljna optimizacija sestave storitev proizvodnje v oblaku po metodi genetskega algoritma, izboljšanega z entropijo oblaka SI 99 Rade Vasiljević, Milomir Gašić, Mile Savković: Vplivni parametri dinamičnega vedenja nosilne konstrukcije portalnega žerjava tipa H SI 100 Cuneyt Uysal, Kamil Arslan, Huseyin Kurt: Numerična analiza lastnosti toka fluida in prenosa toplote v pravokotnih mikrokanalih za nanofluid ZnO-etilenglikol SI 101 Xiaozhou Li, Manlin Zhu, Jiancang Xie: Numerična simulacija regulacije prehodnega tlaka v črpalnem vodovodnem sistemu z izboljšano obvodno cevjo SI 102 Osebne vesti SI 103 Naravna konvekcija pri nenewtonskih fluidih v pravokotni votlini z lokaliziranim virom toplote Afrasiab Raisi Univerza Shahrekord, Tehniška fakulteta, Iran Naravna konvekcija pogosto šteje za glavni mehanizem prenosa toplote v votlinah in ima množico aplikacij v različnih tehničnih sistemih. Fluidi v votlinah so lahko newtonski, pri mnogih naravnih in umetnih sistemih pa izkazujejo nenewtonsko vedenje. Tehnično pomembni fluidi, kot so nanofluidi, raztaljeni polimeri, barve, prehranski izdelki, barvila, organske snovi in lepila, imajo lahko močne nenewtonske lastnosti. Cilj predstavljene študije je preučitev vpliva parametrov, kot so Rayleighovo število, eksponent potenčnega zakona ter položaj in dolžina vira toplote, na stopnjo prenosa toplote ter na temperaturna in tokovna polja. Podan je shematski diagram problema dvodimenzijske votline, napolnjene z nenewtonsko potenčno tekočino, kakor tudi robni pogoji. Privzeto je konstantno Prandtlovo število 100. Tok skozi votlino je stacionaren, dvodimenzionalen in laminaren, učinki sevanja pa so zanemarljivi. Zato so uvedene brezdimenzijske vodilne enačbe. Na podlagi vodilnih enačb je izračunano lokalno in povprečno Nusseltovo število na površini toplotnega vira. Vodilne enačbe z ustreznimi robnimi pogoji so razrešene v formulaciji kontrolnih volumnov s Patankarjevim algoritmom SIMPLE. Konvekcijsko-difuzijski členi so diskretizirani po potenčnem zakonu in sistem je numerično modeliran v FORTRAN-u. Domena za reševanje je dvodimenzijska enakomerna zamaknjena mreža. Razdelan je vpliv Rayleighovega števila in eksponenta potenčnega zakona. S povečevanjem Rayleighovega števila se izboljšuje naravna konvekcija v votlini, s tem pa se povečuje prenos toplote in zmanjšuje temperatura vira toplote. Pri n < 1 se naravna konvekcija še dodatno izboljša z večanjem Rayleighovega števila. Z zmanjševanjem eksponenta potenčnega zakona (n) se zmanjšuje navidezna viskoznost fluida in se izboljšuje naravna konvekcija v votlini, to pa privede do povečanega prenosa toplote in zmanjšanja temperature vira toplote. Izboljšanje naravne konvekcije zaradi zmanjšanja vrednosti nje izrazitejše pri večjih vrednostih Rayleighovega števila. S povečevanjem dolžine vira toplote se izboljšuje naravna konvekcija v votlini, obenem pa se tudi povečuje temperatura vira toplote in se zmanjšuje povprečno Nusseltovo število. Ko se toplotni vir od leve stene votline približuje sredini spodnje stene, se toplotne razmere v votlini spreminjajo odvisno od tega, ali je eksponent potenčnega zakona večji, manjši ali enak 1. Ključne besede: naravna konvekcija, votlina, vir toplote, nenewtonski, potenčni zakon Zasnova in vrednotenje hierarhičnega algoritma za krmiljenje sistema električnega aktivnega stabilizatorja Zhenxing Kong - Dawei Pi* - Xianhui Wang - Hongliang Wang - Shan Chen Znanstveno-tehniška univerza v Nanjingu, Oddelek za strojništvo, Kitajska Članek predstavlja nov hierarhični algoritem za krmiljenje sistema sprednjega in zadnjega električnega aktivnega stabilizatorja (ASB) na štirikolesnem cestnem vozilu. Njegova naloga je aktivno krmiljenje dinamike stranskega nagibanja in vrtenja okrog navpične osi. Predlagani hierarhični krmilni algoritem vključuje regulatorje na treh ravneh. Na podlagi stabilnosti stranskega nagibanja in udobja pri vožnji je bil izbran referenčni model kota stranskega kota nagibanja, iz katerega izhaja ciljni kot nagiba. Regulator na zgornji ravni izračunava ciljni aktivni protinagibni moment za nadzor nad kotom stranskega nagibanja vozila po metodi drsečega krmiljenja, namenjen pa je izboljšanju stabilnosti stranskega nagibanja vozila. Dinamiki stranskega nagibanja in vrtenja okrog navpične osi sta sklopljeni zaradi nelinearne karakteristike pnevmatik, poleg tega pa se lahko zaradi lateralnega prenašanja obremenitve na sprednji in zadnji osi spreminjajo lastnosti krmiljenja vozila. Regulator na srednji ravni zato skrbi za dinamično porazdeljevanje ciljnega aktivnega protinagibnega momenta med sprednjo in zadnjo osjo in s tem nadzor nad vrtenjem vozila okrog navpične osi. Uporabljen je mehki regulator, ki lahko izboljša stabilnost vrtenja vozila okrog navpične osi s sprotnim spreminjanjem lateralnega prenašanja obremenitve na sprednji in zadnji osi. Na podlagi ciljnega aktivnega protinagibnega momenta na sprednji in zadnji osi se izračunava ciljni tok sprednjih in zadnjih izvršnih členov sistema ASB. Naloga PI-regulatorja na spodnji ravni je krmiljenje toka za izvršne člene ASB, izboljšanje izhodne karakteristike izvršnih členov ASB (izhodni moment in stopnja momenta) ter zagotavljanje točnega aktivnega protinagibnega momenta vozila. Za izvedbo numeričnih simulacij in eksperimentov s strojno opremo v zanki (HIL) so bili v programu MATLAB/Simulink postavljeni model dinamike vozila s 14 prostostnimi stopnjami (DOF), vhodni model ceste, model sprednjih in zadnjih izvršnih členov sistema ASB in model krmilnega sistema ASB. Postavljen je bil eksperimentalni sistem HIL za električni sistem ASB na podlagi realnočasovnega simulacijskega sistema (AutoBox), električne krmilne enote sistema ASB (ASB ECU), električnih izvršnih členov ASB in zaznaval za silo. Opravljene so bile numerične simulacije in eksperimenti HIL za značilne manevre (ničelni vhodi, J-obrat in sinusni val) z namenom vrednotenja učinkovitosti predlaganega krmilnega algoritma. Numerične simulacije in eksperimenti HIL demonstrirajo zanesljivost in natančnost hierarhičnega krmiljenja električnega sistema ABS, ki izpolnjuje zahtevana merila prečnega nagibanja in vrtenja okrog navpične osi. Predlagani krmilni algoritem lahko zmanjša prečni nagibni kot, hitrost ter spekter gostote moči (PSD) vrtenja vozila okrog navpične osi, izboljšuje pa tudi frekvenčni odgovor izvršnih členov sistema ASB za dobro dinamiko prečnega nagibanja vozila. Dinamična porazdelitev protinagibnega navora sistema ABS poleg tega zmanjšuje nihanja hitrosti vrtenja okrog navpične osi ter ohranja stabilnost vozila za dober dinamični odziv vozila. Porazdelitev protinagibnega navora pri manevru s sinusnim valom pa zaradi sklopitve dinamike prečnega nagibanja vozila in vrtenja okrog navpične osi ter zaradi fizičnih omejitev izvršnih členov ASB ne zagotavlja želenih rezultatov krmiljenja. Električni sistem ASB na podlagi predlaganega krmilnega algoritma izboljšuje stabilnost prečnega nagibanja vozila, udobje med vožnjo in splošno stabilnost vodenja. Ključne besede: aktivni stabilizator, hierarhično krmiljenje, drseče krmiljenje, mehko vodenje, numerična simulacija, eksperiment s strojno opremo v zanki SI 98 *Naslov avtorja za dopisovanje: Znanstveno-tehniška univerza v Nanjingu, Oddelek za strojništvo, Nanjing, Jiangsu 210094, Kitajska, pidawei@mail.njust.edu.cn Večciljna optimizacija sestave storitev proizvodnje v oblaku po metodi genetskega algoritma, izboljšanega z entropijo oblaka Yongxiang Li - Xifan Yao* - Jifeng Zhou Fakulteta za strojništvo in avtomobilsko tehniko, Tehniška univerza Južne Kitajske, Kitajska Za optimizacijo razporejanja proizvodnih virov v oblaku in izboljšanje učinkovitosti sestave storitev proizvodnje v oblaku je bil raziskan problem optimizacije sestave storitev proizvodnje v oblaku. V ta namen je predlagan nov pristop: genetski algoritem, izboljšan z entropijo oblaka (CEGA), ki pri optimizaciji upošteva več ciljev: stopnjo primernosti storitev, stopnjo harmonije sestave in kompleksnost sestave storitev. Podane so definicije stopnje primernosti storitev, stopnje harmonije sestave in entropije oblaka s pripadajočimi računskimi metodami. Postavljen je matematični model večciljne optimizacije sestave proizvodnih storitev v oblaku. Oblikovana je matrika stopnje primernosti storitev na podlagi štirih meril: kapacitete, želje, opreme in lokacije. Za popis povezav med storitvami v sestavi so uporabljena orodja mehke matematike, za izračun stopnje harmonije pri zapleteni sestavi storitev je vpeljan količnik interakcije. Uveden je normalen model oblaka za izboljšanje izračunov verjetnosti križanja in verjetnosti mutacij algoritma. Zasnovani so strategija šifriranja, funkcija uspešnosti ter operaciji križanja in mutacije algoritma CEGA. Za verifikacijo predlaganega algoritma CEGA na postavljenem modelu sestave je bila uporabljena naloga izdelave avtomatskega vodenega vozila (AGV). Le-ta je bila razdeljena na šest podnalog: nalogo karoserije w1, nalogo voznega sistema w2, nalogo pomožne opreme w3, nalogo napajalnega sistema w4, nalogo pomožnega krmilnega sistema w5 in nalogo glavnega krmilnega sistema w6. Naloge so bile razporejene v 14 javnih in dve zasebni oblačni storitvi. Vsi parametri so bili nastavljeni in izračunani po predlaganem algoritmu CEGA. Po 48 iteracijah CEGA je bila dosežena optimalna vrednost uspešnosti populacije 99,934, stopnja primernosti storitev 4, stopnja harmonije sestave 11,799, entropija oblaka 8,708, dobavni rok 375 ur, strošek 80.490 evrov, kromosom optimalne rešitve 1001001000100110 ter razdalja med optimalno vrednostjo ciljne funkcije in idealno točko 1,101. Glede na kodirni kromosom optimalne rešitve je naloga w1 razporejena v prvo oblačno storitev množice storitev javnega oblaka Q1, naloga w2 v prvo oblačno storitev množice storitev javnega oblaka Q2, naloga w3 v prvo oblačno storitev množice storitev javnega oblaka Q3, naloga w4 v prvo oblačno storitev množice storitev javnega oblaka Q4, naloga w5 v drugo oblačno storitev množice storitev javnega oblaka Q5, in naloga w6 v prvo oblačno storitev množice storitev zasebnega oblaka Q6. Na istem računalniku s procesorjem Intel core i3-3110 z 2,4 GHz in 4 GB pomnilnika traja reševanje problema po algoritmu CEGA 11,63 s, reševanje standardnega genetskega algoritma (SGA) pa 19,34 s. Rezultati kažejo, da konvergira algoritem CEGA hitreje kot SGA. Glavni prispevki in novosti: (1) Definirani so stopnja primernosti storitev, stopnja harmonije sestave in entropija oblaka, podane so tudi pripadajoče računske metode. S tem je razširjen nabor vplivnih dejavnikov za optimizacijo sestave proizvodnih storitev v oblaku. (2) Postavljen je matematični model za optimizacijo sestave proizvodnih storitev v oblaku, ki upošteva stopnjo primernosti storitev, stopnjo harmonije sestave in entropijo oblaka kot večciljne funkcije z omejitvami časa izvajanja, stroškov izvajanja itn. (3) Z uvedbo teorije modela normalnega oblaka za izboljšanje tradicionalnega genetskega algoritma je predlagan algoritem CEGA za reševanje matematičnega modela sestave proizvodnih storitev v oblaku z več cilji, ki je bil dodatno preverjen glede primernosti in učinkovitosti na primeru izdelave vozila AGV. Optimizacija sestave proizvodnih storitev v oblaku je kompleksen problem s številnimi vplivnimi dejavniki kot so stopnja primernosti storitev, stopnja harmonije sestave in entropija oblaka, ki so preučeni v tem članku. Nadaljnje raziskave se lahko posvetijo ostalim vplivnim dejavnikom in algoritmom. Ključne besede: proizvodnja v oblaku, optimizacija sestave storitev, entropija oblaka, stopnja primernosti storitev, stopnja harmonije sestave, genetski algoritem Vplivni parametri dinamičnega vedenja nosilne konstrukcije portalnega žerjava tipa H Rade Vasiljević1-* - Milomir Gašić2 - Mile Savković2 1 Visoka strokovna tehnična šola, Beograd, Srbija 2 Univerza v Kragujevcu, Fakulteta za strojništvo in gradbeništvo, Srbija Članek obravnava preiskavo parametrov, ki vplivajo na dinamično vedenje nosilne konstrukcije portalnih žerjavov tipa H z veliko nosilnostjo v pogojih premikanja žerjava in nihanja tovora. Uvodoma je predstavljen pomen dinamike konstrukcij. V nadaljevanju so oblikovani in razrešeni prostorski dinamični modeli ročice in portala žerjava. Jedro članka predstavlja modalno analizo ter določitev dinamičnega odziva prostorske konstrukcije po metodi končnih elementov in neposredne integracije. Podan je nov kombinirani pristop k analizi dinamičnega vedenja portalnih vrtljivih žerjavov. Sistem je bil razdeljen na dva podsistema in na podlagi tega sta bila oblikovana novi dinamični model ročice in novi model nosilne konstrukcije (portala) po metodi končnih elementov. Analiziran je vpliv spreminjanja parametrov na dinamični odziv konstrukcije. Izkaže se, da se kot nihanja tovora povečuje s povečevanjem hitrosti gibanja portala. Pri imenski hitrosti žerjava 0,6 m/s je največji kot nihanj tovora 0,0815 rad/s. Žerjav ima dobre dinamične lastnosti, če je njegova prva lastna frekvenca visoka in konstrukcija obravnavanega portalnega žerjava ima v tem smislu ugodno dinamično stabilnost (f1 = 2,65 [Hz]). Matematični model konstrukcije je bil preverjen na modelu po MKE v programski opremi SAP2000®. Analiziran je bil tudi vpliv spreminjanja vrtilnega kota ročice na dinamični odziv konstrukcije. Analiza odmikov elastičnih podpor je pokazala, da ima elastična podpora največji pozitivni odmik pri vrednosti zasuka ročice 45° in v praksi se pri tem najpogosteje privzdigne noga portalnega žerjava. Največji dinamični odmik kritične elastične podpore je znašal 78,93 mm. Po parametrih matematičnega modela konstrukcije portalnega žerjava je bila pridobljena pričakovana vrednost največjega dinamičnega odmika. Konstrukcija obravnavanega portalnega žerjava je torej najbolj občutljiva na položaj ročice 45° in ta položaj je zato najustreznejši za popis dinamičnega stanja konstrukcije. Rezultati raziskave zagotavljajo dobro osnovo in zaključke za konstruiranje portalnih žerjavov in nadaljnje raziskave na tem področju. Pomembne so široke možnosti za praktično uporabo rezultatov iz članka. Vrednosti vplivnih parametrov, ki so opredeljene v tem članku, lahko pomagajo pri omejevanju vse pogostejšega pojava poškodb na portalnih žerjavih. Ker konstrukcijske rešitve žerjavov pogosto niso optimalne, lahko vrednosti teh parametrov prispevajo k optimizaciji nosilne konstrukcije portalnih žerjavov. Rezultate dinamične analize iz tega članka bodo lahko uporabili tako obstoječi kot tudi potencialni proizvajalci tovrstnih žerjavov in si tako izboljšali svoj konkurenčni položaj na trgu. Konstruktorji žerjavov bodo lahko uporabili vrednosti vplivnih parametrov za oblikovanje določenih domnev v začetni fazi snovanja portalnih žerjavov za zanesljivost konstrukcije. Ključne besede: portalni žerjav tipa H, vrtljiva ročica, vzbujanje, analiza dinamičnega vedenja, metoda končnih elementov, dinamični odziv SI 100 *Naslov avtorja za dopisovanje: Visoka strokovna tehnična šola, Zorana Đinđića 152A, 11000 Beograd, Srbia, rade.vasiljevic@visokatehnicka.edu.rs Numerična analiza lastnosti toka fluida in prenosa toplote v pravokotnih mikrokanalih za nanofluid ZnO-etilenglikol Cuneyt Uysal* - Kamil Arslan - Huseyin Kurt Univerza Karabuk, Oddelek za strojništvo, Turčija Članek predstavlja numerično preiskavo lastnosti toka fluida in prenosa toplote v pravokotnih mikrokanalih z različnimi razmerji stranic (a = 1 to 2) za nanofluid ZnO-etilenglikol (EG). Uporabljeni so bili različni volumski deleži nanodelcev (ф = 1 % to 4 %) v nanofluidu ZnO-EG. Tok je obravnavan kot enofazen, trirazsežnosten, stacionaren, nestisljiv, toplotno razvijajoči se laminami tok. Uporabnost nanofluida ZnO-EG kot medija za prenos toplote pri hlajenju elektronskih naprav in vpliv geometrijske konfiguracije kanala sta bila preiskana z numerično analizo. Sistem je oblikovan kot mikrokanal, pritrjen na vročo površino za odvajanje toplote, ki nastaja v poljubni elektronski napravi (npr. v računalniških čipih). Toploti, ki jo ustvarja elektronska naprava, je tako izpostavljena spodnja stena mikrokanala. Za numerične izračune je bila uporabljena metoda končnih volumnov. Vodilne enačbe so razrešene skupaj z robnimi pogoji. Konvektivni členi v vodilnih enačbah so diskretizirani z gorvodno shemo drugega reda, za diskretizacijo tlaka pa je bila uporabljena standardna shema. Za razrešitev sklopitve hitrosti in tlaka je bil uporabljen algoritem SIMPLE. Momentne in energijske enačbe so bile diskretizirane po Green-Gaussovi celični metodi. Pri izračunih niso bile opažene nikakršne težave s konvergenco. Točnost kode je bila validirana s primerjavo rezultatov predstavljene študije in podatkov iz literature, ki bazirajo na eksperimentalnih podatkih za pravokotne mikrokanale v toplotno razvijajočem se laminarnem toku v razmerah konstantnega toplotnega toka. Rezultati se dobro ujemajo s podatki iz literature. Izkazalo se je, da je pri odvajanju toplote s konvektivnim prenosom najučinkovitejši mikrokanal z razmerjem stranic a = 1. Ugotovljeno je bilo tudi to, da se z dodatkom nanodelcev ZnO v čisti EG zmanjša temperaturna razlika med steno in maso nanofluida. S povečevanjem volumskega deleža nanodelcev v nanofluidu ZnO-EG se zato povečuje koeficient konvektivnega prenosa toplote, Nusseltovo število pa nasprotno s povečevanjem volumskega deleža nanodelcev v nanofluidu ZnO-EG pada. Vzrok za zmanjšanje Nusseltovega števila leži v večjem porastu koeficienta prenosa toplote s prevajanjem v primerjavi s porastom koeficienta konvektivnega prenosa toplote, ko se v čisti EG dodajajo nanodelci ZnO. Pri dani potrebni črpalni moči je bila največja vrednost toplotnega upora dosežena pri a = 2, najmanjša pa pri a = 1. Toplotni upor se zmanjšuje z rastjo volumskega deleža nanodelcev v nanofluidu ZnO-EG. Najmanjše vrednosti toplotnega upora pri dani potrebni črpalni moči so bile izračunane pri 4 % deležu ZnO v nanofluidu. Iz rezultatov sledi sklep, da je mogoče z geometrijsko konfiguracijo mikrokanala zelo učinkovito vplivati na lastnosti toka fluida in prenosa toplote. Hidrodinamično in toplotno zmogljivost mikrokanalov je tako mogoče izboljšati le z geometrijsko konfiguracijo sistema. Mikrokanali z vrednostjo a = 1 dajejo boljše rezultate tako glede toka fluida kot glede prenosa toplote. Z uporabo nanofluida ZnO-EG se je izboljšal konvektivni prenos toplote, obenem pa se je tudi povečal tlačni padec. V tem primeru nam vrednost toplotnega upora pri dani črpalni moči omogoča vpogled v prednosti in slabosti uporabe nanofluida. Nanofluid ZnO-EG je v prednosti pred čistim EG zaradi nižje vrednosti toplotnega upora pri dani črpalni moči. Ključne besede: pravokoten mikrokanal, nanofluid ZnO-EG, razmerje stranic, konvektivni prenos toplote, tlačni padec, toplotni upor, črpalna moč Numerična simulacija regulacije prehodnega tlaka v črpalnem vodovodnem sistemu z izboljšano obvodno cevjo Xiaozhou Li- Manlin Zhu* - Jiancang Xie Tehniška univerza Xi'an, Državni laboratorij za ekološko hidrotehniko na sušnih območjih, Kitajska Študija ima tri cilje: 1) Predlagan je ukrep za preprečevanje vodnih udarov z izboljšano obvodno cevjo (IBP), ki odpravlja omejitve konvencionalnih obvodnih cevi (CBP). Pojasnjen je tudi način delovanja IBP; 2) Na podlagi predhodnih študij in metode karakteristik (MOC) je postavljen matematični model robnih pogojev IBP in opisana je rešitev tega matematičnega modela; in 3) Kot primer uporabe zaščitnega ukrepa IBP in za numerično simulacijo učinka preprečevanja vodnega udara na črpalki z zaščitno napravo IBP je uporabljen vodovodni sistem v industrijski coni. Uporaba obvodne cevi v sistemu črpalke je ekonomična in preprosta metoda za preprečevanje vodnih udarov pri regulaciji ekstremnih prehodnih tlakov, ki lahko nastopijo v primeru odpovedi črpalke. Pri tehnični uporabi konvencionalnih obvodnih cevi (CBP) trenutno obstajajo določene omejitve in za njihovo odpravo je bila cev CBP v tej študiji modificirana v izboljšano obvodno cev (IBP). Za razliko od CBP je krmilni ventil v sistemu IBP hidravlični krmilni ventil, ki uporablja krmilnik in brezprekinitveni vir napajanja (UPS) za odpiranje in zapiranje ventila. Modifikacija omogoča natančno upravljanje hidravličnega krmilnega ventila in preprečuje vodne udare. Za preučitev učinkovitosti regulacije prehodnega tlaka v sistemu IBP z uporabo metode karakteristik (MOC) je bil postavljen matematični model kompleksnega sistema, ki vsebuje IBP, nepovratni ventil in vodno črpalko. Matematični model je bil postavljen na naslednjih domnevah: (1) Od trenutka odpovedi črpalke do trenutka, ko pretok na izhodu vodne črpalke pade na nič, je lokalni koeficient upora nepovratnega ventila konstanten. Ko je pretok na izhodu vodne črpalke Q6 < 0, se nepovratni ventil v trenutku zapre. (2) Lokalne tlačne izgube v vozliščih tipa T nad in pod vodno črpalko so zanemarjene. (3) Vozlišča tipa T, ki povezujejo cev z obvodno cevjo nad in pod vodno črpalko, so zelo kratka. Pretočne izgube so zanemarjene in vozlišča so analizirana kot integralna enota. Za primer numerične simulacije prehodnih pojavov po odpovedi črpalke s programsko opremo MATLAB je bil uporabljen vodovodni sistem v industrijski coni. Podana je analiza in primerjava ekstremnih tlačnih sprememb in ovojnice tlaka v cevi v značilnem prerezu vodovodnega sistema s črpalko. Primerjava rezultatov simulacij je pokazala, da lahko ukrep za zaščito pred vodnim udarom z IBP (s spremembo nastavitve hidravličnega krmilnega ventila v trenutkih t1 = 2,5 s, t2=3 s, t3 = 5 s in t4 = 15 s) učinkovito zmanjša največji tlak pri vodnem udaru v vodovodnem sistemu s črpalko. Nekateri problemi, ki so bili razkriti v tej raziskavi, pa zahtevajo še dodatno pozornost: 1) Matematični model IBP ne upošteva dolžine in premera obvodne cevi, vanj pa je vgrajen le lokalni količnik upora ventila obvodne cevi. V naslednjem koraku bi bilo zato treba vključiti tudi dolžino l in premer d obvodne cevi za preučitev vpliva dolžine in premera obvodne cevi na sposobnost preprečevanja vodnih udarov. 2) Numerična simulacija in preračun prehodnega pojava v vodovodu s črpalko sta pokazala, da je rešitev IBP sposobna zaščititi vodovod pred škodo zaradi vodnih udarov. Za določitev okvirnih vrednosti za tehnične aplikacije pa bi bilo treba določiti še občutljivost ključnih parametrov, ki vplivajo na učinek preprečevanja vodnega udara z IBP (dolžina l in premer d obvodne cevi, zakasnitev odpiranja hidravličnega krmilnega ventila obvodne cevi tj, trenutek odpiranja t2, trajanje odpiranja t3 in trenutek zapiranja t4). Ključne besede: izboljšana obvodna cev, hidravlični prehodni pojav, vodovod s črpalko, krmilnik z brezprekinitvenim napajanjem, numerična simulacija *Naslov avtorja za dopisovanje: Tehniška univerza Xi'an, Državni laboratorij za ekološko hidrotehniko na sušnih območjih, NO.5 South Jinhua Road, Xi'an 710048, Shaanxi ,Kitajska, zhuml@xaut.edu.cn DOKTORSKA DISERTACIJA, ZNANSTVENA MAGISTRSKA DELA DOKTORSKA DISERTACIJA Na Fakulteti za strojništvo Univerze v Ljubljani so obranili svojo doktorsko disertacijo: • dne 16. septembra 2016 Miha PIPAN z naslovom: »Karakterizacija in krmilni algoritem pnevmatičnih mišic« (mentor: prof. dr. Niko Herakovič); Krmiljenje pnevmatičnih mišic pri dinamičnih spremembah raztezka zahteva podrobno poznavanje njihovih statičnih in dinamičnih karakteristik ter krmilnih in ostalih elementov, ki sestavljajo celoten sistem. Predvsem je treba podrobno poznati in upoštevati vse parametre, ki definirajo matematični model mišice in vplivajo na njeno dinamično odzivnost ter natančnost pozicioniranja. V nasprotju s konvencionalnimi pnevmatičnimi aktuatorji (batnični valji, brezbatnični valji, rotacijski moduli itd.) imajo pnevmatične mišice nelinearno karakteristiko generirane sile in pomika v odvisnosti od tlaka ter masnega toka. V doktorski disertaciji so zato podrobno obravnavani in analizirani najpomembnejši parametri, ki vplivajo na statične in dinamične karakteristike pnevmatičnih mišic. Nadgrajen in znatno izboljšan je najbolj uveljavljen Chouov statični matematični model pnevmatične mišice ter dinamični model z variabilno nelinearno togostjo ter dušenjem, kar v obstoječih modelih še ni bilo obravnavano. Prav tako je v doktorskem delu izvedena analiza pretočnih karakteristik hitrega preklopnega ventila, krmiljenega s PWM signalom, s pomočjo katere smo izdelali krmilni algoritem hitrih preklopnih ventilov. Velik poudarek je na razvoju adaptivnega krmilnega algoritma (A-PID) za krmiljenje hitrega preklopnega pnevmatičnega ventila za hitro dinamično krmiljenje krčenja pnevmatične mišice z upoštevanjem spreminjanja njenih karakteristik in s tem parametrov, ki vplivajo na odzivnost hitrega preklopnega ventila. Izdelan je linearni aktuator s pnevmatično mišico in hitrima preklopnima ventiloma za testiranje A-PID algoritma. Po podrobni teoretični in eksperimentalni analizi linearnega aktuatorja je narejena primerjalna analiza rezultatov simulacije in eksperimentalno določene odzivnosti aktuatorja. Narejena je tudi primerjava rezultatov razvitega A-PID algoritma za krmiljenje pnevmatične mišice s hitrimi preklopnimi ventili z obstoječimi uveljavljenimi krmilnimi algoritmi in simulacijskimi modeli; • dne 20. septembra 2016 Somayeh AKBARI z naslovom: »Tribochemical mechanisms in boundary and mixed lubrication regime (SLO: Tribokemijski mehanizmi v mejnem in mešanem režimu mazanja)« (mentor: prof. dr. Mitjan Kalin); V doktorski nalogi smo se osredotočili na temeljno razumevanje interakcij med prevleko in aditivom ZDDP. S tega vidika smo najprej raziskali tvorbo termičnega filma ZDDP, pri dobro kontroliranem statičnem termalnem testu tj. brez tribološkega kontakta, na H-DLC in Si-DLC prevlekah. Za referenco smo enako storili tudi na jeklenih površinah. Raziskali smo vpliv različnih temperatur in koncentracij na termično tvorbo ZDDP na vseh površinah. Na podlagi dobljenih rezultatov smo predlagali mehanizme reakcij termičnih filmov aditiva ZDDP na DLC prevleke. Določili smo kemijske interakcije ZDDP s površinami pri statičnih pogojih, kar predstavlja običajen pristop za primerjavo in dopolnitev triboloških rezultatov, z namenom obrazložitve tribokemijskih mehanizmov. Po raziskavi tvorbe termičnih filmov ZDDP na vse obravnavane površine, smo tako kot pri statičnih testih, ovrednotili tvorbo triboloških filmov ZDDP na taistih površinah, ki so nastali pri različnih temperaturah in koncentracijah. Določili smo vpliv teh parametrov na tribološko obnašanje vseh kontaktov z ZDDP v mejnem režimu mazanja. Reaktivnost vseh površin z ZDDP, v primeru termičnih in triboloških testov, je bila določena predvsem s kemično karakterizacijo, za kar je bila uporabljena tehnika infrardeče spektroskopije s Fourierjevo transformacijo z zmanjšano totalno refleksijo (ATR-FTIR), ki v literaturi ni pogosto uporabljena, in rentgensko fotoelektronsko spektroskopijo; • dne 21. septembra 2016 Hubert KOSLER z naslovom: »Sprotna triangulacija in adaptivna daljinska obdelava z vlakenskim laserjem« (mentor: prof. dr. Janez Možina, somentor: izr. prof. dr. Matija Jezeršek); Doktorska disertacija obravnava razvoj robotskega sistema za lasersko daljinsko obdelavo, ki deluje na osnovi sprotnega merjenja oblike površine obdelovanca. Najprej je bil razvit merilnik oblike površine obdelovanca na osnovi laserske triangulacije vključno z algoritmi iskanja in sledenja zvarnega spoja. Vzporedno s tem so bile obravnavane tudi nekatere možnosti uporabe razvitega merilnika za nadzor širine rege in oblike zvarov ter za adaptivno raziglevanje ulitkov. V drugem delu je bil najprej postavljen robotski sistem za varjenje rotorjev z vlakenskim laserjem moči 500 W v kombinaciji s horizontalno rotacijsko zunanjo robotsko osjo. Vodenje laserskega snopa po obdelovancu je bilo v tem primeru izpeljano s konvencionalnim učenjem robota. Nadalje je bila razvita laserska celica za lasersko daljinsko varjenje z močjo 6 kW iz razdalje 600 mm od obdelovanca. Načrtovanje obdelovalne trajektorije poteka v tem primeru vizualno z dvema pilotnima snopoma. Kot sinteza sistemov laserske triangulacije, robotskega adaptivnega vodenja in daljinske obdelave je bil postavljen sistem za adaptivno robotsko daljinsko obdelavo s sprotnim merjenjem oblike površine obdelovanca. Za kalibracijo sistema je bila razvita izvirna metoda na osnovi primerjave z meritvami oblike referenčne površine. Delovna trajektorija laserskega snopa se v tem primeru določi z multikriterijsko detekcijo roba oziroma spoja na osnovi sprotne meritve oblike površine. Natančnost nastavitve položaja gorišča laserskega snopa znaša 70 ^m v ravnini obdelave in 120 ^m v smeri pravokotno na to ravnino. Izvirni algoritmi za avtomatsko določanje obdelovalne trajektorije omogočajo do 30-kratno skrajšanje časa učenja sistema v primerjavi z vizualno metodo. Sistemi laserske triangulacije in daljinske obdelave, ki so bili razviti v okviru disertacije, predstavljajo pomembno obogatitev prodajnega programa slovenskega sistemskega integratola robotskih sistemov; • dne 23. septembra 2016 Vili PEPEL z naslovom: »Nastanek in rast poškodbe duktilnega materiala pri malocikličnem utrujanju« (mentor: prof. dr. Ivan Prebil, somentor: izr. prof. dr. Robert Kunc); Določitev dobe trajanja izdelka je pomembna naloga razvojnih inženirjev. Napovedovanje natančne dobe trajanja privede do nižjih stroškov proizvodnje, lažjih konstrukcij, preprečevanja dolgotrajnih in dragih zastojev in lahko prepreči tudi potencialno življenjsko nevarne situacije v fazi obratovanja. Poleg znanih metod merjenja in določanja nastanka in razvoja makro razpok pri malocikličnem obremenjevanju je uporabljena tudi metoda akustične emisije. Za analizo vhodnih podatkov akustične emisije smo uporabili Metodo empirične dekompozicije (EMD), Metodo valčne transformacije (WT) in Kalmanov filter (KF). Metodo akustične emisije pri malocikličnem obremenjevanju smo za določanje nastanka poškodbe preverili s pomočjo preizkusov. Eksperimentalna preverjanja smo izvedli na materialu nižjo natezno trdnostjo EN 42CrMo4 in materialu z višjo natezno trdnostjo Armox 500T. Predstavljeni so rezultati izvedenih meritev ter prednosti in omejitve uporabljenih neporušitvenih metod za določanje dobe trajanja strojnih delov. S pomočjo akustične emisije se doseže zgodnejša določitev nastanka poškodbe kot pri metodah modula elastičnosti, napetosti in plastične deformacijske energije; • dne 26. septembra 2016 Alexander KUZNETSOV z naslovom: »Optimisation and control of an annular laser beam droplet generation process from a metal wire (Optimizacija in krmiljenje procesa tvorjenja kapljic iz kovinske žice s kolobarjastim laserskim žarkom)« (mentor: prof. dr. Edvard Govekar); V delu so predstavljeni rezultati raziskave laserskega tvorjenja kapljic (LTK) iz kovinske žice z uporabo kolobarjastega laserskega žarka soosnega z žico. Obravnavano je lasersko tvorjenje kapljic na zahtevo (LTKnZ) in proces laserskega tvorjenja zaporedja kapljic z vnaprej določeno frekvenco. Predlagana je bila nova strategija procesa LTKnZ, kjer so vrednosti glavnih parametrov procesa določene glede na zahtevane lastnosti kapljice (premer, temperatura). Za namen optimizacije in krmiljenja procesa je bila izvedena analiza in karakterizacija negotovosti sistema in procesa LTK. Da bi zagotovili ponovljivo LTKnZ, je bila izvedena optimizacija procesnih parametrov z ozirom na dinamiko procesa tvorjenja kapljic in pripadajočo temperaturo kapljice. Pri tem so bili v odvisnosti od lege vratu kapljice glede na gorišče kolobarjastega žarka opaženi dinamsko različni režimi ločevanja viseče kapljice vključujoč; ločitev viseče kapljice s prisotnostjo izbrizgov taline, uspešna ločitev, ponovna združitev ločene kapljice s koncem žice in neuspešna ločitev viseče kapljice. Na primeru Ni žice je bilo pokazano, da iz žice premera 0,6 mm lahko tvorimo kapljice s premerom v območju od 0,85 mm do 1,25 mm z ločljivostjo 50 |im in negotovostjo 15 |im. Temperatura tvorjenih kapljic je odvisna od premera kapljice in se nahaja v intervalu od 1650 °C do 1850 °C. Za tvorjenje zaporedja kapljic je bil zasnovan krmilni sistem, ki ob upoštevanju negotovosti sistema in procesa LTK, za vsako tvorjeno kapljico v zaporedju izbere optimizirane vrednosti procesnih parametrov glede na začetne pogoje konca žice. S predlaganim krmilnim sistemom smo v našem primeru dosegli stabilno tvorjenje zaporedja kapljic z vnaprej določenimi lastnostmi in z izbrano frekvenco tvorjenja kapljic do 3 Hz; • dne 26. septembra 2016 Janez UREVC z naslovom: »Karakterizacija mehanskega odziva humane karotidne arterije v fiziološkem stanju« (mentor: prof. dr. Boris Štok, somentor: prof. dr. Milan Brumen); V prvem delu disertacije je numerično obravnavan mehanski odziv človeške karotidne arterije in vivo z upoštevanjem zaostalih napetosti. Zaostalo napetostno stanje arterije je določeno na osnovi termomehanskega (TMA) postopka, ki smo ga razvili in ga predlagamo v tem delu. Metodologija omogoča določiti mehanski odziv arterije na osnovi njenih geometrijskih podatkov in vivo; torej brez poznavanja začetne breznapetostne oblike arterije. Podatki arterije, na osnovi katerih je metodologija prikazana, so pridobljeni in vivo z uporabo neinvazivnih metod na 28-letnem zdravem moškem osebku. Mehanski odziv žile je obravnavan kot nelinearno elastičen in izotropen, pristop TMA pa je prikazan na obravnavi eno- in dvo-slojne strukture žile. Rezultati izračunov napetostnega stanja kvantitativno sovpadajo s sorodnimi študijami iz literature. Prikazan in analiziran je tudi vpliv zaostalih napetosti na mehanski odziv karotidne arterije pri vstavitvi žilne opornice. Drugi del disertacije se naša na določitev hemodinamskih razmer v predelu razvejišča karotidne arterije. Prikazan je vpliv spremenjene viskoznosti krvi, kot posledica spremenjene deformabilnosti rdečih krvnih celic, na odziv žilne stene oz. na parametre, ki so v povezavi z patogenezo ateroskleroze; • dne 28. septembra 2016 Marko JOVANOVIĆ z naslovom: »Optimizacija modularnih proizvodnih strojev za izdelavo gumarskih izdelkov« (mentor: prof. dr. Ivan Prebil, somentor: prof. dr. Marko Starbek); Proizvajalci pnevmatik, ki želijo v kompleksnem proizvodnem okolju sodobnega tržišča ostati konkurenčni, se morajo hitro prilagajati spremembam, ki nastanejo v proizvodnih procesih. Osnovna problema, ki jih v takšnih razmerah izkažejo tradicionalni proizvodni sistemi, sta počasen odziv v procesu optimizacije proizvodnega procesa in nizka stopnja prilagodljivosti na motnje v sistemu. Cilj naše raziskave je razvoj distribuiranega in adaptivnega krmilnega pristopa, ki temelji na konceptu holonskega krmiljenja. Predstavljeni pristop distribuiranega holonskega krmiljenja proizvodnega sistema za proizvodnjo pnevmatik omogoča: dinamični odziv v primeru pojava optimizacijskih zahtev, zmanjšanje vpliva motenj na produktivnost sistema, dobro izkoriščanje delovnih sredstev in v prihodnosti manjša investicijska vlaganja v delovna sredstva in opremo. Razvili smo okolje za izvedbo distribuiranega holonskega pristopa krmiljenja proizvodnega sistema za proizvodnjo pnevmatik z uporabo funkcijskih blokov, ki so v soglasju s standardom IEC 61499. Razvito virtualno proizvodno okolje omogoča analizo proizvodnega procesa, ponazoritev proizvodnih operacij, krmiljenje virtualnega proizvodnega sistema, upravljanje s parametri proizvodnih operacij in analizo načina delovanja krmilnega sistema še preden bo vgrajen v realni proizvodni sistem. Vrednotenje distribuiranega holonskega pristopa krmiljenja v virtualnem proizvodnem okolju je bilo izvedeno s pomočjo simulacij za različne proizvodne scenarije, ki upoštevajo obratovanje sistema v stabilnem in nestabilnem stanju. Pričakujemo, da bo realna izvedba predstavljene tehnologije v prihodnosti zvišala produktivnost, izkoriščenost delovnih sredstev in robustnost proizvodnega sistema; • dne 28. septembra 2016 Franc RAVNIK z naslovom: »Akustična emisija pri gašenju jekel« (mentor: prof. dr. Janez Grum); Toplotna obdelava pogosto predstavlja zadnjo stopnjo v proizvodnji strojnih delov. Namen nadzora procesa gašenja z izbiro najprimernejših parametrov gašenja je zagotoviti zahtevano trdoto in zaostale napetosti, zlasti na površini strojnega dela. Predvsem v masovni proizvodnji, kjer se stremi po čim boljših mehanskih lastnostih in kakovosti ter tako k zmanjšanju teže in stroškov, je ta še bolj izrazita. Končne mehanske lastnosti so tako odvisne od optimalne izbire parametrov gašenja in nadzora nad samim procesom gašenja. Program doktorskega dela je obsegal poglobljeno raziskavo zvočnih pojavov pri gašenju jekel. V okviru doktorskega dela je bila raziskana možnost povezave zvočnih pojavov s kinematiko omočenja gasilnega sredstva in vročega predmeta z drugimi pojavi, ki potekajo med procesom gašenja. Rezultate poteka procesa gašenja smo povezali z mehanskimi lastnostmi jekel po gašenju. Nakazana je možnost, da iz zajetih zvočnih signalov identificiramo primernost oziroma kakovost procesa gašenja, kar bi omogočilo boljši nadzor kakovosti procesa gašenja jekel. V ta namen je bil izdelan eksperimentalni sistem s hidrofonom, ki je omogočil izvedbo procesa in zajem nastalih zvočnih signalov ter napoved pojavov iz nadaljnje obdelave zajetih signalov. Obdelanih je bilo več jeklenih vzorcev različnih oblik, ki so bili gašeni v gasilnih sredstvih različne sposobnosti hitrosti ohlajanja. Primerjava rezultatov je pokazala, da uporabljeni sistem za zajem signalov omogoča ustrezne rezultate. Raziskana je bila tudi možnost identifikacije zvočnih pojavov, ki nastanejo kot posledica deformacije predmeta/vzorca in nastajanja razpok zaradi previsokih notranjih napetosti. Rezultati kažejo, da ta možnost nakazuje k pristopu največje aplikativne vrednosti. Pridobljeni eksperimentalni rezultati, dobljeni po različnih metodah procesiranja, so kriteriji na podlagi katerih se odločamo o uspešnosti posamezne metode za nadzor procesa gašenja in spremljanje kakovosti izdelkov; • dne 29. septembra 2016 Senad OMEROVIĆ z naslovom: »Vpliv položaja prsnega koša na gibanje in obremenitev vratne hrbtenice v fazi trka« (mentor: prof. dr. Ivan Prebil, somentor: izr. prof. dr. Robert Kunc); Nihajne poškodbe vratne hrbtenice predstavljajo izjemno ekonomsko breme za družbo s simptomi poškodbe, ki lahko v določnih primerih trajajo več kot eno leto. Številne raziskave nihajne poškodbe so bile narejene za normalen položaj sedenja, kar ima za posledico slabo raziskan vpliv različnih leg sedenja udeležencev v prometnih nezgodah na nihajno poškodbo. Zato je bil razvit podroben numerični model človeškega vratu z metodo končnih elementov. Izvedena je bila primerjalna analiza vpliva štirih tipskih načinov sedenja na možnost pojava nihajne poškodbe. Štirje raziskani položaji sedenja so: Normalen položaj, v katerem je hrbet naslonjen na naslon sedeža z glavo v pokončni legi; Trup nagnjen naprej, za približno 10° glede na naslon sedeža; Glava nagnjena naprej, tako da vrat tvori kot okoli 20° glede na lego vratu v Normalni legi in Glava in trup nagnjena naprej za 20° oz. 10°. Primerjalna študija je obsegala analizo deformacije kapsularnega, anteriornega, posteriornega in ligamenta flava ter nivo nefiziološke zakrivljenosti vratu v obliki črke S. Rezultati raziskave so pokazali, da sta položaja sedenja Glava nagnjena naprej in Glava in trup nagnjena naprej najbolj nevarna za poškodbe kapsularnih ligamentov v zgornjem oz. spodnjem delu vratne hrbtenice. Za nefiziološki pojav zakrivljenosti vratne hrbtenice v obliki črke S, sta najbolj nevarna položaja sedenja Glava nagnjena naprej in Trup nagnjen naprej. Predstavljena raziskava je pokazala vpliv položaja sedenja na možnosti poškodbe posameznih mehkih tkiv vratne hrbtenice. Z uporabo razvitega numeričnega modela človeškega vratu pa je možno preučevanje mehanizma nihajne poškodbe s ciljem njenega preprečevanja. • dne 30. septembra 2016 Ivo PRAH z naslovom: »Kalibracija simulacijskih modelov motorjev z notranjim zgorevanjem« (mentor: izr. prof. dr. Tomaž Katrašnik, somentor: prof. dr. Ferdinand Trenc); Delo obravnava postopek računalniško krmiljene kalibracije kombiniranega 0D in 1D termodinamskega simulacijskega modela MNZ s prisilnim polnjenjem. Cilj kalibracije je določitev takšnih vrednosti vhodnih podatkov simulacijskega modela, ki vodijo v čim manjše razlike med simulacijskimi rezultati in izmerjenimi veličinami ob zagotavljanju čim krajšega časa kalibracije. Postopek računalniško krmiljene kalibracije je zasnovan na inovativni interakciji optimizacijskih metod in fizikalnih osnov posameznih podsklopov MNZ, ki so uporabljene pri razdelitvi celotnega simulacijskega modela MNZ na ustrezne podsisteme in upoštevanju fizikalnih zakonitosti pri določitvi parametrov komponent. Inovativna večstopenjska interakcija optimizacijskih metod in fizikalnih osnov tako omogoča v nasprotju z ustaljeno uporabo zgolj optimizacijskih metod, uspešno kalibracijo velikega števila vhodnih podatkov, kar je ključno za uspešno kalibracijo modelov sodobnih motorjev. * Na Fakulteti za strojništvo Univerze v Mariboru sta obranila svojo doktorsko disertacijo: • dne 1. septembra 2016 Martin PODGRAJŠEK z naslovom: »Razvoj modela za napovedovanje zdržljivosti preoblikovalnih orodij za vroče kovanje s trdimi prevlekami« (mentor: prof. dr. Zoran Ren); V doktorski disertaciji je obravnavana tematika razvoja modela za določevanje zdržljivosti orodja za vroče kovanje, zaščitenega s trdo prevleko. Kovaško orodje za delo v vročem je med obratovanjem izpostavljeno visokim cikličnim termomehanskim, kemičnim in tribološkim obremenitvam, ki po določenem času na tribološkem stiku povzročajo nastanek različnih mehanizmov obrabe, poškodb, plastične deformacije in termomehanskega utrujanja materiala. Zdržljivost orodja lahko izboljšamo z oplemenitenjem efektivnih površin gravure na osnovi ustrezne difuzijske plasti in trde prevleke. Pri oplemenitenju površin je vedno bolj v uporabi tehnologija zaščite gravure s postopkom »dupleks«, ki ga sestavlja kemo-termična obdelava ter nanos trde obrabno obstojne prevleke. V okviru računalniških simulacij je bil izdelan numerični model vročega kovanja, na podlagi katerega je bila izdelana analiza glavnih parametrov vročega kovanja. Na osnovi numeričnega modela smo v gravuri kovaškega orodja določili temperaturno in napetostno - deformacijsko polje, ki predstavljata osnovo za izvajanje eksperimentalnega dela ter določitve dobe trajanja kovaškega orodja. Pri eksperimentalnem delu je bila izvedena podrobna metalografska analiza poškodb in obrabe kovaškega orodja, zaščitenega s trdo prevleko. S pomočjo elektronskega mikroskopa smo določili karakteristična področja gravure in analizirali glavne mehanizme poškodb površine orodja. Z uporabo trdomera smo analizirali potek trdote v odvisnosti od globine delovne plasti gravure. Metalografska analiza je na izpostavljenih mestih gravure pokazala izrazite poškodbe površine v obliki luščenja materiala zaradi prisotnosti razpok v difuzijski coni. Kritične poškodbe delovne plasti gravure so nastale že po 6.500 udarcih kladiva. Nanos trde prevleke na kovaškem orodju ni prispeval k izboljšanju zdržljivosti kovaškega orodja. Zdržljivost orodja je bila pri izbrani konstrukciji gravure in pri danih termomehanskih obremenitvah manjša od pričakovane. Za izboljšanje zdržljivosti orodja je potrebno optimizirati konstrukcijo orodja in karakteristike difuzijske plasti. Izvedeni sta bili tudi analizi statične trdnosti in analiza parametrov nizkocikličnega utrujanja orodnega jekla pri povišani temperaturi 150 °C in 500 °C na osnovi testiranj preizkušancev z in brez trde prevleke. Eksperimentalni preizkusi so pokazali, da je v režimu nizkocikličnega utrujanja pri napetostih okoli meje tečenja vpliv PVD površinske prevleke nezanemarljiv. V režimu visokocikličnega utrujanja pri napetostih v elastičnem področju pa ima PVD prevleka pozitiven učinek na utrujenostno obnašanje materiala, kar prispeva k daljši obratovalni dobi v primerjavi z materiali brez prevleke. Ugotovljene utrujenostne parametre je mogoče uporabiti za napovedovanje obratovalne dobe podobnih orodij za vroče kovanje v pogojih nizkocikličnega utrujanja; • dne 28. septembra 2016 Matej BRAČIČ z naslovom: »Površinska obdelava silikona s polisaharidi za razvoj protimikrobnih uroloških katetrov« (mentor: prof. dr. Lidija Fras Zemljič); V doktorskem delu je predstavljena uporaba alternativnih površinskih obdelav na osnovi polisaharidov za izboljšanje protimikrobnih lastnosti in znižanje adhezije biomolekul na silikonskih površinah, ki se uporabljajo za izdelavo uroloških katetrov. Uporaba katetrov je namreč tesno povezana z visokim tveganjem nastanka mikrobioloških infekcij, ki se velikokrat končajo z dolgoročnimi zdravstvenimi težavami. V ta namen so, kot alternativa antibiotikom in obdelavam s kovinskimi ioni, v doktorskem delu bile uporabljene površinske obdelave na osnovi polisaharidov kot so hitozan, karboksimetiliran hitozan in sinergistična formulacija med hialuronsko kislino in naravnim površinsko aktivnim sredstvon na osnovi lizina. Polisaharidi so bili nanešeni na površino silikona iz raztopin in disperzij, ki so bile prvotno okarakterizirane z nevtralizacijskimi titracijami, dinamičnim sipanjem svetlobe in vrstično elektronsko mikroskopijo z namenom pridobitve informacij o njihovem površinskem naboju in velikosti delcev v disperziji v odvisnosti od vrednosti pH medija. Sledila je podrobna študija interakcij med raztopinami in disperzijami polisaharidov z modelnimi silikonskimi površinami, kjer je bil preučen vpliv njihove vrednosti pH, koncentracije soli in različnih metod površinske aktivacije silikona s pomočjo uporabe natančne kremenove mikrotehtnice. Kot modelni sistem so bili uporabljeni tanki silikonski filmi narejeni z raztapljanjem silikona v toluenu in tvorbo filmov s »spin-coating« metodo. Znanje pridobljeno iz študija interakcij modelnih filmov je bilo uporabljeno za površinsko obdelavo realnih silikonskih sistemov kot so silikonske plošče in silikonske cevke. Študiju adsorpcije/desorpcije je sledila podrobna karakterizacija površinske morfologije, površinske kemije ter mehanske in kemijske stabilnosti površinskih obdelav na osnovi polisaharidov za modelne in realne sisteme. V ta namen so bile uporabljene različne mikroskopske in spektroskopske metode, nevtralizacijske titracije ter metode za določanje mehanskih lastnosti. Na zadnje so bile ovrednotene še protimikrobne lastnosti ter lastnosti proti adheziji biomolekul na površini silikona. Protimikrobne analize so bile izvedene na po gramu pozitivne in po gramu negativne mikroorganizme, ki se najpogosteje nahajajo v okuženem urinu, medtem ko so bile lastnosti proti adheziji biomolekul izvedene z adsorpcijo proteina albumina, fobrinogena in lizocima na zgolj modelne funkcionalne filme Rezultati so pokazali, da je homogeno nanešene, kakor tudi stabilne površinske obdelave silikonskih površin mogoče doseči z adsorpcijo polisaharidov iz njihovih disperzij, kjer delci dosegajo velikosti med 200 nm in 300 nm. Takšne disperzije je mogoče pripraviti z obarjanjem polisaharidov; tj. z natančnim uravnavanjem vrednosti pH raztopine hitozana na pH = 6.5, karboksimetil hiztozana na pH = 7 in z mešanjem hialuronske kisline in naravne površinsko aktivne snovi na osnovi lizina pri koncentracijah 2.5x10-4 mol/L in 1.25x10-3 mol/L za hialuronsko kislino in 5.0x10-4 mol/L in 1.2x10-3 mol/L za površinsko aktivno snov.. Z višanjem števila nanesenih slojev polisaharidov je prav tako mogoče povečati maso in debelino slojev, kar neposredno vpliva na povišanje protimikrobnih lastnosti, in se odraža v do 90 odstotnem zaviranju rasti patogenih mikroorganizmov kot so Escherichia coli, Staphylococcus aureus in Pseudomonas aeruginosa.. Pri tem se je kot najboljša izkazala formulacija med hialuronsko kislino in naravno površinsko aktivno snovjo na osnovi lizina, ki zaradi svojega amfoternega značaja izkazuje tudi najboljše lastnosti proti adheziji biomolekul, medtem ko se obe vrsti hitozana izkažeta slabše. Iz dobljenega znanja in rezultatov doktorske naloge je mogoče sklepati, da površinska obdelava silikonskih površin na osnovi polisaharidov uspešna, v kolikor se izvede iz disperzij polisaharidov. Takšni površinski nanosi uspešno zavirajo rast mikroorganizmov, ki povzročajo infekcije pri uporabi katetrov. ZNANSTVENA MAGISTRSKA DELA Na Fakulteti za strojništvo Univerze v Ljubljani so z uspehom zagovarjali svoje magistrsko delo: • dne 6. septembra 2016 David KRALJ z naslovom: »Razvoj sistema za sledenje in analitiko obdelovalnih sistemov« (mentor: prof. dr. Alojzij Sluga); • dne 6. septembra 2016 Katja OZIMIC z naslovom: »Stroškovni vidiki sočasnega osvajanja izdelka« (mentor: prof. dr. Marko Starbek, somentor: izr. prof. dr. Janez Kušar); • dne 7. septembra 2016 Borivoj GRD EN z naslovom: »Reinženiring logističnih procesov« (mentor: prof. dr. Marko Starbek, somentor: izr. prof. dr. Janez Kušar); • dne 7. septembra 2016 Borut POGAČNIK z naslovom: »Optimizacija logističnih procesov pri servisnih posegih na letalu« (mentor: prof. dr. Jožef Duhovnik, somentor: izr. prof. dr. Jože Tavčar); • dne 8. septembra 2016 Andrej FARAZIN z naslovom: »Termoekonomska analiza hibridnega sistema daljinskega hlajenja« (mentor: prof. dr. Alojz Poredoš); • dne 13. septembra 2016 Jožef VALENTINČIČ z naslovom: »Simulacija delovanja težkega orožja pri alternativnem izstrelku in enakem obremenilnem učinku na vitalne dele orožja« (mentor: prof. dr. Franc Kosel); • dne 14. septembra 2016 Samo RESNIK z naslovom: »Načrtovanje in vodenje individualne proizvodnje« (mentor: prof. dr. Marko Starbek); • dne 15. septembra 2016 Andrej PAPEŽ z naslovom: »Analiza konvektivnega prenosa toplote s pomočjo infrardeče termografije« (mentor: prof. dr. Iztok Golobič); • dne 16. septembra 2016 Danjela PRINČIČ z naslovom: »Vpliv polnila in robnih učinkov na izolacijsko sposobnost vakuumskega izolacijskega panela« (mentor: prof. dr. Iztok Golobič); • dne 19. septembra 2016 Iza UKMAR z naslovom: »Vpliv tehničnih zahtev na zmanjšanje rabe energije in na energijsko učinkovitost klimatskih naprav/sistemov v realnem okolju« (mentor: prof. dr. Vincenc Butala); • dne 23. septembra 2016 Denis ORAČ z naslovom: »Optimiranje preoblikovanja termoplastičnih kompozitov« (mentor: izr. prof. dr. Zlatko Kampuš, somentor: izr. prof. dr. Ivan Bajsić); • dne 29. septembra 2016 Andrej WEISS z naslovom: »Vpliv načina ohlajanja zvarnega spoja po varjenju na mehanske lastnosti jekla z oznako P91« (mentor: prof. dr. Janez Tušek); • dne 29. septembra 2016 Iztok PALČIČ z naslovom: »Zmanjšanje zaostalih napetosti v zvarnih spojih z vibriranjem« (mentor: prof. dr. Janez Tušek); • dne 29. septembra 2016 Andrej BIČEK z naslovom: »Identifikacija, vrednotenje in zmanjševanje hrupa sesalne enote« (mentor: prof. dr. Mirko Čudina); • dne 29. septembra 2016 Niko KRANJAC z naslovom: »Vodenje sistema ravnanja z okoljem v podjetju« (mentor: izr. prof. dr. Andrej Senegačnik, somentor: prof. dr. Matija Tuma). * Na Fakulteti za strojništvo Univerze v Mariboru so z uspehom zagovarjali svoje magistrsko delo: • dne 5. septembra 2016 Igor IVANOVSKI z naslovom: »Varnost in zanesljivost geotermalne elektrarne« (mentor: izr. prof. dr. Jure Marn); • dne 5. septembra 2016 Damijan KANDUTI z naslovom: »Eksperimentalna določitev toka zmesi pepela in vode« (mentor: izr. prof. dr. Jure Marn); • dne 15. septembra 2016 Anton NAHTIGAL z naslovom: » Prispevek k razvoju računskega modela za katalitično pretvorbo povratnega toka izpušnih plinov « (mentor: prof. dr. Matjaž Hriberšek); • dne 15. septembra 2016 Matjaž FRAS z naslovom: »Varnost in zanesljivost čistilne naprave« (mentor: izr. prof. dr. Jure Marn); • dne 15. septembra 2016 Bernarda PODGORŠEK KOVAČ z naslovom: »Vpliv suhega zapolnjevanja sadre na koncentracijo delcev PM10 v zunanjem zraku« (mentor: prof. dr. Matjaž Hriberšek); • dne 26. septembra 2016 Alenka AMBROŽ JURGEC z naslovom: »Pouk mehanike s sodobnimi učnimi metodami in ob uporabi učnih okolij 21. stoletja« (mentor: prof. dr. Boris Aberšek); • dne 26. septembra 2016 Mitja SELIČ z naslovom: »Analiza sistemske plošče z difuzorjem v sistemu HVAC« (mentor: prof. dr. Jure Marn); • dne 26. septembra 2016 Boštjan VUZEM z naslovom: »Uravnavanje temperature v rastlinjaku s pomočjo toka tekočine med dvema steklenima ploščama« (mentor: izr. prof. dr. Jure MARN); • dne 29. septembra 2016 Andrej MEKIŠ z naslovom: »Razvoj pilotne naprave za optimizacijo kakovosti vode z uvajanjem ogljikovega dioksida« (mentor: prof. dr. Aleksandra Lobnik); • dne 29. septembra 2016 Mojca PLEVNIK ŽNIDAREC z naslovom: »Ugotavljanje izpolnjevanja pogojev odlaganja komunalnih odpadkov težke frakcije po mehansko biološki obdelavi« (mentor: prof. dr. Niko Samec); • dne 29. septembra 2016 Damijana ZLATOLAS z naslovom: »Priprava in karakterizacija elektroprevodnih tekstilnih nanosov« (mentor: prof. dr. Aleksandra Lobnik). Information for Authors All manuscripts must be in English. Pages should be numbered sequentially. The manuscript should be composed in accordance with the Article Template given above. The maximum length of contributions is 10 pages. Longer contributions will only be accepted if authors provide juastification in a cover letter. For full instructions see the Information for Authors section on the journal's website: http://en.sv-jme.eu . SUBMISSION: Submission to SV-JME is made with the implicit understanding that neither the manuscript nor the essence of its content has been published previously either in whole or in part and that it is not being considered for publication elsewhere. All the listed authors should have agreed on the content and the corresponding (submitting) author is responsible for having ensured that this agreement has been reached. The acceptance of an article is based entirely on its scientific merit, as judged by peer review. Scientific articles comprising simulations only will not be accepted for publication; simulations must be accompanied by experimental results carried out to confirm or deny the accuracy of the simulation. Every manuscript submitted to the SV-JME undergoes a peer-review process. The authors are kindly invited to submit the paper through our web site: http://ojs.sv-jme.eu. The Author is able to track the submission through the editorial process - as well as participate in the copyediting and proofreading of submissions accepted for publication - by logging in, and using the username and password provided. SUBMISSION CONTENT: The typical submission material consists of: - A manuscript (A PDF file, with title, all authors with affiliations, abstract, keywords, highlights, inserted figures and tables and references), - Supplementary files: • a manuscript in a WORD file format • a cover letter (please see instructions for composing the cover letter) • a ZIP file containing figures in high resolution in one of the graphical formats (please see instructions for preparing the figure files) • possible appendicies (optional), cover materials, video materials, etc. Incomplete or improperly prepared submissions will be rejected with explanatory comments provided. In this case we will kindly ask the authors to carefully read the Information for Authors and to resubmit their manuscripts taking into consideration our comments. COVER LETTER INSTRUCTIONS: Please add a cover letter stating the following information about the submitted paper: 1. Paper title, list of authors and their affiliations. 2. Type of paper: original scientific paper (1.01), review scientific paper (1.02) or short scientific paper (1.03). 3. A declaration that neither the manuscript nor the essence of its content has been published in whole or in part previously and that it is not being considered for publication elsewhere. 4. State the value of the paper or its practical, theoretical and scientific implications. What is new in the paper with respect to the state-of-the-art in the published papers? Do not repeat the content of your abstract for this purpose. 5. We kindly ask you to suggest at least two reviewers for your paper and give us their names, their full affiliation and contact information, and their scientific research interest. The suggested reviewers should have at least two relevant references (with an impact factor) to the scientific field concerned; they should not be from the same country as the authors and should have no close connection with the authors. FORMAT OF THE MANUSCRIPT: The manuscript should be composed in accordance with the Article Template. The manuscript should be written in the following format: - A Title that adequately describes the content of the manuscript. - A list of Authors and their affiliations. - An Abstract that should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, as well as the methodology employed. It should summarize the results and state the principal conclusions. - 4 to 6 significant key words should follow the abstract to aid indexing. - 4 to 6 highlights; a short collection of bullet points that convey the core findings and provide readers with a quick textual overview of the article. These four to six bullet points should describe the essence of the research (e.g. results or conclusions) and highlight what is distinctive about it. - An Introduction that should provide a review of recent literature and sufficient background information to allow the results of the article to be understood and evaluated. - A Methods section detailing the theoretical or experimental methods used. - An Experimental section that should provide details of the experimental set-up and the methods used to obtain the results. - A Results section that should clearly and concisely present the data, using figures and tables where appropriate. - A Discussion section that should describe the relationships and generalizations shown by the results and discuss the significance of the results, making comparisons with previously published work. (It may be appropriate to combine the Results and Discussion sections into a single section to improve clarity.) - A Conclusions section that should present one or more conclusions drawn from the results and subsequent discussion and should not duplicate the Abstract. - Acknowledgement (optional) of collaboration or preparation assistance may be included. Please note the source of funding for the research. - Nomenclature (optional). Papers with many symbols should have a nomenclature that defines all symbols with units, inserted above the references. If one is used, it must contain all the symbols used in the manuscript and the definitions should not be repeated in the text. In all cases, identify the symbols used if they are not widely recognized in the profession. Define acronyms in the text, not in the nomenclature. - References must be cited consecutively in the text using square brackets [1] and collected together in a reference list at the end of the manuscript. - Appendix(-icies) if any. SPECIAL NOTES Units: The SI system of units for nomenclature, symbols and abbreviations should be followed closely. Symbols for physical quantities in the text should be written in italics (e.g. v, T, n, etc.). Symbols for units that consist of letters should be in plain text (e.g. ms-1, K, min, mm, etc.). Please also see: http://physics.nist.gov/cuu/pdf/sp811.pdf . Abbreviations should be spelt out in full on first appearance followed by the abbreviation in parentheses, e.g. variable time geometry (VTG). The meaning of symbols and units belonging to symbols should be explained in each case or cited in a nomenclature section at the end of the manuscript before the References. Figures (figures, graphs, illustrations digital images, photographs) must be cited in consecutive numerical order in the text and referred to in both the text and the captions as Fig. 1, Fig. 2, etc. Figures should be prepared without borders and on white grounding and should be sent separately in their original formats. If a figure is composed of several parts, please mark each part with a), b), c), etc. and provide an explanation for each part in Figure caption. The caption should be self-explanatory. Letters and numbers should be readable (Arial or Times New Roman, min 6 pt with equal sizes and fonts in all figures). Graphics (submitted as supplementary files) may be exported in resolution good enough for printing (min. 300 dpi) in any common format, e.g. TIFF, BMP or JPG, PDF and should be named Fig1.jpg, Fig2.tif, etc. However, graphs and line drawings should be prepared as vector images, e.g. CDR, AI. Multi-curve graphs should have individual curves marked with a symbol or otherwise provide distinguishing differences using, for example, different thicknesses or dashing. Tables should carry separate titles and must be numbered in consecutive numerical order in the text and referred to in both the text and the captions as Table 1, Table 2, etc. In addition to the physical quantities, such as t (in italics), the units [s] (normal text) should be added in square brackets. Tables should not duplicate data found elsewhere in the manuscript. Tables should be prepared using a table editor and not inserted as a graphic. REFERENCES: A reference list must be included using the following information as a guide. Only cited text references are to be included. Each reference is to be referred to in the text by a number enclosed in a square bracket (i.e. [3] or [2] to [4] for more references; do not combine more than 3 references, explain each). No reference to the author is necessary. References must be numbered and ordered according to where they are first mentioned in the paper, not alphabetically. All references must be complete and accurate. Please add DOI code when available. Examples follow. Journal Papers: Surname 1, Initials, Surname 2, Initials (year). Title. Journal, volume, number, pages, DOI code. [1] Hackenschmidt, R., Alber-Laukant, B., Rieg, F. (2010). Simulating nonlinear materials under centrifugal forces by using intelligent cross-linked simulations. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, DOI:10.5545/sv-jme.2011.013. Journal titles should not be abbreviated. Note that journal title is set in italics. Books: Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication. [2] Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken. Note that the title of the book is italicized. Chapters in Books: Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages. [3] Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordić, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553576. Proceedings Papers: Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages. [4] Štefanić, N., Martinčević-Mikić, S., Tošanović, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427. Standards: Standard-Code (year). Title. Organisation. Place. [5] ISO/DIS 16000-6.2:2002. Indoor Air — Part 6: Determination of Volatile Organic Compounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva. WWW pages: Surname, Initials or Company name. Title, from http://address, date of access. [6] Rockwell Automation. Arena, from http://www.arenasimulation.com, accessed on 200909-07. EXTENDED ABSTRACT: When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters). The instruction for composing the extended abstract are published on-line: http://www.sv-jme.eu/information-for-authors/ . COPYRIGHT: Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to transfer copyright to SV-JME when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Transfer Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on http://en.sv-jme.eu/. PUBLICATION FEE: Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 240.00 EUR (for articles with maximum of 6 pages), 300.00 EUR (for articles with maximum of 10 pages), plus 30.00 EUR for each additional page. The additional cost for a color page is 90.00 EUR. These fees do not include tax. Strojniški vestnik - Journal of Mechanical Engineering Aškerčeva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu http://www.sv-jme.eu Contents Papers Afrasiab Raisi: Natural Convection of Non-Newtonian Fluids in a Square Cavity with a Localized Heat Source Zhenxing Kong, Dawei Pi, Xianhui Wang, Hongliang Wang, Shan Chen: Design and Evaluation of a Hierarchical Control Algorithm for an Electric Active Stabilizer Bar System Yongxiang Li, Xifan Yao, Jifeng Zhou: Multi-objective Optimization of Cloud Manufacturing Service Composition with Cloud-Entropy Enhanced Genetic Algorithm Rade VasiLjević, MiLomir Gašić, Mile Savković: Parameters Influencing the Dynamic Behaviour of the Carrying Structure of a Type H Portal Crane Cuneyt Uysal, Kamil Arslan, Huseyin Kurt: A Numerical Analysis of Fluid Flow and Heat Transfer Characteristics of ZnO-Ethylene Glycol Nanofluid in Rectangular Microchannels Xiaozhou Li, Manlin Zhu, Jiancang Xie: Numerical Simulation of Transient Pressure Control in a Pumped Water Supply System Using an Improved Bypass Pipe 55 9770039248001