Strojniški vestnik Journal of Mechanical Engineering Strojniški vestnik - Journal of Mechanical Engineering (SV-JME) Aim and Scope The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis. The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue. The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s). Editor in Chief Vincenc Butala University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Technical Editor Pika Škraba University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Founding Editor Bojan Kraut University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerčeva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www.sv-jme.eu Print: Grafex, d.o.o., printed in 310 copies Founders and Publishers University of Ljubljana, Faculty of Mechanical Engineering, Slovenia University of Maribor, Faculty of Mechanical Engineering, Slovenia Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association President of Publishing Council Branko Širok University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Vice-President of Publishing Council Jože Balič University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: Image shows an interpenetrating phase composite with unique combination of the cellular ceramic material, with high strength and stiffness at low fractional densities, and the ductile metallic phase. The infiltration of the melt into pores was achieved by gravity casting. The interfacial reaction products were formed between the matrix and reinforcement, and affected the properties of the composite. Image Courtesy: Matej Steinacher, University of Maribor, Faculty of Mechanical Engineering, Slovenia ISSN 0039-2480 International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, University of Engineering and Technology, Pakistan Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, UL, Faculty of Mechanical Engineering, Slovenia Franci Čuš, UM, Faculty of Mechanical Engineering, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Igor Emri, UL, Faculty of Mechanical Engineering, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Janez Grum, UL, Faculty of Mechanical Engineering, Slovenia Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, UM, Faculty of Mechanical Engineering, Slovenia Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan Julius Kaplunov, Brunel University, West London, UK Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany Jernej Klemenc, UL, Faculty of Mechanical Engineering, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, UL, Faculty of Mechanical Engineering, Slovenia Gorazd Lojen, UM, Faculty of Mechanical Engineering, Slovenia Thomas Lübben, University of Bremen, Germany Janez Možina, UL, Faculty of Mechanical Engineering, Slovenia George K. Nikas, KADMOS Engineering, UK José L. Ocana, Technical University of Madrid, Spain Miroslav Plančak, University of Novi Sad, Serbia Vladimir Popović, University of Belgrade, Faculty of Mech. Eng., Serbia Franci Pušavec, UL, Faculty of Mechanical Engineering, Slovenia Bernd Sauer, University of Kaiserlautern, Germany Rudolph J. Scavuzzo, University of Akron, USA Arkady Voloshin, Lehigh University, Bethlehem, USA 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 62, (2016), number 2 Ljubljana, February 2016 ISSN 0039-2480 Published monthly Papers Matej Steinacher, Borut Žužek, Darja Jenko, Primož Mrvar, Franc Zupanič: Manufacturing and Properties of a Magnesium Interpenetrating Phase Composite 79 Junning Li, Wei Chen, Libo Zhang, Taofeng Wang: An Improved Quasi-Dynamic Analytical Method to Predict Skidding in Roller Bearings under Conditions of Extremely Light Loads and Whirling 86 Franc Cimerman, Matej Jarm, Branko Širok, Bogdan Blagojevič: Taking in Account Measuring Errors of Volume Conversion Devices in Measuring of the Volume of Natural Gas 95 Ali Majd, Ahmad Ahmadi, Alireza Keramat: Investigation of Non-Newtonian Fluid Effects during Transient Flows in a Pipeline 105 Sreten Simović, Vladimir Popović, Milanko Damjanović: Analysis of the Influence of Pavement Irregularities on the Lifespan of a Vehicle's Drive-Wheel Half Shaft 116 Yabing Cheng, Shuaibing Yin, Xiaopeng Wang, Lichi An, Huan Liu: Design and Analysis of Doubleside Meshing and Dual-phase Driving Timing Silent Chain System 127 César Garda-Hernandez, Rafael-Maria Gella-Marin, José-Luis Huertas-Talón, Nikolaos Efkolidis, Panagiotis Kyratsis: WEDM Manufacturing Method for Noncircular Gears, Using CAD/CAM Software 137 Strojniški vestnik - Journal of Mechanical Engineering 62(2016)2, 79-85 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2840 Original Scientific Paper Manufacturing and Properties of a Magnesium Interpenetrating Phase Composite Matej Steinacher1* - Borut Žužek2 - Darja Jenko2 - Primož Mrvar3 - Franc Zupanič1 1 University of Maribor, Faculty of Mechanical Engineering, Slovenia 2 Institute of Metals and Technology, Slovenia 3 University of Ljubljana, Faculty of Natural Sciences and Engineering, Slovenia The manufacturing and properties of the AE44 magnesium alloy reinforced with SiC-Al2O3-SiO2 ceramic foam were studied. The interpenetrating phase composite was manufactured by gravity casting at different preheating temperatures of the ceramic foam. The samples were investigated using optical and electron microscopy, energy dispersive X-ray spectroscopy, X-ray diffraction, transmission electron microscopy and compression testing. The interfacial reaction products (AlSiRE and AlMgSiRE) between the metal phase and ceramics were influenced by the preheating temperatures of ceramic foam and reduced the compression strength significantly Keywords: AE44 magnesium alloy, ceramic foam, interpenetrating phase composite, reaction product, mechanical properties Highlights • The interpenetrating phase composite was produced by gravity casting. • The AE44 Mg-alloy and SiC-Al2O3-SiO2 ceramic foam were used. • The AE44 Mg-alloy strongly reacted with the SiC-Al2O3-SiO2 ceramics. • The main interfacial reaction products were MgO, AlSiRE, and AlMgSiRE. • The hard AlSiRE and AlMgSiRE decreased mechanical properties. 0 INTRODUCTION Metal-matrix composites (MMCs) show improved performances over their matrix alloys. Magnesium matrix composites can offer potential applications within the automobile and aircraft industries. Interpenetrating phase composites (IPCs) usually display superior mechanical properties when compared to conventional MMCs reinforced with particles, intermetallic phases, ceramics, and carbon fibres. A unique combination of cellular ceramic materials, with high mechanical strength and stiffness at low fractional densities, and the ductility of the metallic phase may be considered as a major advantage of metal/ceramic IPCs [1]. They can be produced using various ways [2] and [3]. The infiltration can be achieved by a spontaneous capillary-driven metal infiltration [4], gas pressure-assisted infiltration [5] and [6], or squeeze-casting into a cellular ceramic preform [1], [7] and [8]. Interface behaviour between the matrix and the reinforcement can profoundly affect the properties of the MMCs [9]. The reinforcement type, alloying elements, solidification condition, and heat-treatment of MMCs can affect the local chemical composition and the extent of the interfacial reactions of the MMCs [10]. Several IPCs with Mg-matrix have been reported [1], [7] and [8]. The AZ91, and AZ31 magnesium alloys were used for the matrix, whilst the reinforcing phase was the ceramic foam based on SiO2, SiO2-Al2O3, and SiC-SiO2-C-Si. The mechanical properties (strength and Young's modulus) of the investigated IPCs were higher by up to 50 % when compared to the unreinforced magnesium alloy at room and elevated temperature. During this work, we used the AE44 magnesium alloy for the matrix and ceramic foam consisting of SiC, Al2O3, and SiO2 as the reinforcing phase. The IPC was produced using a simple and low-cost procedure, characterised and mechanically tested. It is to be expected that different manufacturing conditions can profoundly affect the properties of the resulting composite. 1 EXPERIMENTAL An AE44 magnesium alloy and SiC-Al2O3-SiO2 ceramic foam were used for manufacturing the IPC. The composition of the AE44 alloy was determined using ICP-AES (Table 1a). Polyurethane foam of the desired shape (Fig. 1a) was used as a preform for the manufacturing of the ceramic foam (Fig. 1b). Low viscous ceramic slurry (Table 1b) was infiltrated into a polyurethane preform. The excessive slurry was squeezed-out, and the coated preform was dried. Finally, it was heat-treated in order to remove the polyurethane skeleton and to sinter the ceramic powder. The ceramic foam is commercially accessible (ETI d.d.) and is used for filtration of the melt before the casting, because it is stable up to a temperature of 1700 °C. Fig. 1. The cellular material; a) a polyurethane preform and b) a ceramic foam Table 1. The chemical composition of the AE44 alloy and ceramics (wt. %) a) The AE44 magnesium alloy_ Al Mn Zn Si Ce 4.94 0.21 0.03 0.02 1.95 La Nd Pr Mg 1.71 0.48 0.28 Balance b) The ceramics SiC Al2O3 SiO2 Fe2O3 CaO 72.3 18.3 8.8 0.5 0.1 The IPC was made in several steps. A steel mould (Fig. 2) with a properly sized gating system and coated with a boron-nitride was made first. The gating system provided sufficient metal-static pressure that the melt could fill all pores of the ceramic foam. The ceramic foam had been inserted into the mould cavity before the mould was closed with the second part. Then the mould with inserted ceramic foam was preheated to a temperature of 500 °C, 600 °C, or 700 °C. During heating, the AE44 magnesium alloy was induction- melted and heated to a casting temperature of 730 °C. When the mould with inserted ceramic foam attained a suitable temperature, it was taken out of the furnace, placed on the vibration plate, and the melt was gravity cast through the inlet channel. After casting, the insulation cover was placed on the top of the mould, which acted as an insulated feeder. During the casting and solidification, the mould was vibrated so that the air from the pores of the ceramic foam and the mould cavity was removed. The metallic vibration plate transferred the heat from the bottom of the mould very well and thus enabled directional solidification. The AE44 magnesium alloy was also cast into the mould under the same conditions as ICP. Fig. 2. The mould with an inserted ceramic foam For the stress-strain measurements the samples with square cross-section (7 mm x 7 mm x 9.5 mm) were cut from the ICP and AE44 alloy. The ICP's samples contained 10.9 to 25.5 vol. % ceramics. The compression tests were carried out at room temperature (20 °C) and 200 °C by using Gleeble 1500D. Light microscopy (LM) work was done using the Olympus BX61 with the Analysis Materials Research Lab 5.0 software, and the scanning electron microscopy (SEM) in a FEI SIRION NC. The transmission electron microscopy (TEM) was carried out on a FEI Tecnai F20. A TEM specimen was cut out at a specific site using the focussed ion beam (FIB) in a FEI Nova 200. The hardness was determined using a Wilson Instruments Tukon 2100 B (load of 10 N) and micro-hardness using an Agilent Nano Indenter G200 testing machine (a Berkovich diamond indenter, depth limit of 500 nm, strain rate target of 0.05 cycles/s, harmonic displacement target of 2 nm, and a frequency target of 45 Hz). The composition of the ceramics was determined using an X-ray fluorescence (XRF) analyser Niton XL3t GOLDD+ (50 kVp). The X-ray diffraction (XRD) for the ceramics was carried out in a Philips 17-10 using Cu Ka radiation with a scan rate of 1.2 °/min, and for the alloy EA44 in a PANalytical B.V PW3830/40 using Cu Ka radiation with a scan rate of 0.25 °/min. 2 RESULTS AND DISCUSSION 2.1 The Magnesium Alloy Fig. 4. The fracture of the ceramics with sites of EDS analyses The characterisation of the commercial AE44 magnesium alloy using EDS-analysis showed that it consisted of a-Mg, AlnRE3, Al2RE, and Al10RE2Mn7. These phases were also identified using XRD-analysis. Each intermetallic compound contained all rare-earth elements; however, the content of Ce was the highest. Also other authors [11] to [13] found the same phases within microstructure of the AE44 alloy. 2.2 The Ceramic Foam Fig. 1b shows the structure of the SiC-Al2O3-SiO2 ceramic foam with the interconnected primary and mainly closed secondary porosity. The cellular-shape of the primary porosity had a mean-cell diameter of 4.23 mm, which was almost identical to the shapes and sizes of those pores in the polyurethane foam. The secondary porosity had a triangular void in the struts, which may reduce the strength significantly, is typical for reticulated foams. The mean-strut thickness was 0.55 mm. The XRD (Fig. 3) and EDS (Fig. 4, Table 2) results revealed that the ceramics was composed of four compounds: a-Al2O3, a-SiC, ß-SiC, and SiO2. Table 2. The EDS analyses of the ceramics (at. %, Fig. 4) site C O Al Si compound 1 54.6 2.5 / 42.9 SiC 2 5.2 63.6 27.1 4.1 Al2O3 3 7.3 60.7 0.5 31.5 SiO2 10 20 30 40 50 60 70 diffraction angle 20 / ' Fig. 3. The XRD pattern of the ceramics 2.3 Interpenetrating Phase Composite The reaction time between the melt and ceramic foam was equal to the solidification time, which was approximately 5 s, 30 s, and 60 s at the preheating temperatures of 500 °C, 600 °C, and 700 °C, respectively. During this rather short period, the melt completely filled the primary porosity, whilst the secondary porosity was filled only at preheating temperatures of 600 °C and 700 °C. The secondary porosity was probably filled by a combination of melt penetration through the strut walls and the melt infiltration through the holes in the strut walls, which represented direct links between the primary and secondary porosities. Fig. 5 shows that the melt did not just fill ceramic foam but also reacted with it. The infiltration at a preheating temperature of 500 °C resulted in the composite that suffered from partial debonding at the metal/ceramic interface and cold junctions (Fig. 5a). The widths of the cracks between the metal and ceramic skeleton were up to 50 ^m. In contrast, at the preheating temperature of 600 °C (Fig. 5b) and 700 °C (Fig. 5c) obtained a crack and debonding free interfaces but with the interfacial reaction products. These were as a result several reactions between the Mg-melt, SiO2, and Al2O3. The reaction products formed not only at the interfaces but also within the penetrated strut walls. The widths of the interfacial reaction products were more than 200 ^m. EDS analyses (Table 3) revealed the presence of MgO, and two intermetallic compounds with the general formulae AlSiRE and AlMgSiRE (Fig. 6), which were undetermined by XRD analysis because of their small amounts (Fig. 7). metal Fig. 5. The cross-section through a strut of the IPC manufactured at the different preheating temperatures; a) 500 °C, b) 600 °C, and c) 700 °C Table 3. The chemical compositions of the AlSiRE and AlMgSiRE phases as determined using EDS analyses (at. %) compound Al Mg Si Ce La Nd Pr AlSiRE 36.9 / 26.3 20.0 7.7 6.8 2.3 AlMgSiRE 14.1 24.2 33.1 16.6 6.6 3.9 1.5 AlMgSiRE AlSiRE AlMgSiRE 51 t « 0 Л M • (002 11 \ • ) 4(111 , •и лй 11) f (ООО 11 1 #(111) ) • (1111 w(1 t 1 ч # (002 è W 1 1 4 ) • • i 2 1/nm ss ф (012)# (002) • (012)# (ОТО)* (ООО) t (010) • (012)» (002)^ (012)» t k Ф • m 2 1/nm Fig. 6. The reaction products at the metal/ceramic interface; a) SEM, b) TEM micrograph, c) SAED AlSiRE, zone axis [110], d) SAED AlMgSiRE, zone axis [100] + a-Mg o a-AlA s a-SiČ 10 20 30 40 50 60 70 diffraction angle 20 /'0 Fig. 7. The XRD pattern of the interpenetrating phase composite Therefore, TEM investigations were carried out. The corresponding SAED-patterns for AlSiRE and AlMgSiRE are shown in Figs. 6c, d. Both are consistent with tetragonal structures. The AlSiRE phase always precipitated at first. It formed on the MgO film that covered the SiC and Al2O3. The AlSiRE particles later represented the nucleation sites for the AlMgSiRE phase. Thus, AlSiRE was usually partly or even completely surrounded by the AlMgSiRE, and the transformation from AlSiRE to AlMgSiRE took place with a reaction similar to a peritectic reaction. These compounds were presented in detail in reference [14]. The chemical reaction between Mg Fig. 8. Stress-strain diagrams of the IPC and AE44 alloy manufactured at different preheating temperatures; a) and b) 500 °C, c) and d) 600 °C, e) and f) 700 °C and tested at room temperature; a), c) and e) 20 °C and b), d) and f) 200 °C and SiO2, which was added as a binding agent at the manufacturing of the ceramic foam, caused the decomposition of the ceramics at the longest reaction time of 60 s (Fig. 5c). 2.4 Mechanical Properties Fig. 8 shows the stress-strain diagrams of the IPC and AE44 alloys. At both testing temperatures (20 °C and 200 °C) the IPC showed a significant improvement of the yield strength and the Young's modulus to AE4 alloy. But the rupture modulus of the IPCs, manufactured at preheating temperatures of 500 °C and 600 °C, were lower than for the non-reinforced AE44 alloy (Figs. 8a to d), while the IPC manufactured at 700 °C had higher values (Figs. 8e and f). Gibson and Ashby, crushing of the ceramic foams is according to the micro-mechanical model of based on bending of the struts [15]. Thus, in the interpenetrating metal/ceramic phase composites, the metal phase stabilies of the ceramic struts and partially prevents strut bending. This results in a major improvement of the mechanical properties [1]. However, the IPC, made at preheating temperatures of 500 and 600 °C, showed inferior values modulus of the rupture. At 500 °C, the load transfer from the metal matrix to the reinforcing ceramic skeleton was reduced due to air gaps between the infiltrated metal and the ceramic struts. Accordingly, no pronounced reinforcing effect was observed. The firm metal/ ceramic interfaces with interfacial reaction products were at 600 °C. The large interfacial reaction products, which are usually brittle and defective, decreased the compression strength of the IPC. Similarly, the tensile strength of the Ti-MMCs reinforced with SiC fibres decreased with the increase of the reaction zone but the strength was not reduced while the thickness of the interfacial reaction products was up to about 1 ^m [16]. At 700 °C, the ceramic foam disintegrated into the particles and the interfacial reaction products floated into the melt. These hard particles reinforced the matrix, similarly as in the conventional metal-matrix composite reinforced with particles because the mechanical properties were higher compared to the bulk AE44. 2.5 Hardness and Micro-hardness of AlSiRE and AlMgSiRE The measurements were carried out within the matrix and at the interface (Fig. 9a). The hardness of the matrix was 40 HV 1 and of the interface was 88 HV 1, respectively. The micro-hardness of the AlSiRE was 1336 ± 78 HV and of AlMgSiRE was 917 ± 32 HV (Fig. 9b). Thus, both phases were very hard, yet, AlSiRE was much harder than AlMgSiRE. Therefore the formations of these phases decreased in strength. On the other hand these can increase the resistance to cutting and wearing. Fig. 9. The measuring points for; a) hardness and b) micro-hardness 3 CONCLUSIONS According to the results of this work the following conclusions can be drawn: The interpenetrating phase composite can be manufactured by simple gravity casting. The AE44 magnesium alloy strongly reacted with the SiC-Al2O3-SiO2 ceramics. The preheating temperature of the ceramic foam effected the bonding and interfacial reaction. The main reaction products were MgO and the two novel phases AlSiRE and AlMgSiRE. The hard AlSiRE and AlMgSiRE decreased the strength of the interpenetrating phase composite at room and elevated temperatures. At a preheating temperature of 700 °C the ceramic foam disintegrated into the particles, which together with reaction products floated into the melt and reinforced the matrix, because the strength was increased. 4 ACKNOWLEDGEMENT This work was partly financed by the Slovenian Research Agency (ARRS), projects 1000-09-310152 and L2-2269. The authors also wish to thank Mrs. Petra Juvan, ETI, for preparing the ceramics and Mr. Tomaž Martinčič, University of Ljubljana, Faculty of Natural Sciences and Engineering, for manufacturing the composite. 5 REFERENCES [1] Zeschky, J., Lo, J., Höfner, T., Greil, P. (2005). Mg alloy infiltrated Si-O-C ceramic foams. Materials Science and Engineering: A, vol. 403, no. 1-2, p. 215-221, D0l:10.1016/j. msea.2005.04.052. [2] Claussen, N., Urquhart, A.W. (1990). Directed Oxidation of Molten Metals. Pergamon Press, Oxford. [3] Liu, W., Köster, U. (1996). Criteria for formation of interpenetrating oxide/metal-composites by immersing sacrificial oxide preforms in molten metals. Scripta Materialia, vol. 35, no. 1, p. 35-40, D0I:10.1016/1359-6462(96)00084-X. [4] Toy, C., Scott, W.D. (1990). Ceramic-metal composite produced by melt infiltration. Journal of the American Ceramic Society, vol. 73, no. 1, p. 97-101, DOI:10.1111/j.1151-2916.1990. tb05097.x. [5] Prielipp, H., Knechtel, M., Claussen, N., Streiffer, S.K., Müllejans, H., Rühle, M. Rödela, J. (1995). Strength and fracture toughness of aluminum/alumina composites with interpenetrating networks. Materials Science and Engineering: A, vol. 197, no. 1, p. 19-30, DOI:10.1016/0921-5093(94)09771-2. [6] Skirl, S., Hoffman, M., Bowman, K., Wiederhorn, S., Rödel, J. (1998). Thermal expansion behavior and macrostrain of Al2O3/Al composites with interpenetrating networks. Acta Materialia, vol. 46, no. 7, p. 2493-2499, DOI:10.1016/S1359-6454(98)80033-5. [7] Zeschky, J., Goetz-Neunhoeffer, F., Neubauer, J., Lo, S.H.J., Kummer, B., Scheffler, M., Greil P. (2003). Preceramic polymer derived cellular ceramics. Composites Science and Technology, vol. 63, no. 16, p. 2361-2370, DOI:10.1016/ S0266-3538(03)00269-0. [8] Zeschky, J., Lo, S.H.J., Scheffler, M., Hoeppel, H.W., Arnold, M., Greil, P. (2002). Polysiloxane-derived ceramic foam for the reinforcement of Mg alloy. Zeitschrift für Metallkunde, vol. 93, no. 8, p. 812-818, DOI:10.3139/146.020812. [9] Kainer, K.U. (2006). Basics of Metal Matrix Composites. Kainer, K.U. (ed.). Metal Matrix Composites: Custom-made Materials for Automotive and Aerospace Engineering. Wiley-VCH & KGaA, Weinheim, p. 1-52, DOI:10.1002/3527608117. [10] Mordike, B.L., Lukač, P. (2001). Interfaces in magnesium-based composites. Surface and Interface Analysis, vol. 31, no. 7, p. 682-691, DOI:10.1002/sia.1094. [11] Rzychon, T., Kietbus, A. (2008). Microstructure and tensile properties of sand cast and die cast AE44 magnesium alloy. Archives of Metallurgy and Materials, vol. 53, no. 3, p. 901907. [12] Zhu, S.M., Nie, J.F., Gibson, M.A., Easton, M.A., Bakke, P. (2012). Microstructure and creep behavior of high-pressure die-cast magnesium alloy AE44. Metallurgical and Materials Transactions A, vol. 43, p. 4137-4144, DOI:10.1007/s11661-012-1247-9. [13] Rzychon, T., Kietbus, A., Cwajna, J., Mizera, J. (2009). Microstructural stability and creep properties of die casting Mg-4Al-4RE magnesium alloy. Materials Characterization, vol. 60, no. 10, p. 1107-1113, DOI:10.1016/j.matchar.2009.05.014. [14] Steinacher, M., Mrvar, P., Zupanič, F. (2015). Interaction between AE44 magnesium alloy and SiC-Al2O3-SiO2 ceramic foam. Transactions of Nonferrous Metals Society of China, vol. 25, no. 3, p. 1011-1019, DOI:10.1016/S1003-6326(15)63692-5. [15] Gibson, L.J., Ashby, M.F. (1997). Cellular Solids: Structure and Properties, 2nd ed. Cambridge University Press, Cambridge, DOI:10.1017/cbo9781139878326. [16] Onzawa, T., Suzumura, A., Kim, J.H. (1991). Influence of reaction zone thickness on tensile strength for titanium matrix composites reinforced with SiC fiber. Tsai, S.W., Springer G.S. (eds.) Composites: Design, Manufacture and Application, Proceedings of the 8th International Conference on Composite Materials, Honolulu, p. 19J/1-19J/10. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)2, 86-94 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2848 Original Scientific Paper An Improved Quasi-Dynamic Analytical Method to Predict Skidding in Roller Bearings under Conditions of Extremely Light Loads and Whirling Junning Li12* - Wei Chen2 - Libo Zhang2 - Taofeng Wang2 1 Xi'an Technological University, School of Mechatronic Engineering, China 2 Xi'an Jiaotong University, Key Laboratory of Education Ministry for Modern Design and Rotor-Bearing System, China Slip often occurs in high-speed and light-load roller bearings (HSLLRBs) when the frictional drive force is inadequate to overcome the drag forces between the rolling elements and the raceway. Formerly, skidding analysis of HSLLRBs considering bearing whirling based on a simplified method using the Dowson-Higginson empirical model, although the analytical results of the cage slip fraction show significant discrepancies with the experimental data, as for extremely light radial loads. One of the main reasons is the inaccuracy of the evaluation of oil drag forces using the empirical equations. In this study, the elastohydrodynamic lubrication (EHL) method was adopted to calculate the oil film thickness and pressure distribution of HSLLRBs, so as to obtain more accurate oil drag forces. The cage speed and cage slip fraction were obtained by combining the whirl orbits, drag forces, load, kinematic equations and other related equations and then solved using the Newton-Raphson method. The skidding mechanism was investigated in terms of various operating parameters such as whirl orbit radii, radial load and viscosity. The results showed that the cage slip fraction and cage speed oscillate over time because of the whirl, which leads to an increase in the risk of bearing skidding damage. Under the extremely light load and high speed, the slip and influence of the whirl on bearing skidding increases as the whirl radius and radial load increases, while the viscosity shows a reverse trend. Therefore, in order to reduce slip and skidding damage, the whirl radius and radial load should be decreased suitably, while the viscosity should be increased moderately. A comparison between the calculated and experimental results shows that the proposed method is both feasible and valid. Keywords: skid, whirl, roller bearing, squeeze film damper, EHL, operating parameters Highlights • Extremely light load and whirling are considered in the skidding analysis of roller bearings. • An improved quasi-dynamic skidding analysis method coupled with EHL is proposed. • The proposed method proved to be feasible and useful for predicting skidding in HSLLRBs. • The influence of the whirl orbit, radial load and viscosity on HSLLRBs skidding are analyzed. • Suggestions for preventing skidding are summarized. 0 INTRODUCTION Rolling-element bearings are key precision components used for rotor support in nearly all machinery. In high-speed and light-load roller bearings (HSLLRBs), such as the typical main shaft bearings of an aircraft engine, the centrifugal forces on the rollers are considered to play a major role in its mechanics. Under these conditions, the tractive forces between the inner raceway and the rollers are frequently insufficient to overcome the drag on the rolling element assembly, which results in the phenomenon of slip [1]. Extreme slip between the roller and the raceway can cause wear on the rolling contact surfaces and subsequently result in a smearing type of surface damage [2]. Several analytical models have been developed for the prediction of slip in roller bearings under different conditions [1] to [7]. Dowson and Higginson analyzed the effects of film thickness and frictional forces on cage slip and derived equations for calculating the various forces acting on a roller under rigid and elastohydrodynamic lubrication (EHL) regimes [3], but did not include the effect of centrifugal forces. Harris [4] and Harris and Kotzalas [5] proposed an analytical method for predicting skid in high-speed roller bearings; this method yielded good results under conditions of medium and excess loading. However, analytical data on the cage speed show significant discrepancies with the experimental data at extremely light radial loads. Poplawski developed a roller bearing model to estimate the cage slip and cage forces, whose analytical predictions of slip correlated well with his test results [6], however they were not in agreement with Harris's cage speed predictions and experimental results for extremely light radial loads. Since the slip velocities at the rolling element-cage contact are typically large, a constant friction coefficient is used at the rolling element-cage pocket contact [7]. Tu et al. [8] and Shao et al. [9] presented an analytical model to investigate skidding during acceleration of a rolling element bearing, which takes into consideration the contact force and friction force between the rolling elements races and the cage, as well as the gravity and centrifugal force of the rolling elements, and analyzed the effects of a localized surface defect on the vibration response of the cylindrical roller bearing. In addition, based on the computed tomography, Zbrowski and Matecki analyzed the grinding smudges and subsurface defects in roller bearing rings [10]. Zhang [11] and Li et al. [12] analyzed various factors that influence slip and developed skidding analysis software for high-speed roller bearings. Most of the above mentioned theoretical analyses assume that roller bearings are installed directly in the bearing house. In practice, squeeze-film damper bearings (SFDBs) are widely used to inhibit the effects of vibration in a rotor-bearing system [13] and [14]. The inner ring of SFDBs is often mounted on the outer ring of the roller bearing with an interference fit that makes all of them whirl together, and thus inevitably influencing skidding damage of a roller bearing. Severe skidding of HSLLRBs in a whirling squeeze film damper often leads to vibration and failure of rotor-bearing systems or even entire machines. Many studies have focused on HSLLRBs skidding analysis without considering bearing whirl. Based on a simplified method using the Dowson-Higginson empirical model, Li and Chen analyzed the effects of different structure parameters on skidding of high-speed roller bearings in SFDBs considering bearing whirling [12], in which the radial load exceeds 500 N. As for the extremely light radial loads, the analytical data on slip have significant discrepancies with the experimental data. One of the main reasons is the inaccuracy of the evaluation of oil drag forces with the Dowson-Higginson empirical equations, which have a significant effect on the bearing skidding analysis, especially under conditions of extremely light radial loads. In this manuscript, HSLLRBs in a whirling squeeze film damper were taken as an example for skidding analysis. In order to obtain more accurate oil drag forces to assist in the skidding analysis of HSLLRBs, the EHL equations were solved using multigrid methods to obtain the values of the oil film thickness and pressure distribution, and then to acquire the fluid frictional drag forces. The skidding mechanism was investigated systematically in terms of various operating parameters such as whirl orbit radii, radial load and lubricating oil viscosity. 1 SQUEEZE-FILM DAMPER BEARINGS Squeeze-film damper bearings (SFDBs) are widely used to inhibit the effects of vibration in a rotor-bearing system. It has been shown that correctly designed SFDBs are a very effective means for reducing both the amplitude of rotor motion and the force transmitted to the bearing support [13] and [14]. In general, the rotating shaft carries a roller bearing, whose outer ring whirls with the inner ring of the SFDBs in the oil-filled clearance space between the inner and outer rings of the SFDBs. The outer ring of the bearing forms a damper so that rotation is inhibited by an anti-rotation pin [15] and [16]. Therefore, the damping effect of SFDBs mainly relies on an inner ring whirl, which squeezes the oil film in the clearance and forms resistance from the film. The inner ring of SFDBs is often mounted on the outer ring of roller bearing with an interference fit that makes all of them whirl together, thus inevitably influencing skidding damage to the roller bearing. The HSLLRBs of a whirling squeeze film damper is taken as an example for skidding analysis. The contact between roller and raceway is viewed as a rigid contact, the roller purely rolls along the outer ring raceway for the unloaded zone, and the cage normal forces are equal for every roller over all contact areas of the roller bearing. 2 MATHEMATICAL MODELS 2.1 Whirl Orbits Bearing whirl is majorly stimulated by shaft rotation. Therefore, the whirling frequency is relevant to shaft speed. In this manuscript, the authors considered the example of a roller and assumed its whirl orbits to be a circle with radius e with simple harmonic vibration in the X and y directions [17] and [18]. Taking the time of the maximum amplitude as the initial moment, the coordinates of orbit center Ow can be expressed as follows: So, er = e cosa t, ey = esmrnjt. er = -ea cosa,t, e,, = -ею sin m t. (1) (2) (3) (4) 2.2 Fluid Frictional Drag Forces 2.3 EHL Formulas As for the extremely light radial loads, the cage slip shows significant discrepancies with the experimental results. One of the main reasons is the inaccuracy of the evaluation of oil drag forces using the Dowson-Higginson empirical equations, which have a significant effect on the bearing skidding analysis, especially under conditions of extremely light radial loads [19]. The oil drag forces relate to the film thickness and pressure. The fluid frictional drag forces model of the HSLLRBs is shown in Fig. 1. Fig. 1. Model of the fluid frictional drag forces The fluid frictional drag forces T can be expressed as follows [20]: T = -f°h dx + (u2 - ul) — dx. (5) Jx,. 2 dx Jx- h h Here, x and xo are the coordinates of the inlet and outlet. Eq. (5) can be expressed in terms of dimensionless quantities as follows: F 4 x Hb dp ,v rx y.Rfi dx + fx JX X 8leR dx JX Hleb -dX, F = - Xo Hb dp ,v , rx voRon dX + fx x -dX. (6) X 8leR3 dX JX Holeb The solution to the fluid frictional drag forces equations are related to the oil film thickness and pressure distribution, which are the main challenge and the emphasis of this study. Here, the EHL method is adopted to calculate the oil film thickness and pressure distribution of the roller bearing, and then the fluid frictional drag forces can be obtained using Eq. (6) to replace the inaccuracy values obtained using the Dowson-Higginson empirical equations. As for the isothermal EHL in the line contact, the basic equations and their dimensionless forms can be expressed as follows: Reynolds equation: ( dp) d (pH ) dX Г dX )- dX = 0, (7) where U 4W2 ' The Reynolds boundary conditions of the equation are: in the area of the oil inlet X=Xin, P = 0; and in the area of oil outlet X=Xout , P=dP / dX= 0. (1) Film thickness equation: __\_ 2 2k h = H0 + ^— ^ J °и P(S) ln(X - S)2 dS. (8) (2) Viscosity equation varying with pressure (Roelands equation): П = exp{(lnn0 + 9.67)[(1 + 5.1x10-9pHP)z -1]}. (9) (3) Density equation: P = 1 + 0.6 xl0-9 pHP Л (10) l +1.7 x10 pHP (4) Load balancing condition equation: f x-'pdX = П. JX<„ 2 Eqs. (7) to (11) can be expressed in terms of x о dimensionless quantities as follows: X = —, p=—, b Po (11) H = R hR П = По P = P, PH и - Пои E ' R W = ■ w E ' R Ph Here, the EHL equations are solved using multigrid methods to obtain the values of the oil film thickness and pressure distribution, and then to acquire the fluid frictional drag forces using Eq. (6) to assist in the skidding analysis [19] and [20]. 2.4 Kinematics and Mechanical Models of the Roller The roller whirls along with the bearing. The kinematics model of a roller is shown in Fig. 2, and 6 П 2 b the forces acting on the roller at angular location qj are shown in Fig. 3 [12]. Fig. 2. Kinematics model Consider a separate roller as an analysis object; its kinematics model is shown in Fig. 2. Thus, the sliding velocities can be expressed as follows: Vj= \(dm - Db )(Ю-Юс ) - i Dbvmj, (12) Voj= 1(dm + Db )Юс - iDbvmj. (13) Additionally, the fluid entrainment velocities are defined as follows: Uj= 2 (dm - Db )(Ю-Юс ) + i Db^mj-, (14) Uoj= \(dm + Db )g>c + 2 Dbo>a]. (15) The sliding and the fluid entrainment velocities can be expressed in terms of dimensionless quantities as follows: V = Vn^L j ER : v = jo oi er U. = " ER U = Uo^L o ER ' (16) (17) (18) (19) In the loaded zone, the dynamic balance equations are expressed as follows: Qyij + Fa- Qyoj) = mr e, (20) Fj+ Qxj - Qxoi- Foj- Fdj = mr p, (21) Fj+ Foj- fj= 0. (22) In the unloaded zone, the dynamic balance equations are expressed as follows: F — O = m e , a z-syoju r y > F — Q — F = me dju z~-xoju oju r x Foju + fdj = 0. (23) (24) (25) Eq. (20) to (25) can be expressed in terms of dimensionless quantities. In the loaded zone, - R - - mè Q - + (-o-)(F - Q - ) = —-y- R Ш z-'yoj' IE ' R i ' \g/\zZxoj oj dj> iE IR Fi'+ ( R ) Fl- fdj= 0 In the unloaded zone, (26) Fj + Qxj - (biQxoj + F + Fdj ) = Ä, (27) (28) and (4) with the kinematic equations, load equations, Dowson-Higginson formulae, and other related equations and solved using the Newton-Raphson method [4] and [12]. 3 ANALYSIS AND DISCUSSION A single-row cylinder roller bearing is taken as an example. The default parameters are listed in Table 1. The effects of various operating parameters on the skidding of HSLLRBs taking into consideration bearing whirl are investigated. F -Q = -^y- ш Üyoju IE IR F - q - F = dju z~-xoju oju IE ' R„ Fj + fdj = 0. For entire rollers, Ы= 0. j=i (29) (30) (31) (32) Table 1. Details of the bearing and lubrication Qyij and Qyoj denote normal forces transmitted by the raceways to the roller owing to bearing radial load and geometry. The new fluid frictional drag forces equations are adopted here: f.=-f X°Hb dPdx +fxVRn 1 Jx, 8leR dX j- F = - X 8leRi dX Jx, H,leb X Hb dP,v , rX V0R0n 1-х dX + [ Jx X 8leR3 dX Jxi Hleb dX, (33) dX, (34) Qxj = 18.4(1 - D)G(35) In addition, the cage drag normal force and tangential force acting on a roller are expressed as follows: _ Л 2.447П/ Vj ■ (Dbx /2)2 , (37) Fdj= 2J-^-dx, (37) fi 0.5482ц/ a, ■ (Dbx /2)2 , (38) fdj= 2J-/=-dx (38) yfi0 In this study, the EHL method is adopted to calculate the oil film thickness and pressure distribution of the roller bearing, which allows the fluid frictional drag forces to be determined. Eqs. (27) and (30) can be further substituted into Eq. (32), and then Eqs. (28), (31) and (32) form an equation set with Z+1 equations. Then, the cage speed and cage slip fraction can be obtained by combining Eqs. (3) Type d [m] dm [m] G E' [Pa] Pr [kg/m3] NU214 0.07 0.0975 5000 2.28 e11 7800 Cdyn db [m] l [m] Z Pd [m] П0 [Pas] 1.37 e5 0.013 0.013 18 2.5 e-5 0.08 The cage slip fraction Sf can be expressed as follows: Sf = 1--, nm nš = - n. (1 - D ). 2 A dj (39) (40) 3.1 Effect of Different Whirl Radii on Skidding Bearing whirl induces the generation of additional forces in roller-cage and roller-rings pairs. More importantly, the oil film forces magnitudes of the roller-rings and roller-cage to oscillate as well. These forces oscillate over time, which often effect fatigue life and skidding damage to the bearing. Fig. 4 shows that the cage slip fraction and cage speed vary with time as the result of the whirl. In addition, the amplitudes of the cage speed and cage slip fraction increase with an increase in the whirl radius. So the degree of skidding can be reduced by suitably decreasing the whirl radius. 3.2 Effect of Different Radial Loads on Skidding In contrast to the previous results, Fig. 5 shows that the cage slip fraction increases as the radial load increases, but that the cage speed shows a reverse trend. On the other hand, Table 2 shows that the amplitudes of cage speed and cage slip fraction decrease with an increase in the radial load. In other words, the influence of whirl on bearing skidding increases as the radial load increases. Therefore, the degree of skidding can be и 0 a) b) Fig. 4. Skidding analysis under different whirl radii (ni = 5000 r/min; Fr = 200 N; Cdyn/P=675); a) cage speed vs. whirl radius, and b) cage slip fraction vs. whirl radius к T3 1300 128012601240 12201200 - Fr= 100 N Fr= 150 N Fr= 200 N 0.00 0.01 a) 0.02 0.03 Time t [s] 0.04 0.05 45 [ 44 Й О •■ö 43 И 42 . & И 41 i i=1 where E0, Ei and Ej are constant coefficients of the regression estimator, zt is an independent variable and y is a dependent variable.. The coefficients are determined using the method of least squares. For the system's response we look at the function of approximation ferror, for which we demand that it has a minimum: n n 2 ferror =ZZU - y(E0, Ei, Eij; j >i > Eii Zi > Z j>i )) . (13) i=1 j=1 The function of the approximation error will have a minimum if the first derivative of the function ferror with respect to the parameters Ei is equal to 0, but the second derivative of the function is positive. The acquired system of equations Eq. (14) is usually solved via the Gaussian method. dferrm- = 0 = 2 J ( - y(Eo, Ei, zi, E2, z2,...) )), 5E„ i,j=1 = 0 = 2 J (- y( ^ ^ z, ^ Z2,..) )), df™. = 0 = 2 J y( E0, El, Zi, E2, z2,...) ), i,j=1 dE2 dE2 dfer - = 0 = 2 J(y, -y(E0,Ei,zi,E2,z2,...))-E-.(14) dE 2.1 Determining the Equations of the Regression Model with the Help of Experimental Design Theory The equation of the regression models can be determined on the basis of understanding the mathematical model that describes a specific physical process. However, the physical process is not always sufficiently well known, and in such an example we may help ourselves with the experimental design theory from [13] and the graphics. Within the experimental design theory there is also the perfect factor experiment, which is used when we are interested in the influence of two or more factors (k) on at least two levels (r), with several repetitions of the experiments under unchanged conditions (h). The number of required experiments is: к i i = r ■ h. (15) If we assume that the effect of the factors is quadratic, then we need 32 experimental points. The geometrical interpretation is shown in Fig. 5, where the independent variables are marked as zj and z2: У = y(ZP Z2 ). (16) In Fig. 5, the experimental points for which we carry out the measurements are presented. The plan of the experiment 32 also has a central point. On this basis the corresponding regression polynomial that we determine with the method of least squares is equal to: У = E0 + E1Z1 + E2 Z2 + E 12 Z1Z2 + E 11 z,2 + E22 z22. (17) T3=(-U) T8=ro,i) T4=(U) Г1 at \ ' r\ T5=(-1,0) T9H(0,0) T7=(1,0) i,j=i dE T,=(-1,-1) T6=|C0,-1) T2=(1,-1) Fig. 5. The geometric interpretation of the quadratic perfect factorial experiment 32 2.2 The use of regression analysis for estimating the VCD measurement error The functional dependency of the measurement error is defined with Eq. (11) fC = fC (t, p). With a consideration of Eq. (17), we can move this dependency with a regression polynomial, Eq.(18): fC = E0 + Ej + E2p + El2t ■ p + Ellt2 + E22p2. (18) The main requirement of a regression analysis is that the number of available measurements is as numerous as possible. For the testing of type 1 VCDs, Fig. 4, N= 15 measurement points were measured, thus we must determine M = 6 parameters in the approximation equation. The measure for the estimate of suitability for the approximation is also the standard error of the estimate SEE, which takes into account the function of errors and the number of degrees of freedom v=N—M. It is defined with the expression in Eq. (19): error of the estimate SEE, which is measured in units dependent on the variable fC and is very small. SEE =i ferror N - M (19) The easiest way to determine the coefficients in Eq. (18) is with the use of a regression analysis, where we consider the use of variance analysis ANOVA. The F statistic in Table 3 determines whether the variation between the sample means is significant on the level 0.05. An example of such a printout from a regression analysis on the basis of the measurements from Fig. 4 and the use of computer programs SPSS [17] or Excell are shown in Table 3. Table 3. Regression analysis and ANOVA with the use of computer programs SUMMARY OUTPUT Regression Statistics Multiple r 0.98775 r2 0.97565 Adjusted r2 Square 0.96212 SEE 0.01581 Observations 15 ANOVA v SS MS F Significance F Regression 5 0.09008 0.01802 72.116 5.50981E-07 Residual 9 0.00225 0.00025 Total 14 0.09233 Coefficients Standard Error E00 0.604690 0.023 E1 -0.091278 0.015 E2 -0.002149 0.001 E12 0.000600 0.000 E11 0.009292 0.002 E22 0.000022 0.000 The results of the regression analysis are shown in Fig. 6. On the basis of the variance analysis, Table 3 and Fig. 6, we can conclude that the regression curve fits well with the measurements. Therefore, the significance F is very small, the coefficient of multi regression r2 is very large, as is the standard 2.000 3.000 4.000 Absolute pressure p [bar] Fig. 6. A comparison between the measured values and the approximated results in accordance with Eq. (18) at parameters r2 = 0.976 and SEE=0.016 % For the data from the regression analysis that are collected in Table 3 and Fig. 6, we can determine the next statistical parameters that estimate the suitability of the approximations chosen in Table 4. Table 4. Regression analysis and ANOVA Statistic Value SEE 0.0158 % r2 0.97565 Significance 5.51 • 10-7 3 EXPERIMENTAL ANALYSIS OF THE USE OF THE MULTI-REGRESSION MODEL FOR MEASUREMENT ERRORS DURING A VCD TEST For the use of a regression analysis at natural gas VDC we chose eight different models of the VDC from different manufacturers, i.e., M1 to M8. The VDCs M6, M7, and M8 met the requirements of the MID, meaning their maximum-permissible error (MPE) was ±0.5 %. The first three models of the VCDs had an allowed error of ±1 %. Within the selected VCDs we have also used a type 2 VDC (M4) that we calibrated as a type 1VCD. The sample size of the testing was 149, namely; 30 M1, 13 M2, 22 M3, 30 M4, 20 M5, 13 M6, 16 M7, and 5 M8 VCDs. Fig. 7 shows the results of the regression analysis yields for the M7 VDC that was MID approved. Fig. 7 shows very clearly the correspondence between the measurements and the error curve. This holds for all the analysed VCDs that met the requirements of the MID. 0.000 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000 9.000 10.000 11.000 12.000 Absolute pressure p [bar] Fig. 7. Comparison between the measured and approximated data for the M7 VDC in accordance with Eq. (18) at parameters r2 = 0.981 and SEE = 0.011 %% 3.1 Parameters of the Properties for the Regression Equations The most important statistical estimator for assessing the multi-regression model is the standard error of the estimate SEE. For the upper threshold value we took the testing facilities' measurement uncertainty of ±0.035 %. This means that the measurement errors of the VCDs, which have SEE values of SEE = ±0.035 % or less, which we calculated using a regression analysis, fit the regression model more than 90 %. This means that for every VCD that had its calculated SEE value via a multi-regression analysis lower than ±0.035 %, this regression model is suitable for the approximation. The second statistical estimator with which we assessed the approximation of the regression model, is the regression coefficient r2. It is more difficult to set the upper threshold value for the regression coefficient, because r2 is dependent on the selected regression curve. By considering the obtained results we can set a lower threshold value for r2 on 0.8, which means that the measured values lie 80 % on the regression curve. The VCDs with an r2 value lower than 0.8 have a worse fit. How well the regression model fits with the measurements determines the estimator of the significance F. The statistical significance F is calculated within the ANOVA analysis of the variance framework. The upper threshold value for the significance estimator F is 0.05. For all the VCDs approximations whose value is under 0.05, we can claim with a 95 % probability that at least one of their independent variables has a statistical effect on a dependent variable. In Table 5 we have a list of values for the estimators SEE, r2, and the significance F for the best and the worst VCDs that we checked during their testing. The second is the consecutive numbering of the tested VCD. As is clear from Table 5, the SEE values for the VCDs, apart from VCD M5-10, are within the selected threshold SEE values for SEE = ±0.1 %. For the VCD model M5 it turned out that no approximation of this VCD is sufficient for the required SEE values. From this we can conclude that because of the specific characteristics of the VCDs M5, the chosen regression model Eq. (18) is not suitable for approximating its measured quantities. The cause of such a poor approximation is the specific metrological characteristics of the built-in pressure and temperature sensors. Here we must add that the VCDs M5 have been out of production for a long time. The selected threshold SEE value of SEE = ±0.035 % is equal to the allowed share of the measurement uncertainty of the VCD testing facilities, which is why we also studied the other suitable values. We increased the allowed values in increments and observed how many approximated VCDs are suitable for a particular threshold value. The SEE value of SEE = ±0.075 % is of an approximately equal order of magnitude to the measurement uncertainty of the test facility for checking the electronic VCDs (±0.167 %). As it turns out, the results do not differ much from the SEE value of SEE = ±0.1 %. The combined measurement uncertainty of the testing facility and the approximation do not exceed ±0.3 %, which is also the acceptable value for the measurement uncertainties of the newest flow computers. Table 5. Values of statistical estimators of regression analysis for electronic VCDs VCD Statistical estimator (corrector) Standard Error r2 Significance F M1 M1-18 0.015 0.993 0.0000 M1-20 0.098 0.628 0.0706 M2 M2-5 0.031 0.843 0.0020 M2-2 0.095 0.823 0.0000 M3 M3-5 0.014 0.976 0.0000 M3-4 0.099 0.839 0.0022 M4 M4-12 0.012 0.989 0.0000 M4-21 0.089 0.311 0.5686 M5 M5-19 0.039 0.966 0.0000 M5-10 0.154 0.489 0.2248 M6 M6-6 0.018 0.971 0.0000 M6-3 0.044 0.123 0.9281 M7 M7-5 0.007 0.990 0.0000 M7-14 0.034 0.388 0.0000 M8 M8-5 0.008 0.968 0.0000 M8-1 0.027 0.774 0.0090 The data about the number of VCDs that are approximated with Eq. (18) and meet the criteria of the SEE are collected in Table 6. Besides the number, the relative shares of suitable approximations are also presented. The number of VCDs used for the approximation is indicated in the last column of the table. From Table 6 we can see that only 4 of the 149 approximations do not comply with the SEE value of SEE = ±0.1 %. All four approximations are of the model M5 VCDs, for which we predicted that we would need to find a better regression model. Table 6. Threshold values of SEE and the number of suitable VCDs SEE [%] Threshold limit Total < 0.035 < 0.05 < 0.075 < 0.1 number M1 number 13 20 28 30 30 % 43.33 66.67 93.33 100.00 M2 number 1 7 10 13 13 % 7.69 53.85 76.92 100.00 M3 number 13 20 21 22 22 % 59.09 90.91 95.45 100.00 M4 number 17 26 28 30 30 % 56.67 86.67 93.33 100.00 M5 number 0 1 10 16 20 % 0.00 5.00 50.00 80.00 M6 number 9 13 13 13 13 % 96.23 100.00 100.00 100.00 M7 number 16 16 16 16 16 % 100.00 100.00 100.00 100.00 M8 number 5 5 5 5 5 % 100.00 100.00 100.00 100.00 Table 7. SEE threshold values and the number of suitable VCDs excluding M5 VCDs _SEE [%]_ M1, M2, M3, M4, M5, M6, M7, M8 Limit_Total number_% < 0.035 74 57.36 < 0.050_107_82.95 < 0.075_121_93.80 < 0.100_129_100.00 Total 129 100,00 If we exclude all the model M5 VCDs from the analysis, then we only take into account 129 VCDs. The results are listed in Table 7. From these results collected in Tables 6 and 7 we can see that the VCDs with MID approval meet the SEE threshold already at 0.05 %. These results confirm our hypothesis, that it is possible to approximate the measurements of natural-gas VCDs sufficiently well with the multi-regression model in Eq. (18). The exactness of the approximation is dependent on the allowed threshold SEE values. Based on the results of the multi-regression analysis we can claim that the regression model only poorly approximates the measurements from the VCDs of model M5, and we may, as a consequence, exclude them from our analysis. However, it must be understood that our research did not demonstrate the model M5 VCD to be unfit for use, but because of its uncertainties our selected regression model cannot describe it well enough. 3.2 Reducing the Numbers of Measurement Points The optimum choice of measurement points for the VCDs testing in the field could be determined on the basis of [18] and [19]. The number of measurements points was obtained with the aim of reducing the calibration curve's uncertainty. The authors in [19] present the design of a calibration using experimental design techniques. The optimum calibration plan for the measurement chain is identified by suitably elaborating the error-propagation law suggested by the ISO Guide [20]. In our case we want to keep the valid methodology and national procedure during the testing procedure in our laboratory. Our findings were that the suggested regression Eq. (18) fits our measurements exceptionally well, especially for the MID-approved VCDs that are being produced recently. Based on these findings it stands to reason that we need only 9 measurement points and the regression analysis. The number of selected measurement points without repetition comes from the experimental design theory, Fig. 5. A reduced number of measurement points is proposed for the testing of the VCDs in the field. From this we propose the following reduced designated testing points, as depicted in Table 8. Table 8. The designated testing measurement points for 9 measuring points Temperature [°C] p1min Рз p5max t1 34 t2 64 t3 7^ 8^ ^9 Measuring point, Fig. 5 z1min z10 z1max z1min T1 T5 T3 z10 Тб T9 T8 z1max T2 T7 T4 For the example in Fig. 7 we have made an approximation with only 9 measurement points. The Table 9. The comparison of SEE values considering 9 or 15 measurement points fc [%] p t Approximation: 9 points Approximation: 15 points [bar] [°C] fc,apr fc,apr-fc [%] fc,apr fc,apr-fc [%] 0.06 0.98749 -0.01 0.069 0.009 0.066 0.006 -0.07 5.49396 0.00 -0.078 -0.008 -0.080 -0.010 -0.08 9.98705 -0.01 -0.081 -0.001 -0.076 0.004 0.08 0.98749 15.02 0.069 -0.011 0.064 -0.016 -0.08 5.49396 15.00 -0.078 0.002 -0.082 -0.002 -0.09 9.98705 15.00 -0.081 0.009 -0.078 0.012 0.07 0.98749 30.00 0.072 0.002 0.068 -0.002 -0.08 5.49396 30.00 -0.074 0.006 -0.078 0.002 -0.07 9.98705 30.00 -0.078 -0.008 -0.074 -0.004 SEE= ±0.0122 SEE= ±0.0080 difference in the SEE values, if we have approximated it with 9 points, or if we have calculated the coefficients of the estimator on 15 points, is less than 0.005 %. The data for the comparison of regressions is collected in Table 9. 3.3 The Combined Estimate of the Measuring Uncertainty with a Consideration of the Regression Analysis For our estimate of the measurement uncertainty we followed the GUM [20] and [21]. The measurement uncertainty of the measurement errors Ufc) is dependent on the measurement uncertainty of the testing facility for the control of VCDs Ufc)ms, which accounts for 1/3 of the measurement error, which is 0.167 %, and of the measurement uncertainty of the approximation of the measurement errors Ufc)apr. Table 10. The combined measurement uncertainty of the VCDs' measurement error U ( fc ) =V U ( fc ) ^ 2 + U ( fc ) 2 C ' apr (20) If we take into account the SEE because of the approximation, then the estimate of the measurement uncertainty for the correlation factor equals: U( fC) =4(1 / 3MPE)2 + (2 • SEE)2. (21) Taking in consideration our analysis, we can estimate the combined measurement uncertainty of the correction factor on the basis of Eq. (21), with the results collected in Table 10. Based on our estimate (Table 10) our combined measurement uncertainty is less than ±0.3 %, which is for all VCDs lower than the acceptable measurement errors. Modern VCDs must allow the entry of regression equation coefficients, so that the combined measurement uncertainty for standard or norm natural-gas volumes is less than ±0.5 %. SEE(fc) 1/3 MPE 2 SEE U(fc) SEE= ±0.000 % ±0.167 ±0.000 ±0.167 SEE= ±0.035 % ±0.167 ±0.070 ±0.181 SEE= ±0.070 % ±0.167 ±0.140 ±0.218 SEE= ±0.100 % ±0.167 ±0.200 ±0.260 Modern flow meters enable the measurement of natural-gas volumes to better than ±0.5 %. An example of this would be some ultrasonic flow meters, which can reach measurement uncertainties values of less than ±0.2 %. With the proposed measurement method it is possible to reduce the measurement uncertainty of the VCDs to less than ±0.3 %, which would enable a new generation of VCDs. 4 CONCLUSIONS Based on the multi-regression analysis results we were able to come to the following conclusions: • The most important statistical estimator for estimating the suitability of a multi-regression analysis used for any natural-gas VCD is the standard error of the estimate SEE. Additional important statistical parameters are the regression coefficient r2 and the significance F. • The proposed regression model is suitable for use with the control of all VCDs, especially modern MID-approved VCDs. However, the estimator is not suitable for the M5 VCDs that are out of production, because of the specifics of the sensors. • The VCD approximation suitability on the basis of the chosen regression model differs depending on the model of the VCD. Of all the tested VCDs, the regression model approximates VCD M7 with the best fit. • The regression model approximates the VCD measurements by considering the demanded SEE threshold value of SEE = ±0.035 % in 57.4 % of cases. And for an SEE threshold value of SEE = ±0.075 %, which approximately corresponds to the measuring uncertainty of the testing facility (±0.167 %), it suitably approximates for more than 93 % of the VCDs. • The proposed method of measurement represents an upgrade to the national verification procedures for the testing of VCDs. • In the case that we perform a multi-regression analysis of the VCD test on the basis of the reduced number of measurement points (9 points), and the results conform to the requirements, then the probability of an error is less than 25 %. In the case of model M4 the probability of an error was 13.3 %, so we may claim with a certainty of 86.7 % that the results will also conform to a 15-point test. Meanwhile, for VCDs M6, M7, and M8, the prediction is 100 %. A reduced number of measurement points is well suited to the field testing of VCDs. • In a comparison of the regression analysis based on 15 or 9 points, in most cases the SSE values will increase, as will the values of the r2 estimator, but the values of the significance F will, however, decrease in most cases. • Modern VCDs generally allow the entry of the regression-estimator constants from Eq. (18). • A regression analysis is especially valuable for VCD testing, because it can shorten the testing times. • We have demonstrated that we may use a regression analysis based on 9 points instead of 15, which is the keystone of field measurement procedures. • The proposed testing procedure can be automated. 5 REFERENCES [1] EN 12405-1:2005+A2:2010. Gas meters - Conversion devices - Part 1: Volume conversion. European Committee for Standardization, Brussels. [2] PTB - Prüfregeln: Band 20:1993. Elektronische Mangenumwerter für Gas. PTB, Braunschweig, Berlin. [3] Rules on Measurement Instruments (Pravilnik o spremembah in dopolnitvah pravilnika o merilnih instrumentih) (2006). Official Gazette of the Republic of Slovenia, nr. 42/2006. (in Slovene) [4] L135 (2004), L137 (2009). Directive 2004/22/EC of the European Parliament and of the Council of 31 March 2004 on measuring instruments. vol. 47, p. 1-80. [5] Oprešnik, M. (1992). Thermodynamics, University of Ljubljana, Faculty of Mechanical Engineering, 4th ed., Ljubljana (in Slovene). [6] EN 12213-1,2,3:2009. Natural gas - Calculation of compression factor. European Committee for Standardization. Brussels. [7] Buonanno, G., Carotenuto, A., Dell'Isola, M. (1998). The influence of reference condition correction on natural gas flow measurement. Measurement, vol. 23, p. 77-91, D0l:10.1016/ s0263-2241(98)00009-8. [8] Ficco, G., Dell'Isola, M., Vigo, P., Celenza, L. (2015). Uncertainty analysis of energy measurements in natural gas transmission networks. Flow Measurement and Instrumentation, vol. 42, p. 58-68, DOI:10.1016/j.flowmeasinst.2015.01.006. [9] Arpino, F., Dell'Isola, M., Ficco, G., Vigo, P. (2014). Unaccounted for gas in natural gas transmission networks: Prediction model and analysis of the solutions. Journal of Natural Gas Science and Engineering, vol. 17, p. 58-70, DOI:10.1016/j. jngse.2014.01.003. [10] Kramer, R., Mickan, B. (2009). Gas Measurement and Gas Meters Testing in Practise. PTB-Mitteilungen, vol. 119, no. 1, p. 16-22. [11] Blagojevič, B., Cimerman, F. (2013). Knowledge of natural gas chemical composition in Plinovodi d.o.o. transmission system. ZSIS, Ljubljana. (in Slovene) [12] PTB - Prüfregeln: Band 30:2003. Messgeräte für Gas, Hochdruckprüfung von Gaszählern. PTB, Braunschweig und Berlin. [13] Montgomery, D.C. (2013). Design and Analysis of Experiments. 8th Ed., John Wiley & Sons, Inc., Hoboken. [14] Širok, B., Blagojević, B., Bullen, P. (2008). Mineral Wool -Production and properties. Woodhead Publishing Limited, Cambridge, DOI:10.1533/9781845694456. [15] Spiegel, M., Schiller, J., Srinivasan, A. (2013). Schaum's -Probability and Statistics, 4th Ed., McGraw-Hill, New York. [16] Grabec, I., Gradišek, J. (2000). Description of the Random Phenomena. University of Ljubljana, Faculty of Mechanical Engineering, Ljubljana. (in Slovene) [17] IBM Software: SPSS Statistic Desktop version 22.0 for Windows. [18] Betta, G., Dell'Isola, M. (1996). Optimum choise of measurement points for sensor calibration. Measurement, vol 17., no. 2, p. 115-125, DOI:10.1016/0263-2241(96)00019-X. [19] Betta, G., Dell'Isola, M., Frattolillo, A. (2001). Experimental design techniques for optimising measurement chain calibration. Measurement, vol 30, p. 115-127, DOI:10.1016/ S0263-2241(00)00060-9. [20] ISO/IEC Guide 98-3:2008. Uncertainty of measurement - Part 3. Guide to the expression of uncertainty in measurement (GUM: 1995). International Organization for Standardization Geneva. [21] Adunka, F. (2007). Messunsicherheiten: Theorie und Praxis. Auflage 3, Vulkan-Verlag, Essen. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)2, 105-115 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2787 Original Scientific Paper Investigation of Non-Newtonian Fluid Effects during Transient Flows in a Pipeline Ali Majd1 - Ahmad Ahmadi1* - Alireza Keramat2 1 Shahrood University, Civil Engineering Department, Iran 2 Jundi Shapur University of Technology, Engineering Department, Iran A sudden change in the flow rate brings about significant pressure oscillations in the piping system, known as water hammer (fluid hammer). Unsteady flow of a non-Newtonian fluid due to instantaneous valve closure is studied. Power law and Cross models are used to simulate non-Newtonian effects. Firstly, the appropriate governing equations are derived and then, they are solved by a numerical approach. A fourth-order Runge-Kutta scheme is used for the time integration, and the central difference scheme is employed for the spatial derivatives discretization. To verify the proposed mathematical model and numerical solution, a comparison with corresponding experimental results from the literature are made. The results reveal a remarkable deviation in pressure history and velocity profile with respect to the water hammer in Newtonian fluids. The significance of the non-Newtonian fluid behaviour is manifested in terms of the drag reduction and line packing effect observed in the pressure history results. A detailed discussion regarding the fluid viscosity and its shear-stress diagrams are also included. Keywords: transient pipe flow, generalized Newtonian fluid, shear thinning fluids Highlights • Unsteady flow equations developed for non-Newtonian fluids. • The numerical method employed for non-Newtonian transient fluid flow equations. • Flow characteristics in a pipe section have been considered in detailed. • Significant changes observed in drag reduction and line packing effect. 0 INTRODUCTION Transient flows associated with the water hammer phenomenon are commonly encountered in both natural and engineering systems, such as hydraulic systems, oil transportation systems, and human arterial network. Sudden changes in pressurized pipe flow conditions caused by valve closure, pump operation, etc. are routine events. The excitations arising from these transient events can cause significant pressures leading to devastating forces [1]. The flow of non-Newtonian fluids and slurries in pipes occurs in a wide range of practical applications in the processing industries and many natural systems. If the fluid has a significant yield stress, or if its effective viscosity is high, industrially relevant flow rates may occur in the laminar flow regime. The fluids under consideration in this study are shear-thinning non-Newtonian, whose rheology is described by a generalized Newtonian fluid (GNF) model, i.e. the dependency of isotropic viscosity on flow properties [2] and [3]. In the specific fluids of the present work, the viscosity can be described using either the power law or Cross models. The capability of these models has been investigated by several researchers, including such as Pinho and Whitelaw [4], Toms [5] and Bird et al. [2] in experimental and numerical studies. To conduct computations on non-Newtonian fluids, the strain rate has to be evaluated. It requires a two-dimensional analysis to provide the velocity profile of a flow cross section. The two-dimensional analysis and computations of unsteady pressure and velocity profiles during water hammer have been developed by several researchers [6] to [10]. Pezzinga proposed a quasi-two-dimensional model for the unsteady turbulent flow of a pipe network and obtained better results than 1D models [11] and [12]. Vardy and Brown [13] have had significant contributions on non-Newtonian unsteady pipe flows especially for modelling fluids with time-dependent viscosities. More recently, Wahba [14] compared shear-thinning and shear-ticking fluids in response to a water hammer event using the power law model. Herein, unsteady pipe flow of a non-Newtonian fluid is studied. This work may be seen as a new extension to the classic water hammer model in which transients of a Newtonian fluid contained in a straight elastic pipe supported at the valve and along the pipeline with sufficient longitudinal anchors to suppress fluid-structure interaction effects is investigated. Having done this fluid hammer simulation for the power law and Cross models, several alternate works (in terms of studying the other effects [15]) on the transients of these fluids can be offered for future research. They include viscosity of the pipe wall, fluid-structure interaction, column separation, each of which or any of their combinations, e.g. Ahmadi and Keramat [16]; [15], Soares et al. [17], Hadj-Tai'eb and Hadj-Tai'eb [18] in conjunction with the present non-Newtonian behaviour can reveal new aspects of transient flow in possible systems of corresponding significance. Pezzinga et al. analysed transients in pressurized polymeric pipes using a two-dimensional (2D) KelvinVoigt viscoelastic model [19]. Differences between the transients in viscoelastic and elastic pipes are pointed out by considering a 2D model. They showed that viscoelastic models precisely represent a faster decay of pressure oscillations and velocity profiles because of a time-lag between pressure oscillations and retarded circumferential strain. Brunone et al. [20] and Kim [21] considered pressure and energy dissipation and unsteady friction in laminar transient flows. They compare their numerical result with those of experiments. Meniconi et al. also studied rapidly decelerating turbulent pipe flow. They proposed a new approach to estimate energy dissipation and pressure decay [22] and [23]. In the present study, laminar transient non-Newtonian pipe flow is simulated using the power and Cross models. To this aim, the quasi-2D equations of water hammer for non-Newtonian fluids are derived and then they are solved with appropriate numerical solutions based on the finite difference method. Computational results are provided in terms of velocity, shear stress, and viscosity distribution at the flow cross section in the middle of the pipeline. The results reveal that the non-Newtonian fluid effects significantly contribute to cross-sectional flow characteristics. 1 MATHEMATICAL MODELLING 1.1 Governing Equations The continuity equation for transient pipe flow in a cylindrical coordinate system is as follows [24]: +P dp dp — + vr — dt dr f 1 d(rvr ) r dr .YttdP. r de 1 dve I---e r dt dp ~dz dv __z_ dz = 0, (1) dv dv PI —- + vr —-n dt r dr dv дт dr ve dv. + ——- + v.—z r дв dz 1 дтв- дт т + ~ dz = Pgz + + + ^ + —, (2) r дв where т- (i,j=r, в, z) are stress components in the liquid in the corresponding surface and directions. To derive unsteady flow equations for a non-Newtonian fluid some assumptions and simplifications to be stated in the following are made. The flow is quasi-two-dimensional. The term 'quasi' indicates that ve , vr=0 meaning that vz is the only velocity component that varies along radial and axial directions. The convective terms are neglected [25] and [26]. The non-Newtonian fluid behaviour is based on power law and Cross models. The equation of state: dpc2 = pgdH [25] is valid, where H and c are pressure head and wave speed, respectively. With the application of the assumptions above as well as some algebraic simplifications, the continuity Eq. (1) is reduced to: c g (dH+v dH 21 dt z dz I dz 1 d(r dr = 0. (3) Integration over the flow cross section and neglecting convective terms leads Eq. (3) to: ~2 dV d_H+^ = о dt g dz (4) in which V is the average velocity over the flow cross section. Likewise, Eq. (2) with considering the assumptions above yields: dv dv dv dH 1 drx — + vz —~ + vr —— = -g-+--, (5) dt dz dr dz rp dr The convective terms in the above equation may be neglected so that it reduces to: dvz + dH 1 drr dt dz rp dr (6) The shear stress term on the right-hand side of Eq. (6) represents the fluid dynamic forces and is be calculated by the constitutive rheological fluid property. So far, no particular assumption is made for the type of fluid in governing equations; consequently, the above equations are valid for all fluid types. where vr , ve , vz are radial, angular and axial velocity components, p is density and t is time. The momentum equation in a cylindrical coordinate in an axial direction is: 1.2 Non-Newtonian Fluid Equations Non-Newtonian fluids may be classified into three general classes: time independent, time dependent and viscoelastic fluids [3]. Among them, the first branch is investigated in transient flows of the present study. Time independent fluids, which are placed in the inelastic fluids category, are also known as generalized Newtonian fluids (GNF). This category is similar to Newtonian fluids, but the shear stress and rate of deformation tensor are no longer a linear relation anymore. In fact, shear stresses are a nonlinear function of the rate of deformation. This nonlinear function is originated from natural features. On this basis, this category is divided into that with a yield stress and without yield stress. In the no yield stress group, there are two types: pseudo-plastics and dilatant fluids. The simulation of the former is the focus of this article. The viscous fluid flow is defined in terms of the velocity gradient that includes the rate of deformation and spin tensor. The constitutive relation between the shear stress т in Eq. (6) and the shear rate of the fluid t =nY yx ' (7) where n and у are apparent viscosity and shear rate respectively. Power law, Carreau, Cross, Ellis, etc. are different models that exist in the literature for the mathematical modelling of pseudo plastics, each of which has strengths and weaknesses [3]. In this study, the power law and the Cross model are applied. The power law is the simplest with the fewest possible parameters. It is described by the following equation: П= m (yx ) 1. (8) where m and n are two empirical curve fitting parameters and are known as the fluid consistency coefficient and the flow behaviour index respectively. In this equation, if n equals one and m is set to n0 , the Newtonian fluid is achieved. In this study, m is fixed to n0 while several quantities for n are selected. This allows for the investigation of n in the power law model during transient flows. The other model is the Cross model, which has the following description: Y=yF4lTD D =-1 2 du1 du1 dx1 dx1 where Dj is the rate of deformation tensor, and 11D is the second invariant of Dj. This representation of the shear rate in the r, в, z coordinate system with assuming the one-direction flow pattern reduces to [2] and [24]: 7: = = J-4 x - (d )2 - trD2 ) = dvz dr (12) 1.3 Initial Condition The fluid filled pipe is assumed to convey steady state flow before the transient event starts. So, the initial condition corresponds to the steady state flow. The momentum and continuity equations of the steady state flow can be written as: dP -2t0 dz dv _z_ dz R = 0. (13) (14) 1.4 Boundary Condition The transient flow in a reservoir-pipe-valve system is simulated. The quasi 2D analysis consists of three sets of boundaries including reservoir, valve, and internal pipe walls in contact with the flow. At the valve boundary, the velocity distribution is set to zero after the valve closure. A constant pressure head is associated with the reservoir boundary condition. The flow boundaries in contact with the pipe wall have zero velocity. These boundary equations can be written as follows: Vz(r=R) = 0 Vz( valve) = 0 H rezervoir = const. (15) n-n« 1 n1+k Y ) (9) where n and к are two fitting parameters whereas n0 and n^ are the limiting values of the apparent viscosity at low and high shear rates, respectively. In addition, for using shear rate and its independence from the coordinate system, Eq. (10) is applied [2]. 2 NUMERICAL METHOD The proposed unsteady flow equations in the previous section are solved using the finite difference method. A fourth order Runge-Kutta scheme is used to integrate the system of equations in time. Spatial derivatives are discretized using second order central differencing scheme. Second order dissipative terms are added to eliminate numerical oscillations. These terms perform only in the high gradient region, and they are effectively switched off in smooth regions [7]. To start with the numerical implementation, the two Eqs. (4) and (6), are combined to one matrix-form equation with the unknown vector W = {H V*}: dW „ dW ^ -+ В-= C, dt dx (16) where B and C are matrices of equation coefficients. The elements of vector W is sequentially evaluated from the discretized form of Eqs. (4) and (6) based on the Runge-Kutta scheme. The axial velocity profile is firstly evaluated from Eq. (6) so V* = Vz. Then its average is applied in the continuity Eq. (4) to calculate pressure head H, thus herein V* = V . The dissipative terms represented by A(W) are added to Eq. (16) to suppress the artificial numerical fluctuations: dW dW — + B— = C + A(W). dt dx (17) The employed dissipative scheme is based on the Jamson method [27] which is directly added to the basic equations. It is evaluated as follows: A(W) = ^kl/2 (W+i - W)-M-,/2 (W - W-)), (18) £i+i/2 =1 max ( ,ai+i ,ai,ai-i ) (19) in which a is a numerical variable that behaves like a switch, to be on or off on high and low gradients of the unknowns, respectively. In this article, a total variation diminishing (TVD) scheme is adapted to distinguish high values of gradients according to: a = \Wi+1 - 2Wi + (1 -m)WT mW ¥= \WM + 2). + ) ^TVD = \WM - W\ + W - W-11 • (20) (21) (22) If m is equal to zero, the above equation will be reduced to: - = -(23) W+ - W 1 - A a, =-1. (24) ' 1 +1 For calculating the average velocity of flow, the following equation is used, which is numerically computed by the Simpson integration scheme: fR 2nrv dr V = -J » z nRl (25) 3 MODEL VERIFICATION In order to validate the mathematical model and corresponding numerical solution and its implementation, the computations are compared with experimental results. To this end, an experiment done by Holmboe and Rouleau [28] on a reservoir-pipe-valve system with the following characteristics is considered. The pipe made of copper has an inner diameter of 0.025 m and a length of 36.09 m. Pressure signals directly upstream of the valve and at the pipe midpoint are recorded. The operating fluid in the laminar flow condition (Reynolds number = 82) is high-viscosity oil (p = 0.03484 Ns/m2), and the wave speed is measured to be 1324 m/s. The test is initiated by sudden valve closure that causes excision of fluid flow in the valve place and creates oscillations in pressure and velocity propagating along the pipe. Experiment results illustrate pressure values at various times after valve closure at two points along the pipe (valve and midpoint). In Figs. 1 and 2, non-dimensional experimental results taken from experiment [28] are compared with those of numerical results for the fluid pressure in the valve (Fig. 1) and midpoint (Fig. 2). 0 2 Fig. 1. Pressure time-history at the valve According to these comparisons, there is good agreement between the numerical and experimental results, thus validating the proposed mathematical model and numerical implementation. Fig. 2. Pressure time-history at midpoint of pipe The average velocity of fluid flow (computed using Eq. (25)) at the reservoir and midpoint cross sections are provided in Fig. 3 and, as can be seen, the computational results have a similar pattern to those of the conventional one-dimensional solutions. Another comparison is made for the non-dimensional axial velocity profile at the midpoint at Fig. 3. Computational average velocity history at midpoint and reservoir several time sections, being factors of the pipe length over the pressure wave speed (L/c). The velocity profiles and their gradients can be compared with the corresponding computations provided by [7] (Fig. 4). a) b) Fig. 4. Velocity profiles at the pipe midpoint for the laminar water hammer, a) present study, and b) Wahba [7]) aj b) Fig. 5. Pressure time history at the pipe midpoint for the laminar water hammer, a) the calculations of the present study, and b) Brunone et al. study [20] Again, the consistency between the two set of results verifies the implemented computer code. As more evidence for the correctness of the numerical model and its implementation, the results of Brunone et al. [20] were used. It explains that an experiment consists of a reservoir-copper pipe-valve system with 141.07 m length and 0.020 m inner diameter. The other specifications of the system are: Reynolds number is 815, pressure wave speed 1120 m/s, valve closure 0.11 s, and water temperature 17 °C. The modelling results of this study using the current simulation are compared with those of Brunone et al. [20] in Fig. 5. Fig 5b corresponds to the pressure heads obtained using a 1D model with unsteady friction [29] to [31] while Fig 5a is computed using the current 2D model. 4 INVESTIGATION OF NON-NEWTONIAN FLUID EFFECTS To recognise the significance of the non-Newtonian fluid behaviour that is manifested in viscosity variations during a transient flow, a couple of numerical examples are presented and discussed in detail via several figures. A pseudo plastic liquid that behaves as a shear thinning fluid is studied because it is the most common non-Newtonian fluid in applications. The coefficients of the power law model are chosen to be (n = По) and n = 0.8 and 0.6. In the cross model the initial viscosity (n0) is chosen to be equal to that of Newtonian fluid (defined in the previous section for the verifying case), and the ultimate viscosity (n^) equates to 20 % or 50 % of the initial viscosity. The two remaining parameters of this model are assumed to be n = 2/3 and к = 2. Note that this way of allocation of the initial viscosity in the two models enables the computational results to be favourably compared with the corresponding simulations for the Newtonian fluids so as to discriminate the deviations introduced to the flow characteristics as a result of the nonlinear fluid property. The set above of coefficients for the power and Cross models leads the viscosity values to those presented in Figs. 6 and 7 for the various shear rates. As can be seen, the power law model (Fig. 7) considerably suffers from a lack of accuracy in the regions of low shear rate while the Cross model (Fig. 6) is in accordance with the reality that is herein assumed to be the Newtonian constant viscosity. Several fluid properties defined via the power and Cross models are taken into account as the input fluid data for the transient flow analysis. The aim is to investigate transient pressures due to instantaneous downstream valve closure (see the previous section). Considering Joukowsky's pressure increase formula (AH=cV0/g), the value of the transient pressure just after the excitation is directly related to the pressure wave speed and initial velocity (steady state). Fig. 6. Cross model viscosity variations vs. strain rate 0.12 0.1 0.08 0.06 0 04 0.02 ^ \ \ \ \ \ ---power law n=0.8 -power law n=0.6 ......... power law n=0.9 -----Newtonian • ; 10" 10 10' 10' 10 У [1/s] Fig. 7. Power law model viscosity variations vs. strain rate These two quantities are kept unchanged so as to only scrutinize the non-Newtonian fluid properties during a transient flow. The defined various fluids, in turn, develop various pressure gradients and head loss during the steady state flow. The calculated head losses in the mentioned three power law cases are 2.693 cm, 1.343 cm and 0.6658 cm per metre and in the Cross model are 2.693 cm, 1.410 cm and 0.6404 cm per metre. According to the above conditions, the proposed numerical method produces the following pressures at the endpoint (valve) and the midpoint of pipe depicted in Figs. 8 and 9. Note that in all simulations, the flow pattern is laminar. The selected non-Newtonian fluids represent viscosities, though changing over the flow area but always smaller than that of the constant viscosity of Newtonian flow, (see Fig. 8 and Fig. 9). It means that smaller shear stresses develop on the pipe wall that correspond to less unsteady friction and causes to less pressure drop in the subsequent transient periods with respect to that of linear fluids. At the same time, during the first half period, the mentioned viscosity reduction of the shear thinning fluids enhances the fluid flow in the original direction and leads to less flow barrier and pressure gradient causing a smaller transient pressure rise in that time interval. In other words, a reduced packing effect is expected due to the reduced fluid viscosity, and this is in agreement with the computational figures that are provided for the pressure history at the valve and midpoint in Figs 8 and 9, respectively. The Figs 8a and 9a correspond to Cross and Figs 8b and 9b correspond to the power law model. Another manifestation of non-Newtonian fluids is observed in velocity profiles at various time and space sections. This is shown for the aforementioned shear thinning fluids in Figs. 10 and 11. In these figures, the velocity profile in the middle section of the pipeline for the Newtonian and two non-Newtonian fluids are compared. According to Figs. 10 and 11, the velocity distribution significantly changes as a result of the non-Newtonian fluid behaviour, and this change occurs throughout the unsteady fluid flow. The difference is such that with a reduction of fluid viscosity fluctuations, the variations of the velocity profile increase. In other words, the amplitude of the velocity gradient is increased in the flow cross-section. The growth in the velocity gradient in the vicinity of the pipe wall causes a drop in the viscosity value (shear thinning); in turn, this affects the values of fluid velocity and shear stress beside the pipe walls. Furthermore, the viscosity drop causes the maximum relative velocity to occur closer to the pipe walls. The velocity profiles also reveal that the central core area of flow has almost a rigid movement, and it is gradually affected by wall shear stress and viscosity variations of fluid near the pipe wall. In other words, the high values of wall shear stresses tend to penetrate in the core region and this pattern seems to be more progressive with the increase in viscosity variation а) июо b) u,uo Fig. 10. Velocity profiles at the midpoint for the power law fluid model (thick line) vs. Newtonian model (thin line); a) n = 0.6, and b) n = 0A a) "'"o b) « Fig. 11. Velocity profiles at the midpoint for the Cross fluid model (thick line) vs. Newtonian model (thin line); a) n» = 50 % no, and b) n» = 20 % По (with respect to the Newtonian fluid viscosity) in the present shear thinning fluids. This issue is explained in more detail by shear stress and viscosity distributions to be provided in the coming figures. According to the depicted shear stress and viscosity profiles (Figs. 12 to 16), their values in the core area of the pipe cross section remain almost unchanged. In fact, the flow in this region demonstrates a rigid movement (no relative displacement). According to the shear stress profiles at the pipe midpoint, the greater the shear thinning behaviour of the fluid, the less the expansion of the wall shear stresses to the core region of flow. Indeed, the area of rigid flow in the pipe cross section is extended, and the wall effects are more limited to the radial flow boundaries. In the meantime, the behaviours of the different non-Newtonian fluids in the core area of flow are remarkably similar to each other. The fluctuations of viscosity in Figs. 15 and 16 during the fluid hammer can be interpreted in terms of velocity profiles in Figs 10 and 11 and the viscosity variations versus the strain rate in Figs. 6 and 7. According to Figs. 6 and 7, the maximum value of the shear rate corresponds to the least value of viscosity, and this occurs at the pipe wall annulus. This can be found from the first derivative of velocity profiles with respect to the pipe radius. In contrast, the minimum value of the shear rate leads to the largest values of viscosity, which in the power law model 0.5 0.4 03 Q 02 0.1 0 -4 -2 0 2 4 6 T Fig. 12. Shear stress distribution for Newtonian model a) a) bi Fig. 13. Shear stress distribution for power law model; a) n = 0.8, b) n = 0.6 b) Fig. 14. Shear stress distribution for Cross model; a) = 50 % ц0, and b) = 20 % ц0 ai bi Fig. 15. Viscosity distribution in the pipe section for power law model; a) n = 0.8, b) n = 0.6 a) b) Fig. 16. Viscosity distribution in the pipe section for Cross model; a) цгл = 50 % , and b) цгл = 20 % ц0 tends to infinity and in the Cross model is a constant quantity called initial viscosity. According to the velocity profiles, around the central axis of the flow, the shear rate is zero, and it smoothly increases. This trend can be observed in the viscosity distribution, where it shows its largest value. There are more local maxima of viscosity in these figures that correspond to peaks in the velocity distribution. As an example, one can observe in Fig. 16b at t = 4L/c, which its maxima corresponds to the points indicated by arrows in Fig. 11b. The aforementioned figures of viscosity and shear stress can be used to interpret the pressure time history results of non-Newtonian fluids. In fact, the increase in the shear thinning property of a liquid corresponds to viscosity and shear stress variations in the annulus of the pipe cross-section that are closer to the pipe walls. This trend, which is clearly manifested in Figs. 12 to 16, can also result from the velocity profiles in Figs. 10 and 11. This behaviour of the shear-thinning fluids indicates that the region of more energy dissipation is limited to a smaller area, which in turn leads to less energy loss and pressure drops during the transient event. This issue is demonstrated by Figs. 8 and 9, which show that the pressure history of the liquid with the greater shear-thinning property shows smaller pressure drop over time. 5 CONCLUSION In this article, the non-Newtonian fluid effects in unsteady flows have been studied. Based on derived governing equations of transient non-Newtonian flows, a fourth-order Runge-Kutta numerical method has been used for the approximation of time phrases and second-order central difference scheme has been used for discretization in space. Furthermore, second-order dissipation phrases have been used for elimination of numerical fluctuations. In order to validate the proposed mathematical model and numerical solution, computational results have been compared with those of available experimental ones from the literature. The differentiating pattern between Newtonian and non-Newtonian flows which mainly stems from the nonlinear dependency of fluid viscosity on velocity gradient is observed in pressure variation, velocity profile and wall shear stress. Non-Newtonian power law and Cross models verified previously were then investigated through case studies to see the axial velocity profile at various times. Some of the most important results follow. The increase in the shear-thinning property of a liquid corresponds to the viscosity and shear stress variations in the annulus of the pipe cross-section that is closer to the pipe walls. This behaviour of the shear-thinning fluids indicates that the region of greater energy dissipation is limited to a smaller area, which in turn leads to reduced energy loss and pressure drops during the transient event. The pressure history of the liquid with the greater shear-thinning property shows a reduced pressure drop over time. Increasing the shear-thinning property of the non-Newtonian pseudoplastic fluid and, thus, the relative drop in the apparent viscosity decreased the amount of head loss in the pipe; comparatively, the pressure at the valve grows. Furthermore, due to the reduction of the apparent viscosity at the wall, a reduced line-packing effect is observed compared to Newtonian models. The shear-thinning behaviour of non-Newtonian fluids causes the region of high gradient velocities to move towards the pipe walls and the maximum relative velocities occur closer to the radial boundaries, thus leading to severe fluctuations in the cross-sectional velocity profile. The provided 2D computational results at a cross section reveal the significance of the non-slip boundary of the inner pipe wall during flow transients in terms of lags between the mean flow direction and velocities at several radii. The non-Newtonian fluid effect tends to vary the velocity profiles at each time. In the midpoint of the pipe, a semi-rigid movement with negligible relative velocity variations was observed which is illustrated in terms of viscosity variations. According to the velocity profiles, around the central axis of the flow, the shear rate is zero, and it smoothly increases. This trend can be observed in the viscosity distribution where it shows its largest value. There are more local maxima of viscosity in these figures that correspond to peaks in the velocity distribution. 6 REFERENCES [1] Wylie, E.B., Streeter, V.L.A. Suo, L. (1993). Fluid Transients in Systems. Prentice Hall, Upper Saddle River. [2] Bird, R.B., Armstrong, R.C. Hassager, O. (1987). Dynamics of Polymeric Liquids, 2nd ed., vol. 1. Wiley Interscience Publication, Hoboken. [3] Chhabra, R.P., Richardson, J.F. (2011). Non-Newtonian Flow and Applied Rheology: Engineering Applications, 2nd ed. Elsevier Science, Oxford. [4] Pinho, F.T. Whitelaw, J.H. (1990). Flow of non-Newtonian fluids in a pipe. Journal of Non-Newtonian Fluid Mechanics, vol. 34, no. 2. p. 129-144, DOI:10.1016/0377-0257(90)80015-R. [5] Toms, B.A. (1948). Some observation on the flow of linear polymer solutions through straight tubes at large Reynolds numbers. Proceedings of 1st International Congress on Rheology, vol. 2, p. 135-141. [6] Ghidaoui, M.S., McInnis, D.A., Axworthy, D.H., Zhao, M. (2005). A review of water hammer theory and practice. Applied Mechanics Reviews, vol. 58, no. 1, p. 49-76, DOI:10.1115/1.1828050. [7] Wahba, E.M. (2006). Runge-Kutta time-stepping schemes with TVD central differencing for the water hammer equations. International Journal for Numerical Methods in Fluids, vol. 52, no. 5, p. 571-590, DOI:10.1002/fld.1188. [8] Riasi, A., Nourbakhsh, A., Raisee, M. (2009). Unsteady velocity profiles in laminar and turbulent water hammer flows. Journal of Fluids Engineering, vol. 131, no. 12, p. 121202, DOI:10.1115/1.4000557. [9] Brunone, B. Berni, A. (2010). Wall Shear stress in transient turbulent pipe flow by local velocity measurement. Journal of Hydraulic Engineering, vol. 136, no. 10, p. 716-726, DOI:10.1061/(ASCE)HY.1943-7900.0000234. [10] Brunone, B., Karney, B., Mecarelli, M., Ferrante, M. (2000). Velocity profiles and unsteady pipe friction in transient flow. Journal of Water Resources Planning and Management, vol. 126, no. 4, p. 236-244, DOI:10.1061/(ASCE)0733-9496(2000)126:4(236). [11] Pezzinga, G. (1999). Quasi-2D Model for Unsteady Flow in Pipe Networks. Journal of Hydraulic Engineering, vol. 125, no. 7, p. 676-685, DOI:10.1061/(ASCE)0733-9429(1999)125:7(676). [12] Pezzinga, G. (2000) Evaluation of Unsteady Flow Resistances by Quasi-2D or 1D Models. Journal of Hydraulic Engineering, vol. 126, no. 10. p. 778-785, DOI:10.1061/(ASCE)0733-9429(2000)126:10(778). [13] Vardy, A.E., Brown, J.M.B. (2011). Laminar pipe flow with time-dependent viscosity. Journal of Hydroinformatics, vol. 13, no. 4, p. 729-740, DOI:10.2166/hydro.2010.073. [14] Wahba, E.M. (2013). Non-Newtonian fluid hammer in elastic circular pipes: Shear-thinning and shear-thickening effects. Journal of Non-Newtonian Fluid Mechanics, vol. 198, no. 0, p. 24-30, DOI:10.1016/j.jnnfm.2013.04.007. [15] Keramat, A., Tijsseling, A.S., Q. Hou, Ahmadi, A. (2012). Fluid-structure interaction with pipe-wall viscoelasticity during water hammer. Journal of Fluids and Structures, vol. 28, p. 434455, DOI:10.1016/j.jfluidstructs.2011.11.001. [16] Ahmadi, A., Keramat, A. (2010) Investigation of fluid-structure interaction with various types of junction coupling. Journal of Fluids and Structures, vol. 27, no. 7-8, p. 1123-1141, DOI:10.1016/j.jfluidstructs.2010.08.002. [17] Soares, A.K., Covas, D.I.C., Carrigo, N.J.G. (2012). Transient vaporous cavitation in viscoelastic pipes. Journal of Hydraulic Research, vol. 50, no. 2, p. 228-235, DOI:10.1080/00221686 .2012.669143. [18] Hadj-Taieb, L., Hadj-Taieb E. (2009). Numerical simulation of transient flows in viscoelastic pipes with vapour cavitation. International Journal of Modelling and Simulation, vol 29, p. 206-213, DOI:10.2316/Journal.205.2009.2.205-5100. [19] Pezzinga, G., Brunone, B., Cannizzaro, D., Ferrante, M., Meniconi, S., Berni, A. (2014). Two-Dimensional Features of Viscoelastic Models of Pipe Transients. Journal of Hydraulic Engineering, vol. 140, no. 8, p. 04014036, DOI:10.1061/ (ASCE)HY.1943-7900.0000891. [20] Brunone, B., Ferrante, M., Cacciamani, M. (2005). Decay of Pressure and Energy Dissipation in Laminar Transient Flow. Journal of Fluids Engineering, vol. 126, no. 6, p. 7, DOI:10.1115/1.1839926. [21] Kim, S. (2011). Holistic Unsteady-Friction Model for Laminar Transient Flow in Pipeline Systems. Journal of Hydraulic Engineering, vol. 137, no. 12, p. 1649-1658, D0I10.1061/ (ASCE)HY.1943-7900.0000471. [22] Meniconi, S., Duan, H., Brunone, B., Ghidaoui, M., Lee, P., Ferrante, M. (2014). Further Developments in Rapidly Decelerating Turbulent Pipe Flow Modeling. Journal of Hydraulic Engineering, vol. 140, no. 7, p. 04014028, DOI:10.1061/(ASCE)HY.1943-7900.0000880. [23] Meniconi, S., Brunone, B., Ferrante, M., Massari, C. (2014). Energy dissipation and pressure decay during transients in viscoelastic pipes with an in-line valve. Journal of Fluids and Structures, vol. 45, no. 0, p. 235-249, DOI:10.1016/j. jfluidstructs.2013.12.013. [24] Lai, W.M., Knovel, Krempl, E., Rubin, D. (2010). Introduction to Continuum Mechanics. 4th ed., Butterworth-Heinemann, Elsevier, Amsterdam, Boston. [25] Tijsseling, A.S. (1993). Fluid - Structure Interaction in Case of Waterhammer with Cavitation. PhD Thesis, Delft University of Technology, Delft. [26] Tijsseling, A.S. (1996). Fluid-structure interaction in liquid-filled pipe systems: A review. Journal of Fluids and Structures, vol. 10, no. 2, p. 109-146, DOI:10.1006/jfls.1996.0009. [27] Jameson, A., Schmidt, W., Turkel, E. (1981). Numerical Solutions of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-stepping Schemes, AIAA, Fluid and Plasma Dynamics Conference, Palo Alto. [28] Holmboe, E.L., Rouleau, W.T. (1967). The Effect of Viscous Shear on Transients in Liquid Lines. Journal of Basic Engineering, vol. 89, no. 1, p. 174-180, DOI:10.1115/1.3609549. [29] Zielke, W. (1968). Frequency-Dependent Friction in Transient Pipe Flow. Journal of Basic Engineering, vol. 90, no. 1, p. 109115, DOI:10.1115/1.3605049. [30] Bratland, 0. (1986). Frequency-Dependent Friction and Radial Kinetic Energy Variation in Transient Pipe Flow. Proceedings of the 5th International Conference on Pressure Surges, p. 95101. [31] Vardy, A.E., Hwang, K. (1991). A characteristics model of transient friction in pipes. Journal of Hydraulic Research, vol. 29, no. 5, p. 669-684, DOI:10.1080/00221689109498983. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)2, 116-126 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2495 Original Scientific Paper Analysis of the Influence of Pavement Irregularities on the Lifespan of a Vehicle's Drive-Wheel Half Shaft Sreten Simović1* - Vladimir Popović2 - Milanko Damjanović1 1 University of Montenegro, Faculty of Mechanical Engineering, Montenegro 2 University of Belgrade, Faculty of Mechanical Engineering, Serbia As a complex technical system, a vehicle requires a full analysis of its dynamic behaviour. A vehicle is the research subject requiring a detailed analysis of all separate systems both individually and in mutual interaction. One of the technical aspects of vehicle analysis is the lifespan of each element. This paper presents an analysis of the impact of the characteristics of pavement irregularities on the lifespan of a half shaft. In this paper, a mathematical model and a computer simulation of vehicle half-shaft loads have been set up. The importance of this approach is reflected in the use of state-of-the-art computational methods of analysis. The obtained results enable the drawing of clear conclusions about the effects of pavement properties on the lifespan of a vehicle half shaft. Keywords: vehicle, pavement, half shaft, load, lifespan Highlights • The influence of pavement irregularities on a vehicle drive-wheel half shaft is analysed. • A mathematical model of vehicle half-shaft loads has been set up and the influence on lifespan is analysed. • Simulation results of dynamic behaviour were tested and verified through experiments. • Results show that the presented analytical model can be used for the analysis of a vehicle's behaviour when wheel passing over irregularities. • The pavement characteristics significantly affect a vehicle half shaft's stress and lifespan. 0 INTRODUCTION During their service life, vehicles are faced with a wide spectrum of vibrations. These vibrations are primarily triggered by surface irregularities and internal sources, where the forces and moments generated by the contact between a surface and wheels have the greatest influence [1] to [3]. In addition to the demands of efficacy, emissions, and safety, vehicles are simultaneously faced with many demands that are reflected in the achievement of comfort, driveability, low wear, availability and long-term functionality [4]. Besides the complexity of the vehicle structure, dynamic behaviour and interaction between vehicle systems and the environment, most analyses consider only individual parts of vehicle systems [5]. In that sense, the vehicle drivetrain system is frequently the subject of research, which is to be expected given its importance in its interactions with other vehicle systems. Thus, a great deal of reference literature sources shows the analysis of a vehicle drivetrain system with road irregularities as the input in the system or in engine and vehicle models to emulate vehicle drivetrain loads [6] to [8]. There is a noticeable lack of analyses and conclusions on the influences of the vehicle systems' characteristics and pavement characteristics on the lifespan of the vehicle drivetrain system elements. Indeed, the complexity of a vehicle structure and the behaviour of the deformable wheel when it passes over an irregularity is a huge problem in this analysis, because there is no complete and satisfactory theory to explain that process [9]. Regarding the process of force and torque generation, the characteristics of a tyre are also very important. A wheel has a function to protect a vehicle against surface irregularities, which requires a deep analysis to provide the characteristics of a vehicle use and the need for the accomplishment of a long lifespan. Such an analysis must provide a prediction of behaviour according to drive inputs and influences of the surface irregularities and surface properties. The theoretical research frequently argues that the dynamic vehicle behaviour also depends on the characteristics of change of a vertical force that, among other factors, depends on the characteristics of a suspension system [10]. The previously stated indicates the conclusion that the analysis of vehicle systems' dynamic behaviour must include the influence of a wheel as a deformable element, the determination of reliance between the displacement of the unsprung and sprung masses and analysis of the oscillatory system consisting of a wheel, drive axle, suspension, and a vehicle frame [6]. Papers [11] and [12] present an analysis of the various vehicle suspension systems. There are many tyre models, e.g. Magic Formula, Fiala, Delft Tyre 97, etc. [1]. However, the fact is 116 *Corr. Author's Address: University of Montenegro, Faculty of Mechanical Engineering, Džordža Vašingtona bb, Podgorica, Montenegro, sretens@ac.me that when a model is simpler, it will possibly be less accurate, while a more complex model will be less suitable for the vehicle motion simulations. A detailed analysis of the wheel dynamics is usually done with the data obtained from the performed experimental investigations where the wheel-pavement contact is analysed and the conclusions about the particular characteristics were drawn, for example that the main parameters of the tyre influencing the generated force are the vertical and the radial stiffness of a tyre, where the pneumatic damping characteristic may be neglected, [13] to [15]. The need to perform an analysis of a vehicle systems element lifespans is also clearly indicated. The importance of this analysis is proven by fact that a vehicle durability depends on design geometry, material properties, and load, with road roughness as an important parameter influencing the vehicle load. The analysis of the lifespan of vehicle drivetrain system elements in relation to the diversity of pavement characteristics shows that the design of vehicle systems should be adjusted to the vehicle exploitation condition, [10]. Thus, different approaches in the process of solving certain problems in a lifespan analysis are used [16]. In the analysis of lifespan, the real conditions are usually replaced by the values of the estimated forces and stresses acting on the vehicle system elements [17] and [18]. The fatigue and reliability tests may be performed in terms of analysing different roads, analysing interacting influences of some mass and geometrical parameters of the vehicles or analysis of the relative fatigue damage of the drivetrain systems for specific roads [19] to [21]. 1 DETERMINATION OF LOADS AND ESTIMATION OF LIFESPAN The analyses of vehicle elements' lifespan are usually conducted in terms of a service life and an analysis of the load of some vehicle elements. In this regard, an analysis investigating the interactive influence of the vehicle elements and the necessity of their analysing as a unique dynamic system is very important [22]. Due to the complexity of the analysed system, any analysis of the lifespan of the vehicle drivetrain requires a meticulous approach. Therefore, in the process of the development of a vehicle the measurements and simulation data must be combined; for instance, in that process the vehicle performance may be predicted fairly exactly, but a fatigue life estimation, as a result of a more complex analysis, permits only relative statements, and they are always supported by tests [23]. In that sense, [24] presents a system for the continual tracking of working conditions and the estimation of a remaining lifespan of automotive transmission gear wheels, while the model of a drive system and the damage of a gearwheel are presented in [25]. Within the analysis of a lifespan of a driving axle, the time of a vehicle use in some gear ratios and types of road must be taken into account [26] and [27]. Since vehicle motion is characterised by a constant change of wheel-road interactions, i.e. generated forces and stresses, thence determination of a variability, values and laws on the distribution of working and permissible stresses in relation to road conditions represent a major problem in the analysis [28]. The analytic determination of the loads requires a mathematic description of a dynamic behaviour of a drivetrain and the processing of the stresses in elements. When the stress recording is processed, its discretisation is carried out in a manner in which the real load is replaced by a load equivalent to the original load. In real exploitation conditions, the stress changes are multi-cyclic with a changeable amplitude and a mean stress. The elements' strength in such a change of loads is a working strength that can be determined by experimental testing, using the basic strength or in a theoretically-experimental way using the hypothesis about the fatigue accumulation. Working strength is most accurately obtained by direct experimental testing, but this is mostly not acceptable because comprehensive, long-term, and expensive research is necessary. The hypotheses of the material damage accumulation enabling the calculation estimation of an element's lifespan is widely used in the determination of working strength. The first hypotheses about a material damage accumulation were developed by Palmgrin, and Miner expanded them [28]. Palmgrin-Miner's damage hypothesis represents a basis for other hypotheses about the damage accumulation, on the basis of which the total damage is calculated as a sum of individual items of damage due to stresses and assume that a fatigue crack occurs when the interaction function, ar, is equal to 1. The value of interaction function is: = V —; for г. > : ^N 1 (1) where щ is the number of cycles of the stress of the ith level, t, for the entire lifespan, N is the number of changes of stresses on the ith level that an element can endure, j is the total number of the levels of stress, td is the fatigue limit strength, Fig. 1. a r =1 Fig. 1. Wöhler's diagram In the total number of cycles, the changes of stress of all weights to the failure of an element NR the participation of each stress Ti is comparable to the participation of this same stress in a unit spectrum nbi/nb namely: n = Nr , (2) where NR is the total number of changes of working stresses of all levels that an element can endure to a fracture, nbi is the number of changes of stresses of the ith level in a unit spectrum, nb is the total number of changes of stresses of all levels in a unit spectrum. From the equation of a curve of a dynamic endurance, Wöhler's curve, we obtain: T ■ Ni =tD ■ Nd =t1 ■ Nu ( Y N = Nn = N, f \m T (3) (4) where ND is a base cycle number, N1 number of changes of maximum working stress, t1 maximum working stress, Ti working stress on ith level. Using Palmgrin-Miner's damage hypothesis, the component's service life in load cycles, NR, is calculated according to: Nr =- Nd ( Y -; for г. >td (5) For the purpose of the improvement of Palmgrin-Miner's damage hypothesis, several corrected linear hypotheses were developed and the most important are Corten-Dolan's, Haibach's and Sorensen-Kogaev's hypotheses. Corten-Dolan's damage hypothesis adds the elementary damage along the modified fatigue line for all sorts of stresses, and, according to Haibach's damage hypothesis, the influence of stresses with values lower than the fatigue limit strength is assumed according to a fictitious line on Wöhler's diagram without the removal of smaller stresses. Both hypotheses assume that a fatigue crack occurs when the interaction function, ar, is equal to 1, [28]. The comparison of previously mentioned hypotheses is presented on Wöhler's diagram given in Fig. 2. Sorensen-Kogaev's damage hypothesis represents a modification of Palmgrin-Miner's hypothesis where a fatigue crack occurs when a sum of relative damages reaches the value of the calculated interaction function. The impact on the dynamic strength of stresses with values lower than the fatigue limit strength is assumed by the coefficients of a stress spectrum shape by which the calculation accuracy is improved in relation to Palmgrin-Miner's hypothesis. With the use of Sorensen-Kogaev's damage hypothesis, the calculation of the influence of a stresses with values lower than the fatigue limit strength on a working strength of elements is carried out for each specific case. Log N(n) Fig. 2. The comparison of hypotheses The value of interaction function for Sorensen-Kogaev's damage hypothesis is calculated using Eq. (6): T t n t 'Z - k-Zr ктв <"1 T - к ' Z D (6) where к is a constant determining the bottom limit of stresses causing the damage of the material (к = 0.5). Using Sorensen-Kogaev's damage hypothesis, the component's service life in load cycles, NR, is calculated according to: bi b Ur = KTD J b Nr = a • Nr »=i nb / \m T V ld ; ■; for t» > TD. (7) In order to compare the results of use of the mentioned hypothesis, some papers show results in which a reliability evaluation in use of Palmgrin-Miner's and Corten-Dolan's hypothesis leads to higher values of a driving shaft reliability than reliability evaluation using Haibach's and Sorensen-Kogaev's hypothesis, yielding the result for a lifespan that is closer to a real life cycle. The comparison of the results of a reliability evaluation in the use of Haibach's and Sorensen-Kogaev's damage hypothesis shows higher values of driving shaft reliability for road conditions with a higher share of roads being macadam and roads with a lower quality of a pavement surface in use of Haibach's damage hypothesis, [29]. 2 ANALYSIS OF THE INFLUENCE OF PAVEMENT IRREGULARITIES In order to conduct the analysis, relatively simple mathematical models of the drivetrain and suspension systems and the pavement characteristics were created, using the schematic for vehicle systems shown in Fig. 3. Fig. 3. Vehicle systems The simulation of the dynamic behaviour of the vehicle system presented by equivalent mass, momentum, stiffness, and damping was performed through the linearized state-space equations of the quarter car model, as follows: a. the drivetrain: • engine: Ш. +^ .Ah, (8) dm dh Je ' А'Фе = * ЛЧ - [К ' ( - АфоШ,е ) + +c (e -АФout,e ) , (9) where ke and ce are the equivalent stiffness and damping of the engine, Je the momentum of the engine, фе the angle of rotation of the crankshaft, pout,e the angle of rotation of the crankshaft, h the position of the accelerator and Me the output torque, gearbox: Ji +ГТ ■ Ф = M-- те n.g k ■ Ч — -Ф . i Tout.g + c. Фе ■ — -Ф . Tout.g V g (10) where Mpg is gearbox input torque, kg and cg the equivalent stiffness and damping of the gearbox, J1, J2 the momentums of the gearbox presented with two masses, yg the angle of rotation of the input gearbox shaft, yout>g the rotation angle of the output transmission shaft and ig the transmission ratio, • half shaft: Jhs ■ phs = Mp,hs - [khs • (Phs - Souths ) + + Chs -((phs (out,hs )], (11) where MPhhs are the half-shaft input torque, khs and chs are the equivalent stiffness and damping, Jhs the momentum of the half shaft, yhs the rotation angle of the half shaft and yout>hs the rotation angle of the output, differential: r J2 J1 +—2 + i2 J J •h = (p(l) • io )2 ((d) • io ) = Mp,d -\_Mres,f(D + Mres,r(r)] • (12) M s,f (() lo • lf (() Kout,f (l ) Vd lo ' lf (l) -m 'out,/(() +c 'out,/ (l) Vd l0 ' lf(() --m M res,r(r) io ■ iHr] +c 'out.r (r] kout,r(r] Pd out,f (( ) Pd (13) io ■ ir{r] -Pout r(r --P 'out,r(r] (14) g 1 1 where Mpd are the differential input torque, /0>д1 ^ the gear ratios in the differential, kout,m (15) (16) (17) (18) epn,ekv + A ' B ' epn,ekv + A ' ßpn,ekv A ' C ' epn,ekv,ul, (19) epn,ekv epn,ekv0 (for Zp,t,peg 0), (20) where Md is the wheel input torque, kpn, t and cpn, t are the tangential stiffness and damping of the tyre, Jrim, Jpn the momentums of the wheel rim and tyre, yrim, фрп the rotation angle of the wheel rim and tyre, Mres the resistance torque, Fpn,z,ekv the equivalent vertical force on the tyre, epn, ekv the equivalent arm of the vertical force, r0 the outer radius of the tyre, s the distance travelled in 0.01 s and zptpeg the effective value of the irregularity height, A, B, C constants. b. suspension system, Fig. 6: Fig. 6. Suspension system m ■ "z = k -(z — z ) + c -Iz — z I, (21) sp sp s \ unsp sp 1 s \ unsp sp I? X^^J nsp Zunsp (zp,t Zunsp Cr,w (zp,t Zunsp k -(z - z ) + c -('z - z VI • cos S •a, (22) s \ unsp sp J s \ unsp sp f J b Fpn,z,ekv = kr,w ' \_zp,t — zunsp ] + Cr,w ' p,t ~ 'Zunsp ] , (23) F = k -[z t - z 1 + C -['z , -'z 1, (24) pn,z r,w I p,t,peg unsp I r,w I p,t,peg unsp P v ' where a and b are the distances in the vehicle suspension system, S the angle of inclination of the suspension system, ks and cs the equivalent stiffness and damping of the suspension system, msp and munsp the sprung and unsprung masses, krw and crw the radial tyre stiffness and damping, zsp and zunsp the position of the sprung and unsprung masses, Fpnz and Fpnz,ekv the vertical and equivalent vertical force on the wheel and zp, t and z^g tl irregularity height. and zp, t and zptpeg the road irregularity and equivalent c. road irregularities: As the tyre is a deformable element of the system simultaneously in contact with a number of irregularities, it is necessary to take into account the impact of the deformable tyre that is reflected in smoothing irregularities. This characteristic of a tyre is given in the form of the transfer function of the dynamic chain in which the input value is the ordinate of the road surface micro-profile and the output obtained is its mean height at the length of the contact zone between the tyre and the road surface, [6] and [30]: K, (jmf + jm- kpl -V2+ k (25) where kpl is a coefficient that is calculated by the expression: кр! =(0.9 .O)-V (26) where lt is the length of contact between the road surface and tyre, and v is the vehicle speed. The parameter lt is obtained by using the expression: I, = 2 -J 0.1 • Ht •(( - 0.1-Ht ), (27) where Dt and Ht are the outer radius and the height of the tyre profile. The spectral density of the process smoothed by the tyre Speg is obtained by the transfer function: H = Wpeg (HWg (- H^h H Speg H) = К, (( -H ) + 2 • kpl H • S (h), (28) where Sh is the spectral density of the micro profile. Based on the transfer function shown above, the following differential equation is obtained, enabling the calculation of the equivalent height of an irregularity. This is relevant to the calculation of the behaviour of the system under the influence of the irregularities: 'Zp,t,peg + 'Zp,t,peg ' kpl + Zp,t,peg ' kpl kpl ' Zp,t ' (29) where zp t is the accidental height of the irregularity. 3 VALIDATION OF MATHEMATICAL MODEL The validation of a mathematical model requires the comparison of the obtained results with the data about the real system behaviour. In that sense, before the execution of an analysis of a lifespan, a validation of a model via the comparison of results of computer simulation with the results of the experimental testing should be carried out [4]. Regarding the necessary steps in the analysis, a previously performed individual analysis of the steps for passenger cars with front drive shafts and McPherson struts was used [31]. The experiment was performed on a real vehicle, where the engine drive torque and the number of revolutions were maintained at an approximately constant value. The vehicle was moving in a straight line over an ideally flat concrete surface at a constant speed of approximately 20 km/h and over a triangular irregularity with a height of 55 mm and length of 190 mm. The validation was done via the comparison of data obtained from the simulations of a one-quarter car model, performed under the conditions in which the experiment was conducted. The comparison was performed for a change of a spring-damper length, ^unsp zsp, and for the change of drive torque on the half shaft, Md. The parameters used for a computer simulation are shown in Table 1 [31]. The results of the validation are shown in Figs. 7 and 8 which show the results of experimental testing and computer simulations. The diagrams presented in Figs. 7 and 8 clearly indicate that the simulation and the experimentally obtained results regarding the suspension and drivetrain system show a very good compliance, by which it can be concluded that the presented analytical model has an advantage and can be used for the analysis of vehicle behaviour when passing over irregularities. 4 SIMULATION OF DYNAMIC BEHAVIOUR The computer simulation of the dynamic behaviour of the drivetrain and suspension systems under the influence of disturbances caused by pavement Table 1. Simulation parameters System Symbol Value Unit M0 150 Nm dMe / drne -0.5 Nms Drive engine Je 0.15 kgm2 ke 3 kNm/rad ce 0.1 Nms/rad J1, J2 0.02 kgm2 ig 3.6 - Gearbox kg 3 kNm/rad cg 0.1 Nms/rad J1, J2, J3, J4 0.02 kgm2 i0 3.58 - Differential kout,f(1) 3.5 kNm/rad Cout,f(1) 0.1 Nms/rad Jhs 0.0003 kgm2 Half shaft khs 5 kNm/rad Chs 0.1 Nms/rad Jrim 14 kgm2 k "r.w 220 kN/m c - - kpn,t 4.5 kNm/rad cvn.t 2 Nms/rad Drive-wheel Jph 0.3 kgm2 r0 0.293 m e0 0.04 m A 130 1/s B 1.5 - C 3 - msp 280 kg ks 19 kN/m cs 2 kNs/m Suspension m 45 kg system s 7.6 о a 0.27 m b 0.32 m kpl 1.1 1/s Table 2. Pavement characteristics (coefficient of approximation of the auto-correlative function) Road surface pavement Dk A [cm2] 1 A2 a1 a2 ß 1. Asphalt (good) 0.664 1 0 0.13 0 1.05 2. Asphalt (used) 1.21 0.15 0.85 0.05 0.2 0.6 3. Gravel 6.3 0.047 0.9553 0.049 0.213 1.367 4. Stubble (field) 10.63 0.1 0.9 0.2 0.7 1.57 irregularities, for the purpose of an analysis of a lifespan of a vehicle half shaft, was done using the previously presented analytical model. The pavement characteristics are presented as a stationary ergodic 0.7 0.8 Time, t [s] Fig. 7. Change of front axle spring-damper length 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 Time, t [s] Fig. 8. Change in the drive torque on the half shaft random process with an auto-correlative function, Rk(t), defined using the expression: Rk(t) = D "[A ■ ^ ■ co.(ßk ■ t)+A* ■ ^]. (30) Therefore, simulations of the dynamic behaviour of the system when crossing four types of pavement were made, and pavement characteristics are given in Table 2 [30]. The dynamic behaviour was simulated in MATLAB Simulink, in a time corresponding to a vehicle motion for 30 seconds, with a step of random excitation of the pavement irregularity of 0.01 seconds. Reflecting that, Figs. 9 and 10 show the best and worst pavement characteristics in terms of surface irregularities (pavements 1 and 4 respectively), and Figs. 11 to 14 show the change in a torque on the half shaft when a drive-wheel crosses over the four characteristic pavements. 5 ANALYSIS Since the measurement readings of a random change in the load cannot be used directly in the calculation of the lifespan of a vehicle half shaft, the readings were converted into a suitable form. To achieve this, Fig. 9. Changes in the random pavement height profile under the drive-wheel when passing over Pavement 1; (gray thiner line - random input, gray thicker line - smoothed pavement height characteristic) Fig. 11. Changes in torque on the half shaft when drive-wheel passing over Pavement 1 Fig. 10. Changes in the random pavement height profile under the drive-wheel when passing over Pavement 4 (gray thiner line - random input, gray thicker line - smoothed pavement height characteristic) Fig. 12. Changes in torque on the half shaft when drive-wheel passing over Pavement 2 Fig. 13. Changes in torque on the half shaft when drive-wheel passing over Pavement 3 Fig. 14. Changes in torque on the half shaft when drive-wheel passing over Pavement 4 the rain flow method was used (pagoda roof method), [32]. After the processing of the diagrams of the torque change on the half shaft when a drive-wheel passes over a random pavement irregularity (shown in Figs. 11 to 14) which are divided into twenty classes of amplitude range, Figs. 15 to 22 show the histograms of the vehicle half-shaft torsion stress and the histograms of the relative number of occurrences of amplitude load. The case of the stress of the analysis is presented for all the pavements noted in Table 2. In the Figs. 15 to 22, a significant change in the half-shaft torque can be noticed, depending on the pavement characteristics. 123 Fig. 15. Half shaft load, т [N/mm 2] Histogram of the half-shaft torsion torque stress when drive-wheel passing over Pavement 1 Half shaft load, г [N/mm2] Fig. 17. Histogram of the half-shaft torsion torque stress when drive-wheel passes over Pavement 2 no 100 90 о 70 60 0 ^ 50 Ц g 40 1 30 £ 20 10 U 15 Load class [N/mm2] Fig. 16. The relative number of occurrences of the amplitude of the torsion torque of the half shaft when drive-wheel passing over Pavement 1 Load class [N/mm2] Fig. 18. The relative number of occurrences of the amplitude of the torsion torque of the half shaft when drive-wheel passes over Pavement 2 For an analysis of the half-shaft lifespan, a half-shaft material (Č.1530) with the characteristics Td(-1) = 180 N/mm2 and tD(0)=230 N/mm2 was used as an example. The calculation was done according to the Sorensen-Kogaev method. The following values were used in the calculation: stress Tmin = 0.6 X td, base cycle number 5.6 x 106 and Wöhler curve slope exponent 2.6. The half-shaft diameter is 22.8 mm. Through the lifespan analysis and according to the load characteristics obtained by the simulated movement of the vehicle drive-wheel across the four types of pavement, the half-shaft lifespans for Pavements 3 and 4 were obtained. In the conditions of differing pavement types, the lifespan reduction amounts to 88 % in comparison with Pavements 3 and 4. The lifespan values are not calculated for the first two categories of a pavement because the method used for the calculating of a lifespan uses only the stress values exceeding 60 % of td. Based on the results of the analysis, it is clear that the pavement characteristics, in terms of a height of irregularities, significantly affect a vehicle half-shaft stress and lifespan in a way that the lifespan decreases with increasing pavement irregularities. As for the input data for calculating the half-shaft service life and the calculated service life, they are given only as an example for a general conclusion on the stated influences, and not with the aim of dwelling on a specific half shaft. 6 CONCLUSIONS This paper gives a description of the set mathematical model of the dynamic behaviour of a drive and support system in a vehicle, thus enabling the execution of the analysis of the lifespan of the vehicle drive system elements. In addition to the mathematical model, this paper presents the conducted analysis of the impact of surface irregularities on the load of a vehicle drive-wheel half shaft, with a given example of the analysis of lifespan of the half shaft as a part of the drive system. 20 40 60 80 100 Half shaft load, r [N/mm2] Fig. 19. Histogram of the half-shaft torsion torque stress when drive-wheel passes over Pavement 3 20 40 60 80 100 Half shaft load, г [N/mm2] Fig. 21. Histogram of the half-shaft torsion torque stress when drive-wheel passes over Pavement 4 rki.....i.........I........I........I.........I........i.........I........I i m ;........; 160l I ; É||r„J 10 15 20 25 30 35 40 45 Load class [N/mm2] Fig. 20. The relative number of occurrences of the amplitude of the torsion torque of the half shaft when drive-wheel passes over Pavement 3 no Load class [N/mm2] Fig. 22. The relative number of occurrences of the amplitude of the torsion torque of the half shaft when drive-wheel passes over Pavement 4 Based on the above content, it is possible to formulate the following conclusions: • when the wheel passes over pavement with irregularities, forces of an extremely dynamic character are generated, with directions different from those in a vehicle moving on a flat pavement, • the pavement surface irregularities have a significant effect on the interaction of the wheel with the road, • the forces generated at the contact point between the wheel and the pavement depend on the characteristics of pavement irregularities. The results and conclusions concerning the impact of pavement irregularities on the load and a lifespan of the vehicle half shaft require the choice of parameters of a vehicle systems not only in terms of their impact on driving comfort, handling and stability, but also, from the standpoint of impact on load and the lifespan of the drive train elements. This is especially vital for the designs of special vehicles spending their service life primarily on roads with bad characteristics, and this would further enhance the process of the development of design and control for such systems. Based on the experiments and the simulation described in this paper, it can be concluded that the characteristics of a wheel force and a change of a suspension position are very complex in cases in which a vehicle moves along an uneven road surface. The presented model has an advantage of a simple adjustment of a vehicle configuration to different vehicle types and to the different configurations of their systems. 7 REFERENCES [1] Blundell, M., Harty, D. (2004). The Multibody Systems Approach to Vehicle Dynamics. Elsevier ButterworthHeinemann, Oxford 125 [2] Genta, G. (1997). Motor Vehicle Dynamics - Modeling and Simulation. World Scientific Publishing, Singapore, D0l:10.1142/9789812819765_fmatter. [3] Gillespie, T.D. (1992). Fundamentals of Vehicle Dynamics, SAE, Warrendale, DOI:10.4271/R-114. [4] Kiencke, U., Nielsen, L. (2000). Automotive Control Systems -For Engine, Driveline And Vehicle. SAE, Springer-Verlag Berlin, Heidelberg, New York [5] Kulakowski, B.T. (ed.) (1994). Vehicle-Road Interaction. ASTM, Philadelphia, D0l:10.1520/stp1225-eb. [6] Beleckii, A.V. (2008). Modelling of Profile of Road Basement in Tasks of Analysis of Dynamic of Powertrain of Wheeled Machines. Civil Engineering, Road Machines and Technic, Facilities for Civil Engineering. from http://sdm.str-t.ru/ publics/page_6/ accessed on 2014-12-01. (in Russian) [7] Ahlawat, R., Jiang, S., Mendoza, D., Smith, H.M. (2012). On emulating engine and vehicle transient loads for transmission-in-the-loop experiments. Mechatronics, vol. 22, no. 7, 989996, D0l:10.1016/j.mechatronics.2012.07.001. [8] Caruntu, C.F., Lazar, M., Gielen, R.H., van den Bosch, P.P.J. (2013). Lyapunov based predictive control of vehicle drivetrains over CAN. Control Engineering Practice, vol. 21 no. 12, p. 1884-1898, D0l:10.1016/j.conengprac.2012.05.012. [9] Wallenowitz, H. (2006). Automotive Engineering II - Tire Technology - Simulation und Testing. Automotive Institute, Aachen. [10] Dixon, J. (1996). Tires, Suspension and Handling, 2nd edition, SAE, London, D0l:10.4271/R-168. [11] Lajqi, S., Pehan, S. (2012). Designs and optimizations of active and semi-active non-linear suspension systems for a terrain vehicle. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 12, p. 732-743, D0l:10.5545/sv-jme.2012.776. [12] Popović, V., Vasić, B., Petrovič, M., Mitić, S. (2011). System approach to vehicle suspension system control in CAE environment. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 2, p. 100-109, D0l:10.5545/sv-jme.2009.018. [13] Pacejka, H.B. (2006). Tyre and Vehicle Dynamics, 2nd edition. Butterworth-Heinemann, Elsevier, Oxford. [14] Cebon, D. (1993). Interaction between heavy vehicles and roads. SAE Technical Paper 930001, D0l:10.4271/930001. [15] Cole, D.J., Cebon, D. (1996). Truck tyres, suspension design and road damage. Proceedings of the International Rubber Conference, Manchester. [16] Abdullah, S., Choi, J.C. (2006). Using fatigue data editing approach for analysis cycle sequence effects in variable amplitude road loadings. Journal - The Institution of Engineers, vol. 67, no. 2, p. 47-54. [17] Sakai, E., Tominaga, H., Kato, Y., Yamamoto, N., Matsuki, H. (2005). Development of technique for estimating suspension input force during endurance run on rough road. Mitsubishi motors - Technical review, No. 17, p. 45-48, JSAE Paper Number: 20050500. [18] Huhtala, M., Halonen, P., Sikiö, J. (1998). Measurements on dynamic effects of dual and wide base single tyres. COST- 344 Report, Technical Research Centre of Finland, VTT-Yhdyskuntatekniika, Espoo. [19] Sener, A.S. (2012). Determination of vehicle components fatigue life based on FEA method and experimental analysis. International Journal of Electronics, Mechanical and Mechatronics Engineering, vol. 2, no. 1, p. 133-145. [20] Strelkov, M.N. (2007). Estimational investigation of low frequency correlatively connected oscillations of powertrain and suspension of light duty vehicle. Herald IžGTU, no. 2. (in Russian) [21] Lee, S.-H.; Lee, J.-H., Lee; Goo, S.-H.; Cho, Y.-C.; Cho, H.-Y. (2009). An evaluation of relative damage to the powertrain system in tracked vehicles. Sensors, vol. 9, no. 3, 1845-1859, D0l:10.3390/s90301845. [22] Pesterev, A.V., Bergman, L.A., Tan, C.A. (2002). Pothole-induced contact forces in a simple vehicle model. Journal of Sound and Vibration, vol. 256, no. 3, p. 565-572, D0l:10.1006/jsvi.2001.4220. [23] Foulard, S., Rinderknecht, S., Ichchou, M., Perret-Liaudet, J. (2015). Automotive drivetrain model for transmission damage prediction. Mechatronics, vol. 30, p. 27-54, D0l:10.1016/j. mechatronics.2015.06.008. [24] Foulard, S., Ichchou, M., Rinderknecht, S., Perret-Liaudet, J. (2015). Online and real-time monitoring system for remaining service life estimation of automotive transmissions -Application to a manual transmission. Mechatronics, vol. 30, p. 140-157, D0l:10.1016/j.mechatronics.2015.06.013. [25] Willmerding, G., Häckh, J., Berthold, A. (2000). Driving cycle, load and fatigue life predictions based on measured route data. Society of Automotive Engineers, paper 01ATT120. [26] Seherr-Thoss, H. C., Schmelz, F., Aucktor, E. (2006). Universal Joints and Driveshafts, Analysis, Design, Applications. 2nd ed. Springer-Verlag, Berlin, Heidelberg. [27] Lechner, G., Naunheimer, H. (1999). Automotive Transmissions - Fundamentals, Selection, Design and Application. SpringerVerlag Berlin, Heidelberg. [28] Durković, R. (2015). Power Transmission. University of Montenegro, Faculty of Mechanical Engineering, Podgorica. [29] Jovanović, J. (2009). Influence of damage hypothesis on reliability evaluation of vehicle transmission elements. 13th International Research/Expert Conference "Trends in the Development of Machinery and Associated Technology', Hammamet. [30] Hačaturov, A.A., Afanasjev, V.L., Vasiljev, V.S., Goljdin, G.V., Dodonov, B.M., Žigarev, V.P., Koljcov, V.l., Jurik, V.S., Jakovlev, E.I. (1976). Dynamics of System Road-Pneumatic-Vehicle-Operator. Engineering, Moscow. (in Russian) [31] Simović, S. (2011). The Influence of the Suspension System on the Loading and Lifetime of the Vehicle Drivetrain. PhD thesis, Faculty of Engineering, Kragujevac. (in Serbian) [32] Niestony, A. (2009). Determination of fragments of multiaxial service loading strongly influencing the fatigue of machine components. Mechanical Systems and Signal Processing, vol. 23, no. 8, p. 2712-2121, D0l:10.1016/j.ymssp.2009.05.010. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)2, 127-136 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2837 Original Scientific Paper Design and Analysis of Double-side Meshing and Dual-phase Driving Timing Silent Chain System Yabing Cheng* - Shuaibing Yin - Xiaopeng Wang - Lichi An - Huan Liu Jilin University, School of Mechanical Science and Engineering, China Based on the structure of automotive engines and the layout of their timing system, the timing silent chain system including the overall layout, the structure of link plate and sprockets, and chain length are designed. The design method of double-side meshing and dual-phase driving timing silent chain system is presented. A dynamic analysis model is built, and the fluctuation of the tension sprocket, and the vibration of the system, and the rotation speed stability of exhaust camshaft sprocket, and the transmission error of the intake camshaft sprocket and exhaust camshaft sprocket, and the contact force between the link plate and other components are analysed. Furthermore, dynamic characteristics of the system are contrasted with a single-phase timing silent chain system. Analysis results show that the double-side meshing and dual-phase driving timing silent chain system has a distinct advantage in reducing vibration and enhancing driving stability. The system designed in this paper can meet the requirements of engine timing system, which further verifies the scientific validity and effectiveness of the design method of the double-side meshing and dual-phase driving timing silent chain system proposed in this paper. Keywords: design method, double-side meshing, dual-phase driving, timing silent chain system, dynamic characteristics Highlights • Established the geometric model of link plate and sprockets. • Presented systematic design method of double-side meshing and dual-phase driving timing silent chain system. • Designed a double-side meshing and dual-phase driving timing silent chain system. • Studied the dynamic characteristics of the system. • The double-side meshing and dual-phase driving timing silent chain system have better performance than a single-phase driving timing silent chain system. 0 INTRODUCTION The engine timing system is a core component to ensure normal operation of automotive engines. Recently, with the rapid development of the auto industry, timing chains have increasingly been used in engine-timing systems instead of timing belts or gear transmissions. However, only a few enterprises, such as BorgWarner and DID, own most of the key technologies and patents about new type automobile engine timing chains. Because of business secrecy, only a few references are available. Thus, many enterprises encounter technical barriers in the production of high-quality timing chains. Research on the chain-driving system has been conducted by many scholars and experts. Based on the link plate-sprocket-hob meshing mechanism, Meng et al. [1] introduced the choice and calculation method for the chain, analysed the working condition and the main loss of effectiveness form of automobile chain, proposed multi-variant variation design method, and examined the effects of the testing temperature and fluctuant speed on the wear resistance of the new silent chain [2] to [3]. Chintien et al. [4] investigated the relative motion between the sprocket and a link plate from the tight span stage to the seated stage. Xu et al. [5] developed a mathematical model to calculate the dynamic response of a roller chain drive working at a constant or variable speed condition. Troedsson and Vedmar [6], researched the oscillations in a chain drive under different speeds using dynamic analysis. Pereira et al. [7], built chain drive automatic multibody models from a minimal set of data, and aimed to overcome the difficulty of building manually complex models of chain drives. Wang et al. [8] designed a new sprocket tooth profile to reduce the polygonal action and meshing impact. Zhang [9] analysed the transverse vibration of an axially moving silent chain by using multi-body dynamics. Xue et al. [10] to [11], researched the design and manufacture of the involute sprocket through CAD software, and studied the engagement theory of the involute sprocket by the engagement analysis between the silent chain and involute sprocket. Pedersen and Hansen et al. [12] and Pedersen [13] developed a dynamic model of a roller chain drive including the impact of guide bars, the inner friction and polygonal action of roller chains. Cheng et al. [14] proposed a design method of the dual-phase Hy-Vo silent chain transmission and studied the polygon effect through mathematical analysis. Recently, with the increase of individualized needs, the research on timing silent chains has been varied, including meshing mechanism variation, *Corr. Author's Address: Jilin University, School of Mechanical Science and Engineering, Renmin street 5988 , Changchun, China. chengyb@jlu.edu.cn 127 shape variation, and parameter variation. Meshing mechanism variation includes outer meshing, inner-outer compound meshing, and outer meshing with inner-outer compound meshing. Shape variation includes pinhole shape variation, shape variation of the pin cross section, and guide plate shape variation. Parameter variation includes the variation of positioning offset angle, the apothem, the plate tooth profile angle, and sprocket pressure angle. However, all the silent chains researched in the above papers are single-side meshing chains based on a single-phase chain driving system. In this paper, a double-side meshing and dual-phase driving timing silent chain system are designed, the system is established, and the design method of the system is put forward. Through dynamic simulation, the rationality and scientific validity of the design method of the double-side meshing and dual-phase driving timing silent chain system are verified, and the advantages of the system are revealed. 1 ESTABLISHMENT OF THE TIMING SILENT CHAIN SYSTEM STRUCTURE As shown in Fig. 1, the timing silent chain system designed in this paper is composed of a crankshaft sprocket, intake camshaft sprocket, exhaust camshaft sprocket, sprocket tension group, idle gear group A, idle gear group B, and timing chains. The crankshaft sprocket, intake camshaft sprocket, and exhaust camshaft sprocket are dual-phase sprockets. Fig. 1. The structure of timing silent chain system The structure of the dual-phase sprocket is shown in Fig. 2. It is composed of two connected identical sprockets, and the two sprockets have a phase difference of n/z, where z is the tooth number of the sprocket. The tension sprocket group, idle gear group A, and idle gear group B have the same structure, and they all consist of a pair of coaxial single-phase sprockets, which can rotate without interference. The structure of the single-phase sprocket is shown in Fig. 3, and the assembly diagram of the idle sprocket groups is shown in Fig. 4. The crankshaft sprocket is the driving wheel, and has 19 teeth. The intake camshaft sprocket and the exhaust camshaft sprocket are driven wheels, and it has 38 teeth. There are two timing chains with the same length in the system, and each chain meshes with one sprocket of the sprocket groups. This design ensures that the sprocket tension group and idle gear groups can work more efficiently and minimize the vibration of the system. b) Fig. 2. Schematic diagram of the dual-phase sprocket; a) Front view, b) side view Fig. 3. Schematic diagram of the single-phase sprocket The coordinate system is established as shown in Fig. 1. The crankshaft sprocket centre is regarded as the coordinate origin, and the horizontal direction is regarded as the X-axis, and the vertical direction is regarded as the Y-axis. The intake camshaft sprocket and exhaust camshaft sprocket are symmetrical about the Y-axis. The coordinate of the exhaust camshaft sprocket centre is (-75,160), and the coordinate of intake camshaft sprocket centre is (75,160). To ensure the stability of the driving, the centre distance between the driving sprocket and driven sprocket is usually 30 to 50 times that of the chain pitch. In order to reduce engine vibration and noise, sprocket tension group is installed on the loose side, and idle gear group B is installed on the tight side. Fig. 4. Assembly diagram of sprocket tension group and Idle sprocket groups 2 GEOMETRIC MODEL ESTABLISHMENT OF THE DOUBLE-SIDE MESHING AND DUAL-PHASE DRIVING TIMING SILENT CHAIN SYSTEM 2.1 Geometric Model Establishment of the Link Plate Fig. 5 shows the structure of the link plate. The link plate is an inner and outer compound mechanism that can further reduce the fluctuation of the system [15]. Supposing that A is the hole pitch of the link plate, f is the benchmark apothem between the circle centre O and link plate outside linear profile, 2a is the tooth profile angle of the link plate, h is the distance between the centreline and the crotch of the link plate, r is the curvature radius of inner tooth profile, and Oj (xby1) is the central coordinate of it. When the chain is straightened, the medial profile overhang of the link plate is S. The crankshaft sprocket, tension sprockets, and idle gears have the same number of teeth. The basic pitch of the double-side meshing chain is p=8 mm, the gap between the pinhole and the pin is A = 0.04 mm, A = p - A = 7.96 mm, S = 0.217 mm, a = 30°, r = 13.37 mm, h = 1.4 mm. When the link plate is a broad waist link plate, f= 0.40, p = 3.2 mm. The coordinate of the inside tooth profile curvature centre is calculated as: r-f-s ■ a X =---sin ( A + cos A r - f -S (1 y =-cos (A - 60° 60° (1) cos A where X = arctan (x1 /y1) - 60° = 10°. According to Eq. (1), x1 = 9.513 mm, y1 = 3.463 mm can be obtained. At the onset of the meshing, the inside convex flank of the link plate contacts the sprocket first. With the relative rotation of the link plate with the adjacent link plate, the engagement point transfers from the inside flank of the link plate to the outside flank of the adjacent link plate, and the inner-outer compound meshing is completed. Because of the symmetry of the link plate, the meshing process is exactly the same when the two sides of the link plate mesh with the sprockets. This meshing mechanism reduces the meshing impact between link plates and sprockets, relieves vibration, and improves the wear resistance of the system. Fig. 5. The structure of the double-side meshing link plate 2.2 Establishment of the Geometric Model of the Sprockets Fig. 6 shows the meshing diagram of link plates with a dual-phase sprocket. Involute tooth sprockets are designed to reduce the impact between the link plates and sprockets. Moreover, it will lead to a lower impact velocity and a higher wear resistance of the sprockets. Supposing p0 is the initial pitch of the silent chain, f0 is the initial apothem. p2 is the normal pitch of hob, and normal tooth form angle of the hob is a2 which is equal to a. a is the phase difference of the sprocket, dj is the pitch circle diameter of the sprocket, and d is the reference diameter of the sprocket. x is the modification coefficient of the hob, straight line k-k is the reference line of the hob before modification, and straight line e-e is the reference line of the hob after modification. Therefore, x is a minus, and the distance between e-e and k-k is -xm2. As shown in the diagram, dotted lines stand for one phase of the system, while full lines stand for another phase. dR is the diameter of the gauge pin. 2.2.1 The Parameters of Crankshaft Sprocket The tooth number of the crankshaft sprocket is z = 19. When the tooth number of the sprocket is too small, a slight pressure angle of the sprocket will make the Fig. 6. The design system of the double-side meshing and dual-phase driving silent link plate- sprocket-hob tooth root fragile. Thus, the sprocket pressure angle is set as 36°. The normal pitch of the hob is p2 = p0 = 7.455 mm, and the normal tooth form angle of the hob is a2 =a =30°. The pitch of the crankshaft sprocket is p1 = (p2cos a2) / cos a1 = 7.98 mm. The pitch increment is Ap =p-p0 = 0.545 mm. When a=30°, the side-heart distance increment is Af = [cot (n / z)-л/з]Ap /4 = 0.581mm . The initialized side-heart distance is f0 =f- Af= 2.619 mm. Thus, the modification coefficient of the hob is calculated as: ncot a2 z -2 — + 4 2 2 tan— ^sma2 z (2) The measured distance is calculated as: Mr = -pf sin h п л cos--+ dR, 2 z (3) where dR is 5 mm. According to Eqs. (2) and (3), x=-0.938, MR = 50.138 mm, can be obtained. The phase difference of the sprocket is a=n/z = 9.474°. Since the tooth number of the tension sprockets and the tooth number of the idle sprockets are equal to the tooth number of the crankshaft sprocket, their parameters are equal to those of the crankshaft sprocket. 2.2.2 The Parameters of Intake Camshaft Sprocket and Exhaust Camshaft Sprocket The parameters of the intake camshaft and exhaust camshaft are equal. Tooth number is z=38, the sprocket pressure angle is a1 = 31.5°, the normal pitch of hob is p2 =p0 = 7.876 mm, and the normal tooth form angle of the hob is a2 =a =30°. The pitch of camshaft sprocket is p1 = (p2 cos a2) / cos a1 = 8 mm. The pitch increment is Ap = 0.124 mm. When a=30° the side-heart distance increment is Af = [cot (n / z ) - л/3 ] Ap /4 = 0.320 mm, f0 = 2.520 mm, x=-0.703, MR = 100.293 mm, a=n/z = 4.737°. 2.3 Position Calculation of Tension Sprocket Group and Idle Gear Groups As shown in Fig. 7, the loose side concave distance C2 is C1, and the tight side concave distance C2 is C1. When the centre distance between the sprockets is small, and the transmission is horizontal or approximately horizontal, C1 = (5% to 8%)a, C2 = (2 % to 5 %)a is. When the centre distance between the sprockets is large, and the transmission is vertical or nearly vertical, C1 = (8 % to 12 %)a, C2 = (5 % to 8 %) a [1]. O (x,y) is the centre of the crankshaft sprocket, 01 (x1,y1) is the centre of the tension sprocket group, 02 (x2,y2) is the centre of the idle sprocket group B, O3 (x3, y3) is the centre of the exhaust camshaft sprocket, O4 (x4,y4) is the centre of the intake camshaft sprocket, and O5 (x5,y5) is the centre of the idle sprocket group A. Pitch circle radiuses of x the sprockets are r, r1, r2, r3, r4, r5 respectively, and r = r1 = r2 = r5 = 24.438 mm, r3 = r4 = 48.438 mm. Fig. 7. The design diagram of the system layout 2.3.1 Position Calculation of Tension Sprocket Group As shown in Fig. 7, the coordinate of O is (0, 0), and the coordinate of O1 is (-75, 160). It can be ascertained that the centre distance between the crankshaft sprocket and the exhaust camshaft sprocket is a1 = yj(x3 - x)2 + (y3 - y)2 = 176.706 mm, and C1 = (8 % to 12 %)a1. According to Fig. 7, an equation set can be obtained: b=^/((-x) +~by\—yf d =V(x3 -xi)2 +(Уз -yi)2 sin ßi r3 - r (4) cos ß2 = ai + b12 - di2 2a1b1 ß4 =ß2 =ß3 -ßi C1 = r1-b1 sin ß4+r According to Eq. (4), O1 (-75, 50) and C1 = 17.538 mm can be obtained to satisfy the equation set. 2.3.2 Position Calculation of Idle Sprocket Group B As shown in Fig. 7, the coordinate of O is (0, 0), and the coordinate of O4 is (75, 160). It can be ascertained that the centre distance between crankshaft sprocket sprocket = 176.706 mm , is and and intake_camshaft a2 =V(X4 - X)2 + ( - 7)2 C2 = (5 % to 8 %)a2. According to Fig. 7, an equation set can be obtained: b2 =^j(x2 - X)2 + ( - y )2 d2 =\I(X4 - X2 )2 +(У4 - У2 )2 • a r4 - r sin ß5 = —- (5) C0s ß6 = fl22 + - 2 2 '2 d2 2a2b2 ßs = ßl =ß6 -ß5 C2 = r2-b2 sinßs+r According to Eq. (5), O2 (75, 50) and C2 = 12.880 mm C2 = 12.880 mm can be obtained to satisfy the equation set 29.45°. 2.3.3 Position Calculation of Idle Sprocket Group A As shown in Fig. 7, the centre distance between the exhaust camshaft sprocket and intake camshaft sprocket is a3 =^(x4 - x3 )2 + (y4 - y3 )2 = 150 mm . Owing to the effect of gravity, the chain between the exhaust camshaft sprocket and intake camshaft sprocket hangs down badly. Therefore, C3 should be large enough. The coordinate of point O5 can be set as (0, 195), and C3 = 37.74 mm can be obtained. 2.4 The Length Calculation of the Silent Chain As shown in Fig. 8, the length of the silent chain is composed of the lengths of six arcs and the lengths of six straight lines. According to Fig. 8, an equation set can be obtained: fb +r )2 f2 =>/-(r, + Г, )2 ß9 =ßl0 -ßl. ßl2 ßl0 = arccos f j 2 12 2 Л b. + d.2 - a.z 2b.d. (6) ■■ arccos r + r. V bi J ßl2 = arccos d. а 2 a ß,3 Fig. 8. The diagram of chain length calculation According to Eq. (6), ß9 = 23.88°, f = 70.98 mm, f2 = 82.67 mm, can be obtained. Similarly, the values of f3,f4,f>,f6,fi3,fi4,fi5,fi6, fi7 can be obtained. Since r=r1 = r2 = r5, r3 = r4, it can be ascertained that the length of the chain can be calculated as: S = (ß9 + ßl3 + ßl4 + ßiv)x r + (ßl5 + ßie)x r3 + +( fi + f + f + f4 + f5 + fe)- (7) Moreover, the link number is calculated as: L = S / P. (8) According to Eqs. (7) and (8), L = 96.1 can be obtained. Considering the gaps between adjacent links, L should be larger than the exact value and be an even number. By adjusting the position of the sprocket tension group and the initial phase of all the sprockets, it can be ascertained that the link number of the system is L = 98. 3 DYNAMIC SIMULATION OF THE TIMING SILENT CHAIN SYSTEM As shown in Fig. 9, based on the design parameters above, a dynamic model is established based on RecurDyn dynamic simulation software in which the decoupling algorithm of implicit numerical integration is used. In the system, a rigid body contact model is used, and all the components in the system are rigid bodies whose stiffness is close to infinity. The dynamic model is simplified in this paper, taking no account of the chain moving along the sprocket axial direction and the friction of the system. The timing silent chain is simplified to 1X1 form, and the pin length is more than 2 times the thickness of the link plate. As for a 4-cylinder and 4-stroke engine, the load torque on the exhaust camshaft sprocket is a sine function of the rotation angle of the camshaft sprocket, and it goes through four cycles while the exhaust camshaft sprocket takes a turn. According to the principles of the engine valve system, intake camshaft sprocket is also a sine function only with a phase difference about exhaust camshaft sprocket, and the phase difference is 1/4 of the period T. Therefore. the functions between load torque, and the rotation angle of camshaft sprockets can be established. It can be ascertained that the period of the load torque is T = 2n/4 = n/2, the frequency is m = 2n/T=4, the phase difference between exhaust camshaft sprocket and exhaust camshaft sprocket is ф = T/4 = n/8. Moreover, for the convenience of research, amplitude A and initial torque B can be valued as A = 10 and B=2. Therefore, it can be obtained that the load torque exerted on the exhaust camshaft sprocket is M1 = (10 x sin(46> +n/8) + 2), and the load torque exerted on the intake camshaft sprocket is M2 = (10 x sin(46>) + 2). As shown in Fig. 10, curve 1 represents M1 and curve 2 represents M2, where в is the rotation angle of the camshaft sprocket. Idle gear Fig. 9. The simulation model of double-side meshing and dualphase driving timing silent chain system After simulation, the result can be obtained through Plot Result in RecurDyn. The fluctuation of tension sprocket, the vibration of the chain system, the rotation speed stability of exhaust camshaft sprocket, the transmission ratio of the intake camshaft sprocket and exhaust camshaft sprocket, and contact force between link plate and other components are analysed when the crankshaft sprocket rotation speed is n = 2000 r/min. As shown in Fig. 11, the input driving on the crankshaft sprocket is step (time, 0, 0, 0.05, -66.7*n). Furthermore, the simulation results of the dual-phase silent chain system and single-phase silent chain system are compared. Tension sprockets move along a straight line, and the angle between the line and X axis is 36.2°. Thus, the fluctuation of tension sprockets can be obtained by measuring the fluctuation along the X axis. Due to the existence of the assembly gap, the fluctuation is serious during the initial stage of the running. Take the period of 0.11 s to 0.18 s to study the fluctuation when the system is approximately steady, and the fluctuation curve is shown in Fig. 12. 0 45 90 135 180 225 270 315 360 0[°] Fig. 10. Load torque on the camshaft sprockets -250- 0 0.09 0.18 0.27 0.36 t M Fig. 11. Input driving on the crankshaft sprocket 3.1 Fluctuation Analysis of Tension Sprocket As shown in Fig. 12, I is the fluctuation of singlephase silent chain system, and II is the fluctuation of the dual-phase silent chain system. The minimum value and maximum value of curve I is A (-60.131) and B (-60.058). The minimum value and maximum value of curve II are C (-60.119) and D (-60.072). yA, yB, yC, yD are respectively the y-direction ordinates of point A, B, C, and D. Thus, the maximum fluctuation quantity can be obtained: d = (yB -yA) / cos 36.2° = 0.090 mm, d2 = (yD -yC) / cos 36.2° = 0.058 mm. ■60.04 -60.141_ 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.1S t[s] Fig. 12. Tension sprocket fluctuation along X-axis It can be observed that the fluctuation of the dual-phase silent chain system is small. Moreover, the fluctuation of the dual-phase silent chain system is smaller than that of the single-phase silent chain system. Furthermore, it is obvious that the fluctuation of curve I is more drastic than the fluctuation of curve II. So it can be concluded that the transmission of double-side meshing and the dual-phase driving system is steady and credible. 3.2 Link Vibration Analysis Because of the influence of the chain polygon effect, the impact between the link plates and the sprockets, and the gravity of the chain, transverse vibration, and longitudinal vibration are unavoidable in the process of transmission. In this paper, the link vibration is analysed through researching link trajectory and the rotational speed of a link plate. In Fig. 13, the link trajectory is shown. In Fig. 14, the rotational speed of a link plate is shown, where T1 and T2 are two periods of the running. 200 150 J 100 > 50 0 -50 -150 -100 -50 0 50 100 150 X [mm] Fig. 13. The link trajectory From Fig. 13, it can be seen that the link trajectory is smooth and steady, and the vibration amplitude is very small. From Fig. 14, it can be seen that at the beginning of the running, the fluctuation of the rotational speed is bigger and bigger, which is caused by the assembly gap of the system and the effect of the sprocket tension group. During link plate meshing with sprockets, the fluctuation is small, which indicates, according to [16], that the meshing is smooth. In summary, double-side meshing and dualphase driving timing silent chain system can work smoothly and has little vibration, which can enhance reliability and decrease the noise of the system. ^ = Дю / ю 2, (9) Fig. 14. The link plate rotational speed 3.3 Rotational Speed Analysis of Exhaust Camshaft Sprocket Since the rotational speed of the crankshaft sprocket is constant, the stability of the exhaust camshaft sprocket becomes an important reliability indicator of a silent chain system. The theoretical rotational speed of the exhaust camshaft sprocket is ю2 = ю1 / n = 104.72 rad/s, where ю1 is the rotational speed of the crankshaft sprocket, and n is the transmission ratio of the exhaust camshaft sprocket and the crankshaft sprocket. As shown in Fig. 15, I is the rotational speed curve of the single-phase silent chain system, and II is the rotational speed curve of a dual-phase silent chain system, and III is the theoretical rotational speed curve. Fig. 15. The rotational speed of the exhaust camshaft sprocket Rotational speed stability index ц is always used to measure the stability of sprocket rotation. Moreover, ц is calculated as: where Дю is the fluctuation quantity of the rotational speed. From Fig. 15, it can be seen that the tendencies of curve I and II are roughly the same. According to Eq. (9), it can be obtained that the rotational speed stability index of a single-phase silent chain system is 0.66 %, the rotational speed stability index of a dual-phase silent chain system is 0.38 %>, and both of them can meet the design requirements of the system, but the dual-phase silent chain system has better performance than the single-phase silent chain system. 3.4 Transmission Error Analysis Because of the influence of the chain polygon effect and impact load, the transmission error between exhaust camshaft sprocket and intake camshaft sprocket is unavoidable. As shown in Fig. 16, I stands for the transmission ratio of a single-phase silent chain system, II stands for the transmission ratio of a dual-phase silent chain system, and III stands for the theoretical transmission ratio. Since the intake camshaft sprocket and the exhaust camshaft sprocket have same tooth number, the theoretical transmission ratio is a constant one. Fig. 16. The transmission ratio of intake camshaft sprocket and exhaust camshaft sprocket From Fig. 16 it can be seen that the tendencies of curve I and curve II fluctuate around straight line III and they are similar. The transmission error of single-phase silent chain system is 0.80 %>, while the transmission error of the dual-phase silent chain system is 0.43 %>. It can be concluded that the transmission errors of the two systems are small, and both of them can meet the design requirements of the system. However, the dual-phase silent chain system can work more reliably than the single-phase silent chain system can. 3.5 Contact Force Analysis The contact force between link plate and other components is shown in Fig. 17. Stage I represents the contact force when the link plate contacts the crankshaft sprocket. Stage II represents the contact force when the link plate contacts the sprocket tension group. Stagelll represents the contact force when the link plate contacts the exhaust camshaft sprocket. Staged represents the contact force when the link plate contacts the idle gear group A. Stage V represents the contact force when the link plate contacts the intake camshaft sprocket. Staged represents the contact force when the link plate contacts the idle gear group B. Fig. 17. The link contact force From Fig. 17, it can be seen that the contact force between the link plate and other components is small, which indicates that the double-side meshing and dual-phase driving timing silent chain system can reduce action force and friction between link plates and sprockets and prolong the life of the system. From Stages I, III and V, it can be seen that the contact force is large at the beginning and end of the meshing. This is because that at the beginning of the meshing between the link plate and sprocket, the majority of the load is imposed on the link plate. With the progress of the meshing, the adjacent link plate starts to engage with the sprocket and undertakes an increasing load. Thus, the contact force becomes smaller and smaller. At the end of the meshing, with the former link plate dropping out of the engagement, an increasing load will be imposed on the next link plate. 4 CONCLUSIONS (1) Based on engine structure, a new type of doubleside meshing and dual-phase driving timing silent chain system is designed. Furthermore, a systematic design method of the system is proposed, which lays the theoretical foundation for the design of the engine timing silent chain system. (2) The fluctuation of the tension sprocket, the vibration of the chain system, the rotation speed stability of the exhaust camshaft sprocket, the transmission error of the intake camshaft sprocket and exhaust camshaft sprocket, and the contact force between link plate and other components are analysed. The analysis results show that the system designed in this paper can meet the performance requirements of the engine timing silent chain system. Furthermore, it verifies the scientific validity and the validity of the design method proposed in this paper. (3) By comparing the single-phase timing silent chain system, the double-side meshing and dualphase driving timing silent chain system has less fluctuation of the tension sprocket, and a steadier rotation speed of exhaust camshaft sprocket, and a smaller transmission error of the intake camshaft sprocket and exhaust camshaft sprocket. 5 ACKNOWLEDGEMENTS The research presented in this paper was supported by Jilin Province Science and Technology Development Plan Item (No. 20150204075GX), and National Natural Science Foundation of China (No. 51305154). 6 REFERENCES [1] Meng, F., Li, B., Lu, X., Wu, L., Wang, H. (2009). Design method of timing chain system for automotive engine. Journal of Harbin Institute of Technology, vol. 41, no. 5, p. 121-124. (in Chinese) [2] Meng, F., Qu, S., Dong, C. (2012). Multi-variation design method of Hy-Vo silent chain. Transactions of the Chinese Society of Agricultural Machinery, vol. 43, no. 2, p. 230-234, DOI: 10.6041/j.issn.1000-1298.2012.02.044. (in Chinese) [3] Meng, F., Feng, Z., Li, C., Yin, D. (2004). Experimental research on the wear mechanisms and temperature and speed characteristics of a new silent chain. Tribology, vol. 24, no. 6, p. 560-563. (in Chinese) [4] Huang, C., Kosasih, L., Huang, C.C. (2005). The tooth contact analysis of round pin jointed silent chains. ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, vol. 5, p. 605613, DOI:10.1115/DETC2005-84065. [5] Xu, L., Yang, Y., Chang, Z., Liu, J. (2010). Dynamic modeling of a roller chain drive system considering the flexibility of input shaft. Chinese Journal of Mechanical Engineering, vol. 23, no. 3, p. 367-374, DOI:10.3901/CJME.2010.03.367. [6] Troedsson, I., Vedmar, L. (1999). A dynamic analysis of the oscillations in a chain drive. Journal of Mechanical Design, vol. 123, no. 3, p. 395-401, D0l:10.1115/1.1374196. [7] Pereira, C.M., Ambròsio, J.A., Ramalho, A.L. (2010). A methodology for the generation of planar models for multibody chain drives. Multibody System Dynamics, vol. 24, no. 3, p. 303-324, D0I:10.1007/s11044-010-9207-x. [8] Wang, Y., Ji, D.S., Zhan, K. (2013). Modified sprocket tooth profile of roller chain drives. Mechanism and Machine Theory, vol. 70, p. 380-393, DOI:10.1016/j. mechmachtheory.2013.08.006. [9] Zhang, W. (2005). An analysis on vibration of moving silent chain by multi-body dynamics. ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, vol. 6, p. 1707-1712, D0I:10.1115/DETC2005-85311. [10] Xue, Y., Wang, Y., Wang, X. (2006). CAD and machining simulation of involute sprockets. Journal of Wuhan University of Technology, vol. 28, no. 5, p. 105-107. (in Chinese) [11] Xue, Y., Wang, Y., Wang, X. (2007). Engagement theory of silent chain mechanism with involute tooth profile. Journal of Jiangsu University (Natural Science Edition), vol. 28, no. 2, p. 104-107. (in Chinese) [12] Feng, Z., Meng, F., Li, C. (2005). Meshing mechanism and simulation analysis of a new silent chain. Journal of Shanghai JiaotongUniversity, vol. 39, no. 9, p. 1427-1430. (in Chinese) [13] Pedersen, S.L., Hansen, J.M., Ambròsio, J.A.C. (2004). A roller chain drive model including contact with guide-bars. Multibody System Dynamics, vol. 12, no. 3, p. 285-301, D0I:10.1023/ B:MUB0.0000049131.77305.d8. [14] Pedersen, S.L. (2005). Model of contact between rollers and sprockets in chain-drive systems. Archive of Applied Mechanics, vol. 74, no. 7, p. 489-508, DOI:10.1007/s00419-004-0363-4. [15] Cheng, Y., Wang, Y., Li, L., Yin, S., An, L., Wang, X. (2015). Design method of dual phase Hy-Vo silent chain transmission system. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 4, p. 237-244, DOI:10.5545/sv-jme.2014.2318. [16] Cali, M., Sequenzia, G., Oliveri, S.M., Fatuzzo, G. (2015). Meshing angles evaluation of silent chain drive by numerical analysis and experimental test. Meccanica, p. 1-15, DOI:10.1007/S11012-015-0230-0. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)2, 137-144 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2994 Original Scientific Paper WEDM Manufacturing Method for Noncircular Gears, Using CAD/CAM Software César Garcia-Hernandez1* - Rafael-Maria Gella-Marin1 - José-Luis Huertas-Talón1 -Nikolaos Efkolidis1 - Panagiotis Kyratsis2 1 University of Zaragoza, Department of Design and Manufacturing Engineering, Spain 2 Technological Education Institution of Western Macedonia, Department of Mechanical Engineering & Industrial Design, Greece Noncircular gears are used in several technological applications in order to enhance the performance of different mechanical instruments (flow meters, bikes, internal combustion engines, etc.), in order to unify speed in assembly lines and in research. Noncircular gears are typically manufactured by shaving: milling each tooth or by generation. This requires controlling the geometric and kinematic variables in the process. In this research, a method to manufacture elliptical and oval gears using wire electro-discharge machining (WEDM) is presented. This is a continuous procedure, and its performance is not inferior to the previously mentioned methods. Mathematical models for manufacturing elliptical and oval gears are presented, simulations are carried out, and this method is implemented in a WEDM machine, obtaining two pairs of elliptical and oval gears. This method could be useful in the manufacturing of injection moulds or custom-made metallic gears. Finally, a discussion using bibliographic references is presented about the surface finish and the consequences of using WEDM in comparison to other shaving methods which do not involve a material phase change. Keywords: WEDM, noncircular gears, worksheet, CAD/CAM Highlights • Wire Electro Discharge Machining is applied to noncircular gear manufacturing. • The mathematical models for manufacturing elliptical and oval gears are included. • The algorithm of the mathematical models is implemented using MatlabTM. • Simulations are carried out, and the method is implemented in a real WEDM machine. • Two pairs of elliptical and oval gears are obtained applying this method. 0 INTRODUCTION The manufacture of noncircular gears (in which elliptical and oval gears are included) is an interesting research topic from several aspects. It has practical commercial applications, since these gears are used in several mechanisms [1], such as in flow meters, chain transmission systems, transmissions of bike plates [2], oil pumps and even in chain belts used in combustion engines [3] (Fig. 1). Noncircular gears are also studied as a way of developing different types of gears that can meet movement specifications [4] or as a way of experimenting with several machining methods [5]. They have also been used to transmit movement in several different ways, as we can see in different museums [6]. Similarly, this field has been addressed in research. In [7] a literature review is carried out, analysing the analytical methods used to obtain the equations which define the pitch circumference. Furthermore, the methods to obtain these gears are *Corr. Author's Address: University of Zaragoza, Dept. of Design and Manufacturing Engeenering, Campus Rio Ebro, C/ Maria de Luna, 3 - 50018 - Zaragoza, Spain, garcia-hernandez.cesar@unizar.es revised, and an experimental analysis is carried out of the transmission ratio relating to speed and torque. Fig. 1. Flowmeter and oval gears It is important to point out the kinematic differences between noncircular gears. This research focuses on designing and manufacturing elliptical and oval gears. As it has been widely described [8], elliptical gears rotate around their focus. However, this kinematic relation is not the same in oval gears. According to Bloomfield [9], two oval gears mesh when they rotate around their centre. For this reason, this research addresses the development procedures in order to design and manufacture elliptical gears and oval gears. Even though two different methods are addressed to develop such gears, a common flowchart has been developed to organize the process. It is described in Fig. 2. Fig. 2. Flow chart for design and manufacturing According to this flowchart, the designing process is arranged using Matlab™ and an Excel™ macro. Starting from the initial input data, Matlab™ is used to calculate the semi-axis, the pitch ellipse/ oval, and the initial points distributed on it from where the involute arcs define the start of the teeth. Then, a macro is used to calculate the involute arcs that define the teeth profiles. These points define the gear, and they are exported to general purpose CAD software [10] in order to obtain the 3D gear model; afterwards, they are post-processed in a CAM software to obtain the milling trajectories required to manufacture the gear. Elliptical and oval gears may be similar, but their behaviour is not the same, not only from a geometrical aspect but also from a kinematic point of view as described below. In this article, a Cartesian coordinate system is used in the elliptical gear case and a polar coordinate system to develop the oval one. Being a and b, the semi-major and the semi-minor axes respectively, elliptical gears rotate with respect to their focuses and oval gears rotate around their centres. These differences, along with the equations describing the pitch ellipse and the pitch oval, are represented in Figs. 3 and 4, as shown below. Elliptical gears rotating with respect to their focuses (Fig. 3): 2 2 x y 1 --ъ — = 1. 2 J.2 a b (1) Fig. 3. Elliptical gears Oval gears rotating with respect to their focuses (Fig. 4): 2ab (a + b )-(a - b }cos 26 (2) Fig. 4. Oval gears 1 METHODS The procedure presented in this paper approaches the design by approximating and calculating the teeth as in the spur gears, but taking into account the local curvature of the ellipse [11] and the oval. In order to obtain the teeth of the gears, the radius of curvature of Eqs. (3) and (4) can be applied. It is expressed in Cartesian and polar coordinates respectively by [12]: P = 1 + dy d x d2 y d x2 (3) P=- r2 + d r de 2 , o' dr r + 2 ' d e d2r de2 (4) For the ellipse Eq. (1), applying Eq. (3), the radius of curvature in Cartesian coordinates is: P = Г 4 2 ! z,4 21 /2 lay + b X I 3 2 2 2 2 As an example, the extreme values of the radius of curvature, for a = 100 and b = 60, are represented in Fig. 5. Fig. 5. Extreme radii of curvature of an ellipse using Eq. (1) The radius of curvature obtained for the oval Eq. (2), as a result of applying Eq. (4), is a complex expression. It is important to note that the oval Eq. (2) presents points with zero slope. When a <2b its aspect is similar to an ellipse. For greater values of the parameter a, six points with zero slope are obtained, Fig. 6a. In Fig. 6b, the maximum and the minimum radii of curvature are represented for a = 100 and b = 60. Fig. 6. Different ovals; a) with a > 2b, and b) with a = 100 and b = 60 In the intersection point of the symmetry axis of the tooth with the ellipse/oval, the constructive radius of this tooth is considered the radius of curvature of the pitch ellipse/oval in order to generate the tooth, Fig. 7. The following considerations must also be taken into account. The tooth thickness measured in the pitch ellipse/ oval must be the same for all the teeth; otherwise, the teeth would not mesh with the hole between the flanks of consecutive teeth. As the final thickness is always lower than the nominal one, there is a bigger space between consecutive teeth. Fig. 7. The tooth of the ellipse/oval corresponds to the tooth of a circle with the same curvature In a pair of gears, each one has z teeth, and a tooth pitch measured on the ellipse/oval of p=n m (where m is the module), which is the double of the tooth thickness. If the tooth pitch in the pitch ellipse/oval is nm, then the perimeter is ж-m-z, this last value being the pitch ellipse/oval length. In this case, the calculation starts using the pitch ellipse/oval length, obtaining the values of a and b through another condition, such as the eccentricity. In this way, if a gear with module m = 2.5, teeth number z = 53 and eccentricity e = 0.85 is manufactured, the ellipse/oval length is l = к-m-z = n-2.5-53 = 416.2610. In order to obtain simpler values of a and b, another non-conventional module can be applied. Here, two different methods are discussed to obtain the semi-axis a and b. They are not the same because different approaches are used with an ellipse and an oval, as follows. 1.1 Semi-axis a and b, calculation in an Ellipse The ellipse length is obtained applying this integral: a - l = 4^dx2 + dy2. (6) 0 Solving Eq. (6) leads to an elliptic integral of the second order. After evaluating it with Matlab™ and equalizing it to the value of the ellipse perimeter, it is End I v v I End Fig. 8. Numerical method to obtain semi-axis a and b, solving Eq. (7) possible to solve it and to obtain the value of the semi-axis a and b. 1.2 Semi-axis a and b, Calculation in an Oval The oval length is obtained applying the integral: = 4 f + r2 d в. (7) It is not possible to solve this equation directly with Matlab™; numerical methods can be applied in order to evaluate the integral and equalize it to the value of the oval perimeter with an error smaller than the value set by the gear designer. The numerical method to obtain the value of the semi-axis a and b is presented in Fig. 8. Table 1. Uncertainties in the algorithm l Oval perimeter calculated with the input data_ a Major semi-axis b Minor semi-axis E Eccentricity Error between the calculated perimeter by the numerical method and by the input data er Oval polar equation p Oval perimeter calculated by integration Der Error increment between iterations Da Value to decrement or increment the semi-major axis initial value between iterations The uncertainties in the algorithm are summarized in Table 1. In this algorithm, the initial values for the major and the minor semi-axis a and b are considered. The perimeter value is calculated using input data. Then, using the algorithm, the initial value of the major semi-axis a is increased or decreased until the difference between the perimeter calculated using these values and the perimeter obtained from the input data is smaller than a given error. 1.3 Tooth Distribution A starting point for the pitch ellipse/oval (Fig. 9) is chosen, in this case (a, 0), and the points of the pitch ellipse or the pitch oval are calculated. The calculation of the teeth distribution starts with a Matlab™ program and ends with an Excel™ macro. The Matlab™ program determines the pitch ellipse points with a chordal error lower than a given value, Fig. 6. In the environment of point A(x1,y1), point B(x2, y2) is obtained satisfying the condition that the secant AB maintains a maximum distance d to the curve. This distance is obtained by substituting the coordinates T(x3, y3) in the normal equation of the line. a • x3 + b • y3 + c a + h2 (8) The sign of the distance d must also be kept in mind, which indicates the semiplane where the point l 0 r taken as a reference for that distance is placed. In this case, the square of the distance has been evaluated in order not to take into consideration the distance sign. 100—, X \ ) 40 -1 20 -1 io -г 0 -6 О -4 0 -2 0 2 0 4 0 6 0 8 о ь 30 1 .0 1> 100- Fig. 9. Pitch ellipse Fig. 10. Obtaining the ellipse points (a • x3 + b • y3 + c )2 (a2+b2) (9) The point T is the point of tangency of the parallel line to the secant segment AB with the function graph. Therefore, the first derivate in T is equal to the tangent of the a angle. yx=x = У2 - y Xo Xi (10) In this way, the two conditions which will allow solving the two uncertainties (x2,y2) are found. It is not necessary to analyse the sign of the tangent, because the trajectory is known. In order to solve this system of two equations with numerical methods, MATLAB™ makes possible the application of the following command [13]: maple(fsolve({equ1,equ2},{var1,var2},{var1=v1initial.. v1final,var2=v2initial..v2final})) Knowing the tooth thickness, the ellipse arc segments of the calculated length are distributed on the ellipse and the radius of curvature of these points are determined (Fig. 11). Fig. 11. Points to position the flanks and other characteristic points In order to obtain the initial points of the involute arcs, that the ellipse/oval length obtained by Eqs. (6) and (7) lead to elliptical integrals must be taken into account. For this reason, these calculations have been carried out in MATLAB™, but once the points of the pitch ellipse have been obtained, the separation value of the half tooth pitch is calculated, finding the points that have a separation of this value between them. The following equation is used: l = 1- x )2+(y+i- y )2 • (11) When the length is greater than the tooth thickness on the pitch ellipse/oval, an additional point is interpolated, and this additional point becomes a part of the series of points of the linear approximation (Fig. 12). -10 -9 -8 -7 -6 -5 -4 -3 -2 Fig. 12. Starting flank points calculated on the pitch ellipse/oval The absolute error, which is obtained on a perimeter of 416.2589, is smaller than 0.0021 for a chordal error below 0.001. On these points, the involute arcs corresponding to the radius of the local curvature have been placed (Fig. 13). These arcs start from the base ellipse/oval (equivalent to the base circumference of the local radius of curvature), and have been extended beyond 2 the exterior ellipse/oval. The real tooth will be located between the interior and the exterior ellipse/oval. However, it must be taken into account that not all the teeth can reach the bottom of the calculated involute profile, because when there is a low number of local teeth, where the curvature is smaller, the inner locally equivalent diameter is bigger than the base one. This involute profile, Fig. 13b, has been calculated as in reference [14], but using an Excel™ macro, as it is possible to export the Excel™ data to a general purpose CAD software. In this case, the command "curve by table" of Solid Edge has been used. This process can be automated with a Visual Basic macro [15]. This command also makes it possible to export the points of the addendum, dedendum, and pitch ellipses/ovals. Fig. 13. Points of the different involutes according to the radius of curvature 2 EXPERIMENTAL MANUFACTURING PROCESS Once all these curves have been exported to Solid Edge, the correspondent extrusions and the rounding of the teeth fillet are carried out (Figs. 14 and 15). Flank Involute Fig. 15. Extrusions and fillets All these operations could have been carried out directly with Excel™ or Matlab™, calculating the position of the rounding points for each tooth, as described in [14]. The continuity of the points of the dedendum, addendum, and the flank points, which correspond to a specific tooth profile, provide the path of the WEDM. A wire radius correction plus a gap is applied to these points. However, in this case, general purpose CAD software was used in the process and, later, CAM software was applied. As an example, EdgeCAM™ [16] has been used to make a pair of gears with an eccentricity of 0.75 and 80 teeth are shown. Fig. 16. Points post-processed in a CAM software Fig. 17. Two manufactured gears; a) elliptical; b) oval Fig. 14. Calculated points imported in Solid Edge 3 RESULTS AND DISCUSSION The method presented in this work does not use machines dedicated to a specific manufacturing process, such gear-milling machines, using universalpurpose programmable machines. The technology applied involves a phase shift to obtain the chip removal, which can originate some problems in certain materials. In previous research [17] and [18], a miniature brass gear was manufactured with WEDM, and the advantages of using this method were analysed. Those papers described the obtained results and analysed the superficial parameters of microhardness and microstructure, concluding that this method is appropriate to obtain good quality levels (N5). In the case of steel, surface quality and fatigue studies were carried out using high alloy steel for manufacturing tools, applying WEDM [19], demonstrating that an adequate selection of parameters is needed for an optimum WEDM finishing. Fig. 18. Validation of the algorithm with a coordinate measuring machine In different research [20], the same conclusions were reached with Inconel 718. The superficial finishing was very good, but the conditions of intensity, pulse and discharge time must be taken into consideration in order to obtain a superficial microstructure that does not damage the future integrity of the manufactured part. However, there are previous works [21] in which it was shown that WEDM involves a roughness increase and a decrease of fatigue resistance, in comparison to test pieces manufactured with turning processes. The discrepancy of these results was explained by the difference of twelve years between the machines used for the experimental works ([21] vs. [22]). During this time, manufacturers applied technological improvements to their machines [22]. This is congruent with the conclusions obtained in [23] for Ti-6Al-4V alloy. The algorithm was validated by measuring the eroded teeth with a coordinate measuring machine (Fig. 18). These gears had a tolerance value better than 5 дт, between the maximum and minimum profiles of the manufactured involute, so the manufactured gear offers optimum quality. 4 CONCLUSIONS Even though the application field of elliptical/oval gears is not very extensive, it can be very convenient when a non-constant torque must be transmitted, as well as for research. Simplifying the gear manufacturing process without specific machinery, makes the processes of research, design, and manufacture more accessible to a larger number of people. The example described in this study has been generalized to obtain machining paths, using CAM software. The use of CAM allows manufacturing with several methods: WEDM, water jet cut, laser cut, milling, etc. In this way, a good method is presented so the developer can control the gear design and also the manufacturing methods which ensure a better gear finishing. Although MATLAB™ has been used to solve the equations, it is also possible to apply Excel™, implementing numerical methods with a macro. Finally, with the analysed references about the obtained superficial features in the process, it can be concluded that, with the available machinery and an optimum selection of the conditions which define the process, e.g. intensity, pulse time, etc. (paying attention to the manufacturer advice), the obtained metrological, physical and chemical properties are similar to any other manufacturing process. 5 REFERENCES [1] Hebbale K, Li D, Zhou J, Duan C, Kao C.K., Samie, F., Lee, C., Gonzales, R. (2014). Study of a non-circular gear infinitely variable transmission. ASME 2014 Dynamic Systems and Control Conference, paper no. DSCC2014-6083, p. V003T49A003, D0l:10.1115/dscc2014-6083. [2] Bini, R.R., Dagnese, F. (2012). Noncircular chainrings and pedal to crank interface in cycling: a literature review. Brazilian Journal of Kinantrophometry and Human Performance, vol. 14, no. 4, p. 470-482, D0l:10.5007/1980-0037.2012v14n4p470. [3] Volkswagen (2015). from www.volkspage.net/technik/ssp/ ssp/SSP_337.pdf, accessed on 2015-04-06. [4] Litvin, F.L., Gonzalez-Perez, I., Fuentes, A., Hayasaka, K. (2008). Design and investigation of gear drives with non-circular gears applied for speed variation and generation of functions. Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 45-48, p. 3783-3802, DOI:10.1016/j.cma.2008.03.001. [5] Li, J., Wu, X., Mao, S. (2007). Numerical computing method of noncircular gear tooth profiles generated by shaper cutters. The International Journal of Advanced Manufacturing Technology, vol. 33, no. 11-12, p. 1098-1105, DOI:10.1007/ s00170-006-0560-0. [6] Efstathiou, K., Basiakoulis, A., Efstathiou, M., Anastasiou, M., Seiradakis, J.H. (2012). Determination of the gears geometrical parameters necessary for the construction of an operational model of the Antikythera Mechanism. Mechanism and Machine Theory, vol. 52, p. 219-231, DOI:10.1016/j. mechmachtheory.2012.01.020. [7] Quintero Riaza H.F. (2014). from http://www.tdx.cat/ handle/10803/6418, accessed on 2014-04-06. [8] Litvin, F.L., Fuentes, A. (2004). Gear Geometry and Applied Theory. Cambridge University Press, Cambridge, DOI:10.1017/ CBO9780511547126. [9] Bloomfield, B. (1960). Noncircular Gears. Product Engineering Special Report; p. 159-165. [10] Haba, S.A., Oancea, G. (2015). Digital manufacturing of air-cooled single-cylinder engine block. The International Journal of Advanced Manufacturing Technology, vol. 80, no. 5, p. 747759, DOI:10.1007/s00170-015-7038-x. [11] Fuentes-Aznar, A., Gonzalez-Perez, I., Hayasaka, K. (2009). Noncircular gears: Design and Generation. Cambridge University Press, Cambridge, p. 36-37 and p. 47-52. [12] Mataix, C.M. (1967). Anälisis algébraico e infinitesimal. Editorial Nuevas Graficas, Madrid. (in Spanish) [13] del Rio, J.A.I., Cabezas, J.M.R. (2007). Métodos Numéricos: Teorfa, problemas y präcticas con Matlab™ (eng: Numerical Methods: Theory, problems and practice with Matlab ™). Ed. Piramide, Madrid. (in Spanish) [14] Talón, J.L.H., Ortega, J.C.C., Gómez, C.L., Sancho, E.R., Olmos, E.F. (2010). Manufacture of a spur tooth gear in Ti-6Al-4V alloy by electrical discharge. Computer-Aided Design, vol. 42, no. 3, p. 221-230, DOI:10.1016/j.cad.2009.11.001. [15] Siemens. (2015). Solid Edge "Programming User's Guide" Version 1.0 and 2.0, Plano. [16] Vero Software (2015). from www.edgecam.com, accessed on 2015-04-06. [17] Gupta, K., Jain, N.K. (2014). Comparative study of wire-EDM and hobbing for manufacturing high-quality miniature gears. Materials and Manufacturing Processes, vol. 29, no. 11-12, p. 1470-1476, DOI:10.1080/10426914.2014.941865. [18] Gupta, K., Jain, N.K. (2014). On surface integrity of miniature spur gears manufactured by wire electrical discharge machining. The International Journal of Advanced Manufacturing Technology, vol. 72, no. 9-12, p. 1735-1745, DOI:10.1007/s00170-014-5772-0. [19] Ghanem, F., Fredj, N.B., Sidhom, H., Braham, C. (2011). Effects of finishing processes on the fatigue life improvements of electro-machined surfaces of tool steel. The International Journal of Advanced Manufacturing Technology, vol. 52, no. 5-8, p. 583-595, DOI:10.1007/s00170-010-2751-y. [20] Li, L., Guo, Y.B., Wei, X.T., Li, W. (2013) Surface integrity characteristics in wire-EDM of Inconel 718 at different discharge energy. Procedia CIRP, vol. 6, p. 220-225, DOI:10.1016/j.procir.2013.03.046. [21] Ferrer, C., Carcel, A., Pascual-Guillamón, M., Pérez-Puig, M.A., Segovia, F., Razzaq, K.A. (2003). Investigación de las causas que modifican la resistencia a fatiga en flexión rotativa del proceso de mecanización de electroerosión por hilo en las aleaciones de aluminio AA-2000 (eng. Research of the causes that modify the fatigue resistance in rotated flexion of the process by wire electro discharge machining in the aluminum alloys AA-2000). Anales de mecänica de la fractura (eng. Annals of Fracture Mechanics), vol. 20, p. 201-206. (in Spanish) [22] Ona Electroerosión (2015). from www.ona-electroerosion. com, accessed on 2015-04-06. [23] Mower, T.M. (2014). Degradation of titanium 6Al-4V fatigue strength due to electrical discharge machining. International Journal of Fatigue, vol. 64, p. 84-96, DOI:10.1016/j. ijfatigue.2014.02.018. Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 62, (2016), številka 2 Ljubljana, februar 2016 ISSN 0039-2480 Izhaja mesečno Razširjeni povzetki (extended abstracts) Matej Steinacher, Borut Žužek, Darja Jenko, Primož Mrvar, Franc Zupanič: Izdelava in lastnosti kompozita z infiltrirano magnezijevo zlitino SI 13 Junning Li, Wei Chen, Libo Zhang, Taofeng Wang: Izboljšana kvazidinamična analitična metoda za napovedovanje zdrsavanja kotalnih ležajev v pogojih izjemno majhnih obremenitev in opletanja SI 14 Franc Cimerman, Matej Jarm, Branko Širok, Bogdan Blagojevič: Upoštevanje merilnih pogreškov elektronskih korektorjev pri merjenju prostornin zemeljskega plina SI 15 Ali Majd, Ahmad Ahmadi, Alireza Keramat: Preiskava vplivov nenewtonske tekočine med prehodnimi tokovi v cevovodu SI 16 Sreten Simović, Vladimir Popović, Milanko Damjanović: Analiza vpliva nepravilnosti na vozišču na življenjsko dobo pogonskih polgredi vozila SI 17 Yabing Cheng, Shuaibing Yin, Xiaopeng Wang, Lichi An, Huan Liu: Zasnova in analiza sistema dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige SI 18 César Garria-Hernandez, Rafael-Maria Gella-Marin, José-Luis Huertas-Talón, Nikolaos Efkolidis, Panagiotis Kyratsis: Postopek žične elektroerozijske obdelave nekrožnih zobnikov s programsko opremo CAD/CAM SI 19 Osebne vesti SI 20 Izdelava in lastnosti kompozita z infiltrirano magnezijevo zlitino Matej Steinacher1* - Borut Žužek2 - Darja Jenko2 - Primož Mrvar3 - Franc Zupanič1 1 Univerza v Mariboru, Fakulteta za strojništvo, Slovenija 2 Inštitut za kovinske materiale in tehnologije, Slovenija 3 Univerza v Ljubljani, Naravoslovnotehniška fakulteta, Slovenija V delu je obravnavan kovinsko-keramični material z infiltrirano magnezijevo zlitino (kompozit IPC; ang. interpenetrating phase composite). Preiskovani kompozit IPC je bil izdelan s postopkom gravitacijskega kokilnega litja. Metalografsko smo ga preiskali s svetlobno mikroskopijo, vrstično in presevno elektronsko mikroskopijo, energijsko disperzijsko spektroskopijo in rentgensko difrakcijo, prav tako pa smo opredelili njegove mehanske lastnosti. Osnova kompozita, magnezijeva zlitina AE44, ki je vsebovala 4,94 mas. % Al in 4,42 mas. % kovin redkih zemelj (RE), je bila sestavljena iz primarnih kristalov večkomponentne trdne raztopine a-Mg in intermetalnih faz Al RE Al2RE in Al10RE2Mn7. Utrjevalna sestavina kompozita, keramična pena, je bila sestavljena iz a-Al2O3, a-SiC, ß-SiC in SiO2. Keramična pena je imela odprto primarno in večinoma zaprto sekundarno poroznost. Za tehnologijo izdelave smo uporabili gravitacijsko kokilno litje, kjer smo z ustreznim ulivnim sistemom in opredeljenimi livnimi parametri ter vibriranjem kokile z vstavljeno keramično peno med litjem in strjevanjem izdelali kompozit IPC. Pri tem je magnezijeva zlitina AE44 zapolnila primarne pore, v sekundarne pore pa se je delno infiltrirala, delno pa penetrirala skozi mostičke keramične pene. V mejnih območjih med zlitino AE44 in keramično peno se je pojavila močna reakcija, ki je vplivala na mikrostrukturo nastalega kompozita IPC, zato je bilo največ dela usmerjenega v natančno opredelitev mehanizmov in kinetike kemijskih reakcij v mejnih območjih. Študije mejnih območij med AE44 in keramično peno so bile narejene po različnih temperaturah predgretja kokile z vstavljeno keramično peno, ki so znašale 500 °C, 600 °C in 700 °C. Glavni reakcijski produkti v mejnih območjih med zlitino AE44 in keramično peno ter v penetriranih mostičkih keramične pene so bili MgO, AlSiRE in AlMgSiRE. Najprej je nastal MgO z redukcijo SiO2 in Al2O3 z magnezijem. Nato je na MgO nastala faza AlSiRE, na kateri je kasneje epitaksialno kristalizirala faza AlMgSiRE, katere delež se je z daljšim reakcijskim časom povečeval. Z mikrokemijsko analizo je bila opredeljena njuna kemijska sestava, dognano pa je bilo tudi, da imata fazi tetragonalno kristalno zgradbo s povsem jasno medsebojno kristalografsko orientacijo. Ugotovljeno je bilo, da fazi AlSiRE in AlMgSiRE ne ustrezata nobeni znani fazi. Kompozit IPC je imel v vseh preiskanih stanjih večjo napetost tečenja in modul elastičnosti ter manjšo tlačno trdnost kot zlitina AE44, razen pri predgretju na temperaturo 700 °C. Faza AlMgSiRE je bila manj trda kot faza AlSiRE. Ključne besede: magnezijeva zlitina AE44, keramična pena, kompozit z infiltrirano kovinsko osnovo, gravitacijsko kokilno litje, reakcijski produkt, mehanske lastnosti. *Naslov avtorja za dopisovanje: Univerza v Mariboru, Fakulteta za strojništvo, Smetanova ulica 17, 2000 Maribor, Slovenija, matej.steinacher@gmail.com SI 13 Izboljšana kvazidinamična analitična metoda za napovedovanje zdrsavanja kotalnih ležajev v pogojih izjemno majhnih obremenitev in opletanja Junning Li12* - Wei Chen2 - Libo Zhang2 - Taofeng Wang2 1 Tehniška univerza Xi'an; Univerza Xi'an Jiaotong, Kitajska 2 Univerza Xi'an Jiaotong, Kitajska Zdrsavanje se pogosto pojavlja pri visokohitrostnih in malo obremenjenih kotalnih ležajih (HSLLRB), ko pogonska torna sila ni dovolj velika, da bi premagala upor med kotalnimi elementi in tečino. Opletanje lahko povzroči obrabo in predčasno odpoved ležajev, še posebej pri značilnih ležajih glavnih gredi letalskih motorjev. Ta študija predlaga izboljšano kvazidinamično analitično metodo za napovedovanje zdrsavanja kotalnih ležajev v pogojih izjemno majhnih obremenitev in opletanja, s ciljem ugotavljanja primernih ukrepov proti zdrsavanju ležajev HSLLRB. Predhodna analiza zdrsavanja ležajev HSLLRB ob upoštevanju opletanja bazira na poenostavljenem empiričnem modelu Dowson-Higginson, saj se analitični podatki o deležu zdrsavanja kletke signifikantno razlikujejo od eksperimentalnih podatkov za izjemno majhne radialne obremenitve. Eden glavnih razlogov tiči v netočnem vrednotenju sil upora olja z empiričnimi enačbami. V predstavljeni študiji je bil ta problem razrešen z novo metodo kvazidinamične analize zdrsavanja v povezavi z EHL. V raziskavi je bila za izračun debeline oljnega filma in porazdelitve tlaka v kotalnem ležaju pri izjemno majhnih obremenitvah uporabljena metoda EHL. Sile upora olja so bile uporabljene kot zamenjava za netočne rešitve empiričnih enačb po Dowson-Higginsonu. Nato sta bila ugotovljena hitrost kletke in delež zdrsavanja kletke s kombiniranjem orbit opletanja, sil upora, obremenitev, kinematičnih enačb in drugih povezanih enačb. Za reševanje je bila uporabljena Newton-Raphsonova metoda. Končno je bil preučen mehanizem zdrsavanja z ozirom na razne obratovalne parametre, kot so polmeri orbite opletanja, radialna obremenitev in viskoznost mazalnega olja. Raziskovalno področje, ki ga obravnava članek, je tribologija. 1. Delež zdrsavanja in hitrost kletke se spreminjata v času zaradi opletanja, s tem pa se poveča tveganje poškodb zaradi zdrsavanja ležaja. Stopnja zdrsavanja je neposredno odvisna od polmera opletanja in jo je mogoče zmanjšati z ustreznim zmanjšanjem polmera opletanja. 2. Delež zdrsavanja kletke se povečuje z rastjo radialnih obremenitev pri izjemno majhnih obremenitvah in velikih hitrostih, trend pri hitrosti kletke pa je ravno obraten. Vpliv opletanja na zdrsavanje ležaja se po drugi strani zmanjšuje z večanjem radialnih obremenitev. Stopnjo zdrsavanja je zato mogoče zmanjšati z ustreznim zmanjšanjem radialne obremenitve v pogojih izjemno majhnih radialnih obremenitev. 3. Delež zdrsavanja kletke se zmanjšuje s povečevanjem viskoznosti pri izjemno majhnih obremenitvah in velikih hitrostih, trend hitrosti kletke pa je ravno obraten. Opletanje ima po drugi strani manjši vpliv na zdrsavanje ležaja, ko se povečuje viskoznost. Stopnjo zdrsavanja je torej mogoče zmanjšati z ustreznim povečanjem viskoznosti v pogojih izjemno majhnih radialnih obremenitev. Pri analizi v tej študiji niso bile upoštevane izgube, npr. zaradi vrtinčenja. Izgube lahko vplivajo na analizo zdrsavanja kotalnih ležajev, zato jim bo treba posvetiti pozornost pri prihodnjih raziskavah. 1. Predstavljena analiza zdrsavanja kotalnega ležaja upošteva izjemno majhne obremenitve in opletanje. 2. Predlagana je izboljšana metoda za kvazidinamično analizo zdrsavanja v povezavi z EHL. 3. Predlagana metoda je primerna in uporabna za napovedovanje zdrsavanja pri ležajih HSLLRB. 4. Analiziran je vpliv orbit opletanja, radialnih obremenitev in viskoznosti na zdrsavanje HSLLRB. 5. Zbrani so tudi predlogi ukrepov za preprečevanje zdrsavanja. Raziskava ima teoretično in praktično vrednost v podpori analizam mehanizmov odpovedi in vibracij pri ležajih HSLLRB. Ključne besede: zdrsavanje, opletanje, kotalni ležaj, mazani elementi za viskozno blaženje, EHL, obratovalni parametri SI 14 *Naslov avtorja za dopisovanje: Tehniška univerza Xi'an; Univerza Xi'an Jiaotong, 4 Jinhua Road,, Xi'an, Kitajska, junningli@outlook.com Upoštevanje merilnih pogreškov elektronskih korektorjev pri merjenju prostornin zemeljskega plina Franc Cimerman1* - Matej Jarm1 - Branko Širok2 - Bogdan Blagojevič1 1 Plinovodi d.o.o., Slovenija 2 Univerza v Ljubljani, Fakulteta za strojništvo, Slovenia Za merjenje standardnih prostornin zemeljskega plina uporabljamo naprave za pretvarjanje prostornin (VCD) ali korektorje. Standard EN 12405-1 predpisuje zahteve in preizkušanja za konstruiranje, delovanje, varnost in skladnost korektorjev v povezavi s plinomeri. Korektorji se uporabljajo za pretvarjanje prostornin zemeljskega plina, ki jih izmerijo plinomeri pri merjenem tlaku in merjeni temperaturi ter upoštevajo sestavo zemeljskega plina, na referenčne pogoje. Korektorji sestojijo iz računske enote in temperaturnega zaznavala, ali pa iz računske enote, temperaturnega in tlačnega zaznavala. Na podlagi standarda EN 12405-1 ločimo korektorje tipa 1 in tipa 2. Korektorji tipa 1 so kompaktni ali nerazstavljivi merilniki, medtem ko so korektorji tipa 2 razstavljivi. Merilne pogreške za korektorje tipa 1 določimo v treh različnih temperaturnih in v petih različnih tlačnih stanjih, medtem ko merilne pogreške korektorjev tipa 2 preverjamo za vsak sestavni element posebej. Kljub temu bo z novim standardom za računalnike pretoka (tudi korektorje tipa 2), ki bo izšel v letu 2016, dana možnost, da lahko merilne pogreške preverjamo enako, kot jih preverjamo za korektorje tipa 1. Z namenom, da bi čim bolj natančno določili zaznano količino zemeljskega plina, je v članku prikazana možnost upoštevanja merilnih pogreškov, ki jih določimo po veljavnem nacionalnem postopku overjanja korektorjev. Na podlagi evropske direktive MID (Measuring Instrument Directive 2004/22/EC) je največji dopustni merilni pogrešek za korektorje tipa 1 in tipa 2 ±0,5 %. Veljavni nacionalni postopki overjanja korektorjev so enaki postopkom prvih overitev, ki se bodo po zahtevah MID direktive od aprila 2016 izvajali pri izdelovalcih korektorjev. Cilj prispevka je pokazati, da lahko za različne tipe elektronskih korektorjev prostornin zemeljskega plina njihove merilne pogreške, ki so bili določeni pri kontrolah korektorjev, aproksimiramo z regresijsko enačbo po metodi najmanjših kvadratov in ugotovimo značilnost modela. Zato je v prispevku na podlagi načrtovanja eksperimentov in multiregresijske analize podana in ocenjena povezanost med merilnimi pogreški in neposredno merjenimi veličinami tlaki in temperaturami pri kontroli korektorjev. Pri eksperimentalni analizi smo upoštevali osem različnih tipov korektorjev. V analizi so bili obravnavani korektorji, ki že imajo odobritev tipa po MID, in starejši korektorji, ki imajo nacionalno ali evropsko odobritev tipa. Kljub temu je bilo naše vodilo, da morajo vsi korektorji izpolnjevati kriterij dopustnega merilnega pogreška ±0,5 %. Na podlagi analize rezultatov so v članku predlagane mere za določitev primernosti regresijskega modela za posamezen korektor. Najbolj pomembna statistična cenilka za ocenjevanje primernosti uporabe multiregresijske analize pri kontroli elektronskih korektorjev zemeljskega plina je standardni pogrešek aproksimacije SEE. Pomembni statistični cenilki sta tudi regresijski koeficient r2 in signifikantnost F. V kolikor opravljamo multiregresijsko analizo kontrole korektorjev na podlagi zmanjšanega števila merilnih točk in upoštevamo zahteve Pravilnika o merilnih instrumentih za korektorje, je lahko zmanjšano število merilnih točk pri kontroli zelo primerno pri vmesnih kontrolah korektorjev na terenu. Na podlagi predlaganega načina upoštevanja merilnih pogreškov v članku lahko za poljuben elektronski korektor tipa 1 in tipa 2 izpolnimo zahteve MID glede merilnih pogreškov, ocenjena merilna negotovost pa poleg zahteve 1/3 dopustnih merilnih pogreškov upošteva še merilno negotovost zaradi aproksimacije predlagane metode. Zato lahko opisana metoda upoštevanja merilnih pogreškov predstavlja možnost za nadgraditev nacionalne procedure overjanja korektorjev, ki je primerna za oba tipa korektorjev. Ključne besede: naprave za pretvarjanje prostornin ali korektorji, merilni postopki, multiregresija, merilna negotovost Preiskava vplivov nenewtonske tekočine med prehodnimi tokovi v cevovodu Ali Majd1 - Ahmad Ahmadi1* - Alireza Keramat2 1 Tehniška univerza Shahrood, Oddelek za gradbeništvo, Iran 2 Tehniška univerza Jundi Shapur, Oddelek za strojništvo, Iran Nenadna sprememba pretoka v cevovodu povzroči velika tlačna nihanja, ki so znana pod imenom vodni (hidravlični) udar. V predstavljeni študiji je bil preučen laminami prehodni tok nenewtonske tekočine, ki se pojavi zaradi hipnega zapiranja ventila. Za simulacijo nenewtonskih vplivov so bili uporabljeni Crossovi modeli in potenčni zakon. Delo podaja razširitev klasičnega modela vodnega udara, ki se uporablja za preučevanje prehodnih pojavov v newtonski tekočini v ravni elastični cevi, podprti pri ventilu in z zadostnim številom podpor vzdolž cevovoda za preprečitev interakcij med tekočino in konstrukcijo. V ta namen so bile izpeljane kvazidvodimenzionalne enačbe vodnega udara za nenewtonske tekočine in razrešene s primernim numeričnim postopkom na podlagi metode končnih razlik. Za integracijo sistema enačb v času je bila uporabljena shema Runge-Kutta četrtega reda. Prostorski odvodi so bili diskretizirani s centralno diferenčno shemo. Za odpravo numeričnih nihanj so bili dodani disipativni členi drugega reda. Ti členi učinkujejo samo v območju velikih gradientov, v gladkih območjih pa so praktično izključeni. Za validacijo matematičnega modela ter pripadajoče numerične rešitve in njene implementacije so bili rezultati izračunov primerjani z razpoložljivimi rezultati eksperimentov. Med rezultati izračunov so bile hitrosti, strižne napetosti in porazdelitev viskoznosti po pretočnem prerezu na sredini cevovoda. Rezultati so razkrili, da vplivi nenewtonske tekočine signifikantno prispevajo k lastnostim toka po prerezu. Za podčrtanje pomena vedenja nenewtonske tekočine, ki se odraža v spremembah viskoznosti med prehodnim tokom, je predstavljenih in podrobno opisanih več numeričnih primerov, podprtih s slikami. Najprej je preučena psevdoplastična tekočina, ki se obnaša kot tekočina s strižno odvisnim upadom viskoznosti in je najbolj pogosta nenewtonska tekočina v praktični uporabi. Vhodni podatek o tekočini za analizo prehodnega toka so bile lastnosti tekočine, opredeljene s potenčnim zakonom in s Crossovimi modeli. Cilj je bil preučiti prehodne tlake zaradi hipnega zapiranja ventila na dolvodni strani. Pri spremembah tlaka, hitrostnem profilu in strižnih napetostih na steni so se pokazali vzorci razlik med newtonskimi in nenewtonskimi tokovi, ki izhajajo predvsem iz nelinearne odvisnosti med viskoznostjo tekočine in hitrostnim gradientom. Povečanje lastnosti strižno odvisnega upada viskoznosti tekočine, ki ustreza spremembam viskoznosti in strižnih napetosti v obroču prereza cevi, je bilo ugotovljeno bližje stenam cevi. Takšno vedenje tekočin s strižno odvisnim upadom viskoznosti kaže, da je območje povečane disipacije energije omejeno, zato sta manjša tudi izguba energije in tlačni padec med prehodnim pojavom. Iz časovne vrste tlaka v tekočini je bilo ugotovljeno, da je padec tlaka v tekočini z večjim strižno odvisnim upadom viskoznosti v določenem časovnem obdobju manjši. S povečanjem strižno odvisnega upada viskoznosti nenewtonske psevdoplastične tekočine in s tem relativnega padca navidezne viskoznosti so se zmanjšale tlačne izgube v cevi. Zaradi zmanjšanja navidezne viskoznosti ob steni je manjši tudi učinek povečanja linijske gostote v primerjavi z newtonskimi modeli. Sledi sklep, da strižno odvisni upad viskoznosti nenewtonskih tekočin povzroči premik območja velikih gradientov hitrosti proti stenam cevi in da so največje relativne hitrosti bližje radialnim mejam, zaradi česar prihaja do močnih fluktuacij hitrostnega profila po prerezu. Ključne besede: prehodni tok v cevi, posplošena newtonska tekočina, tekočine s strižno odvisnim upadom viskoznosti, profil hitrosti, časovna zgodovina tlaka, porazdelitev viskoznosti Analiza vpliva nepravilnosti na vozišču na življenjsko dobo pogonskih polgredi vozila Sreten Simović1* - Vladimir Popović2 - Milanko Damjanović1 1 Univerza Črne gore, Fakulteta za strojništvo, Podgorica, Črna gora 2 Univerza v Beogradu, Fakulteta za strojništvo, Beograd, Srbija Vozilo kot kompleksen tehnični sistem, ki je izpostavljen širokemu spektru vibracij, zahteva celovito analizo dinamičnega vedenja vozila, vseh ločenih sistemov in, posledično, analizo življenjske dobe elementov vozila. Ker te vibracije povzročajo predvsem nepravilnosti na vozišču in notranji viri, imajo največji vpliv nanje sile in momenti, ki nastanejo na stiku med površino in kolesi, ter vedenje deformabilnega kolesa. Analiza dinamičnega vedenja sistemov vozila mora zato običajno vključevati tudi vpliv kolesa kot deformabilnega kolesa, določitev relacij med odmiki nevzmetenih in vzmetenih mas, ter analizo nihajnega sistema, sestavljenega iz kolesa, pogonske gredi, vzmetenja in karoserije vozila. Zaradi opaznega pomanjkanja analiz in ugotovitev glede vplivov lastnosti sistemov vozila in vozišča na življenjsko dobo elementov pogonskega sklopa vozila, se lahko preizkusi utrujanja in zanesljivosti izvedejo z analizo utrujenostnih poškodb posameznih elementov pogonskega sklopa vozila pri uporabi vozila na različnih vrstah cest. Analiza vpliva nepravilnosti na vozišču na življenjsko dobo polgredi vključuje določitev matematičnega modela ter računalniško simulacijo obremenitev sistemov vozila in polgredi. Simulacija dinamičnega vedenja vozila je bila opravljena s pomočjo lineariziranih enačb prostora stanj za četrtinski model vozila, ki je vključeval pogonski sklop, sistem vzmetenja in nepravilnosti na cesti. Matematični model je bil validiran z uporabo eksperimentalnih podatkov, pridobljenih z individualno analizo potniških vozil s sprednjim pogonom in MacPhersonovo vzmetno nogo, pri čemer je bil eksperiment opravljen na realnem vozilu, ki se je premikalo v ravni črti in s konstantno hitrostjo po idealni ravni betonski površini in čez nepravilnost trikotne oblike. Rezultati so pokazali, da so sile na kolo in spremembe položaja vzmetenja v primeru, ko se vozilo premika po neravni površini vozišča, zelo kompleksne. Ko vozilo prevozi nepravilnosti na vozišču, nastanejo sile, ki so izjemno dinamične narave, smer delovanja pa je drugačna kot pri vožnji vozila po ravnem vozišču. Nepravilnosti na površini vozišča pomembno vplivajo na interakcije med kolesom in cesto, sile v točki stika med kolesom in voziščem pa so odvisne od lastnosti nepravilnosti na vozišču. Iz tega izhaja sklep, da je treba pri izbiranju parametrov sistemov vozila upoštevati njihov vpliv ne le na udobje pri vožnji, upravljanje in stabilnost, temveč tudi na obremenitve in življenjsko dobo elementov pogonskega sklopa. Lastnosti vozišča signifikantno vplivajo na napetosti in življenjsko dobo polgredi vozila tako, da se življenjska doba s povečevanjem nepravilnosti na vozišču skrajšuje. Analiza življenjske dobe je pokazala, da je treba zasnovo sistemov vozila prilagoditi pogojem uporabe vozila. Pomen takšnega pristopa se odraža v možnosti uporabe sodobnih računalniških metod analize in v dejstvu, da je analizo mogoče opraviti z razmeroma preprostimi matematičnimi modeli pogonskega sklopa, sistemov vzmetenja in lastnosti vozišča. Model omogoča analizo s preprosto prilagoditvijo konfiguracij vozila različnim tipom vozil in konfiguracijam njihovih sistemov. Ključne besede: vozilo, vozišče, nepravilnost, vzmetenje, polgred, obremenitev, življenjska doba *Naslov avtorja za dopisovanje: Univerza Črne gore, Fakulteta za strojništvo, Džordža Vašingtona bb, Podgorica, Črna gora, sretens@ac.me SI 17 Zasnova in analiza sistema dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige Yabing Cheng* - Shuaibing Yin - Xiaopeng Wang - Lichi An - Huan Liu Univerza Jilin, Fakulteta za tehniške vede in strojništvo, Kitajska Z namenom izboljšanja zmogljivosti sistema tihe krmilne verige sta bila zasnovana nova vrsta dvostranskega ubiranja in sistem dvofaznega pogona tihe krmilne verige. Predlagana je metoda za snovanje takšnih sistemov. V raziskavi dinamičnih lastnosti sistema dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige je bila preverjena racionalnost in znanstvenost metode snovanja omenjenega sistema dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige. Predstavljena je metoda snovanja sistema dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige. Rezultat snovanja je geometrijski model oblike členov tihe verige in verižnikov, izračunan pa je tudi položaj skupin napenjalnih in prostih verižnikov ter dolžina tihe verige. Na ta način je vzpostavljen sistem dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige. S primerjalno analizo sistema dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige ter enofaznega pogona verige so bile preučene dinamične lastnosti sistema (vključno s stopnjo fluktuacij napenjalnega verižnika, vibracijami členov, vrtilno hitrostjo verižnika na izpušni odmični gredi, napakami prenosa in kontaktno silo). Za podroben opis delovanja omenjenega sistema sta podana shematska diagrama dvofaznega in enofaznega verižnika ter sestavni risbi skupin napenjalnih in prostih verižnikov. Članek opisuje zasnovo sistema vključno s členi verige, verižniki in razporeditvijo sestavnih delov. Na podlagi simulacije sistema so bile preučene dinamične lastnosti za preverjanje znanstvenosti in učinkovitosti sistema. Momenti, s katerimi sta obremenjena verižnika sesalne in izpušne odmične gredi, so bili izpeljani ob upoštevanju načel delovanja sistema motornih ventilov. Stabilnost vrtenja verižnikov je bila določena s pomočjo formule za indeks stabilnosti vrtilne hitrosti. Na podlagi konstrukcije motorja je bila zasnovana nova vrsta dvostranskega ubiranja in sistem dvofaznega pogona tihe krmilne verige, predlagana pa je tudi metoda za sistematično snovanje takšnih sistemov. Preučene so dinamične lastnosti sistema in rezultati analize kažejo, da predlagani sistem izpolnjuje zahteve za delovanje tihe krmilne verige na motorju. Potrjena je tudi znanstvenost in veljavnost predlagane metode snovanja. Rezultati primerjave z enofaznim pogonom tihe krmilne verige kažejo na večjo zmogljivost sistema dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige. Članek predstavlja metodo snovanja sistema dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige ter preučuje dinamične lastnosti sistema. Opraviti bi bilo mogoče tudi analizo togo-fleksibilno sklopljenega modela. V prihodnjih raziskavah bi bilo treba sistem dvostranskega ubiranja in dvofaznega pogona tihe krmilne verige preučiti tudi eksperimentalno. Zasnovana je nova vrsta dvostranskega ubiranja in sistem dvofaznega pogona tihe krmilne verige, predlagana pa je tudi metoda za sistematično snovanje takšnih sistemov. Z uporabo takšne krmilne verige na avtomobilskem motorju se zmanjšajo napake v prenosu in vibracije sistema, delovanje krmilnega sistema pa postane zanesljivejše. Ključne besede: Metoda snovanja, dvostransko ubiranje, dvofazni pogon, sistem tihe krmilne verige, napenjalni verižnik, dinamične lastnosti SI 18 *Naslov avtorja za dopisovanje: Univerza Jilin, Fakulteta za tehniške vede in strojništvo, Renmin street 5988 , Changchun, Kitajska, chengyb@jlu.edu.cn Postopek žične elektroerozijske obdelave nekrožnih zobnikov s programsko opremo CAD/CAM César Garcia-Hernandez1* - Rafael-Maria Gella-Marin1 - José-Luis Huertas-Talón1 -Nikolaos Efkolidis1 - Panagiotis Kyratsis2 1 Univerza v Zaragozi, Oddelek za konstruiranje in proizvodno strojništvo, Španija 2 Ustanova za tehnično izobraževanje Zahodne Makedonije, Oddelek za strojništvo in industrijski dizajn, Grčija Pri konvencionalni izdelavi kovinskih zobnikov se uporabljajo postopki za oblikovanje profila zob ali preoblikovalna orodja. V zadnjih letih je bilo veliko pozornosti posvečene iskanju univerzalnih postopkov za izdelavo zobnikov, ki ne zahtevajo specialnih strojev. Nekrožni zobniki se uporabljajo v različnih tehnoloških aplikacijah za izboljšanje zmogljivosti raznih mehanskih instrumentov ali poenotenje hitrosti na montažnih linijah, kakor tudi v raziskovalne namene. Postopek izdelave tovrstnih zobnikov prinaša tudi poseben izziv prilagajanja konvencionalnih metod snovanja za preoblikovanje ali rezkanje. Te metode je mogoče uporabiti tudi pri nenamenski opremi za žično elektroerozijsko obdelavo (WEDM) ali rezkanje ter za ustvarjanje parametričnih 3D-modelov. Predstavljena je metoda za izdelavo eliptičnih in ovalnih zobnikov s postopki žične elektroerozijske obdelave. Gre za neprekinjen postopek, katerega zmogljivost ni nič manjša kot pri konvencionalnih metodah. Zobniki nekrožne oblike se običajno izdelujejo z odvzemanjem materiala, torej z rezkanjem vsakega zoba posebej, ali z oblikovanjem profila. To zahteva nadzor nad geometrijskimi in kinematičnimi spremenljivkami procesa. Razviti so bili matematični modeli za izdelavo eliptičnih in ovalnih zobnikov, opravljene so bile simulacije in metoda je bila implementirana na stroju WEDM za izdelavo dveh parov eliptičnih in ovalnih zobnikov. Metoda bi bila lahko uporabna tudi pri izdelavi orodij za brizganje plastike ali kovinskih zobnikov po meri. Najprej je bil ustvarjen 3D-model in nato je bila s pomočjo programske opreme CAM pripravljena pot žice. S tem sta bila izpolnjena dva cilja - t. j. modeliranje in sama izdelava. Metoda, predstavljena v članku, dokazuje, da lahko univerzalni numerično krmiljeni stroji uspešno nadomestijo specialne stroje za določen proizvodni proces, kot so stroji za rezkanje zobnikov. Algoritem je bil validiran z meritvijo erodiranega zoba na koordinatnem merilnem stroju. Toleranca med minimalnim in maksimalnim profilom evolvente pri izdelanih zobnikih je bila boljša od 5 ^m in izdelani zobniki tako zagotavljajo optimalno kakovost. Članek podaja primerjavo kakovosti površine in drugih parametrov pri WEDM in pri drugih postopkih obdelave z odvzemanjem materiala, kjer ne prihaja do faznih premen materiala. Sledi končni sklep, da je z razpoložljivimi stroji in optimalno izbiro pogojev procesa, kot so intenziteta, čas impulzov itd. (ob upoštevanju nasvetov proizvajalca), mogoče doseči merilnotehnične, fizikalne in kemijske lastnosti, ki so podobne kot pri ostalih proizvodnih procesih. dstranjevanje odrezkov pri uporabljeni tehnologiji je povezano s fazno premeno, ki lahko pri določenih materialih povzroči težave. Čeprav eliptični/ovalni zobniki še niso močno razširjeni v praktični uporabi, so zelo prikladni za prenašanje nekonstantnih momentov, kakor tudi za raziskave. S poenostavitvijo postopka izdelave zobnikov, ki ne zahteva specialnih strojev, postanejo procesi raziskav, razvoja in izdelave dostopnejši večjemu številu ljudi. Postopek WEDM je primeren za izdelavo nekrožnih zobnikov in predstavljeni so potrebni matematični modeli. Za izvedbo simulacij je bil implementiran algoritem po matematičnih modelih. Sledila je praktična obdelava dveh parov eliptičnih in ovalnih zobnikov na stroju za WEDM kot alternativi za konvencionalno rezkanje. Ključne besede: obdelava, žična elektroerozijska obdelava, zobniki nekrožne oblike, delovni list, CAD, CAM *Corr. Author's Address: Univerza v Zaragozi, Oddelek za konstruiranje in proizvodno strojništvo, Campus Rio Ebro, C/ Maria de Luna, 3 - 50018 - Zaragoza, Španija, garcia-hernandez.cesar@unizar.es DOKTORSKA DISERTACIJA Na Fakulteti za strojništvo Univerze v Ljubljani je obranil svojo doktorsko disertacijo: • dne 7. januarja 2016 Mitja FRANKO z naslovom: »Napovedovanje zanesljivosti naključno obremenjenih strojnih elementov« (mentor: prof. dr. Marko Nagode, somentor: prof. dr. Matija Fajdiga); Napoved zanesljivosti strojnega elementa, ki je obremenjen z naključnimi obremenitvami, predstavlja kompleksen problem. V doktorski raziskavi je predstavljen postopek, ki omogoča napoved zanesljivosti naključno obremenjenega strojnega elementa po vozliščih modela končnih elementov. Postopek temelji na statičnem modelu zanesljivosti, pri katerem je treba poznati gostoti porazdelitve verjetnosti obremenitvenega kolektiva ter zdržljivosti na vsakem napetostnem nivoju in gostoto porazdelitve verjetnosti ekvivalentnih amplitudnih napetosti. Gostota porazdelitve verjetnosti ekvivalentnih amplitudnih napetosti je pridobljena z razvito statistično transformacijo, s katero se informacija o porazdelitveni funkciji amplitud in srednjih vrednosti napetosti rainflow matrike ohranja. Information for Authors All manuscripts must be in English. Pages should be numbered sequentially. The manuscript should be composed in accordance with the Article Template given above. The maximum length of contributions is 10 pages. Longer contributions will only be accepted if authors provide juastification in a cover letter. For full instructions see the Information for Authors section on the journal's website: http://en.sv-jme.eu . SUBMISSION: Submission to SV-JME is made with the implicit understanding that neither the manuscript nor the essence of its content has been published previously either in whole or in part and that it is not being considered for publication elsewhere. All the listed authors should have agreed on the content and the corresponding (submitting) author is responsible for having ensured that this agreement has been reached. The acceptance of an article is based entirely on its scientific merit, as judged by peer review. Scientific articles comprising simulations only will not be accepted for publication; simulations must be accompanied by experimental results carried out to confirm or deny the accuracy of the simulation. Every manuscript submitted to the SV-JME undergoes a peer-review process. The authors are kindly invited to submit the paper through our web site: http://ojs.sv-jme.eu. The Author is able to track the submission through the editorial process - as well as participate in the copyediting and proofreading of submissions accepted for publication - by logging in, and using the username and password provided. SUBMISSION CONTENT: The typical submission material consists of: - A manuscript (A PDF file, with title, all authors with affiliations, abstract, keywords, highlights, inserted figures and tables and references), - Supplementary files: • a manuscript in a WORD file format • a cover letter (please see instructions for composing the cover letter) • a ZIP file containing figures in high resolution in one of the graphical formats (please see instructions for preparing the figure files) • possible appendicies (optional), cover materials, video materials, etc. Incomplete or improperly prepared submissions will be rejected with explanatory comments provided. In this case we will kindly ask the authors to carefully read the Information for Authors and to resubmit their manuscripts taking into consideration our comments. COVER LETTER INSTRUCTIONS: Please add a cover letter stating the following information about the submitted paper: 1. Paper title, list of authors and their affiliations. 2. Type of paper: original scientific paper (1.01), review scientific paper (1.02) or short scientific paper (1.03). 3. A declaration that neither the manuscript nor the essence of its content has been published in whole or in part previously and that it is not being considered for publication elsewhere. 4. State the value of the paper or its practical, theoretical and scientific implications. What is new in the paper with respect to the state-of-the-art in the published papers? Do not repeat the content of your abstract for this purpose. 5. We kindly ask you to suggest at least two reviewers for your paper and give us their names, their full affiliation and contact information, and their scientific research interest. The suggested reviewers should have at least two relevant references (with an impact factor) to the scientific field concerned; they should not be from the same country as the authors and should have no close connection with the authors. FORMAT OF THE MANUSCRIPT: The manuscript should be composed in accordance with the Article Template. The manuscript should be written in the following format: - A Title that adequately describes the content of the manuscript. - A list of Authors and their affiliations. - An Abstract that should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, as well as the methodology employed. It should summarize the results and state the principal conclusions. - 4 to 6 significant key words should follow the abstract to aid indexing. - 4 to 6 highlights; a short collection of bullet points that convey the core findings and provide readers with a quick textual overview of the article. These four to six bullet points should describe the essence of the research (e.g. results or conclusions) and highlight what is distinctive about it. - An Introduction that should provide a review of recent literature and sufficient background information to allow the results of the article to be understood and evaluated. - A Methods section detailing the theoretical or experimental methods used. - An Experimental section that should provide details of the experimental set-up and the methods used to obtain the results. - A Results section that should clearly and concisely present the data, using figures and tables where appropriate. - A Discussion section that should describe the relationships and generalizations shown by the results and discuss the significance of the results, making comparisons with previously published work. (It may be appropriate to combine the Results and Discussion sections into a single section to improve clarity.) - A Conclusions section that should present one or more conclusions drawn from the results and subsequent discussion and should not duplicate the Abstract. - Acknowledgement (optional) of collaboration or preparation assistance may be included. Please note the source of funding for the research. - Nomenclature (optional). Papers with many symbols should have a nomenclature that defines all symbols with units, inserted above the references. If one is used, it must contain all the symbols used in the manuscript and the definitions should not be repeated in the text. In all cases, identify the symbols used if they are not widely recognized in the profession. Define acronyms in the text, not in the nomenclature. - References must be cited consecutively in the text using square brackets [1] and collected together in a reference list at the end of the manuscript. - Appendix(-icies) if any. SPECIAL NOTES Units: The SI system of units for nomenclature, symbols and abbreviations should be followed closely. Symbols for physical quantities in the text should be written in italics (e.g. v, T, n, etc.). Symbols for units that consist of letters should be in plain text (e.g. ms-1, K, min, mm, etc.). Please also see: http://physics.nist.gov/cuu/pdf/sp811.pdf . Abbreviations should be spelt out in full on first appearance followed by the abbreviation in parentheses, e.g. variable time geometry (VTG). The meaning of symbols and units belonging to symbols should be explained in each case or cited in a nomenclature section at the end of the manuscript before the References. Figures (figures, graphs, illustrations digital images, photographs) must be cited in consecutive numerical order in the text and referred to in both the text and the captions as Fig. 1, Fig. 2, etc. Figures should be prepared without borders and on white grounding and should be sent separately in their original formats. If a figure is composed of several parts, please mark each part with a), b), c), etc. and provide an explanation for each part in Figure caption. The caption should be self-explanatory. Letters and numbers should be readable (Arial or Times New Roman, min 6 pt with equal sizes and fonts in all figures). Graphics (submitted as supplementary files) may be exported in resolution good enough for printing (min. 300 dpi) in any common format, e.g. TIFF, BMP or JPG, PDF and should be named Fig1.jpg, Fig2.tif, etc. However, graphs and line drawings should be prepared as vector images, e.g. CDR, AI. Multi-curve graphs should have individual curves marked with a symbol or otherwise provide distinguishing differences using, for example, different thicknesses or dashing. Tables should carry separate titles and must be numbered in consecutive numerical order in the text and referred to in both the text and the captions as Table 1, Table 2, etc. In addition to the physical quantities, such as t (in italics), the units [s] (normal text) should be added in square brackets. Tables should not duplicate data found elsewhere in the manuscript. Tables should be prepared using a table editor and not inserted as a graphic. REFERENCES: A reference list must be included using the following information as a guide. Only cited text references are to be included. Each reference is to be referred to in the text by a number enclosed in a square bracket (i.e. [3] or [2] to [4] for more references; do not combine more than 3 references, explain each). No reference to the author is necessary. References must be numbered and ordered according to where they are first mentioned in the paper, not alphabetically. All references must be complete and accurate. Please add DOI code when available. Examples follow. Journal Papers: Surname 1, Initials, Surname 2, Initials (year). Title. Journal, volume, number, pages, DOI code. [1] Hackenschmidt, R., Alber-Laukant, B., Rieg, F. (2010). Simulating nonlinear materials under centrifugal forces by using intelligent cross-linked simulations. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, D0I:10.5545/sv-jme.2011.013. Journal titles should not be abbreviated. Note that journal title is set in italics. Books: Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication. [2] Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken. Note that the title of the book is italicized. Chapters in Books: Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages. [3] Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordić, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553576. Proceedings Papers: Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages. [4] Štefanić, N., Martinčević-Mikić, S., Tošanović, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427. Standards: Standard-Code (year). Title. Organisation. Place. [5] ISO/DIS 16000-6.2:2002. Indoor Air — Part 6: Determination of Volatile Organic Compounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva. WWW pages: Surname, Initials or Company name. Title, from http://address, date of access. [6] Rockwell Automation. Arena, from http://www.arenasimulation.com, accessed on 200909-07. EXTENDED ABSTRACT: When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters). The instruction for composing the extended abstract are published on-line: http://www.sv-jme.eu/information-for-authors/ . COPYRIGHT: Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to transfer copyright to SV-JME when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Transfer Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on http://en.sv-jme.eu/. PUBLICATION FEE: Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 240.00 EUR (for articles with maximum of 6 pages), 300.00 EUR (for articles with maximum of 10 pages), plus 30.00 EUR for each additional page. The additional cost for a color page is 90.00 EUR. These fees do not include tax. Strojniški vestnik - Journal of Mechanical Engineering Aškerčeva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu http://www.sv-jme.eu Contents Papers Matej Steinacher, Borut Žužek, Darja Jenko, Primož Mrvar, Franc Zupanič: Manufacturing and Properties of a Magnesium Interpenetrating Phase Composite Junning Li, Wei Chen, Libo Zhang, Taofeng Wang: An Improved Quasi-Dynamic Analytical Method to Predict Skidding in Roller Bearings under Conditions of Extremely Light Loads and Whirling Franc Cimerman, Matej Jarm, Branko Širok, Bogdan Blagojevič: Taking in Account Measuring Errors of Volume Conversion Devices in Measuring of the Volume of Natural Gas Ali Majd, Ahmad Ahmadi, Alireza Keramat: Investigation of Non-Newtonian Fluid Effects during Transient Flows in a Pipeline Sreten Simović, Vladimir Popović, Milanko Damjanović: Analysis of the Influence of Pavement Irregularities on the Lifespan of a Vehicle's Drive-Wheel Half Shaft Yabing Cheng, Shuaibing Yin, Xiaopeng Wang, Lichi An, Huan Liu: Design and Analysis of Double-side Meshing and Dual-phase Driving Timing Silent Chain System César Garda-Hernandez, Rafael-Maria Gella-Marin, José-Luis Huertas-Talón, Nikolaos Efkolidis, Panagiotis Kyratsis: WEDM Manufacturing Method for Noncircular Gears, Using CAD/CAM Software 79 9770039248001