Journal of Mechanical Engineering no. 9 year 2022 volume 68 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škerceva 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 Mihael Sekavcnik 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: Using smaller and smaller rulers the coast length limits to infinity. If this logic is applied to the fractal heat sink geometry, infinite cooling capacity should be obtained using fractals with mathematically infinite surface area. On the figures the formation of Koch snowflake geometry from flat line to 1024 elements on the 5th step is shown. Using colours, the temperature field of cooling air flowing from the left-hand side is plotted. One should note the similarity of geometry and temperature development up to the last figure, where the changes become minor. Image Courtesy: Matjaž Ramšak, University of Maribor, Faculty of Mechanical Engineering, Slovenia © The Author, CC BY 4.0 Int. 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 Lübben, University of Bremen, Germany George K. Nikas, KADMOS Engineering, UK Tomaž Pepelnjak, University of Ljubljana, Slovenia Vladimir Popovic, 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 Vasic, 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. ISSN 0039-2480, ISSN 2536-2948 (online) We would like to thank the reviewers who have taken part in the peer-review process. © 2022 with Authors. The journal is subsidized by Slovenian Research Agency. 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. Strojniški vestnik - Journal of Mechanical Engineering is available on https://www.sv-jme.eu. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9 Contents Contents Strojniški vestnik - Journal of Mechanical Engineering volume 68, (2022), number 9 L jubljana, September 2022 ISSN 0039-2480 Published monthly Papers Matjaž Ramšak: Fractal Geometry as an Effective Heat Sink 517 Youyu Liu, Liteng Ma, Siyang Yang, Liang Yuan, Bo Chen: MTPA- and MSM-based Vibration Transfer of 6-DOF Manipulator for Anchor Drilling 52 Kuat Kombayev, Murat Muzdybayev, Alfiya Muzdybayeva, Dinara Myrzabekova, ojciech ieleba, Tadeusz Leniewski: FunctionalSurface Layer Strengthening and ear Resistance Increasing of a Low Carbon Steel by Electrolytic-Plasma Processing 542 Mateusz rzochal, Stanisaw Adamczak, Grzegorz Piotrowicz, Sylwester nuk: Industrial Experimental Research as a Contribution to the Development of an Experimental Model of Rolling Bearing Vibrations 552 Yanbing Ni, enliang Lu, Shilei Jia, Chenghao Lu, Ling hang, Yang en: Limit-protection Method for the orkspace of a Parallel Power Head 560 iaolin Huang, Chengzhe ang, Jiaheng ang, Nengguo ei: Nonlinear Vibration Analysis of Functionally Graded Porous Plates Reinforced by Graphene Platelets on Nonlinear Elastic Foundations 571 Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528 Received for review: 2022-01-24 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-04-25 DOI:10.5545/sv-jme.2022.28 Original Scientific Paper Accepted for publication: 2022-06-13 Fractal Geometry as an Effective Heat Sink Matjaž Ramšak1 1Faculty of Mechanical Engineering, University of Maribor, Slovenia "HowlongisthecoastofBritain?"wasthequestionstatedbyMandelbrot.Using smallerand smallerrulersthecoastlengthlimitstoinfnity.Ifthislogicis applied to the fractal heat sink geometry, infnite cooling capacity should be obtained using fractals with mathematically infnite surface area. The aim of this article is to checkthis idea using Richardson extrapolation of numerical simulation results varying the fractal element length from one to zero. As expected, theextrapolated cooling capacityhasanoninfnite limit.The presentedfractalheatsink geometryis non-competitivetothestraightfns. Keywords: fractal heat sink, light-emitting diode and central processing unit cooling, conjugate heat transfer, laminar fow, boundaryelement method,Kochsnowfake Highlights • Infnitefractalheattransfer area leadstozero Nusseltnumber.Their product limitsto fnitevalueofheatfow. • Richardsonextrapolationofnumericalsimulation results confrm that. • Fractalfowpatternfollowsfractal geometry. • Fractal heat sink could not compete with straight fns. • Conjugate heat transfer is simulated using boundaryelement method inhouse code. 0 INTRODUCTION The convective heat transfer per unit time Q.between a heat sink and the fuid is computed as Q.= hA˜T ,(1) where h is the heat transfer coeffcient, A the heat transfer area, and ˜T the temperature difference between the solid surface and fuid free stream. If the heat transfer area tends to infnity, the heat transfer power should tend to infnity too, presuming that the heat transfer coeffcient and ˜T are fnite nonzero values. The aim of this research is to test this assumption using computational fuid dynamics (CFD).Thetestis performedusingasequenceofKoch snowfakefractal formation starting fromastraight line to some small fnite value, as shown in the last fgure in the article. Using Richardson extrapolation [1], the results are extrapolated to the infnite surface area. If the result of this study confrms the idea, a very effective coolingdevice shouldbegained using fractal geometry. The fractal geometry is encountered in nature in many shapes and purposes. For example, the blood veins with capillaries in lungs for heat and mass exchange. It is well known that the nature evolution solutions are superior to engineering ones in many areas. If this is true for heat sinks, then a fractal heat sink must exist to be better than simple one with the straight fns which is used nowadays in most cases. As itisapopularslogan: thereisa roomforimprovement. The geometry of heat sinks is subject of many recent articles [2–6]. The fractals are not an unknown topic in this subject. In most papers, fractal geometry is used to characterise the surface roughness [7], where the correlation between the critical heat fow and the fractal surface roughness of the pressure tubes from the cool heavy water moderator is investigated. A3D laminarfowinamicrochannelis analysed numerically in [8], where near wall swirls are obtained of size 1 µm. In the minority group of papers, fractal geometry is used to describe the material porosity, for example [9–12]. Elaborate review article of fractal heat exchangers is written in [13] where the list of similar articles is presented. Almost all of them are numerical simulations. The numerical feasibility study of fractal heat sink is published in [14], where fractal like heat sink is proposed for cooling of electronic device. The conjugate heat transfer in a fractal tree like channels network heat sink is studied numerically and experimentally [15]. A conclusion is made that a fractal heat sink has lower pressure drop, more uniform temperature feld distribution, and higher coeffcient of performance than that of the traditional helical channel net heat sink. In the latest review article [16] on optimization design of heat sinks the fractal geometry is mentioned only in one cited article [17]. A simple heat sink is simulated consisting of a fractal split microchannel up to second iteration formation only. The challenge of high cooling power in electronic devices and its dissipation using bionic Y-shaped fractals is the subject of very recent work of He et al. [18]. TheKoch snowfake fractal geometry, again only up to the second iteration formation, is found in [19], where it is applied for micro mixer baffes geometry. The laminar fow is solved up to Re = 100, coincidentally the same as in our work. Table 1. Results sensitivity on mesh densityfor Re = 1and solid fuid heat conductivityratio k = 1 Mesh data Mesh name coarse medium fne fnest Number of Nodes [*1000] 47 92 231 457 Nodes on Koch shortest element 3 3 5 7 Nodes on Outlet 90 180 226 300 Resultsfor isothermal computation Number of iterations 940 880 955 940 CPU [h] (serial run) 14 36 100 351 AVG(Interface vorticity) -2.872 -2.779 -2.682 -2.662 Error to coarser mesh [%] - 3.35 3.62 0.75 Num. acc. estimation [1] [%] - - 1.43 0.50 Resultsfor thermal computation AVG(Interface temp.) 0.6536 0.6533 0.6529 0.6528 Error to coarser mesh [%] - 0.05 0.06 0.02 Num. acc. estimation [1] [%] - - 0.04 0.02 Nusselt number 0.5173 0.5206 0.5211 0.5218 Error to coarser mesh [%] - 0.62 0.10 0.13 Num. acc. estimation [1] [%] - - 0.08 0.16 The conjugate heat transfer in this work is computed using in house code based on mixed boundary elements and subdomain boundary element method (BEM). The idea of mixed boundary elements [20] is to split the function and fux approximation using continuous interpolation polynomials for function and discontinuous for function derivative in normal direction to the boundary element. In this manner the problem of undefned normal direction on the corner fux nodal points is elegantly avoided. The main advantage of subdomain technique is sparse matrix in comparison to the classic BEM, where only the boundary of computational domain must be discretised. The subdomain technique in its limit version by [21], where each subdomain is consisted of three or four boundary elements as triangle or quadrilateral subdomain, resulted in extremely sparse system matrix like in fnite element method (FEM). The interface boundary conditions between mixed elements of subdomains lead to overdetermined matrix, which is solved using fast iterative least squares method, [22]. The code has been successfully used and validated for the conjugate heat transfer Benchmark revision [23]. The paper is organised as follows. The Problem defnition is stated after the Introduction. In the next section Results and Discussion, the main subsection is the last one titled The infnite heat fow idea. Prior to this, various tests are performed for the shortest simulated fractal length (1/3)5 = 0.004 in the ffth iteration consisting of 45 = 1024 elements: mesh sensitivity, Reynolds dependencyfor isothermal solution, infuence of solid/fuid thermal conductivity ratio and infuence of Reynolds number value on thermal solution. The paper fnishes with concluding remarks on enhanced heat sinks using fractal geometry. 1 PROBLEM DEFINITION The geometry of the Koch snowfake represents the frst cut out of the heat sink of the high heat density source, such as light-emitting diode (LED) or central processing unit (CPU) processor, as shown in Fig. 1. The bottom of the solid wall is heated to a constant temperature. The cooling fuid is fowing into the domain with constantvelocity. The zero gradient outlet boundary condition is prescribed. It is not physically adequate sincethefow feldisnotfullydeveloped,but the obtained fow feld is as expected at the outfow region and serves well for the numerical example aim. The buoyancy effect is neglected. This could be justifed by small fuid solid temperature difference or orienting the gravity in the third dimension not infuencing the fow feld in the 2D computational cross-section shown in Fig. 1. The problem is nondimensionalized as follows. Length quantities are nondimensionalized by the length of the computational domain in the mean fow direction. The velocity is nondimensionalized by the inlet velocity. The Reynolds number is computed as Re = Velocity@inlet ×Length/˜. The steady laminar fow of air is presumed as a cooling fuid. The Prandtl number is set to 0.71. The bottom dimensionless temperature is 1.0 and inlet temperature is zero. In this manner the temperature difference ˜T is defned to be 1in Eq. (11). Computing the steady state solution, the fuid solid thermal diffusivity ratio is the last solution parameter investigated in next sections. The non-dimensional form of governing equations for a 2D incompressible laminar fow are written using the nondimensional stream function vorticity formulation of Navier-Stokes equations. Stream function equation °is .2°.2° .(2) .x2 + .y2 = -. Inlet:˜v =(1,0)Outlet: , °˜ °n =0, ˜ =y, °.°. °T T =0. =0. °n 0.9 Interface: ˜v =(0,0), =0, °˜ ˜ °n =0, Left/Right solid wall: °T °n =0. Ts =Tf , °Ts °Tf ks =kf 0.1 . °n °n .˜ Lower wall: T =1. Fig. 1. The geometryof computational domain and boundaryconditions Table 2. Results sensitivity on solid fuid conductivityratio k for Re =1 Solid fuid conductivity ratio k 1 10 100 1000 104 105 AVG(Interface temperature) 0.6529 0.9175 0.9900 0.9990 0.9999 1.0000 Change to prior -lower k [%] - 28.84 7.32 0.90 0.09 0.01 AVG(Outlet temperature) 0.4030 0.5236 0.5526 0.5562 0.5565 0.5565 Change to prior -lower k [%] - 23.03 5.25 0.65 0.05 0.00 AVG(Nusselt number) 0.5211 1.1167 1.3659 1.3994 1.4033 1.4036 Change to prior -lower k [%] - 53.33 18.25 2.39 0.28 0.02 Heat fow .Q [W] 0.0309 0.0663 0.0811 0.0831 0.0833 0.0833 Change to prior -lower k[%] - 53.33 18.25 2.39 0.28 0.02 Heat transfer coef. h [W/m2K] 0.0073 0.0157 0.0192 0.0197 0.0198 0.0198 Change to prior -lower k[%] - 53.33 18.25 2.39 0.28 0.02 Table 3. Results sensitivityonReynoldsnumberfor k =10 where T is non-dimensional fuid temperature. Energy equation within the solid region is °T °t ˜.s °= .f 1 ˜2T ,RePr (5) where .s and .f are diffusivities for the solid and fuid regions respectively. For details on equations derivation [23]. The mechanism of heat conduction andits backgroundis clearlyexplainedinLiuetal.[9]. The interface boundary conditions on the wall Reynolds number Re 1 10 100 AVG(Interface vorticity) -2.682 -2.777 -3.629 AVG(Interface temperature) 0.9175 0.9054 0.8471 AVG(Outlet temperature) 0.5236 0.4766 0.3331 Nusselt number 1.1167 1.2231 1.7366 Heat fow .Q [W] 0.0663 0.0726 0.1031 Heat transfer coef. h [W/m2K] 0.0157 0.0172 0.0245 Vorticity equation.is between solid and fuid are written as temperature equality °.°(vx.) °(vy.) 1 ++= ˜2.,(3) Tf =Ts ,(6) °t °x °y Re and heat fux equality as where vx is thevelocityin x direction computed as vx = °˜/°y and vy as vy =-°˜/°x. The energy equation ˜°T °˜°T ° kf =-ks ,(7) within the fuid region is °n °n f ,interface s,inter f ace °T °(vxT ) °(vyT ) 1 where the n is unit normal direction to the fuid solid ++= ˜2T ,(4) interface. The local Nusselt number is defned as the °t °x °y RePr CoarseN = 47,000 FinestN = 457,000 Fig. 2. Meshesused; Increasingthemeshdensitythenumberof outletnodesare90,180,226and300whileontheshortestfractal elementoftheKoch boundaryare3,3,5and7nodes Table 4. Average values at the Koch fuid solid interface for fractal geometry sequence; the fractal element length is denoted using l and the complete interface length using L; the Richardson extrapolation [1] to infnite boundaryis presented in the line l =0including extrapolation accuracy estimation interval. In the last line, the best guestimate values are stated l L ˜Koch TKoch Nu .Q h 1.000 1.000 -18.58 0.9999 8.158 0.1149 0.1149 0.333 1.333 -15.52 0.9998 6.071 0.1140 0.0855 0.111 1.778 -8.958 0.9998 4.980 0.1247 0.0701 0.037 2.370 -6.397 0.9998 3.954 0.1320 0.0557 0.012 3.160 -4.578 0.9998 3.051 0.1358 0.0430 0.004 4.214 -3.629 0.9998 2.332 0.1384 0.0328 0.0 ˜ -2.518 0.9998 -0.899 0.1446 -0.0127 acc. est. ± 1.0 0.0000 1.6 0.0081 0.0219 guestimate -2.518 0.9998 0 0.1446 0 temperature derivativeat the fuid side of the fuid solid interface as °T Nu =-.(8) °n ˜ f ,inter f ace The average Nusselt number Nu is the integral value computed as 1 °L Nu = Nu dl ,(9) L 0 where L is the interface wall length. The heat fow Q .is computed using its defnition as ..°T . Q =-Akf =(L ·1)k f Nu ,(10) °n f whereAisthe actual interface area computedas L ·1. The unity length is defned on the 3rd dimension. Heat transfer coeffcient h is computed as . Q h = (11) A°T where °T is one in this case. In this paper only a steady solution is computed. The steady solution is solid thermal diffusivity .s independent since the time derivative is zero, see Eq. (5). In this manner the solid fuid thermal conductivity ratio k defned as k =ks/kf is the only solid material Fig. 3. The Reynolds number dependency of the stream function and vorticity contour plots parameter arising from interface boundary condition Eq. (7), infuencing the solution. All governing equations can be written in the same general form and solved using practically the same multidomain BEM solver applying different boundary conditions for each governing equation. The solver is explained in a detail in [24]. The validation of the developed multidomain BEM solver [23] where the benchmark solution of conjugate heat transfer of backwardfacing step problemis computed. 2 RESULTS AND DISCUSSION The basic aim of the article is to check the assumption about the infnite heat fow for infnite lengthoffractal solid-fuid interface. Before this main numerical test, a few necessary tests are performed using the fnest fractal geometry which is numerically the most cumbersome to solve. 2.1 Mesh Sensitivity Study The aim of this test is to choose the appropriate mesh density and numerical solution accuracy estimation using standard procedure described in [1]. Four mesh densities were used: coarse, medium, fne and fnest, see Fig. 2. In Tab. 1 the results are shown for the selected case Re = 1and solid fuid conductivity ratio k = 1. Three integral values are selected as numerical solution accuracy indicators: average (AVG) value of vorticity, temperature and Nusselt value on the solid fuid interface. The AVG (Interface vorticity) Fig. 4. Geometryandfow (stream function)self similarityatReynoldsnumbers1and100 numerical solution accuracy is the worst among all, being 0.50 %, which is to be expected since the interface vorticity is the most nonlinear and therefore diffcult to solve. The thermal solutions are less mesh sensitive, resulting in maximal 0.02% error for average temperatures and 0.16% error for the Nusselt number. Obviously and intuitively, the temperature profle over the interface is less mesh dependent than the vorticity one. Based on the CPU consumption and basic aim, the fne mesh was chosen as a default mesh for all further computations resulting in less than 1 % numerical solution accuracy. Similar numerical accuracyis published in the work [23] where using the same BEM code the benchmark conjugate heat transfer problem was solved. 2.2 Reynolds Number Dependencyof Isothermal Solution The contour plots are presented in Fig. 3 for stream function and vorticity feld at Re = 1 and Re = 100. Close to the Koch snowfake boundary, the streamlines in Figs. 3 are almost symmetrical with respect to the upstream and downstream sides for Re = 1. For a higher Re number the symmetry feature is altered, since the recirculation zone appears. In contrast to Re = 1, in general, the bright blue colour representing positive stream value change to dark blue colour representing negative stream value on the downstream side. The dark blue zones represent the clockwise and light blue counter clockwise rotation vortices. The symmetry issues are more evident in the smallest cavities. This feature is also visible on vorticity contour plots as a dark and bright red colour, representing the zero-value vorticity contour. The fow self-similarity is discussed in the next paragraph. Streamlines are shown in Fig. 4, where two successive zooms are enlarged on the right-hand side of the full-scale plot. Zooms 1 and 2 were selected in such a manner that the geometric self-similarity is evident.Thefow patternself-similarityisobvioustoo, especially for Re = 1. The fow pattern is discussed in manyaspects in our latest work completely dedicated to the isothermal solution [25]. 2.3 Infuence of Solid Fluid Conductivity Ratio k The real solid fuid conductivity ratio k values are 16000, 9000 and 2300 for Cu, Al and steel heat sink material, respectively. The test values are chosen to be in the range from 1 to 105. The Nusselt Fig. 5. Infuenceof solid fuid conductivityratio k on temperature feldfor Re = 1 Fig. 6. InfuenceofReynoldsnumberon temperature feldfor k = 10 values are not changing signifcantly (only 2.39 %) when increasing k over 100, see Tab. 2 and Fig. 5, indicating no signifcant changes. Increasing k, solid gradients decrease obtaining an almost constant solid temperature. Replacing the steel with Al the cooling heat fow increases for 0.3%in this problem confguration. The neglecting improvementof 0.02% is obtained replacing Al with Cu. Interesting temperature contours are obtained for k = 1, equalling the heat conductivity in fuid and solid, see Fig. 5. In this case (Re = 1), the solid domain infuence on temperature nearlyvanishes, since the conduction equals convection, obtaining large temperature gradients in the solid domain. The Biot number (Bi = Length ·h/ks) analysis follows. The characteristic Length is defned in this case as Length = A/L where the A is fnite heat sink cross section approximately 0.1·1 and L the length of the fractal cooling surface area. Increasing L to the infnity the characteristic length and Bi limits to zero. Additionally, increasing k and consequently ks values the Bi number also limits to zero value. Both parameters clearly indicating very low Bi values and consequently the uniform solid temperature as already mentioned. Using results in Tab. 2 the maximum Bi number value is Bi = 0.0002 using L = 4.214, Length = 0.0237, h = 0.0073 and k = 1, confrming that heat conduction in solid prevails heat convection to fuid. 2.4 Infuence of Reynolds Number on the Thermal Solution For this test, the fuid solid conductivity ratio k value is fxed to 10 in order to obtain temperature changes in the solid. In this academic case the heat sink is made using isolative material such as wood or plastic. Increasing the Re number, the heat convection prevails over diffusion, resulting in nearly linear growth of the Nusselt number values, see Tab. 3 and Fig. 6. One should expect that increasing Re and consequently the cooling heat rate, the outlet temperature should increase too. Wrong, the outlet temperature decreases since the mass fow 0.15 0.145 0.14 0.135 hQ . 0.13 0.125 0.12 0.115 0.11 Fractal element length(l) Fractal element length(l) 0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 Fractal element length(l) Fractal element length(l) Fig. 7. Graphical presentation of tabulated results in Tab. 4; the shaded area is accuracy estimation of extrapolated result value computed using [1] increases. The proper thermodynamic conclusion in this case would be: the enthalpyof outfow, computed as mass fow and temperature product, increases and exergy decreases. Increasing Re number the interface temperature decreases too, indicating higher temperature gradients and consequently, higher heat fux in the solid domain. 2.5 The Infnite Heat Flow Idea The infnite heat fow ideais testedby formingKoch snowfake fractal geometry starting from the fat heat sink geometry denoted by l = 1, where l is the fractal element length. The next geometry iteration is obtainedbydividing each fractal elementbyfactor3 and creating4new elements using theKoch snowfake formation procedure, see Fig. 8. In this manner, the Koch boundary length is increased by factor 4/3 at each formation iteration giving an infnite limit. The fnal boundary length in this research is (4/3)5=4.214 using 5th iteration formation and fractal element size l =(1/3)5 =0.004. Using Richardson extrapolation [1],the numerical resultvalues limitis computedusing l as mesh size parameter which tends to zero. The accuracyestimation is also a part of the extrapolation procedure results. The infnite heat fow assumption is tested using Re =100 and solid fuid conductivity ratio k =104 which matches the Aluminium heat sink material approximately. In the Figs. 8, 9 and 10 the contour plots of stream function, vorticity and temperature are plotted respectively for each fractal element length l. It is interesting and natural to expect that the fow felds have the same fractal nature as the geometry has in the Koch snowfake formation procedure. This is more evident in the upwind side of the heat sink in comparison to the downwind side where the recirculation zone is present. The quantitative results are presented in Tab. 4 and graphically in Fig. 7. Observing the heat fow Q .dynamics the smooth response is obtained generally. The only exception is in the frst iteration between l =1.0andl =0.333. The frst questionwouldbehow it is possible, that the cooling heat fowis slightly lower for 0.8% usinga single rib(l =0.333) comparing to the fat surface(l =1.0). While the surface area Fig. 8. Stream function contoursfor theKoch snowfakeformation procedure denotedbyfractal element length l;the fow feld have the same fractal nature as the geometry of ribbed surface A is increased, the heat transfer coeffcient h is signifcantly decreased resulting in slightly lower heat fow which is their product Q .= h ·A ·˜T. In this case the cooling rib is more of a fuid fow and thermal obstacle than cooling enhanced as expected intuitively. Additionally, the recirculation zone performs thermalisolation increasing the thermal boundary layerthicknessindownwind area comparing to the fat surface case, see Fig. 10. The next objectivity of discussion is to verify that the fve iterations of the Koch snowfake formation are enough for the testing aim. Comparing results in the last two fgures, namely l = 0.012 and l = 0.004, the macro fow solution does not change any more signifcantly. Decreasing l further approaching the roughness sizevalue, the fowwould change onlyvery close to the boundary until the viscous forces would damp the smallest swirls in cavities limiting to the hydraulic smooth surface fow. If the dimension of the heat sinkwouldbe1 cm, then the shortest fractal element length l = 0.004 would be 40 µm, which is equivalent to the N11 Roughness Grade Number which is obtained using sand cast or hot roll manufacturing of heat sink, [26]. Finally, the minimal edge dimension of l = 40 µm is still big enough to be in the continuum mechanics having a Knudsen number value of 600. 3 CONCLUSIONS The infnite cooling capacity idea is very naive. The numerical experiment annulated the idea of course. Decreasing the fractal length l to zero the main conclusions are: • The area of interface surface converges to infnity. • Nu and h converged to zero setting the convective heat fow to zero (bearing in mind the extrapolation error). • The Q .is increasing monotonically to a specifc fnite value: heat diffusion. Fig. 9. Vorticity function contoursfor theKoch snowfakeformation procedure denotedbyfractal element lengthl In this manner the resulting heat transfer Eq. (11) could be written as ˜° Q .= lim (h ˜0)(A ˜˜)(°T = 1) ,(12) l˜0 as stated in Eq. (10). Since the Nusselt number represents the ratio between convection and diffusion, setting the Nu ˜0annulated heat convection leaving the diffusion the only heat transfer mechanism in the solid fuid interface, as it is thefact. Thefact is also, that each real surface has a roughness, might be in the fractal manner or not, and that the heat convection from solidsurfaceto fuidis achievedbyheatdiffusiononly at the fuid solid interface. The numerical experiment in this article confrms thisfact. The fractal geometry heat sink as an effective heat sink? No. This kind of fractal heat sink is non-competitive to the simple straight fn heat sink. This is clearly revealed by almost constant fractal fn temperature for aluminium – air confguration inTab. 2and4indicating the optimal fn should be slenderer generating larger temperature changes. 4 ACKNOWLEDGEMENTS Authors acknowledge the fnancial support from the Slovenian Research Agency (research core funding No. P2-0196). 5 REFERENCES [1] AIAA (1994). Editorial policy statement on numerical accuracy and experimental uncertainty, AIAA Journal, vol. 32, no. 1., p. 3-8, DOI:10.2514/3.48281. [2] Mocnik, U., Blagojevic, B., Muhic, S. (2020). Numerical analysis with experimental validation of single-phase fluid flow in a dimple pattern heat exchanger channel. Strojniški vestnik - Journal of Mechanical Engineering, vol. 66, no. 9, p. 544­553, DOI:10.5545/sv-jme.2020.6776. [3] Yadav, S., Pandey, K.M. (2018). A parametric thermal analysis of triangular fins for improved heat transfer in forced convection. Strojniški vestnik - Journal of Mechanical Engineering, vol. 64, no. 6, p. 401-411, DOI:10.5545/sv­jme.2017.5085. [4] Al-Kouz, W.G., Kiwan, S., Alkhalidi, A., Sari, M., Alshare, A. (2018). Numerical study of heat transfer enhancement for l = 1.0 l = 0.333 l = 0.111 l = 0.037 l = 0.012 l = 0.004 low-pressure flows in a square cavity with two fins attached to the hot wall using Al2O3-air nanofluid. Strojniški vestnik - Journal of Mechanical Engineering, vol. 64, no. 1, p. 26-36, DOI:10.5545/sv-jme.2017.4989. [5] Carnogurská, M., Príhoda, M., Lázár, M., Jasminská, N., Gallik, R., Kubik, M. (2016). Measuring selected parameters of polypropylene fibre heat exchangers. Strojniški vestnik - Journal of Mechanical Engineering, vol. 62, no. 6, p. 381-388, DOI:10.5545/sv-jme.2015.3202. [6] Yayla, S. (2013). Flow characteristic of staggered multiple slotted tubes in the passage of a fin tube heat exchanger. Strojniški vestnik - Journal of Mechanical Engineering, vol. 59, no. 7-8, p. 462-472, DOI:10.5545/sv-jme.2012.902. [7] Fong, R.W.L., McRae, G.A., Coleman, C.E., Nitheanandan, T., Sanderson, D.B. (2001). Correlation between the critical heat flux and the fractal surface roughness of zirconium alloy tubes. Journal of Enhanced Heat Transfer, vol. 8, no. 2, p. 137-146, DOI:10.1615/JEnhHeatTransf.v8.i2.60. [8] Chen, Y., Zhang, C., Shi, M., Peterson, G.P. (2009). Role of surface roughness characterized by fractal geometry on laminar flow in microchannels. Physical Review E, vol. 80, art. ID 026301, DOI:10.1103/PhysRevE.80.026301. [9] Liu, J., Xue, Y., Zhang, Q., Wang, H., Wang, S. (2022). Coupled thermo-hydro-mechanical modelling for geothermal doublet system with 3D fractal fracture. Applied Thermal Engineering, vol. 200, art. ID 117716, DOI:10.1016/j. applthermaleng.2021.117716. [10] Cai, W., Chen, W., Wang, F. (2017). Three-dimensional hausdorff derivative diffusion model for isotropic/anisotropic fractal porous media. Thermal Science, vol. 22, no. suppl. 1, p. s1-s6, DOI:10.2298/TSCI170630265C. [11] Gao, Y., Gao, F., Dong, G., Yan, W., Yang, X. (2019). The mechanical properties and fractal characteristics of the coal under temperature-gas-confining pressure. Thermal Science, vol. 23, no. suppl. 3, p. 789-798, DOI:10.2298/ TSCI180610094G. [12] Wang, Q.-L., Li, Z.-B., Kong, H.-Y., He, J.-H. (2015). Fractal analysis of polar bear hairs. Thermal Science, vol. 19, no. suppl. 1, p. 143-144, DOI:10.2298/TSCI15S1S43W. [13] Huang, Z., Hwang, Y., Aute, V., Radermacher, R. (2016). Review of Fractal Heat Exchangers. International Refrigeration and Air Conditioning Conference. Paper 1725. [14] Ikegami, S., Umetani, K., Hiraki, E., Sakai, S., Higashino, S. (2018). Feasibility study of fractal-fin heat sink for improving cooling performance of switching power converters. IEEE International Telecommunications Energy Conference, p. 1-8, DOI:10.1109/INTLEC.2018.8612377. [15] Xia, C., Fu, J., Lai, J., Yao, X., Chen, Z. (2015). Conjugate heat transfer in fractal tree-like channels network heat sink for high-speed motorized spindle cooling. Applied Thermal Engineering, vol. 90, p. 1032-1042, DOI:10.1016/j. applthermaleng.2015.07.024. [16] Ahmed, H.E., Salman, B.H., Kherbeet, A.Sh., Ahmed, M.I. (2018). Optimization of thermal design of heat sinks: A review. International Journal of Heat and Mass Transfer, vol. 118, p. 129-153, DOI:10.1016/j.ijheatmasstransfer.2017.10.099. [17] Xu, S., Wang, W., Fang, K., Wong, C.-N. (2015). Heat transfer performance of a fractal silicon microchannel heat sink subjected to pulsation flow. International Journal of Heat and Mass Transfer, vol. 81, p. 33-40, DOI:10.1016/j. ijheatmasstransfer.2014.10.002. [18] He, Z., Yan, Y., Zhao, T., Zhang, L., Zhang, Z. (2021). Multi-objective optimization and multi-factors analysis of the thermal/hydraulic performance of the bionic y-shaped fractal heat sink. Applied Thermal Engineering, vol. 195, art. ID 117157, DOI:10.1016/j.applthermaleng.2021.117157. [19] Zhang, S., Chen, X., Wu, Z., Zheng, Y. (2019). Numerical study on stagger koch fractal baffles micromixer. International Journal of Heat and Mass Transfer, vol. 133, p. 1065-1073, DOI:10.1016/j.ijheatmasstransfer.2019.01.009. [20] Ramšak, M., Škerget, L. (1999). Mixed boundary elements for laminar flows. International Journal for Numerical Methods in Fluids, vol. 31, no. 5, p. 861-877, DOI:10.1002/(SICI)1097­0363(19991115)31:5<861::AID-FLD900>3.0.CO;2-S. [21] Žagar, I., Škerget, L. (1995). The numerical simulation of non-linear separation columns by boundary-domain integral formulation. Computers & Chemical Engineering, vol. 19, Suppl. 1, p. 785-790, DOI:10.1016/0098-1354(95)87130-6. [22] Ramšak, M., Škerget, L. (2000). Iterative least squares methods for solving over-determined matrix for mixed boundary elements. ZAMM Berlin, vol. 80, no. 3, p. 657. [23] Ramšak, M. (2015). Conjugate heat transfer of backward-facing step flow: A benchmark problem revisited. International Journal of Heat and Mass Transfer, vol. 84, p. 791-799, DOI:10.1016/j.ijheatmasstransfer.2015.01.067. [24] Ramšak, M., Škerget, L. (2014). A highly efficient multidomain bem for multimillion subdomains. Engineering Analysis with Boundary Elements, vol. 43, p. 76-85, DOI:10.1016/j. enganabound.2014.03.009. [25] Ramšak, M. (2019). Multidomain bem for laminar flow in complex fractal geometry. Engineering Analysis with Boundary Elements, vol. 101, p. 310-317, DOI:10.1016/j. enganabound.2019.01.003. [26] DeGarmo, E.P., Black, J.T., Kohser, R.A. (1979). Materials and processes in manufacturing, Macmillan, New York. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 529-541 Received for review: 2022-01-17 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-07-07 DOI:10.5545/sv-jme.2022.15 Original Scientific Paper Accepted for publication: 2022-08-26 MTPA- and MSM-based Vibration Transfer of 6-D O F Manipulator for Anchor D rilling Youyu Liu1,3 – Liteng Ma1,3 – Siyang Yang1,3 – Liang Yuan2 – Bo Chen1,3 1 Key Laboratory of Advanced Perception and Intelligent Control of High-end Euipment, Ministry of Education, China 2 Technology Department, uhu Yongyu Automobile Industry Co., China 3 School of Mechanical Engineering, Anhui Polytechnic University, China An anchor drilling for a coal mine support system can liberate an operator from heavy work, but will cause serious vibration, which will be transmitted to the pedestal from the roof bolter along a manipulator. Based on the multi-level transfer path analysis (MTPA) and modal superposition method (MSM), a vibration transfer model for the subsystem composed of the joints of a manipulator with six degrees of freedom (DOF) was established. Moreover, its frequency response function matrix was also built. The 6-DOF excitation of the roof bolter was deduced. The exciting force on the roof bolter transmitted to the pedestal along the 6-DOF manipulator was analysed with a force Jacobian matrix, to identify the external loading on the pedestal. A case in engineering practice shows that the amplitude of each DOF of the pedestal from large to small is as follows: bending vibration (component 1), longitudinal vibration, torsional vibration, bending vibration (component 2), rotational vibration around z-axis, rotational vibration around y-axis. The pedestal is mainly in the form of bending vibration. The theory of vibration transfer along the 6-DOF manipulator for anchor drilling proposed in this article can provide a theoretical foundation for the development of vibration-damping techniques and the design of absorbers. Keywords: manipulator, multi-level transfer path analysis, modal superposition method, vibration transfer, force Jacobian matrix Highlights • Based on MTPA and MSM, a mathematical model of vibration transfer of a 6-DOF manipulator for anchor drilling is established. • The external loading of the response point of the manipulator pedestal is analysed by using the force Jacobian matrix. • The vibration responses on each DOF of the pedestal and the resonance frequency are obtained. • The case studied in an engineering practice shows that the pedestal is mainly in the form of bending vibration. 0 INTRODUCTION Roof bolters are key mechanical euipment for a coal mine supporting system. In the past, drilling was done manually. orking in an area with high concentrations of dust for a long time, workersphysical and mental health will be seriously threatened. At present, the development of a coal mine tunnel support tends to be automatic and intelligent. The manual labour of roof bolters has been gradually replaced by mechanical clamping [1]. An operator can control the roof bolter to drill automatically by human-computer interaction, which can liberate the operator from heavy work, and improve the stability and safety of the coal mine support. Due to the comprehensive excitation of different geotechnical parameters, axial thrust, torue and other factors, there are complex vibrations on drill strings during construction. The main forms include bending vibration, longitudinal vibration, and torsional vibration, which interact with each other to form a nonlinear coupled vibration [2]. The excitation vibration of each degree of freedom (DOF) of the roof bolter is transmitted to the pedestal along a manipulator, which causes the manipulator to vibrate violently, shorten its service life, and then affect the support effect. Transfer path analysis (TPA) is a tool to study vibration transfer [3] and [4]. There are several methods, such as operational transfer path analysis (OTPA) [5] and [6], global transmissibility direct transmissibility (GTDT) [7] and [8], inverse substructure TPA (ITPA) [9], multi-level transfer path analysis (MTPA) [10] and [11], and so on. Lee and Lee [5] proposed the OTPA method using an emerging deep neural network model, which can successfully predict the path contributions using only operational responses. Yoshida and Tanaka [6] attempted to calculate the vibration mode contribution by modifying OTPA, and then considered the relationship between the principal component and the vibration mode, as well as the associated the principal components with the vibration modes of a test structure. High contributing vibration modes to the response point have been found. It is easily disturbed by factors such as excitation coupling and noise employing this method when calculating the transfer matrix. ang [7] developed further the prediction capabilities of the GTDT method, which can predict a new response using measured variables of an original system, even though operational forces are unknown. Guasch [8] addressed some issues concerning the prediction capabilities of the GTDT method when blocking transfer paths in a mechanical system and outlined differences with the more standard force TPA. ang [9] developed the SDD method further by considering the mass effect of resilient links, which can identify decoupled transfer functions accurately, whilst eliminating the mass effect of resilient links. However, its manoeuvrability is poor for a serial system with many substructures. Gao [11] used MTPA to find the critical paths of seat jitter caused by dynamic unbalance excitation of the drive shaft. The key technology of this method is to identify the external excitation loading, which has good operability for series system. For the vibration problem of pedestal from manipulator caused by the excitation of a roof bolter, a response amplitude matrix in pedestal is established by the modal superposition method (MSM) in this article. According to the excitation of the roof bolter, the external loading of the response point of the manipulator pedestal is analysed using the force Jacobian matrix. The 6-DOF freuency response function of each subsystem of the manipulator is derived by MTPA, and then the freuency response matrix is constructed, which can solve the problem of Transfer path analysis with low accuracy and poor operability. It will provide a theoretical foundation for the development of vibration damping techniues and the design of absorbers. 1 MULTI-LEVEL TRANSFER PATH ANALYSIS To reduce the influence of non-important factors while analysing the vibration transfer of manipulator for anchor drilling, some simplifications are made as follows. 1) Each linkage of the manipulator is euivalent to a bar with uniform mass; 2) some transfer mechanisms, such as belt driving and harmonic decelerator in the manipulator, are euivalent to linear massless springs; 3) the modal parameters of the manipulator are linear, namely, the output caused by any combined input are eual to the combination of respective outputs; 4) it satisfies the assumption of time-invariance [12], namely, the dynamic properties of the system are not vary with time. According to the above simplifications, a vibration transfer model of the manipulator for anchor drilling is established, as shown in Fig. 1. This is a multi-input and multi-output system, in which the six joints of the manipulator are connected in series. The vibration source is the roof bolter that provides excitation; the vibration receiver is the pedestal. The manipulator can be divided into six subsystems by the rotating joint. The external loading of the J6-subsystem is the excitation from the roof bolter; that of other subsystems is the output of the previous subsystem. The output (response point) of the subsystem is the input (exciting point) of the next subsystem. Nevertheless, the output of the J1­subsystem acts on the pedestal. According to different excitations of spatial degrees of freedom, each subsystem has m inputs and n outputs. According to MTPA, the transfer function of manipulator for anchor drilling can be expressed by the product of the transfer functions of all subsystems [11]. H ˜HH °°°H .(1) JJJ 126 The vibration response of the excitation from the roof bolter transmitted to the pedestal is expressed as follows. S= HF. (2) 2 MODAL SUPERPOSITION METHOD 2.1 Response Amplitude Matrix The dynamic euation of the manipulator for anchor drilling is as follows. Fig. 1. Vibration Transfer model of the manipulator .. . n F  MS+CS+KS=F. (3) oa opa Sp ˜ .(12) a ˜1 2 °Ma Ka /Ma . Ma  ..1 According to the linear hypothesis [12], the Ka 2j a .KK .. aa  .ˆ displacement column vectors of each subsystem of the manipulator can be expressed as the linear addition of Their freuencyresponse function (FRF) is as E. each order of modal shapes. (13). S ˜°n qa . a .(4) a ˜1 Sp H Substituting E. (4) into E. (3), then, po ˜ ˜ Fo  oa n  ° n ° n ° n (13) pa M C K q .. q . F . . . . . ˜ (5) q . .  1 j 1 ˜ TT ˜ M ..nq .. ˜ .ˆ˜ C ..nq . ˜ .According to the linear superposition assumption b aab aa a °1a °1 [12], when F ˜°FF . FN .T, the response n 12 TT ˆ˜ K ..q ˜ .°˜ F .. (6) b aab a °1amplitude of each point of the system is as E. (14). a a a a a a ˜1 ˜1˜1   a a a a ° ...ˆ M  KM 2M  / K a a a ... .  a 2 Both ends of E. (5) are multiplied by . b T, then, a a K K  a a According to the orthogonality of the main modal S 1 HH H F 1 ° ° °. ° . . shape [13], E. (7) can be obtained from Es. (5) and N 11121 .... . ... ˜ ........ .... .... (6). S HH H F 2 ˜ ° N S ˜ . (14) 221222 . Mab , a ° . ˜ ˜˜.˜ ˜ ˜ TT ˆˆ. ˆˆ. b a S HH H FN ° . .ˆ . ˆ ˆ Cab °, . (7) a ˜ ˜ MCK ° N N N NN 12 b a ° Kab , a ˆˆ. ˆˆ. ˜ ˜ . T ab Substituting E. (7) into E. (6), then, Since the excitation of the roof bolter includes all DOF .. . TT Mq Cq ˜˜Kq °. F °. F.(8) in space, and each subsystem has m(6) inputs and n (6) aaaaaaa b outputs, the freuency response of each subsystem of The excitation loading and displacement response the manipulator is a 66th order matrix, as follows. in E. (8) are expressed in complex form as follows. b a 0, 2.2 Parameter Identification JJJ °.i °.i °. i ˆ jt jt ˆ. ° .... ˜˜ Ff qQ a . . H H ˜H e () 111216 ˆ  JJJ °.i °.i °. i H H ˜H e ˆ  a H 212226 ˜ .(15) ˆ J i °°.° Substituting E. () into E. (8), ˆ   . H °. H °J.˜H °J. ii i ˆ . J Tf . 616266 Qa ˜ (10) a . °.2M .j .C .K a aa The parameters of the J1 to J6 freuency response curves are identified in the freuency domain [14] by Substituting Es. (10) and () into E. (4), STcan be obtained as E. (11). Tjt the rational polynomial method [15]. Its mathematical ˆ. . f e aa model is a rational formula of freuency response .(11) T S ˜.an ˜12 °ˆM .jC K a ˆa . function, as follows. a It is assumed that the system has two points: o and p, substitute Es. (11) and () into E. (4), the response amplitude of the point p can be expressed as 5432 .x ..x ..x ..x ..x .. 123456 H J..(16) follows. ˜° mni 5432 ˆx .ˆx .ˆx .ˆx .ˆx .ˆ 123456 3 EXTERNAL EXCITATION LOADING E. (1) is expressed in matrix form as follows. FczD ° FcyD . ˆ 3.1 Excitation from Roof Bolter Fc ˜ 0 Fcz Fcy 0 . (20) .. .. 22 Force and moment on roof bolter: Fg, Fa, F, F, and c M. The direction of Fg, Fa and Fis along the shaft d of the roof bolter, and their vector expressions is as follows. The excitation from the roof bolter at the tip of the manipulator is as E. (21). .F g ˜.°Fg 00000 .. F z ˜ˆFz 00000 .. (17) F a ˜ˆFa 00000 .   The direction of Md is along the shaft of the roof bolter, and its vector expression is as follows. Md ˜°000 Md 00 .. (18) hile the lateral displacement of the roof bolter is greater than the distance of them, the roof bolter will collide with rock-soil. The collision force in the and y axes is as following [16], respectively. ˆ.ˆ°kvt .. °sgn vt .. vt .. ..˜ Fcz . 0 else . . . (1) . .ˆ.°kwt .. °sgn wt .. wt .. .  F ˜ .cy . ..0 else  .FD .FcyD . cz F˜°.Fa .Fg .Fz Fcz Fcy Md .(21) 22 ˆ 3.2 Excitation to Pedestal Force and torue on pedestal: Fg2, 1. Each joint of the manipulator can rotate independently. To accurately describe the mechanical properties of the excitation transmitted to the pedestal through the joints of the manipulator, a force Jacobian matrix of the manipulator is introduced [17]. The transfer relationship between the excitation and the joint generalized driving force is as follows [18]. ˜˜°J..q T F (22) The torues of each joint of the manipulator is as follows (E. (23)).  J J dcs a000 ac 44536   ˜ . Fx y JJ  . ˆ ˆ .ˆ ss 46  asc a3 000 1 x x 124 dccc 45656 ......... ......... ˜˜˜˜˜ ......... ......... .................. asc  3  F JJ  ss 46 000 2 1y 2y 4  dccc 45656 Fz 3 ˜° ˜ 6 where J ˜°d .cccc °ss °ssc .ac ac °ds sc c cs .; ci ˜cos(; si °) 1x 2ˆ23.45646.2356.°.22323423..45646°i )˜sin(i ; J ˜°d °cccc .ss ssc .ac .ac °ds sc c °cs .; s.°.) 1y 2.ˆ23.45646..2356..22323423..454623˜sin(2; 563 ˜d ccs °sc °ac °ac .ds ss c ˜cos(.°.); J ......; 23231z 223452352232342345 J ˜asc °dccc °ss asccc °ss csc ; 2x 3564.45646..2.ˆ3.45646..356. J ˜°asc °d .°ccc °ss ..as .°ccc °ss ..csc ˜ac °dcs °a ..scs °cc .; 2y 3564456462ˆ.345646356.; J 2z 36445234536 °˜° F ˜F sincos.; FF sincos.; F ˜F cos°. xyz ° 1z 2z ,(23)  M   ss 46 csc 2356 scc 456 cs 46 scc 456 cs sc 465 56 s 60 4 x 23  sccc 456  M  ss 46 scs csc 2356456 scc cc 46456 cc 46 ss 5 c 60 5 y 23  sccc 456 M  scs 2345 ss 45c 0 cc 23545ss 51 6 z 4 CASE IN ENGINEERING PRACTICE 4.1 Essential Parameters A 6-DOF manipulator for anchor drilling in a coal mine in Huainan, China, is taken as the research object. The size of the two-wing drill adopted is .32 mm; the length of drill string is 86 mm. The drilling object is sandstone, and its mechanical parameters [19] are as follows: . 2600 kg/m3; R 38 MPa; R 0.34 MPa; E 12 GPa; µ 0.25; F 6000 N; Md130 N; 0 mm; k 10 Nm–1; self-weight of the roof bolter is 40 kg; self-weight of the manipulator is 550 kg. Moreover, the parameters of the linkages of the manipulator are shown in Table 1 [20]. The three-dimensional model of the manipulator is shown in Fig. 2. Some finite element models of the manipulator are established, as shown in Fig. 3. The rotating joints of the manipulator are divided into some subsystems, and its exciting points and response points are determined. Based on MTPAand MSM, the computation flow chart of the vibration transfer of 6-DOF manipulator for anchor drilling is shown Fig. 4. ma Table 1. Parameters of the linkages of the manipulator Linkages Angle variable ai –1 [m] d i [m] Angle range [rad] 1 0 0.56 –3.14 to 3.14 2 0.90 0 –2.27 to 1.22 3 0.16 0 –1.40 to 3.05 4 0 1.01 –6.28 to 6.28 5 0 0 –2.09 to 2.09 6 0 0.2 –6.28 to 6.28 Fig. 3. Grid division of subsystems and exciting points, response points; a) J1, b) J2, c) J3, d) J4, e) J5, f) J6 4.2 Frequency Response Curves The excitation of the roof bolter is high-freuency vibration; the freuency range in practice is 0 Hz to 200 Hz [21]. Substituting Es. (1) and (21) into E. (13), the freuency responses of subsystems are analysed using ABAQUS [22]. The acceleration freuency response curves of J1 to J6 are shown in Figs. 5 to 10. According to Figs. 5 to 10, there are resonance peaks [23] in the freuency response curves of each subsystem of the manipulator for anchor drilling in the Fig. 4. Computation flow chart Fig. 5. Frequency response curves of J1 ; a) translational acceleration, and b) rotational acceleration Fig. 6. Frequency response curves of J2; a) translational acceleration, and b) rotational acceleration Fig. 7. Frequency responses curve of J3; a) translational acceleration, and b) rotational acceleration Fig. 8. Frequency response curves of J4; a) translational acceleration, and b) rotational acceleration Fig. 9. Frequency response curves of J5; a) translational acceleration, and b) rotational acceleration freuency range of 0 Hz to 200 Hz. The freuencies Table 2. Frequencies corresponding to peak values corresponding to the peak values of the 6-DOF Joints DOF Frequencies freuency response curves is shown in Table 2. RRR x yx y J1 191.3 191.3 192.2 197.3 197.4 197.2 4.3 FRF Matrixes of Subsystems J2 88.95 90.32 90.32 91.55 91.34 91.24 J3 42.19 42.19 91.77 86.2 82.9 88.1 The freuency response curves of each subsystem are J4 137.5 135.2 137.9 135.3 137.8 137.8 imported into MATLAB, which are fitted according to J5 178.8 187.4 178.8 187.4 178.7 178.8 E. (16) with the Curve Fitting Tool toolbox [24]. J6 91.27 91.34 91.54 81.66 73.24 95.54 The coefficients of each element in the freuency response matrixes are shown in Tables 3 to 8. Table 3. Coefficients of H (J1) mn Coefficients H 11(J1) H 12(J1) H 13(J1) H 14(J1) H 15(J1) H 16(J1) a1 1.053×10–5 –2.359×10–4 0 –0.04122 –2.474 –3.082 a2 18 7.91×10–4 5.53×10–4 9.311 569.6 708.6 a3 –14.27 8.967 –2.796×10–3 –99.29 –8613 –1.05×104 a4 –6.116 –9.886 2.15 –988.8 18.74 –328.7 a5 3.172 –0.2581 –3.32 –155.4 58.47 7.649 a6 1.583 2.063 1.303 –22.72 11.5 3.809 ß1 0 0 0 1 1 1 ß2 1 0 0 –398.7 –398.8 –398.8 ß3 –0.793 1 0 3.98×104 3.983×104 3.982×104 ß4 –0.3398 –1.105 1 1979 376.4 1047 ß5 0.1762 –2.802×10–2 –1.519 204.9 –200.2 –50.2 ß6 0.08792 0.2307 0.5861 30.14 –35.51 –10.31 Table 4. Coefficients of H (J2) mn Coefficients H 11(J2) H 12(J2) H 13(J2) H 14(J2) H 15(J2) H 16(J2) a1 14.37 0 0 0 0 0 a2 –5269 18.7 1345 0 0.01195 8.488×10–3 a3 7.419×105 –4285 –2.042×105 0 –1.265 –1.264 a4 –4.749×107 3.445×105 2.55×106 93.75 –30.33 28.8 a5 1.163×109 –1.361×107 4.565×108 –1.865×104 5454 1803 a6 4.135×107 3.364×108 –4.691×108 9.654×105 –2.66×104 –9660 ß1 1 0 1 0 0 0 ß2 –367.8 1 26.8 1 0 0 ß3 5.199×104 –221.8 –1.898×104 –406 1 1 ß4 –3.341×106 1.735×104 –2.815×105 6.166×104 –185.9 –185.8 ß5 8.21×107 –6.762×105 9.554×107 –4.15×106 9099 9093 ß6 3.09×106 1.654×107 –9.949×107 1.045×108 –4.014×104 –4.077×104 Table 5. Coefficients of H (J3) mn Coefficients H 11(J3) H 12(J3) H 13(J3) H 14(J3) H 15(J3) H 16(J3) a1 1.396×10–2 2.342 0.4121 1.18×10–3 8.639×10–7 0.000374 a2 –3.895 –276.7 –83.52 3.508 5.586 27.74 a3 416.3 1.218×104 5613 –720 –1106 –5552 a4 –1.863×104 –2.488×105 –1.427×105 3.802×104 5.885×104 3.039×105 a5 2.958×105 2.244×106 1.487×106 –2.085×105 –4.029×105 –2.49×106 a6 961 –1.232×106 –10.48 2.724×105 1.164×106 5.792×106 ß1 1 1 1 0 0 0 ß2 –183.2 –155.3 –242.3 1 1 1 ß3 1.311×104 9218 2.06×104 –194.5 –191.7 –193.6 ß4 –4.286×105 –2.496×105 –7.338×105 1.006×104 9927 1.028×104 ß5 5.332×106 2.666×106 1.026×107 –5.507×104 –6.78×104 –8.394×104 ß6 1.613×104 –1.849×106 –74.75 7.189×104 1.955×105 1.951×105 4.4 Torques of Joints symmetrically with the change of joint angle of the manipulator. hen .2is at the ultimate angle of –2.27 Substituting the data in Table 1 into E. (23), the rad and .4 is at that of –6.28 rad, t1 is only –220 values of t1 change with ., . and . are obtained by Nm as .1 and .3 change. hile .3 . (–1.24 1.13) i MATLAB, as shown in Fig. 11. t1 is distributed rad, t1 shows a trend of decay, when .3 . (1.13 Table 6. Coefficients of H (J4) mn Coefficients H 11(J4) H 12(J4) H 13(J4) H 14(J4) H 15(J4) H 16(J4) a1 0 0 0 0 0.9853 –2.179 a2 1.452 10.87 3.019 1.283 –174.3 1054 a3 –542.2 –3655 –729.4 15.19 6824 –1.6×105 a4 8.655×104 3.114×105 4.697×104 –9.851×104 8.145×104 8.051×106 a5 –7.894×106 –6.236×105 –2.925×105 1.105×107 –6.495×105 –1.18×107 a6 3.419×108 3.929×104 3.765×105 –3.454×104 6.07×105 5.426×104 ß1 0 0 0 1 1 1 ß2 1 1 1 –624.9 –256.5 –627.4 ß3 –538.6 –341.5 –277.6 1.458×105 1.548×104 1.606×105 ß4 1.113×105 2.958×104 2.022×104 –1.506×107 1.606×105 –1.892×107 ß5 –1.044×107 –5.924×104 –1.247×105 5.81×108 –1.351×106 8.4×108 ß6 3.746×108 3828 1.529×105 –1.813×106 1.274×106 –3.868×106 Table 7. Coefficients of H (J5) mn Coefficients H 11(J5) H 12(J5) H 13(J5) H 14(J5) H 15(J5) H 16(J5) a1 0 0 0 2.245 9.23×10–3 0.0141 a2 1.223×10–2 4.005 12.8 –661.9 –3.182 –5.048 a3 –4.295 –1464 –3564 3.512×104 312.4 481 a4 388.1 1.366×105 6.619×104 2.645×106 –6312 –3513 a5 –31.56 –4.933×105 3.108×107 –8.582×106 5.686×104 –3.849×104 a6 46.17 3.554×105 –2.391×107 2.1×106 436 1.514×105 ß1 0 0 0 1 0 0 ß2 1 1 1 –23.44 1 1 ß3 –346.4 –366.1 696.9 –8.757×104 –364.4 –369.2 ß4 2.834×104 3.423×104 –3.441×105 1.099×107 3.424×104 3.642×104 ß5 3.039×105 –1.237×105 3.376×107 –4.608×107 –1.834×105 –4.209×105 ß6 6662 8.921×104 –2.461×107 4.1×107 –3.878×104 1.139×106 Table 8. Coefficients of H (J6) mn Coefficients H 11(J6) H 12(J6) H 13(J6) H 14(J6) H 15(J6) H 16(J6) a1 2.102×10–3 –3.807 0.06361 –7.357×10–6 –1.485×10–5 4.141 a2 –0.1981 1476 –14.32 3.304×10–3 18.62 –834.1 a3 –153.7 –1.796×105 1087 13.72 –3857 4.244×104 a4 2.772×104 7.604×106 –3.055×104 –2855 2.287×105 7347 a5 –1.185×106 –5.151×107 2.653×105 1.522×105 –2.966×106 –450.4 a6 1.48×107 4.774×107 –4.187×105 –5.962×105 8.696×106 –211.3 ß1 010001 ß2 1 –413.3 1 0 1 –195.1 ß3 –354.4 6.999×104 –162.5 1 –200.1 9678 ß4 4.741×104 –5.504×106 4823 –202.7 1.16×104 –1.156×104 ß5 –2.837×106 1.66×108 1.585×105 1.072×104 –1.497×105 5793 ß6 6.444×107 –1.577×108 –3.222×105 –4.198×104 4.382×105 1667 2.30) rad, t1 shows a steady trend; and .3 . (2.30 distributed obviously; and t1 max 117 Nm. The 3.04) rad, t1 showed a slight upward trend. hen .6 maximum of t1 applied to the pedestal is 117 Nm, is at the ultimate angle of –6.28 rad, t1 is distributed which is transmitted along the 6-DOF manipulator symmetrically with the change of .5; and t1 max 6018 with any position and posture in space. Nm. As . and f change, t1 is symmetrically a) t1 change with .1, .2, b) t1 change with .3 .4, c) t1 change with .5, .6, d) t1 change with ., f 4.5 Vibration Response of the Pedestal Substituting the freuency response matrix of each subsystem into E. (1) and substituting Es. (1) and (24) into E. (2), the vibration responses on each DOF of the pedestal are shown in Fig. 12, in which the positive and negative values of the vibration response only represent the direction. As shown in Fig. 12a, the response of the longitudinal vibration (S ) of the pedestal reaches the x maximum, being 1.6510–2m, when the freuency is 45 Hz. hile the freuencies are 0 Hz and 180 Hz, the responses of the two components (Sy and S ) of the bending vibration of the pedestal reach the maximum, being 2.1210–2m and 8.0610–3m, respectively. At this time, the freuencies corresponding to the peaks are integer multiples of each other, so the phenomenon of resonance will happen. As shown in Fig. 12b, the response of the torsional vibration (S ) of the pedestal reaches the maximum, R x being 10–3 m, when the freuency is 10 Hz. The two components (Syand S ) of rotational vibration RR RRR around y and axes sharply increase in the freuency interval 166-181 Hz and 151-180 Hz, and reach the maximum, being 2.0810–4 m and 2.4610–4 m respectively, at the freuencies of 180 Hz and 181 Hz. At this time, the freuencies corresponding to the peaks are very similar, which is easy to induce a resonance. The amplitude of each DOF of the pedestal is from large to small is as follows: bending vibration (component 1) Sy, longitudinal vibration Sx , torsional vibration S , bending vibration (component R x 2) S , rotational vibration around -axis SR , rotational vibration around y-axis Sy R . 5 CONCLUSIONS (1) Based on MTPA and MSM, a mathematical model of vibration transfer of 6-DOF manipulator for anchor drilling is established. The freuency response matrix of subsystems is derived under multi-DOF excitations. hen the external excitation loading is determined, the freuency response function of each DOF at the response point can be calculated by the mathematical model, which is universal for series systems. (2) ithin the freuency range (0 Hz to 200 Hz) of the excitation of the roof bolter, the corresponding freuency ranges of the peak values in the 6-DOF vibration direction of each subsystem are 11.3, 17.4 Hz, 88.5, 1.55 Hz, 42.1, 1.77 Hz, 135.2, 137.8 Hz, 178.8, 187.4 Hz and 73.24, 5.54 Hz, respectively. In the above freuency ranges, the vibration response will be the largest, and the resonance can be avoided by changing the excitation freuency of the roof bolter. (3) A case in engineering practice shows that the amplitude of each DOF of the pedestal is from large to small is as follows: bending vibration (component 1) S (2.1210–2 m) at 0 Hz, y longitudinal vibration S (1.6510–2 m) at 45 x Hz, torsional vibration S (10–3m) at 10 Hz, R x bending vibration (component 2) S (8.0610–3 m) at 180 Hz, rotational vibration around -axis S (2.4610–4m) at 180 Hz, rotational vibration R around y-axis S (2.0810–4 m) at 180 Hz. yR Obviously, the pedestal is mainly in the form of bending vibration. (4) The case also shows that a resonance will occur among the two components of the bending vibration at the freuencies of 0 Hz and 180 Hz; a resonance among the two components of rotational vibration around the y and axes is highly likely to occur at the freuencies of 180 Hz and 181 Hz. The theory of vibration Transfer along the 6-DOF manipulator for anchor drilling proposed in this article can provide a theoretical foundation for the development of vibration damping techniues and the design of absorbers. 6 NOMENCLATURES a the ath order modal shape T the b th order modal shape b .oa modal shape of the oth DOF of the ath modal vector .pa modal shape of the pth DOF of the ath modal vector q modal participation factors a o exciting point p response point n modal order M modal mass coefficient a Ca modal damping coefficient K modal stiffness coefficient a M mass matrix C damping matrix K stiffness matrix S vibration displacement S .vibration speed .. S vibration acceleration J(q )T force Jacobian matrix J joint J elements in Jacobian matrix d distance between two adjacent linkages along the i common axis m@ a–1 common perpendicular length between joint i –1 i and joint i m@ moving DOF of -axis x moving DOF of x -axis y moving DOF of y-axis rotational DOF of -axis R xR rotational DOF of x -axis yR rotational DOF of y-axis Sy bending vibration (component 1) m@ S longitudinal vibration m@ x S bending vibration (component 2) m@ S rotational vibration around x -axis (torsional R x vibration) m@ S yR rotational vibration around y-axis m@ S rotational vibration around -axis m@ R F external loading N@ Fy force projected on the y-axis N@ F loading at o point N@ o F Ji excitation force transmitted to joints N@ F force projected on the -axis N@ F force projected on the x -axis N@ x Fg2 gravity of manipulator N@ Fg gravity N@ Fa axial thrust N@ Fdisturbing force N@ Fimpact force N@ c Md torue M torue projected on the x -axis Nm@ x My torue projected on the y-axis Nm@ My torue projected on the -axis Nm@ v bending vibration displacement (component 2) m@ w bending vibration displacement (component 1) m@ µ Poissons ratio D outer diameter of drill string m@ . material density kgm-3@ R compressive strength MPa@ R tensile strength MPa@ m E elastic modulus GPa@ k elastic coefficient of impact force Nm-1@ . angle between force and -axis rad@ f angle between force and x -axis rad@ .joint angle rad@ i clearance between drill string and rock-soil mm@ .a damping ratio 2 torue matrix 1 torue from the J1-axis ai molecular coefficients, i 1 to 6 ßi denominator coefficients, i 1 to 6 7 ACKNOWLEDGEMENTS This work was supported in part by Special Fund for Collaborative Innovation of Anhui Polytechnic University & Jiujiang District under Grant No. 2021CYTB3, and by and by Natural Science Research Project of Higher Education of Anhui Province of China under Grant No. KJ2020A0357. 8 REFERENCES [1] Kang, Y.S., Liu, Q.S., Xi, H.L., Gong, G.Q. (2018). Improved compound support system for coal mine tunnels in densely faulted zones: A case study of China’s Huainan coal field. Engineering Geology, vol. 240, p. 10-20, DOI:10.1016/j. enggeo.2018.04.006. [2] Liu, C.C., Zheng, X.G., Arif, A., Xu, M.B. (2020). Measurement and analysis of penetration rate and vibration on a roof bolter for identification rock interface of roadway roof. Energy Sources, Part A: Recovery, Utilization, and Environmental Effects, vol. 42, no. 22, p. 2751-2763, DOI:10.1080/155670 36.2019.1618987. [3] Shu, J.C., He E.M. (2020). Review on vibration transfer path analysis methods in frequency domain. Mechanical Science and Technology for Aerospace Engineering, vol. 39, no. 11, p. 1647-1655, DOI:10.13433/j.cnki.1003-8728.20200270. [4] van der Seijs, M. V., de Klerk, D., Rixen, D.J. (2016). General framework for transfer path analysis: History, theory and classification of techniques. Mechanical Systems and Signal Processing, vol. 68-69, p. 217-244, DOI:10.1016/j. ymssp.2015.08.004. [5] Lee, D.H., Lee, J.W. (2020). Operational transfer path analysis based on deep neural network: Numerical validation. Journal of Mechanical Science and Technology, vol. 34, p. 1023­1033, DOI:10.1007/s12206-020-0205-5. [6] Yoshida, J.J., Tanaka, K.K. (2016). Contribution analysis of vibration mode utilizing operational TPA. Mechanical Engineering Journal, vol. 3, no. 1, p. 574-589, DOI:10.1299/ mej.15-00589. [7] Wang, Z.W., Zhu, P., Zhao, J.X. (2017). Response prediction techniques and case studies of a path blocking system based on global transmissibility direct transmissibility method. Journal of Sound and Vibration, vol. 388, p. 363-388, DOI:10.1016/j.jsv.2016.10.020. [8] Guasch, O. (2009). Direct transfer functions and path blocking in a discrete mechanical system. Journal of Sound and Vibration, vol. 321, no. 3-5, p. 854-874, DOI:10.1016/j. jsv.2008.10.006. [9] Wang, Z.W., Peng, Z.K., Liu C., Shi, X. (2019). Virtual decoupling of mechanical systems considering the mass effect of resilient links: Theoretical and numerical studies. Mechanical Systems and Signal Processing, vol. 123, p. 443-454, DOI:10.1016/j. ymssp.2019.01.028. [10] Sakhaei, B., Durali, M. (2014). Vibration transfer path analysis and path ranking for NVH optimization of a vehicle interior. Shock and Vibration, vol. 2014, art. ID 697450, DOI:10.1155/2014/697450. [11] Gao F., Li Y.O., Zhang J. (2015). The control strategy rseearch and engineering application of seat vibration at high speed of fr vehicle. Proceedings of the 2015 Annual Meeting of the Chinese Society of Automotive Engineering, p.192-195. [12] Shah, V-.N., Bohm, G.J., Nahavandi, A.N. (1979). Modal superposition method for computationally economical nonlinear structural analysis. Journal of Pressure Vessel Technology, vol. 101, no. 2, p. 134-141, DOI:10.1115/1.3454612. [13] Fu, Z.F., He, J.M. (2001). Modal Analysis. Elsevier, Woburn. [14] Karaagaçli T., Özgüven H.N. (2020). A frequency domain nonparametric identification method for nonlinear structures: Describing surface method. Mechanical Systems and Signal Processing, vol. 144, art. ID. 106872, DOI:10.1016/j. ymssp.2020.106872. [15] Dong, X.G., Guo, X.Y., Wang, Y.X. (2020). Local polynomial method for frequency response function identification. Systems Science & Control Engineering, vol. 8, no. 1, p. 534­540, DOI:10.1080/21642583.2020.1833784. [16] Piovan, M.T., Sampaio, R. (2006). Non linear model for coupled vibrations of drill-strings. Mecanica Computacional, p. 1751-1766. [17] Turkovic, K., Car, M., Orsag, M.(2019). End-effector force estimation method for an unmanned aerial manipulator. 2019 Workshop on Research, Education and Development of Unmanned Aerial Systems, p. 96-99, DOI:10.1109/ REDUAS47371.2019.8999670. [18] Ding, W.H., Deng, H., Zhang, Y., Ren, Y.Q. (2012). Optimum design of the jaw clamping mechanism of forging manipulators based on force transmissibility. Applied Mechanics and Materials, vol. 157-158, p. 737-742, DOI:10.4028/www. scientific.net/AMM.157-158.737. [19] Liu, S.W., Fu, M.X., Zhang, H. (2016). Vibration mechanism and characteristics analysis of drill rod when drilling roof bolt hole. Journal of China University of Mining & Technology, vol. 45, no. 5, p. 893-900, DOI:10.13247/j.cnki.jcumt.000562. [20] EFORT (2019). Robot Products, from https://www.efort.com. cn/product/prodetail/44.html, accessed on 2019-05-09. [21] Mendil C., Kidouche M., Doghmane M.Z., Benammar,, S., Tee, K.F. (2021). Rock-bit interaction effects on high-frequency stick-slip vibration severity in rotary drilling systems. Multidiscipline Modeling in Materials and Structures, vol. 17, no. 5, p. 1007-1023, DOI:10.1108/MMMS-10-2020-0256. [22] Weeger O., Wever U., Simeon B. (2014). Nonlinear frequency response analysis of structural vibrations. Computational Mechanics, vol. 54, p. 1477-1495, DOI:10.1007/s00466-014­1070-9. [23] Li, C.S., Luo, S.H., Ma, W.H. (2013). Vibration isolation performance analysis of a HXN3 diesel locomotive cab based on frequency response functions. Journal of Vibration and Shock, vol. 32, no. 19, p. 210-215, DOI:10.13465/j.cnki. jvs.2013.19.025. [24] Grauer, J.A. (2015). An interactive MATLAB program for fitting transfer functions to frequency responses. AIAA Scitech 2021 Forum, p, 1-14, DOI:10.2514/6.2021-1426. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 542-551 Received for review: 2022-04-05 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-08-02 DOI:10.5545/sv-jme.2022.147 Original Scientific Paper Accepted for publication: 2022-08-10 Functional Surface L aye r Strengthening and We ar R esistance Increasing of a L ow Carbon Steel by E lectrolyt ic-Plasma Processing Kuat Kombayev1 – Murat Muzdybayev1 – Alfiya Muzdybayeva1 – Dinara Myrzabekova1 – ojciech ieleba2 – Tadeusz Leniewski2, 1 East Kazakhstan State Technical University, Technological Machines and Transport Department, Republic of Kazakhstan 2 rocaw University of Science and Technology, Faculty of Mechanical Engineering, Poland The modified technology for strengthening the surface layer of low-carbon steel for machine components by electroplasma processing is proposed. This technology is to be used as an alternative carburizing method with subsequent hardening. The developed technology of strengthening by the proposed method is based on the author’s invention. The parameters of this process are given, resulting in a thickness of the reinforced surface layer of 1000 µm to 1700 µm. An increase in microhardness of 1.5 to 2 times (compared to the initial state) was observed. With hardening, chemical modification of the material’s surface layer occurs. The microstructure of the treated surface of the steel samples is characterised by a dark modified surface layer: a fine needle-like structure of martensitic origin under the dark layer transforms into an initial perlite-ferritic structure. The advantage of strengthening based on the electrolytic plasma processing consists of low energy consumption at high hardening speeds and the possibility of local processing of surface areas, especially large parts of complex shapes. In addition, the proposed surface treatment method using electrolytic plasma processing (EPP) not only achieves a smooth surface but also improves the service qualities of the components, specifically wear resistance. Keywords: low-carbon compound steel, strengthening, electrolytic plasma processing Highlights • The author’s technology of strengthening the surface layer of low-carbon compound steel for details of machines by electrolytic plasma processing is developed. • This technology is proposed as an alternative method of carburizing with subsequent hardening; the total thickness of the strengthened surface layer is 1000 µm to 1700 µm, and the microhardness of the strengthened area of a martensitic structure is 6500 MPa to 7200 MPa. • The microstructure of the processed surface of steel samples is characterized by a dark modified surface layer. A fine needle-like structure of martensitic origin under the dark layer changes into the initial pearlite-ferritic structure. • It is stated that in the structure of 18CrNi3Mo-Sh steel sample after electrolytic plasma processing, particles of residual cementite Fe3C, a-phase based on Fe, Cr0.6 Fe1.4 - phase and Fe2.7 Mo0.8 Ni0.1 - phase appear, which confirms availability of martensite hardening. • The advantage of strengthening based on the electrolytic plasma processing consists of low energy consumption at high hardening speeds and the possibility of local processing of surface areas, especially large parts of complex shapes. 0 INTRODUCTION Mechanical engineering development is largely related to solving the tasks of increasing machines¶ reliability. To increase the operational reliability of machines, it is necessary to increase the endurance and durability of their parts [1] and [2] and reduce the risk of their failures [3] and [4]. This is achieved by various methods, including strengthening the operating surfaces of component parts and creating new technological processes of highly wear-resistant materials manufacturing [5] to [7]. Toughening the reuirements concerning the microstructure and properties of the surface layers of machines¶ component parts reuires new methods for their modification and strengthening [8]. Techniues involving the influence of concentrated energy flows on the surface of component parts appeared to be widely used [9]. The most promising technology for surface strengthening of component parts is electrolytic plasma processing (EPP) [10]. EPP is an electrochemical heat treatment process and mass transfer (coating) on a component part (product) surface using a plasma jet. A plasma jet is a partially or fully ionized gas in combination with electric discharged phenomena at the boundary of the main electrode-water electrolyte solution at high capacities up to 1000 V[11]. ith a gradual increase in the applied constant voltage, salt electrolysis occurs. According to Ohms law, the current increases (section from 0 to A, Fig. 1). This area is characterized by scaling up the current with voltage increasing. Doing this, the electrolyte temperature also increases. This is a conseuence of the current passing through the electrolyte. hen a certain voltage value is reached (from 100 V to 176 V), the electrolyte starts boiling on the surface of the cathode, and an active release of bubbles near the surface occurs (bubble boiling). At bubble boiling, the temperature of the component part is close to the steam point. hen bubble boiling occurs around the active electrode, large current pulsations can be observed. Their amplitude is significantly reduced when the component part is heated above 470 C. Fig. 1. Volt-ampere characteristic of electrolyte plasma processing Over a further increase in rectified voltage, film boiling appears (point A, Fig. 1), which is characterized by the disappearance of bubble boiling and a sharp drop in current because the resulting gas-steam jacket has a higher electrical resistance than a liuid electrolyte (section from Ato B, Fig. 1). Since the gas-steam jacket is less electrically conductive, the main voltage drop occurs precisely where more heat is generated. The low-temperature plasma of a specific blue glow around the component part is formed due to a gas-steam jacket forming an electric current passing through it. The brighter the blue colour of plasma burning, the more ions it contains, including ion modifiers. An abnormal discharge is observed when increasing the voltage [12]. The principle of the technological process of strengthening based on the method of EPP is as follows. The cathode is made of 18CrNi3Mo-Sh steel (which corresponds to the standard 18HN3MA- TU 3-850-80) for manufacturing roller cones. It is immersed in the electrolyte (10 water solution Na2CO3) to 4 mm to 6 mm depth. The anode is stainless 10 CrNiTi 18-10 (corresponding to the standard 1218H10T GOST 54-75 Steel). It is in the form of a disk with a diameter of 50 mm and a thickness of 2 mm; 4 mm holes are drilled in the disk. Plasma arises between the cathode and liuid electrolyte. Negative ions ease excess electrons when passing through the holes of the anode of stainless 10 CrNiTi 18-10 steel. The cations are entrained by the hydrodynamic flow of electrolyte and recombine at the cathode (the surface of the sample material of the drill bit). The conversion of electrical energy into heat occurs mainly in the plasma layer on the heated cathode surface. hen being processed, cathode surfaces are cemented with carbon ions and hardened in an electrolyte in a short time. The modification of the surfaces of machines¶ component parts with the EPP method contributes to the compositional restructuring of a surface layer that increases strengthening properties and wear resistance due to physical input (ions of high-temperature plasma, electrical discharge). The EPP method is currently used when modifying the surfaces of component parts of bent shafts, iron cylinders of diesel engines, circular saws [13], and others. The main advantages of the EPP method are the following: a complex profile strengthening, internal surfaces and cavities; no need for special preparation of surfaces before strengthening; ecological safety (use of special treatment facilities is not reuired). Analysis of existing technologies for low-carbon compound steel products hardening suggested that electrolytic plasma strengthening could be the most appropriate technology for the thermal hardening of component parts. The aim of this research is to identify experimentally the values of parameters for the technological process of electrolyte-plasma hardening of 18CrNi3Mo-Sh steel. The article presents the research methodology description, the obtained results, their discussion, conclusions, and acknowledgements. 1 METHODS AND EXPERIMENTAL APPROACH A multiple-factor experiment has been fully implemented to select optimal electrolyte plasma processing conditions. Experiments of this type can localize the area necessary and present all possible independent variables¶ values. Heating a sample of a component part in depth is associated with an inhomogeneous temperature distribution in the volume of the material of a component part [14]. A drill bit (Fig. 2) provided by JSC Vostokmashzavod (VKM JSC) was selected to study the effectiveness of various methods of strengthening. Contact durability, abrasive and shock-abrasive wear resistance of the drill bits component parts are provided by gas carburizing followed by hardening. Deformation, low carbon steel cracking, and high labour and energy intensity are the disadvantages of that heat processing. In order to eliminate the disadvantages of the above treatment method, for EPP, the following samples from drill bits component parts were used: 18CrNi3Mo-Sh (0.16 to 0.18 ; 3.3 Ni; 0. Cr; 0.51 Mo; 0.44 Mn; 0.34 Si; 0.05 Al; 0.008 S; 0.012 P; 0.015 N; 0.01 O; 0.01 H) NSS 4543-71 [13]. Low-carbon, compound, heat-resistant steel of electroslag remelting can be used at –70 to 450 . Carbon, chromium, molybdenum, manganese, and silicon are alloying constituents, increasing strengthening properties. Nickel provides strength and good toughness; molybdenum provides heat resistance to steel. Samples of 10 mm 10 mm 20 mm were cut out of roller cones with a diamond blade with a thickness of 1 mm, which was immersed in cooling fluid. The samples were made from roller cones in the initial condition and compared after thermal processing on the JSC VKM base. Two series of samples (3 repetitions in each) hardened according to the method given by JSC VKM and processed using the EPP method were examined. Fig. 2. Roller cone At low cutting speeds (n 350 rpm) and low load (m 250 g), the sample does not experience significant deformations and thermal effects. The samples were first polished with chromium dioxide cement and then etched with 5 alcoholic nitric acid solution to do a metallographic microanalysis. Application investigations and mechanical testing were carried out in the Regional University Engineering Laboratory Irgetas and the research and manufacturing complex Machine industry of D. Serikbayev East Kazakhstan State Technical University [15]. Metallographic analysis was carried out using an Axioscop-2MAT microscope euipped with a Sony digital camera. Qualitative and uantitative phase analysis of the steel samples¶ structure was carried out using an Pert PRO -ray diffractometer of the PAN analytical firm, applying Cu-Ka radiation. The electron images were made with the following euipment: a raster electron microscope JSM-630LV JEOL Ltd. (Japan) with an energy dispersive microanalysis system INCA Energy Penta FET 3 OFORD Instruments Analytical Limited. Microhardness measurements were carried out using a PMT-3 device, euipped with a diamond pyramid, with a load on the 2H indenter by 450-76 GOST. The wear resistance of the samples was estimated by the loss of mass per unit time when the sample was being tested by attrition on the abrasive disc when sliding rubbing without lubricant [16]. To measure the mass of the samples, an electronic weighing unit of the VL­120 model with an accuracy of 0.1 mg was used. 2 RESULTS The influence of technological and electrical parameter changes during electrolyte (10 water solution of Na2CO3 calcined soda) on the strengthening uality was studied. It was experimentally proved that the considered factors of EPP modes affect the uality of the steel surface, which is strengthened. It should be noted that the main technological parameters, such as the thickness of strengthening, microhardness, and wear resistance of steel depend on the heating time, hardening time, the number of strengthening cycles by the EPP method and the electric current voltage [14]. Fig. 3 shows the initial microstructure of the material of the sample of a strengthened component part of 18CrNi3Mo-Sh steel. It is seen that a coarse-grained pearlite-ferritic microstructure characterizes the material of the sample. Fig. 3. The initial structure of the sample material of a component part of 18CrNi3Mo-Sh steel In the course of the study, a series of searching experiments were carried out, involving EPP testing for sample processing to identify the logical modes of the process. The ranges of the numerical values of the parameters studied are presented in Table 1. Table 1. Technological parameters for hardening of 18CrNi3Mo-Sh steel by the EPP method Terms in the program of the experiment Physical factors Factor levels min max The heating time of a work X1 piece from the plasma 1 15 temperature [s] X2 The hardening time in the flow of electrolyte [s] 1 10 X3 DC voltage [V] 180 260 X4 Number of cycles [-] 20 40 It was decided to switch to a three-factor experiment model and settle the number of processing cycles at 30 to conduct a comparative analysis of technological modes¶ influence on the material of the component part strengthening. To systematize the assessment of the EPP modes¶ effect on the uality indicators of steel of drill bits strengthening, depending on such main parameters as heating time T heati ng from ionized plasma and hardening time T harde ni ng in the electrolyte flow, taking into account the direct current (DC) voltage U between the electrodes, we converted micrographs of material microstructure into a tabular form (Table 2). In one case, the images located along the edges lack the modes for optimal hardening of the material of a sample (left). In another case, the samplematrix burns out (right), transforming the hardening process in EPP into the process of ion transfer of metal from the anode to the strengthened cathode. Table 2. EPP modes while steel drill bit strengthening Modes EPP processing Parameters ABC Heating time theati ng [s] 2 415 Hardening time tharde ni ng [s]2 4 15 DC voltage U [V] 190 200 250 Fig. 4 shows the microsections of electrolyte-plasma processed steel to the depth of the hardened layer (cross-section). A comparison of microsections (cross-sectional cuts to the deep from left to right), presented in Fig. 4, makes it possible to note the features of the microstructure of the material obtained at different modes of EPPprocessing: when processing at mode A, the influence of high-temperature plasma is observed (Fig. 4a, the dark area from left to right). The layer of coarse needle-like martensite smoothly transfers into the initial structure as it deepens into the material. hen processing at mode B, there is also an area of high-temperature plasma exposure. Its depth is comparatively larger than at mode A and is approximately 200 m deep (Fig. 4b, dark area from left to right). The availability of fine needle-like martensite with a smooth transition to the initial structure of the sample material is noted. a) b) c) Fig. 4. Micro sections of the cross-section (to the depth of material) of electrolytic plasma processed steel After 15 s of EPP processing (at mode C), a deposited layer of anode material (steel 10 CrNiTi 18-10) on the microsection occurs. The sample substrates border and martensitic transformation are clearly visible (Fig. 4c, to the depth from left to right). Thus, images of microsections reveal: mode Alacks parameter values of EPPprocessing for the optimal hardening of the sample material (Fig. 4a), mode B burns out the material of a sample while EPP processing (Fig. 4b), transforming the hardening process into the process of ion transfer of metal from the anode to the strengthened cathode. It is stated that electrolyte-plasma processing reuires low power (within 3 kh) when rapid motion occurs. The obtained data are used in [17]. In actual practice, a thermal exchange is a complex process. To facilitate its study and simplify the correlations obtained, the concept of elementary types of thermal exchange is introduced: experimentally determined thermocouples made by a natusral thermal junction. Two thermocouples are installed in two layers of the sample workpiece at a depth of h1 and h2 from the heated surface (Fig. 5). Thermocouples are installed to control the temperature mode of the process to prevent the materials melting from being strengthened. Fig. 5. The layout of thermal junctions in the surface layer of the sample when measuring the temperature in the process of electrolytic plasma heating A 10 mm 10 mm 20 mm steel sample acts as a cathode. The temperature is measured by a thermocouple mounted in the cathode body. The thermocouple is a thermoelectric converter of the type, conditional designation () according to GOST 3044-84. The measurement range is from minus 200 to plus 1000 . Admission class 1. The limits of permissible deviations are 0,004. when a measurement is in the range above plus 375 to plus 1350 . The sensitivity of the thermocouple is 40 V/ . The range of measured temperatures corresponds to the temperature of phase transformations from plus 840 to plus 860 . Heating temperature is the main parameter of phase transformations for 18CrNi3Mo-Sh steel and is 860 C [18]. At EPP, the bulk temperature increases to the hardening temperature of 860 C, and the temperature of the ionized plasma overheats the surface. Its temperature exceeds 1100 C. hen ionized plasma is excited (plasma temperature ranging from 6000 K to 30000 K), a gas-vapor layer appears on the sample surface due to electrolyte dissociation [19]. The gas-vapor layer prevents the ingress of electrolyte on the overheated surface. It reduces the cooling rate, eliminating the formation of thermal (hardening) cracks. As a result, the operating durability of steel increases. Thus, the main factors determining the uality of steel strengthening at EPP have been experimentally stated. They are heating time, hardening time, and electric current voltage. A mathematical model is developed to describe the change in the key parameter of the technological process of steel strengthening at the EPP method: the heating temperature T . The following regression euation expresses the logarithmic dependence of the temperature T on the main factors: ln() T Ca ln( ) .Db ln() ˜°°t °°U heating .°°ln( tcooling ), (1) Ed where . is steel heating time s, theati ng heating time s, tc ooling cooling time in the electrolyte flow s, U voltage V, and C 4.5, D 4.8, E 18, a 2, b 1, c 1 the coefficients for E. (1) were found by using logarithm in the Deductor Studio Academic software. Then E. (1), which describes the dependence of the heating temperature on the heating time, cooling time and voltage, was transformed to E. (2): T ˜°t 2.°.°.U 18t 45.heating 48cooling ,(2) where . is steel heating time s, theati ng heating time s, tc ooling cooling time in the electrolyte flow s, and U voltage V [20]. Experimentally determined optimal modes of steel hardening by the EPP method (theati ng 4 s, tc ooling 4 s, U 200 V) have a good correlation with the established dependence, see E. (2). A raster elemental analysis of the processed surface (Fig. 6) revealed chemical modification of the surface layer of the metal happening along with hardening during the electrolytic plasma heating of the sample [21]. Fig. 6. Raster elemental analysis of the 18CrNi3Mo-Sh steel surface after EPP strengthening Table 3. Qualitative and quantitative analysis of the sample processed [23] No. of spectrum Spectrum 1 Spectrum 2 Spectrum 3 C 0.56 0.69 0.71 Na – 0.38 – Si – – 0.31 Cr 0.66 0.66 0.54 Mn 0.48 – 0.54 Fe 95.73 95.75 95.48 Ni 2.57 2.52 2.42 Total 100.00 100.00 100.00 Temper increasing (Table 3), relative to the initial state, is because charged carbon ions [22], saturating the surface of the sample, are generated in the plasma layer of electrical gas discharge from the water solution of soda ash Na2CO3 when electric current flows. In water solution Na2CO3, such ions as 2Na, CO32-, -, occur. hen passing through the holes in the anode, anions ease away excess electrons. The Table 4. The phase composition of 18CrNi3Mo-Sh steel samples a) b) Fig. 7. X-ray diffraction pattern of 18CrNi3Mo-Sh steel; .) in the initial state, b) after EPP strengthening cations are entrained by the hydrodynamic electrolyte flow and recombine on the cathode (the products surface). All the results in Table 3 are presented in weight fractions and expressed in percentage terms. -ray structural analysis of 18CrNi3Mo-Sh steel samples in the state, as received (Fig. 7a), and after EPP strengthening according to mode B (Fig. 7b) revealed the a phase line, based on Fe, line Cr0.6 Fe1.4 phase, and also the line Fe2.7 Mo0.8 Ni0.1 phases. Processing type Phase composition 2Theta [deg] d [A] h k l l [%] 44.677 2.0267 1 1 0 100 a-phase 65.028 1.4331 2 0 0 11.5 82.344 1.1701 2 1 1 17.4 63.845 1.456730 3 0 2 17 Initial Cr0.6 Fe1.4 67.048 1.3947 2 0 5 13 61.739 1.5013 2 1 3 21 After EPP 96.981 1.0286 4 0 0 6 Fe2.7 Mo0.8 Ni0.1 98.149 1.0195 4 0 1 1 99.472 1.0094 2 2 4 20 44.750 2.0235 1 1 0 10.8 Fe3C 65.108 82.444 1.4315 1.1689 2 2 0 1 0 1 10 37 99.064 1.0125 2 2 0 14 Functional Surface Layer Strengthening and Wear Resistance Increasing of a Low Carbon Steel by Electrolytic-Plasma Processing 547 3 DISCUSSION After EPP, an expansion of intensity acceleration and diffraction lines can be observed (Fig. 7b) relative to the initial state (Fig. 7a). This indicates a stress state due to thermal exposure. The phase composition of samples of 18CrNi3Mo steel is presented in Table 4. As is known from [24], while hardening, the martensitic (A M) transformation is not fully completed, and the decay products are left in steel. Lines after EPP, a phase of Fe3C residual cementite and a-phase based on Fe, line Cr0.6 Fe1.4 phase, line Fe2.7Mo0.8Ni0.1 phase indicate martensitic hardening. According to the Kurdjumow-Sachs theory [25], a martensitic crystal appears on the forming shear area. Stresses play a very important role. Sources of stress are heat gradients along the cross-section, the heterogeneity of chemical composition, structural imperfections, different spatial crystal orientations, the specific volumes of austenite and martensite, and various coefficients of phase linear extension. There is a fine-grained lamellar structure of a fine needle-like martensitic class on the cross-section cut of the electrolytic plasma-processed sample (Fig. 8). The microstructure to the processed surface (left) is characterized by a dark layer, which is the structural phase transformation formed under the influence of an ionized high plasma temperature [23]. The methods and methodology relating to Fig. 8 are mentioned in the patent specification [24]. Fig. 8. The microstructure of the cross section of 18CrNi3Mo-Sh steel after EPP strengthening (100 times increase) It is established that electrolytic plasma processing makes it possible to obtain a strengthened layer of thickness from 1000 m to 1700 m (Fig. ). Microhardness at the cross-section of the strengthened area (martensitic structure) ranged from 6500 MPa to 7200 MPa. The results of microhardness measurements on the surface of a heat-treated sample in a plant (VKM JSC) averaged 7000 MPa (Fig. ). Microhardness decreases to the initial state and averages 3000 MPa when moving away from the processed surface [26]. The microrelief on the worn EPP surface of the steel sample after wear testing (Fig. 10a) indicates the least abrasive wear. ith fine-grained martensite, a strong surface structure is formed that cannot be broken even when worn (during the experiment) by severe forms of abrasive wear, thereby preventing large areas of destruction. Fig. 9. The value of the microhardness of 18CrNi3Mo-Sh steel after EPP strengthening The worst wear can be observed on the initial untreated steel sample (Fig. 10c). The pearlite-ferritic structure is easily abraded by abrasion, confirming the low wear resistance of thermally unprocessed steel. A steel sample that was thermally treated at JSC VKM wears out less than the initial sample (Fig. 10b). However, in terms of wear resistance, it is significantly inferior to the sample strengthened with the EPP method. Asample after EPPstrengthening has better wear resistance. As can be seen from Table 5, the steel resistance to abrasion increased significantly after EPP strengthening. Table 5. Value of wear of 18CrNi3Mo-Sh steel samples Processing type Hv [MPa] Wear [mg/h] EPP strengthening 6817 54.4 Thermo processed at JSC “VKMZ” 6298 100.4 State as received 3271 150.0 for abrasive wear resistance, 100 times increase; .) after EPP strengthening, b) thermo processed at VKMZ JSC, and c) initial state Based on the data obtained, it can be concluded that electrolytic plasma processing significantly increases the wear resistance of drill bit component parts at low labour intensity, which reduces the cost of the product as a whole [27]. The practical significance of the obtained results of scientific research in the field of working surface hardening by the EPP method means that the developed technology applies to a wider range of parts. In particular, they include the vast majority of the articulated elements of machines that work as friction pairs. Aspecial feature of the operating modes of articulated joints is the intense relative movement of the mated working surfaces, high contact load, unsatisfactory lubrication conditions of the friction surfaces due to the extrusion of the lubricant, the ingress of corrosive media into the mating zone, as well as abrasive particles. Such conditions are typical for loading and delivering vehicles when working in underground mines and open pit mining. The experimental research was conducted to ensure the operability of parts on the example of the rotary mechanism hinge assembly of the underground loader Caterpillar R1300G confirmed the possibility of effective use of the developed technology to strengthen the surface parts based on the EPP method. Thus, it has been stated by experiments that controlling EPP modes makes it possible to affect the uality of steel strengthening and obtain wear-resistant surfaces and significantly improve EPP productivity. This fact suggests the controllability of the EPP technological process and the possibility of practical implementation of the developed technology in the production 4 CONCLUSIONS The effect of electrolytic plasma processing on phase transformations and structural changes and the properties of low-carbon and compound steels has been theoretically and experimentally studied. The microstructure of the processed surface of steel samples is characterized by a dark modified surface layer. A fine needle-like structure of martensitic origin under the dark layer changes into the initial pearlite-ferritic structure. The total thickness of the strengthened surface layer is 1000 m to 1700 m. The microhardness of the strengthened area of a martensitic structure is 6500 MPa to 7200 MPa. The structure of 18CrNi3Mo-Sh steel sample in the initial state is composed of particles of a–phase based on Fe, Cr0.6 Fe1.4 phase, as well as Fe2.7 Mo0.8 Ni0.1 phase after electrolytic plasma processing in the structure of 18CrNi3Mo-Sh steel samples there appear particles of residual cementite Fe3C, a-phase based on Fe, Cr0.6 Fe1.4 phase, Fe2.7 Mo0.8 Ni0.1 phase, that testifies availability of martensite hardening. The most effective surface strengthening method of drill bit parts is electrolytic plasma processing. The advantages of this method are low energy consumption at high speeds of hardening, the possibility of local surface processing of parts of complex shapes operating under conditions of intense loads and the simplicity of the process. 5 ACKNOWLEDGEMENTS This research was funded by the Science Committee of the Ministry of Education and Science of the Republic of Kazakhstan within the framework of the grant project IRN AR0058518. 6 REFERENCES [1] Rembeza, A.I. (ed.) (1986). Reliability and Efficiency in Technology. Methodology. Organization. Terminology. Maszynostroienie, Moscow. (in Russian) [2] Sokolski, P., Smolnicki, T. (2021). A method for monitoring the technical condition of large-scale bearing nodes in the bodies of machines operating for extended periods of time. Energies, vol. 14, no. 20, art. ID 6637, DOI:10.3390/en14206637. [3] Henli, A.J., Kumamoto, H. (1984). Reliability of Technical Systems and Risk Assessment. Translated from English. Syromjatnikova, V.S., Demina, G.S. (eds.): Machine Industry, Moscow. (in Russian) [4] Gui, W., Lin, J., Hao, G., Qu, Y., Liang, Y., Zhang, H. (2017). Electrolytic plasma processing-an innovative treatment for surface modification of 304 stainless steel. Scientific Reports, vol. 7, art. ID 308, DOI:10.1038/s41598-017-00204-w. [5] Davis, J.R. (ed.) (2002). Surface Hardening of Steels: Understanding the Basics. ASM International, Materials Park. [6] Nestler, K., Bttger-Hiller, F., Adamitzki, W., Glowa, G., Zeidler, H., Schubert, A. (2016). Plasma electrolytic polishing - an overview of applied technologies and current challenges to extend the polishable material range. Procedia CIRP, vol. 42, p. 503-507, DOI:10.1016/j.procir.2016.02.240. [7] Jumbad, V.R., Chel, A., Verma, U., Kaushik, G. (2020). Application of electrolytic plasma process in surface improvement of metals: a review. Letters in Applied NanoBioScience, vol. 9, no. 3, p. 1249-1262, DOI:10.33263/ LIANBS93.12491262. [8] Lachowicz, M., Lesniewski, T., & Lachowicz, M. (2022). Effect of Dual-stage Ageing and RRA Treatment on the Three-body Abrasive Wear of the AW7075 Alloy. Strojniški vestnik - Journal of Mechanical Engineering, vol. 68, no. 7-8, p. 493-505, DOI:10.5545/sv-jme.2022.142. [9] Janicki, D., Musztyfaga, M. (2016). Direct Diode Laser Cladding of Inconel 625/WC Composite Coatings. Strojniški vestnik - Journal of Mechanical Engineering, vol. 62, no. 6, p. 363-372, DOI:10.5545/sv-jme.2015.3194. [10] Rakhadilov, B.K., Kozhanova, R.S., Popova, N.A., Nugumanova, A.B., Kassymov, A.B. (2020). Structural-phase transformations in 0.34C-1Cr-1Ni-1Mo-Fe steel during plasma electrolytic hardening. Material Science - Poland, vol. 38, no. 4, p. 699­706, DOI:10.2478/msp-2020-0073. [11] Kulikov, I.S., Vashhenko, S.V., Kamenev, A.J. (2010). Electrolyte-Plasma Processing of Materials. Bielaruskaja Navuka, Minsk. (in Russian) [12] Pogrebnjak, A.D., Kul’ment’eva, O.P., Kobzev, A.P., Tyurin, Yu.N., Golovenko, S.I., Bojko, A.G. (2003). Mass transfer and alloying processes during electrolytic plasma processing of cast iron. Letters to Technical Physics Journal, vol. 29, no. 8, p. 8-15, DOI:10.1134/1.1573301. [13] Rakhadilov, B., Zhurerova, L., Pavlov, A. (2016). Method of electrolyte-plasma surface hardening of 65G and 20GL low-alloy steels samples. IOP Conference Series: Materials Science and Engineering, vol. 142, art. ID 012028, DOI:10.1088/1757­899X/142/1/012028. [14] Rakhadilov, B.K., Kurbanbekov, Sh.R., Kylyshkhanov, M.K., Kenesbekov, A.B. (2018). Changing the structure and phase states and the microhardness of the R6M5 steel surface layer after electrolytic plasma nitriding. Eurasian Journal of Physics and Functional Materials, vol. 2, no. 3, p. 259-266, DOI:10.29317/EJPFM.2018020307. [15] Gorelik, S.S., Skakov, J.A., Rastorguev, L.N. (2002). X-ray and electron-optical analysis, (edition 4, updated and revised). Moscow Institute of Steel and Alloys, Moscow. (in Russian) [16] Geller, J.A., Rahshtat, A.G. (1989). Materials Sciences, (edition 6, updated and revised). Metallurgy, Moscow. (in Russian) [17] Skakov, M., Zhurerova. L., Scheffler, .. (2012). Way of hardening surface coating of details from steel 30CrMnSi in electrolytic plasma. Key Engineering Materials, vol. 531­532, p. 178-181, DOI:10.4028/www.scientific.net/KEM.531­532.178. [18] Skakov, ..K., Uazyrkhanova, G.K., Popova, N.A., S.heffler, M. (2012). Influence of heat treatment and deformation on the phase-structural state of steel 30CrMnSiA. Key Engineering Materials, vol. 531-532, p.13-17, DOI:10.4028/www.scientific. net/KEM.531-532.13. [19] Skakov, M., Rakhadilov, B., Scheffler, M., Batyrbekov, E. (2015). Microstructure and tribological properties of electrolyte plasma nitrided high-speed steel. Materials Testing, vol. 4, no. 57, p. 360-364, DOI:10.3139/120.110709. [20] Kombaev, K.K., Kylyshkanov, M.K., Lopuhov, J. (2009). Influence of electrolytic plasma processing of 18CN3..-ShSh steel on the surface microstructure and hardness. Journal of the Siberian Federal University, Engineering and Technology, vol. 2, no. 4, p. 394-399, handle/2311/1562, accessed on 2022-08-29. (in Russian) [21] Kozha, E., Smagulov, D.U., Ahmetova, G.E., Kombaev, K.K. (2017). Laboratory installation for electrolytic-plasma treatment of steel. News of the National Academy of Science of the Republic of Kazakhstan, vol. 4, no. 424, p. 219-224. [22] Tjulenin, A.N., Tjurin, Ju.N., Gradnev, A.I. (1988). Hysteresis of current-voltage characteristic of power supplies. MITOM. .oscow, no. 1, p. 9-12. [23] Kombaev, K.K., Kozha, E., Smagulov, D.U., Sadeh, B. (2016). Structural phase transitions of low-carbon alloy steels during electrolytic-plasma processing. Proceedings of the 2nd International Conference on Artificial and Industrial Engineering, p. 491-495, DOI:10.2991/aiie-16.2016.114. [24] Kombaev, K.K., Skakov, M.K., Putrenko, N.F. (2010). Patent #649, Republic of Kazakhstan. Test Stand for Testing Materials for Friction and Wear. Patent of the State Treasury Enterprise “Eastern Kazakhstan State Technical University named after D. Serikbayev “Ministry of Education and Science of the Republic of Kazakhstan” No. 2010 / 016.2. Kazakhstan, Nur-Sultan. [25] Kylyshkanov, M.K., Kombaev, K.K., Pogrebnyak, A.D. (2009). Patent #23178, of Republic of Kazahstan. Method of Electrolytic Plasma Strengthening of Drill Bit Component Parts, Nur-Sultan. [26] Kozha, E., Smagulov, D.U., Akhmetova, G.E., Kombaev, K.K. (2017). Laboratory installation for electrolytic plasma treatment of steel. News of National Academy of Sciences of the Republic of Kazakhstan, Almaty, vol. 4, no. 424, p. 219­225. [27] Rahadilov, B.K., Satbayeva Z., Baizhan, D. (2019). Effect of electrolytic-plasma surface strengthening on the structure and properties of steel 40KHN. Proceedings 28th International Conference on Mettallurgy and Materials, p. 950-955, DOI:10.37904/metal.2019.739. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 552-559 Received for review: 2022-05-04 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-07-05 DOI:10.5545/sv-jme.2022.184 Original Scientific Paper Accepted for publication: 2022-08-17 Industrial Exp erimental R esearch as a Contribution to the D evelopment of an Exp erimental Model of R olling Bearing Vibrations Mateusz rzochal1, – Stanisaw Adamczak1 – Grzegorz Piotrowicz2 – Sylwester nuk2 1 Kielce University of Technology, Poland 2 Fabryka Lozysk Tocznych-Krasnik S.A., Poland The influence of rolling bearing contamination on the level of generated vibrations is ignored in many cases. The influence is clearly visible in the theoretical and experimental works on the mathematical modelling of rolling bearing vibration. The authors of the most recent works in this field use substantial simplifications, often ignoring the presence of solid particles between the rolling bearing mating surfaces. However, bearing in mind the continuous improvement of the rolling bearings production processes, a series of exploratory tests in industrial conditions should be carried out, taking into account the factors affecting the bearing vibrations that are crucial in terms of quality control of the manufactured bearings. The analysis presented in this paper show that the cleanliness of the bearings is a critical factor that determines whether the rolling bearings meet quality requirements. Keywords: rolling bearings, vibrations, bearings cleanliness, quality control Highlights • Bearing contamination is an important factor in deciding whether a manufactured bearing is positively evaluated during final vibration inspection. • Ignoring the cleanliness condition in the mathematical models of bearings or experimental tests is not justifiable. • Bearing manufacturers focused on the development of precision bearings should, in addition to investing in the manufacturing process, develop the process of ensuring technical purity in parallel. • It is necessary to conduct detailed research on the analysis of the influence of the degree of rolling bearing contamination on its vibration and to study the degree of detectability of this type of incompatibility by applied instruments and measurement methods in the bearing industry. 0 INTRODUCTION Various mathematical models of rolling bearing vibrations can be found in the available publications. These models are oriented toward various applications important from the point of view of the model the authors present. The theoretical studies [1] carried out thus far on the mathematical models of rolling bearings enabled us to observe that the models are constructed with far-reaching simplifications and refer only to a few factors affecting the bearing level of vibrations [2] and [3]. There is no mathematical model that includes the vast majority of real factors, and the contamination of bearings and grease is one of the most commonly ignored factors in vibration analysis [4]. In addition, most analytical models focus on a bearing with an obvious defect that may appear only as a result of a prolonged operation or inappropriate conditions of operation [5] and [6]. Therefore, the models proposed in the latest literature on the subject are not universal and cannot be applied in industrial practice [7] and [8]. The analysis of the recently published studies shows that even if contamination is included in the experimental studies, they are mainly focused on grease contamination as a result of bearing operational use [9] and [10]. Therefore, the possible contamination in the process of production and its impact on the result of uality control is definitely a separate problem [11]. The construction of a model aimed at typical industrial applications is important because there are no reference bearings (which define the correct value of the vibration level of a bearing with specific characteristics), and even more so because there are no exact reference systems with which we could compare the results obtained on individual industrialsystems for vibration measurement that differ (e.g., in design) [12] and [13]. Therefore, efforts should be made to create an experimental model for brand-new bearings. To create such a model, a series of tests should be performed, including all factors influencing the vibration level, with the use of multicriteria statistics. The desired model should consider not only the imperfections affecting discrete vibrations (deviations in shape, size and position, excessive waviness, micro waviness and roughness) but also and above all the presence of solid particles in the bearing, which is an inherent result of the production process [14] and [15]. The basic parameters determining the uality of rolling bearings include durability (operating time in given conditions), rotation accuracy (the amount of the momentary deviation of the bearing from the operating position), rolling friction moment (mainly determining the efficiency of the device the bearing is installed in, which affects its energy efficiency), noise level (important from the everyday operation), as well as the level of vibrations (which affects the behaviour of the entire structure) [16] and [17]. The mutual proportionality of the vibration level and other indicators defining the bearing uality allows, in a simplified form, to make the bearing uality dependent on the vibration level itself (slight vibrations:high uality; significant vibrations:low uality) [18] to [20]. This uality criterion means that every bearing (without exception) is subjected to vibration control in bearing production plants. Other parameters are checked randomly or at the customers reuest [21]. It can also be noted that among the applied diagnostic methods, the measurement of vibrations gives the best results and is used much more often than, for example, acoustic- or temperature-based uality control [22]. Therefore, the importance studying the vibrations of rolling bearings is justified. order to understand the position of cleanliness control in the overall bearing acceptance inspection, a flow chart was prepared, as shown in Fig. 1. 1 VIBRATION MEASUREMENT ON THE PRODUCTION LINE Two types of measuring devices are used to measure the level of vibrations generated by bearings in industrial conditions: automatic and semi-automatic measuring devices. Shown in Fig. 2a, automatic measuring devices are generally used in automatic assembly lines, while semi-automatic devices are used to control small series of dismantlable bearings (Fig. 2b). The main components of such devices are a measuring unit with a vibration sensor, a measuring spindle (hydrodynamic, hydrostatic or pneumatic) and its drive, a pressure unit of the tested bearing, a bearing feeding system to the measuring stand (manual or automatic), a supporting structure with shock absorbers, and a control system. The tested rolling bearing is mounted on a shaft rotating in the spindle (the spindle has a universal socket for shafts adapted to different types of bearings). The correct vibration measurement reuires an axial load of the tested bearing. The load is introduced by pressing the outer ring with a force depending on the type of bearing. Usually, the clamps are adjustable, replaceable or, within certain limits, are universal to match the size of the tested bearing. Radial vibrations of the running and loaded bearing are recorded by an electrodynamic vibration velocity sensor with a direct contact with the stationary outer ring. The sensor is mounted in a holder and can be moved along two axes. In the variant for measuring vibrations on a production line shown in Fig. 2b, the bearings move in turn through a trough where they are placed on the measuring shaft by a pneumatic pressure. In the case of the manual device shown in Fig. 2b, the operator places the bearing on the shaft by applying the pusher. The above-mentioned measuring devices allow to measure the level of vibrations in three freuency bands: Low (50 Hz to 300 Hz), Medium (300 Hz to 1800 Hz) and High (1800 Hz to 10 000 Hz). The freuency bands listed above are related to the constant spindle rotational speed of 1800 min–1 at which all tests are carried out. This rotational speed is the reference speed and is determined by the standards for industrial measurement of bearing vibration levels (i.e., ISO 15242-1:2015 [23], ISO 15242-2:2015 [24], ISO 15242-3:2015 [25]). During the measurement, the effective value of the vibration velocity (RMS) and the maximum value (peak) are determined for each freuency band. The vibration level is expressed in different units, depending on the clients reuirements: dB, m/s, m/s, A (anderon) or (the conventional unit of one of the clients). a) b) Fig. 2. Industrial vibration measurement devices; a) in an automatic production cycle directly on the production line, and b) in a semi-automatic/manual production process In the automatic cycle, the bearing is measured twice, and the measurement is taken on both sides of the bearing; in the semi-automatic/manual cycle, the bearing is measured four times, with two measurements on each side of the bearing, every 0 in relation to the axis of measurement of the bearing. During semi-automatic/manual measurement, the bearing is additionally checked for sounds it emits. The bearing should emit a hum appropriate for a given bearing. Audible disturbances such as rattling, whistling, or crackling caused by contamination, damage to the race of the inner ring, outer ring, or rolling elements are unacceptable. If the device is not euipped with a vibration spectrum analyser, oscilloscopes are additionally used for measurements. The image of the vibration amplitude observed on the oscilloscope should be uniform and should not show peaks and irregularities of waves higher than one half of the baseband width for standard bearings and one third of the baseband width for bearings with a reduced level of vibrations (i.e., electric bearings (C66)). These bearings (C66) are 100 vibration controlled, while for standard bearings, the measurement freuency is 3 for each assembly lot but not less than 10 pieces. The RMS value is displayed on the screen of the device in numerical form and additionally visualized in the form of bar graphs. For the bearing being measured, the RMS and peak´ limits for each band are entered before starting the measurements. The comparison of the four values obtained with the values set as the limit is decisivefor the recognition of the bearing as meeting the uality reuirements. Fig. 3 shows two different measurement results. Fig. 3a indicates the good uality of the tested rolling bearing: the vibration level does not exceed the assumed limits, and all bars are green). Fig. 3b shows the result of insufficient uality of the tested rolling bearing: the assumed limits are exceeded in the medium and high range (the second and third bars turn red). The vibration level generated by the bearings is measured on open and empty bearings (without grease and seals). Such bearings are protected only with a thin film of oil. If the exposed bearing meets the uality reuirements regarding the vibration level, then the closed (i.e., lubricated and sealed) bearing will meet such conditions all the more (the grease dampens vibrations in the bearing). The customer may wish to measure the vibration level in a closed bearing. In this case, the drive of the measuring device must be more powerful, and the contact pressure must be greater. The increase of these parameters (higher power of the drive and higher bearing pressure force) may accelerate the wear of the measuring device. This risk increases with the diameter and the weight of the tested bearing. The result of the vibration level test proves its technical perfection. By making this measurement, we obtain information about potential defects of individual bearing components (outer ring, inner ring, rolling elements, cage) and possibly about impurities a) in the bearing. The results of the vibration level measurement are subject to continuous analysis and, depending on the result, the uality controllers react in accordance with the procedures in place. hether the defective bearings with exceeded vibration levels are isolated cases or there is an error in the process of production should be investigated. Dirty bearings are rewashed. Bearings with exceeded vibration levels that do not show signs of significant damage, are reclassified as lower-class bearings and still offered for sale. Bearings with exceeded vibration levels that show symptoms of significant damage are dismantled. Bearing components to be dismantled: inner rings, outer rings, and rolling elements are subject to microscopic inspection. The results of this inspection have an impact on the next steps. If the problem is isolated, the production procedures should be analysed, but if more pieces are affected, the root cause of the problem should be identified and resolved. Defective bearing components are subject to regeneration if possible, and if not, they are scrapped. The detected deficiencies have an impact on the improvement of the technological process. Each improvement in the technological process upgrades the parameters of the manufactured bearing. All the improvements automatically influence the process of design and production of new, more and more accurate devices for measuring the vibration level. 2 MAIN CAUSES OF NON-CONFORMITIES FOUND BY THE QUALITY CONTROL Contamination of the lubricating layer with solid particles is an important parameter affecting the durability of rolling bearings. The degree of life reduction caused by solid particles in the lubricating layer depends on the type, size, hardness, and uantity of particles, as well as the thickness of the lubricant layer (viscosity) and the size of the bearing. Durability tests, however, are destructive and time-consuming tests. The level of vibration generated by a rolling bearing is a critical parameter and is inversely proportional to service life. Bearing vibration tests are much less time-consuming and do not damage the bearing. However, both in theoretical models of bearing vibrations and in experimental tests, the cleanliness of the bearing interior is largely ignored. This is due to the fact that, firstly, the contamination is very difficult to simulate mathematically, and, secondly, the experimental studies are focused on high-amplitude components in the low and medium freuency range initiated by damage or surface geometry. To discover the most important causes of rejects, a six-month observation of the production process of rolling bearings of ten different types was performed. A total of 46,811 bearings were tested. Products that were rejected from the production process during the uality control involving the measurement of vibrations were thoroughly analysed. The causes of the defects of a given product were grouped into four categories: Defects in the rolling element related to an incorrect shape of one or more rolling elements, a significant surface defect, or an unsuitable size. Defects of the inner or outer ring related to the inadeuate geometric structure of their surfaces or shape errors, less often to a local defect. Contamination was understood as residues left in the process of production. In most cases, after thorough washing, the bearing met the uality reuirements. Other defects not related to the rolling element rings and contamination (i.e., cage defect or exceeding the limit of one of the three freuency bands) not related to any of the potential causes. The test results are presented in Table 1. Table 1. Percentage of specific causes of rejects of rolling bearings from the production line Reject rate [%] Type Ball defect Outer/ inner ring defect Contamination Other 6016Z 53 16 31 0 63007-2RS 41 18 35 6 6311Z 64 21 12 3 6313P63 80 8 12 0 6317-2RS 52 18 20 9 6208Z 60 5 13 21 6211Z 84 16 0 0 6217-2RS 40 24 20 16 6219Z 66 21 12 2 6405 35 36 28 1 The data presented in Table 1 show that the main cause of the rejection due to the level of vibration was a ball defect. This is a special situation because this type of defect is a critical defect that, if occurs, disualifies the bearing even though it does not exceed the limits on any of the analysed freuency bands. The reason for this situation was that the observation coincided with the delivery of defective balls deliveries from one of the suppliers. Therefore, further evaluation of the bearing vibration level was focused mainly on the detection of defective balls. After the defect was discovered, a complaint procedure was initiated in order to find the root cause, and the collected materials served as evidence confirming the scale of the problem. In addition to the defects of the balls in the observed production process, the contamination of the bearings was also a significant problem that significantly affected the level of vibration and made it difficult to assess the bearing, to capture single defects on the rings or the aforementioned defects on the balls. Therefore, when conducting research aimed at analyzing the vibrations of rolling bearings in the industrial environment, it is impossible to ignore the contamination, which constitutes a significant cause of all discrepancies. It is very important, especially for high-precision bearings designed for high-speed applications for which the evaluation of the vibration level is an important parameter determining the uality of the bearings. 3 TESTING CLEANLINESS OF BEARINGS IN INDUSTRIAL CONDITIONS As shown by the presented long-term tests on several types of bearings, contamination is a significant cause of discrepancies discovered during uality control in bearing production plants. Bearings that are suspected of not meeting the cleanliness conditions set for them are re-washed in accordance with the procedure in place. Bearing samples, both contaminated and clean, are sent to the laboratory. There, the degree of pollution and its structure is determined. The test is performed using weight analysis and microscopic analysis. The weight analysis consists of washing the bearing parts in extraction naphtha using ultrasound. The liuid after washing and rinsing the bearing parts is filtered through a previously soaked (in extraction naphtha), dried, and weighed nylon membrane with a mesh side of 0.8 m. The filtrate with the bearing residue is dried and then weighed to the nearest 0.1 mg. The difference in weight of the filter after and before filtration represents the weight of the contaminations. The resulting contaminant mass is compared to the reuirements. Table 2. The result of the weight analysis of the contaminants found in the bearings Contaminants above Bearing Mass of contaminants Bearing 300 µm Quantity type, [mg/kg] number (max. Dimension) quantity Requirements Result Requirements Result 1. 6208Z 6.6 0 max max 2. 6208Z 4.8 1 (420) 5 mg/kg 300 µm 3. 6208Z 3.7 0 After weight analysis, the filters are examined using a microscope. Microscopic analysis involves the computerized counting of particle sizes divided into five ranges: 25 m to 50 m, 50 m to 100 m, 100 m to 150 m, 150 m to 200 m, and 200 m to 300 m. Contaminants below 25 m are disregarded, while those above 300 m are not acceptable. The counting of contaminants on the surface of the dried nylon membrane is carried out at a zoom of 100. The test allows a specific cleanliness class to be assigned according to the standards described in ISO 16232:2018 [26]. Examples of the results of the weight analysis are given in Table 2. Table 3 presents an example of the result of the microscopic analysis of the tested bearings contaminants, while Fig. 4 shows the filters and the captured foreign bodies. Table 3. The result of the microscopic analysis of the contaminants found in the bearings Size ranges of Amount of metallic contaminants [pcs] contaminants 25-50-100-150-200->300 (max. class [µm] 50 100 150 200 300 dimension) 12177815 4 4 0 8 Bearing 6208Z no. 36230114 0 0 7 ofc - out of class In order to ensure the appropriate cleanliness class of the bearings, modern devices with extensive filtration systems are increasingly designed for washing individual bearing elements, as well as assembled bearings. Thus, thanks to accurate systems used for measuring the level of vibration generated by bearings, we are able to determine the cause of the increased vibration level. e can also determine the reason for the problem and which stages of the technological process of bearing production can be improved. 4 CONCLUSIONS Contamination of the lubricating layer with solid particles is a critical problem that reduces the life of rolling bearings. The degree of reduction in service life caused by solid particles in the lubricant layer depends on the type, size, hardness, and uantity of the particles, the thickness of the lubricant layer (viscosity) and the size of the bearing [27] and [28]. However, durability tests are destructive tests and are very time-consuming. The level of vibration generated by a rolling bearing is a critical parameter that is inversely proportional to service life. Bearing vibration tests, however, are significantly less time-consuming and do not damage the bearing. However, both in theoretical models of bearing vibration and in experimental studies, bearing cleanliness is largely neglected. This is due to the fact that, firstly, the contamination is very difficult to simulate mathematically, and, secondly, the experimental studies are generally focused on high-amplitude components in the low and medium freuency range due to the damage or the geometric structure of the surface. Fig. 4. Examples of microscopic photos of the filters showing contamination of the rolling bearings The most important conclusion from the descriptions given in the paper is that the contamination in bearings is a very important factor that determines the positive assessment of the manufactured bearing during the final control of the vibration level. Therefore, ignoring this factor in mathematical models or experimental research is unjustified, especially since the industry trends in downsizing and striving for ever-higher operating performance make rolling bearings operate in increasingly difficult conditions with greater precision of rotation, high operating speed, minimizing noise, vibrations, and heat. Therefore, bearing companies focused on the development of precision bearings, in addition to investing in the production process, should simultaneously develop the process of ensuring technical cleanliness. Providing and maintaining perfect cleanliness in the production process is a real challenge. In companies producing parts for the automotive industry, the technical cleanliness of components in accordance with the reuirements of ISO 16232 and VDA1 is one of the most important issues. Due to the reuirements, not only the production process and the washing process affect the final level of technical cleanliness of the components that are delivered to the final assembly. Of key importance will also be the level of cleanliness of the rooms in which the bearing production and assembly operations are performed, as well as the packaging and measures to prevent re-contamination. In the future, the authors of this paper plan to conduct detailed research on the analysis of the impact of the degree of contamination of the rolling bearing on its vibrations and to test the degree of detection of this type of non-compliance by the instruments and measurement methods used in the bearing industry; as a result of the planned tests, it will be possible to predict the importance of bearing cleanliness to limit the level of vibrations. In turn, the planned research will be part of the preliminary research for broader studies related to an attempt to develop a model based on multi-criteria statistics that could be used to forecast vibrations of a newly manufactured bearing in the industrial environment. A model built on the basis of real data collected in a rolling bearings production plantwill undoubtedly have great practical significance. 8 REFERENCES [1] Wrzochal, M., Adamczak, S. (2020). The problems of mathematical modelling of rolling bearing vibrations. Bulletin of the Polish Academy of Sciences: Technical Sciences, vol. 68, no. 6, p. 1363-1372, DOI:10.24425/bpasts.2020.135398. [2] Shi, P., Su, X., Han, D. (2016). Nonlinear dynamic model and vibration response of faulty outer and inner race rolling element bearings. Journal of Vibroengineering, vol. 18, no. 6, p. 3654-3667, DOI:10.21595/jve.2016.16944. [3] Deák, K., Mankovits, T., Kocsis, I. (2017). Optimal wavelet selection for the size estimation of manufacturing defects of tapered roller bearings with vibration measurement using Shannon entropy criteria. Strojniški vestnik - Journal of Mechanical Engineering, vol. 63, no. 1, p. 3-14, DOI:10.5545/ sv-jme.2016.3989. [4] Liu, J., Shao, Y. (2015). Vibration modelling of nonuniform surface waviness in a lubricated roller bearing. Journal of Vibration and Control, vol. 23, no.7, p. 1115-1132, DOI:10.1177/107754631 5589675. [5] Do, V.T., Nguyen, L.C. (2016). Adaptive empirical mode decomposition for bearing fault detection. Strojniški vestnik - Journal of Mechanical Engineering, vol. 62, no. 5, p. 281-290, DOI:10.5545/sv-jme.2015.3079. [6] Kong, F., Huang, W., Jiang, Y., Wang, W., Zhao, X. (2019). Research on effect of damping variation on vibration response of defective bearings. Advances in Mechanical Engineering, vol. 11, no. 3, p. 1-12, DOI:10.1177/1687814 019827733. [7] Soni, A., Pateriya, A., Dewada, S. (2018). Study the effects of solid and liquid contamination in ball bearing through vibration analysis. International Research Journal of Engineering and Technology, vol. 5, no. 5, p. 4397-4401. [8] Jamaludin, N., Mba, D., Bannister, R.H. (2001). Condition monitoring of slow-speed rolling element bearings using stress waves. Proceedings of the Institution of Mechanical Engineers, Part E: Journal of Process Mechanical Engineering, vol. 215, no. 4, p. 245-271, DOI:10.1177/095440890121500401. [9] Malla, C., Panigrahi, I. (2019). Review of condition monitoring of rolling element bearing using vibration analysis and other techniques. Journal of Vibration Engineering & Technologies, vol. 7, p. 407-414, DOI:10.1007/s42417-019-00119-y. [10] Dadouche, A., Conlon M.J. (2016). Operational performance of textured journal bearings lubricated with a contaminated fluid. Tribology International, vol. 93, p. 377-389, DOI:10.1016/j. triboint.2015.09.022. [11] Maru, M.M., Castillo, R.S., Padovese, L.R. (2007). Study of solid contamination in ball bearings through vibration and wear analyses. Tribology International, vol. 40, no. 3, p. 433­440, DOI:10.1016/j.triboint.2006.04.007. [12] Shah, D.S., Patel V.N. (2014) A review of dynamic modeling and fault identifications methods for rolling element bearing. Procedia Technology, vol. 14, p. 447-456, DOI:10.1016/j. protcy.2014.08.057. [13] Boškoski, P., Petrovcic, J., Musizza, B., Juricic, Ð. (2010). Detection of lubrication starved bearings in electrical motors by means of vibration analysis. Tribology International, vol. 43, no. 9, p. 1683-1692, DOI:10.1016/j.triboint.2010.03.018. [14] Xiong, Q., Zhang, W., Xu Y., Peng, Y., Deng, P. (2018). Alpha-stable distribution and multifractal detrended fluctuation analysis-based fault diagnosis method application for axle box bearings. Shock and Vibration, vol. 2018, art. ID 1737219, DOI:10.1155/2018/1737219. [15] Ouadine, A., Mjahed, M., Ayad H., El Kari, A. (2018). Aircraft air compressor bearing diagnosis using discriminant analysis and cooperative genetic algorithm and neural network approaches. Applied Sciences, vol. 8, no. 11, art. ID 2243, DOI:10.3390/ app8112243. [16] Kahaei, M.H., Torbatian, M., Poshtan, J. (2007). Bearing-fault detection using the meyer-wavelet-packets algorithm. Strojniški vestnik - Journal of Mechanical Engineering, vol. 53, no. 3, p. 186-192. [17] Li, J., Li, M., Zhang, J., Jiang, G. (2019). Frequency-shift multiscale noise tuning stochastic resonance method for fault diagnosis of generator bearing in wind turbine. Measurement, vol. 133, p. 421-432, DOI:10.1016/j. measurement.2018.10.054. [18] Adams, M.L. (2018). Bearings: Bearings Basic Concepts and Design Applications. 1st ed., CRC Press, Boca Raton, p. 92­123, p. 212-231, DOI:10.1201/b22177. [19] Harris, T.A., Kotzalas, M.N. (2006). Rolling Bearing Analysis Essential Concepts of Bearing Technology, 5th ed., CRC Press, Boca Raton, p. 329-348, DOI:10.1201/9781420006582. [20] Taylor, J.I. (2003). The Vibration Analysis Handbook, 2nd ed. VCI, p. 167-221. [21] Pan, H., He, X., Tang, S., Meng, F. (2018). An improved bearing fault diagnosis method using one-dimensional CNN and LSTM. Strojniški vestnik - Journal of Mechanical Engineering, vol. 64, no. 7-8, p. 443-452, DOI:10.5545/sv-jme.2018.5249. [22] Delgado-Arredondo, P.A., Moríñigo-Sotelo, D., Osornio-Rios, R.A., Avina-Cervantes, G.J., Rostro-Gonzales, H., Romero-Troncoso R.J. (2017). Methodology for fault detection in induction motors via sound and vibration signals. Mechanical Systems and Signal Processing, vol. 83, p. 568-589, DOI:10.1016/j.ymssp.2016.06.032. [23] ISO 15242-1:2015. Rolling bearings - Measuring methods for vibration - Part 1: Fundamentals. International Organization for Standardization, Geneve. [24] ISO 15242-2:2015. Rolling bearings - Measuring methods for vibration - Part 2: Radial ball bearings with cylindrical bore and outside surface, International Organization for Standardization, Geneve. [25] ISO 15242-3:2006. Rolling bearings - Measuring methods for vibration - Part 3: Radial spherical and tapered roller bearings with cylindrical bore and outside surface. International Organization for Standardization, Geneve. [26] ISO 16232:2018. Road vehicles - Cleanliness of components and systems. International Organization for Standardization, Geneve. [27] Espejel, G. M., Gabelli A., Ioannides S. (2010). Lubrication and Contamination effects on bearing life. Evolution, vol. 2, p. 25-30, from https://evolution.skf.com/wp-content/ uploads/2013/02/EN_lubrication0.pdf, accessed on 2022­07-28. [28] Espejel, G. M., Gabelli A., Ioannides S. (2010). Lubrication and Contamination effects on bearing life, part 2. Evolution, vol. 3, p. 25-31, from https://evolution.skf.com/wp-content/ uploads/2014/10/EN_lubrication_part_2.pdf, accessed on 2022-07-28. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 560-570 Received for review: 2022-02-21 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-07-01 DOI:10.5545/sv-jme.2022.69 Original Scientific Paper Accepted for publication: 2022-07-11 L imit-protection Method for the Wor kspace of a Parallel Pow er H ead Yanbing Ni1,– enliang Lu2 – Shilei Jia2 – Chenghao Lu2 – Ling hang2 – Yang en1 1 Tianjin University, School of Mechanical Engineering, China 2 Key Laboratory of Mechanism Theory and Euipment Design of Ministry of Education, China The end pose of a one-translation two-rotation (1T2R) parallel mechanism is a mapping of the servo motor action in the joint space. Because it is difficult to obtain information about the end attitude state, we have designed and implemented a simplified algorithm for determining the attitude of such a mechanism. The kinematic inverse solution of the robot and the modelling analysis of the workspace are carried out. From this, it is deduced that the length transformation of the three branch chains of the mechanism reflects the position and attitude of the end-motion platform. Based on this algorithm, the limit protection of a parallel power head under arbitrary configuration is realized. The correctness of the calculation method is verified by simulation. Finally, based on the software and hardware conditions of an existing control system, experimental verification is carried out. The experimental results show that the simplified algorithm can implement limit protection for this type of machine. Keywords: parallel power head, modelling, workspace analysis, position and pose judgment, limit-protection Highlights • The kinematics of the 1T2R head are analysed, and the workspace of the mechanism is obtained. • A new simplified algorithm for position and attitude judgment is proposed. • The new simplified algorithm of position and attitude judgment is combined with the mechanism control system to realize a fast limit. • The algorithm is verified with simulation and experimental data. 0 INTRODUCTION Due to its compact structure, high stiffness to mass ratio, high precision, and good dynamic characteristics [1] and [2], the parallel power head has been of interest to both academia and the manufacturing industry and has been widely used in high-end manufacturing fields, such as large aircraft structural parts [3]. The working space is the key index for evaluating the performance of parallel power heads. Unlike tandem robots, which are driven by a series of connecting rods and rotating joints in series, parallel power heads are connected by at least two independent kinematic chains. Although the working space is relatively small, their structure is precise and compact with high repetitive positioning accuracy [4] and [5]. orkspaces can be partitioned in a number of ways, depending on performance reuirements and selected parameters [6]. The factors influencing the size and shape of the workspace include the constraints of the length of the branch chain, the rotation angle, and the size of the revolute joint [7] and [8]. The methods used for the analysis of the workspace include the geometric, numerical [5], and discretization methods [9] and [10]. Gui et al. [11] proposed a reliability mathematical model based on random probability and presented a measurement and calculation method for the evaluation of the reliability level of mechanism motion. Shao et al. [12] analysed the new spatial-planer parallel mechanism using geometric methods and verified that it has good accuracy and efficiency. Kaloorazi et al. [13] used the structural geometry method to determine the maximum non-singularity workspace of a 3-degree of freedom (3-DOF) parallel mechanism, and Huang et al. [9] obtained the workspace of the Stewart-Gough manipulator by a discrete method. The scientific community has also made many achievements in numerical method research. hang [14] et al. used the fast search method to calculate the workspace volume and took it as the optimization objective function. In addition, Gao and hang [15] and hu et al. [16] used the modified boundary search method to obtain the workspace of a parallel mechanism. Farzaneh-Kaloorazi et al. [17] used an interval analysis approach for the barrier-free working space of a parallel mechanism. Majid [18] et al. analysed the workspace of a three-prismatic­prismatic-spheric-revolute (3-PPSR) manipulator using the numerical method and found that three regions in the workspace corresponded to the postures of a type of manipulator. In order to make more reasonable and effective use of the existing working space and ensure the safe and reliable operation of the parallel power head, the method of trajectory planning and space limitation is often adopted [19] and [20]. Khoukhi et al. [21] proposed a multi-objective dynamic trajectory planning method for the parallel mechanism under the constraints of task, workspace, and mechanism with the help of the discrete augmented Lagrangian techniue. Reveles et al. [22] proposed a kinematic redundant parallel robot joint trajectory planning method using a feasibility map, which plans the joint trajectory while avoiding parallel singularities through the graphical evaluation of the robot pose related to four working modes. Dash et al. [23] proposed a numerical path planning method to avoid the singularity of the mechanism in the accessible workspace of a parallel robot by clustering the singularities and using the local routing method based on Grassmann line geometry to avoid the singularities. Ibal et al. [24] discussed two complex control strategies: computational torue control (CTC) and variable structure control (VSC) and improved the trajectory-tracking performance of the robot. Manzoor et al. [25] integrated the functions of mechanical computer-aided design (CAD) and robot CAD into the same platform and achieved the accurate control of the robot through various three-dimensional models in the platform. Alam et al. [26] considered two different methods based on sliding mode control (SMC) to achieve the nonlinear control of an elastic joint robot. This control method enables the robot to obtain a locally stable closed-loop system. Because the pose of the end of the one-translation two-rotation (1T2R) parallel mechanism in the operation space is a nonlinear mapping of the motion of the servo motor in the joint space, the modelling and calculation process is complicated, the motion controller takes a long time. It is therefore impossible to realize the real-time operation to prevent exceeding the limit during the operation of the mechanism. The purpose of this paper is to propose a simplified algorithm for position and attitude judgment, which can greatly reduce the calculation amount of limit, and realize the limit protection of the 1T2R parallel power head through a combination of software and hardware. Firstly, the inverse kinematics analysis was carried out, and the mapping relationship between the terminal pose of the mechanism and each input value was constructedbased on the inverse kinematics model. Then, according to the mechanism structure, scale parameters, range of motion of each pair, interference and other constraints, the working space of 1T2R power head was determined. By analysing and summarizing the motion rules of 1T2R power head, a simplified algorithm for judging the position and pose of 1T2R power head was obtained, which was used to realize the limit protection of the 1T2R head. Finally, the algorithm was verified via experimentation. 1 KINEMATIC ANALYSIS The position inverse solution of the 1T2R head is to solve the rod length uantity of each branch chain motion joint by knowing the positional parameters of the end reference point of the tool providing the theoretical model for mechanism error analysis and control. 1.1 Machine Tool Profile As shown in Fig. 1, a 1T2R mechanism is a parallel mechanism with one translational and two rotational degrees of freedom. A1T2R power head is composed of a moving platform, a static platform, and three RPS branch chains, in which R, P, and S represent revolute joints, active prismatic pairs and spherical joints, respectively. One end of each RPS branch chain is connected to the moving platform through a spherical joint, and the other end is connected to the static platform through a revolute joint. The motorized spindle is installed on the moving platform and the active prismatic pair is driven by a servo motor. Fig. 1. Model of 1T2R spindle head 1.2 Inverse Kinematic Solution The power head of the 1T2R mechanism diagram is shown in Fig. 2, A and B represent the centres of the spherical joint and the revolute joint, respectively; A 1A 2A 3 and B 1B 2B 3 are euilateral triangles; points O ' and O are the geometric centres of the two triangles; the moving platform is represented by A 1A 2A 3; the static platform is represented by the plane of B 1B 2B 3; point P is the tip at the end of the mechanism; e is the distance from P to the plane of the moving platform; a, b are the circumcircle radii of the moving and static platform, respectively; l represents the initial branch length. ii ' The fixed and moving coordinate systems O xy ' and O x y' are established at the centre O and O of ' the fixed and moving platforms, respectively. In the initial position, the axis and axis are perpendicular to the planes A 1A 2A 3 andB 1B 2B 3, respectively, as ' shown in Fig. 2. The x axis and x axis are along the ˜°˜°' ˜˜˜˜˜˜˜˜ AAdirection and BB, respectively. The axis 3232 and y axis are determined according to the right-hand rule. a and ß represent rotation about the x axis and y axis, respectively. y ' rP ˜°xP yPzP .T.(3) The following vector euation can be obtained: r ˜b °q wa ei ˜,,, .°w 123(4) P iiii where a, bis the position vector of A , B ; aRa0; ii iiii bi ˜b°cos.i sin.i 0.T; ai 0 is the measurement of A in the connected coordinate system of moving i platform, ai ˜a °cos.sin.0.T; fi is the 0ii ˜°˜˜˜˜° ˜˜˜ structural angle of 'i and OBi x axis OArelative to and x axis respectively, ˜i .2°.i .1ˆ 3.° 2; q and i = ' w are the length of the branch chain and the unit i vector, respectively. Since the revolute joint will restrict the branch chain from moving along the x axis, dot product x i axis direction unit vector uix =° cos ˜i sin ˜i T0.at both ends of E. (4): .r P a ˜°w e .u T.ˆ0.i ix (5) Solve E. (5) to obtain x P ˜ a sin°2ˆ2 ..°1cos...e sinsinˆ., a y ˜c os°2ˆ..°1c os...e cossˆiin., Fig. 2. Mechanism sketch of 1T2R mechanism P 2 ' ' ' The attitude transformation matrix of O x y' .˜.ˆ.(6) in the connected body coordinate system of a moving From E. (6) platform relative to O xy in the static platform coordinate system is R : R˜Rot°z ,..Rot°x ,ˆ.Rot°z ,.. q = r ˜b °a ˜e w = r ˜b ˜s i Pii P ii .ˆ.ˆ . (7) , r ˜b ˜s P ii w = i .,, 123 i q i  cc .scs ....ˆ.ss ...ˆ.cs .scc .ˆ sc ccs ....ˆ..cs ... . ss scc   .ˆ. .ˆ ss ccc = where sis the vector of the geometric centre of i the spherical pairs pointing to the tool reference point ˆ.ˆ.ˆ ' ' ' =uvw,(1) where s sin, c cos, u, v, w represent the measurement of three coordinate axes of O x y' in the static platform coordinate system O xy , respectively, and w can be used to represent the tool attitude vector. ., ., . are the precession angle, nutation angle and spin angle, respectively, which are related to a and ß:   ˆ... in the static coordinate system, and s is the vector of the geometric centre of the spherical pairs pointing to the tool reference point in the moving coordinate system. s ˜RRss ,˜°0ae .T,(8) ii cos .i °  ..sin .i ° 2  2 .0 R i ˜sin i ° . cos 2 .0, 2 ..° i 0 01 . °ˆ . ˜ . i ˜,, () 12 3 . . .. ° ˜ .˜° .  (2) cossin sin atan2 , 2 WORKSPACE ANALYSIS sin sin As shown in Fig. 2, the position vector of the The constraints of the 1T2R head include its structure, arccos coscos ..  tool point P at the end of the mechanism in the static scale parameters (a, b , l , e), motion range of each platform coordinate system O xy is: prismatic pair and interference. orkspace analysis Ni, Y. – Lu, W. – Jia, S. – Lu, C. – Zhang, L. – Wen, Y. Fig. 3. Interference position of 1T2R spindle head is used to determine the set of all reachable space position points of the tool reference point under the above constraints. 2.1 Constraint Analysis The 1T2R power head is restricted by its own mechanical structure and other conditions and can only work within a certain space range. After analysing the structure of the 1T2R head, the constraint conditions are as listed in Table 1. Table 1. Restrictions of 1T2R mechanism Constraint type Constraint value Active prismatic pair length constraint 0.4 m = q =0.915 m i Rotation angle constraint of revolute joint –12° = . = 3° i Angle constraint of spherical joint |ai = 45°, |ßi = 90° Interference between principal axis and d= 0.15 m i branched chain Interference between tool point and table = 1.250 m P q min and q are the maximum and minimum ii max lengths of the active prismatic pairs, respectively; .min i and .are the maximum and minimum values of i max each revolute joint angle, respectively; amin and ß ii max are the maximum values of angles a and ß of each spherical joint, respectively; dis the linear distance i from point S of the spindle end to the branch chain; di min is the minimum allowable actual linear distance between the spindle end and the branch chain. As shown in Fig. 3, d can be obtained as follows: i ˜i = wr °. w r wwb i °.P .e .l .i ˆ .(10) i BS i 2.2 Workspace Description A hierarchical processing approach can be used to divide the workspace into multiple subspaces. The boundary region of the subspace is then determined by means of the uasi-spherical coordinate search method [27], and the envelope surface and the stereogram of the workspace are described, defining the accessible working space for 1T2R heads as shown in Fig. 4. It can be seen from the calculation results that 1T2R head can realize the attitude space range with a maximum of . . 0 40, . . 0 360. As we can see from the reachable workspace, not all regions can achieve the maximum nutation angle. Fig. 4. Reachable workspace of 1T2R 3 POSE DETERMINATION AND SIMPLIFICATION ALGORITHM By analysing and summarizing the motion rules of 1T2R head, a simplified algorithm for judging the position and pose of 1T2R parallel power head is derived. This algorithm can realize the fast response of the limit protection function under the existing hardware configuration and is used to realize the limit protection of the 1T2R head. To facilitate calculation, the relationship between the coordinates of the tool reference point in the and the coordinates of the geometric centre point of the the static platform to the geometric centre point of the O' moving platform in the direction of is: moving platform. zz e°cos.(11) ˜O .. 3.2 Pose Constraint Condition Thus, the input parameters of the inverse kinematics model of the 1T2R head are replaced by the coordinates of the geometric centre point O' of the moving platform in the direction, the process angle . and the nutation angle .. Q is defined as the kinematic chain length of the branch chain movement; Q' is the difference between the kinematic chain lengths of any two branch chains. 3.1 Position Constraint Condition As shown in Fig. 2, let the vector from the geometric centre point O of the static platform to the geometric centre point O' of the moving platform be o. rcan be P expressed as r P o e (12) ˜°w. Substituting E. (12) into E. (3), we obtain: q w ˜a °b .o i ˜123. (13) ii ii ,, By summing the above euation in order of i , and taking into account the geometric relation of the prototype a1 a2 a3 b1 b2 b3 0, we obtain: q w +q w +q w = 3 o .(14) 112233 According to E. (13), the left side of the above euation is eual to a1 a2 a3 b1 b2 b3 0; therefore: q 1w 1 + q 2w 2 + q 3w 3 ˜ .(15) 3 o ° a 1.b 1 ° a 2.b 2 ° a 3.b 3 Because a1– b1 a2– b2 a3– b3 has a maximum value of 6a sin(.max/ 2) at the maximum nutation angle, the following can be obtained by sorting: a 1˜b 1 ° a 2˜b 2 ° a 3˜b 3 ..6a sin.ˆmax 2..(16) . max By combining E. (15) and E. (15), and substituting q 1w 1 q 2w 2 q 3w 3 Q , the following is obtained: 3 o ˜˜3o Q °6asin..max 2..(17) It can be seen from E. (17) that, theoretically, the sum of kinematic chains of the 1T2R power head branch prismatic pair is always approximately three times the distance from the geometric centre point of A similar calculation is used for the difference of length between any two branch chains. The length difference between branched chains 1 and 2 is used as an example, which can be determined from E. (13): q w ˜q w °a ˜a ˜b .b .(18) 11221212 The difference between the two sides of a triangle is less than the third: q w 1 ˜q w 12 2 °q w 1 ˜q w 12 .2 (1) Substitute E. (23) into, and we obtain: q w ˜q w 1122 °a ˜a ˜b .b ,1212 (20) because a1 aRV1, a2 aRV2, b1 b V1, b2 b V2, Vi = cosfsinf0T, i 1, 2. ii According to the actual parameters of the mechanism, a = b is substituted, and then: aabb ˜˜° . 1212 .a RE)T (˜ .21a .˜cos...1˜2ˆˆ R˜E ,(21) where, E is identity matrix, T cosf1–cosf2 sinf1–sinf2 0T. According to E. (1), the matrix R is related to the precession angle ., nutation angle . and spin angle ., and further: a ˜a ˜b °b .a ..ˆ 21˜cos23..f ..,...(22) 1212 Substituting E. (22) into E. (21), we obtain: q w ˜ q 2w 2 °21.. a .˜cos23..f .ˆ,...(23) 11 According to E. (26), we then obtain: Q ˜ °a ..s. 21co23ˆ..f ..,...(24) It can be seen from E. (24) that,theoretically, the upper bound of the length difference between any two branching chains of the 1T2R head is a function of . and .. 3.3 Simulation The form of motion in which the tool makes an arc trajectory parallel to the plane of the static platform is relatively simple and easy to describe. Therefore, under this motion, the above conclusions are verified by combining the existing 1T2R parallel mechanism data. hen the values of and . remain unchanged O' at . .0 360, the length of 1T2R head branch chain prismatic pair is calculated as shown in Fig. 5. Since the leg length data q 1, q 2 and q 3 satisfy the symmetric three-phase sine uantity. The values of O and . can be changed and the leg length data summed to obtain Table 2. The differenceof the leg length data is taken to obtain Table 3. Table 2. Summation of 1T2R head branches’ length Fig. 5. Data distribution type of 1T2R branches’ length . 3 . 2. . 1. . . . 0. max/min 1.8720 Q Q max min 1.879424 1.876457 1.874194 1.873482 1.872384 1.872296 1.872018 1.872016 1.872000 1.872000 1.879424 1.872000 1.9695 Q Q max min 1.976447 1.973780 1.971561 1.970920 1.969863 1.969783 1.969517 1.969515 1.969500 1.969500 1.976447 1.969500 2.0670 Q Q max min 2.073528 2.071116 2.068943 2.068363 2.067343 2.067271 2.067016 2.067015 2.067000 2.067000 2.073528 2.067000 2.1645 Q Q max min 2.170656 2.168465 2.166338 2.165811 2.164826 2.164760 2.164516 2.164514 2.164500 2.164500 2.170656 2.164500 2.2620 Q Q max min 2.267824 2.265824 2.263744 2.263262 2.262310 2.262250 2.262015 2.262013 2.262000 2.262000 2.267824 2.262000 Table 3. Subtraction of any two 1T2R head branches’ length O ' . 3 . 2. . 1. . . . 0. 0.6240 Q Q ' max ' min 0.271960 -0.271960 0.209725 -0.209725 0.140941 -0.140941 0.067737 -0.067737 0 0 0.6565 Q Q ' max ' min 0.272005 -0.272005 0.209745 -0.209745 0.140944 -0.140944 0.067737 -0.067737 0 0 0.6890 Q Q ' max ' min 0.272056 -0.272056 0.209763 -0.209763 0.140947 -0.140947 0.067737 -0.067737 0 0 0.7215 Q Q ' max ' min 0.272098 -0.272098 0.209779 -0.209779 0.140950 -0.140950 0.067737 -0.067737 0 0 0.7540 Q Q ' max ' min 0.272134 -0.272134 0.209793 -0.209793 0.140952 -0.140952 0.067737 -0.067737 0 0 Average value Q Q ' max ' min 0.272051 -0.272051 0.209761 -0.209761 0.140947 -0.140947 0.067737 -0.067737 0 0 According to the data of Fig. 5 and Table 2, the where, Q represents the lower bound of the length following relationship can be fit: sum of the branch chain movement, Q ˜3 zO , ° Q represents the upper bound of the length sum of Q ==QQ , (25) the branch chain movement Q ˜3z °0012251z .0015190.. . O .O . From the data in Table 2 and the E. (25), it can be seen that the length of the 1T2R head branch chain prismatic pair is always approximately eual to three times of , which satisfies the E. (21). Therefore, O the motion characteristic law of this form of motion can be extended. By the same token, through Fig. 5 and Table 3, Q' and . can be fitted as follows: Q ˜ °Q ˜, (26) where, Q ' represents the upper bound of the length difference Q' of the auxiliary leg of the branch chain movement, Q ˜°...0003748. 0006997. It can be seen from E. (26) that in the case of the same nutation angle, the upper bound of the difference between the length of any two kinematic chains is almost the same, with only an error of orders of magnitude. This satisfies E. (24). Therefore, the motion characteristic law of this form of motion can be extended. The above is the analysis of the motion characteristics of a 1T2R head and the simplified algorithm of pose judgment. The aim of the algorithm is to further divide the working space of the mechanism through the information of the three branch chains, so that the judgment of the attitude of the mechanism can be obtained without establishing a complex mapping model of the servo motor action. ithin its workspace, the sum length of the 1T2R head branch chain prismatic pair is always approximately eual to three times the coordinate value of the geometric centre of the moving platform. The absolute value of the difference between the lengths of any two branch chains is always less than a specific value related to the nutation angle. Thus, the limit protection of 1T2R head can be realized simply, reliably, and efficiently. 4 LIMIT PROTECTION IMPLEMENTATION PROCESS Based on the above analysis results of the movement characteristics of a 1T2R power head, the simplified algorithm of position and pose judgment, the limit protection method is designed in combination with the actual mechanical structure of the 1T2R power head servo feed system and its control system. This method will adopt two methods: proactive limit and preventive limit. The active limits are all areas of the 1T2R power head dexterous workspace except the boundary. The preventive limits are to prevent the mechanism from exceeding the workspace boundary, i.e., the boundary area of the 1T2R power head dexterous workspace. Real-time data of the prismatic pair length of the branch chain can be obtained through the feedback signal of the encoder of the servo motor. The PLC program in the motion controller PMAC is used to sum the leg lengths of the three branches and calculate the difference between the lengths of any two branches. The results of the operation are fed into the motion controller register. Then, according to the above pose judgment algorithm, whether the limit area is exceeded can be determined. The specific realization method of the 1T2R power head limit is shown in Fig. 6. Fig. 6. Schematic diagram of limit protection process For the active limit, according to Table 2, Q min 1.872000 is determined to be the lower trigger condition value of the active limit, and Q 2.267824 is determined to be the upper trigger max condition value of the active limit. For the preventive limit position, according to the data shown in Table 3, it is determined that Q' 1.272134 is the upper max limit trigger condition value of the preventive limit position, and Q' min –0.272134 is the lower limit trigger conditionvalue of the preventive limit position. 5 EXPERIMENTAL ANALYSIS The experimental platform is shown in Fig. 7, and the 1T2R parallel power head is driven by three chains. Fig. 8 shows the servo feed control system. The real-time length of three branches can be obtained through the code disk of the branch servo motor. The trajectory is an approximately circular trajectory parallel to the plane of the static platform, as shown in Fig. . According to the scale parameters (a, b , l , e) mentioned above, after consulting the machine tool operation manual, it is found that collision interference occurs easily when the nutation angle is 40, so the nutation angle . is selected as 30. The motion position of the tool reference point P at the end of the mechanism in the direction is 1.1687 m, the nutation angle . is 30, and the precession angle . is varied from 0 to 360. The length data of each chain obtained during the experiment are shown in Table 4. Fig. 10 shows the variation rule of the length of the branch chain corresponding to the trajectory. Fig. 9. Schematic diagram of trajectory reachable workspace Fig. 11a shows the numerical distribution law after the summation of the length data of the motion track. It can be seen from the distribution of the summation result that the sum of the length is distributed within the interval range of 2.2514 m to 2.2521 m and has periodic data fluctuation with a small amplitude. As shown in Fig. 11b, in comparing the length of the data points and the result data with the threshold value of the over-limit setting, it is found that none of the data points in the trajectory exceed the limit. t [s] q q q q q q q q q 0.0000 0.0221 0.0443 0.0664 0.0885 0.1107 0.1328 0.1549 1 [m] 0.6269 0.627 0.627 0.627 0.627 0.627 0.627 0.6271 2 [m] 0.8121 0.8121 0.8118 0.8114 0.811 0.8106 0.8101 0.8097 3 [m] 0.8123 0.8128 0.8132 0.8136 0.814 0.8144 0.8149 0.8153 t [s] 0.1771 0.1992 0.2214 0.2435 0.2656 0.2878 0.3099 … 1 [m] 0.6271 0.6271 0.6271 0.6271 0.6271 0.6272 0.6272 … 2 [m] 0.8093 0.8089 0.8084 0.8080 0.8076 0.8071 0.8067 … 3 [m] 0.8157 0.8161 0.8165 0.8169 0.8173 0.8177 0.8181 … t [s] 35.9922 36.0143 36.0365 36.0586 36.0807 36.1029 36.1250 36.1471 1 [m] 0.627 0.627 0.627 0.627 0.627 0.627 0.627 0.627 2 [m] 0.8129 0.8126 0.8125 0.8125 0.8125 0.8125 0.8125 0.8125 3 [m] 0.8121 0.8124 0.8125 0.8125 0.8125 0.8125 0.8125 0.8125 distribution form of the difference results of different branch leg lengths is the same except for the phase differences. As shown in Fig. 12b, by comparing the result data of the leg length difference of data points with the threshold value of the over-limit setting, it can be seen that none of the data points in the trajectory exceed the limit. Through the above simulation verification of the simplification algorithm of pose judgment, it is proved that the simplification algorithm of pose judgment can realize real-time judgment of the terminal pose state of the mechanism during the action process of a 1T2R power head. The result of pose judgment of this algorithm is found to be accurate. The active limit and preventive limit set according to the simplified algorithm of position and pose judgment can realize the limit protection function of 1T2R mechanism accurately and reliably and ensure the safe and reliable operation of the mechanism. a) numerical distribution of the summation, and b) extra-limit judgment of summation 6 CONCLUSION (1) This paper presents a simplified algorithm for the determination of the position and pose of 1T2R power head. Compared with the pose judgment method using the inverse solution model, it is estimated to be .42 faster with no need to upgrade the original system. The simplified algorithm is easy to implement, efficient and reliable. (2) The real-time limit protection of a 1T2R power head can be realized simply and reliably with the help of the simplification algorithm of position and pose judgment. Based on the software and hardware conditions of the existing numerical control system, this method can effectively solve the real-time limit protection problem of parallel machine tools by combining software and hardware and improve the operation safety of such topological machine tools. 7 ACKNOWLEDGEMENTS The authors would like to acknowledge the financial supported from the National Natural Science Foundation of China (Grant No.51575385) and National Key Research and Development Program of China (Grant No.201YFA070004). The authors also acknowledge Key Laboratory of Mechanism Theory and Euipment Design of Ministry of Education and School of Mechanical Engineering the in Tianjin University for providing the experimental environment. 8 REFERENCES [1] Guo, J.Z., Wang, D., Fan, R., Chen, W.Y. (2017). Design and workspace analysis of a 3-degree-of-freedom parallel swivel head with large tilting capacity. Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, vol. 231, no. 10, p. 1838-1849, DOI:10.1177/0954405415607782. [2] Ni, Y., Zhou, H., Shao, C., Li, J. (2019). Research on the error averaging effect in a rolling guide pair. Chinese Journal of Mechanical Engineering, vol. 32, art. ID. 72, DOI:10.1186/ s10033-019-0386-y. [3] Cheng, G., Gu, W., Yu, J.-l., Tang, P. (2011). Overall structure calibration of 3-ucr parallel manipulator based on quaternion method. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 10, p. 719-729, DOI:10.5545/sv-jme.2010.167. [4] Ni, Y., Zhang, Y., Sun, K., Wang, H., Sun, Y. (2016). Interpolation control algorithm for a three-RPS parallel spindle head. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems & Control Engineering, vol. 230, no. 7, p. 661-671, DOI:10.1177/0959651816645687. [5] Baizid, K., Cukovic, S., Iqbal, J., Yousnadj, A., Chellali, R., Meddahi, A., Devedžic, G., Ghionea, I. (2016). IRoSim: Industrial robotics simulation design planning and optimization platform based on CAD and knowledgeware technologies. Robotics and Computer-Integrated Manufacturing, vol. 42, p. 121-134, DOI:10.1016/j.rcim.2016.06.003. [6] Kumar, V. (1992). Characterization of workspaces of parallel manipulators. Journal of Mechanical Design, vol. 114, no. 3, p. 368-375, DOI:10.1115/1.2926562. [7] Xu, P., Cheung, C.-F., Li, B., Ho, L.-T., Zhang, J.-F. (2017). Kinematics analysis of a hybrid manipulator for computer controlled ultra-precision freeform polishing. Robotics and Computer-Integrated Manufacturing, vol. 44, p. 44-56, DOI:10.1016/j.rcim.2016.08.003. [8] Macho, E., Altuzarra, O., Amezua, E., Hernandez, A. (2009). Obtaining configuration space and singularity maps for parallel manipulators. Mechanism & Machine Theory, vol. 44, no. 11, p. 2110-2125, DOI:10.1016/j.mechmachtheory.2009.06.003. [9] Huang, C.K., Tsai, K.Y. (2015). A general method to determine compatible orientation workspaces for different types of 6-DOF parallel manipulators. Mechanism & Machine Theory, vol. 85, no. p. 129-146, DOI:10.1016/j.mechmachtheory.2014.11.011. [10] Zhang, D., Gao, Z. (2012). Optimal kinematic calibration of parallel manipulators with pseudoerror theory and cooperative coevolutionary network. IEEE Transactions on Industrial Electronics, vol. 59, no. 8, p. 3221-3231, DOI:10.1109/ TIE.2011.2166229. [11] Cui, G., Zhang, H., Zhang, D., Xu, F. (2015). Analysis of the kinematic accuracy reliability of a 3-DOF parallel robot manipulator. International Journal of Advanced Robotic Systems, vol. 12, no. 2, p. 15-26, DOI:10.5772/60056. [12] Shao, J., Chen, W., Fu, X. (2015). Position, singularity and workspace analysis of 3-PSR-O spatial parallel manipulator. Chinese Journal of Mechanical Engineering, vol. 28, p. 437­450, DOI:10.3901/CJME.2015.0122.018. [13] Kaloorazi, M.H.F., Masouleh, M.T., Caro, S. (2015). Determination of the maximal singularity-free workspace of 3-DOF parallel mechanisms with a constructive geometric approach. Mechanism & Machine Theory, vol. 84, p. 25-36, DOI:10.1016/j.mechmachtheory.2014.10.003. [14] Zhang, D., Wei, B. (2017). Interactions and optimizations analysis between stiffness and workspace of 3-UPU robotic mechanism. Measurement Science Review, vol. 17, no. 2, p. 83-92, DOI:10.1515/msr-2017-0011. [15] Gao, Z., Zhang, D. (2011). Workspace representation and optimization of a novel parallel mechanism with three-degrees-of-freedom. Sustainability, vol. 3, no. 11, p. 2217­2228, DOI:10.3390/su3112217. [16] Zhu, C., Guan, L., Han, J., Wang, L. (2009). A new resolution of workspace problem of parallel machine tool. IEEE International Conference on Automation and Logistics, p. 1002-1007, DOI:10.1109/ICAL.2009.5262566. [17] FarzanehKaloorazi, M.H., Masouleh, M.T., Caro, S. (2017). Collision-free workspace of parallel mechanisms based on an interval analysis approach. Robotica, vol. 35, no. 8, p. 1747­1760, DOI:10.1017/S0263574716000497. [18] Majid, M.Z.A., Huang, Z., Yao, Y.L. (2000). Workspace analysis of a six-degrees of freedom, three-prismatic-prismatic­spheric-revolute parallel manipulator. International Journal of Advanced Manufacturing Technology, vol. 16, p. 441-449, DOI:10.1007/s001700050176. [19] Rouhani, E., Nategh, M.J. (2015). Workspace, dexterity and dimensional optimization of microhexapod. Assembly Automation, vol. 35, no. 4, p. 341-347, DOI:10.1108/AA-03­2015-020. [20] Wan, J., Yao, J., Zhang, L., Wu, H. (2018). A weighted gradient projection method for inverse kinematics of redundant manipulators considering multiple performance criteria. Strojniški vestnik - Journal of Mechanical Engineering, vol. 64, no. 7-8, p. 475-487, DOI:10.5545/sv-jme.2017.5182. [21] Khoukhi, A., Baron, L., Balazinski, M. (2009). Constrained multi-objective trajectory planning of parallel kinematic machines. Robotics and Computer-Integrated Manufacturing, vol. 25, no. 4-5, p. 756-769, DOI:10.1016/j.rcim.2008.09.002. [22] Reveles, R.D., Pamanes, G.J.A., Wenger, P. (2016). Trajectory planning of kinematically redundant parallel manipulators by using multiple working modes. Mechanism and Machine Theory, vol. 98, p. 216-230, DOI:10.1016/j. mechmachtheory.2015.09.011. [23] Dash, A.K., Chen, I.-M., Yeo, S.H., Yang, G. (2005). Workspace generation and planning singularity-free path for parallel manipulators. Mechanism and Machine Theory, vol. 40, no. 7, p. 776-805, DOI:10.1016/j.mechmachtheory.2005.01.001. [24] Iqbal, J., Ullah, M.I., Khan, A.A., Irfan, M. (2015). Towards sophisticated control of robotic manipulators: An experimental study on a pseudo-industrial arm. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 7, p. 465-470, DOI:10.5545/sv-jme.2015.2511. [25] Manzoor, S., Islam, R.U., Khalid, A., Samad, A., Iqbal, J. (2014). An open-source multi-DOF articulated robotic educational platform for autonomous object manipulation. Robotics and Computer-Integrated Manufacturing, vol. 30, no. 3, p. 351­362, DOI:10.1016/j.rcim.2013.11.003. [26] Alam, W., Mehmood, A., Ali, K., Javaid, U., Alharbi, S., Iqbal, J. (2018). Nonlinear control of a flexible joint robotic manipulator with experimental validation. Strojniški vestnik - Journal of Mechanical Engineering, vol. 64, no. 1, p. 47-55, DOI:10.5545/sv-jme.2017.4786. [27] Mahmoodi, A., Sayadi, A.. (2015). Six-Dimensional space expression of workspace of six-DoF parallel manipulators using hyper spherical coordinates (HSC). Advanced Robotics, vol. 29, no. 23, p. 1527-1537, DOI:10.1080/01691864.2015 .1076345. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582 Received for review: 2022-07-08 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-09-02 DOI:10.5545/sv-jme.2022.274 Original Scientific Paper Accepted for publication: 2022-09-07 Nonlinear Vibration Analys is of Functionally G raded Porous Plates R einforced by G raphene Platelets on Nonlinear Elastic Foundations iaolin Huang – Chengzhe ang – Jiaheng ang – Nengguo ei Guilin University of Electronic Technology, School of Architecture and Transportation Engineering, China This paper presents a nonlinear vibration analysis of functionally graded graphene platelet (GPL) reinforced plates on nonlinear elastic foundations. Uniformly or non-uniformly distributed internal pores were present in the plates. Based on the modified Halpin-Tsai micromechanics model and the extended rule of mixture, the material properties were evaluated. The governing equations, coupled with the effect of nonlinear foundations, were derived by using the higher-order plate theory and general von Kármán-type equations. A two-step perturbation technique was employed to obtain the nonlinear frequency and transient response. After the present method was verified, the effects of pores, GPLs, and elastic foundations were investigated in detail. A new finding is that the influence of the porosity coefficient on the natural frequency and dynamic response is relevant to foundation parameters. Moreover, the influence of the nonlinear foundation parameter can be negligible. Keywords: functionally graded porous nanocomposites, graphene platelets, pores, nonlinear elastic foundation, nonlinear vibration, transient response Highlights • The material properties model of functionally graded (FG) graphene platelets reinforced plates is modified. • A two-step perturbation technique for the nonlinear vibration of porous plates on nonlinear elastic foundations is presented. • Some interesting conclusions about the effects of pores and nonlinear elastic foundations are drawn. 0 INTRODUCTION Due to the excellent load-carrying capacity with stronger bonding between the matrix and carbonaceous nanofillers [1] and [2], nano graphene platelets (GPL) are now increasingly used in many engineering fields, including aerospace, automobile, and civil engineering. It is necessary to study the dynamic characteristics of GPL-reinforced composite structures. In recent years, few results about the dynamic behaviour for functionally graded (FG) plates reinforced with GPL have been reported [3]. Song et al. [4] employed the first-order shear deformation plate theory to study the free and forced vibrations of FG multilayer GPL/polymer plates. Their results illustrated that a small amount of GPL can lead to higher natural freuency and lower dynamic response. By using the finite element method, hao et al. [5] studied investigated the bending and free vibration behaviours of FG trapezoidal plates. Their results also illustrated that the transient deflection was decreased by using a small number of graphene platelets. Gholami and Ansari [6] investigated the nonlinear harmonically excited vibration of FG graphene-reinforced composite plates. They found that the nonlinear buckling was increased with the rise of GPLs weight fraction. Moreover, u et al. [7] studied the dynamic stability of FG nanocomposite plates subjected to thermal and mechanical loads. Recently, a novel kind of graphene-reinforced nanocomposite was made by dispersing GPL into metal foams [8] to [10]. Due to the excellent mechanical and thermal properties, the novel composite materials have attracted some attention. Kittipornchai et al. [11] proposed a micromechanical model to estimate the typical properties of FG graphene reinforced nanocomposite, in which both GPLand internal pores are uniformly dispersed within each layer. Later, the model was employed to study the dynamic behaviour of graphene-reinforced nano-composite plates [12] and [13]. In these papers, the effect of internal pores has been discussed in detail, and the results showed that both pore volume fraction and the distribution of pores can affect the dynamic characteristics of porous structures. In practical engineering, plates are sometimes rested on elastic foundations. Thus, a suitable model is needed to evaluate foundation interaction. The simplest is the inkler model, in which a series of springs are used to calculate the tensile and compressed forces of elastic foundations. This model was laterdeveloped into the Pasternak model, in which a shear spring is added to estimate the shear forces of the inkler foundations. Thereafter, the Pasternak model has been widely used to study the static and dynamic behaviours of plates on elastic foundations [14] to [16]. To estimate the nonlinear interaction between foundations and plates, some researchers recently presented a nonlinear foundation model. For example, Nath et al. [17], Civalek [18] and Najafi et al. [19] studied the nonlinear dynamic behaviour of composite plates by using a nonlinear foundation model. Their results revealed that the nonlinear foundation parameter had a distinct influence on the dynamic characteristics. As mentioned above, the number of studies focused on the dynamic behaviour of FG porous plates reinforced with GPL is still rather scarce. According to the authors¶ knowledge, no previous work has been done to study the nonlinear dynamic behaviour of FG graphene platelets reinforced porous plates on nonlinear elastic foundations. Hence, this paper attempts to study the nonlinear vibration of the plates. Amodified material properties model is proposed, and the effects of internal pores and a nonlinear elastic foundation are discussed. 1 A POROUS NANO-COMPOSITE PLATE As depicted in Fig. 1, a functionally graded graphene-reinforced porous plate (length a, width b , thickness h) on a nonlinear elastic foundation is taken into account. The origin of the coordinate system (X , Y , Z ) is located at one corner of the middle plane of the plate. The -axis is perpendicular to the X -Y plane and points upwards. Three types of porosity distributions are considered (Fig. 2). P-1 indicates that the largest-size pores are distributed in the middle. In contrast, P-2 indicates that the largest-size pores are on the bottom and top surfaces. The symbol of even distribution is P-3. Unlike the material properties model presented by Kittipornchai et al. [11], the present model is based on the volume fraction of pores, which are assumed to be: e 0cosZh P °,1° ˜/˜.  * .˜Zh P ,2° VZ ˜°.e .1cos/°..(1) P0ˆ.˜.  1.,˜P .3°  * In the above euation, e0, e0(0 e0(e0*) 1 and a denote the porosity coefficients for P– 1, P– 2, and P– 3 distributions, respectively. Yangs elastic modulus E (Z ) and shear elastic modulus G (Z ) of the porous plate can be expressed by using the rule of mixture: EZ ˜E 0(°VP , ()1)GZ ˜G 0(°VP ,(2) ()1) where E 0 and G 0 are the corresponding variations of the graphene-reinforced nanocomposites without internal pores. Fig. 1. A functionally graded GPL-reinforced porous plate on a nonlinear elastic foundation a) b) c) Fig. 2. Three types of porosity distributions; a) P –1, b) P –2, and c) P –3 Because the mass density .(Z ) cannot be calculated by using the rule of mixture [11], it is assumed to be Based on the rule of mixture, Poissons ratio ., ˜0..1.e 1cos.°Z /h ˆ,.P .1ˆ and mass density .0 can be calculated as ˜().˜.1.e (1.cos°Z /h , Z 0.1*.ˆ.P .2ˆ ,(3) ˜().˜V .˜1.V ), Z ( GPLGPL m GPL ˜,.P .3ˆ 0'°.°V .°(1.V ),(8) 0GPLGPL m GPL where .is the mass density of GPL, and .m and .m GPL where .0 is the mass density of the nanocomposites *are the mass density and Poissons ratio of the matrix. without pores. Coefficients e1, e1and a' are the mass The shear modulus G 0 can be obtained by density coefficients for the three types of porosity distributions. E 0G ˜.() 0 For open-cell metal foams, the relationship 21.. .°0 between Yangs elastic modulus and mass density is Ez ()°.Z 2 ().* [12]. Hence, the coefficients e1, e ˜..1 E . 0.0ˆ and a' can be estimated as 1   Z  Z ° ..ˆ°..ˆ ,  P ˜ 1,  ˜ . ˜  e 1 e 0 1 .. .. cos cos h h 1 Fig. 3. Three types of GPL dispersion patterns; a) G -1, b) G -2, and c) G -3  Z  Z ° ..ˆ°..ˆ ° ..ˆ°..ˆ  ˜ 2  ,(4) P ˜˜ . ˜˜ * * e 1 e 0 1  .. 1 .. 1 .. .. cos cos h h , P ˜ 3  .  As shown in Fig. 3, three types of GPLdispersion   patterns (G-1, G-2, G-3) are considered. V is GPL expressed as It is assumed that the masses of all plates with .vi [1°cos(Z /h )].G °1. GPL different porosity distributions are euivalent. Hence, 1 . . .. V cos(Z /)]h G ˜ ° (10) z () v i 2 e0* and a can be calculated by , 1 . v . G °3. ˆi 3 1 1  h h /2  /2  Z Z °..ˆ°..ˆ where .i 1, .i 2 and .i 3 are the maximum value of   dZ dZ ˜˜ .˜ *  e 0 e 0 (1 .... coscos h h V GPL, i (i 1, 2, 3) indicate the three types of porosity distributions. .i 1, .i 2 and .i 3 can be reckoned as h /2 h /2 ˜ ˜ , (5) 1 h /2h /2  Z  °..ˆ   . ˜ dZ e 0 .. cos follows: h h /2 ˜ h ˜ /2 Based on the Halpin-Tsai micro-mechanical VT ˜WGPL .m , GPL (11) mode, the effective Youngs elastic modulus E 0can be WGPL .m °.GPL .WGPL .GPL expressed by ˜ Z () ° Z H  1 Z . dZ ..ˆ . h /2 . cosvi  .1°V ..1°V . LLGPL WWGPL 1 ˜ *h /2 3 5 E m , (6) E E ˜ ° ˆ ˆ  m 0 V LGPL . V WGPL . ˜ Z () ˜° Z ..ˆ. h h 81 81 . . . . () /2/2 VGPL T .  1 .  dZ dZ ,(12) . v i cos ˜˜ H ** h /2 h /2 where E and V GPL are Youngs elastic modulus of m matrix and GPL volume fraction, respectively. The coefficients ., ., . and . are defined by LWLW ˜ () Z h /2 .  dZ v i 3 ˜ *h /2 T in which VGPL and W GPL are the total volume and 2a 2b GPL GPL ˜L .h ,˜L .h ,weight fractions of GPL, respectively. GPL GPL E /E .1E /E .1 GPL m GPL m 2 FORMULATIONS ˜L .,˜W .,(7) E /E .°E /E .° GPL mL GPL mW 2.1 Governing Equations in which E , a, b and hare Youngs GPL GPL GPL GPL modulus, average length, width, and thickness of As noted by Civalek [18], the effect of nonlinear plate-graphene platelets, respectively. foundation interaction on the dynamic response of plates on elastic foundations must not be neglected. Therefore, the following three-parameter nonlinear foundation model is adopted: 444 ..W .W .W .3 R KWK °..2.,(13) ˜.KW f1242243 .X XY .Y .. ˆ where K , K , and K are inkler, Pasternak, and 12 3 nonlinear foundation parameters. According to Reddys higher-order thick plate theory [20], the displacements of the thick composite plate are assumed to be E .E .. . Q ˜Q ˜2,Q ˜2, 112212 1°.1°. E .. ... Q ˜Q ˜0,Q ˜Q ˜Q ˜.(17) 1626445566 2(1 ..) The in-plane forces Ni , bending moment Mi , higher-order bending moment Pi , shear forces Qi , and higher-order shear forces Ri are h /2 ˜NMP °..ˆ˜1 ZZ 3 °dZ ,( i .126, , , ,, ,,) iii i .h /2 h /2 2 ˜QR , °..ˆ˜1, Z ° dZ , 22 4 .h /2   2  .ˆ. .. Z h W X W Y  .ˆ. .. 4 h /2 ˜QR °.˜, ° , ˆ12 11 5 . . UXYt (,,) Z  ˜ ° . °  Z d Z . (18) u 1 , 11 3 h /2  . By using the Hamilton principle, the euations of .. 2 3   2  .ˆ. .. Z h .ˆ. .. 4 VXY (, , t Z  ° . ° ˜  u ) , motion for the plate can be derived as y y 2 3  N ˜ N ˜ W ˜ .. IU (14) 4 WXYt (,,) ˜ u I I ° . ° . 6 Y ˜ X ˜ 1 , 3 , 114 h h 43 X ˜ 2 in which U , V , W , .1 and . 2 are the displacements N ˜ N ˜ W ˜ .. IV . .. 2 I I 6 ° . ° . 2 and rotations of a point (X , Y ) on the mid-plane. , 124 Y ˜ ˜ X Q Y ˜ ˜˜˜ X 2 The von Karman strains associated with the displacement field in E. (14) can be stated as ..°.32U 1W 4Z W 11.Z °.....˜ .,2Y ..2W ,2Y  , .ˆ..X ...ˆ..32Y ..W .2Y 2X 3h . .4Z °2.Z .2Y h 3 °ˆ...2.W 4Z .2Y h 2X .ˆ.. 2W Y ..ˆ.. °. ˜ 42 1X . V 12Y 2....˜ .˜ 0,3 ˜˜W X ˜W N 1 ...˜˜X W ˆ..˜ ˜ W Y ° ° N ° Q 1 26 ...˜ Y ˜ 4 2 ˆ.. ... ˜ R 1X ˜˜ R 2Y ˆ.. ˜˜ ˜˜ Y N ° ° N  ° 62 ˜ Y 222 .PP ˆ 162 X h 2 ° P ° ° 4 R ˜ ˜ ˜ ˜ . . 2 f °° YY 3h °X 2 °Y 22 . . 2.. 2.. . 2°W °W .. ˆ .. .4ˆ qI1WI 7.2.3. ˜ ˜ . .. . . , h °X 2°Y 2  , W . 4 Z W . 2 . . ˜  . . ˆ  . 1 51 h X . X . 2 .ˆ. .. M ˜ M ˜ P 1 ˜ P 6 ˜ 4 h 43 2 I Q1 .. W R 1˜ ° . ° . ° 61 ° .° . W .. Z 2 . Y W .. X W .. . VX .. .. hY U .. 3 X ˜  Y ˜ h Y X ˜ ˜ 2 ˜ . . 12 6 X . Y . .. IU h 4 3  M ˜ .. 3 ..  . °° . Z . 3 I  ° . 4 , ˆ . . (15) 1 2 x 25 X ˜ 2 2 X . Y . .. XY 2 .ˆ. .. M ˜ P 6 ˜ ˜ Y X P 2 ˜ 4 h 43 ° I Q2 .. W R 2 ˜ . ° . ° 61 Based on Hooks law, the relationship between Y ˜ X ˜ h ˜ 22 stresses and strains can be expressed as .. IV 4 h  .. 3 3 , In E. (1), the constants Ij and Ij were given by Reddy [20]. The superposed dots indicate the .. I (1)  ° . .ˆ Q Q ˜˜˜ ..... ° 0 . ˆ .ˆ y 25 Y ˜ 21112 .......... 1 1 .... ° ........ .. QQ 0  2 12222 . ° 00Q 6 6 66 differentiation with respect to time. ˆ..  .... . Q 44 0 . 0Q 55 ˆ...... ° ° ˜˜ ... The strain compatibility euation is 4 5 ˆ.. 4 ,(16) 5 202020222 ˜˜˜.˜W .2˜W ˜W ° . . . ..22 ˜˜ (20) 612 2 in which ˜Y ˜X ˜˜ XY ˜˜ XY X Y 2 ˆ  The in-plane forces Ni can be expressed by stress function FXYt (,,): 222 °F °F °F N ˜2,N ˜2,N ˜..(21) 126 °° °Y °X XY By substituting Es. (18) and (21) into Es. (1) and (20), the governing euations of nonlinear vibration for the plate can be derived as follows: ().()lF R lW ˜l ()˜l .°()° 1112113214f .. ... ... .. 12 .lWFlW I (,)°()°° I °q , 1788 .X .Y lF °l .°l .˜lW .˜1(, ()()()()lWW ), 2122123224 2 .. .W .. l (W l )°()˜l ()°lF I ..().°I ., 31321332349101 .X .. ..().lW ()˜l ()° l ()°lF I .W °I ... ,(22) 41421432449102 .Y where the constants Ij ( j 8, , 10), linear operators lij and nonlinear operator l were given by Shen [21] and Huang and heng [22]. The four edges of the plate are assumed to be simply supported. The boundary conditions are expressed as X ˜0,a:W ˜°˜M ˜P ˜N ˜0, 1116 Y ˜0,b:W ˜°˜M ˜P ˜N ˜0.(23) 2226 2.2 Solution Procedure To solve the nonlinear euations, Es. (22) and (23), we first introduce the following dimensionless parameters: °X °YZa W x ˜,y ˜,z ˜,.˜,W ˜/, ****14 a bhb[DD AA ] 11221122 F (,..)a 12 F ˜,(,)˜, .. ** 12****/ [DD ]/12°[DD AA ]14 112211221122ˆa 4a 2 , .424 **** °°Ka DDAA .311221122 (,KK )˜,K ˜, 123 *4* (,)° DKK D 111211 qa 4°t E .˜,˜ m ,(24) q 4*****14 DDDAA ]/a  °[ 1111221122m where the stiffness constants A ij, Bi j, and D ij are defined in the standard way [20]. The dimensionless form of the euation (22) can be rewritten as ()˜lF KWK W .2 lW l (.)˜l (.)°ˆ()°˜ 1112113211412 32... 1... 2 .. °KW .ˆ.lWF °l()°ˆ°ˆ.° (,)lW , 311722q .x .y 12 ()ˆ.(,), lF ()°ˆl (.)°ˆl (.)˜ˆlW .˜lWW 2122212232 24242 2 .. ()°()˜l ()°ˆ() .ˆ.W °ˆ,(25) lW l lF.. 1 3132133213445 .x where the constants .(i= 1, 2, ..., 5), the dimensionless operators l ij and l can be seen in the previous work [22]. The dimensionless form of the boundary conditions in E. (23) are also rewritten as: i x˜0,°:W ˜.˜M ˜P ˜N ˜0, 1116 y ˜0,°:W ˜.˜M ˜P ˜N ˜0.(26) 2226 Atwo-step perturbation techniue is employed to solve the nonlinear governing euations, E. (25). As the essence of this procedure, the asymptotic solution is supposed to be ˜i ˜ Wx(,,,y ˜°)..°wi (xy ,,),˜ i .1 ˜i ˜ (,,,).i (,,), Fxy ˜°.°f xy ˜ i .1°i ˜ . (,,,xy ˜°)..°.(,,), xy ˜ x 1i i .1 .y (,,,xy ˜˜°)..°i .2i (xy ,,)˜˜, i .1 ˜i ˜ ˆq (,,,xy ˜°)..°ˆi (xy ,,).(27) ˜ i .0 In E. (27), the time parameter t. (t. et) is used to improve the perturbation procedure. Substituting E. (27) into Es. (25) and (26), then solving the perturbation euations step by step, the displacements W , ., .y and stress function F can be obtained. The x dimensionless transverse load . can be derived as q ˜xy °.[gw °.g .w .. °mx sin (,,).()()]sinny q 411431 .(.w ())(°2g cocs2mx .g cos2ny ) 1441442 333 42[w 1()](sinmx sin)ny .o (). .g .°.(28) In E. (28), it should be noted that t. is replaced by t. Multiplying E. (28) by (sin mx sin ny ) and integrating over the plate area, the following nonlinear ordinary differential euation can be obtained: 2 d (˜w 1)23. g .g (w ).g (˜w ).g (˜w )..q (),(2) ˜° 1213141 d °2 in which . .. ˜°().4..˜q (,,)xy °sinmx sinny dxdy .(30) 00 q .2 The nonlinear ordinary euation, E. (2) can be Exam ple 2. The dynamic response of a FG solved by using the Runge-Kutta iteration Scheme GPL reinforced plate is investigated in this example. [23]. For the free vibration problem (.q (t) 0), the The dispersion pattern of GPLis G – 2. The plate is approximate nonlinear freuency can be derived as subjected to the explore load: / 2 .9gg 10g .12 , (31) qXYt ,) (, ˆ... ˜ p p , (32) ˜NL ° A . 2 .  Pm ° .t .t . t / t g 41 (1 ), 0 424144 g 12gg ˆ434143 where AW h is the dimensionless vibration = / amplitude. If Amax0, the dimensionless natural freuency is ˜L °g /g . 4143 3 RESULTS AND DISCUSSION In the section, the several dimensionless parameters are used as follows: 4242 Ka Ka Kah ww 3 k ˜,k ˜,k ˜, 123 DD D mm m Eh 3.a 2 . . mLm m D ˜,ˆ˜ ,ˆ. ˜.h . 2 m 121(°.)hE LE m mm . tt where the peak pulse pm is 500 kPa and the loading time tp is 0.01 s. The material properties of the matrix are E 3.0 GPa, .m 120 kg/m and .m 0.34. m The corresponding parameters of GPLs are E G PL PP 1.01 TPa, .GL 1062.5 kg/m and .GL 0.186. The geometric parameters of GPLs and the plate are a b 0.45, h 0.045 m, aG PL 2.5 m, bGPL 1.5 m and h1.5 nm. The curves of central transient GPL deflection versus time are depicted in Fig. 4. It is found that the present results are agreement with those given by Song et al. [4]. 0.06 . t 0, 0 and p 3.1 Comparison Studies To verify the accuracy and effectiveness of the present method, two numerical examples are presented in this subsection. Exam ple 1. The dimensionless fundamental freuencies of GPL-reinforced porous plates resting on elastic foundations are calculated and listed in Table 1. The material properties and dimensions of GPL are EP 1.01 TPa, .P 1062.5 kg/m3, GL GL .G PL 0.186, aG PL 2.5 m, bGPL 1.5 m, hG PL 1.5 nm. The material properties of the Matrix are E m 200 GPa, .m 808 kg/m3, .m 0.31. The geometrical parameters of the plate h 0.05 m, a b 1.0 m. The GPL weight fraction and porosity coefficient are W 5 , e0 0.4. It can GPL be observed that the present results are close to those given by Gao et al. [13]. The maximum error is about 2.3 . This is because Gao et al. [13] employed the classic plate theory and differential uadrature method to calculate the fundamental freuencies, which is different from the present method. t ms@ Fig. 4. Comparison of the central transient deflection for an FG GPL-reinforced plate 3.2 Parametric Studies In what follows, the effects of material properties and foundation parameters are investigated in detail. The material and geometric parameters are the same as those in Example 1. Unless specially stated, the weight fraction W and pore coefficient e0 are 5 GPL and 0.2. The dimensionless foundation parameters (k 1, k 2, k 3) are (50, 50, 50). Table 1. Comparison of dimensionless fundamental frequencies .. for GPL reinforced plates resting on Winkler-Pasternak elastic foundations (k 1, k 2) (100, 100) Method Ref. [13] present Error [%] P –1, G –1 P –2, G –1 0.0637 0.0591 0.0652 0.0603 2.3 2.0 P –1, G –3 0.0745 0.0761 2.1 P –2, G –3 0.0696 0.0710 2.0 576 Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G. Tables 2 to 4 list the fundamental natural freuencies of the plate with different GPLdispersion pattern, porosity distribution, weight fractions W , GPL porosity coefficients e0 and foundation parameters (k 1, k 2). It can be seen that the natural freuency is increased with the rising weight fraction W and GPL foundation parameters (k 1, k 2). The freuency for G–1 is higher than those for G–2 and G–3. This illustrates that G–1 can strengthen the plate more effectively than the other two patterns. Also, it can be seen that the freuency for P–1 is lower than those for P–2 and P–3. The fact demonstrated that P–1 can weaken the plate more seriously. In past studies [11] to [13], a conclusion was drawn that the natural freuency monotonously decreases with the increasing porosity coefficient. However, in the present case, the conclusion was correct only on P–2. On P–1 and P–3, the effect of the porosity coefficient Table 2. Dimensionless fundamental frequencies for porosity distribution P –1 (k 1, k 2) =(0, 0) (k 1, k 2) =(50, 0) (k 1, k 2) =(50, 50) GPL e0 [wt%] [wt%] [wt%] 2.0 5.0 8.0 2.0 5.0 8.0 2.0 5.0 8.0 0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 G –1 0.2 8.970 12.149 14.623 9.280 12.422 14.889 14.055 16.952 19.400 0.4 8.928 12.076 14.530 9.265 12.376 14.820 14.378 17.247 19.683 0.0 7.359 9.072 10.500 7.709 9.411 10.842 12.783 14.580 16.172 G –2 0.2 7.311 8.997 10.404 7.687 9.363 10.773 13.056 14.848 16.440 0.4 7.283 8.946 10.337 7.672 9.345 10.739 13.415 15.212 16.812 0.0 8.014 10.364 12.260 8.336 10.663 12.554 13.172 15.420 17.371 G –3 0.2 7.982 10.323 12.210 8.328 10.643 12.526 13.444 15.690 17.643 0.4 7.973 10.310 12.198 8.348 10.559 12.440 13.803 16.056 18.021 Table 3. Dimensionless fundamental frequencies for porosity distribution P –2 (k 1, k 2) =(0, 0) (k 1, k 2) =(50, 0) (k 1, k 2) =(50, 50) GPL e0 [wt%] [wt%] [wt%] 2.0 5.0 8.0 2.0 5.0 8.0 2.0 5.0 8.0 0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 G –1 0.2 8.283 11.328 13.526 8.617 11.534 13.822 13.629 16.313 18.595 0.4 7.397 10.035 12.088 7.800 10.394 12.435 13.484 15.889 17.959 0.0 7.359 9.072 10.500 7.709 9.411 10.842 12.783 14.580 16.172 G –2 0.2 6.775 8.382 9.719 7.180 8.774 10.112 12.765 14.486 16.017 0.4 6.110 7.606 8.845 6.593 8.072 9.312 12.720 14.468 15.942 0.0 8.014 10.364 12.260 8.336 10.663 12.554 13.172 15.420 17.371 G –3 0.2 7.337 9.488 11.223 7.712 9.836 11.566 13.072 15.155 16.977 0.4 6.553 8.473 10.024 7.004 8.894 10.438 13.037 14.945 16.630 Table 4. Dimensionless fundamental frequencies for porosity distribution P 3 (k 1, k 2) =(0, 0) (k 1, k 2) =(50, 0) (k 1, k 2) =(50, 50) GPL e0 [wt%] [wt%] [wt%] 2.0 5.0 8.0 2.0 5.0 8.0 2.0 5.0 8.0 0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 G –1 0.2 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 0.4 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 G –2 0.2 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 0.4 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202 G –3 0.2 8.721 11.821 14.233 9.039 12.103 14.506 13.898 16.719 19.317 0.4 8.372 11.348 13.664 8.731 11.666 13.972 14.041 16.747 19.354 1.25 1.20 1.20 1.15 1.15 1.05 1.05 1.00 1.00 W /h maxW /h a) a) max 1.25 1.25 //....NL LNL L 1.10 /.. NL L 1.10 0.0 0.2 0.4 0.6 0.8 1.0 1.05 1.05 1.00 1.00 W /h b) maxW /h b) max Fig. 5. Effect of GPL dispersion pattern on the frequency ratio; Fig. 7. Effect of GPL weight fractions on the frequency ratio; a) P 1,and b) P 2 a) G 1/P-3, and b) G 2/P 3 1.20 1.25 1.20 1.20 1.15 /..NL L 1.15 1.10 1.10 0.0 0.2 0.4 0.6 0.8 1.0 1.20 1.15 /.. NL L 1.15 //.... NL LNL L 1.10 1.10 1.05 1.05 1.00 1.00 W /h W /h max a) a) max 1.25 1.25 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.20 1.15 1.10 1.20 /.. NL L 1.15 1.10 1.05 1.05 1.00 1.00 W /h W /h maxmax b) b) Fig. 6. Effect of porosity distribution on the frequency ratio, Fig. 8. Effect of porosity coefficient on the frequency ratio; a) G 1, and b) G 2 a) G 3/P 1, and b) G 3/P 2 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 is relevant to the values of foundation parameters (k 1, 6 k 2). If foundation parameters are (50, 0) or (0, 0), the G-1 G-2 natural freuency is decreased. In contrast, the natural G-3 freuency for (50, 50) is increased. 3 [kN.m/m@ W mm@ Figs. 5 to 8 reveal the influences of GPL dispersion pattern, porosity distribution, GPL weight fraction W GPL , and pore coefficient e0on the nonlinear 0 to linear freuency ratio ./ .. As can be observed, NL L the freuency ratio for G – 2 / P– 2 is higher than those for other GPL and pore distributions. The freuency -3 ratio is decreased as pore coefficient e0 rises but rose 0 5 10 15 t ms@ a) with the increase of GPL fraction W . GPL 80 G-1 1.5 60 G-2 k 1,k 2,k 3) (0,0,0) (((( G-3 1.4 40 20 k 1,k 2,k 3) (50,0,0) k 1,k 2,k 3) (50,50,0) k 1,k 2,k 3) (50,50,50) //.... NL LNL L 1.3 x M 0 1.2 -20 1.1 -40 0 5 10 15 1.0 t ms@ b) 0.0 0.2 0.4 0.6 0.8 1.0 Fig. 11. Effect of GPL dispersion patterns on the central transient W /h maxresponses for P 3; a) dynamic deflection, and Fig. 9. Effect of elastic foundation on the frequency ratio b) dynamic bending moment for G 3/P 3 5.0 P-1 1.25 P-2 P-3 k 30 1.20 2.5 k 3 50 k 3 100 W mm@ 1.15 0.0 1.10 1.05 -2.5 0 5 10 15 1.00 t ms@ a) 0.0 0.2 0.4 0.6 0.8 1.0 80 W /h P-1 P-2 max Fig. 10. Effect foundation parameter k3 on the frequency ratio P-3 [kN.m/m@ for G 3/P 3 The effect of nonlinear elastic foundation 40 x parameters on freuency ratio is shown in Figs. and 0 M 10. The two figures demonstrated that the freuency ratio reduces with the increasing parameters k 1 and k 2. However, the ratio rises as the nonlinear parameter -40 k 3 increases. 0 5 10 15 The curves of central transient deflections and b) t ms@ bending moments for different distributions of GPLs Fig. 12. Effect of porosity distributions on the transient responses and pores are depicted in Figs. 11 and 12. It is found for G 3; a) dynamic deflection, and b) dynamic bending moment Nonlinear Vibration Analysis of Functionally Graded Porous Plates Reinforced by Graphene Platelets on Nonlinear Elastic Foundations 579 that the maximum dynamic deflection for G – 2 / P– 2 4 is the largest among all patterns and distributions. e00 e0 0.2 e0 0.4 In contrast, the amplitude of the dynamic bending moment for G – 2 / P– 2 is the smallest. 2 The effects of GPL weight fraction W and GPL pore coefficient e0 on the transient responses is illustrated in Figs. 13 and 14. It is discerned that W mm@ 0 the rise of W reduces the amplitude of transient GPL deflection, but increases the amplitude of bending moment. The maximum deflection increases by about 8 as the porosity coefficient e0rises from 0.0 to 0.4. -2 0 5 10 15 Therefore, a conclusion may be made that the effect of a) t ms@ the porosity coefficient on the dynamic response can 80 be negligible. e00 e0 0.2 e0 40 0.4 W GPL 0 W GPL 2 3 W GPL 5 [kN.m/m@ x M [kN.m/m@ W mm@ x 0 M 0 -40 0 5 10 15 t ms@ b) -3 Fig. 14. Effect of porosity coefficient on the transient responses for 0.000 0.005 0.010 0.015 t ms@ P 1 / G 1; a) dynamic deflection, and b) dynamic bending moment a) 80 W GPL 0 8 (((( k 1,k 2,k 3) (0,0,0) 2 5 4 W GPL k 1,k 2,k 3) (50,0,0) W GPL k 1,k 2,k 3) (50,50,0) 40 k 1,k 2,k 3) (50,50,50) mm@ 0 W 0 -40 0 5 10 15 t ms@ b) Fig. 13. Effect of GPL weight fraction on the transient responses a) -4 0 200 5 t ms@ 10 15 (((( k 1,k 2,k 3) (0,0,0) 3) (50,50,0) for P 1 / G 1; a) dynamic deflection, and b) dynamic bending moment 150 k ,k ,k ) (50,0,0) 123 k k k 1,2, k 1,k 2,k 3) (50,50,50) 100 The effect of foundation parameters (k 1, k 2, k 3) on dynamic responses is presented in Fig. 15. As expected, inkler and Pasternak elastic foundation parameters reduce the dynamic responses. The dynamic response for Pasternak elastic foundations M [kN.m/m@ x 50 0 -50 (50, 50, 0) is very close to those for the nonlinear elastic foundations (50, 50, 50). Hence, a conclusion -100 0 5 10 15 may be drawn that the effect of nonlinear foundation t ms@ b) parameter k 3 on dynamic responses may be neglected, Fig. 15. Effect of elastic foundation on the transient responses for which is different from the statement mentioned by P 1 / G 1; a) dynamic deflection, and b) dynamic bending moment 580 Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G. Civalek [18] and Najafi et al. [19]. They deemed that the nonlinear foundation parameter k 3 had a significant effect on the dynamic responses of laminated and FGM plates. 4 CONCLUSIONS The present work presents a reliable and effective method to investigate the nonlinear free and forced vibrations of functionally graded GPL reinforced porous plates on nonlinear elastic foundations. The effects of internal pores, GPLs, and nonlinear elastic foundations are discussed in detail. Some interesting conclusions can be drawn from the numerical results: 1. Both GPL dispersion patterns and porosity distribution can affect the nonlinear vibrations and responses for porous plates. Furthermore, the effect of a GPL dispersion pattern is more significant than that of a porosity distribution. 2. The GPL weight fraction and foundation parameters k 1 and k 2 increase the natural freuency but decrease the nonlinear to linear freuency ratio and transient deflection. 3. The increase of porosity coefficient does not always lead to the rise of natural freuency and transient responses. 4. Nonlinear foundation parameters have insignificant effects on the nonlinear to linear freuency ratio and transient response. 5 ACKNOWLEDGEMENTS The authors are thankful for the financial support of the National Natural Science Foundation of China 12162010 and the Natural Science Foundation of Guangxi 2021GNSFAA 220087. 6 REFERENCES [1] Yusufoglu, E., Avey, M. (2021). Nonlinear dynamic behavior of hyperbolic paraboloidal shells reinforced by carbon nanotubes with various distributions. Journal of Applied and Computational Mechanics, vol. 7, p. 913-921, DOI:10.22055/ JACM.2021.36043.2783. [2] Avey, M., Yusufoglu, E. (2020). On the solution of large-amplitude vibration of carbon nanotube-based double-curved shallow shells. Mathematical Methods in the Applied Sciences, vol. 8, p. 1-13, DOI:10.1002/mma.6820. [3] Zhao, S.Y., Zhao, Z., Yang, Z.C., Ke, L.L., Kitipornchai, S., Yang, Y. (2020). Functionally graded graphene reinforced composite structures: A review. Engineering Structures, vol. 210, art. ID 110339, DOI:10.1016/j.engstruct.2020.110339. [4] Song, M., Kittipornchai S., Yang, J., (2017). Free and forced vibrations of functionally graded polymer composite plates reinforced with graphene nanoplatelets. Composite Structures, vol. 159, p. 579-588, DOI:10.1016/j.compstruct.2016.09.070. [5] Zhao, Z., Feng, C., Wang, Y., Yang, J. (2017). Bending and vibration of functionally graded trapezoidal nanocomposite plates reinforced with graphene nanoplatelets (GPLs). Composite Structures, vol. 180, p. 799-808, DOI:10.1016/j. compstruct.2017.08.044. [6] Gholami R., Ansari R. (2018). Nonlinear harmonically excited vibration of third-order shear deformable functionally graded graphene platelet- reinforced composite rectangular plates. Engineering and Structures, vol. 156, p. 197-209, DOI:10.1016/j.engstruct.2017.11.019. [7] Wu, H., Yang, J., Kittipornchai, S. (2018). Parametric instability of thermo-mechanically loaded functionally graded graphene reinforced nanaocomposite plates. International Journal of Mechanical Science, vol. 135, p. 431-440, DOI:10.1016/j. ijmecsci.2017.11.039. [8] Bartolucci, S.F., Paras, J., Rafiee, M.A., Rafiee, J., Lee, S., Kapoor, D. (2011). Graphene-aluminum Nanocomposites. Material Science and Engineering: A, vol. 528, no. 27, p. 7933-7937, DOI:10.1016/j.msea.2011.07.043. [9] Rashad, M., Pan, F., Tang, A., Asif, M. (2014). Effect of graphene nanoplatelets addition on mechanical properties of pure aluminum using a semi-powder method. Progress of Natural Science: Materials International, vol. 24, p. 101-108, DOI:10.1016/j.pnsc.2014.03.012. [10] Wu, H.L, Yang, J., Kitipornchai, S. (2020). Mechanical analysis of functionally graded porous structures: A review. International Journal of Structural Stability and Dynamics, vol. 20, no. 13, art. ID 2041015, DOI:10.1142/S0219455420410151. [11] Kittipornchai, S., Chen, D., Yang, J. (2017). Free vibration and elastic buckling of functionally graded porous beams reinforced with graphene platelets. Material Design, vol. 116, p. 656-665, DOI:10.1016/j.matdes.2016.12.061. [12] Yang, J., Chen, D., Kittipornchai, S. (2018). Buckling and free vibration analysis of functionally graded graphene reinforced porous nanocomposite plates on Chebyshev-Ritz method. Composite Structures, vol. 193, p. 281-294, DOI:10.1016/j. compstruct.2018.03.090. [13] Gao, K., Gao, W., Chen, D., Yang, J. (2018). Nonlinear vibration of functionally graded graphene platelets reinforced porous nanocomposite plates resting on elastic foundation. Composite Structures, vol. 204, p. 831-846, DOI:10.1016/j. compstruct.2018.08.013. [14] Saidi, A.R., Bahaadini, R., Majidi-Mozafari, K. (2019). On vibration and stability of porous plates reinforced by graphene platelets under aerodynamic loading. Composite Part B: Engineering, vol. 163, p. 778-799, DOI:10.1016/j. compositesb.2019.01.074. [15] Duc, N.D., Lee, J., Nguyen-Thoi, T., Thang, P.T. (2017). Static response and vibration of functionally graded carbon nanotube-reinforced composite rectangular plates resting on Winkler-Pasternak elastic foundations. Aerospace Science and Technology, vol. 68, p. 391-402, DOI:10.1016/j. ast.2017.05.032. [16] Haciyev, V.C., Sofiyev, A.H., Kuruoglu, N. (2018). Free bending vibration analysis of thin bidirectionally exponentially graded orthortropic rectangular plates resting on two-parameter elastic foundations. Composite Structures, vol. 184, p. 372­377, DOI:10.1016/j.compstruct.2017.10.014. [17] Nath, Y., Prithviraju, M., Mufti, A.A. (2006). Nonlinear statics and dynamics of antisymmetric composite laminated square plates supported on nonlinear elastic subgrade. Communication in Nonlinear Science and Numerical Simulation, vol. 11, p. 340-354, DOI:10.1016/j. cnsns.2004.11.003. [18] Civalek, Ö. (2013). Nonlinear dynamic response of laminated plates resting on nonlinear elastic foundations by the discrete singular convolution-differential quadrature coupled approaches. Composite Part B: Engineering, vol. 50, p. 171­179, DOI:10.1016/j.compositesb.2013.01.027. [19] Najafi, F., Shojaeefard, M.H., Googarchin, H.S. (2016). Low-velocity impact response of functionally graded plates with nonlinear elastic foundation in thermal fields. Composite Part B: Engineering, vol. 10, p. 123-140, DOI:10.1016/j. compositesb.2016.09.070. [20] Reddy, J.N. (1984). A refined nonlinear theory of plates with transverse shear deformation. International Journal of Solids and Structures, vol. 20, p. 881-896, DOI:10.1016/0020­7683(84)90056-8. [21] Shen, H. (1997). Kármán-type equations for a higher-order shear deformation plate theory and its use in the thermal postbuckling analysis. Applied Mathematics and Mechanics, vol. 18, p. 1137-1152, DOI:10.1007/BF00713716. [22] Huang, X.-L., Zheng, J.-J. (2003). Nonlinear vibration and dynamic response of simply supported shear deformable laminated plates on elastic foundations. Engineering Structures, vol. 25, p. 1107-1119, DOI:10.1016/S0141­0296(03)00064-6. [23] Pearson, C.E. (1986). Numerical Methods in Engineering & Science. Van Nostrand Reinhold Company, New York. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9 Vsebina Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 68, (2022), številka 9 L jubljana, september 2022 ISSN 0039-2480 aa eeo R az širjeni povz etki (extended abstracts) Matjaž Ramšak: Fraktalna geometrija reber za ucinkovit hladilnik SI 6 Youyu Liu, Liteng Ma, Siyang Yang, Liang Yuan, Bo Chen: Uporaba metod MTPAin MSM za analizo prenosa vibracij po manipulatorju za vrtanje sidrnih vrtin s 6 prostostnimi stopnjami SI 70 Kuat Kombayev, Murat Muzdybayev, Alfiya Muzdybayeva, Dinara Myrzabekova, ojciech ieleba, Tadeusz Leniewski: Utrjevanje funkcijskih površinskih slojev in izboljševanje obrabne obstojnosti strojnih elementov iz maloogljicnega jekla s plazemsko elektrolizo SI 71 Mateusz rzochal, Stanisaw Adamczak, Grzegorz Piotrowicz, Sylwester nuk: Eksperimentalna raziskava v industriji kot prispevek k razvoju eksperimentalnega modela vibracij kotalnih ležajev SI 72 Yanbing Ni, enliang Lu, Shilei Jia, Chenghao Lu, Ling hang, Yang en: Varovalna metoda za omejevanje delovnega prostora paralelne delovne glave SI 73 iaolin Huang, Chengzhe ang, Jiaheng ang, Nengguo ei: Analiza nelinearnih vibracij funkcionalno gradientnih poroznih plošc, ojacenih s plošcicami grafena, na nelinearnih elasticnih temeljih SI 74 Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, SI 69 Prejeto v recenzijo: 2022-01-24 © 2022 Avtor. Prejeto popravljeno: 2022-04-25 Odobreno za objavo: 2022-06-13 ratala eoetria reer a ioit ladiliN Matjaž Ramšak Univerza v Mariboru, Fakulteta za strojništvo, Slovenija Fraktalna geometrija se v naravi srecuje v vec oblikah in za razlicne namene. Omenimo tok krvi v razvejanem žilnem tokokrogu od glavnih žil do kapilar. V tehniškem žargonu je namen prekrvavitve tkiv prenos mase in tudi prenos toplote. nano je, da mnogo evolucijskih rešitev prekaša današnje inženirske na mnogih podrocjih. e to velja tudi za prenosnike toplote, bi morala obstajati fraktalna geometrija hladilnih reber, ki bi bila ucinkovitejša od preprostih ravnih reber, ki so danes najpogosteje uporabljena za hlajenje naprav, kot pravi priljubljen slogan: obstaja prostor za izboljšave. Kako dolga je obala Britanije je vprašanje, ki ga je zastavil pionir fraktalne matematike Mandelbrot. uporabo vedno krajšega merila se dolžina obale povecuje v neskoncnost. e to logiko uporabimo za fraktalno oblikovana hladilna rebra, dobimo neskoncno hladilno površino (A ) in ob nic razlicni prestopnosti toplote (h) bi morali dobili neskoncno toplotno moc ( Q . ), ki je produkt obeh Q . hA T , pri cemer je T konstantna temperaturna razlika. Cilj prispevka je preveriti omenjeno idejo z uporabo Richardsonove ekstrapolacije rezultatov numericnih simulacij prestopa toplote iz zaporedja Kochovih snežink na tekocino pri variaciji dolžine osnovnega fraktalnega elementa od 1 do 0. a simulacijo prenosa toplote smo uporabili lasten program za Racunalniško dinamiko tekocin na osnovi Metode robnih elementov (BEM). Katedra za energetiko in procesno strojništvo na Fakulteti za strojništvo ima že 40 letno tradicijo razvoja BEM za prenosne pojave v trdninah in tekocinah. V tem prispevku uporabljamo BEM z mešanimi elementi in s tehniko podobmocij v limitnem režimu, kjer je vsako podobmocje sestavljajo 3 robni elementi za trikotno podobmocje in 4 za štirikotno. Na tak nacin se izognemo polni nesimetricni matriki, ki je glavna hiba BEM in dobimo prazno matriko kot pri metodi koncnih elementov. Razviti programom smo uspešno validirali in objavili revizijo Benchmark primera vezanega prenosa toplote v trdnini in tekocini. Ideja o neskoncni toplotni moci fraktalnih hladilnih reber je seveda naivna, kar je potrdil tudi pricujoci numericni eksperiment. a neskoncno velikost fraktalne površine (A ) smo izracunali vrednost toplotne prestopnosti enako nic (h 0). Limita toplotne moci, ki je produkt velikosti površine in toplotne prestopnosti, je koncna vrednost ( Q . konstanta). Slednje je splošno znano dejstvo. Dejstvo je tudi, da je vsaka površina hrapava najsi je fraktalne ali kakšne druge oblike in da je prestop toplote iz trdne stene na tekocino pravzaprav difuzija toplote. amišljeni fraktalni hladilnik v obliki Kochove snežinke ni konkurencen ravnim hladilnim rebrom. Nadaljnje raziskave v tej smeri nimajo pravega smisla. le eede ratalo ladilo telo E i laee oirai reo tolote laiari to etoda roi eleeto ocoa eia Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, SI 70 Prejeto v recenzijo: 2022-01-17 © 2022 Avtorji. Prejeto popravljeno: 2022-07-07 Odobreno za objavo: 2022-08-26 oraa etod i a aalio reoa iraci po manipulatorju z a vrtanje sidrnih vrtin s 6 prostostnimi stopnjami Youyu Liu1,3 – Liteng Ma1,3, – Siyang Yang1,3 – Liang Yuan2 – Bo Chen1,3 1Državni laboratorij za senzoriko in pametno upravljanje visokotehnološke opreme, Ministrstvo za izobraževanje, Kitajska 2 uhu Yongyu avtomobilska industrija, Tehnološki oddelek, Kitajska 3Politehnika Anhui, Šola za strojništvo, Kitajska Strojno vrtanje sidrnih vrtin za podporne sisteme v premogovnikih odpravlja težko rocno delo, povzroca pa mocne vibracije, ki se od stroja za sidranje po manipulatorju prenašajo na podstavek. lanek obravnava mehanizem prenosa vibracij po manipulatorju med vrtanjem sidrnih vrtin. TPAje eno od orodij za preucevanje prenosa vibracij, kot so tudi OTPA, GTDTin ITPA. Obstajajo dolocene omejitve, npr. zaradi motecih dejavnikov, kot sta sklapljanje vzbujanj in šum, manevrirna sposobnost pa je slaba zaradi serijskega sistema z mnogimi podstrukturami. Na podlagi vecstopenjske analize prenosne poti (MTPA) in metode modalne superpozicije (MSM) je bil postavljen model prenosa vibracij za podsistem, sestavljen iz zgibov manipulatorja s šestimi prostostnimi stopnjami (DOF). Dolocena je bila matrika njegovega frekvencnega odziva. Izpeljano je bilo vzbujanje stroja za sidranje v šestih prostostnih stopnjah. Vzbujalna sila stroja za sidranje, ki se po manipulatorju prenaša na podstavek, je bila analizirana z Jacobijevo matriko sil in na ta nacin so bile dolocene zunanje obremenitve podstavka. Primer iz inženirske prakse kaže, da si amplitude prostostnih stopenj podstavka sledijo v naslednjem vrstnem redu od vecjih proti manjšim: upogibne vibracije (komponenta 1), vzdolžne vibracije, torzijske vibracije, upogibne vibracije (komponenta 2), rotacijske vibracije okrog osi z in rotacijske vibracije okrog osi y. Na podstavku se pojavljajo predvsem upogibne vibracije. (1) Na podlagi metod MTPA in MSM je bil postavljen matematicni model prenosa vibracij po manipulatorju pri vrtanju sidrnih vrtin. Izpeljana je bila matrika frekvencnega odziva podsistemov na vzbujanje v vec prostostnih stopnjah. Na podlagi dolocitve zunanjih vzbujalnih obremenitev je mogoce izracunati funkcijo frekvencnega odziva vsake prostostne stopnje v tocki odziva po matematicnem modelu, ki je univerzalen za serijske sisteme. (2) notraj frekvencnega obmocja vzbujanja stroja za sidranje (0 Hz do 200 Hz) se nahajajo frekvencna obmocja z najvecjimi vrednostmi vibracij v šestih prostostnih stopnjah posameznih podsistemov 11,3, 17,4 Hz, 88,5, 1,55 Hz, 42,1, 1,77 Hz, 135,2, 137,8 Hz, 178,8, 187,4 Hz in 73,24, 5,54 Hz. Vibracijski odziv je najvecji v omenjenih frekvencnih obmocjih, resonanci pa se je mogoce izogniti s spremembo vzbujalne frekvence stroja za sidranje. (3) Primer iz inženirske prakse kaže, da si amplitude prostostnih stopenj podstavka sledijo v naslednjem vrstnem redu od vecjih proti manjšim: upogibne vibracije (komponenta 1) (2,1210-2 m) pri 0 Hz, vzdolžne vibracije (1,6510-2 m) pri 45 Hz, torzijske vibracije (10-3 m) pri 10 Hz, upogibne vibracije (komponenta 2) (8,0610-3 m) pri 180 Hz, rotacijske vibracije okrog osi (2,4610-4 m) pri 180 Hz, rotacijske vibracije okrog osi y (2,0810-4 m) pri 180 Hz. Na podstavku se pojavljajo predvsem upogibne vibracije. (4) Izkazalo se je, da resonanca nastopi med komponentama upogibnih vibracij pri frekvencah 0 Hz in 180 Hz, resonanca med komponentama rotacijskih vibracij okrog osi y in pa je najbolj verjetna med frekvencama 180 Hz in 181 Hz. Predlagana teorija prenosa vibracij po manipulatorju za vrtanje sidrnih vrtin s šestimi prostostnimi stopnjami bo lahko teoreticna osnova za razvoj tehnik blaženja vibracij in konstruiranje blažilnikov. Na podlagi metod MTPA in MSM je bil postavljen matematicni model prenosa vibracij po manipulatorju za vrtanje sidrnih vrtin s šestimi prostostnimi stopnjami. unanja obremenitev v tocki odziva na podstavku manipulatorja je bila analizirana z Jacobijevo matriko sil. Dolocena sta bila vibracijski odziv po vsaki prostostni stopnji podstavka in resonancna frekvenca. Študija primera iz inženirske prakse je pokazala, da je podstavek izpostavljen predvsem upogibnim vibracijam. le eede ailator tro a idrae etoea aalia reoe oti etoda odale superpoz icije, prenos vibracij, Jacobijeva matrika sil SI 70 *Naslov avtorja za dopisovanje: Wuhu Yongyu avtomobilska industrija, Tehnološki oddelek, Wuhu, Kitajska, ahpulyy@163.com Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, SI 71 Prejeto v recenzijo: 2022-04-05 © 2022 Avtorji. Prejeto popravljeno: 2022-08-02 Odobreno za objavo: 2022-08-10 treae cii orii loe in iz boljševanje obrabne obstojnosti strojnih elementov i alooliea ela laeo eletrolio Kuat Kombayev1 - Murat Muzdybayev1 - Alfiya Muzdybayeva1 - Dinara Myrzabekova1 - ojciech ieleba2 - Tadeusz Leniewski2 1 Državna tehniška univerza vzhodnega Kazahstana, Oddelek za tehnološko opremo in transport, Republika Kazahstan 2 nanstveno-tehniška univerza v roclawu, Fakulteta za strojništvo, Poljska Plazemska elektroliza (EPP) je najucinkovitejši postopek za utrjevanje delov vrtalne opreme. Vpricujoci raziskavi je bila opravljena eksperimentalna identifikacija parametrov tehnološkega procesa utrjevanja jekla 18CrNi3Mo­Sh s plazemsko elektrolizo. Tehnologija utrjevanja površinskega sloja strojnih elementov iz maloogljicnega jekla je možna alternativa ogljicenju in kaljenju (postopek, ki se uporablja pri kotalnih dletih). Predlagana metoda za utrjevanje površinskega sloja delov iz maloogljicnega jekla vkljucuje segrevanje dela do temperature 30 C do 40 C, ki mu sledi naogljicenje površine z ogljikovimi ioni do globine 2 mm in nato kaljenje na temperaturi 800 C do 820 C. a postopek je znacilno plazemsko segrevanje površine dela, ki je potopljen v tekocem elektrolitu do globine 4 mm do 6 mm. a preverjanje ucinkovitosti utrjevanja po postopku proizvajalca kotalnih dlet JSC VKM in po metodi EPP sta bili obdelani dve seriji preizkušancev. Rastrska elementna analiza obdelane katodne površine je pokazala, da EPPpoleg utrjevanja povzroci tudi kemicno modifikacijo površinskega sloja komponente. EPPomogoca izdelavo utrjenih slojev debeline 1000 m do 1700 m. Ugotovljeno je bilo 1,5- do 2-kratno povecanje mikrotrdote v primerjavi z zacetnim stanjem. a mikrostrukturo obdelane površine jeklenih preizkušancev je znacilen temen modificirani površinski sloj. Fina iglicasta struktura martenzitnega izvora pod temnim slojem se spremeni v zacetno perlitno-feritno mikrostrukturo. Prednost utrjevanja s plazemsko elektrolizo je v majhni porabi energije, visoki hitrosti utrjevanja ter zmožnosti lokalne obdelave površin, zlasti pri vecjih delih zahtevnih oblik. Predlagana metoda za obdelavo površin po postopku EPP ne zagotavlja le gladkih površin, temvec tudi izboljša delovne lastnosti komponent, npr. obrabno obstojnost. Ta je bila ocenjena z izgubo mase preizkušanca na enoto casa v obrabnih preizkusih z drsenjem po abrazivnem disku brez maziva. Mikrorelief obrabljene površine jeklenega preizkušanca je pokazal najmanjšo abrazivno obrabo pri delu, obdelanem po postopku EPP. Pri finozrnatem martenzitu se na površini oblikuje mocna mikrostruktura, ki je ni unicila niti mocna abrazivna obraba (med eksperimenti) in ki tako preprecuje razsežnejše poškodbe. Jekleni preizkušanec, ki je bil toplotno obdelan po postopku JSC VKM, je inferioren v primerjavi s preizkušancem, utrjenim po postopku EPP. Eksperimenti so tako potrdili, da je z reguliranjem postopka EPP mogoce vplivati na kakovost utrjevanja jekel, pridobiti obrabno obstojne površine in znatno izboljšati produktivnost EPP. Tehnološki postopek EPP je obvladljiv ter primeren za prakticno uvedbo razvite tehnologije v proizvodnji. le eede aloolio elo treae laea eletrolia relee ieiri ori obrabna obstojnost Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, SI 72 Prejeto v recenzijo: 2022-05-04 © 2022 Avtorji. Prejeto popravljeno: 2022-07-05 Odobreno za objavo: 2022-08-17 Eksperimentalna raz iskava v industriji kot prispevek k raz voju eksperimentalnega modela iraci otali leaeY Mateusz rzochal1, Stanisaw Adamczak1, Grzegorz Piotrowicz2, Sylwester nuk2 1Tehniška univerza v Kielcah, Poljska 2Polish Bearing Factory, Kranik SA, Poljska V odsotnosti fizicnih simulatorjev za kotalne ležaje obstaja realna potreba po razvoju matematicnih modelov za cim tocnejšo aproksimacijo delovanja kotalnih ležajev. a dosedanje teoreticne študije z matematicnimi modeli kotalnih ležajev so znacilne velike poenostavitve in obravnava le nekaterih dejavnikov, ki vplivajo na raven vibracij ležajev. Trenutno ni matematicnega modela, ki bi upošteval vecino realnih dejavnikov, pri analizi vibracij pa je najveckrat prezrta kontaminacija ležajev in masti. Poleg tega vecina analiticnih modelov obravnava ležaje z izrazitimi napakami, ki bi nastale le v primeru dolgega obratovanja ali neprimernih obratovalnih pogojev. Tovrstni modeli v razpoložljivi svetovni literaturi zato niso univerzalni in niso primerni za uporabo v industrijski praksi. Potreben je eksperimentalni model za popolnoma nove ležaje. Teoreticni modeli vibracij ležajev in eksperimentalni preizkusi obicajno prezrejo cistoco v notranjosti ležajev. Matematicna simulacija kontaminacije je prvic zelo težavna, drugic pa eksperimentalne študije obravnavajo visokoamplitudne komponente v nizko- in srednjefrekvencnem obmocju, ki so posledica poškodb ali geometrije površin. Da bi razkrili glavne vzroke za izmet, je bil šest mesecev opazovan proces proizvodnje desetih razlicnih vrst kotalnih ležajev. Preizkušenih je bilo skupaj 46.811 ležajev. Vsi izdelki, ki so bili izloceni po meritvah vibracij v okviru kontroli kakovosti, so bili skrbno pregledani. Vzroki napak na izdelkih so bili razvršceni v štiri kategorije. a pripravo eksperimentalnega modela kotalnega ležaja je treba opraviti vrsto preizkusov, ki upoštevajo vse dejavnike vpliva na stopnjo vibracij, in vkljuciti tudi veckriterijsko statisticno analizo. Model mora upoštevati nepopolnosti, ki vplivajo na diskretne vibracije (odstopanja oblike, velikosti in položaja, cezmerna valovitost, mikrovalovitost in hrapavost), ter predvsem prisotnost trdnih delcev v ležaju, ki je neizogibno povezana s proizvodnim procesom. eprav so bile študije opravljene na zelo veliki skupini razlicnih ležajev, specificnih vzrokov kontaminacije ne opisujejo. Njihovo poznavanje bi lahko omogocilo izboljševanje tehnoloških procesov izdelave kotalnih ležajev. Avtorji clanka nameravajo v prihodnje opraviti podrobnejše raziskave in analizo vpliva stopnje kontaminacije kotalnih ležajev na vibracije ter preveriti stopnjo zaznavanja tovrstnih neskladnosti z merilno opremo in metodami, ki jih uporabljajo v industriji ležajev. nacrtovanimi preizkusi bo mogoce napovedati vpliv cistoce ležajev na omejevanje vibracij. Neupoštevanje cistoce v matematicnih modelih in eksperimentalnih raziskavah ni upraviceno. Proizvajalci ležajev, ki se ukvarjajo z razvojem preciznih ležajev, morajo ob naložbah v proizvodnjo razviti tudi procese za zagotavljanje tehnicne cistoce. Na koncno stopnjo tehnicne cistoce komponent, ki se vgrajujejo v koncne sestave, ne vplivata le proizvodni proces in pranje. Kljucnega pomena so tudi stopnja cistoce v prostorih, kjer potekajo proizvodne in montažne operacije, ter embalaža in ukrepi za preprecitev ponovne kontaminacije. le eede otali leai iracie itoa leae otrola aooti Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, SI 73 Prejeto v recenzijo: 2022-02-21 © 2022 Avtorji. Prejeto popravljeno: 2022-07-01 Odobreno za objavo: 2022-07-11 Varovalna metoda z a omejevanje delovnega prostora paralelne delovne glave Yanbing Ni1,– enliang Lu2 – Shilei Jia2 – Chenghao Lu2 – Ling hang2 – Yang en1 1 Univerza v Tianjinu, Fakulteta za strojništvo, Kitajska 2 Državni laboratorij za teorijo mehanizmov in konstrukcijo, Ministrstvo za šolstvo, Kitajska Paralelna delovna glava 1T2R je topološko zasnovana kot paralelni mehanizem z eno translacijsko in dvema rotacijskima (1T2R) prostostnima stopnjama. Lega konca mehanizma znotraj delovnega prostora je dolocena z nelinearno preslikavo delovanja servomotorja v prostor sklepov. Informacije o legi konca mehanizma je mogoce pridobiti samo s kompleksnim modeliranjem in izracuni. Kljucni problem, ki ga je treba nujno razrešiti za tocno in realnocasovno omejevanje gibanj tovrstnih obdelovalnih strojev, je torej zmanjšanje obsega obdelave in racunanja pri ocenjevanju lege konca mehanizma. V clanku sta predstavljeni zasnova in implementacija poenostavljenega algoritma za ocenjevanje lege omenjenih mehanizmov. Najprej je sistematicno predstavljena paralelna delovna glava 1T2R, temu pa sledi predstavitev kinematicnega modela paralelnega mehanizma. Sledita analiza položaja ter oblikovanje preslikave med vrednostmi vhodov in gibanjem konca mehanizma z inverzno kinematicno rešitvijo. Na osnovi konstrukcije mehanizma, parametrov skaliranja, obmocja gibanja ter interferenc posameznih kinematicnih parov in drugih omejitev sta doloceni množica vseh dosegljivih prostorskih položajev referencne tocke orodja na koncu mehanizma in množica dosegljivih orientacij orodja v dani referencni tocki. analizo zakona gibanja glave 1T2R je nato izpeljan poenostavljen algoritem za oceno njenega položaja in orientacije. Položaj referencne tocke na orodju glave je ocenjen s pomocjo vsote dolžin treh nog verige. Orientacija glave je ocenjena s pomocjo razlik v dolžini nog katerih koli dveh od treh vej. Kombinacija obeh zakonov omogoca hitro oceno položaja in orientacije glave z nezahtevnimi izracuni ter hiter odziv varovalne omejevalne funkcije z obstojeco strojno opremo za zašcito pred prekoracitvijo delovnega obmocja glave. Pravilnost predlaganega poenostavljenega algoritma so potrdili tudi rezultati simulacij v nadaljevanju. Poenostavljeni algoritem za ocenjevanje položaja in orientacije je bil uporabljen v kombinaciji z mehansko konstrukcijo in krmilnim sistemom servopodajanja glave pri zasnovi dveh metod omejevanja – aktivne in preventivne. Eksperimenti so koncno dokazali, da je s poenostavljenim algoritmom za oceno položaja in lege izvedljivo realnocasovno ocenjevanje položaja in lege konca mehanizma med obratovanjem paralelne delovne glave 1T2R. Algoritem daje tocne ocene položaja in lege. Na osnovi rezultatov kinematicne analize in analize nacina gibanja mehanizma je bil oblikovan predlog algoritma za oceno lege paralelne delovne glave 1T2R. Ta v primerjavi z metodo za ocenjevanje lege po modelu inverznega položaja prihrani ,42 racunskega casa. a razliko od metod za ocenjevanje lege, ki uporabljajo zunanje senzorje, ne zahteva nadgradnje originalnega sistema ter je ucinkovit, zanesljiv in enostavno izvedljiv. Poenostavljeni algoritem za oceno položaja in lege tako omogoca preprosto in zanesljivo realnocasovno omejevanje položaja glave. Algoritem trenutno pokriva celotni dosegljivi prostor mehanizma v obmocju omejevanja in tako izboljšuje zmožnosti hitrega odziva omejevalne funkcije za varnost obratovanja obdelovalnih strojev s tako topologijo. le eede aralela deloa laa odelirae aalia deloea rotora ocea oloaa i lee varovanje z omejevanjem, inverz na kinematika Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, SI 74 Prejeto v recenzijo: 2022-07-08 © 2022 Avtorji. Prejeto popravljeno: 2022-09-02 Odobreno za objavo: 2022-09-07 Analiz a nelinearnih vibracij funkcionalno gradientnih oroi lo oaei loicai raea a elieari elatii teeliK iaolin Huang – Chengzhe ang – Jiaheng ang – Nengguo ei Univerza za elektroniko v Guilinu, Šola za arhitekturo in transport, Kitajska Pricujoci clanek obravnava nelinearne vibracije plošc. Predstavljen je model na podlagi modificiranega mikromehanskega modela Halpin-Tsai in razširjenega zakona zmesi, namenjen vrednotenju lastnosti nanokompozitov z notranjimi porami, ojacenih z grafenom. Predmet obravnave je funkcionalno gradientna porozna plošca na nelinearnih elasticnih temeljih, ojacena z grafenom. Upoštevane so tri vrste porazdelitve por, tako enakomerne kot neenakomerne. V študiji je predstavljen modificirani model, ki upošteva volumski delež por. Na podlagi modificiranega modela so bile ocenjene dejanske materialne lastnosti porozne plošce. Izpeljane so enacbe nelinearnih vibracij plošce v okviru teorije plošc višjega reda in splošnih Krmnovih enacb, ki upoštevajo vpliv sklopitve z nelinearnimi elasticnimi temelji. a dolocitev nelinearnega frekvencnega in dinamicnega odziva sta bili uporabljeni tehnika dvostopenjskih perturbacij in Galerkinova metoda. Podrobno so preuceni vplivi por, plošcic grafena in elasticnih temeljev na lastno frekvenco porozne plošce. V predhodnih študijah je bilo ugotovljeno, da se lastna frekvenca monotono znižuje z rastjo koeficienta poroznosti. Izkazalo se je, da pravilo ne velja za dani primer. Koeficient poroznosti pri porazdelitvah por P-1 in P-3 vpliva na vrednosti parametrov temeljev. Pri vrednostih parametrov temeljev (50, 0) in (0, 0) se lastna frekvenca zmanjša, za lastno frekvenco pri parametrih (50, 50) pa velja nasprotno. Preucen je bil tudi vpliv prehodnega odklona. Najvecji odklon se ob dvigu koeficienta poroznosti z 0,0 na 0,4 poveca za približno 8 Sledi sklep, da je vpliv koeficienta poroznosti na dinamicni odziv zanemarljiv. Reševanje nelinearnih vodilnih enacb gibanja porozne plošce je zelo težavno, zato je bilo mogoce pridobiti samo asimptotske rešitve tretjega reda. Avtorji nameravajo v prihodnjih raziskavah poiskati asimptotske rešitve višjega reda. V clanku je predstavljen modificiran model za vrednotenje materialnih lastnosti funkcionalno gradientnih nanokompozitov, ojacenih z grafenom. Podanih je vec zanimivih ugotovitev v zvezi z vplivom por in nelinearnih elasticnih temeljev. Rezultati bodo lahko uporabni pri projektiranju podobnih konstrukcij v praksi. le eede cioalo radieti oroi aoooiti loice raea ore elieari elatii temelji, nelinearne vibracije, prehodni odgovor 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 suggested length of contributions is 10 to 20 pages. 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 pub­lished 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 inter­est. 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 high­light what is distinctive about it. - An Introduction that should provide a review of recent literature and sufficient back­ground 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 pre­viously 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 includ­ed. 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 con­tain 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 col­lected together in a reference list at the end of the manuscript. - Appendix(-icies) if any. 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 vest-nik - 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. 553­576. 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 Com­pounds 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 2009­09-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 publish the paper uder CC-BY 4.0 Int. or CC-BY-NC 4.0 Int. when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright 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 https://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. SPECIAL NOTES Units: The SI system of units for nomenclature, symbols and abbreviations should be Strojniški vestnik -Journal of Mechanical Engineering followed closely. Symbols for physical quantities in the text should be written in italics (e.g. Ašker°eva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu http://www.sv-jme.eu Contents Papers 517 Matjaž Ramšak: Fractal Geometry as an Effective Heat Sink 529 Youyu Liu, Liteng Ma, Siyang Yang, Liang Yuan, Bo Chen: MTPA- and MSM-based Vibration Transfer of 6-DOF Manipulator for Anchor Drilling 542 Kuat Kombayev, Murat Muzdybayev, Alfiya Muzdybayeva, Dinara Myrzabekova, Wojciech Wieleba, Tadeusz Lesniewski: Functional Surface Layer Strengthening and Wear Resistance Increasing of a Low Carbon Steel by Electrolytic-Plasma Processing 552 Mateusz Wrzochal, Stanislaw Adamczak, Grzegorz Piotrowicz, Sylwester Wnuk: Industrial Experimental Research as a Contribution to the Development of an Experimental Model of Rolling Bearing Vibrations 560 Yanbing Ni, Wenliang Lu, Shilei Jia, Chenghao Lu, Ling Zhang, Yang Wen: Limit-protection Method for the Workspace of a Parallel Power Head 571 Xiaolin Huang, Chengzhe Wang, Jiaheng Wang, Nengguo Wei: Nonlinear Vibration Analysis of Functionally Graded Porous Plates Reinforced by Graphene Platelets on Nonlinear Elastic Foundations