_ Sinceri 955 r 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 380 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 Vice-President of Publishing Council Jože Balič University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: Tunnels are strategic structures that enable adequate, safer, and faster transport connections. The article deals with the effect of wind with its characteristics on the longitudinal ventilation of a road tunnel. The influence of wind gusts is researched using CFD simulations. Results present the importance of taking wind gustiness into account when designing and managing the tunnel ventilation as a segment of road tunnel safety. Courtesy: Sandi Radovan, DARS d.d. ISSN 0039-2480 International Editorial Board Kamil Arslan, Karabuk University, Turkey 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 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 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. © 2015 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 61, (2015), number 7-8 Ljubljana, July-August 2015 ISSN 0039-2480 Published monthly Papers Aleš Suban, Stojan Petelin, Peter Vidmar: Effect of Gusty Wind on Road Tunnel Safety 421 Qiang Cheng, Chengpeng Zhan, Zhifeng Liu, Yongsheng Zhao, Peihua Gu: Sensitivity-based Multidisciplinary Optimal Design of a Hydrostatic Rotary Table with Particle Swarm Optimization 432 Vlada Sokolović, Goran Dikić, Goran Marković, Rade Stančić, Nebojsa Lukić: INS/GPS Navigation System Based on MEMS Technologies 448 Regita Bendikiene, Gintaras Keturakis, Tilmute Pilkaite, Edmundas Pupelis: Wear Behaviour and Cutting Performance of Surfaced Inserts for Wood Machining 459 Jamshed Iqbal, Muhammad Imran Ullah, Abdul Attayyab Khan, Muhammad Irfan: Towards Sophisticated Control of Robotic Manipulators: An Experimental Study on a Pseudo-Industrial Arm 465 So Young Hwang, Naksoo Kim, Cheol-soo Lee: Numerical Investigation of the Effect of Process Parameters during Aluminium Wheel Flow-Forming 471 Adam Bureček, Lumfr Hružik, Martin Vašina: Determination of Undissolved Air Content in Oil by Means of Compression Method 477 Strojniški vestnik - Journal of Mechanical Engineering 61(2015)7-8, 421-431 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2433 Original Scientific Paper Effect of Gusty Wind on Road Tunnel Safety Aleš Suban* - Stojan Petelin - Peter Vidmar University of Ljubljana, Faculty of Maritime Studies and Transport, Slovenia This article deals with the effect of wind with its characteristics on the longitudinal ventilation of a road tunnel as a key segment of fire safety. Using CFD simulations, the influence of wind gusts, which differ and is more difficult to define than the influence of a constant wind on tunnel ventilation, is researched in more detail. On the basis of the simulation results, which have been validated with real measures, the findings are presented regarding the importance of taking into account the characteristics of non-stationary wind to ensure the adequate safety of tunnel users in the event of a fire, as well as of fire-fighting units during intervention. Keywords: road tunnel, Kastelec tunnel, longitudinal ventilation, wind gusts, the Bora wind, CFD simulations Highlights • The Bora wind characteristics that have an effect on road tunnel ventilation are defined. • Composition of CFD model in FDS by taking into account the wind characteristics is presented. For case study the Kastelec motorway tunnel is used. • Results present the importance of taking wind gustiness into account when designing and managing the tunnel ventilation as a segment of road tunnel safety. • For model validation real measurements of gusty wind Bora and air speed in Kastelec tunnel are used. 0 INTRODUCTION Tunnels are strategic structures that, predominantly in hilly terrain, enable adequate, safer, and faster transport connections. All newly constructed tunnels in Slovenia, as in the rest of the European Union (EU), that are longer than 500 m must be constructed in accordance with the EU Directive 2004/53/EC, which determines the minimal safety requirements for tunnels in the Trans-European Road Network. The directive specifies the minimal requirements for tunnels, which they must ensure the control of emissions that are created due to the vehicles, as well as control of the extraction of smoke and heat in the event of a fire. Mechanical ventilation must provide the adequate extraction of smoke and heat, visibility, decrease the spreading of the fire, as well as provide adequate conditions for fire-fighters, so that they are able to gain access to extinguish any fires and to rescue people from the tunnel [1] and [2]. Ventilation of unidirectional motorway tunnels is predominantly implemented longitudinally in the direction of the traffic flow, while some longer tunnels utilise combined ventilation [2] and [3]. Winds that blow on the tunnel's portal either help or hinder mechanical ventilation. When the wind is blowing in the opposite direction to the direction of the ventilation it can cause problems with ensuring adequate air velocities to adequately ventilate the tunnel, as it opposes the thrust of the ventilators. The effect of the wind on the safety of the road tunnel must be taken into account when designing the ventilation, as well as during its subsequent use. A number of tunnels in the Slovenian motorway network face the problem of wind effect; these are the Kastelec tunnel, the Dekani tunnel, and the Markovec tunnel. They latter lie in very windy areas where powerful gusts from the Bora wind blow, in addition to which, the direction of the portals is in the same direction as the wind. The Bora wind, due to its high velocity gusts, effects the ventilation of the tunnels [1]. When designing and finding solutions for the ventilation of tunnels that are affected by winds, it is necessary to take into account the wind's characteristics. In the ensuing article it shall be proven that, gustiness as a characteristic of the wind, should be taken into consideration when designing the ventilation. This paper uses the Kastelec tunnel as a case study where the issue of strong gusty wind which affects the tunnel's portal is dealt with. The ventilation in the Kastelec tunnel is longitudinal and is under effect of gusty Bora wind with speeds up to 25 m/s, which covers approximately 90 % of windy days in the year [1]. In rare cases, in approximately 8 % of days in the year, the Bora wind gusts reach speeds of up to 30 m/s or more [1]. In these cases, traffic flow through the tunnel is foreseen to be halted, as it's predicted that the effect of the wind on the ventilation is such that smoke and heat could not be properly extracted in the event of fire. *Corr. Author's Address: University of Ljubljana, Faculty of Maritime Studies and Transport, Pot pomorščakov 4, SI-6320 Portorož, Slovenia, a.suban@gmail.com 421 0.1 Adequate Ventilation in Road Tunnels as a Safety Segment Adequate ventilation that extracts smoke and heat in the event of a fire is essential for ensuring the safety of road tunnel users [3] and [4]. We need to ask ourselves, how can we determine that the tunnel ventilation is adequate in the event of a fire, when is effected by wind? In literature [4] and [5] is stated that the ventilation must extract the smoke and heat in the event of a fire to ensure the safe conditions for the users to exit the tunnel as well as to allow effective fire-fighting operations [6]. The ventilation must fulfil this requirement even under various meteorological influences [1] and [6]. Which is why, the tunnels that are affected by wind must have increased capabilities as well as adequate ventilation management so that the effect of the wind is eliminated. The release of heat and the creation of smoke depend on a fire's intensity. If the longitudinal ventilation velocity in unidirectional road tunnels in the direction of the traffic flow is too low, lower than that of the spreading smoke, the smoke will travel against the direction of the ventilation flow, also known as back-layering [2] to [5] and [7]. Back-layering endangers both the tunnel users and fire-fighters. Smoke contains numerous toxic combustion products, has a high temperature, and reduces visibility. It hinders, and in some cases makes successful interventions by fire-fighters impossible [6]. The critical velocity uc of the air in the tunnel is defined as the minimal wind velocity required to prevent the back-layering of smoke [2], [3] and [7]. The critical velocity depends on the geometry of the tunnel [2] and [7], and the size of the fire (HRR) [3]. Generally, the critical velocity in tunnels is approximately up to 2 m/s for smaller vehicle fires or initial fires (up to 10 MW), 2 to 3 m/s for fires of lorries carrying cargo (30 to 50 MW) and does not exceed 4 m/s for fires of tankers carrying fuel (above 100 MW) [1] and [3]. On the basis of this, we can establish that longitudinal tunnel ventilation ensures adequate safety for tunnel users and fire-fighter interventions if the ventilation speed is greater than the critical velocity of a certain intensity of fire in order to prevent backlayering [2] to [7]. As an example, if the possibility of an approximately 30 MW fire is foreseen, the tunnel ventilation must ensure an air velocity of above 3 m/s even under windy conditions. 1 WIND CHARACTERISTICS The effect of wind that has a constant velocity is not the same as that of gusty wind. It is easier to predict the tunnel conditions with wind with a constant velocity while also being easier to study its effect. The velocity hardly changes over time, making it possible to determine the pressure that the wind establishes on the tunnel's portal, which is a measure of the pressure that the mechanical ventilation must provide. On the other hand, the velocity of gusty wind varies considerably over time. At the peak of a gust, the velocity can reach 3 to 4 times the average, while a period of almost still air can immediately follow. The prevalent winds above the moderate geographical latitudes and in the altitudes over Europe are westerly winds. From a general westerly direction, the flow veers to the north and south, while large, closed spiralling winds occasionally break away. These are more common in the lower levels of the atmosphere and are known as cyclones and anticyclones. In relation to the weather conditions elsewhere in Europe, the winds in Slovenia are relatively weak, and when they are strong, they are both temporally and spatially restricted [8]. Regional winds are winds whose area of occurrence is in the range of approximately 100 km in size and are a consequence of general winds. In Slovenia, the Alpine and Dinaric mountain barriers have an effect on them. Close to the ground there are three distinct wind regimes: sirocco (a southerly wind), the Bora wind, and Karavanke foehn [8]. According to characteristics, winds are classified into steady, gusty and spiralling winds (whirlwinds) [8]. Steady winds blow at a constant velocity, which means they do not have gusts. The winds form and blow in flat areas where there are fewer obstacles and friction. They have smaller velocity oscillations during their constant average velocity, which, if present, are smaller and short-lasting with speeds around the average. Examples of such types of winds in Slovenia are the sirocco and thermal winds of the mistral, burin, etc. [8]. Whirlwinds are a rare phenomenon, occurring predominantly during severe storms. They are a weaker form of tornado. The wind direction is a spiral around the centre due to large differences in pressure, the velocities reached are therefore very high, while the surface area is small. In Slovenia, whirlwinds are very rare phenomenon, mainly occurring on the surface of the sea and lasting only a very short time [8]. 1.1 Gusty Winds A gust of wind is defined as a fast, strong, and abrupt rush of wind for a short time period which is usually followed by a period of still air [8] and [9]. Gusts are predominantly short-lived, lasting 20 to 50 seconds. A good example of gusty wind is the Bora wind which blows in the areas of SW Slovenia. The characteristics of the Bora wind are described below to get a sense of the gusty winds that are the subject of this study. Gusty winds form mainly due to the difference between air pressures in more hilly areas. The equalisation of these pressures leads to gusts forming. The Bora wind is cold, dry, and very gusty. A strong Bora wind is a regional phenomenon that is most prominent on slopes and just beneath them, as well as where converge and gather speed [8]. The Bora wind usually blows in a stable atmosphere (cold air pushes through under the warm air) which means that the local Bora wind speeds can greatly increase due to the convergence of winds between landforms, while there are also sheltered areas [8]. The Bora wind is dangerous due to its powerful gusts that effect traffic, structures and vegetation. The development of the Bora wind to its full intensity is relatively fast. In most cases, the Bora wind can develop from low intensity initial individual gusts to full intensity in 20 minutes on average; exceptional cases are also possible, where the Bora wind reaches its highest values in its initial gusts [9]. The average duration of the Bora wind is one day, but cases have been seen where it lasts for 12 days. One must be careful when using the average velocity of gusty wind like the Bora. An important characteristic of the Bora wind is its gustiness; measurements with a precise anemometer show almost still air between individual gusts [9]. Due to these periods of still air, the average velocity is low, however its gusts can reach three to four times its average speed [9], which is also confirmed by the measurements [10] and [11]. This is why one must take into account the high velocity of a gust when calculating the effects of gusty wind. From experience and measurements done in the north Adriatic, strong gusts of the Bora wind repeat in periods from 3 to 4 minutes (sometimes up to 11 minutes) [9]; measurements taken in the Karst region of Slovenia [10] prove the same. On the basis of analyses of Bora measurements [9] to [11] and its varying gust intensities (velocities), a simple and important Bora wind pattern could be created for research presented in this article (Fig. 1). The form of the gust is clearly visible from the pattern (the full line). A gust can emerge relatively quickly and then subside, it can also occur gradually and subside quickly, or occur gradually and subside gradually. Until the next gust, a period of almost still air occurs, or a period with very low wind velocities. Regardless of the form of the gust, we can conclude that is occurs quickly and subsides quickly. I 10 3 о 0 200 400 600 Time [s] — different shapes and duration of the wind gusts with usuai time gap of 3 min to 4 min — other possible time gaps between wind gusts Fig. 1. Gust patterns of the Bora wind The form itself is not an essential characteristic of the effect of a gust of wind. The characteristics that significantly affect the environment, traffic and, in our case, road tunnel ventilation are: • Gust intensity - the highest velocity reached by the wind during a gust. The Bora wind can reach speeds of up to 50 m/s in exposed areas [9] and [11]. • The duration of a gust - the length of time the gust lasts. Bora wind gust last between 20 s and 60 s [9] and [10]. The various gust durations are shown in Fig. 1 (the full line). The most frequent duration of gusts is around 30 s to 40 s, while shorter gusts are more frequent than longer ones [9]. • Period between individual gusts - the interval between gusts. Typical periods of stronger Bora gusts range from 3 min and 11 min [9] and [10]. 2 COMPOSITION OF THE CFD MODEL Using computer generated models, the tunnel with its characteristics and ventilation can be modelled. The boundary and initial conditions are used to model the effect of wind on the tunnel's ventilation. The use of simulation tools allows us to study the problem of the effect of winds on tunnel ventilation in more detail, taking into account that the person has to construct the model correctly. Researches described in this article define the significance of the problem, the approach to solving it, as well as the characteristic and specifics that are necessary to be considered. Simulation models that can be used to study and analyse the field of road tunnel ventilation can be roughly divided into zone models, one-dimensional (1D) models and three-dimensional (3D) models, to which apply the principle of computational fluid dynamics (CFD). CFD models are most suitable for carrying out studies regarding the effects of winds and other meteorological phenomenon on tunnel ventilation. Models that are based on CFD are being used more and more often. The numerical solving of the system of partial differential equations, which describe the physical behaviour of the flow give a better understanding of the underlying laws of fluid flow. Models can be based on various methods: the relatively simple RANS (Reynolds-averaged Navier-Stokes) and TRANS (transient RANS) models appropriate for very large systems like tunnels [12]; the large eddy simulation (LES) method that utilises turbulent models with which to deal with larger systems with good accuracy, so the method is very appropriate for the tunnels [13]; the most accurate model that is direct numerical simulation (DNS) which is an exceptionally demanding method to solve and it is only used for very small systems and not so much for the tunnels. The fire dynamics simulator (FDS) program which was used for research in this article is based on the LES method. The software was conceived at the American National Institute for Standards and Technology (NIST). It comes with the FDS -Smokeview (FDS-SMV) companion visualisation program that enables the graphic display of geometries and calculated parameters from the FDS program. FDS is a suitable program to explore the presented problem of the effect of winds and to carry out simulations. It enables CFD simulations, the input model with boundary pressure conditions with which the gust of wind can be modelled, and includes a combustion model to describe the fire. The Governing equations for FDS model are for conservation of mass (Eq. (1)), momentum (Eq. (2)) and energy (Eq. (3)) for a Newtonian fluid and are presented in [14] in detail: dp w - ''' + V • pu = m b dt д — (pu ) + V-puu + Vp = pg + f + V - t ij (1) (2) д Dp — (phs ) + V • phsu = D + f - q'b -V • q ' +e, (3) where m'b is the production rate of species by evaporating droplets or particles, t is time, p is density and u is the velocity vector (Eq. (1)). The term uu in momentum equation (Eq. (2)) is a dyadic tensor, g is gravity vector, fb represents external forces (i.e. drag exerted by liquid droplets) and т is the viscosity stress tensor. In energy conservation equation (Eq. (3)), hs is sensible enthalpy, q' is heat release rate per unit volume from a chemical reaction, q" is the energy transferred to the evaporating droplets, qb represents the conductive and radiative heat fluxes and e is the dissipation rate. Turbulence is treated by means of the Smagorinsky form of Large Eddy Simulation (LES) as a default mode of operation [14]. FDS uses rectangular grid elements and the core algorithm is an explicit predictor-corrector scheme that is second order accurate in space and time [14]. Hydrodynamic model solves numerically an approximate form of the Navier-Stokes equations, appropriate for low speed and thermally-driven flow. The approximation involves the filtering out of acoustic waves [14]. The literature [14] also contains a detailed description of the numerical solution of the momentum and pressure equations and for that, it is sufficient to consider the momentum equation written as [14]: ^ + f + VH = 0. dt (4) The pressure equation is obtained by taking the divergence of the momentum equation [14]: V2H = —(V-u )-V-F. dt ' (5) the For h=q 22+ /2 /p [14]. External pressure pexl outflow, pressure is defined by where q = | u | and p is set to pext on an open vent can be set by the user (DYNAMIC_PRESSURE; by default it is set to 0). For the inflow, when fluid enters the domain at an open boundary condition, it is assumed that the Bernoulli equation is valid and that the fluid on the boundary accelerates from the state along a streamline [14]. The FDS software is often used to solve problems in tunnels like smoke spread and air velocity in case of fire [15] and [16], reconstruction of tunnel fires to evaluate the conditions likely to be found in this type of fires [17], clarification of some fire phenomena, for example pulsation of a tunnel fire [18], and for numerical simulation of real scale fire tests [19]. FDS is validated by the developers of the software [20]. Various other researchers have also validated the software for different convection flows [21], for fires in buildings [22] as well as for tunnel fires [23]. When a simulation CFD model starts being constructed, it is important that the model is constructed gradually. Beginning with simple inputs, it then slowly evolves into a more complex model that describes the real situation. The optimum between accuracy, calculation time, and costs must be found. By taking into account the wind characteristics, it is possible to contribute to the planning of tunnel ventilation and increasing the safety level. By using CFD simulation tools, we are able to perform analyses in which the analyser must take certain rules into account in order to ensure the most realistic results. As already mentioned in the introduction, the research of wind issue is most easily presented using an actual case study of a tunnel, where we present the importance of taking the wind characteristics into account. In the case study it is shown what was considered when carrying out the simulations. Characteristics of the Kastelec tunnel will be used: the geometry, ventilation, and wind characteristics of the Bora wind that blows at the location of the tunnel. The tunnel is suitable as a case study, as the geometry of the tunnel is like the majority of newly built motorway tunnels in Slovenia, while similarities can be found to many other motorway tunnels in the EU. The same applies to the supervision, managing, and machinery of the tunnel. 2.1 Tunnel Geometry and Fire Modelling In the simulation model, the geometry must be properly designed (a good approximation of the real situation), adequate initial and boundary conditions must be set. The FDS enables the modelling of a fire, where it is possible to describe using a suitable combustion model or as a source of heat and mass. The second approach is usually more simple and suitable for the field dealt with. As we focus on the effect of winds on ventilation, the actual development of the fire is not significant, as the formation of smoke and heat release is greatest with a fully developed fire. If the ventilation adequately extracts the smoke and heat of a fully developed fire, it will also do so with a less intense fire. This is why we model the intensity of a fully developed fire in the tunnel. The intensity of the fire is predominantly determined using the heat release rate (HRR), measured in kW or MW. In the model, we determine the location of the fire and set it as the heat and substance source of certain intensity HRR of a fully developed fire. This is done by determining the HRR per unit area. The Kastelec tunnel is constructed as a two-lane twin-tube tunnel with pavements, passageways, and emergency lay-bays. Including the portals, the left tube - direction Koper - Ljubljana (KP-LJ) - which is used in this research is 2320 m long. The width of an individual tunnel tube is 9.60 m and the height is 6.5 m. Both tubes are connected with passageways every 400 m. By taking into account the horizontal course of the track, the maximum cross slope of the road is determined at 2.5 %. The geometry described is created in the simulation model. In simulation, a section of tunnel 400 m long (x-axis), 9 m wide (y-axis) and 6.5 m high (z-axis) is used (Fig. 2). There were 64800 cells used in numerical study. Mesh distribution in x-axis is 400 cells, in y-axis 18 cells and z-axis 9 cells (Fig. 2). By using a section of the tunnel, we can lower the number of cells, which drastically cuts the calculation time. Regarding the size of the section, it must not be too small, and it must enable the fluid flow to develop properly. If not, the results will not be accurate enough. Vehicles that may be present in the tunnel are not considered, as their type and number cannot be foreseen. Fig. 2. Mesh distribution and computational domain geometry Grid independence study was done for numerical study with finer cells throughout the domain (108000 cells) and the results for the ventilation speed in tunnel tube were almost the same through both simulations, only at few points the maximum differences were less than 1.2 %. This is a satisfactory outcome for the reduction of computational time and does not affect the findings of the final results. In the tunnel, forced longitudinal ventilation with axial fans is foreseen for the right and left tunnel tube. Since the direction of the prevailing wind - the Bora wind - coincides with the longitudinal position of the tunnel tubes, its effect is favourable on the tunnel tube leading from LJ to KP and unfavourable in the direction from KP to LJ. To maintain longitudinal ventilation in the direction of traffic in the KP-LJ tunnel pipe considering the wind velocity of 25 m/s in the opposite direction, 14 fans are required for normal operation, which are positioned in pairs every 100 m between fans and 150 m from portals (Fig. 3). Every jet fan provides 1100 N of thrust in the main direction. In exceptional cases, wind velocities of more than 30 m/s (108 km/h) can occur, which equate to 8 % of windy days in a year, where the closure of the tunnel is foreseen. It should be mentioned that in the ventilation project, the wind pressure on the portal at 25 m/s is foreseen to be constant; however the Bora wind blows in gusts. All ventilators are foreseen to operate at full power during periods of maximum wind velocities to overcome its effects. LJ-KP <^BORA WIND ( Ш Ш Ш Ш KP-LJi ш CD □ □ о о 1 □ ш iZZi Fan 1,2 3,4 5,6 7, S i, 10 11, 12 13, 14 Fig. 3. Layout of fans in the Kastelec tunnel 2.2 Boundary Conditions and Accounting for the Wind Boundary conditions are defined as pressures, temperatures, densities, velocities, and surface characteristics. The beginning and end of the tunnel geometry is defined as an open boundary condition, which does not represent an obstacle in the flow, and the remaining surfaces as walls. For air (gas) movement in the tube pressure boundary conditions are used: on the one side, pressure created by the fans, on the other side, pressure caused by a gust of wind. Pressure loss that occurs along the length of the tunnel is taken into account. Boundary condition 1 is the pressure exerted on the starting point of the section marked as p1 and is the pressure caused by all the fans in the tunnel, less the losses in the tunnel tube. Boundary condition 2 is the pressure on the final point of the geometry, marked with p2 and is the pressure caused by the wind (Bora) blowing in the opposite direction to the ventilation, less the losses in the tunnel tube (Fig. 4). Fig. 4. The tunnel section with pressure boundary conditions Equation for calculating the pressure for boundary condition 1: Pi = Pa + nf • Pf- pia (6) where nf is the number of fans,pf is the pressure of an individual fan, whereby it is considered that all fans provide the same pressure, the equation is presented in [3]; ploss are pressure losses in the tube calculated using the Darcy equation. The surface roughness (friction coefficient) is accounted for in the losses, which also takes into consideration the various construction obstacles present in the tunnel. Stationary vehicles are not taken into consideration as their type and number are difficult to predict. Although, stationary vehicles would represent an obstacle in the flow of both the mechanical ventilation and the wind, which would mean that it would affect both boundary conditions. Atmospheric pressure pa needs to be considered in the event it differs at both ends of the tunnel tube. This primarily occurs in longer tunnels. If the atmospheric pressure is equal at both ends of the tunnel, it can be left out of the calculation of both boundary conditions. Equation for calculating the pressure for boundary condition 2: p2 " pa + pw - pwloss , (7) where pw is dynamic pressure caused by wind, pwloss is the loss of pressure in the tube during various wind velocities at the portal. Since gusty wind is being considered, the speed on the portal is not constant but varies over time, which must be accounted also for losses. The same is valid for atmospheric pressure pa as is for boundary condition p1. The dynamic wind pressure pw on the portal is calculated using the Eq. (8): A, = ZPUw, (8) where Z is the coefficient of local loss regarding the shape of the entrance of the tunnel tube (portal), the values range between 0.5 and 1.0 ([3] and [24]) during winds that work directly in the direction of the tunnel, perpendicularly to the portal. In cases where the wind is not directly perpendicular to the portal, the coefficient can be greater than 1. p is the air density, uw is the wind velocity at the tunnel portal. Very similar equations for boundary conditions are presented in literature [25], only this ones are moderated for the system dealt with. The maximum gust velocity of the Bora wind used in the simulation was 30 m/s. In view of the gust characteristic (Fig. 1), an input model was constructed for boundary condition 2 (p2). The velocities were recalculated into pressures which were then entered into the simulation input model (Fig. 5). Calculated and used pressures for boundary condition 1 in simulation was 46 Pa, and for boundary condition 2 at maximum gust speed 30 m/s was 59 Pa. It is important to know that the duration of gust's maximum speed is just a few seconds. For the example we take a gust of 40 s which was used in one of the simulation scenarios: the maximum speed emerges quickly in 1 s, lasts for 10 s, then in next 18 s drops to 83 % of maximum speed, in another 10 s drops to 66 % and after that drops quickly in 2 s to the low wind velocity of 27 % of maximum speed until the next gust appears. For longer gusts, percentages were the same, just to each segment of the time additional few seconds were added. Fig. 5. Example of Input model for boundary conditions px and p2 In scenario with two wind gusts of high speed 3 TUNNEL VENTILATION AND EFFECT OF GUSTY WIND 3.1 Simulation Scenarios Using numerous simulations of the effect of wind of tunnel ventilation, the field was carefully studied. Many combinations were simulated, using various wind velocities, different amounts of gusts, various forms of gusts, gust durations, periods between gusts, and various fire intensities. Models using constant winds were constructed for comparison. It immediately became clear that the gust's form does not have a significant function, which is why all gusts were modelled in accordance with the p2 boundary condition in Fig. 5. On the basis of all the scenarios, a few key ones were chosen and presented. All scenarios also simulate a fire of 10 MW, except Scenario 6 in which the fire size is 15 MW. The key wind scenarios were the following: • Scenario 1: 2 gusts of wind at 30 m/s, duration of gust 40 seconds, period between gusts 3 minutes (usual Bora wind pattern). • Scenario 2: 3 gusts of wind at 30 m/s, duration of gust 40 seconds, period between gusts 1 minute (very rare Bora wind pattern). • Scenario 3: 2 gusts of wind at 30 m/s, duration of gust 50 seconds, period between gusts 1 minute (rare Bora wind pattern). • Scenario 4: 2 gusts of wind, 1st at 30 m/s, 2nd at 27 m/s, duration of gust 50 seconds, period between gusts 1 minute (usual Bora wind pattern). • Scenario 5: effect of constant wind at 21 m/s (comparison). • Scenario 6: 2 gusts of wind at 30 m/s, duration of gust 50 seconds, period between gusts 1 minute, fire HRR is 15 MW (comparison between fire sizes). 3.2 Results To determine the effect of wind on the tunnel's ventilation, six scenarios were constructed which were then analysed using FDS simulation tools. The most important characteristics of gusty wind were divided in the scenarios, such as the period between gusts, length of the gust and its intensity. To compare effects, scenario with a constant wind was also constructed. The results showed that the effect of gusty wind differs from the effect of constant wind. Results has shown that even in the event when the pressure of the gust on the portal is somewhat greater than the mechanical ventilation pressure, the extraction of smoke and heat from the tunnel may be suitable because of shortness of the gust. This is a consequence of the air flow inertia due to mechanical ventilation. In Fig. 6, the speed of velocities in the tunnel tube in various scenarios is shown, as well as for critical velocities uc for particular fire intensities. In the most common Bora wind pattern (Scenario 1), the speed of the tunnel ventilation remains above 4 m/s, which ensures the adequate extraction of smoke and heat for modelled fire. In the rare occurrences of more than one strong consecutive gusts (Scenario 2) between which the period is short (1 minute), the interaction of gusts can be seen as well as a greater drop in ventilation velocity. Despite this, the extraction of smoke and heat is appropriate for lower to medium intensity fires. The gusts are predominantly short-lived (to 40 s), in rare cases they can last longer and vary between 50 s and 60 s. Scenario 3 describes the occurrence of longer gusts of Bora wind, lasting 50 s with a period of 1 minute. In this case, speed in tunnel tube drops at lowest point of all the scenarios, but still remains above 3 m/s, which ensures the adequate extraction of smoke and heat at 10 MW fire. Real cases of two or three consecutive gusts of exactly the same intensity are very rare with the Bora wind. That is why a Time [s] -» Scenario 1 Scenario 2 Scenario 3 —Scenario 4 Scenario 5 Fig. 6. The results of various wind scenarios for a 10 MW fire comparison of a more realistic situation was carried out (Scenario 4), where the consecutive gusts differed slightly regarding velocity. In Scenario 4, a model was constructed of the effect of two longer gusts (50 s) with a period of 1 minute, and with varying speeds. The second gust was 3 m/s slower that the first. From the result in Fig. 6, there is a visible difference of the effect on the speed in the tunnel in comparison with Scenario 3 in which speed of both wind gusts is equal. In Scenario 5, constant wind speeds that affect the portal is taken into consideration. The velocity in the tunnel does not fluctuate when affected by constant wind on the portal. Therefore, it is visible in Fig. 5 that even at a much lower velocity of 21 m/s, the critical velocity of the 10 MW fire is almost ____________ / ^ \ // i / ** / * // 1 /л r \ AS / / / / / / ' \ h / \ U l / ' \ /' uc: > 100 MW f W Л ji j f uc\ 30 to 50 MW / V/ / uc: < 10 MW Time [5] -Scenario 3 — -Scenario 6 Fig. 7. The difference between 10 MW and 15 MW fire reached. At higher speed than 21 m/s, the tunnel can be established to not provide adequate smoke and heat extraction, as the velocity in the tunnel drops below the critical velocity of a 10 MW fire, that is under 2 m/s. In Scenario 6, the wind characteristics are the same like in Scenario 3, but fire has a larger HRR of 15 MW. Results are seen in Fig. 7. The speed in Scenario 6 drops slightly under critical velocity of 15 MW fire for a short period of time, which causes a very short back-layering. The entire occurrence in tunnel lasts 70 s and is quickly eliminated by the mechanical ventilation (Fig. 8). The smoke remains stratified during this period, while the velocity in the lower part of the tunnel where the users, and later firefighters, are located adequately extracts the smoke and heat (Fig. 8). With comparison between Scenario 3 and Scenario 6 we can see that is important to determine to what size of the fire the ventilation fulfil the safety requirement for a tunnel under various wind influences. Time 364 s Fig. 8. Short-term back-layering in Scenario 6 4 VALIDATION OF THE RESULTS WITH MEASUREMENTS Verification of the constructed theoretical scenarios was created using real measurements. The Motorway Company of the Republic of Slovenia (DARS) , which is the manager of the Slovenian motorway network, carries out meteorological measurements in the area of the Kastelec tunnel, more specifically on the Smelavc viaduct [11]. In the Kastelec tunnel itself, air velocity meters are fitted in the tunnel pipe [26]. The wind speed sensor (Boschung Mecatronic WG/WR) accuracy is ±3 % in measuring range from 0.5 m/s to 50 m/s, the tunnel tube air speed sensor (SICK Flowsic 200) typical accuracy is ±0.1 m/s in measuring range -20 m/s to +20 m/s. A theoretical CFD validation model was constructed on the basis of real measurements [11] and [26]. The CFD model characteristics were the same as in the tunnel at the time of measuring: wind in gusts of 30 m/s, the tunnel closed for all traffic (empty tunnel tube), the mechanical ventilation is completely shut down. From the measurements [11] it is not clear how long an individual gust lasts, which is why three simulations were constructed with gust durations of 40 s (Sim 1), 50 s (Sim 2) and 60 s (Sim 3), which are the most common lengths of Bora wind gusts [9] and [10]. The results are presented in Table 1. The difference between the theoretical CFD models and real measurements from the tunnel is around 5 %, while almost being the same in certain cases (simulation 2). It is very difficult to postulate all parameters (gust characteristics, coefficients, losses, etc.) of a theoretical model that affect the air velocity in the tunnel. As the difference in results of the simulations and real measurements is about 5 %, the theoretical CFD models can be considered a very close approximation to a realistic situation. The other CFD models used in the research scenarios are based on the same assumptions, which is why we can claim that the deviations from a real situation are around 5 % too. Table 1. Comparison of results from FDS simulations and real measurements [23] and [26] - model validation The speed of air movement in the tunnel tube under the influence of 30 m/s gusty wind on the portal FDS Kastelec tunnel Difference result measurements in speed Difference [m/s] [m/s] [m/s] [%] Sim 1: 40 s -5.20 -5.47 0.27 4.9 Sim 2: 50 s -5.48 -5.47 0.01 0.2 Sim 3: 60 s -5.76 -5.47 0.29 5.3 5 CONCLUSION Road tunnels are strategic structures for the transport of people and goods. Fire is very hazardous for tunnel users, and likewise for fire-fighters who attempt to extinguish the fire. A key element of safety for users and for a safe intervention is the adequate ventilation of tunnels, which must ensure the extraction of smoke and heat. When designing the ventilation system it is necessary to consider various factors, one of the more important factors being meteorological effects. A key meteorological factor that affects ventilation is wind. Wind, which blows on the portal of unidirectional road tunnels with longitudinal ventilation in the opposite direction of the traffic counteract the ventilation. Mechanical ventilation must provide a suitable air velocity in the tunnel to overcome friction and wind pressure exerted on the portal. It is important that the mechanical ventilation provides air velocity in the tunnel that is greater than the fire's critical velocity, which ranges from 1.5 m/s for fires up to 10 MW, and up to 4 m/s for fires over 100 MW. This article establishes the importance of taking into account the gustiness of winds when determining the adequate safety of road tunnels. As a case study, research of the Kastelec tunnel was carried out, in an area where the gusty wind of the Bora blows. Research has shown that the effect of gusty wind, such as the example of the Bora wind, on the speed in a tunnel differs from that of a constant wind. It is necessary to consider a gust's intensity, duration, and period of repetition. The velocity at the tunnel portal during the gust fluctuated depending on its curve. This effect also transfers to the movement of air in the tunnel. As the occurrence of the maximum speed of the gust is only short-term, followed by a period of low velocity wind, the mechanical ventilation air flow inertia prevents a greater drop in speed in the tunnel. By comparing the effect of constant wind, a significant difference between the wind velocities at which the tunnel air flow velocity is still adequate is determined. In all wind scenarios where the gusts were at 30 m/s, the ventilation velocity in the tunnel was adequate for the determined 10 MW fire and provided the extraction of smoke from the tunnel. Theoretical CFD models of the effect of wind have been validated with real wind measurements in the area of the Kastelec tunnel [23] and measurements from the tunnel tube [26]. Verification indicated up to a 5 % deviation in the results. On the basis of this, designed CFD models can be categorised as a good approximation of reality. 6 ACKNOWLEDGEMENT Authors thank DARS d.d. The Motorway Company of the Republic of Slovenia, especially Mr. Marko Korošec for providing measurement [11] and [26] and Mr. Boris Milič who helped with other information about motorway tunnels in Slovenia. 7 REFERENCES [1] Petelin, S., Vidmar, P., Perkovič, M., Paliska, D., Klasek, Z., Filli, B., Polman, Z., Čufer, A., Malovrh, M., Gerzevič, R., Resman, N., Rom, J. (2005). Fire Safety Field Research Paper: Success of an Intervention in Tunnels. Republic of Slovenia, Ministry of Defence. (in Slovene) [2] Drakulič, M. (2006). The Role of the Longitudinal Ventilation System in Road Tunnel Fire Accidents. Doctoral dissertation, Faculty of Mechanical Engineering and Naval Architecture, Zagreb. (in Croatian) [3] Tunnel Study Centre. (2003). Ventilation. Road Directorate, Pariz. (in French) [4] Modic, J. (2003). A Model of a Tunnel and a Simulation of Ventilation in the Case of Fire. Strojniški vestnik - Journal of Mechanical Engineering, vol. 49, no. 9, p. 458-468. [5] Rhodes, N. (2011). Operational Strategies for Emergency Ventilation. World Road Association - PIARC, Report no. 2011R2. [6] Suban, A. (2012). Analyses to Support Intervention in the Event of a Fire in Road Tunnel. M.Sc. Thesis, Faculty for Maritime Studies and Transport, University of Ljubljana, Portorož. (in Slovene) [7] Wu, Y. (2010). The critical velocity and the Fire Development. Fourth International Symposium on Tunnel Safety and Security Frankfurt Conference Proceedings, p. 407-417. [8] Rakovec, J., Žagar, M., Bertalanič, R., Cedilnik, J., Gregorič, G., Skok, G., Žagar, N. (2009). Winds in Slovenia. ZRC, Ljubljana. (in Slovene) [9] Petkovšek, Z. (1987). Main bora gusts - a model explanation. Geofizika, vol. 4, p. 41-50. [10] Perkovič, M., Batista, M., Najdovski, D., Vidmar, P., Luin, B.. (2011). Measuring the effect of crosswinds on dynamics of road vehicles, focussing on the Bora analyses. ICTS Portorož Conference proceedings. [11] DARS (2014). Measurements of Bora wind on Smelavc viaduct (45°34'52.32"N, 13°54'34.60"E) - March 4 and March 5, 2014. Company database archive, Ljubljana. [12] Muhasilovic, M., Duhovnik, J. (2012). CFD-Based Investigation of the Response of Mechanical Ventilation in the Case of Tunnel-Fire. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 3, p. 183-190, D0l:10.5545/sv-jme.2009.091. [13] Rahmani A., Carlotti P., Gay, B., Buffat, M. (2004). Simulating fires in tunnels using large eddy simulation. International Conference Tunnel Safety and Ventilation, Graz Conference Proceedings, p. 111-118. [14] McGrattan, K., Hostikka, S., Floyd, J., Baum, H., Rehm, R., Mell, W., McDermott, R. (2010). Fire Dynamics Simulator (Version 5) Technical Reference Guide, Volume 1: Mathematical model. NIST/VTT, Department of Commerce, Washington. [15] Weisenpacher, P., Halada, L., Glasa, J. (2011). Computer simulation of fire in a tunnel using parallel version of FDS. MCS7 Conference Proceedings. [16] Yan-ling, L., Xin-feng, L., Ping, L. (2009). FDS analysis of smoke spreading in the tunnel with flue under fire. Journal of Chemistry & Chemical Engineering, vol. 3, no. 4, p. 39-45. [17] Huczek, J., Mintz., T., Bajwa, C., Das, K., Axler, K. (2010). FDS simulation of the newhall pass tunnel fire. 4th International Symposium on Tunnel Safety and Security, Conference Proceedings. [18] Kim, M.E. (2006). A Study on Pulsation in Runehamar Tunnel Fire Tests with Forced Longitudinal Ventilation. M.Sc. Thesis, Worcester Polytechnic Institute, Worcester. [19] Li, Y.Z., Ingason, H., Lönnermark, A. (2012). Numerical simulation of Runehamar tunnel fire tests. 6th International Conference Tunnel Safety and Ventilation. Conference Proceedings, p. 203-210. [20] McGrattan, K., Hostikka, S., McDermott, R., Floyd, J., Weinschenk, C., Overholt, K. (2014). Fire Dynamics Simulator Technical Reference Guide Volume 3: Validation. NIST/VTT, Department of Commerce, Washington. [21] Smardz, P. (2006). Validation of Fire Dynamics Simulator (FDS) for Forced and Natural Convection Flows. M.Sc. Thesis, University of Ulster, Ulster. [22] Nielsen, J.G. (2013). Validation Study of Fire Dynamics Simulator. M.Sc. Thesis, Department of Energy Technology, Aalborg University, Aalborg. [23] Vidmar, P. (2007). Deterministic Model of Fire in Tunnel. Ph.D. Thesis, Faculty for Maritime Studies and Transport, University of Ljubljana, Portorož. (in Slovene) [24] Kraut, B., Puhar, J., Stropnik, J. (2001). Kraut's Engineering Handbook, 13th edition. Littera picta, Ljubljana. [25] Vidmar, P., Petelin, S. (2003). An analysis of a fire resulting from a traffic accident. Strojniški vestnik - Journal of Mechanical Engineering, vol. 49, no. 5, p. 254-266. [26] DARS (2014). Measurements in Kastelec tunnel tube - March 4 and March 5, 2014. Company database archive, Ljubljana. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)7-8, 432-447 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2478 Original Scientific Paper Sensitivity-based Multidisciplinary Optimal Design of a Hydrostatic Rotary Table with Particle Swarm Optimization Qiang Cheng12 - Chengpeng Zhan1 - Zhifeng Liu1* - Yongsheng Zhao1 - Peihua Gu3 1 Beijing University of Technology, China 2 Huazhong University of Science and Technology, Digital Manufacturing Equipment and Technology Key National Laboratories, China 3 Department of Mechatronics Engineering, Shantou University, China In five-axis machine tools, a rotary table is often used as a means for providing rotational motion and supporting the workpiece. Its rigidity, precision and carrying capacity is directly related to the machining ability and the accuracy of the NC machine tool. Traditional rotary table design is normally performed by teams, each with expertise in a specific discipline, which causes excessive iterations and cannot provide users with products of reliable working performance and bearing capacity. To achieve an optimal design with less cost and better performance, this paper considers the mutual interaction of hydrostatics and structure disciplines involved in the design of hydrostatic rotary tables, and a sensitivity-based multidisciplinary optimal design procedure of a hydrostatic rotary table is proposed. Sensitivity analysis is adopted to identify the key design parameters that have a major influence on the performance of rotary tables to improve the convergence of optimization process. The constrained multi-objective optimization problem is solved by using a particle swarm optimization approach. A hydrostatic rotary table of a five-axis heavy duty machine tool is selected as an illustration example. The results show that the proposed method can realize the multidisciplinary optimization resulting in a rotary table of good rigidity and bearing capacity. Keywords: hydrostatic rotary table, sensitivity analysis, NC machine tool, particle swarm optimization, multidisciplinary optimal design Highlights • The mutual interaction of hydrostatics and structure disciplines is considered in the design of a hydrostatic rotary table. • A sensitivity-based multidisciplinary optimization design method of a hydrostatic rotary table is proposed. • Sensitivity analysis is adopted to identify the key design parameters. • The constrained multi-objective optimization problem is solved by using a particle swarm optimization approach. • A hydrostatic rotary table of a five-axis heavy duty machine tool is selected as an illustration example. 0 INTRODUCTION Numerical control (NC) machine tools are one of the most important components in modern manufacturing facilities, and high-performance machines are required. The accuracy of the machine tool motion has a significant influence on the quality of the machining operations and the machined parts. The industrial demand for manufacturing geometrically complex parts often calls for multi-axis machine tools to have a tool orientation capability [1]. Five-axis machine tools are becoming increasingly popular and can be found in a large number of manufacturing applications [2]. However, such tools are usually more complicated and less rigid in structure in comparison with traditional three-axis machine tools, which leads to lower machining accuracy [1]. In five-axis machine tools, a rotary table (or rotary-tilting table) is often used as a means for providing rotational motion and is now widely used in the machine shop. As a consequence, a resultant volumetric error, i.e. the relative error between the cutting tool and the machined part, is more complicated than three-axis machine tools. Hence, the reduction of the error to improve the accuracy of the machine tool is crucial. In a five-axis machine tool, each element contributes some degree of inaccuracy due to their manufacturing and assembling limitations, whereas machine structure stiffness, machine foundation, machine control, operating systems, and environmental conditions further add some inaccuracy. The errors/ inaccuracies can be reduced with the structural improvement of the machine tool through better design, manufacturing and assembly practices [3]. In heavy duty 5-axis machine tools, a hydrostatic rotary table handles supporting and rotating the workpiece, and the location of a rotary axis constitutes a significant error source [4]. Therefore, the accuracy of the rotary table is crucial for part manufacturing with multi-axis machine tools [1]. Its rigidity, precision and carrying capacity is directly related to the machining ability and the accuracy of the NC machine tool. Due to the lack of effective means of analysis and experiment, the designed rotary table cannot 432 *Corr. Author's Address: Beijing University of Technology, No.100, Pingleyuan, Chaoyang District, Beijing, China, lzfeng1@126.com provide users with reliable working performance and bearing capacity, which directly affects the machining accuracy, and the bearing grounding will cause serious economic loss to users. Therefore, how a hydrostatic rotary table is designed for a heavy duty five-axis machine tool is very important for machine manufacturers, which requires determining the relations between the design variables that refer to different disciplines and their effects on various performance objectives. Traditional engineering design is normally performed by teams, each with expertise in a specific discipline, such as hydrostatics or structures. Each team uses the experience and judgment of its members to develop a workable design, usually sequentially. This causes contradictory results from different design teams; therefore, multiple revisions throughout the design process may result [5]. These excessive iterations will clearly increase the cost and time of the design process. Starting in the 20th century, engineers decided to use multidisciplinary design optimization (MDO) approaches to solve similar problems [5]. MDO is widely studied and applied in both academia and industry, and can reduce the time required to execute the design process. By using MDO methods, designers may quickly and efficiently conduct alternative design points over a wide range of parameters. MDO has become essential for solving complex engineering design problems. MDO processes allow an evaluation of the constraints for multiple disciplines from the early stages of the design; the expense of making approximations or corrections is thus reduced. Although the multidisciplinary optimization has been successfully applied in the aircraft industry and resulted in more reliable and better products, it is rarely used in the design of machine tools. To make a high-performance machine, not only can the kinematic functions of mechanisms used in the machine be optimized but also the structure and even the controllability of the mechanisms must be optimized. Therefore, the MDO technique can be used to improve the overall performance of machines. Based on these reasons, the aim of this paper is to present a multidisciplinary optimization method for a hydrostatic rotary table based on sensitivity analysis. The main contribution of this paper is two-fold. Firstly, the proposed multidisciplinary optimization method takes into comprehensive consideration the hydrostatics and structure disciplinary characteristics, which can reduce the iterative modification caused by sequential design by experts in different disciplines. Secondly, because a hydrostatic rotary table is a complex product and has many design parameters, sensitivity analysis is introduced to identify the key design parameters that have significant influence on the performance; therefore, the optimization is realized with quick convergence. The remainder of this paper is organized as follows. In next section, the background reviews of sensitivity-based MDO are given. The third section presents multidisciplinary characteristics of modelling in both the hydrostatics and structure disciplines. Subsequently, MDO based on sensitivity analysis is implemented in Section 4. Section 5 demonstrates the proposed method with a designed hydrostatic rotary table. The final section contains the conclusions. 1 BACKGROUNDS 1.1 Multidisciplinary Design Optimization MDO is essential to the design and operation of a complex system, because it simultaneously takes into consideration all relevant disciplines to find the global optimum that is superior to a solution from a sequence of local optimizations in individual disciplines. MDO can be traced back to Schmit [6] and Haftka [7], who extended their experience in structural optimization to include other disciplines. The intention was to address these challenges and, in particular, the coupling within design hierarchies and between disciplines [8]. The research area of MDO has been intensively investigated over the previous decades [8], and the focus of MDO has shifted dramatically, as faculties and researchers are finding new ways to use MDO methods and tools on a wide array of problems [9] and [10]. There have been many advances to capture, represent, and propagate couplings in analysis, design, and organizations, yet the design of complex engineered systems continues to be challenging. MDO can enhance system design by exploiting synergies among different disciplines. However, there are two major challenges in applying MDO: organizational and computational complexities. The organizational complexities mean that a simultaneous consideration of multiple disciplines may increase the difficulty of data origination and coordination between different disciplines or different computer-aided engineering (CAE) software. The increased complexity of the optimization system inevitably causes increased amounts of computing, and the convergence problem increases the computational complexities [11]. To address these two challenges, one of the research focuses in MDO has been on optimization procedure [9] and [10]. Optimization procedures can be categorized into two types: single-level and multi-level approaches. Single level approaches employ a system optimizer for the whole problem, which is straightforward to understand and easy to implement. Multi-level approaches utilize decomposition strategies to allow disciplinary autonomy in design and optimization while managing interdisciplinary consistency via system coordination [12]. One of the first applications of MDO was aircraft wing design, where aerodynamics, structures, and controls are three strongly coupled disciplines [13]. Since then, the application of MDO has been extended to complete aircraft [14] and a wide range of other engineering systems, such as bridges [15], buildings [16], automobiles [17] and [18], ships [19], and spacecraft [20]. Yifei et al. [21] proposed the framework of multidisciplinary energy-saving optimization design for bridge cranes, and realized metal structures level, transmission design level, and electrical system design as well as the optimization design of bridge cranes. There have been many summaries of MDO since the 1990s. In a collection of articles Kroom [22] provided a comprehensive overview of MDO, including a description of both monolithic and distributed architectures. Sobieszczanski-Sobieski and Haftka [23] presented a detailed summary of the MDO literature. Because one of the most important considerations when implementing MDO is how to organize the discipline analysis models, approximation models, and optimization software in concert with the problem formulation, a combination of problem formulation and organizational strategy is referred to as an MDO architecture. Martins and Lambe [24] provided a survey of all the architectures that had been presented at length in the literature. In engineering design, to achieve high reliability and safety in complex and coupled multidisciplinary systems, reliability-based multidisciplinary design optimization (RBMDO) has received increasing attention. Since the 1990s, the consideration of the effect of uncertainty has been one of the focus areas in engineering design [25] and [26]. RBMDO can improve the system design by exploiting the synergistic effects between coupled disciplines by interdisciplinary collaboration, and can also enhance the reliability by taking uncertainties into consideration in the design phase [27]. There are successful applications of RBMDO clearly demonstrating its efficacy [28]. Yao et al. [27] summarized two categories of RBMDO procedures: the single level procedure, and the decomposition and coordination-based procedure, which were mostly developed under random uncertainties with probability theory [8]. 1.2 Hydrodynamics Design and Analysis The design of hydrostatic rotary tables has been studied by researchers in recent years, mostly concentrating on the hydrostatic part. Traditional optimization research work is mostly based on single objective, such as pump power, friction factor, bearing-capacity [29] to [31], and many researchers have used this method to design the hydrostatic bearings [32] to [35]. Lin studied the influence of factors including surface roughness and inertia, and optimized the bearing [34]. Zhao et al. used Isight software to optimize hydrostatic guideways with multiple pockets for a heavy duty CNC vertical turning mill [36]. Solmaz compared single and multi-objective optimization solutions of hydrostatic radial bearings and thrust bearings [37]. Several researchers have intensively studied the optimization of hydrostatic bearings. As for the rotary table, because it has a huge volume and many rib plates, and because there are many design parameters, it is difficult to realize the multidisciplinary optimization quickly and with good convergence. Fortunately, sensitivity analysis can evaluate the variation in dynamic model outputs with respect to variation in model parameters. Therefore, sensitivity analysis needs to be taken into account in order to weigh all of the parameters to get a more accurate optimization results. 1.3 Sensitivity Analysis Sensitivity analysis (SA) can be used to identify the effect of system parameter uncertainty variation on system responses and to identify the most critical parameters [38] and [39]. It is one approach to identifying and quantifying the relationships between input and output uncertainties [40], and can evaluate the variation in dynamic model outputs with respect to variation in model parameters. Therefore, SA can be used to perform uncertainty analysis, estimate model parameters, analyse experimental data, guide future data collection efforts, and suggest the accuracy to which the parameters must be estimated [41]. For a review on methods for SA, see Saltelli et al. [42], and Helton and Davis [43]. SA is divided into the local sensitivity analysis (LSA) and global sensitivity analysis (GSA) [44]. LSA, emphasizing the effect of small parameter variations on model responses, is used to determine model response changes with an individual parameter. GSA is applied to understand how the model response varies with the model parameters to determine interaction strengths among the parameters, such as Fourier amplitude sensitivity test, regression-based methods, Sobol method, and McKay's oneway ANOVA method [45]. The global SA examines the global response (averaged over the variation of all the parameters) of model output(s) by exploring a finite (or even an infinite) region. The local SA, easier to implement, can only inspect one point at a time, and the sensitivity index of a specific parameter is dependent on the central values of the other parameters. There are several numerical methods for the calculation of LSA, e.g. finite differences [44], iterative approximation with directional derivatives [46] or the direct decoupled method [47], but the calculated values should be identical within the numerical accuracy of the method used. Common GSA techniques include correlation and regression modelling, variance decomposition analysis [48], factorial screening [49], and partitioning based generalized sensitivity analysis. Saltelli et al. provided an extensive list of other techniques that have also been found useful in this context [44]. 2 MULTIDISCIPLINARY CHARACTERISTICS MODELING 2.1 Structure of Hydrostatic Rotary Table According to engineering experience, a heavy hydrostatic rotary table whose diameter is more than 5 metres will have double support circles, as shown in Fig. 1. Generally, a hydrostatic rotary table consists of five parts: the oil supply system, the drive system, the countertop, the base system, and the support system. The countertop and support system are the most important parts of the rotary table, which are needed to focus on during the design process. The support system includes supporting oil pads, preloaded oil pads and the radial bearing. The supporting oil pad is a circular step pad, and a rotary table has 24 supporting pads in total. They are arranged in two supporting circles and the number of pads in the second supporting circle is twice of that of the first supporting circle. The preloaded oil pad is annular step recess pad and provides a pre-pressure that can enhance the stiffness of the turntable. The radial bearing has four recesses and is mounted in the centre of the turntable. The supporting system has twenty-nine pads (include 24 supporting pads, 1 preload pad and 4 radial pads) which need a constant flow of oil, but it is expensive Fig. 1. The structure of hydrostatic rotary table Fig. 2. The structure of countertop and unnecessary to assemble twenty-nine constant displacements. The hydrostatic oil supply is divided into many portions by multi-point gear oil separators, which are always connected by a shaft. Therefore, only a small number of oil pumps can supply all the oil pads. For multi-pad hydrostatic bearings, compensation is needed, because of the issues of mutual influence between the oil pads. Restrictors (e.g. orifice, capillary, cylindrical-spool valves and so on) are widely used as compensation devices as described in detail in references [29] to [31]. Constant flow is another effective compensation method via the use of multi-point gear oil separators and so every pad has the same flow rate, and there is little interaction between the oil pads. All the equal gear pumps are connected by one shaft, which forces them to rotate at the same speed and produce the same flow rate. A detailed compensation theory of constant flow is described on page 17 of reference [30]. 2.2 Characteristic Modelling of Structural Discipline The rotary table countertop structure diagram is shown in Fig. 2, and it is mainly composed of a rib structure. The size parameters directly determine the performance of the countertop. Ideally, the rotary table should have a minimum mass and deformation while having the best static and dynamic performance. The decrease of mass, structure deformation and the increase of natural frequency can not only reduce costs but also enhance the accuracy of the rotary table. Therefore, the optimization goals for the structure discipline are to minimize the mass and deformation and to maximize the natural frequency. The value of countertop's mass, deformation and natural frequency are determined by all the design parameters. Then, the mass function, deformation function, and modal frequency function can be written as follows: freq = f ( ax, a2, a3, a4, a5, a6, a7, RL, Rs, Dy, Dr, b), (1) mass = fš(a1,a2,a3,a4,a5,a6,a7,RL,Rs,Dy,Dr,b), (2) deform = f Ц,a2,a3,a4,a5,a6,a7,RL,Rs,Dy,Dr,b). (3) RT is a design parameter of the countertop and has some influence on its mass, deformation, and natural frequency. In Eqs. (1) to (3), RT is not considered not because of its lowest influence but that it is mainly decided by the machining capacity of the NC machine tool. In other words, RT is decided by the designer of machine tool not by the designer of the rotary table. Therefore, when the design of the rotary table is started toRT has already been determined as a constant value. The objective function in structure part can be written as: f = max[ freq, -mass, -deform], (4) where the minus sign in Eq. (4) means the optimization goals is minimum. 2.3 Characteristic Modelling of Hydrodynamics Discipline 2.3.1 Establishment of Reynolds Equations In this study, it is assumed that thin film lubrication theory is applicable, and the flow in bearing is isothermal, laminar and axisymmetric. The Reynolds-type equation and the radial fluid flow equations can be given as (The specific derivations are shown in Appendix): 1 d ( rhh dp ) _dh r dr 12n dr dt Q(r )= - nrk3 dp 6n dr (5) (6) where r is the radius, h is the film thickness, p is the film pressure, t is the time, and n is the lubricant viscosity. 2.3.2 Calculation of Supporting Oil Pad As shown in Fig. 3, the supporting pad of the rotary table is a circular recess pad, R1 and R2 are the inner radius and outer radius of the pad respectively, Q0 is the flow rate supply to the pad. Fig. 3. The structure of the supporting oil pad By solving Eq. (5) with boundary conditions: p, ( r=R1 ) = p0i ; p, ( r=R2 ) = 0 , the film pressure profile p i (r) can be obtained. Then, by substituting p, (r) into Eq. (6), recess pressure p0i can be obtained. Finally, the load-carrying capacity of the bearing is calculated by integrating the film pressure. F = dh 3n[2Q0 -n(R22 + ](R22 -R2) dt 2h (7) where the detailed derivations of p0i , pi (r) and Fi are shown in Appendix. In addition, the stiffness, the damping coefficients and the pump power are obtained as follows: K Si = Csi = N.. = 9Qn(R - R2) h 3щ(R - -Ri4) 2 h 6n ln(R )Qo2 Ri (8) nh 2.3.3 The Calculation of the Preloaded Oil Pad wherep0y andpy (r) are the recess pressure and the film pressure profile of a preloaded oil pad, respectively. Then the stiffness, the damping coefficients and the pump power are obtained as follows: R R 9Qn((r2c4 -Я2з)Н^) + (R2C1 -R2c2)ln(-RCL)) K =-rRr-Rc^~, (10) ) y RR ЛС 2Л C 4 3цл ( К - ^2 + Кз - ) RC 2 RC 4 C = У -(RCl - RC2 + RC3 - Rf74 ) 2hi\n(Rc^ ) (13) RC 2 RC 4 NTy = 6^ln(^ )ln( ^ )Q2 nh in( IRc4Rc2 ) (12) ЛС 3Л C1 where the detailed derivations of p0y , py (r) and Fy are shown in Appendix. The preloaded oil pad is an annular recess pad as shown in Fig. 4, RC1, RC2, RC3, RC4 are the structural parameters of the annular recess pad. Fig. 4. The structure of preloaded oil pad The load-carrying capacity of the bearing can be calculated by integrating the hydrostatic film pressure: 77 3n F =---X y , R^R 2h3 Inf C 'C3 ) y RC2RC4 x[2ft (RC4 -RC3)ln(^ +(RC1 -RC2)ln(^ C 2 C 4 RR ^ (R4 r4 + r4 r4 )1п(дC1 C3) (RC1 - RC2 + RC3 - RC4) 1n(R R " C 2 C 4 (D2 D2 i D2 D2 \2 -(RC1 - RC2 + RC3 - RC4 dh (9) 3.3.4 Optimization Objective Function in Hydrostatics The support system consists of supporting oil pads, preloaded oil pad and radial bearing. So the stiffness, damping coefficients and the pump power can be calculated as follows: K7 = V K +У K f + K Z si s] y i=l nI ]=1 Kt = У K si ( Rl sincri ))2 + У Ks](Rs sin(^] ))2 i=i nl j =l Cz = y C +У C - C Z I si s] y . (13) i=l nl ]=l с, = У Csi(Rl sincri))2 + У Cs](Rs sinfa]))2 i=1 nl ]=1 N =У NTl +У N]+ NT i=l ]=l Forces on the turntable are balanced in the initial state, so we have Eq. (14). nI n 2 Z Fi + Z F j - Fy - G = 0. (14) i=1 j=1 Turntable manufacturers tend to determine the film thickness of the supporting oil pads and the preloaded oil pad at first and then determine the flow + rate of the oil pad, so in this work the flow rate of preload oil pad Q1 is calculated by Eq. (14). The objective function can be written as: f = max[ Kz, K t, C z , C,, - NT (15) Therefore, the total objective function in hydrostatics discipline can be written as: f = min[f1 , f2]. 3 MDO BASED ON SENSITIVITY ANALYSIS (16) 3.1 Sensitivity Analysis Different design parameters of the rotary table have different influences on its performance. Because there are many design parameters, it is difficult to realize optimization with a good convergence. Therefore, it is feasible to select the key parameters that have a significant influence on performance to implement optimization. The sensitivity is the gradient of a concern target to the design parameters of the rotary table, and SA can help identify the key parameters [43]. Here, SA can be divided into two parts: the first part is about the rotary table structure SA; the second part is about hydrostatic part SA. The derivations of objective function fb f2 to the parameters x1 and x2 can be found respectively. Where, Xi = [«i, a2, a3, a4, a5, a6, a7, RL, R., Dy, Dr, b] X2 = [ RS' RL' R1' R2' RC1' RC 2' RC 3' RC4,n' Q0 ' n1, К 0' hy 0]. So, the sensitive matrix can be written as: si[ / , /гУ- dfreq( xj ) ddeform( xj ) dKZ da. da. dR, dfreq( Xj ) ddeform( Xj ) dKZ db db Sh„, dNT dR, dNT dhy0 (17) As for the first part, a parameterized FEM model is established to calculate SIf1) in ANSYS software, and the sensitivity of the second part is calculated in MATLAB according to Eq. (13) and Eq. (17). The partial differential equations in Eq. (17) can be solved by numerical method. For example, dfreq(x1)/da 1 can be calculated by converting it into a differential equation \freq(a1h)-freq(a1l)] / (a1h - a1), in which a1h and a1{ are upper boundaries and lower boundaries of parameter a1 respectively, and all design parameters have no change except a1 in the calculation offreq(a1h) andfreq(a1l). In this way, Eq. (17) becomes Eq. (18): freq{ a h )-freq( a,, ) SI[ f , f2 ] = ddeform{ a1h )-ddeform(a1l ) KZ (RSh ) — KZ (RSl ) Да, ddeform{ bh )-ddeform{ bl ) Д Ms Kz ( ^h ) — Kz (wi ) Aw NT ( Rsh ) — Nz ( Rsl ) ARs Nt ( Wh ) — Nz ( wi ) Aw . (18) 3.2 MDO Strategy The MDO strategy starts from the design problem itself and is aimed at calculation structure and information organization. It is necessary to decompose the coupling relationship for multidisciplinary analysis and optimization. Information organization of MDO can make virtual design between different disciplines and different CAE analysis software possible and practical. With current computing resources, techniques of multidisciplinary optimization can be integrated effectively with multi-objective optimization algorithms to search for optimal designs; the detailed process is illustrated in Fig. 5. Isight is a generic software framework for the integration, automation, and optimization of design processes [50]. In this optimization, commercial software, such as CATIA, ANSYS, and MATLAB, are integrated with Isight so that they can input design parameters and output analysis results through a unified software platform The optimization implementation is completed in five steps. Step 1: start Isight software; Step 2: set up in Isight software including connect MATLAB, CATIA, ANSYS with Isight, drag optimization component to the task, and set up boundary conditions, variables, objectives and other optimization parameters; Step 3: Run Isight and optimization component; Step 4: Iterations are implemented until the maximum iteration steps are attained; Step 5: Output the results and the optimization is completed. Fig. 5. Optimization implementation framework In this study, particle swarm optimization (PSO) is selected as the optimization algorithm to optimize the turntable. PSO was proposed by Kennedy and Eberhart [51], Shi, and Eberhart [52] and Kennedy [53]; it is a stochastic optimization method based on swarm intelligence theory. It has been successfully applied in continuous optimization problems such as neural network training [54], voltage stability control [55], distribution route selection [56] and the optimization of cutting parameters [57]. PSO mimics the social behaviour of animal groups such as flocks of birds or fish shoals [58]. The process of finding an optimal design point is similar to the food foraging activity of animals. During the searching process, an animal can obtain maximum global optimization results via group co-operation. In PSO, particles represent potential solutions of the problem, every particle associated with two parameters: the position xid and velocity Vid in dimension d. When PSO is used to search for the best solution of a problem, each particle's movement is influenced by its local best known position, but is also guided toward the best known positions in the search-space, which are updated as better positions are found by other particles. This is expected to move the swarm toward the best solutions [51]. Specifically, position xi, d and velocity Vi: d are updated in (+ 1)th iteration by the equations as follows: Vid (t +1) = wVid (t) + +c lrl(pbestXìd (t ) - Xid (t )) + +c2 r( gbestXd (t ) - Xi,d (t )), (19) Fig. 6. The implementation procedure of PSO Xìtd (t +1) = Xìtd-i(t ) + Vid-i(t ). (20) Here i is particle's index and d represents the dh design parameters, in other words, if there are np particles and nd design parameters need to be optimized, V and X are matrices of np rows and nd columns. pbestXid is particle i's personal best experience; gbestXd is the group best experience found by all the particles so far; c1 and c2 are the acceleration coefficients; r1 and r2 are two random numbers that generated with the uniform distribution in the range of [0, 1]; and w is the inertia weight that is used to balance the global/local searches of particles. The steps of PSO adopted in this research are showed in Fig. 6. First, we need to set up initial values and boundary conditions for all design parameters, maximum iteration times, inertia weight w and the number of particles np. Then the maximum fly velocity Vmax should be set up. Thirdly, iterations are processed to obtain the best position of particles by four steps. Step 1: calculate velocity Vid and position Xid by Eq. (19) and Eq. (20); Step 2: use new Xi, d to calculate new fitness values, where fitness values are KZ, Kt, CZ, Ct, Nt and freq, mass, deform and they are calculated in MATLAB and ANSYS respectively; Step 3: compare fitness values with pbestX,d and gbestXd then update them according to Eq. (21). Step 4: Finally, determine whether the program meets the termination conditions; if yes, the iterations will stop; if not, return to Step 1 and continue iterating. if : adpi,d (t +1) > pbestXi,d (t) ^ pbestX i d (t +1) = adpi, d (t +1); if : adpid (t +1) < pbestX i d (t ) ^ pbestX t d (t +1) = pbestX i d (t ); , , (21) if :max(adp=b.^,d (t +1)) > gbestXd (t) ^ gbestXd(t +1) = max(adpt=1~ np ,d(t +1)); if :max(adpi=x„„p4(t +1)) < gbeslXd(t) ^ gbestX d (t +1) = gbestX d (t ). 3.3 Multidisciplinary Design Optimization The parameterized FEM of the rotary table countertops is established in CATIA to analyse its sensitivity. The sensitivity of hydrostatic part will be computed according to Eq. (13) and Eq. (16). The normalized results are shown in Fig. 7. According to Fig. 7a to c, among all the design parameters of the countertop, a1, a3, a5, a7 have major influence on both mass and natural frequency, which means designer should take them seriously and a) b) Fig. 7. d) Analysis results; a) sensitivity of mass, b) sensitivity of natural frequency, c) sensitivity of deformation, and d) sensitivity of f2 c carefully and decide their values. As for the hydrostatic part, film thickness parameters (like: hs0, hy0) have the biggest influence on the performance, and in the next place are the flow rate (lQ0) and the structural parameters of pads (R1,R2,RC1,RC2,RC3,RC4). The influence of other parameters such as q, n1, w is very small and can be ignored. Therefore, Q0, R1, R2, RC1, RC2, RC3, RC4 can be selected as key parameters to do optimization analysis. The total objective function and constraints can be listed as follows: Objective: f= min [fi, f, ], subject to: 0.64 m X wj Ш i š p Aß щ2 •eq3 ▲ i ■ X di sform 1 w к X К ■ ■ *xx X > i i 0.75 0.85 0.95 1.05 normalized mass 0 ü 1.6 £ 41-4 -a aj JS 1.2 "л 1 . о с à ▲ 1 F * ф 4 & > i > XX v И ft ♦ X *Kz UKt L Cz P 1 X Ct 1 II ■ ■ —1 1— 0.85 b) 0.95 1.05 normalized NT 1.15 Fig. 8. The pareto optimization results; a) of structure part, and b) of hydrostatic part were performed on the turntable bearing after assembly. The first test was a sweep test, as shown in Fig. 10, done by exciting countertop of the rotary table about 20 minutes with a vibration exciter, and the sweep range was 5 Hz to 100 Hz; speed was 0.08 Hz/s. Finally, we obtained the first, second and third natural frequencies of the countertop. Furthermore, FEM models of the initial and optimized countertops were established in ANSYS to perform a modal analysis. Test and modal analysis results are shown in Table 2. The natural frequency of the countertop after optimization is 33.79 Hz, and the experimental value is 37.25 Hz. The error between them is not large; this means that the simulation results are accurate. Therefore, we can draw that the dynamic performance of the rotary table has been significantly improved. Table 2. Natural frequencies of rotary table First Second Third order order order Simulation values for initial design parameters [Hz] 31.04 55.38 64.82 Simulation values for optimized design parameters [Hz] 33.79 60.93 68.69 Experimental values [Hz] 37.25 64.09 75.58 Fig. 9. The optimized hydrostatic rotary table Fig. 10. Sweep test setup Secondly, the bearing capacity of the rotary table was tested, as shown in Fig. 11. The average film thickness of each oil pad was measured by a dial indicator when the rotary table was under different load states. The load varied from 0 t, 150 t, 280 t to 410 t. Fig. 12 shows the measurement results when the load is 410 t. According to prior experience, the film thickness should be greater than 0.08 mm when under the maximum load of 410 t. Therefore, it can be assumed that the designed rotary table can meet the requirement. The error between experimental and optimized simulation curve is less than 20 %, which is within a reasonable range of error. The optimized simulation curve is slightly better than the former film thickness curve, which indicates that the stiffness of the film has gained a certain improvement. Therefore, the above results indicate that this optimization method is effective and feasible. а 0.16 0.15 0.14 0.13 0.12 0.11 0.1 0.09 0.08 -♦-simulati oli be fore optimiz —Ш— nimiilnti atioii nn ifì РГ OJ itimi z Uli ail ation vi operin- leiital value 0 150 load 280 IO3 [kg] Fig. 12. Film thickness variation of rotary table 410 Fig. 11. Bearing capacity experiment setup 5 CONCLUSION A heavy duty hydrostatic rotary table is often used as a means for providing rotational motion and handles supporting and rotating the workpiece in heavy-duty five-axis machine tools. Its rigidity, precision and carrying capacity are directly related to the machining ability and the accuracy of the NC machine tool. Due to the lack of effective means of analysis and experiment, the traditionally designed rotary table cannot provide users with reliable working performance and bearing capacity, and the rail grinding risk will cause serious economic loss to users. Therefore, properly designing a hydrostatic rotary table of a heavy-duty five-axis machine tool is very important for machine manufacturers, which aims to reduce its mass, supply pump power and save energy on the premise of good performance. The traditional design of a hydrostatic rotary table is normally performed by different teams sequentially, with expertise in a specific discipline. This will cause contradictions and result in excessive iterations, higher costs or longer design process time. In this paper, a sensitivity-based multidisciplinary optimization design method of a hydrostatic rotary table is proposed. The method takes hydrostatics and the structure's disciplinary characteristics into consideration, and a comprehensive optimization model is established. To achieve the optimization goal, PSO is introduced in Isight with the integration of CATIA, ANSYS, and MATLAB. In order to optimize with good convergence, SA is adopted to identify the key design parameters that have significant influence on the performance of a rotary table. Characteristics of this method are summarized as follows. Compared with the current methods in manufacturing, the proposed multidisciplinary optimization method establishes an optimization model with the integration of hydrostatics and structure disciplinary characteristics together, which can reduce the iterative modification caused by sequentially design by experts in different disciplines. A hydrostatic rotary table is a complex product and has many design parameters, which increases the difficulty of optimization. SA is introduced to identify the key design parameters that significantly influence the performance of hydrostatic rotary tables; therefore, the optimization convergence is improved. Despite the progress, it is important to note the limitation of the method needed to be further addressed to perfect the current work. In this work, the thermal effect is not considered, and the temperature is still cannot be ignored in the design process. Therefore, the development of an effective approach to integrate the thermodynamics characteristic into the optimization process is another focus for future research. 6 ACKNOWLEDGEMENT The authors are most grateful to Beijing Nova Program (xxjh2015106), Jinghua Talents Programme of Beijing University of Technology, the Leading Talent Project of Guangdong Province, and Open Project of State Key Lab of Digital Manufacturing Equipment & Technology (Huazhong University of Science and Technology), Shantou Light Industry Equipment Research Institute of science and technology Correspondent Station (2013B090900008) for supporting the research presented in this paper. 7 REFERENCES [1] Suh, S.H., Lee, E.S., Jung, S.Y. (1998). Error modelling and measurement for the rotary table of five-axis machine tools. The International Journal of Advanced Manufacturing Technology, vol. 14, no. 9, p. 656-663, D0l:10.1007/ BF01192286. [2] Zhang, Y., Yang, J. Zhang, K. (2013). Geometric error measurement and compensation for the rotary table of five-axis machine tool with double ballbar. The International Journal of Advanced Manufacturing Technology, vol. 65, no. 1-4, p. 275-281, D0I:10.1007/s00170-012-4166-4. [3] Mekid, S., Ogedengbe, T. (2010). A review of machine tool accuracy enhancement through error compensation in serial and parallel kinematic machines. International Journal of Precision Technology, vol. 1, no. 3-4, p. 251-286, DOI:10.1504/IJPTECH.2010.031657. [4] Florussen, G.H.J., Spaan, H .A.M. (2012). Dynamic R-test for rotary tables on 5-axes machine tools. Procedia CIRP, vol. 1, p. 536-539, D0I:10.1016/j.procir.2012.04.095. [5] Jodei, J., Ebrahimi, M., Roshanian, J. (2009). Multidisciplinary design optimization of a small solid propellant launch vehicle using system sensitivity analysis. Structural and Multidisciplinary Optimization, vol. 38, no. 1, p. 93-100, D0I:10.1007/s00158-008-0260-5. [6] Schmit, L.A. (1981). Structural synthesis-its genesis and development. AIAA Journal, vol. 19, no. 10, p. 1249-1263, D0I:10.2514/3.7859. [7] Haftka, R.T. (1977). Optimization of flexible wing structures subject to strength and induced drag constraints. AIAA Journal, vol. 15, no. 8, p. 1101-1106, D0I:10.2514/3.7400. [8] Lu, S., Kim, H.M. (2010). A regularized inexact penalty decomposition algorithm for multidisciplinary design optimization problems with complementarity constraints. Journal of Mechanical Design, vol. 132, no. 4, p. 1-12, D0I:10.2514/3.7400. [9] Sobieszczanski-Sobieski, J., Haftka, R.T. (1997). Multidisciplinary aerospace design optimization: survey of recent developments. Structural Optimization, vol. 14, no. 1, p. 1-23, D0I:10.1007/BF01197554. [10] Agte, J., De Weck, O., Sobieszczanski-Sobieski, J., Arendsen, P., Morris, A., Spieck, M. (2010). MDO: assessment and direction for advancement-an opinion of one international group. Structural and Multidisciplinary Optimization, vol. 40, no. 1-6, p. 17-33, D0I:10.1007/s00158-009-0381-5. [11] Bi, Z., Wang, L., Wu, C., Yang, G., Zhang, D. (2013). Multidisciplinary Design Optimization in Engineering. Mathematical Problems in Engineering, vol. 2013, article ID 351097, D0I:10.1155/2013/351097. [12] Yao, W., Chen, X., Ouyang, Q., van Tooren, M. (2012). A surrogate based multistage-multilevel optimization procedure for multidisciplinary design optimization. Structural and Multidisciplinary Optimization, vol. 45, no. 4, p. 559-574, D0I:10.1007/s00158-011-0714-z. [13] Ning, A., Kroo, I. (2010). Multidisciplinary considerations in the design of wings and wing tip devices. Journal of Aircraft, vol. 47, no. 2, p. 534-543, D0I:10.2514/1.41833. [14] Alonso, J.J., Colonno, M.R. (2012). Multidisciplinary optimization with applications to sonic-boom minimization. Annual Review of Fluid Mechanics, vol. 44, p. 505-526, D0l:10.1146/annurev-fluid-120710-101133. [15] Balling, R., Rawlings, M.R. (2000). Collaborative optimization with disciplinary conceptual design. Structural and Multidisciplinary Optimization, vol. 20, no. 3, p. 232-241, DOI:10.1007/S001580050151. [16] Geyer, P. (2009). Component-oriented decomposition for multidisciplinary design optimization in building design. Advanced Engineering Informatics, vol. 23, no. 1, p. 12-31, D0l:10.1016/j.aei.2008.06.008. [17] McAllister, C.D., Simpson, T.W. (2003). Multidisciplinary robust design optimization of an internal combustion engine. Journal of Mechanical Design, vol. 125, no. 1, p. 124-130, D0l:10.1115/1.1543978. [18] Jixin, W., Mingyao, Y., Yonghai, Y. (2011). Global optimization of lateral performance for two-post ROPS based on the Kriging model and genetic algorithm. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 10, p. 760-767, D0l:10.5545/sv-jme.2010.246. [19] Peri, D., Campana, E.F. (2003). Multidisciplinary design optimization of a naval surface combatant. Journal of Ship Research, vol. 47, no. 1, p. 1-12. [20] Cai, G., Fang, J., Zheng, Y., Tong, X., Chen, J., Wang, J. (2010). Optimization of system parameters for liquid rocket engines with gas-generator cycles. Journal of Propulsion and Power, vol. 26, no. 1, p. 113-119, D0l:10.2514/1.40649. [21] Yifei, T., Wei, Y., Zhen, Y., Dongbo, L., Xiangdong, L. (2013). Research on multidisciplinary optimization design of bridge crane. Mathematical Problems in Engineering, vol. 2013, article ID 763545, D0l:10.1155/2013/763545. [22] Kroo, I.M. (1997). MDO for Large-Scale Design. In Alexandrov, N.M., Hussaini, M.Y. (eds.) Multidisciplinary Design Optimization: State-of-the-Art, SIAM, Philadelphia, p. 22-44. [23] Sobieszczanski-Sobieski, J., Haftka, R.T. (1997). Multidisciplinary aerospace design optimization: survey of recent developments. Structural Optimization, vol. 14, no. 1, p. 1-23, D0l:10.1007/BF01197554. [24] Martins, J.R., Lambe, A.B. (2013). Multidisciplinary design optimization: a survey of architectures. AIAA Journal, vol. 51, no. 9, p. 2049-2075, D0l:10.2514/1.J051895. [25] Zhang, X., Huang, H.Z. (2010). Sequential optimization and reliability assessment for multidisciplinary design optimization under aleatory and epistemic uncertainties. Structural and Multidisciplinary Optimization, vol. 40, no. 1-6, p. 165-175, D0l:10.1007/s00158-008-0348-y. [26] Du, X., Guo, J., Beeram, H. (2008). Sequential optimization and reliability assessment for multidisciplinary systems design. Structural and Multidisciplinary Optimization, vol. 35, no. 2, p. 117-130, D0l:10.1007/s00158-007-0121-7. [27] Yao, W., Chen, X., Luo, W., van Tooren, M., Guo, J. (2011). Review of uncertainty-based multidisciplinary design optimization methods for aerospace vehicles. Progress in Aerospace Sciences, vol. 47, no. 6, p. 450-479, D0l:10.1016/j. paerosci.2011.05.001. [28] Zeeshan Q, Yunfeng D, Rafique AF, et al. (2010). Multidisciplinary robust design and optimization of multistage boost phase interceptor. 51st AIAA/.ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Orlando, D0l:10.2514/6.2010-2920. [29] Rowe, W.B. (2013). Hydrostatic and Hybrid Bearing Design. Elsevier, London. [30] Bassani, R., Piccigallo, B. (1992). Hydrostatic Lubrication. TribologySeries, vol. 22, Elsevier, Amsterdam. [31] Hamrock, B.J., Schmid, S.R., Jacobson, B.O. (2004). Fundamentals of Fluid Film Lubrication. CRC Press, Boca Raton, D0l:1201/9780203021187. [32] Ting, L.L., Mayer, J.E. (1971). The effects of temperature and inertia on hydrostatic thrust bearing performance. Journal of Tribology, vol. 93, no. 2, p. 307-312, D0l:10.1115/1.3451572. [33] Deb, K., Goyal, M. (1997). Optimizing Engineering Designs Using a Combined Genetic Search. Proceedings of the 6th International Conference on Genetic Algorithms, p. 521-528. [34] Lin, J.R. (2000). Surface roughness effect on the dynamic stiffness and damping characteristics of compensated hydrostatic thrust bearings. International Journal of Machine Tools and Manufacture, vol. 40, no. 11, p. 1671-1689, D0l:10.1016/S0890-6955(00)00012-2. [35] Bakker, O.J., Van Ostayen, R.A.J. (2010). Recess depth optimization for rotating, annular, and circular recess hydrostatic thrust bearings. Journal of Tribology, vol. 132, no. 1, p. 11-13, D0l:10.1115/1.4000545. [36] Zhao Ming, Huang Zhengdong, Li Bing, Chen Liping. (2007). Hydrostatic pressure calculation and optimization for design of beam and slide-rest guideways in heavy duty CNC vertical turning mill. Chinese Journal of Mechanical Engineering (English Edition), vol. 20, no. 5, p. 16-22, D0I:10.3901/ CJME.2007.05.016. [37] Solmaz, E., Öztürk, F. (2006). Optimisation of hydrostatic journal bearings with parameter variations based on thermodynamic effects. Industrial Lubrication and Tribology, vol. 58, no. 2, p. 118-122, D0l:10.1108/00368790610651530. [38] Karkee, M., Steward, B.L. (2010). Local and global sensitivity analysis of a tractor and single axle grain cart dynamic system model. Biosystems Engineering, vol. 106, no. 4, p. 352-366, D0l:10.1016/j.biosystemseng.2010.04.006. [39] Bo, L., Kyu, P.N. (2013). Sensitivity analysis for identifying the critical productivity factors of container terminals. Strojniški vestnik - Journal of Mechanical Engineering, vol. 59, no. 9, p. 536-546, DOI:10.5545/sv-jme.2012.931. [40] Xu, C., Gertner, G. (2007). Extending a global sensitivity analysis technique to models with correlated parameters. Computational Statistics & Data Analysis, vol. 51, no. 12, p. 5579-5590, D0l:10.1016/j.csda.2007.04.003. [41] Rodriguez-Fernandez, M., Banga, J.R. (2009). Global sensitivity analysis of a biochemical pathway model. 2nd International Workshop on Practical Applications of Computational Biology and Bioinformatics, Advances in Soft Computing, vol. 49, p. 233-242, Springer, Berlin, Heidelberg, D0l:10.1007/978-3-540-85861-4_28. [42] Saltelli, A., Ratto, M., Tarantola, S., Campolongo, F., Commission, E. (2006). Sensitivity analysis practices: Strategies for model-based inference. Reliability Engineering & System Safety, vol. 91, no. 10-11, p. 1109-1125, D0l:10.1016/j.ress.2005.11.014. [43] Helton, J.C., Davis, F.J. (2003). Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliability Engineering & System Safety, vol. 81, no. 1, p. 23-69, DOI:10.1016/S0951-8320(03)00058-9. [44] Saltelli, A., Chan, K., Scott, E.M. (eds.) (2000). Sensitivity Analysis, vol. 1, Wiley, New York. [45] Sudret, B. (2008). Global sensitivity analysis using polynomial chaos expansions. Reliability Engineering & System Safety, vol. 93, no. 7, p. 964-979, DOI:10.1016/j.ress.2007.04.002. [46] Maly, T., Petzold, L.R. (1996). Numerical methods and software for sensitivity analysis of differential-algebraic systems. Applied Numerical Mathematics, vol. 20, no. 1-2, p. 57-79, DOI:10.1016/0168-9274(95)00117-4. [47] Leis, J.R., Kramer, M.A. (1988). Algorithm 658: ODESSA-an ordinary differential equation solver with explicit simultaneous sensitivity analysis. ACM Transactions on Mathematical Software, vol. 14, no. 1, p. 61-67, DOI:10.1145/42288.214371. [48] Saltelli, A., Tarantola, S., Chan, K.P.-S. (1999). A quantitative model-independent method for global sensitivity analysis of model output. Technometrics, vol. 41, no. 1, p. 39-56, DOI:10.1080/00401706.1999.10485594. [49] Morris, M.D. (1991). Factorial sampling plans for preliminary computational experiments. Technometrics, vol. 33, no. 2, p. 161-174, DOI:10.1080/00401706.1991.10484804. [50] Koch, P.N., Evans, J.P., Powell, D. (2002). Interdigitation for effective design space exploration using iSIGHT. Structural and Multidisciplinary Optimization, vol. 23, no. 2, p. 111-126, DOI:10.1007/s00158-002-0171-9. [51] Kennedy, J., Eberhart, R. (1995). Particle Swarm Optimization. Proceedings of IEEE International Conference on Neural Networks, Vol. 4, p. 1942-1948, DOI:10.1109/ ICNN.1995.488968. [52] Shi, Y., Eberhart, R.C. (1998). A modified particle swarm optimizer. Proceedings of IEEE International Conference on Evolutionary Computation, p. 69-73, DOI:10.1109/ icec.1998.699146. [53] Kennedy, J. (1997). The particle swarm: social adaptation of knowledge. Proceedings of IEEE International Conference on Evolutionary Computation, p. 303-308, DOI:10.1109/ ICEC.1997.592326. [54] Van den Bergh, F., Engelbrecht, A.P. (2000). Cooperative learning in neural networks using particle swarm optimizers. South African Computer Journal, vol. 26, p. 8490. [55] Fukuyama, Y., Nakanishi, Y. (1999). A particle swarm optimization for reactive power and voltage control considering voltage stability. Proceedings of IEEE international conference on intelligent system applications to power systems, Rio de Janeiro. [56] Wu, Z. (2014). Optimization of distribution route selection based on particle swarm algorithm. International Journal of Simulation Modelling, vol. 13, no. 2, p. 230-242, DOI:10.2507/ IJSIMM13(2)CO9. [57] Župerl, U., Čuš, F., Gečevska, V. (2007). Optimization of the characteristic parameters in milling using the pso evaluation technique. Strojniški vestnik - Journal of Mechanical Engineering, vol. 53, no. 6, p. 354-368. [58] Yeh, W.C., Lin, Y.C., Chung, Y.Y., Chih, M. (2010). A particle swarm optimization approach based on Monte Carlo simulation for solving the complex network reliability problem. IEEE Transactions on Reliability, vol. 59, no. 1, p. 212-221, DOI:10.1109/TR.2009.2035796. 8 NOMENCLATURE CSl [Ns/m] damping coefficients of ith supporting pads Ct [Nrad/s] incline damping coefficients of the turntable Cy [Ns/m] damping coefficients of preloaded pads Cz [Ns/m] axial damping coefficients of the turntable f [-] objective function of all the part fl [-] objective function for countertop f2 [-] objective function for hydrostatic part F [N] load-carrying capacity of ith supporting pads Fy [N] load-carrying capacity of preloaded pads G [N] weight of the countertop h. [mm] film thickness of ith supporting pads hy [mm] film thickness of preloaded pads Ks, [N/m] stiffness of ith supporting pads Kt [Nm/rad] incline stiffness of the turntable Ky [N/m] stiffness of preloaded pads Kz [N/m] axial stiffness of the turntable NT [W] overall pump power NTl [W] pump power of i'th supporting pads NTy [W] pump power of preloaded pads P0i [Pa] recess pressure of i'th supporting pads P0y [Pa] recess pressure of preloaded pads Pi(r) [Pa] film pressure of i'th supporting pads Ply(r) [Pa] film pressure of inner land in preloaded pads P2y(r) [Pa] ilm pressure of outer land in preloaded pads Qo [m3/s] flow rate which supply to every supporting pads Qi [m3/s] flow rate which supply to preloaded pads 9 APPENDIX In this study, it is assumed that thin film lubrication theory is applicable and the flow in bearing is isothermal, laminar and axisymmetric. So the N-S equations can be simplified as: 1 d(rur ) + d(vj = 0, r dr dp dz = d\ dr dz2 d-P = 0. & (A1) (A2) (A3) Film pressure has no change in z direction according to Eq. (A3) above, so integrating both sides of Eq. (A2) with boundary conditions: z=h, ur = 0; z = 0, ur = 0; at the surface of the pad, then the radial velocity is found to be: z ( z - h) dp 2ц dr (A4) Substituting Eq. (A4) into Eq. (A1) and integrating it with boundary conditions: z=h, ur = 0; vz=dh / dt; z = 0, ur = 0, vz = 0, then the Reynolds equation is obtained as: 1 d ( rh dp ) _dh r dr 12n dr dt (A5) Integrating Eq. (A4) at z e (0, h) and ре (0,2n) the flow rate is: Q(r) = J J urrdzdq> - nrhh dp вц dr (A6) where r is radius, h is film thickness, p is film pressure, t is time, and n is lubricant viscosity. For supporting pad (shown in Fig. 3) h ^ h (i = 1,2,3,~ n); p ^ pt. Integrating both sides of Eq. (A5) two times for r, then we have 12n r2 dh Pi (r) = -f(Aln(r) + B + — -L), hi 4 dt (A7) where A and B are unknown constants. The boundary conditions for circular recess pads can be provided in Eq. (A8): (when r = Rj, pi = p0i [when r = R2,pi = 0 Substituting the boundary conditions into Eq. (A7), we can gotten that: (A8) 1 ln( •A) 12n hi R - R2 dh ^ —4 If) (A9) R) 1 (h ln(R,) - R2 ln(R,) - R2ln(R) dhL) ( 12n Pm - 4 di' Substituting Eq. (A9) into Eq. (A7), then the pressure distribution can obtained: r ln( - ) Pi (r ) = Po i-— ■ ■n< f) 3n(—2ln(-) - ->(—)) — — ) Щ ln(-) ' K, (A10) Assuming fluid incompressible, the lubricant flowing out from ith pad is can be divided in two parts; one is Q0 flow from the pump, another is caused by squeeze velocity dh and the flow rate of this part is: -nR —. So the flow dt dh continuity equation is Q(R2 ) = Q0 - nR —- . Substituting dt Eq. (A9) and Eq. (A10) into the flow continuity equation, then recess pressure can be calculated: p = 6ПкыRl) - nR - R2) dhL 0 nh) ( R hi dt (A11) Then, the load-carrying capacity of the bearing can be calculated by integrating the mean steady hydrostatic film pressure: 2 FF = kRp0i + 2n j rpt (r)dr - 3nm-n( R + R1) f ](R,2 - R?) _dt_ 2h3 ' (A12) In addition, the stiffness, the damping coefficients and the pump power are obtained as follows: K F = 9Qm(R22 -R2) Ksi dh 'dhu h CSi=- dF]_ 3 R24 - R,4) 5( ^ ) dt M 2h 6nln3Qo2 (A13) Ri R, Nn= Poi (-£ = 0)Qo =- .3 dt л h For annular recess pad (shown in Fig. 4), its boundary conditions are: 0 0 A 2 B p = py ; h = hy when r = RC1, py = 0; when r = RC2, py = p0y when r = RC3, py = p0y; when r = RC4, py = 0 . (A14) Qi -n(R2C4 -R2C1)^ = -Q(Rcl) + Q(Rc4) dt Substituting Eq. (A6), (A15) and (A17) into flow continuity equation: dt Qi -n(R2C4 -R21)-y = -Q(RC1) + Q(RC4) then recess pressure can be calculated: When r g (RC3, RC4), the pressure distribution is same as Eq. (A10). By replacing R1, R2, p0i and h, with RC3 , RC4 , p0y and hy , p2y(r) can be written as: Pl у(r) = Po r ln( — ) RC 4 _ ln( Rc3 ) RC 4 3n Po, 2ßi ln( R4) + RC2 RC3 +П1Ш( r^ - о+- io nh in( Rc,Rc 3 ) yRR . (A18) The load-carrying capacity of the bearing can be calculated by integrating the hydrostatic film pressure. 3nr +Hv + 2 3n(R.4ln(—) - ЛСзН—)) R, Кы^3) R dhy (A15) When r g (RCl, RC2), substituting the boundary conditions py(RC1)=0, py(RC2) = p0y into Eq. (A7), we can get that: A = _1_( hy p - RC2 RCI 'У ) ln(R2) -ln(RCI)(12n y ^ ) ri, - Rci h 4 dt B = - 1 ( hyin(Rci) ln(RCI) - ln( RC 2) 12n RC,ln(Rci) -RCiln(Rc,) dhy Po У (A16) 4 dt Substituting Eq. (A16) into Eq. (A7) ,we have: Fy =n(RlC3 - RlCl)P0y + 2П J rPly(r)dr + 2П J rPly(r)dr 3n 2h3ln( RciRc3 ) У J? J? RC 2 RC [2Qi R (Rc4 - RC3)ln^-C1) + RC2 +(R2I - R?2)ln(^) V RC +П с 4 RR ^ (RC1 - RC2 + RC3 - RC4) ln( „Cl 'C3 ) -(RCl - RC2 + RC3 - RC4 )2 dhy_ dt . (A19) Then the stiffness, the damping coefficients and the pump power are obtained as follows: K,. R R 9Qn((Rc4 -RC3)in(RCL) + (Rži -R^HR3)) Rc 2 Rc h4]n(RiRc3 ) y DD , (A20) r ln( — ) Pl у (r ) = Po у— ln( ) R„ r 2 3n(RCiM—) - R^M—)) 3r]r R h-3 ln( R2 ) ^) ^. (A17) dt Ъцл f RR \ (RC1 - RC2 + RC3 - C1 nC3 ) \-(RCl - RC2 + R<23 - RC4)2 2h3 ln( RciRc3 ) y RR Лс ^C 4 (A21) = Po 6^ln( )ln( )ßf nhl ln(Rc4Rc2) y RR Лс ^Cl (A22) R у + у Strojniški vestnik - Journal of Mechanical Engineering 61(2015)7-8, 448-458 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2014.2372 Original Scientific Paper INS/GPS Navigation System Based on MEMS Technologies Vlada Sokolović1* - Goran Dikić1 - Goran Marković2 - Rade Stančić1 - Nebojsa Lukić3 1 University of Defence in Belgrade, Military Academy, Serbia 2 University of Belgrade, School of Electrical Engineering, Serbia 3 Serbian Armed Forces, Military Test Centre, Serbia Integrated navigation systems are widely used in all aspects of life, i.e. in various civilian and military applications. In the previous two decades, most of the attention in this area has been directed toward the use of low-cost sensors, which are characterized by the low accuracy. Therefore, many algorithms have been developed to implement inertial sensor errors damping, which are based on the integration of inertial sensors with the external sources of information (GPS receivers, magnetometers, barometers, etc.) or mathematical models developed for the self-damping of inertial sensors errors. This paper proposes a specific method of gyroscope drift compensation by using a PI controller based on the magnetometer measurements, as well as a method of errors compensation in the horizontal channel of navigation system by using the adaptive control signals. The experimental results obtained for the proposed solutions presented in this paper confirmed the possibility of the successful application of these solutions in the real-world environment. Keywords: barometer, global positioning system, inertial navigation, integration, Kalman filter, magnetometer, MEMS Highlights • Usage of magnetometer and accelerometers for gyroscopes drift estimation. • Specific method of gyroscope drift compensation by using PI controller. • Error compensation in the horizontal channel of navigation system by using adaptive control signals. • Analysis of experimental results and validation of the proposed solutions. 0 INTRODUCTION In accordance with current global trends, specifically the demand to achieve the highest possible accuracy of navigation devices while simultaneously minimizing their costs, low-cost sensors, such as sensors made in micro-electro-mechanical systems (MEMS) technology, have increasing practical application. The largest contribution to the total error of inertial navigation systems (INS) is from the errors of inertial sensors [1]. In comparison to high accuracy sensors, the errors of low-cost sensors have a nonlinear nature that is very difficult to model and to alleviate. Consequently, there is a large error in determining the position and the attitude of the object when these low-cost sensors are used [1]. Bearing in mind that the computation of the new velocity and position values by the navigation algorithm are obtained based on the previous computations, it is obvious that these errors are cumulative and that their values increase rapidly in time, especially for a long-term navigation process. Therefore, it is necessary to perform a periodic correction of the position and velocity values by using an additional navigation system or other independently acquired external measurements [1]. By increasing the accuracy of calculated navigation solution developed by using the low-cost sensors, the construction costs of the INS equipment would decrease. As a result, the relatively weak demand for their usage, which is a key problem in current research in the field of navigation systems based on MEMS technology, should increase. The Global Positioning System (GPS) is most commonly used as an external source of information for the correction of INS navigation parameters [1]. However, GPS accuracy depends on the spatial arrangement of the GPS satellites relative to the observed receiver, the occurrence of multipath signal propagation, the specific influence of the ionosphere on the GPS signal propagation and the presence of noise in the GPS receiver [2] to [4]. In order to improve the accuracy of the integrated navigation solution, especially for the periods when the GPS receiver does not provide an output with the navigation information, other sensors, such as magnetometers or barometric altimeters (composed of a temperature and a pressure sensor), can be used in addition to the integrated INS/GPS solution. The previously mentioned sensors are also characterized by the specific errors that can be modelled and damped. In this paper, the results of the research conducted with the goal to develop a model of an integrated INS/GPS/Barometer/Magnetometer navigation system based on a loosely-coupled (LC) integration method is presented. The validity of this research 448 *Corr. Author's Address: University of Defense in Belgrade, Pavla Jurišića Šturma 33, Serbia, vlada.sokolovic@va.mod.gov.rs lies in a continuous development of methods and algorithms for sensors integration, which are applied in accordance with the advances in low-cost sensor technology, as the one of the main reasons for further research in this field. The aim of this research was to develop a model that will provide the unified processing of all collected measurements acquired by using the INS, GPS, barometer and magnetometer, in order to obtain more accurate information on the object's spatial position, velocity and attitude. The practical goal of this study was to enable the development of the real-world navigation system and to create the suitable conditions for the further research in the field of multi-sensor integrated navigation. The applied methods of modelling, which include the dominant application of mathematical models, are used to determine the errors models for individual sensors and to enable successful implementation of the observed navigation system model in the MATLAB environment. The performed analysis has shown the positive contributions of magnetometer usage in order to determine the attitude of the object, as well as of baro-altimeter usage in order to determine the vertical component of velocity and altitude. Throughout this study, random processes are described as first-order Gauss-Markov or the additive white noise processes. Furthermore, the errors that occur due to quantization, averaging, rounding of the measured values and conversion between different data types are ignored. The accuracy of the navigation solution is established based on the control points (CP) defined for the test trajectory. These referent CPs are precisely determined by the differential GPS. A verification of the proposed solutions is performed through the experiment on the land vehicle. Therefore, the proposed solutions are tested based on the data collected during this experiment and with the adequate software implemented in MATLAB. In this work, several results represent the original contribution in the field of research, of which the most important are the following: • By using the proportional-integral (PI) controller, the algorithm of the magnetometer and gyroscope triad integration has improved the determination of attitude of the object; • A method for error damping of the horizontal channel of integrated navigation system, based on the adaptive error damping coefficients, is defined; • The verified results certainly represent a notable improvement in the given field and, therefore, proposed solutions are suitable for practical application in the real-world navigation systems. In the second section, the attitude determination with gyroscope drift compensation based on magnetometer and accelerometer triad measurements by utilizing the PI controller is described. In the third section, the error compensation in the horizontal and vertical channel of a "strap-down" INS (SINS) is presented. The general scheme of the proposed INS/ GPS/BAR/MAG integration and basic theory of the Kalman filter (KF) with an adaptive control signal, which is used for the integration, is detailed in the fourth section. The fifth section contains a description and the results of the experimental testing of MEMS sensors as well as of the experimental tests that have been realized by using a land vehicle. Due to the nonlinear nature of the dynamic systems in practice, there is no single solution to the navigation algorithm. Thus, in the present research, several methods of the estimation of states have been applied in the integrated navigation systems. Bian et al. proposed an adaptive navigation solution, in which a covariance matrix and a matrix of Kalman gains are adapted by the maximum probability function [5]. Stancic and Graovac analysed an INS/ GPS navigation system, in which the integration is realized with the extended KF (EKF) with the control signal [6]. Noureldin et al. [7] used neural network instead of KF for the INS/GPS integration. The damping of errors in the vertical channel of INS during the outage of GPS information is usually realized with the baro-altimeter, as described by Readdy and Saraswat [8]. In this work, the attenuation of vertical channel errors is performed by the vertical channel error-damping loop with the external source of altitude from GPS or baro-altimeter measurements. Sokolovic et al. [9], used the adaptive EKF control signal for errors damping in the vertical channel of the integrated navigation system. The magnetometer is the most commonly used device for object attitude determination. For attitude error damping, Madgwick et al. [10] used the negative gradient between the axis of the gyroscopes and the magnetometer. Sokolovic et al. used gain coefficients for gyros drift-damping error compensation [11]. Moreover, the magnetometer was used for the initial attitude determination and for the INS heading angle correction, e.g. in Hao et al. [12], Zhao and Wang [13] and Han and Wang [14]. 2 ATTITUDE DETERMINATION In this section, the main idea regarding obtaining information about the attitude from two sources is described. The first source is gyroscopes and the second source is a magnetometer and accelerometers as external sensors. Since gyroscope measurements are subjected to certain known errors, the external sensors are used in order to estimate the gyroscope measurement error by utilizing PI controller approach. If m, denotes the components of angular rotation velocity of the body coordinate system relative to the inertial space, then the quaternion innovation at the time к + 1 is given as [15], 4 +1 2 qk (aT = qkAq' (1) where Г is a sampling time, and Aq is updating quaternion. In order to define the horizontal and vertical reference plane by measurements of accelerometers and magnetometers, an objective function is applied here [10] and [16], f (q N, d n ) = q N (2) where dN is some referent vector described in frame N, qN is quaternion that describes the orientation of the frame N relative to frame B and q^* is the conjugate of qN. In our case, N is the navigation frame and B is the body (sensor) frame. We assumed that the sensor frame and body frame coincided. A complete form of the objective function can be written as: 2dx(2 - q - ?42) + 2dy (?1?4 + ?2?3) + 2dz (?2?4 - ?1?з) f faN > d N ) = 2dx (?2?3 - ?1?4) + 2dy (1 - 42 - q2) + 2dz (?1?2 - ?3?4) 2dx(?1?3 + ?2?4) + 2dy (?3?4 - ?1?2) + 2dz (-2 - 4/2 - lì) (3) where dx, dy, and dz are elements of the vector dN, and the indexes (x, y, z) denote the axis of the navigation coordinate system. The reference vector of gravity coincides with the vertical axis of the navigation frame, so the reference vector dN is defined by using the normalized value of the vector of gravity and is given as: d s = [0 0 0 1]. (4) Estimated reference directions of gravity, based on Eq. (3) and Eq. (4), can be written as: v = - qq 2(qq2 -q3q4) 2 2 2 2 q - q2- qs + q4 (5) The reference vector of the magnetic field is calculated based on Eq. (6), where m is the normalized vector of geomagnetic field measured by magnetometer, h w = 2m °.5 - q2 - ql qq - qq q2 q2 + qq ЧгЧз+qq 0.5 - q22- ql q3q2- qq q2 q2- q qз q3 q2 + q qi 0.5 - q22 - q23 . (6) The geomagnetic field can be considered to have components in one horizontal axis and the vertical axis as; b = h and the reference vector dN can be defined as, d = [0 bx 0 bz (7) (8) Based on Eqs. (3) and (8), the estimated reference directions of magnetic flux are given as: w = 2bx(°.5 -q2 - q2)+2bz(qq -ш) 2bx (q q- qqO + 2bz (Mi + Ъ1л) 2bx (Мз + ?4) + 2bz (0.5 - ^ - ?32) (9) The gyroscope measurement error or the deviation of the measurement vector from the reference vectors is calculated as, e = (vx)a + (wx)m, (10) where (vx) and (wx) denote the skew-symmetric matrices of gravity, and magnetic flux estimated directions and a is the normalized vector of acceleration. The gyroscope drift and noise corrections are performed based on the estimated value of the error in the current iteration: о = ° + Kp • e + Ki • e^ (11) where ю denote the measurements of gyroscopes, ra c denote the estimated (corrected) values of gyroscopes measurements, Kp is the gyro drift damping coefficient, K is the gyro measurement noise damping coefficient, and eint is the integral of error in current iteration. The coefficients Kp = 0.01, and K = 0.005 represent the coefficients of the PI controller and are determined based on the Gauss-Smith method. The PI controller is chosen in order to allow the removal of coarse gyroscope drift and to enable the faster response to the sudden appearance of measurements error independently on the previous measurements. Corrected values of gyroscopes measurements are used to determine the updating quaternion and are based on Eq. (1) to calculate quaternion in the current iteration. Based on these values, the final quaternion, transformation matrix Cp, is formed. Index P represents the platform frame. The platform frame is an image of the navigation frame, which is established by an on-board computer. For the ideal sensors, the platform frame coincides with navigation frame N. In this paper, a north-east-down (NED) frame is used as a navigation frame. The fine filtering is realized by the EKF, which performs the estimation of the attitude errors that are used to form the matrix CNP . Thus, the final solution of the attitude is determined by the expression: CB = CP CB (12) where CNB is the resultant matrix, whose elements are used to determine the attitude of the object body relative to the navigation coordinate system, i.e. roll, pitch and yaw angles, marked as ф, в, щ. The functional diagram of the strap-down INS, with the coordinate systems that are used in this work, is shown in Fig. 1. Fig. 1. The functional diagram of SINS In Fig. 1, denotes gyroscopes measurements in the body frame and fB denotes specific forces, measured by the accelerometers, in the body frame. 3 INS ERRORS DAMPING The proposed solution, specifically the error compensation for the horizontal and vertical channels of INS, is achieved with the introduction of external information, e.g. the additional measurements obtained with the other sensors. In this section, the algorithm for errors damping for the horizontal and vertical channels of INS is presented. The single channel error diagram for INS with error damping by introducing external information from GPS is shown in Fig. 2. In this figure Kj and K2 are coefficients, Фе is attitude error in east direction (difference between gyro platform and navigation frames), Фе(0) is initial misalignment error between platform and navigation frames, SVn is velocity error, Re is the radius of the Earth, SRn is position error, SaBn are accelerometers "bias" and SgBe is gyro drift. Fig. 2. The North Channel error diagram of INS, with error damping by the GPS measurements usage Based on the diagram in Fig. 2, the single channel error model is defined as: SVn = gфe - K1SVn +SaBn SV Фе =- K2SVn + e Re 2 n (13) By obtaining the derivatives of the first equation in Eq. (13) (derivatives of SVn ), and introducing the second equation of Eq. (13) into the first equation of Eq. (13), we obtain: 8Vn + K1Sfrn + (m* + K2 g )SVn = - g • SgBe + SaBn. (14) The left side of Eq. (14) describes a second-order oscillator. The selection of optimal coefficients K1 and K2 is based on a compromise between the value of the static error and the system bandwidth. A selection of high values of K1 and K2 provides a low static error value, but in this case, the system is characterized by a wide bandwidth in order to allow transmission of high-frequency noise components [14]. Therefore, the selection of optimal coefficients, K1 and K2, should be performed as a compromise between these two opposite requirements. As the MEMS sensors inherently produce high static errors, the values of K1 and K2, can be very high. Such behaviour may cause the instability of the EKF; therefore, we have to separate these coefficients. The first pair of coefficients is marked as K1INS and K2IN5 and are used for attenuation of INS errors, while the second pair of coefficients, K1KF and K2KF, are used in the EKF in order to form control signals, as given in: K = -KinKF sinh(5 Vn ), < = -KieKF sinh(5 Ve ), < = -K2nKF Sinh(S Ve ), ul = -K2nKF Smh(SVn ), (15) where Kw = Kw = 0.01, and ^kf = ^2eKF = 0.007. The control signals at the horizontal channel, which are shown in Fig. 3, are defined by using the maximum expected ^velocity difference 8Vn (in our case 1.1 m/s), i.e. S Ve for the eastern channel. During the period of GPS data outage, the control signals are formed on the basis of the moving average values computed for the last 16 measurements of velocity, Vn and Ve, in order to obtain average values Vnma and Vema. Fig. 3. The control signals at the horizontal channel of the proposed INS solution In existing solutions [6] and [11], these values were used for the auto-correction of the horizontal channel errors, while in the solution proposed here , these values are used in order to compute the velocity vector of movement in the horizontal plane VHma =-\lV2ma - Vema , which can be decomposed into the north and east components of velocity by using the yaw angle. This way, we can determine: INS SVn = VniNs - VHma C0S(V ) , and INS SVe = VeiNs - VHma Sm(l^) . In order to achieve more effective error damping in case of GPS data outage, the higher values of coefficients K are necessary. For this reason the adaptation of these coefficients by using the nonlinear hyperbolic functions sinh() is performed and the values of coefficients K are selected to reach the maximum expected value of SVne . The function sinh() is defined for a wide domain of input values, with the property that as the input signal error values rise the function response increases slowly for small inputs (almost with a linear rise), while it increases more quickly for the higher signal error values. The adaptation of the coefficients in the vertical channel is adopted as described in [9], and the model of the vertical channel is given as: Si = SVz - C sinh(Sh - ShBAR), SVz = SAz - C2 sinh(Sh-ShBAR )-Sa + 2a2sSh, (16) where SVz is the velocity error, 8h is the INS height error, 8Az is the vertical acceleration error, 8hBAR is the baro-altimeter error, Sa is the steady state error, ms is Shuler's frequency, and Cb C2 are control coefficients pair for vertical channel error. The high values of these coefficients can cause system instability and, therefore, the separation of these coefficients is introduced, with C1INS, C2INS being coefficients for the INS correction, and C1KF, C2KF, being coefficients for the control signals in EKF. The coefficient values that were used in this paper are C1INS = 0.04, C2INS = 0.8, C1KF = 0.035 and C2KF = 3.2x10-4. When there are no valid GPS measurement data, the EKF works in the prediction mode, based on the covariance matrices Q (system noise) and R (measurement noise). 4 INTEGRATED NAVIGATION SYSTEM The block scheme of the INS/GPS/BAR/MAG integration is shown in Fig. 4. In accordance with a solution for INS errors damping, the errors model, used in the presented integrated scheme for INS errors estimation, is given in [1]. The variables given in Fig. 4 are defined as: SVn, SVe, SVd are linear velocity errors NED components, Sy, SA, Sh are position errors components, фп, фе, фй present orientation errors NED components for calculated platform, ю*are gyro drifts, w(0 is Gaussian white noise, while Bn, Be, Bd present NED components of the accelerometer biases. Based on the estimation error model, the correction of the velocity components in the NED coordinate system, is defined as: Vc = V™ -SVn,Vec = = vr -SjfV = VdNS -SVd, (17) where SVn, SV e, S V d are estimates of velocity errors defined at the output of EKF. The position correction is given as: =Vms -S(p,Ìc = ÌNS-SÌ,h = hNS-Sir, (18) where S(p SX, S h are estimates of position errors for the LLh (latitude, longitude, and height) coordinates. The state vector in EKF, in our case, is given as, x = [8гф 8rx 8rh 8V„ 8Ve 8V? Sep 8в 8V Bn Be B? af a? a? ]. (19) The states in the EKF with control signals can be represented as: i + = Ф, i+-1 + u k + K k (z k - H k Фк i+-1 - H k Uk ), (20) where u is matrix defined as: Uk-1 = diag ([° 0 Uh K K U K Ue °]) > (21) while Kk is the Kalman gain matrix, and Ф is the state transition matrix. The control signals are given as: «I = K1nKF Sinh(SVn ), Ue = -KleKF Sinh(SVe ^ «I = -KnKF SÌnh(5Ve ), « = -K2eKF Siuh(5Vn ), uh = -C1KF sinh(5h),uvd = -C2KF sinh(5h). (22) Fig. 4. Block scheme of the observed INS/GPS/BAR/MAG Integration The Kalman gains are time varying, with the steady state defined for the equilibrium between the process noise variance accumulation and the measurement noise variance. These steady-state gains in the Kalman filter are the same as the optimum K and C coefficients used for error damping in the INS. It should be noted that the adaptation of these coefficients is performed only when GPS data are not available. In the case of GPS drop-out, EKF is used to perform prediction with the Kalman gain matrix formed as zero matrix (K = 0). 5 EXPERIMENTAL RESULTS The integrated navigation system used in this study is formed by using the low cost INS MEMS sensors: the MPU-60X0, which combines a 3-axis gyroscope and a 3-axis accelerometer, "Gms-u1LP" GPS receiver (L1 C/A code, updating frequency 10 Hz) and the MS5611-01BA03 as a barometric pressure sensor. We assumed the non-orthogonality of the INS sensors as defined in product specifications. The stochastic noise characteristic of these sensors was determined by the Allan variance method based on a three-hour measurements. These data are used for the sensor calibration and for the covariance matrices used in EKF initialization. The results of those tests are given in Table 1. Experimental testing of the integrated navigation system was performed by using a land vehicle. The device was installed on the center of mass of the vehicle and the antenna of the GPS receiver was placed on the top of the vehicle, wherein the applied software solution was used to adapt to the position differences between the GPS antenna and the inertial sensors. The vehicle was moving along the predetermined trajectory, in the urban environment, which contains ups and downs, and which is defined with twenty CP. Table 1. Sensors error parameterization Acc Bias [m/s2] White noise [m/s/s1/2] Random walk [m/s2/s1/2] X -11.5x10-2 1.77x10-2 6.48x10-4 y -1.66x10-2 1.86x10-2 5.73x10-4 z 9.8x10-1 1.87x10-2 1.94x10-3 Gyro Drift [rad/s] White noise [m/s/s1/2] Random walk [m/s2/s1/2] X 9.28x10-3 11.3x10-3 3.7x10-4 y 2.8x10-3 10.5x10-3 3.04x10-4 z 1.02x10-4 10.2x10-3 6.07x10-5 Baro rc [s] White noise [m/s/s1/2] Random walk [m/s2/s1/2] 200 9.31x10-3 13.9x10-2 The covariance matrix R, is formed based on the 12 hours long GPS measurements at the known position and is given as: Ru = 5.7296e-8; R2,2 = 4.9e-8; R3,3 = 1.55; R4,4 = 0.01; R5 5 = 0.01; R66 = 0.01. Covariance matrix Q is determined based on the sensors' noise characteristics and is given as: Q33 =1e-4; Q44 = 4.77e-5; Q55 = 2.75e-6; Q66 = 7.3e-3; Q7j = 1.5e-3; Q8,8 =1.5e-3; Q9>9 = 1e-3; Q^q = 1.2e-3; Q1U1 = 1.2e-3; Q^,^ = 1e-4; QB,13 = 1e-4; QH14 = 1e-4; Q15,15 = 1e-4; The other elements of the matrices R and Q are zero. The paper presents the test results obtained for a period of 300 seconds. All system settings were set in stationary conditions. In Fig. 5, the accelerometers' row data acquired during the test are shown, as well as the results after filtering performed based on the parameters of the stochastic sensors. Samples [ Fig. 5. Accelerometers measurements obtained during the test In Fig. 6, the results of the gyroscopes measurements, acquired during the test are shown, with the results of these measurements after filtration. In these figures, in both cases, the large influence of noise is evident, especially for the low values of the sensors measurements. This fact is particularly noticeable on the vertical axis of gyroscopes triad in Fig. 6. Fig. 6. Gyroscope measurements obtained during the test The normalized values of the magnetometer measurements after calibration, in the horizontal plane, are shown in Fig. 7. As the graph describes a full circle, it is clear that the hard iron measurement errors are removed when the calibration method is applied. Deviations of the graph from the exact circle are a consequence of the soft iron errors, which have not been entirely removed. Fig. 7. Normalized magnetometer measurements in the horizontal plane, after the calibration, acquired during the test The results of pressure and temperature sensor measurements during the test are shown in Fig. 8. In this figure, the black lines represent the filtered measurements. Bearing in mind the noise characteristics of the pressure sensor, shown in Table 1, there is no obvious strong presence of noise in Fig. 8, because the dominant barometer error is random walk noise. The step graph of temperature row data, shown in the lower graph in Fig. 8, is a consequence of sensor resolution of 0.1 °C. Fig. 8. Measurements of the barometer and temperature sensors, before and after filtration A profile of the trajectory along which the test was performed is shown in Fig. 9. The CPs are selected in such a way to allow full control of determining the position of the vehicle during its movement, e.g. at critical points on the road, such as the changes of ups and downs and curves. 1'7|(^486 20.488 20.49 20.492 20.494 20.496 20.498 20.5 Longitude [deg] Fig. 9. The trajectory of vehicle movement in the horizontal plane, with marked control points Horizontal and vertical profile of vehicle trajectory in NED coordinates, relative to the GPS measurements are given in Figs. 10 and 11. During these tests, there was no interruption of the GPS information. The vertical profile of the trajectory is particularly interesting because of the large deviations of GPS measurements relative to the CPs, which resulted with the large covariance measurement error of the altitude. Fig. 10. The north and east coordinates, of the vehicle trajectory in NED coordinate system (when GPS is always on) The vehicle velocity in the NED coordinates is shown in Fig. 12. As the updating frequency of GPS data is 10 Hz, the likelihood of a large navigation error accumulation between the GPS measurements is very low (negligible). Furthermore, with the well-defined EKF parameters, the measurement noise is removed, as it is evident in Fig. 12. Fig. 11. The altitude profile of vehicle trajectory, (when GPS is always on) In the upper graph in Fig. 13, the heading angle of the vehicle is shown. The angle was determined based on full integration of gyros, accelerometers, and the magnetometer, where the angle change is in accordance with the trajectory of the vehicle. The deviation of yaw angle for different cases of sensor integration is shown in the lower graph in Fig. 8. All these deviations are reduced to a small value by the use of additional external measurements, especially those of magnetometer. о Fig. 12 Velocity profiles of the vehicle in NED coordinate system, (GPS on) The particularly large positive contribution of gyroscopes and magnetometers integration can be seen through the determination of the pitch angle, as given in Fig. 14. When observing the results of measurements relative to the specific angles defined with CPs, it is clear that the best results are obtained when the gyroscopes drift suppression is done by using the magnetometers and accelerometers. The CPs can only approximately determine the pitch angle; therefore, it is not possible to accurately determine the measurement error. 200 ~~ Gyro/Acc/Mag -200 0 50 100 150 200 250 300 S Fig. 13. The flow of the yaw angle during the experimental testing 20 CP 50 100 150 200 250 300 t [s] Fig. 14. The flow of the pitch angle during the experimental testing The improvement of the roll angle accuracy is evident in Fig. 15, which shows the transition of the roll angle along the experimental trajectory. By observing the roll and pitch angle curves, it is evident that there are large divergences in cases in which the angles are determined only based on the gyroscope measurements, which is the consequence of the gyro sensors imperfections. A partial gyroscope drift compensation is achieved by the introduction of accelerometer measurements into the orientation algorithm. However, accelerometers do not produce perfect data and therefore gyroscopes drift cannot be absolutely removed. The largest deviations of the pitch and roll angles occur when there are changes of direction and velocity of the vehicle, which is in accordance with Figs. 10 to 12. The best possible results are achieved with the integration of all three sensor types. 30 ---Gyro 50 100 150 200 250 300 t M Fig. 15. The flow of the roll angle during the experimental testing The graphics of the vehicle position errors in the horizontal and vertical plane for the period of GPS information absence (artificially introduced between 180 s to 210 s) for a period of 30 seconds are shown in Fig. 16. The self-damping algorithm of INS measurements is realized as described in Section 3. The rapid increasing of the position error is typical for the MEMS sensors. The curves in Figure 16 show the case when there are linear acceleration of the vehicle, in all three directions, during the manoeuvre. The maximum, mean and root mean square error (rms) of measurements are shown in Tab. 2. The case in which there are linear and angular accelerations of the vehicle, during the GPS outage, exist in a period between 150 s and 160 s. In this case, there is a large error in determination the position and velocity of the vehicle. The errors of the position and velocity in this case are also shown in Table 2. For comprehensive and unequivocal monitoring of navigation parameters, the velocity deviations, during the GPS data outages, are also displayed for the time interval between 180 s and 210 s in Fig. 17. The errors of vehicle velocity are also shown in Table 2. If we look at the graphics in Figs. 16 and 17, and the results that are shown in Table 2, it can be concluded that the proposed solutions for the object attitude, position and velocity errors damping allow the successful realization of the real multi-sensor integrated navigation system, which is based on MEMS sensors. I e а z £ lH. a [Ü t „ " E 1 GPS off 30s....... GPS. off 30»....... GPS off 30s t M Fig. 16. Vehicle position errors during the GPS outage Fig. 17. The velocity estimation in NED coordinates, during the GPS outage Table 2. Summary of the experimental results Parameter ■ Error (GPS off 10 s) Error (GPS off 30 s) max mean rms max mean rms N [m] 8.209 4.021 4.809 5.141 1.801 2.320 E [m] 9.301 2.231 3.491 13.894 4.840 6.430 D [m] 0.071 0.024 0.029 0.109 0.035 0.042 Vn [m/s] 1.158 0.042 0.746 1.231 0.497 0.602 Ve [m/s] 1.561 0.031 0.823 1.759 0.882 1 Vd [m/s] 0.123 0.009 0.002 0.135 0.009 0.019 6 CONCLUSIONS The results of the comprehensive analysis, observed in the accordance with the objectives of the research, clearly show that the proposed solutions for the error damping justify the initial assumptions and objectives of the conducted research. It can be also concluded that the usage of the magnetometer for gyro drift compensation through the PI controller usage contributes to the improvement in the attitude determination of the vehicle. Finally, we have shown that the proposed solution for the self-error damping in the horizontal channel produces reliable results in the user position and velocity determination process, although only the low-cost sensors are used. 7 ACKNOWLEDGEMENT This work was done within the research project of the Ministry of Science and Technological Development of Serbia III47029. 8 REFERENCES [1] Salychev, O.S. (2012). MEMS-based Inertial Navigation: Expectations and Reality The Bauman Moscow State Technical University, Moscow [2] Sokolović, V. (2011). Analyse of signal acquisition in GPS software receiver. Vojnotehnički glasnik / Military Technical Courier, vol. 59, no. 1, p. 81-95. (in Serbian) [3] Kuščer, L., Diaci, J. (2013). Measurement Uncertainty Assessment in Remote Object Geolocation. Strojniški vestnik - Journal of Mechanical Engineering, vol. 59, no. 1, p. 32-40, D0l:10.5545/sv-jme.2012.642. [4] Gülal, E., Dindar, A.A., Akpinar, B., Tiryakioglu, I., Aykut, N.O., Erdogan, H. (2015). Analysis and management of GNSS reference station data. Tehnični vjestnik - Technical Gazette, vol. 22, no. 2, p. 407-414, DOI:10.17559/tv-20140717125413. [5] Bian, H., Jin, Z., Tian, W. (2006). IAE-adaptive Kalman filter for INS/GPS integrated navigation system. Journal of Systems Engineering and Electronics, vol. 17, no. 3, p. 502- 508, DOI:10.1016/S1004-4132(06)60086-8. [6] Stancic, R., Graovac, S. (2010). The integration of strap-down INS and GPS based on adaptive error damping. Robotics and Autonomous Systems, vol. 58, no. 10 p. 1117-1129, DOI:10.1016/j.robot.2010.06.004. [7] Noureldin, A., Karamat, T.B., Eberts, M.D., El-Shafie, A. (2009). Performance enhancement of MEMS-based INS/ GPS integration for low-cost navigation applications. IEEE Transactions on Vehicular Technology, vol. 58, no. 3, p. 10771096, DOI:10.1109/TVT.2008.926076. [8] Readdy, G.S., Saraswat, V.K. (2013). Advanced navigation system for aircraft applications. Defence Science Journal, vol. 63, no. 2, p. 131-137, DOI:10.14429/dsj.63.4254. [9] Sokolovic, V., Dikic, G., Stancic, R. (2014). Adaptive error damping in the vertical channel of the INS/GPS/Baro-altimeter integrated navigation system. Scientific Technical Review, vol. 64, no. 2, p. 14-20. [10] Madgwick, S.O.H., Vaidyanathan, R., Harrison, A.J.L. (2010). An Efficient Orientation Filter for IMU and MARG Sensor Arrays. Report, University of Bristol, Department of Mechanical Engineering, Bristol. [11] Sokolovic, V., Dikic, G., Stancic, R. (2013). Integration of INS, GPS, magnetometer and barometer for improving accuracy navigation of the vehicle. Defence Science Journal, vol. 65, no. 5, p. 451-455, D0l:10.14429/dsj.63.4534. [12] Hao, Y., Zhang, Z., Xia, Q. (2010). Research on data fusion for SINS / GPS / magnetometer integrated navigation based on modified CDKF. IEEE International Conference in Progress in Informatics and Computing, p. 1215-1219. [13] Zhao, H., Wang, Z. (2012). Motion measurement using inertial sensors, ultrasonic sensors, and magnetometers with extended Kalman filter for data fusion. Sensors Journal, vol. 12, no. 5, p. 943-953, D0I:10.1109/JSEN.2011.2166066. [14] Han, S., Wang, J. (2011). A novel method to integrate IMU and magnetometers in attitude and heading reference systems. Journal of Navigation, vol. 64, no. 4, p. 727-738, DOI:10.1017/ S0373463311000233. [15] Ahmed, M., Cuk, D. (2005). Strapdown attitude algorithms using quaternion matrix and random inputs. Scientific Technical Review, vol. 55, no. 1, p. 3-13. [16] Kuipers, J.B. (1999). Quaternions and Rotation Sequences: A Primer with Applications to Orbits, Aerospace and Virtual Reality, Princeton University Press, Princeton. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)7-8, 459-464 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2581 Original Scientific Paper Wear Behaviour and Cutting Performance of Surfaced Inserts for Wood Machining Regita Bendikiene1* - Gintaras Keturakis1 - Tilmute Pilkaite1 - Edmundas Pupelis2 1 Kaunas University of Technology, Faculty of Mechanical Engineering and Design, Lithuania 2 Kaunas University of Technology, Laboratory Centre, Lithuania Machining, in general, and wood machining, in particular, are complex to explain and many factors influence the process. Natural wood is a nonhomogeneous biological material, so each species of wood requires different machining conditions and tools. Understanding the properties of wood and choosing the proper cutting tools and machining parameters can improve the quality of wooden products, increase production efficiency and improve machining. The objective of this study was to find the most suitable cutting tools for the machining of oak wood grown in Lithuania. This paper describes tests of two commercial cutting inserts made of high speed tool steel (CT01M-LA2 and 8X6FT) and two experimental inserts (EI) made by surfacing using a submerged arc welding (SAW) technique and a mixture of alloying elements spread on the surface under industrial flux. The results obtained from the milling tests demonstrated the suitability of the suggested surfaced inserts for machining oak wood. All tool wear parameters, such as cutting edge radius, edge recession, nose width, and cutting power, have been evaluated and compared. The cutting edge radius of EI after 3200 m of cutting length was ~ 5.00 /um, 3 to 4 times lower than the wear of standard tools; similarly, the edge recession of the surfaced inserts showed 2 to 3 times lower results. The results of this study indicate that the cutting power increases as the feed per cutter increases. These findings suggest that surfaced inserts can replace the commercial inserts used for wood machining. Keywords: wood milling, surfacing, edge recession, tool wear, oak Highlights • Cutting tool inserts were made by surfacing using a submerged arc welding (SAW) technique. • Wear behaviour of the presented inserts was compared with commercial inserts. • Tests were performed on samples of oak wood. • Surfaced inserts showed better wear behaviour and cutting performance. • The suggested technology can be used for manufacturing wood machining tools. 0 INTRODUCTION The literature has emphasized the importance of both wood cutting technology and the natural properties of wood. Machining of different species of natural wood (oak, pine, birch, etc.) and timber products requires the usage of high quality cutting tools. Each type of material requires different machining conditions, so it is not possible to use one type of cutting tool to achieve the best efficiency. A number of authors have claimed that high quality tools are necessary due to the particular properties of the natural materials machined, such as the possible significant anisotropy of the structure and the cutting behaviour of the material in different directions (along or across to the fibres), large sizes (for the initial processing from the logs), unexpected changes in the structure and a sudden rise in the internal stresses of the material machined, the appearance of hard and brittle particles and changes in density in a cross section, etc. [1] and [2]. Cutting tools employed in the field of wood machining do not allow the full potential of the machines to be fulfilled [3]. Several reports have shown that the main problems in wood processing arise due to differences in the physical and chemical structure of wood and metal. Firstly wood has relatively good machinability allowing high cutting and feed speed; but it contains some water, making it very corrosive [3] and [4]. Secondly, the natural defects (knots, wane, particles of hard mineral contamination, etc.) present in the timber can initiate blunting of the cutting edges; therefore, very hard but brittle materials are not suitable as cutting tools in this case. Hardened steels, high speed tool steels, stellites, tetrahedral amorphous carbon, composites of titanium carbides and polycrystalline diamond wood cutting tools are currently used in the wood industry; among them, the most common are cemented carbides, because of their good wear resistance and relatively low cost compared to diamond based tools [1] and [5]. All of the studies reviewed here support the hypothesis that cutting tools manufactured from high speed tool steel or surfaced with high wear resistance coatings can replace expensive tools made of sintered carbides or sintered polycrystalline diamond [6]. Finally, machining of wood is performed under a very high working speed and extremely sharp *Corr. Author's Address: Kaunas University of Technology, Studentu str. 56, Kaunas, Lithuania, regita.bendikiene@ktu.lt 459 cutting edges are needed. Thus far, previous studies have reported that the main wear mechanism is the erosion of the cutting tool material; hence, coatings for woodworking tools should be very hard, adherent and smooth and exhibit high wear resistance [5] and [7]. Another problem is the rounding of the cutting edge radius during the cutting process. Even a small increase in the edge radius leads to increased tool wear and worse surface quality of the final wood product. Thus this paper presents the results of wear behaviour of two standard cutterhead inserts, made of high speed tool steels for wood cutting, supplemented by test results of experimental inserts (EI) surfaced using submerged arc welding (SAW) technique and a properly chosen chemical composition of alloying flux. 1 METHODS The base material for surfacing was cheap plain carbon steel (C 0.14 % to 0.22 %; Si 0.12 % to 0.13 %, Mn 0.4 % to 0.65 %, S< 0.05 %, P < 0.04%) provided as an 8 mm thick plate. The surfacing process was performed on 40 mm x 100 mm samples in a single pass using the SAW technique with alloying materials mixture (~ 6 g) spread on the surface under the flux. The chemical composition of the materials mixture is presented in Table 1. A single 1.2 mm diameter electrode low carbon wire (C < 0.1 %; Si < 0.03 %, Mn 0.35 % to 0.6 %, Cr < 0.15 %, Ni < 0.3%) was used for the surfacing. The SAW was carried out with an automatic welding device (torch MIG/MAG EN 500 78), with welding parameters: welding current 180 A to 200 A, voltage 22 V to 24 V, travel speed 14.4 m/h, and the wire feed rate 25.2 m/h. Table 1. Chemical composition of the spread materials mixture Composition of materials mixture [wt.%] "SIC Cr W Fe-70%Mn Flux EI 40 10 40 10 AMS1* * LST EN 10204:2004 SiO2 and MnO > 50 % A blended powder of materials was spread on the surface of the base metal and fused using a metal arc. Additional standard flux AMS1 (GOST 9087-81; SiO2 38 % to 44 %, MnO 38 % to 44 %, CaF2 6 % to 9 %, CaO < 6.5 %, MgO < 2.5 %, Al2O3 < 5 %, Fe2O3 < 2 %, S < 0.15 %, P < 0.15 %) was used to shield and to protect the welding area. The presence of chromium in the surfaced layer affects the formation of the retained austenite, thus slowing the decomposition of austenite, since chromium provides some corrosion resistance. Silicon carbide (SiC) was used as a deoxidizer in the welding flux. Deoxidizers react with oxygen at the welding temperature and significantly decrease the quantity of oxides in the bead, thus increasing the quality of the weld. Adding silicon into the flux improves the metal mass transfer coefficient and the form of the weld, as well as modifying the slag [7]. One testing lot of surfaced inserts were heated to 1100 °C afterwards and hammered in order to examine the influence of plastic deformation on the wear properties of the tool. During hot hammering, the face of the surfaced samples was plastically deformed and flattened to the level of the base metal and, as a result, the time of insert machining was reduced (no additional cutting needed). The second positive outcome of the smithing process was the self-hardening of inserts in the air, i.e. tempering following hot plastic deformation. The hardness of surfaced inserts reached 60 HRC, while the surfaced and hammered reached 50 HRC. After tempering at 500 °C, the hardness values were changed to 57 HRC and 55 HRC, respectively. The mechanical behaviour of surfaced experimental inserts and commercial inserts was assessed in terms of hardness and wear properties. Table 2 shows the chemical composition of commercial inserts. Hardness measurements of the layers were accomplished on the wrought and heat-treated (tempered) inserts using Rockwell tester TK - 2 at a load of 1470 N using a diamond indenter. Table 2. Chemical composition of commercial inserts Steel Quantity of elements [wt %] C Cr Mo Ni Ti V W Co 1* 0.80 5.87 0.91 0.10 0.20 2** 1.00 3.83 0.21 0.12 2.60 12.55 0.12 * 8X6FT (GOST 6567-75) ** Freud CT01M-LA2 The most important characteristics (Fig. 1) selected to define the wear behaviour of inserts were: cutting edge radius p [mm], edge recession Aa, [mm], nose width b [mm], and cutting power P [W] [8]. The actual values of the edge rounding radius were assessed using a lead imprint method with a Nikon Eclipse E200 optical microscope and Lumenera Infinity 1 digital camera. Infinity Analyze Release 5.0.2 software was used to analyse and evaluate the obtained results with an accuracy of ± 2 mm.The experimental results were subjected to statistical analysis. Cutting power P was determined by measuring the available power and taking out the idle motion power. Available and idle motion power were measured with an accuracy of ± 10 W on a universal power tester K506. Table 3. Characteristics of inserts Fig. 1. The main geometrical parameters of the cutting tool: cutting angle 8, sharpness angle ß, rake angle y, and clearance angle а The values of the cutting edge radius, edge recession, nose width, and cutting power were recorded and measured at intervals of cutting length L: (0; 50; 100; 150; 200; 400; 800; 1200; 1600; 2400, and 3200) m. Each value at every specified cutting length was an average of 5 tests. 2 EXPERIMENTAL Experimental inserts were strait sharpened and their edges were converged according general grinding procedures for inserts. The sequence of inserts manufacturing is presented in the Fig. 2. a) b) c) Fig. 2. Sequence of experimental inserts preparation: a - blank; b - surfaced; c - insert Characteristic EI1 * EI2 Ihrdnr-- lirC Not tempered 50 60 nardness, hrc After tempering 55 57 Dimensions [mm] 60x30x3.55 Sharpness angle ß [degree] 40 Weight of insert m [g] 45.69 45.19 Roughness of rake face Ra [^m] 0.135 0.152 Roughness of clearance face Ra [^m] 0.083 0.066 * hammered after surfacing. Ten wood test samples were prepared from oak wood (Quercus robur) grown in Lithuanian (Table 4) with dimensions of 1000 mm x 100 mm x 45 mm. Special care was taken to select samples as free as possible of knots or other defects. The characteristics of experimental inserts are presented in Table 3. Hereafter accurate dimensions of inserts were ensured by measuring inserts using electronic callipers with an accuracy of ± 0.001 mm. The surface roughness tester, profilotemer Mahr MarSurf PS 1, was used to evaluate the roughness of the rake face and clearance face [9]. Table 4. Physical properties of Quercus robur Average Average number Average width Average moisture of annual rings of annual ring density content ю [%] per 1 cm [mm] [kg/m3] 10.2 3.00 3.33 690 The average moisture content was estimated using a Gann Hydrometer Compact A electronic moisture tester with an accuracy of ± 1 %. The number of annual rings per 1 cm was determined by counting the rings in the end section perpendicular to the wood fibers [8]. Samples were weighted on electronic scales (accuracy ± 0.01 g) for determination of density. Average ambient temperature of wood samples testing was 18 ± 2 °C, while relative air humidity was 60 ± 5 %. Table 5. Milling test conditions Name Values EI1 EI2 Cutting speed nc [m/s] 31 Feed per insertf [mm] 1.00 0.50 Feed speed Vf [m/min] 6 3 Depth of milling h [mm] 2 Milling width b [mm] 45 Diameter of cutterhead d, [mm] 103 Number of inserts z [unit] 1 Cutting angle 8 [degree] 60 The wear performance of experimental and standard inserts was carried out on a typical industrial thickness planer (SR3-6) with a face milling cutterhead using oak samples as the workpiece. The milling was conducted according longitudinal milling, with vectors of cutting speed nc and feeding speed nf. Milling conditions were the same for each of the tested inserts and are shown in Table 5. The inserts made of different steel grades and surfaced layers were replaced for each test. The cutterhead was designed to have two cutting edges, to avoid the imbalance that can appear with one insert, however while two experimental inserts were mounted [10], only one was tested. An indirect method was used to change the thickness of chips through the feed per one insert f = 0.50 mm and 1.00 mm. The cutting speed was constant for all tested samples nc = 31 m/s. The rotating frequency of the cutterhead measured with an accuracy of 10 min1 using a Tachometer Stroboscope SC-5 was n = 5790 min1. 3 RESULTS AND DISCUSSION The wear measurement was based on the determination of the edge recession (cutting edge radius, edge recession, and nose width) after each defined cutting length (effective cutting path of the blade). The wear test results of two experimental and two commercial inserts showed that there are three phases, which characterize the evolution of the insert wear recession: running (intensive wear), linear wear (stable), and vital wear (tool failure). It can be stated that for each type of insert, the running period was up to 800 m of cutting length due to intensive wear [2]. The results, as shown in Fig. 3, indicate that further wear evolution of inserts was relatively linear or stable. 0 800 1600 2400 3200 Cutting length [m] Fig. 3. Variation of the cutting edge radius for the feed per insert (fz = 1.0 mm) In Fig. 3 there is a clear trend towards an increase in the inserts cutting edge radius with increasing cutting length. The cutting edge radius of the surfaced and subsequently hammered insert EI1 showed the lowest wear evolution when compared with commercial inserts: 5.22 цт. The results of the cutting edge test are in line with those of the previous test, as the tendency of wear of the EI2 inserts was the same: 5.35 цт. The cutting edge radius of the industrial insert made of CT01M-LA2 was 11.2 цт, while the maximum wear values were achieved on 8X6FT inserts with 19.8 цт. Previous studies have attempted to explain why tools with a cutting edge radius of more than 25 цm cannot be used for machining any more [5]. The smallest edge recession was noticed when testing surfaced insert EI2 (Fig. 4). о 0 J---- 0 800 1600 2400 3200 Cutting length [m] Fig. 4. Edge recession for the feed per insert (fz = 0.5 mm) 1 ^ 'pC г J"=-- —<> —•—Ell □ 0 800 1600 2400 3200 Cutting length [m] Fig. 5. Nose width growth of surfaced inserts There was no significant difference between edge recession of surfaced inserts EI1 and EI2 (10.66 цm and 10.75 цm respectively), while the edge recession of industrial inserts was ~ 20 to 35 цm. Overall, these results indicate that the edge recession of experimental tools was 2 to 3 times lower. As shown in Fig. 5, the nose width test was used to analyse tool wear. The milling test was interrupted at defined intervals of cutting length in the same way as for previous tests. Average tool wear or blunting can be defined as the difference between the resultant nose width and the initial nose width. Cutting tool temperature is another important factor affecting tool wear in wood machining, because the hardness, toughness, and chemical properties of tool material degrade when the tool's temperature increases [11]. Continuous plastic deformation and shear during chip formation generates thermal energy and friction, which appear on the rake and clearance face of the tool, at the same time there is also friction between the sample and the back face of the tool. The heat generated is transferred to the cutting tool and work sample. This heat has a negative effect on the quality and accuracy of the machined products and on the main parameters of cutting: cutting speed, depth of cut, blunting and cutting power. Consequently, the cutting power of all inserts was tested over the whole cutting length (Figs. 6 and 7 and Table 6). Fig. 6. Cutting power of inserts forfz = 1.0 mm The values of cutting power up to 400 m of cutting length were high because of crumbling of the top of blade's cutting edge [12]. At this stage of machining the cutting edge radius grew rapidly as well. Linear or stable cutting power intensity was observed for all inserts at cutting lengths from 1200 m to 3200 m. Wear by crumbling of inserts blades was displaced by a plastic wear phase. Surfaced inserts EI1 and EI2 showed very similar results for cutting power, therefore they can be used for machining of wood. Table 6. Cutting power of inserts (W) (fz = 0.5 mm) Cutting length [m] 50 100 150 200 400 800 1600 3200 EI2 200 220 200 220 200 162 180 180 1 * 218 198 175 183 175 200 220 216 2** 165 150 165 132 140 132 139 166 1 * 8X6FT; 2 ** CT01M-LA2 Cutting 1 Fig. 7. Cutting power of inserts for fz = 0.5 mm 3 CONCLUSIONS The lowest cutting edge radius occurred on surfaced and additionally plastically deformed insert EI1; the cutting edge radius after 3200 m of cutting length was 5.22 цт, while for unhammered EI2 it was 5.35 цт. In summary, these results show 3 to 4 times lower wear than commercial tools. Better wear performance was achieved by testing the edge recession of surfaced inserts EI1 and EI2. In summary, edge recession of the suggested experimental tools was 2 to 3 times lower. The most obvious finding to emerge from the analysis is that the relatively hard coatings (55 to 57 HRC) surfaced on soft plain carbon steel can replace some commercial inserts made of high speed tool steels for use in oak wood machining, thus reducing friction and wear of the wood cutting tool. 4 REFERENCES [1] Warcholinski, B., Gilewicz, A., Ratajski, J. (2011). Cr2N/ CrN multilayer coatings for wood machining tools. Tribology International, vol. 44, no. 9, p. 1076-1082, D0l:10.1016/j. triboint.2011.05.004. [2] Vobroucek, J. (2015). The influence of milling tool geometry on the quality of the machined surface. Procedia Engineering, 25th DAAAM International Symposium on Intelligent Manufacturing and Automation, vol. 100, p. 1556-1561, D0I:10.1016/j.proeng.2015.01.528. [3] Faga, M.J., Settineri, L. (2006). Innovative anti-wear coatings on cutting tools wood machining. Surface and Coatings Technology, vol. 201, no. 6, p. 3002-3007, DOI:10.1016/j. surfcoat.2006.06.013. [4] Djouadi, M.A., Beer, P., Marchal, R., Solokowska, A., Lambertin, M., Precht, W., Nouveau, C. (1999). Antiabrasive coatings: application for wood processing. Surface and Coating Technology, vol. 116-119, p. 508-516, D0l:10.1016/S0257-8972(99)00236-4. [5] Endler, I., Bartsch, K., Leonhardt, A., Scheibe, H.J., Ziegele, H., Fuchs, I., Raatz, Ch. (1999). Preparation and wear behaviour of woodworking tools coated with superhard layers. Diamond and Related Materials, vol. 8, no. 2-5, p. 834-839, D0l:10.1016/S0925-9635(98)00359-8. [6] Gilewicz, A., Warcholinski, B., Myslinski, P., Szymanski, W. (2011). Anti-wear multilayer coatings based on chromium nitride for wood machining tools. Wear, vol. 270, no. 1-2, p. 32-38, DOI:10.1016/j.wear.2010.09.002. [7] Ambroza, P., Bockus, S., Kavaliauskiene, L. (2013). Formation of build up layers microstructure by arc automatic overlay welding using secondary raw material powders. Archives of Metallurgy and Material, vol. 58, no. 2, p. 549-553, DOI:10.2478/amm-2013-0034. [8] Keturakis, G., Lisauskas, V. (2010). Influence of the sharpness angle on the initial wear of the wood milling knives. Materials Science (Medžiagotyra), vol. 16, no. 3, p. 205-209. [9] Keturakis, G., Juodeikiene, I. (2007). Inverstigation of Milled wood surface roughness. Materials Science (Medžiagotyra), vol. 13, no. 1, p. 47-51. [10] Aknouche, H., Outahyon, A., Nouveau, C., Marchal, R., Zerizer, A., Butaud, J.C. (2009). Tool wear effect on cutting forces: In routing process of Aleppo pine wood. Journal of Materials Processing Technology, vol. 209, no. 6, p. 2918-2922, DOI:10.1016/j.jmatprotec.2008.06.062. [11] Darmawan, W., Gottlober, Ch., Oertel, M., Wagenfuhr, A., Fischer, R. (2011). Performance of helical edge milling cutters in planning wood. European Journal of Wood and Wood Products, vol. 69, no. 4, p. 565-572, DOI:10.1007/s00107-010-0517-8. [12] Horman, I., Busuladzic, I., Azemovic, E. (2014). Temperature influence on wear characteristics and blunting of the tool in continuous wood cutting process. Procedia Engineering, 24th DAAAM International Symposium on Intelligent Manufacturing and Automation, vol. 69, p. 133-140, DOI:10.1016/j.proeng.2014.02.213. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)7-8, 465-470 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2511 Original Scientific Paper Towards Sophisticated Control of Robotic Manipulators: An Experimental Study on a Pseudo-Industrial Arm Jamshed Iqbal12* - Muhammad Imran Ullah2 - Abdul Attayyab Khan3 - Muhammad Irfan4 1 National University of Computer and Emerging Sciences, Islamabad, Department of Electrical Engineering, Pakistan 2 COMSATS Institute of Information Technology, Islamabad, Department of Electrical Engineering, Pakistan 3 University of Genova, DIBRIS, Italy 4 University of Engineering and Technology, Taxila, Department of Electrical Engineering, Pakistan Robotic manipulators have reshaped industrial processes. The scientific community has witnessed an ever increasing trend in robots deployed to accomplish various tasks in industry The complex nature and constrained requirements of robots may demand non-trivial control approaches. This paper deals with the design, simulation and hardware realization of two sophisticated control strategies: computed torque control (CTC) and variable structure control (VSC) on a pseudo-industrial manipulator with six degrees of freedom (DOF). Based on the derived kinematic and dynamic models of the robot, control laws have been formulated, which are then subjected to various test inputs in a simulation to characterize the tracking performance. The simulation results were then validated by implementing control laws on a custom-developed pseudo-industrial autonomous articulated robotic educational platform (AUTAREP). The experimental results show the effectiveness of the control strategies to track a desired trajectory. Keywords: robot control, manipulator, robust laws, industrial robots Highlights • This paper presents advanced strategies to control a highly non-linear system like a multiple Degree Of Freedom (DOF) robotic arm. • The strategies include computed torque control (CTC) and variable structure control (VSC). • Design parameters using both strategies have been investigated in a simulation. • The strategies were carried out on a custom-developed autonomous articulated robotic educational platform (AUTAREP). • Trajectory tracking results showed that the derived laws can effectively track the desired reference input for both strategies. 0 INTRODUCTION Today, robots are being deployed to accomplish tasks having strict requirements of accuracy, precision, repeatability, mass production and quality. A major breakthrough was reported with the advent of feedback control systems and self-correcting mechanisms. The development of multiple degree of freedom (multi-DOF) manipulators has contributed significantly to modern robots. Industrial robots are primarily multi-DOF anthropomorphic manipulators. The past few years have seen a large rise in the use of industrial robots and this trend will likely continue, as highlighted in [1]. The performance of a robotic manipulator is characterized by a well-defined control approach [2]. Classical or trivial control strategies are usually based on linear control laws while modern approaches are nonlinear in nature. The approach to controlling a multi-DOF manipulator must be robust enough to cope with the effects of inherent nonlinearities and coupling in the robot dynamics [3] to [5]. Classical approaches suffer from various issues, which can be avoided by merging these with advanced control strategies or using an advanced strategy in standalone. Control of multi-DOF robotic manipulators is a vital research area today. However, most of the reported research is either limited to the implementation of linear control approaches or the simulation of sophisticated control strategies. In contrast, the present work investigates advanced approaches like computed torque control (CTC) and variable structure control (VSC) from a simulation viewpoint as well as its physical realization on a custom-developed autonomous articulated robotic educational platform (AUTAREP). 1 CONTROL DESIGN The study of manipulators for diversified applications has highlighted the need for sophisticated algorithms for their control and trajectory planning. The objective in the design of robotic manipulators is to control both the position and the orientation of the tool in a 3D workspace. The scientific community has reported both classical and robust strategies for robot control. *Corr. Author's Address: Department of Electrical Engineering, National University of Computer and Emerging Sciences, A.K Brohi Road, H-11/4, Islamabad, Pakistan iqbal.jam@gmail.com With regard to classical approaches, Iqbal et al. have proposed proportional integral derivative (PID) controllers for mobile robots [6] and multi-DOF serial link robotic exoskeletons [7] and [8]. The role of PID control in industrial automation has been presented in [9], which formulates a nonlinear PID control law to ensure global asymptotic stability (GAS). Classical approaches, when combined with modern control strategies, improve transient response in uncertain scenarios as highlighted in [10]. Combining VSC with PID and adaptive control strategies, Jingmei et al. have improved precision in trajectory tracking of a robotic manipulator [11]. Here, chattering has been reduced via increased system response time. Tahir and Jaimoukha have proposed a model predictive robust controller [12] for linear discrete-time systems subjected to polytopic constraints and bounded disturbances. The proposed control approach is novel in that the outer controller incorporates a state-feedback structure where feedback gains are considered as decision variables in online optimization [13]. PID has been the main workhorse in the industrial sector. However, researchers have recently shown an active interest in the development and applications of nonlinear control methodologies applied to robotic manipulators. A comprehensive review of control strategies for manipulators is reported in [3]. The overall control problem consists of kinematic and dynamic modeling, followed by the design of a control law. The kinematics of the AUTAREP manipulator has been derived in [14] using Denavit-Hartenberg (DH) representation while the system dynamics has been modeled in [15] using Euler Lagrange equation (Eq. (1)). T = M(q,q)q + V(q,q) + G(q) + f (q), (1) where M (q,q) is a 4*4 inertia matrix, V (q,q ), G (q ) and f (q ) are 4*1 vectors of Coriolis centrifugal force, Gravitational force and Frictional force, respectively. t is the 4*1 torque vector applied to the joints of the robot and q, (j and q are 4*1 vectors for angular position, velocity and acceleration respectively. This paper deals with two modern control strategies namely CTC and VSC. 1.1 CTC CTC is a special type of feedback linearization technique with symmetric, constant and positive definite controller gains. CTC can be utilized effectively in cases where there are known nonlinear dynamic parameters and uncertainties. For the CTC law, the expression for the manipulator system in Eq. (1) can be written as Eq. (2). t = M(q,q)(qd -2U-X2e) + V(q,q) + G(q), (2) where the vector q = [q1 q2 q3 q4]T corresponds to the first four joints of the manipulator. qd , qd and qd are the desired joint angle position vector (qd = [qd1 qd2 qd3 qd4]T) and its 1st and 2nd derivatives respectively. e=q - qd is the error signal with e as its 1st derivative. t is the required control output and is represented by t = T т2 т3 t4]t. The gain matrix k=diag X2 A3 A4} can be used to alter the system dynamics. 1.2 VSC VSC has the potential to eliminate uncertainties and disturbances present in the system. A switching surface is designed and the main task of the controller is to drive the system states to this surface. The system then remains on the switching surface to reduce disturbances and modeling uncertainties. VSC is a robust control strategy using high frequency switching control to alter the dynamics of the nonlinear system. The VSC law for the robotic manipulator, using the sliding manifold S = è + Ce, is given in Eq. (3). T = M(q,q)(qd -Cè)+V(q,q)+G(q)--K sgn(Cè + è), (3) where the matrices C=diag {c1 c2 c3 c4} and K = diag {k1 k2 k3 k4} are the switching gain constant and sliding surface constant, respectively, and can be changed to alter the system dynamics. 2 SIMULATION The controller s-functions were developed in a MATLAB/Simulink based on the derived dynamic model of the manipulator. Various desired trajectories including step, sinusoidal and ramp have been applied to analyze the robustness and effectiveness of the proposed control laws. The overall effect of the plant has been investigated by choosing different values for the system constants, i.e. k for CTC and C, K for VSC. In the case of CTC, varying k revealed that the performance of the system improves by increasing the value of k. For example, when considering the second joint (shoulder), a simulated step response is used to investigate the effect of the assigned X2 values. Keeping Àu A3 and Л4 constant (unity), the step response of the shoulder joint for X2 = 1,2, 3, 4 is illustrated in Fig. la with the corresponding torques plotted in Fig. 1b. No overshoot is observed in any case. However, a significant difference in rise time and settling time is noticeable. When X2 is increased from 1 to 4, the rise time is reduced from 7.13 s to 1.97 s while the settling time is reduced from 5.20 s to 1.86 s at ±5% of the desired joint angle. The initial magnitude of torque (25.9 Nm) is required to keep the shoulder joint at its initial position against the gravitational force. The final magnitude of torque is 12.9 Nm. - 0.8 0.6 0.4 0.2 -—Lamda : — ~~ a) ■*-1 .amda2 2 mua_ 1 / 1 )esircd 4 6 Time [s] 10 b) ,amda2 [ ,amda2 4 V'" •*— I.amda2 1 .amda2 1 g -10 § -15 ^ -20 .2 <"-25 -30 ~ 0 2 4 6 8 10 Time [s] Fig. 1. Shoulder joint a) CTC step response for different values of K2 b) corresponding torques Where all joints are moving simultaneously, the same value (X = 4) has been set for each joint of the manipulator. The corresponding step responses and plots of torque are shown in Fig. 2. With a slight difference in the behaviour of various joints, it is inferred from the results that the various joints exhibit different torque requirements. The shoulder joint requires final torque of 5.4 Nm due to the movement of other joints. The initial and final torque requirements of the base joint are 0 Nm due to the zero gravitational effect. The same requirements in case of the elbow joint are 5.9 Nm and 3.1 Nm, respectively. Having few gravitational effects in wrist joint, initial and final torque requirements are 0 Nm and 0.2 Nm, respectively. The equivalent control is designed such that the states of the system are on the sliding manifold, which is defined as S = è + Ce where C must satisfy the Hurwitz condition, i.e. C > 0. Considering first joint (base) for the sake of brevity, the Lyapunov function, l1 = 0.5 si restricts k1 > 0 The corresponding entries of C and K matrices cannot be selected independently without violating the above conditions [16] due to the joints' coupling effects. i ■ 0.8- '0.6 I 0.4 0.2- Shoulde Waist a) f/*— Wrist ///*—Klbow *—f— 1 )csircd 4 6 Time [s] 10 Юг -10 -20 -30, «—Waist 1 ^ Wrist ......... b) -./•*—l-lbow —Shoulder Time [s] 10 15 Fig. 2. Joints moving simultaneously, a) CTC step response b) corresponding torques With c1 =4, the step response for different values of k1 is illustrated in Fig. 3. Comparing the responses, it is clear that the waveform corresponding to k1 = 8 exhibits relatively less rise time and settling time. Thus, the optimum response is achieved with k1 = 8. 2 4 6 Time [s] 8 10 Fig. 3. VSC step response of the base joint for different values of k1 With tuned values of the C and K matrices, the step response of various joints moving simultaneously is shown in Fig. 4a. It can be inferred from the plot that the elbow joint and wrist joint, even after reaching their desired position, are not stable until the shoulder 0 joint reaches its final position and is stabilized. This is in accordance with the coupling effects of the joints. Looking at the applied torque plot shown in Fig. 4b, the switching effect is observable. - 0.8 0.6 0.4 0.2 10 ? 0 £ ao -10 a--20 < -30, a) —/-Wrist /—llbmv , it -Waist -Desired <— Shoulder )/ 2 4 6 Time [s] 10 Waist у Wrist S b) --------------- Elbow --^—Shoulder ; 2 4 6 8 Time [s] Fig. 4. Joints moving simultaneously a) VSC step response, b) corresponding torques 10 different gains (Ц) is then observed for base, shoulder, elbow and wrist joints with each joint moved at a time. Fig. 6 illustrates the step responses of the shoulder joint and elbow joint with 20° set as the desired angle. Fig. 5. AUTAREP - A custom-developed pseudo-industrial framework CTC simulation and hardware implementation both show that the increase in the gain-constant results in better performance by improving the system's response. In contrast to the simulation, the hardware results show that for same value of gain constant, each joint exhibit slightly different response. This is due to the fact that each joint's motor produces a different torque and speed. 3 EXPERIMENTAL SETUP AND RESULTS To validate the simulation results, experiments have been conducted on an indigenously developed AUTAREP illustrated in Fig. 5. It is a mini-industrial open-source, novel and complete robotic system that finds potential in training interns, imparting mechatronics concepts to engineering students and validating advanced algorithms for trajectory generation and control, object manipulation and grasping, path planning, etc. [17]. The mechanical system of the platform is built around a 6-DOF serial robotic manipulator. The arm's geometrical configuration resembles that of the human arm. Six precise DC servo motors actuate the robot while the sensory system is comprised of encoders and force sensing resistors (FSR) in addition to an on-board camera. The primary features of AUTAREP are mentioned in [15]. The designed electronic system mainly consists of an embedded controller DSPIC33F and 6A/50V rated custom-designed motor drivers. The hardware and software architectures of the platform are detailed in [17]. First, the CTC law is implemented on the manipulator. The trajectory tracking performance for 0.4 0.3 0.2 0.1 a) Desired S j*-f-Lambda 9 / H-Lambda 6 r 4 6 Time [s] 10 0.4- 0.3 0.2- 0.1 b) Desired n-Lambda n<—Lambda= =9 =6 1 / Time [s] 10 Fig. 6. CTC step response of the a) shoulder joint, b) elbow joint After simple trajectory tracking experiments, the coupling effect was also studied in order to determine the VSC law on the physical manipulator. Joints were moved simultaneously as well as independently to observe this effect. The results are illustrated in Fig. 7. It can be inferred from Fig. 7a that the shoulder joint exhibits overshoot for a longer period of time and settles to the reference angle after settling of the elbow joint. The coupling effect on the shoulder joint causes it to move faster. VSC simulation and hardware implementation both confirm a strong coupling effect. In contrast to the simulation, the elbow joint in the hardware implementation does not come to rest until the shoulder joint is stabilized at its destination. This is quite expected due to the fact that a link near the fixture or base is less exposed to mechanical stability. 0.4 '0.3 SP 0.2 0.1 Desired ГК-.....:............. a) H--У— Elbow li*—-f- Shoulder /'/ L—Wais 4 6 Time [s] 10 0.4 '0.3 SP 0.2 0.1 Desired , C, ..... b) / \h—A Vaist / if-Shoulder hit-Elbow // / 4 6 Time [s] 10 Fig. 7. Step response of base, shoulder and elbow joints when the joints are actuated: a) at the same time, b) one at a time 4 CONCLUSIONS This paper presents the design, simulation and hardware realization of CTC and VSC strategies. The simulation results have been verified through experimental implementation on a physical platform. Trajectory tracking results showed that the derived laws can effectively track the desired reference input for both non-linear control methods. The coupling effects present in the joints are less visible in the simulation but are more prominent in the hardware implementation. Future work will include a task dependent performance comparison of robust control strategies on multi-DOF robotic manipulators. 5 REFERENCES [1] World Robotics (Executive Summary). (2014). International Federation of Robotics, from: http://www.worldrobotics.org/ uploads/media/Executive_Summary_WR_2014_02.pdf, accesed at 2015-06-11, p. 11-24. [2] Caglar, C., Ali, K., Selcuk, M., Sadettin, K., Hakan, Y. (2014). An enhanced control technique for the elimination of residual vibrations in flexible-joint manipulators. Strojniški vestnik -Journal of Mechanical Engineering, vol. 60, no. 9, p. 592-599, D0I:10.5545/sv-jme.2014.1698. [3] Ajwad, S.A., Ullah M.I., Khelifa, B., Iqbal, J. (2014). A comprehensive state-of-the-art on control of industrial articulated robots. Journal of Balkan Tribological Association, vol. 20, no. 4, p. 499-521. [4] Siciliano, B., Khatib, O. (2008). Springer Handbook of Robotics. Springer Science & Business Media, DOI:10.1007/978-3-540-30301-5. [5] Ajwad, S.A., Iqbal, J., Ullah, M.I., Mehmood, A. (2015). A systematic review of current and emergent manipulator control approaches. Yang, S. (ed.) Frontiers of Mechanical Engineering, Higher Education Press, Beijing and SpringerVerlag Berlin Heidelberg, D0I:10.1007/s11465-015-0335-0. [6] Iqbal, J., Saad, M.R., Malik, A., Tahir, A.M. (2014). State estimation technique for a planetary robotic rover. Revista Facultad de Ingenieria, vol. 73, p. 58-68. [7] Iqbal, J., Tsagarakis, N.G., Caldwell, D.G. (2014). Human hand compatible underactuated exoskeleton robotic system. IET Electronic Letters, vol. 50, no. 7, p. 494-496, D0I:10.1049/ el.2014.0508. [8] Iqbal, J., Khan, A.H., Tsagarakis, N.G., Caldwell, D.G. (2014). A novel exoskeleton robotic system for hand rehabilitation - Conceptualization to prototyping, Biocybernetics and Biomedical Engineering, vol. 34, no. 2, p. 79-89, D0I:10.1016/j.bbe.2014.01.003. [9] Antonio, Y., Victor, S., Javier, M.V. (2011). Global asymptotic stability of the classical PID controller by considering saturation effects in industrial robots. International Journal of Advance Robotic Systems, vol. 8, p. 34-42, D0I:10.5772/45688. [10] Kadir, Z.A., Mazlan, S.A., Zamzuri, H., Hudha, K., Amer, N.H. (2015). Adaptive fuzzy-PI control for active front steering system of armored vehicle: Outer loop control design for firing on the move system. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 3, p. 187-195, D0I:10.5545/sv-jme.2014.2210. [11] Jingmei, Z., Haiyang, H., Bo, K. (2012). Studies of adaptive control methods based on VSC for trajectory tracking of robotic manipulators. Proceedings of International Conference on Robotics and Biomimetics, p. 429-434. [12] Tahir, F., Jaimoukha, I.M. (2013). Robust feedback model predictive control of constrained uncertain systems. Journal of Process Control, vol. 23, no. 2, p. 189-200, D0I:10.1016/j. jprocont.2012.08.003. [13] Tahir, F., Jaimoukha, I.M. (2013).Causal state-feedback parameterizations in robust model predictive control. Automatica, vol. 49, no. 9, p. 2675-2682, D0I:10.1016/j. automatica.2013.06.015. [14] Iqbal, J., Islam, R.U., Khan, A.H. (2012). Modeling and analysis of a 6 DOF robotic arm manipulator. Canadian Journal on Electrical and Electronics Engineering, vol. 3, no. 6, p. 300306. [15] Manzoor, S., Islam, R.U., Khalid, A., Samad, A., Iqbal, J. (2014). An open-source multi-DOF articulated robotic educational platform for autonomous object manipulation. Robotics and Computer Integrated Manufacturing, vol. 30, no. 3, p. 351362, DOI: 10.1016/j.rcim.2013.11.003. [16] Shafiei, S.E. (ed.) (2010). Sliding mode control of robot manipulators via intelligent approaches. Advanced Strategies for Robot Manipulators, Intech Open Access Publisher, p. 135172, D0I:10.5772/10193. [17] Iqbal, U., Samad, A., Nissa, Z., Iqbal, J. (2014). Embedded control system for AUTAREP - a novel autonomous articulated robotic educational platform. Tehnički vjesnik -Technical Gazette, vol. 21, no. 6, p. 1255-1261. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)7-8, 471-476 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2014.2180 Original Scientific Paper Numerical Investigation of the Effect of Process Parameters during Aluminium Wheel Flow-Forming So Young Hwang - Naksoo Kim - *Cheol-soo Lee Sogang University, Department of Mechanical Engineering, Korea Aluminium wheel flow-forming is an incremental forming technique that lengthens the material in the axial direction and thins it in the radial direction. There are many advantages to wheel flow-forming, including reducing the weight of the product, which leads to improved fuel economy. However, in terms of geometric accuracy, it is difficult to manufacture with high precision. In this study, process conditions, such as headstock load, feed rate, and radial compression, as well as the geometry of the rolls, are investigated to improve the precision of the product. To do this, a simulation was performed with SHAPE-RR, a type of FEM software. The result of the FE analysis was verified by comparing the bulging amount of the wheel. The material property of aluminium was determined through a hot compression test. As a result, geometric accuracy was improved by increasing the feed rate of the rolls. Moreover, applying a lower compression amount in the radial direction for the third roll increased the geometric precision. Keywords: aluminium wheel flow forming, backflow defect, geometric precision, compression amount Highlights • Process condition study for wheel flow forming is done using Finite Element method. • Validation of the simulation is done by comparing the bulged amount according to the headstock load. • Optimal radius tip and size of the rolls are revealed. • Effect of the compression amount for each roll is investigated. 0 INTRODUCTION Ecofriendly, lightweight, and high-efficiency are the most important key words in manufacturing fields, especially the automobile industries. There have been many studies on lightweight automotive design [1] to [4]. For this, reducing the weight of the wheels can be a solution. Both steel wheels and aluminium alloy wheels for passenger cars are used. The aluminium wheel's weight can be one-third less than that of a steel wheel. If the weight below the spring of the automobile is reduced by 1 kg, the total weight of the car can be reduced 10 to 15 kg [5]. This can lead to higher efficiency and more environmentally friendly automobiles. Aluminium wheels can be manufactured by casting or by casting followed by flow forming. Friction stir welding can also be used for the production of aluminium wheels, and many studies have been conducted on this technique [6] to [10]. If flow forming is included in the manufacturing process, the weight can be reduced by 10 to 18% because the mechanical characteristics are improved. During aluminium flow forming, three side rolls are used. The first roll starts to operate when the process begins and the second roll begins to operate a few seconds later, followed by the third roll. The three rolls operate at the same time during the process except at the beginning and the end. However, because of its complicated shape, there is no systematic design approach to determine the optimal process condition to increase geometric accuracy. Therefore, in the field, manufacturers rely on trial-and-error methods to find the best process condition. Many studies have been performed to understand the flow-forming process. Wong et al. reviewed the principles and developments of flow forming and insisted on the great potential of manufacturing complex shapes [11]. Yao and Murata experimentally investigated the process parameters, including feed rate, radial force, and surface roughness, for aluminium tube flow forming [12]. Molladavoudi and Djavanroodi experimentally investigated the effects of thickness reduction on the mechanical properties and spinning accuracy [13]. Lee and Lu demonstrated the relationship between the load of roller and mandrel to the strain rate and surface roughness. They also confirmed that it is beneficial to have larger rollers since the surface finishing is better than with smaller rollers [14]. Davidson et al. found that the depth of cut in turning operation is an important process parameter affecting the percentage of elongation, by using the Taguchi method [15]. While many experimental studies on flow forming have been done, there also have been studies using the finite element method (FEM) to investigate the flow forming process. Xu et al. used a 3D rigid *Corr. Author's Address: Department of Mechanical Engineering, Sogang University, 35 Baekbeom-ro, Mapo-gu, Seoul, 121-742 South Korea, nskim@sogang.ac.kr 471 plastic model to establish an FEM model of tube spinning, and the model was verified experimentally [16]. Xia et al. studied the effect of the main forming parameters, such as offset amount, feed rate, and path direction, on the spinning forces, theoretically and experimentally, using a finite element (FE) simulation [17]. Research exists on the defects that occur during spinning or flow forming, but it remains at the level of simulating and verifying the phenomenon by comparing experimental and analytical results. However, in the field, it is important to achieve the precise feature that is desired. Therefore, the purpose of this study is to understand the effect of the process parameters of headstock load, compression amount. feed rate, and roll shape on the geometric accuracy. To do this, a sensitivity analysis was performed using FE analysis to determine the effect of the process variables on the geometric accuracy. The SHAPE-RR program, which is exclusively used for the rotational manufacturing process, was used for FE analysis. To verify the validity of this analysis, the bulging amounts from the analytical and experiment results were compared, which was the verification method used by Mohebbi and Akbarzadeh to enhance the reliability of the FE analysis [18]. which affects the precision of the product. Therefore, the optimal headstock load must be studied. Axial force Fig. 1. Mechanism of backflow defect 2 WHEEL FLOW-FORMING MODEL To perform a sensitivity analysis of the wheel flow forming process, the SHAPE-RR program, which was developed based on the FEM, was used. Prior to the computation analysis, the material properties of aluminium had to be determined, and the validity of the analysis program had to be verified. 1 BACKFLOW DEFECT 2.1 Material Characteristics According to Cha et al., the so-called backflow defect occurs during the wheel flow-forming process [19]. Therefore, backflow must be reduced to enhance the geometric accuracy. The backflow defect is a bulging phenomenon located on the opposite side of the forming direction. As shown in Fig. 1, when the roller applies a load to a part where the diameter rapidly varies, the resultant load can be divided into radial and axial forces. The radial force reduces the thickness of the material, and the axial force elongates in the axial direction. When the roller is used to reduce the thickness, a material point is moved in the opposite direction of the roller movement to cancel out the axial force. Moreover, when the load of the headstock is enormous, backflow is generated, which creates the bulging defect, as shown in Fig. 1. According to Cha et al. [19], the amount of bulging is related to the headstock load. If the headstock load is excessive, the reaction force on the headstock is increased, which leads to a greater amount of bulging. However, if the headstock load is insufficient, the wheel will lift up from the mandrel, To describe the material flow during wheel flow forming, the Johnson-Cook flow stress model shown in Eq. (1) was used. In the equation, spl is the effective plastic strain, s pl the effective plastic strain rate, and Tmdt is the melting temperature. This is suitable for considering the strain rate and temperature terms because wheel flow forming is a bulk deformation process with time-dependency. Uf. = (a + B(spl )й )(1 + C ln tpl )(1 - T*m), (1) T * = T - T T - T melt ro (2) The Johnson-Cook model uses A, B, n, C and m as material coefficients. Friction and these coefficients were obtained through a hot compression test. The dimensions of the cylindrical specimen are a length of 40 mm and a diameter of 15 mm. To obtain the terms of strain rate and temperature, five experimental cases were set. The initial temperature of the aluminium specimen for the flow forming process was 360 °C. Therefore, the experimental conditions of temperature 4 6 8 Stroke [mm] Fig. 2. Comparison of experiment and simulation for 330 °C -10 mm/min 4 6 8 Stroke [mm] Fig. 3. Comparison of experiment and simulation for 360 °C -10 mm/min 4 6 8 Stroke [mm] Fig. 4. Comparison of experiment and simulation for 390 °C -10 mm/min were considered to be 390 °C and 330 °C. Each case was repeated three times. After obtaining the load-stroke curve from the experiment, it was averaged for each case. The coefficients and the friction that minimized the load-stroke curve between the experiment and the simulation were obtained. The coefficients are shown in Table 1, and the friction coefficient is 0.2. The load-stroke curves obtained from the experiment and the simulation are compared for each case in Figs. 2, 3, and 4. Table 1. Material coefficients of the Johnson-Cook model A B n C m T 1 melt 154 MPa 63 MPa 0.476 0.015 1.0 887K 2.2 Simulation Model The wheel flow forming simulation model is shown in Fig. 5. Three-dimensional eight-node solid elements were used in the model. The second and third rollers were placed 125° from the first roller, which put them 115° from each other. roller #3 Fig. 5. Full FE model of wheel flow forming a) b) Fig. 6. Numerical model for the cross-section of the wheel; a) actual shape, and b) simplified finite element mesh To increase the computational efficiency, a complicated shape in the corner of the work piece, shown in Fig. 6a, was removed. A mesh test was also performed to determine the maximum size of the mesh. The longitudinal mesh length on the side was studied, and with the same forming condition, the geometrical difference was compared, and 4 mm was chosen. The cross-section of the initial work piece and the boundary condition are shown in Fig. 6b. The boundary condition is given as 500 rpm. 2.3 Simulation Validation To validate the reliability of the simulation, the amount of bulging according to headstock load was compared with the experimental model. The obtained material properties were applied to the simulation. At each load condition, three experiments were conducted to ensure the accuracy of the results. The initial work piece and the flow-formed wheel are shown in Fig. 7. Fig. 7. Initial work piece and manufactured wheel 30000 60000 90000 Headstock load [N] Fig. 8. Comparison of experiment and simulation of bulging amount As shown in Fig. 8, the amount of bulging tends to increase as the load increases. This is because the load in the longitudinal direction is related to the bulging mechanism, as explained in Section 1. The amount of bulge converges at the end because of the geometric limitations. As shown in Fig. 8, the averaged bulging amount of the experimental model and the simulation is similar, with a 6"% difference between them. 3 GEOMETRIC ACCURACY Manufacturing the product with high geometric accuracy is the most important factor during the wheel flow-forming process. If it is accurate, the cutting process afterwards will be reduced. To perform a sensitivity analysis, simulation conditions have to be determined. The location of the rolls and the work piece during the process is shown in Fig. 5. Three side rolls were placed, as shown in Fig. 5. The gap between the mandrel and the work piece after the forming process was calculated as the cross-section area between them (Fig. 9). During wheel flow forming, the work piece is lifted up from the mandrel, which creates a gap. This gap has to be minimized to increase the geometric accuracy. First, it is important to investigate the effect of headstock load on the geometric accuracy. As mentioned in the prior chapter, it is important to determine the suitable headstock load. If the load is too high, the bulging amount will increase. However, if the headstock does not have enough force to settle the work piece down onto the mandrel, the work piece is lifted up from the mandrel. Therefore, according to the headstock load, the gap between the mandrel's top surface and the work piece is calculated. As a result, with a higher load, the gap is reduced (Fig. 10). However, a headstock load over 60,000 N creates bulging of over 0.5 mm, which is the limitation. Furthermore, the gap must be less than 400 mm2 for the wheel to be a non-defective product. If the gap exceeds 400 mm2, then the wheel cannot become used because of the geometry tolerance. Therefore, 30,000 N is the most appropriate headstock load. Second, the shapes of the side-rolls were investigated. Specifically, the diameter of the roll and the radius of the roll-tip is considered. The level of each parameter is shown in Table 2. In this table, r1 indicates the radius of the tip for the first roll, and r2 indicates the radius of the tip for the second roll. The tip radius of the third roll was not considered because it is the finishing roller. The upper bound was determined according to the limitations of geometry. To run an efficient amount for the experiment, the Taguchi method was used. Three variables with three levels fit the L9(34) orthogonal array. Therefore, a total of nine cases were simulated. For each case, the gap difference was calculated as shown in Table 3, and this was used to obtain the objective function as shown in Eq. (3): y = 1093 - 110.5 r1 - 34 r1 + 2 R + 3.3 r12 + 2.4 r22 . (3) By optimizing the objective function using the response surface method [20], the optimum values for each parameter to maximize the geometric accuracy were r1 = 20 mm, r2 = 12.7 mm, and R = 176 mm. The results indicate that the first roller must have the largest tip radius. The role of the first roll is to push down the material so that it can be flow-formed close to the mandrel. If the tip radius is large, it will be easier for the roller to push down the material, which makes it easier for the second and third rolls to shape the wheel. Fig. 9. Gap between mandrel and material 30000 60000 90000 Headstock load [N] Fig. 10. Gap according to headstock load In contrast, the compression amount of each roll was considered. As explained above, three side-rolls were used during the process, and the radial thickness that had to be reduced was 10 mm. Therefore, the amount of compression for each roll was adjustable, so the compression amounts of the rolls were investigated. Parameter d2 indicates the compression amount in the radial direction of the second roll and d3 represents this for the third roll. Automatically, the compression amount of the first roll is the rest of the 10 mm. Table 4 shows the results. As the compression amount is small for the third roll, it has a better geometric precision. This is because the third roll is the last roll, which affects the final shape of the work piece. Therefore, as much as the first roll compresses, it is easier for the second and third rolls to form precise geometries. Table 2. Roll-shape parameters and levels r1 [mm] r2 [mm] R [mm] Level 1 10 10 120 Level 2 15 15 180 Level 3 20 20 240 Table 3. Taguchi method matrix r1 [mm] r2 [mm] R [mm] Error Gap [mm2] 1 10 10 120 1 315 2 10 15 180 2 347 3 10 20 240 3 401 4 15 10 180 3 209 5 15 15 240 1 199 6 15 20 120 2 469 7 20 10 240 2 258 8 20 15 120 3 332 9 20 20 180 1 463 Table 4. Geometric accuracy according to compression amount d2 [mm] d3 [mm] Gap [mm2] Case 1 3 2 191.7 Case 2 4 2 238.5 Case 3 3 1 130.2 Table 5. Geometric accuracy according to feed rate Feed rate [mm/rev] Gap [mm2] 1 1.7 mm/rev 146.9 2 1.6 mm/rev 149.5 3 1.5 mm/rev 165.6 4 1.4 mm/rev 166.1 5 1.3 mm/rev 201.9 Finally, the feed rate was considered. The feed rate refers to the speed in the longitudinal direction while rotating. Therefore, 1.5 mm/rev means that the roller moves downward at 1.5 mm per revolution, and this is related to the forming process time. For each feed rate considered, the result is shown in Table 5. Case 1 is when the forming process was done in 14 seconds. As shown in Fig. 11, the gap, or the form error, is reduced as the feed rate increases. However, after 1.6 mm/rev, the gap between the mandrel and the work piece converges to 145 mm2. 4 CONCLUSIONS In this research study, the process conditions that can affect geometric precision were studied. The first factor that was considered was the headstock load, which is the main cause of bulging defects. As this is increased, the bulging amount increases. However, the lifted gap between the work piece and the mandrel decreases as the headstock loads increase. Considering that both factors are equally important, 30,000 N is the most proper headstock load for the presented case study. Furthermore, for the first roll, it is better to have a large tip radius. Due to geometric limitations, 20 mm is the largest tip radius. In contrast, the optimal tip radius of the second roll is 12.7 mm. In terms of compression amount, it is better for the first roll to compress more and to reduce the compression amount for the second and third rolls for geometric precision. The gap can be reduced down to 130.2 mm2. For the feed rate, which is related to process time, it is better for it to be fast. When the feed rate is over 1.6 mm/ rev, the gap is less than 150 mm2. These results only consider the geometric precision. Finally, this study can be extended to considering the lifespan of the rolls. 4 ACKNOWLEDGMENTS This work was supported by Sogang University (Grant no. 201210031) and the National Research Foundation of Korea (NRF) grant, funded by the Korean Government (NRF-2010-0023152). 5 REFERENCES [1] Belingardi, G., Beyene, A.T., Koricho, E.G., Martorana, B. (2015). Alternative lightweight materials and component manufacturing technologies for vehicle frontal bumper beam. Composite Structures, vol. 120, p. 483-495, D0l:10.1016/j. compstruct.2014.10.007. [2] Jiang, J., Wang, Y., Chen, G., Liu, J., Li, Y., Luo, S. (2012). Comparison of mechanical properties and microstructure of AZ91D alloy motorcycle wheels formed by die casting and double control forming. Materials and Design, vol. 40, p. 541549, DOI:10.1016/j.matdes.2012.04.029. [3] Kobold, D., Pepelnjak, T., Gantar, G., Kuzman, K. (2010). Analysis of deformation characteristics of magnesium AZ80 wrought alloy under hot conditions. Strojniški vestnik- - Journal of Mechanical Engineering, vol. 56, no. 12, p. 823-832. [4] Hirsch, J., Al-Samman, T. (2013). Superior light metals by texture engineering: Optimized aluminum and magnesium alloys for automotive applications. Acta Materialia, vol. 61, no. 3, p. 818-843, DOI:10.1016/j.actamat.2012.10.044. [5] Jambor, A., Beyer, M. (1997). New cars-new materials. Material & Design, vol. 18, pp. 203-209, DOI:10.1016/S0261-3069(97)00049-6. [6] Besharati-Givi, M.K., Asadi, P. (eds.) (2014) Advance in Friction Stir Welding and Processing. Woodhead Publishing, Cambridge. [7] Kallee, S., Nicholas, D. (1998). Application of friction stir welding to lightweight vehicles. SAE Technical Paper, 982362, DOI:10.4271/982362. [8] Podržaj, P., Jerman, B., Klobčar, D. (2015). Welding defects at friction stir welding. Metalurgija, vol. 54, no. 2, p. 387-389. [9] Klobčar, D., Kosec, L., Pepelnjak, T., Tušek, J. (2012). Microstructure and mechanical properties of friction stir welded AlMg4.5Mn alloy. Engineering Review, vol. 32, no. 2, p. 104-110. [10] Herakovic, N., Simic, M., Trdic, F., Skvarc, J. (2011). A machinevision system for automated quality control of welded rings. Machine Vision and Applications, vol. 22, no. 6, p. 967-981, DOI:10.1007/s00138-010-0293-9. [11] Wong, C.C., Dean, T.A., Lin, J. (2003). Review of spinning, shear forming and flow forming processes. Journal of Machine Tools and Manufacture, vol. 43, no. 14, p. 1419-1435, DOI:10.1016/S0890-6955(03)00172-X. [12] Yao, J., Murata, M. (2002). An experimental study on paraxial spinning of one tube end. Journal of Materials Processing Technology, vol. 128, no. 1-3, p. 324-329, DOI:10.1016/ S0924-0136(02)00473-9. [13] Molladavoudi, H.R., Djavanroodi, F. (2011). Experimental study of thickness reduction effects on mechanical properties and spinning accuracy of aluminum 7075-o during flow forming. International Journal of Manufacturing and Technology, vol. 52, no. 9-12, p. 949-957, DOI:10.1007/s00170-010-2782-4. [14] Lee, K.S., Lu, L. (2001). A study on the flow forming of cylindrical tubes. Journal of Materials Processing Technology, vol. 113, no. 1-3, p. 739-742, DOI:10.1016/S0924-0136(01)00585-4. [15] Davidson, M.J., Balasubramanian. K., Tagore, G.R.N. (2008). Experimental investigation on flow-forming of AA6061 alloy—A Taguchi approach. Journal of Materials Processing Technology, vol. 200, no. 1-3, p. 283-287, DOI:10.1016/j. jmatprotec.2007.09.026. [16] Xu, Y., Zhang, S.H., Li, P., Yang, K., Shan, D.B., Lu, Y. (2001). 3D rigid-plastic FEM numerical simulation on tube spinning. Journal of Materials Processing Technology, vol. 113, no. 1-3, p. 710-713, DOI:10.1016/S0924-0136(01)00644-6. [17] Xia, Q.X., Xie, Sh., W., Huo, Y.L., Ruan, F. (2008). Finite element simulation and experimental investigation on the forming forces of 3D non-axisymmetrical tubes spinning. Journal of Materials Processing Technology, vol. 206, no. 1-3, p. 500508, DOI:10.1016/j.jmatprotec.2007.12.066. [18] Mohebbi, M.S., Akbarzadeh, A. (2010). Experimental study and FEM analysis of redundant strains in flow forming of tubes. Journal of Materials Processing Technology, vol. 210, no. 2, p. 389-395, DOI:10.1016/j.jmatprotec.2009.09.028. [19] Cha, W.G., Hwang, S.Y., Kim, N., Lee, C.S. (2014). Analysis of mechanism of backflow defect of the aluminum wheel flow forming. Journal of Precision Engineering and Manufacturing, vol. 15, no. 6, p. 1075-1080, DOI:10.1007/s12541-014-0439-1. [20] Noordin, M.Y., Venkatesh, V.C., Sharif, S., Elting, S., Abdullah, A. (2004). Application of response surface methodology in describing the performance of coated carbide tools when turning AISI 1045 Steel. Journal of Materials Processing Technology, vol. 145, no. 1, p. 46-58, DOI:10.1016/S0924-0136(03)00861-6. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)7-8, 477-485 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2471 Original Scientific Paper Determination of Undissolved Air Content in Oil by Means of a Compression Method Adam Bureček* - Lumn Hružik - Martin Vašina VŠB - Technical University of Ostrava, Department of Hydromechanics and Hydraulic Equipment, Czech Republic This article describes a combination of experimental and mathematical methods for the determination of undissolved air content in hydraulic oil. The experimental part consists of the determination of the oil bulk modulus, considering the influence of undissolved air by means of a volume compression method in a steel pipe. A multiphase model of an oil/undissolved air mixture is subsequently defined using Matlab SimHydraulics software. The multiphase model permits the volume compression of oil and air bubbles independently of each other. Furthermore, time dependencies of pressures are mathematically simulated during the compression of the multiphase mixture of oil and undissolved air for different concentrations of the latter. The undissolved air content is determined by comparing the mathematically simulated and experimentally measured time dependencies of pressure increases. Keywords: oil/air mixture, bulk modulus, undissolved air content, hydraulic system Highlights • Experimental determination of secant bulk modulus and tangent bulk modulus of oil by means of a compression method. • Multiphase mathematical model of compressibility of oil/undissolved air mixture. • Mathematical simulation and measurement of time dependence of pressure during compression of oil/undissolved air mixture in steel pipe. • Determination of undissolved air content in oil by comparing the mathematical model with the measurement. • The undissolved air concentrations in the measured hydraulic system were determined in the range of 0.22 % to 0.49 %. 0 INTRODUCTION Basic properties of liquids are described by their density, viscosity and compressibility, which can be expressed by bulk modulus [1] and [2], resistance to deformation, or capacity. Liquid with a content of undissolved air is considered to be a mixture. The bulk modulus of the liquid/undissolved air mixture is significantly influenced by the concentration of undissolved air in the mixture. The mixture bulk modulus generally increases with increasing liquid pressure and decreasing temperature [3]. It is possible to determine the mixture's bulk modulus by different experimental methods, e.g. by acoustic [4] and [5], capacity [6], piezoelectric impedance [7] or volume [8] methods. Hydraulic oil, which is the most frequently used energy carrier in hydraulic systems, is the investigated liquid in this paper. The air content in hydraulic oil is typically in two states, i.e. in dissolved and undissolved states. In the dissolved (i.e. diffused) state, air in hydraulic oil is in the form of oxygen and nitrogen molecules that are mixed with oil molecules. The contents of other gases in air are negligible in comparison to the oxygen and nitrogen volumes. In the case of the undissolved state, the oxygen and nitrogen molecules are clumped together. For this reason, air bubbles are created. A volume of released and dissolved air in oil is given by Henry's law [9]. Oil with air bubbles creates an oil/undissolved air mixture. The mixture is characterized mainly by a higher compressibility, which corresponds to the relevant bulk modulus (sometimes referred to as effective bulk modulus). Therefore, the effective bulk modulus of an oil/air mixture includes the influence of undissolved air [4], [10] and [11]. The compressibility of the oil/ undissolved air mixture has a negative influence on the static and dynamic properties of hydraulic systems [12]. For this reason, it is necessary to eliminate the air content in this mixture. For a more accurate definition of a mathematical multiphase model of the mixture, it is necessary to define not only the oil bulk modulus, but also the undissolved air content, which is very difficult to measure. There are different methods for the experimental determination of the undissolved air content. It is possible to measure the bubble size distribution in liquid, e.g. by image analysis [13], drift flux analysis [14], as well as by acoustical [15], optical [16] and electro-resistivity [17] methods. However, a given method of the bubble size measurement is generally applicable for a certain bubble size range [18]. The aim of the paper is to describe a specially developed method for the determination of the undissolved air content in hydraulic oil on the basis *Corr. Author's Address: VŠB - Technical University of Ostrava, Department of Hydromechanics and Hydraulic Equipment, 17. listopadu 15/2172, Ostrava - Poruba 708 33, Czech Republic, adam.burecek@vsb.cz of a comparison of experimental measurements and mathematical modelling. 1 THEORETICAL BACKGROUNDS 1.1 Dissolved Air in Oil Air in the dissolved state presents a chemical bond of oxygen and nitrogen molecules to oil molecules. Petroleum oils will generally dissolve 8.5 % ± 0.5 % by volume of air at atmospheric pressure and room temperature [19]. For pressures higher than atmospheric levels, absorption follows Henry's law, which is defined as [20]: H = Coi C (1) where H is the dimensionless Henry's constant, Cair is the solute concentration in air, and CoU is the solute concentration in oil. The amount of dissolved air increases with increasing liquid pressure [21]. In the case of a disturbance from the equilibrium state as a consequence of pressure or temperature changes, air molecules are released and air bubbles are generated. For this reason, an oil/air mixture is created, or the air can be further dissolved in oil. This process is time-dependent. The time of the air release in oil is much shorter in comparison to the time of its dissolution. 1.2 Undissolved Air in Oil and Its Influence on Bulk Modulus of Oil/Air Mixture Dissolved air in oil is most frequently released to the undissolved state (i.e. to air bubbles) at pressure and temperature changes. Air bubbles can also enter the oil through different leaks. The oil/undissolved air mixture thus obtained has different properties in comparison to oil without air bubbles. It is possible to partly reduce the formation of bubbles at high operating pressures or by means of de-aeration devices, i.e. so-called gas separators. The bulk modulus, which is important in a mathematical model, is a significant property of the oil/undissolved air mixture. It is desirable to achieve high values of the bulk modulus in practice. The bulk modulus of the oil/undissolved air mixture is affected by many factors, e.g. by the pressure, temperature and volume of undissolved air. The amount of undissolved air has the greatest influence on the bulk modulus of the oil/ undissolved air mixture, mainly at low pressures. In this case, the air is much more compressible. The oil/undissolved air mixture can be defined as a multiphase mixture in mathematical models. The bulk modulus of the multiphase mixture is changed according to pressure and the amount of undissolved air. Then, the bulk modulus KPM of the multiphase oil/ air mixture is given by the equation [2], [22] and [23]: 1 + a kpm = ko ■ ( V1" Pa Pa + P 1 + a ■ KO 1 In Pa (2) ■(Pa + P ) (n+1)ln where KO is the bulk modulus of oil without air content, a = Va/VO is the relative air content in oil at atmospheric pressure, Va is the air volume at atmospheric pressure, VO is the oil volume at atmospheric pressure, n is the isentropic coefficient (n = 1.4), p is the working pressure, and pa is the atmospheric pressure. The above-mentioned equation allows the inclusion of the compression of the oil volume and undissolved air independently of each other. The oil bulk modulus KO is defined by the following equation [1] and [2]: K = Vo Ap AVo ' (3) where Ap is the pressure difference, and AVO is the difference of oil volume before and after compression. The bulk modulus of the oil/undissolved air mixture can also be taken into account in a mathematical model in a simplified manner. It is possible to assume a single-phase mixture, in which the compressibility of air bubbles is included in the bulk modulus KM of the oil/undissolved air mixture. In this case, the bulk modulus is defined by a constant for the given working pressure p. There is a certain inaccuracy in the mathematical model, especially at lower pressures. The thermodynamic effect, which is affected by the compression speed, occurs during the compression of the oil/undissolved air mixture. The isothermal effect proceeds at slow compression. In contrast, the isentropic effect is typical for rapid progressive compression. There are four different bulk modulus types of the oil/undissolved air mixture. From one standpoint, there are the secant bulk modulus and the tangent bulk modulus of the mixture [24] and [25] (see Fig. 1). Furthermore, each of these can be further divided into isothermal and isentropic moduli [26]. Fig. 1. Determination of secant (S) bulk modulus and tangent (T) bulk modulus of oil/undissolved air mixture The secant bulk modulus KM,S of the mixture is defined by the formula [10] and [25]: K = V M,S r M Ap (4) where VM is the volume of oil/undissolved air mixture, and AVM is the volume difference of an oil/undissolved air mixture before and after compression. The tangent bulk modulus KM,T of the mixture is expressed by the equation [10] and [25]: K = V ^M ,T r M d p dVT Fig. 2. Schematic diagram of experimental hydraulic circuit for determination of bulk modulus of oil/undissolved air mixture (5) Fig. 3. View of experimental hydraulic circuit for determination of bulk modulus of oil/undissolved air mixture 2 EXPERIMENTAL MEASUREMENT OF INVESTIGATED MIXTURE 2.1 Description of the Experimental Equipment The schematic diagram of the experimental equipment is shown in Fig. 2. The equipment consists of the hydraulic pump HP, the check valve CV, the relief valve RV, the steel pipe P, the seat valve SV, the reservoir R, the measuring equipment M 5050, the measuring point MP, and the pressure sensor PS. The M 5050 measuring equipment allows scanning, display and recording of measuring data from sensors that are used in hydraulics (e.g. from pressure, temperature and flow sensors). The hydraulic pump HP represents a flow source of the hydraulic system. If the seat valve SV at the pipe end (see Fig. 3) is open, hydraulic oil flows through the pipe P and the valve into the reservoir R. The seat valve is subsequently closed and, therefore, the flow through the valve is interrupted. Nevertheless, the pump HP supplies further liquid to the pipe P. Therefore, oil pressure is increased and a mixture of oil and air bubbles is compressed inside the steel pipe P. If the pressure of the mixture is increased to the value (i.e. p = 200 bar), which is adjusted by the relief valve, the relief valve RV is opened and subsequently the mixture of oil and air bubbles flows from the pump HP through the relief valve RV into the reservoir R. At the same time, the oil/undissolved air mixture in the pipe is compressed under the pressure that is adjusted by the relief valve RV. The pressure increase in the pipe was measured depending on the time. The pressure of the oil/undissolved air mixture inside the pipe P was recorded by the pressure sensor PS and the M 5050 measuring equipment. An example of the time dependence of the pressure p is shown in Fig. 4. The time interval At of the pressure scanning was set to 1 ms in this case. Measuring data were processed using Hydrowin software. The parameters of the steel pipe P are as follows: outside diameter DP = 0.03 m, inside diameter dP = 0.022 m, wall thickness sP = 0.004 m, length lP = 1.88 m, Young's modulus of elasticity EP = 2.1*10n Pa and Poisson ratio vP = 0.3. 250 200 150 100 50 0 S ce jO,