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: Abografika, printed in 300 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 Mitjan Kalin University of Ljubljana, Faculty of Mechanical Engineering, Slovenia International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, University of Engineering and Technology, Pakistan Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, UL, Faculty of Mechanical Engineering, Slovenia Franci Čuš, UM, Faculty of Mechanical Engineering, Slovenia Janez Diaci, UL, 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 Lubben, University of Bremen, Germany George K. Nikas, KADMOS Engineering, UK José L. Ocana, Technical University of Madrid, Spain 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 Branko Vasič, University of Belgrade, Faculty of Mechanical Eng., Serbia Arkady Voloshin, Lehigh University, Bethlehem, USA Vice-President of Publishing Council Jože Balič University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: Prototype of the Pure Hydro-Elastic Actuator (PHEA): a hydraulic force control system using high volumetric expansion hoses. Image courtesy: LASHIP/UFSC - Laboratory of Hydraulic and Pneumatic Systems, Brazil ISSN 0039-2480, ISSN 2536-2948 (online) 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://www.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. © 2018 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. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)10 Contents Contents Strojniški vestnik - Journal of Mechanical Engineering volume 64, (2018), number 10 Ljubljana, October 2018 ISSN 0039-2480 Published monthly Papers Job Angel Ledezma Pérez, Edson Roberto De Pieri, Victor Juliano De Negri: Force Control of Hydraulic Actuators using Additional Hydraulic Compliance 579 Simon Marčič, Sebastijan Seme, Julijan Jan Salamunic, Zdravko Praunseis, Stanislav Pehan: Heating by Flat Plate Collector Combined with a Heat Pump 590 Jun Liu, Weilong Liu, Xiaofeng Wang: Research of Crack Detection and Delay Crack Propagation on a Nonlinear Rotor 601 Riad Ramadani, Marko Kegl, Jožef Predan, Aleš Belšak, Stanislav Pehan: Influence of Cellular Lattice Body Structure on Gear Vibration Induced by Meshing 611 Fabrício Alves de Almeida, Guilherme Ferreira Gomes, Rachel Campos Sabioni, José Henrique de Freitas Gomes, Vinícius Renó de Paula, Anderson Paulo de Paiva, Sebastiäo Carlos da Costa: A Gage Study Applied in Shear Test to Identify Variation Causes from a Resistance Spot Welding Measurement System 621 Qihui Yu, Xin Tan, Jianguo Meng, Fei Shang, Xueqing Hao: Performance Characteristics Research on Power System of a Four-Valve Compressed Air Engine 632 Daniel Nasulea, Gheorghe Oancea: Integrating a New Software Tool Used for Tool Path Generation in the Numerical Simulation of Incremental Forming Processes 643 Strojniški vestnik - Journal of Mechanical Engineering 64(2018)10, 579-589 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2018.5339 Original Scientific Paper Received for review: 2018-03-07 Received revised form: 2018-05-29 Accepted for publication: 2018-07-02 Force Control of Hydraulic Actuators using Additional Hydraulic Compliance Job Angel Ledezma Pérez - Edson Roberto De Pieri - Victor Juliano De Negri* Federal University of Santa Catarina, Brazil This paper presents an approach for improving the performance of hydraulic force control systems by adding hydraulic compliance. Typically, force control systems require a flexible coupling system, such as a spring between the actuator and the load, to achieve a non-oscillatory response. To avoid the use of springs, this study proposes the use of special hoses to add compliance to the system hydraulically. With this approach, a direct connection between the actuator and the load is feasible, simplifying the mechanical assembly and saving physical space in the axial direction. An analytical model is proposed to estimate the appropriate hydraulic stiffness value and to select a commercial hose based on catalogue data. Moreover, a robust controller based on quantitative feedback theory is designed to improve the stability, performance and disturbance rejection of the force control system. The system performance is demonstrated by dynamic simulation and experimental results. Keywords: hydraulic force control, hydraulic compliance, high expansion volumetric hose, QFT-based control Highlights • A hydraulic force control system with added hydraulic compliance is presented. • Passive compliance is provided with the use of hydraulic hoses instead of mechanical springs. • An analytical procedure for selecting and sizing hydraulic hoses is proposed. • The presented achievements are validated by dynamic simulation and experimental results. 0 INTRODUCTION Several industrial applications require controlled force for handling materials or machining tasks in which the end effector is in contact with the piece or equipment. In most cases, position/trajectory control is required at the same time as contact force control. Assembly of mechanical components, polishing and grinding of complex pieces, milling, and endurance tests are some examples for which there are demands for force-controlled actuators in substitution to manual work [1] and [2]. As discussed in [3], the problem of force control is more complex than the problem of position control, which often makes its practical implementation more difficult. In contrast, to force control systems, the dynamics of the load is not included in the closed-loop in the case of positioning systems, and thus good performance and disturbance rejection are more easily achieved. One way to overcome the inherent problems associated with force control systems is by increasing the compliance. System compliance can be controlled or modified in active or passive ways [4] and [5]. Active compliance is obtained with force-based feedback control (software), allowing more flexibility to perform tasks for which compliance is required. Typically, this approach requires additional measuring for control and collision detection [6]. In contrast, passive compliance is typically obtained by introducing mechanical elastic components (hardware) between the actuator and the environment. This approach enables energy storage, and it is appropriate for safety applications [7]. The series elastic actuator (SEA) represents one well-known example of a passive compliance actuation system [8]. This kind of actuator has been applied to robotics and biomechanics due to its simplicity and good performance in force control [9] and [10]. In an SEA, a spring is placed between the actuator and the load, providing good force fidelity, low output impedance, and shock tolerance capability. However, the system bandwidth is reduced due to the use of this compliant element [8]. In addition to the use of a spring to increase the system compliance and improve the system performance, researchers have employed a spring as a simplified representation of the load. Good force responses have been obtained using control techniques such as nonlinear PD (proportional-derivative), fuzzy PID (proportional-integral-derivative), QFT (quantitative feedback theory), and backstepping [11] to [14]. References [15] and [16] included additional load motion measurements in order to cancel the disturbance characteristic and obtain better force tracking. In contrast, the fluid compressibility in hydraulic actuators introduces a spring effect that is *Corr. Author's Address: Federal University of Santa Catarina, Mechanical Engeeniring Department, LASHIP, Campus Universitario, Trindade, 88040-900, Florianopolis, SC, Brazil, victor.de.negri@ufsc.br characterized by the hydraulic stiffness. This can be shaped by using hydraulic capacitive components such as flexible hoses or accumulators. The use of accumulators connected to the cylinder chambers has been described in conference papers and patents [17] to [19]. However, these publications do not take into account the mathematical modelling to determine the hydraulic compliance needed to achieve the force control requirements and select the hydraulic capacitive component. In this context, the present paper introduces a pure hydro-elastic actuator (PHEA) that uses only hydraulic components to reduce the transmission stiffness. Mathematical modelling and analysis are carried out, and a procedure for sizing the required hydraulic components in order to increase the compliance is determined. Instead of accumulators, the use of hoses with high volumetric expansion (HVE) is investigated. Considering the effects of compliance addition and the possible parametric variations of the system, this paper considers the use of a robust controller based on the quantitative feedback theory (QFT) technique [20]. Furthermore, the load is assumed to be continuously in contact with the actuator, and the load-displacement acts as a disturbance over the force control system. The following section provides a brief description of the hydraulic force control system. The nonlinear and linear modelling is presented in Section 2 followed by the proposed method for determining the hydraulic stiffness. The force controller based on the QFT technique is then designed, and simulation and experimental results are discussed. The last section presents the conclusions. 1 THE FORCE CONTROL SYSTEM A hydraulic actuation system is very stiff due to the low fluid compressibility. Direct contact between the cylinder and the load would cause the entire system to be stiffer, and even a small control signal sent to a servovalve could generate large force variations due to the high open-loop gain. Reference [21] emphasizes the advantages of including a spring in a force control system. In [8], the procedure of the spring selection for the SEA was improved by using two types of linear actuators: an electromechanical and a hydraulic actuator. In both cases, the introduction of a compliant element between the actuator and the load decreases the output impedance. Additionally, the measurement of the spring compression is used to calculate indirectly the force applied over the load (Fig. 1a). component a) b) Fig. 1. Compliance addition in a hydraulic force system using a) spring and b) HVE hoses The use of the spring limits the rate of the force applied by the hydraulic actuator, ensuring some degree of isolation between the hydraulic system and the movement of the load (environment) and maintaining a stable and robust applied force. The use of a hydraulic compliant component instead of a spring is proposed. Fig. 1b shows an option using high volumetric expansion (HVE) hoses between the servovalve and cylinder. The use of accumulators instead of hoses could also be possible for similar purposes. 2 MATHEMATICAL MODELING Two mathematical models are explained in this section: a nonlinear model, which is used for simulation in order to obtain force responses; and a linear model, used for the controller design and hose sizing. Both models are based on Fig. 1b. 2.1 Nonlinear Modelling of the Hydraulic System 2.1.1 Servovalve Modelling In the design of the hydraulic circuit, a symmetrical servovalve is used (Fig. 2) [22]. The dynamic relationship between the input control signal (Uc) and the spool displacement, represented by an equivalent voltage (Ucsp), can be approximated by a second-order function: U = 1 dUCp dU Csp dt2 dt U Csp: (1) where ^v is the natural frequency of the valve and represents the damping ratio of the valve. The flow rate through the valve, including the effects of internal leakage can be described by the following equations [23], where a control signal range of -10 V to +10 V is assumed: for UCsv > 0: IvA = IvB = K U Csp UC U Csp K U v ucn K ■K,„. yjps "pa -kvinp4pä—pt,(2) VPb - Pt - KvinpVPs - Pb > (3) for UcSp < 0: ( IvA =■ K„ Csp I tfC -Kv. \/Pa - Pt +Kvinpy[pS~ vinp ( Pa , (4) IvB =- K U Csp U K -jps- Pb +Kvinp4pb~ Pt , (5) where qvA and qvB are the flow rates at ports A and B, respectively; pA andpB are the pressures in lines A and B, respectively; pS and pT represent the supply and the reservoir pressures, respectively; UCn is the nominal control signal. Kvp is the partial flow coefficient of the valve; and Kvinp is the internal partial leakage coefficient. iva l \1m Pa T t pB èé .....a......mm *-> i I »-> iw Ps1 1Pt Fig. 2. Flow rates and pressures at the servovalve 2.1.2 Cylinder Modelling Considering Fig. 3 and applying the continuity equation for each cylinder chamber yields: and q - A dip+VdPA , dxp VB dpB vB B dt ße dt (6) (7) where Aa and AB are the piston areas in chambers A and B, respectively; xp represents the piston displacement, and is the effective bulk modulus. Fig. 3. Variables and parameters of the cylinder The displacement of the rod can be modelled using Newton's second law: F - F - F = M- d2. dt2 (8) where FH = pAAA - pBAB is the hydraulic force, Ffr represents the friction force, Fe is the applied force, and M is the piston mass. In turn, the force applied over the load can be represented as: Fe = KS (p - ) (9) such that KS is the load cell stiffness and xL represents the load displacement. 2.2 Linear Modelling of the Hydraulic System 2.2.1 Servovalve Modelling The flow rate through the valve (Eqs. (2) to (5)), linearized for any UCspi andpLi, can be represented as: IvC = KuPCsp - KC¡PL, (10) where: KqUi = Kv UC yj(ps -pt -sgn(ucsp)pli)2, (11) and K„; = K,„ U vp I Csp; I V8 uCa V Ps - Pt - sgn(uCsp, )Pl, V8 1 Ps - Pt - sgn(uCsp,)Pl, 1 Ps - Pt +sgn(uCsp,)Pl, (12) such that KqUi is the flow-voltage gain at operating point i; Kci the flow-pressure coefficient at operating point i; pL = pA - pB the load pressure; and qvC the control flow rate considered as the average of the qvA and qvB. p 1 + + 2.1.2 Cylinder Modelling Since the total variation of volume is the same for the two chamber volumes, it is possible to substitute VA and VB in Eqs. (6) and (7) with the parametric uncertainty V. This approach is suitable for the use of the quantitative feedback theory (QFT). Therefore, the hydraulic stiffness (KH) can be generalized as: KH = Pe A ße A 2ße AU V V V A B ' (13) and the linear model of the symmetrical cylinder results: . dXp A dpL = 4, —- + , * ^ dt kh dt (14) where Au = Aa = AB is the useful area of the piston. The friction force is a nonlinear function of velocity and can be described by a variable viscous friction coefficient (f) [24]. This coefficient can be considered as an uncertain parameter, and this assumption allows the motion equation to become linear, yielding: dx d2 xp AuPl - fv dp - F = M —p. (15) dt dt 2.3 Open Loop Transfer Function For the open-loop (OL) transfer function, the values of KqUi and Kci are considered to be evaluated at null operating point (i = 0), where UCsp0 = 0, qvC0 = 0, and PL0 = 0. Combining Eqs. (1), (9), (10) and (14) with (15), the OL transfer function that relates the output force (Fe) to the input control signal (UC) and the load displacement (xL) is: Fe(s) = C0UC(s) a5s5 + a4 s4 + a3s3 + a2 s2 + a1s + a0 KS (b4s4 + b3s3 + b2s2 + bjs + b0 )sXL(s) (16) where: a0 =®nv2KcoKHKS> a1 =®nv2Au2(KH + KS) + Kc0KH(2Îv®nvKS + ®nv2/v), a2 =24®nvAu2(KH + Ks) + fflnv2Au2/v + KcoKH(2iv®nv/v + ®nv2M + Ks), a3 =2^v®nvAu2/v + ®nv2Au2M + KÄ^^nM + /) + Au2(Kh + Ks), a4 =2tv®n-Au2M + Au2/v+ KcoKhM, a5 =Au2M, bo =®nv2Au2KH + ®nv2/vKcoKH, b =24®nv4u2KH + fflnv2Au2/v + K^h^®,/ + ®nv2M), b2 =2^v®nvAu2/v + ®nv2Au2M + Au2KH + Kc0KH(2{v®nvM+/v), b3 =2Zv®nvAu2M + A/ + KcKM, b4 =Au2M, c0 =AuKHKqU0KS® 3 SELECTION OF THE HYDRAULIC HOSE In this section, a procedure for sizing the hydraulic hose required to add hydraulic compliance to the system is proposed. It consists of a few analytical expressions that allow the selection of commercial HVE hoses in a simple way. To demonstrate its effectiveness, dynamic responses using the nonlinear model as well as experimental results are reported in Section 5. 3.1 Selecting the Desired Hydraulic Stiffness Based on Eq. (16), and expressing the input control signal UC(s) as G(s)(Fref (s) - Fe(s)), where G(s) and Fref(s) represent the controller transfer function and the reference force, respectively, the closed-loop (CL) transfer function results in: G (s) N( s) F( S) = D( s) + G (s) Nj( s) N2 (s) F-ef (s) -XL( s), (17) D(s) + G(s) Nl(s) where: N(s) = C0, N2(s) = KS (b4s4 + b3s3 + b2s2 + b1s + b0) s, D(s) = a5s5 + a4s4 + a3s3 + a2s2 + a1s + a0. The first part of Eq. (17) determines the system performance in response to a reference force. The second part of Eq. (17) describes how the system behaves when a disturbance input XL(s) occurs. Therefore, it is related to the output impedance and can be considered as a measure of the system's capability against external disturbances. The hydraulic stiffness KH affects both parts of Eq. (17) as reported below. Therefore, the hose selection will be based on a trade-off between performance and disturbance rejection. 3.1.1 Output Impedance The output impedance can be analysed considering the transfer function Fe(s)/XL(s) presented in the second part of Eq. (17). The presence of a zero at the origin introduces an implicit disturbance rejection. However, based on the initial value theorem, an initial overshoot equal to the equivalent stiffness (KHKS/(KH+KS)) times a disturbance amplitude is expected. Therefore, increasing the rejection capability against load movement requires that KH is as low as possible. Nevertheless, this value must be greater than or equal to the hydraulic stiffness required to achieve the required performance, which is discussed next. 3.1.2 Desired Tracking Control Ratio An approximation of the first part of Eq. (16) expressed by first and second order terms can be written as: f V V / A V J A Fe(s) Uc( s) KqU0 Keq / Au s + K co K eq/ \ s + dls + d0 j (18) where: do = (kh + ks)/m, di = f/m, fiï 2 eo ei K eq KsKh/(Ks+Kh), representing the equivalent stiffness. The first term is equivalent to a linear modelling of the open-loop system neglecting mass, friction and valve dynamics. The middle term is related to Newton's second law including the effect of the hydraulic stiffness. Finally, the last term represents the servovalve dynamics. It is noteworthy that, using the parameters presented in the Appendix, Eq. (18) is an almost exact approximation of Eq. (16), with a maximum difference between curves of 0.2 dB in magnitude and 0.003° in phase. The open loop poles of Eq. (18) are shown in Table 1. for different hydraulic stiffness values for which KHn corresponds to the nominal stiffness when using rigid pipes (instead of hoses). The real pole related to the hydraulic subsystem is -Kc0Keq/Au2, and it is dominant over the other two complex poles. This fact becomes more noticeable as the hydraulic stiffness is reduced. The damped natural frequency of the mechanical subsystem is loosely influenced by the hydraulic stiffness variation. Due to the dominance of the real pole, the first term in Eq. (16) can be substituted by the first order term in Eq. (18) and the corresponding closed-loop transfer function is: Fe(s) = G(s)(KqU0Keq / Au) Fref (s) S + ( Kc0 Keq / Au2) + G(s)(KqU0 Keq / Au ) Table 1. Pole locations for different values of KH Hydraulic Mechanical Servovalve dynamics dynamics dynamics Kh [N/m] KqUoKe,/ A d0 e0 (s + Kco Keq/ A ) (s2 + dls + d0 ) (s2 + es + e0 ) KHn = 1.5x107 -2.963 -11.6 ± 2.93 x103i -989 ± 479i 0.1 KHn -0.336 -11.4 ± 2.75 x103i -989 ± 479i 0.01 KHn -0.034 -11.4 ± 2.73x103i -989 ± 479i A desired tracking control ratio for the system can be specified by: Fe(s)/ Ffs) = Kss / (zds +1), (20) where KSS is the steady-state gain, and rd is the desired time constant. Comparing Eqs. (19) and (20) and assuming a proportional controller, a correlation between the required hydraulic stiffness and the proportional gain (Kp) can be obtained: Kh = Ks A KT (Kc0 + AuKpKql]0) - (21) 3.2 Selecting a Commercial Hydraulic Hose The hose diameter can be calculated based on the maximum flow rate (qvmax) and the recommended fluid velocity in the hose (voil) according to: Dho =yj41vmax1 nVoil . (22) Hose manufacturers recommend using fluid velocities of 2 m/s to 4 m/s for pressure lines [25] and [26]. Given that Au = Aa = AB, the maximum and minimum values for hydraulic stiffness (KH) can be expressed by: K ße Au2 (Au L + 2Vh0 ) Vh0 (Au L + Vho ) > at xp = 0 or xp = L, and K 4ße A L at xp - —, (23) (24) A L + 2Vho where Vho represents the volume of oil trapped in the hose coupled to each cylinder chamber, and L is the 0 piston stroke. For simplification, henceforth Eq. (24) will be used for the calculation of KH. The effective bulk modulus can be represented as: r PoPhoSS ße = ßoßho ßo +ßho ßo + r ßh (25) where ßo is the fluid bulk modulus, ßho represents the dynamic hose bulk modulus, ßhoSS is the static hose bulk modulus, and is the ratio ßho/ßhoss- According to [27], the hose bulk modulus changes during the system operation, i.e., dynamic bulk modulus (ßho) is higher than static bulk modulus (ßhoSS) due to a phenomena called dynamic hardening. Those authors state that the rß is about 4 to 5 for nylon braid hoses, which is the case of the HVE hoses used in this paper. According to the experimental results presented in [27], the static hose bulk modulus can be expressed by: AioSS _ V dp oho y (26) V Oho where AVho represents the change in the hose volume, V0ho is the initial volume of the hose, and p is the working pressure. This equation can be rewritten in order to calculate fihoSS as a function of hose catalogue data, yielding: (Ah0 + E)_( + 4E) AioSS _ " 4s (27) where E represents the volumetric expansion of the hose, defined as AVho/Lho, being Lho the hose length; Dho is the hose diameter, and s corresponds to dE/dp. E is a parameter determined at the working pressure (p), and s can be assumed constant for a specific hose, corresponding to the slope of straight lines in graphs of E versus p provided by manufacturers [25] and [26]. Combining Eqs. (24), (25) and (27) yields: Lo = K>2 + 4E) - 2AL. (28) nDho2Kh ( + n (02 + 4E nDh Eq. (28) shows that to obtain smaller hydraulic stiffness Kh, longer hose Lho is required. Furthermore, the greater the hose's volumetric expansion, the lower the hose length to achieve the same hydraulic stiffness. If necessary, the process of hose selection using Eqs. (21) and (28) can be iterative to avoid obtaining long hoses, which may not be suitable for implementation. For example, if a calculated hydraulic stiffness results in a long hose, the proportional gain can be decreased in Eq. (21), resulting in a higher KH and smaller Lho. Therefore, the final value of hydraulic stiffness is obtained from a trade-off between KH, Kp, and Lho. In contrast, Eq. (28) can also be used to calculate Kh based on a pre-selected hose with the required Dho (from Eq. (22)) and the desired Lho. After that, the proportional gain to achieve the required closed loop response can be calculated by isolating it in Eq. (21). The use of a proportional controller in this design stage allows obtaining a very good approximation for the required hydraulic stiffness. A specific controller design and the system dynamic analysis using a nonlinear model can then be carried out as shown in Sections 5 and 6. 3.3 Example of Hose Selection Consider a Pure Hydro-Elastic Actuator (PHEA), as shown in Fig. 1b, able to apply forces up to 9000 N. The desired tracking control ratio is specified according to Eq. (20) with a time constant rd equal to 50 ms. The system parameters are according to the Appendix, with KqU0 and Kc0, calculated at null operating based on Eqs. (11) and (12). Considering the maximum flow rate equal to the test rig pump supply (1.6X10-4 m3/s) and fluid velocity in the line equal to 2 m/s, a commercial hose diameter of 12.7 mm (1/2 in) was selected using Eq. (22). Analysing HVE hose catalogues, the EATON Synflex ® 3130-08 was selected [25]. The hose volumetric expansion is 1.56x10-5 m3/m @ 7x106 Pa (4.7 cc/ft @ 1000 psi), resulting in a static bulk modulus 0#hoss) of 6.36x107 Pa. Based on Eq. (26), an rp value equal to 5 was considered. As discussed before, the system performance must be determined by a trade-off between KH, Kp, and Lho. In this study, a Lho equal to 1.5 m was specified and using Eq. (28) the resulting KH is 2.24x106 N/m and, based on Eq. (21), Kp is 2.7 x10-4 The hydraulic force control system according to the configuration shown in Fig. 1b was assembled in a test rig as shown in Fig. 4. The system comprises two hoses interconnecting each cylinder chamber with the servovalve. A load cell is fixed at the cylinder rod, and it will be in contact with a metal block attached to the test rig frame. The parameter values of this experimental setup are presented in the Appendix. For the QFT controller design and dynamic simulation, KS equal to 2x107 N/m was used, corresponding to the d i: equivalent stiffness resulting from the load cell and environment compliances. Fig. 4. Hydraulic test bench used for experiments 4 CONTROLLER DESIGN USING LINEAR QFT The QFT allows designing a linear controller for a nonlinear system, assuming it as linear under the effect of disturbances and parametric variations [20]. In this study, the controller design was carried out using the QFT frequency domain control design toolbox for use with Matlab [28]. The objective of the QFT-based control design is to synthesize a prefilter (F(s)) and a controller (G(s)) such that the force responses of the system always fall within a predefined time-domain tolerance described by upper and lower limits. Typically, these limits are based on second-order specifications, such that the time constant value (rd = 50 ms) used at the hose selection procedure (Section 3.1.2) is converted on a settling time (?si% = 200 ms). In this study, the upper and lower limits were defined as: • Upper Limit (Bu(t)): settling time (tsi%) of 0.1 s and maximum overshoot of 1 %. • Lower limit (BlW): settling time (tsi%) of 0.3 s without overshoot. Translating the specifications from the time domain to the frequency domain yields: W = 2 2116 s" + 73.6s + 2116 A 75 s + 75 and TlS = 400 s2 + 40s + 400 A s + 50 50 (29) (30) The addition of a zero and a pole in Eqs. (29) and (30), respectively, relax the requirements for the controller design in the frequency domain without affecting the responses in the time domain. Table 2 shows the parametric uncertainties assumed in the plant model described by Eq. (16). Table 2. Parameter uncertainties at Eq. (16) Description Nominal value Range fv [Ns/m] 100 100 to 3x104 Kqu [m3/(s-V)] 5.4x10-5 3.23x10-5 to 5.61x10-5 Kci [m3/(s-Pa)] 6.42x10-13 6.42x10-13 to 7x10-11 Kh [N/m] 2.24x106 2.24x106 to 4.11x106 ^nv [rad/s] 1099 879 to 1319 & 0.9 0.72 to 1.08 To define the uncertainties, the following considerations were assumed: • fv varies from a minimum value related to the Coulomb friction to a maximum value representing the stiction, which was obtained experimentally, • the range of parametric uncertainty for Kqu was calculated using Eq. (11), assuming no load (Pli=0) and the operating point at the maximum power (PLi=2ps/3). The nominal value was obtained through experiments according to ISO 10770-1 [29], • the range for Kci was calculated using Eq. (12). The nominal value was obtained based on the results of the internal leakage measurement test according to ISO 10770-1 [29], • Kh varies from the minimum to maximum hydraulic stiffness calculated by Eqs. (23) and (24), • Finally, a ±20% variation around the nominal values extracted from the catalogue is assumed for ^nv and tv. The other parameter values used for the controller design are shown in the Appendix, with the exception of KS that was changed in order to represent the equivalent stiffness instead of just the load cell characteristic. Based on the measurement of the environment deflection, the resulting Ks was 2*107 N/m. For the QFT boundary generation, robust stability, disturbance rejection and reference tracking criteria are considered in this study. The robust stability bounds are defined through [20]: TiM P(jö)G (jö) 1 + P(jö)G(jö) 1 + L(jö) < Ôi(w>) = 1.3, ©gQ15 (31) where P(]oS) represents the plant, G(j^) represents the controller, L(]w) is known as the loop transmission function, and is a constant constraint calculated assuming a gain margin of 5 dB. The subset of analysis frequencies ^ is {0.01, 0.1, 1, 5, 10, 50, 100, 150} rad/s. For the disturbance rejection bounds, the constraint function represents the disturbance control ratio: Tj = Y (jw) 1 D(jw) 1 + L(jW) (33) where R(joS) is related to the reference signal, F(j®) represents the prefilter, and the constraint functions ^5sup(®) and ¿5infp) are the transfer functions TrU(^) and ^(s) defined in equations (29) and (30), respectively. The subset of analysis frequencies Q5 is {0.01, 0.1, 1, 5, 10, 50, 100, 150} rad/s. After the bound generation and intersection, the loop shaping process is carried out to obtain the controller function (Fig. 5) [20]. -270 -225 -180 -135 -90 Open-loop phase [deg] Fig. 5. Loop shaping of the controller The resulting controller is: G(s) = (0.001s + 0.008)/s. (34) The synthesis of the prefilter is similar to that of the controller, resulting in: F (s) = 910 /(s2 + 83s + 910). (35) The hydraulic force control system, including the prefilter and controller, was implemented in Simulink, using the block diagram shown in Fig. 6. QFT-based control Hydraulic System Described in Section 4 Described in Section 2.1 QFT QFT Prefilter Controller u Ref F(s) G(s) ^vAB ^ F Valve Cylinder Output Force Fig. 6. Block diagram of the hydraulic force control system 5 EXPERIMENTAL AND SIMULATION RESULTS Experimental tests and simulations were carried out using the selected HVE hose according to Section 3.3. A force-tracking reference was defined as several force steps of different magnitudes (Fig. 7a). To show that the model used for simulation has a good matching with the experimental results, an enlarged plot can be observed in Fig. 7b. The hydraulic force (FH), chamber pressures (pA and pB), and the control signal sent to the valve (Uc) also demonstrate the good representativeness of the nonlinear dynamic model. In both simulation and experiment, the QFT-based prefilter and controller were used. The systems parameters are those presented in the Appendix. The advantages of using HVE hoses for adding compliance to the system for improving the force control performance can be observed when comparing with the system using rigid pipes. The controller G(s) and prefilter F(s) were designed for the case in which rigid pipes are used, assuming the same performance specifications applied to the system with hoses, resulting: and G(s) = (0.0005s + 0.002) / s, (36) F (s) = 484 / (s2 + 44s + 484). (37) Fig. 8 shows the system force simulation responses using pipes and hoses. A 4 kN step reference force was applied at 2 s. The disturbance input was a filtered step, used to reproduce a load movement disturbance, with 5 mm of magnitude at 3 s. 13 14 15 16 17 18 19 20 21 22 23 b) Time [s] Fig. 7. a) Comparison between experimental and simulated force tracking responses using HVE hoses, and b) enlarged plot of region A and hydraulic force, chamber pressures and control signal curves Time Fig. 8. Comparison between force responses using pipes or HVE hoses The step response using rigid pipes was more oscillatory in comparison to the system with HVE hoses. The resulting overshoots were 3.87 % and 2.15 % and settling times (?S1%) equal to 0.33 s and 0.20 s, respectively. The settling time of the system with hoses is between the time domain specifications presented in section 4; however, the overshot is higher than 1%. The system with pipes was not able to strictly achieve the specifications using the designed QFT controller. It is necessary to observe that the QFT controller was designed based on the linear model considering the parameter uncertainties. Therefore, some deviation from the specifications can be expected. The main drawback of a force control system is the disturbance rejection. As can be seen in Fig. 8, by reducing hydraulic stiffness through the use of HVE hoses, the output impedance is decreased, such that the load displacement causes small changes in the applied force. The control signal has a higher amplitude and is more oscillatory when using pipes, as expected, since there is not a hydraulic or mechanical compliant element. 6 CONCLUSIONS A study on a hydraulic force control system using high volumetric expansion hoses to increase the hydraulic compliance was carried out. The major achievements of this study are to demonstrate the possibility of adding compliance hydraulically in order to perform the force control without using a spring and to establish an analytical procedure for selecting the hydraulic hoses necessary to achieve this goal. The inclusion of more hydraulic capacitive components allows a reduction in the effective bulk modulus and makes the system more compliant. Based on the presented results, it was concluded that the force response performance of the system using hydraulically compliant components is improved in comparison to a system with rigid pipes. Moreover, the use of the QFT-based linear controller, applied in a nonlinear system, resulted in suitable performance and disturbance rejection capabilities. Since the reduction in the hydraulic stiffness makes the system more compliant and conservative, the complexity of the controller is also reduced giving the perspective of the use of other types of controllers, including the classical PID. Future work should include developing a selection procedure for accumulators with a focus on force control and series hydraulic damping analysis aimed at improving the force control performance of a PHEA-based system. 6 ACKNOWLEDGEMENTS Special thanks to CNPq by the financial support and to M.S. Craig Borghesani from Terasoft, Inc. for providing the QFT Toolbox for this study. 7 REFERENCES [1] Wang, J., Zhang, G., Zhang, H., Fuhlbrigge, T. (2008). Force control technologies for new robotic applications. IEEE International Conference on Technologies for Practical Robot Applications, p. 143-149, D0l:10.1109/TEPRA.2008.4686689. [2] Parzer, H., Gattringer, H., Müller, A., Naderer, R. (2016). Robot force/position control combined with ILC for repetitive high speed applications. International Conference on Robotics in Alpe-Adria Danube Region, p. 12-19, D0I:10.1007/978-3-319-49058-8_2. [3] Alleyne, A., Liu, R. (1999). On the limitations of force tracking control for hydraulic servosystems. Journal of Dynamic Systems, Measurement, and Control, vol. 121, no. 2, p. 184190, D0I:10.1115/1.2802453. [4] Whitney, D. (1985). Historical perspective and state of the art in robot force control. IEEE International Conference on Robotics and Automation, p. 262-268, D0I:10.1109/ R0B0T.1985.1087266. [5] Calanca, A., Muradore, R., Fiorini, P. (2016). A review of algorithms for compliant control of stiff and fixed-compliance robots. IEEE/ASME Transactions on Mechatronics, vol. 21, no. 2, p. 613-624, D0I:10.1109/TMECH.2015.2465849. [6] Udai, A.D., Hayat, A.A., Saha, S.K. (2014). Parallel active/ passive force control of industrial robots with joint compliance. IEEE/RSJ International Conference on Intelligent Robots and Systems, p. 4511-4516, D0I:10.1109/IR0S.2014.6943201. [7] Boaventura, T. (2013). Hydraulic Compliance Control of the Quadruped Robot HyQ. Ph.D. Thesis, University of Genova, Genova. [8] Robinson, D.W. (2000). Design and Analysis of Series Elasticity in Closed-loop Actuator Force Control. Ph.D. Thesis, Massachusetts Institute of Technology, Massachusetts. [9] Lee, Y.-F., Chu, C.-Y., Xu, J.-Y., Lan, C.-C. (2016). A humanoid robotic wrist with two-dimensional series elastic actuation for accurate force/torque interaction. IEEE/ASME Transactions on Mechatronics, vol. 21, no. 3, p. 1315-1325, D0I:10.1109/ TMECH.2016.2530746. [10] Losey, D.P., Erwin, A., McDonald, C.G., Sergi, F., O'Malley, M.K. (2016). A time domain approach to control of series elastic actuators: adaptive torque and passivity-based impedance control. IEEE/ASME Transactions on Mechatronics, vol. 21, no. 4, p. 2085-2096, D0I:10.1109/TMECH.2016.2557727. [11] Xu, Y., Hollerbach, J.M., Ma, D. (1995). A Nonlinear PD controller for force and control transient control. IEEE Control Systems, vol. 15, no. 1, p. 15-21, D0I:10.1109/37.341859. [12] Alleyne, A., Liu, R. (2000). A simplified approach to force control for electro-hydraulic systems. Control Engineering Practice, vol. 8, no. 12, p. 1347-1356, D0I:10.1016/S0967-0661(00)00081-2. [13] Ahn, K.K., Truong, D.Q., Thanh, T.Q., Lee, B.R. (2008). Online self-tuning fuzzy porportional-integral-derivative control for hydraulic load simulator. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, vol. 222, no. 2, p. 81-95, D0l:10.1243/09596518JSCE484. [14] Nakkarat, P., Kuntanapreeda, S. (2009). Observer-based backstepping force control of an electrohydraulic actuator. Control Engineering Practice, vol. 17, no. 8, p. 895-902, D0I:10.1016/j.conengprac.2009.02.011. [15] Plummer, A.R. (2007). Robust electrohydraulic force control. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, vol. 221, no. 4, p. 717-731, D0I:10.1243/09596518JSCE370. [16] Lamming, C.P.G., Plummer, A.R., Hillis, A.J. (2010). Analysis of robust electrohydraulic force control. 7th International Fluid Power Conference, p. 1-12. [17] Wells, D.L., Iversen, E.K., Davis, C.C., Jacobsen, S.C. (1990). An investigation of hydraulic actuator performance tradeoffs using a generic model. IEEE International Conference on Robotics and Automation, p. 2168-2173, D0I:10.1109/ R0B0T.1990.126325. [18] Petersen, N.R. (2002). Loading Assembly Having a Soft Actuator. U.S. Patent 6.457.369 B1, World intellectual property organization, Geneva. [19] Zoppi, M. (2013). Method for Adapting Stiffness in a Variable Stiffness Actuator. U.S. Patent 2013/0047596 A1, World intellectual property organization, Geneva. [20] Houpis, C.H., Rasmussen, S.J. (1999). Quantitative Feedback Theory: Fundamentals and Applications. Marcel Dekker, New York. [21] Pratt, G.A., Williamson, M.M., Dillworth, P., Pratt, J., Wright, A. (1997). Stiffness isn't everything. 0. Khatib, J. K. Salisbury (eds.), Experimental Robotics IV, Springer, Berlin, p. 253-262, D0I:10.1007/BFb0035216. [22] M00G (2007). 760 Series - Servovalves - ISO 10372 Size 04. MOOG. East Aurora, New York Catalog CDL6335 Rev G 500213 807. [23] Szpak, R., Ramos Filho, J.R.B., De Negri, V.J. (2010). Theoretical and experimental study of the matching between proportional valves and symmetric and asymmetric cylinders. 7th International Fluid Power, p. 155-166. [24] Gomes, S.C.P., da Rosa, V.S. (2003). A new approach to compensate friction in robotic actuators. IEEE International Conference on Robotics and Automation, p. 622-627, D0I:10.1109/R0B0T.2003.1241663. [25] EATON Co. (2008). Synflex Hose and Fittings Master Catalog, Eden Prairie, Catalog E-H00V-MC002-M. [26] Parker Hannifin Corp. (2012). Parflex® Thermoplastic & Fluoropolymer Products - Hose, Tubing, Fittings & Accessories, Ravenna, Catalog CAT 4660. [27] Johnston, D.N., Way, T.M., Cone, K.M. (2010). Measured dynamic properties of flexible hoses. Journal of Vibration and Acoustics, vol. 132, no. 2, p. 021011-1-021011-8, D0I:10.1115/1.4000774. [28] Borghesani, C., Chait, Y., Yaniv, 0. (2003). The QFT Frequency Domain Control Design Toolbox: For use with Matlab. User's Guide. Terasoft. [29] IS0 10770-1:2009. Hydraulic fluid power - Electrically modulated hydraulic control valves - Part 1: Test methods for four-port directional flow-control valves. International Organization for Standardization, Geneva. 8 APPENDIX Table 3. Technical data used for simulations Hydraulic power unit pS = 7x106 Pa pT = 0 Pa ße = 1.2x109 Pa Model: BOSCH REXROTH CGT3 MS2 50/22-500/Z1X/B1 Double rod and double L = 0.5 m acting Au = 1.5833x10-3 m2 cylinder m = 13.1 kg Model: MOOG 760 C263-A 4/3 closed center servovalve Kvp*= 3x10-7 m3/(s-Pa05)@70x105 Pa Kvinp*= 2.405x10-9 m3/(s-Pao.5)@70 x105 Pa Un = 10 VDC wnv = 1099 rad/s @ ±40 % of opening valve ¿v = 0.9 @ ± 40 % of opening valve High Model: EATON Synflex ® 3130-08 volumetric Dho = 12.7 x 10-3 m expansion E = 1.56 x 10-7 m3/m @ 70 x 105 Pa hose Lho = 1.5 m Pipes Dp: 10.5 x 10-3 m Lp: 0.5 m Model: HBM U2AD1-1t Load cell Maximum capacity: 9800 N KS = 9.8 x 107 N/m ' Obtained through experimental tests Strojniški vestnik - Journal of Mechanical Engineering 64(2018)10, 590-600 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2018.5371 Original Scientific Paper Received for review: 2018-03-26 Received revised form: 2018-06-08 Accepted for publication: 2018-06-29 Heating by Flat Plate Collector Combined with a Heat Pump Simon Marčič1 - Sebastijan Seme1 - Julijan Jan Salamunic1 - Zdravko Praunseis1 - Stanislav Pehan2 1 University of Maribor, Faculty of Energy Technology, Slovenia 2 University of Maribor, Faculty of Mechanical Engineering, Slovenia The aim of the research was to establish the efficiency of the combination flat plate collector-heat pump versus ambient air-heat pump. The extensive measurements were carried out and analyzed of the flat plate collector alone and combination flat plate collector-heat pump. In the beginning, the heat pump introduction into the serial combination with flat plate collectors was analyzed. The appropriate test rig with measurement equipment was designed and built. The test rig base consisted of the flat plate collector, heat pump, and 3001 boiler that was a replacement for a house's heating. Detailed temperature measurements were performed during rainy, sunny and cloudy weather in summer and winter. To estimate the efficiency of heating, the consumed electric energy versus heat is the universal criterion that is derived from the Coefficient of performance (COP), which is calculated on the basis of measured temperatures. The comparison of the calculated COP for the combination flat plate collector-heat pump and the air-water heat pump shows that the largest increase of the COP of the combination flat plate collector-heat pump in comparison with the air-water heat pump is on sunny winter days, and that it amounts to 38 %. In rainy and cloudy weather, the COP of the flat plate collector-heat pump combination is lower than the air-water heat pump COP The presented analyses clearly indicated that use of the flat plate collector is limited. Keywords: flat plate collector, heat pump, the Coefficient of performance, test rig Highlights • The test rig for combination flat plate collector-heat pump was built to prove the extended use of flat plate collectors. • The heating tests were carried out at various weather conditions and at different sun illuminations in summer and winter days. • On the basis of the temperature measurements, the Coefficients of performance of the combination flat plate collector-heat pump and ambient air-heat pump were calculated and compared. • According to research, the combination flat plate collector-heat pump makes sense mostly on a sunny winter day. 0 INTRODUCTION Recently, many houses have been equipped with solar collectors in order to ensure the sanitary hot water and heat. Nowadays, the heating systems based on heat extraction from the ambient air have become more efficient and, in particular, cheaper. The question arises as to whether the systems with solar collectors should be really entirely written-off, or whether they could, nevertheless, be upgraded sensibly and effectively. The basic goal of the presented research work was to find a way to raise the heating capacity of a heating system that is based on solar collectors [1] and [2]. The second goal was to present the theoretical reasons that would convince the owners and the investors to preserve, and even upgrade, the existing house heating systems with solar collectors to provide a longer heating season. The key component of active solar heating systems are the solar collectors that gather the sun's energy, transform this into heat, and then transfer that heat to a fluid, usually water. Mainly two types of solar collectors are used presently: Flat plate collectors [3] and [4] and evacuated tube collectors [5] to [7]. Flat plate collectors usually consist of a glass cover and metal plate absorber with the integrated pipes and housing. The solar energy absorbed by the metal plate absorber is transferred into the liquid flowing within the pipes. The glazing on the upper surface and the insulated housing bottom limit the thermal losses. In the flat plate collector, in the gap between the glazing and the metal plate absorber, the air is present. From the metal plate absorber, this air is warmed, then heats the glazing and becomes the source of losses. To avoid this harmful phenomenon, the so-called evacuated tube collectors have been introduced. They consist of glass evacuated tubes, through which the absorption pipes are installed [8] and [9]. The evacuated tube collectors are more effective but, in comparison with heat pumps, both types of collectors alone in the long-term are losing heat. For domestic heating and domestic hot water gathering a heat pump alone [10] to [12] is used increasingly. They are much more effective than simple electrical resistance heaters using the same amount of electricity [13] to [14]. Heat pumps are designed to move the thermal energy in the opposite direction of spontaneous heat transfer by absorbing heat (anergy) from a cold space and releasing it to a warmer one. To 590 *Corr. Author's Address: University of Maribor, Faculty of Mechanical Engineering, Smetanova 17, Maribor, Slovenia, stanislav.pehan@um.si run the compressor, a heat pump uses a small amount of external electric energy (exergy). The source of anergy can be any surrounding media; air, water etc. To describe the ratio of the heat flow (anergy plus exergy) and external electric energy (exergy), the term Coefficient of performance (COP), is used. Theoretically, the COP increases when the temperature difference between the cold space and the warmer one decreases. That fact encourages the heating designer to choose low-temperature underfloor house heating rather than high-temperature radiators [15] and [16]. Conventionally, all types of collectors and heat pumps are used separately. In the paper, it is a suggestion to join both systems in a serial connection. Heated water that is coming from the collector is led directly into the input side of a heat pump to gain the higher temperatures. Hypothetically, this serial connection should ensure higher temperatures of heating as each heating system works independently. To verify the effectiveness of the serial connection flat plate collector - heat pump, a convenient test rig was built on a reduced scale. Measurements of suitable quantities were carried out and analyzed. The test rig consists of a flat plate collector, water-water heat pump, and boiler. The boiler is a replacement for house heating. The water in the boiler accumulates the heat that is produced over the day. The tests were carried under various weather conditions in the summer and winter seasons. Inlet and outlet temperatures from the flat plate collector and heat pump were measured, as well as the ambient air temperature and the water temperature in the boiler. The paper comprises the testing of the flat plate collector alone and combination flat plate collector-heat pump. The purpose of this test is to establish the efficiency of each system. Comparisons of both systems were carried out. 1 HEATING BY FLAT PLATE COLLECTOR The dedicated test facility is shown in Fig. 1. On the building roof, several flat plate collectors are positioned oriented to the south. The flat plate collector positioned on the left side of the roof that is tested has an active absorber area of two square meters. The boiler, capacity 300 liters, is located in the building. The flat plate collector heats the boiler water directly by adequate heat exchanger built in the boiler. Basic flat plate collector characteristics are follows: 2 square meter active surface, copper absorber with special selective layer, crystal prism glass, housing made from anodized aluminum, rock wool isolation, four pipe connections and steam flap. The measurement schema is shown in Fig. 2. Fig. 1. Test facility Heat exchanger Fig. 2. Measurement schema of heating boiler water by flat plate collector Measurements, Figs. 3 and 4, show that in cloudy and rainy weather, the efficiency of the flat plate collector is low. In cloudy weather, the water temperature in the boiler T5 is mostly lower than the temperature of the ambient air during the measurement process. In cloudy weather, the contribution of diffusion light to effective heating is low. This phenomenon is in accordance with the measurements of the illumination, which was, in perfectly clear weather, from 90,000 Lux to 95,000 Lux, and in the cloudy from 3,000 Lux to 5,000 Lux. 2 HEATING BY SERIAL CONNECTION FLAT PLATE COLLECTOR-HEAT PUMP Flat plate collector is in serial connection with the heat pump water-water, 1.5 kW under standard measurement conditions, Fig. 5. A comparison of heating with a flat plate collector and a combination flat plate collector-heat pump shows the obvious advantage of the latter, Fig. 6. With the system, even in the most unfavorable weather Fig. 3. Temperatures in the boiler by heating with flat plate collector in July 2015 Fig. 4. Temperatures in the boiler by heating with flat plate collector in August 2015 Heat pump Fig. 5. The setup of the boiler water heating by the combination flat collector-heat pump conditions, high temperatures can be achieved in the about 25 °C, and the combination is from 18 °C to boiler. On a rainy day, for example, only flat plate about 39 °C. collectors heat the water in a boiler from 19 °C to A 50 40 30 Temperature 20 [°C] (Measured at . „ 15:40 p.m.) 0 49 888KB , Wat er ten perat ure in B888H 39 29 23 T- 25 18 Part ly clo udy 15 — VR ainy \ — Rainy L ^ainy Partly cloudy 30. 31. 01. 02. 03. Day in July and August 2015 Ambient temperatures Partly cloudy Rainy Fig. 6. Temperature in the boiler as a result of heating with flat plate collector-heat pump Fig. 7. Temperatures in the heating system, of the flat collector and the heat pump, on a rainy summer day 3 FLAT PLATE COLLECTOR-HEAT PUMP VERSUS AMBIENT AIR-HEAT PUMP To understand the functioning of the flat plate collector-heat pump, it is worthwhile to know more precisely the generation of heat energy over the entire current day. Therefore, temperature measurements were carried out at the characteristic points of the heating system, every 10 minutes, from 9:30 to 15:30. The temperature of the water at the inlet and at the outlet from the flat plate collector and at the inlet and at the outlet of the heat pump, the boiler temperature, and the ambient air temperature are measured. In the following figures, Figs. 7 to 11, the measurements of temperatures are presented in three summer and two winter days, in order to see the operation of the heating system in extreme cases. From Fig. 7 it follows that, on a rainy summer day, the ambient air temperature is higher than the temperatures in the flat plate collector, as a result of the low illumination. The contribution of the diffusion light is too small to warm the water at a higher temperature than the ambient air. At the beginning of the measurement, at 9:30, the temperatures T1 and T2 are higher than the ambient air temperature because the system is stationary and there is no water flow. In the stationary state, diffusion light radiation can heat the water above the temperature of the ambient air. At the beginning of the measurement, Fig 8, when there is no water flow in the system, the inlet and the outlet temperatures T and T2 are higher than Fig. 8. Temperatures in the heating system, of a flat plate collector and heat pump, on a partly cloudy summer day Fig. 9. Temperatures in the heating system of a flat plate collector and heat pump, on a sunny summer day the ambient air by14 °C. By noon the cloudy weather prevailed, after that the weather became mostly clear. By 12:00 the outlet water temperature T2 is equal to or two degrees Celsius higher than the ambient air. After 12:00, when the weather became mostly clear, the temperatures T1 and T2 exceeded the ambient air temperature. The largest surplus is 7 °C at 15:15. Fig. 9 shows the course of temperatures in the misty sunny weather when the illumination did not exceed 55,000 Lux. After the circulation pumps switched on, the outlet water temperature T2 was equated with the temperature of the ambient air and, till 14:30, remained practically the same. After the sunlight dispersed the mist, the temperatures T1 and T2 became higher than the ambient air Tam by about 5 °C. Fig. 10 shows the course of temperature and brightness on December 29th, 2015. Illumination varied between 49,000 Lux and 63,000 Lux. In the stationary state, when there is no flow, water (a mixture of glycol and water) in the flat plate collector was warmed up to 18 °C higher temperature than the ambient air. After the circulation pumps were switched on, the temperatures T1 and T2 were equalized to the temperature Tam. As the illumination increases and, above all, because of the rise of the Sun above the horizon, the incidence angle of light towards the flat plate collector increases, so more light falls on the flat plate collector. The highest difference between the outlet water temperature T2 and the ambient air temperature Tam was 9 °C, at 12:30, when the Sun Fig. 10. The ambient illuminance and the temperatures in the heating system of a flat plate collector and heat pump, on a sunny winter day Fig. 11. The ambient illuminance and the temperatures in the heating system, of a flat plate collector and heat pump, on a cloudy winter day was in the highest position, and then the temperature difference began to decrease, and it was, at 15:30, only 2 °C. At the beginning of the measurements, Fig. 11, at 9:30, when the circulation pump was stationary and the illumination was low, the temperatures T1, T2, and Tam were the same. The inlet T3 and the outlet temperature T4 of the heat pump were the same as the water temperature in the boiler T5. After the circulation pumps were switched on, the outlet water temperature from the flat plate collectors T2 is reduced by an average amount of 10 °C. Lowering the temperature was due to the lower illumination. The flat plate collector received the heat from the diffusion light, the portion of which was low, and from the heat transfer, which was low also due to the good insulation. On the basis of the temperature measurements and exergetic efficiency given by the heat pump manufacturer, COP is calculated using the following equation: Fig. 12. The COP and energies in the heating system with the heat pump on a rainy summer day Fig. 13. The COP and energies in the heating system with the heat pump on a partly cloudy summer day COP = T -*- A T - T 4 2 (1) where T4 is the heat pump outlet water temperature and T2 is the temperature the of the heating anergy. In the case of the air-water heat pump, the COP is calculated by the following equation COP = - T4 - Tan ■Š, (2) where Tam stands for ambient temperature. The purpose of the flat plate collector-heat pump combination testing is to establish under what weather conditions this combination is better than a conventional air-water heat pump. Heat generated by both systems is calculated using the equation Q = m ■ Cw -(T5 - T5int ), (3) where Q is the generated heat, m is the mass of water in the boiler, cw is specific water heat, T5 is the actual temperature of the water in the boiler, and T5int is the initial water temperature in the boiler (18 °C). Eq. 4 Fig. 14. The COP and energies in the heating system with the heat pump on a sunny summer day Fig. 15. The COP and energies in the heating system with the heat pump on a sunny winter day (3) is used to calculate the actually generated heat accumulated in the boiler. Subsequently, by means of the COPs calculated by Eqs. (1) and (2) the heating exergy or electric energy consumed by both heat pumps is calculated. Fig. 12 shows the calculated COPs for flat plate collector-heat pump combination and air-water heat pump. Given that, in rainy summer weather, the outlet temperature from the flat plate collector T2 is lower than the ambient air temperature Tam, the COPs for the flat plate collector-heat pump combination are lower than the COP of the air-water heat pump. As a result, the heating exergy of the flat plate collector-heat pump combination is higher than the air-water heat pump heating exergy. In other words, to heat the water to the same temperature T5, the flat plate collector-heat pump combination consumes more electric energy than the air-water heat pump. If the weather is changeable or partly cloudy, Fig. 13 shows that the COP of the flat plate collector-heat pump combination is better than of the air-water heat pump by 3 % at 11:30 and by 33 % at 15:30. Moreover, the consumption of electric energy or exergy of the flat plate collector-heat pump combination is correspondingly lower than that of the air-water heat pump. If it is a hazy and sweltering sunny day, Fig. 14, the COP of the flat plate collector-heat pump combination is, on average, higher by 4 % between 9:30 and 14:30 than of the air-water heat pump, but at 15:00 it rises to 19 %. The consumption of the electric energy of the flat plate collector-heat pump combination decreases Fig. 16. The COP and energies in the heating system with the heat pump on a cloudy winter day in the same portion. The increase of COP and lower electric energy consumption result from the clearer weather. The sun dispersed the mist covering the sky until then. The greatest improvement of the COP of the flat plate collector-heat pump combination in comparison with the air-water heat pump occurs in sunny winter weather, Fig. 15, and amounts to as much as 38 % at 12:30. The smallest improvement is at 15:30 and amounts to 17 %. Such improvements result from a low ambient air temperatures Tam and, consequently, from the low proportion of the heating energy in the heating thermal energy. In winter months when the ambient air temperatures Tam are low, the air-water heat pump must overcome wider temperature differences to achieve the same temperature of the heated water as on summer days, resulting in low COPs. Likewise, the consumption of electric energy of the flat plate collector-heat pump combination is reduced in comparison with the air-water heat pump. On cloudy winter days when the illumination is low, Fig. 16, the COP of the flat plate collector-heat pump combination is worse than the COP of the air-water heat pump. The flat plate collector fails to heat the water that is circulating in the collector-evaporator system, Fig. 5, of the heat pump of the flat plate collector-heat pump combination, to a temperature higher than the ambient air temperature Tam. The flat plate collector receives heat from diffusion light, which is very low. Consequently, the heat transfer by convection from the ambient air to the absorption flat plate collector is very low also. The electric energy consumption of the flat plate collector-heat pump combination is also increased accordingly in comparison with the air-water heat pump. 4 RESULTS AND DISCUSSION On sunny summer days, the solar collector heated the boiler water to a temperature from 34 °C to 36 °C. In cloudy and rainy weather, the boiler water temperature was always lower than the ambient air temperature and amounted to 19 °C in rainy weather and to a temperature from 23 °C to 27 °C in cloudy weather. In sunny weather, water is heated to 51 °C by using the flat plate collector-heat pump combination, up to 29 °C in partly cloudy weather, and to 19 °C in rainy weather. In sunny and partly cloudy weather, the flat plate collector-heat pump combination achieves considerably higher temperatures than the solar collectors. Figs. 7 to 11 show the trend of characteristic temperatures between 9:30 and 15:30 on different summer and winter days and under various weather conditions. The ambient air temperature Tam and water outlet temperature from collectors T2 are the most important temperatures. In sunny weather, T2 is always higher than Tam. The biggest difference between these two temperatures is on a sunny December day, and it amounts to 9 °C. In rainy and cloudy weather, temperature T2 is lower than Tam. The biggest undercooling is in January and amounts to 11 °C. The comparison of the calculated COPs for the flat plate collector-heat pump combination and the heat pump air-water system shows that the largest increase of the COP of the flat plate collector-heat pump combination in comparison with the air-water heat pump is on sunny winter days and that it amounts to 38 %, Fig. 15. In rainy and cloudy weather, the COPs of the flat plate collector-heat pump combination are lower than the air-water heat pump COPs. The consumption of electric energy follows the COPs strictly, Table 1. Table 1. Energy - obtained and consumed Summer Winter Sunny Partly cloudy Rainy Sunny Cloudy Energy in boiler [kWh] 11.498 11.149 8.362 4.878 2.787 Heat pump, energy from ambient air [kWh] 1.543 1.585 1.190 1.223 0.588 Average COP 7.5 7.0 7.0 4.0 4.7 Heat pump, energy from flat plate collector [kWh] 1.501 1.416 1.772 1.004 0.827 Average COP 7.7 7.9 4.7 4.9 3.4 5 CONCLUSIONS The main novelty of the work is the clarification about the usefulness of the flat plate collectors for most effective and competitive domestic heating. In short, these devices are for write off. This radical statement is based on quite comprehensive measurements of flat plate collectors functioning at different weather condition in summer and winter as the heat provider, that bringing the sun energy in form of heat to the heat pump. The competitors of them have been the ambient air, as it is. Normally, as the majority of the authors are dealing with, the bivalent system of flat plate panel and the heat pump was designed. In order to improve the efficiency of the heating as best that can be done, in this research, the serial connection of flat plate collector and the heat pump is made and tested. Due to the allegedly higher temperature input into the heat pump, from the serial connection hypothetically the evident improvements of the heating has been expected. Instead of dealing with the real house heating system, in the experiment the daily collected heat was gathered by using the water boiler, the capacity of 300 liters, that proved to be the very useful and practical idea. The most obvious advantage was the relatively simple test rig design, which does not include any unnecessary systems. An additional advantage was the simplicity of energy measurement that is carried out by using the thermometers only. All the rest was done by simple computations, which could be easy to control. The major disadvantage of these simple experiments represents the absolute values of COPs. In comparison with the generally expected values, COPs obtained in the presented tests, were suspiciously high. This phenomenon was not due to the excellence of the system components, more due to the relatively low-temperature differences of the boiler water. Even the experts should be warned not to rely on one individual data of the heating system but to try to evaluate the system in whole. It is shown that at rare distinctively sunny days the COPs reach unexpected high values. Considering the fact, that sunny day is the exception rather than the rule, in all seasons, mostly partly cloudy weather prevails. In such transient weather conditions, the heat pump, that works after definition with the simplest air/fluid heat exchanger, is the best solution as it is proved. 6 REFERENCES [1] Pandey, K.M., Chaurasiva, R. (2017). A review on analysis and development of solar flat plate collector. Renewable and Sustainable Energy Reviews, vol. 67, p. 641-650, D0l:10.1016/j.rser.2016.09.078. [2] Tian, Y., Zhao, C.Y. (2013). A review of solar collectors and thermal energy storage in solar thermal applications. Applied Energy, vol. 104, p. 538-553, D0I:10.1016/j. apenergy.2012.11.051. [3] Powell, K.M., Edgar, T.F. (2012). Modeling and control of a solar thermal power plant with thermal energy storage. Chemical Engineering Science, vol. 71, no. 27, p. 138-145, D0I:10.1016/j.ces.2011.12.009. [4] Dubey, S., Tiwari, G.N. (2009). Analysis of PV/T flat plate water collectors connected in series. Solar Energy, vol. 83, no. 9, p. 1485-1498, D0I:10.1016/j.solener.2009.04.002. [5] Yang, J., Jiang, Q., Hou, J., Luo, C. (2015). A study on thermal performance of a novel all-glass evacuated tube solar collector manifold header with an inserted tube. International Journal of Photoenergy, vol. 2015, Article ID 409517, D0I:10.1155/2015/409517. [6] Kaya, H., Arslan, K., Eltugral, N. (2018). Experimental investigation of thermal performance of an evacuated U-Tube solar collector with ZnO/Etylene glycol-pure water nanofluids. Renewable Energy, vol. 122, p. 329-338, D0I:10.1016/j. renene.2018.01.115. [7] Sabiha, M.A., Saidur, R., Mekhilef, S., Mahian, 0. (2015). Progress and latest developments of evacuated tube solar collectors. Renewable and Sustainable Energy Reviews, vol. 51, p. 1038-1054, D0I:10.1016/j.rser.2015.07.016. [8] Morrison, G.L., Tran, N.H., McKenzie, D.R., Onley, I.C., Harding, G.L., Collins, R.E. (1984). Long-term performance of evacuated tubular solar water heaters in Sydney, Australia. Solar Energy, vol. 32, no. 6, p. 785-791, D0I:10.1016/0038-092X(84)90253-6. [9] Morrison, G.L., Budihardjo, I., Behnia, M. (2005). Measurement and simulation of flow rate in a water-in-glass evacuated tube solar water heater. Solar Energy, vol. 78, no. 2, p. 257-267, D0I:10.1016/j.solener.2004.09.005. [10] Chua, K.J., Chou, S.K., Yang, M.W. (2010). Advances in heat pump systems: A review. Applied Energy, vol. 87, no. 12, p. 3611-3624, D0I:10.1016/j.apenergy.2010.06.014. [11] Oprešnik, M. (1969). The influence of the different working means properties on the exergy loss due to the compressor heat pump processes irreversibility. (Vpliv lastnosti različnih delovnih sredstev na izgube eksergije zaradi nepovračljivosti v kompresorju pri procesih toplotne črpalke). Strojniški vestnik - Journal of Mechanical Engineering, vol. 15, no. 4-5, p. 113116. (in Slovene) [12] Sakellari, D. Lundqvist, P. (2001). Towards energy-use optimization of a domestic heating system based on a heat pump. Strojniški vestnik - Journal of Mechanical Engineering, vol. 47, no. 8, p. 506-511. [13] Haložan, H. (2000). Heat-Pumping Technologies. Strojniški vestnik - Journal of Mechanical Engineering, vol. 46, no. 7, p. 445-453. [14] Hepbasli, A., Kalinci, Y. (2009). A review of heat pump water heating systems. Renewable and Sustainable Energy Review, vol. 13, no. 6-7, p.1211-1229, D0I:10.1016/j. rser.2008.08.002. [15] Smith, D.C., Elmore, A.C. (2018). The observed effects of changes in groundwater flow on a borehole heat exchanger of a large scale ground coupled heat pump system. Geothermic, vol. 74, p. 240-246, D0I:10.1016/j.geothermics.2018.03.008. [16] Jeong, J., Hong, T., Kim, J., Chae, M., Ji, C. (2018). Multi-criteria analysis of a self-consumption strategy for building sectors focused on ground source heat pump systems. Journal of Cleaner Production, vol. 186, p. 68-80, D0I:10.1016/j. jclepro.2018.03.121. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)10, 601-610 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2018.5242 Original Scientific Paper Received for review: 2018-01-26 Received revised form: 2018-05-31 Accepted for publication: 2018-06-13 Research of Crack Detection and Delay Crack Propagation on a Nonlinear Rotor Jun Liu12* - Weilong Liu1 - Xiaofeng Wang12 1 Tianjin Key Laboratory of the Design and Intelligent Control of the Advanced Mechanical System, China 2 National Demonstration Center for Experimental Mechanical and Electrical Engineering Education, China The existence and propagation of a crack on an engine's nonlinear rotor can easily cause a serious accident, so it is quite essential to detect a crack in the rotor system through analysis of the rotor's vibration signals and to delay the crack propagation. The vibration signals of the cracked rotor show the non-periodic and the nonlinearity within a transient response. The fault characteristics of the crack can be easily hidden in other vibrational components and are difficult to extract from transient signals. In the paper, the Hilbert-Huang marginal spectrum (HHMS) method is proposed to realize the extraction of crack's fault characteristics in a nonlinear rotor system. In addition, the crack's opening and closing degrees at any rotating position are also put forward representing the crack propagation, and some tactics of the delaying crack propagation are proposed based on the effect on the crack's opening and closing degree caused by the change of different rotor parameters. The simulation results show that the HHMS method can detect a crack early at about five per cent of the depth of the rotor's diameter, and the methods of delaying crack propagation can effectively delay propagation. Experimental results verify the effectiveness of the proposed HHMS method. Keywords: cracked rotor, transient response, Hilbert-Huang marginal spectrum, crack detection, delaying crack propagation Highlights • The crack fault can be effectively extracted by the HHMS method from complex vibration signals. • The effectiveness of the proposed HHMS method was verified by the experimental and simulation results. • Increasing the rotational acceleration can delay the crack propagation. • When the angle 3 and the support of the rotor are adjusted effectively, the delay crack propagation can be achieved, such as reducing the symmetric nonlinear term 3(0) and increasing the value of the asymmetric nonlinear term ec(1) and es<1). 0 INTRODUCTION Engines are essential to the aircraft, and the rotor system is the key component of the engine. The rotor system usually works in difficult environments with multiple coupling fields, which include high temperatures, high pressures, high loads, and periodic alternating loads. Such an environment makes the rotor system generate micro-cracks and micro-defects, and the defects could cause severe accidents when the crack propagates. In recent decades, researchers have identified cracks through the analysis of vibration signals of a rotor system, and this research has supplied some significant achievements. The Fourier transform [1] and [2] can be used to transfer the periodic signal from the time domain to the frequency domain, and engineers use the frequency characteristics of the detected signals to analyze any faults. However, it is only able to deal with periodic signals and can not process non-periodic signals. The short-time Fourier transform (STFT) [3] and [4] can obtain the time-varying spectrum of the detected signal by imposing a time window on the detected signal and then sliding the time window to do the Fourier transform with data of the time window. However, a contradiction of the time domain resolution and frequency domain resolution exists; therefore, the resolution accuracy could be affected. The Winger-Ville distribution [5] and [6] has a good time-frequency clustering, and it is a quadratic time-frequency method. When it is used to analyse the multi-component signals [7], it will produce serious cross-item interference and a false frequency that does not exist in the actual signal. The wavelet transform [8] and [9] can express local characteristics in the time-frequency domain. Because it needs to select different parameters of the wavelet function according to different detected signals, the analysis results depend heavily on the selecting parameters of the wavelet function. The Hilbert-Huang transform (HHT) [10] and [11] is a novel signal analysis method on the time-frequency domain, and the marginal spectrum of the signal can be obtained by integrating the Hilbert spectrum [12] to [14], which can accurately reflect the change of the signal amplitude on the frequency domain. The HHT has already been applied to voice recognition, bearing fault detection, and other fields. An aircraft's engine's rotor often works under the transient conditions when the aircraft takes off or lands, performs stunts, or follows a climbing-diving flight model [15]. The vibration signal of the cracked rotor system exhibits the non-periodic and nonlinear *Corr. Author's Address: Tianjin University of Technology, No. 391 Binshuixi Road, Tianjin 300384, China, liujunjp@tjut.edu.cn 601 transitional characteristics under the transient response, and a weak fault characteristic is easily hidden in other transition components. Based on the above problems, the HHMS method is proposed to realize the extraction of the crack's fault on a cracked rotor. The simulation results show that the proposed method can effectively extract the crack's fault characteristics, and it can overcome above the faults of other methods. Based on the novel model of the crack's opening and closing, which can flexibly describe the crack's opening and closing degree of the rotor at any rotational position, the strategy of delaying crack propagation is discussed by analysing the crack's opening and closing degree with different rotors' parameters. Experimental results have also verified the effectiveness of the HHMS method. 1 MODEL DESCRIPTION 1.1 Rotor Model The Jeffcott rotor model with a lateral crack is shown in Fig. 1, and the disc is mounted on the centre of the elastic symmetrical shaft with a certain eccentric mass, and the crack is close to the disc. The horizontal rotor is supported at one end by a double-row of self-aligning ball bearing and at the other end by a single-row of deep-groove ball bearing. Because the single-row of deep-groove ball-bearing radial tilts have some rotational constraints, the elastic restoring force to be applied to the rotor will show nonlinearities, and the vibrations of the rotor will also appear to have nonlinear vibration characteristics. The differential equations of motion of the nonlinear rotor system under the transient response are shown below. mx + cx + kx+Nx = mem2 cos(m0t+—Xt2 + $0) 1 2 +meXsin(m0t+—Xt + $0)+mg my + cy + ky + Ny = mem2 sin(w0t+—Xt2 + $0) 1 +meXcos(w0t+—Xt +$0) where m is the disc mass, c is the damping coefficient, k is the stiffness coefficient, m is the rotational speed, X is the rotational acceleration, and e is the static unbalance of the disk. Nx and Ny are the nonlinear elastic restoring forces [16], corresponding to the x-and y-directions of the nonlinear term, respectively. 1.2 Stiffness Model of Cracked Rotor The relationship among the parameters is shown in Fig. 2. Due to the presence of the crack breath behaviour in a cracked rotor system, the rotor stiffness is changed during per revolution cycle. The stiffness model obtained by the neutral axis theory can accurately describe the crack's opening and closing process [17] and [18], and the mathematical equation [19] of the stiffness of the crack direction is shown as follows. Fig. 2. Schematic diagram of coordinate and parameter of a cracked rotor Fig. 1. Jeffcott rotor model with a lateral crack kz {a,p) / k = • 1 a sin(4a) n 4n 3 a + sin(4a) + p sin(4n) + cotp 4 2n 8n 2n 8n n 1, 7 a sin(4a) p sin(4p) cotp n n ---+ a < p <--a 2 2 /4 • 4 x n n (cos a - sin p), y -a continues to increase, the crack changes from the semi-open and the semi-closed state to the fully closed state. 2 HILBERT-HUANG MARGINAL SPECTRUM Under the transient response of a rotor system, the vibration signals of a nonlinear cracked rotor show the non-periodic change and nonlinearity. Because the weak fault characteristic is easily hidden in other transitional components, the fault characteristic of a crack is difficult to extract from the vibration signals. The HHMS method is proposed to process the complex non-periodic and transient vibration signals. The frequency components in the signals of the transient response can be decomposed, and the fault characteristic can be extracted from the complex vibration signals by using the HHMS method. 2.1 Empirical Mode Decomposition The complex non-stationary signal x(t) can be represented a sum form of the intrinsic modal function (IMF) ci and a residual component rn by using the empirical mode decomposition (EMD) [7] and [21]. 2.2 Hilbert Spectrum After using the EMD, the Hilbert transform (HT) is used to compute each IMFs The computed equation is as follows. *)=- P r ^ (7) n J—w t — t Through the HT processing, the analytical signal can be reconstructed. The expression is shown as follows. z (t) = c(t) + jy(t) = a(t ), (8) where a(t) = (c2 + y2)1/2 and d(t) = tan-1(y / c). Where the instantaneous frequency and amplitude of each IMFs can be obtained, so the original signal can be expressed as follows. x(t) = Re ±at (t )e' ^ ")d'. (9) i=1 Because mi(t) is the instantaneous frequency, the residual function can be omitted. The Hilbert spectrum can be computed to use the following equation. H t) = Re ¿a (t )e' " (')d'. (10) i=1 Some information can be seen from the above equation. The a,(t) and œ,(t) are functions of time variables and can form the three-dimensional spectrum picture with the time, the frequency and the amplitude. 2.3 Marginal Spectrum The marginal spectrum [22] can be obtained by integrating the Hilbert spectrum along the time variable, and the computational expression is shown below. h(a) = fQH (a, t )dt. (11) The marginal spectrum is the sum of the amplitude (energy) of a certain frequency in the signal, and it can also be understood that the sum of the amplitude of the frequency can be computed when the frequency appears at the whole processing. The frequency does not necessarily occur at the whole-time domain and does not necessarily only occur at a certain time. It may appear at several times with the different or same amplitude. The marginal spectrum realizes the distribution of the total amplitude for all frequencies in the signal, and it can accurately reflect the actual frequency components in the signal. 3 NUMERICAL SIMULATION The effectiveness of the proposed HHMS method for the nonlinear cracked rotor system is investigated with numerical simulation. Using Eq. (4), the transient response of the cracked rotor is obtained based on the Runge-kutta method. The response curves can be calculated by using Eq. (7) to Eq. (11). The computational parameters are as follows. Obtaining suitable values c=0.1 and e = 0.1, the angle between crack direction and the eccentric ft=0, the symmetric nonlinear coefficient ffV = 0.1 and the asymmetric nonlinear coefficient ec(1)=e/1-1 = 0.1. The resonance curves and the crack's opening and closing degree are discussed with the influence of different rotor parameters. 3.1 Crack Detection Based on Transient Response The transitional vibration signals of a rotor with a crack about five per cent of the depth [16] and [23] of the rotor's diameter are simulated and analysed using the wavelet method and the HHMS method, respectively. The simulation results of the two methods are compared with each other. There is only the major a) b) d) Fig. 4. Simulation results; a) and b) using wavelet method; c) and d) using HHMS method c) harmonic vibration component in the spectrum of the normal rotor. However, there is not only the major harmonic vibration component, but also the super harmonic vibration component in the spectrum of the cracked rotor. The vibration component of the super harmonic can be used to determine the existence of a crack in a rotor system. The simulation results are shown in Figs. 4 and 5 based on the different rotational accelerations and crack angles. The abscissa is the vibration frequency component and the ordinate is the amplitude R. 3.1.1 Wavelet Method Figs. 4a and b are the spectrum pictures of a normal rotor and a crack angle a = n/6 under the rotational acceleration A=0.01 to be calculated by the wavelet method. When the rotor does not have a crack, the calculation result is shown in Fig. 4a. There is the major harmonic vibration component, the sub-harmonic of order 1/2 vibration component, and the super-harmonic vibration component. When the rotor appears a crack and the crack angle a = n/6, the calculation result is shown in Fig. 4b. The vibration components in the figure are similar to the result shown in Fig. 4a, and the corresponding amplitude does not change greatly. At each harmonic resonance, some clutter components exist, which are easily ignored. A misjudgment, whether the rotor has a crack or not, is easy to be made from above results by computing with the wavelet method. 3.1.2 HHMS Method The analysis results with the HHMS method are shown in Figs. 4c and d. When the rotor does not have a crack, only the major harmonic vibration component exists, shown in Fig. 4c. When the rotor has a crack, and the crack angle a = n/6, the corresponding calculation results shown in Fig. 4d under the rotational acceleration A=0.01. There is the sub-harmonic of order 1/2 vibration component, the major harmonic vibration components and the super-harmonic vibration component to be caused by the existence of a crack. It can be seen that the results of the HHMS can be used to clearly determine whether the crack of the rotor system exists or not. b) d) Fig. 5. Simulation results; a) and b) using wavelet method; c) and d) using HHMS method 3.1.3 Crack Detection Based on Acceleration Effect Increasing the rotational acceleration A to 0.1 and keeping other parameters unchanged, the calculation results are shown in Fig. 5 by using the wavelet method and the HHMS method, respectively. When the rotor does not have a crack, all the vibrational components are reduced by the wavelet method from comparing Fig. 5a with Fig. 4a. Comparing Fig. 5c with Fig. 4c to be computed by the HHMS method, the major harmonic vibration component is also reduced. However, the sub-harmonic of order 1/2 vibration component exists owing to the nonlinearity. When the rotor has a crack and the crack angle a = n/6, the major harmonic vibration component is increased to compare Fig. 5b with Fig. 4b and the other vibration components are reduced. Comparing Fig. 5a with Fig. 5b, the vibration components are similar, and the amplitudes of components show little change. Therefore, changing the acceleration offers little help in detecting the rotor's crack with the wavelet method. Comparing the results shown in Fig. 5d with Fig. 4d, the major harmonic vibration component reduced due to the large rotational acceleration. However, the sub-harmonic of order 1/2 vibration component is almost reduced to zero, and the super harmonic vibration component increases remarkably. Therefore, increasing the rotational acceleration A is beneficial to detect a crack with the HHMS method. 3.2 Delay Crack Propagation Based on Rotor Parameters Effect Keeping the other rotor parameters unchanged, the crack's opening and closing degree is investigated by using the three types of rotor parameters. When the appropriate value of rotor parameters can reduce the crack's opening and closing degree, the delay crack propagation could be achieved. It is very important that we can adjust some parameter values to reduce the crack's opening and closing degree in an actual rotor system for the delay of crack propagation. The simulation results are shown in Fig. 6. The abscissa is time t, and the ordinate is the crack's opening and closing degree p. 3.2.1 Delay Crack Propagation Based on Acceleration Effect The calculation results of the crack's opening and closing degree with different rotational accelerations X are equal to 0.001, 0.01, and 0.1, and are shown in Figs. 6a, b, and c. With the increase of the rotational acceleration X, the crack's opening and closing degree decreases significantly. Its maximum value decreases from 0.027 shown in Fig. 6a to 0.016 shown in Fig. 6c, and the crack's opening and closing time is reduced from 1600 s to 100 s. The crack's opening and closing number reduces remarkably, shown in Figs 6a, b, and c. Therefore, when the rotor already has a crack, the acceleration value could appropriately increase under the transient running to reduce the crack's opening and closing degree to realize the delay of crack propagation. 3.2.2 Delay Crack Propagation Based on Angle p Effect In [20], the influence of the angle ft was discussed with the rotor vibration. This research has obtained a similar conclusion through simulation results. The crack's opening and closing degree is significantly reduced with the increase of the angle ft. When the angle ft changes to a small angle, the eccentric of the shaft is located in the same direction as the angle of the crack. The vibration amplitude of the rotor appears a large amplitude, so the crack's opening and closing degree is significantly increased. Therefore, when the rotor already has a crack, the delay of crack propagation could be achieved by adding a small mass on the disc to adjust the value of the angle ft, shown in Fig. 2. 3.2.3 Delay Crack Propagation Based on Nonlinear Parameters Effect The calculation results of the crack's opening and closing degrees with the value of the symmetric nonlinear parameter ft(0), is equal to 0.05, 0.15, and 0.25, shown in Fig. 7. When the symmetric nonlinear parameter ft(0) increases from 0.05 to 0.15, the corresponding crack's opening and closing degrees at the same time are reduced to some extent. When it increases from 0.15 to 0.25, the corresponding crack's opening and closing degree reduces at the forepart and increases at the rear part (t = 160 s later). The crack's opening and closing number significantly increases a) bj Time t [s] cj Fig. 6. Crack's opening and closing degree with acceleration changing; a) A = 0.001; b) A = 0.01; and c) A = 0.1 a) Time t [s] b) Time 1 [s] cj Fig. 7. Crack's opening and closing degree with nonlinear parameter p<°) changing; a) p<°) = 0.05; b) p<°) = 0.15; and c) p<°) = 0.25 b) Time t [s] c) Fig. 8. Crack's opening and closing degree with nonlinear parameters ec<1) and es(1) changing; a) ecm = esm = 0.1; b) ep = esm = 0.2; and c) ecm = esm = 0.3 a) from 140 s to 160 s; that is, a small is benefit to the delay crack propagation of a rotor system. Therefore, when a rotor already has a crack, the rotor support could be adjusted to decrease the value of the symmetric nonlinear parameter ft0) to realize the delay crack propagation. The calculation results of the crack's opening and closing degree with the value of the asymmetric nonlinear parameters ec(1) and e/1), are equal to 0.1, 0.2, and 0.3, and are shown in Fig. 8a, b, and c. The characteristics, seen from the figures, are that the crack's opening and closing degree corresponding to the same time is reduced and the crack's opening and closing number is also reduced with the increase of the asymmetric nonlinear parameters ec(1) and es(1). Therefore, it can be concluded that we can adjust the rotor support to increase the value of the asymmetric nonlinear parameters ec(1) and es(1) to realize the delay of crack propagation. 4 EXPERIMENT 4.1 Experimental Equipment The experimental equipment is shown in Fig. 9. The disc is located on the centre of the elastic rotor, the left end of the rotor is supported by a self-aligning double-row ball bearing, and the right end is supported by a single-row deep-groove ball bearing. The length of the rotor is 505 mm; the diameter of the rotor is 10 mm. The disc diameter and thickness are 100 mm and 20 mm, respectively. The rotor with a crack is located to the left end of 290 mm, the crack length is 20 mm, and the depth is 1 mm. The crack is machined by the EDM machine, and the cut part was re-fixed into the rotor to form the crack. In the left side of the rotor sets the pulse encoder to detect the rotor rotational speed, its displacements in x- and y- directions are measured with the position sensor that were set at the right side of the disc. The experimental data is collected by computer and the Donghua software, the initial position is set at the rotational speed m = 300 rpm. Fig. 9. Experimental setup 4.2 Experimental Result 4.2.1 Without Crack on the Rotor First, the experiment was carried out by using a normal rotor. The rotational acceleration was set to 500 rpm/min and the experimental range was set from the rotational speed of 300 rpm to the rotational speed of 2000 rpm. The experimental resonance curves were tested by the above setup and processed by using the wavelet and the HHMS method. The experimental results are shown in Fig. 10. The vibration component of the rotor system is mainly 150 Hz. The vibration characteristics were lost to use the wavelet method shown in Fig. 10a. The experimental results show that the HHMS method is the better extraction than the wavelet method on the transient response. 4.2.2 Cracked Rotor The experiment was carried by using the above-mentioned cracked rotor with a rotational acceleration X = 500 rpm/min and the experimental range set from a) b) Fig. 10. Experimental result without a crack; a) using wavelet method; and b) using HHMS method 0 100 200 300 400 500 0 100 200 300 400 500 aj Frequency/[Hz] y Frequency/[Hz] Fig. 11. Experimental result with a crack; a) using wavelet method; and b) using HHMS method the rotational speed of 300 rpm to the rotational speed of 2000 rpm. The experimental resonance curves were also processed by using the wavelet and the HHMS method. The experimental results are shown in Fig. 11. The vibration components of the rotor system mainly have 150 Hz, 75 Hz and 100 Hz, respectively. In particular, the super harmonic vibration component is clearly revealed at 75 Hz, because the crack exists in the rotor system. In addition, the sub-harmonic of order 1/2 vibration component also appeared at the rotational position of 280 Hz due to the nonlinearity shown in Fig. 11b. Experimental results showed that the HHMS method is very effective for the fault detection of the transient rotor system. 5 CONCLUSIONS In this research, the HHMS method is proposed for detecting a crack within the nonlinear rotor under the transient response. The simulation results show that it can effectively judge the existence of the early crack in a rotor system, and the crack which has about five percent of the depth of the rotor's diameter can be clearly judged by using HHMS. The experimental results also verify the effectiveness of the proposed HHMS method for the crack fault extraction from the transient response. This paper also proposes a novel method of describing the crack's opening and closing degree and discusses the delaying crack propagation with the different rotor parameters. The research results are proposed as follows. 1. The experimental and simulation results verify the effectiveness of the proposed HHMS method for the crack's fault extraction from the transient response, and the HHMS can accurately reflect the actual frequency components in the transient signal because of its integral effect. 2. When the aircraft takes off or landings, performs stunts, or makes a climbing-diving flight, the delay crack propagation can be achieved by increasing the rotational acceleration and deceleration. 3. When the angle ft and the support of the rotor are adjusted effectively, the delay crack propagation can be achieved, such as reducing the symmetric nonlinear term ft(0) and increasing the value of the asymmetric nonlinear term ec(1) and es(1). 6 ACKNOWLEDGEMENTS This work is supported by the Tianjin Natural Science Foundation of China (17JCZDJC38500 and 17JCYBJC18800) and the National Key Research and Development Program of China (2017YFB1303304). 7 NOMENCLATURES a depth of a crack, [m] c damping coefficient, [Nsm-1] C intrinsic modal function, [m] e static unbalance of the disk, [m] h thickness of the disk, [m] k stiffness coefficients, [Nm-1] l length of a rotor, [m] m mass of disc, [kg] Nx , Ny nonlinear elastic restoring forces, [N] rn residual component, [m] x, y displacements in x- and y-axis directions, [m] x(t) complex non-stationary signal, [m] a one-half angle of the crack [rad] ft the angle between crack direction and the eccentric [rad] ft0) symmetric nonlinear coefficient, [-] ô dimensionless factor, [-] ec(1), e/1-1 asymmetric nonlinear coefficient, [-] 0 inclination angle, [rad] k rotational acceleration, [rads-2] f initial phase angle, [rad] m rotational speed, [rads-1] 8 REFERENCES [1] Koc, A., Bartan, B., Gundogdu, E., Cukur, T., Ozaktas, H.M. (2017). Sparse representation of two-and three-dimensional images with fractional Fourier, Hartley, linear canonical, and Haar wavelet transforms. Expert Systems with Applications, vol. 77, p. 247-255, D0l:10.1016/j.eswa.2017.01.046. [2] Yao, J., Tang, B., Zhao, J. (2016). Improved discrete Fourier transform algorithm for harmonic analysis of rotor system. Measurement, vol. 83, p. 57-71, D0I:10.1016/j. measurement.2016.01.028. [3] Chandra, N.H., Sekhar, A.S. (2016). Fault detection in rotor bearing systems using time frequency techniques. Mechanical Systems and Signal Processing, vol. 72-73, p. 105-133, D0I:10.1016/j.ymssp.2015.11.013. [4] Bashar, S.K., Bhuiyan, M.I.H. (2016). Classification of motor imagery movements using multivariate empirical mode decomposition and short time Fourier transform based hybrid method. Engineering Science and Technology, an International Journal, vol. 19, no. 3, p. 1457-1464, D0I:10.1016/j. jestch.2016.04.009. [5] Dash, S.K., Rao, G.S. (2016). Arrhythmia Detection Using Wigner-Ville Distribution Based Neural Network. Procedia Computer Science, vol. 85, p. 806-811, D0I:10.1016/j. procs.2016.05.269. [6] Cai, C.L., Xiong, H.L., Ling, Y. (2008). Blind processing method for frequency hopping signals based on Winger-Ville rearrangement and stationary wavelet. Journal of Dalian Maritime University, vol. 34, no. 4, p. 115-118. [7] Wang, Y., Markert, R., Xiang, J., Zheng, W. (2015). Research on variational mode decomposition and its application in detecting rub-impact fault of the rotor system. Mechanical Systems and Signal Processing, vol. 60-61, p. 243-251, D0I:10.1016/j.ymssp.2015.02.020. [8] Ren, Z., Zhou, S., E, C., Gong, M., Li, B., Wen, B. (2015). Crack fault diagnosis of rotor systems using wavelet transforms. Computers & Electrical Engineering, vol. 45, p. 33-41, D0I:10.1016/j.compeleceng.2015.04.010. [9] Wang, D., Tsui, K.-L. (2017). Dynamic Bayesian wavelet transform: New methodology for extraction of repetitive transients. Mechanical Systems and Signal Processing, vol. 88, p. 137-144, D0I:10.1016/j.ymssp.2016.11.003. [10] Huang, Q., Jiang, D., Hong, L. (2009). Application of Hilbert-Huang transform method on fault diagnosis for wind turbine rotor. Key Engineering Materials, vol. 413-414, p. 159-166, D0I:10.4028/www.scientific.net/KEM.413-414.159. [11] Yan, J., Lu, L. (2014). Improved Hilbert-Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis. Signal Processing, vol. 98, no. 4, p. 74-87, D0I:10.1016/j. sigpro.2013.11.012. [12] Li, H., Zhang, Y., Zheng, H. (2009). Hilbert-Huang transform and marginal spectrum for detection and diagnosis of localized defects in roller bearings. Journal of Mechanical Science and Technology, vol. 23, no. 2, p. 291-301, D0I:10.1007/s12206-008-1110-5. [13] Wang, C., Lu, J. (2010). Fault diagnosis of diesel engine based on HHT marginal spectrum. Journal of Vibration Measurement & Diagnosis, vol. 34, p. 1265-1272. [14] Fu, K., Qu, J., Chai, Y., Zou, T. (2015). Hilbert marginal spectrum analysis for automatic seizure detection in EEG signals. Biomedical Signal Processing & Control, vol. 18, p. 179-185, D0I:10.1016/j.bspc.2015.01.002. [15] Hou, L., Chen, Y., Lu, Z., Li, Z. (2015). Bifurcation analysis for 2:1 and 3:1 super-harmonic resonances of an aircraft cracked rotor system due to maneuver load. Nonlinear Dynamics, vol. 81, no. 1-2, p. 531-547, D0I:10.1007/s11071-015-2009-1. [16] Ishida, Y., Yamamoto, T. (2013). Linear and Nonlinear Rotordynamics: A Modern Treatment with Applications. Wiley-VCH. [17] Lu, Z., Hou, L., Chen, Y., Sun, C. (2016). Nonlinear response analysis for a dual-rotor system with a breathing transverse crack in the hollow shaft. Nonlinear Dynamics, vol. 83, no. 1-2, p. 169-185, D0I:10.1007/s11071-015-2317-5. [18] Cavalini, A.A.Jr., Sanches, L., Bachschmid, N., Steffen, V.Jr. (2016). Crack identification for rotating machines based on a nonlinear approach. Mechanical Systems and Signal Processing, vol. 79, p. 72-85, D0I:10.1016/j. ymssp.2016.02.041. [19] Wang, Z., Lin, W., Wen, B. (2010). Analysis on the stiffness of the rotor system with a switching crack. Journal of Vibration and Shock, vol. 29, no. 9, p. 69-72. [20] Liu, J., Fan, Y., Wang, L.F., Chen, J., Wang, X., Liu, Z. (2016). Research of the delaying crack propagation on an aero-engine rotor. International Journal of Mechatronics & Automation, vol. 5, no. 4, p. 190-200, D0I:10.1504/IJMA.2016.084217. [21] Lei, Y., Lin, J., He, Z., Zuo, M.J. (2013). A review on empirical mode decomposition in fault diagnosis of rotating machinery. Mechanical Systems and Signal Processing, vol. 35, no. 1-2, p. 108-126, D0I:10.1016/j.ymssp.2012.09.015. [22] Huang, N.E., Shen, Z., Long, S.R., Wu, M.C., Shih, H.H., Zheng, Q., Yen, N.C., Chi, C.T., Liu, H.H. (1998). The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings Mathematical Physical & Engineering Sciences, vol. 454, no. 1971, p. 903-995, D0I:10.1098/rspa.1998.0193. [23] Ishida, Y. (2008). Cracked rotors: Industrial machine case histories and nonlinear effects shown by simple Jeffcott rotor. Mechanical Systems & Signal Processing, vol. 22, no. 4, p. 805-817, D0I:10.1016/j.ymssp.2007.11.005. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)10, 611-620 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2018.5349 Original Scientific Paper Received for review: 2018-03-12 Received revised form: 2018-06-21 Accepted for publication: 2018-07-02 Influence of Cellular Lattice Body Structure on Gear Vibration Induced by Meshing Riad Ramadani1* - Marko Kegl2 - Jožef Predan2 - Aleš Belšak2 - Stanislav Pehan2 1 University of Prishtina, Faculty of Mechanical Engineering, Kosovo 2 University of Maribor, Faculty of Mechanical Engineering, Slovenia This paper discusses the influence of a gear body structure on gear vibrations induced by meshing. For this purpose the spur gear body was designed as cellular lattice structure. In order to reduce the stress levels as much as possible and to remove stress concentrations, the lattice structure was optimized by engaging a topology optimizer. The obtained lattice structure was expected to have a positive influence on vibrations reduction due to longer pressure waves travelling paths and several path direction changes. To verify this experimentally, the spur gear was produced from titanium alloy Ti-6Al-4V ELI by using selective laser melting technique. Furthermore, a new precise closed loop test rig was designed and produced to measure experimentally vibrations caused by rotating and lubricated gear pairs. Vibrations input data were obtained by measuring accelerations on the housing of the test rig. The signals were analyzed in frequency and time-frequency domains. Experimental results confirm that the cellular lattice structure of the gear body, especially if the voids are filled with a polymer, has a positive effect on reduction of vibrations induced by meshing of engaged gears. Keywords: gear vibration, cellular lattice structure, topology optimization, test rig, signal analysis Highlights • The body of a spur gear was designed as a cellular lattice structure. • A spur gear with optimized lattice body structure was produced using selective laser melting technique. • A closed-loop test rig for testing gear pairs was developed and produced. • The time signals of acceleration were analyzed in frequency and time-frequency domain. • Experimental measurement of acceleration has verified that cellular lattice structure and polymer inserted into it, can effectively reduce gear vibration. 0 INTRODUCTION Gears are primary components in many transmission applications. Therefore, the reduction of noise and vibration, generated by engaged gears, is of great importance. Gear vibrations are induced by external and internal influences. Typical external sources are time variable driving torque and speed. On the other hand, the internal sources originate from tooth meshing [1] and [2]. This process is plagued by stiffness variation (due to variable number of teeth engaged) and interacting teeth impacts [3] and [4] due to deformation, wear, and imperfections. In this way each toothed gear ring becomes an internal generator of vibrations. Gear vibrations are transmitted to the housing by shafts and bearings, which creates annoying noise and decreases the quality of transmission operation. Consequently, any measure that may potentially reduce vibrations is worth of being investigated. Earlier research work on gear vibration reduction was focused mainly on tooth profile modification [5] and [6]. More recently, this work was enriched by introducing sophisticated optimization procedures with the aim to improve tooth flanks [7] to [9]. Progress and improvements were achieved, but regardless of any profile or flank modification, tooth meshing still remained a generator of vibrations. An interesting and successful vibration reduction was achieved by using a gear pair with phasing [10]. Such gears run smoothly but they are comparatively heavier as plain gears due to the complicated inner structure. Another approach, taken by some authors [11] and [12], was an attempt trying to reduce gear vibrations by improving the function of the gear body. For example, a convenient powder material was used to fill the bores in the gear body. To some extent this was successful because this hinders the pressure wave propagation and induces some damping. Till now the aim of majority of modifications related to the gear body was to reduce the weight of the gear. Weight reduction resulted in thinner walls and introduction of holes. These measures reduced the weight of the gear but did not have any significant impact on vibration generation and emission. From the dynamic point of view, a weaker and slimmer gear body might even make the situation worse [13]. So far, all investigations related to gear vibration reduction were based on the assumption that the *Corr. Author's Address: University of Prishtina, Faculty of Mechanical Engineering, Bregu i diellit, 10000 Prishtina, Kosovo, riad.ramadani@uni-pr.edu 611 gear has to be produced by known conventional technologies. However, in recent years new additive manufacturing (AM) technologies came into foreground and the impressive development of these technologies opened new potentials and approaches to structural design. This was accompanied by a respectful development of numerical optimization methods, focused on the potential benefits offered by the new AM technologies. This opened up a completely new perspective of how to address the design of structural parts. The aim of this paper is to engage the newly emerged production technologies and optimization software in order to address gear body design. The objective is to design a light-weight gear body that will exhibit low generation of vibrations and will of course fulfil adequate load-carrying, stiffness, and strength requirements. Additionally, the gear body should contain sufficient empty space that could potentially be filled by a suitable polymer with the aim to achieve additional vibration reduction. In this work the gear body is designed as a cellular lattice structure. Such a structure prolongs the pressure wave traveling paths and the traveling direction is often changed. Longer paths dissipate more energy due to internal friction and direction changes also influence positively energy dissipation. It is therefore expected, that such a design could exhibit good vibration reduction properties. A loaded lattice structure can easily exhibit stress concentrations, especially at lattice joints. This should be avoided by all means because stress concentrations could easily induce crack initiation, which could lead to structural failure. To remove stress concentrations and to lower the stress levels to a minimum, a topology optimizer CAESS ProTOp was engaged. This high-performance topology optimizer enables to configure numerically a solid part as a lattice structure. The configured design is immediately ready for optimization, where material is redistributed so that stress concentrations are removed and stress levels reduced as much as possible. An optimized lattice structure can typically be produced only by AM technologies. Therefore, in this work the obtained optimal spur gear was produced by selective laser melting (SLM) technique. This gear was then tested on a closed loop test rig. Although, various rigs have been used to test the dynamic behavior of gears so far [14] and [15], in this work a special test rig was developed in order to enable precise and simultaneous measurements of displacement, sound pressure, and acceleration during gear operation. The acceleration was measured on the housing and the signals taken from the accelerometer were analyzed in frequency and time-frequency domains; choosing these domains is typically the most useful approach in gear vibration analysis [16] to [18]. The recorded signals were processed within LabView. This paper is organized as follows. In Section 1 gear body design and optimization is described briefly. Section 2 describes the developed test rig and measurements setup. Sections 3 is devoted to frequency and time-frequency analysis. 1 DESIGN AND OPTIMIZATION OF GEAR BODY STRUCTURE In this work a spur gear corresponding to the data shown in Table 1, was addressed. The diameter of the fixing hole was 25 mm and the adopted nominal torque was 200 Nm. The numerical finite element analysis (FEA)-ready model of the gear with a full solid body was prepared in SIMULIA Abaqus [19]. A linear elastic and isotropic material was selected. Table 1. Basic data of a spur gear Data Value Number of teeth 34 Module 2.5 mm Pressure angle 20° Width 10 mm Pitch diameter 85 mm Base diameter 79.874 mm Root diameter 78.558 mm Tip diameter 90 mm Module of elasticity 206,000 MPa Poisson's ratio 0.3 The gear was meshed with linear hexahedral elements of type C3D8, which are eight-node isoparametric elements with 24 displacement degrees of freedom. Hexahedral elements were used because the generated mesh exhibited symmetries which influence positively the symmetries of the final design. Due to the planned topology optimization, the generated mesh needed to be finer than typically engaged for usual FEA. Consequently, the number of elements was relatively high; the mesh contained a total of 2,918,412 elements. To influence positively the design symmetries and cyclic lattice pattern, a total of 68 static load cases were defined: one load case per one side of each tooth. The applied load was calculated from the nominal torque, which was converted into pressure of p = 618 MPa and applied to a tooth flank. The loading Fig. 1. Von Misses stress; a) of initial lattice structure, and b) of optimized lattice structure point was the outer single meshing point. The inner surface of the hub was fixed by adequate displacement constraints. This solid body FEA model was then imported into CAESS ProTOp [20] for further design configuration and optimization. According to [21] topology optimization of a lattice structure modeled by solid finite elements differs from conventional topology optimization in one crucial step: configuring of the optimization domain. This configuration step enforces the chosen lattice structure on a usual solid domain meshed by solid finite elements and prevents the optimizer to add material outside of the domain of the lattice structure. To reconfigure the solid gear body into a cellular lattice structure, two different lattice cell types have been chosen; a cube diagonal lattice cell and a plane diagonal lattice cell. The lattice structure has been defined in the cylindrical coordinate system. The inner part of the gear body was designed by using a cube diagonal lattice cell, while for the outer part of the gear body a plane diagonal lattice cell was used (Fig. 1). A structure configured in ProTOp is immediately ready for topology optimization. The process is iterative and each iteration requires one full finite element analysis per load case and one computation of topology parameters. In our case each iteration required 68 full finite element (FE) analyses. A detailed explanation for cell parameters used to design and optimize the gear body structure can be found in [22]. The stresses corresponding to the initial configured design with a volume part of 41.2 % exhibited peak levels about 590 MPa (Fig. 1; gear tooth region), while the stresses within the lattice structure reached maximum levels about 400 MPa (Fig. 1; yellow color). The structure was then optimized and the volume part of the free domain was increased from 41.2 % up to 49.3 %. Surprisingly this also reduced stresses within the tooth region, which now peaked to about 550 MPa. More importantly, however, are the stress levels within the optimized lattice structure. These stresses dropped dramatically with maximum stress levels reaching now only about 150 MPa. The optimized structure exhibited drastically reduced stress levels and thus represents a good design from the mechanical point of view. From the technological point of view, however, it can practically not be produced by any of the conventional technologies. Therefore, the optimized spur gear was produced by using a titanium alloy Ti-6Al-4V ELI powder and the AM technology called the SLM technique. SLM allows a part to be built additively, layer by layer, from an adequate provided powder metal [23]. Lattice structure _ ^ Grinded tooth flank Fig. 2. AM-produced optimized titanium gear The parameters that have been used to produce the gear were: laser power 75 W, scan speed 600 mm/s, and layer thickness 35 ^m. The teeth flank and the inner surface of the hub were grinded in the same tolerances as the gear with a solid body. The AM-produced titanium gear is shown in Fig. 2. 2 DEVELOPMENT OF GEAR TEST RIG A completely new closed-loop test rig for gear vibration measurements was developed and manufactured in order to meet all internal requirements and cover all planned activities. 2.1 Design and Operation The test rig consists of two heavy steel blocks, which are firmly connected with two connecting shafts; these components form a very rigid frame. Heavy-duty rotating components are firmly positioned into the steel blocks and the test gears are mounted from the outside. In this way, it is not necessary to disassemble the bearing system to replace the tested gear pair. The axial distance between the shafts is 90 mm, Fig 3. Torque Steel applying device block 'i iV Connecting shaft Electromotor Locking <•' block" Test gear Fig. 3. Layout of the gear test rig The closed-loop consist mainly of two shafts connected by two gear pairs. One gear pair just transmits the load torque from one shaft to the other; the other gear pair is the tested one. One shaft includes a clutch and the other shaft includes a torsional spring. The load torque is applied with a plain digital torque wrench while the clutch is open and one side of the clutch is firmly tightened by the locking block. At the desired torque level, the clutch is closed and the locking block and torque applying device removed from the test rig. The test rig is run by an electric motor and the RPM is controlled by a plain frequency regulator. To enable monitoring of lubrication conditions a special partly transparent and easy removable housing was made. Depending on the equipment this system enables simultaneous displacements, sound pressure, and acceleration measurements. 2.2 Experimental Setup To measure the vibrations, the National Instruments NI PXI 4472 system and the PCB accelerometer (model 3 56A3 2) were used. The accelerometer was placed on the steel block. The signals taken from the accelerometer were transmitted to the data acquisition system and then to the computer for signal data analysis. The test rig was placed on a firm table and a rubber layer was placed between the table and the steel blocks to isolated the test rig from any external vibrations. The experimental setup is shown in Fig. 4. Fig. 4. Experimental setup The vibration measurements were performed with 3 different gear pairs as follows: • S-S gear pair is a reference pair consisting of two usual solid body gears; • S-L gear pair consists of a usual solid body gear and of a gear with optimized lattice structure body; • S-LP gear pair is the same as S-L gear, except that the void regions between lattices are filled with a polymer. The S-L gear pair is shown in Fig. 5. Fig. 5. The S-L gear pair mounted on the test rig 3 ACCELERATION FREQUENCY AND AMPLITUDE ANALYSIS In order to estimate the dynamic behavior of various gear pairs, the measured acceleration signals were analyzed. For this purpose, the original accelerations given in the time domain were transformed to the frequency and time-frequency domains. 3.1 Frequency Domain Analysis Conventionally, this is achieved by using the Fourier transformation by which a time-dependent quantity x(t) can be represented in the frequency domain by its Fourier transform (FT) Xf) defined as [24]: X(t) = | x(t) • e~j2nftdt, (1) —» where f is the frequency to analyze and j is the imaginary number [24]. For a computer implementation, this function is replaced by the discrete Fourier transform (DFT), which can be obtained by numerical integration of Eq. (1) as [24]: N -1 X(fk) = j x(t,) • e-J2f (t+ -t,), k = 1,..., N, (2) i=0 where f is the frequency of the kth harmonic. For computational efficiency reasons, the DFT is typically computed by using the fast Fourier transformation (FFT). This procedure delivers amplitudes at discrete frequencies of the underlying harmonics. In the case of gear pairs the gear mesh frequency is known to be the fundamental frequency or the first harmonic. The dynamic behavior of engaged gears is characterized from the distribution of exposed frequencies (the ones with significant amplitudes). To estimate the dynamic behavior of tested gear pairs, time signals from an accelerometer placed on the housing were acquired. The applied torque was T = 60 Nm and the number of rotations was n = 310 rpm. This means that the fundamental frequency was fi = 176 Hz. The frequencies were analyzed up to 7 kHz, while the measured time signal was 0.2 s long. The time signals are converted to frequency and time-frequency domains by using LabView. The acceleration time signal of the reference S-S gear pair is shown in Fig. 6, whereas its frequency spectrum is shown in Fig. 7. From the frequency spectrum it is evident that the dominating amplitude appeared at the 5th harmonics f = 5f = 880 Hz. The other frequencies with notable amplitudes were those corresponding to the 1st, 3rd, and 8th harmonics. The frequency components after 8th harmonics are observable but relatively low. The acceleration time signal for the S-L gear pair with lattice structure is presented in Fig. 8, whereas its frequency spectrum is given in Fig. 9. 0.450.40.350.3.1 0.251 °"2" ^ 0.150.10.05- il UJ /j'VvkA ..... 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 Frequency [Hz] Fig. 7. Frequency spectrum of acceleration for the reference s-s gear pair One can see from the frequency spectrum that the dominating amplitude is now the one corresponding to the 6th harmonics f = 6/1 = 1056 Hz. An only somewhat lower amplitude is now observed for the 3rd harmonics. Still somewhat lower but nicely exposed are the 7th, 8th, and 9th harmonics. The acceleration time signal for the S-LP gear pair with lattice structure and filled with polymer is shown in Fig. 10; its frequency spectrum is shown in Fig. 11. One can see that the frequencies of the 6th and 9th harmonics are the dominating ones. A somewhat lower but nicely notable frequencies are those of the 4th, and 5th harmonics. As shown in Fig. 7. the largest amplitude was observed for the S-S gear with a solid body; a peak of approximately 0.42 g appeared at the frequency of 880 Hz. For the S-L gear pair with lattice structure (Fig. 9), the maximal amplitude was substantially lower; a peak of around 0.26 g appeared at the frequency of 1056 Hz. A further reduction of the maximal amplitude value was achieved by filling the void regions within the lattice structure with a polymer (Fig. 11). The S-LP gear pair exhibits a maximum of 0.165 g at the frequency of 1056 Hz. By comparing the frequency spectrums of all tested gear pairs (Fig. 7, 9, and 11), one can see that the dominant frequency amplitude is the highest and notable frequency components are relatively widely spread for the solid body gear. A lattice structure of the gear body obviously lowers the maximum amplitude and the amplitudes of all notable harmonics are more equilibrated. This trend is even more observable for the lattice-polymer gear, whose spectrum also indicates an obvious reduction of vibration energy. This example indicates that a lattice body structure can potentially do a good job in reducing gear vibrations. It looks also that the inserted polymer is quite efficient in increasing the structural damping, which also helps a lot. 3.2 "Time-Frequency Domain Analysis In the spectrum of the frequency domain, it is possible to determine, which frequencies were present in the signal but it is not possible to know the time when those frequencies appeared. The time-frequency analysis introduces the time variable and makes it possible to determine how those frequencies of non-stationary signals change with time, and what the levels of signal energy are [25]. Time-frequency analysis is a three-dimensional time, frequency, and amplitude representation of a signal. The short time Fourier transformation (STFT) is a linear time-frequency transformation and it is used for analyzing non-stationary signals. The basic idea of STFT is to divide a signal into short time segments, and then _jdj U A -WAJL J, . „JU-t ..,,,„, X......... 1 . 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 Frequency [Hz] Fig. 9. Frequency spectrum of acceleration for the S-L gear pair 0.450.40.35- ■3 °3" J 0.25- < 0.150.10.05- li ! ^LL LIA * 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 Frequency [Hz] Fig. 11. Frequency spectrum of acceleration for the S-LP gear pair the Fourier analysis is carried out for each segment separately. It is assumed that signals are stationary within each segment. The STFT is obtained in the following way: for any signal x(t), let's suppose that w(t - t) is a window function centered at the time t, to produce a segmented signal x(t) w(t -t), where t is a time variable. By applying the FFT to each of such segmented signals, the STFT is obtained as [26]: STFTx (t, f) = | x(t) • w(t -1 )e— fT. (3) —<» The time window type and positioning must be selected properly in order to get a good compromise between good time and frequency resolution of the result. For reference S-S gear pair with a solid body, the time-frequency spectrogram is presented in Fig. 12. Although the result is smeared over time and frequency, one can roughly estimate that the dominating 5th harmonic of acceleration is present in about 62 % of the time interval (Fig. 12; red color). One can also see that the 5th harmonic is not stable, but exhibits an occasional vanishing and pulsating character. For the S-L gear pair with lattice structure the spectrogram is given in Fig. 13. One can see that the dominating 6th harmonic of acceleration is present in about 45 % of the time (Fig. 13; red color). The 3rd harmonic with slightly lower amplitude can also be observed; however, its time presence is significantly lower. One can say that the 6th harmonic is semi-stable with a pulsating character. For the S-LP gear pair with polymer filled lattice structure the spectrogram is shown in Fig. 14. One can see that the dominating 6th harmonic of acceleration is present in about 33 % of the time (Fig. 14; red color). The amplitudes are also observable lower than in case of other gear pairs. The 6th harmonic is semi-stable with a pulsating character. When the time-frequency spectrograms of all tested gear pairs are compared, it is obvious that the notable frequency components of the S-S gear pair are longer present in time. This presence decreases for the S-L and even more for the S-LP gear pairs. A similar trend can be observed for the amplitudes, with the S-LP gear pair exhibiting by far the lowest amplitudes. Regarding the signal stability one can say that the lattice structure caused a more stable response in the sense that the dominating frequencies are more evenly present during the observed time interval. At this point it might be worth noting that besides of reduced vibrations, the optimized lattice body gear also exhibits a substantially reduced mass. Namely, the volume of the gear with a solid body is about 51130 mm3, while the volume of the gear with the optimized lattice structure is about 32570 mm3. Thus, Fig. 14. Time-frequency spectrogram of acceleration for the S-LP gear pair the volume of the used material is reduced by around 0.40137 kg, whereas the weight of the gear with the one third. The weight of the gear with a solid body is optimized lattice structure is 0.25567 kg. 4 CONCLUSION A new approach to gear vibration reduction was presented, based on the gear body structure modification. The body structure was designed as a cellular lattice structure that was optimized for stress level and stress concentrations by employing topology optimization. A spur gear with optimized lattice structure was produced by using Selective Laser Melting technique. To measure the dynamic behavior of such gears, a totally new closed-loop test rig has been developed and produced. The test rig enables to simultaneously measure displacement, sound pressure and acceleration, while the gears rotate. According to experimental results of accelerations for the tested gears, it was proven that the lattice structure can effectively reduce gear vibration analyzed in frequency spectrum and time-frequency spectrogram. A lattice structure provides additionally a nice opportunity to be filled with an adequate polymer. According to the results, this may further reduce the vibrations. The main disadvantage of a lattice structure is that currently it is a rather expensive choice since traditional manufacturing is not possible and consequently the new AM technologies have to be engaged. Additionally, the teeth flank of a printed gear must be grinded which requires an extra production operation. Another, less obvious, point of concern is that currently printed materials are not very resistant against fatigue crack initiation and consequent structural failure. To mitigate this drawback, a printed part has to be optimized by all means in order to reduce the stresses as much as possible. 5 ACKNOWLEDGEMENT The article was written under favourable circumstances of Faculty of Mechanical Engineering, University of Maribor, Slovenia. The authors want to thank Erasmus Mundus JoinEU-SEE PENTA scholarship program, Slovenian Research Agency (research core funding No. P2-0063) as well as CAESS company for the research support. 6 REFERENCES [1] Mohammed, O.D., Rantatalo, M., Aidanpaa, J.O. (2013). Improving mesh stiffness calculation of cracked gears for the purpose of vibration-based fault analysis. Engineering Failure Analysis, vol. 34, p. 235-251, D0l:10.1016/j. engfailanal.2013.08.008. [2] Saxena, A., Chouksey, M., Parey, A. (2017). Effect of mesh stiffness of healthy and cracked gear tooth on modal and frequency response characteristics of geared rotor system. Mechanism and Machine Theory, vol. 107, p. 261-273, D0I:10.1016/j.mechmachtheory.2016.10.006. [3] Agemi, F.M., Ognjanovic, M. (2004). Gear vibration in supercritical mesh-frequency range. FME Transaction, vol. 32, p. 87-94. [4] Ognjanovic, M., Agemi, F. (2010). Gear vibrations in supercritical mesh-frequency range caused by teeth impacts. Stojniski vestnik - Journal of Mechanical Engineering, vol. 56, no. 10, p. 653-662. [5] Kahraman, A., Blankenship, G.W. (1999). Effect of involute tip relief on dynamic response of spur gear pairs. Journal of Mechanical Design, vol. 121, no. 2, p. 313-315, D0I:10.1115/1.2829460. [6] Liu, G., Parker, R.G. (2008). Dynamic modeling and analysis of tooth profile modification for multimesh gear vibration. Journal of Mechanical Design, vol. 130, no. 12, p. 121402-13, D0I:10.1115/1.2976803. [7] Faggioni, M., Samani, F.S., Bertacchi, G., Pellicano, F. (2011). Dynamic optimization of spur gears. Mechanism and Machine Theory, vol. 46, no. 4, p. 544-557, D0I:10.1016/j. mechmachtheory.2010.11.005. [8] Ma, H., Pang, X., Feng, R., Wen, B. (2016). Evaluation of optimum profile modification curves of profile shifted spur gears based on vibration responses. Mechanical Systems and Signal Processing, vol. 70-71, p. 1131-1149, D0I:10.1016/j. ymssp.2015.09.019. [9] Ghosh, S.S., Chakraborty, G. (2016). On optimal tooth profile modification for reduction of vibration and noise in spur gear pairs. Mechanism and Machine Theory, vol. 105, p. 145-163, D0I:10.1016/j.mechmachtheory.2016.06.008. [10] Cheon, G.-J. (2010). Numerical study on reducing the vibration of spur gear pairs with phasing. Journal of Sound and Vibration, vol. 329, no. 19, p. 3915-3927, D0I:10.1016/j. jsv.2010.04.005. [11] Xiao, W., Li, J., Wang, S., Fang, X. (2016). Study on vibration suppression based on particle damping in centrifugal field of gear transmission. Journal of Sound and Vibration, vol. 366, p. 62-80, D0I:10.1016/j.jsv.2015.12.014. [12] Xiao, W., Huang, Y., Jiang, H., Jin, L. (2016). Effect of powder material on vibration reduction of gear system in centrifugal field. Powder Technology, vol. 294, p. 146-158, D0I:10.1016/j. powtec.2016.01.038. [13] Li, S. (2008). Experimental investigation and FEM analysis of resonance frequency behavior of three-dimensional, thin-walled spur gears with a power-circulating test rig. Mechanism and Machine Theory, vol. 43, no. 8, p. 934-963, D0I:10.1016/j.mechmachtheory.2007.07.009. [14] Arun A.P., Senthil kumar, A.P., Giriraj, B., Faizur rahaman, A. (2014). Gear test rig - A review. International Journal of Mechanical & Mechatronics Engineering, vol. 14, no. 5, p. 1626. [15] Akerblom, M. (1999). Gear test rig for noise and vibration testing of cylindrical gears. Proceedings OST Symposium on Machine Design, p. 183-199. [16] Belsak, A., Flasker, J. (2006). Method for detecting fatigue crack in gears. Theoretical and Applied Fracture Mechanics, vol. 46, no. 2, p. 105-113, D0l:10.1016/j.tafmec.2006.07.002. [17] Belsak, A., Flasker, J. (2008). Vibration analysis to determine the condition of gear unit. Stojniski vestnik - Journal of Mechanical Engineering, vol. 54, no. 1, p. 11-24. [18] Mohammed, O.D., Rantatalo, M. (2016). Dynamic response and time-frequency analysis for gear tooth crack detection. Mechanical Systems and Signal Processing, vol. 66-67, p. 612-624, D0I:10.1016/j.ymssp.2015.05.015. [19] Abaqus/CAE 6.14 User's guide. (2014). Dassault Systèmes Corp., Providence. [20] CAESS ProTOp, from http://caess.eu/, accessed on 2018-0309. [21] Harl, B., Predan, J., Gubeljak, N., Kegl, M. (2017). On configuration-based optimal design of load-carrying lightweight parts. International Journal of Simulation Modelling, vol. 16, no. 2, p. 219-228, D0I:10.2507/ijsimm16(2)3.369. [22] Ramadani, R., Belsak, A., Kegl, M., Predan, J., Pehan, S. (2018). Topology optimization based design of lightweight and low vibration gear bodies. International Journal of Simulation Modelling, vol. 17, no. 1, p. 92-104, D0l:10.2507/ ijsimm17(1)419. [23] Olakanmi, E,O., Cochrane, R.F., Dalgarno, K.W. (2015). A review on selective laser sintering/melting (SLS/SLM) of aluminium alloy powders: Processing, microstructure, and properties. Progress in Materials Science, vol. 74, p. 401-477, D0I:10.1016/j.pmatsci.2015.03.002. [24] Brigham, E.O. (1988). The Fast Fourier Transform and Its Applications. Prentice Hall, Michigan. [25] Boashash, B. (2003). Time Frequency Analysis. Elsevier Science, Amsterdam. [26] Feng, Z., Liang, M., Chu, F. (2013). Recent advances in time-frequency analysis methods for machinery fault diagnosis: A review with application examples. Mechanical Systems and Signal Processing, vol. 38, no. 1, p. 165-205, D0I:10.1016/j. ymssp.2013.01.017. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)10, 621-631 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2018.5235 Original Scientific Paper Received for review: 2018-01-23 Received revised form: 2018-05-23 Accepted for publication: 2018-06-13 A Gage Study Applied in Shear Test to Identify Variation Causes from a Resistance Spot Welding Measurement System Fabricio Alves de Almeida1* - Guilherme Ferreira Gomes2 - Rachel Campos Sabioni3 -José Henrique de Freitas Gomes1 - Vinicius Reno de Paula1 - Anderson Paulo de Paiva1 - Sebastiâo Carlos da Costa1 federal University of Itajuba - Institute of Industrial Engineering and Management, Brazil 2Federal University of Itajuba - Mechanical Engineering Institute, Brazil 3Sorbonne Universités, Université de Technologie de Compiègne - Department of Mechanical Engineering, France Resistance welding processes, especially spot welding, have wide applicability in the industry especially in the automotive sector, due to its fast execution and the non-use of consumables. In addition, the search for quality improvement of the final product is incessant and, in a capable process, there should be no error related to the measurements. In this study, the NGR&R was used by the ANOVA method to identify the variation components of the measurement system in the shear test for two quality characteristics: tensile-shear strength and ultimate strain. The experiments were conducted on a hot dip galvanized steel by using design of experiments to select the parts in order to represent the real amplitude of the process. From the results it was possible to verify that one of the destructive test machines used in this study has a strong variability, evidencing that some adjustments and improvements are necessary in the coupling of the specimens for steels with coating layers (such as galvanized steel, which has a layer of zinc). Keywords: spot welding, measurement system analysis, shear test, NGR&R, ANOVA Highlights • A nested GR&R study applied in shear test for a resistance spot welding process. • Design of experiments to select the parts in order to represent the real amplitude of the process. • The analysis of variance method to identify variation causes for two tensile machines and two quality characteristics: tensile-shear strength and ultimate strain. • The results showed that Machine 1 presents greater contribution on the system variability, with measurement results outside the control, as well as a lower degree of repetitiveness than Machine 2. 0 INTRODUCTION Resistance spot welding (RSW) is a structural joint technique widely used in the automotive sector [1]. RSW is highlighted among welding processes due to its features that favor the industry such as agile operation, which is easily suitable for automatic processes, simple handling, diverse applications and low cost [2] to [4]. Because of its wide applicability and importance, new methodologies for parameter adjustment have been applied to RSW improvement, contributing to the process control and capability. Among the available methods for verification of the weld point, there is the shear test, which is characterized by the application of opposing forces causing stress in a sliding movement for a given sample. Since this type of test allows to evaluate the quality of welded point, it is being increasingly used, as described by Feng et al. [5], Zhang et al. [6], Martin et al. [7], Shan et al. [8], Chen et al. [9], Manladan et al. [10]. The search for quality improvements has been leading industries to improve their efficiency. However, devoting improvements only to the process may not contribute to make it better, as the variability can also be caused by the measurement system. Therefore, it is necessary to verify the measurement system variability in industrial processes, such as RSW. There are several methods for controlling and monitoring quality in the RSW process, such as: expulsion detection in materials [1] and [11]; strength estimation based on sonic emission [12]; welding current analysis on weld strength [13]; temperature measurement [14]; electrode displacement [15] and [16]; and other types of control (i.e. electrical variables, ultrasound transmission and acoustic emission) [17]. However, the control approaches must be verified through specific tests, which illustrate the mechanical characteristics necessary for their capability evaluation, such as the shear test. The shear test is characterized as a destructive test that evaluates the mechanical strength of the weld point in relation to shear stresses. Destructive tests are performed from time to time, by sampling, being widely employed in the automotive sector. *Corr. Author's Address: Federal University of Itajuba, 1303 BPS Avenue, Itajuba, Brazil, fabricio.alvesdealmeida@gmail.com 621 Thereunto, quantitative methods are used to verify the process quality, in which the measuring device must be validated before data collection [18]. The variability in destructive test results may arise from the measurement system itself as well as from the manufacturing process [19] to [21]. In this case, the measurement error must be avoided in experimental procedures. On the other hand, in quality methodologies such as Six Sigma, before analyzing the process, it is necessary to verify the capability of the measurement system (MS). One of the techniques used to evaluate the variation components of MS, according to Peruchi et al. [22], is gage repeatability and reproducibility (GR&R), in which the MS variability is quantified and analyzed through an analysis of variance (ANOVA). The repeatability (Fig. 1a) is characterized by the variation within the system under fixed and already defined conditions of measurements (part, environment, operator, instrument among others), i.e., the variation acquired in a measuring equipment used several times by an operator, which is based on a single part [23] to [27]. In this way, the reproducibility (Fig. 1b) is characterized by the average variation between evaluated systems, being the variation found in the mean of different operators using the same equipment to perform the measurement of a single part [23] and [28] to [31]. ANOVA is a statistical method applied together with GR&R, as described by: Shi et al. [32]; Deshpande et al. [33]; Zhua et al. [34]; García & del Río [35]; Woodall & Borror [20]; Johnson et al. [36]. However, there are few applications focused on the RSW process, reinforcing the potential contribution of this work. There are several studies in scientific literature related to the measurement system analysis on RSW processes, such as: Wan et al. [37]; Degidi et al. [38]; Wang et al. [39]; Xia et al. [40]; Simoncic, & Podrzaj [15]; Al-Jader et al. [41]; Lei et al. [42]; Lai et al. [43]; Podrzaj et al. [11]. Conduct the measurement system analysis by comparing tensile machines for destructive tests (such as the shear test) on RSW process is original, since published studies that have performed this type of analysis in this process do not present a similar evaluation. In addition, considering the accessed references in this current work, it is possible to confirm the relevance of applying the galvanized steel in this process. Given the importance that a measurement system presents in industrial processes (especially in RSW), this paper aims to perform a gage repeatability and reproducibility for nested design (NGR&R) applied to the shear test to verify variation causes in a RSW process measurement system, using univariate method ANOVA, comparing two tensile testing machines. The limits of the process parameters were defined from preliminary tests, in order to guarantee desired failure modes. Design of experiments (DOE) methodology was used, being a statistical strategy to model experiments [44], thus that the characteristics of the parts represented a real amplitude of the RSW process. 260 340 420 500 580 660 740 X Fig. 1. a) Repeatability and b) reproducibility on a measurement system This article is organized as follows: a bibliographic review on the RSW process and GR&R study. Subsequently, section 3 describes the materials and methods used in the development of this work. Experimental and statistical results are presented in section 4. Finally, section 5 presents pertinent and relevant conclusions about the subject. ? 1 Stage 2 Stage 3 Stage 4 Stage 5 Fig. 2. The spot welding cycle showing the five main phases of the process 1 THEORETICAL BACKGROUND 1.1 The Process of RSW Commonly used in large-scale manufacture, the RSW process consists in joining two metal parts through the fusion of the metal, overlaid by two electrodes that generate sufficient force and heat at the weld point, during the passage of an electric current [45]. The welding cycle for RSW presents a series of stages, which are described below: • Stage 1: the electrodes intercept the parts to be welded, providing a certain force (F) on them and ensuring a good settlement; • Stage 2: still under pressure, the electric current (Iw) passes through the system initiating the weld point formation; • Stage 3: after the point being established, the electric current is interrupted, but the mechanical pressure generated by the electrodes is maintained on the pieces until the point solidification; • Stage 4: the exerted force (F) ceases; • Stage 5: the electrodes stop intercepting the parts. The sequence of the welding process is presented in Fig. 2. RSW is controlled by three parameters: welding current, welding time and electrode pressure. These control parameters are presented in several scientific researches which use RSW, such as: Wan et al. [46], Zhang et al. [47], Pouranvari [48], Ighodaro et al. [49], Fan et al. [50], Moos & Vezzetti [51], Amaral et al. [52], Florea et al. [53], Podrzaj & Simoncic [54], Podrzaj et al. [55]. Knowing controllable and uncontrollable process factors, as well as adequately configuring these parameters, helps to ensure a good welded product, since these factors impact on its geometry and final quality [56]. There are several types of tests to evaluate the quality and performance of RSW [57], such as destructive tests, which are usually performed by sampling. Such methods are described by the American Welding Society (AWS) [58], in order to guarantee and monitor the quality of welding nugget, as shown in Fig. 3. It should be noted that in the shear test, the specimens are fixed in specific equipment for tensile tests, in which opposing forces are applied until their rupture, according to the scheme shown in Fig. 3. In this test it is possible to collect the responses of tensile-shear strength (SS) and ultimate strain (US), which may indicate the energy absorbed (toughness) by the specimen. Fig. 3. Shear test on standard specimens after welding 1.2 GR&R Study The total process variation comprises the sum of variations derived from measurement system and part-to-part. The GR&R seeks to estimate how much of the total process variation is caused by the measurement system, and determining how significant it is, compared to the part-to-part variation [59]. Two cases of control can be considered for MS, product and process: • Product: associated with binary decisions, for approved and unapproved products, under 100 % sample inspections or inspections, in which GR&R aims to estimate the tolerance of the product, without verifying the process; • Process: associated to decisions directed to the adequacy of the measurement system to the process control. Searching to stabilize (and understand) the natural process variability in order to make this appropriate. The process can be divided into a certain number of categories and may be represented from a variability statistics of the MS named number of distinct categories (ndc) [60] and [61]. This paper presents a literature review about this subject. The number of distinct categories must be greater or equal to 5. Table 1 presents the guidelines for the MS acceptance [23]. When a single characteristic is evaluated in GR&R studies, a single response variable is analyzed in order to verify the measurement system capability. Besides the traditional univariate approaches, Wang & Chien [31] say that ANOVA, among the known statistical methods, is the best one. ANOVA method stands out in relation to the mean and amplitude (M&A), because it estimates the variance more accurately. In addition, the ANOVA method presents more information about the data than the M&A approach. Table 1. Classification criteria for the measurement system Measurement System GR&R [%] Acceptable < 10 Marginal 10 to 30 Unacceptable > 30 1.2.1 Analysis of Variance (ANOVA) The variability in measurements of univariate cases can be originated from operator mistakes, measuring instrument variation or even from the product itself. Thus, for a complete GR&R study, it is mandatory to follow the model described in Eq. (1) [24] and [26]. Y = p + a, +ßj + (aß ).. +e ijk i = 1,2,..., p j = 1,2,..., o . (1) k = 1,2,..., r In the Eq. (1), Y refers to the response variable, j to the values mean, a, ~N(0, oa) is the random variable for each part, fy ~N(0, op) is the random variable for operator, afyj ~N(0, oa[j) refers to the interaction and sijk ~N(0, os) is the estimated error term. When operators cannot measure all parts, which is a common feature in destructive tests, a NGR&R must be used [62]. This type of design does not present an interaction term between the factors. The variance components for the model are described in Eq. (2). Y = Vn +ßj +a(ßl(J) + £jk i = 1,2,..., p j = 1,2,..., o k = 1,2,..., r (2) where iN is the average values from the nested design, fyj ~N(0, ofy) and a(fi)i(j) ~N(0, o^) are the random and independent variables for operator and for parts nested within operators, respectively. The eijk ~ N(0, o6) is the estimated error term. Still in the Eq. (2), p is the number of parts, o the number of operators and r the number of replicas. The variation components of a NGR&R study with no significant interaction are estimated as in Table 2. Table 2. Variance components for nested design a2 process ^2 a(ß) 2 MSa{ß) - MSs aa(ß) r _2 repeatability a] = MSs a2 reproducibility MSß - MSa(ß) aß pr NGR&R _ _2 repeatability + reproducibility a2 _ _ 2 °total _ a process + NGR&R Table 2 describes the mean squares for operators, parts within operators and the error term. In a nested design, the %R&R and ndc are the two main indicators commonly used to measure and evaluate the MS [60]: • %R&R is the percentage statistics of repeatability & reproducibility (R&R), which measures the MS standard deviation against the total standard deviation, represented by Eq. (3). ndc, also known as the signal-noise index (SNR), measures the variability of the MS. Eq. (4). % NGR & R = •100 [%], N - ndc = 2al \ ®NGR&R = 1.41- ar (3) (4) 2 EXPERIMENTAL PROCEDURE AND PARAMETERS DEFINITION All test specimens used in this study were performed on a stationary classification machine (Presol Transweld® brand, model TWPRV50) with rated power of 50 kVA, AC and maximum current of 6 kA, as shown in Fig. 4. In addition, a chromium-zirconium electrode (Group A, class 2) was used for welding the specimens (0.10 % to 0.15 % C, 0.3 % to 0.6 % Mn, 0.005 % Al, < 0.03 % P, < 0.05 % S, 40 g/m2 to 50 g/m2 Zn) with 0.8 mm thickness. The dimensional specifications of Specimens were made according to the AWS [58] standard, as shown in Fig. 5. failure mode does not occur at the weld point. Table 3 shows the defined values (maximum and minimum). The electrode pressure was set at 2 bar. Table 3. Control factors and respective levels Setup Levels unit -1 + 1 T 1 preheating Cycles 5 11 1preheating % kA 66 74 T 1 welding Cycles 7 17 1welding % kA 75 83 From the limits specified in Table 3, the DOE statistical technique was used to generate the fractional factorial design (FFD) as shown in Table 4, in order to have parts with different characteristics as well as representing the process amplitude. Table 4. Experimental matrix Setup Run T 1 preheating [cycles] 1preheating [% kA] T 1 welding [cycles] 1welding [% kA] 1 11 74 7 83 2 11 66 17 83 3 5 74 7 83 4 11 74 7 75 5 5 66 17 83 6 11 66 17 75 7 5 74 17 83 8 11 74 17 83 Fig. 4. TWPRV50 Presol Transweld® machine Fig. 5. Dimensions of the specimen for the shear test The main welding parameters were established based on preliminary tests, thus the minimum parameter limits ensure that the "interfaciaF type of Fig. 6. Equipment used in the shear test: a) EMIC® (Machine 1) and b) Instron® (Machine 2) Forty-eight welds were performed for shear test. Two operators (or simply the tensile machines) were considered in the NGR&R study: An EMIC® DL2000 (Fig. 6a) with an axial force of 30 kN and an Instron® hydraulic servomotor model 8801 (Fig. V ®total J 6b) with an axial force of 100 kN. Table 5 shows the 48 measurements with two critical-to-quality characteristics (CTQ). 3 RESULTS AND DISCUSSION From the data collected in the shear test, two quality characteristics will be analyzed separately, which are tensile-shear strength (TSS) and ultimate strain (US). Table 5. Measurement of quality responses (CTQ) i k - Machine 1 Machine 2 TSS [N] US [mm] TSS [N] US [mm] 1 1 3217.58 1.47 3737.14 0.3 1 2 2783.24 0.28 3139.77 0.22 1 3 3247.12 0.35 3105.44 0.2 2 1 4779.46 1.45 5030.16 1.14 2 2 5498.73 2.67 5458.82 1.21 2 3 5326.73 1.29 5450.11 1.31 3 1 2096.99 0.17 2100.57 0.09 3 2 1747.78 0.14 2098.47 0.09 3 3 1848.55 0.15 2720.33 0.14 4 1 1228.31 0.09 2539.37 0.11 4 2 2199.49 0.22 2379.70 0.11 4 3 1388.15 0.15 2864.78 0.14 5 1 4484.11 1.16 5580.03 1.35 5 2 4395.51 1.60 4275.08 0.81 5 3 4482.38 1.15 4364.20 0.91 6 1 2413.19 0.22 4329.65 0.9 6 2 2593.87 0.27 3067.30 0.15 6 3 2011.86 0.14 2813.81 0.13 7 1 5385.80 1.97 4983.90 1.16 7 2 5378.85 6.63 4931.34 1.13 7 3 4935.82 1.31 5145.93 1.15 8 1 4656.11 1.20 5838.03 1.37 8 2 4607.46 1.22 4692.68 0.95 8 3 4711.71 1.20 4453.84 0.82 3.1 Results for Tensile-shear Strength In destructive tests, the TSS vector, which holds the set of original responses of the shear test, does not present an interaction term, so it can be represented by Eq. (2). The analysis of variance for the TSS characteristic is found in Table 6. From the ANOVA, the hypothesis of parts being equal must be rejected, but it is not possible to reject the null hypothesis of different operators replicate the same measurement to a specific part, since /»-values are equal to 0.561 and 0.000, respectively. From Eqs. (3) and (4), it is possible to measure the square roots of variances, the %R&R and ndc indicators for the TSS response. Considering the results from Table 7, %R&R presented a value of 29.91 %, being considered marginal, given the conditions established in Table 1. The number of distinct categories (ndc) identified by the system was classified as unacceptable, according to the AIAG [23] recommendations. Table 6. ANOVA for TSS results Source DF SS MS F P Operators 1 1952795 1952795 0.354 0.561 Part (Operators) 14 77233406 5516672 31.5346 0.000 Repeatability 32 5598088 174940 - - Total 47 84784289 - - - Table 7. Variance components contribution of TSS Source G % Contribution ggr&r 418.26 29.91 Grepeatability 418.26 29.91 Greproducibility 0.00 0.00 Gpart-to-part 1334.38 95.42 Gj 1398.40 100 ndc 4 In addition, Fig. 7 shows the result of measurements for the TSS characteristic. It is possible to verify the dispersion between the measured values and the average of these values. For this characteristic, it is possible to verify that the two machines (operators) present different values of average for each part, properly representing the process amplitude. Fig. 7. Measurement results for the TSS characteristic 3.2 Result for Ultimate Strain The same procedure performed for the TSS characteristic was repeated to verify the study for the US quality characteristic. From the analysis of variance for the US characteristic (Table 8) it is possible to verify that it does not reject the null hypothesis of operators replicating an equal measurement for the same part (p-value equal to 0.309). However, the null hypothesis of parts being equal must be rejected, confirming that the choice of parts represents the process amplitude. Table 8. ANOVA for US results Source DF SS MS F P Operators 1 2.3464 2.34638 1.11418 0.309 Part (Operators) 14 29.4829 2.10592 3.41443 0.002 Repeatability 32 19.7367 0.61677 - - Total 47 51.566 - - - Table 9 presents the results concerning the variation components for the US characteristic calculated from Eqs. (3) and (4). Based on the outcomes, it is possible to verify that the %R&R indicator presented a value of 74.4 %, being deemed as unacceptable under the established conditions from Table 1. In addition, the number of distinct categories (ndc) identified by the system presented a value equal to 1, which is considerably unacceptable according to the AIAG [23] recommendations. Table 9. Variance components contribution of US Source G % Contribution ggr&r 0.7917 74.70 Grepeatability 0.7854 74.1 Greproducibility 0.1001 9.44 Gpart-to-part 0.7046 66.48 Gj 1.0598 00 ndc 1 After the results were obtained, the control R chart (Fig. 8a) and boxplot (Fig. 8b) were plotted in order to represent the performance of the two operators to measure the US characteristic. The control R chart highlights the out-of-control point on the Machine 1 for the part 7, which indicates that there was not a considerable measurement repetitiveness of this part in the Machine 1. In addition to that, the boxplot reinforces that there is a measurement problem on the Machine 1, where it is possible to identify the presence of the outlier. In order to identify an outlier in part 7 for the Machine 1, a new analysis was performed for the US characteristic. Given the parameter configuration x*uncoded = [7; 74; 17; 83], a new experiment was carried out under the same conditions, aiming to replace the outlier and perform a new diagnosis (TSS = 5219.7 N; US = 1.72 mm). The results of the new NGR&R study, for US characteristics, showed that the %R&R was 48.31 %, after outlier removal, reducing the variability of the measurement system. However, it is still classified as unacceptable according to AIAG recommended criteria [23]. It is important to highlight that for the TSS there was no statistically significant difference. & c ra ® 1.5 g >3 0.0 7 6 5 ? 4 JE S 3 2 1 0 b) outlier __i ~ ~--=Pl— Machine 1 Machine 2 Operators Fig. 8. a) Control R Chart for US characteristic, and b) boxplot for US characteristic E Jl 1.5 Fig. 9. Measurement dispersions to the US after outlier removal Fig. 9 highlights the metric dispersions for the US characteristic after the removal of an outlier in part 7, evidencing the metrics behavior of each part. Because of outlier removal, the upper limit control (UCL) presented a new range, going from Machine 2 6.0 4.5 a value UCL = 1.755 to UCL = 1.005. Lower limit control (LCL) remains the same. Thus, the control R chart now has two out of control points in part 1 and 2 for Machine 1, as seen in Fig. 10. 3 1.0. R=0.390 LCL=0 1234567812345678 Parts Fig. 10. Control R Chart for US characteristic without outiler 3.3Analysis of the Results After verifying the results from the gage study by using the ANOVA method, it was found that, although one of the characteristics presents a result of variation classified as acceptable, there is a great variability attributed to the Machine 1 (EMIC® machine), presenting less repeatability for the measurements, as well as an out-of-control point for part 7. After detecting this variation cause, the new study showed that the %R&R decreased, but still remained unacceptable for MS. This result can be explained since the Machine 1 did not exhibit an adequate behavior during the tests, as it led to specimens slips, implying in a non-faithful reading of the CTQ even with the proper preparation usually performed for this procedure. It can be inferred, from this analysis, that Machine 1 needs adjustments and improvements, especially with respect to the parts fitting in the machine, since its coupling presented slips during the test due to the metal part composition, which presents a layer of zinc (galvanized steel). Such an improvement could hold the coupling of the part during the shear test in order to not compromise future diagnostics in tests performed by it. Fig. 11 shows the boxplot, which illustrates the form, central tendency and variability of the sample analyzed for the TSS characteristic (Fig. 11a) and US (Fig. 11b). Fig. 12, in turn, presents the result from the gage run chart for the two characteristics analyzed, in which is possible to verify that the test outcomes from Machine 2 present greater homogeneity than the results obtained by the Machine 1. 4 CONCLUSION Machine 1 Machine 2 Operators Machine 1 Machine 2 Operators Fig. 11. Boxplot for a) TSS and b) US a) St TÜ1 Operators Operators Machine 1 Machine 2 Operators Fig. 12. Gage Run Chart for a) TSS and b) US In order to verify the assigned variability of the measuring instrument for the resistance spot welding process. A MSA for the shear test, by using different tensile machines, was performed in this paper, in Machine 1 Machine 2 „ 1.5 UCL = 1.005 2000 2000 Panel variable Pa order to verify the reliability of the outcomes obtained in the RSW process. From the shear test analysis, it was possible to verify that the Machine 1 was responsible for the greater contribution on the system variability, presenting measurement results outside the control, as well as a lower degree of repetitiveness than Machine 2. In addition to that, Machine 1 was not well adjusted, since some specimens have slipped during the tests. This result evidences that some improvements in the parts coupling, which have a coating (such as galvanized steel) are necessary in order to avoid future slips and, therefore, to favor more reliable results, without compromising future diagnoses in tests performed by this type of equipment. Regarding the specimens, it was verified that part 7 presented greater variability in measurements, especially for Machine 1, which may indicate a generalized measurement error for this part. In addition, after identifying this outlier, the new study showed a decrease in the variability of the measurement system to the US characteristic, being able to identify new out of control points in parts 1 and 2. 5 ACKNOWLEDMENTS The authors would like to acknowledge the support from CNPq, CAPES and FAPEMIG. 6 REFERENCES [1] Gomes, G.F., Vieville, P., Durrenberger, L. (2017) Dynamic behavior investigation of spot welding machines and its influence on weld current range by modal analysis. Journal of the Brazilian Society of Mechanical Sciences and Engineering, vol. 39, no. 3, p. 765-773, DOI:10.1007/s40430-016-0580-0. [2] Zhang, Y., Li, Y., Luo, Z., Yuan, T., Bi, J., Wang, Z.M., Wang, Z.P., Chao, Y.J. (2016). Feasibility study of dissimilar joining of aluminum alloy 5052 to pure copper via thermo-compensated resistance spot welding. Materials & Design, vol. 106, p. 235246, DOI:10.1016/ j.matdes.2016.05.117. [3] Bi, J., Song, J., Wei, Q., Zhang, Y., Li, Y., Luo, Z. (2016). Characteristics of shunting in resistance spot welding for dissimilar unequal-thickness aluminum alloys under large thickness ratio. Materials & Design, vol. 101, p. 226-235. DOI:10.1016/j.matdes.2016.04.023. [4] Li, Y.B., Li, Y.T., Shen, Q., Lin, Z.Q. (2013). Magnetically assisted resistance spot welding of dual-phase steel. Welding Journal, vol. 92, no. 4, p. 124-132. [5] Feng, Y., Li, Y., Luo, Z., Ling, Z., Wang, Z. (2016). Resistance spot welding of Mg to electro-galvanized steel with hot-dip galvanized steel interlayer. Journal of Materials Processing Technology, vol. 236, p. 114-122, DOI:10.1016/j. jmatprotec.2016.05.015. [6] Zhang, Y., Luo, Z., Li, Y., Liu, Z., Huang, Z. (2015). Microstructure characterization and tensile properties of Mg/ Al dissimilar joints manufactured by thermo-compensated resistance spot welding with Zn interlayer. Materials & Design, vol. 75, p. 166-173, DOI:10.1016/j.matdes.2015.03.030. [7] Martin, 0., De Tiedra, P., San-Juan, M. (2017). Combined effect of resistance spot welding and precipitation hardening on tensile shear load bearing capacity of A286 superalloy. Materials Science and Engineering: A, vol. 688, p. 309-314, DOI:10.1016/j.msea.2017.02.015. [8] Shan, H., Zhang, Y., Li, Y., Luo, Z. (2017). Dissimilar joining of AZ31B magnesium alloy and pure copper via thermo-compensated resistance spot welding. Journal of Manufacturing Processes, vol. 30, p. 570-581, DOI:10.1016/j. jmapro.2017.10.022. [9] Chen, N., Wang, H. P., Carlson, B. E., Sigler, D. R., Wang, M. (2017). Fracture mechanisms of Al/steel resistance spot welds in lap shear test. Journal of Materials Processing Technology, vol. 243, p. 347-354, DOI:10.1016/j.jmatprotec.2016.12.015. [10] Manladan, S.M., Yusof, F., Ramesh, S., Zhang, Y., Luo, Z., Ling, Z. (2017). Microstructure and mechanical properties of resistance spot welded in welding-brazing mode and resistance element welded magnesium alloy/austenitic stainless steel joints. Journal of Materials Processing Technology, vol. 250, p. 45-54, DOI:10.1016/j.jmatprotec.2017.07.006. [11] Podrzaj, P., Polajnar, I., Diaci, J., Kariz, Z. (2004). Expulsion detection system for resistance spot welding based on a neural network. Measurement Science and Technology, vol. 15, no. 3, p. 592-598, DOI:10.1088/0957-0233/15/3/011. [12] Podrzaj, P., Polajnar, I., Diaci, J., Kariz, Z. (2005). Estimating the strength of resistance spot welds based on sonic emission. Science and Technology of Welding and Joining, vol. 10, no. 4, p. 399-405, DOI:10.1179/174329305X44107. [13] Podrzaj, P., Polajnar, I., Diaci, J., Kariz, Z. (2006). Influence of welding current shape on expulsion and weld strength of resistance spot welds. Science and Technology of Welding and Joining, vol. 11, no. 3, p. 250-254, DOI:10.1179/174329306X101391. [14] Podrzaj, P., Simoncic, S. (2013). Resistance spot welding control based on the temperature measurement. Science and Technology of Welding and Joining, vol. 18, no. 7, p. 551-557, DOI:10.1179/1362171813Y.0000000131. [15] Simoncic, S., Podrzaj, P. (2012). Image-based electrode tip displacement in resistance spot welding. Measurement Science and Technology, vol. 23, no. 6, p. 1-7, DOI:10.1088/0957-0233/23/6/065401. [16] Simoncic, S., Podrzaj, P. (2014). Resistance spot weld strength estimation based on electrode tip displacement/velocity curve obtained by image processing. Science and Technology of Welding and Joining, vol. 19, no. 6, p. 468-475, DOI:10.1179/ 1362171814Y.0000000212. [17] Podrzaj, P., Polajnar, I., Diaci, J., Kariz, Z. (2008). Overview of resistance spot welding control. Science and Technology of Welding and Joining, vol. 13, no. 3, p. 215-224, DOI:10.1179/174329308X283893. [18] Majeske, K.D. (2008). Approval criteria for multivariate measurement systems. Journal of Quality Technology, vol. 40, no. 2, p. 140-153, DOI:10.1080/00224065.2008.11917721. [19] Pereira, R.B.D., Peruchi, R.S., de Paiva, A.P., da Costa, S.C., Ferreira, J.R. (2016). Combining Scott-Knott and GR&R methods to identify special causes of variation. Measurement, vol. 82, p. 135-144, D0l:10.1016/j.measurement.2015.12.033. [20] Woodall, W.H., Borror, C.M. (2008). Some relationships between gage R&R criteria. Quality and Reliability Engineering International, vol. 24, no. 1, p. 99-106, D0I:10.1002/qre.870. [21] Senol, S. (2004). Measurement system analysis using designed experiments with minimum a-p Risks and n. Measurement, vol. 36, no. 2, p. 131-141, D0I:10.1016/j. measurement.2004.05.001. [22] Peruchi, R.S., Paiva, A.P., Balestrassi, P.P., Ferreira, J.R., Sawhney, R. (2014). Weighted approach for multivariate analysis of variance in measurement system analysis. Precision Engineering, vol. 38, no. 3, p. 651-658, D0I:10.1016/j.precisioneng.2014.03.001. [23] Automotive Industry Action Group - AIAG (2010). Measurement systems analysis - Reference Manual. 4th ed. Chrysler, Ford, General Motors Supplier Quality Requirements Task Force. [24] Al-Refaie, A., Bata, N. (2010). Evaluating measurement and process capabilities by GR&R with four quality measures. Measurement, vol. 43, no. 6, p. 842-851, D0I:10.1016/ j.measurement.2010.02.016. [25] Awad, M., Erdmann, T.P., Shanshal, Y., Barth, B. (2009). A measurement system analysis approach for hard-to-repeat events. Quality Engineering, vol. 21, no. 3, p. 300-305, D0I:10.1080/08982110902852344. [26] Burdick, R.K., Borror, C.M., Montgomery, D.C. (2003). A review of methods for measurement systems capability analysis. Journal of Quality Technology, vol. 35, no. 4, p. 342-354, D0I:10.1080/00224065.2003.11980232. [27] Wu, C.W., Pearn, W.L., Kotz, S. (2009). An overview of theory and practice on process capability indices for quality assurance. International Journal of Production Economics, vol. 117, no. 2, p. 338-359, D0I:10.1016/j.ijpe.2008.11.008. [28] Erdmann, T.P., Does, R.J., Bisgaard, S. (2009). Quality quandaries*: A gage R&R study in a hospital. Quality Engineering, vol. 22, no. 1, p. 46-53, D0I:10.1080/08982110903412924. [29] Knowles, G., Vickers, G., Anthony, J. (2003). Implementing evaluation of the measurement process in an automotive manufacturer: a case study. Quality and Reliability Engineering International, vol. 19, no. 5, p. 397-410, D0I:10.1002/qre.533. [30] Polini, W., Turchetta, S. (2004). Test protocol for micro-geometric wear of sintered diamond tools. Wear, vol. 257, no. 3-4, p. 246-256, D0I:10.1016/j.wear.2003.12.008. [31] Wang, F.K., Chien, T.W. (2010). Process-oriented basis representation for a multivariate gauge study. Computers & Industrial Engineering, vol. 58, no. 1, p. 143-150, D0I:10.1016/j.cie.2009.10.001. [32] Shi, L., Chen, W., Lu, L.F. (2014). An approach for simple linear profile gauge R&R studies. Discrete Dynamics in Nature and Society vol. 2014, D0I:10.1155/2014/816980. [33] Deshpande, A.A., Ramya, A., Vishweshwar, V., Deshpande, G.R., Roy, A.K. (2014). Applications of gage reproducibility & repeatability (GRR): understanding and quantifying the effect of variations from different sources on a robust process development. Organic Process Research & Development, vol. 18, no. 12, p. 1614-1621, D0I:10.1021/op5002935. [34] Zhu, X., Zhao, Z., Wang, L., Zhang, L. (2014). A new method to measure fat content In coconut milk based on Y-type optic fiber system. Optik-lnternational Journal for Light and Electron Optics, vol. 125, no. 20, p. 6172-6178, D0I:10.1016/ j.ijleo.2014.06.115. [35] Garcia, A.C., del Rio, A.G. (2013). Number of distinct data categories and gage repeatability and reproducibility. A double (but single) requirement. Measurement, vol. 46, no. 8, p. 2514-2518, D0I:10.1016/j.measurement.2013.04.065. [36] Johnson, J.A., Widener, S., Gitlow, H., Popovich, E. (2006). A "Six Sigma"© black belt case study: GEP Box's paper helicopter experiment part A. Quality Engineering, vol. 18, no. 4, p. 413-430, D0I:10.1080/08982110600875894. [37] Wan, X., Wang, Y., Zhao, D., Huang, Y., Yin, Z. (2017). Weld quality monitoring research in small scale resistance spot welding by dynamic resistance and neural network. Measurement, vol. 99, p. 120-127, D0I:10.1016/j. measurement.2016.12.010. [38] Degidi, M., Caligiana, G., Francia, D., Liverani, A., Olmi, G., Tornabene, F. (2016). Strain gauge analysis of implant-supported, screw-retained metal frameworks: Comparison between different manufacturing technologies. Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, vol. 230, no. 9, p. 840-846, D0I:10.1177/0954411916653623. [39] Wang, L., Hou, Y., Zhang, H., Zhao, J., Xi, T., Qi, X., Li, Y. (2016). A new measurement method for the dynamic resistance signal during the resistance spot welding process. Measurement Science and Technology, vol. 27, no. 9, p. 1-12, D0I:10.1088/0957-0233/27/9/095009. [40] Xia, Y.J., Zhang, Z.D., Xia, Z.X., Zhu, S.L., Zhang, R. (2015). A precision analogue integrator system for heavy current measurement in MFDC resistance spot welding. Measurement Science and Technology, vol. 27, no. 2, p. 1-11, D0I:10.1088/0957-0233/27/2/025104. [41] Al-Jader, M.A., Cullen, J.D., Shaw, A., Al-Shamma'a, A.I. (2011). Theoretical and practical investigation into sustainable metal joining process for the automotive industry. Journal of Physics: Conference Series, vol. 307, no. 1, p. 012044, D0I:10.1088/1742-6596/307/1/012044. [42] Lei, Z., Kang, H.-T., Reyes, G. (2010). Full field strain measurement of resistant spot welds using 3D image correlation systems. Experimental Mechanics, vol. 50, no. 1, p. 111-116, D0I:10.1007/s11340-008-9186-5. [43] Lai, X.M., Luo, A.H., Zhang, Y.S., Chen, G.L. (2009). Optimal design of electrode cooling system for resistance spot welding with the response surface method. The International Journal of Advanced Manufacturing Technology, vol. 41, no. 3-4, p. 226-233, D0I:10.1007/s00170-008-1478-5. [44] Almeida, F.A., Gomes, G.F., De Paula, V.R., Correa, J.E., Paiva, A.P., Gomes, J.H.F., Turrioni, J.B. (2018). A weighted mean square error approach to the robust optimization of the surface roughness in an AISI 12L14 free-machining steel-turning process. Strojniški vestnik - Journal of Mechanical Engineering, vol. 64, no. 3, p. 147-156, D0I:10.5545/sv-jme.2017.4901. [45] Zhang, H., Senkara, J. (2011). Resistance Welding: Fundamentals and Applications. CRC press, Boca Raton. [46] Wan, Z., Wang, H.P., Chen, N., Wang, M., Carlson, B.E. (2017). Characterization of intermetallic compound at the Interfaces of Al-steel resistance spot welds. Journal of Materials Processing Technology, vol. 242, p. 12-23, DOI:10.1016/j. jmatprotec.2016.11.017. [47] Zhang, Y., Shan, H., Li, Y., Guo, J., Luo, Z., Ma, C.Y. (2017). Joining aluminum alloy 5052 sheets via novel hybrid resistance spot clinching process. Materials & Design, vol. 118, p. 36-43, DOI:10.1016/j.matdes.2017.01.017. [48] Pouranvari, M. (2017). Fracture toughness of martensitic stainless steel resistance spot welds. Materials Science and Engineering: A, vol. 680, p. 97-107, DOI:10.1016/j. msea.2016.10.088. [49] Ighodaro, O.L., Biro, E., Zhou, Y.N. (2016). Comparative effects of Al-Si and galvannealed coatings on the properties of resistance spot welded hot stamping steel joints. Journal of Materials Processing Technology, vol. 236, p. 64-72, DOI:10.1016/j.jmatprotec.2016.03.021. [50] Fan, Q., Xu, G., Gu, X. (2016). Expulsion characterization of stainless steel resistance spot welding based on dynamic resistance signal. Journal of Materials Processing Technology, vol. 236, p. 235-240, DOI:10.1016/j.jmatprotec.2016.05.026. [51] Moos, S., Vezzetti, E. (2015). Resistance spot welding process simulation for variational analysis on compliant assemblies. Journal of Manufacturing Systems, vol. 37, p. 44-71, DOI:10.1016/j.jmsy.2015.09.006. [52] Amaral, F.F., Almeida, F.A., Costa, S.C., Leme, R.C., Paiva, A.P. (2018). Application of the Response Surface Methodology for Optimization of the Resistance Spot Welding Process in AISI 1006 Galvanized Steel. Soldagem & Inspeqao. DOI:10.1590/0104-9224/SI2302.02. [53] Florea, R.S., Bammann, D.J., Yeldell, A., Solanki, K.N., Hammi, Y. (2013). Welding parameters influence on fatigue life and microstructure in resistance spot welding of 6061-T6 aluminum alloy. Materials & Design, vol. 45, p. 456-465, DOI:10.1016/j.matdes.2012.08.053. [54] Podrzaj, P., Simoncic, S. (2011). Resistance spot welding control based on fuzzy logic. The International Journal of Advanced Manufacturing Technology, vol. 52, no. 9-12, p. 959-967, DOI:10.1007/s00170-010-2794-0. [55] Podrzaj, P., Jerman, B., Simoncic, S. (2016). Poor fit-up condition in resistance spot welding. Journal of Materials Processing Technology, vol. 230, p. 21-25, DOI:10.1016/j. jmatprotec.2015.11.009. [56] Zhou, M., Zhang, H., Hu, S. J. (2003). Relationships between quality and attributes of spot welds. Welding Journal, vol. 82, no. 4, p. 72-S-77-S. [57] Darwish, S.M., Al-Dekhial, S.D. (1999). Micro-hardness of spot welded (B.S. 1050) commercial aluminium as correlated with welding variables and strength attributes. Journal of Materials Processing Technology, vol. 91, no. 1-3, p. 43-51, DOI:10.1016/S0924-0136(98)00414-2. [58] AWS/SAE (2002). Recommended Practices for Test Methods for Evaluating the Resistance Spot Welding Behavior of Automotive Sheet Steel Materials. American Welding Society, Florida. [59] Shiau, Y.R. (2001). Decision support for off-line gage evaluation and improving on-line gage usage. Journal of Manufacturing Systems, vol. 19, no. 5, p. 318-331, DOI:10.1016/S0278-6125(01)89004-X. [60] Montgomery, D.C., Burdick, R.K., Lawson, C.A., Molnau, W.E., Zenzen, F., Jennings, C.L., Shah, H.K., Sebert, D.M., Bowser, M.D., Holcomb, D.R. (2005). A University-based Six Sigma Program. Quality and Reliability Engineering International, vol. 21, no. 3, p. 243-248, DOI:10.1002/qre.631. [61] White, T.K., Borror, C.M. (2011). Two-dimensional guidelines for measurement system indices. Quality and Reliability Engineering International, vol. 27, no. 4, p. 479-487, DOI:10.1002/qre.1144. [62] Aquila, G., Peruchi, R.S., Rotela, P.J., Rocha, L. C.S., Queiroz, A.R., Pamplona, E.O., Balestrassi, P.P. (2018). Analysis of the wind average speed in different Brazilian states using the nested GR&R measurement system. Measurement, vol. 115, p. 217-222, DOI:10.1016/j.measurement.2017.10.048. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)10, 632-642 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2018.5397 Original Scientific Paper Received for review: 2018-04-03 Received revised form: 2018-06-01 Accepted for publication: 2018-07-03 Performance Characteristics Research on Power System of a Four-Valve Compressed Air Engine Qihui Yu* - Xin Tan - Jianguo Meng - Fei Shang - Xueqing Hao Inner Mongolia University of Science & Technology, School of Mechanical Engineering, China To alleviate air pollution and eliminate pollutants from internal combustion engines, many researchers have paid more attention to compressed air engines (CAEs). However, the low energy utilization efficiency and low output power of a CAE hinder its development. In this paper, a new kind of CAE, which has four valves is proposed, which reduces the kinetic energy of the flowing air. Through analysis of the CAEs, a mathematical model of the working processes was developed. Using the simulation method, the performance of the new kind of CAE was obtained. To verify the theoretical model, a prototype of a four-valve of single cylinder piston type CAE was designed, and the experiment was carried out. Parameter influence analysis and optimization were studied to improve the four-valve CAE performance. Results show that, first, the output torque and output power of the four-valves of CAE are higher than of the two-valves of CAE. When the intake pressure is set at 0.6 MPa and 0.7 MPa, and the rotation speed is 450 rpm, the output torque of the four-valve of CAE is about 2.7 and 2.2 times that of the two-valve of CAE. The output torque and power differences between the two kinds of CAE ascend sharply with increased rotation speed. Second, the energy efficiency and output power are impacted significantly by the intake valve close angle and exhaust passage diameter. Thirdly, the optimization of the four-valve CAE can significantly enhance the output power compared with two-valve CAE. This research can be referred to in the design of the CAEs. Keywords: four-valve compressed air engine, energy efficiency, output power, optimisation Highlights • A design of a four-valve compressed air engine to obtain high output power and efficiency is presented. • A prototype of a four-valve single cylinder piston type compressed air engine is designed, and experiments are carried out. • Parameter influence analysis and optimization are studied to improve the four-valve compressed air engine's performance. • The proposed four-valve compressed air engine can enhance significantly the output power in comparison to two-valve compressed air engine. 0 INTRODUCTION One major source of air pollutants is generated from internal combustion engines (ICE) for vehicles [1]. It has been reported that the pollutants produced by ICE are responsible for more than 30 % of air pollutants. To alleviate air pollution and eliminate pollutants from ICE, various alternatives, including battery power driven engines [2], hydrogen fuel cell-based engines [3], liquid nitrogen fuel power driven engines, and hybrid energy driven engines [4] and [5] have been developed. However, looking at these various alternatives, each option has its limitations. For example, battery-operated engines suffer from long charging time, high initial cost, heavy metal pollution, and periodical replacement cost [6]. Distribution and storage are difficult for hydrogen fuel cell driven engines, and building hydrogen generator stations is expensive [7]. Liquid nitrogen vehicles are difficult to use given their high cost. The main power source of hybrid energy driven engines remain fossil fuel, so this kind of engine does not meet zero pollution vehicle requirements [8]. In these situations, the technology of the compressed air engine (CAE) and its utilization in vehicles has become an attractive option, owing to its potential of zero pollution. However, in early tests, the vehicle running on compressed air alone was limited to the storage capacity of the tanks. To obtain high energy density and store more energy in limited space, the storage pressure reaches several hundred bars. Once the compressed air is released from the tanks, its pressure falls to 30 bar as it passes through the pressure reducer. Even at high pressure, the compressed air carries much less energy than other transportation energy sources. Therefore, the CAE applications will need to overcome limitations in compressed air energy density. To enhance the output power of CAE and improve the energy efficiency, numerous studies have considered the high efficiency and high output power energy conversion system for CAE. Shen et al. built an air-powered motorcycle prototype with a fuzzy logic speed controller; the efficiency is above 70 % [1]. Liu et al. investigated a new rotary intake and exhaust system for a piston type CAE and obtained 2.15 kW output power at 13 bar [9]. To overcome leakage under high air supply pressure, a new pressure-compensated intake valve was designed by Yu et al. [10]. Zhang et *Corr. Author's Address: Inner Mongolia University of Science & Technology, School of Mechanical Engineering, Arding Road No. 7, Baotou, Inner Mongolia, China, yqhhxq@163.com al. studied a single-valve reciprocating expander and pointed out that the pressure loss is mainly caused by the air intake process [11]. Xu et al. designed a heat exchange system of the CAE to prevent ice blocking and improve energy efficiency [12].Yu et al. analysed the influences of some parameters and optimized the energy efficiency and output power with improved NSGA-II [13]. The output power and energy efficiency are significantly influenced by the open and close duration, angle and pressure of the intake and exhaust valves [14].To apply of CAEs as the main power system in automobiles, modifications should be made on a compressed air distribution system of CAE. In this paper, a four-valve CAE is proposed to obtain more energy during the working process. The challenge of this work lies in the energy efficiency. Thus, a multi-objective optimal method is applied to improve energy efficiency. The study aims to obtain a CAE whose energy efficiency is above 20 %, and output power is over 5.5 kW. 1 ENERGY EFFICIENCY OF THE TRADITIONAL CAE The energy of compressed air could be converted to mechanical work during the state change process of compressed air. The working process of a single cylinder CAE could be divided into three stages in Fig. 1. The working process has been illustrated in the previous literature. Intake Expand Fig. 1. Working cycle of CAE For CAE, the ideal scenario for performing the maximum level of mechanical work is described as three ideal processes: quasistatic intake processes, isothermal expand process, and quasistatic exhaust process. According to the literature [15], the work done in these processes can be calculated by Eq. (1). V. = pV ln p, Pa where Wt is the maximum level of mechanical work, ps is the pressure of supply air, Vs is the volume of supply air, and pa is the environment pressure. To determine energy losses and investigate CAE internal energy distributions, the air power is introduced in this paper, which is described by Eq. (2). P = p G ln^. Pa (2) According to the exact definition of air power, the kinetic power of the flowing air, isothermal expansion power and exhaust power are included in the air power calculation. The kinetic energy, expansion, energy and exhaust energy can be expressed with the following equations. „ d i 1 Pk = — I — mv Pe = dt i 2 rv2 JV pv dV V K dt Pd = paGa ln p . (3) (4) (5) In the above equations, Pk, PE and Pd are the kinetic power, expansion power and exhaust power, respectively, V1 and V2 are the volumes of beginning expansion and finishing expansion respectively, pd is the pressure of exhaust air. Fig. 2 displays the ratio of Pk in intake air power. When the average air velocity is below 100 m/s, the pressure of air is above 0.3 MPa, and kinetic energy accounts for less than 5 % of the available energy. When the average air velocity is above 200 m/s, kinetic energy accounts for more than 18 % of the available energy. Thus, when the effective cross-sectional area is narrow, and the ratio of the input pressure and output pressure is below the critical pressure ratio, the kinetic energy should be considered. Fig. 3 shows expansion energy loss in air power. The bigger the polytropic exponent is, the larger the expansion loss is. Therefore, achieving quasi-isothermal expansion is an important measure to improve energy efficiency and output power. It is obvious that the higher the exhaust pressure is, the greater the exhaust energy loss is. Theoretically, energy loss occurs in an irreversible process. Therefore, irreversible processes in CAE systems will reduce the air power. The irreversible process will influence the energy efficiency of a CAE. 2 a 40 _ 35 £ 30 £25 ° 20 o pS 10 ii V -v=340 m/s ■ ' \ ----v=200 m/s 11 \ M \ ! 1 \ ---------v=100 m/s ----v=50 ni's \ I \ t \ |! \ J v ■i. i V . -r — - -----_ ----- 0.1 0.5 0.9 1.3 1.7 2.1 P . [MPa] Fig. 2. The ratio of kinetic energy in air power Polytropic exponent Fig. 3. The ratio of expansion energy loss in air power In this paper, the efficiency of a CAE is defined as the ratio of the output shaft energy to the input air power of CAE as follows [16]. 2 THE NEW KIND OF CAE WITH FOUR VALVES To have precise, continuous valve timing control and fast response, the flow of air moving into and out of the cylinder is controlled by using a simple cam mechanism. The effective area of intake and the exhaust port is determined by cam lift. According to the ratio of the pressure of the upstream side pu to the pressure of downstream pressure pd, the flow equation for the flow through a restriction can be written as follows [17]: dm dty AePu$ (6) where 0 = 2 R i k + 1 k+l k-l Pd < 0.528 2 k+l " 2k l f Pd jK l Pd j k Pd (k- 1) R l Pu j l Pu j Pu where m is the rotation speed of the crankshaft. The air mass flowing into the cylinder is inversely proportional to the rotation speed of the crankshaft. Therefore, the higher the speed of the crankshaft, the lower the output torque. To improve the performance under high speed, increasing the effective area is an effective method. To increase the effective area, a four-valve CAE is proposed. The schematic drawing of the proposed four-valve mechanical structure is illustrated in Fig. 4. Rotations of a crankshaft of CAE are transmitted to both inlet and exhaust cams by using a timing belt and timing pulleys. The cam lift follows the cam profile when appropriate preload springs are used [18]. Fig 4. The schematic drawing of the four-valve mechanical structure The four-valve CAE, which has been applied in the field of traditional internal combustion engines, has a bigger flow effective area and provides more power than the two-valves CAE. Adding more valves improves the mass flow of intake and exhaust compressed air, thereby enhancing output power at high speed. The principle of the four-valve CAE is same as the two-valve CAE, which is introduced in the literature [16]. k u 3 THERMODYNAMIC MODEL AND EXPERIMENT To model a pneumatic system, the law of conservation of energy, the equation of continuity and the equation of state of the ideal gas is fundamental [19]. 3.1 The Volume The volume of the cylinder is variable due to the movement of the piston. From the piston-crank geometry, the piston displacement, y, is expressed by y = L + r - L cos ß- r cosç, where ß = sin r sm^ L ' The volume of the cylinder is given by V = nD 2 y (8) (9) (10) transfer coefficient which is a constant value,^ is the heat transfer area. 3.5 The Performance The output power and energy efficiency are important evaluate indicators to engine. The output power can be expressed as: ^ (14) P =-e 9550 T = IT where i id r f p , (15) (16) where Pe is the output power, Ti is the instantaneous torque, Tid is the indicated torque, Tr is the reciprocating torque, and Tf is the friction torque. Details of torque calculation are presented in the literature [20]. 3.2 The Mass According to the law of mass conservation, air mass in the cylinder in one cycle can be given as: — = p(G. -G ), (11) dt \ in out ] where subscript in and out represent the fact that compressed air flows into and out the cylinder. 3.3 The Pressure Pressure changes of the air inside the cylinder can be calculated by deriving the state equation of ideal gases: dp 1 dt V R0G + mRpdV dt dt (12) 3.4 The Energy The temperature differential can be calculated by solving the energy equation. For CAE system, intake process and exhaust process are the open system, expand process is a closed system. For ideal air, accordingly, the energy balance is written by, C d° h A C m— = h A v dt e w (-o)+( ROG dV -ROG - p-, out dt c o - c o)r. pa v J it (13) where Cv is the specific heat at constant volume, Cp is the specific heat at constant pressure, he is the heat 3.6 Performance of the Four-Valve CAE Compared with the two-valve CAE, the four-valve CAE can achieve a low Mach number, which can reduce the kinetic energy of the flowing air, and the four-valve mechanism can increase the mass flow rate, which can enhance the output power. To obtain the characteristic of four-valve CAE, the initial values of the four-valve CAE mechanism are shown in Table 1. For simplification, the intake valve and exhaust valve have the same structural parameter. MATLAB/ SIMULINK software is used for modelling the simulation. Fig. 5 depicts the output power and energy efficiency characteristics of the multi-valve CAE and the two-valve CAE at different rotational speeds when the supply pressure is equal to 2 MPa. Table 1. Four-valve of CAE specifications Type Value Bore 85 mm Stroke 88 mm Displacement 0.5 L Intake port diameter 16 mm Exhaust port diameter 16 mm Compression ratio 10 Minimum diameter of passage 12 mm Number of intake valve 2 Number of exhaust valve 2 As can be seen from Fig. 5, the curve of the energy efficiency of the four-valve CAE is slightly larger than that of the two-valve CAE, and the output power t 4 of the four-valve CAE is higher than that of the two-valve CAE by approximately 0.6 kW, which is about 1.1 to 2 times that of the two-valve CAE. Especially when the rotational speed is 600 rpm, the highest output power is 3.6 kW. a) b) Fig. 5. Efficiency and power of CAE in different rotational speed; a) curves of efficiency, b) curves of power The output power and energy efficiency characteristics of the four-valve CAE have been studied by setting the rotational speed of the four-valve CAE at 600 rpm, under the supply pressure of the four-valve, at 1 MPa to 3 MPa. Fig. 6 illustrates the relationship between the performances and supply pressure of the four-valve CAE. With an increase in the supply pressure, the output power of the four-valve CAE increases stably. In the first instance, the energy efficiency increases with an increase in the supply pressure; when the supply pressure is larger than 1.8 MPa, the energy efficiency slowly drops. The primary reason is that the indicated pressure inside the cylinder increases with an increase in the supply pressure. The higher the indicated pressure inside the cylinder is, the larger output power is. However, the residual pressure inside the cylinder also increases when the supply pressure is larger than 1.8 MPa under the exhaust process, which leads to exhaust energy loss. 2.2 Perssure , [MPa] Fig. 6. Efficiency and power of four-valve CAE in different supply pressure 3.7 Experiment and Analysis To verify the above theoretical model, we designed a multi-valve single-cylinder piston type CAE model and prototype, which are shown in Fig. 7. The effective areas of intake and exhaust valves are shown in Fig. 8. Basic parameters are shown in Table 2. Fig. 7. The model and prototype of multi-valves of CAE The friction coefficient of the ring assembly is important to calculate friction torque. The correct oil can keep an engine running smoothly. Because the principles of the CAE are different from ICE, it is necessary to find the correct lubricating oil to keep the CAE running smoothly and effectively. According to the literature [21], the lower the temperature is, the higher the viscosity is. The working temperature range of CAE is 220 K to 293 K, so we adopted refrigerator oil as a lubricating oil; its viscosities value are 15 mm2/s and 3.1 mm2/s when the temperatures are 313 K and 373 K, respectively. 180 270 Crank angle ,[°] Fig. 8. Effective areas of intake valve and exhaust valve Table 2. Single cylinder CAE specifications Type Value Displacement 0.5 L Number of valves 4 Opening angle of intake valve 0° Closing angle of intake valve 90° Maximum cam lift of intake valve 3.5 mm Opening angle of exhaust valve 180° Closing angle of exhaust valve 360° Maximum cam lift of exhaust valve 5.72 mm Intake port diameter 22 mm Exhaust port diameter 22 mm Minimum diameter of passage 14 mm To compare performance, the output power and torque of two-valves of CAE and multi-valves of CAE, the experiment is carried out under same working conditions. The experimental apparatus and the procedure have been introduced in the literature [10]. The intake pressure is set at 0.6 MPa and 0.7 MPa. Experiments of the output torque and output power of the two kinds of CAEs were done, and the results are shown in Fig. 9. Fig. 9a illustrates the relationship between the output torque and the rotation speed of the two kinds of CAEs. As can be seen, the curves of the output torque of the two kinds of CAEs decline sharply with the rotational speed. Moreover, the output torque of the four-valve of CAE is much higher than that of the two-valve of CAE in the same intake pressure and rotational speed. When the intake pressure is set at 0.6 MPa, and the rotational speed is 450 rpm, the output torque of the four-valve CAE is about 2.7 times that of the two-valve CAE. When the intake pressure is set at 0.7 MPa, and the rotational speed is 450 rpm, the output torque of the four-valve CAE is about 2.2 times that of the two-valve CAE. Moreover, the output torque difference of the two kinds of CAE ascends sharply with increasing rotational speed. a) 0.6MPa four-valve 0.6MPa two-valve ©— 0.7MPa four-valve G— 0.7MPa two-valve 600 800 1000 1200 y Speed, [rpm] Fig. 9. The output torque and power of the CAEs in different intake pressure; a) torque vs speed and b) power vs speed Fig. 9b illustrates the relationship between the output power and the rotational speed of the two kinds of CAEs. As can be seen, the curves of the output power of the two kinds of CAEs decline sharply with the rotational speed. Moreover, the output power of the four-valve CAE is much higher than that of the two-valve CAE in the same intake pressure and rotational speed. The output power difference of two kinds of CAEs ascends sharply with increased rotational speed. 4 PARAMETER INFLUENCE ANALYSIS In conventional mechanical valve trains, intake and exhaust valves lift are determined by the cam profile. In this paper, the cam profile is described by a five polynomial curve which is given by: h (a) = C„ + C ap + C aq + C ar + C ^ ' 0 d a r s where C0 h„ (17) (18) (19) (20) where x is the cam angle, xB is the half wrap angle, superscript p, q, r, s are plural. For simplification, we suppose the following equations [22]: p = 2n q = 2n +l r = 2n + 2l s = 2n + 4l where, n, l are positive integers. Cp, Cq, Cr, Cs can be calculated by boundary conditions [23] which are functions of n, l. In this paper, intake cam profile variable nt and lt are equal to 3 and 6, respectively, and exhaust cam profile variable ne and le are equal to 6 and 3, respectively. In an actual four-valve CAE system, the compressed air charged from the buffer container flows into the cylinder through the regulator, and the air mass that flows into the cylinder is limited by the effective sectional area of the regulator. The values of the initial parameters are shown in Table 2. The rotation speed is 600 rpm, and the supply pressure is 3 MPa. Except the above parameters, according to the thermodynamic model, the output power and energy efficiency are influenced by cylinder clearance distance, LC, is the distance between the bottom of the valve and the top of the piston when the piston reaches to the top dead centre, intake cam maximum lift, IVL, intake valve close angle, IVC, intake valve head diameter, IHD, exhaust cam maximum lift, EVL, exhaust valve close angle, EVC, exhaust passage diameter, EPD, exhaust valve head diameter, EHD. The initial values of the above parameters are shown in Table 3. Table 3. The initial values of the parameters Parameter LC IVL IVC IHD Value 7 mm 3.2 mm 90° 16 mm Parameter EVL EVC EPD EHD Value 5.72 mm 360° 12 mm 16 mm According to the thermodynamic model above, each parameter can be changed for comparison while all other parameters are kept constant, and the simulation results are varying each parameter, as illustrated in Fig. 10. It can be seen from Fig. 10 that: a) The energy efficiency and output power are impacted significantly by the intake valve close angle and exhaust passage diameter. The energy efficiency decreases from 28 % to 20 % with the intake valve close angle increasing from 70° to 130°. In the first instance, the output power increases with an increase in the intake valve close angle; when the intake valve close angle is larger than 110°, the output power begins to decrease. When the exhaust passage diameter increases from 10 mm to 16 mm, the energy efficiency increases from 18.2 % to 25 %, and the output power increases from 4 kW to 6.7 kW. b) The energy efficiency and output power are slightly affected by intake cam maximum lift, exhaust cam maximum lift and exhaust valve close angle. c) Other parameters have little influence on the output power and energy efficiency. 5 OPTIMUM PERFORMANCE PARAMETERS According to the above analysis, the output power and energy efficiency of the four-valve CAE are influenced by intake valve close angle, exhaust passage diameter, intake cam maximum lift, exhaust cam maximum lift and exhaust valve close angle. Therefore, in this section, the above five parameters will be optimized. Meanwhile, to keep the four-valve CAE running smoothly, the average temperature inside the cylinder is used as another objective function. The coupling between parameters has no explicit expressions, according to the literature [15]. The independent optimal of individual parameters shows that the output power, energy efficiency, and temperature of a four-valve CAE cannot achieve X a = x B a) b) d) g) h) 17 18 19 EHD, [mm] Fig. 10. Influence of parameters on the efficiency and power optimum simultaneously. The ranges of parameters are shown in Table 4. Table 4. Parameters range Parameters Range Intake cam maximum lift (IVL) 2 mm to 4 mm Exhaust cam maximum lift (EVL) 4 mm to 7 mm Intake valve close angle (IVC) 60° to 120° Exhaust valve close angle (EVC) 330° to 360° Exhaust passage diameter (EPD) 9 mm to 12 mm Constraint functions such as the maximum acceleration of the valve train components and the minimum curvature radius are considered in the optimization procedure. In this paper, an orthogonal design [24] and [25] was used to obtain representative combinations of parameters. The optimal parameter combination for the expected indicators was obtained using grey relation analysis [26] and [27]. A part of the orthogonal array is formed by the parameters in Table 5. Table 5. The orthogonal array of the four-valve CAE's parameters (partial) No. IVL [mm] EVL [mm] IVC [°] EVC [°] EPD [mm] 1 2 4 60 330 9 2 2 4.75 75 337.5 9.75 3 2 5.5 90 345 10.5 4 2 6.25 105 352.5 11.25 5 2 7 120 360 12 6 2.5 4 75 345 11.25 7 2.5 4.75 90 352.5 12 8 2.5 5.5 105 360 9 9 2.5 6.25 120 330 9.75 10 2.5 7 60 337.5 10.5 Details of the optimization procedure are presented in [10]. Through simulating calculations with the virtual prototype, the output power, energy efficiency and temperature are shown in Table 6, and the corresponding combination of parameters is shown in Table 7. Table 6. Optimal performance Indicators Pe [kW] n [%] 0 [K] range 4 to 6 20 to 40 210 to 250 value 5.95 23.1 249.9 Table 7. Optimal combination of the parameters Parameter IVL [mm] EVL [mm] IVC [°] EVC [°] EPD [mm] Value 3.5 4 105 337.5 12 6 CONCLUSIONS In this paper, energy loss during the working process was obtained by using air power. To improve output power and energy efficiency, a four-valve CAE was introduced. The prototype was designed to experiment and to compare the performance of the output power and torque of two-valve CAE and four-valve CAE; the experiment was carried out under same working conditions. The eight factors related to the performance of the multi-valves were discussed, and some key parameters were optimized by using orthogonal design and grey relation analysis. The conclusions are summarized as follows: 1. The output torque and output power of the four-valve CAE are higher than of the two-valve CAE. When the intake pressure is set at 0.6 MPa, and the rotation speed is 450 rpm, the output torque of the four-valve CAE is about 2.7 times that of the two-valve CAE. When the intake pressure is set at 0.7 MPa, and the rotation speed is 450 rpm, the output torque of the four-valves of CAE is about 2.2 times that of the two-valve CAE. Moreover, the output torque and power differences of the two kinds of CAE ascend sharply with increasing rotational speed. 2. The energy efficiency and output power are impacted significantly by the intake valve close angle and exhaust passage diameter. When the intake valve close angle increases from 70° to 130°, the energy efficiency decreases from 28 % to 20 %. Moreover, in the first instance, the output power increases with an increase in the intake valve close angle; when the intake valve close angle is larger than 110°, the output power begins to decrease. 3. When the inlet pressure is 3 MPa, the rotational speed is 600 rpm, intake valve lift is 3.5 mm, exhaust valve lift is 4mm, intake valve close angle is 105°, exhaust valve close angle is 337.5 and exhaust passage diameter is 12 mm, optimal output power equals 5.95 kW, the working efficiency equals 23.1 %, and the average temperature inside cylinder equals 249.9 K. The optimization of the multi-valves CAE can significantly enhance the output power compared with a two-valve CAE. Moreover, the air power efficiency of the four-valve CAE is slightly higher than that of the two-valve CAE. 7 ACKNOWLEDGMENTS 9 REFERENCES The research work presented in this paper is financially supported by the science and technology innovation Financial Grant (2017CXYD-2) from the Finance Department of Inner Mongolia, Grant 61765012 of the National Natural Science Foundation of China, and Grant NJZZ18139 of the project of science foundation of the educational department of Inner Mongolia. 8 NOMENCLATURE A effective area, [m2] Aw heat transfer area, [m2] D diameter of piston, [m] Cp specific heat at constant pressure, [J/(kg K)] Cv specific heat at constant volume, [J/(kg K)] h cam profile displacement, [m] he heat transfer coefficient, [J/(m2K)] G mass flow, [kg/s] L length of connecting rod, [m] M mass of the piston, rings, pin and small end of the connecting rod, [kg] m mass, [kg] P power, [w] G volume flow, [m3/s] p pressure, [Pa] p air density, [kg/m3] V volume, [m3] v velocity, [m/s] W work done per cycle, [J] P the angle of link and crank, [rad] y piston displacement, [mm] 0 temperature, [K] k specific heat ratio, 9 crank angle, [rad] Y polytropic exponent, n efficiency, [%] ra crank speed, [rad/s] Subscripts a ideal 5 supply u upstream d downstream in inflow out outflow [1] Shen, Y.-T., Hwang, Y.-R. (2009). Design and implementation of air-powered motorcycles. Applied Energy, vol 86, no. 7-8, p. 1105-1110, D0l:10.1016/j.apenergy.2008.06.008. [2] Tesar, D. (2015). Open architecture vehicles of the future. Mechanism and Machine Theory, vol. 89, p. 107-127, D0I:10.1016/j.mechmachtheory.2014.11.007. [3] Hwang, H.T., Varma, A. (2014). Hydrogen storage for fuel cell vehicles. Current Opinion in Chemical Engineering, vol. 5, p. 42-48, D0I:10.1016/j.coche.2014.04.004. [4] Huang, K.D., Tzeng, S.-C. (2005). Development of a hybrid pneumatic-power vehicle. Applied Energy, vol. 80, no. 1, p. 4759, D0I:10.1016/j.apenergy.2004.02.006. [5] Huang, K.D., Tzeng, S.-C., Ma, W.-P., Chang, W.-C. (2005). Hybrid pneumatic-power system which recycles exhaust gas of an internal-combustion engine. Applied Energy, vol. 82, no. 2, p. 117-132, D0I:10.1016/j.apenergy.2004.10.006. [6] Chan, C.C. (2007). The state of the art of electric, hybrid, and fuel cell vehicles. Proceedings of the IEEE, vol. 95, no. 4, p. 704-718, D0I:10.1109/JPR0C.2007.892489. [7] Offer, G.J., Howey, D., Contestable, M., Clague, R., Brandon, N.P. (2010). Comparative analysis of battery electric, hydrogen fuel cell and hybrid vehicles in a future sustainable road transport system. Energy Policy, vol. 38, no. 1, p. 24-29, D0I:10.1016/j.enpol.2009.08.040. [8] Singh, B.R., Singh, O. (2008). Development of a vaned-type novel air turbine. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 222, no. 12, p. 2419-2426, D0I:10.1243/09544062JMES993. [9] Liu, C.-M., You, J.-J.,, Sung, C.-K., Huang, C.-Y. (2015). Modified intake and exhaust system for piston-type compressed air engines. Energy, vol. 90, part 1, p. 516-524, D0I:10.1016/j. energy.2015.07.085. [10] Yu, Q., Cai M., Shi, Y., Xu, Q. (2015). Optimization study on a single-cylinder compressed air engine. Chinese Journal of Mechanical Engineering, vol. 28, no. 6, p. 1285-1292, D0I:10.3901/CJME.2015.0520.072. [11] Zhang, X., Xu, Y., Jian Xu, J,, Xue, H., Chen, H. (2016). Study of a single-valve reciprocating expander. Journal of Energy Institute, vol. 39, no. 3, p. 400-413, D0I:10.1016/j. joei.2015.02.013. [12] Xu, Q., Cai, M., Shi, Y. (2014). Dynamic heat transfer model for temperature drop analysis and heat exchange system design of the air-powered engine system. Energy, vol. 68, p. 877-885, DOI:10.1016/j.energy. 2014.02.102. [13] Yu, Q., Cai, M., Shi, Y., Fan, Z. (2014). Optimization of the energy efficiency of a piston compressed air engine. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 6, p. 395-406, D0I:10.5545/sv-jme.2013.1383. [14] Xu Q., Shi Y., Yu, Q., Cai, M. (2014). Virtual prototype modeling and performance analysis of the air-powered engine. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 228, no. 14, p. 2642-2651, D0I:10.1177/0954406214520818. [15] Cai, M., Kawashima, K., Kagawa, T. (2006). Power Assessment of flowing compressed air. Transactions of the ASME, Journal of Fluids Engineering, vol. 128, no. 2, p. 402-405, DOI:10.1115/1.2170129. [16] Yu, Q., Cai, M., Shi, Y., Yuan, C. (2015). Dlmenslonless study on efficiency and speed characteristic of a compressed air engine. Journal of Energy Resources Technology, vol. 134, no. 4, p. 044501, DOI:10.1115/1.4029867. [17] Pipan, M., Herakovič, N. (2016). Volume flow characterization of PWM-controlled fast-switching pneumatic valves. Strojniški vestnik-Journal of Mechanical Engineering, vol. 62, no. 9, p. 543-550, DOI:10.5545/sv-jme.2016.3531. [18] Nagaya, K. (1989). Nonlinear transient response of a high speed driven value system and stresses in valve spring for internal combustion engines. Transaction of ASME, Journal of Vibration, Acoustics, Stress, and Reliability in Design, vol. 111, no. 3, p. 264-271, DOI:10.1115/1.3269851. [19] Shi, Y., Cai, M. (2011). Working characteristics of two kinds of air-driven boosters. Energy Conversion and Management, vol. 52, no. 12, p. 3399-3407, DOI:10.1016/j. enconman.2011.07.008. [20] Zweiri, Y.H., Whidborne, J.F., Seneviratne, L.D. (1999). Dynamic simulation of a single-cylinder diesel engine including dynamometer modelling and friction. Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, vol. 213, no. 4, p. 391-402, DOI:10.1243/0954407991526955. [21] Seeton, C.J. (2006). Viscosity-temperature correction for liquids. Tribology Letters, vol. 22, no. 1, p. 67-78, DOI:10.1007/ s11249-006-9071-2. [22] Cui, Y., Song, Y., Dai, Z., Deng, K. (2010). Power oriented cam profile optimization for gasoline engine. Chinese Internal Combustion Engine Engineering, vol. 31, no. 4, p. 21-24. (in Chinese) [23] Shi, Y., Xiao, J.S., Liu C.X., Cui, D.Z., (2002). Sequential unconstrained minimization technique and augmented Lagrange method for optimal design of engine inlet valve cam. Journal of Wuhan University of Technology, vol. 26, no. 3, p. 77-80. (in Chinese) [24] Chen, J.L., Au, K.C., Wong, Y.S., Tam, N.F.Y. (2010). Using orthogonal design to determine optimal conditions for biodegradation of phenanthrene in mangrove sediment slurry. Journal of Hazardous Materials, vol. 176, no. 1-3, p. 666-671, DOI:10.1016/j.jhazmat.2009.11.083. [25] Liao, Y., Lian, Z., Feng, J., Yuan, H., Zhao, R. (2018). Effects of multiple factors on water hammer induced by a large flow directional valve. Strojniški vestnik - Journal of Mechanical Engineering, vol. 64, no. 5, p. 329-338, DOI:10.5545/sv-jme.2017.5109. [26] Xia, X.T., Li, J.H. (2011). Influence of vacuum on friction torque of rolling bearing for space applications based on optimum grey relational space. Journal of Grey Systems, vol. 23, no. 4, p. 327-334. [27] Wei, G.W. (2010). GRA method for multiple attribute decision making with incomplete weight information in intuitionistic fuzzy setting. Knowledge-Based Systems, vol. 23, no. 3, p. 243-247, DOI:10.1016/j.knosys.2010.01.003. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)10, 643-651 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2018.5475 Original Scientific Paper Integrating a New Software Tool Used for Tool Path Generation in the Numerical Simulation of Incremental Forming Processes Daniel Nasulea - Gheorghe Oancea* Transilvania University of Brasov, Department of Manufacturing Engineering, Romania In the incremental sheet forming simulation, finite element modelling typically was used to anticipate material behaviour, predict deformation forces, thickness reduction or to identify other information. Because the forming tool motions are difficult to be implemented in finite element method (FEM) software systems, and because of the large number of points that describe the tool path, in order to reduce data preparation time, this paper presents the implementation of a new software tool conceived by the authors in the process of the numerical simulation of incremental sheet forming. The software tool uses a CNC file in G-code format to reveal the interpolation point coordinates of the tool motions and the positioning time in a specific ANSYS format. Usually, a tool path for an incremental sheet forming process consists of thousands of interpolation points and is very difficult and time-consuming to implement manually into a FEM software system. The new software tool solves this issue, at least for ANSYS software, being swift and easy to use. The paper also presents how the software tool is integrated and validated in a case study of an incremental sheet forming process simulation concerning different part configurations. Keywords: incremental sheet forming, FEM, numerical simulation, tool path, software tool Highlights • The paper states the necessity of reducing the FEM model preparation time for ISF processes by reducing the difficulty of describing the forming tool motions for numerical simulation applications. • A new stand-alone software tool is developed to convert a CNC G-code file into a format specific for ANSYS, in order to facilitate swift tool path description in ISF processes. • The new software tool can be easily used by any user without additional software skills. • The software tool is integrated into an ISF numerical simulation process. • A case study proves the utility of the developed software tool. 0 INTRODUCTION In recent years, flexible manufacturing processes have gained increased attention from researchers and producers. Incremental sheet forming (ISF) is a "die-less" and flexible manufacturing process from the cold-pressing industry, used for customized products, rapid prototyping, and small batch production series of sheet metal parts [1] and [2]. The ISF working principle presented in Fig. 1 entails the sheet metal blank mounted on a fixing device consisting of a blank holder with or without a backing plate to support the part flange, and a clamping plate which is placed over the blank, fixing it most often by the tightening of several screws. The pack of plates is lifted to a greater height than the part depth, using various configurations of auxiliary components to ensure enough space below the blank to enable the part forming. Plastic forming is performed locally, usually by a simple shape tool with a hemispherical head [3], which is rotated by its axis, following a numerically controlled path. The part is divided into a number of layers from top to bottom, and at each layer the tool follows the CNC path, which is, in fact, the outer contour of the part. The distance between two consecutive layers is called an incremental step down (Az). When one layer is completed, the tool descends by an incremental step and follows the part contour again. This operational stage is repeated until the final part is obtained [2] and [3]. The ISF process undergoes final shape Fig. 1. ISF working principle [2], adapted from [3] *Corr. Author's Address: Transilvania University of Brasov, Department of Manufacturing Engineering, Mihai Viteazul 5, Brasov, Romania, gh.oancea@unitbv.ro 643 continuous development, and for its improvement, in addition to experimental research, more modern techniques for the design of experiments (DOE) or finite element methods (FEM) are used [1] and [4]. FEM analysis has increased in confidence and can provide a first check on what works in the real world and what does not [5]. Various topics are simulated, starting with the most common mechanical issues up to thermal analysis, vibration frequency, forming simulations [6] to [8], fluid flow or injection moulding of plastic parts. Indeed, recently numerous specialists and users around the world have come to understand and implement FEM analysis. These make FEM suitable to be used and to improve many aspects: save time and money in a product lifecycle, simulation and validation of parts behaviour, force prediction, geometric deviations or parts stiffness. 1 FEM IMPLEMENTATION IN ISF PROCESSES ISF processes have multiple applications in different industries [2] and [3], i.e. instrument and body panels in the aerospace industry, door and hood panels in the automotive industry, low volume or unique sheet metal products, products in medical industry: denture plates, ankle supports, metal helmets, and so on. This is one of the reasons that FEM is frequently used in ISF process simulation. Usually software systems like ABAQUS, LS-DYNA, or ANSYS are used to analyse or simulate the process and to study [2]; strains, forces, thickness reduction [9], the influence of tool parameters and path on deformation behaviour [10], springback amplitude and tool path correction [11], or the wrinkling phenomena of thin-walled parts manufactured with ISF or other processes [12] and [13]. The main steps for implementing FEM analysis for an incremental forming process are [2] and [14] design of the part and fixing device and devising the computer-aided design (CAD) models, generating the tool movements or the tool path, designing finite element assembly in FEM software, defining materials and apply boundary conditions, solving the model, and obtaining the results. One of the main problems in FEM simulation in ISF processes is the considerable CPU time, an issue mentioned by many researchers. This is because of many nonlinear problems [14], because of tool trajectory complexity that consists of thousands of integration points, and because of a fine mesh or an adaptive meshing that has to be devised in the contact area between the tool and the sheet metal blank to obtain results as accurate as possible. Many solutions have been found for reducing this computational time. Robert et al. [15] used the incremental deformation theory of plasticity as a solution to improve the central processing unit (CPU) time in the ABAQUS FEM simulation of a cup made from aluminium sheet of 1.5 mm thickness. Computational time can also be reduced by using explicit instead of an implicit solver [14]. Two other solutions used very often entail applying time scaling or mass scaling to the prepared FEM model, without affecting the results to a significant extent [16]. Benedetti et al. [17] applied mass scaling in SPIF numerical simulation, using ABAQUS explicit code, in order to achieve a user-specified time step size, reducing analysis computational cost and time. Naranjo et al. [18] concluded that a medium size blank mesh used in the FEM simulation of an incremental forming process in ANSYS, is in good agreement with the experimentally obtained values, while at the same time, by using this medium size mesh, CPU time was reduced. A simplified numerical approach was developed by Ben Ayed et al. [19] to simulate the ISF process. This approach called ISF-SAM is a simplified finite element model that succeeds to reduce the computation time by more than 50 % in comparison to the time required when using ABAQUS implicit code, with sufficiently accurate results prediction. A decoupling algorithm that separates the finite element model into two zones, the elastic and elasto-plastic deformation zones is proposed by Sebastiani et al. [20]. To reduce the computation time, during simulation this algorithm concentrates on the small elasto-plastic contact zone between tool and sheet, simplifying in this way the model solving. Bambach [21] applied a FEM simulation onto the asymmetrical incremental sheet forming, using an adaptive remeshing strategy based on a multi-mesh method combined with subcycling. It was concluded that by using these methods the calculation time is reduced by up to 80 %. Zhang et al. [22] proposed a new selective element fission approach (SEF) to improve CPU time, reducing the unnecessary calculations of the area outside the plastic deformation zone. These methods entail different mesh densities for blank areas: a fine mesh in the tool-sheet contact area to improve the accuracy of results, and a coarse mesh outside the contact area to reduce unnecessary calculation. It was concluded that the SEF approach could save up to 74 % in terms of computational time with satisfactory accuracy for ISF processes. Therefore, many researchers aim to improve the computational time for FEM simulations of incremental sheet forming processes; in literature, many other implemented solutions can be found for this issue. A time-consuming phase in FEM implementation in ISF processes is not only the computational time but also the FEM model preparation time. In model preparation, because in ISF processes the tool path is not a simple one, the most difficult and time-consuming phase is to specify the tool movements. The described tool path should be the same as for the manufacturing process. Usually, tool path generation is achieved by computer-aided manufacturing (CAM) software systems, which yield a computer numerical control (CNC) file, where the trajectories are described using G code format. However, FEM software systems usually do not allow introducing the CNC file or the G-codes directly as input, in order to describe the tool movements. Thus, solutions are required to describe the tool trajectories for incremental forming processes in FEM simulation. The input file of the tool trajectory for numerical simulation should be defined in terms of two parameters: tool points coordinates for each X, Y and Z direction, and time position. In simple cases, the tool path point coordinates can be manually processed by means of simple spreadsheets [14], [16] and [23]. For complex shapes, the point coordinates have to be extracted from the CAM programmed tool path or from the CNC G-code file manually or by other complex methods. Another research team [24] used LS-DYNA for the simulation of the incremental forming process, dedicated FEM software, and in order to specify the input of tool movement, started with the CNC program from which, using in a first stage a programming software called Python, an ordered file was extracted containing an array of five columns with the coordinates of tool positions for linear interpolations and two circle centre coordinates for circular interpolations [24]. This ordered file has been processed using Matlab software, in order to generate the final output for LS-DYNA simulation, the file which contains the tool position versus time data. Using more than one software tool to transfer data from a CNC to a FEM format, as used in reference [24] can be difficult, errors can occur, and the users certainly need special knowledge to work with such software systems. In this paper, the authors propose to solve this generation issue with a new software tool that can be integrated into ISF numerical simulation. From the user viewpoint, this software tool is a standalone application and generates an input file for the ANSYS software system, which describes the tool movements in terms of tool points coordinates by all directions X Y Z and the tool movement time in each interpolation point. 2 METHOD DESCRIPTION The FEM model preparation stage for an ISF process is a complex and time-consuming operation, especially because of the enormous working volume required to describe the tool motions in FEM software. This stage can be improved in terms of preparation time reduction. Thus, a software application named tool motion points generator (TMPG) has been developed by the authors and integrated into a numerical simulation of an ISF process. The objective of this software tool is to convert a CNC G-code file obtained from a CAM software system into a file which contains the necessary input data in ANSYS format for ISF simulation. This software tool allows users, with no specialized software knowledge, to perform the quick conversion of a CNC G-code file into a text file that contains the tool motion points coordinates and time position necessary as inputs in ANSYS, for the description of the complete tool trajectories in ISF processes. It reduces the FEM model preparation time with a beneficial impact on the total implementation time of a numerical simulation for an ISF process. The TMPG software tool flow is presented in the block diagram in Fig. 2. Fig. 2. Software tool block diagram 2.1 Software Tool Description Many programming languages and systems allow the development of new customized software tools, such as LISP, Python, C++, Delphi, Matlab, Visual Basic and so on [24] to [26]. The software tool was conceived in the Delphi environment which uses the Object Pascal programming language. The graphical user interface (GUI) of the TMPG application shown in Fig. 3 is divided into four separate functional zones by their role in the software tool. Those areas which are different windows are highlighted with black frames. In the left upper corner area, coordinates of the interpolation points extracted from the CNC G-code file are displayed. Another text window displays the time values for each interpolation tool position. The area where the final output is displayed is the right upper corner. Thus, in this area, four successive arrays are written, one for each parameter, X Y, Z coordinates and position time for all points. The fourth area is used for launching commands, is placed below the display zone and consists of eight buttons and an input box for the feed rate. Each button from the GUI associated to the TMPG software tool has a specific role. When "Extract coordinates from CNC file" is pressed, first a warning message is displayed which warns the user that the CNC input file has to include in the first row with tool motion G-codes, all three coordinates X, Y, Z, otherwise the initial tool position is unknown. Once this aspect has been verified, the CNC file to be processed can be selected from a dialogue box. The X Y, and Z coordinates for each interpolation point will be further displayed in the first zone of the TMPG interface. All coordinates are counted, and their numbers are displayed below each correspondent window. These numbers have to be equal; otherwise, errors occur in the final file. Taking into consideration the extracted points coordinates, the time for each point can be calculated by pressing the "Time calculation" button. This process uses the feed rate specified by the user and the distance between the previous and the current point in order to calculate the specific time for each tool position. The distance between two points A(xa, ya, za) and B(xb, yb, zb) is calculated using Eq. (1) [27]. AB = ylixb - X ) + (( - ya) + ( - za)) (1) The obtained tool movement times are displayed in the "Calculated time" window. The button "Generate ANSYS format" generates and displays in the right upper area all four arrays in a specific format, needed as input into ANSYS for tool motion description. While running the software application allows saving additional files as follows: only the points coordinates, points coordinates and calculated time in intermediary files. The final ANSYS format output is also generated. The "Restart" has to be used for the processing of a new CNC file. Fig. 3. TMPG software tool interface 646 Nasulea, D. - Oancea, G. 2.2 TMPG Input and Output Files Considering that TMPG aims to generate the tool movements for an incremental forming process in a format that can be specified in ANSYS software, the software tool needs an input file. This input is a CNC G-code generated by a CAM software system. Fig. 4 presents a sequence of a G-code structure for a CNC file generated by using CATIA "Advanced Machining" workbench. In this example, the file extension is "NCCode". There can be linear or circular interpolations. For circular interpolations, the trajectory is also described as a sum of very small linear movements, which at the end describe a circular motion. Thus, first, the CNC file has to be generated by using a CAM software system, to ensure an input data for the TMPG software tool. The CNC input file will be processed by using the developed TMPG software tool, the final result being the obtained ANSYS output file. The data that can be accepted by ANSYS are a set of lines or arrays in a specific format for each position point parameter X, Y and Z and also for the tool position time. Fig. 5 shows the text format for each parameter that will be specified as input in the ANSYS software system. OlOOO * INTELLI6ENTMANUFACTORYSOITWAREWWW.IMS-SOITWARE.COM * * IMSPOST VERSION : 7.4R * ' USER VERSION :1 * N1 G49 G54 G17 G80 G40 G90 G23 G94 G01 G98 (TOOL DATA : T1 END MILL D 8} N2T1 M6 N3S500 M3 N4G0X9.29Y-2.15Z5 N5G43 Z.78 HI N6G1 X9.96Y-1.04Z.66 F300. N7X10.55Y.61 Z.51 N618X2.78 Y2.78 N619Z2. N620 M30. Fig. 4. CNC G-code structure a) ParameterX *SET, UUX(1,1,1), XI value *SET. UUX(n, 1,1), Xn value Parameter Z *SET,UUZ(1,1,1), Z1 value *SET, UUZfn, 1,1), Zn value b) d) Parameter Y *SET. UUY(1,1,1), Y1 value *SET. UUY(n, 1,1), Yn value Parameter Time *SET.TTimed, 1,1), Tl value *SET, TTime(n, 1,1) , In value Fig. 5. Text format for ANSYS input; a) for the X-axis; b) for the Y-axis; c) for the Z-axis; d) for the tool position time The symbols used in these files have the following meanings: • UUX, UUY, UUZ and TTime are the names for each set of lines or arrays; • n is the number of interpolation points in the CNC file; • X1..Xn, Y1...Yn, Z1...Zn, T1...Tn - are the values for the X Y, Z coordinates and time for each interpolation point. In ANSYS parametric design language (APDL) there are two methods by that the available commands in ANSYS software can be used. The first method is to use the GUI to access the commands, and the second method is to use the "command prompt" to run the specific commands by writing the codes that are the abbreviation of these commands [28]. For specifying the generated parameters for tool movements, the proper way is to use the "command prompt" to input the generated data by using the developed TMPG software tool. Thus, the TMPG software tool generates four arrays, one for each parameter X Y Z and Time. All these arrays are automatically written by the developed software tool in a single text file which can be introduced as input for tool movements, or simply copy the text file to ANSYS "command prompt". 3 INTEGRATING AND VALIDATING THE NEW SOFTWARE TOOL IN AN ISF PROCESS -A CASE STUDY The objective of the case study is to integrate the new software tool into an ISF process, to verify if various dimensional configurations of a frustum of the cone can be manufactured without fracture and, at the same time, to predict the values of the developed forces. These force values will be used to determine whether a Victor Vcentre-55 CNC milling machine can be used for experimental research in incremental sheet forming. In incremental forming, the forming tool applies force on the sheet blank, thus producing a local plastic deformation of the contact surface between tool and sheet. The three components of the force are the radial force Fx, the tangential force Fy and the axial force Fz that typically is the largest [2] and [29]. Further determined in this case study whether this milling machine is capable of applying enough force in X, Y and Z directions, for an incremental deformation of the experimental sheet metal parts manufactured from DC05 deep drawing steel of 1mm thickness. To check if the part configurations can be manufactured and to estimate the values of the forces, a numerical simulation in ANSYS/LS-DYNA FEM software was performed, using the TMPG software tool to obtain, c) in a short time, the time tool movements and points coordinates. For a better overview of the process behaviour and the force variation, six numerical simulations were performed for different geometrical configurations of the experimental part presented in Fig. 6. Table 1. Axial and radial forces for each part configuration Fig. 6. Experimental frustum of a cone Different parts variants with different wall angles $ (40°, 50° and 60°) and height values (35 mm, 45 mm and 50 mm) are analysed, manufactured from deep drawing sheet metal blank DC05 of 1 mm thickness. For each variant, two step depth Az were used: 0.5 mm and 1 mm. The FEM models for simulation were prepared in ANSYS APDL and solved using an LS-DYNA explicit solver [30] and [31]. Fig. 7 presents the main steps for the implementation of FEM simulation in the incremental forming process. The 3D model of the part was designed in CATIA V5 CAD workbench, and then the tool trajectories were generated also in CATIA V5 Advanced machining workbench, from where the CNC file was obtained. Then, the CNC file was processed by using the developed TMPG software tool to obtain the needed ANSYS file. Thus in a few minutes, the application generated all four arrays with a minimum effort, which serve as input for ANSYS to describe the tool motion. Furthermore, in ANSYS, the materials were defined, the boundary conditions were applied, and an output K type file "Ansys_output.k" was obtained. This K file serves as input for the LS-DYNA explicit solver that will run and solve the simulation [30] and [31]. As an example of the obtained final result, Fig. 8 presents the plots for the frustum of a cone with a wall angle $ = 60° and H = 50 mm formed by using a step depth Az = 1 mm. The results for all the six studied parts were summarized in Table 1 that features the maximum values for the axial and radial forces. Because the frustum of a cone is an axis-symmetric part, the maximum force in the X direction is equal to the maximum force value in the Y direction. For this reason, the radial force is measured only in one direction, in this case in the X direction. Only the largest values of the forces from all six simulations results are necessary to check if the selected CNC Maximum Maximum Part configuration axial force radial force [N] [N] $ = 40°, H= 35 mm, Az = 0.5 mm 2062.2 711.17 $ = 40°, H= 35 mm, Az = 1 mm 2189.3 786.44 $ = 50°, H= 45 mm, Az = 0.5 mm 2516.0 1104.8 $ = 50°, H= 45 mm, Az = 1 mm 2760.7 1196.7 $ = 60°, H= 50 mm, Az = 0.5 mm 2650.2 1443.2 $ = 60°, H= 50 mm, Az = 1 mm 2931.9 1624.1 Fig. 7. Main steps for the implementation of FEM simulation a) b) J 1 Fig. 8. Obtained results for $ = 60°, H = 50 mm and Az = 1 mm part configuration; a) resulted FEM model; b) axial force plot; c) radial force plot Fig. 9. Setup used for ISF process implementation; a) Victor Vcentre-55 CNC milling machine; b) sheet metal fixing device; c) elements for ISF implementation, adapted from [33] Fig. 10. Experimental parts; a) FEM model and manufactured part for $ = 40°, H = 35 mm, Az = 0.5 mm configuration; b) FEM model and manufactured part for $ = 50°, H = 45 mm, Az = 0.5 mm configuration; c) FEM model and manufactured part for $ = 60°, H = 50 mm, Az = 1 mm configuration milling machine can be used for the ISF process. From Table 1, it follows that in order to be deformed the part configuration with wall angle 60", height 50 mm and step depth 1 mm requires the largest forces both axial and radial. These values are 2931.9 N for the axial force and 1624.1 N for the radial force. The CNC milling machine Victor Vcentre-55 is capable of developing, based on its X, Y, Z axis ball screws transmission, over 22000 N force along the X and Y axes and over 33500 N force along the Z axis [32]. To avoid damaging the CNC machine mechanisms over time, the manufacturer recommends that in ISF processes, the developed forces do not exceed half of the above-specified values. In this case, the maximum radial and axial forces resulted from FEM simulation are smaller than the recommendation of the milling machine manufacturer. Fig. 9 shows the CNC milling machine, the sheet metal blank fixing device and a short overview of the setup used to implement the ISF process [33]. With regard to parts forming, the numerical simulation prediction that all the part configurations can be manufactured without fracture was validated by machining the parts. Fig. 10 presents both the numerically simulated and the manufactured parts for three configurations using a step depth of 0.5 mm. The conclusion of the case study is that the new software tool is successfully integrated into an ISF process, and the selected part configurations were manufactured without fracture in accordance with the numerical simulation prediction. The considered CNC milling machine can be used safely for machining the presented parts within future experimental research in the field of incremental sheet metal forming. 4 CONCLUSIONS In the numerical simulation of incremental forming processes, the preparation of the FEM models is a time-consuming stage because the FEM software systems do not allow introducing the CAM programmed tool paths directly as inputs. The paper presents an original software tool named TMPG which ensures an easy way to obtain the tool path points coordinates and their position time that have to be introduced into ANSYS software for the description of the tool motion. This software tool allows any user with no specialized software knowledge, the quick processing of the tool path CNC file, finalized with the generation of a text file that contains the information needed to describe the same tool trajectories in ANSYS software. The TMPG utility was proved in a case study of simulating an ISF process of a different frustum of cone geometrical configurations, highlighting the importance of the software tool in the FEM model preparation process. The special utility of the new software tool developed by the authors consists in the fact that the ANSYS file for tool movement description is obtained in a very short time. By using the TMPG application, the FEM model preparation time is reduced substantially, thus avoiding the manual processing or other complex software tools that require additional knowledge and special skills. 5 REFERENCES [1] Mulay, A., Ben S., Ismail S., Kocanda A. (2017). Experimental investigations into the effects of SPIF forming conditions on surface roughness and formability by design of experiments. Journal of the Brazilian Society of Mechanical Sciences and Engineering, vol. 39, no. 10, p. 3997-4010, D0l:10.1007/ s40430-016-0703-7. [2] Nasulea, D., Oancea, G. (2017). Incremental deformation: A literature review. 8th International Conference on Manufacturing Science and Education, p. 03017, D0I:10.1051/matecconf/201712103017. [3] Jeswiet, J., Geiger, M., Engel, U., Kleiner, M., Schikorra, M., Duflou, J., Neugebauer, R., Bariani, P., Bruschi, S. (2008). Metal forming progress since 2000. CIRP Journal of Manufacturing Science and Technology, vol. 1, no. 1, p. 2-17, D0I:10.1016/j.cirpj.2008.06.005. [4] Sharma, V., Gohil, A., Modi, B. (2017). Experimental investigation of single point incremental forming of aluminum sheet in groove test. Applied Mechanics and Materials, vol. 867, p. 177-183, D0I:10.4028/www.scientific.net/ AMM.867.177. [5] Felippa, C. (2015) Introduction to Finite Element Methods, Appendix O index: The Origins of the Finite Element Method. from: http://www. colorado. edu/engineering/CAS/courses. d/ IFEM.d/Home.html, accessed on 2018-26-04. [6] Pepelnjak, T., Milutinovic, M., Plančak, M., Vilotic, D., Randjelovic, S., Movrin D. (2016). The influence of extrusion ratio on contact stresses and die elastic deformations in the case of cold backward extrusion. Strojniški vestnik -Journal of Mechanical Engineering, vol. 62, no. 1, p. 41-50, D0I:10.5545/sv-jme.2015.3051. [7] Deng, W.J., Zhang, J.Y., Liu, L.W., He, D., Xia, W. (2017). Simulation analysis of a new chips recycling process termed forming extrusion cutting. International Journal of Simulation Modelling, vol. 16, no. 4, p. 694-706, D0I:10.2507/ IJSIMM16(4)C016. [8] Trzepiecinski, T., Lemu, H.G., Fejkiel, R. (2017). Numerical simulation of effect of friction directionality on forming of anisotropic sheets. International Journal of Simulation Modelling, vol. 16, no. 4, p. 590-602, D0I:10.2507/ IJSIMM16(4)3.392. [9] Blaga, A., Oleksik, V. (2013). A study of the influence of the forming strategy on main strains, thickness reduction and forces in single point incremental forming process. Advances in Materials Science and Engineering, art. id. 382635, DOI:10.1155/2013/382635. [10] Gatea, S., Ou, H., McCartney, G. (2016). Review on the influence of process parameters in incremental sheet forming. The International Journal of Advanced Manufacturing Technology, vol. 87, no. 1-4, p. 479-499, DOI:10.1007/s00170-016-8426-6. [11] Fu, Z., Mo, J., Han, F., Gong, P. (2013). Tool path correction algorithm for single-point incremental forming of sheet metal. The International Journal of Advanced Manufacturing Technology, vol. 64, no. 9-12, p. 1239-1248, DOI:10.1007/ s00170-012-4082-7. [12] Lemeš, S., Zaimovic-Uzunovic, N. (2008). Using buckling analysis to predict wrinkling in incremental sheet metal forming. Strojniški vestnik - Journal of Mechanical Engineering, vol. 54, no. 2, p. 115-121. [13] Liu, N., Yang, H., Li, H., Yan, S. (2016). Plastic wrinkling prediction in thin-walled part forming process: A review. Chinese Journal of Aeronautics, vol. 29, no. 1, p. 1-14, DOI:10.1016/j.cja.2015.09.004. [14] Nimbalkar, D.H., Nandedkar, V.M. (2013). Review of incremental forming of sheet metal components. International Journal of Engineering Research and Application, vol. 3, no. 5, p. 39-51. [15] Robert, C., Dal Santo, P., Delameziere, A., Potiron, A., Batoz, JL. (2008). On some computational aspects for incremental sheet metal forming simulations. International Journal of Material Forming, vol. 1, p. 1195-1198, DOI:10.1007/s12289-008-0155-4. [16] Senthil, R., Gnanavelbabu, A. (2014). Numerical analysis on formability of AZ61A magnesium alloy by incremental forming. Procedia Engineering, vol. 97, p. 1975-1982, DOI:10.1016/j. proeng.2014.12.352. [17] Benedetti, M., Fontanari, V., Monelli, B., Tassan, M. (2015). Single-point incremental forming of sheet metals: Experimental study and numerical simulation. Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, vol. 231, no. 2, p. 301-312, DOI:10.1177/0954405415612351. [18] Naranjo, J., Miguel, V., Martinez-Martinez, A., Gomez-Lopez, L.M., Manjabacas, M.C., Coello, J. (2015). Analysis and simulation of single point incremental forming by ANSYS. Procedia Engineering, vol. 132, p. 1104-1111, DOI:10.1016/j. proeng.2015.12.602. [19] Ben Ayed, L., Robert, C., Delameziere, A., Nouari, M., Batoz, J.L. (2014). Simplified numerical approach for incremental sheet metal forming process. Engineering Structures, vol. 6263, p. 75-86, DOI:10.1016/j.engstruct.2014.01.033. [20] Sebastiani, G., Brosius, A., Tekkaya, A.E., Homberg, W., Kleiner, M. (2007). Decoupled simulation method for incremental sheet metal forming. AIP Conference Proceedings, vol. 908, p. 1501-1506, DOI:10.1063/1.2741021. [21] Bambach, M. (2016). Fast simulation of Incremental sheet metal forming by adaptive remeshing and subcycling. International Journal of Material Forming, vol. 9, no. 3, p. 353360, DOI:10.1007/s12289-014-1204-9. [22] Zhang, M.H., Lu, B., Chen, J., Long, H., Ou, H. (2015). Selective element fission approach for fast FEM simulation of incremental sheet forming based on dual-mesh system. The International Journal of Advanced Manufacturing Technology, vol. 78, no. 5-8, p. 1147-1160, DOI:10.1007/s00170-014-6723-5. [23] Salah, B.M.E., Hrairi, M. (2011). Process simulation and quality evaluation of incremental sheet forming. IIUM Engineering Journal, vol. 12, no. 3, p. 185-196. [24] Suresh, K., Khan, A., Regalla, S.P. (2013). Tool path definition for numerical simulation of single point incremental forming. Procedia Engineering, vol. 64, p. 536-545, DOI:10.1016/j. proeng.2013.09.128. [25] Gusel, L., Rudolf, R., Brezocnik, M. (2015). Genetic based approach to predicting the elongation of drawn alloy. International Journal of Simulation Modelling, vol. 14, no. 1, p. 39-47, DOI:10.2507/IJSIMM14(1)4.277. [26] Aguado, S., Velazquez, J., Samper, D., Santolaria, J. (2016). Modelling of computer-assisted machine tool volumetric verification process. International Journal of Simulation Modelling, vol. 15, no. 3, p. 497-510, DOI:10.2507/ IJSIMM15(3)9.353. [27] Distance formula - MATHEMATICS (2018). from https:// www.britannica.com/science/distance-formula, accessed on 2018-12-06. [28] ANSYS Documentation (2018). from http://www.ansys.stuba. sk/html/elem_55/chapter1/ES1-1.htm, accessed on 201826-04. [29] Behera, A. K., De Sousa, R.A., Ingarao, G., Oleksik, V. (2017). Single point incremental forming: An assessment of the progress and technology trends from 2005 to 2015. Journal of Manufacturing Processes, vol. 27, p. 37-62, DOI:10.1016/j. jmapro.2017.03.014. [30] Oleksik, V. (2016). Numerical study about the influence of wall angle about main strains, thickness reduction and forces on single point incremental forming process. ACTA Universitatis Cibiniensis - Technical Series, vol. 68, p. 1-6, DOI:10.1515/ aucts-2016-0001. [31] Oleksik, V. (2005). Theoretical and Experimental Research on Incremental Forming of Indigenous Metal Sheets, PhD thesis. Lucian Blaga University of Sibiu, Sibiu. (in Romanian) [32] Balls screws - technical information (2018) from http://pdf. directindustry.com/pdf/hiwin/ballscrews/14370-67561. html, accessed on 2018-26-04. [33] Nasulea, D., Oancea, G. Design and manufacturing of a fixing device for incremental sheet forming process. 22nd International Conference on Innovative Manufacturing Engineering and Energy, vol. 178, art. no. 2004, DOI:10.1051/ matecconf/201817802004. Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 64, (2018), številka 10 Ljubljana, oktober 2018 ISSN 0039-2480 Izhaja mesečno Razširjeni povzetki (extended abstracts) Job Angel Ledezma Pérez, Edson Roberto De Pieri, Victor Juliano De Negri: Regulacija sile hidravličnih aktuatorjev z dodatno hidravlično podajnostjo SI 83 Simon Marčič, Sebastijan Seme, Julijan Jan Salamunic, Zdravko Praunseis, Stanislav Pehan: Gretje s ploščatimi sončnimi kolektorji v kombinaciji s toplotno črpalko SI 84 Jun Liu, Weilong Liu, Xiaofeng Wang: Raziskava odkrivanja in upočasnitve napredovanja razpok v nelinearnih rotorjih SI 85 Riad Ramadani, Marko Kegl, Jožef Predan, Aleš Belšak, Stanislav Pehan: Vpliv celične strukture v telesu zobnika na vibracije, ki jih povzroča ubiranje SI 86 Fabrício Alves de Almeida, Guilherme Ferreira Gomes, Rachel Campos Sabioni, José Henrique de Freitas Gomes, Vinícius Renó de Paula, Anderson Paulo de Paiva, Sebastiâo Carlos da Costa: Preverjanje sposobnosti merilnega sistema za identifikacijo virov variabilnosti pri strižnem preizkusu uporovnih točkovnih zvarov SI 87 Qihui Yu, Xin Tan, Jianguo Meng, Fei Shang, Xueqing Hao: Raziskava zmogljivosti štiriventilskega motorja na stisnjen zrak SI 88 Daniel Nasulea, Gheorghe Oancea: Integracija novega programskega orodja za ustvarjanje poti orodja v numerične simulacije procesa inkrementalnega preoblikovanja SI 89 Regulacija sile hidravličnih aktuatorjev z dodatno hidravlično podajnostjo Job Angel Ledezma Pérez - Edson Roberto De Pieri - Victor Juliano De Negri* Zvezna univerza v Santa Catarini, Brazilija V mnogih industrijskih aplikacijah obstaja potreba po nadzorovanju sile, npr. pri rokovanju z materialom ali pri obdelavah, kjer je vrh robota v stiku z opremo. Problem regulacije sile je kompleksnejši od problema regulacije položaja in praktična izvedba je pogosto težavnejša, saj je dinamika bremena del zaprtozančnega sistema. Eden od načinov za premagovanje problemov, povezanih s sistemi za regulacijo sile, je s povečanjem podajnosti, ki jo je mogoče aktivno ali pasivno upravljati oz. spreminjati. Pasivni podajni sistemi za nadzorovano silo potrebujejo za odzivanje brez oscilacij fleksibilno povezavo med aktuatorjem in bremenom, npr. z vzmetjo. Ta pristop omogoča shranjevanje energije in se uporablja pri varnostno kritičnih aplikacijah. Stisljivost fluida v hidravličnih aktuatorjih zagotavlja učinek vzmeti, ki ga je mogoče opisati s hidravlično togostjo. Izkoristiti ga je mogoče s kapacitivnimi hidravličnimi komponentami, kot so gibke cevi ali zbiralniki. Da bi se izognili uporabi vzmeti, je v članku predstavljen pravi hidroelastični aktuator (PHEA), ki uporablja za zmanjšanje togosti prenosa samo hidravlične komponente. S takim pristopom je mogoče zagotoviti neposredno povezavo med aktuatorjem in bremenom, tako pa se poenostavi mehanska konstrukcija in prihrani prostor v aksialni smeri. Uveden je analitični model za ocenjevanje ustrezne vrednosti hidravlične togosti in na osnovi tega je podan postopek za izbiro komercialne hidravlične gibke cevi po katalogu. Sposobnost upiranja premikom bremena zahteva čim manjšo hidravlično togost, vrednost pa mora biti vsaj večja ali enaka hidravlični togosti, ki je potrebna za zahtevano funkcijo. Zmogljivost sistema je opredeljena s kompromisom med hidravlično togostjo, proporcionalnim ojačenjem regulatorja in dolžino gibke cevi. Predstavljen je nelinearni model, ki je bil validiran z eksperimentalnimi rezultati s preizkuševališča. Zasnovan je robustni regulator na osnovi kvantitativne teorije povratnih zank, ki zagotavlja izboljšano stabilnost, delovanje in odpornost sistema za regulacijo sile proti motnjam. V primerjavi s sistemi togih cevi se je izkazalo izboljšanje odziva sistema s hidravlično podajnimi komponentami na stopničaste spremembo sile. Zmanjšanje hidravlične togosti z uporabo gibkih cevi z velikim volumskim raztezkom prinese zmanjšanje izhodne impedance in odmiki bremena povzročijo manjše spremembe sile. Ker je sistem zaradi zmanjšanja hidravlične togosti bolj podajen in konzervativen, se zmanjša kompleksnost regulatorja in tako bi bilo mogoče uporabiti tudi druge vrste regulatorjev, vključno s klasičnim PID. Glavni prispevek študije je torej v predstavitvi možnosti za hidravlično dodajanje podajnosti in s tem za regulacijo sil med aktuatorjem in zunanjim bremenom brez vzmeti, ter v določitvi analitičnega postopka za izbiro hidravličnih gibkih cevi, ki so potrebne za doseganje tega cilja. V prihodnjih raziskavah bi bilo treba razviti postopke izbire zbiralnikov s poudarkom na regulaciji sile in analizi hidravličnega blaženja za izboljšanje učinkovitosti regulacije sile pri sistemih s pravim hidroelastičnim aktuatorjem (PHEA). Ključne besede: hidravlična regulacija sile, hidravlična podajnost, gibka cev z velikim volumskim raztezkom, regulacija QFT, stisljivost fluidov, dinamično modeliranje *Naslov avtorja za dopisovanje: Zvezna univerza v Santa Catarini, Oddelek za strojništvo, LASHIP, Campus Universitario, Trindade, 88040-900, Florianopolis, SC, Brazilija, victor.de.negri@ufsc.br SI 83 Gretje s ploščatimi sončnimi kolektorji v kombinaciji s toplotno črpalko Simon Marčič1 - Sebastijan Seme1 - Julijan Jan Salamunic1 - Zdravko Praunseis1 - Stanislav Pehan2* 1 Univerza v Mariboru, Fakulteta za energetiko, Slovenija 2 Univerza v Mariboru, Fakulteta za strojništvo, Slovenija Namen predstavljene raziskave je bil ugotoviti učinkovitost ploščatih sončnih kolektorjev v kombinaciji s toplotno črpalko. Problem, ki nastopi pri gretju stanovanjskih stavb in gretju sanitarne vode je, da ni povsem jasno, kaj sploh narediti s sončnimi kolektorji, ki so bili na objekte postavljeni pred nekaj leti, torej pred časom, ko so se pojavile učinkovite toplotne črpalke. Pred lastniki in uporabniki nastopi dilema, ali je sploh smiselno še naprej uporabljati in vzdrževati sistem kolektorjev, pojavi se dilema ali se naj takšen sistem nadgradi s toplotno črpalko kapljevina-kapljevina ali pa kolektorje sneti in uporabiti zgolj toplotno črpalko zrak-kapljevina. Prikazana je osnovna teorija, ki omogoča ocenjevanje sistemov za gretje, predvsem pa so izvedeni eksperimenti ter narejene obširne analize, da bi dobili odgovor na vprašanje o smiselnosti nadaljnje uporabe ploščatih sončnih kolektorjev. Narejeno je bilo priročno namensko preizkuševališče, postavljeno na terenu, ki je sestavljeno iz sončnega kolektorja, toplotnih črpalk in 300-litrskega hranilnika toplotne energije. Preizkuševališče je opremljeno z ustreznimi tipali in sistemom za zajemanje podatkov. Da bi dobili dovolj argumentirane informacije, so meritve in analize podatkov potekale v različnih letnih časih in ob skrajno različnih vremenskih razmerah. Meritve so potekale preko celega dne. Presledek med zajemanji podatkov je bil okoli 10 minut. Narisani so ustrezni diagrami, ki predstavljajo realne situacije, ki se dogajajo v praksi. Na temelju izmerjenih veličin; temperatur, svetilnosti, orientacije sončnih kolektorjev, vremenske situacije in rabljene električne energije, so bili narejeni izračuni toplotnih števil, sledila pa je primerjalna analiza toka energij. Analiza rezultatov meritev je pripeljala do zaključka, da so ploščati sončni kolektorji, ne glede, kako jih vežemo v sistem gretja s toplotno črpalko, kot samostojni grelci, vzporedno ali zaporedno vezani, pravzaprav neuporabni. Sončni kolektorji so učinkovitejši od same toplotne črpalke zrak-tekočina zgolj v primeru izrazito sončnega poletnega dne. Izračuni grelnih števil so v bistvu zavajajoči, vsi dvomi o dominaciji toplotnih črpalk zrak-tekočina pa se razblinijo, ko se izpostavijo energijski tokovi in ko iščemo odgovore o ekonomičnosti gretja. Novosti v članku predstavljajo meritve in analize več različnih modernih sistemov in njihovih kombinacij za gretje stanovanj in sanitarne vode z namenom, da se nedvoumno izpostavi njihova učinkovitost glede rabe energije. Analize so podprte s teoretičnim ozadjem, ki omogoča logično pojasnjevanje zaključkov, ki sledijo iz meritev. Raziskava je obširna, a se ne spušča v tehnične podrobnosti posameznih ogrevalnih sistemov, ker je vsa pozornost osredotočena na primerjalno energijsko učinkovitost posameznih principov delovanja. Ključne besede: ploščati sončni kolektorji, toplotne črpalke, ogrevanje, toplotno število, preizkuševališče, učinkovitost gretja SI 84 *Naslov avtorja za dopisovanje:Univerza v Mariboru, Fakulteta za strojništvo, Smetanova 17, Maribor, Slovenija, stanislav.pehan@um.si Raziskava odkrivanja in upočasnitve napredovanja razpok v nelinearnih rotorjih Jun Liu12* - Weilong Liu1 - Xiaofeng Wang12 1 Državni laboratorij v Tianjinu za konstruiranje in inteligentno upravljanje naprednih mehanskih sistemov, Kitajska 2 Nacionalno demonstracijsko središče za izobraževanje na področju eksperimentalne mehanike in elektrotehnike, Kitajska Pojav in napredovanje razpok v nelinearnih rotorjih motorjev lahko privede do hudih nesreč, zato je zelo pomembno, da je z analizo vibracij mogoče zaznati razpoko v sistemu rotorja ter upočasniti njeno napredovanje. Signali vibracij razpokanega rotorja v prehodnem odgovoru so neperiodični in nelinearni. Značilnosti razpoke se zlahka skrijejo med drugimi komponentami vibracij in jih je težko izločiti iz prehodnih signalov. Zato je za izločanje informacij o razpokah v rotorjih predlagana Hilbert-Huangova metoda marginalnega spektra (HHMS). Za postavitev matematičnega modela in reševanje diferencialnih enačb z numerično integracijo so bili uporabljeni teorija dinamike rotorjev, nelinearna dinamika, enačbe gibanja sistema rotorja z razpoko ter model stopnje odpiranja in zapiranja razpoke. Predlagana je uporaba metode HHMS za izločitev značilnosti napake iz signala vibracij sistema nelinearnega rotorja. Predstavljena je tudi primerjava uporabljene metode in valčne metode. V raziskavi je bil uporabljen pristop z izločitvijo značilnosti razpoke iz prehodnega odgovora na osnovi teorije dinamike rotorjev. Za izločitev značilnosti napake iz signala vibracij nelinearnega rotorja je uporabljena metoda HHMS, upočasnitev napredovanja razpoke pa je uresničljiva s prilagoditvijo raznih parametrov rotorja. Uspešnost metode je bila sistematično preverjena s simulacijami in eksperimenti. Raziskava je bila omejena na prečne razpoke in na upočasnitev napredovanja razpok v sistemu rotorja. V raziskavi je predlagana uporaba metode HHMS na prehodnem odgovoru za odkrivanje razpok v nelinearnih rotorjih. Rezultati simulacije so pokazali, da je predlagani pristop primeren za ocenjevanje prisotnosti zgodnjih razpok v sistemu rotorja. HHMS poleg tega omogoča nedvoumno ugotovitev obstoja razpok z globino, ki dosega približno 5 odstotkov premera rotorja. V članku je predstavljena nova metoda za opisovanje stopnje odpiranja in zapiranja razpok, obravnavana pa je tudi upočasnitev napredovanja razpok z različnimi parametri rotorja. Sledijo rezultati raziskave. 1. Rezultati eksperimentov in simulacij potrjujejo učinkovitost predlagane metode HHMS za izločanje značilnosti razpok iz prehodnih odgovorov. HHMS s svojim integracijskim učinkom točno odraža dejanske frekvenčne komponente v prehodnem signalu. 2. Med vzletom ali pristankom letala, izvajanjem akrobatskih prvin, vzpenjanjem in strmoglavim letom je mogoče napredovanje razpoke upočasniti s povečanjem vrtilnega pospeška in pojemka. 3. Z ustrezno nastavitvijo kota p in podpore rotorja je mogoče upočasniti napredovanje razpoke, npr. z zmanjšanjem simetričnega nelinearnega člena p(0) in s povečanjem vrednosti asimetričnega nelinearnega člena ec(1) in es(1). Metoda HHMS se je v raziskavi izkazala za zanesljivo pri odkrivanju razpok v nelinearnih rotorjih na osnovi prehodnega odgovora. Trenutno je na voljo več metod za identifikacijo in diagnosticiranje razpok v sistemu rotorja, le redke pa so primerne za odkrivanje razpok v nelinearnih rotorjih na osnovi prehodnega odgovora in lahko izkoristijo mehanizme razpokanega rotorja za upočasnitev napredovanja razpoke. Sistem rotorja mora biti za možnost zanesljivega odkrivanja razpok razmeroma dobro uravnotežen. Z metodo HHMS je mogoče učinkovito izločiti signale razpok iz kompleksnih signalov vibracij. Metoda HHMS je v primerjavi z valčno metodo občutljivejša in točnejša metoda za izločanje značilnosti iz prehodnih signalov vibracij sistemov rotorjev z razpoko. Eksperimenti so potrdili rezultate teoretičnih simulacij. Metoda je primerna tudi za industrijsko rabo. Ključne besede: razpokan rotor, prehodni odgovor, Hilbert-Huangov marginalni spekter, valčna metoda, odkrivanje razpok, upočasnitev napredovanja razpok *Naslov avtorja za dopisovanje: Državni laboratorij v Tianjinu za konstruiranje in inteligentno upravljanje naprednih mehanskih sistemov, Kitajska, liujunjp@tjut.edu.cn SI 85 Vpliv celične strukture v telesu zobnika na vibracije, ki jih povzroča ubiranje Riad Ramadani1* - Marko Kegl2 - Jožef Predan2 - Aleš Belšak2 - Stanislav Pehan2 1 Univerza v Prištini, Fakulteta za strojništvo, Kosovo 2 Univerza v Mariboru, Fakulteta za strojništvo, Slovenija Članek obravnava metodo za omejevanje vibracij in hrupa zobnikov. Problem, ki nastopa pri vsakovrstnih zobniških prenosnikih, so vibracije, ki nastopajo zaradi samega principa prenosa moči med zobatima kolesoma. Ker se zaradi ubiranja zob spreminja vzmetna konstanta torzijskega sistema, je ta sprememba že sam po sebi razlog za vibracije in hrup, pa četudi se prenaša povsem konstanten moment. Glavna metoda, ki zmanjša prenos vibracij na okolje, je spremenjena notranja konstrukcija zobniškega telesa. Namesto polnega zobniškega telesa je to telo sestavljeno iz celic, nekakšne palične strukture. Topologija celične strukture je optimizirana tako, da se bistveno zmanjša masa celotnega zobnika, napetosti v strukturi pa se dvignejo, a še vedno tudi na mestih prehodov ostanejo daleč pod nivojem kritičnih. Takšna struktura postane v primerjavi s polnim telesom za velikostni razred bolj elastična. Koncentrirane lokalne obremenitve, ki nastanejo na zobnih bokih, se najprej razpršijo po telesu zobnika, nato pa skozi strukturo potujejo po zapletenih poteh, kjer nastopa več raznovrstnih deformacij, pri čemer se energija vibracij zmanjšuje, kar vodi do njihovega dušenja. Postavljena je bila hipoteza, da bi k zmanjšanju vibracij pripomogla polimerna matrica, ki bi jo naknadno vtisnili v odprto celično strukturo zobniškega telesa. Hipoteze so bile preverjene tako, da so se s 3D-laserskimi printerji izdelali vzorčni zobniki s celično strukturo, nato pa so se njihove lastnosti primerjale z referenčnimi zobniki na preizkuševališču. Preizkuševališče je bilo narejeno tako, da so zobniki ubirali v zaprti zanki, ker se je na ta način lahko zelo udobno in natančno kontroliralo obremenitve, ki so jih prenašali. Zaradi relativne majhnosti, enostavnosti in kompaktnosti preizkuševališča je bilo možno v laboratorijskih razmerah realizirati obratovanje referenčnih in zobnikov s celično strukturo telesa v ustreznem ohišju, kjer so bili kontrolirani pogoji obratovanja tako glede mazanja, vrtilne frekvence, momenta kot tudi glede temperature. Na zobnike, ki so bili predmet opazovanja, so bili nameščeni mini senzorji oziroma merilni lističi, mini mikrofoni ter miniaturna tipala pospeškov, in sicer enkrat tik ob izvoru vibracij, na zobnem vencu zelo blizu zobnega korena, ter drugič na pestu zobnika oziroma v njegovi neposredni bližini. Med obema setoma senzorjev je bilo le zobniško telo, zato se je na ta način lahko meril vpliv strukture telesa na prenos vibracij. Ohišje testnega gonila je bilo opremljeno s prosojno steno, tako da je bilo mogoče stanje zobnikov tudi med preizkušanjem spremljati s pogledom. Tipala so bila povezana z mikro-drsniki z ustrezno merilno opremo in analizatorji signalov, tako da so bili narisani karakteristični diagrami, ki so vsebovali elementarne podatke o tem, kako se s časom spreminjajo amplitude in frekvence nihanj, pospeški in hrup, ki se prenaša preko okoliškega zraka. Analiza podatkov, ki so bili dobljeni z eksperimenti, je potrdila hipotezo, da celična struktura zobniškega telesa odločilno vpliva na zmanjšanje nihanj, ki so posledica delovanja zobnikov. Analiza je pokazala, da pa je treba celično strukturo skrbno optimizirati, saj se s takšnim pristopom izognemo lokalnim koncentracijam napetosti, ki bi lahko bistveno poslabšale trdnostno strukturo celotnega sistema. Posebno dodano vrednost v obliki močnega dušenja vibracij predstavlja vtisnjena polimerna matrica. Predlagan pristop za zmanjšanje vibracij zobnikov je novost v svetu po teoretični, znanstveni in praktični plati tako zaradi uporabljene optimizacije topologije celične strukture kakor tudi zaradi integracije polimerne matrice. Ključne besede: vibracije, zobniki, celična struktura, optimizacija topologije, preizkuševališče, analiza signalov SI 86 *Naslov avtorja za dopisovanje: Univerza v Prištini, Fakulteta za strojništvo, Bregu i diellit, 10000 Prishtina, Kosovo, riad.ramadani@uni-pr.edu Preverjanje sposobnosti merilnega sistema za identifikacijo virov variabilnosti pri strižnem preizkusu uporovnih točkovnih zvarov Fabricio Alves de Almeida1* - Guilherme Ferreira Gomes2 - Rachel Campos Sabioni3 -José Henrique de Freitas Gomes1 - Vinicius Reno de Paula1 - Anderson Paulo de Paiva1 - Sebastiâo Carlos da Costa1 1 Zvezna univerza v Itajubi, Inštitut za industrijski inženiring in upravljanje, Brazilija 2 Zvezna univerza v Itajubi, Strojniški inštitut, Brazilija ' Univerza Sorbona, Tehniška univerza v Compiègnu, Oddelek za strojništvo, Francija V raziskavi je bila opravljena gnezdena študija ponovljivosti in obnovljivosti (NGR&R) merilnega sistema z analizo variance (ANOVA) za identifikacijo komponent variabilnosti pri strižnem preizkusu z delovanjem natezne sile, s katerim se ugotavljata dve lastnosti: strižna in natezna trdnost. Eksperimenti so bili opravljeni z vroče pocinkanim jeklom po načrtu DOE, v katerem izbrani preizkušanci predstavljajo realno amplitudo procesa. Uporovno točkovno varjenje (UTV) je tehnika spajanja konstrukcijskih elementov, ki je razširjena v avtomobilski industriji zaradi svoje prilagodljivosti, primernosti za avtomatizacijo, preprostega rokovanja, raznovrstnih možnosti in nizkih stroškov. Zaradi razširjenosti in pomena tehnologije UTV v industriji so bile razvite nove metodologije za prilagajanje parametrov, ki prispevajo k nadzoru in sposobnosti procesa. Preizkušanec je v strižnem preizkusu izpostavljen delovanju sil v nasprotnih smereh. Težnje po boljši kakovosti vodijo industrijska podjetja k nenehnemu izboljševanju učinkovitosti. Tudi kadar je vsa pozornost posvečena izboljšavam procesa, ta zato ni nujno boljši, saj so lahko viri variabilnosti tudi v merilnem sistemu. Obstaja torej potreba po verifikaciji variabilnosti merilnih sistemov v industrijskih procesih, kot je UTV Za nadzorovanje in spremljanje kakovosti rezultatov UTV je na voljo več preizkusov za določanje mehanskih lastnosti, kot je strižni preizkus. Gre za porušno preiskavo, s katero se vrednoti mehanska trdnost zvarne točke. Porušne preiskave v avtomobilski industriji se izvajajo z občasnim vzorčenjem. Kakovost procesov se verificira s kvantitativnimi metodami in variabilnost v rezultatih porušnih preiskav lahko izhaja tako iz proizvodnega procesa kakor tudi iz samega merilnega sistema. V tem primeru se je treba pri eksperimentalnih postopkih izogniti merilnim napakam. Ponovljivost je določena z variabilnostjo sistema pri vnaprej določenih pogojih merjenja (del, okolje, operater, merilna oprema) enega samega preizkušanca. Obnovljivost je določena s povprečno variabilnostjo merilnega sistema med različnimi operaterji, ki merijo isti preizkušanec z isto opremo. Predstavljena analiza merilnega sistema s primerjavo strojev za natezni preizkus (stroja 1 in 2), na katerih se izvajajo porušne preiskave (denimo strižni preizkus) rezultatov procesa UTV, je izviren prispevek, saj objavljena poročila o tovrstnih analizah procesa UTV ne vključujejo podobnih vrednotenj. Poleg tega je na osnovi pregledane literature mogoče potrditi tudi pomen uporabe pocinkanega jekla v tem procesu. Glede na pomen merilnih sistemov v industrijskih procesih (še posebej pri UTV) je v članku predstavljena študija NGR&R merilnega sistema za strižni preizkus, s katero so bili preverjeni viri variabilnosti merilnega sistema v procesu UTV, in sicer s primerjavo dveh strojev za natezni preizkus po metodi ANOVA. Omejitve procesnih parametrov so bile za zagotavljanje želenih načinov odpovedi opredeljene na podlagi preliminarnih preizkusov. Uporabljena je bila metoda načrtovanja eksperimentov (DOE), tako da so lastnosti uporabljenih preizkušancev predstavljale realno amplitudo procesa UTV. Ugotovljeno je bilo, da daje stroj št. 1 večji prispevek k variabilnosti sistema in da njegovi merilni rezultati niso pod nadzorom, kakor tudi da ima ta stroj manjšo stopnjo ponovljivosti kot stroj št. 2. Poleg tega stroj št. 1 ni bil dobro nastavljen, saj je več preizkušancev med testi zdrsnilo. Iz rezultatov sledi sklep, da so potrebne določene izboljšave pri vpenjanju preizkušancev s prevlekami (kot je pocinkana pločevina), s katerimi bodo odpravljeni zdrsi, rezultati preizkusov s tovrstno opremo pa bodo zanesljivejši. Ugotovljena je bila tudi večja variabilnost meritev pri preizkušancu št. 7, še posebej na stroju št. 1. Nova študija po identifikaciji odstopajočega dela je pokazala zmanjšanje variabilnosti sistema za merjenje natezne trdnosti, novi točki zunaj nadzora pa sta bili ugotovljeni pri preizkušancih št. 1 in 2. Ključne besede: točkovno varjenje, analiza merilnega sistema, strižni preizkus, NGR&R, ANOVA, DOE, variabilnost *Naslov avtorja za dopisovanje: Zvezna univerza v Itajubi, Inštitut za ind. inženiring in upravljanje, 1303 BPS Av., Itajuba, Brazilija. fabricio.alvesdealmeida@gmail.com SI 87 Raziskava zmogljivosti štiriventilskega motorja na stisnjen zrak Qihui Yu* - Xin Tan - Jianguo Meng - Fei Shang - Xueqing Hao Znanstvenotehniška univerza Notranje Mongolije, Šola za strojništvo, Kitajska Motorji z notranjim zgorevanjem (MNZ) v transportnih vozilih so eden glavnih virov onesnaževanja zraka. Onesnaževala iz MNZ po nekaterih poročilih predstavljajo do 30 odstotkov vseh onesnaževal v zraku. Z namenom zmanjševanja teh vplivov se porajajo številne alternative za MNZ, med njimi tudi motorji, ki jih poganja energija iz baterij. Zanimivo možnost za transportna vozila predstavlja tehnologija motorjev na stisnjen zrak (MSZ), ki lahko obratujejo praktično brez onesnaževanja. Vozila, ki jih poganja stisnjen zrak, so bila v prvih preizkusih omejena z zmogljivostjo rezervoarjev. Da bi lahko začeli uporabljati MSZ v vlogi glavnih avtomobilskih pogonov, bodo potrebne še določene spremembe na sistemu za distribucijo stisnjenega zraka. V članku so opisani štiriventilski MSZ, ki lahko izkoristijo več energije v delovnem procesu. Problem, ki ga obravnava študija, je izboljšanje energijskega izkoristka, za reševanje pa je bila uporabljena metoda večkriterijske optimizacije. Študija zajema MSZ z energijskim izkoristkom nad 20 % in izhodno močjo nad 5,5 kW. V članku je opisan izkoristek tradicionalnih MSZ. V izračunu moči motorja je upoštevana kinetična energija toka zraka, energija izotermne ekspanzije in moč v izpuhu. Kadar je efektivni prerez majhen ter je razmerje med vhodnim in izhodnim tlakom pod kritično vrednostjo, je treba upoštevati kinetično energijo. Učinkovit ukrep za izboljšanje zmogljivosti pri visokih hitrostih je povečanje efektivnega prereza in v ta namen je predlagan štiriventilski MSZ. Vrtilno gibanje ročične gredi MSZ se prenaša na sesalno in izpušno odmično gred prek krmilnega jermena in jermenic. Za analizo energijskega izkoristka in izhodne moči batnega MSZ je bil razvit matematični model batnega MSZ na osnovi energijske enačbe, zakona o ohranitvi mase, enačbe pretoka, enačbe stanj, enačbe strukture in enačbe gibanja. Zmogljivost štiriventilskega MSZ je bila določena s simulacijami, za eksperimentalno preverjanje teoretičnega modela pa je bil zasnovan prototip štiriventilskega MSZ. Vpliv parametrov je bil analiziran s pomočjo matematičnega modela. Za optimalen energijski izkoristek in moč so bili s pomočjo ortogonalnega načrta in sive relacijske analize optimizirani parametri kot zapiranja sesalnega ventila, premer izpušnega prehoda, največji dvig sesalnega odmikala, največji dvig izpušnega odmikala in kot zapiranja izpušnega ventila. Rezultati so zbrani v nadaljevanju: 1) Izhodni navor in izhodna moč štiriventilskega MSZ sta večja kot pri dvoventilskem MSZ. Ko je sesalni tlak nastavljen na 0,6 MPa in je vrtilna hitrost 450 vrt./min, je izhodni navor štiriventilskega MSZ približno 2,7-krat večji kot pri dvoventilskem MSZ. Ko je sesalni tlak nastavljen na 0,7 MPa in je vrtilna hitrost 450 vrt./ min, je izhodni navor štiriventilskega MSZ približno 2,2-krat večji kot pri dvoventilskem MSZ. Razlika v izhodnem navoru in moči obeh izvedb se strmo povečuje z višanjem vrtilne hitrosti. 2) Na energijski izkoristek in na izhodno moč pomembno vplivata kot zapiranja sesalnega ventila in premer izpušnega prehoda. Ko se kot zapiranja sesalnega ventila poveča s 70° na 130°, se energijski izkoristek zmanjša z 28 % na 20 %. V prvem primeru se izhodna moč povečuje s povečevanjem kota zapiranja sesalnega ventila do 110°, nato pa se začne izhodna moč spet zmanjševati. 3) Pri vhodnem tlaku 3 MPa, vrtilni hitrosti 600 vrt./min, dvigu sesalnega ventila 3,5 mm, dvigu izpušnega ventila 4 mm, kotu zapiranja sesalnega ventila 105°, kotu zapiranja izpušnega ventila 337,5° in premeru izpušnega prehoda 12 mm je optimalna izhodna moč 5,95 kW, delovni izkoristek je 23,1 % in povprečna temperatura v cilindru je 249,9 K. V predstavljenem delu je podan predlog novega MSZ, ki zmanjšuje kinetično energijo zračnega toka. Rezultati raziskave bodo uporabni pri projektiranju batnih MSZ in podajajo metodo za optimizacijo energijskega izkoristka. V prihodnjih raziskavah bo razvit prototip avtomobila s štiriventilskim MSZ. Izdelan bo regulator hitrosti in vozilo bo preizkušeno na cesti. Ključne besede: štiriventilski motor na stisnjen zrak, energijski izkoristek, izhodna moč, prototip, ortogonalni načrt, siva relacijska analiza SI 88 *Naslov avtorja za dopisovanje: Znanstvenotehniška univerza Notranje Mongolije, Šola za strojništvo, Arding Road No. 7, Baotou, Notranja Mongolija, Kitajska, yqhhxq@163.com Integracija novega programskega orodja za ustvarjanje poti orodja v numerične simulacije procesa inkrementalnega preoblikovanja Daniel Nasulea - Gheorghe Oancea* Transilvanska univerza v Brasovu, Oddelek za proizvodno strojništvo, Romunija Znano je, da se v simulacijah inkrementalnega preoblikovanja pločevine za napovedovanje vedenja materiala, deformacijskih sil, zmanjšanja debeline in drugih parametrov običajno uporablja modeliranje z metodo končnih elementov (FEM). Gibanje preoblikovalnega orodja je težko implementirati v programskih paketih za FEM in je povezano z velikim številom točk, ki opisujejo pot orodja. Da bi se skrajšal čas priprave podatkov, je v članku predstavljena uporaba novega programskega orodja pri numeričnih simulacijah inkrementalnega preoblikovanja pločevine. Avtorji so novo programsko orodje TMPG (angl. Tool Motion Points Generator) razvili v okolju Delphi, ki uporablja programski jezik objektni Pascal. Programsko orodje prevede G-kodo za CNC-obdelavo iz programskega sistema CAM v datoteko, ki vsebuje potrebne vhodne podatke za simulacijo ISF v formatu ANSYS. Programsko orodje uporabnikom omogoča hitro pretvorbo datoteke z G-kodo v besedilno datoteko, ki vsebuje koordinate točk gibanja orodja in časovne koordinate za ANSYS, s čimer je mogoče popisati celotno trajektorijo orodja v procesih ISF. Gre za samostojno aplikacijo z grafičnim uporabniškim vmesnikom, ki jo lahko uporablja kdorkoli in ne zahteva dodatnih računalniških znanj. Uporabnost orodja TMPG je bila preizkušena v študiji primera simulacije procesa ISF za različne geometrije prisekanega stožca. Novo orodje lahko odigra pomembno vlogo v procesu priprave modela FEM in njegovi integraciji v proces ISF. Cilj študije primera je bila simulacija procesa ISF in preverjanje, ali je mogoče kakovostno izdelati prisekane stožce v različnih konfiguracijah in obenem napovedati vrednosti nastalih sil. Te so bile uporabljene pri odločanju o primernosti CNC-rezkalnega stroja Victor Vcentre-55 za eksperimentalne raziskave procesa inkrementalnega preoblikovanja pločevine. Analiziranih je bilo šest različnih delov z različnimi vlečnimi koti $ (40°, 50° in 60°) in različnimi višinami (35 mm, 45 mm in 50 mm), izdelanih iz 1 mm debele pločevine za globoki vlek DC05. Iz rezultatov numerične simulacije je sledila napoved, da bi bilo mogoče vse različice dela izdelati brez loma. Napoved je bila v celoti potrjena na CNC-rezkalnem stroju Victor Vcentre-55. Posebna prednost novega programskega orodja je v tem, da lahko zelo hitro pripravi datoteke ANSYS za opis gibanja orodja. Z aplikacijo TMPG se občutno skrajša čas priprave modela FEM in odpade potreba po ročni obdelavi oz. po zahtevnih programskih orodjih, ki zahtevajo dodatna znanja in posebne veščine. Ključne besede: inkrementalno preoblikovanje kovine, FEM, numerična simulacija, pot orodja, programsko orodje *Naslov avtorja za dopisovanje: Transilvanska univerza v Brasovu, Oddelek za proizvodno strojništvo, Mihai Viteazul 5, Brasov, Romunija, gh.oancea@unitbv.ro SI 89 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 vestnih - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, DOI:10.5545/sv-jme.2011.013. Journal titles should not be abbreviated. Note that journal title is set in italics. Books: Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication. [2] Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken. Note that the title of the book is italicized. Chapters in Books: Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages. [3] Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordic, 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] Štefanic, N., Martinčevic-Mikic, S., Tošanovic, 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-ime.eu/mformation-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-ime.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 320 EUR (for articles with maximum of 6 pages), 400 EUR (for articles with maximum of 10 pages), plus 40 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 Job Angel Ledezma Pérez, Edson Roberto De Pieri, Victor Juliano De Negri: Force Control of Hydraulic Actuators using Additional Hydraulic Compliance Simon Marčič, Sebastijan Seme, JuLijan Jan Salamunic, Zdravko Praunseis, Stanislav Pehan: Heating by Flat Plate Collector Combined with a Heat Pump Jun Liu, Weilong Liu, Xiaofeng Wang: Research of Crack Detection and Delay Crack Propagation on a Nonlinear Rotor Riad Ramadani, Marko Kegl, Jožef Predan, Aleš Belšak, Stanislav Pehan: Influence of Cellular Lattice Body Structure on Gear Vibration Induced by Meshing Fabrício Alves de Almeida, Guilherme Ferreira Gomes, Rachel Campos Sabioni, José Henrique de Freitas Gomes, Vinícius Renó de Paula, Anderson Paulo de Paiva, Sebastiao Carlos da Costa: A Gage Study Applied in Shear Test to Identify Variation Causes from a Resistance Spot Welding Measurement System Qihui Yu, Xin Tan, Jianguo Meng, Fei Shang, Xueqing Hao: Performance Characteristics Research on Power System of a Four-Valve Compressed Air Engine Daniel Nasulea, Gheorghe Oancea: Integrating of a New Software Tool Used for Tool Path Generation in the Numerical Simulation of Incremental Forming Process 579 590 601 621 632 643 9770039248001