Strojniški vestnik Journal of Mechanical Engineering no. 6 year 2021 volume 67 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: Demat d.o.o., printed in 250 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 Vice-President of Publishing Council Bojan Dolšak University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: The cover image shows the real object and virtual simulation of a pair of line gears. The teeth of the line gears are designed to be twisted three-dimensional cantilevers to provide conjugated meshing curves, and therefore the line gears can complete precise transmission in limited space through the continuous point contact of these conjugated meshing curves. Courtesy: Guangxi University, College of Mechanical Engineering, Guangxi Key Laboratory of Manufacturing System & Advanced Manufacturing Technology, China ISSN 0039-2480, ISSN 2536-2948 (online) International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, King Fahd U. of Petroleum & Minerals, Saudi Arabia Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, University of Ljubljana, Slovenia Filippo Cianetti, University of Perugia, Italy Janez Diaci, University of Ljubljana, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Igor Emri, University of Ljubljana, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, University of Maribor, 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, University of Ljubljana, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, University of Ljubljana, Slovenia Gorazd Lojen, University of Maribor, Slovenia Darko Lovrec, University of Maribor, Slovenia Thomas Lubben, University of Bremen, Germany George K. Nikas, KADMOS Engineering, UK Tomaž Pepelnjak, University of Ljubljana, Slovenia Vladimir Popovič, University of Belgrade, Serbia Franci Pušavec, University of Ljubljana, Slovenia Mohammad Reza Safaei, Florida International University, USA Marco Sortino, University of Udine, Italy Branko Vasič, University of Belgrade, Serbia Arkady Voloshin, Lehigh University, Bethlehem, USA General information Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue). Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor. Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http:// 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 peer-review process. © 2021 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 https://www.sv-jme.eu. Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6 Contents Contents Strojniški vestnik - Journal of Mechanical Engineering volume 67, (2021), number 6 Ljubljana, June 2021 ISSN 0039-2480 Published monthly Papers Jiang Ding, Aiping Deng, Liwei Liu, Mengen Lu: Dynamic Analysis of Line Gear Pair Based on Numerical Manifold Method 275 Radomir Bokic, Jovan Vladic, Milan Kljajin, Vesna Jovanovic, Goran Markovic, Mirko Karakasic: Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators 287 Ignas Sokolnikas, K^stutis Ciuprinskas, Jolanta Ciuprinskiene: Minimization of the Lifecycle Cost of a Rotary Heat Exchanger Used in Building Ventilation Systems in Cold Climates 302 Anmol Bhatia, Reeta Wattal: Process Parameters Optimization for Maximizing Tensile Strength in Friction Stir-Welded Carbon Steel 311 Boopathi Sampath, Sureshkumar Myilsamy: Experimental Investigation of a Cryogenically Cooled Oxygen-mist Near-dry Wire-cut Electrical Discharge Machining Process 322 Gautham Velayudhan, Prabhu Raja Venugopal, Ebron Shaji Gnanasigamony Thankareathenam, Mohanraj Selvakumar, Thyla Pudukarai Ramaswamy: Reliability-Based Design Optimization of Pump Penetration Shell Accounting for Material and Geometric Non-Linearity 331 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 © 2021 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2021.7118 Original Scientific Paper Received for review: 2021-02-01 Received revised form: 2021-04-28 Accepted for publication: 2021-05-14 Dynamic Analysis of Line Gear Pair Based on Numerical Manifold Method Jiang Ding1-2-* - Aiping Deng1 - Liwei Liu1 - Mengen Lu1 !Guangxi University, College of Mechanical Engineering, China 2Guangxi University, Guangxi Key Laboratory of Manufacturing System & Advanced Manufacturing Technology, China As a kind of tiny gear based on space curve meshing theory, the line gear is very suitable for miniaturized machines due to its compact size and low weight. However, the line gear usually suffers from serious vibration problems since its line teeth are designed as twisted three-dimensional cantilevers to provide conjugated meshing curves. A dynamic model of the line gear pair is established in this paper using the numerical manifold method (NMM) to alleviate its vibration conditions, which can simultaneously provide mathematical and physical covers. The displacement function is first derived for the line teeth, and the dynamic equations of the manifold element are acquired. After inspecting the reasons that cause meshing excitation, the dynamic response of the line teeth is attained in all three orthogonal directions. The attained dynamic response shows that the vibration in the axial gear direction is more significant than that in the curvature direction. Furthermore, the vibration differential equations of the line teeth are solved through a detailed example, and the relationship between the design parameters and the natural frequency is revealed. The vibration characteristics of the first four order of the line gear are revealed through the method of NMM and compared with the result that is carried out through the commercial finite element method (FEM). The comparison shows that NMM can efficiently relieve the vibration problems of the line gear. Keywords: line gear, dynamic response, vibration, numerical manifold method Highlights • The line gear is designed to be applied in miniaturized machines. • The line teeth of the line gear are twisted three-dimensional cantilevers. • A dynamic model of the line gear is established through the numerical manifold method. • The numerical manifold method provides mathematical and physical covers simultaneously. 0 INTRODUCTION As core parts to transmit motion and moment in miniaturized machines, tiny gear-boxes and the mini-gears inside them have raised significant concern. These days, many new gears have been proposed [1] to [3]. Among them, the line gear, as a tiny gear based on space curve meshing theory, has shown great potential in precise transmission in limited space. The line gear has been used in several kinds of small gear transmissions [4], but its line teeth usually suffer from fatal vibration problems and meshing transmission failure during the meshing process [5]. Many traditional gears also face similar vibration problems. Yi Yang et al. [6] investigated the non-linear dynamic response of a spur gear based on periodic mesh stiffness to improve the dynamic characteristics, and Zong Meng et al. [7] studied the vibration response and analysed fault characteristics of gears. To alleviate the vibration conditions, Belingardi et al. [8] made a dynamic analysis of a gear transmission system for an electric vehicle through a multibody approach, and Marco Cirelli et al. [9] presented a novel implementation of a specific multibody model through the tip relief micro-modification on spur gears. In addition, the method of simulating gear pair dynamic response is a potential research direction. Cirelli et al. [10] proposed a refined methodology to simulate the non-linear dynamic response of spur gears through the multibody model based on a penalty contact formulation and considers teeth. Ebrahimi and Eberhard [11] established a rigid-elastic modelling of meshing gear wheels to investigate the effects of multi-tooth contact, as well as backlash and left and right-hand side contact of the meshing teeth. The calculation of gear meshing stiffness is indispensable in the dynamic analysis of gear pairs, and is also one of the research directions for scholars. Cooley et al. [12] studied the calculation methods of gear meshing stiffness, and compared the average slope method and local slope method. Luo et al. [13] proposed a tooth tip modelling method based on defect ratio and independent of gear shape to better calculate the time-varying meshing stiffness of the gear. However, due to the irregular shape of the line teeth, the dynamic study on the traditional gears cannot be applied to alleviate the vibration conditions of the line gear [5]. As shown in Fig. 1, the line teeth can be considered multiple twisted three-dimensional space curved cantilever beams axisymmetrically fixed *Corr. Author's Address: Guangxi University, College of Mechanical Engineering, China, jding@gxu.edu.cn 275 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 around the end of the wheel cylinder. The vibration problems of the line gear are mainly caused by these cantilever beams [5]. For now, many theoretical models of three-dimensional cantilever beams have been studied. Wang and Li [14] calculated the natural frequency of functionally graded cantilever beam and analysed the effects of material gradient parameters, the length-depth ratio of cantilever beam and boundary conditions on the vibration. Zhang et al. [15] investigated the natural frequency and mode of pyramidal Timoshenko beams, and analysed the effects of mechanical properties, rotational inertia and shear deformation on the natural frequency of beams with different taper ratios. Zhao and Wu [16] established the motion coupling equation of a rotating three-dimensional cantilever beam, and showed the effects of Coriolis term and steady-state axial deformation on coupling vibration. It is worth noting that stress and deformation analysis of the cantilever beam is an important research direction. For instance, Zhou et al. [17] investigated the 3D dynamics of a rotary functional gradient cantilever beam and revealed the couplings among the axial, flap-wise, and chordwise deformations. However, the previous models mainly focus on a 3D cantilever of normal shapes, rather than twisted three-dimensional space curved cantilever beams in the line gears. In this paper, we propose the numerical manifold method (NMM) to relieve the dynamic problems of the line gear. The NMM is a numerical calculation method based on the concept of manifold elements, and has been widely used in solving continuous linear elasticity problems [18], crack problems [19], and continuous stress-strain field problems [20]. The research results show that NMM can improve computational accuracy and convergence compared with the finite element method [21]. Wei-bin WEN et al. [22] presented a new NMM based on quarticuniform B-spline interpolation, which has high interpolation accuracy and rapid convergence. Compared to other methods, the NMM can provide mathematical covers and physical covers simultaneously. The mathematical cover is not limited by the physical cover since unknown variables contained in mathematical covers are no longer corresponding to the node displacement. Therefore, the NMM can reduce the sensitivity of unit distortion in the line gears and efficiently obtain numerical results with decent accuracy. In this paper, the dynamic equations of manifold element are established through the NMM, and the vibration differential equation of the cylindrical helical line tooth is derived from the stress analysis of the line tooth micro-element. The effects of the main design 276 Ding, J. - Deng, A.P. - parameters of line gears on their natural frequency are analysed. The dynamic response of the line gear pair and the natural vibration modes of the first four orders are revealed, and the comparison between the natural frequency result computed with the commercial finite element method (FEM) and the NMM is carried out to verify further that the dynamic model established by the NMM. Fig. 1. Space curved cantilever beams on a pair of line gears 1 DYNAMIC MODEL The usual interpolation functions for NMM include polynomial function, trigonometric function, and B-spline function. Polynomials and trigonometric functions are widely used due to their convenient calculation, but they are inaccurate and too sensitive for the mesh distortion when describing the displacement continuity in a pair of line gears. In contrast, the B-spline interpolation function possesses advantages of continuity, local support, local control and modification, and can be employed to solve mechanical problems effectively [23]. Therefore, the B-spline function is adopted as the interpolation function of NMM in this paper for the line teeth, which are sensitive to the unit deformation. 1.1 Displacement Function of Manifold Element The line gear meshing model and the meshing coordinate systems of line gears are shown in Fig. 2. The driving gear is a cylindrical helical line gear. The rotating speed of driving gear is donated as mj, and the rotating speed of driven gear as m2. oJ - xj yj zj and o'J- x\y\z\ are the fixed coordinate system and the rotating coordinate system of the driving gear, respectively. 02 - /2 ¿2 and o'2- x'ly'iz\ are the fixed Liu, L.W. - Lu, M.G. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 coordinate system and the rotating coordinate system of the driven gear, respectively. The angle of the centre axis is (n- 0). The pair of interaction meshing forces between the driving gear and the driven gear are donated as F and F', respectively. a) b) Fig. 2. Simulation diagrams of line gears: a) meshing model; b) meshing coordinate systems The B-spline interpolation function has different kinds of definitions in mathematics. To ensure the coming programming practical, the B-spline basis function in this paper is explicitly defined by recursion and expressed in polynomial form. The B-spline basis function is defined as: . 1 xf < x.,, b,„( x) = ■' ' ¡+1 Ba (x) = „ others x - x Xi+k x< 'i+k+1 -vi+1 B.,k-i (x) , — B+1,k-1 (x ) (1) where [xv, x+i] is defined as the B-spline node interval i', i is the number of B-spline, and k is the order of B-spline function. To mesh the line gear model, the three-dimensional B-spline interpolation function needs to be established. Through Eq. (1), B-spline functions in x, y and z directions can be obtained, respectively. Three-dimensional B-spline interpolation functions can be obtained as follows: k1 k2 k3 N = L L L Bt_k. + niA( x) B n=0 n =0 n =0 =i, (y)Bl-fe + n ,k (z) j-k2 +^2,k2 y* ' i-k^ + n,k (2) where kj, k2 and k3 are the order of B-spline function in x, y and z directions respectively, and i, j and l are the number of B-spline in x, y and z directions respectively. Eq. (2) indicates that the three-dimensional B-spline function satisfies the requirement of weight function; thus, the B-spine function can be taken as the covering weight function of the three-dimensional manifold element. The k order B-spline polynomial function with non-repeated nodes has k - 1 order continuity and coordination in the whole element, which ensures the solution accuracy. The B-spline basis function has high order continuity and coordination, and obtains decent interpolation precision. Each weight function should correspond to the local mathematical cover Ui-k1+n1,;-k2+n2, i-k3+%, and each manifold element corresponds to at least one local cover function. The whole displacement function of manifold element e is: U(e) = Ne x(sD(e)) _ j(e)D(e) (3) where Ne, S, T(e) and D(e) are the interpolation function matrix of the element e, the order matrix of the local cover function, the covering matrix of the element e and the degree of freedom matrix of the element e, respectively. Also, matrix D(e) contains time variable, t. T(e) and D(e) are expressed as: TW = [N1S N2S D W _ = [d- N ,S], dj ] (4) (5) The cover form of element e is shown in Fig. 3. The overlap part of each local cover function and the physical domain is manifold element e. The global cover function on the solution domain can be obtained through the cover function and weight function on each manifold element. X Dynamic Analysis of Line Gear Pair Based on Numerical Manifold Method 277 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 M D(e)+ C D(e) + K D(e) = F , p. p. p. p' Fig. 3. Relation between manifold element and local mathematical cover 1.2 Dynamic Equation of Manifold Element The vibration characteristic of the manifold element e is investigated without loss of generality. The linear elastic strain energy matrix Ve is: V = 2 iiil>' T 2 (6) where oe and £e are the stress matrix and strain matrix of the manifold element e, respectively. Q is the integral domain. oe and £e are expressed as: £e = QU(e) = QT(e)D(e) = Be D(e), Ge = Eee = EQT(e)D(e) = EB eD(e), (7) (8) where Q, E and B are the differential operator matrix, the elastic matrix and the element strain transformation matrix, respectively. The strain energy matrix Ve and kinetic energy matrix We of the manifold element e can be obtained as follows: Ve = I[D(e) J JJJ[be ]TEBedxdydz ■ D1 2 n 3U(e)( x, y, z, t) (9) W -2JJJ dt dxdydz [D(e) Y JJJ[T(e)]T pT(e)dxdydz • D(e). (10) Through the Lagrange equation of the second kind, the differential equation of motion is: _d rsw dt 1 3D 3W dV ^ --+ — = F. 3D 3D (11) The viscous damping theory is adopted to take into account the damping influence, and the dynamic equation of manifold element is: (e) — (12) where Me, Ce and Ke are the element mass matrix, the element damping matrix and the element stiffness matrix, respectively. They are expressed as: ~\T Me =JJJ[T(e)] pT(e)dxdydz n Ce =JJJ[T(e) ] cT(e)dxdydz . n Ke = JJJ[B(e) ]T EB(e) dxdydz (13) To acquire the final global dynamic equation expediently, a transformation matrix is defined as C'(e). The global generalized mass matrix M, the global generalized damping matrix C, and the global generalized stiffness matrix K are expressed as: M = J [C'(e)]r M eC(e} e-1 c = jr [ew]r cec,(e) . e=\ K = J [C'(e)]r K eC{e) (14) The global dynamic equation of manifold element can be obtained as follows: Mi) + CI) + KD = F. (15) 2 ANALYTICAL SOLUTION OF EQUATION 2.1 Initialization and Calculation Method of Manifold Element Equation It is notable that D, D and D in the dynamics equation of manifold element are generalized unknowns and no longer correspond to the displacement, velocity, and acceleration in the solution domain. In the NMM, the nodes of the manifold element only reflect the cover information and irrelevant to the actual physical nodes. As a result, the corresponding initial values D0, D0 and D0 must be derived by transformation. The total number of manifold elements is N'x*N'y* N'z, where N'x, N'y and N'z are the number of discrete interpolation elements in x, y and Z directions, respectively. The final calculation nodes are: 278 Ding, J. - Deng, A.P. - Liu, L.W. - Lu, M.G. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 x1-k 1 - x2-kj - " ' - x0 - "' - xN' +k '+1 yi-k2 - y2-k2 - " - yo- ■■■- y*ry+k-+i. z1-k3 - z2-k3 - " - z0 - " - zN' +k'+1 k nodes are defined out of the physical cover to cover the entire solution domain. These nodes can help to clarify all the weight functions and meanwhile ensure that the interpolation functions near the node satisfy the coordination requirements. The initial displacement and the initial velocity of the given node are Uo = U (xt, yj, z^ and U = U (xt, yj, z^. Do and D0 can be obtained as follows: 2.2 Natural Frequency and Mode of the Line Tooth In Fig. 2b, F represents the meshing load on the line tooth. The meshing load F can be converted into equivalent load f on the constant section microelement of the curve length ds. It is assumed that the torsion deformation and shear deformation of the line tooth micro-element are ignored. Without considering the central axial stress of the line tooth, the microelement is subjected to the equivalent external force f, the shear force q, the bending moment m' and the inertial force pAds• du2(s,t) / dt2. The force analysis is shown in Fig. 4. where T = Do = VUo, D o = VU o. T/l, Ji,L G1, Ji NfxN[xNfT Uo = U((yvzi) " U((,y^y,zN[) (16) (17) (18) (19) In Eq. (18), G/, j, L is the transformation matrix. And through Eq. (15), D 0 is attained as follows: Do = M-1 (F - KD0 - CD„). (20) It is worth noting that the matrices in the dynamic equations remain symmetric and positive, so the time integral method can be used for calculation. Considering the great increase of freedom degrees in the manifold elements, the Newmark method with good numerical dissipation is adopted. According to Eq. (15), the dynamic equation of manifold element at time t+A t can be written as MD t+At + CD t+Ai- ■KDt+At = Ft+At. (21) The basic formulas are: D t+A = D, +[(1 -7' )D t +7D t+Ai ]Ai, (22) D,+Ai = D, + D, At- - -p | D, + pD,+A, At2. (23) When the conditions /> 0.5 and P'> (0.5 + yf / 4 are satisfied, the algorithm is stable with good numerical dissipation unconditionally, and the dynamic response function of the line teeth can be obtained. Fig. 4. Force analysis of the micro-element at any position of the line tooth The force balance equation and the moment balance equation are established as follows: 82u (s, t) 8q (s, t) , , pAds—+ ' ds ■ cos2-f (s,t) = 0 dt2 8s 5m'(s, t) 8s ds -q(s,t)■ ds ,(24) f( A i ds , 82u(s,t)(s,t) ds2 +f (s, t )cos 2--pA--—P—- cos 2-= 0 J V ' 2 H dt1 2 u (s, t ) = G (s )■ H (t )■ h (t) where u(s, t), £(s) and X are the displacement vector function, the mode function, and the helical angle, respectively, and h(t) is the unit vector. The second-order component of the moment equation is omitted, and the motion differential equation of the cylindrical helical line tooth for the undamped free vibration is obtained as follows: pA d2u (s, t) d4u (j, t) 1 EJ ~ dt2 dt4 -cos A = 0. (25) Through the separated variable method, Eqs. (26) and (27) can be attained, d2 H (t) dt2 + m2 H (t) = 0, (26) Dynamic Analysis of Line Gear Pair Based on Numerical Manifold Method 279 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 d4G( s) ds4 ~fi4G (s) = 0, where = _ œ2p EJ cos X The solution of the Eq. (26) is: H(t) = CH sin(fflt + ç), (27) (28) (29) where CH and y are undetermined coefficients determined by initial conditions. Supposing that G(s) is eas, and the solution of the Eq. (27) is: G(s) = Ol sin fis + O2 cos fis +O3 sinh fis + O4 cosh fis. (30) The boundary conditions are that one end of the line tooth is fixed and the other is free. That means: EJ G (s) = 0 d2G (s ) dG (s ds - = 0 (s = 0 ), (31) rd3G (s ) - = 0 EJ-^ = 0 (s = s,\. (32) ds2 ds3 V " Substituting the boundary conditions into Eq. (30), the frequency equation is: cos fist ■ cosh fist +1 = 0. (33) Thus, the natural frequencies of each order are obtained as follows: fi,s, | EJ cos A «"I s H—• (34) where fi1st = 1.875, fist = 0.5(2;- 1), (i = 2, 3, ..., n), and s = fV x\t) + y '(*) + z)dt. (35) 1 i The mode functions corresponding to each order of natural frequency are: G (s) = ^ (sin fiist - sinh fiist) + cos fiist - cosh fis,, (36) where cos fiist + cosh fiist sin fi,s, + sinh fi,s, (i = 1,2,...). (37) 3 VIBRATION ANALYSIS OF THE LINE TEETH 3.1 Meshing Excitation Analysis of Line Gear The dynamic excitation of line gear pair in meshing transmission can be divided into three types: meshing stiffness excitation, error excitation, and meshing impact excitation. ® Meshing stiffness excitation The contact ratio of line gear in the meshing process is a periodic function. The change in the number of meshing gear pairs will cause the change in the gear angular velocity. Due to the inhomogeneity of the velocity, the vibration of the line teeth is generated. When the contact ratio of line gears changes periodically during the meshing process, the meshing stiffness of gear pairs alters in the same period. The product of the meshing stiffness ks of the gear pair and the relative displacement As't of the meshing line teeth in the normal direction is defined as the meshing stiffness excitation force, Fk(t). Fk(t) = ks -As',, (38) ks =\_km + kv ■ sin (wmt + yk)], (39) km = k\k'2f (k \ + k\), (40) kv = km (s-1), (41) As', = 1 -0.5(kn -k2), (42) where km, k'1 and k'2 are the average stiffness of gear pairs, the single-tooth stiffness of the driving line tooth and the driven line tooth, respectively, kR and k2 are the correction coefficients of the centre curve of the driving line gear and the driven line gear, respectively, and e is the coincidence degree of the line gears. © Error excitation The meshing error of the line gears is caused by machining error and installation error. The machining and manufacturing methods of line gears can be divided into three types: profiling, stereo lithograph apparatus (SLA), and selective laser melting (SLM). The comparison of these three manufacturing methods is shown in Table 1. Table 1. Comparison of manufacturing methods of the line gear [24] to [26] Shape Dimensional Surface roughness, accuracy accuracy [mm] Ra [urn] Profiling poor <0.02 5 SLA good <0.02 20 SLM good <0.03 25 As shown in Table 1, error excitations of the line teeth are different, depending on the methods by which they were manufactured. The error excitation is mainly caused by the deviation of tooth shape and the surface roughness of the teeth. The deviation of 280 Ding, J. - Deng, A.P. - Liu, L.W. - Lu, M.G. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 tooth shape comes from the discrepancy between the theoretical meshing curve and the actual meshing curve and leads to the vibration during the meshing. The surface roughness of line teeth can meet the requirements by secondary machining. Therefore, the vibration caused by roughness can be ignored. The error excitation is defined as: F (t) = ks-(1 -kc ), (43) where kc is the correction coefficient of the meshing curve. The approach shock force at any time t in the At can be defined as: F = Md ■AUJ At (46) where Md and AUd are the mass and the velocity variation of the single meshing line tooth. Due to the inertia, the recess shock force will reduce to a negligible level, so only the approach shock in the meshing impact excitation needs to be considered. ® Meshing impact excitation Due to the size error and pitch error of the line gear, the actual meshing point of engaging-in and engaging-out usually deviate from the theoretical meshing point of engaging-in and engaging-out in the meshing process, which disturbs the rotating speed of the meshing line teeth. The vibrations at the beginning and the end of gear meshing are defined as the approach shock and recess shock, respectively. When the line teeth fail to rotate in time into the meshing state at the meshing point of engaging-in due to the pitch deviation, the approach shock will occur. Meanwhile, when the meshing line teeth fail to separate in time at the meshing point of engaging-out due to tooth shape error or pitch deviation, and the latter pair of the meshing line teeth may rotate into a meshing state at the meshing point of engaging-in. In this case, the former pair of line teeth has to change their velocity to maintain their continuous motion transition, which will cause the recess shock. The difference between the theoretical time TO and the actual time t0 at the meshing point of engaging-in is defined as the action time of the approach shock. At = T -t0l T = 20 2n (44) (45) 3.2 Example In the meshing process of the line gear, meshing stiffness excitation and error excitation are the periodic excitation and the meshing impact excitation is the transient excitation. As the dynamic response of the line teeth mainly belongs to the steady-state response, only meshing stiffness excitation and error excitation are needed to be concerned. Supposing the angular velocity mm of the driving gear is n [rad s-1], and the excitation function is: F (t) = Ft (t) + F, (t) = F [1 + (e- l)sin(fflmt + q>k)], (47) F = kk^r [2" 0.5( + k2)"K], (48) s = Açn / (2^) (49) where A^1 and yk are the rotation angle of the single meshing line tooth and the initial phase. The centre curve equations of the driving line tooth are: xCi = m cos tc yCi = m sin tc , zc = ntc + nn -n< t <-0.475n. (50) According to the theory of line gear meshing, the centre curve equations of the driven line tooth can be obtained as follows: n +1„ ■ n + tc y 2c =sm—— ( i- cos6- nt + nn - ( m -cosQ - nt + nn - ( z2c = —m - sin6 — kmd d 2 2k, \ . knd . t + n sin 6 + a + sin —- 2 i. k„ d + d d \ 2k sin6 + a k d t +n cos - (51) kmd d nt + nn ----1-- c 2 2k„ k„ = V2 2 m + n k= n V2 2 m + n — n< t <—0.45n n,a>. x2 = cos 12 12 12 171 Dynamic Analysis of Line Gear Pair Based on Numerical Manifold Method 281 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 where il2 is the transmission ratio, d is the included angle of the gear surfaces, nx is the tooth number of the driving gear, m is the base circle radius of the driving gear, n is the thread pitch parameter, kld is the clearance coefficient of driving line gear, a is the centre distance, d is the diameter of line teeth, and n is the number of the driven gear. The main parameters of the cylindrical helical driving line gear and the corresponding driven line gear are shown in Table 2 as an example. Table 2. Main parameters of the line gears in the example i\2 0 [°] m [mm] n [mm] kid 2.5 150 6 30 20 1.5 k!2 kn n\ d [mm] a [mm] kc 0.95 0.9 15 1.5 120 0.95 The Newmark parameters / and fi' are selected as 0.5 and 0.25, respectively. Let the damping c be zero, and the dynamic response function of the driving tooth is: U ( x, y, z, t ) = 2F z (sin ai (x + m) + sin ay ) + (2/ -l) • sin a,. PAs i=i (2i +l)-®2 [l-(©/®)2] •z © . sin ©t--sin ©f (52) where at = (2i- 1) / 2 s The displacement curves of undamped forced vibration of the driving tooth in x, y and z directions can be drawn as follows: Fig. 5. Undamped forced vibration response diagram of the driving tooth in x, y and z directions From Eq. (33), the natural frequency function is attained: © = d ißiSt) E cos X 4m2 (l + tan2 X)M2\ P ßtst = 1.875, ßtst = k, (i = 2,3,...,n) (53) where p is the material density, m is the base circle radius, X is the helix angle, At is the line tooth length parameter, and E is the elasticity modulus of the material. Eq. (53) indicates a one-to-one correspondence mathematical relationship between the natural frequency and the design parameters of the line gears. The parameter values of the line gears have an applicable range shown in Table 3, and the function relationship is shown in Fig. 6. Table 3. Applicable ranges of parameters affecting the natural frequency d [mm] p [g-cm-3] A [°] m [cm] At [-] 1 to 10 1 to 10 25 to 25 2 to 10 n/3 to 3n a) b) =i 282 Ding, J. - Deng, A.P. - Liu, L.W. - Lu, M.G. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 d) e) Fig. 6. Function relationship between the natural frequency and main parameters of line gears; a) material density, b) line gear base radius, c) helix angle of line tooth, d) line tooth diameter; and e) line tooth length parameter As shown in Fig. 6, the natural frequency of the line tooth is proportional to the diameter of the line tooth but inversely proportional to the helix angle approximately. The natural frequency is strictly inversely proportional to material density, base circle radius and line tooth length. The transmission error can be caused by various types of errors, for instance, installation error and tooth profile error. The line gear meshing stress nephogram and the transmission error curve are shown in Fig. 7. The transmission error results computed with FEM and NMM are essentially consistent. We use an aluminium alloy as the manufacturing material. The elasticity modulus E of the material is 72 GPa. The material density p is 2.8 gcm-3. The Poisson's ratio fi is 0.3. The natural mode of the line tooth can be obtained from the natural frequency. The natural modes of the first four orders and their specific descriptions are shown in Table 4. a) 0 0.2 0.4 0.6 0.8 1 b) Rotation angle of driving gear, angle of initial phase, [rad] $ rotation angle of driving gear, [rad] y transmission error, [°] m excitation frequency, [Hz] 7 REFERENCES [1] Zheng, F., Han, X., Lin, H., Zhang, M., Zhang, W. (2018). Design and manufacture of new type of non-circular cylindrical gear generated by face-milling method. Mechanism and Machine Theory, vol. 122, p. 326-346, D0l:10.1016/j. mechmachtheory.2018.01.007. [2] Chen, Y., Lv, Y., Ding, J., Chen, Z. (2013). Fundamental design equations for space curve meshing skew gear mechanism. Mechanism and Machine Theory, vol. 70, p. 175-188, D0I:10.1016/j.mechmachtheory.2013.07.004. [3] Hai, H., Yin, Y. Q., De Xin, T., Hong, X. L. (2013). Multi-objective optimization design of hard gear-face point-line meshing gear. Applied Mechanics and Materials, vol. 271-272, p. 10271031, DOI:10.4028/www.scientific.net/AMM.271-272.1027. [4] Chen, Y., Huang, H., Lv, Y. (2016). A variable-ratio line gear mechanism. Mechanism and Machine Theory, vol. 98, p. 151163, D0I:10.1016/j.mechmachtheory.2015.12.005. [5] Chen, Y., Liang, S., Ding, J. (2014). The equal bending strength design of space curve meshing wheel. Journal of Mechanical Design, vol. 136, no. 6, art. ID 061001, DOI:10.1115/1.4027160. [6] Yi, Y., Liyan, C., Hang, L., Yiping, D. (2019). Non-linear dynamic response of a spur gear pair based on the modeling of periodic mesh stiffness and static transmission error. Applied Mathematical Modelling, vol. 72, p. 444-469, DOI:10.1016/j. apm.2019.03.026. [7] Meng, Z., Shi, G., Wang, F. (2020). Vibration response and fault characteristics analysis of gear based on time-varying mesh stiffness. Mechanism and Machine Theory, vol. 148, art. ID 103786, DOI:10.1016/j.mechmachtheory.2020.103786. [8] Belingardi, G., Cuffaro, V., Cura, F. (2018). Multibody approach for the dynamic analysis of gears transmission for an electric vehicle. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 232, no. 1, p. 57-65, DOI:10.1177/0954406216674981. [9] Cirelli, M., Giannini, O., Valentini, P.P. Pennestri, E. (2020). Influence of tip relief in spur gears dynamic using multibody models with movable teeth. Mechanism and Machine Theory, vol. 152, art. ID 103948, DOI:10.1016/j. mechmachtheory.2020.103948. [10] Cirelli, M., Valentini, P.P., Pennestri, E. (2019). A study of the non-linear dynamic response of spur gear using a multibody contact based model with flexible teeth. Journal of Sound and Vibration, vol. 445, p. 148-167, D0I:10.1016/j. jsv.2019.01.019. [11] Ebrahimi, S., Eberhard, P. (2006). Rigid-elastic modeling of meshing gear wheels in multibody systems. Multibody System Dynamics, vol. 16, p. 55-71, DOI:10.1007/s11044-006-9021-7. [12] Cooley, C.G., Liu, C., Dai, X., Parker, R.G. (2016). Gear tooth mesh stiffness: A comparison of calculation approaches. Mechanism and Machine Theory, vol. 105, p. 540-553, DOI:10.1016/j.mechmachtheory.2016.07.021. [13] Luo, Y., Baddour, N., Liang, M. (2019). A shape-independent approach to modelling gear tooth spalls for time varying mesh stiffness evaluation of a spur gear pair. Mechanical Systems and Signal Processing, vol. 120, p. 836-852, DOI:10.1016/j. ymssp.2018.11.008. [14] Wang, X., Li, S. (2016). Free vibration analysis of functionally graded material beams based on Levinson beam theory. Applied Mathematics and Mechanics-English Edition, vol. 37, p. 861-878, DOI:10.1007/s10483-016-2094-9. [15] Zhang, J., Ge, R., Zhang, L. (2019). Transverse free vibration analysis of a tapered Timoshenko beam on visco-Pasternak foundations using the interpolating matrix method. Earthquake Engineering and Engineering Vibration, vol. 18, p. 567-578, DOI:10.1007/s11803-019-0522-9. [16] Zhao, G., Wu, Z. (2017). Coupling vibration analysis of rotating three-dimensional cantilever beam. Computers & Structures, vol. 179, p. 64-74, D0I:10.1016/j.compstruc.2016.10.024. Dynamic Analysis of Line Gear Pair Based on Numerical Manifold Method 285 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 275-286 [17] Zhou, D., Fang, J., Wang, H., Zhang, X. (2019). Three-dimensional dynamics analysis of rotating functionally gradient beams based on Timoshenko beam theory. International Journal of Applied Mechanics, vol. 11, no. 4, art ID 19500404, D0l:10.1142/S1758825119500406. [18] Yang, Y., Sun, G., Cai, K., Zheng, H. (2018). A high order numerical manifold method and its application to linear elastic continuous and fracture problems. Science China Technological Sciences, vol. 61, p. 346-358, D0I:10.1007/ s11431-016-9070-8. [19] Yu, C., Liu, F., Xu, Y. (2018). An h-adaptive numerical manifold method for solid mechanics problems. Science China Technological Sciences, vol. 61, p. 923-933, D0I:10.1007/ s11431-017-9143-9. [20] Yang, Y., Sun, G., Zheng, H. (2020). A high-order numerical manifold method with continuous stress/strain field. Applied Mathematical Modelling, vol. 78, p. 576-600, D0I:10.1016/j. apm.2019.09.034. [21] Zhang, Z., Zhang, X., Lue, W. (2010). Numerical method based on compatible manifold element for thin plate bending. Chinese Journal of Mechanical Engineering, vol. 23, no. 1, p. 100-109, D0I:10.3901/CJME.2010.01.100. [22] Wen, W., Jian, K., Luo, S. (2013). 2D numerical manifold method based on quartic uniform B-spline interpolation and its application in thin plate bending. Applied Mathematics and Mechanics-English Edition, vol. 34, p. 1017-1030, D0l:10.1007/s10483-013-1724-x. [23] Vadlamani, S., Arun, C.O. (2019). Construction of beam elements considering von Karman non-linear strains using B-spline wavelet on the interval. Applied Mathematical Modelling, vol. 68, p. 675-695, D0I:10.1016/j. apm.2018.11.042. [24] Chen, Y., Sun, L., Wang, D., Yang, Y., Ding, J. (2010). Investigation into the process of selective laser melting rapid prototyping manufacturing for space-curve-meshing-wheel. Advanced Materials Research, vol. 135, p. 122-127, D0I:10.4028/www.scientific.net/AMR.135.122. [25] Kruth, J.P., Froyen, L., Van Vaerenbergh, J., Mercelis, P., Rombouts, M., Lauwers, B. (2004). Selective laser melting of iron-based powder. Journal of Materials Processing Technology, vol. 149, no. 1-3, p. 616-622, D0I:10.1016/j. jmatprotec.2003.11.051. [26] Wu, W., Yang, Y., Lai, K. (2007). Process analysis of rapid prototyping with selective laser melting. Journal of South China University of Technology Natural Science Edition vol. 35, p. 22-27, D0I:10.1002/jrs.1570. (in Chinese) 286 Ding, J. - Deng, A.P. - Liu, L.W. - Lu, M.G. Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 © 2021 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2020.7179 Original Scientific Paper Received for review: 2021-03-21 Received revised form: 2021-05-06 Accepted for publication: 2021-05-07 Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators Radomir Dokic1* - Jovan Vladic1 - Milan Kljajin2 - Vesna Jovanovic3 - Goran Markovic4 - Mirko Karakašic5 1 University of Novi Sad, Faculty of Technical Sciences, Serbia 2 University North, University Center Varaždin, Croatia 3 University of Niš, Faculty of Mechanical Engineering, Serbia 4 University of Kragujevac, Faculty of Mechanical and Civil Engineering in Kraljevo, Serbia 5 University of Slavonski Brod, Mechanical Engineering Faculty, Croatia Modelling the dynamic behaviour of elevators with high lifting velocities (contemporary elevators in building construction and mine elevators) is a complex task and an important step in the design process and creating conditions for safe and reliable exploitation of these machines. Due to high heights and lifting velocities, the standard procedures for dynamic exploitation are not adequate. The study presents the method of forming a dynamic model to analyse nonstationary vibrations of a rope with time-varying length with nonholonomic boundary conditions in the position where the rope is connected with the cabin (cage) and in the upcoming point of its winding onto the pulley (drum). A unique method was applied to identify the basic parameters of the dynamic model (stiffness and damping) based on experimental measures for a concrete elevator. Due to the verification of this procedure, the experiment was conducted on a mine elevator in RTB Bor, Serbia. Using the obtained computer-experimental results, the simulations of the dynamic behaviour of an empty and loaded cage were shown. In addition, the study shows the specific method as the basis for forming a control program that would enable the decrease in vertical vibrations during an elevator starting and braking mode. Keywords: high-speed elevators, dynamic analysis, a rope with time-varying length, mechanical characteristics of steel ropes, longitudinal oscillations, control program Highlights • The complexity of the dynamic analysis of elevators is because these are systems for lifting (lowering) people and load to great heights (depths) with high velocities and variable parameters. • Determining the critical hoisting velocity of the elevator car can be performed in the function with mechanical characteristics of ropes, such as the elasticity modulus and damping and loads in steel ropes. • Based on the theory of free harmonic damping oscillations, the mechanical characteristics of steel ropes can be determined through the oscillation diagrams obtained by measurement. • By defining the basis for the driving mechanism control program, it is possible to provide minimum dynamic loads of elevators based on adequate models and simulations of their operation in real conditions. 0 INTRODUCTION Mines with underground exploitation and operating levels up to 2500 m, along with the rising number of exceptionally tall buildings, with heights reaching up to 850 m nowadays - such situations require electric elevators with specific characteristics, whose velocity reaches up to 20 m/s, and load capacity up to over ten tons. The elevator quality is estimated according to several important indicators. Safety, comfort, and reliability are especially important features [1]. These indicators depend on, first of all, vibrations occurring while the elevator is in motion. Vibrations are a consequence of driving parameters, inertial characteristics and elasticity of the binding elevator elements. Fig. 1 shows the elements with the greatest impact. An elevator can be divided into two basic parts according to the dynamic impact on the vibration values. The first part is a driving mechanism (engine, reductor, brake, and couplings), while the second part is a cabin lifting system, mostly made of steel ropes for lifting the counterweights on one end and cabin (cage) on the other end and their guide rails. The driving mechanism comprises elements that are much more rigid (q, c2, Fig. 1) and have a smaller mass than the cabin lifting system (c, Fig. 1), which in turn causes the oscillations in smaller amplitudes and higher frequencies. As well as that, bearing in mind that the oscillations are indirectly transferred to the cabin (cage) via ropes, it can be deduced that the lifting system has a much bigger influence on the comfort during the motion than driving mechanism elements does. Deep shaft mines require special mining ropes to hoist personnel and materials safely and efficiently. They are made of round wires that must be either bright or galvanized. The values of the rope safety *Corr. Author's Address: University of Novi Sad, Faculty of Technical Sciences, Trg D. Obradovica 6, 21000 Novi Sad, Serbia, djokic@uns.ac.rs 287 Strajniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 factors for hoisting in mines depend on shaft depth and rope number and are higher in the case of personnel hoisting. According to ISO standards for the area of mine elevators, safety factors vary between 4 and 8 for new ropes and between 3.6 and 6.4 for ropes to be discarded, depending on shaft depth. The highest value of 13 is required by special regulations for hoisting people to depths of up to 600 m. Fig. 1. Elements with the biggest Impact on vibration occurrence In elevators While the elevator is in operation mode, the hoist ropes increase and decrease their free length, so the parameters, such as rope stiffness and damping, are constantly changing [2] to [4]. In high-speed elevators, dynamic instability may occur during lifting (reducing free length) due to increased relative deformation. This instability seriously impacts the safety of the passengers. Since classic models are based on elastic body (rope) oscillations with constant dynamic parameters (mass, stiffness and damping), it is necessary to form dynamic models that will enable the analysis and definition of the dynamic behaviour of elevators with variable parameters [5] and [6]. Because special attention has to be paid to the accuracy of installing and making of cabin guide rails and counterweight in high-speed elevators, the following conclusion arises. Without the addition of external influences, it can be concluded that longitudinal oscillations are dominant, as opposed to transversal oscillations [7] and [8]. Taking into consideration that up until this moment, the problems in driving mechanism vibrations were the subject of a great number of scientific and research papers, with standard analyses as the most frequently applied method, it seems logical that the main focus of dynamic research of elevators should be pointed towards innovative methods for analysing the longitudinal oscillations with variable parameters [9]. 1 DYNAMIC MODELS FOR THE ANALYSIS OF ELEVATOR LONGITUDINAL OSCILLATIONS 1.1 Standard Models Many researchers are interested in studying longitudinal oscillations, and their studies have been based on the general theory about the application of oscillation of elastic bars with constant parameters (mass, stiffness, and damping). These are the so-called standard models [8] and [10]. Those models are acceptable for analysing elevators with low lifting velocities and heights. Figs. 2a and b show models with one and two degrees of freedom and a rope of constant lengths, represented here as Hook's, i.e., Calvin's body. A certain improvement has been made with the analysis of high-lift elevators (> 35 m) and low velocities (till 3 m/s) by using the model represented in Fig. 2c. The model represents a bar of a constant length with an equally spread mass q (kg/m), i.e., it is a model of an elastic body with an unlimited number of degrees of freedom and a concentrated mass at the bottom end as the boundary condition. Based on the analysis which was shown in detail in [11] to [14], in the case when the free rope length is small compared to the cabin mass, it is possible to significantly simplify the dynamic model analysis. Fig. 2. Standard dynamic models; a) with one degree of freedom; b) with two degrees of freedom; c) a "heavy" constant length bar Fig. 3 shows an oscillation diagram for the first three harmonics. Due to the very small oscillation amplitudes in higher harmonics, their influence can be neglected. Thus, the total oscillation process with infinite degrees of freedom, whose total oscillation 288 Bokic, R. - Vladic, J. - Kljajin, M. - Jovanovic, V. - Markovic, G. - Karakasic, M. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 form is shown in Fig. 3 with a dashed line (a), can be replaced with a straight line (a), with satisfactory accuracy. In other words, it is replaced with a system with one degree of freedom and constant dilatation (s) along the free rope end. d a ,b x Fig. 3. Oscillation shapes (forms) of the first three harmonics (a, b and c), and the summary oscillation form for a = 0.1 In accordance with the above, regarding high-lift elevators (> 35 m) and low velocities (to 3 m/s), it is roughly possible to create load oscillation models with one degree of freedom, with a "heavy" spring, which was studied in general literature. Also, it is necessary to replace the total mass (of both load and rope) with an equivalent mass Me = M+ (1/3) • qL, reduced in the cabin place, [15] and [16]. 1.2 Dynamic Elevator Models with Dynamic Variable Parameters Fig. 4a shows the most common solutions of the lifting systems for high-speed elevators with a driving pulley, while the corresponding dynamic model is shown in Fig. 4b. a) b) Fig. 4. Elevator models; a) Koepe system, b) a dynamic model of a high-speed elevator In order to secure comfort during the motion, control programs are used in contemporary elevators. They define the circumferential velocity of the pulley. Thus, they also define the cabin motion velocity (a kinematic condition), as opposed to the previous period when the motion velocity depended a great deal on the driving electromotor's mechanical characteristics and brake system (the dynamic equilibrium condition). In the earlier periods, replacing one-speed engines with two-speed engines was observed as a significant improvement. This improved the motion comfort in braking instances and aided the accuracy of stopping the cabin. As for the process of a regular elevator, in cases in which there is no slipping of the steel rope on the driving pulley and when the driving characteristic is represented via a rope velocity at the meeting point of the rope and pulley, the elevator model can be simplified and represented in the form shown in Fig. 5. Upon observing just the upcoming rope end, the model can be represented as a system with an unlimited number of degrees of freedom; at one end it is rolling onto the pulley at a v(t) velocity, while on the other it is burdened with concentrated mass. Due to the variable rope length during the motion, the stiffness (c= EA/L) changes. This is a characteristic of parametric oscillations and contributes to the possible occurrence of resonance. To this end, it is necessary to complete certain steps in the analysis of dynamic behaviour. The critical lifting velocity, during which the unstable motion occurs, i.e., the rope strain is increased when its free length is reduced, needs to be determined. jo Fig. 5. Elevator model with a rope of a variable length with boundary conditions The deformation of an arbitrary cross-section represents the function of the position x, and the time t, i.e.: u = f (x, t). (1) Upon observing the equilibrium of the elementary part (dx), it can be deduced that: Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators 289 Strajniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 q ■ dx d2u(x, t) g dt2 = -S ( x, t ) + S ( x, t ) dS q ■ dx H--dx + q ■ dx +--a. dx g (2) By representing the rope as Calvin's model, where the influence of internal friction can be taken into consideration via the so-called rope resistance force factor (bf), the dependence of the inner force in the rope on its deformation can be noted in this form: 5 ( x, t ) = E • A •— dx u( x, t ) + bf du ( x, t ) • dt (3) If the Eq. (2) is divided with (q-dX) / g and a replacement for S(x,i), this is obtained: 82u(x, t) g ■ E ■ A 8 St2 q +g ± a. — | u(x, t) + bf 8u( x, t) 8t (4) In Eqs. (2) and (3), (E) is the rope elasticity modulus whose magnitude depends on the elasticity modulus of wires (Er = 2.1 • 105 MPa) and the construction of the rope [17] and [18], whose value can be twice smaller (stranded wire ropes). By using the differential equation, (Eq. (4)) and the equilibrium condition of moments on the driving pulley, it is possible to form a system of equations that describes a dynamic equilibrium on the driving pulley in the case of a model shown in Fig. 4b in this form [12] and [19]: q d2u1(x,t) g dt2 = E • A • dx2 +q-|i ± -g u1 ( x, t ) + bf q d2u2(x,t) 82 ---;-= E ■ A--- g 8t 8x u2 ( x, t ) + bf du1 ( x, t ) dt 8u2 ( x, t) 28t (5) +q ■(I ± - M = —. E • ^ • — (6) i .-q dx O U (l1, t ) - u2 (l2, t ) + bf— [ul (l1, t ) - w2 (l2, t )] dt a • i - J1 • R • (7) 1.3 Boundary Conditions In order to solve a partial differential equation, (Eq. (7)), i.e., an equation system (Eqs. (5) to (7)), it is necessary to define the boundary conditions in the incoming point of the rope to the pulley in the meeting point of hoist ropes and the elevator cabin. Boundary conditions at point C, where the rope makes the first contact with the pulley, are dependent on whether the winding is with or without rope slipping. Fig. 6a shows the distribution of the forces on the wrap angle of the pulley (a) for a quasistatic case of elevator operation (without the influence of dynamic forces). The regular elevator operations must not allow for the rope to slip in the whole wrap angle (a), which is regulated with a safety degree against slipping, i.e., with the existence of a suitable angle (aM) with the so-called relative abeyance of the rope on a driving pulley. It should be noted that, in quasistatic conditions, the zone with the elastic slipping of the rope on a driving pulley (aK), Fig. 6a, always occurs on the descending side. However, due to the oscillation, the force in the incoming end changes, so slipping is probable in the pulley's incoming zone. The change of force in the wound rope part can only be maintained if this change is smaller than the adhesive force making it possible. Fig. 6d shows different cases of force distribution over the wound rope length as a function of elevator velocity. With high-speed elevators, it should be expected that the force change in the wound part of the rope is smaller than the friction force change in the incoming zone (curve v3, Fig. 6d). Basically, in this case, it can be accepted that there is no slipping in the rope and rope pulley meeting point. The problem is analysed in great detail in [12] and briefly summarized in [20]. Fig. 6. Boundary conditions; a) force distribution on the wrap angle of the pulley, b) driving pulley without slipping, c) cage (cabin), d) different cases of the force distribution 2 2 290 Bokic, R. - Vladic, J. - Kljajin, M. - Jovanovic, V. - Markovic, G. - Karakasic, M. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 The boundary condition at the meeting place of the rope and the pulley without slipping, Fig. 6b (point C), is presented in this form: , * =1 t du (l, t) ( dl o dx \ dt dt. (8) The boundary condition in the meeting point of the rope and elevator cabin (mine elevator cage) or counterweight is given in this form: Q = E - A-Af u ( L, t ) + bf•dUÎM dx { f dt ( Q d2u(L, t) g ' dt2 \ - F, (9) where is Ff FN • sign(v) friction force between sliding guide shoes and guide rails. Fig. 7. Deformation of a rope 1.3.1 Estimation of Critical Velocity The non-integrated boundary condition (Eq. (8)) prevents the solution of the partial differential equation (Eq. (5)) so that the solution can be sought through the formation of integral equations, which contain both the differential equations and the corresponding boundary conditions. Simplifying the mathematical model, friction within the boundary condition is omitted in Eq. (9). In Fig. 7, a weightless fibre is shown, loaded at point (5) by the force Sj, where the displacement of the rope point (without sliding over the drive pulley) in the region s x , for s < x (10) and displacement u(x) = K(x,s)S. For the case with several acting forces with acting points at x=Sj there is a displacement: i( x) = ^ K ( x, st ) • S. (11) For the load case with evenly distributed rope load (rope's weight), there is the expression: rL i(x) = J K(x, s) ■ q ■ ds. (12) By applying this procedure to the differential equation (Eq. (5)), replacing the floating argument with (s), multiplying by the function K(x,s,l) and performing the necessary mathematical transformations, the rope deformation in the form is obtained [20]: u( x, t ) = - £ K ( x, s, I ) • q{s) • I ^ - g dt2 ds r' MhH^t)dt + E • A\LK(x,s,l) pjy Jl du dx -b•E•A• ds dt ds d 2u (L, t) dxdt K ( x, L, l ). (13) Using the method of particular integrals, taking into account only the first form of oscillations of the rope and that the length of the free rope (L-l) changes "slowly" with time, and if a new variable v(x,t) is introduced, so that: u(x, t) = v(x, t) + u(l, t) , the deformation of rope free side in the phase of load-lifting will be: v( x, t ) = ( x -1 ) • h0 x -1 L - ln Ï4 E • A Q- L -1 q • (2 • L + x +1) 2 • cos (j; m(i )dt ) i+^ (14) where: ho =~ q a E ■ A ' g ' Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators 291 =i Strajniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 œ\l ) = Q g ■ E ■ A q ■ ( L -1 ) ' 3 ( L -1 ) m = - b ■ f ' m2(l)-dt. 2 Jo v ' Analysing Eq. (7), it can be concluded that during the reduction of the free length of the rope (load-lifting), its deformation can increase under certain conditions, which causes a permanent increase in the dynamic load of the hoist rope. Such a phenomenon, which is usually described as unstable lifting, will not occur if the lifting velocity is less than the critical velocity, which is defined by the expression: = 2 • b • E • A q • ( l -1) 3 ■ = 2 • b Q E • A ■ ~M~ ' (15) where Me =1 Q + q • ( L -1 ) 3 Ig is the reduced mass of the cabin and the rope part as a "heavy spring". The level of critical velocity depends on the damping characteristic of the wire rope and the elevator-s basic characteristics. When it comes to elevators with high load capacity, high velocity and lifting height (high-speed elevators and mine elevators), the lifting velocity can exceed the critical velocity, which is why it is necessary to check the stability of the lifting process already in the design phase. It should be taken into account that Eq. (15) does not contain the friction in guide rails, which adds to the oscillation reduction, and simultaneously increases the critical velocity value. 2 EXPERIMENTAL DETERMINATION OF DYNAMIC PARAMETERS Experimental studies shown in this paper represent the sequel to the research in [13]. The measuring was done on a mine elevator in RTB Bor, Serbia; the maximum projected cage lifting velocity was 16 m/s, and the projected lifting height was 523 m in the first phase and 763 m in the second phase of mining. The driving mechanism is shown in Fig. 8. The other significant characteristics of this exploitation machine can be seen in Tab. 1. A mine shaft is like a round cross-section with a 10 m diameter. The transfer of the driving moment to the hoist ropes is done through friction (Koepe system) from the grooved drum with a 2.5 m diameter. Fig. 8. Driving machine of a mine elevator and the connecting tools for connecting the hoist ropes and the cage Table 1. Mine elevator technical characteristics Driving electric motor Power: 1500/2860 kW, number of revolutions: 122.2 rpm, torque: 117.2/233.4 kNm_ Mass: 13 t (it includes the cage and the connecting tools for connecting the hoist ropes), Fig. 8_ Cage (cabin) Counterweight Mass: 21 t z = 6 pieces, d = 27 mm (150 wires per cross-section). Lang's lay ropes, three right-hand Lang's lay ropes, and three left-hand Lang's lay ropes. Breaking force: 561 kN, tensile strength of the wires: 1700 MPa. Unit mass: 3.02 kg/m. Hoist ropes z = 2 pieces, d = 50 mm (222 wires per cross-section). Connected with the cage via rotating hooks to prevent unwinding. Unit mass: 9.64 kg/m. Compensation ropes 2.1 Equipment Used and the Method for Measuring Data Acquisition The following measuring equipment was used for the mine elevator experiment: • Universal 8-channel measuring amplifier (2 pieces), QUANTUM X MX480B, pos. 10a and 10b (Fig. 9), • A computer for storing the measuring signals, pos. 13 (Fig. 9), • Antennas for the wireless transfer of the measuring signal (2 pieces) NanoStation loco NS2L, made by IBIQUITI NETWORKS; 5 km range with the antennas optical visibility, pos. 12a and 12b (Fig. 9), • Incremental encoder - the number of revolutions sensor AINS 41, made by Meyer Industrie-Electronic GmbH - MEYLE; measuring range up 292 Bokic, R. - Vladic, J. - Kljajin, M. - Jovanovic, V. - Markovic, G. - Karakasic, M. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 to 6000 rpm and with an allowed axial/radial load at the output shaft of 30/20 N, pos. 9 (Fig. 9), Optical measuring device for the number of revolutions ROS type, made by Monarch Instrument; measuring range varies from 1 to 250 000 rpm, pos. 1 (Fig. 9), Accelerometer HBM B12, pos. 2 (Fig. 9), Strain gauges HBM LY41-6/120, pos. 3 to 8 (Fig. 9), Software for acquisition and processing of measuring signals, HBM catmanEasy-AP. Fig. 9. Schematic layout of the measuring places The experiment was done according to the principle of summary data acquisition via a wireless data transfer by using the NanoStation loco NS2L antennas (pos. 12a and 12b, Fig. 9). In this way, a complete synchronization of the collected data on the cage roof and machine room was achieved. The drum number of revolutions was measured with an AINS 41 incremental encoder (pos. 9, Fig. 9). It was set on the adaptive carrier with a 127 mm diameter measuring wheel, which was directly leaned onto the rim of the brake disc, Fig. 10a. The encoder measuring signal was led to the amplifier (pos. 10b, Fig. 9), and then via a LAN switch (pos. 11b, Fig. 9) and the antenna (pos. 12b, Fig. 9) wirelessly connected to the other radio antenna positioned on the cage (pos. 12a, Fig. 9). The measuring signal was stored in a common computer file (pos. 13, Fig. 9) placed on the cage. It was done via another LAN switch (pos. 11a, Fig. 9). The lifting and lowering of the velocity of the cage was measured via a roller guide, whose change in the number of revolutions was registered with an optical sensor (pos. 1, Fig. 9), positioned on a small magnet table, Fig. 10b. The change in the cage's acceleration during the measuring process was registered by an HBM B12 accelerometer (pos. 2, Fig. 9), set on the connecting tools of the cage through a magnetic holder, Fig. 11a. The force changes in the hoist ropes were monitored by measuring the deformation of the connecting tools (pos. 3 to 8, Fig. 9). The deformations were measured on each connecting tool (out of 6 in total) by using the HBM LY41-6/120 strain gauges, Fig. 11b. ajD^^^^^^^^M^B b) Fig. 10. Measuring places; a) the incremental encoder connected with the measuring wheel positioned on the rim of the brake disc, and b) the optical sensor (number of revolutions sensor) positioned on the cage aj^B^HKB b)l Fig. 11. Measurements on the cage connecting tools; a) cage acceleration, and b) deformations All measuring signals were led to the other eight-channel measuring amplifier (pos. 10a, Fig. 9), and they were recorded on the computer (pos. 13, Fig. 9) via a LAN switch (pos. 11a, Fig. 9). 2.2 Protocols and Results of the Experiment In Fig. 12, a schematic layout of a mine shaft elevator is shown with labelled levels, whose positions were set with altitude, using plus and minus. The ground level was at an altitude of +436 m, while the machine room was placed at +465 m. Experiment results and determination of the dynamic model parameters in this paper are shown in four characteristic motion cases (lifting and lowering) of the cage, with and without load. A significantly higher number of measurements was completed for different cases (different positions for starting and stopping the cage motion etc.), and the results are presented in [11]. Here, the results of the lowering and lifting of an empty cage and a loaded cage (traction machine) are presented. According to those, the system oscillation parameters are defined. The characteristics of experimental cases are as follows: I. Lowering an empty cage from position +104 m to position +56 m (48 m). The action of Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators 293 Strajniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 Fig. 12. Schematic layout of a vertical mine shaft in "Jama" mine of RTB Bor II. III. IV stopping the cage was performed by suddenly switching off the driving electromotor (at a high velocity, Fig. 13a). Lifting of an empty cage from position +56 m to a position +199 m (143 m). Stopping the cage was done just like in the case above, with the sudden switching off the power to the driving electromotor, at a high velocity, Fig. 13b. Lifting a loaded cage (the traction machine for wagons, mass ~9.35 t) from position a -71 m, to position +52 m (123 m). Stopping the cage was done at a low speed, Fig. 13c. Lowering the loaded cage (the traction machine for wagons, mass ~9.35 t) from position +52 m to -71 m (123 m). Stopping the cage was done at a low speed, Fig. 13 d. The results of the conducted experimental research can be seen in the following figures, which are showing the changes in the winding velocities of hoist ropes that were wound onto the drum; cage acceleration and changes in the forces in the hoist ropes, i.e., on the elements for connecting the ropes to the cage. The results are valid for all four chosen examples of motion. The diagram shows the character of the changes of the said parameters. It shows the cage oscillation amplitudes that occur after stopping the driving machine. This part of the diagram (after stopping the driving machine) is used to set the mechanical characteristics of steel ropes (harmonic oscillations). Fig. 14. Diagram with the results obtained during lowering an empty cage and sudden stopping, at a "high" velocity (I motion case) Fig. 13. Mine elevator parameters that are relevant for the analysis; a) lowering an empty cage, b) lifting of an empty cage, c), lifting a loaded cage, d) lowering the loaded cage Fig. 15. Diagram with the results obtained during lifting the empty cage and sudden stopping at a "high" velocity (II motion case) 294 Bokic, R. - Vladic, J. - Kljajin, M. - Jovanovic, V. - Markovic, G. - Karakasic, M. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 Fig. 16. Diagram with the results obtained during the loaded cage lifting process (III motion case) Fig. 17. Diagram with the results obtained during loaded cage lowering process (IV motion case) 2.3 Determination of Dynamic Parameters Based on Measurement Results The stiffness and damping parameters can be determined, and the elasticity modulus as mechanical characteristics of a steel rope can be indirectly determined too, by using the above diagrams and determining the oscillation period, i.e., the frequency and logarithm decrement of the damping. In [13], and especially in [17], the said characteristics are described in detail, so this study only mentions the important parts for explaining the results obtained with new measurements. Stiffness and elasticity modulus. Rope stiffness is defined with an expression [13]: E • A c =-. (16) L -1 (t) The stiffness depends on the change in the rope's free length, but it also depends on the elasticity modulus, usually taken as a constant parameter in dynamic behaviour analyses. However, due to the complex construction of a steel rope (consisting of wires, strands, and a core), its change depending on the strain magnitude and loading character (loading-unloading) cannot be neglected. According to [18], a distinction can be made between the so-called tangent elasticity modulus (Et) and an average (medium) elasticity modulus between two stresses (ES), Fig. 18. The tangent elasticity modulus represents a theoretical angle of the inclination of the tangent to the c-e curve for the current stress value (ctz). The constant variability in the elasticity modulus magnitude can be observed, and significant differences occur in loading and unloading the rope. The average (median) value of elasticity modulus can be set for boundary working stresses - Es Slower, ^upper). Dilatation [ Fig. 18. Tangent and average rope elasticity modulus Damping. In many practical systems, oscillation energy is gradually transformed into heat or sound, the effect mostly known as damping. Although the quantity of such energy is relatively small, it is important to dedicate some attention to damping to predict the oscillatory system's response, such as the elevator systems considered in this paper. Damping at the elevator is complex, and it happens because of the rope's inner friction (viscous and hysteretic damping), Coulomb's friction on guide rails and damping due to the airflow around the cabin (cage) in the elevator shaft [13] and [15]. Since the presented experiment results (Tab. 2) refer to the moment of stopping the cage (v= 0), damping as a consequence of airflow around the cage can be neglected. Fig. 19. An oscillating system with Coulomb's friction; a) motion with the load, b) a model with Coulomb's damping The friction on the elevator guide rails creates the damping force (Coulomb's friction), which is constant in its magnitude, but it is of the opposite direction compared to the motion of the oscillating load. Since this type of elevator requires special attention to be paid to the guiding accuracy and reducing the guide Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators 295 Strajniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 rails friction, in the example with centrical cabin load, the friction forces can be neglected in relation to the total load. Nevertheless, the influence can be of importance for cabin oscillation analyses [21]. A typical diagram of an oscillator with Coulomb's friction is presented in Fig. 19a. The declining of amplitudes, as opposed to inner friction oscillation, is linear in its nature [15] and [17]. As it was explained in [13], the total force of Coulomb's friction on guide rails is: F = n ■M' Fn (17) For the most part, the literature presents the inner friction in the rope as viscous friction, like in homogenous bodies, Figs. 2a and b. Nevertheless, while the rope is being deformed, the energy dissipation also appears due to the friction between wires and strands, which slide against one another during the deformation process. This creates the hysteresis loop (Fig. 20a). The effect causes the damping known as hysteretic or structural damping. The loss of energy per unit of material volume in one loading and unloading cycle is equal to the area closed by the hysteresis loop. Experiments have confirmed that the energy loss per cycle is approximately proportional to the square of the oscillating amplitude. A model with hysteretic damping is presented in Fig. 20b. Similarly to the above, the rope's inner friction (br) is a combination of viscous and hysteretic damping, Fig. 20c. Fig. 20. Hysteretic damping; a) hysteresis loop, b) a model with hysteretic damping, and c) an equivalent model with internal rope damping The total damping of the considered system is an overall influence of the damping due to the rope's internal friction and the damping caused by the friction in guide rails, so it can be represented by the parallel connections in the oscillator model, Fig. 21b. This is seen in the diagrams shown in Figs. 14 to 17, which show the cumulative influence of rope's internal friction and the guide rails' friction is generally of a viscous damping character. The size of the impact caused by rope damping and damping due to guide rails friction can be defined via the "overlapping" of the diagram measurement results and simulation results, as presented in Fig. 22. Based on the theory about free harmonic damping oscillations (Fig. 21a), it is concluded that for determining the dynamic parameters, one should take the part of the cage oscillation diagram after the driving machine is stopped. Fig. 21. Oscillating system with viscous integral damping; a) damped oscillations b) a lifting system model 160 165 Time [si Fig. 22. The impact of inner and Coulomb's friction on oscillating systems when b=5161 Ns/m (Table 2) By measuring the amplitudes and oscillating periods of free damping oscillations (changes in cage acceleration), the damping coefficient (S) can be determined via a logarithmic decrement, and based on this, a resistance force coefficient (b): D = ln^- = iln-X- = S-T D, T b = 2-8-M„ (18) (19) Determining the stiffness coefficient (c) and elasticity modulus (E) of hoist ropes can be performed 296 Bokic, R. - Vladic, J. - Kljajin, M. - Jovanovic, V. - Markovic, G. - Karakasic, M. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 Table 2. Values obtained by measuring and parameters of hoist ropes L [m] Me [kg] S [s-1] m [rad/s] c [N/m] b [Ns/m] E [MPa] S/m I (Fig. 13a) | 409 18 797 0.33 4.67 409 163 12 222 98 440 0.070 II (Fig. 13b) Î 266 20 525 0.21 5.69 663 329 8718 103 792 0.037 III (Fig. 13c) Î 413 28 099 0.09 4.26 510 803 5161 124 095 0.022 IV (Fig. 13d) | 536 26 613 0.18 3.88 401 206 9639 126 498 0.047 with the measured oscillation parameters (oscillating period, oscillating amplitudes, etc.): c = M, -m E = - A (20) (21) Fig. 23 is an illustration of determining the said mechanical characteristics of hoist ropes based on the cage acceleration changes diagram (after the driving machine was stopped) in the III motion case. Time [s] Fig. 23. The measured values of cage acceleration for III motion case Based on the diagram in Figs. 14 to 17 and the expressions (Eq. (18) to (21)), mechanical characteristics of steel ropes are defined, as shown in Table 2. Upon analysing the measured results, it can be concluded that the data for elasticity modulus are in agreement with the data from the literature [18]. The higher values of elasticity modulus for loaded cabins are seen as a consequence of the rope constructed by putting the lays of wires in strands and strands into a rope. This confirms the validity of the applied procedure, enabling us to define the real (exploitation) values in mine elevators. The damping coefficient values, for which there is no significant comparative data, are not constant in size but appear to be different in the analysed experiments. Regarding the limitations in conducting the experiments in working conditions, it was concluded that the experiment should be prepared in a laboratory setting in the future. 2.4 Computer Simulations of Mine Elevator Dynamic Behaviour and Result Correlation In order to simulate the mine elevator operation process, the study has set the change in hoisting rope velocity at the point of winding onto the driving drum as the change (diagram), obtained by direct measuring on the driving machine via an incremental encoder (pos. 9, Fig. 9). An oscillation diagram was separated with the purpose of verifying the dynamic model during the whole period while the cage was moving, as well as after the complete stop in the driving machine moving. The acceleration diagram in the III motion case is shown, "overlapping" the diagram with the measured results, Fig. 24. Judging by the diagram in Fig. 24, it can be concluded that the dynamic model with dynamic parameters, determined by the measuring with satisfactory accuracy, shows the real behaviour of a mine elevator. As a result, this makes computer simulations possible, which will analyse the real loads of the exploitation facility [11] and [22]. The described method is a new approach that will ensure the analysis of dynamic behaviour in systems for vertical hoisting, which can be found in exploitation. Time [s] Fig. 24. Cage acceleration diagrams in the motion case III (Fig. 13c) After comparing the results of numerical analysis with measurement results, a conclusion arises that it would be useful to have some form of numerical value for result correlation. In the literature, this value is known as Pearson's correlation coefficient. This Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators 297 Strajniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 study did not determine the correlation coefficient for individual motions, so that it could be taken as one of the ideas for future research. As an illustration of the possibility of dynamic analysis, the following section shows diagrams that represent changes in individual values. They were created in specialized software for dynamic analysis using the experimental data shown in Table 2. Fig. 25. Diagram showing the changes in the cage position, velocity, and acceleration, for I motion case Time [s] Fig. 26. Comparative diagrams showing the changes in the elongation, stiffness coefficient, and resistance force coefficient of hoist ropes, for I motion case Fig. 29. Diagram showing the changes in cage position, velocity, and acceleration, for III motion case ----- ......... w S-.........* -rope elongation (?/) ---stiffness coefficient (c) resistance force coefficient (/>) -ppw*- 5 5 z Time [s] Fig. 30. Comparative diagrams showing the changes in the elongation, stiffness coefficient, and resistance force coefficient of hoist ropes for III motion case o -20 S "40 c -60 o -so ¿-100 -120 -140 05 ~ 0 —-0.5 X..... -cage velocity (v) ---cage acceleration (a) \ T "V \ |r ■ — 1 Time [s] Fig. 31. Diagram showing the changes in cage position, velocity, and acceleration for IVmotion case 5 ^4 I rr.:. ---- — V \ 'WiÎ{»»l / : \ WMïi11': / \ / —cage v elocity (v) xeleration (a) osition (A) ----- — cage p 40 50 Time [s] 2 i 'iT "p 0 — -i! cd -2 jy o ■2 o -4 Fig. 27. Diagram showing the changes in cage position, velocity, and acceleration, for II motion case Time [s] Fig. 28. Comparative diagrams showing the changes in the elongation, stiffness coefficient, and resistance force coefficient of hoist ropes, for II motion case Time [s] Fig. 32. Comparative diagrams showing the changes in the elongation, stiffness coefficient, and resistance force coefficient of hoist ropes for IV motion case 3 AN ANALYSIS OF THE HOISTING VELOCITY IMPACT AS THE BASIS FOR CONTROL SYSTEM The material of this paper confirms that there are different analyses of individual parameters for vertical hoist systems. One part is of special interest, and that is the research on the changes in hoist velocities, i.e., defining the control system [23]. The above method can be used for determining an optimal form 298 Bokic, R. - Vladic, J. - Kljajin, M. - Jovanovic, V. - Markovic, G. - Karakasic, M. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 for changing the hoist velocity, to reach maximum capacity, i.e., the shortest time for the motion, along with securing a high level of motion, especially in modern passenger elevators, whose velocities reach up to 20 m/s. As stated above, the dynamic parameters during the elevator or exploitation facility design process can only be determined with limited accuracy, so the exploitation deviation can be significant. The method shown in this paper can be used as a basis, which shows that it is possible to apply straightforward methods and form an optimal elevator facility. The application of the illustrated method makes it possible to determine real parameters and perform dynamic analysis even after the installation process. Therefore, the method makes it possible to adjust and correct the exploitation facility and improve the previously defined optimal operating conditions. This particularly concerns defining the appropriate control program, which ensures the shortest motion time with the maximum level of comfort, i.e., with the control of the dynamic loads. The following figures present some of the possible ways of impacting driving kinematic characteristics through a control program. Fig. 33 shows how the changes in velocity form characteristics can impact the cabin (cage) oscillating amplitudes, as well as the motion comfort. As it appears, it is important to note that by controlling the drive in the acceleration period, and especially by choosing the right moment to switch from acceleration to stationary velocity, the system can be significantly "relieved", i.e., the motion comfort can be bettered (the amplitudes can be quite decreased), pos. 4 in Fig. 33. Installing an acceleration transducer (accelerometer) on the points where ropes connect with the cabin (cage) has acceleration change as a response at any moment. This is vital with high-speed elevators, in cases in which nominal velocity is not reached between two neighbouring floors, but the acceleration is followed by braking, Fig. 34. The right control is key for braking at the most convenient moment, as shown by line (2) in the 5 10 15 0 4 8 12 Time [s] Time [s] a) b) Fig. 33. The influence of driving system kinematic characteristics on the motion comfort; a) change of velocity, b) change of acceleration Fig. 34. The layout of the influence of controlling the magnitude of oscillating system amplitudes; a) change of velocity, b) change of acceleration Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators 299 Strajniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 velocity change diagram, Fig. 34a, which causes much smaller oscillating amplitudes, Fig. 34b. The distance the cabin (cage) has passed and the accuracy of station landing must be taken into account when creating a control system. 4 CONCLUSION The dynamic analysis and scientific approach shown in this paper are of special importance taking into account the nature of the systems; they are meant for lifting (lowering) people and load to great heights (depths) at velocities up to 20 m/s. The complexity of the dynamic analysis can also be seen because these are oscillating processes with variable parameters and boundary conditions at the incoming/outgoing points of the steel rope and the driving pulley. Furthermore, the stiffness and elasticity modulus of steel ropes are not constant in magnitude but rather change depending on the changes in cage position, as well as construction, stress, and rope exploitation time. The dynamic model that is suitable for the analysis can be formed by observing a particular system and by simplifying (i.e., omitting) the small values of higher order. The presented approach on forming the differential motion equation for the incoming end of the rope onto the driving pulley helps determine the critical hoisting velocity in the function with mechanical characteristics (elasticity modulus and damping) and stress (load) in steel ropes. Moreover, this method also makes it possible to execute computer simulations and dynamic behaviour analysis during the motion by using the appropriate specialized software for the dynamic analysis of mechanical systems. The paper presents an experimental method that can be used to determine dynamic parameters such as stiffness, elasticity modulus, and damping in exploited elevator steel ropes. Specific elasticity modulus values presented in Table 2 reveal the noticeable dependence on the load, i.e., rope stress. In addition, the damping coefficient in ropes is not a constant size, but it depends on the position of a cage. It can be deduced that ropes experience a combination of viscous damping and hysteretic damping, which should be further investigated in these systems. By forming adequate dynamic models, determining real values for dynamic parameters, and simulations of dynamic behaviour for real facilities, it is possible to define the basis for a program that would control the driving mechanism. That would ensure minimal dynamic loads, which is especially relevant for comfort during the motion, while simultaneously securing the optimal motion time, i.e., the efficiency of high-speed passenger elevators, especially in transient operation regimes. As for mine elevators, strict safety conditions are required in order for the facility to operate, especially regarding elevators for the transport of people. It is necessary to secure special conditions at the moment of power transfers and friction motion to avoid the slipping of the rope as a whole along the driving pulley. Finally, this is vital from the point of view of defining suitable boundary conditions for a dynamic model at the points where the ropes get on or off the driving pulley. 5 NOMENCLATURES uj, U2 elastic deformation of the rope on the incoming and outgoing rope end, [mm] E rope elasticity modulus, [MPa] A rope's cross-section area, [mm2] a acceleration of the driving mechanism, [m/s2] Mm driving motor torque, [Nm] i gear ratio, [-] n driving mechanism efficiency, [-] J moment of inertia of rotational masses, reduced to the shaft of a driving pulley, [kgm2] R driving pulley radius, [m] q rope weight per meter, [kg/m] l wound rope length, [m] v (v = dl/dt) circumferential velocity, [m/s] Ff friction force between sliding guide shoes and guide rails, [N] l0 length of winded rope at time t = 0, [m] nv (nv = 4) number of the roller guide groups, [-] nt (nt = 3) number of rollers in the guide group, [-] p rolling resistance of the roller guide on the guide rail, [-] FN force of the roller guide pressure on the guide rail for centric load, which depends on the spring tightening during installation, [N] Xi (also, xi+j and xi+l) measured oscillating amplitudes, [mm] T measured oscillating period of free damping oscillations, [s] Me reduced oscillating mass (Me = M+ (q • L(t)) / 3), [kg] M total mass hanging on hoist ropes (loaded cage and compensation ropes), [kg] L(t) hoist ropes' free length (L(t) = L - J v(t)dt), [m] v(t) circumferential velocity of the pulley (drum), [m/s] 300 Bokic, R. - Vladic, J. - Kljajin, M. - Jovanovic, V. - Markovic, G. - Karakasic, M. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 287-301 m circular frequency of free oscillations ((» = J(b2 +V ), [rad/s] a> damping oscillation frequency (co = 2n/T), [rad/s] L hoist ropes length, [m] S damping coefficient, [s-1] c ropes stiffness coefficient, [N/m] b resistance force coefficient, [Ns/m] 6 REFERENCES [1] Peng, Q., Xu, P., Yuan, H., Ma, H., Xue, J., He, Z., Li, S. (2020). Analysis of Vibration Monitoring Data of Flexible Suspension Lifting Structure Based on Time-Varying Theory. Sensors, vol. 20, no. 22, 6586, D0l:10.3390/s20226586. [2] Zhang, Q., Yang, Y., Hou, T., Zhang, R. (2019). Dynamic analysis of high-speed traction elevator and traction car-rope time-varying system. Noise & Vibration Worldwide, vol. 50, no. 2, p. 37-45, D0I:10.1177/0957456519827929. [3] Bao, J., Zhang, P., Zhu, C. (2011). Modeling of Rope Longitudinal Vibration on Flexible Hoisting System with Time-Varying Length. Applied Mechanics and Materials, vol. 130-134, p. 2783-2788, D0I:10.4028/www.scientific.net/ AMM.130-134.2783. [4] Zhang, Y., Pota, H.R., Agrawal, S.K. (2002). Modification of residual vibrations in elevators with time-varying cable lengths. Proceedings of the American Control Conference Anchorage, vol. 6, p. 4962-4966, D0I:10.1109/ACC.2002.1025450. [5] Arrasate, X., Kaczmarczyk, S., Almandoz, G., Abete, J.M., Isasa, I. (2014). The modelling, simulation and experimental testing of the dynamic responses of an elevator system. Mechanical Systems and Signal Processing, vol. 42, no. 1-2, p. 258-282, D0I:10.1016/j.ymssp.2013.05.021. [6] Herrera, I., Su, H., Kaczmarczyk, S. (2010). Investigation into the damping and stiffness characteristics of an elevator car system. Applied Mechanics and Materials, vol. 24-25, p. 7782, D0I:10.4028/www.scientific.net/AMM.24-25.77. [7] Kimura, H., Ito, H., Nakagawa, T. (2007). Vibration analysis of elevator rope (Forced vibration of rope with time-varying length). Journal of Environment and Engineering, vol. 2, no. 1, p. 87-96, D0I:10.1299/jee.2.87. [8] Guo, Y., Zhang, D., Chen, K., Feng, C., Ge, S. (2018). Longitudinal dynamic characteristics of steel wire rope in a friction hoisting system and its coupling effect with friction transmission. Tribology International, vol. 119, p. 731-743, D0I:10.1016/j.triboint.2017.12.014. [9] Watanabe, S., Okawa, T. (2018). Vertical vibration of elevator compensating sheave due to brake activation of traction machine. Journal of Physics: Conference Series, vol. 1048, no. 012012, D0I:10.1088/1742-6596/1048/1/012012. [10] Li, C., Hua, C., Qin, J., Zhu, Z. (2019). Research on the Dynamic Characteristics of High-Speed Elevator System. Advances in Engineering Research, vol. 181, p. 105-109, D0I:10.2991/ ice2me-19.2019.23. [11] Bokic, R. (2016). Dynamics researching and development of vertical transport machines using numerical-experimental procedures. PhD thesis, University of Novi Sad, Faculty of Technical Sciences, Novi Sad. [12] Vladic, J. (1982). Contribution to the determination of the safety degree against slipping of a dynamically loaded rope in the power transmission system by a driving pulley. MSc thesis, University of Novi Sad, Faculty of Technical Sciences, Novi Sad. [13] Vladic, J., Jovanovic, M., Bokic, R., Kljajin, M., Karakašic, M. (2015). Theoretical and experimental analysis of elevator dynamic characteristics. Tehnički vjesnik / Technical Gazette, vol. 22, no. 4, p. 1011-1020, D0I:10.17559/TV-20150107175453. [14] Vladic, J., Bokic, R., Kljajin, M., Karakašic, M. (2011). Modelling and simulations of elevator dynamic behaviour. Tehnički vjesnik / Technical Gazette, vol. 18, no. 3, p. 423434. [15] Rao, S. S. (2011). Mechanical vibrations. Prentice Hall, Upper Saddle River. [16] Sinha, A. (2010). Vibration of Mechanical Systems. Cambridge University Press, Cambridge. [17] Vladic, J., Bokic, R., Jovanovic, M. (2019). Steel wire rope and computational-experimental procedures for the analysis of specific transport machines. University of Novi Sad, Faculty of Technical Sciences, Novi Sad. [18] Feyrer, K. (2007). Wire Ropes. University of Stuttgart, Stuttgart. [19] Goroshko, O.A. (2007). Evolution of the dynamic theory of hoist ropes. International Applied Mechanics, vol. 43, no. 1, p. 64-67, D0I:10.1007/s10778-007-0007-9. [20] Vladic, J., Malešev, P., Šostakov, R., Brkljač, N. (2008). Dynamic analysis of the load lifting mechanisms. Strojniški Vestnik - Journal of Mechanical Engineering, vol. 54, no. 10, p. 655-661. [21] Zhang, X., Li, H., Meng, G. (2008). Effect of Friction on the Slide Guide in an Elevator System. Journal of Physics: Conference Series, vol. 96, no. 012074, D0I:10.1088/1742-6596/96/1/012074. [22] Vladic, J., Bokic, R. (2017). Characteristics of mathematical methods and specialized software systems for dynamic analysis of elevators and mining elevators. IMK-14 - Research and development in Heavy Machinery, vol. 23, no. 2, p. 31-38, D0I:10.5937/imk1702031v. [23] Knezevic, B.Z., Blanusa, B., Marcetic, D.P. (2017). A synergistic method for vibration suppression of an elevator mechatronic system. Journal of Sound and Vibration, vol. 406, p. 29-50, D0I:10.1016/j.jsv.2017.06.006. Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators 301 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 © 2021 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2021.7168 Original Scientific Paper Received for review: 2021-03-17 Received revised form: 2021-05-05 Accepted for publication: 2021-05-19 Minimization of the Lifecycle Cost of a Rotary Heat Exchanger Used in Building Ventilation Systems in Cold Climates Ignas Sokolnikas1 - K^stutis Čiuprinskas2-* - Jolanta Čiuprinskiene2 i Salda UAB, Lithuania 2 Vilnius Gediminas Technical University, Lithuania This article presents an analysis of rotary heat exchangers (RHE) used as heat recovery units in building ventilation systems in cold climates. Usually heat exchangers with the highest heat transfer efficiency are the preferable option for this purpose. However, such exchangers usually have the highest media pressure drop, thus requiring the highest amount of energy for media transportation. In this study, the problem is solved by analysing the lifecycle cost (LCC) of the RHE including both the recovered heat and the electricity consumed in the fans of the air handling unit (AHU). The purpose of the investigation was to determine the optimal set of geometrical characteristics such as the exchanger's length, foil thickness, the height and width of the air channel. Two hundred and seventy different combinations were examined using analytical dependencies and ANSYS simulations. The results are compared with experimental data obtained earlier at the KOMFOVENT laboratory. The results show that the best overall energy efficiency is obtained in heat exchangers that do not offer the best heat recovery efficiency, and LCC differences in the same climatic and economic conditions can go as high as 31 %, mainly due to the geometrical parameters of the heat exchanger. Keywords: rotary heat exchanger, heat recovery, ventilation system, temperature efficiency, pressure loss, ANSYS, lifecycle cost Highlights • The energy efficiency of a heat exchanger should comprise not only the efficiency of heat recovery but also the energy needs for media transportation through the heat exchanger. • Heat exchangers with the highest heat recovery efficiency usually have the highest media pressure drop and, as such, the highest energy requirements for media transportation. • During the lifecycle of the heat exchanger, the best overall energy efficiency is obtained in heat exchangers with sub-par heat recovery efficiency. • During the lifecycle of the heat exchanger, the differences in the overall energy efficiency can reach 31 % due to the geometrical parameters of the heat exchanger • The best energy performance was obtained with the rotary heat exchanger that was the longest and had the smallest thickness of foil. 0 INTRODUCTION Nowadays, mechanical ventilation systems are an indispensable part of building engineering systems. Ventilation systems consume about 10 % of energy in commercial buildings [1]. Ventilation systems also increase heating and cooling system operating costs, which constitute 25 % and 9 % of building energy, respectively [1]. The main element of a ventilation system that saves the majority of thermal energy is the heat exchanger located in the air handling unit (AHU). Rotary and plate heat exchangers are most commonly used for this purpose. In Lithuania, as well as in other cold climate countries, rotary heat exchangers (RHE) are the preferred type due to lower frosting in low outdoor temperatures, humidity recovery, compactness, and no drainage requirement. Rotary heat exchangers usually are classified into three types: condensation, enthalpy, and sorption heat exchangers. In this research, a condensation-type rotary heat exchanger is investigated. The initial aim of the literature review was to discover some constructional recommendations in pursuance of maximal overall energy efficiency of the RHE, including heat transfer efficiency and the energy required to overcome the aerodynamic resistance of the RHE. As expected, no such specific recommendations were found. The majority of studies focus on the heat and moisture transfer efficiency of heat exchangers and ignore the pressure loss. The primary goal of research [2] was to develop a mathematical and numerical model to predict the energy wheel's effectiveness and to take into account the condensation and evaporation processes. Paper [3] presents the modelling results of the recovery unit working in real conditions with regard to changes in temperature, air velocity, and rotor speed. In [4], the heat and mass transfer in a desiccant wheel was modelled into a set of linear differential equations under the linearization assumptions on the temperature 302 *Corr. Author's Address: Vilnius Gediminas Technical University, Lithuania, kestutis.ciuprinskas@vilniustech.lt Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 and humidity profiles and the psychrometric relation. In [5], the optimum operational conditions of the rotary regenerator were obtained using genetic algorithm optimization technique subject to a list of constraints. The objective function in the optimization technique was the thermal effectiveness, while the design parameters (decision variables, the optimum operating conditions of the rotary regenerator were obtained using a genetic algorithm. The objective function in the optimization technique was the thermal effectiveness, while the design parameters (decision variables) were volumetric flow rates of cold and hot air streams, matrix rotational speed, and the exchanger's frontal area (heat transfer surface area). The purpose of the study in [6] was to propose a mathematical model for the heat and mass transfer characteristics (moisture removal capacity and the moisture removal regeneration) of a zeolite-coated heat exchanger. The lack of such recommendations has also been observed by RHE manufacturers [7]. In a study of geometrical parameters on heat and mass transfer processes [8], pressure drop over the RHE is used as an indicator of the frosting progress only; energy consumed by fans is not included in the evaluation. In some research projects, energy consumed in AHU fans is taken into account [9] and [10] and even is referred to as a "significant part" of the total primary energy used [10]; however, it still omitted from the overall energy efficiency evaluation. AHU fan energy is included in power usage effectiveness (PUE) used as an energy efficiency indicator of data centre facilities, and the study in [11] is quite informative in this context; however, it only gives a distinct presentation of electrical power, but not energy use, and no energy increase due to the use of the RHE is given. Regardless of the different field of application, one of the best formulations of the problem is given in an RHE optimization study [12]. Recovered heat and energy used in fans are defined as "conflicting objectives", and a solution is given as a Pareto optimum. In the field of air conditioning, the connection between heat transfer efficiency and pressure loss is probably given the most attention in-depth in two studies [13] and [14]. In [13], the optimal geometrical parameters of the RHE are determined using the Pareto front concept; in [14] the boundaries of a more generalized porosity factor are defined for the typical operational conditions of the RHE. Unfortunately, neither of the studies analyse the ratio between the recovered heat and the consumed electrical energy. Simplified analytical models are used in most of the sources reviewed, though it is well known that solving non-linear partial differential equations is generally needed in the case of heat and mass transfer. In [4], the advantages of analytical versus computational methods are presented, and an analytical solution using the linearization of differential equations of heat and mass transfer in a desiccant wheel is proposed. In contrast, some authors state that it is impossible to accurately calculate the temperature efficiency of an RHE by applying the analytical method, because the process of heat exchange is not stationary [15] and [16]. The application of numerical models produces results that are closer to experimental results yet are not always very reliable [2] and [17]. As expected, there is no clear dependence of temperature efficiency or pressure loss on one parameter. Due to the aforementioned reasons, there are no distinct guidelines for designing an RHE, or they are too abstract. Although some sources [13], [14] and [18] suggest that the RHE with the highest temperature efficiency ratio would not be expected to be the most energy-efficient due to the high-pressure losses of such a heat exchanger. Therefore, this research focuses on the optimal ratio of temperature efficiency to pressure losses. The objective of this research was to find the optimal set of geometrical RHE parameters, taking into consideration the heat recovered by the heat exchanger and the energy used for air transportation through the exchanger itself (i.e., the difference in the amount of electricity used in the fans with and without the RHE). These energy amounts are assessed during the lifecycle of the RHE. Additional outcomes from this research are the range of variation of lifecycle costs (LCC) and the possibility to use the results obtained in formulating the selection criteria of the RHE. 1 METHODS The computational fluid dynamics (CFD) model was adopted to determine the main variable of the problem: the temperature efficiency ratio of the RHE. The pressure drop was initially calculated using the same CFD model, which was later replaced with an analytical model due to its shorter calculation time. The CFD and analytical calculation results were validated by experimental tests results. The total amount of heat recovered, electricity consumed as well as energy and materials costs, were defined for Lithuanian weather and economic conditions. Minimization of the Lifecycle Cost of a Rotary Heat Exchanger Used in Building Ventilation Systems in Cold Climates 303 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 The set of the most significant parameters, such as the foil wave height and wavelength, the length of the RHE, the air velocity through the RHE, the foil thickness and the rotational speed of the RHE was determined on the basis of the literature review. where £2 is the temperature of the supplied air [°C]; t21 is the outdoor air temperature [°C]; t11 is the indoor air temperature [°C]. 1.2 CFD Model 1.1 Geometrical Model Constructively, an RHE is a drum (cylinder) of aluminium foil arranged in concentric rings. Both in the CFD model and in the experimental tests, the wave height and the wavelength of the RHE are determined using the Eurovent certification methodology [7]. The wavelength (b) was calculated using Eq. (1); the measurement scheme is shown in Fig. 1. Fig. 1. Measuring the wavelength of the rotary heat exchanger b = D ■ arcsin(c / D) (1) where b is the calculated wavelength; D is the diameter of the rotor [mm]; c is the length of the measured segment [mm]; n is the number of the waves in the measured segment. The wave height was calculated by measuring the height of 10 to 20 waves, which was then divided by the number of waves. The wave height is shown below (Fig. 2): SZESZI' Fig. 2. Measuring the wave height (a) of the rotary heat exchanger The temperature efficiency of the rotary heat exchanger was calculated using Eq. (2) [7]: ¿11 ¿21 (2) The temperature efficiency of the RHE is determined using the ANSYS FLUENT 18.0 finite volume-based software [19]. The CFD model of an RHE consists of the RHE itself and four airflow sections. The RHE model is divided into a finite volume network consisting of 2,039 nodes and 8,087 finite volumes. The finite volume network is built from tetrahedrons. The model and the finite volume network thereof are shown in Fig. 3. Fig. 3. The finite volume network of the rotary heat exchanger model To describe the conditions of the problem, the 1-epsilon viscosity model is selected. The change of the physical parameters of air over time is described in accordance with the polynomial functions. To describe the physical parameters of aluminium, the same data as in experimental tests are used (see below). Then the boundary conditions (the outdoor and indoor air temperatures, as well as air velocities) are determined. The characteristics of the rotary heat exchanger are introduced by applying the function of porous material. In this case, the porosity of the material c, the ratio between the surface area and volume A/V and the heat transfer rate hT are introduced [13]. A a = - AM + Ac A =- 2 a' = a - 2s, b' = b - 2s, (3) (4) (5) (6) (7) 304 Sokolnikas, I. - Čiuprinskas, K. - Čiuprinskiene, J. n Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 2b 2 , X2 3 + I -I ( b 1 ( a 1 \rna ) 1 - 1 4 -1 -ft 1 — ( 2b Y 12 ) 12 ) „ 4- . ft a, PL A _ v ~{AC + Am)L ' Nu ■ k hT = - D, (8) (9) (10) Nu = 1.1791 Deq = a 1 + 2.77011 ^J- 3.1901^-^7 +1.99751 - 0.4966i^ 1.0542 - 0.4670 — -0.1180 ( a b I b (11) +0.17941 - 0.0436 (â- (12) where a is the porosity of the RHE; 5 is the foil thickness, [m]; a' is the internal wave height of the channel, [m]; b' is the internal width of the channel, [m]; A m is the total area of foil cross-section, [m2]; Ac is the total area of air channel cross-sections, [m2]; P is the total perimeter of channels, [m2]; L is the length of the RHE, [m]; A is the surface area of foil, [m2]; V is the volume of the RHE, [m3]; hT is the heat transfer coefficient, [Wm-2K-1]; Nu is the Nusselt number; ka is the thermal conductivity of air, [Wm 'K-']; Deq is the equivalent diameter of the channel, [m]. First, the problem is approached in a stationary mode until the air velocity and temperature stabilizes while the rotor is not running. Then the non-stationary mode is activated, and the non-stationary problem of heat transfer is solved. The selected time-step is 0.25 s, while the maximum number of iterations during this time-step is 10. The process is halted when the variation in the temperature of the supplied air drops below 0.3 °C. This condition is reached within 150 s to 200 s of the process simulation, while calculations take 30 min to 40 min. 1.3 Calculation of Pressure Loss Unlike temperature efficiency, the pressure loss of an RHE can be calculated using the analytical method. Compared to numerical modelling methods, this method saves a significant amount of calculation time. The calculations were carried out using the following equations [13]: AP = 4 2 pu2 + 4/ pu2, uD Re = f a (13) (14) (15) fRe = 9.5687 1 + 0.07721 ^1 + 0.86191 a -0.83141^-1 - 0.2907 {Ç\ - 0.0338 M (16) where £,c is the factor for air expansion and compression at the entrance and exit of the rotor (0.2); p is the density of air under standard conditions (1.2 kgm-3); is the air velocity in the channel; fis the Fanning friction factor; Re is the Reynolds number; H is the kinematic viscosity [m2s-1]; vf is the air face velocity through the RHE, [ms1]. 1.4 Calculation of Life Cycle Costs The LCC are assumed as the sum of costs of the RHE, electricity consumed in the supply and exhaust air fans, and the heat consumed in the air heating coil, displaced after the RHE in the supply air stream. The lifetime of the RHE is assumed to be equal to the typical lifetime of AHU, which is 10 years. The cost of the RHE is assumed as the sum of two parts: the cost of manufacturing and the cost of aluminium foil. The manufacturing cost was fixed at 345 EUR, because the same size of all RHEs investigated. The cost of aluminium foil is variable because of the different lengths and porosity of the RHE. The price of aluminium foil was assumed to be 8.62 EUR per kg. The volume of aluminium foil consumed is calculated under Eq. (17): V = 0.25nD2 L (l-a). (17) The electrical power and energy consumed in the fans are calculated using the following equations: P = IO-3 qv Apfan / %, Eei = {{,sup + Pel,ex = = 2■ Pe, -10• 365-12• 5/7. (18) (19) where qv is the volumetric airflow; Apfan is the pressure difference of the fan, calculated by adding the 200 2 u = 2 2 2 Minimization of the Lifecycle Cost of a Rotary Heat Exchanger Used in Building Ventilation Systems in Cold Climates 305 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 Pa of the ventilation system aerodynamic resistance to the resistance of the RHE obtained from CFD calculations; ne is the efficiency of the fan assembly, assumed to be constant and equal to 0.625; top is the annual time of the fan's operation. It is assumed that the ventilation system is used 12 hours per each working day for 10 years. The heat consumed in the air heat coil is calculated using the following equations and assumptions: Ph = q Pcpmax tsa - ts,RHE - Atfan) , (20) W = toa + VsHE {L ~ *oa) - (21) where cp is the specific heat of the air (1 kJ kg-1K-1); tsa is the temperature of the supply air (20 °C); tea is the temperature of the exhaust air (20 °C); toa is the temperature of the outdoor air; tsRHE is the temperature of the air after the RHE; Atfan is the increase of the air temperature due to the supply air fan operation (1 °C); nRHE is temperature efficiency of the RHE obtained from CFD calculations. The heat energy calculations were performed for each hour of the ventilation system uptime. The outdoor air temperatures were obtained using RETscreen software [20]. As used for the purposes of the LCC calculations, the price of electricity was 0.099 EUR per kWh, and the price of the heat 0.0463 EUR per kWh. 1.5 Validation of the Models In order to compare the experimental and calculation results, tests with three different condensation rotary heat exchangers were performed at a laboratory located at the factory of KOMFOVENT, a company that specializes in the production of ventilation equipment. Over the course of the testing procedure, the temperature efficiency and pressure loss of the RHE was measured at air speeds between 0.5 ms-1 and 4 ms-1. The length of every RHE tested was 200 mm while the diameter was 1,000 mm, the rotational speed was 12 revolutions per minute and the thickness of the foil was 0.065 mm. The characteristics of aluminium foil used were as follows: density of 2,730 kgm-3, specific heat capacity of 900 Jkg-1K-1, thermal conductivity coefficient of 175 Wm-1K-1. The experiments were carried out under the following conditions: outdoor air temperature of 2 °C, relative humidity of 80 %, indoor air temperature of 22 °C, relative humidity of 45 %. As a result, the methodology for calculating the temperature efficiency and pressure loss of a rotary heat exchanger was validated. The comparison of the results for three different cases is provided below; Table 1 shows the following parameters of the heat exchangers under comparison: the air velocity during the test, the wave height, and the wavelength of the foil. As shown in Table 1, the calculation results are rather similar to the experimental results. In terms of temperature efficiency, the most significant relative difference compared to the experimental result was 2.54 %. As for pressure loss, the maximum relative difference was even smaller, only 1.19 %. Minor differences between experimental and calculation results most likely occur due to calculation errors for the wavelength and wave height of the rotary heat exchanger since these parameters vary across the entire area of the rotor. In addition to these aspects, there might be other reasons that cause the aforementioned differences: the accuracy of the instruments used to measure temperature, the pressure and airflow rate and the excessive sparsity of the finite volume network. However, a 2.54 % margin of calculation error for temperature efficiency and a 1.19 % margin of calculation error for pressure loss of a rotary heat exchanger constitutes a very good result. This confirms that appropriate methods were selected to determine the temperature efficiency and pressure loss of a rotary heat exchanger. 2 RESEARCH SCOPE In the search for the optimal set of RHE parameters, establishing the scope of research is key since there Table 1. Results of experiments and calculations _Parameters_Experiment_Calculations_Absolute difference_Relative difference v, a, b, Temperature Pressure Temperature Pressure Temperature Pressure Temperature Pressure [ms-1] [mm] [mm] efficiency, [%] loss, [Pa] efficiency, [%] loss, [Pa] efficiency, [%] loss, [Pa] efficiency, [%] loss, [%] 2.5 1.69 3.85 79.36 129.8 78.61 129.29 0.75 0.47 0.95 0.36 2.5 1.39 2.60 85.17 234.0 83.01 235.22 2.16_-1.22_2.54_-0.52 1.5 1.39 4.04 84.82 102.5 82.77 103.70 2.05 -1.22 2.42 -1.19 306 Sokolnikas, I. - Čiuprinskas, K. - Čiuprinskiene, J. Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 are many variables, and the number of possible combinations thus grows exponentially. The time for calculating temperature efficiency using the ANSYS FLUENT software was reduced as much as possible. The primary calculations took more than 24 hours; however, the time required for the calculations was reduced from 30 min to 40 min (computer specifications: i5-5200U CPU @ 2.2 GHz, 16 GB RAM). This was achieved by shortening the airflow ducts from 3 m to 1 m, reducing the number of finite volumes from 105,867 to 8,087, increasing the time-step from 0.05 s to 0.25 s and expanding the conditions for the solution from At = 0.1 °C to At = 0.3 °C. The consistency between the calculation results and the experimental data shows (see Table 1) that these changes are acceptable. The literature review has shown that once the rotational speed of the rotary heat exchanger reaches a certain value, there is virtually no effect on the temperature efficiency [2], [5] and [13]. The results of a large number of studies show that this value does not exceed 12 revolutions per minute. Therefore, in this study, taking into consideration the duration of calculations, the rotational speed of all rotary heat exchangers was pegged at 12 revolutions per minute. During this study, in search for the optimal set of RHE parameters, air velocity was considered a constant, primarily due to the fact that in almost all the cases of preliminary calculations, the best results were achieved at the lowest air velocity. However, this imposes the largest dimensions of the AHU, which quickly becomes unacceptable from the practical point of view. The analysis of the heat exchangers tested at the KOMFOVENT laboratory showed that the highest temperature efficiency is achieved when the air velocity is about 1.5 ms-1. Therefore, in this study, the air velocity in all analysed rotary heat exchangers was the same: 1.5 ms-1. Having established the rotational speed and air velocity of the RHE as immutable values, this research focuses on the length of the RHE, the length and width of the foil wave and foil thickness. A matrix of parameter values was made considering the time required for model calculation. These values are provided in Table 2. Table 2. Variable parameters of rotary heat exchangers Parameter Values of the parameter L, [mm] 200 300 400 - - - 5, [mm] 0.06 0.08 0.10 - - - a, [mm] 1.4 1.5 1.6 1.7 1.8 - b, [mm] 2.5 3.0 3.5 4.0 4.5 5.0 Two hundred and seventy different combinations of rotary heat exchanger parameters were simulated and compared in total. The heat exchangers were compared on the basis of LCC. 3 RESULTS AND DISCUSSION The LCC of different RHE cases are provided below (Fig. 4). A heat exchanger with a length of 400 mm, foil thickness of 0.06 mm, wave height of 1.8 mm, wavelength of 5 mm, pressure loss of 116.5 Pa and Case number Fig. 4. The life-cycle costs of different rotary heat exchanger options (each vertical line represents a separate case simulated) Minimization of the Lifecycle Cost of a Rotary Heat Exchanger Used in Building Ventilation Systems in Cold Climates 307 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 temperature efficiency of 88.2 % produces the lowest lifecycle costs. The lifecycle costs of this kind of heat exchanger are 3,596 EUR over 10 years, i.e., 31 % lower than the heat exchanger with the highest lifecycle costs of 5,204 EUR over 10 years. As seen in Fig. 5, there is a certain, visually quite obvious correlation (R = 0.50) between the LCC and temperature efficiency. It can be stated that the lowest LCC is typical of heat exchangers with a temperature efficiency between 85 % and 90 %. The figure also shows that both an excessively low and excessively high-temperature efficiency negatively affects the LCC. In this case, when the temperature efficiency is very high, the pressure loss is also significant, leading to increased LCC. Fig. 6 shows the dependence of the LCC on pressure loss. A quite obvious correlation (R = 0.52) is evident here as well. In this case, to reduce the 5400 S200 5000 4800 ^ 4600 en 3 Si 4400 (J u ^ 4200 4000 3800 3600 3400 * _m_# • m / ► / • • • t m • »v •• i v- \ L" / • / • • vV v •vi? J !. » • / • T •• i • H • • • w • • • • f »v hJ ■ * __ri R! = 0,4969 • /.v S3 & MS** rj» • • ••vv • 75 77 79 81 91 93 83 85 87 89 Temperature efficiency [%] Fig. 5. The dependence of the LCC on temperature efficiency (each dot represents a separate case simulated) 95 5400 5200 5000 4800 4600 Si 4400 4200 4000 3800 3500 3400 • t v V/.1. • . ™ 4* M__* ■ ..... a i ■ v • • • • r *V ' - m * W 50 100 150 200 i k 250 300 350 ; [Pa] Pressure I Fig. 6. The dependence of the LCC on pressure loss (each dot represents a separate case simulated) 308 Sokolnikas, I. - Čiuprinskas, K. - Čiuprinskiene, J. Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 LCC, the pressure loss of the heat exchanger should be within 100 Pa and 200 Pa. This chart, much like the one above, shows that both an excessively low and excessively high-temperature efficiency leads to increased LCC. Pressure loss that is too low leads to lower efficiency, which prevents more significant energy savings, while high-pressure loss leads to increased electricity costs. In order to determine the optimal parameters of rotary heat exchangers, all heat exchangers with the LCC below the considered limit of EUR 3,800 were analysed (i.e., a total of 22). The length of the first five heat exchangers with the lowest LCC is 400 mm. It should also be noted that the foil thickness of all these heat exchangers is 0.06 mm, which further tightens the limits of the pressure loss that should be taken into consideration when producing a heat exchanger that is the most appealing in terms of the LCC, i.e., 101 Pa to 179 Pa. In this case, the temperature efficiency of heat exchangers is between 85 % and 90 %. Table 3. Optimal parameters of heat exchangers The length of a RHE [mm] The number of heat exchangers of respective length [-] b/a [-] s [mm] a [mm] b [mm] 200 2 1.73 0.060 1.45 2.50 300 11 1.86 0.060 1.71 3.18 400 9 2.59 0.060 1.72 4.44 Table 3 lists the average parameters of heat exchangers. It becomes evident that the most typical lengths of heat exchangers with the lowest LCC are 300 mm and 400 mm (11 and 9 units, respectively). Only two heat exchangers are 200 mm long. The table also shows that the foil of these heat exchangers is the thinnest out of those analysed in this study: 0.06 mm. The average ratio of wavelength to wave height, in this case, is 1.73, 1.86, and 2.59 for the heat exchangers with the length of 200 mm, 300 mm, and 400 mm, respectively. The average wave height of heat exchangers that are the most appealing in terms of the LCC is 1.45 mm (for 200 mm heat exchangers), 1.71 mm (for 300 mm heat exchangers), and 1.72 mm (for 400 mm heat exchangers). The average wavelength of heat exchangers is 2.50 mm for 200 mm, 3.18 mm for 300 mm, and 4.44 mm for 400 mm. Neither the wave height and wavelength nor the ratio of the two have been observed to have any impact. This can be explained by the narrow gauge of the air channels. While they are narrow enough to keep the physical proprieties of the air stream homogenous in the same cross-section, the shape of the air channel's cross-section makes no difference. Obviously, the total foil surface area exposed to the air has a much greater influence. This parameter affects the RHE's temperature efficiency and its pressure drop. 5 CONCLUSIONS 1. The literature review reveals a lack of specific recommendations for the constructional parameters of rotary heat exchangers. Most research projects focus on simulating heat and mass transfer processes and on improving the efficiency of these processes. There is also a lack of discussion about the criteria of rotary heat exchanger optimization. Most cases only refer to the efficiency of heat and moisture recovery as a criterion. 2. Good coherence of the calculation results and laboratory tests made at the KONFOVENT laboratory proves the suitability of the CFD model for the calculation of the temperature efficiency of rotary heat exchanger described. These tests also prove the suitability of the analytical model presented in [13] for the calculation of pressure drop of a rotary heat exchanger. 3. The lifecycle costs analysis of 270 different RHE variants calculated using the CFD and analytical models has shown that the best results are achieved when the length of the rotary heat exchanger is maximal (400 mm in the cases analysed), foil thickness is minimal (0.06 mm in the cases analysed), and pressure loss is between 100 Pa and 180 Pa. The temperature efficiency of an RHE in cold climate conditions, in this case, is expected to be between 85 % and 90 %. 4. This research confirms that it does not take the best heat recovery efficiency to achieve the best energy efficiency of a rotary heat exchanger due to its pressure loss. The difference between the heat exchangers with the lowest LCC and the heat exchangers with the highest LCC can reach 31 % in the same weather and economic conditions; this difference depends mainly on the geometrical characteristics of the rotary heat exchanger. 8 REFERENCES [1] EIA (2012). Commercial Buildings Energy Consumption Survey: Energy Usage Summary., from https://www.eia.gov/ consumption/commercial/reports/2012/energyusage/, accessed on 2021-01-22. [2] Al-Ghamdi, A.S. (2006). Analysis of Air-to-Air Rotary Energy Wheels, PhD Thesis, Ohio University, Ohio. Minimization of the Lifecycle Cost of a Rotary Heat Exchanger Used in Building Ventilation Systems in Cold Climates 309 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 302-310 [3] Grzebielec, A., Rusowicz, A., Rucinski, A. (2014). Analysis of the performance of the rotary heat exchanger in the real ventilation systems, The 9th International Conference Environmental Engineering, D0I:10.3846/enviro.2014.259. [4] Kang, H., Lee, G., Lee, D.-Y. (2015). Explicit analytic solution for heat and mass transfer in a desiccant wheel using a simplified model. Energy, 93, p. 2559-2567, D0I:10.1016/j. energy.2015.10.091. [5] Sanaye, S., Jafari, S., Ghaebi, H. (2008). Optimum operational conditions of a rotary regenerator using genetic algorithm. Energy and Buildings, vol. 40, no. 9, p. 1637-1642, D0I:10.1016/j.enbuild.2008.02.025. [6] Park, B., Lee, S. (2020). Investigation on heat and mass transfer characteristics for a zeolite-coated heat exchanger using comparatively low-temperature energy: heating humidification mode and cooling dehumidification mode. Indoor and Built Environment, p. 1-17, in press, D0I:10.1177/1420326X20942291. [7] Eurovent (2021). Air to Air Rotary Heat Exchangers. Certification Programme, Eurovent, Paris. [8] Kanas, P., Jedlikowski, A., Anisimov, S. (2019). the influence of geometrical parameters on heat and mass transfer processes in rotary heat exchangers. SN Applied Sciences, vol. 1, art ID 526, D0I:10.1007/s42452-019-0540-2. [9] Merzkirch, A., Maas, S., Scholzen, F., Waldmann, D. (2016). Field tests of centralized and decentralized ventilation units in residential buildings - specific fan power, heat recovery efficiency, shortcuts and volume flow unbalances, Energy and Buildings, vol. 116, p. 376-383, D0I:10.1016/j. enbuild.2015.12.008. [10] Dodoo, A. (2020). Primary energy and economic implications of ventilation heat recovery for a multi-family building in a nordic climate. Journal of Building Engineering, vol. 31, art. ID 101391, D0I:10.1016/j.jobe.2020.101391. [11] Mok, S., Kumar, S., Hutchins, R.R., Joshi, Y.K. (2016). Impact of a rotary regenerative heat exchanger on energy efficiency of an air cooled data center. 15th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems, p. 1182-1190, D0I:10.1109/ITHERM.2016.7517682. [12] Wang, L., Bu, Y., Chen, X., Wei, X., LI, D., Che, D. (2018). Multi-objective design optimization of rotating regenerative air preheater using genetic algorithm. Proceedings of the ASME Power Conference with ASME 12th International Conference on Energy Sustainability and ASME Nuclear Forum. DOI:10.1115/POWER2018-7417. [13] De Antonellis, S., Intini, M., Joppolo, C.M., Leone, C. (2014). Design optimization of heat wheels for energy recovery in HVAC systems. Energies, vol. 7, no. 11, p. 7348-7367, D0l:10.3390/en7117348. [14] Mioralli, P.C., Da Silva, M.A.B., Avallone, E., Palota, P.H., Natividade, P.S.G. (2019). Computational analysis for good thermal exchange and low pressure drop in regenerative air preheaters. Journal of Energy Research and Reviews, vol. 2, no. 4, p. 1-15, D0I:10.9734/jenrr/2019/v2i430083. [15] Fathieh, F., Besant, R.W., Evitts, R.W., Simonson, C.J. (2015). Determination of air-to-air heat wheel sensible effectiveness using temperature step change data. International Journal of Heat and Mass Transfer, vol. 87, p. 312-326, D0I:10.1016/j. ijheatmasstransfer.2015.04.028. [16] Yilmaz, T., Buyukalaca, O. (2003). Design of regenerative heat exchangers. Heat Transfer Engineering, vol. 24, no. 4, p. 3238, DOI:10.1080/01457630304034. [17] Zendehboudi, A. (2016). Implementation of GA-LSSVM modelling approach for estimating the performance of solid desiccant wheels. Energy Conversion and Management, vol. 127, p. 245-255, D0I:10.1016/j.enconman.2016.08.070. [18] Rashidi, S., Kashefi, M.H., Kim, K.C., Samimi-Abianeh, O. (2019). Potentials of porous materials for energy management in heat exchangers - A comprehensive review. Applied Energy, vol. 243, p. 206-232, D0I:10.1016/j.apenergy.2019.03.200. [19] Ansys Fluent | Fluid Simulation Software (2021). from https:// www.ansys.com/products/fluids/ansys-fluent, accessed on 2021-05-03. [20] Government of Canada, (2010). Retscreen. Clean Energy Management Software, from https://www.nrcan.gc.ca/maps-tools-and-publications/tools/modelling-tools/retscreen/7465, accessed on 2021-05-03. 310 Sokolnikas, I. - Čiuprinskas, K. - Čiuprinskiene, J. Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 © 2021 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2021.7203 Original Scientific Paper Received for review: 2021-03-31 Received revised form: 2021-05-16 Accepted for publication: 2021-05-21 Process Parameters Optimization for Maximizing Tensile Strength in Friction Stir-Welded Carbon Steel Anmol Shada1,2^ - Reeta Wattal2 JThe NorthCap University, Department of Mechanical Engineering, India 2Delhi Technological University, Department of Mechanical Engineering, India The present study focuses on improving the ultimate tensile strength of friction stir welded carbon steel (AISI1018). The effect of the process parameters (welding speed, tool RPM, and shoulder diameter) on the response parameters (ultimate tensile strength, percentage elongation and percentage reduction in area) were studied. Response surface methodology was used to develop the mathematical model for response parameters, and the adequacy of the model was checked using analysis of variance (ANOVA). The welding speed and tool RPM were found to affect the ultimate tensile strength significantly. The percentage elongation was affected only by welding speed. The percentage reduction in the area was affected by welding speed and shoulder diameter. The microstructure and microhardness of the weld have been studied and reported in the study. Keywords: friction stir welded steel, ANOVA, response surface methodology, ultimate tensile strength, microstructure and microhardness Highlights • In the present study, experimental analysis was done on friction stir-welded AISI1018 carbon steel joints to determine the effect of process parameters on response parameters (i.e., ultimate tensile strength, percentage elongation and percentage reduction in area) by developing mathematical models using the response surface methodology. • The developed regression models proved to be adequate for determining the output response on friction stir-welded steel joints at a 95 % confidence interval. • The surface plots are also presented in the study, which illustrates the individual effect of process parameters on the mechanical properties of the friction stir-welded steel joint. • The microstructure and microhardness of the friction stir-welded steel joint under optimized condition are also presented, and the observable pattern is presented in detail. 0 INTRODUCTION Steel is the most popular choice of engineers because it has a good combination of mechanical properties and is the most widely used material for industrial applications. Steel is generally welded with a fusion-welding process that leads to various problems, including porosity, hot cracking, hydrogen embrittlement, microstructural changes in the heat-affected zone, etc. [1]. The reason for the occurrence of these defects is mainly due to the inherent melting of base metal. In order to overcome these defects, the friction stir-welding process is considered the best choice as it does not involve the melting of workpiece; thus, the chance of defects arising due to melting are eliminated. The process was patented by The Welding Institute, U.K. [2]. The setup consists of a non-consumable tool having a shoulder through which a small pin is extended. The plates that are to be welded are rigidly fixed, and the rotating tool is plunged within the joint of plates to the extent the shoulder touches the plate surface, resulting in the generation of heat. This localized heating results in softening the material around the pin. The material moves from the front to the back of the pin because of the rotational and translational movements of the tool; in this manner, a solid-state joint is produced [3]. The initial work on friction stir welding was done on aluminium and its alloys, as they have a low melting point [4] to [6]. Many aluminium alloys like Al 7075, Al 6063 [7] to [9] were the popular choice of the researchers. Then, the microstructural examinations were carried out on those alloys [10] and [11]. Later, the process parameter optimization studies were reported on friction stir-welded aluminium and its alloys [12] and [13]; other studies including corrosion analysis and tensile strength optimization were also reported for aluminium [14] to [16] and magnesium alloys [17] but nowadays the research focus is on high-temperature materials. A feasibility study related to friction stir welding of mild steel was initially presented by Linert et al. [18], in which defect-free welds were successfully produced, but excessive tool wear was reported; the paper emphasized the need to work on the geometry of the tool in order to eliminate the issue of tool wear. The study demonstrated that the friction stir-welding of mild steel is feasible *Corr. Author's Address: Delhi Technological University, Department of Mechanical Engineering, Delhi, India, anmolbhatla@ncuindia.edu 311 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 without the loss of tensile properties. Ueji et al. [19] butt-welded ultra-fine plain low carbon steel with the help of a tungsten carbide tool. The paper reported the microstructural characterization and microhardness of the joint, but no mechanical testing was performed on the weld. Thereafter, a study on friction stir-welded 1012 and 1035 carbon steels was conducted by Fujii et al. [20] using a tungsten carbide tool; the paper presented the effect of welding speed on the ultimate tensile strength of the weld nugget; mechanical properties and microstructural examination were also presented. An effort was then made to weld high carbon steel using FSW by Cui et al. [21]; the study was conducted on AISI 1070 material using a tungsten carbide tool. Metallurgical examinations were reported in the paper with an enhanced joint strength than the base metal. The effect on friction stir-welding conditions on the properties and microstructure was reported by Miles et al. [22]; the paper concluded that desired mechanical properties can be achieved by adjusting the process parameters. Later, a study was conducted by Ghosh et al. [23] in which the effect of process parameters was analysed on friction stir-welded medium carbon steel; the paper concluded that the strength of the stir zone was higher than the parent metal, and HAZ was reported to be the weakest section. Also, the ductility of the weld was found to be lower than the base metal. Lakshminarayanan et al. [24] conducted friction stir welding on 5 mm thick AISI 1018 plates using a tungsten carbide tool and reported an increased tensile strength value as compared to the base metal. The microstructural examination was also reported in the paper in which smaller pearlitic structures were observed in the weld region. Ghosh et al. [25] performed experimental analysis on lap-welded friction stir-welded M190 steel sheets. The variable parameters were kept as rotational and transverse speeds with the range of 600 mm/s to 1200 mm/s and 0.85 mm/s to 3.39 mm/s, respectively. The microstructural examination revealed three regions: the nugget and two heat affected zones (HAZ-1 and HAZ-2) in which HAZ-1 was reported the weakest region constituting ferrite pearlite structures, while HAZ-2 contains martensite along with the increased value of hardness. The authors also reported that the best process parameter window to obtain joint efficiency of more than 60 % was 1000 rpm tool rotational speed and 1.69 to 3.39 transverse speed. Yabuuchi et al. [26] analysed the effect of rotational speed on the microstructure and tensile strength of friction stir-welded oxide dispersion-strengthened steel. The investigations were carried out on the speed ranging from 250 rpm to 400 rpm. The study revealed that the hardness was reduced when welding was carried out at all the speed because of the occurrence of recrystallisation. The change in ultimate tensile strength was also reported similar to the change in hardness, meaning that the joint strength was also reduced by grain growth. Karami et al. [27] analysed the effect of the tool rotational speed and transverse speed on microstructure and tensile behaviour of friction stir welded mild steel. A tunnel defect was observed in the specimens with lower welding speed or higher rotational speeds as the amount of heat input was not sufficient to provide enough flowability. The microstructural study indicated that the stir zone or HAZ had an austenitic structure and, on cooling, transforms to ferritic and pearlitic structures. The study also revealed that all the specimens exhibit high yield strength and a lower value of uniform elongation because of the fine microstructure in stir zone and heat affected zone of the welded specimen. The microstructural and tensile behaviour of friction stir-welded trip steel were studied by Mironov et al. [28], who reported that the thermal effect of welding gave rise to thermal softening of the material in HAZ and that martensite formation occurred in the stir zone. A premature failure was also reported during the tensile test, because of the non-homogeneous strain distribution that occurred because of microstructural changes. Mahmoudiniya et al. [29] applied a friction stir-welding technique to weld DP700 steel at different rotational speeds 600 rpm, 800 rpm, and 1000 rpm. The authors evaluated the microstructural and mechanical properties of the weld; they reported the occurrence of ferrite bands in the stir zone at 600 rpm, which lead to reduced strength. For the condition of 800 rpm and 1000 rpm, joint failure occurred at softened HAZ. The maximum joint strength was reported for the rotational speed of 800 rpm. Lee et al. [30] studied the friction stir welded medium Mn steel using a tungsten carbide tool. The paper reported an increase in ultimate tensile strength and total elongation because of the increased strain hardening rate. Although few studies were conducted to analyse the effect of process parameters on microstructure and mechanical properties for friction stir-welded steels, there is a gap in understanding the mechanical properties and concerned microstructural characterization of the friction stir-welded mild steel. Therefore, an attempt is made to establish empirical relationships to determine the effect of tensile strength, percentage elongation and percentage reduction in area for a friction stir-welded AISI 1018 312 Bhatia, A. - Wattal, R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 joint. Analysis of variance was also incorporated to check the adequacy of the developed model. Thus, the objective of the present study is to investigate the effects of process parameters on output response parameters such as tensile strength, percentage elongation and percentage reduction in area. RSM is also incorporated to optimize the process parameters for tensile strength. 1 MATERIAL AND METHODS A butt-welded joint was produced on friction stir-welding machine, which is developed by R.V. Machine Tools, India. Table 1 shows the specifications of the machine. The material chosen for the experimental analysis was AISI 1018, with a thickness of 3 mm. The steel plate was 200 mm in length, 40 mm in breadth, and 3 mm thick. The composition and mechanical properties of the material is illustrated in Tables 2a and b. The friction stir-welding was done with a tungsten carbide tool with 7 % cobalt. The sectional view of the friction stir-welded plate with the start and stop position of the tool is illustrated in Fig. 1. The composition of the tool used in the present study was made with EDAX analysis and is presented in Fig. 2a. Fig. 1. Photograph of friction stir-welded joint The dimensions of the tool used are presented in Fig. 2b. The profile of the tool pin used is presented in Fig. 2c. The length of the pin was kept as 2.6 mm so that its proper penetration takes place in the plates. A 1.5-degree tool tilt angle was kept throughout the experimentation, as shown in Fig. 2d. From the literature review, the process parameters chosen for the present study were welding speed, tool rpm, and shoulder diameter. These parameters are the major contributors in generating heat and affecting the tensile strength. The selected process parameters and their levels are tabulated in Table 3. The range of the selected parameters was determined by conducting several trials of input parameters in which obtaining the defect-free weld was the key factor. This was carried out by changing one of factors from minimum to maximum, keeping the other factors as constant. The observations from the pilot experimentation are presented in Table 4. Table 1. Specifications of friction stir welding machine used for experimentation Spindle motor 11 kW, 1440 rpm, 440 V 3-Phase, AC Drive, Flange mounting (STARK) Spindle speed 1440 rpm (max) Spindle housing tilting Angle (-5° to +5°) Spindle pulley type Timing Pulley Z-axis stroke 300 mm Z-axis thrust 3922.6 N (min)- 39226.6 N (max) _(Adjustable in steps)_ Z-axis feed rate 0 mm/min to 2000 mm/min X-axis stroke 600 mm Table 600 mm x 400 mm 'T' slot- 18 mm x 3 mm Tool Holder ISO 40 Taper-Side lock holder Z-axis feed force feedback Load cells of capacity 5000 kgs, with least count 1 kg is provided Controller PC based control system Table 2a. Chemical composition of base metal (results by spectrometry) Material C [wt%] S [wt%] Mn [wt%] P [wt%] Fe [wt%] AISI 1018 Steel 0.18 0.03 0.65 0.04 Bal Table 2b. Mechanical properties of base metal U.T.S. U.T.S. for reduced Elongation Yield strength Reduction [MPa] section [MPa] [%] [MPa] in area [%] 416 435 21.8 280 50.77 Table 3. Chosen process parameters with their levels Parameter chosen Notation Units - Levels 1 2 3 Welding speed V mm/min 60 110 210 Tool rpm N mm/min 430 550 750 Shoulder diameter D mm 15 18 20 2 EXPERIMENTATION The L9 orthogonal array was used as a design matrix for the conduction of experiments. Transverse tensile test specimens were prepared as per the E-8 standard of ASTM [31]. Two tensile tests were conducted on the welds: transverse tensile testing and reduced section tensile testing. The former is used to locate the region of lowest strength in weldment and to determine the strength in that region; however, when the strength of weld nugget is to be determined then the later one that is reduced section tensile test is used in which area of the cross-section is deliberately reduced at Process Parameters Optimization for Maximizing Tensile Strength in Friction Stir-Welded Carbon Steel 313 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 0 2 4 6 keV B 10 12 Fig. 2. a) Results obtained from EDAX analysis of tool, b) dimensions of the tool, c) tool profile, and d) 1.5-degree tool tilt angle used in generating the weld samples Table 4. Observations made while selecting the range of process parameters Variable parameter Input parameters Macrostructure Observable remarks Welding speed Welding speed < 60 mm/min Tool RPM = 430 mm/min Shoulder diameter = 20 mm Root stuck to the surface because of increased contact time of hot metal and backing plate. Welding speed > 210 mm/min Tool RPM = 430 mm/min Shoulder diameter = 20 mm Groove defect was observed on the weld surface. This is because of the hindrance between the frictional heat and the material flow. Tool RPM Tool RPM < 430 mm/min Welding speed = 60 mm/min Shoulder diameter = 20 mm Excessive degradation of the tool pin was observed. This is due to the breakage of tool pin because of excessive heat generation. Tool RPM >750 mm/min Welding speed = 60 mm/min Shoulder diameter = 20 mm Tunnel defect with excessive flash was observed. The excessive flash occurred because of high rotational speed. Shoulder diameter Shoulder diameter < 15 mm Welding speed = 60 mm/min Tool RPM = 430 mm/min Surface defect was observed. This is because the area between the shoulder pin and shoulder diameter was reduced, which ultimately reduces the flow of material. Shoulder diameter > 20 mm Welding speed = 60 mm/min Tool RPM = 430 mm/min Tunnel defect was observed on the welded specimen. This is because excessive heat was generated and material flow was inadequate. 314 Bhatia, A. - Wattal, R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 weld nugget; this ensures that the transverse tensile test specimen fractures at the weld nugget region although weld nugget may not be the weakest region of the weldment. Three tensile test specimens were prepared for each welding condition. The samples for tensile test, microhardness, and microstructure evaluation were extracted from the welded specimen, as illustrated in Fig. 3 a. The dimensions of the transverse tensile test specimen and reduced section tensile test specimen are shown in Fig. 3b and Fig. 3 c, respectively. The tensile testing of the specimen was done with a Zwick/Roell Z250 computer-controlled machine. The speed of testing was restricted between 1.5 mm/ min to 15 mm/min according to the E-8 standard of ASTM; the average of the results of ultimate percentage elongation and percentage reduction in area for three specimens are recorded and tabulated in Table 5. The microstructure of the FSW joint for the optimized condition was observed with a scanning electron microscope. The specimen extracted for microstructural evaluation was further polished and etched with a 3 % Nital solution. Microhardness testing was done as per ASTM E-92. Vickers's hardness at a 100-gram load was taken over the wide region of the welded specimen under optimized conditions. The fractured tensile tested specimen was also observed under a scanning electron microscope to observe the fractured surface morphology. 3 RESULTS AND DISCUSSIONS 3.1 Developing Mathematical Model Tensile strength is a function of process parameters such as welding speed (V), tool rpm (N) and shoulder diameter (D). As the number of levels of the above process parameters was three, so up to second-order terms have been considered in the regression equation. The second-order regression equation used to represent the effect is represented by Eq. (1). Effect = Constant + b1{V) + b2(V)2 + b3(N)2 + b4(N)2 +b5(D + be(D)2, (1) where b1, b2, b3, b5 and bg are the regression coefficients. During the regression analysis, some of the coefficients may come out to be zero, meaning that those parameters do not affect the particular mechanical property. Standard t-test and regression Fig. 3. a) Test coupon for taking out samples for mechanical and metallurgical testing, b) dimensions of transverse tensile test specimen, and c) dimensions of reduced section tensile test specimen Table 5. Run order and corresponding output responses Run Coded values Original values For transverse tensile test specimen For reduced section specimen order V N D V N D UTS [MPa] UTS [MPa] Elongation [%] Reduction in area [%] 1 1 1 1 60 430 15 407.13 538.73 20.49 44.74 2 1 2 2 60 550 18 430.00 460.70 20.62 45.72 3 1 3 3 60 750 20 414.10 526.60 18.84 42.48 4 2 1 2 110 430 18 408.83 533.56 22.00 58.82 5 2 2 3 110 550 20 410.23 525.96 23.26 52.30 6 2 3 1 110 750 15 408.53 549.76 23.82 47.93 7 3 1 3 210 430 20 408.93 556.80 25.53 60.34 8 3 2 1 210 550 15 410.60 525.16 24.29 64.55 9 3 3 2 210 750 18 419.56 564.36 24.09 65.41 Process Parameters Optimization for Maximizing Tensile Strength in Friction Stir-Welded Carbon Steel 315 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 analysis was performed using MINITAB software. F-test was also performed on the data at a 95 % confidence level to prove that the data fit the proposed model. The strength of the weld nugget could only be determined by using the reduced section tensile test specimen. The final mathematical model developed for reduced section specimen (i.e., ultimate tensile strength [MPa], elongation [%] and reduction in area [%]) is presented below: UTS = 624.1 + 20.05(V) - 162.2(N) + 41(A2), (2) Elongation = 18.32 + 2.326(V), (3) Reduction in area = 19.86 + 9.56(V) + 18.02(D) - 4.59 (D2). (4) Analysis of variance was done to identify the significant factors that affect the ultimate tensile strength, percentage elongation, and percentage reduction in area. The results of ANOVA are presented in Table 6. 3.2 Checking the Adequacy of Developed Model To check the adequacy of the developed model, the coefficient of determination (RE2) and the square root of the coefficient of determination (RE) were determined. The closer is the value of RE to one, the better is the model. The sum of the square of residual (SS residual), regression model sum of squares (SS model), mean square of residual (MS residual), and the degree of freedom (Df model and Df residual) are used for the calculation of the coefficient of determination (RE2). Adjusted RE2 that adjusts the number of terms in a model was also calculated. The correlation parameters for the developed model for the ultimate tensile strength, percentage elongation, and percentage reduction in area were calculated and are presented in Table 7. The 'p' value for all the developed models is less than 0.05, which proves that the model is adequate. Table 6a. ANOVA table for Ultimate tensile strength of reduced section specimen Effect Error /-value /-value Significant Constant 48.0293 12.5563 0.0000 Yes Welding speed 8.20525 3.7432 0.0406 Yes (Welding speed)2 Pooled Tool RPM 53.5256 -3.2989 0.0128 Yes (Tool RPM)2 13.2450 3.23889 0.0208 Yes Shoulder diameter Pooled No (Shoulder diameter)2 Pooled No Table 6b. ANOVA table for percentage elongation of the weld Effect Error /-value /-value Significant Constant 0.8655 21.88 0.00000 Yes Welding speed 0.4021 5.334 0.00092 Yes (Welding speed)2 Pooled No Tool RPM Pooled No (Tool RPM)2 Pooled No Shoulder diameter Pooled No (Shoulder diameter) ,2 Pooled No Table 6c. ANOVA table for percentage reduction in area of the weld Effect Error /-value /-value Significant Constant 6.993 2.5673 0.028756 Yes Welding speed 1.028 8.7654 0.000294 Yes (Welding speed)2 Pooled No Tool RPM Pooled No (Tool RPM)2 Pooled No Shoulder diameter 7.567 2.65431 0.041513 Yes (Shoulder diameter) i2 1.889 -2.65435 0.034563 Yes Table 7. Correlation parameters of the developed model Correlation parameter For UTS For percentage For percentage Elongation reduction in area Multiple Re 0.8938 0.97423 0.97231 Multiple RE2 0.8038 0.8518 0.93831 Adjusted RE2 0.6860 0.8024 0.90132 SS model 5801.13 32.7336 590.971 Df model 1 1 3 MS model 1933.71 31.9065 196.990 SS residual 1530.286 7.1123 38.870 Df residual 5 7 5 MS residual 308.6411 1.0349 7.7741 F 6.83 29.861 25.34 P 0.032 0.00094 0.002 3.3 Optimization of Welding Parameters In order to optimize the welding parameters, response surface methodology is used. Figs. 4 to 6 represent the surface plots and contour plots for UTS, percentage elongation, and percentage reduction in area for the different combination of process parameters. 3.3.1 Effect of Process Parameters on UTS of Welds The temperature generated during welding is inversely proportional to welding speed. Thus, recrystallisation takes place at a lower temperature when welding is carried out at a higher speed, which restricts the grain growth after recrystallisation, which leads to a finer grain structure. Hence welding speed significantly 316 Bhatia, A. - Wattal, R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 affects the ultimate tensile strength, and the strength of the weld nugget increases with speed, but only at a maximum limit of welding speed. Hence, in order to achieve the maximum tensile strength, welding speed should be kept as high as possible. The friction temperature also increases with the increase in tool RPM. This high temperature leads to softening of the material; therefore, the value of flow stress will become lower. Higher strain rates are also achieved by increasing tool RPM. When the tool RPM is increased initially, this factor is non-dominant as the increase in temperature largely governs the size of the recrystallized grain. On further increasing the tool RPM, the increase in the strain rate takes place, and the strength starts to increase, which can be seen clearly from Fig. 4. Fig. 4. a) Surface plot showing variation of UTS with the process parameters, and b) contour plot of UTS showing variations with process parameters Increasing the tool shoulder diameter increases the amount of heat generated, but the larger tool diameter takes away the heat due to conduction. Therefore, due to the larger tool shoulder, the heat is spread over the larger area. Thus, this factor has not affected the ultimate tensile strength of the weld. 3.3.2 Effect of Process Parameters on Percentage Elongation The percentage elongation depends upon the strain hardening behaviour of the material that is, in turn, dependent on the grain size. It is evident from Fig. 5 that the percentage elongation increases with the welding speed. This is because of the reduction in grain size due to the higher strain rates. The effect of Tool RPM in increasing the strain rate and increasing the temperature during welding seem to have compensated each other within the range; therefore, it is independent of variation of tool RPM. Also due to the larger tool diameter, heat is spread all over the larger area; thus, it does not affect the percentage elongation. Fig. 5. a) Surface plot showing variation of percentage elongation with the process parameters, and b) contour plot of percentage elongation showing variations with process parameters 3.3.3 Effect of Process Parameters on Percentage Reduction in Area As seen from Fig. 6, the percentage reduction in area increases with the welding speed. The reason for this is the reduction in grain size due to the higher strain rate. The effects of tool RPM in increasing the strain rate and increasing the temperature during welding seem to have compensated for each other within the range in which process parameters were varied. A Process Parameters Optimization for Maximizing Tensile Strength in Friction Stir-Welded Carbon Steel 317 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 slight effect of tool shoulder diameter is also observed on the percentage reduction in area. 2.0 2.5 Welding speed Fig. 6. a) Surface plot showing variation of percentage reduction in area with the process parameters, and b) contour plot of percentage reduction in area showing variations with process parameters 3.4 Results of the Confirmation Experiment A predicted maximum tensile strength is 564.37 MPa for the values of welding speed 210 mm/min, Tool RPM as 750 mm/min and tool shoulder diameter as 20 mm. The value of UTS as determined experimentally under the same conditions was 553.5 MPa, which shows the consistency of the model. The results of the confirmation experiment are illustrated in Table 8. 3.5 Observations on Microstructure Fig. 7 illustrates the microstructure of the weld obtained under optimized conditions at different locations. The location from which the micrographs have been taken is shown in Fig. 7. The base metal exhibits ferrites and small grains of pearlites, as shown in Fig 7a. The optical micrograph consists of regions containing the fine-grained heat-affected zone, interface between heat affected zone and stir zone, subregion inside the weld nugget as shown by Figs. 7b to d. The reason behind this ferrite and pearlite growth is due to relatively slow cooling rates. Fig. 7d represents the formation of refined grain structure in the stir zone. The microstructural analysis reveals that heat affected zone consists of ferrite and pearlites. The ferrites are represented by the light-etched constituents, and pearlite is represented by the dark-etched constituent. The refinement of grains was also observed in the stir zone as a result of which refined grain structures are observed in this zone. The boundary between the transformational zone is successfully identified, which shows the clear difference of recrystallisation. Fig. 7. Optical micrograph of various regions of the friction stir welded steel; a) base metal (region A), b) heat affected zone (region B), c) interface between the heat affected zone and stir zone (region C), and d) stir zone (region D) Table 8. Results of confirmation experiment Process parameters Ultimate tensile strength Error [%] Welding speed (V) [mm/min] Tool RPM (N) [mm/min] Shoulder diameter (D) [mm] Predicted value Experimental value [MPa] [MPa] 210 750 20 564.37 553.5 1.86 318 Bhatia, A. - Wattal, R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 3.6 Observation on Microhardness of the Weld The plan followed for making indentations to conduct the microhardness testing is shown in Fig. 8. The specimen was polished on transverse section of the weld for testing purpose. Indentations were made 1 mm apart, covering the entire region of weldment on the transverse section. A total of 20 mm distance was covered during the indentation scheme, which was sufficient to cover the weld nugget, HAZ and a portion of the unaffected base material. The results of the microhardness are shown in Fig. 9. The hardness taken in the stir zone shows a constituently high value. The microhardness values in the weld nugget are higher as compared to the base metal. Although the hardness values vary across the weld nugget, they never fall below the hardness value of base metal. The increase in hardness is primarily due to the fine grain structure of the weld nugget. ? Fig. 8. Scheme followed for making indentations for conducting microhardness test Fig. 9. Microhardness profile (transverse section) of the specimen welded under optimized conditions 3.7 Fractured Surface Fig. 10 illustrates the microstructure of fractured tensile test specimen in which dimpled surfaces are predominantly seen over the entire area. The mode of failure appears to be ductile fracture, which is further confirmed by appearance of sheared surfaces and deformed craters. ? . . • . .v - , ■ ..v," mmmy ■ ■ . Fig. 10. SEM image of fractured joint region 4 CONCLUSIONS The major conclusions that can be drawn from the present study are: 1. Mathematical models were developed for the ultimate tensile strength, percentage elongation and percentage reduction in area, using regression analysis and ANOVA was used to predict the tensile strength of the joint at a 95 % confidence interval. 2. The weld nugget was reported to be the strongest region, and the unaffected base metal was reported to be the weakest region. 3. The confirmation experiment for ultimate tensile strength depicts that the developed model was adequate. 4. The ultimate tensile strength was significantly affected by the welding speed and tool RPM. The percentage elongation of the weld was affected by welding speed only, and the percentage reduction in the area was significantly affected by welding speed and tool shoulder diameter. 5. The microstructure of the welded joint for the optimized conditions was observed. Pearlite and ferrites were observed in the heat-affected zone, and refined pearlitic structures were observed in the stir zone. 6. The microhardness values in the weld nugget are higher as compared to the base metal. The increase in hardness is primarily due to the fine grain structure of the weld nugget. 7. The fracture morphology of the tensile test specimen was also carried out and revealed ductile fracture because the dimpled regions were observed over the entire area. Process Parameters Optimization for Maximizing Tensile Strength in Friction Stir-Welded Carbon Steel 319 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 5 ACKNOWLEDGEMENTS The authors would like to thank DTU, Delhi for providing the Experimental setup and Central Research Facility, NCU, Gurugram for providing the necessary testing facilities required for the present study. 6 REFERENCES [1] AWS. (1996). AWS Handbook, 8th ed., vol.3, American Welding Society, Miami. [2] Thomas, W.M., Nicholas, E.D., Needham, J.C., Murch, M.G., Temple-Smith, P., Dawes, C.J. (1995). Friction Welding, Patent no. 5460317. United States Patent and Trademark Office (USPTO), Washington. [3] Mishra, R.S., Ma, Z.Y. (2005). Friction stir welding and processing. Materials Science and Engineering: R: Reports, vol. 50, no. 1-2, p. 1-78, D0l:10.1016/j.mser.2005.07.001. [4] Hassan, Kh.A.A., Norman, A.F., Price, D.A., Prangnell, P.B. (2003). Stability of nugget zone grain structures in high strength Al-alloy friction stir welds during solution treatment. Acta Materialia, vol. 51, no. 7, p. 1923-1936, D0I:10.1016/ s1359-6454(02)00598-0. [5] Heinz, B., Skrotzki, B. (2002). Characterization of a friction-stir-welded aluminum alloy 6013. Metallurgical and Materials Transactions B, vol. 33, p. 489-498, D0I:10.1007/s11663-002-0059-5. [6] Liu, H.J., Fujii, H., Maeda, M., Nogi, K. (2003). Tensile properties and fracture locations of friction-stir-welded joints of 2017-T351 aluminum alloy. Journal of Materials Processing Technology, vol. 142, no. 3, p. 692-696, D0I:10.1016/s0924-0136(03)00806-9. [7] Mahoney, M.W., Rhodes, C.G., Flintoff, J.G., Bingel, W.H., Spurling, R.A. (1998). Properties of friction-stir-welded 7075 T651 aluminum. Metallurgical And Materials Transactions A, vol. 29, p. 1955-1964, D0I:10.1007/s11661-998-0021-5. [8] Sato, Y.S., Kokawa, H. (2001). Distribution of tensile property and microstructure in friction stir weld of 6063 aluminum. Metallurgical and Materials Transactions A, vol. 32, p. 30233031, D0I:10.1007/s11661-001-0177-8. [9] Sato, Y.S., Kokawa, H., Enomoto, M., Jogan, S., Hashimoto, T. (1999). Precipitation sequence in friction stir weld of 6063 aluminum during aging. Metallurgical and Materials Transactions A, vol. 30, p. 3125-3130, D0I:10.1007/s11661-999-0223-5. [10] Sato, Y.S., Urata, M., Kokawa, H., Ikeda, K. (2003). Hall-petch relationship in friction stir welds of equal channel angular-pressed aluminium alloys. Materials Science and Engineering: A, vol. 354, no. 1-2, p. 298-305, D0I:10.1016/s0921-5093(03)00008-x. [11] Su, J.-Q., Nelson, T.W., Mishra, R., Mahoney, M. (2003). Microstructural investigation of friction stir welded 7050-T651 aluminium. Acta Materialia, vol. 51, no. 3, p. 713-729, D0I:10.1016/s1359-6454(02)00449-4. [12] Lakshminarayanan, A.K., Balasubramanian, V. (2008). Process parameters optimization for friction stir welding of RDE-40 aluminium alloy using Taguchi technique. Transactions of Nonferrous Metals Society of China, vol. 18, no. 3, p. 548554, DOI:10.1016/s1003-6326(08)60096-5. [13] Vijayan, S., Raju, R., Rao, S.R.K. (2010). Multiobjective optimization of friction stir welding process parameters on aluminum alloy AA 5083 using Taguchi-based grey relation analysis. Materials and Manufacturing Processes, vol. 25, no. 11, p. 1206-1212, D0I:10.1080/10426910903536782. [14] Elatharasan, G., Senthil Kumar, V.S. (2014). Corrosion analysis of friction stir-welded AA 7075 aluminium alloy. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 1, p. 29-34, D0I:10.5545/sv-jme.2012.711. [15] Kadaganchi, R., Gankidi, M.R., Gokhale, H. (2015). Optimization of process parameters of aluminum alloy AA 2014-T6 friction stir welds by response surface methodology. Defence Technology, vol. 11, no. 3, p. 209-219, D0I:10.1016/j. dt.2015.03.003. [16] Babu, N., Karunakaran, N., Balasubramanian, V. (2017). A study to estimate the tensile strength of friction stir welded AA 5059 aluminium alloy joints. The International Journal of Advanced Manufacturing Technology, vol. 93, p. 1-9, D0I:10.1007/s00170-015-7391-9. [17] Satheesh, C., Sevvel, P., Senthil Kumar, R. (2020). Experimental identification of optimized process parameters for FSW of Az91C mg alloy using quadratic regression models. Strojniški vestnik - Journal of Mechanical Engineering, vol. 66, no. 12, p. 736-751, D0I:10.5545/sv-jme.2020.6929. [18] Lienert, T.J., Stellwag, W.L., Jr., Grimmett, B.B., Warke, R.W. (2003). Friction stir welding studies on mild steel. Supplement To The Welding Journal, vol. 82, p. 1s-9s. [19] Ueji, R., Fujii, H., Cui, L., Nishioka, A., Kunishige, K., Nogi, K. (2006). Friction stir welding of ultrafine grained plain low-carbon steel formed by the martensite process. Materials Science and Engineering: A, vol. 423, no. 1-2, p. 324-330, D0I:10.1016/j.msea.2006.02.038. [20] Fujii, H., Cui, L., Tsuji, N., Maeda, M., Nakata, K., Nogi, K. (2006). Friction stir welding of carbon steels. Materials Science and Engineering: A, vol. 429, no. 1-2, p. 50-57, D0I:10.1016/j.msea.2006.04.118. [21] Cui, L., Fujii, H., Tsuji, N., Nogi, K. (2007). Friction stir welding of a high carbon steel. Scripta Materialia, vol. 56, no. 7, p. 637-640, D0I:10.1016/j.scriptamat.2006.12.004. [22] Miles, M.P., Nelson, T.W., Steel, R., Olsen, E., Gallagher, M. (2009). Effect of friction stir welding conditions on properties and microstructures of high strength automotive steel. Science and Technology of Welding and Joining, vol. 14, no. 3, p. 228-232, D0I:10.1179/136217108x388633. [23] Ghosh, M., Hussain, M., Kumar Gupta, R. (2012). Effect of welding parameters on microstructure and mechanical properties of friction stir welded plain carbon steel. ISIJ International, vol. 52, no. 3, p. 477-482, D0I:10.2355/ isijinternational.52.477. [24] Lakshminarayanan, A.K., Balasubramanian, V., Salahuddin, M. (2010). Microstructure, tensile and impact toughness properties of friction stir welded mild steel. Journal of Iron and Steel Research International, vol. 17, p. 68-74, D0I:10.1016/ s1006-706x(10)60186-0. 320 Bhatia, A. - Wattal, R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 311-321 [25] Ghosh, M., Kumar, K., Mishra, R.S. (2012). Process optimization for friction-stir-welded martensitic steel. Metallurgical and Materials Transactions A, vol. 43, p. 19661975, DOI:10.1007/s11661-012-1084-x. [26] Yabuuchi, K., Tsuda, N., Kimura, A., Morisada, Y., Fujii, H., Serizawa, H., Nogami, S., Hasegawa, A., Nagasaka, T. (2014). Effects of tool rotation speed on the mechanical properties and microstructure of friction stir welded ODS steel. Materials Science and Engineering: A, vol. 595, p. 291-296, DOI:10.1016/j.msea.2013.12.022. [27] Karami, S., Jafarian, H., Eivani, A.R., Kheirandish, S. (2016). Engineering tensile properties by controlling welding parameters and microstructure in a mild steel processed by friction stir welding. Materials Science and Engineering: A, vol. 670, p. 68-74, DOI:10.1016/j.msea.2016.06.008. [28] Mironov, S., Sato, Y.S., Yoneyama, S., Kokawa, H., Fujii, H.T., Hirano, S. (2018). Microstructure and tensile behavior of friction-stir welded trip steel. Materials Science and Engineering: A, vol. 717, p. 26-33, DOI:10.1016/j. msea.2018.01.053. [29] Mahmoudiniya, M., Kokabi, A.H., Kheirandish, S., Kestens, L.A.I. (2018). Microstructure and mechanical properties of friction stir welded ferrite-martensite DP700 steel. Materials Science and Engineering: A, vol. 737, p. 213-222, DOI:10.1016/j.msea.2018.09.013. [30] Lee, S.-J., Park, T.M., Nam, J.-H., Choi, W.S., Sun, Y., Fujii, H., Han, J. (2019). The unexpected stress-strain response of medium Mn steel after friction stir welding. Materials Science and Engineering: A, vol. 744, p. 340-348, DOI:10.1016/j. msea.2018.12.041. [31] Standard test methods for tension testing of metallic materials, in E 8M-04,2004, AST M International: West Conshohocken, PA 19428-2959, p. 1-6. Process Parameters Optimization for Maximizing Tensile Strength in Friction Stir-Welded Carbon Steel 321 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 © 2021 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2021.7161 Original Scientific Paper Received for review: 2021-03-13 Received revised form: 2021-05-28 Accepted for publication: 2021-06-18 Experimental Investigation of a CryogenicaUy Cooled Oxygen-mist Near-dry Wire-cut Electrical Discharge Machining Process Boopathi Sampath1* - Sureshkumar Myilsamy2 1Muthayammal Engineering College, Department of Mechanical Engineering, India 2Bannari Amman Institute of Technology, Department of Mechanical Engineering, India In this paper, a novel method of cryogenically cooled (low-temperature nitrogen gas) wire tool is used during the oxygen-mist near-dry wire-cut electrical discharge machining (NDWEDM) process to cut Inconel 718 alloy material. The current, pulse-width, pulse-interval, and flow rate are the controllable variables for response characteristics, such as the material removal rate (MRR) and wire wear ratio (WWR). The Box-Behnken method is applied to design the experiments to collect the observations from experiments. The mathematical models for each response were developed using significant individual, interaction, and quadratic terms by the sequential sum of the square test. The response surfaces were developed. It was revealed from the analysis that 52.92 % of current, 24.63 % of Pulse-width, 12.81 % of pulse- interval and 5.75 % of flow rate contributed to MRR, while 14.89 % of current, 9.75 % of pulse-width, 62.20 % of pulse-interval, and 5.44 % of flow rate contributed to WWR. The pulse-width has more contribution on MRR due to the long period of spark between the wire and work materials. It was also observed that the pulse-interval has more effect on WWR due to the more ideal period (high spark-pause-time) between two consecutive high-temperature sparks in the wire tool. The wear of the wire tool has been analysed using scanning electron microscopy (SEM) photographs. The desirability principles were first applied to obtain multi-objective solutions with a combination of process parameters to achieve the optimal values of both responses. The predicted combination of results has been validated by data that were collected from confirmation experiments. Keywords: cryogenically cooled, oxygen-mist, near-dry, wire-cut EDM, MRR, WWR, Box-Behnken method Highlights • A novel method of cryogenically cooled (low-temperature nitrogen gas) wire electrode tool in the oxygen-mist near-dry wire-cut electrical discharge machining (NDWEDM) process was experimented with to cut the Inconel 718 alloy material. • Wire wear ratio and material removal rate of cryogenically cooled near-dry WEDM process was first investigated in this research. • It was revealed that the pulse-interval has more effect on wire wear ratio due to a more ideal period (high spark-pause-time) between two consecutive high-temperature spark in the wire tool. • The wear of the wire tool has been analysed using Scanning Electron Microscopy test photographs. • The desirability principles were first applied to obtain the multi-objective solutions to optimize both responses. 0 INTRODUCTION In an unconventional machining process, the relationships between manufacturing parameters and environmental impact are developed to analyse the material removal mechanics, tool change, minimum rejections in production, and the effects of cutting-fluid flow [1]. The environmental impact of machining processes should be analysed for minimizing environmental impacts by the modification of existing technology and the development of new manufacturing methods [2]. In these aspects, research into the modification of EDM and WEDM processes was developed to make a trade-off between machining performance and machining pollutions [3]. The analytical relationship of EDM processes was developed to reveal the wear of the tool and workpiece, the dielectric fluid flows, and toxicity and flammability [3] and [4]. The inferences of cooling electrode wear and surface roughness of the workpiece have been investigated by changing process parameters, such as voltage, pulse-width, current, and pulse-interval [5] and [6]. The parametric analysis of dry electric discharge machining of mild steel was investigated, and response models were developed using response surface methodology [7]. Generally, the machining performance of the dry EDM process is very low compared to the conventional process. It was revealed that the tool wear of the dry EDM process is significantly reduced by cryogenic cooling of the electrode and workpiece [8]. The mechanism of the gas-liquid-powder mixture in EDM was investigated in both dry and near-dry processes to improve material erosion. The cryogenically tested brass wire produces a 22.55 % greater material removal rate (MRR) compared to an untreated brass wire [9]. A parametric study was performed using a molybdenum wire tool and tool steel workpiece with an air dielectric to investigate the influence of air-mist pressure, voltage, pulse-duration, pulse-width, and current on the MRR and Ra using the Taguchi technique [5] and [10]. Later, the oxygen gas near-dry WEDM experiments were conducted using Taguchi's L27 orthogonal array and multi-objective artificial bee colony (MOABC) 322 *Corr. Author's Address: Department of Mechanical Engineering, Muthayammal Engineering College, Rasipuram-637408, Namakkal (Dt), India, boopasangee@gmail.com. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 algorithm [11] and [12]. Cryo-treated wire electrodes, liquid nitrogen, and zinc-coated brass wires were investigated to reduce the tool wear rate for the green environment [4], [13] to [15]. Recently, the electrical conductivity of cryo-cooled molybdenum wire has been increased by cryogenic treatments [16]. However, no synchronized cryo-treated gas-liquid mixer was found in the near-dry EDM and WEDM process. In this research, the data from 29 observations from cryogenically cooled oxygen-mist near-dry WEDM experiments were collected. The collected data are used to predict the data for better results of the near-dry WEDM process. The significant process parameters were identified to improve the quality of cutting processes. The desirability approach is used to convert the single-objective problem into a multi-objective optimization problem. The WEDM machine manufacturer and operators utilize the optimum results to set the optimal process parameters for the best machining performances. 1 EXPERIMENTATION 1.1 Experimental Setup The cryo-cooled wire electrode setup was developed in numerically controlled (NC) wire-cut electrical discharge machining. Liquid nitrogen was stored in a dewar flask to maintain the cryogenic temperature. The molybdenum wire was cooled on both sides of the electrode movements. The oxygen and dielectric fluid mixture are used as a working medium in the reciprocating WEDM machine. Based on trial experiments, the input parameters and their significance levels are identified. The billet size of 718 is 50 mm x 50 mm x 5 mm Inconel 718 is used as work material for near-dry WEDM machining processes. The experimental setup of cryo-cooled near-dry WEDM is shown in Fig. 1. Sub-123 K of liquid nitrogen was stored. The 253 K temperature and 10 g/s mass flow rate of N2 were used to cool the molybdenum wire during the cutting process. The Table 1. Parameter and Machine Setting level oxygen and the dielectric fluid mixture were used as a dielectric medium in reciprocating the WEDM machine. The MRR in [mm3/min] can be calculated by the volume of materials removed concerning time using Eqs. (1) and (2). Kerf = wire diameter x (2 x sparking gap), thickness x Kerf x length of cut MRR = - time (1) (2) The wire wear ratio (WWR) has been measured from the loss of wire materials during the cutting process concerning the time, and the initial weight of wire is to be taken before machining process [17] (Eq. (3)). WWR = weight loss of wire (3) initial weight of wire Based on exploratory experiments, the input parameters and their significance levels are identified. The levels of each process variable are tabulated in Table 1. The near-dry WEDM experiments are conducted [18], and the MRR and WWR values observed through the experiments are shown in Table 2. Fig. 1. Cryogenically cooled oxygen-mist near-dry WEDM experimental set Description Symbol Parameter Units Low High Mean C Current A 3 5 4 PW Pulse-width ^s 15 25 20 PI Pulse-interval ^s 45 75 60 F Flow rate ml/min 10 20 15 Dielectric medium Oxygen gas mixed with water Wire treatment Cryogenic Nitrogen gas during the machining process Output parameters MRR in [mm3/min], and WWR Experimental Investigation of a Cryogenically Cooled Oxygen-mist Near-dry Wire-cut Electrical Discharge Machining Process 323 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 Table 2. Design of Experiments and observations using Box-Behnken method Exp. No. C PW PI F MRR WWR 1 3 20 75 15 6.13 0.464 2 5 20 75 15 9.45 0.626 3 4 15 60 20 8.38 0.643 4 4 20 45 10 9.28 0.994 5 4 25 60 10 9.57 0.921 6 4 20 60 15 9.12 0.783 7 4 20 60 15 9.09 0.724 8 4 15 75 15 7.12 0.424 9 5 20 60 20 10.9 0.829 10 3 25 60 15 8.51 0.864 11 3 20 60 10 6.58 0.745 12 4 25 45 15 10.91 0.991 13 4 20 75 10 7.78 0.635 14 3 20 45 15 7.76 0.848 15 4 15 60 10 7.38 0.767 16 4 15 45 15 8.67 0.921 17 4 20 75 20 8.74 0.478 18 5 20 45 15 10.87 1.050 19 5 20 60 10 9.44 0.934 20 4 25 75 15 9.38 0.670 21 5 25 60 15 10.77 0.872 22 4 20 45 20 10.29 0.905 23 5 15 60 15 9.47 0.913 24 4 20 60 15 9.12 0.778 25 4 20 60 15 9.11 0.78 26 4 25 60 20 10.63 0.807 27 4 20 60 15 8.77 0.772 28 3 15 60 15 6.02 0.502 29 3 20 60 20 7.24 0.621 1.2 Design of Experiments The Box-Behnken method is used to conduct the experiments by design expert software. The Box-Behnken design is a self-determining quadratic design, which does not contain the partial factorial design. The designs have restricted capability related to the central composite designs. Five central points are repeated to avoid bios errors. The design uses 8 trails from (2 x 4 = 8) two levels of k parameters, 16 trails of two factorial design (24 = 8), and five repeated central points to calculate lack-of-fit. Twenty-nine sets of experiments were conducted, and the observed responses are tabulated in Table 2. Based on the analysis of the variance test, the significant individual and interaction and quadratic terms were identified [5]. Insignificant terms are eliminated from the model. If the response is 'f(x)', the independent variables are x1, x2.....xn, and the response model is developed by following general Eq. (4). f ( x) = P0 -YLjXj+. h, ¿=i j and i < j, (4) where k is the number of process variables; P0, P, and Pj are the model coefficients; $ is the statistical error, which represents variability by other noises. 2 RESULT ANALYSIS AND DISCUSSIONS The sequential sum of the square test was used to select the optimum model for the analysis. Initially, the linear model is selected [12]. However, the model is not significant due to the coefficient of determination (R2) 324 Sampath, B. - Myilsamy, S. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 value of MRR and WWR is very low (MRR R2 =0.072 and WWR R2 = 0.81). The two factors interaction model (2FI) also did not fit with the solution due to the R2 value of both models is minimum (MRR R2 = 0.608 and WWR R2 = 0.075). Then, the quadratic model of MRR was selected due to the R2 value of 0.998 of both responses. The insignificant terms of the models were eliminated from quadratic models. The lack of fit of the model is insignificant and thus acceptable. The cubic models of both responses were not selected due to more allies' terms in the models. The lack of fit tests of MRR and WWR are shown in Tables 3 and 5, respectively. The analysis of the variations of MRR and WWR concerning process parameters are shown in Tables 4 and 6, respectively. The regression models of MRR and WWR were developed to identify the inference of process variables as shown in Eqs. (5) and (6) respectively. MRR = -8.241 + 5.335 x C + 0.45 x PW - 0.065 x PI -0.0575 x F - 0.0595 x Cx PW +3.5x10-3 x Cx PI + 0.04x Cx F -0.425 x C2, (5) R2 = 99.40 %, Adjusted R2 = 99.14 %, and Predicted R2 = 98.18 %. WWR = -0.529 + 0.501x C + 0.061x PW -6.161x10-3 x PI +1.727 x 10-3 x F -0.0201x Cx PW + 5.867 x10-4 x PWx PI -2.267x10-4 xPIx F-1.298x 10-4 xPI2, (6) R2 = 99.40 % , Adjusted R2 = 99.16 %, and Predicted R2 = 99.04 %. The insignificant terms are eliminated from MRR and WWR regression models to improve the predicted R2 and adjusted R2. The regression models are used to plot the response surface between response and process variables. The influences of the interaction effects of process parameters have been studied using the response surfaces. The response surface of MRR for flow rate and current is shown in Fig. 2. It is also significantly enhanced by increasing the flushing flow rate due to the quick disposal of debris from the cutting zone. Fig. 3 shows that the response surface of MRR by the pulse-width vs current. The MRR is improved by the increase in spark current between work materials and wire due to high spark strength [19]. It was observed that the large increase in MRR is seen from low pulse-width to high value by enhancing spark strength, as shown in Fig. 4 [12]. However, the MRR is increased by reducing pulse-interval due to an increase in spark ideal time. 1000 3 Fig. 2. Response surface for MRR concerning flow rate and current 15 3 Fig. 3. Response surface for MRR concerning pulse-width and current Fig. 4. Response surface for MRR concerning current and pulse-interval Experimental Investigation of a Cryogenically Cooled Oxygen-mist Near-dry Wire-cut Electrical Discharge Machining Process 325 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 Table 3. Lack of fit test for material removal rate model Source Sum of squares Degree of freedom Mean sum square F-value /»-value > F Comments Linear 2.040 20 0.102 4.383 0.081 Insignificant 2FI 1.513 14 0.108 4.644 0.075 Insignificant Quadratic 0.215 10 0.022 0.924 0.584 Suggested Cubic 0.036 2 0.018 0.775 0.520 Aliased Pure Error 0.093 4 0.023 - - - Table 4. Data analysis of material removal rate Source Sum of squares Degree of freedom Mean square F-value /-value > F Remarks Contribution [%] Model 54.491 8 6.811 404.036 < 0.0001 Significant 99.39 Current (C) 29.016 1 29.016 1721.184 < 0.0001 Significant 52.92 Pulse-width (PW 13.504 1 13.504 801.052 < 0.0001 Significant 24.63 Pulse-interval (PI) 7.023 1 7.023 416.571 < 0.0001 Significant 12.81 Flow rate(F 3.152 1 3.152 186.962 < 0.0001 Significant 5.75 C x PW 0.354 1 0.354 21.000 0.0002 Significant 0.65 C x PI 0.011 1 0.011 0.654 0.4282* Not significant 0.02 C x F 0.160 1 0.160 9.491 0.0059 Significant 0.29 C 2 1.271 1 1.271 75.369 < 0.0001 Significant 2.32 Residual 0.337 20 0.017 - - - - Lack of fit 0.244 16 0.015 0.656 0.7580 Not significant - Pure error 0.093 4 0.023 - - - - Cor total 54.8281 28 - - - - - Table 5. Lack of fit test for WWR model Source Sum of squares Degree of freedom Mean square F-value /-value > F Comments Linear 0.06 20.00 0.00 4.78 0.070 Insignificant 2FI 0.01 14.00 0.00 0.92 0.601 Insignificant Quadratic 0.00 10.00 0.00 0.11 0.998 Suggested Cubic 0.00 2.00 0.00 0.08 0.927 Aliased Pure Error 0.00 4.00 0.00 - - - Table 6. Data analysis of wire wear ratio Source Sum of squares Degree of freedom Mean square F-value /-value > F Remarks Contribution [%] Model 0.7747 8.0000 0.0968 411.8384 < 0.0001 Significant 99.40 Current (C) 0.1160 1.0000 0.1160 493.4670 < 0.0001 Significant 14.89 Pulse-width (PW) 0.0760 1.0000 0.0760 323.2220 < 0.0001 Significant 9.75 Pulse-interval (PI) 0.4848 1.0000 0.4848 2061.8105 < 0.0001 Significant 62.20 Flow rate(F) 0.0424 1.0000 0.0424 180.1662 < 0.0001 Significant 5.44 C x PW 0.0406 1.0000 0.0406 172.6734 < 0.0001 Significant 5.21 PW x PI 0.0077 1.0000 0.0077 32.9337 < 0.0001 Significant 0.99 PI x F 0.0012 1.0000 0.0012 4.9162 0.0384 Significant 0.15 PI2 0.0060 1.0000 0.0060 25.5181 < 0.0001 Significant 0.77 Residual 0.0047 20.0000 0.0002 - 0.60 Lack of fit 0.0023 16.0000 0.0001 0.2360 0.9842 Not significant - Pure error 0.0024 4.0000 0.0006 - - - - Cor total 0.7794 28.0000 - - - - - The minimization of WWR is one of the goals of pulse-interval and pulse-width due to fine and soft this study. The WWR is minimum at the low value of spark in the cutting zone, as shown in Fig. 5. While 326 Sampath, B. - Myilsamy, S. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 the increase in pulse-width, the spark strength is significantly increased with material removal rate and reduced the WWR. The increase in pulse-interval is increasing the WWR due to the high spark pause time [5]. Fig. 6 shows the interaction effects of pulse-width and current on WWR. While increasing the spark current, the WWR value is exploiting due to heave spark intensity [20] and [21]. Similarly, the WWR value is also slowly increased due to the growing MRR by fast flushing debris, as shown in Fig. 7. Fig. 5. Response surface for WWR concerning pulse-width and pulse-interval -I^OO 3 00 Fig. 6. Response surface for WWR concerning pulse-width and current The wire wear ratio has been analysed with scanning electron microscopy (SEM) photograph, as shown in Fig. 8. It was observed that the wear on the wire is linearly along the axis due to the reciprocating wire longitudinally. The crater of materials in the wire is high at the point 'P' due to the sudden supply of current to the wire. The path 'QR' is sparking a portion of the wire, which is in direct with the workpiece. The high temperature along the 'QR' path due to frequent changes of power supply between the wire and workpiece causes to increase the wear ratio. The path 'ST' is a non-spark portion of the wire tool, which has the minimum crater of wear. 70« TS 00 Fig. 7. Response surface for WWR concerning flow rate and pulse-interval Fig. 8. SEM photograph of wire electrode wear 3 MULTI-OBJECTIVE OPTIMIZATION A ND VALIDATION OF PREDICTION In this stage, the combination of experimental parameters and their response based on the standard ranges defined for the responses are predicted [22]. Based on the desirability function, techniques were used to predict the best results of both responses. This process detects a point that improves the desirability function [23]. For validating the developed models, some solutions were selected randomly. The optimized predicted values of the output responses were compared to experimentally obtained values [24]. The combination of variables that presents the overall optimum desirability (99 %) of response and contour plots is displayed in Fig. 8. Considering all the quality attributes and using the optimization method with Experimental Investigation of a Cryogenically Cooled Oxygen-mist Near-dry Wire-cut Electrical Discharge Machining Process 327 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 Table 7. Predicted Results from desirability and validation by confirmation experiments No__PW[is] P/fcs] F [ml/min] MRR [mm3/min] WWR_Description 1_5_25_75_20_10.91_059_Desirability Approach 2_5_25_75_20_H00_0.62_Near-dry WEDM 3 5 25 75 2 litres/min 17.42 0.83 Conventional WEDM which the quality parameters were put into standard ranges (MRR and WWR), formulation one consisting of 25 ps pulse-width, 75 ps pulse-interval, 5 A current, and 20 ml/min flow rate was selected as having the maximum desirability. The combined optimizing responses are predicted as 10.91 mm3/min of material removal rate, and 0.59 of WWR. Confirmation tests were conducted to validate the predicted optimum process parameters for best MRR and WWR, as shown in Table 7. The multi-optimization results were confirmed by the mean observed values of the test. The 5 A current, pulse-width 25 ps, pulse-interval 75 ^s, and 20 ml/min flow rate gives 11 mm3/min of MRR and 0.62 % of WWR. The desirability principles were applied to obtain the multi-objective solutions to optimize both responses. It is very useful for machine operators to select the best process parameters for minimum WWR and maximum MRR. These results are used by WEDM machine manufacturers to fix the best (default) setting for cryogenically cooled oxygen-mist near-dry wire-cut electrical discharge machining (NDWEDM) process. 4 COMPARISON OF CRYOGENICALLY COOLED NEAR-DRY AND CONVENTIONAL WEDM The predicted best combinations of input parameters were considered for the comparative analysis. In Fig. 9, the MRR and WWR of cryogenically cooled near-dry WEDM are compared with the conventional process. The range of input parameters of both processes are current 5 A, pulse-width 25 ^s, and pulse-interval 75 ^s. The flow rate of oxygen-mist of near-dry WEDM is 20 ml/min, and the water flow rate of conventional WEDM is 2 l/min. As per the literature, the near-dry WEDM is an effective eco-friendly process while comparing with conventional WEDM [4], [13] to [15]. The MRR of the near-dry WEDM process was comparatively lower than the conventional process due to the significant flush out of removed material from the workpiece. The WWR of the near-dry process is lower than conventional WEDM because the heat dissipation from the wire is improved by increasing the thermal conductivity of the cryogenic cooling wire [16]. The minimum WWR promotes the life of reusable wire electrode. 5 CONCLUSIONS In this research, data were collected from the experiments of oxygen-mist cryo-cooled wire near-dry WEDM were carried out to maintain enough temperature in the cutting zone to cut the Inconel 718 alloy material. The electrical conductivity of cryogenically cooled molybdenum wire was significantly increased in the near-dry WEDM process. It was observed that pulse-width, pulse-interval, current, and flow rate are significant parameters on the material removal rate and wire wear ratio. The current and pulse-width are the most important factors for material removal rate. Increasing the current and pulse-width, the material removal rate Fig. 9. Comparisons of machining performances between cryogenically cooled near-dry and conventional WEDM 328 Sampath, B. - Myilsamy, S. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 was increasing and same time decreasing the WWR. While increasing pulse-interval from 45 ps to 75 ^s, both the material removal rate and WWR were decreased. The WWR has been significantly reduced at the dielectric fluid flow rate of 20 ml/min. It was observed from the SEM photograph that the wire wear rate is linear to the axis, and the wear rate is high along the spark path of the wire tool. Thus, the optimum process parameters were obtained for the cryo-cooled molybdenum wire electrode used in the oxygen-mist near-dry wire-cut electrical discharge machining process. The desirability technique was used to find the best solution among the conflict behaviour of MRR and WWR. It was observed that the behaviour of the parameters against response variables is the same as with the response surface method. The current 5 A, pulse-width 25 ps, pulse-interval 75 ps, and 1 ml/min flow rate give 10.97 mm3/min of MRR and 0.59 of WWR. The multi-optimization results were validated by the mean observation value of confirmation experiments. These results are used by manufacturers and operators to fix the best (default) setting for the cryogenically cooled oxygen-mist NDWEDM process. In near-dry WEDM, the cryogenically cooling of wire significantly contributed to reducing WWR than to MRR. While comparing conventional WEDM, the WWR of near-dry WEDM was very lower. The lowest WWR increases the life of reusable molybdenum wire tools. As per the literature, the near-dry WEDM is an effective eco-friendly process while comparing with conventional WEDM. However, the MRR of near-dry WEDM is lower than that of the conventional WEDM process. 6 ACKNOWLEDGEMENTS This research was conducted in the Bannariamman Institute of Technology, Erode, Tamilnadu, India. for the data collection, and interpretation of results to submit the article for publication. 7 NOMENCLATURES PW pulse-width, [ps] PI pulse-interval, [^s] C spark current, [A] F flowrate, [ml/s] MRR material removal rate [mm3/min] WWR wire wear ratio [-] fX) Function of 'X x1. x2.....xn independent variables fio. fii fij model coefficients $ statistical error R2 coefficient of determination FI factor interaction F-value Values from statistical "F" table 8 REFERENCES [1] Shabgard, M.R., Seyedzavvar, M., Oliaei, S.N.B. (2011). Influence of input parameters on the characteristics of the EDM process. Strojniski vestnik - Journal of Mechanical Engineering, vol. 57, no. 9, p. 689-696, D0I:10.5545/sv-jme.2011.035. [2] Munoz, A.A., Sheng, P. (1995). An analytical approach for determining the environmental impact of machining processes. Journal of Materials Processing Technology, vol. 53, no. 3-4, p. 736-758, D0I:10.1016/0924-0136(94)01764-R. [3] Yeo, S.H., Tan, H.C., New, A.K. (1998). Assessment of waste streams in electric-discharge machining for environmental impact analysis. Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, vol. 212, no. 5, p. 393-400, D0I:10.1243/0954405981515996. [4] Hewidy, M.S., El-Taweel, T.A., El-Safty, M.F. (2005). Modelling the machining parameters of wire electrical discharge machining of Inconel 601 using RSM. Journal of Materials Processing Technology, vol. 169, no. 2, p. 328-336, D0I:10.1016/j.jmatprotec.2005.04.078. [5] Boopathi, S., Sivakumar, K. (2016). Optimal parameter prediction of oxygen-mist near-dry Wire-cut EDM. International Journal of Manufacturing Technology and Management, vol. 30, no. 3-4, p. 164-178, D0I:10.1504/IJMTM.2016.077812. [6] Abdulkareem, S., Khan, A.A., Konneh, M. (2009). Reducing electrode wear ratio using cryogenic cooling during electrical discharge machining. International Journal of Advanced Manufacturing Technology, vol. 45, no. 11-12, p. 1146-1151, D0I:10.1007/s00170-009-2060-5. [7] Saha, S.K., Choudhury, S.K. (2009). Experimental investigation and empirical modeling of the dry electric discharge machining process. International Journal of Machine Tools and Manufacture, vol. 49, no. 3-4, p. 297-308, D0I:10.1016/j. ijmachtools.2008.10.012. [8] Jiang, Z., Zhang, H., Sutherland, J.W. (2012). Development of an environmental performance assessment method for manufacturing process plans. International Journal of Advanced Manufacturing Technology, vol. 58, no. 5-8, p. 783790, D0I:10.1007/s00170-011-3410-7. [9] Bai, X., Zhang, Q.-H., Yang, T.-Y., Zhang, J.-H. (2013). Research on material removal rate of powder mixed near dry electrical discharge machining. International Journal of Advanced Manufacturing Technology, vol. 68, no. 5-8, p. 1757-1766, D0I:10.1007/s00170-013-4973-2. [10] Boopathi, S., Sivakumar, K. (2013). Experimental investigation and parameter optimization of near-dry wire-cut electrical discharge machining using multi-objective evolutionary algorithm. International Journal of Advanced Manufacturing Technology, vol. 67, p. 2639-2655, D0I:10.1007/s00170-012-4680-4. 329 Experimental Investigation of a Cryogenically Cooled Oxygen-mist Near-dry Wire-cut Electrical Discharge Machining Process 460 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 322-330 [11] Boopathi, S. (2012). Experimental comparative study of near-dry wire-cut electrical discharge machining (WEDM). European Journal of Scientific Research, vol. 75, no. 4, p. 472-481. [12] Boopathi, S., Sivakumar, K. (2014). Study of water assisted dry wire-cut electrical discharge machining. Indian Journal of Engineering and Materials Sciences, vol. 21, no. 1, p. 75-82. [13] Rajyalakshmi, G., Venkata Ramaiah, P. (2013). Multiple process parameter optimization of wire electrical discharge machining on Inconel 825 using Taguchi grey relational analysis. International Journal of Advanced Manufacturing Technology, vol. 69, p. 1249-1262, D0I:10.1007/s00170-013-5081-z. [14] Garg, M.P., Kumar, A., Sahu, C.K. (2017). Mathematical modeling and analysis of WEDM machining parameters of nickel-based super alloy using response surface methodology. Sadhana - Academy Proceedings in Engineering Sciences, vol. 42, no. 6, p. 981-1005, D0I:10.1007/s12046-017-0647-3. [15] Nayak, B.B., Mahapatra, S.S. (2016). Optimization of WEDM process parameters using deep cryo-treated Inconel 718 as work material. Engineering Science and Technology, an International Journal, vol. 19, no. 1, p. 161-170, D0I:10.1016/j.jestch.2015.06.009. [16] Myilsamy, S., Boopathi, S., Yuvaraj, D. (2021). A study on cryogenically treated molybdenum wire electrode. Materials Today: Proceedings, vol. 45, p. 8130-8135, D0I:10.1016/j. matpr.2021.02.049. [17] Shather, S.K., Mohammed, M.T. (2018). Investigate WEDM process parameters on wire wear ratio, material removal rate and surface roughness of steel 1012 AISI. Engineering and Technology Journal, vol. 36, no. 3, p. 256-261, D0I:10.30684/ etj.36.3A.3. [18] Boopathi, S. (2019). Experimental investigation and parameter analysis of LPG refrigeration system using Taguchi method. SN Applied Sciences, vol. 1, no. 8, p. 892, D0l:10.1007/s42452-019-0925-2. [19] Boopathi, S., Sivakumar, K., Kalidas, R. (2012). Parametric study of dry WEDM using Taguchi method. International Journal of Engineering Research and Development, vol. 2, no. 4, p. 63-68. [20] Pan dey, H., Dhakar, K., Dvivedi, A., Kumar, P. (2015). Parametric investigation and optimization of near-dry electrical discharge machining. Journal of Scientific and Industrial Research, vol. 74, no. 9, p. 508-511. [21] Singh, J., Sharma, R.K. (2017). Green EDM strategies to minimize environmental impact and improve process efficiency. Journal for Manufacturing Science and Production, vol. 16, p 1-18, D0I:10.1515/jmsp-2016-0034. [22] Srivastava, M., Maheshwari, S., Kundra, T., Rathee, S. (2017). Multi-response optimization of fused deposition modelling process parameters of ABS using response surface methodology (RSM)-based desirability analysis. Materials Today: Proceedings, vol. 4, no. 2, p. 1972-1977, D0I:10.1016/j.matpr.2017.02.043. [23] Senthil, S.M., Parameshwaran, R., Ragu Nathan, S., Bhuvanesh Kumar, M., Deepandurai, K. (2020). A multi-objective optimization of the friction stir welding process using RSM-based-desirability function approach for joining aluminum alloy 6063-T6 pipes. Structural and Multidisciplinary Optimization, vol. 62, p. 1117-1133, D0I:10.1007/s00158-020-02542-2. [24] Vahdati, M., Moradi, M., Shamsborhan, M. (2020). Modeling and optimization of the yield strength and tensile strength of Al7075 butt joint produced by FSW and SFSW using RSM and desirability function method. Transactions of the Indian Institute of Metals, vol. 73, p. 2587-2600, D0I:10.1007/ s12666-020-02075-8. 330 Sampath, B. - Myilsamy, S. Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 © 2021 Journal of Mechanical Engineering. All rights reserved. DOI:10.5545/sv-jme.2021.7104A Short Scientific Paper Received for review: 2021-01-22 Received revised form: 2021-04-04 Accepted for publication: 2021-05-14 Reliability-Based Design Optimization of Pump Penetration Shell Accounting for Material and Geometric Non-Linearity Gautham Velayudhan* - Prabhu Raja Venugopal - Ebron Shaji Gnanasigamony Thankareathenam - Mohanraj Selvakumar - Thyla Pudukarai Ramaswamy PSG College of Technology, Department of Mechanical Engineering, India The roof slab of a nuclear reactor supports all the components and sub-systems. It needs to resist the seismic loads in accordance with load-carrying criteria. The static stress analysis of the reactor roof slab reveals that high-stress concentration was present in the pump penetration shell (PPS) that supports the primary sodium pump. This paper presents an assessment of collapse load and the optimization of the pump penetration shell through the reliability approach, accounting for material non-linearity, geometrical non-linearity and randomness in loading. In addition, the load-carrying capacity of PPS was determined considering two different materials: IS2062 and A48P2. The design of experiments (DoE) was formulated considering the flange angle and flange thickness as parameters. An empirical model for load function was formulated from the results of the collapse load obtained for various combinations of design parameters. The above function was used to perform the reliability-based geometry optimization of the PPS of the roof slab. Keywords: nuclear reactor; buckling; optimization; reliability; limit load; genetic algorithm Highlights • The reliability-based design of the pump penetration shell is adapted to account for material properties that are uncertain. • The methodology of incorporating the reliability index for reliability-based design has been studied, and a genetic algorithm optimization technique is adopted for performing design optimization of the pump penetration shell. • The prediction of collapse load of a single-layer PPS reveals that the PPS made of IS2062 has a higher load-carrying capacity than the PPS made of A48P2. • The PPS made of IS2062 has less influence of effect of mode shapes and imperfection on the collapse load compared with the PPS made of material A48P2. • The load-carrying capacity of the PPS increasing with the increase of angle and thickness of the flange. 0 INTRODUCTION Considering the geometrical configuration of the roof slab of a nuclear reactor, a large box-type structure with many penetrations made of carbon steel, posed many difficulties during manufacturing, particularly due to lamellar tearing. Alternatively, taking advantage of high-load carrying capacity with possibly minimum thickness, a dome-shaped roof slab was conceived comprised of a conical shell connected to the vertical shell through a short torus portion [1]. The structural integrity of the roof slab is significant for the safe and reliable functioning of the reactor [2]. Structural integrity refers to the condition of a structural system and implies that the structure and its components remain intact over the intended lifetime of the structure. The roof slab is subjected to fatigue loading due to the reaction taking place inside the main vessel and the intermittent flow of sodium [1]. The fatigue phenomenon is a complex and progressive failure; one possibility of occurrence could be from a small weld defect. The failure begins with the formation of micro-cracks and crack propagation, leading to the failure of the structure. Reliability is an aspect of engineering uncertainty. It is the probability that an item will perform a required function without failure under stated conditions for the stated period [3] to [7]. Durability is a particular aspect of reliability related to the ability of an item to withstand the effects of time or operating cycles, etc. The main objective of any structural design is to ensure the safety and economy of the structure operating under a given environment. Hence, the demand should not exceed the capacity of the structure. Capacity (C > Demand (D) The above-cited condition should be satisfied so that the structure's safety is ensured for the intended purpose for which the structure is built. The deterministic approach is based on the premise that a given problem can be stated in the form of a question or a set of questions to which there is an explicit and distinct answer. A probabilistic approach is based on the concept that several or varied outcomes of a situation are possible to this approach in which uncertainty is recognized. Probabilistic modelling aims at the study of a range of outcomes for the given set of input data. Accordingly, the description of a PSG College of Technology, Peelamedu, Coimbatore, India, 1701rm01@psgtech.ac.in 331 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 physical situation or system includes the randomness of data and other kinds of uncertainties. Fig. 1. Probabilistic approach The probability of failure is the common area shared by both the demand and capacity curves as shown in Fig. 1. The reason for this is that C is always less than D in the shaded area, hence ensuring failure. Considering the probability density function for the normal distribution of capacity and demand curve, the measure of reliability is defined with a term known as the reliability index (#), as shown in Fig. 1, as proposed by Cornell et al. [8]. It is the shortest distance measured from the origin when the capacity and demand are expressed. The random variables are characterized by their first two moments only, i.e., mean and variance. This is achieved by first-order approximation of the non-linear function; therefore, this method is termed a first-order second moment (FOSM) method. Pump Fig. 2. Dome shaped roof slab model The pump penetration shell is one of the critical components of the dome-shaped roof slab (Fig. 2). Hence, it is essential to analyse the load-carrying capacity of the roof slab to ensure structural integrity. The primary objective of the work is to perform reliability-based design optimization of the pump penetration shell (PPS) using a genetic algorithm and compare the load-bearing capacity of PPS made of materials IS2062 and A48P2. It is essential to optimize the design parameters of the PPS for maximum load-carrying capacity with a reduced volume of PPS. Prior to performing the non-linear analysis of the roof slab, a benchmarking study that involves the prediction of collapse load of a simpler structure is carried out. 1 BENCHMARK STUDY ON PREDICTION OF COLLAPSE LOAD FOR SIMPLE STRUCTURE In order to ascertain the methodology to be adopted for performing the non-linear analysis of roof slab, a simpler problem proven experimentally by Zhou et al. [9] is taken up for benchmark study. The problem involves the determination of the limit load of the steel pyramid-to-tube socket connection subjected to uniform compression. Ansys Workbench finite element software is used. The geometric and meshed models of the steel structure are shown in Figs. 3 and 4, respectively. The material properties of the structure are given in Table 1. t Fig. 3. Geometric model of structure Table 1. Material properties of the structure [9] Component Tube Socket Steel Pyramid Steel type STKR400 SS400 Measured plate thickness [mm] 4.2 4.2 Young's modulus [N/mm2] 207893 204076 Yield strength [N/mm2] 378 348 Ultimate strength [N/mm2] 454 433 Uniform elongation [%] 17 17 Initially, the critical load was predicted by performing linear buckling analysis by fixing the bottom end of the pyramid and applying a uniform compressive load at the top, without considering the non-linear behaviour of the structure [10] to [13]. The critical buckling load of the structure is found to be 1867 kN and the corresponding displacement contour is shown in Fig. 5. Fig. 4. Meshed model 332 Gautham, V. - Prabhu Raja, V. - Ebron Shaji, G.T. - Mohanraj, S. - Thyla, P.R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 lÂJWO-^Kifl ■ 1 thiM^hi^militWr'IH B IHn *IMTIlW 111b UH4 OHM jran (.111)1 - ^^^^ ri__ii: ___xarq .-ni Fig. 5. Buckling mode of steel structure (Mode 1) Table 2. Comparison of numerical and benchmark results Description Numerical Benchmark Deviation result result [9] [%] Buckling load [kN] 1867 1910 2.24 The results obtained from numerical analysis are found to be in close agreement with experimental results by Zhou et al. [9], as shown in Table 2. Hence, the above methodology can be extended for the numerical prediction of the collapse load of the PPS. Since the present work deals with optimizing the pump penetration shell for higher load-carrying capacity, a benchmark study that involves the application of reliability approach and genetic algorithm is carried out prior to performing the reliability-based design optimization of the reactor roof slab. 2 A BENCHMARK STUDY OF RELIABILITY-BASED DESIGN OPTIMIZATION A cantilever beam subjected to bending loads is considered for performing reliability-based design optimization. In general, design parameters, such as material properties and geometrical dimensions, as well as boundary and loading conditions, are prefixed while designing any structure. However, in the case of reliability-based design optimization, the limit state function that relates the capacity (R) and demand (S) of the function should be derived as a function 'g of the variables. When g is greater than zero, then the structure is considered to be a safe structure. The limit state function is given by: g = R - S. (1) In reliability-based design, the limit state function g is a random variable because of the uncertainties involved in load and yield strength. The reliability index ft is a direct measure of the reliability of the system and hence a larger value of ft represents higher reliability. For structural applications, the reliability index is generally taken as 3, for which the probability of failure is 0.00135 [14]. A cantilever beam shown in Fig. 6 is taken up for benchmark study for which Wang et al. [15] have prescribed the results of the reliability-based design optimization. The width, height, and thickness of the beam of a length (I) of 3048 mm, are designated as w, h and t, respectively. The loads acting at the free end of the cantilever beam, along x and /-axes are designated as Px and P/, respectively. The moments of inertia of the cantilever beam about x-axis and /-axis are represented as Ixx and Iy/ respectively. In the limit state function, R and S represent the yield strength of material and the maximum bending stress, respectively. The mean and standard deviation of the load acting in x-axis are 2224 N (P) and 667 N (sx) (respectively while the corresponding values for load acting in y-axis are 4448 N (P/) and 222 N (s/). The mean and standard deviations of the material yield strength are 276 MPa (R) and 14 MPa (sR), respectively. The objective is to design the cantilever beam for a minimum weight with the given constraints: ft >3 ; 50 mm < w, h < 254 mm; 2.5 mm < t < 13 mm. Since the system involves more than two random variables, the FOSM method is used to obtain the Fig. 6. Loads acting on the cantilever beam Reliability-Based Design Optimization of Pump Penetration Shell Accounting for Material and Geometric Non-Linearity 333 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 mean (pg) and standard deviation (og) of the limit state function. rlwP lhPy^ V 2lyy 24 Ol = ,dR S Sg ySP y ( „ Y KdPv y y y 5 2 dR 1 dg _ lw VP* = dg _ _lh_ dPy ~ 2Ia ag = ^K. I 5 2 ydR 1 R ( Y a„ ( -, \ V y / (2) (3) (4) (5) (6) (7) (8) The above equations are incorporated into the MATLAB optimization tool, and the objective function is minimized using a genetic algorithm [16] to [18]. The optimization results are compared with that of Wang et al. [15] and are shown in Table 3. must be considered during structural design. The objective is to minimize the weight of the PPS, which is achieved by minimizing the volume of the structure. The volume of the PPS is given by the sum of the volumes of the cylindrical part and the conical part. V and V2 are the volumes of the hollow cylinder and truncated hollow cone, respectively. R1 and R2 are the outer and inner radii of the bottom surface of the cone, r1 and r2 are the outer and inner radii of the cylinder, hy and hc are the heights of the cylinder and the cone respectively. V ={nhy )[r,2 - r22 ], (9) where, r2 = 88; r1 = r2 + t, V = nhc [( + Rr + r,2)-(R22 + R2r2 + r?)], (10) where, R2 = 128; R = R2 + t, v = v + v2. (11) The constraints imposed on the PPS for reliability-based design are given as follows, ft >3; 0° < 6 < 30°; 1 mm < t < 3 mm. For the limit state function (g function) given in Eq. (1), R and S represent collapse load and applied load or operating load (10 kN), respectively. The calculation of collapse load for PPS is discussed in later sections. 3.1 Determination of Collapse Load for PPS 2 2 2 2 2 S y Table 3. Comparison of results based on reliability-based optimization Area [mm2] Thickness Width Height (objective t w h function) [mm] [mm] [mm] Based on genetic algorithm 1770 2.54 169 184 Results by Wang et al. [15] 1772 2.54 179 175 It is inferred from the above table that the results obtained from that predicted by Wang et al. [15] are in good agreement. Hence, the validated methodology for reliability-based optimization is extended for the pump penetration shell, which is described in the next section. 3 RELIABILITY-BASED DESIGN OPTIMIZATION OF PUMP PENETRATION SHELL In order to ensure the reliability of the PPS, the variability in the flange thickness and flange angle 334 The elasto-plastic analysis of PPS is carried out considering two different materials, namely IS2062 and A48P2, which are commonly used in the fabrication of the prototypes of the reactor roof slab. The collapse load behaviour of the PPS made of IS2062 and A48P2 will be compared in this section. The PPS is idealized as a shell body, and it is scaled down to 1:10 ratio from the actual reactor roof slab to facilitate experimentation. The geometry of PPS shown in Fig. 7 is modelled and meshed with SHELL 181 elements. The mesh convergence study was carried out to idealize the element size. Prior to determining the collapse load, it is necessary to determine the buckling mode, which has a considerable influence on the failure of the structure. From experimental investigation on static loading, it was observed that buckling occurred in the conical portion of the structure with imperfection. Hence, Modes 2, 9, 16, 23, and 30, which correspond to buckling in the conical portion of the structure are considered for numerical analysis, both for linear buckling analysis and elasto-plastic analysis. The Gautham, V. - Prabhu Raja, V. - Ebron Shaji, G.T. - Mohanraj, S. - Thyla, P.R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 mode shapes of PPS made of IS2062 are shown in Fig. 8. Similar mode shapes were also observed for the PPS made of A48P2. Subsequently, elasto-plastic analysis was performed, and the collapse load of PPS was predicted by using the elastic slope method twice [19] and [20]. It is inferred from Fig. 9 that there is no considerable influence of mode shape on the collapse load of PPS made of IS2062. In contrast, the influence of mode shape on the collapse load of PPS made of A48P2 is considerable, as shown in Fig. 10. It is found that the collapse load is the lowest corresponding to Fig. 7. Geometric model of pump penetration shell a) Mode 2 b) Mode 9 c) Mode 16 d) Mode 23 Fig. 8. Mode shapes of PPS made of IS2062 e) Mode 30 0.5 0.75 Deformation [mm] Fig. 9. Effect of mode shape on collapse load of PPS made of IS2062 0.5 0.75 Deformation [mm] Fig. 10. Effect of mode shape on collapse load of PPS made of A48P2 Reliability-Based Design Optimization of Pump Penetration Shell Accounting for Material and Geometric Non-Linearity 335 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 Mode 23. It is observed that Mode 23 is characterized by multiple lobes, whereas only two lobes are present in Mode 30. The presence of multiple smaller lobes will offer lower resistance to collapse and vice-versa. Hence, Mode 30 for A48P2 has a higher collapse load than Mode 23. Thus, the above mode is considered for further study on the effect of imperfection on collapse load. The surface irregularities are considered to be imperfections in the PPS. It is found from the analysis that the imperfection below 50 % of shell thickness does not influence the collapse load of the structure. Any imperfections greater than 90 % of shell thickness are not acceptable as per manufacturing standards. Hence, the imperfection levels considered for this analysis are 50 %, 70 %, and 90 % of the thickness of the PPS. It is inferred from Fig. 11 that there is no significant influence of imperfection on collapse load of PPS made of IS2062. However, the influence of imperfection on the collapse load of PPS made of A48P2 is considerable, as shown in Fig. 12. It is inferred from Fig. 9 to 12 that the influence of mode shape and geometrical imperfection on load-deformation characteristics with material A48P2 is minimal. In contrast, the above influence is pronounced with material IS2062, which is attributed to the variation in slope in the non-linear zone of the stress-strain curve. It is also found that the collapse load is the lowest, corresponding to an imperfection level of 90 %, and thus the above imperfection level is considered. In order to determine the collapse load for the generated design points, the buckling mode shape (Mode 23) and the geometrical imperfection level of 90 % are considered for further investigation. 3.2 Experimental Validation of Collapse Load on PPS An experimental investigation was done on the pump penetration shell. The PPS prototype model taken for experimentation has the diameters of the cylindrical part and conical skirt as 180 mm and 264 mm, respectively, with 3 mm thickness as shown in Fig. 7. The PPS was mounted in the fixture, and load was gradually applied on PPS using UTM until it attained permanent deformation. The geometrical shapes of PPS before and after collapse are shown in Figs. 13 and 14, respectively. 0.5 0.75 Deformation [mm] Fig. 11. Effect of imperfection on collapse load of PPS made of IS2062 0.5 0.75 Deformation [mm] Fig. 12. Effect of imperfection on collapse load of PPS made of A48P2 336 Gautham, V. - Prabhu Raja, V. - Ebron Shaji, G.T. - Mohanraj, S. - Thyla, P.R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 Since, the numerically predicted Mode 23 matches with the experimentally observed mode of collapse in the conical region of the PPS specimen, a comparison of the load-deformation characteristics corresponding to Mode 23 is presented in Fig. 15. Considering the numerically predicted load-displacement curve, the collapse load of the structure is found to be 122 kN by using the twice elastic slope method. In contrast, the experimental collapse was identified as 129 kN by the snap-through point in the experimental curve. Hence the deviation in collapse load found from numerical and experimental methods is 6 %. Furthermore, the validated numerical methodology is extended for predicting the collapse load for varying parameters of flange angle and flange thickness and subsequently for reliability-based design optimization of the PPS. 3.3 Design of Experiments The design of experiments is a systematic method to determine the factors affecting a process and the output of that process. This information is needed to manage process inputs in order to optimize the output. Response surface methodology (RSM) is used to formulate the design of experiments for developing a regression equation to predict the collapse load of PPS, which accounts for the relationship between the flange angle and flange thickness. The flange angle is limited to 30° and flange thickness to 3 mm due to manufacturing limitations of PPS. The minimum, medium, and maximum values of the two design parameters are given in Table 4. Table 4. Design optimization parameters Parameter Minimum Medium Maximum (-1) (0) (1) Angle, 0 [deg] 0 15 30 Thickness, t [mm] 1 2 3 The load-deflection curves for all randomly generated experiments are obtained by elasto-plastic Reliability-Based Design Optimization of Pump Penetration Shell Accounting for Material and Geometric Non-Linearity 337 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 analysis considering Mode 23 and imperfection level of 90 %. Mode 23 is considered since the collapse load corresponding to the above mode is the minimum compared to the other modes under consideration. Using the twice elastic slope method, the collapse load is estimated for the design combinations as shown in Table 5. Table 5. Generated design points and corresponding collapse load of PPS Angle Thickness Collapse load [N] [deg] [mm] IS2062 A48P2 0 0 42,500 37,220 -1 5,500 5,500 0 1 71,250 67,220 1 55,000 46,070 1 1 120,000 107,500 0 23,570 20,350 -1 1 33,880 34,000 1 0 95,000 80,000 0 0 42,500 37,220 0 0 42,500 37,220 0 0 42,500 37,220 -1 0 17,140 17,140 Fig. 16. Contour plot of collapse load [N] vs. thickness [mm], angle [deg] of material IS2062 10 IS 20 26 30 Angle [deg] Fig. 17. Contour plot of collapse load [N] vs. thickness [mm], angle [deg] of material A48P2 It is observed that the load-carrying capacity of PPS is increasing with an increase in flange angle and flange thickness, and also follows the same pattern in the contour plot as shown in Figs. 16 and 17 for the PPS made of IS2062 and A48P2, respectively. Considering the flange angle to be constant, it is found that every 0.1 mm increment of flange thickness will contribute to around a 5 % increase in collapse load. Similarly, keeping the thickness constant, every 1° increment in flange angle will contribute to around a 4 % increase in collapse load for PPS made of IS2062 andA48P2. 3.4 Reliability-Based Design Optimization of PPS The random variables considered for design optimization are the flange angle and flange thickness (Table 5), which follow normal distribution within a range of 15°±15° and 2 mm±1 mm respectively. The equations for the load function obtained based on the results of collapse load for the materials IS2062 and A48P2 are given by, R = 43519 + (355800) + (23510?) +(10005<92) + (1345?2) + (91550?), (12) R = 38114 + (308030) + (22008?) +(6649<92) + (3434?2) + (94200?), (13) where 6 and t are the flange angle and flange thickness, respectively. Eqs. (14) and (15) are the limit state equations for the materials IS2062 and A48P2, respectively. £ = [43519 + (355800)+ (23510?) + (1OOO502) + (1345?2) + (91550?) - S, (14) £ = [38114 + (308030) + (22008?) + (664902) + (3434?2) + (94200?) - 5, (15) 3; 0° < 6 < 30°; 1 mm < t < 3 mm), a MATLAB program is written to obtain optimum parameters of PPS that will maximize the collapse load with the minimal possible volume. Eqs. (8) to (11) account for volume, while Eqs. (14) and (16) 2 2 338 Gautham, V. - Prabhu Raja, V. - Ebron Shaji, G.T. - Mohanraj, S. - Thyla, P.R. Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 50 0.25 0.5 0.75 1 1.25 Deformation [mm] Fig. 18. Comparison of collapse load of optimized geometry of PPS pertain to collapse load. The program is executed using a MATLAB optimization tool based on a genetic algorithm. The output of optimization is the optimum values of the flange angle and flange thickness for the PPS as given in Table 6. Table 6. Optimized results for PPS Material IS2062 Material A48P2 Angle [deg] Thickness [mm] Angle [deg] Thickness [mm] 19 16 14 21 It is inferred from Fig. 18 that the collapse load of PPS is 4 and 4.5 times more than the operating load of 10 kN considering IS2062 and A48P2, respectively. The PPS made of IS2062 has a better load-carrying capacity than that made of material A48P2. The flange angle is found to contribute significantly to the load carrying capacity of PPS. Even though the thickness of PPS made of IS2062 is comparatively lesser than that made of A48P2, the collapse load is higher. 4 CONCLUSIONS The reliability-based design optimization of pump penetration based on non-linear static analysis has been carried out, and the resultant optimum values are obtained for the given objective and constraints. The material non-linearity and the geometric imperfection have been taken into account. A genetic algorithm-based optimization technique is adopted for the reliability-based design of the pump penetration shell. A study on the influence of buckling mode in predicting the collapse load reveals that the maximum deviations in collapse load are 1.3 % and 7.6 %, respectively, for PPS made of IS2062 and A48P2. It is inferred that the collapse load of PPS made of IS2062 is not dependent on the buckling mode shape when compared with the PPS made of A48P2. Based on the optimized results, the PPS made of IS2062 is found to have an 11 % higher load-carrying capacity with 21 % less material volume than the PPS made of A48P2. There is a significant improvement in the load-carrying capacity of the optimized geometry of PPS by 54 % when compared to PPS with 0° flange angle made of IS2062 and A48P2. 5 ACKNOWLEDGENENTS The research described in this paper was financially supported by the Board of Research in Nuclear Science (BRNS), Mumbai, India (36(2)/14/39/2014-BRNS/1179). The authors are indebted to Board of Research in Nuclear Science (BRNS), Mumbai and Indra Gandhi Centre for Atomic Research (IGCAR), Kalpakkam, Tamil Nadu for their support throughout the investigation of this project work. 6 NOMENCLATURES g limit state function, [-] R capacity, [-] S demand, [-] l length of the beam, [mm] w width of the beam, [mm] h height of the beam, [mm] t thickness of the beam, thickness of the flange, [mm] Ixx moment of inertia about x-axis, [mm4] Iyy moment of inertia about y-axis, [mm4] Reliability-Based Design Optimization of Pump Penetration Shell Accounting for Material and Geometric Non-Linearity 339 Strojniski vestnik - Journal of Mechanical Engineering 67(2021)6, 331-340 Px load acting on x-axis, [N] Py load acting on y-axis, [N] sR standard deviation of yield strength, [MPa] sx standard deviation of load acting in x-axis, [N] sy standard deviation of load acting in y-axis, [N] H mean, [-] a standard deviation, [-] P reliability index, [-] V1 volume of hollow cylinder, [mm3] V2 volume of truncated hollow cone, [mm3] R1 outer radius of hollow cylinder, [mm] R2 radius of hollow cylinder, [mm] r1 outer radius of truncated hollow cone, [mm] /2 radius of truncated hollow cone, [mm] hy height of hollow cylinder, [mm] hc height of truncated hollow cone [mm] 6 flange angle, [deg] S6 standard deviation of flange angle, [deg] St standard deviation of flange thickness, [mm] 7 REFERENCES [1] Chellapandi, P., Puthiyavinayagam, P., Balasubramaniyan, V., Raghupathy, S., RajanBabu, V., Chetal, S.C., Baldev Raj (2011). Development of innovative reactor assembly components towards commercialization of future FBRs. Energy Procedia, vol. 7, p. 359-366, D0I:10.1016/j.egypro.2011.06.047. [2] Prabhu Raja, V. Ramu, M., Thyla, P.R., Sriramachandra Aithal, Rajan Babu, V., Chellapandi, P. (2016). Structural design optimization of roof slab of a pool type sodium cooled fast reactor. Advances in Engineering Software, vol. 102, p. 97104, D0I:10.1016/j.advengsoft.2016.09.006. [3] Aven, T., Jensen, U. (2013). Stochastic Models in Reliability. Springer, New York. [4] Miura, D., Ishida, Y., Miyasaka, T., Aoki, H. Shinya, A. (2020). Reliability of different bending test methods for dental press ceramics. Materials, vol. 13, no. 22, art ID. 5162, D0I:10.3390/ma13225162. [5] Ormon, S.W., Cassady, C.R., Greenwood, A.G. (2002). Reliability prediction models to support conceptual design. IEEE Transactions on Reliability, vol. 51, no. 2, p. 151-157, D0I:10.1109/TR.2002.1011519. [6] Xin, T., Zhao, J., Cui, C., Duan, Y. (2020). A non-probabilistic time-variant method for structural reliability analysis. Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, vol. 234, no. 5, p. 664-675, D0I:10.1177/1748006X20928196. [7] Zhao, X., Huang, X., Sun, J. (2020). Reliability modeling and maintenance optimization for the two-unit system with preset self-repairing mechanism. Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, vol. 234, no. 2, p. 221-234, D0I:10.1177/1748006X19890739. [8] Cornell, C.A. (1969). A Probability-Based Structural Code. Journal of American Concrete Institute, vol. 66, no. 12, p. 974985. [9] Zhou, Z., Nishida, A., Kuwamura, H. (2011). Applicability of finite element method to collapse analysis of steel connection under compression. Progress in Nuclear Science and Technology, vol. 2, p. 481-485, D0I:10.15669/pnst.2.481. [10] Hornung, U., Saal, H. (2002). Buckling loads of tank shells with imperfections. International Journal of Non-Linear Mechanics, vol. 37, no. 4-5, p. 605-621, D0I:10.1016/S0020-7462(01)00087-7. [11] Rust, W., Schweizerhof, K. (2003). Finite element limit load analysis of thin-walled structures by ANSYS (implicit), LS-DYNA (explicit) and in combination. Thin-Walled Structures, vol. 41, no. 2-3, p. 227-244, D0I:10.1016/S0263-8231(02)00089-7. [12] Schenk, C.A., Schueller, G.I. (2003). Buckling analysis of cylindrical shells with random geometric imperfections. International Journal of Non-Linear Mechanics, vol. 38, no. 7, p. 1119-1132, D0I:10.1016/S0020-7462(02)00057-4. [13] Zhang, J., Zhu, B.Y., Wang, F., Tang, W.X., Wang, W.B., Zhang, M. (2017). Buckling of prolate egg-shaped domes under hydrostatic external pressure. Thin-Walled Structures, vol. 119, p. 296-303, D0I:10.1016/j.tws.2017.06.022. [14] Wang, S., Li, Q., Savage, G.J. (2015). Reliability-based robust design optimization of structures considering uncertainty in design variables. Mathematical Problems in Engineering, vol. 2015, art ID 280940, D0I:10.1155/2015/280940. [15] Wang, W., Wu, J.Y., Antonio, S., Lust, R. V. (1997). Deterministic design, reliability based-design and robust design. MSC Aerospace Users' Conference Proceedings, p. 1-11. [16] Jenkins, W.M. (1991). Towards structural optimization via the genetic algorithm. Computers and Structures, vol. 40, no. 5, p. 1321-1327, D0I:10.1016/0045-7949(91)90402-8. [17] Rani, P., Mahapatra, G.S. (2019). A neuro-particle swarm optimization logistic model fitting algorithm for software reliability analysis. Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, vol. 233, no. 6, p. 958-971, D0I:10.1177/1748006X19844784. [18] Wang, Z., Huang, H.-Z., Liu, Y. (2010). A unified framework for integrated optimization under uncertainty. Journal of Mechanical Design. Transactions of the ASME, vol. 132, no. 5, S. 0510081-0510088, D0I:10.1115/1.4001526. [19] Muscat, M., Mackenzie, D., Hamilton, R. (2003). A work criterion for plastic collapse. International Journal of Pressure Vessels and Piping, vol. 80, no. 1, p. 49-58, D0I:10.1016/ S0308-0161(02)00105-9. [20] Paik, J.K., Seo, J.K. (2009). Nonlinear finite element method models for ultimate strength analysis of steel stiffened-plate structures under combined biaxial compression and lateral pressure actions - Part I: Plate elements. Thin-Walled Structures, vol. 47, no. 8-9, p. 1008-1017, D0I:10.1016/j. tws.2008.08.005. 340 Gautham, V. - Prabhu Raja, V. - Ebron Shaji, G.T. - Mohanraj, S. - Thyla, P.R. Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6 Vsebina Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 67, (2021), številka 6 Ljubljana, junij 2021 ISSN 0039-2480 Izhaja mesečno Razširjeni povzetki (extended abstracts) Jiang Ding, Aiping Deng, Liwei Liu, Mengen Lu: Analiza dinamike linijske zobniške dvojice po numerični mnogoterostni metodi SI 37 Radomir Bokic, Jovan Vladic, Milan Kljajin, Vesna Jovanovic, Goran Markovic, Mirko Karakašic: Modeliranje dinamike, eksperimentalna identifikacija in računalniške simulacije nestacionarnih vibracijpri hitrih dvigalih SI 38 Ignas Sokolnikas, K^stutis Čiuprinskas, Jolanta Čiuprinskiene: Minimizacija stroškov življenjskega cikla rotacijskih prenosnikov toplote za stavbne prezračevalne sisteme v hladnih podnebjih SI 39 Anmol Bhatia, Reeta Wattal: Optimizacija parametrov procesa za povečanje natezne trdnosti ogljikovega jekla, torno varjenega z gnetenjem SI 40 Boopathi Sampath, Sureshkumar Myilsamy: Eksperimentalna raziskava procesa skoraj suhe žične elektroerozijske obdelave v kisikovi meglici s kriogenim hlajenjem žice SI 41 Gautham Velayudhan, Prabhu Raja Venugopal, Ebron Shaji Gnanasigamony Thankareathenam, Mohanraj Selvakumar, Thyla Pudukarai Ramaswamy: Optimizacija konstrukcije prebojnega skoznjika črpalke za zanesljivost ob upoštevanju materiala in geometrijskih nelinearnosti SI 42 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)6, SI 37 © 2021 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2021-02-01 Prejeto popravljeno: 2021-04-28 Odobreno za objavo: 2021-05-14 Analiza dinamike linijske zobniške dvojice po numerični mnogoterostni metodi Jiang Ding1'2'* - Aiping Deng1 - Liwei Liu1 - Mengen Lu1 1 Univerza regije Guangxi, Kolidž za strojništvo, Kitajska 2Univerza regije Guangxi, Državni laboratorij regije Guangxi za proizvodne sisteme in napredno proizvodno tehnologijo, Kitajska Linijski zobniki na podlagi teorije prostorskih ubirnic so drobni in lahki zobniki, ki so zelo primerni za miniaturne stroje. Zanje pa je značilen resen problem vibracij, saj morajo biti ti zobniki za doseganje skladnosti ubirnic oblikovani kot upognjene tridimenzionalne konzole. Za izboljšanje vibracijskih razmer je v članku predstavljen dinamični model linijske zobniške dvojice na osnovi numerične mnogoterostne metode (NMM). NMM je metoda, ki uporablja koncept mnogoterostnega elementa za potrebe matematičnih in fizikalnih modelov. V primerjavi s tradicionalno metodo končnih elementov (MKE) je točnejša in manj občutljiva na deformacije elementov. Najprej je bila izpeljana funkcija premika za linijske zobnike po metodi NMM, nato pa so bile z Lagrangeovo enačbo določene dinamične enačbe mnogoterostnega elementa. Vibracije linijske zobniške dvojice so lahko posledica notranjih in zunanjih dejavnikov. Zaradi poenostavitve analize je privzeto, da linijska zobniška dvojica izpolnjuje vse zunanje pogoje, članek pa obravnava samo notranje dejavnike. Vpliv dušenja je prezrt, glavni notranji povzročitelj težav z vibracijami linijske zobniške dvojice pa je vzbujanje z ubiranjem. Po analizi vzbujanja je bil določen dinamični odziv linijskih zobnikov v vseh treh ortogonalnih smereh. Iz dinamičnega odziva je razvidno, da so vibracije v aksialni smeri zobnika močnejše kot v smeri krivin. Diferencialne enačbe vibracij linijskega zobnika so bile razrešene na podrobnem primeru in razkrita je bila odvisnost med konstrukcijskimi parametri in lastno frekvenco. Nato so bile preučene lastnosti vibracij ob ubiranju zobniške dvojice. Za zmanjšanje vibracij so bili prilagojeni konstrukcijski parametri linijskega zobnika in zmanjšana amplituda vzbujanja. Vibracijske lastnosti linijskega zobnika prvih štirih redov so bile ugotovljene po metodi NMM in nato primerjane z lastno frekvenco, določeno po MKE. Opravljena je bila tudi primerjava napak pri prenosu linijskih zobnikov, izračunanih z metodama NMM in MKE. Izkazalo se je, da je metoda NMM primerna za odpravo težav z vibracijami teh zobnikov. V članku je predstavljena uporaba mnogoterostnega elementa za zmanjšanje vibracij linijskih zobnikov, kakor tudi točnejši teoretični rezultati. Predstavljena teorija NMM bo lahko osnova za raziskave vibracij linijskih zobnikov, uporabna pa bo tudi na sorodnih področjih, npr. pri optimizaciji profila zob, analizi kontaktnih napetosti in analizi upogibnih napetosti v korenu zob. Teorija NMM zagotavlja osnovo za razširitev praktične uporabnosti linijskih zobnikov, s tem pa podpira razvoj industrije zobniških prenosnikov. Ključne besede: linijski zobnik, dinamični odziv, vibracije, numerična mnogoterostna metoda, zobniški prenosnik, drobni zobniki, upognjena tridimenzionalna konzola, metoda končnih elementov *Naslov avtorja za dopisovanje: Univerza regije Guangxi, Kolidž za strojništvo, Kitajska, jding@gxu.edu.cn SI 37 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29 © 2021 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2021-01-26 Prejeto popravljeno: 2021-04-14 Odobreno za objavo: 2021-04-23 Modeliranje dinamike, eksperimentalna identifikacija in računalniške simulacije nestacionarnih vibracij pri hitrih dvigalih Radomir Dokič1* - Jovan Vladic1 - Milan Kljajin2 - Vesna Jovanovic3 - Goran Markovic4 - Mirko Karakašic5 JUniverza v Novem Sadu, Fakulteta tehniških znanosti, Srbija 2Vseučilišče Sever, Univerzitetno središče v Varaždinu, Hrvaška 3Univerza v Nišu, Fakulteta za strojništvo, Srbija 4Univerza v Kragujevcu, Fakulteta za strojništvo in gradbeništvo v Kraljevu, Srbija 5Univerza v Slavonskem Brodu, Fakulteta za strojništvo, Hrvaška Študija predstavlja metodo za oblikovanje dinamičnega modela za analizo nestacionarnih vibracij vrvi s časovno spremenljivo dolžino v anholonomnih robnih pogojih. Obravnavan je odsek vrvi med kabino (kletko) in točko navijanja na vrvenici (bobnu). Predstavljena dinamična analiza in znanstveni pristop sta še posebej pomembna ob upoštevanju narave teh sistemov, ki so namenjeni dviganju (spuščanju) ljudi in bremen na velike višine (globine) s hitrostmi do 20 m/s. Predstavljeni pristop z oblikovanjem diferencialne enačbe za gibanje konca vrvi, ki se navija na vrvenico, pomaga pri določanju kritične hitrosti dvigovanja v odvisnosti od mehanskih lastnosti (modula elastičnosti in dušenja) ter napetosti (obremenitev) v jeklenih vrveh. Za identifikacijo osnovnih parametrov dinamičnega modela jeklenih vrvi, kot so togost, modul elastičnosti in dušenje, je bila uporabljena edinstvena metoda na osnovi eksperimentalnih meritev na konkretnem rudniškem dvigalu. Koeficient dušenja vrvi ni konstanten, ampak je odvisen od položaja kletke. Sklepati je mogoče, da v vrveh deluje kombinacija viskoznega in histereznega dušenja, ki jo bo treba še dodatno raziskati. Med delovanjem dvigal se nenehno povečuje in zmanjšuje prosta dolžina dvižnih vrvi, posledično pa se tudi spreminjata togost in dušenje vrvi. Pri hitrih dvigalih lahko povečanje relativne deformacije privede do dinamičnih nestabilnosti med dvigovanjem (krajšanjem proste dolžine). Te nestabilnosti resno vplivajo na varnost potnikov. Glede na to, da so problemi zaradi vibracij pogonskih mehanizmov predmet večjega števila znanstvenih člankov, v katerih so največkrat uporabljene standardne analize, je logično, da morajo biti v središču pozornosti raziskave dinamike dvigal inovativne metode za analizo longitudinalnih nihanj s spremenljivimi parametri. Eksperimentalni del raziskav je bil opravljen na rudniškem dvigalu s predvideno najvišjo hitrostjo dviganja 16 m/s in dvižno višino 523 m v prvi fazi oz. 763 m v drugi fazi eksploatacije rudnika. Eksperimenti so bili opravljeni z brezžičnim zbiranjem podatkov. S tem je bila dosežena popolna sinhronizacija podatkov, zbranih na strehi kletke in v strojnici. Rezultati eksperimentov in parametri dinamičnega modela so prikazani za štiri različne primere gibanja (dviganja in spuščanja) prazne oz. naložene kletke. Na osnovi rezultatov meritev so bile opredeljene mehanske lastnosti jeklenih vrvi. Iz analize rezultatov sledi sklep, da se podatki o modulu elastičnosti ujemajo s podatki iz literature. To potrjuje veljavnost uporabljenega postopka ter omogoča opredelitev realnih (obratovalnih) vrednosti za rudniška dvigala. Predlagana metoda predstavlja nov pristop, ki bo omogočil analizo dinamičnega vedenja realnih sistemov za vertikalno dviganje. Na osnovi eksperimentalnih rezultatov so bile opravljene simulacije dinamike prazne in naložene kletke. Študija predstavlja tudi posebno metodo za oblikovanje krmilnega programa, ki bi zmanjšal vertikalne vibracije med zagonom in zaviranjem dvigala. Pri hitrih dvigalih je pomembno raziskati spremembe v hitrostih dviganja in opredeliti ustrezno krmiljenje za primere, ko med sosednjima nadstropjema ni dosežena nominalna hitrost, pospeševanju pa takoj sledi zaviranje. Pravilno krmiljenje je ključno za zaviranje v najprimernejšem trenutku, ki zagotavlja bistveno manjšo amplitudo nihanj. Ključne besede: hitra dvigala, dinamična analiza, vrv s časovno spremenljivo dolžino, mehanske lastnosti jeklenih vrvi, vzdolžna nihanja, krmilni program SI 38 *Naslov avtorja za dopisovanje: Univerza v Novem Sadu, Fakulteta tehniških znanosti, Trg D. Obradovica 6, 21000 Novi Sad, Srbija, djokic@uns.ac.rs Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29 © 2021 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2021-01-26 Prejeto popravljeno: 2021-04-14 Odobreno za objavo: 2021-04-23 Minimizacija stroškov življenjskega cikla rotacijskih prenosnikov toplote za stavbne prezračevalne sisteme v hladnih podnebjih Ignas Sokolnikas1 - K^stutis Čiuprinskas2* - Jolanta Čiuprinskiene2 i Salda UAB, Litva 2 Tehniška univerza Gediminas v Vilniusu, Litva Prvotni cilj raziskave je bil določitev konstrukcijskih priporočil za doseganje maksimalne energijske učinkovitosti rotacijskih prenosnikov toplote (RHE) ob upoštevanju izkoristka prenosa toplote in energije, ki je potrebna za premagovanje aerodinamičnega upora pri RHE. Pregled literature je razkril pomanjkanje specifičnih priporočil glede konstrukcijskih parametrov RHE. Večina raziskovalnih projektov je osredotočenih le na simulacije procesov prenosa toplote in snovi ter na izboljšanje njihovega izkoristka. Manjka tudi razprava o merilih za optimizacijo rotacijskih prenosnikov toplote. V večini primerov se obravnava zgolj rekuperacija toplote in vlage, medtem ko so tlačne izgube prezrte. V nekaterih raziskovalnih projektih je bila sicer upoštevana energija za pogon ventilatorjev klimatov, ki je tam okarakterizirana kot »pomemben del« celotne porabe energije, kljub temu pa ni bila vključena v analizo celotne energijske učinkovitosti. Raziskave kažejo, da ni jasnih odvisnosti temperaturnega izkoristka ali tlačnih izgub od enega samega parametra. Prav tako ni smernic za načrtovanje RHE, ali pa so te preveč abstraktne. Cilj raziskave je bil zato opredeljen na novo kot iskanje optimalne množice geometrijskih parametrov RHE ob upoštevanju vračanja toplote v prenosniku ter energije, ki se porabi za transport zraka skozi prenosnik. Narejena je bila ocena količin teh energij v življenjskem ciklu RHE. Dodatna rezultata raziskave sta tudi razpon variabilnosti stroškov življenjskega cikla ter možnost uporabe rezultatov v oblikovanju meril za izbiro RHE. Za določitev temperaturnega izkoristka RHE kot glavne spremenljivke je bil privzet model na osnovi računalniške dinamike fluidov (CFD). Tlačni padec je bil sprva določen na istem modelu CFD, pozneje pa je bil ta zamenjan z analitičnim modelom, ki je skrajšal računski čas. Rezultati modela CFD in analitičnih izračunov so bili eksperimentalno validirani. Skupna količina vrnjene energije, poraba energije ter stroški elektrike in materiala so bili ocenjeni na podlagi podnebja in cen v Litvi. Najpomembnejši parametri, kot so višina in dolžina valov, dolžina RHE, hitrost zraka v RHE, debelina folije in vrtilna hitrost RFE, so bili določeni na osnovi pregleda literature. Stroški življenjskega cikla so vsota stroškov RHE, porabe električne energije ventilatorjev na odvodni in odvodni strani ter poraba toplote na grelniku dovodnega zraka za RHE. Za življenjsko dobo RHE je bila privzeta značilna življenjska doba klimata, ki znaša 10 let. Analiza stroškov življenjskega cikla za 270 različic RHE, izračunanih na osnovi modela CFD in analitičnih modelov, je razkrila najboljše rezultate pri maksimalni dolžini rotacijskega prenosnika toplote (v tem primeru 400 mm), minimalni debelini folije (v tem primeru 0,06 mm) in tlačnem padcu med 100 Pa in 180 Pa. Pričakovani temperaturni izkoristek RHE v hladnih klimatskih razmerah je med 85 % in 90 %. Dobro ujemanje rezultatov izračunov in laboratorijskih preizkusov dokazuje primernost modela CFD za izračun temperaturnega izkoristka rotacijskih prenosnikov toplote. Preskusi dokazujejo tudi primernost analitičnega modela za izračun tlačnega padca na rotacijskem prenosniku toplote. Raziskava je potrdila, da za najboljšo energijsko učinkovitost RHE ni nujen najboljši izkoristek vračanja toplote, razlog za to pa je v tlačnem padcu. Razlika med prenosniki toplote z najnižjimi in z najvišjimi LCC lahko pri enakem podnebju in cenah na trgu doseže 31 %, odvisna pa je predvsem od geometrijskih lastnosti rotacijskega prenosnika toplote. Ključne besede: rotacijski prenosnik toplote, rekuperacija toplote, prezračevalni sistem, temperaturni izkoristek, tlačni padec, ANSYS, stroški življenjskega cikla *Naslov avtorja za dopisovanje: Tehniška univerza Gediminas v Vilniusu, Litva, kestutis.ciuprinskas@vilniustech.lt SI 39 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29 © 2021 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2021-01-26 Prejeto popravljeno: 2021-04-14 Odobreno za objavo: 2021-04-23 Optimizacija parametrov procesa za povečanje natezne trdnosti ogljikovega jekla, torno varjenega z gnetenjem Anmol Bhatia1,2,* - Reeta Wattal2 1 Univerza NorthCap, Oddelek za strojništvo, Indija 2 Tehniška univerza v Delhiju, Oddelek za strojništvo, Indija Jeklo je eden najbolj razširjenih materialov v industriji. Inženirji ga radi izbirajo zaradi dobre kombinacije mehanskih lastnosti in nizke cene. Jeklo se običajno vari s talilnimi postopki, ki pa so pogosto povezani z raznimi težavami, kot so poroznost, toplo pokanje, vodikova krhkost, mikrostrukturne spremembe v toplotno vplivani coni itd. Večina napak pri talilnem varjenju jekel nastane kot posledica taljenja osnovnega materiala. Torno varjenje z gnetenjem je proces, ki odpravlja večino teh napak, saj se obdelovanec med varjenjem ne tali. V predstavljeni študiji je bil opravljen poskus izdelave zvarnih spojev brez napak in z izboljšanimi mehanskimi lastnostmi, in sicer z uporabo orodja iz volframovega karbida. Vhodni parametri imajo pomembno vlogo pri doseganju izboljšanih lastnosti. V dostopni literaturi je nekaj raziskav na temo vpliva parametrov procesa na mikrostrukturo in mehanske lastnosti jekel, torno varjenih z gnetenjem, obstajajo pa še vrzeli v razumevanju mehanskih lastnosti in karakterizaciji mikrostrukture mehkih jekel, torno varjenih z gnetenjem. V pričujoči študiji je bil zato opravljen poskus določitve empiričnih odvisnosti, ki določajo vpliv vhodnih parametrov na natezno trdnost, raztezek pri pretrgu in zoženost pri pretrgu za torno varjen spoj iz materiala AISI 1018. V središču pozornosti je bilo izboljšanje porušitvene natezne trdnosti ogljikovega jekla, torno varjenega z gnetenjem. Obravnavan je bil vpliv parametrov procesa (hitrost varjenja, št. vrtljajev orodja in premer rame) na parametre odziva (porušitvena natezna trdnost, raztezek pri pretrgu in zoženost pri pretrgu). Za razvoj matematičnega modela parametrov odziva je bila uporabljena metoda odzivnih površin, ustreznost modela pa je bila preverjena z analizo variance (ANOVA) za 95 % interval zaupanja. Ugotovljeno je bilo, da hitrost varjenja in število vrtljajev orodja pomembno vplivata na porušitveno natezno trdnost. Raztezek pri pretrgu je bil odvisen samo od hitrosti varjenja. Zoženost pri pretrgu je bila odvisna od hitrosti varjenja in od premera rame. Največja trdnost je bila ugotovljena v zvarni leči, najšibkejši pa je nedotaknjeni osnovni material. Opravljena je bila tudi karakterizacija mikrostrukture zvarnega spoja, ustvarjenega v optimalnih pogojih. V toplotno vplivani coni je bila ugotovljeno formiranje perlita in ferita, v območju mešanja pa so bile opažene rafinirane perlitne strukture. Vrednost mikrotrdote v zvarni leči je bila višja kot v osnovnem materialu zaradi formiranja finozrnatih struktur v območju zvarne leče. Določena je bila tudi morfologija loma preizkušancev za natezni preizkus. Prisotnost jamičastih predelov po celotni površini razkriva duktilni lom. Članek predstavlja poskus optimizacije parametrov procesa za doseganje maksimalne natezne trdnosti zvarnega spoja pri jeklu AISI 1018. Porušitvena natezna trdnost v optimalnih pogojih je bila določena matematično in nato preverjena eksperimentalno. Učinkovitost spajanja pri vseh možnih kombinacijah vhodnih dejavnikov presega 90 odstotkov. Ta vrednost je v primerjavi z rezultati predhodnih študij presenetljivo visoka. Ključne besede: torno varjenje jekla z gnetenjem, ANOVA, metoda odzivnih površin, porušitvena natezna trdnost, morfologija loma, mikrostruktura in mikrotrdota SI 40 *Naslov avtorja za dopisovanje: Univerza NorthCap, Oddelek za strojništvo, Gurugram, Indija, anmolbhatia@ncuindia.edu Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29 © 2021 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2021-01-26 Prejeto popravljeno: 2021-04-14 Odobreno za objavo: 2021-04-23 Eksperimentalna raziskava procesa skoraj suhe žične elektroerozijske obdelave v kisikovi meglici s kriogenim hlajenjem žice Boopathi Sampath1 * - Sureshkumar Myilsamy2 tehniški kolidž Muthayammal, Oddelek za strojništvo, Indija 2Inštitut za tehnologijo Bannari Amman, Oddelek za strojništvo, Indija Članek obravnava postopek skoraj suhega elektroerozijskega rezanja zlitine Inconel 718 v kisikovi meglici s kriogeno hlajeno žično elektrodo (s plinastim dušikom) in napoved optimalnih nastavitev parametrov za najboljšo stopnjo odvzema materiala (MRR) in stopnjo obrabe žice (WWR). Nadzorovane spremenljivke tok, širina impulza, interval impulza in pretok vplivajo na odzivne veličine, kot sta stopnja odvzema materiala (MRR) in stopnja obrabe žice (WWR). Eksperimenti so bili zasnovani po metodi Box-Behnken. Najprej so bila uporabljena načela zaželenosti za napoved optimalnih vrednosti parametrov za najboljšo stopnjo odvzema materiala (MRR) in stopnje obrabe žice (WWR). Napovedi so bile validirane s podatki, zbranimi v potrditvenih eksperimentih. Za preučitev vpliva posameznih dejavnikov, interakcij in kvadratnih členov po metodi Box-Behnken so bili razviti matematični modeli in odzivne površine. Z večciljno optimizacijo je bila določena kombinacija procesnih parametrov za optimalne vrednosti obeh odzivov po načelu zaželenosti. Vrednosti WWR in MRR za predlagano skoraj suho WEDM v kisikovi meglici so bile primerjane s tistimi v konvencionalnem procesu. Analiza odzivnih površin je razkrila naslednje prispevke k MRR: 52,92 % tok, 24,63 % širina impulza, 12,81 % interval impulza in 5,75 % pretok. Prispevki k WWR pa so: 14,89 % tok, 9,75 % širina impulza, 62,20 % interval impulza in 5,44 % pretok. Širina impulza ima večji prispevek k MRR zaradi daljšega trajanja iskre med žico in obdelovanim materialom. Ugotovljeno je bilo tudi to, da ima interval impulza večji vpliv na WWR zaradi bolj idealne periode (dolg čas med iskrama) med dvema zaporednima visokotemperaturnima iskrama. Primerjava s konvencionalnim postopkom WEDM je pokazala nižjo vrednost WWR pri skoraj suhi WEDM. Kriogeno hlajenje žice pri skoraj suhi WEDM pomembno vpliva k zmanjšanju WWR. Literatura navaja, da je postopek skoraj suhe WEDM prijaznejši do okolja kot konvencionalna WEDM. Vrednost MRR pa je pri skoraj suhi WEDM nižja kot pri konvencionalnem procesu WEDM. V raziskavi sta bili najprej eksperimentalno določeni stopnja obrabe žice (WWR) in stopnja odvzema materiala (MMR) pri skoraj suhi WEDM v kisikovi meglici s kriogeno hlajeno žico. Obraba žice je bila analizirana z vrstično elektronsko mikroskopijo (SEM). Pri vrednosti toka 5 A, širini impulza 25 ^s, intervalu impulza 75 ^s in pretoku 1 ml znaša vrednost MRR 10,97 mm3/min, vrednost WWR pa 0,59 %. Rezultati večciljne optimizacije bodo uporabni za proizvajalce in operaterje pri določanju najboljših (privzetih) nastavitev za proces NDWEDM v kisikovi meglici s kriogenim hlajenjem žice. Ključne besede: kriogeno hlajenje, kisikova meglica, skoraj suho, žična elektroerozija, MRR, WWR, metoda Box-Behnken *Naslov avtorja za dopisovanje: Tehniški kolidž Muthayammal, Oddelek za strojništvo, Rasipuram-637408, Namakkal (Dt), Indija, boopasangee@gmail.com. SI 41 Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29 © 2021 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2021-01-26 Prejeto popravljeno: 2021-04-14 Odobreno za objavo: 2021-04-23 Optimizacija konstrukcije prebojnega skoznjika črpalke za zanesljivost ob upoštevanju materiala in geometrijskih nelinearnosti Gautham Velayudhan* - Prabhu Raja Venugopal - Ebron Shaji Gnanasigamony Thankareathenam -Mohanraj Selvakumar - Thyla Pudukarai Ramaswamy Oddelek za strojništvo, Tehniški kolidž PSG, Indija Strešna plošča jedrskega reaktorja mora nositi več komponent in podsistemov. Prevzemati mora vse seizmične obremenitve v skladu z zahtevami nosilnosti. Analiza statičnih napetosti v strešni plošči reaktorja je razkrila visoko koncentracijo napetosti v prebojnem skoznjiku črpalke (PPS), ki nosi glavno črpalko za natrij. PPS je torej ena od kritičnih komponent strešne plošče. Glavni cilj procesa konstruiranja take komponente je zagotovitev varnosti in ekonomičnosti delovanja v danem okolju z visoko zanesljivostjo. Zanesljivost je dejavnik upravljanja negotovosti in predstavlja verjetnost, da bo element v danih pogojih brez odpovedi opravljal svojo funkcijo v načrtovanem časovnem obdobju. Članek obravnava optimizacijo konstrukcije prebojnega skoznjika črpalke za zanesljivost ob upoštevanju materiala in geometrijskih nelinearnosti. Debelina in kot prirobnice PPS kot konstrukcijska parametra sta bila optimizirana za maksimalno nosilnost ob manjši količini materiala. Najprej sta bili opravljeni dve primerjalni študiji, prva za napoved porušne obremenitve preprostejše konstrukcije ter druga na osnovi zanesljivosti in genetskega algoritma za optimizacijo konzolnega nosilca na osnovi zanesljivosti. Za zanesljivost PPS je bila v analizah obravnavana variabilnost debeline in kota prirobnice. Opravljena je bila elasto-plastomehanska analiza PPS po metodi končnih elementov za materiala IS2062 in A48P2, ki se pogosto uporabljata pri izdelavi prototipov strešnih plošč reaktorjev. Narejena je bila primerjava porušnih obremenitev PPS iz omenjenih materialov ob upoštevanju geometrijskih nepopolnosti. Eksperimenti so pokazali, da prihaja do uklona v koničnem delu konstrukcije z nepopolnostmi. V numeričnih analizah linearnega uklona in v elasto-plastomehanskih analizah so bila zato obravnavana stanja, ki ustrezajo uklonu koničnega dela konstrukcije. Numerična analiza je zajela oblike, ki izkazujejo opisano vedenje. Validirana numerična metoda je bila razširjena za napovedovanje porušne obremenitve pri variabilnih vrednostih kota in debeline prirobnice in nato za optimizacijo konstrukcije PPS za zanesljivost. Zasnova eksperimentov je bila oblikovana po metodi odzivnih površin (RSM) za razvoj regresijske enačbe, ki napoveduje porušno obremenitev PPS ob upoštevanju odvisnosti med kotom in debelino prirobnice. Porušna obremenitev in prostornina sta bili uporabljeni kot vhodni funkciji za optimizacijo na osnovi zanesljivosti z omejitvami. Izdelan je bil program v okolju MATLAB za določanje optimalnih parametrov PPS, ki zagotavljajo maksimalno porušno obremenitev ob minimalni količini materiala. Program je bil izvršen z optimizacijskim orodjem MATLAB na osnovi genetskega algoritma. Porušna obremenitev PPS na osnovi optimiziranih parametrov za IS2062 in A48P2 je od 4- do 4,5-krat večja od obratovalne obremenitve. Ugotovljeno je bilo, da porušna obremenitev za PPS iz materiala IS2062 ni odvisna od uklonske oblike. Kot prirobnice pomembno prispeva k nosilnosti PPS. Čeprav je debelina PPS iz materiala IS2062 manjša od debeline PPS iz materiala A48P2, je porušna obremenitev prvega PPS večja. Po rezultatih optimizacije ima PPS iz materiala IS2062 za 11 % višjo nosilnost in za 21 % manj materiala kot PPS iz materiala A48P2. Z optimizacijo geometrije PPS je bilo pri teh materialih doseženo 54-odstotno izboljšanje nosilnosti v primerjavi s PPS s kotom prirobnice 0°. Uporabljeno metodologijo bo mogoče razširiti za optimizacijo konstrukcije drugih kritičnih delov v jedrski industriji za zanesljivost. Ključne besede: jedrski reaktor, uklon, optimizacija, zanesljivost, mejna obremenitev, genetski algoritem SI 42 *Naslov avtorja za dopisovanje: Tehniški kolidž PSG, Peelamedu, Coimbatore, Indija, 1701rm01@psgtech.ac.in Guide 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 12 pages (approx. 5000 words). Longer contributions will only be accepted if authors provide justification 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. One corresponding author should be provided. 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. Kordič, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553576. Proceedings Papers: Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages. [4] Štefanič, N., Martinčevič-Mikič, S., Tošanovič, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427. Standards: Standard-Code (year). Title. Organisation. Place. [5] ISO/DIS 16000-6.2:2002. Indoor Air - Part 6: Determination of Volatile Organic Compounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva. WWW pages: Surname, Initials or Company name. Title, from http://address, date of access. [6] Rockwell Automation. Arena, from http://www.arenasimulation.com, accessed on 200909-07. EXTENDED ABSTRACT: When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters or approx. 600 words). The instruction for composing the extended abstract are published on-line: http://www.sv-jme. eu/information-for-authors/ . COPYRIGHT: Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to transfer copyright to SV-JME when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Transfer Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on http://en.sv-jme.eu/. PUBLICATION FEE: Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 380 EUR (for articles with maximum of 6 pages), 470 EUR (for articles with maximum of 10 pages), plus 50 EUR for each additional page. The additional cost for a color page is 90.00 EUR (only for a journal hard copy; optional upon author's request). 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 Jiang Ding, Aiping Deng, Liwei Liu, Mengen Lu: Dynamic Analysis of Line Gear Pair Based on Numerical Manifold Method Radomir Bokic, Jovan VLadic, Milan KLjajin, Vesna Jovanovic, Goran Markovic, Mirko Karakasic: Dynamic Modelling, Experimental Identification and Computer Simulations of Non-Stationary Vibration in High-Speed Elevators Ignas SokoLnikas, Kestutis Ciuprinskas, JoLanta Ciuprinskiene: Minimization of the Lifecycle Cost of a Rotary Heat Exchanger Used in Building Ventilation Systems in Cold Climates AnmoL Bhatia, Reeta WattaL: Process Parameters Optimization for Maximizing Tensile Strength in Friction Stir-Welded Carbon Steel Boopathi Sampath, Sureshkumar MyiLsamy: Experimental Investigation of a Cryogenically Cooled Oxygen-mist Near-dry Wire-cut Electrical Discharge Machining Process Gautham VeLayudhan, Prabhu Raja VenugopaL, Ebron Shaji Gnanasigamony Thankareathenam, Mohanraj SeLvakumar, ThyLa Pudukarai Ramaswamy: Reliability-Based Design Optimization of Pump Penetration Shell Accounting for Material and Geometric Non-Linearity 275 287 9770039248001