Since 1955 „ > Strojniški vestnik Journal of Mechanical Engineering Strojniški vestnik - Journal of Mechanical Engineering (SV-JME) Aim and Scope The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis. The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue. The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s). Editor in Chief Vincenc Butala University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Technical Editor Pika Škraba University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Founding Editor Bojan Kraut University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerčeva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www.sv-jme.eu Print: Koštomaj printing office, printed in 275 copies Founders and Publishers University of Ljubljana, Faculty of Mechanical Engineering, Slovenia University of Maribor, Faculty of Mechanical Engineering, Slovenia Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association President of Publishing Council Mitjan Kalin University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Vice-President of Publishing Council Bojan Dolšak University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: Once a year, a hydropower plant general refit takes place, and the turbine trash rack is dismantled, inspected and cleaned. The right-hand Figure shows the trash rack full of permanently wedged debris after one year of operation, while the assembled clean trash rack is presented on the left. The research presented in the journal is focused on the prediction of instantaneous trash rack clogging, its influence on head losses, and analysing different cleaning strategies within a one year period. Image courtesy: University of Maribor, Faculty of Mechanical Engineering, Heat Machines and Measurements Laboratory, Slovenia & DEM d.o.o, Maribor, Slovenia ISSN 0039-2480, ISSN 2536-2948 (online) International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, King Fahd U. of Petroleum & Minerals, Saudi Arabia Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, University of Ljubljana, Slovenia Filippo Cianetti, University of Perugia, Italy Janez Diaci, University of Ljubljana, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Jožef Duhovnik, University of Ljubljana, Slovenia Igor Emri, University of Ljubljana, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Janez Grum, University of Ljubljana, Slovenia Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, University of Maribor, Slovenia Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan Julius Kaplunov, Brunel University, West London, UK Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany Jernej Klemenc, University of Ljubljana, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, University of Ljubljana, Slovenia Gorazd Lojen, University of Maribor, Slovenia Darko Lovrec, University of Maribor, Slovenia Thomas Lubben, University of Bremen, Germany Jure Marn, University of Maribor, Slovenia George K. Nikas, KADMOS Engineering, UK Tomaž Pepelnjak, University of Ljubljana, Slovenia Vladimir Popovič, University of Belgrade, Serbia Franci Pušavec, University of Ljubljana, Slovenia Mohammad Reza Safaei, Florida International University, USA Marco Sortino, University of Udine, Italy Branko Vasič, University of Belgrade, Serbia Arkady Voloshin, Lehigh University, Bethlehem, USA General information Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue). Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor. Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http:// www.sv-jme.eu. You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content. We would like to thank the reviewers who have taken part in the peer-review process. © 2020 Strojniški vestnik - Journal of Mechanical Engineering. All rights reserved. SV-JME is indexed / abstracted in: SCI-Expanded, Compendex, Inspec, ProQuest-CSA, SCOPUS, TEMA. The list of the remaining bases, in which SV-JME is indexed, is available on the website. The journal is subsidized by Slovenian Research Agency. Strojniški vestnik - Journal of Mechanical Engineering is available on https://www.sv-jme.eu. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2 Contents Contents Strojniški vestnik - Journal of Mechanical Engineering volume 66, (2020), number 2 Ljubljana, February 2020 ISSN 0039-2480 Published monthly Papers Kamil Urbanowicz, Huan-Feng Duan, Anton Bergant: Transient Liquid Flow in Plastic Pipes 77 Guangjian Wang, Li Su, Shuaidong Zou: Uneven Load Contact Dynamic Modelling and Transmission Error Analysis of a 2K-V Reducer with Eccentricity Excitation 91 Tomasz Kozior, Al Mamun, Marah Trabelsi, Lilia Sabantina, Andrea Ehrmann: Quality of the Surface Texture and Mechanical Properties of FDM Printed Samples after Thermal and Chemical Treatment 105 Peilin Zhang, Bo Li, Xuewen Wang, Chaoyang Liu, Wenjie Bi, Haozhou Ma: The Loading Characteristics of Bulk Coal in the Middle Trough and Its Influence on Rigid Body Parts 114 Mihaela Milesan, Constantin Cristinel Girdu, Liviu Cirtina, Constanta Radulescu: Mathematical Modelling Study of Hardox400 Steel Parts' Roughness and Hardness, Cut with CO2 Laser 127 Aleš Hribernik: Evaluation of Clogged Hydropower Plant Trash Rack Losses 142 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2019.6324 Original Scientific Paper Received for review: 2019-09-11 Received revised form: 2019-11-07 Accepted for publication: 2019-12-03 Transient Liquid Flow in Plastic Pipes Kamil Urbanowicz1,* - Huan-Feng Duan2 - Anton Bergant3,4 1 West Pomeranian University of Technology, Department of Mechanical Engineering and Mechatronics, Poland 2 The Hong Kong Polytechnic University, Department of Civil and Environment Engineering, China 3 Litostroj Power d.o.o., Slovenia 4 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Nowadays, plastic pipes play a key role in fluid conveyance, for example, in urban water supply systems. Dynamic analysis of designed water pipe networks requires the use of numerical methods that allow for solving basic equations that describe transient flows occurring in plastic pipes. In this paper, the main equations were formulated with pressure (p) and velocity (v) as the principal variables. A novel simplified retarded strain solution is used to properly model the pipe-wall viscoelasticity effect during transient flow process. Unsteady friction is calculated with the use of a filtered weighting function (with three exponential terms). The proposed numerical method enables the efficient modelling of transient flow in plastic pressurized pipes. Extensive simulations are performed and compared with experimental results known from three different European research centres (London, Cassino, and Lyon), with the aim of demonstrating the impacts of plastic pipe properties and frequency-dependent hydraulic resistance on transient pipe flow oscillations. The effectiveness of the proposed method for determining the creep functions of plastic pipes is also examined and discussed. Keywords: plastic pipes, viscoelasticity, water hammer, transient flow, method of characteristics, unsteady friction Highlights • A proper mathematical model is derived based on the principal variables of pressure and velocity. • An efficient numerical solution of transient flow model in plastic pipes is used. • A simplified, effective numerical convolution integral describing the retarded strain is presented. • A detailed qualitative and quantitative analysis for the results of comparative tests is carried out. • A discussion on the influence of unsteady friction on transient flows is presented. 0 INTRODUCTION The work of Tison, written in 1958, [1] has been considered to be foundational regarding the dynamic behaviour of the viscoelastic pipe wall and its influence on the transient flow behaviour associated with a rapid valve opening. In 1962, Hardung [2] described the physical phenomena that come into play in the arterial system as a consequence of heart action. He studied the influence of internal wall friction and liquid viscosity on damping of pressure wave velocity. From the linearized dynamic equations of viscous incompressible fluid flow in viscoelastic tubes, Martin [3] calculated frequency response of specific tube flow system in 1969. Lorenz and Zeller [4] analysed wave propagation in viscoelastic pipes, to find the analytical solution based on an analogy with a piston problem. Their main conclusion was that the shock waves could not be formed in viscoelastic pipes like those in elastic ones. Instead, a fully dispersed wave is formed, which has a width of several times of the tube diameter. In these early scientific studies, the main focuses was in analytical analysis. A milestone to modern modelling was attributed to deriving general equations of the transient flow in pipes with viscoelastic properties by Rieutord and *Corr. Author's Address: Department of Mechanical Engineering and Mechatronics, West Piastow 19, Szczecin 70-310, Poland, kamil.urbanowicz@zut.edu.pl Blanchard in 1972 [5]. In 1977, Güney, in his PhD thesis [6], presented the first working numerical method (based on finite difference method) that included a laminar unsteady friction [7] and was able to calculate the pressures and mean velocities of flow in plastic pipes during water hammer event including column separation. For cavitation modelling, Güney successfully applied the simplest discrete vapour cavity model. Meanwhile, Rieutord and Blanchard [8] presented a theoretical study in a one-dimensional form on the effect of viscoelastic material properties on the pipe-conveying process. Their analysis results indicated that, in addition to exponential attenuation of wavefronts, the disturbance is transformed into a diffusion front. Limmer and Meißner [9] incorporated the complex creep compliance into the transient flow models and then derived the wave speed and damping factor for an oscillating pressure wave propagating in a thin-walled viscoelastic pipeline. Franke and Seyler [10] utilized an impulse response method to calculate water hammer. In Italy, in the mid-1980s, the idea of reducing the unsteady flow oscillations in pressure pipelines by inserting an in-line polymer section, which is characterized with low pressure wave speed, was examined by Ghilardi and Paoletti [11]. This ranian University of Technology Szczecin, 77 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 solution was then further studied theoretically and experimentally by Pezzinga and Scandura [12]. The authors used a Kelvin-Voigt model with one element. The coefficients of creep function were calibrated based on experimental results. The additional plastic pipe in their tests significantly reduced unsteady-flow oscillations. Tijsseling [13] presented an extended literature review about FSI (fluid-structure interaction) in steel and plastic pipes. At the beginning of this millennium, Covas et al. investigated transient flows in a laboratory high-density polyethylene (HDPE) pipe with a length of 277 m [14] to [16]. A new simplified model for simulating viscoelastic transient pipe flows was implemented into the numerical scheme of the method of characteristics (MOC), which allows a fast calculation of pressure runs in viscoelastic conduits. Thereafter, Duan et al. [17] argued that an inaccurate representation of unsteady friction or turbulence in one-dimensional (1D) model may lead to non-physical calibration of viscoelastic parameters and that the influence mechanism of pipe wall viscoelasticity is totally different from friction effects (steady and unsteady) in terms of both amplitude damping and phase shifting of transient pressures. Their results have been validated by the experimental test data of Brunone and Berni [18]. Similar results regarding the contribution of unsteady friction and viscoelasticity were also obtained by Bergant et al. [19]. The fluid-structure interaction (FSI) in plastic pipes was studied by Keramat et al. [20]. The authors found that damping in transient flow may be attributed to friction effects (steady and unsteady), valve resistance, gas (air) in liquid (dissolved and trapped), pipe-wall viscoelasticity, as well as fluid-structure interaction and other non-elastic behaviours of supports. As many of these effects are unknown and have not been properly included in the transient solver, the KelvinVoigt model not only represents viscoelasticity but also all the other effects. As a result, the calibrated results of viscoelastic parameters for such complex situations (i.e., various influence factors) may become non-physical. Pezzinga et al. [21] used a microgenetic algorithm to calibrate creep compliance function coefficient. Calibration of Kelvin-Voigt models with 2, 3, 5, and 7 parameters, respectively, proved the substantial independence of the elastic modulus and the dependence of the retardation time on the pipe period (i.e., the pipe length). In other recent studies, e.g., [22] to [24], it is concluded that the one-element Kelvin-Voigt model is reasonably accurate to achieve satisfactory numerical simulation for typical transient flows in common plastic pipes. Kodura [25] showed that the characteristic of butterfly valve closure has a significant influence on water hammer in polyethylene (PE) pipes for the valve closure time larger than 25 % of the wave characteristic period. For shorter closure time, good agreement could be obtained between the model and experimental test in terms of maximum pressures. More recently, Ferrante and Capponi [26] introduced a generalized Maxwell model, based on fractional derivatives, to better represent the viscoelastic behaviour of plastic pipes during transient flow process. Through different tests for HDPE and oriented polyvinyl chloride (PVC-O), they concluded that the proposed model performs slightly better than the well-known generalized Kelvin-Voigt model. Today, not only pipes but also many hydraulic devices can be made of different plastics, as needed [27]. As a result, various previous studies on unsteady flows in plastic pipes have demonstrated the importance and necessity of this research subject. Note that in plastic pipes during water hammer process, the same phenomena may occur that are known for the flow in metal pipes, such as cavitation, frequency-dependent friction and fluid-structure interaction (FSI). However, in plastic pipes, the main factor responsible for transient damping is the time-dependent viscoelastic material behaviour, which is the main difference from the elastic pipes studied previously. The objective of this work is to present a mathematical model and its effective numerical solution method that allows the simulation of the phenomenon of transient liquid flows in plastic pipes with sufficient accuracy and efficiency. The model will consider two important factors: unsteady friction of pipe flows and viscoelasticity of retarded strain. These two factors described by convolutional integrals require effective solutions to accelerate the calculation process. The integral describing unsteady friction was effectively solved by Urbanowicz [28], using the suggestions of Vardy and Brown [29]; in contrast, the integral for the viscoelastic nature of retarded strain has been previously solved in a complicated manner [14] and [30]. In this study, it will be solved in a much simpler way, using Schohl's effective solution [31], to further improve its effectiveness. An attempt is made to assess, in both qualitative and quantitative manners, the impacts of pipe-wall viscoelasticity and unsteady friction on transient flows in plastic pipes. The situations of laminar, transitional, and turbulent flows in pipes are tested and analysed. 78 Urbanowicz, K. - Duan, H.-F. - Bergant, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 1 MATHEMATICAL MODEL The hyperbolic partial differential equations of motion and continuity, which are described by the basic parameters representing flow: p pressure and v averaged velocity in the cross-sectional area of the pipeline, including the pipe wall friction in a horizontal pipe, can be defined in the following form [5]: 1 dp dv „'r dP(u ) - + — + a I —— • w pc d' dx d' J ('-u) du = 0 dv dp 32« 16« P- + — + -V + -T- d' dx D2 .dv, (1) (u) D2 0 d' • wu du = 0 Among the different numerical methods enabling resolving the system of above equations, particular attention should be paid to the MOC scheme, which perfectly interprets the essence of physical phenomena of transient flow, and at the same time is characterized by fast convergence, ease of incorporating various boundary conditions, together with the high accuracy of calculation results. The detailed numerical solution may be found in recent papers [30], [32] and [33]. 1.1 Simplified Numerical Solution of Retarded Strain In the case of flow in plastic pipes, the equation of continuity has the same form as for metal pipes. A polymer pipeline does not respond according to Hook's law when it is subjected to a certain instantaneous stress a. It consists of an immediate-elastic response and a retarded-viscous response. Thus the strain can be decomposed into a sum of instantaneous-elastic strain ee and a retarded strain er [6] and [16]. The partial time derivative of retarded strain is a convolution integral of pressure and derivative of the creep function J that describes viscoelastic behaviour of the pipe-wall material as: de ) dt s r dp(u ) d — —— • w,,, )du. 2 J dt J(>-u) (2) An efficient numerical solution scheme to the above type of convolution integral has been developed by Schohl [31]: ds (t+At) dt z ) S J,. —e • +--- 2 At 1 - e Ut+A,) [P(t+At) P(t)) (3) where yt is the time-dependent recursion coefficient, of which the starting value (before transient) is equal to zero. Re-ordering the equation yields: -p^tx-t{x'P(,)-rt,)Y')' (4) r(+At) dt » -— " J where: Y. =— e R and X. =— ' 2 ' 2 At 1 - e Let us now define the time-dependent term G^ in Eq. (4): G«=Z (XiP(. )~yK')Y< ). (5) Then, the final form of the simplified equation that describes the derivative of retarded strain is as follows, de„ r(t+At ) dt P(t+At)X - G(t)• (6) Comparing with an existing solution that was introduced by Covas et al. [14] and [16], this novel solution simplifies the numerical algorithm for analysis of transient events in engineering polymer pipes. It may be very useful for the various pipe engineering systems: design of new pressure systems based on plastic pipes; the modernization of existing plastic pressurized pipe systems; and design strategies to protect systems against large transient loads. 1.2 Numerical Solution of Wall Shear Stress The pipe wall friction effect plays a significant role in modelling transient events in elastic pipes. In viscoelastic pipes, some authors note the importance and impact of unsteady friction on transient flow process [6], [16], [17] and [33], while the vast majority try to simplify the model to a minimum and model unsteady friction by calibrating coefficients in the creep function. The action of the latter means that the calibrated function is no longer a material function describing the strength properties of the pipe material. This function, in a simplified way, describes together all the phenomena that accompany the transient flow (including but not limited to: unsteady friction resistance, the impact of FSI and the phenomenon of energy dissipation resulting from retarded strain of pipe walls) [20]. In this work, due to the fact that only experimentally obtained creep functions are used, there is a need to model resistance in a conventional way, i.e., by including the unsteady component Tu contained in the formula for describing the wall shear stress. Urbanowicz [28] has shown that there is a lack 0 i=i Transient Liquid Flow in Plastic Pipes 79 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 of compliance of the effective friction weighting function with original classic weighting function for small time scales. The corrected final effective solution of convolution integral representing the additional contribution due to unsteadiness becomes [28]: Pft. (+Át) I (+Át) I (+Át) 4Mf D h 8 M nB. [v(,+Át)-v(t)]■ [l-n]C¡ [v(, )- v(,-Át )] (7) where n = f ' J 0 (u ~)du f o'wejr(u )du B = ^ [1 - A ], nAt y,(t+M) A,. = e- C = AB. The correction factor n in the above solution is a quotient of the integral of the classic weighting function wc]as (for laminar flow, Zielke's function [7], and for turbulent flow, Vardy-Brown's function [34]) and the integral of effective weighting function wei. (approximated by a finite sum of exponential terms): w eff. =h me (9) which is just used to approximate the conventional solution. It is a reminder that integration range of both integrals are the same from 0 to first dimensionless time At. The m{ and n{ coefficients are functions of the dimensionless time step At. Their values can be determined with use of the analytical formulas presented in [35]. The range of applicability of the filtered weighting function is from At up to 103 At. 2 COMPARISON OF NUMERICAL SOLUTION WITH EXPERIMENTAL RESULTS To assess the effects of pipe-wall viscoelasticity and unsteady friction on the transient flows, extensive comparative studies have been carried out by researchers. Experimental results obtained by the researchers from European scientific centres (London (UK), Cassino (Italy) and Lyon (France)) are first used to validate the numerical model described in the previous chapter of the paper. For convenience, the comparisons of all the results are presented in dimensionless form. Then the pressure changes observed experimentally (EXV) and simulated by using unsteady (Sv>uF) and quasi-steady hydraulic resistance (Sqsf), and the time are determined by the following formulas: P = P ( )- Pf _ P ( )- Pf and / = t / (4L / c). (10) AP Pcvo For quantitative analysis, it is essential to define a minimum of two coefficients, whose role is to mathematically determine the compliance of the simulated pressure histories with respect to the experimental tests. To the authors' knowledge, only few studies in the literature [24], [28], [36] and [37] in the field of water hammer flows are concerned on such quantitative analysis. Urbanowicz's analysis carried out in [28] was based on pressure histories in dimensional form. However, that analysis method may not be applicable to achieve the purpose of this study. Specifically, for laminar flows, the pressure pulsation from the value of the final pressure to which the transient state ends (equal to the pressure inside the reservoir) is small. Also, the pressure drop along the length of the pipe becomes lower as the Reynolds number is smaller. This is particularly noticeable in laminar flows with very small Reynolds numbers. In preliminary tests, the values of quantitative parameters determining the pressure compliance were calculated using the method presented in [28], which however resulted in very small values for laminar runs. Hence, in this work, to avoid the influence of the Reynolds number on this analysis, the quantitative analysis has been carried out in a dimensionless form. For this purpose, a subprogram was defined to search for maximum and minimum values and their occurrence times (calculated from the beginning of the analysed transient state). For demonstration, Fig. 1 illustrates the procedure of "collecting" maximal and minimal pressures, and their relevant times (circles). Fig. 1. Selecting the maximal and minimal pressures and their times i=i 80 Urbanowicz, K. - Duan, H.-F. - Bergant, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 On the one hand, the coefficient determining the compliance of the maximum and minimum simulated pressures is calculated as follows: IL ep= 100%, (11) k where psi is simulated maximal or minimal dimensionless pressure, and pe. is experimentally recorded maximal or minimal dimensionless pressure. On the other hand, the coefficient that determines the time compliance of subsequent simulated amplitudes was calculated as follows: t s,i t e T Et =- te k - 1 •100%, (12) where tsj is the simulated dimensionless time ^of maximal or minimal dimensionless pressure and te,i is the experimentally recorded dimensionless time of maximal or minimal dimensionless pressure. The degree of simulation compatibility and accuracy increases with decreasing values of the above coefficients (more accurate, matching of the simulated pressure runs with respect to the known experimental runs). Tables 1, 2 and 3 (presented later in this paper) summarize the quantitative analysis results of the Ep and Et parameters for all test cases in this work. The wave speeds were estimated from the empirically observed durations of the first pressure amplitude. Then the instantaneous creep coefficient J can be estimated by the formula: J = 1 1 pzc k;: (13) where S = — £ is expanded pipe constraint e coefficient. For Covas and Evangelista test stands, the pipe constraint coefficient is [3] to [5]: * = f(1 ) + ((14) and for Güney test stand [6] and [7]: 1 + f 1 I2 + V ^ -vp fi -1 D D D (15) In all experiments, the pressure inside the supply reservoir increased with different extents during the water hammer event. The largest increase was observed in the experimental runs carried out by the Evangelista et al. [38] in Casinno (Italy). A relatively smaller increase was recorded during tests carried out by Covas [16] and [39] in London (United Kingdom). Finally, the smallest increase in pressure, inside the supply reservoir, was noted in tests carried out by Guney [6] and [40] in Lyon (France). The main reason for the increase in pressure might be attributed to the poor synchronization of the rapid closure of the valve with the shut-off of the working pump. In the pre-transient state the reservoir pressure prt= 0 is the sum of the theoretical pressure loss over the pipe length Api and the pressure value recorded in the valve cross-section pvt= o: Pr ,t=0 = Pv,t=0 + fLp v2 D 2 = Pv, -APL (16) where Darcy-Weisbach friction factor for laminar flow is /= 64 / Re, and for turbulent flow (assuming smooth pipe walls) f = ^'^i6*4. With use of the above-calculated reservoir pressure at the start of transients prt = o , and known final pressure inside the reservoir ppinai, the value of the pressure increase in the cross-section at the reservoir Apr can be estimated approximately by: so that APr = PFinal ~ Pr ,t P Final =APr + APl + Pvt=0 (17) (18) For all Covas cases with known overpressures, there is no information about valve pressure pvt = o in steady flow prior to transient. Based on the residual information available from their works [16] and [39], the values of pKi= o for Covas test stand were assumed to achieve the simulations (see Table 1). For all test systems, the information needed to conduct simulations (pipe length, inner diameter, wall thickness, the Poisson's ratio, etc.) is collected in Tables 1, 2 and 3. The creep function used for Covas test stand simulations (water temperature T= 20 °C), which was needed to calculate the retarded strain, was archived from experimental test [16], by giving the following creep and retardation time coefficients: J0 = 0.674-10-? Pa-1, J1 = 0.1394-10-> Pa-1, J2 = 0.0062-10-9 Pa-1, J3 = 0.1148-10-9 Pa-1, J4 = 0.3425 10-9 Pa-1, J5 = 0.0928-10-9 Pa-1 and Rt0 = 0 s, Rt1 = 0.05 s, Rt2 = 0.5 s, Rt3 = 1.5 s, Rt4 = 5 s, Rt5 = 10 s. The parameters of creep compliance functions vary strongly with temperature. Especially for Evangelista et al. test stand (water temperature T= 15 °C), the above function has to be scaled. The scaling factor was assumed, as in work [41] to 2 Transient Liquid Flow in Plastic Pipes 81 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 be 0.75 which gives the following values for the simulations: J0 = 0.513-10-9 Pa-1, J = 0.1046-10-9 Pa-i, J2 = 0.0046-10-9 Pa-1, J3 = 0.0861-10-9 Pa-1, J4 = 0.2569-10-9 Pa-1, J5 = 0.0696-10-9 Pa-1 and Rt1 = 0.05 s, Rt2 = 0.5 s, Rt3 = 1.5 s, Rt4 = 5 s, Rt5 = 10 s. The creep function for Guney test stand was experimentally defined for different temperatures. Original values received by Guney are required to be filtered according to the filtration procedure proposed by Urbanowicz and Firkowski [33]. The obtained Table 1. Covas test stand details and quantitative results Case vo [m/s] Reo [-] pv,t= 0 [Pa] Api [Pa] APr [Pa] Unsteady friction (UF) [%] Quasi-steady friction (QSF) [%] AEP AEt aep,uf AEt,UF AEp,QSF AEt,QSF [%] [%] C01 0.028 1417 4.735e5 0.01e4 0.14e4 5.45 0.45 32.68 1.70 27.23 1.25 C02 0.124 6274 4.405e5 0.15e4 0.45e4 2.03 0.97 21.43 1.66 19.40 0.69 C03 0.249 12599 4.405e5 0.50e4 1.00e4 2.77 1.20 15.86 1.87 13.09 0.67 C04 0.373 18874 4.405e5 1.01e4 1.39e4 2.01 0.45 14.94 1.53 12.93 1.08 C05 0.500 25300 4.289e5 1.65e4 2.00e4 3.93 0.29 14.40 1.20 10.47 0.91 C06 0.746 37748 3.425e5 3.39e4 3.11e4 3.04 0.58 12.05 1.30 9.01 0.72 C07 0.870 44022 3.425e5 4.43e4 3.30e4 3.41 0.46 11.21 1.33 7.80 0.87 C08 0.995 50347 4.735e5 0.01e4 0.14e4 3.70 0.56 11.58 1.17 7.88 0.61 L = 271.7 m; D= 0.0506 m; e = 0.0063 m; vp= 0.46 ; Ç = 1.065; S = 8.551; T = 20 °C; p = 998.1 kg/m3; v = 10-6 m2/s; c= 400 m/s Table 2. Evangelista test stand details and quantitative results Case vo [m/s] Reo [-] Pv,t= 0 [Pa] apl [Pa] APr [Pa] Unsteady friction (UF) [%] Quasi-steady friction (QSF) [%] AE AEt [%] aep,uf AEt,UF AEp,QSF AEt,QSF [%] E01 0.070 2702 4.198e5 496.5 3950 7.80 3.28 34.98 4.34 27.18 1.06 E02 0.098 3783 4.183e5 894.5 5160 4.47 3.15 33.09 4.36 28.62 1.21 E03 0.164 6330 4.123e5 2202 8500 3.67 2.76 27.60 4.19 23.93 1.43 E04 0.335 12930 3.992e5 7686 16600 2.80 2.78 20.68 4.10 17.88 1.32 E05 0.657 25358 3.668e5 2.5e4 16800 5.15 2.86 14.94 4.37 9.79 1.51 E06 0.989 38172 3.235e5 5.11e4 48400 5.16 2.75 12.42 4.19 7.26 1.44 E07 1.315 50754 2.700e5 8.414e4 67860 5.16 2.22 10.60 3.76 5.44 1.54 E08 1.645 63491 2.045e5 1.245e5 93000 5.78 2.39 8.61 3.88 2.83 1.49 L = 203.3 m; D= 0.044 m; e=0.003 m; vp = 0.46 ; Ç = 0.937; S = 13.745; T= 15 °C; p = 999.1 kg/m3; v= 1.1410-6 m2/s; c = 365 m/s Table 3. Güney test stand details and quantitative results Case vo [m/s] Reo [-] Pv,t= 0 [Pa] Api [Pa] APr [Pa] Unsteady friction (UF) [%] Quasi-steady friction (QSF) [%] AE AEt [%] aep,uf AEt,UF AEp,QSF AEt,QSF [%] G01 0.490 17422 101325 3423 4750 3.54 3.73 14.28 5.94 10.74 2.21 G02 0.550 25650 101325 3906 5000 8.98 2.77 19.57 4.15 10.59 1.38 G03 0.570 30245 101325 4019 10000 8.81 1.36 19.73 1.87 10.92 0.51 G04 0.550 31646 101325 3695 10000 3.40 2.35 11.91 3.54 8.51 1.19 G05 0.560 34513 101325 3743 8000 9.89 2.56 19.39 3.77 9.50 1.21 G06 0.820 50536 101325 7296 10000 8.84 2.32 16.82 2.94 7.98 0.62 L = 43.1 m; D= 0.0416 m; e=0.0042 m; vp = 0.38 ; Ç = 0.97; S = 9.61 Table 4. Values of temperature dependent Güney's research coefficients Case T[°C] c [m/s] p [kg/m3] v [m2/s] Ki [Pa] Ai [s] Jo [Pa-1] Ji [Pa-1] J2 [Pa-1] Rti [s] Rt2 [s] G01 13.8 310 999.3 1.1710-6 2.14109 0.0043 1.035 10-9 0.637 10-9 0.87110-9 0.0166 1.747 G02 25 265 997.1 0.892 10-6 2.24109 0.0051 1.440 10-9 1.046 10-9 1.237 10-9 0.0222 1.864 G03 31 247 995.3 0.784 10-6 2.27109 0.0055 1.668 10-9 1.397 10-9 1.628 1 0-9 0.0221 1.822 G04 35 240 994.1 0.723 10-6 2.285 109 0.0056 1.772 1 0-9 1.797 10-9 2.349 10-9 0.0265 2.392 G05 G06 38.5 220 992.6 0.675 1 0-6 2.295 109 0.0061 2.12110-9 2.097 10-9 3.570 10-9 0.0347 3.077 82 Urbanowicz, K. - Duan, H.-F. - Bergant, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 values and other temperature-dependent parameters details are collected in Table 4. In all MOC-based simulations, a constant number N= 32 of reaches was used. Such a pipe division gave the following time steps: for Covas test stand At= 0.0212 s, for Evangelista test stand At = 0.0174 s and for Guney test stand, detailed values can be found in Table 4. A rapid-closing valve (tcl < 2L/c) at the downstream end of the pipes was applied to all scenarios (Covas, Evangelista and Guney test stands) to generate water hammer flows in the system. In all of the three test stands, two distinct flow cases were selected (marked with gray in Tables 1 to 3) and investigated in detail. 2.1 Research Based on Covas Experimental Results The results of the selected simulation cases C01 (laminar flow with initial: velocity vo = 0.028 m/s and Reynolds number Re0 = 1410) and C05 (turbulent flow with initial: velocity V) = 0.5 m/s and Reynolds number Re0 = 25300) in relation to the Covas experimental data are plotted in Figs. 2 and 3, respectively. It can be seen from Fig. 2a that the simulation runs realized for laminar flow (C01), taking into account in a simplified "step way" that the recorded pressure increases inside the supply reservoir, the model results present a high accuracy. Fig. 2b presents the comparison of the same experimental data with simulation results, but with implementing only the quasi-steady hydraulic resistance t = Tq in the simulation. It can be observed that in this case the modelled damping (amplitude attenuation) was much smaller, which results in overstated simulation runs (Svqsf). Moreover, comparison in Fig. 2 reveals that a change in the method of modelling the wall shear stress could greatly affect the shift of simulated pressure runs - the use of unsteady friction delays the moment of occurrence of subsequent amplitudes. This means that unsteady friction indirectly affects the change (during the transients) of the pressure wave speed. a) b) Fig. 2. Comparison of simulated and experimental results for Covas C01 (laminar flow); a) unsteady friction model, and b) quasi-steady friction model a) 0.6 0.4 0.2 £ 0 -0,2 -0.4 -0.6 1 b) if-] Fig. 3. Comparison of simulated and experimental results for Covas C05 (turbulent flow); a) unsteady friction model, and b) quasi-steady friction model Transient Liquid Flow in Plastic Pipes 83 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 From the analysis of Fig. 3 obtained for the C05 for turbulent flow scenario, most of the results are similar to those just discussed for laminar case C01. Of the interesting noticeable differences between Figs. 2a and 3a showing the normalized Covas cases analysed in detail in this work (laminar C01 and turbulent C05), the shape of the first amplitude deserves additional attention. In the laminar flow case (C01), it can be found that the maximum value occurred in the initial period, i.e., just after the valve was closed, there was a significant pressure drop at the peak of the first amplitude (see Fig. 4). A similar tendency was also noted from the observation of experimental results. For the turbulent flow case in Fig. 3, it can be seen that there was only a slight decrease in pressure to the minimum value at the top of first amplitude, which maintained a constant value for a long time (normalized times from 0.17 up to 0.35), and, in the final phase, there was a slight increase in pressure. In Fig. 3, the pressure runs at the first amplitude top is almost flat. For clarity, an enlarged plot for this first amplitude is shown in Fig. 4 for a detailed comparison. It is suspected that the above change in pressure behaviour at the top of the first amplitudes is the result of a change in the relationship between two simulated phenomena: friction and retarded strain (both described by convolutional integrals). Fig. 4. Auxiliary drawing - the top of the first pressure amplitude From the quantitative results collected in Table 1 regarding experimental tests carried out by Covas, it can be concluded that: • pressure runs simulated with applying unsteady friction presented a much better accuracy than those in which the resistances were modelled with a quasi-steady component only. This is because for all test cases (C01 up to C08), an inclusion of unsteady friction may reduce the values of Ep and Et coefficients; • the average values of Ep and Et parameters in the test runs, taking into account unsteady resistances, were EpUFmieanC = 3.29 % and Et,UF,meanC = 0.62 %, respectively. However, for test runs simulated with quasi-steady resistances only, EpQSFmeanC = 16.77 % and EtQSFmeanC = 1.47 %, respectively; • the highest Ep value (means the worst simulation fit) was obtained for the laminar flow case. For the model with unsteady friction, the errors were: EpurF,c01 = 5.45 % and EtuF,C01 = 0.45 %, while in the model with quasi-steady resistance: Ep,QSF,C01 = 32.68 % and E^Fm = 1.7 %. The possible reason of these poor results can be attributed to a large pressure signal noise (before filtering) noticed here for laminar results (the same problem also affected other tests with a low initial flow rate (Re < 6500) carried out by Evangelista et al.; see Table 2). 2.2 Research Based on the Experimental Results of Evangelista et al. Two cases for the Evangelista test stand are analysed herein, in which the main difference is the initial velocity of the flowing liquid (Reynolds number). In the first case (E01) (transitional flow case), the initial flow velocity was v0 , E01 = 0.070 m/s, and in the second analysed case (E08) (turbulent case) it was V0E08 = 1.64 5 m/s. The initial velocity significantly affects not only the Darcy friction factor and unsteady friction calculation but also the pressure pulsations appearing after the rapid valve closures. Meanwhile, what follows from equations presented in the first chapter, the simulated retarded strains are also affected (which are usually treated by researchers as a source of additional damping of flow pulsations). Experimental results (especially for low flow rates) were characterized by high noise. For this purpose, the MatLab "filtfilt" filter was used, which performs zero-phase digital filtering by processing the input data in both the forward and reverse directions. After filtering the data in the forward direction, the "filtfilt"function reverses the filtered sequence and runs it back through the filter. The result has the following characteristics: zero phase distortion; a filter transfer function equal to the squared magnitude of the original filter transfer function. The results of selected simulation tests (E01 and E08) are presented in Figs. 5 and 6, respectively. From Fig. 5a, it can be seen that the degree of matching of simulation results (using unsteady friction) to experimental results seems not perfect. This difference may be attributed to the following factors. 84 Urbanowicz, K. - Duan, H.-F. - Bergant, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 First, the Reynolds number Re = 2702 indicates that it is a transitional flow between laminar and turbulent. In these flows, according to experimental research related estimation of friction coefficients [42], there is a sharp change in the value of this coefficient. In the simulation, the friction coefficient and the unsteady hydraulic resistances were calculated in a simplified manner, assuming the laminar nature of the flow (f= 64/Re and m, n representing the weighting function were used for laminar flow). Secondly, the experimental results were initially characterized by high noise, as well as visible effects of reflected waves (which may be the result of narrowing at pressure sensor connection points, sharp knees, intense a) b) Fig. 5. Comparison of simulated and experimental results for Evangelista E01 (transitional flow); a) unsteady friction model, b) quasi-steady friction model a) b) Fig. 6. Comparison of simulated and experimental results for Evangelista E08 (turbulent flow); a) unsteady friction model, b) quasi-steady friction model a) b) Fig. 7. Enlargement of: a) top of the first amplitude, and b) first and second valley of the course of transient pressure changes Transient Liquid Flow in Plastic Pipes 85 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 vibration of the pipeline, etc.) on pressure courses, referring to the sharp increase at the middle of the first amplitude (Fig. 7a) as well as peculiar increase visible on the first and during subsequent valleys of the pressure course (Fig. 7b). The quality of matching of the simulation and experimental runs may be considered satisfactory. From the results in Table 2, describing the quantitative results of Evangelista scenarios, it is evident that transient test runs modelled with additional unsteady friction could match better with experiment than those modelled using the quasi-steady friction only. Along with the increase in the Reynolds number (from case E02 represented with the Reynolds number Re= 3783), a significant decrease in the significance of unsteady friction on the maximum and minimums values of pressure histories was noticed in this system, and evidenced by the decrease of such difference of pressure compliance parameters: AEp = Ep,QSF - Ep Up. (19) For example, in the case E08 in which the highest flow rate was recorded, the difference AEp was only A¿p,£08 = 2.83 %. However, the unsteady friction should not be neglected even for large Reynolds numbers, because it affects the time compliance of simulated transient runs, by referring to the similar difference value of: AEt = Et,QsF - Et,uF , (20) for different Reynolds numbers as shown in Tables 1 and 2. 2.3 Research Based on Güney Experimental Results For Güney test stand six cases of turbulent flow were analysed, in which the main difference was the water temperature. Only two of these tested cases are selected for detailed discussion in this paper. It is known that temperature may significantly affect the mechanical properties of the plastic pipe material. With such influence, a large change in the creep function may occur during the transient tests (Table 4). The creep function derivative is the main function affecting the modeled damping of pressure oscillation, and in this simulation it was filtered according to filtration procedure proposed by Urbanowicz and Firkowski [33]. The numerical simulation results of the experimental tests are provided in Figs. 8 and 9. From the results of the test G01 for relatively low temperature T = 13.8 °C, it is noticed that the model including unsteady friction (Fig. 8a) may match with 86 Urbanowicz, K. - Duan high accuracy the top of the first amplitude (Fig. 10), and even the peak of pressure on the top of this first amplitude could be reflected in the simulation process. This confirms the acceptable accuracy for the used mathematical model. As concluded earlier in this work, the results with better matching accuracy were obtained after applying unsteady friction (Fig. 8a) than using the quasi-steady resistance only (Fig. 8b). From the results in Fig. 9, it is also observed that an increase in temperature in the same system may cause a decrease in the number of appearing amplitudes in the same time interval. For example, when the temperature was 13.8 °C within two seconds after the rapid closing of the valve, there were three pressure amplitudes; while for the temperature of 38.5 °C, only two amplitudes are noticed at the same time. This is mainly due to the following two factors. The first is the change of the wave speed c, from CG01 = 310 m/s for temperature T = 13.8 °C to c^06 = 220 m/s for temperature T = 38.5 °C. The second factor may be attributed to the change of the creep function coefficients J and t{ (see Table 4). The creep parameters describing the properties of the material may strongly depend on the prevailing flow temperature (as seen in Fig. 11 illustrating the effect of temperature on creep functions). In the test case with relatively high temperature (Fig. 9), a large divergence of the modelled peak was observed just after the valve was closed (Fig. 10). It is difficult to clearly indicate the key reason for this. In this example, it is also difficult to confirm the positive effect of unsteady friction by analysing these results only (Fig. 9). Normalized pressure runs show only a slight advantage of the unsteady model over the quasi-steady one. Therefore, it is necessary to conduct a quantitative analysis on such results so as to determine the key factors affecting the modelling accuracy. From the results in Table 3 that were obtained for quantitative analysis for the Guney experimental system, the average value of the obtained coefficients illustrating model compliance was much higher (Ep,uF,meanG = 7.24 % and EtUFmeanG = 2.52 %) than for the other two quantified systems. This may be due to two aspects as follows. Firstly, these results were tested and collected a relatively long time ago (i.e., in 1977), which forced the use of pressure sensors with much larger ranges of measurement errors than those currently used. Secondly, the process of digitizing these results, carried out using the charts attached to the PhD thesis [6], could have also contributed to the erroneous reading of real experimental results. These charts were often characterized by casual information on the value of the ordinate axes. Even though these , H.-F. - Bergant, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 a) 1 L-J b) c l-J Fig. 8. Comparison of simulated and experimental results for Guney G01 (T = 13.8 Re = 17422); a) unsteady friction model, b) quasi-steady friction model 0.6 0,4 0.2 ^ 0 -0,2 -0.4 -0.6 Lf* « U \l Sv.QSF - - EX V r\ I \ L Ï /A 1 \ r \ \ I V / ^ If V / \ /' V A J ^ 0 12 3 aj * [") b) *H Fig. 9. Comparison of simulated and experimental results for Guney G06 (T = 38.5 °C, Re = 50536); a) unsteady friction model, b) quasi-steady friction model 0 0.1 0,2 03 0.4 0.5 0,6 i H Fig. 10. Enlargement of top of the first amplitudes 10"4 10"3 10"2 10"1 10° 101 102 103 t[ s] Fig. 11. Courses of the creep functions used in the simulations differences in this system, it can still be confirmed that the use of unsteady friction may result in a significant decrease of quantitative parameters, and thus the substantial increase of the simulation accuracy. Based on all quantitative analyses of the results for these three considered test systems, it can be concluded that the simulations performed using the experimentally defined creep functions by including additional unsteady friction effect may provide more accurate modelling results for transient flows in plastic pipes. However, it is also revealed that the lack of experimental data regarding the temperature effect on viscoelastic creep function is usually an issue that requires solutions in practical applications, because Transient Liquid Flow in Plastic Pipes 15 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 from this study the temperature is found to be an important factor affecting transient flow behaviour in plastic pipes. Furthermore, the proposed parameters of quantitative analysis in Eq. (11) and Eq. (12) in this study can be useful measures for characterizing and evaluating the importance of unsteady friction effect to transient flows with the existence of the pipe-wall viscoelasticity effect in plastic pipes. 3 DEFINING THE EFFECTIVENESS OF NEW PROPOSED METHOD Further verification of the effectiveness of a computer program is performed herein for the newly developed simple mathematical formulas in this study. As an illustrative example, the laminar Covas C01 run was chosen (details of this case can be found in Table 1), in which the effect of discretization was examined. Computation times told obtained using the original calculation algorithm [30] (complicated solution of retarded strain convolution integral, and unfiltered long exponential weighting function consisting 26 terms is used), were compared with calculation times obtained with the help of new procedures discussed in this work tnpw Table 5. Calculation details N [-] 16 32 64 128 256 512 told [s] 0.30 1.82 11.65 157 1152 10272 -new [s] 0.11 0.52 3.79 66 487 4377 Nodes 469 x17 940 x33 1882 x65 3766 x129 7535 x257 15073 x513 At [e-6] 66.3 33.2 16.6 8.29 4.15 2.07 The tested reservoir-pipe-valve (R-P-V) system was divided into a predefined number of reaches N in accordance with Table 5. An increase in discretization corresponds to the corresponding intensity of the calculation (growing number of computational nodes). The calculation time was measured using a MatLab functions "tic" and "toe". The "tic" command was introduced in the software after specifying the boundary and initial conditions, just before the loop function in which the transient values are calculated. The "toe" command was introduced in the computer program just after the loop function. Simulations were repeated using selected discretization for both calculation algorithms with at least five times. The table summarizes the minimal calculation times obtained. Calculations were carried out on a typical laptop for personal use. It can be seen from Table 5 that the new procedure reduces the calculation time by more than two times (or even three times), in the case of the discretization commonly used in practice, i.e., N< 64. 4 CONCLUSIONS The main results and findings of this paper can be summarized as follows: (1) The comparative studies have demonstrated the significant role of unsteady friction effects in plastic pipes during rapid transient events. These effects could usually be masked in many previous research studies during their creep function calibrations. (2) There is a strong need to conduct detailed experimental research on polymers used for pressure pipes, aiming to determine the exact temperature-dependent creep functions. Such tests cannot be carried out only on individual dynamical mechanical thermal analyser (DMTA) devices, because the frequency of forcing deformations of the tested material may turn out to be too low compared to those factors that affect pipes in real conditions, i.e., systems in which water hammer occurs. A detailed analysis of these creep function courses is required immediately after applying the step load. (3) A new experimental research is needed in straight and horizontal R-P-V (reservoir-pipe-valve) system, in which the local resistances can be minimized. This is because local resistances can significantly affect the flow, as shown by Stosiak et al. [43]. To this end, it is necessary to design pressure sensors working exactly inside the shut-off valve and to use the optimal rounding radius of the plastic pipe in the input cross-section located on the wall of the pressure reservoir. 5 ACKNOWLEDGEMENTS The authors would like to thank Professor Angelo Leopardi from University of Cassino and Southern Lazio for providing the experimental results to one of the tested scenarios. The co-author Duan would like to thank the Hong Kong Research Grants Council (RGC) for the support on his research work (15201017 and 15200719). Finally, Bergant gratefully acknowledges the support of the Slovenian Research Agency (ARRS) conducted through the research project L2-1825 and the programme P2-0162. 88 Urbanowicz, K. - Duan, H.-F. - Bergant, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 6 NOMENCLATURES c pressure wave speed, [m/s] D inner diameter of the pipe, [m] e thickness of the pipe wall, [m] E Young modulus of pipe material, [Pa] f Darcy-Weisbach friction factor, [-] J(T) compliance creep function, [Pa-1] J instantaneous creep coefficient, [Pa-1] J creep-compliance of the ith spring of the KelvinVoigt model, [Pa-1] K bulk modulus of liquid, [Pa] N number of reaches, [-] p liquid pressure, [Pa] Rt retardation times, [s] t time, [s] tcl valve closing time, [s] T temperature, [°C] u time used in convolution integral, [s] v liquid flow velocity, [m/s] W(t_u) friction weighting function,[-] wj(t-u) creep weighting function, [Pa-1s-1] x distance along the pipe axis, [m] se instantaneous elastic strain, [-] sr retarded strain, [-] H dynamic viscosity of liquid, [Pa-s] v kinematic viscosity, [m2/s] £ pipe support coefficient, [-] S expanded pipe support coefficient, [-] p liquid density, [kg/m3] a circumferential stress, [Pa] t wall shear stress, [Pa] At constant time step in the MOC, [s] 7 REFERENCES [1] Tison, G. (1958). Le mouvement non permanent succédant à l'ouverture d'une vanne sur une conduite en polyethylene (The unsteady flow produced by the opening of a valve on a polyethylene pipe). Becetel, communication N°3, Ghent University, Ghent. (in French) [2] Hardung, V. (1962). Propagation of pulse waves in visco-elastic tubings. Handbook of Physiology, (Field, J., ed.) p. 107135, American Physiological Society, Washington. [3] Martin, T. (1969). Pulsierende Strömung in visko-elastischen Leitungssystemen. Ingenieur-Archiv, vol. 37, no. 5, p. 315325, D0l:10.1007/BF00532579. (in German) [4] Lorenz, J., Zeller, H. (1972) An analogous treatment of waves propagation in liquid filled elastic tubes and gaz filled rigid tubes. Proceedings of the 1st International Conference on Pressure Surges, Paper B5. [5] Rieutord, E., Blanchard, A. (1972). Influence d'un comportement viscoélastique de la conduite dans le phénomène du coup de bélier (Influence of a viscoelastic pipe behavior in the phenomenon of water hammer). Reports of the Academy of Sciences, Paris, vol. 274, p. 1963-1966. (in French) [6] Güney, M.S. (1977). Contribution à l'étude du phénomène de coup de bélier en conduite viscoélastique (Contribution to the study of the water hammer phenomenon in viscoelastic pipe), (Phd Thesis). Claude Bernard University, Lyon. (in French) [7] Zielke, W. (1968). Frequency-dependent friction in transient pipe flow. Journal of Basic Engineering, vol. 90, no. 1, p. 109115, D0l:10.1115/1.3605049. [8] Rieutord, E., Blanchard, A. (1979). Pulsating viscoelastic pipe flow - water-hammer. Journal of Hydraulic Research, vol. 17, no. 3, , p. 217-229, D0I:10.1080/00221687909499585. [9] Limmer, J., Meißner, E. (1974). Druckwellen in visko-elstischen Leitungen. In Mitteilungen des Institutes für Hydraulik und Gewässerkunde, Technische Universität München, München, vol. 14. (in German) [10] Franke, P.G., Seyler, F. (1983) Computation of unsteady pipe flow with respect to visco-elastic material properties. Journal of Hydraulic Research, vol. 21, no. 5, p. 345-353, D0I:10.1080/00221688309499456. [11] Ghilardi, P., Paoletti, A. (1987). Viscoelastic parameters for the simulation of hydraulic transients in polymeric pipes. Excerpta, vol. 2, p. 51-62. [12] Pezzinga, G., Scandura, P. (1995). Unsteady flow in installations with polymeric additional pipe. Journal of Hydraulic Engineering, vol. 121, no. 11, p. 802-811, D0I:10.1061/(ASCE)0733-9429(1995)121:11(802). [13] Tijsseling, A.S. (1996). Fluid-structure interaction in liquid-filled pipe systems: a review. Journal of Fluids and Structures, vol. 10, no. 2, p. 109-146, D0I:10.1006/jfls.1996.0009. [14] Covas, D., Stoianov, I., Ramos, H., Graham, N., Maksimovic, C., Butler, D. (2004). Water hammer in pressurized polyethylene pipes: conceptual model and experimental analysis. Urban Water Journal, vol. 1, no. 2, p. 177-197, D0I:10.1080/15730 620412331289977. [15] Covas, D., Stoianov, I., Ramos, H., Graham, N., Maksimovic, C. (2004). The dynamic effect of pipe-wall viscoelasticity in hydraulic transients. Part I - experimental analysis and creep characterization. Journal of Hydraulic Research, vol. 42, no. 5, p. 516-530, D0I:10.1080/00221686.2004.9641221. [16] Covas, D., Stoianov, I., Ramos, H., Graham, N., Maksimovic, C. (2005). The dynamic effect of pipe-wall viscoelasticity in hydraulic transients. Part II - model development, calibration and verification. Journal of Hydraulic Research, vol. 43, no. 1, p. 56-70, D0I:10.1080/00221680509500111. [17] Duan, H.-F., Ghidaoui, M., Lee, P.J., Tung, Y.-K. (2010). Unsteady friction and visco-elasticity in pipe fluid transients. Journal of Hydraulic Research, vol. 48, no. 3, p. 354-362, D0I:10.1080/00221681003726247. [18] Brunone, B., Berni, A. (2010). Wall shear stress in transient turbulent pipe flow by local velocity measurement. ASCE Journal of Hydraulic Engineering, vol. 136, no. 10, p. 716-726, D0I:10.1080/00221681003726247. [19] Bergant, A., Hou, Q., Keramat, A., Tijsseling, A.S. (2013). Waterhammer tests in a long PVC pipeline with short steel end sections. Journal of Hydraulic Structures, vol. 1, no. 1, p. 2334, D0I:10.22055/JHS.2013.10069. Transient Liquid Flow in Plastic Pipes 89 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 [20] Keramat, A., Tijsseling, A.S., Hou, Q., Ahmadi, A. (2012). Fluid-structure interaction with pipe-wall viscoelasticity during waterhammer. Journal of Fluids and Structures, vol. 28, p. 434-455, DOI:10.1016/j.jfluidstructs.2011.11.001. [21] Pezzinga, G., Brunone, B., Meniconi, S. (2016). Relevance of pipe period on Kelvin-Voigt viscoelastic parameters: 1D and 2D inverse transient analysis. Journal of Hydraulic Engineering, vol. 142, no. 12, DOI:10.1061/(ASCE)HY.1943-7900.0001216. [22] Meniconi, S., Brunone, B., Ferrante, M., Massari, C. (2014). Energy dissipation and pressure decay during transients in viscoelastic pipes with an in-line valve. Journal of Fluids and Structures, vol. 45, p. 235-249, DOI:10.1016/j. jfluidstructs.2013.12.013. [23] Weinerowska-Bords, K. (2015). Alternative approach to convolution term of viscoelasticity in equations of unsteady pipe flow. Journal of Fluids Engineering, vol. 137, no. 5, art. ID: 054501, DOI:10.1115/1.4029573. [24] Ferrante, M., Capponi, C. (2017). Comparison of viscoelastic models with a different number of parameters for transient simulations. Journal of Hydroinformatics, vol. 20, no. 1, p. 1-17, DOI:10.2166/hydro.2017.116. [25] Kodura, A. (2016). An analysis of the impact of valve closure time on the course of water hammer. Archives of HydroEngineering and Environmental Mechanics, vol. 63, no. 1, p. 35-45, DOI:10.1515/heem-2016-0003. [26] Ferrante, M., Capponi, C. (2017). Viscoelastic models for the simulation of transients in polymeric pipes. Journal of Hydraulic Research, vol. 55, no. 5, p. 599-612, DOI:10.1080/ 00221686.2017.1354935. [27] Stryczek, J., Banas, M., Krawczyk, J., Marciniak, L., Stryczek, P. (2017). The fluid power elements and systems made of plastics. Procedia Engineering, vol. 176, p. 600-609, DOI:10.1016/j.proeng.2017.02.303. [28] Urbanowicz, K. (2018). Fast and accurate modelling of frictional transient pipe flow. Zeitschrift für Angewandte Mathematik und Mechanik, vol. 98, no. 5, p. 802-823, DOI:10.1002/zamm.201600246. [29] Vardy, A.E., Brown, J.M.B. (2010). Evaluation of unsteady wall shear stress by Zielke's method. Journal of Hydraulic Engineering, vol. 136, p. 453-456, DOI:10.1061/(ASCE) HY.1943-7900.0000192. [30] Urbanowicz, K., Firkowski, M., Zarzycki, Z. (2016). Modelling water hammer in viscoelastic pipelines: short brief. Journal of Physics: Conference Series, vol. 760, no. 1, art. ID: 012037, DOI:10.1088/1742-6596/760/1/012037. [31] Schohl, G.A. (1993). Improved approximate method for simulating frequency - dependent friction in transient laminar flow. Journal of Fluids Engineering, vol. 115, no. 3, p. 420424, DOI:10.1115/1.2910155. [32] Urbanowicz, K., Flrkowskl, M. (2018) Modelling water hammer with quasi-steady and unsteady friction in viscoelastic pipelines. Dynamical Systems in Applications, Awrejcewicz J. (ed.) Springer Proceedings in Mathematics & Statistics, vol. 249, p. 385-399, Springer, Cham, DOI:10.1007/978-3-319-96601-4_35. [33] Urbanowicz, K., Firkowski, M. (2018). Effect of creep compliance derivative in modeling water hammer in viscoelastic pipes. Proceedings of 13th International Conference Pressure Surges, p. 305-324. [34] Vardy, A.E., Brown, J.M.B. (2003). Transient turbulent friction in smooth pipe flows. Journal of Sound and Vibration, vol. 259, no. 5, p. 1011-1036, DOI:10.1006/jsvi.2002.5160. [35] Urbanowicz, K. (2017). Analytical expressions for effective weighting functions used during simulations of water hammer. Journal of Theoretical and Applied Mechanics, vol. 55, no. 3, p. 1029-1040, DOI:10.15632/jtam-pl.55.3.1029. [36] Adamkowski, A., Lewandowski, M. (2006). Experimental examination of unsteady friction models for transient pipe flow simulation. Journal of Fluids Engineering, vol. 128, no. 6, p. 1351-1363, DOI:10.1115/1.2354521. [37] Bertaglia, G., loriatti, M., Valiani, A., Dumbser, M., Caleffi, V. (2018). Numerical methods for hydraulic transients in viscoelastic pipes. Journal of Fluids and Structures, vol. 81, p. 230254, DOI:10.1016/j.jfluidstructs.2018.05.004. [38] Evangelista, S., Leopardi, A., Pignatelli, R., de Marinis, G. (2015). Hydraulic transients in viscoelastic branched pipelines. Journal of Hydraulic Engineering, vol. 141, no. 8, art. lD 04015016, DOI:10.1061/(ASCE)HY.1943-7900.0001030. [39] Covas, D. (2003). Inverse Transient Analysis for Leak Detection and Calibration of Water Pipe Systems Modelling Special Dynamic Effects. PhD Thesis, University of London, London. [40] Gally, M., Guney, M., Rieutord, E. (1979). An investigation of pressure transients in viscoelastic pipes. ASME Journal of Fluids Engineering, vol. 101, no. 4, p. 495-499, DOI:10.1115/1.3449017. [41] Firkowski, M., Urbanowicz, K., Duan, H.F. (2019). Simulation of unsteady flow in viscoelastic pipes. SN Applied Sciences, vol. 1, no. 519, DOI:10.1007/s42452-019-0524-2. [42] Swanson, C.J., Julian, B., Ihas, G.G., Donnelly, R.J. (2002). Pipe flow measurements over a wide range of Reynolds numbers using liquid helium and various gases. Journal of Fluid Mechanics, vol. 461, p. 51-60, DOI:10.1017/ S0022112002008595. [43] Stosiak, M., Zawislak, M., Nishta, B. (2018). Studies of resistances of natural liquid flow in helical and curved pipes. Polish Maritime Research, vol. 25, no. 3, p. 123-130, DOI:10.2478/pomr-2018-0103. 90 Urbanowicz, K. - Duan, H.-F. - Bergant, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 91-104 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2019.6298 Original Scientific Paper Received for review: 2019-08-23 Received revised form: 2019-11-13 Accepted for publication: 2019-12-09 Uneven Load Contact Dynamic Modelling and Transmission Error Analysis of a 2K-V Reducer with Eccentricity Excitation Guangjian Wang - Li Su1,3 - Shuaidong Zou1,3 1 Chongqing University, State Key Laboratory of Mechanical Transmissions, China 2 Chongqing University, School of Automotive Engineering, China 3 Chongqing University, College of Mechanical Engineering, China In this paper, the formula of non-load transmission error (TE) and time-varying backlash is derived through the incremental meshing line method based on the eccentricity error model of 2K-Vgear pair with small tooth difference. Considering the influence of uneven load contact between two gear pairs with small tooth difference, the static transmission error (STE) of a 2K-V gear pair with small tooth difference is deduced. This work helps to predict the TE caused by eccentricity error and provide inference to compensate backlash. Also, we establish the nonlinear dynamic model of 2K-V small tooth difference under eccentricity excitation, in which non-linear factors such as non-uniform load, time-varying backlash, transmission error and contact stiffness are considered, and validated the theoretical model by comparing the dynamic transmission error in simulation and in theory. Furthermore, the transmission error of the 2K-V gear reducer was tested, and the test results further verified the uneven load contact phenomenon of the gear pair with small tooth difference in the meshing process. At the same time, the experimental transmission error frequency spectrum is consistent with the theoretical results, indicating that the transmission accuracy of the 2K-V reducer is mainly caused by the eccentricity error of the gear pair with small tooth difference. Finally, we optimized the initial eccentricity phase and predicted the minimum transmission error, providing theoretical guidance for the assembly of 2K-V small tooth difference gear pairs. Keywords: 2K-V reducer, nonlinear dynamics, transmission error, uneven load contact, eccentricity error Highlights • The formula of the non-load transmission error and time-varying backlash based on the eccentricity error model of a 2K-Vgear pair with small tooth difference is derived. • The static transmission error of a 2K-V gear pair with small tooth difference involving the influence of uneven load contact is deduced. • The nonlinear dynamic model of a 2K-V small tooth difference under eccentricity excitation is established. • The simulation and experiment further validate the dynamic model and theory formula. 0 INTRODUCTION At present, the quality, volume, reliability, and cost of 2K-V reducer have been required higher performance in the aerospace, robotics, and other fields. The 2K-V reducer which is applicable for precision transmission (high-speed input, high-torque output as well as limited volume and mass) has the advantages of large speed ratio, high efficiency, compact structure and low processing cost. In the precision gear transmission system, the eccentricity error is the primary source of the time-varying transmission error in the major cycle and the periodical essential-varying backlash. Therefore, it is significant to study the transmission accuracy of the 2K-V reducer with eccentricity excitation [1]. In terms of the research on the transmission accuracy of the 2K-V reducer, Blanche and Yang [2] and Yang and Blanche [3] studied the transmission accuracy of the cycloidal mechanism with small tooth difference based on the geometry method and deduced the transmission ratio error that was caused by machining error and assembly error - calculation formula. Hidaka et al. [4], using the spring-mass equivalent model, analysed the transmission accuracy of 2K-V small tooth difference gear pairs. In order to increase modelling realism, Zhang et al. [5] established a translational-rotational coupled dynamic model that was linearly time-invariant to predict the vibration modes and natural frequencies. Xu et al. [6] also developed a dynamic model, which integrated gearbox body flexible supporting stiffness, to investigate how gearbox body flexibility influenced on the dynamic response and load sharing among the planetary gears. Considering many factors, manufacturing error, bearing clearance and other factors, Han et al. [7] to [9] established a 2K-V nonlinear dynamic model that was based on D'Alembert's principle and analysed the dynamic transmission error. Besides Zhou and Jia [10] established a 2K-V rigid-flexible coupling dynamic model. Then he carried out theoretical derivation and simulation analysis on transmission accuracy in which the flexible deformation of the planetary wheel with small tooth difference was involved. As for uneven *Corr. Author's Address: Chongqing University, The State Key Laboratory of Mechanical Transmission, Chongqing, 400044, China, gjwang@cqu.edu.cn 91 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 91-104 load sharing, Iglesias et al. [11] pointed out that the unavoidable manufacturing and installation errors would result in the uneven load on the planetary wheel. Then he reduced overall planetary wheel system transmission accuracy. However, Iglesias did not point out the calculation formula of transmission error with the eccentricity error and the uneven load. Mo et al. [12] used the strain measure method to calculate the load-sharing coefficient of the wind turbine gearbox verifying the proposed load sharing calculation model. Wang et al. [13] proposed a novel axial modification method based on composite modification curve with indefinite parameters to improve the contact condition of gear pairs and enhance the meshing performance and bearing capacity. Due to the characteristic that the planetary wheel arrangement in a 2K-V reducer at a symmetrical 180-degree angle, in the transmission process the 2K-V two small tooth difference gear pairs influenced by the eccentricity error may be in different types of contact: with uneven load, snapping off or reversely in contact, etc. Existing literature mainly involves the load-sharing coefficient [14] to [16] itself. However, there are few studies of transmission precision that involve uniform load contact. In addition to 2K-V, there is other research [17] about the nonlinear dynamic characteristics of curve-face gear drive that considered meshing frequency and eccentricity. In this paper, in which the eccentricity error caused by the machining assembly error of 2K-V small tooth difference gear pairs is considered, according to the non-uniform load contact condition, non-load transmission error, time-varying backlash and static transmission error are calculated by using meshing line increment method. Meanwhile. the nonlinear dynamic model of the small tooth difference gear pairs is established by the lumped parameter method, in order to analyse the dynamic transmission accuracy in which the non-uniform load contact is considered, and the correctness of the theoretical model is verified by simulation. Then, in order to improve the precision and dynamic performance of small tooth difference gear drive, the initial eccentricity phase is optimized according to objective function: the minimum transmission error that provides theoretical guidance for the assembly of 2K-V small tooth difference gear pairs. Finally, the 2K-V reducer experimental platform, in order to gain actual transmission error curve that verifies the non-uniform load contact phenomenon of small tooth difference gear pairs in the transmission process, is established. 1 METHODS The meshing line increment method and the lumped parameter method are used to calculate theory static transmission error and dynamic transmission error, respectively. Also, the uneven load contact model is established by the lumped parameter method and solved by the fourth-order five-stage Runge-Kuta numerical integral. Then, the Recurdyn software is used to establish a virtual prototype and verify the theoretical model, as well as optimizing eccentricity initial phase. Finally, the 2K-V transmission error test bench is built to obtain the actual 2K-V transmission error. Combining the above data, the regulation of transmission error is analysed. 2 EXPERIMENTAL 2.1 Calculation of Static Transmission Error and Backlash of 2K-V Gear Pair with Small Tooth Difference Involved Eccentricity Error 2.1.1 Non-Load Transmission Error (NLTE) and Backlash Calculation with Eccentricity Excitation Referring to Figs. 1 and 2, when the input shaft rotates clockwise, the output shaft rotates clockwise (w0j), and the small tooth difference gear pair make counter clockwise revolutions (wg). For the result from machining and assembly error of parts (small tooth difference gear pair, crank shaft, bearing and bearing hole, etc.), the eccentricity error will cause the rotation eccentricity error and revolution rotation eccentricity error of small tooth difference gear pairs. Fig. 1. 2K-V reducer transmission principle: 1 sun-gear, 2 planet of first stage, 3 crankshaft, 4 planet of second stage, 5 ring-gear, and 6 output wheel 92 Wang, G. - Su, L. - Zou, S. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 91-104 a) b) Fig. 2. a) 2K-V reducer transmission principle, and b) increment of meshing line As shown in Fig. 2, let epi (i = 1, 2) defines the planet's eccentricity vector about the rotation centre, 0pi (i = 1, 2) defines the planet's eccentricity phase about the rotation centre, eg defines the planet's eccentricity vector about the revolution centre, 0g defines the planet's eccentricity phase about the revolution centre, eot defines the output wheel's eccentricity vector about the rotation centre, and 0ot defines the output wheel's eccentricity phase about the rotation centre. Then the increment of meshing line [18] and back meshing line caused by eccentricity error of small tooth difference gear pairs at any time is described by: Ai(t) = ePi sin[0pi(t) -a] + eg sin[0g(t) -a] +eot sin[0ot (t) -a] - ep1sin[^p1(t0) -a] -e sm[(9g (t0) -a] - eot sm[0„t(to) -a], (1) A2(t) = ep2 sin[0p2(t) -a] + eg sin[0g(t) -a + n] +eo, sin^t(t) -a+n]-ep2 sin[0p2(t0)-a] -e ¡smM(to)-a+n]-eot sm^,(to)-a+Kl (2) am (t) = epi sin[0pi (t) + a ] + eg sin[0g (t) + a ] +eo> sin[0„t(t) + a] - episin[6pi(to) + a] -eg sin[0g(t0) + a] -eot sin^(t0) + a], (3) An(t) = ep2 sin[0p2(t) + a] + eg sin^(t) + a + n ] +eot Sin[0<,( (t) + a+n] - ep2 Sin[^p2 (t0) + a] -eg sin[0g(to) + a+n] - eot sin^, (to) + a+n], (4) where t and to represent any time and initial time respectively; superscript b refers to the tooth back. Then, according to the meshing line increment and tooth back meshing line increment, the non-load transmission error and tooth back transmission error (subscript i = 1, 2) at any time are derived as: A,, (t) 180 • 60 NLTEi (t) = - NLTEbi(t) = Abi (t) 180 • 60 (5) (6) where rp refers to the base circle radius of ring-gear; the unit of transmission error is arcmin. According to the calculated transmission errors and tooth back transmission error of the gear pairs without load, the varying backlash at any time is given by: Bt (t) = NLTEbi(t) - NLTEi(t) + B(t0). (7) Referring to Fig. 3, due to the eccentricity error excitation, the 2K-V small tooth difference gear pairs will occur the non-uniform load contact phenomenon, but in the case of no-load, only one pair of gear pairs will contact in theory. Hence, the ultimate non-load transmission error of 2K-V is the non-load transmission error of the gear pair, which is actually in contact, i.e., the ultimate transmission error is the maximum of the non-load transmission error of the two small tooth difference gear pairs, and the ultimate time-varying backlash is the minimum of actual time-varying backlash. Therefore, the non-load transmission error and time-varying backlash of the 2K-V small tooth difference gear pairs are described as: NLTE(t) = max( NLTEl (t), NLTE2 (t)), (8) Uneven Load Contact Dynamic Modelling and Transmission Error Analysis of a 2K-V Reducer with Eccentricity Excitation 93 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 91-104 B(t) = min[B1 (t ), B2 (t ) - NLTEl (t ) + NLTE2 (t )] if NLTE1(t) > NLTE2(t) min[ B1 (t ) - NLTE2 (t ) + NLTE1 (t ), B2 (t )] ' if NLTE1(t) < NLTE2(t) (9) A, (t ) = A,, (t) = Al(t) + A 2(t ) T 2• rp • k' STE (t ) = ( A.(t ) + A 2(t ) T 2 • rp • fc' 180•60 rp n (13) (14) 2.1.2 Static Transmission Error Calculation with Eccentricity Excitation Referring to the definition in [19], not only the above motion error, which is caused by the eccentricity error, but also the tooth deformation, which is caused by load torque, is attributed to the static transmission error of 2K-V small tooth difference gear pairs under eccentricity excitation. Assuming that the two small tooth difference gear pairs are all in contact with load torque, based on the static stress balance condition, the contact forces of the gear pair 1 and gear pair 2 are given by: F + F contactl contact 2 (10) and Fcon,ac,1 - Fcontact2 = k(A1 (t) - A2 (t)) F = (A,(t ) - A2(t )) contactl 2r 2 F = (Ai(t) -A2(t ))' contact 2 2r 2 (10) Eq. (10) shows that if \A1(t)-A2(t)| >, there rp ■k is only one gear pair that is in contact and the two gear pairs contact force can be derived as: F„ = Fcontact 2 = 0,if A,(t )-A 2(t ) > 2r F = 0F contact1 ' contact 2 Ti ^ • k = if A 2 (t )-Ai(t ) 2rp rp ■k (15) Therefore, if \A1(t)-A2(t)| , the actual rp ■ k meshing line increment and static transmission error of the 2K-V small tooth difference gear pairs with load torque are calculated as: A, (t) = max(A,(t), A2(t)) --(16) STE (t ) = [max( A1 (t ), A 2 (t )) - rp ■k T , 180•60 rp •k rpn (17) In summary, the static transmission error of the 2K-V small tooth difference gear pair with load torque is found as the following equation: where Fcontact1 stands for contact force of gear pair 1; FCOntac& represents contact force of gear pair 2; k refers to the gear meshing stiffness. When the two gear pairs are all in contact at the same time, the contact force of the gear pairs 1 and 2 must be greater than 0 and the inequality, i.e., T if |Aj(Y) — A2(t)| <—— , the actual meshing line rp • k increments of gear pair 1 and 2 with load torque are derived as: A„ (t ) = Aj(t ) - F Ai(t ) + A2(t ) T 2 • rp • k A 2i (t ) = A 2(t ) - F„ At(t ) + A 2 (t ) T 2 • rp • k (11) (12) Hence, when the two small tooth difference gears are all in contact, the actual meshing line increment and static transmission error with load torque are described as: STE (t) = [max(Aj (t ), A 2(t )) - T1 180 • 60 rp •k if\Aj(t) -A2(t)| rp •k ( Aj +A2 T1 180•60 (18) 2 • rp • Ik if\ Aj(t ) -A 2 (t )|- T1 rp •k When the load torque is 0, it can be seen from Eq. (10) that the static transmission error of the 2K-V small tooth difference gear pairs with the eccentricity excitation is equal to the non-load transmission error, that is, the non-load transmission error of 2K-V small tooth difference gear pair is the static transmission error. Referring to Eqs. (10) and (15), due to the static force balance and the existence of eccentricity error, the contact forces of gear pair 1 and gear pair 2 are uneven, and even one of the gear pairs may occur mesh apart during the transmission process. 94 Wang, G. - Su, L. - Zou, S. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 91-104 2.2 Dynamic Modelling and Dynamic Transmission Error Calculation with Eccentricity Excitation of 2K-V Small Tooth Difference Gear Pairs 2.2.1 Dynamic Modelling The vibration displacement on the meshing line direction of the small tooth difference gear pair is defined as: x = rp°i - r^v ' r ■ n where A(t) = e, Ai(t) = ei, ^(t) = e2, B(t) —--= b . 180■60 In the non-load static transmission process, the vibration displacement (x) of the gear pairs is the meshing line increment corresponding to the non-load transmission. As long as the driving surface of the small tooth difference gear pair is set as a reference, referring to Fig. 3, due to velocity, load, time-varying backlash, damping and time-varying stiffness, when the vibration displacement (x) changes up or down in the meshing line increment corresponding to the transmission error without load, the following situation may occur: Case 1. If x- e> | e1 - e2 | , the driving surfaces of the two small tooth difference gear pairs are in contact; Case 2. If 0 < x- e< | e1 - e2 | and b> 0, one of the two small tooth difference gear pairs contacts with the driving surface, and another gear pair doesn't contact; Case 3. If 0 < x- e< | e1 - e2 | and b< 0, one of the two small tooth difference gear pairs contacts with the driving surface, and another gear pair back contact; Case 4. If -b < x- e < 0, none of the small tooth difference gear pairs contact; Case 5, If x- e< 0 and -| e1 - e2 | < x- e + b< 0, one of the two small tooth difference gear pairs do not contact, and another gear pair back contact; Case 6. If x- e + b< -| e1 - e2 |, all of the small tooth difference gear pairs back contact. According to the lumped parameter method, the dynamic equation of the gear pair is established as follows: m (x - e) - C (x - e) - f (x - e) = -, (19) r p where me stands for equivalent mass; Cg refers to meshing damping of gear pair; f(x- e) represents nonlinear function which is derived as follow: case 1 : k(x - e) + k(x - e -\e1 - e21) case 2: k ( x - e) case 3 : k(x - e - b / 2) case 4: 0 case 5: k ( x - e + b) case 6 : k(x - e + b) + k(x - e + b + \e1 - e21) f ( x - e) = where k is the time-varying stiffness. Expanding k by Fourier series at mesh frequency, it is described by: k(t) = k0 + cos(imt + fa)], (20) Fig. 3. Contact of the small tooth difference gear pairs without load; ej > e2, the gear pair 1 contacts, and b) ej < e2, the gear pair 2 contacts where k0 refers to the average meshing stiffness; a stands for the excitation frequency of gear pair; kj is /h order stiffness amplitude of the meshing gear pair; ^ represents /h order initial phase of the meshing gear pair. 2.2.2 Dynamic Modelling Example Taking a 2K-V small tooth difference reducer as an example to calculate its transmission error, its main parameters are shown in Table 1 and Table 2, in which the eccentricity error of the second-stage small tooth difference gear pair is set responding to the allowable error in the six-level machining accuracy, the eccentricity phase is set arbitrarily, and the delicate transmission error caused by the first-stage gear pair is ignored. As shown in Fig. 4, referring to Eq. (5), the non-load transmission errors with eccentricity excitation of 2K-V small tooth difference gear pairs can be obtained. Define the difference between the maximum and minimum of the transmission error as the peak-to-peak value, and the maximum of the absolute transmission error as the amplitude. As shown in Fig. 4a, the transmission error without load has a peak-to-peak value of 2.06 arcmin, an amplitude of 1.89 arcmin and an average value of -0.86 arcmin. Correspondingly, referring to Fig. 4b, the transmission error without load has a peak-to-peak i=i Uneven Load Contact Dynamic Modelling and Transmission Error Analysis of a 2K-V Reducer with Eccentricity Excitation 95 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 91-104 Table 1. Main parameters of the first-stage in 2K-V reducer Table 2. Main parameters of the second-stage in 2K-V reducer Parameter Unit Num. value Tooth number of sun-gear - 12 Modification coefficient of sun-gear - 0.5 Tooth number of planet - 36 Modification coefficient of planet - -0.5 Addendum coefficient - 1 Centre distance mm 36 Face width mm 7 Contact ratio - 1.433 Normal module mm 1.5 Pressure angle degrees 20 Clearance coefficient - 0.25 Input velocity rpm 300 value of 2.89 arcmin, an amplitude of 1.74 arcmin and an average value of -0.29 arcmin. Furthermore, the non-load transmission error curves in Fig. 4 all have a major cycle of 15.8 s. In order to facilitate the concentrated analysis and observation, the first 3 s of Parameter Unit Num. value Tooth number of planet - 50 Modification coefficient of planet - -0.1060 Tooth number of ring-gear - 52 Modification coefficient of ring-gear - 0.1049 Normal module mm 2.5 Addendum coefficient - 0.6 Clearance coefficient - 0.25 Pressure angle degrees 35.266 Face width mm 10 Centre distance mm 2.8773 Contact ratio - 1.1 Eccentricity error of planet 1 mm 0.01 Eccentricity phase of planet 1 degrees 90 Eccentricity error of planet 2 mm 0.02 Eccentricity phase of planet 2 degrees 180 Eccentricity error of revolution mm 0.01 Eccentricity phase of revolution degrees -135 Eccentricity error of output wheel mm 0.01 Eccentricity phase of output wheel degrees 90 Fig. 4. Transmission error curve of small tooth difference gear pairs without load; a) gear pair 1, and b) gear pair 2 Fig Time fs] 5. Transmission error of small tooth difference gear pair; a) transmission error without load, and b) static transmission error 96 Wang, G. - Su, L. - Zou, S. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 91-104 0 5 ■2.5 DTE mean value w —t— STE mean value ■ - X ■ - I . 1 . t * I 1 250 300 Fig. 6. lJ 2.0 2,5 3,0 0 50 100 150 200 Time [s] Load torque [Nm) Transmission error of small tooth difference gear pair; a) dynamic transmission error, and b) comparison between DTE and STE _ DTE with ÍONptt sinusoidal 1j Thermal treatment temperature [°C] Fig. 4. Roughness surface parameters for heat-treated samples; a) Rz, and b) Ra and Rq a) h] Thermal treatment temperature Fig. 5. Waviness surface parameters for heat-treated samples; a) Wz, and b) Wa and Wq This is different for the waviness parameters depicted in Fig. 5. Here a clear difference between the untreated samples and the ones thermally treated at temperatures of 80 °C or higher is visible, with all three waviness values Wz, Wa, and Wq decreasing significantly to approximately half of the original value. Comparing Figs. 4 and 5 shows that temperature treatment indeed has an influence on the surface of FDM printed objects from PLA, although only the macro-scale (the waviness) is modified, while the micro-scale (the roughness) stays unaltered. Next, the influence of an acetone treatment was investigated. Fig. 6 depicts the roughness values Rz, Ra and Rq of untreated samples in comparison with acetone treatment in gaseous form. While no significant changes are visible for the shorter exposition durations, an increase of the roughness for the longest treatment is visible, which is nevertheless not significant. However, these measurement results fit well to the observation that the surface of these samples seems to be slightly less shiny than before the treatment. a) Acetone gas T □ Ra - M 2.5 60 l)j Treatment duration [min] Fig. 6. Roughness surface parameters for acetone-treated samples (gas); a) Rz, and b) Ra and Rq 108 Kozior, T. - Mamun, A. - Trabelsi, M. - Sabantina, L. - Ehrmann, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 105-113 a) „ 1.0 0.5 DH 0.0 Acetone fluid | Ra CPRq 1 2.5 60 y Treatment duration [min] Fig. 7. Roughness surface parameters for acetone-treated samples (liquid); a) Rz, and b) Ra and Rq 250 200 ^ 150 f 100 CD 50 j Acetone gas - . [□□Wz_\. 1 2.5 60 a) Treatment duration [min] 12 10 « 6 > 4 m ~ 0 Acetone gas □ Wa I Mi 0 1 2.5 60 y Treatment duration [min] Fig. 8. Waviness surface parameters for acetone treated samples (gas); a) Wz, an b) Wa and Wq The impact of treatment with fluid acetone on the surface roughness is shown in Fig. 7. Interestingly, here even the shortest time shows a significant change of Rz, while the other values remain similar. Is should be mentioned that after a treatment duration of 60 min, the original value of the untreated specimen is approached again. While the treatment with acetone in fluid or gaseous form was disadvantageous - or in the best case neutral - in terms of surface roughness, Figs. 8 and 9 show that especially Wz is significantly reduced by an acetone treatment, while the other two waviness values are slightly reduced. This finding underlines that on larger scales, the surface gets more even, while the shine is reduced if the surface is not polished. Apparently, the planned application decides whether an acetone treatment makes sense as an after-treatment, or whether it is counterproductive. 250 200 150 100 CD 50 t Acetone fluid] - nniVi 1 2.5 60 a) Treatment duration [min] 0 1 2.5 60 y Treatment duration [min] Fig. 9. Waviness surface parameters for acetone treated samples (liquid); a) Wz, an b) Wa and Wq Generally, quantitatively analysing the results of metrological measurements, it can be stated that the thermal and chemical treatment with acetone has a slight effect on surface roughness parameters; however, the influence of treatments on the surface waviness of samples is clearly noticeable. The surface Quality of the Surface Texture and Mechanical Properties of FDM Printed Samples after Thermal and Chemical Treatment 109 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 105-113 texture parameters measured after these treatments are concluded in Table 3. Table 3. Surface texture parameters in [jm] No. Ra Rz Rq Wa Wz Wq 1 0.42 8.8 0.59 6.9 240 10.5 2 0.27 8.0 0.38 6.4 215 9.5 3 0.28 6.6 0.37 5.5 199 8.3 4 0.29 8.0 0.40 3.9 120 5.5 5 0.20 6.3 0.29 3.4 103 5.0 6 0.40 7.9 0.56 5.4 194 8.0 7 0.24 7.3 0.32 3.6 132 5.1 8 0.28 8.9 0.40 3.9 126 5.4 9 0.16 5.9 0.23 3.2 107 4.1 10 0.20 7.5 0.27 3.4 114 4.5 11 0.20 7.8 0.28 3.8 127 5.3 12 0.31 8.4 0.46 3.4 131 4.4 13 0.31 7.0 0.45 3.2 128 4.3 14 0.17 7.9 0.25 4.2 128 5.6 15 0.27 8.1 0.38 3.7 123 4.6 16 0.27 7.0 0.38 3.2 133 4.1 17 0.27 7.4 0.38 3.1 134 4.2 18 0.19 6.1 0.26 3.6 102 4.6 9 IQ 0.14 8.5 0.24 3.4 74 4.3 20g 0.33 8.0 0.46 5.9 216 8.3 21g 0.63 11.6 0.93 4.2 91 5.1 19f 0.41 17.9 0.58 4.9 118 6.2 20f 0.54 13.5 0.75 4.2 96 5.6 21f 0.37 9.0 0.48 5.4 75 6.7 Table 4. Mechanical properties of samples: tensile strength Rm and elongation A No. Rm [MPa] A [%] No. Rm [MPa] A [%] 1 35.532 10.594 12 45.734 6.952 2 45.962 8.119 13 44.958 6.813 3 45.177 6.916 14 46.411 6.742 4 46.347 8.647 15 47.262 7.130 5 44.381 8.184 16 39.265 5.704 6 41.884 8.886 17 45.138 5.453 7 39.381 8.298 18 48.035 6.599 8 41.652 7.331 19 g/f 31.579 6.865 9 47.019 8.004 20 g/f 25.68 4.937 10 45.648 8.107 21 g/f 8.369 31.769 11 43.147 6.836 Analysing the results of tensile strength measurements of all types of samples, as depicted in Table 4, it can be seen that in almost all cases there is no clear effect of heat and chemical treatment on the tensile strength of the samples and on elongation. The exception is the tensile strength of the sample that was acetone treated for a period of 60 minutes. In most cases, the tensile strength ranges from 35 MPa to 48 MPa, and the elongation is approximately 3 mm (5 % to 10 %). Tensile strengths of approximately 40 MPa to 50 MPa are typically for PLA 3D printed with a nozzle diameter of 0.4 mm and rectilinear infill. The maximum values of approximately 50 MPa are given by the tensile strength of the pure filament, which is in the same range, while small deviations from a perfect printer setup lead to air holes in the specimens and thus to reduced tensile strength [23]. For sample 21 (acetone 60 min), however, the tensile strength is much lower (8.369 MPa), however, the elongation is almost 12 mm (31.769 %). This is interesting because acetone treatment has positively influenced the waviness parameters. The increased elongation of sample 21 means that acetone treatment allows changing the physicochemical properties of samples and increasing the flexibility of the produced models, which in justified cases is a beneficial phenomenon, despite the lower tensile strength. This is also visible in Fig. 10, showing the clear difference between sample 21 (lilac curve) and all other curves. Obviously, this sample is severely elongated, while the other samples, including Sample 1, which breaks at an elongation of 4 mm, show a brittle break. Fig. 10. Stress-strain curves of all samples under investigation in this study To investigate the reason for this difference, FTIR measurements on the samples with and without after-treatment were performed, which are depicted in Fig. 11. In all cases, we see the typical peaks of the PLA backbone ester group, i.e., 1748 cm-1 (-C=O stretching), 1182 cm-1 (-C-O-C stretching), 1133 cm-1 and 1079 cm-1 (both -C-O-stretching) (all marked by solid lines) [24]. Additional typical PLA peaks can be found at 1452 cm-1 (-CH3 bending), 1268 cm-1 (-C=O bending), 1038 cm-1 (-OH bending) as well as 872 cm-1 and 754 cm-1 (-C-C- stretching of amorphous and crystalline 110 Kozior, T. - Mamun, A. - Trabelsi, M. - Sabantina, L. - Ehrmann, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 105-113 Fig. 11. FTIR measurements of untreated, thermally and acetone-treated specimens; the graphs are vertically offset for clarity phases, respectively) (all marked by dashed lines) [25]. Finally, the small peaks around 1358 cm-1 and 1384 cm-1 can be attributed to -CH vibration and symmetric -CH3 bending, respectively (dotted lines) ([26] and [27]). Differences between the untreated material and the material treated with liquid acetone for 60 min are not visible in the peaks at 872 cm-1 and 754 cm-1, indicating no change of the crystallinity near the surface where acetone-induced changes could be expected to be strongest [28]. Next, a possible modification of the molecular weight is investigated. An increase of the molecular weight by solvent treatment is indeed possible since low-molecular fractions may be dissolved, while higher-molecular fractions stay unaltered [29]. If the molecular weight is altered, changes should mostly be visible in the amount of end functional groups and ester bonds since low-molecular-weight PLA has more carboxylic acid and alcohol groups per mass than higher-molecular-weight PLA. On the other hand, higher-molecular weight PLA has more ester bonds due to enlarging the polymeric chain. This means that the peaks around 1748 cm-1 (carboxylic acid) should be higher for lower-molecular weight PLA, which is not the case for the acetone-treated samples. In [27] the high-molecular-weight PLA also shows an additional small peak around 1715 cm-1 on the shoulder of the large carboxylic acid peak, which is well visible here for most acetone treated specimens, especially for a duration of 60 min (marked by an arrow). These findings might indicate an increase in the molecular weight due to an acetone treatment. On the other hand, the peak at 1268 cm-1 is significantly higher in our measurement of the original PLA, while it nearly vanishes for the acetone-treated specimens. Vanishing of this carbonyl ester group shows a reduced amount of these ester links, i.e. a lower molecular weight. Similarly, the peak at 1210 cm-1 - indicating the -OC-O stretching of carboxylic groups - is much more pronounced for the acetone-treated samples than for the original one, which is again a sign of a lower molecular weight. Since these findings are inconsistent, a possible modification of the molecular weight due to the acetone treatment cannot be verified or falsified by the FTIR measurements. Generally, it can be assumed that treating the specimens with a solvent partly destroys the molecular structure of PLA and thus weakens the molecular chain interaction [30], changing the original brittle behaviour into a more ductile one. 4 CONCLUSIONS After an analysis of the presented research results, the following general conclusions can be formulated: The geometric structure of the surface after thermal and chemical treatment is different from that immediately after printing. It has been shown that this effect is insignificant with respect to roughness parameters; however, it is clearly visible in the case of surface waviness parameters, which in many cases are incorrectly ignored. Tensile strength when using thermal treatment did not change significantly; however, in the case of chemical treatment with acetone, a noticeable decrease in strength intensified by the elongation of the acetone effect can be seen. At the same time, a large increase 111 Quality of the Surface Texture and Mechanical Properties of FDM Printed Samples after Thermal and Chemical Treatment 39 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 105-113 in elongation was observed, especially for samples subjected to 60 minutes of acetone treatment, which is a positive phenomenon in many applications. Tests also showed that the interaction with liquid acetone for a duration longer than one hour causes the samples to break down with significant deformation. The FTIR measurements are inconsistent; a possible modification of the average molecular weight due to the acetone treatment cannot be verified or falsified and should be investigated in a follow-up study using much longer solvent treatment durations. 5 ACKNOWLEDGEMENTS This research was funded by the DAAD, funding program number 57440917 - Research Grant. The APC is funded by the Open Access Publication Fund of Bielefeld University of Applied Sciences and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 414001623. 6 REFERENCES [1] Chua, C.K., Leong, K.F., Lim, C.S. (2003). Rapid Prototyping: Principles and Applications, 2nd Edition. World Scientific Publishing, Singapore. [2] Novakova-Marcincinova, L. (2012). Application of fused deposition modeling technology in 3D printing rapid prototyping area. Manufacturing and Industrial Engineering, vol. 11, no. 4, p. 35-37. [3] Noorani, R. (2005) Rapid Prototyping: Principles and Applications. John Wiley & Sons, New Jersey. [4] Ben-Ner, A., Siemsen, E. (2017). Decentralization and localization of production: the organizational and economic consequences of additive manufacturing (3D printing). California Management Review, vol. 59, no. 2, p. 5-23, DOI:10.1177/0008125617695284. [5] Kousiatza, C., Tzatzis, D., Karalekas, D. (2019). In-situ characterization of 3D printed continuous fiber reinforced composites: A methodological study using fiber Bragg grating sensors. Composites Science and Technology, vol. 174, p. 134-141, DOI:10.1016/j.compscitech.2019.02.008. [6] Li, Q.S., Zhao, W., Li, Y.W., Yang, W.W., Wang, G. (2019). Flexural Properties and Fracture Behavior of CF/PEEK in Orthogonal Building Orientation by FDM: Microstructure and Mechanism. Polymers, vol. 11, no. 4, p. 656, DOI:10.3390/ polym11040656. [7] Kozior, T., Trabelsi, M., Mamun, A., Sabantina, L., Ehrmann, A. (2019). Stabilization of electrospun nanofiber mats used for filters by 3D printing. Polymers, vol. 11, no. 10, p. 1618, DOI:10.3390/polym11101618. [8] Adamczak, S., Zmarzty, P., Kozior, T., Gogolewski, D. (2017). Analysis of the dimensional accuracy of casting models manufactured by fused deposition modeling technology. Engineering Mechanics 2017, 23rd International Conference, p. 66-69, DOI:10.2507/29th.daaam.proceedings.123. [9] Rokicki, P., Budzik, G., Kubiak, K., Bernaczek, J., Dzlubek, T., Magniszewski, M., Nowotnik, A., Sieniawski, J., Matysiak, H., Cygan, R., Trojan, A. (2014). Rapid prototyping in manufacturing of core models of aircraft engine blades. Aircraft Engineering and Aerospace Technology, vol. 86, no. 4, p. 323-327, DOI:10.1108/AEAT-10-2012-0192. [10] ISO 10993-1:2010. Biological Evaluation of Medical Devices - Part 1: Evaluation and Testing within a Risk Management Process. International Standard Organization, Geneva. [11] Kozior, T., Dopke, C., Grimmelsmann, N., Juhász Junger, I., Ehrmann, A. (2018). Influence of fabric pretreatment on adhesion of three-dimensional printed material on textile substrates. Advances in Mechanical Engineering, vol. 10, no. 8, DOI:10.1177/1687814018792316. [12] Mamun, A., Trabelsi, M., Klocker M., Sabantina, L., GroSerhode, C., Blachowicz, T., Grotsch, G., CorneliSen, C., Streitenberger, A., Ehrmann, A. (2019). Electrospun nanofiber mats with embedded non-sintered TiO2 for dye-sensitized solar cells (DSSCs). Fibers, vol. 7, no. 7, p. 60, DOI:10.3390/ fib7070060. [13] Kundera, Cz., Kozior, T. (2014). Research of the elastic properties of bellows made in SLS technology. Advanced Materials Research, vol. 874, p. 77-81, DOI:10.4028/www. scientific.net/AMR.874.77. [14] Kundera, Cz., Martsynkowskyy, V., Gudkov, S., Kozior, T. (2017). Effect of rheological parameters of elastomeric ring materials on dynamic of face seals. Procedia Engineering, vol. 177, p. 307-313, DOI:10.1016/j.proeng.2017.02.230. [15] Krolczyk, G., Raos, P., Legutko, S. (2014). Experimental analysis of surface roughness and surface texture of machined and fused deposition modeling parts. Tehnički Vjesnik - Technical Gazette, vol. 21, no. 1, p. 217-221. [16] Kozior, T., Kundera, Cz., (2017). Evaluation of the influence of parameters of FDM technology on the selected mechanical properties of models. Procedia Engineering, vol. 192, p. 463468, DOI:10.1016/j.proeng.2017.06.080. [17] Bochnia, J., Blasiak, S. (2018). Fractional relaxation model of materials obtained with selective laser sintering technology. Rapid Prototyping Journal, vol. 25, no. 1, p. 76-86, DOI:10.1108/RPJ-11-2017-0236. [18] Garg, A., Bhattacharya, A., Batish, A. (2016). On surface finish and dimensional accuracy of FDM parts after cold vapor treatment. Materials and Manufacturing Processes, vol. 31, no. 4, p. 522-529, DOI:10.1080/10426914.2015.1070425. [19] Rao, A.S., Dharap, M.A., Venkatesh, J.V.L., Ojha, D. (2012). Investigation of post processing techniques to reduce the surface roughness of fused deposition modeled parts. International Journal of Mechanical Engineering and Technology, vol. 3, no. 3, p. 531-544. [20] Beniak, J., Križan, P., Šooš, L., Matuš, M. (2018). Roughness and compressive strength of FDM 3D printed specimens affected by acetone vapour treatment. IOP Conference Series: Materials Science and Engineering, vol. 297, DOI:10.1088/1757-899X/297/1/012018. [21] Havenga, S.P., de Beer, D.J., van Tonder, P.J., Campbell, R.I. (2018). The effect of acetone as a post-production finishing technique on entry-level material extrusion part quality. The 112 Kozior, T. - Mamun, A. - Trabelsi, M. - Sabantina, L. - Ehrmann, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 105-113 South African Journal of Industrial Engineering, vol. 29, no. 4, DOI:10.7166/29-4-1934. [22] Wach, R.A., Wolszczak, P., Adamus-Wtodarczyk, A. (2018). Enhancement of mechanical properties of FDM-PLA parts via thermal annealing. Macromolecular Materials and Engineering, vol. 303, no. 9, D0l:10.1002/mame.201800169. [23] Fafenrot, S., Grimmelsmann, N., Wortmann, M., Ehrmann, A. (2017) Three-dimensional (3D) printing of polymer-metal hybrid materials by fused deposition modeling. Materials, vol. 10, no. 10, p. 1199, DOI:10.3390/ma10101199. [24] Suarez-Franco, J.L., Vázquez-Vázquez, F.C., Pozoz-Guillen, A., Montesinos, J.J., Alvarez-Fregoso, O., Alvarez-Perez, M.A. (2018). Influence of diameter of fiber membrane scaffolds on the biocompatibility of hPDL mesenchymal stromal cells. Dental Materials Journal, vol. 37, no. 3, p. 465-473, D0I:10.4012/dmj.2016-329. [25] Rocca-Smith, J.R., Lagorce-Tachon, A., laconelli, C., Bellat, J.P., Marcuzzo, E., Sensidoni, A., Piasente, F., Debeaufort, F., Karbowiak, T. (2017). How high pressure CO2 impacts PLA film properties. Express Polymer Letters, vol. 11, no. 4, p. 320333, D0I:10.3144/expresspolymlett.2017.31. [26] Silva Murakami, L.M., Baracho Azevedo, J., Diniz, M.F., Silva, L.M., de Cássia Lazzarini Dutra, R. (2018). Characterization of additives in NR formulations by TLC-IR (UATR). Polímeros, vol. 28, no. 3, p. 205-214, D0I:10.1590/0104-1428.06317. [27] Palacio, J., Orozco, V.H., Lopez, B.L. (2011). Effect of the molecular weight on the physicochemical properties of poy(lactic acid) nanoparticles and on the amount of ovalbumin adsorption. Journal of the Brazilian Chemical Society, vol. 22, no. 12, p. 2304-2311, D0I:10.1590/S0103-50532011001200010. [28] Goncalves, C.M.B., Coutinho, J.A.P., Marrucho, I.M. (2010). Optical Properties. Poly(lactic acid): Synthesis, Structures, Properties, Processing, and Applications. Auras, R.; Lim, L.-T.; Selke, S.E.M.; Tsuji, H. (eds.), John Wiley & Sons, Hoboken, p. 97-112, D0I:10.1002/9780470649848.ch8. [29] Renbenberger, R., Cayla, A., Villmow, T., Jehnichen, D., Campagne, C., Rochery, M., Devaux, E., Potschke, P. (2011). Multifilament fibres of poly(s-caprolactone)/poly(lactic acid) blends with multiwalled carbon nanotubes as sensor materials for ethyl acetate and acetone. Sensors and Actuators B: Chemical, vol. 160, no. 1, p. 22-31, D0I:10.1016/j. snb.2011.07.004. [30] Jin, Y.F., Wan, Y., Zhang, B., Liu, Z.Q. (2017). Modeling of the chemical finishing process for polylactic acid parts in fused deposition modeling and investigation of its tensile properties. Journal of Materials Processing Technology, vol. 240, p. 233239, D0I:10.1016/j.jmatprotec.2016.10.003. Quality of the Surface Texture and Mechanical Properties of FDM Printed Samples after Thermal and Chemical Treatment 113 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2019.6267 Original Scientific Paper Received for review: 2019-07-26 Received revised form: 2019-11-06 Accepted for publication: 2019-11-15 The Loading Characteristics of Bulk Coal in the Middle Trough and Its Influence on Rigid Body Parts Peilin Zhang1,2 - Bo Li1,2 - Xuewen Wang1,2,* - Chaoyang Liu1,2 - Wenjie Bi1,2 - Haozhou Ma1,2 1 Shanxi Key Laboratory of Fully Mechanized Coal Mining Equipment, China 2 Taiyuan University of Technology, College of Mechanical and Vehicle Engineering, China The working reliability of a scraper conveyor is related to the bulk coal load characteristics and the interaction between rigid body parts and bulk coal. In the existing relevant research, a load of coal on rigid body parts is equivalent to constant central acting on some points, ignoring time-varying and uneven distribution of the load and new load generated by the interaction, which results in some load characteristics and related faults being obscured. In this paper, the rigid-discrete coupling simulation model of scraper conveyor is established based on the discrete element method and the multi-body dynamic method, and the overall compressive force gradient distribution of bulk coal in the middle trough, the violent interaction between bulk coal and rigid body parts, and the clamping stagnation of scrapers are studied. The results show that the bulk coal compressive force has an obvious distribution gradient in the chute, and in the chute between two scrapers, the bulk compressive force in the first half is about 19 % of that in the second half in the y-direction. The coal particles stuck between the scrapers, the vertical chainrings, and the middle plates violently interact with the rigid body parts, and the coal particles average compressive force is 59 times that of other coal particles. The scraper influenced by bulk coal is easily stuck at the connection of two chutes, and the tension of the chain before the scraper increases sharply, and the maximum tension of the chainrings is about three times that of the maximum contact force between the scraper and the chute. Keywords: scraper conveyor, discrete element method, multi-body dynamics, compressive force distribution gradient, violent interaction Highlights • The rigid-discrete coupling model of a scraper conveyor can simulate the interaction between the bulk coal and the rigid parts. • The compressive force distribution of coal in the middle trough has an obvious gradient, and the load characteristics of the rigid body parts applied by bulk coal can be analysed through the bulk coal compressive force distribution. • Coal particles are easily stuck between scraper, vertical chainring and middle plate, which affects the force and movement of the scraper and the chain as well as the wear of middle plate. • It is easy for the scraper to get stuck at the connection of two chutes, which causes the tension of the chain to increase sharply. 0 INTRODUCTION A scraper conveyor in a fully mechanized mining face is the coal transport equipment, together with shearer and hydraulic support, that comprises a fully mechanized mining face [1] and [2]. The reliability and service life of scraper conveyor have a great influence on the coal mining efficiency and the economic efficiency of coal mining. The working environment of a scraper conveyor is very harsh, which makes it prone to failure. Once a fault occurs, the normal excavation work will be interrupted, and the underground maintenance is very difficult, which wastes a great deal of manpower, material, and financial resources. Therefore, the study of scraper conveyor failure mechanism, efforts to improve its work reliability, and extending the working life are the long-standing pursuits of many researchers. There are primarily three failure modes for the relevant parts of scraper conveyor: the chain break, the wear of middle plate, and of sprocket drum. Researchers have mainly studied the motion and the force characteristics of related parts, the wear 114 *Corr. Author's Address: Taiyuan University < mechanism, and wear resistance measures of the middle plate and sprocket drum. Based on the finite element concept, the physical and mathematical models of scraper conveyor are established to study dynamic characteristics of the chain drive system [3]. The polygon effect of the chain drive has an important impact on the dynamic performance of the scraper conveyor [4] and [5]. The load on the scraper and chain is uneven and unpredictable, and the load under different working conditions and the load change at any point affect the dynamic characteristics of the whole drive system [6] to [9]. The contact force change characteristics between the chain ring, scraper and sprocket, flat ring and vertical ring, as well as the stress and strain of these parts, are studied [10] and [11]. The chain tension also affects the motion and force characteristics of the chain drive system, and real-time monitoring and adjusting the tension is also one significant research direction [12] and [13]. The driving force of the scraper conveyor can affect the dynamic load of the chain drive. Optimizing the control algorithm of the driving system can improve the dynamic behaviour of the scraper conveyor [14] ology, College of Mechanical and Vehicle Engineering, Taiyuan, China, wxuew@163.com Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 and [15]. The dynamic load characteristics, the wear and the wear resistance improvement of sprocket drum are critical research directions [16] to [18]. Many studies have been done on the medium plate wear mechanism, the factors affecting the wear and the improvement of the wear resistance [19] to [21]. Most of these studies assume the load imposed by coal bulk to the relevant parts as a constant value, ignoring the time-varying, uneven distribution of coal bulk load and the possible new load generated by the interaction of coal bulk and parts, resulting in the concealment of some force characteristics and failure mechanism of the relevant parts. The discrete element method (DEM) is a numerical calculation method proposed by Cundall et al. to simulate bulk materials, which has been widely used in studies related to bulk materials [22] and [23]. Wang et al. studied the size distribution, velocity distribution, and transport efficiency of coal particles in the middle trough with DEM, but the rigid body parts in their model could only move along a straight line [24] and [25]. DEM provides a method to simulate bulk material, but the geometric model in DEM commercial software EDEM can only realize simple motion, so it is necessary to combine it with other dynamic software to study the mechanical system with complex motion. Multi-body dynamics (MBD) is a subject that studies the motion rules of multi-body systems; it can be used to establish models to analyse the force and motion rules for complex multi-body systems [26] and [27]. In the related studies on scraper conveyors, MBD is mostly used to study the motion and force characteristics of the chain drive system [11] and [28]. MBD provides a powerful tool for studying the force characteristics among geometries in a multi-body system. This paper focuses on the overall compressive force distribution gradient of coal bulk in the middle trough, the violent interaction between coal bulk and rigid body parts, and the influence of bulk on rigid body parts, especially on the scraper and chain. The overall compressive force gradient distribution reveals the load distribution characteristics of bulk coal, the violent interaction between bulk coal and rigid parts produces new loads, and the influence of bulk coal on rigid parts reveals the principle of some faults. This research is organized as follows. In Section 1, the coupling model is introduced and its reliability is verified; in Section 2, the simulation results with discussions are presented; in Section 3, the conclusions are given. 1 RESEARCH METHOD 1.1 Scraper Conveyor Model The type of scraper conveyor in this study is SGZ880/800, and its specific parameters are shown in Table 1. The complete virtual prototype of scraper conveyor is shown in Fig. 1a, and the simplified model of scraper conveyor in this study is shown in Fig. 1b. The geometric model is modelled with Siemens Unigraphics NX 10.0 (UG 10.0). The longest scraper conveyor used in longwall system can reach to 500 m, and the number of links in the chain is hundreds of thousands [6]. If a complete scraper conveyor model is adopted to establish the simulation model, the number of rigid bodies and kinematic pairs will be greatly increased, and the number of coal particles participating in the simulation will also be increased dozens of times, which can cause simulation on an ordinary computer to be time-consuming or not feasible at all. Therefore, the simulation Table 1. Specific parameters of the scraper conveyor Type Power [kW] Chain dxt [mmxmm] Chute LxBxH [mmxmmxmm] SGZ880/800 2x400 34x126 1500x880x330 Tait drive Spillplate b) Spillplate Head sprocket drum Tail sprocket drum Fig. 1. Virtual prototype of scraper conveyor; a) the complete; and b) the simplified The Loading Characteristics of Bulk Coal in the Middle Trough and Its Influence on Rigid Body Parts 115 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 model is simplified as necessary. The simplified model eliminates the gearboxes, drive motors, and hydrodynamic coupling in the head and tail drives, and only the head and tail sprocket drums, three chutes, scrapers, link chains, and so on are retained. At the same time, adding a bottom chute, scrapers and chains return to tail from the head through the bottom chute. In the simplified model, the scraper and chain move smoothly in the middle trough driven by the sprockets, which meets the research requirements. 1.2 Simulation Model and Its Parameters 1.2.1 The Coupled Model The process of conveying bulk coal by scraper conveyor involves bulk coal and rigid parts of the scraper conveyor (all parts of the scraper conveyor in this research are rigid bodies). DEM can simulate bulk coal, but the discrete element software cannot simulate the complex motion of rigid body parts, so multi-body dynamics software can be used to simulate the complex motion. Therefore, discrete element software and multi-body dynamics software are used to simulate the transport of coal bulk in the scraper conveyor together. The coupling model is established by using multi-body dynamics software Recursive Dynamic (RecurDyn) and discrete element software EDEM. The simplified virtual prototype model established in UG is imported into RecurDyn, where the dynamics model is established, and then the virtual prototype is exported from RecurDyn to a file in WALL format, and then the file is imported into EDEM as the geometric model to establish the discrete element model. In the process of coupling simulation, the two models run simultaneously, and the calculated data are shared through the WALL file in real time. EDEM calculates the forces and motions of bulk coal, Recurdyn calculates the forces and motions of rigid body parts, and the computational results of the two models are shared in real time to be used as the boundary condition of the next time step, as shown in Fig. 2. This is how the coupling model simulates the interaction between bulk and rigid parts. 1.2.2 Discrete Element Model (1) Contact model The contact model defines the calculation method of forces between particles, particle and geometry. The simulation model needs to define the contact model between particles and between particle and geometry, respectively. Both contact models select the Hertz-Mindlin (no slip) model, which is the basic contact model in EDEM software, which can accurately and efficiently calculate the forces. In the model, normal force is based on the Hertz contact theory [29], and the tangential force is based on the research results of Mindlin-Deresiewicz [30]. Let two spherical particles with radius Ri and R2 make elastic contact, and four forces are generated between the two spheres: normal force, normal damping force, tangential force and tangential damping force. Normal force Fn between particles can be obtained by Eq. (1). 4 ,t ,\i/ Fn = 4 E (*) a (1) where E* [Pa], R* [m] and a [m] are equivalent Young's modulus, equivalent particle radius, and normal overlap, respectively. The normal damping force F/ can be obtained by Eq. (2). Ff = -2 J-BJS m v; (2) where ft is a coefficient associated with the coefficient of restitution, Sn [Nm-1], m* [kg] and vrnel [ms-1] are normal stiffness, equivalent mass and normal components of relative velocity, respectively. Tangential force Ft can be obtained by Eq. (3). F = -S, 8, (3) where St [Nm-1] and S [m] are tangential stiffness and tangential overlap, respectively. Tangential damping force Ff can be obtained by Eq. (4). Ff =-2. ¡6MSm v (4) where ft is a coefficient associated with the coefficient of restitution, St [Nm-1] and v"1 [ms-1] are tangential stiffness and the tangential component of the relative velocity, respectively. Please refer to the literature for the calculation of specific parameters in the above formulas [31]. Fig. 2. Schematic diagram of coupling model 116 Zhang, P. - Li, B. - Wang, X. - Liu, C. - Bi, W. - Ma, H. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 (2) Particle model and particle factory The irregular shape of coal particles can be roughly classified into three categories: flat, block, and cone. In EDEM, coal particle model is composed of pellets. A new particle can be formed by defining the radius and relative position of each pellet in the coordinate system. The maximum length of the coal particle model in three coordinate directions is used to represent the particle size of coal particles. Fig. 3 shows the shape of the coal particle model established in EDEM, a) is block, b) is flat and c) is cone. Coal particles are produced by the particle factory above the middle trough near the tail sprocket. The particle factory is set as a cube area that the length, width, and height are 850 mm, 1000 mm, and 200 mm, respectively. The particle factory can produce coal particles at a certain rate and define the initial state of coal particle, and its position is fixed. Fig. 3. The shape of coal particles; a) block; b) flat; and c) cone (3) Related parameters The parameters in the discrete element model mainly involve the intrinsic parameters of bulk coal and rigid materials and the contact parameters of the contact model. The parameters in the simulation are from the existing research results of this research group, and the specific parameters are shown in Tables 2 and 3 [32]. Table 2. The intrinsic parameters of bulk and rigid materials Material Shear modulus [GPa] Poisson's ratio Density [kg-m-3] Coal 0.47 0.3 1229 Steel 80 0.3 7850 Table 3. Recovery coefficient and friction coefficient Recovery coefficient Coal - coal 0.64 Coal - steel 0.65 Static friction Coal - coal 0.329 Coal - steel 0.46 Coefficient of kinetic friction Coal - coal 0.036 Coal - steel 0.032 1.2.3 Dynamic Model RecurDyn software adopts the motion equation theory of relative coordinate systems and the complete recursion algorithm, which is very suitable for solving large-scale multi-body system dynamics problems. And compared with other dynamic software, RecurDyn can truly, efficiently and stably handle contact collision problem. The scraper chain consists of a large number of chain links and the dynamic model is very complex, and RecurDyn can effectively improve the calculation speed. The interaction between the chainrings, the chainring and the chain wheel and so on is expressed by contact pair, which can effectively simulate the movement of scraper conveyor rigid body parts. In the dynamic model of scraper conveyor, the kinematic pairs of rigid body parts are shown in Fig. 4. The action relationship between parts in the dotted frame in the figure is represented by contact pair, including vertical chain and flat chain, chain and scraper, chain and sprocket, and so on. The head and tail sprockets are fixedly connected to the rack (ground) through the revolute pair; the chute is fixedly connected to the ground with a fixed pair. The driving force of sprocket is defined by the revolute pair; two sprockets drive the chain and scrapers simultaneously. The head sprocket drives the chain and scrapers in the forward stroke, and the tail sprocket drives in the return stroke. Rei/olute\ / Fixed Rack Ground Fig. 4. The kinematic pairs between rigid body parts 1.3 Verification Test of Discrete Element Model To ensure the reliability of simulation results, it is necessary to verify that the coupling model can truly simulate the coal transportation process. If the verification test is carried out directly on the scraper conveyor corresponding to the simulation model, the requirement for the test equipment is high, and the cost is substantial. Since the discrete element model and the multi-body dynamics model that constitutes The Loading Characteristics of Bulk Coal in the Middle Trough and Its Influence on Rigid Body Parts 117 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 the coupling model are independent of each other, it is only necessary to verify that the two models meet the research requirements. Discrete element model is the core of the coupling model, the multi-body dynamics model of scraper conveyor can realize the complex motion of rigid body parts and meets these research requirements. so the emphasis is to verify the reliability of the discrete element model. The core parts of a discrete element model are the contact model, the particle model, and the relevant parameters of bulk material and rigid body, etc. Moreover, the geometric model is not very important. As long as these models are appropriate and key parameters are reliable; then the discrete element model can truly simulate bulk coal. The validation test is to verify the reliability of these models and parameters. The idea of verification testing is using the same contact model, particle model, and related parameters of design simulation and check test on a turntable testbed to prove that the simulation model can accurately simulate the bulk coal by comparing the results of simulation and test. The check test was carried out on the turntable testbed in Fig. 5a. The turntable can rotate around its central axis, and the sample above the turntable is fixed on the rack by a three-dimensional force sensor that can measure the force exerted on the sample by the coal. The bottom surface of the sample is 8 mm away from the surface of the turntable, so there is no friction between the sample and the turntable, and the coal particles are not easily stuck between them. Lay 1 kg bulk coal in the turntable which rotates counter clockwise for 1 cycle at n/2 rad/s. In the test, the coal particles with a size of 2 mm to 4 mm are screened. Furthermore, the particle sizes are randomly distributed between 2 mm to 4 mm in simulation. Rotating with the turntable, the bulk coal exerts thrust on the sample in the circumferential tangential direction, and the sample will leave a groove on the bulk coal. The thrust and the groove section curve cut by the plane perpendicular to the turntable upper surface and across the turntable centre are selected as the test results, which represent the force of bulk coal and bulk coal fluidity, respectively. To avoid contingency, three positions are selected to measure the groove section curves. If the initial position of the sample is 0° and one revolution of the turntable is 360°, the groove section curves at the positions of 90°, 180°, and 270° are respectively taken. Fig. 5b is the schematic diagram of the simulation, and the geometric model is composed of the turntable and sample. The thrust in the simulation is obtained from the WALL file of the sample in RecurDyn. By establishing a two-dimensional coordinate system, coordinate values of the groove section curves at some equidistant points were measured, and the groove section curves can be obtained through curve fitting. Take the intersection point of the vertical centreline of the sample and the plane of the turntable upper-end face as the coordinate origin, take the radial direction of the turntable at the origin as the X-axis, and the vertical centreline of the sample as the Y-axis, and establish a two-dimensional coordinate system, as shown in Fig. 5b. Take a point at 8 mm intervals along both sides of the X-axis, measure the y-coordinate values of the groove section curve at these points, and measure 11 points of data together with the origin, and fit the groove section curve from these data points. The groove section curves in simulation and test are obtained by taking points, measuring data and fitting curves. Fig. 5. Verification test equipment and simulation diagram; a) verification test equipment; and b) simulation diagram Fig. 6a is the diagram of the physical test, and Fig. 6b is the shape of the groove section at the 90° position in the simulation. Fig. 6c is the circumferential tangential thrust curve of the sample exerted by bulk coal. The blue curve is the simulation curve, and the red curve is the test curve. There are some mutation 118 Zhang, P. - Li, B. - Wang, X. - Liu, C. - Bi, W. - Ma, H. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 1.5 2 2.5 Time [s] - simulation fitted cu rve -test fitted curve simulation data points test data points R=0.9979 -25 -40 -32 -24 -16 -B e) 0 B 16 24 32 40 x [mm] b) d) f) Fig. 6. Simulation and test results: a) test physical diagram; b) the shape of the groove section at 90° in simulation; c) the thrust curve; d) 90° groove section curve; e) 180° groove section curve; f) 270° groove section curve points in the blue curve and the mutation generated by the extrusion between the particles flowing under the sample and the sample. In the test, due to the low sampling frequency of the sensor, the sudden change of force cannot be captured or presented, and the force curve reflects the thrust without mutation and corresponds to the relatively smooth part in the simulation curve. Therefore, the method of average value is used to compare the closeness of the two curves. The average value of the red curve is 0.94 N, the average value of the blue curve is 0.82 N, and the ratio of the simulation value to the test value is 87 %. The bulk material force in simulation is close to the real working condition. Fig. 6d to f are groove section curves at 90 180 270 respectively. The small red circles and blue triangles in the figures are data points of groove section curves measured respectively in simulation and test. The smooth spline fitting tool in MATLAB is used to fit these data points to get groove contour curves. All of these fit curves have R2 of 0.9901. The correlation coefficients R of the two fitted curves in simulation and test at the positions of 90°, 180° and 270° are 0.9896, 0.9879 and 0.9881, respectively, which indicates that the fluidity of bulk coal simulated by the simulation model is very similar to that in the test. The results of above comparisons indicate that the discrete element model can realistically simulate the bulk coal and meet the research requirements. c The Loading Characteristics of Bulk Coal in the Middle Trough and Its Influence on Rigid Body Parts 119 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 1.4 The Simulation The load distribution of bulk coal is related to the bulk coal compressive force distribution. The force inside the bulk is transmitted through the compressive force. The greater the bulk coal compressive force is, the greater the load on rigid body parts is, and the greater the friction force applied to rigid body parts is. Firstly, the compressive force distribution gradient of bulk coal in the middle trough is analysed to obtain a general understanding of the bulk coal load distribution. Secondly, the interaction between bulk coal and rigid body parts generates new loads, and the large load generated by violent action significantly affects the normal work of scraper conveyor. Find the coal particles that violently act with rigid body parts, and then analyse the violent action process. Finally, some possible faults of rigid parts are analysed. Fig. 7. The rotation speed of sprockets The bulk coal in the simulation is composed of three shapes of particles. The particle factory produces 220 kg bulk coal every second, and the mass of each shape particle accounts for 1/3 of the total mass. Particle size is selected as 40 mm for three shapes of particles and the total particle quantity is unlimited. Coal particles fall to the middle trough at a certain initial velocity, and the components of the velocity are vx = -0.2 m/s, vy = 0.4 m/s, vz = -1 m/s. The initial tension of the scraper chain can be set by adjusting the distance of two sprockets, and the initial tension is about 20,000 N. The scraper conveyor is driven by the sprockets at the same time as same rotation speed, the rotation speed is shown in Fig. 7, and the total simulation time is 6 s. The process of producing coal at the tail end of the scraper conveyor and transporting it to the head end at a speed of about 1.02 m/s is simulated. 2 RESULTS AND DISCUSSION 2.1 Compressive Force Distribution Gradient of Bulk Coal in Middle Trough The total simulation time is 6 s, the whole middle trough is full of bulk coal at 3.3 s, and the transport of bulk coal in the middle trough is relatively stable at 5 s, so choose the time of 5 s to draw the bulk coal compressive force cloud diagram. Fig. 8a is the bulk coal compressive force cloud diagram cut by the plane perpendicular to the x-axis. The plane perpendicular to the x-axis cuts the middle trough along e-e line in Fig. 8c. And the e-e line is in the middle of the chute inner edge and chain, 375 mm from the spillplate edge. Arrow V indicates the motion direction of the scraper and chain. In the figure, the red bulk coal is mainly located in the chute, which is framed with black lines, while the bulk coal above the chute is mainly blue-green. Therefore, the compressive force of bulk coal in the chute is significantly greater than that above the chute. Fig. 8b is the bulk coal compressive force cloud diagram cut by the plane perpendicular to the z-axis. That plane cuts the middle trough along the f-f line in Fig. 8c. The f-f line is 120 mm away from the upper edge of chute, and the maximum compressive force in the cloud diagram increases to 45 N. The area between the two scrapers in the chute middle is divided into six cubic sub-areas marked as A, B, C, D, E and F respectively, and the cross-section shape of the six cubic sub-areas is shown in Fig. 8b and c. In Fig. 8b, A, B and C are mainly red, D, E and F are mainly blue, and B is less red than A and C, while E is more blue than D and F. Therefore, in the y-direction, the bulk coal compressive force in A, B and C is higher than that in D, E and F, respectively, in the x-direction, the bulk coal compressive force in B is lower than that in A and C, and in E is lower than in D and F. The average compressive force F of the bulk coal in the six areas at 5 s is shown in Table 4. Fd / Fa = 21 %, Fe / Fb = 16 %, FF / FC = 19 %, and the average value of the three ratios is 19%. FB / FA = 79 %, FB / FC = 81 %, Fe/ Fd = 59 %, Fe / FF = 70 %, and the average value of the four ratios is 72 %. Therefore, in the area between two scrapers, in the y-direction, the bulk coal compressive force in the first half is about 19% of that in the second half, and in the x-direction, the bulk coal compressive force between the two chains is about 72 % of that between the chain and the chute edge. Fig. 8c is the bulk coal compressive force cloud diagram cut by the plane perpendicular to the y-axis. The plane perpendicular to the /-axis cuts the middle trough along the g-g line in Fig. 8b. The g-g line is in the 120 Zhang, P. - Li, B. - Wang, X. - Liu, C. - Bi, W. - Ma, H. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 c) Fig. 8. Bulk compressive force cloud diagram at 5 s; a) cut by the plane perpendicular to the x-axis; b) cut by the plane perpendicular to the z-axis; c) cut by the plane perpendicular to the y-axis centre of the chain-link adjacent to the scraper, 126 mm away from the centreline of the adjacent scraper, and the maximum value in the cloud diagram is 45 N. This picture offers a new perspective. Table 4. Bulk average compressive force in the 6 areas at 5 s Area A B C D E F Force, F [N] 87 68 84 18 11 16 The bulk coal compressive force in the chute is greater than that above chute. The bulk coal inside chute is squeezed by scrapers, middle plates, and chutes, and moves forward under the force. The bulk coal above chute moves forward under the friction force of bulk coal in the chute, so the compressive force is obviously smaller than that in the chute. The scraper pushes bulk coal move in the chute, the thrust 121 The Loading Characteristics of Bulk Coal in the Middle Trough and Its Influence on Rigid Body Parts 49 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 T iw : 5 3 Ve]ocity(m/&) Fig. 9. Velocity cloud diagram of bulk coal In Z direction at 5 s is transmitted forward through the bulk, and the rear bulk coal pushes the front, so the rear bulk coal bears greater force. In addition, the middle plate and chute exert backward friction force on coal bulk, and the bulk coal near the former scraper moves relative to the backward scraper, which makes the bulk coal near the former scraper gradually loose, and the compressive force is reduced. while the force near the latter scraper becomes bigger due to aggregation. Fig. 9 shows the bulk coal Z component velocity cloud diagram at 5 s and arrow V indicates the motion direction of the scraper and chain. In the back of each scraper, bulk coal flows downward, while in the front of the scraper, bulk coal basically remains unchanged. This is because the bulk behind each scraper will become increasingly loosed, the bulk coal above will flow down to fill the gap left, the shape of the bulk will be concave down. The distance between chain and chute edge is larger than that between two chains. And the more obvious the effect of bulk coal extrusion is, the greater the compressive force is. Secondly, the influence of chain movement on the bulk coal between two chains weakens the bulk compressive force. 2.2 The Violent Interaction between Bulk Coal and Rigid Body Parts Increase the maximum compressive force value to draw the bulk coal compressive force cloud diagram. The particles with extreme compressive force are all Compressive Force (N) 122 Fig. 10. Bulk compressive force cloud diagram; a) 5 s; b) 5.5 s; and c) 6 s Zhang, P. - Li, B. - Wang, X. - Liu, C. - Bi, W. - Ma, H. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 squeezed by scrapers or chainrings, which indicates that violent interaction happens. The cloud diagram seen from the chute bottom upward can clearly show coal particles squeezed by scrapers, chainrings, and chutes. The maximum value is set to 2000 N, chutes and middle plates are hidden, and chains and scrapers are retained. The bulk coal compressive force cloud diagram of 5 s, 5.5 s, and 6 s are drawn respectively, as shown in Fig. 10. The red particles circled in yellow indicate that these particles are squeezed. And the red particles are mainly located in two places: the bottom of scrapers and vertical chainrings, where the coal particles are squeezed between the scrapers, the chainrings and the middle plates. The compressive force value of 12 coal particles marked in yellow in Fig 10a is shown in Table 5. Ignoring the particle 1 and 4 that the force exceeds 10,000 N, the average compressive force of the rest 10 particles is 3,536 N. However, the average compressive force of all coal particles in the whole chute is only 59 N at 5 s, and the average compressive force of the 10 particles is 59 times as much as that of all particles in the chute. Draw the position relation diagram of extreme compressive force coal particles with scrapers, chainrings and middle plates at 5 s, as shown in Fig. 11. The local enlarged picture of the corresponding particles below is shown above in the figure. The particles 1, 3, 5, and 7 are squeezed between scrapers and middle plates; the remaining particles are squeezed between vertical rings and middle plates, which are almost all squeezed at the bottom of the front halfring of the vertical rings. Most of the scrapers and chainrings in the picture do not contact with middle plates, except for a few chainrings on the left. The left chainrings that just enter the chute are less affected by bulk coal, while the remaining chainrings and scrapers fully interact with the bulk coal and leave the middle plates, which makes the particles easy to squeeze. The severe squeezing of scrapers, chainrings and middle plates with coal particles will affect the movement of scrapers and chainrings, but more importantly, it will cause serious wear of middle plates. Studies and practices have shown that the middle plate under two chains formed grooves due to wear, and the wear was significantly greater than other parts of the middle plate [20]. The above analysis reveals the causes of wear, which can be used to improve the shape of chainrings or other measures to reduce the possibility of coal particles being squeezed, so as to reduce the wear of the middle plate. 2.3 The Clamping Condition of Scrapers Observing the movement of scrapers in the middle trough, a scraper (marked as scraper I) operates abnormally at 4.4 s and 5.9 s that one side of the scraper suddenly stop and the other side still advance; then, the scraper returns to normal operation. Both times, the scraper is located at the connection of two chutes. Scraper I is located at the connection of Chute 1 and Chute 2 at 4.4 s, and is located at the connection of Chute 2 and Chute 3 at 5.9 s. Fig. 12 is the posture diagram of scraper I at these two moments, Figs. 12a and c shows the position relation of the stuck side of the scraper and the previous chute, Figs. 12b and d are the bulk coal compressive force cloud diagram at these two moments seen from the chutes' bottom. The compressive force values are set in accordance with Table 5. Compressive force value of 12 extreme particles at 5 s Particle 1 2 3 4 5 6 7 8 9 10 11 12 Force [N] 103646 2855 2631 18524 5035 1835 7603 1841 4212 3867 1865 3614 Fig. 11. The situation of extreme compressive force coal particles at 5 s The Loading Characteristics of Bulk Coal in the Middle Trough and Its Influence on Rigid Body Parts 123 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 114-126 Fig. 10, the chutes and the middle plates behind are hidden, and the middle plate previous is translucent. The part of the red circle in Fig. 12a shows the scraper is stuck at the top edge of the next chute at 4.4 s, and the local enlarged picture shows the scraper is embedded in the chute. In Fig. 12b, the coal particles at the blocked end of the scraper are not subjected to extreme force, indicating that the scraper is not stuck on the chute through coal particles, but directly stuck on the chute. Figs. 12c and d are the position relation diagram of scraper I and chutes at 5.9 s, which is similar to that at 4.4 s and will not be analysed in detail. J&l R - RP\■ dx. (6) The integral is approximated by the area of a rectangle with the base xj _xi_1 and the height f(£) at the intermediate point ^ e xt - xi-1, the function / being continuous. The above-presented results from the theorem of the defined integral average: \f (x )dx = (b - a )f (£). (7) We have proposed a study of the roughness according to some input variables, the prediction and the correlation formula for Ra. As well, the hardness and the cutting width were studied when cutting with CO2 laser parts of a steel material, less well-known in scientific articles, called HARDOX400 (H400). Laser-cutting technology is used in industrial engineering to manufacture hard steel parts. It is a high productivity and profitability process. The laser-cutting process, roughness optimization according to the laser input parameters and quality evaluation of the cut surfaces were studied in a number of papers. The influence of the laser-cutting parameters on Ra was studied by Mesko et al. [16] for S235JR steel, by Rao and Murthy [17] for stainless steel, by Adarsha Kumar [18] for AISI 4340 steel. Riveiro et al. [19] studied the interaction between the assistant gas and the part as well as the role of the nozzle on Ra for AISI 4340 steel. Spena [20] approached the influence of the cutting parameters on steel sheets used for car bodies. Hirano and Fabbro [21] studied the quality of thick steel sheets according to laser wavelength. Zlamal et al. [22] assessed the quality of S235JR steel, 10 and 15 mm thick, after the laser cutting. Lazov et al. [23] used the regression method to establish the cutting quality influenced by the cutting parameters. Wardhana et al. [24] observed that speed has the greatest influence on Ra for 316L stainless steel. Sharifi si Akbari [25] studied laser cutting for Al6061T6 alloy to determine the effects of the process parameters on the cutting zone temperature and cut edge quality. We propose an analysis of the laser-cutting technology as follows: after irradiating steel material with light from a CO2 laser, the incident photons collide with the valence electrons of the constituent atoms, they easily leave the electronic configuration, becoming free electrons, giving birth to the electronic cloud that moves through steel due to the conducting electrons that strike the metal network ions. At the place where the laser beam reaches the steel plate, the temperature is maximum in the centre, the order is 104 to 106 degrees Celsius, so that the ions of the network begin to vibrate very strongly due to the received heat, and the motion is transmitted to the neighbouring ions, which in turn start to vibrate giving rise to a damped oscillatory motion. In the steel sheet, the heat due to the vibration of the conducting ions and electrons having a wavy behaviour appears, the heat released by laser melting and the heat released by the burning reaction of Fe with O2 is transmitted through the plate through the thermal conduction phenomenon. Of particular importance for laser-material interaction is the characteristics of the laser beam. For the CO2 laser, the spot diameter used in the experiment is 0.2 0 D 130 Mile^an, M. - Gírdu, C.C. - Círtína, L. - Radulescu, C. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 mm. The spot has a spherical shape with the focal spot of microns placed in the centre. Laser cutting is the phenomenon of the thermal processing of metallic materials (HARDOX steel sheets) under the influence of a laser beam. The laser pulse lifetime is of seconds in continuous CW regime. Laser cutting is done by melting the material and removing it with O2 jet. With the penetration of the laser in the H400 steel, the laser is pierced through a circular hole with a diameter of 0.3 mm to 1.25 mm and then cutting in a straight profile in the plate with the cutting width of 0.35 mm to 0.84 mm until it penetrates the contour of the part. In the stationary piercing regime, the material was pulled. The role of the assistant gas is to push the melt to the bottom of the sheet, to increase the volume of the piercing channel and to remove the melt on the upper and lower sides of the steel sheet. On the steel sheet surface, melt deposition surrounding the channel is noted. The assistant gas jet has kinetic energy that forces the melt and slag out of the cutting area. The incident intensity I0 of the laser when cutting has a value in the range 105 to 106. By thermal heating, the material is brought to the melting temperature, forming due to the surface tension a melt in the form of an ellipsoid. The volume of the portion heated to the melting temperature Tm forms the liquid melt. The melt in the liquid state is heated until it reaches the laminar flow state at the melting temperature of Fe due to the laser source and the oxygen assistant gas. Removal of the cut material is done in a liquid (molten) state, in addition to the laser radiation, the oxidation reaction: Fe + O ^ FeO + Q takes place, which is an exo-energetic reaction with the release of Fe oxides. The evaporation energy of the melt is equal to the energy of the melt minus the binding energy of the metal electrons. The cutting process is characterized by the penetration of the material and the occurrence of the oxidation reaction. The mechanism of liquid melt removal by CO2 laser consists of the absorption of photons by the metal surface, transforming the light energy into caloric energy, local melting of the surface with melt formation, removing the melt when the temperature in the centre reaches the melting and vaporization temperature, with the help of compressed oxygen gas, which eliminates particles from the cutting area through an incandescent flow. The cutting has the following constituent subprocesses: sublimation of a small amount of material that is heated to the vaporization temperature Tv, the vapours being evacuated by the compressed gas, elimination of a quantity of the molten material in the channel by a stream of incandescent drops at the melting temperature Tm by the cutting gas, oxidation of Fe in the presence of oxygen giving rise to fine solidified droplets on the lower edge of the cut along the contour of the part separated by the plate, the reaction being more intense towards the lower edge, thus resulting in a higher cutting speed. The role of oxygen is to remove the molten material, trigger and maintain the Fe oxidation reaction, but also to protect the focusing lens L. Purity of O2 is 99.97 %. It directly influences combustion and cutting speed. The cutting surface of the material is formed by a melted and resolidified layer due to the cooling of the piece. In laser cutting, the distance between the nozzle and the part is important because the laser beam and CO2 gas are co-axially transmitted to the workpiece. The parameters used in the cutting are: the distance between the nozzle and the cutting area: 0.8 mm, the focal position: -0.5 mm inside the piece, the radius of the correction instrument: 0.3 mm. The reflection degree of the H400 steel is low, which makes it easier to cut. The low thermal conductivity delays the diffusion of heat through the material, which helps the local melting made by the laser. Using the laser installation, the accuracy is better than in any plasma and oxygen cutting, the laser cutting time is lower, and the quality of the cut surface is significantly higher, the roughness close to the roughness obtained by chip removing process, the profitability is higher than in the case of cutting with plasma or oxygen. An analysis of the surfaces reveals uniform oscillations of level, not detectable by the human eye appear, following the evacuation of the melt from the slit. Incandescent liquid drops can give rise to striations and pores on the surfaces affected by the laser. After laser cutting, the material decomposes into a liquid phase and a vapour phase. All these compete for a much shorter processing time compared to plasma or oxygen cutting, small melting mass, reduced amount of lost material. In this way, technological additions for the part design and manufacturing can be substantially reduced, the use of the semifinished part is more profitable, and the loss of material being incomparably smaller. The pollution of the environment is reduced due to the small slot, short cutting time, material saving, reduced maintenance factors, and a much longer service life of the machine compared to the other two processes. One of the benefits of this operation is that the laser cutting does not change the chemical composition. HAZ is much diminished in comparison to the other two procedures. This cutting process does not influence the flatness, rectilinearity, or shape of Mathematical Modelling Study of Hardox400 Steel Parts' Roughness and Hardness, Cut with CO2 Laser 131 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 the semi-finished product. The number of products obtained per hour with laser technology is much higher than in other processes. The laser beam has a circular section; depending on the diameter D before focusing, it can be defocused with ±Af Defocusing has effects on the material melting and, for steel, it is recommended that this one be inside the material. With the focus inside the material, a decrease of the intensity at the surface of the piece results, the size of the focal spot will increase resulting in the increase of the interaction time between the laser and the material, and thus also the Ra value. In the case of thick steel cutting, for a good cutting, it is necessary that the area of parallelism of the beam z be small, that is, the maximum intensity I0 being kept constant, the melting depth zm being correlated with the diameter of the spot d. The laser spot plays an important role in laser cutting, because beneath it a heated portion of a convex plane shape is produced at the beginning, then, by heating, the plane-convex or biconvex melt and the heated convex meniscus shaped surface under it are formed, and finally the melt and heated area are elongated by heating, resulting in vaporization. The increase of the laser power has the effect of increasing the laser intensity, thus the melting occurs, and if the intensity increases further, the vaporization of the material results. The intensity of the laser beam depends on the radiation energy, and for linear contour cutting a linear energy E and a melting energy Ev are necessary. If the cutting speed of the laser increases, the time of interaction between the laser and the material decreases, which implies smaller dimensions of the molten and heated area, less material melted, reduced time of melt elimination with the help of cutting gas, thus leading to better cutting and melting efficiency. The melt is ejected partially on the sidewalls of the cut surfaces. By eliminating the melt its mass decreases and the roughness is better at the bottom than at the upper area, because by the movement of the melt and of the hot drops of liquid, a series of formations are resulting on the cut surface: striations, irregular or inhomogeneous surface, and the oxide film will cover all this. The number, size, and density of the striations increase with the thickness of the steel. A variable describing the interaction of laser radiation with the processed steel part is the focused radiation energy, defined by the mathematical relation: 4P E =-(8) tt • f -e2' w where P is a laser power, tt interaction time, f convergence distance of the lens, and 6 beam divergence. The cutting efficiency is defined as the ratio between the surface cut in a second and the energy developed by the laser in 1 second (the surface cut with the power of 1 J): v-_g P ' (9) where v is laser speed, g material thickness, and P laser power. The cost parameter defines the cutting speed obtained with a certain laser power: C = - CUT p ' (10) Laser-cutting efficiency is defined as the ratio between cutting efficiency and cost parameter: T = P k, ■ d (11) The efficiency is good if the cutting speed / laser power ratio is lower, i.e., low speed cutting and high laser power. The temperature at the surface of the material can be calculated by the expression: T = P k, ■ d (12) where kt is a thermal conductivity of the material, and d diameter of the focused beam or the laser spot. The vaporized mass can be easily deduced as the mass of the molten material from the slot minus the molten mass of oxides found as slag in the waste tray below the grill of the installation: mv = P ■ Kopi, - Moxizi = P ■ Kerf ' g 'l - M kerf (13) where Vtopit is a volume of the melt, l length of the contour of the cut piece, g material thickness, Kkerf average width of the cut, and p density of the melt. For laser thermal cutting, we used a steel plate made of H400, 10 mm thick, with a length of 475 mm and a width of265 mm, which has a content of: 0.19 % carbon, 1.25 % manganese, 0.70 % chromium, 0.38 % molybdenum, 0.50 % silicon. Initially, we prepared a reduced factorial experiment, involving Taguchi's method for the H400 steel plate. Thus, we made the plan of cutting 45 pieces (Fig. 2). A piece was built with three linear sides and a circular one for CO2 laser cutting. Geometrically, a piece is inscribed in a 40 mm side square. Currently, no research available in the literature has been found that describes the dependence of the roughness on the influence parameters of the process 132 Mile^an, M. - Gírdu, C.C. - Círtína, L. - Radulescu, C. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 Fig. 2. a) Sample plate, b) H400 plate, and c) cut parts when cutting Hardox laser plates. In this context, the objective of the research is to determine in what way the roughness of the surface processed by laser cutting is influenced by the input parameters of the process and what the correlation formula between them is. Thus, we aim to obtain a mathematical formula, the dependence of the roughness on the influencing parameters when at least two parameters act simultaneously, and one is kept constant, using the ANOVA regression model. In a preliminary research, we identified the limit values of the selected parameters for which cutting can exist (at too low laser power, the beam does not pierce the 10 mm thickness H400 plate). The optimum cutting speeds, the assistant gas pressure, and other cutting parameters are specified in Table 1. Table 1. Description of the cutting parameters General parameters Value Unit of measure Group STAHL Material STW22 - Plate thickness 10 mm Cutting mode Contin. regime CW Max. laser power 4.400 W Lens with focal length 190.5 mm Nozzle type NK1515 Focal position -0.5 mm Nozzle diameter 1.5 mm The nozzle-piece distance 0.8 mm Assistant gas 02 - Cutting speed I 1.200 mm/min Cutting speed II 1.400 mm/min Cutting speed III 1.600 mm/min In the experiment we have as varied parameters: laser power (P), pressure of the assistant gas (p), laser cutting speed (v) with the values as it follows: power, P1 = 4100 W, P2 = 4200 W, P3 = 4300 W; pressure of the assistant gas, p = 0.35 bar, p2 = 0.45 bar, p3 = 0.55 bar, and cutting speed, vj = 1200 mm/s, v2 = 1400 mm/s, v3 = 1600 mm/s (Table 2). The references from the reduced factorial plan were replicated four times in order to carry out the statistical research observations regarding the study of Ra roughness. Table 2. Reduced factorial experimental plan No. Part Power [W] Pressure [bar] Speed [mm/min] 1 4100 0.35 1200 2 4100 0.45 1400 3 4100 0.55 1600 4 4200 0.35 1400 5 4200 0.45 1600 6 4200 0.55 1200 7 4300 0.35 1600 8 4300 0.45 1200 9 4300 0.55 1400 The investigation of the technological cutting process is focused on modelling to investigate the process accurately. Experimental modelling establishes the correlation between the input parameters, called the influence factors, and the values measured by different procedures using the measuring instruments. The experiment used a factorial array/ matrix that uses the combination of variable levels, i.e., 33. The minimum level has the logical value -1, the average value 0, and the maximum level has the logical value +1. In order to find the coefficients of a polynomial that contains the influencing factors and the products between them, which signify the interactions between these factors, a number of 32 = 9 tests were performed. The optimization algorithm starts running after establishing the extreme values of the parameters presented in Table 3. Table 3. Input parameters and their extreme values Parameters Minimum Maximum Laser power [W] 4100 4300 Speed [mm/min] 1200 1600 Gas pressure [bar] 0.35 0.55 Focal position [mm] -0.5 -0.5 Mathematical Modelling Study of Hardox400 Steel Parts' Roughness and Hardness, Cut with CO2 Laser 133 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 A series consists of three references: S1: 1, 2, 3; S2: 4, 5, 6; S3: 7, 8, 9. In the case of the first factorial series, the power is kept constant at the value of 4100 W and the pressure changed by +0.10 bar from the value of 0.35 bar, and the speed increased by 200 mm/min starting with the value of 1200 mm/min. The second series is obtained according to a certain rule by modifying and increasing the power by 100 W to the value of 4200 W, the pressure changed by +0.10 bar, from the value 0.35 bar to 0.45 bar and 0.55 bar, and the speed modified with 200 mm/min starting with the value of 1400 mm/min. The third series is obtained according to a certain rule by modifying and increasing the power by 100 W to the value of 4300 W, the pressure changed by +0.10 bar from the value 0.35 bar to 0.45 bar and 0.55 bar, and the speed changed by 200 mm/ min, starting with the value of 1600 mm/min. The three factorial series are replicated four times to the same values according to the above-described model. Previously, the authors studied the laser-steel material interaction, considering that when the laser beam enters the material, it enters a diameter of d = 0.2 mm. The wavelength of CO2 laser radiation is X = 10.6 ^m. The beam diameter before focusing is D = 20 mm. The beam divergence angle is calculated by the Rayleigh formula: law, where R is the reflection coefficient and A is the absorption coefficient. 0 = 1.22—. D (14) Calculating with the above values, we obtain 0 = 0.64 mrad. Radiation power equals the ratio of radiation energy E to laser duration t: p=E. T Radiation intensity in focal spot: I = P = 0 S S T E nd2 4 4E nxd2 (15) (16) We can describe the interaction between radiation and material as follows, in Fig. 3. The reflected beam is RI0 and depends on the reflection coefficient of the material. The beam that penetrates the material has an intensity dependent on the penetration depth and decreases after the gradient dI / dz as z increases: I = (1 - R )Ie (17) The intensity of the radiation that has entered the material decreases according to an exponential Fig. 3. Interaction between laser radiation and substance It must be known that surface roughness influences the absorption of laser radiation. 2 RESULTS The surfaces of the laser-cut H400 steel parts are subjected to a heat treatment. In the place where the laser radiation meets the surface of the material, the temperature is maximum. Iron is heated to the result of the Fea-magnetic 20 °C to 750 °C, Fea-non-magnetic 750 °C to 913°C, Fe Y 913 °C to 1390 °C, at 1200 °C iron ignites in the presence of O2 in an oxidation reaction, Fe 5 1390 °C to 1537°C, liquid iron 1537 °C to 3000 °C, vaporized iron >3000 °C. After the laser action stops, iron is suddenly cooled. The heating temperatures do not coincide with the cooling temperatures; the difference between them gives us the thermal hysteresis. By fast cooling, the steel surface is hardened. It can be improved by metal alloying. In the experiment, we initially located on the test plate the point where the cut piece has the best quality given by the combination of parameters: laser power, cutting speed and gas pressure. Based on this determination, the experiment was subsequently designed.Of course, the piece interacts with the geometry of the conical beam which determines the geometry of the cutting slot. One aspect to note for the beam is that it does not have the same intensity at all points, the distribution being given by the Gaussian curve. The incident intensity I0 is calculated in Table 4. The quality of the material is also given by the chemical composition, a higher concentration of carbon and manganese ensures a surface with low roughness, because carbon and manganese burn well in the presence of oxygen. In contrast, we cannot say the same thing about silicon, it melts harder and produces irregular surfaces (craters, pores. etc.). In order to obtain the roughness of the surface, the liquid melt plays a decisive role. Its 134 Mile^an, M. - Gírdu, C.C. - Círtína, L. - Radulescu, C. T Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 shape is given by the surface tension. In the centre of the melt, the temperature is maximum, and it drops to the edges due to the temperature gradient. Table 4. Radiation intensity I0 for H400 parts Parameter Exp. 1 Exp. 2 Exp. 3 Laser durat. CW (r), [s] ri = 14 s r2 = 13 s r3 = 12 s Power, [W] 4100 4200 4300 Radiat. ener. (E), [J] 57.400 54.600 51.600 Radiat. intensity I0, [W/cm2] 9.76x105 10.76x105 14.33x105 It follows that in the droplet, a displacement of liquid metal occurs from the focal point towards the edge, which gives us the concave shape of the droplet. This is removed from the channel with the assistant gas jet blowing and breaking the drop into two fronts containing liquid metal droplets that are ejected through the bottom edge of the plate. The molten droplets give rise to striations, and the craters are formed due to the oxidation reaction and appear on the oxidized surface. By cooling, the solidification of the slit surfaces takes place, resulting in a relief with irregularities. giving rise to the appearance of roughness. From the results of the cutting experiment (Table 5), it is found that the cutting width increases at the bottom, this shows the conicity of the piece at each reference W > Ws. This result is important because the molten mass increases at the bottom, so the number of molten drops also increases due to the oxidation chemical reaction. The best cutting values in the initial experiment are at a minimum power of 4100 W and maximum values of pressure and speed, Ws = 0.6220 mm. The same result is obtained at part no. 39. For the measurement of the parts, an electronic micrometer was used. Hardness was measured with a Krautkramer MIC 2V durometer. The part hardness is given by the high concentration of carbon, which are higher concentration results in higher hardness. For the HARDOX material, hardness varies from the initial value of 29 HRC. The hardness of the piece in the initial experiment is maximum for piece 1, which is 30.80 HRC, so at the minimum values of the power which is 4100 W, the pressure of 0.35 bar and the speed of 1200 mm/s. Subsequently, after heating the material and transmitting heat to the plate, the hardness randomly varies with smaller values obtained compared to part 1, which does not comply with a calculation formula. The lowest hardness value is 23.90 HRC and occurs at maximum power, maximum speed and minimum pressure. The best hardness is obtained for parts 11 and 27 of 31.60 HRC made at low power, medium pressure and speed, respectively at maximum power and pressure and average speed. The roughness was measured in the laboratory of Constantin Brancusi University with a Mytutoyo rugosimeter. The minimum value of Ra 1.7 ^m, is obtained for the 45 pieces at the average value of pressure and maximum cutting speed when the power is kept constant at the value of 4200 W. The graph in Fig. 4 shows a variation of the hardness according to Ra for the proposed experiment. The maximum hardness values are obtained for Ra in the range of 2 |im to 3 |im. From Fig. 5 it follows that at higher values of Ra we obtain lower values of density. The measured HRC points can be approximated by a sine function. Table 6 presents the statistical data obtained by running the program referring at the average roughness and the standard deviation of the roughness as well as the value of the limit coefficients (upper and lower). Table 6 shows Design 33 with 3 levels of laser power variations, assistant gas pressure, cutting speed, Table 5. Values of the measured parameters; the width of the cut up-down, hardness, roughness (replica no. 4) Experiment No. No. crt Ref. Laser power [W] Gas press. O2 [bar] Cutting speed [mm/min] Kerf-[mm] Kerf-W, [mm] Hardness down part [HRC] Ra M 37 4100 0.35 1200 0.5755 0.6880 28.00 3.26 38 4100 0.45 1400 0.6575 0.6750 30.10 2.87 39 4100 0.55 1600 0.6720 0.6765 27.60 2.56 40 4200 0.35 1400 0.5595 0.6870 28.10 2.09 Replica 4 41 4200 0.45 1600 0.6495 0.6500 28.20 3.46 42 4200 0.55 1200 0.6430 0.7615 28.10 4.05 43 4300 0.35 1600 0.5755 0.6865 28.30 3.75 44 4300 0.45 1200 0.7015 0.7295 29.90 2.56 45 4300 0.55 1400 0.6895 0.6980 27.30 3.63 Mathematical Modelling Study of Hardox400 Steel Parts' Roughness and Hardness, Cut with CO2 Laser 135 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 1 block with 45 operated parts where the roughness Ra was measured. The DOE took into account an experiment reduced to 9 observations, replicated at the same points four times. The central point of the experiment was chosen after testing on the H400 test plate at P = 4200 W, p = 0.65 bar, v = 1400 mm/min 4 4.5 A. [foilJ Fig. 4. Variation of hardness depending of the roughness Ra for the initial experiment The minimum and maximum levels were chosen for the input parameters with ±100 W for laser power PM = 4300 W, Pm = 4100 W, ±0.10 bar for the cutting gas pressure pM = 0.55 bar, pm = 0.35 bar, the speed of laser cutting ±200 mm/min, vM = 600 mm/min, vm = 1200 mm/min. These values varied according to the following rule: The power was maintained constant at 4100 W and the pressure was run at the minimum, average, maximum value (0.35 bar, 0.45 bar, 0.55 bar) and laser-cutting speed at the value minimum, average, maximum (1200 mm/min, 1400 mm/min, 1600 mm/min). The power was increased to the constant value of 4200 W and the pressure was run to the minimum, average, maximum value (0.35 bar, 0.45 bar, 0.55 bar), respectively: the speed was increased by 200 mm/min from the value 1400 mm/min, 1600 mm/min, 1200 mm/min. The power was increased to the constant value of 4300 W and the pressure was run to the minimum, average, maximum value (0.35 bar, 0.45 bar, 0.55 bar), respectively the speed started from 1600 mm/min, after which it ran at 1200 mm/ min and at 1400 mm/min. The average roughness was calculated for 5 parts run under observation conditions 1, technologically obtained under the same conditions, and was found: As {f 0»} Fig. 5. Variation of hardness depending of the roughness Ra for Replica 4 Similarly, it is calculated for the other series resulting in a column the average of 3.4142 for the roughness of the 45 pieces. For the standard deviation, the classical formula for n = 5 is applied, and carrying out the square sum SS, a standard error of 0.3812 results for series 1. Calculating the roughness average in the column, summing all the values and averaging Table 6. Statistical data referring to the independent variables and dependant variable: mean and standard deviation (mean square deviation) Design: 3 3-level factors, 1 Blocks, 45 Runs (Rugosity)_ Power Pressure Speed Rugosity - Means Rugosity - Rugosity - -95,% - +95,% - [W] [bar] [mm/min] [pm] Std. Dev. N Cnf.Limt Cnf.Limt 1 4100 0.35 1200 3.2780 0.3812 5 2.8046 3.7513 2 4100 0.45 1400 2.9920 0.5495 5 2.3096 3.6743 3 4100 0.55 1600 3.4580 0.7853 5 2.4828 4.4331 4 4200 0.35 1400 3.1180 0.8537 5 2.0579 4.1780 5 4200 0.45 1600 2.8520 0.8366 5 1.8131 3.8908 6 4200 0.55 1200 4.3360 0.8347 5 3.2994 5.3725 7 4300 0.35 1600 3.5120 0.5962 5 2.7717 4.2522 8 4300 0.45 1200 3.8420 0.8521 5 2.7838 4.9001 9 4300 0.55 1400 3.3400 0.2440 5 3.0369 3.6430 All Runs 3.4142 0.7624 45 3.1851 3.6432 140 M/7e§an, M. - Girdu, C.C. - Cfrtfna, L. - Racfu/escu, C. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 the roughness values for all the parts. we find the value 3.4142 micrometres for the whole block. In the column of standard deviation for roughness, we have the standard error values for each series 1 to 9 resulting in an average square error of 0.7624. N = 45 represents the rank of the block. R = 3-48 + 3'81 + + + 326 = 3.2780. (18) a 5 We can calculate the confidence limit, taking into account the confidence interval Z = 1.96 The calculated in the table data, respectively the average roughness, Z, standard deviation, the mean of the standard deviation with n = 5, will be taken into account. With these values, we can easily calculate the t-student distribution for eliminating aberrant values. The Ra values distributed as a Gauss bell are found in the interval t and -t. Error distribution, probability p is made in the range -95 % and -t, respectively t and +95 %. 3 DISCUSSION Response surfaces (RSM) are study methods that indicate results in an accessible form. RSM presents a 3D graph of the measured value Ra based on two influence factors. With the help of the Statistica.7 software, the correlation formula between the influencing factors was established. The interaction between the influencing factors describes the hierarchy of the interactions between them. There is statistical software in which the influence factors in real values are used in the approximation of the correlation formula (SI), resulting in the comparison. The magnitude of the effect of a parameter is given by the polynomial coefficient module. 1. Prediction: Vary power and speed, pressure remains constant, a. For surface plot (fitted response): The graph shows that at high cutting speeds 1550 mm/s to 1600 mm/s low values of roughness Ra < 3 ^m are obtained when the laser power is 4080 W to 4120 W. and its high values are obtained when the speeds are low, so when the laser stays longer in the material by strongly heating the plate and the cut pieces. In this case, it can be observed that the lowest values of the roughness are obtained at low powers. At low speeds and laser power near maximum, the roughness Ra increases over 4 ^m. The prediction indicates the colour distribution according to the roughness values for the contour plot and shows that the variation of the laser speed and power is inversely Fig. 6. The influence of laser power and pressure on the roughness of the piece when the pressure is constant for surface plot (fitted response) 1850 16M> 1SSO two 145C 14» ose 1300 USD 1100 IIS0 IIOO 4120 4140 4IW 4)80 4200 1120 4240 42ft) J280 HOC 4320 Power ¡W] •4 ■3.5 :3 Fig. 7. The influence of laser power and pressure on the roughness of the piece when the pressure is constant for contour plot (fitted response) Fig. 8. Machining the surface roughness against power and speed with the linear model proportional. Ra increases linearly with the decreasing of cutting speed and the increasing of laser power. The graph shows the measured values of Ra using the Mathematical Modelling Study of Hardox400 Steel Parts' Roughness and Hardness, Cut with CO2 Laser 137 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 SSM method through a fixed response. Using quadric SSM, we indicate Ra which is simultaneously under the influence factors speed and laser power. b. For contour plot (fitted response). 2. The correlation formula of Ra when the factors of influence power and speed vary, and the pressure remains constant. Correlation formula for the linear model (where X is speed, and Y power): Ra =-1.4414 - 0.0014X + 0.00167. (19) Ra is small at high cutting speed and low power value, Ra < 3 ^m. Ra is maximum for low speed and high laser power. The formula shows the increase of Ra under the more pronounced influence of the laser power. which has the coefficient of the higher influence factor. The result given by linear model is in accordance with the prediction. The result indicates that the predominant factor is power. Fig. 9. Machining the surface roughness against power and speed using the quadric model Ra = -75.0092 - 0.0023X + 0.03717 + 9.9083E - -6X•X-6.375E-6X7-3.1667E-677. (20) For expressing Ra as a function of two variables X, Y (laser power and speed), we used the quadratic polynomial model containing all the lower-order hierarchical terms, terms to be added to the model involving interactions, quadratic terms with regression coefficients. The equation defines the response surface, which is how Ra depends on X and Y, which is a surface in a multidimensional graph. In the equation above, E is a mathematical constant. DOE avoids problems caused by correlated predictors. The coefficients can be obtained for calculating effects and interactions. The graph shows the curvature due to each variable X and Y for the designed experiment. The quadric model is a complete model completing the linear model that does not fit the data. Both models are mounted and analysed in the same way. The modulus of the correlation coefficient means that the order of influence is: 1. laser power, and 2. speed. Eq. (20) indicates a pronounced dependence when the factors of influence speed and pressure interact with each other, but also of the dependence given by the quadratic coefficients for speed and power respectively. The graph shown by the correlation formula indicates that the best roughness Ra is obtained at average speed values, when the speed decreases between 1500 mm/min to 1400 mm/min and the laser power increases from 4100 W to 4140 W. Ra changes increases according to Gauss bell maintaining the speed between 1400 mm/min and 1500 mm/min but increasing the laser power from 4200 W to 4300 W using the quadric model. The influence due to the linear coefficients in the quadratic polynomial formula of Ra is very small in comparison to the quadratic coefficients due to the interaction with each one, respectively with the coefficient due to the interaction between the influence factors (speed and power). The coefficient of the linear power term is greater than the speed. The coefficient of the model indicates a decrease of the roughness Ra with a value of 75.0092. The result consists of identifying correlations that include the variations of the influencing factors that approximate the values of the measured parameters obtained from the experiment and which contain factorial series, giving rise to a mathematical model of prediction and establishing the correlation function. It is found that the prediction is in accordance with the correlation formulas of Ra. 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 X Category Boundary Fig. 10. Histogram of predicted values. 3 3 level factors, 1 Blocks, 45 Runs; MS Residual =0.48123 140 M/7e§an, M. - Girdu, C.C. - Cfrtfna, L. - Racfu/escu, C. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 Analysing the histogram, the prediction of the roughness values Ra is shown in Fig. 10 with the arithmetic mean of 3.4 and a standard deviation of 0.7624. Xcategory represents the measured value of Ra. The histogram of the frequencies derives from the values of the pieces n = 45, determinations with a Ra of minimum 2.8 and maximum 4.4, the distribution of these values being slightly asymmetrical, this being due to the cutting defects (striations, pores etc.). There are six intervals, with the following frequencies of appearance of the predicted Ra values: in the range (2.8 |m to 3.0 |m) there are n1 = 5 values; in the range (3.0 |im to 3.2 |im) there are n2 = 5 values; in the range (3.2 |im to 3.4 |im) there are n3 = 8 values, in the range (3.4 |m to 3.6 |m) there are n4 = 8 values, in the range (3.8 |m to 4.0 |m) there are n5 = 5 values and in the range (4.2 |m to 4.4 |m) there are n6 = 2 values. Most Ra values are grouped in the interval (2.8 to 3.6), close to the average. 2.5 2.0 1.5 1.0 0.5 0.0 -0.5 -1.0 -15 -2.0 -2.5 -3.0 *§ 1 1 ^ is f i 1 ! 1 /0 « .95 .75 .55 .¡S .15 .05 -10 -U -10 -05 00 0.5 Residual Fig. 11. The graph of the linear regression of the normal values expectations as function of errors From the graph, one can observe that not all error data points are on the line, most error points are. The total sum of squares that are not on the line, explained by the linear regression model, is the residual sum of squares. In the case of the experiment, the average of the sum of residual squares is MSR = 0.48123. Table 7. Final statistics A SSEI /2 P 0.719563 17.19110 0.347263 0.555671 Table 7 shows the statistical coefficients of the errors for the measurement of the dependent variable, Ra. From Table 7, it follows that the standard error mean is 0.48123. Significance test x2 Pearson indicates the association of columns and lines of a table with two entries: cutting speed and laser power. For our case, the significance of the association is 0.347263. Test x2 shows whether or not the two variables are independent. Finally, a comparison can be made between the significance test x2 (calculated) and the critical test x2, if there is a connection between speed and laser power. So, we can see that 0.347263 < 0.48123, H0 (the null hypothesis) results, i.e. there is no significant link between the speed and the laser power, independent variables B and C. I coefficient is a measure of prediction error decreasing when the independent variables are known, power and speed. I = 0.719563, a value between 0 and 1, so the standard error can be estimated depending on the standard deviation and the root of the block rank n, where the coefficient I is better the higher its value. The statistics are asymmetric, where the dependence variable is Ra. The coefficient I called Goodman and Krushkal measures the reduction of the Ra error in the analysis of crossings for the laser power and cutting speed factors. SSEI = 17.19110 indicates the sum of the squares of errors. Test p indicates that the probability of errors is 0.555671, which is the probability of making an error. As p > 0.05 shows that the statistical connection is insignificant, the power and the speed do not have links between them, i.e., they are independent, so the probability of making an error is quite high. 4 CONCLUSIONS In this article the roughness Ra of the part surface was studied, in the case of laser cutting of an H400 steel sheet. CO2 laser-cutting experiments were performed with the following input cutting parameters: laser power 4100 W to 4300 W, oxygen pressure of 0.35 bar to 0.55 bar and cutting speed of 1200 mm/min to 1600 mm/min. The output parameter was the roughness of the surfaces resulted in the cutting process. Experimental data were statistically analysed by using Anova from Statistics 7.0, in order to establish the linear and quadratic correlation formulas between the influence factors: power, speed, and constant pressure on the one hand and roughness on the other and to determine the predictive model, as well. The main conclusions of the research are: • Following the data processing, the quadratic predictive model shows that at high speed and at low laser power values, low values of roughness, Mathematical Modelling Study of Hardox400 Steel Parts' Roughness and Hardness, Cut with CO2 Laser 127 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 below 3.0 |.un. are obtained, while in the linear predictive model, even if the power increases, the roughness is maintained at same value if speed does not vary. Using the quadric SSM model, the correlation between Ra and the two factors of influence (speed and laser power) is obtained. Thus, it can be concluded that at medium speed and low power, the roughness is improved (low values), and at low speeds and high laser power results high values for Ra. Using the linear SSM model, we obtain the correlation between Ra and the two factors, speed and laser power, proving that the roughness Ra has low values at high speeds and low power, and high values at low speed and high power. It is also deduced from the cutting experiments of the H400 material as well as from the formulas for determining the roughness A',, that laser power is the most sensitive input parameter. Increasing laser power results in an increase of the laser focusing energy; therefore, the heat received by the material has the effect of increasing Ra and, for this reason, the minimum value of the laser power is recommended in order to obtain a low roughness value. • Ra presents a relatively good agreement between the regression model and the predictive model. The representation of the quadric model indicates the variation of the investigated factors of influence, speed and power according to Gauss's function. Ra is minimum at average values of speed and at the lowest power. One explanation would be that the number of drops increases towards the lower part of the H400 steel plate, due to the oxidation reaction, resulting from the increase of the molten mass. Ra increases with the elimination of hot drops that give rise to an irregular surface. As the power increases, Ra increases, which is a coefficient of influence that establishes the hierarchy in the correlation formulas. The result of the interaction between the influence factors delimits the degree of roughness of the surface indicated by the intensity of the colour green. From the separate graphical analysis of the quadric model, we can see that Ra varies directly proportional to the laser power and exponentially with the speed. Outside the bell, the roughness increases exponentially with speed for parts that do not require further processing. By elongating the dependence curve between power and speed, the peak of the Gauss bell moves to higher values of roughness. This means that the temperature inside the slot has increased. Therefore, it is recommended to use minimum power and medium speed for the critical parts of the piece with the potential to be damaged. For a better understanding of the relation between power and speed, we propose equating the linear relationship with the quadratic one in order to observe this dependence, which must be an exponential function. If the exponent of the relation grows, the peak of the bell moves towards a larger Ra. The area of the figure below the graph shows the level of roughness of the surface. The values used in the experiment can be considered optimal for the cutting machine operator and the process engineer when they target the roughness of the cut part surfaces. For example, to achieve a cutting roughness under 3 |im and for 10 mm thickness of material, the laser power value should be 4100 W while the cutting speed should be 1400 mm/min, at a low, constant, gas pressure. 5 ACKNOWLEDGEMENTS Thank you to the manager Radu NUTU and the laser machine operator Liviu from Brasov for their technical support in conducting this C02 laser experiment. 6 REFERENCES [1] Adelmann, B, Hellmann, R. (2011). Fast laser cutting optimization algorithm. Physics Procedia, vol. 12 p. 591-598, D0l:10.1016/j.phpro.2011.03.075 [2] Lutey, A.H.A., Ascarl, A, Fortunato, A, Romoll, L. (2018). Longpulse quasl-CW laser cutting of metals. International Journal of Advanced Manufacturing Technology, vol.94, p 155-162, D0l:10.1007/s00170-017-0913-x [3] Ivarson, A, Powell, J, Slltanen, J. (2015). Influence of alloying elements on the laser cutting process. Physics Procedia, vol. 78, p 84 88, D0l:10.1016/j.phpro.2015.11.020 [4] Pocornl, J, Powell, J, Frostevarg, J, Kaplan, A.F.H. (2018). Dynamic laser piercing of thick section metals. Optics and Lasers In Engineering, vol. 100, p. 82-89, D0l:10.1016/j. optlaseng.2017.07.012 [5] Shulyatyev, Victor B.; Orlshlch, Anatoly M, (2018), Mlcrocraters and surface quality In laser oxygen cutting of thick steel sheets. Journal of Laser Applications, vol. 30, no. 2, art ID 022003, D0l:10.2351/1.5008798 [6] Thombansen, U., Hermanns, T., Stoyanov, S. (2014). Setup and maintenance of manufacturing quality In C02 laser cutting. Procedia CIRP, vol. 20, p. 98-102, D0l:10.1016/j. procir.2014.05.037 [7] Zhang, Y.-L, Lei, J.-H. (2017). Prediction of laser cutting roughness In Intelligent manufacturing mode based on ANFIS. 140 M/7e§an, M. - Girdu, C.C. - Cfrtfna, L. - Racfu/escu, C. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 127-141 Procedia Engineering, vol. 174, p. 82-89, D0l:10.1016/j. proeng.2017.01.152. [8] Zhang, K., Guo, X., Sun, L., Meng, X., Xing, Y. (2019). Fabrication of coated tool with femtosecond laser pretreatment and its cutting performance in dry machining SLM-produced stainless steel. Journal of Manufacturing Processes, vol. 42, p. 28-40, D0I:10.1016/j.jmapro.2019.04.009. [9] Wang, T., Li, Y., Liu, J.,Qin, L., Wang, N., Zhang, L., Wang, H., LI, Z. (2019). Milling force and surface topography of Ti-6Al-4V titanium alloy cladded by the laser. Surface Review and Letters, vol. 26, no. 5, art ID 1850185, D0I:10.1142/ S0218625X18501858. [10] Yi, J.H., Kang, J.W., Wang, T.J., Wang, X., Hu, Y.Y., Feng, T., Wu, P.Y. (2019). Effect of laser energy density on the microstructure, mechanical properties, and deformation of Inconel 718 samples fabricated by selective laser melting. Journal of Alloys and Compounds, vol. 786, p. 481-488, D0I:10.1016/j.jallcom.2019.01.377. [11] Hatala, M., Duplák, J., Dupláková, D., Botki, F. (2019). Effect of traverse speed on surface roughness parameters after laser cutting of non-alloy structural steel. Technology Education Management Informatics Journal, vol. 8, no. 2, p. 402-408, D0I:10.18421/TEM82-12. [12] Kim, D., Lee, S., Park, B. H., Kang, S. (2019). Analysis of the effects of supersonic assist gas for laser cutting using normal shock theory. Transactions of the Korean Society of Mechanical Engineers B, vol. 43, no. 4, p. 231-239, D0I:10.3795/KSME-B.2019.43.4.231. [13] Feng, Y., Hung, T.-P., Lu, Y.-T., Lin, Y.-F., Hsu, F.-C., Lin, C.-F., Lu, Y.-C., Lu, X., Liang., S.Y. (2019). Surface roughness modeling in laser-assisted end milling of Inconel 718. Machining Science and Technology, vol. 23, no. 4, p. 650-668, D0I:10.1080/10910344.2019.1575407. [14] Masoudi, S., Mirabdolahi, M., Dayyani, M., Jafarian, F., Vafadar, A., Dorali, R.M. (2019). Development of an intelligent model to optimize heat-affected zone, kerf, and roughness in 309 stainless steel plasma cutting by using experimental results. Materials and Manufacturing Processes, vol. 34, no. 3, p. 345-356, D0I:10.1080/10426914.2018.1532579. [15] Dragu, D., Popescu, I., Sturzu, A. (1980). Tolerances and Technical Measurements, Didactica and Pedagogica Publishing House Bucharest. (in Romanian) [16] Mesko, J., Zrak, A., Nigrovic, R., Nikolic, R. R. (2018). The effect of selected technological parameters of laser cutting on the cut surface roughness. Tehnicki vjesnik - Technical Gazette, vol. 25, no. 4, p. 997-1003, DOI:10.17559/TV-20160609171348. [17] Rao, V.K., Murthy, P.B.G.S.N. (2018). Modeling and optimization of tool vibration and surface roughness in boring of steel using RSM, ANN and SVM. Journal of Intelligent Manufacturing, vol. 29, no. 7, 1533-1543, D0I:10.1007/ s10845-016-1197-y. [18] Adarsha Kumar, K., Ratnam, C., Venkata, Rao, K., Murthy, B.S.N. (2019). Experimental studies of machining parameters on surface roughness, flank wear, cutting forces and work piece vibration in boring of AISI 4340 steels: modelling and optimization approach. SN Applied Sciences, vol. 1, art ID 26, DOI:10.1007/s42452-018-0026-7. [19] Riveiro, A., Quintero, F., Boutinguiza, M., del Val, J.,Comesana, R., Lusquinos, F., Pou, J. (2019). Laser cutting: A review on the influence of assist gas. Materials, vol. 12, no. 1, art. ID 157, D0I:10.3390/ma12010157. [20] Spena, P.R. (2017). CO2 laser cutting of hot stamping boron steel sheets. Metals, vol. 7, no. 11, art ID 456, D0I:10.3390/ met7110456. [21] Hirano, K., Fabbro, R. (2012). Possible explanations for different surface quality in laser cutting with 1 and 10 |jm beams. Journal of Laser Applications, vol. 24, art ID 012006, D0I:10.2351/1.3672477. [22] Zlamal, T., Malotova, S., Petru, J., Brytan, Z., Musil, V. (2018). The evaluation of the surface quality after laser cutting. Innovative Technologies in Engineering Production, vol. 244, art ID 02009, D0I:10.1051/matecconf/201824402009. [23] Lazov, L., Nikolic, V, Jovic, S., Milovancevic, M., Deneva, H., Teirumenieka, E., Arsic, N. (2018). Evaluation of laser cutting process with auxiliary gas pressure by soft computing approach. Infrared Physics & Technology, vol. 91, p. 137-141, D0I:10.1016/j.infrared.2018.04.007. [24] Wardhana, B.S., Anam, K., Ogana, R.M., Kurniawan, A. (2019). Laser cutting parameters effect on 316L stainless steel surface. IOP Conference Series: Materials Science and Engineering, vol. 494, art ID 012041, D0I:10.1088/1757-899X/494/1/012041. [25] Sharifi, M., Akbari, M. (2019). Experimental investigation of the effect of process parameters on cutting region temperature and cutting edge quality in laser cutting of AL6061T6 alloy. Optik, vol. 184, p. 457-463, D0I:10.1016/j.ijleo.2019.04.105. Mathematical Modelling Study of Hardox400 Steel Parts' Roughness and Hardness, Cut with CO2 Laser 69 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2019.6313 Short Scientific Paper Received for review: 2019-09-03 Received revised form: 2019-11-18 Accepted for publication: 2019-12-19 Evaluation of Clogged Hydropower Plant Trash Rack Losses Aleš Hribernik* University of Maribor, Faculty of Mechanical Engineering, Slovenia A trash rack is applied in front of the turbine to restrict the entrance of significantly sized material present in the water. It obstructs the free flow, and produces energy-losses by generating eddies induced partially by the trash rack bars and partially by the debris collected on it. While the additional static forces due to debris accumulation are considered in the trash rack design process, the debris caused energy losses taking place during plant operation are usually neglected, although a rather simple model was developed to account for them. However, the long term application of this model demands an extensive set of trash rack clogging data and, therefore, no such application has been documented so far. Thus, an analysis was performed to acquire the debris accumulation intensity with time and to evaluate the extra energy losses they caused. Data on one year operation of a hydropower plant aggregate i.e. flow rate and trash rack head losses measured at 15 minute intervals, were acquired and used to build a rack clogging model. Using this model, it was possible to distinguish clearly between debris and rack structure caused head losses, and to analyse different cleaning strategies. It was shown that cumulative debris contributed to almost one half of head losses although the rack was cleaned frequently. This shows clearly that debris caused head losses may not be neglected, moreover, debris removal has to be planned carefully and carried out efficiently. Analyses of acquired data confirmed that incomplete debris removal increased head losses by 18 %, and proved how important regular and thorough rake cleaning is. Moreover, it was found out that the actually applied periodical rake cleaning was not optimal, and that the circumstances required cleaning strategy performed much better. It resulted in similar head losses, while the number of rack cleanings was reduced by 60 %. Keywords: hydropower plant, trash rack, debris, head losses Highlights • A trash rack clogging model was developed using the time dependent debris collection rate acquired from a measured head loss-time series. • Debris caused head losses contribute substantially to overall head losses of a trash rack, even when the rack is cleaned frequently. • Different trash rack cleaning strategies were proposed and analysed from both the technical and economical points of view. • A circumstances required cleaning strategy clearly surpasses the commonly applied periodically performed trash rack cleaning by reducing number of rack cleanings significantly at the same head losses. 0 INTRODUCTION The application of alternative, and especially renewable, energy sources is growing from year to year. It is stimulated by raising the prices of conventional energy sources, the striving of individual countries to reduce their dependence on imported energy sources, and implementing the Kyoto Protocol Directives for the reduction of global emissions from greenhouse gases [1]. Hydropower has a very high potential regarding renewable energy sources, and its share in total electricity production from renewable energy sources exceeds 90 % in many countries [2]; however, the remaining potential of high power Hydro Power Plants (HPPs) is very limited, and the possible investors are now focused on small units and other renewable energy sources. Any efficiency improvement of the existing HPP operation also counts to this category, since it increases energy production without any further impact on the environment and with minimal investment cost. Improving the cleaning strategy of HPP trash racks is a simple example of such measures. 142 *Corr. Author's Address: University of Maribor, Trash racks are a vital part of each HPP. They are installed in front of the turbine in order to restrict the entrance of any material present in the water (like drifting debris, ice or trash), which can damage vital parts of a HPP. A substantial amount of debris, ice and trash drifting in a river can damage vital parts of a HPP. Trash racks are, therefore, used to restrict the entrance of significantly sized material present in the water. The trash rack obstructs the free flow and produces energy-losses by generating large-scale flow structures or eddies, induced partially by the trash rack bars and partially by the debris collected on them [3]. The latter can be reduced significantly when the trash rack is cleaned regularly. Debris is usually removed from a rack by raking. Mechanised rakes are used on large HPPs, while much cheaper manually operated rakes are usually applied on medium and mini size HPPs [4]. Thus, rakes" operation can be fully automated, or they can be operated manually. Application of rakes does not disturb the plant operation much, and racks may be cleaned as frequently as needed, which is definitely the fact when the mechanised racks are used. However, when the manually operated rakes are Faculty of Mechanical Engineering, Smetanova 17, Maribor, Slovenia, ales.hribernik@um.si Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 used a team of workers is needed, and their frequent interventions may increase the operation cost of the HPP significantly. It is, therefore, very important to find an optimal strategy for trash rack cleaning in order to keep the cleaning costs low to moderate on one hand, and not to increase too much the energy losses caused by the collected debris on the other hand. An extensive number of references on trash rack design and trash rack energy losses can be found in the literature [5] to [14], however, the authors mostly neglected the debris caused losses, and only Meusburger [7] considered them in his work. This lack of investigation on debris caused losses is surprising, especially nowadays, when the migratory fish protection law dictates new, fish friendly trash rack design which reduces the clear spacing between the rack bars extensively (values of 10 mm to 30 mm are recommended [15]), as well as lowers the rack inclination below 30° relative to the horizontal [16]. These measures, especially the clear spacing reduction, increase the trash rack losses [11] and [17], as well as the danger of clogging. Meusburger [7] developed a rather simple model to evaluate clogged trash rack losses. However, the application of his model demands an extensive set of trash rack clogging data and, therefore, no such application has been documented so far. Thus, an analysis was performed to acquire the debris accumulation intensity with time and to evaluate any extra energy losses they caused. In the first step, an analysis was performed to evaluate extra energy losses caused by the debris collected on the trash rack over several months" operations. Our goal was to distinguish clearly between the contribution of the rack itself and collected debris to the total energy losses of a clogged trash rack. Data on one year's operation of the selected HPP unit, i.e. flow rate and trash rack head losses measured at 15 minute intervals, were used to calculate cumulative energy losses. Since the aggregate was under general refit at the beginning of the year and its trash rack was inspected, repaired and cleaned manually, measured head losses of the trash rack during the first days of operation corresponded to clean, debris free operation, and were applied to build a (clean) rack energy losses model. When applied later in the year, this model predicted clean rack head losses, and the debris contributed head losses were estimated by deducting them from the measured total head losses. It was shown that cumulative debris contributed to 48 % of total trash rack head losses, despite the fact that the trash rack was cleaned frequently. This shows clearly that debris caused head losses may not be omitted, moreover, they may cause an important reduction in electricity production when the debris removal is not managed correctly. In order to manage debris removal properly, one has to compare different trash rack cleaning strategies, and any trash rack clogging simulation model may be of great support. Such a model should predict the realistic debris accumulation, i.e. trash rack clogging with time under any debris removal frequency, and may allow total head loss calculation. The application of correct, i.e. realistic, debris accumulation rate is crucial here; it should mimic actual debris flow in the river. Debris flow intensity is not constant, and may change significantly within both short and long-time scales, i.e. from day to day and seasonally. It can be obtained from the long-term measured data on HPP operation, as already discussed. As mentioned, it was possible to separate the head losses caused by the debris from those induced by the trash rack structure. This resulted in a trash rack clogging rate model which was used to predict the rate of debris accumulation on the trash rack. It was, therefore, very simple to predict the instantaneous area blocked by the accumulated debris and corresponding extra trash rack losses, and check several scenarios of trash rack cleaning. The weaknesses of the existing cleaning strategy were shown clearly, and the possibility was presented how to reduce the total number of trash rack cleanings, i.e. plant operational cost, while keeping the energy losses, i.e. amount of produced electricity, unchanged. 1 METHODS A trash rack obstructs the free flow and produces unwanted energy-losses, generated partially by the trash rack bars and partially by the accumulated debris. When clean, the trash rack losses are the smallest, and they increase with the amount of trash collected on the rack. Thus, it is possible to distinguish clearly between the rack's bars and debris caused losses when the clean trash rack losses are known. These can be measured or predicted by one of the empirical equations. 1.1 Trash Rack Losses' Prediction Energy losses caused by the trash rack are categorised commonly as a head-loss, and may be calculated from: 2 Ah = , (1) 2 g Evaluation of Clogged Hydropower Plant Trash Rack Losses 143 Strajniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 where £ is the so-called trash rack head-loss coefficient which has to be determined experimentally or numerically. Fig. 1. Trash rack Kirschmer [5] was one of the first who investigated trash rack losses experimentally. He studied different bar thicknesses s and shapes (having a different bar shape coefficient K), assembled with different clear spacing b and inclined under several angles 6 (see Fig. 1 for the details). He proposed the following universal formula: £= K\b 14 sin0, (2) which, till today, was improved further by many researchers. Levin [6] introduced blockage factor p (Eq. (3)), defined as the ratio of area blocked by the vertical bars Ar>s and the horizontal spacing and supporting elements Aah and the total area of the trash rack field Arf in order to account for the influence of the transversal elements: P = - An? + A a A (3) while Meusburger [7] took into account clogging due to the accumulated debris, which further reduced the flow area by Aad: P = ■ An? + A.tt + A. A (4) He carried out an extensive set of model tests at the ETH Zurich, Switzerland. Different trash rack designs (bar spacing) were investigated and 16 different types of clogging considered. He proposed a rather simple equation: £= K 1 - p sin0, (5) which correlated well with the experimental data for almost any type of clogging and satisfied the boundary condition £ ^ ® when p ^ 1. There are several other formulae available for estimating head loss through trash racks. Josiah et. al [8] focused on a circular bar trash rack commonly used in Sri Lanka. Osborne [9] developed a model for rectangular bars, with simplified blockage factor p = 6/(6+s) neglecting supporting elements, which was then improved further by Clark et. al [10] by considering bar shape. The latter influences head losses significantly, which may be reduced by more than 50 % when streamlined profile cross-section bars are used instead of rectangular bars according to Tsikata et. al [14]. Raynal et al. in [11], [12] and [13] investigated fish-friendly trash racks with inclined [11], angled [12] and streamwise oriented [13] trash racks having significantly reduced clear spacing 6. They made a comprehensive comparison of different equations, and stated that the head losses calculated with the original Kirschmer equation (Eq. (2)) were too low in all cases, while the equations proposed in [7], [9], and [10] produced good results for larger spaces between bars (6/s > 2), which are characteristic for the most of already installed trash racks today. According to this, the Eq. (5) proposed in [7] and already adapted to account for the influence of debris accumulated on the trash rack, was applied in our study. 1.2 Trash Rack Losses' Measurement Trash rack losses may be measured accurately in a laboratory, where so-called model-tests are performed, and head-loss is calculated from the energy conservation equation: v2 v2 h + — = h + — + Ah + Ah,, 2g 2 2g ' (6) where hi and /?2 are the measured upstream and downstream water depth, and Ah; is the measured channel head-loss. Both flow velocities vi and V2 are determined from the measured flow rate, taking into account actual flow areas Ai and A2, respectively. Although the tests are performed under laboratory conditions, the common head-loss measurement uncertainty is up to 3 mm, as reported in [11]. Field measurements, reported here, were carried out on HPP Vuhred, which is one of the eight HPPs on the Drava river in Slovenia owned by the company DEM d.o.o. managing the operation of all eight HPPs working in a chain mode. HPP Vuhred is a medium size impoundment type facility with 13 km long reservoir containing 10 million m3 of water, of which 3 144 Hribernik, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 22 % can be used for generating electric power. HPP Vuhred makes use of a 17.4 m available head, and reaches a net capacity of 72 MW at the installed flow rate of 550 m3/s. The dam structure contains three turbine piers placed between four spillways and left and right bank buildings. Vertical Kaplan turbines are built into the turbine piers with generators overhead. Special equipment for controlling individual units is installed in each turbine pier. HPP unit 1 was instrumented with two level transmitters, as shown in Fig. 2, in order to measure flow rate and trash rack loss simultaneously. The entrance of the inflow channel is 12 m high and 14 m wide. It is divided by a supporting pier, thus, two equal trash racks are applied, each at one side of the pier. The entrance of the inflow channel is 12 m high and 14 m wide. It is divided by a supporting pier, thus, two equal trash racks are Fig. 2. Schematic of trash rack head losses' measurement applied, each at one side of the pier. Trash racks were assembled from rectangular bars 12 mm thick (5) and 130 mm deep (t). Steel rods were inserted through the bars at a regular distance of 700 mm. Circular spacers with 30 mm diameter were installed around the rods to obtain 100 mm clear spacing (b) between the bars. The rack is attached to the sides of the flume and supported by three horizontal beams. The trash rack is submersed, and there is no free water surface behind it to apply non-intrusive water level measurement such as the laser scaning reported in [18]. Two horizontal boreholes were, therefore, drilled through the supporting pier, one at each side of the trash racks. Another two inclined borings were made to reach the horizontal boreholes from the upper plateau of the HPP. Submersible level transmitters were inserted into the inclined borings, for water depth measurements h and hi, respectively, which is presented schematically in Fig. 2. Two temperature compensated Hydrobar I sensors with long-term stability less than 0.1 % from the adjusted range, produced by Klay-Instruments [15], were used in our case. A level transmitter measures hydrostatic pressure at the selected depth h0. The venting tube in the centre of the transmitter cable makes reference to the atmospheric pressure, which means that barometric changes do not cause any shift. As reported in [11], the common head-loss measurement uncertainty was up to 3 mm under laboratory conditions. However, under field conditions, i.e. when the measurements are carried out on an actual HPP unit, measurement uncertainty may increase significantly. In order to evaluate it, both type A and type B measurement uncertainty have to be considered. Type A standard measurement uncertainty is defined by the standard deviation of acquired results under similar conditions. In our case, this means under the same flow conditions, which are dictated by the flow rate and debris accumulation. It was, therefore, decided to use one week's data, acquired immediately after the general refit when the trash rack was almost perfectly clean. All the data within 120 m3/s < Q < 130 m3/s interval (the interval having the highest frequency over one year's operation) were used to calculate the standard deviation of head loss. This was 1.77 mm, and represented type A standard measurement uncertainty. The level transmitter accuracy specified by the producer is 0.1 % from the adjusted range. The latter was set to 4 m (approx. 0.4 bar), thus, the type B measurement uncertainty amounted to 2.31 mm. Combined standard measurement uncertainty was, thus, 2.91 mm and the expanded one approximately 5.8 mm at P = 95 %. Evaluation of Clogged Hydropower Plant Trash Rack Losses 143 Strajniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 According to Eq. (6), the channel head-loss has to be deducted from the total measured head loss, in order to obtain the trash rack caused head loss. While under laboratory conditions the trash rack is simply removed and the channel head loss is measured, but this is not possible (allowed) under field tests of an HPP unit. Thus, the channel loss has to be estimated. In [16], the following equation is suggested: 2 g (7) where ^ is the head-loss coefficient which shall be assumed to vary from 0.1 for gradual contractions to 0.5 for abrupt contractions. The former (^ = 0.1) takes place in our case, since the distance between the two hydrostatic pressure readings is only two metres, and the channel cross-section is reduced by less than 5 %. This way, the channel head loss was estimated to 0.7 mm at nominal (maximal) flow rate 185 m3/s, which was much less than the head loss measurement uncertainty. The channel head loss was, therefore, simply neglected. Fig. 3. Measured flow rate and trash rack loss and computed head loss coefficient Instantaneous water flow rate was determined from the turbine operation diagram with the accuracy 1 % from the nominal flow rate. Simultaneously, the signals from the level transmitters were acquired via a computer every 15 minutes, and saved, together with the water flow rate, to the computer's hard disk. It was, therefore, possible to apply Eq. (6) and calculate total head-loss at 15 minutes intervals. Fig. 3 shows the characteristic head-loss and flow rate time-history acquired during a 3 day period. There are three long operation periods and four short periods when the aggregate stood still and flow rate was 0. Significantly high variations in the flow rate were observed during the aggregate's operation and, therefore, the variations of head loss, which changes with the second power of velocity (Eq. (1)), are even higher. Their frequency and amplitude agree well with the flow rate variations, and prove that the sensors were chosen correctly and that the measurements were performed adequately. Also shown, is a computed head loss coefficient which follows from Eq. (1). Its mean value is approximately 0.79, while its fluctuations within ± 0.1 interval (a = 0.04) correspond to measurement uncertainty. 1.3 Clogged Trash Rack Losses Comparison between the theoretical head loss using Kirschmer's development and laboratory tests have found that, for a clean rack, the theory underestimated the head loss by a factor of 1.75 to 2 [4]. This factor, which is increased greatly when the rack begins to become clogged with debris, was found to be as high as 4 with 50 % clogging [4]. As already mentioned in Section 1.2, Meusburger [7] carried out an extensive set of model tests and proposed Eq. (4) to account for the flow area blocked by the accumulated debris Ad Using Eqs. (5) and (4) it is, therefore, possible to calculate blockage factor: P = ■ I K ■ sin0 1 + I K ■ sin0 (8) when the head loss coefficient is known, and then estimate Ad aad - p ■ arf ((rs + aah ). (9) If the Eqs. (8) and (9) are applied for a set of experimental data acquired during a longer time period, one can obtain the debris accumulation variation with time. However, it is important to start the computation at the moment when the trash rack is 100 % clean (Aad = 0 m2), which is only after the general refit, otherwise the history of debris accumulation is not known, and any debris wedged permanently between the bars which cannot be removed by standard rack cleaning, may deteriorate results significantly. Thus, it was decided to skip the results measured before the general refit (October till February) and use the rest of the experimental data acquired within the 8 month operation period between February 21 and October 1, starting immediately after the general refit of the aggregate. During the refit, the trash rack was dismounted for any necessary repairs and cleaning. After the refit, the trash rack operated in a clean state for a period long enough to acquire data on its operation and adopt Eq. (5) for correct energy loss prediction of the clean trash rack (Aad = 144 Hribernik, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 0 m2). One week of data on the flow rate and head losses measured immediately after the refit were used in our case. The blockage area of a clean trash rack (Eq. (4) at Aad = 0 m2) was adjusted step-wise until the measured head loss coefficient and the one predicted by Eq. (5) fit together well (R2 > 0.999), and the correct total blockage area of a clean trash rack was acquired, corresponding only to the bars and horizontal spacing elements (A^ + Aah). This was then applied in Eq. (9) to calculate the blockage area caused by the debris Aad for any measured point. Fig. 4 shows the variation of Aad within the observed period. The results are highly scattered, the first two weeks' standard deviation is approximately 1 m2. Therefore, moving average was applied to represent the data more clearly, as well as an approximation line was added to point out some important features, like periodic trash rake cleaning. For the first two weeks the trash rake was operating clean. It was winter time, when the river, as well as the debris concentration, is low. The first rack clogging is observed in early spring, when melting snow increases the flow in lowland tributaries, and pushes some of the collected debris into the river. Frequent cleaning, every end of the week (see saw like shape of Aad), keeps the amount of the collected debris low. Much more intense rack clogging takes place in mid spring and early summer, when frequent rain and thunderstorms spill debris collected in the surrounding forests during the winter into the river. Moreover, the water level is high during this period, due to melting snow in the Alps, with a lot of drifting debris originating from over-flooded river banks. The debris collection rate is high during this period, and, although cleaned frequently, the rack clogging is substantial, even more, the rack cleaning is aggravated, and some debris remain on the rack after cleaning. The blockage factor p increases up to 45 % and almost never decreases below 30 % during this period, which increases head losses substantially and reduces electricity production. In mid-August, an intervention rack cleaning took place, and, after that, the cleaning was shifted from Friday to Tuesday and the cleaning period was increased to two weeks in September, which is seen clearly in Fig. 4. 2 RESULTS AND DISCUSION 2.1 Debris Caused Head Losses Hribernik [17] and [18] showed how to predict debris caused head losses when a clean trash rack head loss model is known. A simpler approach to account for debris caused head losses is to apply Eqs. (4) and (5) and set the Aad to zero, and deduct this way calculated clean trash rack losses from the measured ones. A result is shown in Fig. 5, where the cumulative one year energy losses and the variation of the flow rate are presented. The losses are flow rate dependent, thus, the increase in cumulative losses is the highest in the autumn high water season and, at the same time, due to the high concentration of drifting trash such as dead algae and leaves, the influence of the collected debris on cumulative losses is the highest too. During one year of operation, debris causes up to 135 MWh of electricity losses, which is 48 % of all losses, and 60 % of all these losses take place in the relatively short three-month-long autumn high water period. Fig. 5 also shows that the aggregate was not operating between January 25 and February 21 when the refit took place. 2.2 Trash Rack Clogging Trash rack clogging is a random process which is not easy to predict. However, as already shown in Section 1.3, an analysis of trash rack clogging is possible if one calculates the blockage area caused by the debris 21-Feb 21-Mar 18-Apr 16-May 13-Jun 11-Jul 08-Aug 05-Sep 03-0ct Fig. 4. Variation of Aad and blockage factor p between February 21 and October 1 Evaluation of Clogged Hydropower Plant Trash Rack Losses 143 Strajniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 Ol-Oct 12-Nov 24-Dec 04-Feb 18-Mar 29-Apr 10-Jun 22-Jul 02-Sep Fig. 5. Energy losses caused by a clean trash rack and by collected debris Aad for any measured point of the HPP aggregate head loss-time history. Moreover, the rate of growth of the area blocked by the accumulated debris dAAr/dt can be estimated simply by differentiating Aad; time history shown in Fig. 4. Numerical differentiation was applied in our case. Since a too small time step may cause high variations due to the head loss measurement uncertainty, the time interval used was one week, which agreed with the trash rack cleaning period; trash rack cleaning took place every Friday. The results are presented in Fig. 6. Again, only the results after the aggregate general refit are shown, when the rack is clean and Aad = 0 m2 apply. In the first 2 weeks after the refit, the debris concentration in the river was low, and the debris accumulation rate equals 0. This is the winter period of low water when the concentration of debris in water is almost zero. In the middle of March, the snow starts to melt in the lowlands and moderately increased water spills the debris from the river bank into the river. Melting snow in nearby hills prolongs this process into April and, at the beginning of May, the highest spike in debris accumulation took place. It is mid-spring, when frequent rain and thunderstorms spill debris collected 27-Dec 07-Mar 16-May 25-Jul 03-0ct 2.0e-03 | '' 1 1 | 1 1 1 1 | 1 1 1 1 | 1 1 1 1 - | l.Se-03 — 1.0e-03 % £ 0,Se-03 0.0e-00 0 10 20 30 40 Week Fig. 6. Week to week change of average rate of growth of Aad in surrounding forests during the winter into the river and its tributaries. During the summer months June and July, the heavy rains become rare, however, the water flow remains high due to snow melting in the Alps and the debris accumulation is still very vivid with the average blockage rate over 5-10-4 m2/min. It then reduces a bit in August and September. 2.3 Debris Accumulation Simulation A simple model was developed which applies the experimentally obtained rate of growth of Aad (area blocked by the accumulated debris) presented in Fig. 6 to predict accumulation of debris between the successive debris removals. It is possible to predict the instantaneous Aad simply by integrating its rate of growth in time, and starting with Aad = Aad,0 (Aad,0 = 0 m2 if all the debris have been removed successfully during cleaning) each time the debris was removed from the trash rack. As presented in Fig. 4, the rake cleaning effectiveness was not always 100 % and, from May on, some of the debris remained on the rake after cleaning. Thus Aad,0 was not simply set to 0, but it was allowed to be adjusted in order to reflect actual cleaning efficiency. A simple example is presented next to explain the application of the Aad growth model in the 9 month period from March till October 2015. Four quasi-constant flow rates were assumed, with three intervals of zero flow (aggregate standstill). The trash rack cleaning interval was set to 1 week and 3 weeks, respectively, and 100 % debris removal was assumed (Aad,0 = 0 m2). Results are presented in Fig. 7. The growth of head losses proportional to the 1.5 power of blockage factor, as well as Aad (see Eqs. (4) and (5)), may be observed in both cases. Head losses drop to 0.017 m when the trash rack is cleaned, and then the growth, with increasing amount of accumulated debris, takes place 144 Hribernik, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 until the next cleaning. The most intense growth of head losses is observed in May and June, when the rate of growth of Aad is the highest (see Fig. 6; weeks 17 to 21), and the maximum head losses of 0.089 m (3 weeks' cleaning interval) are reached. Comparison of the 1 week and 3 weeks interval cleaning strategy shows small differences in head loss within the first two and a half months, while the rest of the time the differences were extremely high, and the average head loss was at least two times higher, with some peaks more than five times higher when the interval between the successive trash rack cleanings was set to three weeks. This suggests that the constant time period cleaning strategy may not be an appropriate solution, while an unnecessarily high number of cleanings is needed to keep the head losses low (see Fig. 7, 1 week cleaning strategy) since, otherwise, the head losses increase enormously (see Fig. 7, 3 weeks cleaning strategy). Therefore, it is better to perform cleaning as circumstances require by determining an optimal upper limit of head loss above which the trash rack should be cleaned. The results are more realistic when the actual flow rate is applied, as shown in Fig. 8. The trash rack cleaning interval was set to 2 weeks and 98 % debris removal efficiency was assumed (Aad,0 = 1.2 m2). Although the Aad growth between rack cleaning intervals is continuous, the head loss is fluctuating highly due to the fluctuating flow rate. Moving average line was, therefore, added, just to improve the representability. Immediately after the rake cleaning, the head loss is between 0.02 m and 0.03 m on average, while a highly blocked rack produces on average, up to 0.1 m head loss with peak value 0.15 m. As we can see, the model is very flexible, and allows very quick prediction of head losses under selected operation conditions, moreover, the actual operating conditions may be reproduced well (Fig. 9), which confirms model accuracy under common operational conditions. R2 for the results shown in Fig. 9 is higher than 0.97, while overall R2 is approximately 0.93. Different rake cleaning strategies were examined by the developed model The influence of debris removing efficiency was analysed in the first approach. As already mentioned, the rack cleaning was 100 % efficient only in the first two months after the general refit, and a substantial amount of the debris remained on the rack after its cleaning in the following months. 14-Feb 14-Mar 11-Apr 09-May 06-Jun 04-Jul 01-Aug 29-Aug 26-Sep Fig. 7. Debris accumulation prediction at 4 quasi-constant flow rates (115, m3/s 125 m3/s, 140 m3/s and 100 m3/s) Fig. 8. Debris accumulation prediction at actual flow rate conditions Evaluation of Clogged Hydropower Plant Trash Rack Losses 143 Strajniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 Fig. 9. Comparison of predicted and measured head losses Week Fig. 10. Trash rack cleaning efficiency The trash rack cleaning efficiency may be defined as: n _ arf (( + aah + aad,o) (10) ARF - ((RS + AAH ) When all debris are removed from the rack, then Ajd>,0 = 0 m2 and the rack cleaning efficiency is 100 %, otherwise it is lower. Fig. 10 shows the rack cleaning efficiency within an 8 month period after the general aggregate refit. The concentration of the debris in the river, as well as the rack clogging, was low to moderate in the first 10 weeks of operation, and the debris removal from the rack was 100 % efficient. The debris removal efficiency reduces with increasing debris concentration and rack clogging. Its average value was approximately 95 %, with eventual drop off below 90 %. In order to evaluate the effect of imperfect rack cleaning, the simulation of head loss was performed with Aad,0 set to 0 for all performed rack cleanings. Comparison shows that the ideal rack cleaning reduces cumulative energy losses within the 8 months' period by 18 %, and proves how important regular and thorough rake cleaning is in order to keep head losses low to moderate. Of course, it is possible to analyse any other cleaning scenario. One of the possibilities is to apply rake cleaning every time the rake blockage factor exceeds a certain value. The clean rake blockage factor is approximately 25 %, thus the simulations were carried out with upper blockage factor values pmax within the range 30 % and 60 %. Rack cleaning efficiency was set to 100 % in all cases. Results are presented in Fig. 11. Both energy losses and number of cleanings necessary to keep the blockage factor below pmax are highly pmax dependent. The energy losses increase while the number of cleanings reduces exponentially with pmax. Thus, almost thirty rake cleanings are necessary when pmax is set to 0.3, where the energy losses are the lowest, while already at pmax = 0.35, the number of cleanings reduces by 50 %. However, the energy losses increase by only 25 %. At pmax = 0.5, the energy losses double, while the number of cleanings reduces to 5. Any further pmax increase does not reduce the number of cleanings significantly, while the energy losses increase enormously. Fig. 11. Trash rack losses and number of cleanings necessary to keep the blockage factor below pmax From the economical point of view, both energy losses and number of cleanings should be kept low, especially when manually operated rakes are used, which is a common practice on medium size HPPs, and is also applied on the observed HPP where a team of two workers operates the rakes. The total cost of their interventions can be transformed into the electric energy equivalent and added to the energy losses for any pmax set. Optimal pmax may then be found simply by searching the minimum sum of energy losses and energy equivalent of the performed rake cleanings. Fig. 12 shows three possible examples where the energy equivalent was set to 10 MWh, 20 MWh and 30 MWh, respectively, per one rake cleaning. In all three cases, the optimal pmax is between 0.4 and 0.5, where higher pmax values correspond to higher cleaning costs, and vice versa. However, the differences in total energy equivalent loss are very small within this interval for any of the three examples. It is advisable, therefore, to operate with lower pmax values, in order to keep the maximum amount of debris accumulated 144 Hribernik, A. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 on the trash rack, as well as the mechanical load implied on the trash rack, smaller, without worsening the rack operation's economy too much. gS 2000 1600 1200 800 400 Cleaning cost equivalent: A 10 MWh/cleaning ♦ 20 MWh/cleaning □ 30 MWh/cleaning • □ y ♦ s □ a aA?§§ 8 0.2 0,3 0,4 0.5 0.6 0,7 0,8 0.9 1.0 Paw H Fig. 12. Optimal trash rack cleaning strategy According to Fig. 11, it is necessary to perform 9 rack cleanings if pmax is set to 0.4, which is an almost 3 times lower value than the number of cleanings actually carried out between February 21 and October 1 (see Figs. 4 and 8). Cumulative energy losses of this period are estimated to 133 MWh, and 37 % of these losses are caused by the debris, and, as already stated, 18 %, i.e. one half of them, can be avoided if the debris removal was always 100 %. It is possible to achieve the same result with only eleven 100 % effective rake cleanings, which should always take place when the blockage factor exceeds 0.375 (pmax = 0.375). Comparison of the proposed and actually carried out rack cleaning processes is presented in Fig. 13. Variation of Aad shows that any rack cleaning before April 18 may be skipped, while, later on, only approximately every second rack cleaning is necessary; however, their timing should be circumstances dependent (pmax = 0.375) and not periodical, as it was in the actual case. This way, the total number of rack cleanings reduces from 25 to 11, while total cumulative energy losses remain the same. By keeping the blockage factor below 0.375, the maximal amount of collected debris, as well as mechanical load of the trash rack, is kept low. As shown in Fig. 13, the maximal rack area blocked by the debris is 10 m2 when pmax is set to 0.375, while it exceeded 16 m2 in the actual case. 3 CONCLUSIONS Trash rack head losses of a 20 MW HPP aggregate were obtained and analysed experimentally. Using a rack clogging model, it was possible to distinguish clearly between the rack structure and collected debris caused losses. The latter contribute 48 % of total head losses, although the trash rack was cleaned frequently. Periodical - once a week cleaning was performed. The most severe rack clogging was observed in the spring and autumn high water periods, when the rake blockage exceeded 50 %. This worsens the aggregate economy and increases mechanical load. An attempt was, therefore, made to improve the rake cleaning strategy. Debris accumulation and the rack clogging model were applied and different cleaning strategies were examined. It was found that commonly applied periodical cleaning did not suit actual rack clogging intensity. More than 50 % of the performed debris removals were not necessary, while a one week cleaning period was too long during high water conditions. A circumstances required cleaning strategy was, therefore, suggested. Maximal allowed blockage factor was set to 0.375, and the cleaning was performed any time this value was exceeded. This way, the total number of rack cleaning was reduced from 25 to 11, while total cumulative energy losses remained the same. Moreover, the maximal amount of collected debris, as well as the imposed mechanical load of trash rack, was 40 % lower. 21-Feb 21-Mar 18-Apr 16-May 13-Jun 11-Jul 08-Aug OS-Sep 03-0ct Fig. 13. Comparison of proposed and actually carried out rack cleaning processes Evaluation of Clogged Hydropower Plant Trash Rack Losses 143 Strajniski vestnik - Journal of Mechanical Engineering 66(2020)2, 142-152 4 ACKNOWLEDGMENTS The research described in this paper was supported financially by the hydropower plant electricity production company DEM Dravske Elektrarne Maribor d.o.o. Its support is greatly appreciated. 5 REFERENCES [1] Markovic-Hribernik, T., Murks, A. (2007). Slovenia's climate policy efforts: CO2 tax and implementation of EU ETS. Climate Policy, vol. 7, no. 2, p. 139-155, DOI:10.1080/14693062.200 7.9685643. [2] Rakic, N.Z., Gordic, D.R., Sustersic, V.M., Josijevic, M.M., Babic, M.J. (2018). Renewable electricity in Western Balkans: Support policies and current state. Thermal Science, vol. 22, no. 6A, p. 2281-2296, DOI:10.2298/TSCI180512169R. [3] Hribernik, A., Fike, M., Markovic-Hribernik, T. (2013). Economical optimization of a trash rack design for a hydropower plant. Proceedings 17th International Research Conference TMT, p. 457-460. [4] Bradley, J., Richards, D., Bahner, C. (2005). Debris Control Structures - Evaluation and Countermeasures. Department of Transportation: Federal Highway Administration, Salem. [5] Kirschmer, O. (1926). Untersuchungen über den Gefällsverlust an Rechen, Mitteilungen des hydraulischen Instituts der TH München, vol. 1., Munich, DOI:10.1515/9783486752953-003. [6] Levin, L. (1968). Formulaire des conduites forcées oléoducs et conduits d'aération [Handbook of Pipes, Pipelines and Ventilation], Dunod, Paris. [7] Meusburger, H. (2002). Energieverluste an Einlaufrechen von Flusskraftwerken, PhD thesis, Bau-Ing., ETH-Zürich, Zürich. [8] Josiah, N.R., Tissera, H.P.S., Pathirana, K.P.P. (2016). An experimental investigation of head loss through trash racks in conveyance systems. Engineer: Journal of the Institution of Engineers, Sri Lanka, vol. 49, no. 1, p. 1-8, DOI:10.4038/ engineer.v49i1.6913. [9] Osborn, J.F. (1968). Rectangular-bar trash rack and baffle head losses. Journal of the Power Division, vol. 94, no. 2, p. 111-123. [10] Clark, S.P., Tsikata, J.M., Haresign, M. (2010). Experimental study of energy loss through submerged trash racks. Journal of Hydraulic Research, vol. 48, no. 1, p. 113-118, DOI:10.1080/00221680903566026. [11] Raynal, S., Courret, D., Chatellier, L., Larinier, M., David, L. (2013). An experimental study on fish-friendly trashracks - Part 1. Inclined trashracks. Journal of Hydraulic Research, vol. 51, no. 1, p. 56-66, DOI:10.1080/00221686.2012.753646. [12] Raynal, S., Chatellier, L., Courret, D., Larinier, M., David, L. (2013). An experimental study on fish-friendly trashracks -Part 2. Angled trashracks. Journal of Hydraulic Research, vol. 51, no. 1, p. 67-75, DOI:10.1080/00221686.2012.753647. [13] Raynal, S., Chatellier, L., Courret, D., Larinier, M. David, L. (2014). Streamwise bars in fish-friendly angled trashracks. Journal of Hydraulic Research, vol. 52, no. 3, p. 426-431, DOI:10.1080/00221686.2013.879540. [14] Tsikata, J.M., Tachie, M.F., Katopodis, C. (2014). Open-channel turbulent flow through bar racks. Journal of Hydraulic Research, vol. 52, no. 5, p. 630-643, DOI:10.1080/00221686 .2014.928805. [15] Böttcher, H., Gabl, R., Aufleger, M. (2019). Experimental hydraulic investigation of angled fish protection systems-Comparison of circular bars and cables. Water, vol. 11, no. 5, art. ID 1056, DOI:10.3390/w11051056. [16] Tomanova, S., Courret, D., Alric, A., De Oliveira, E., Lagarrigue, T., Tetard, S. (2018). Protecting efficiently sea-migrating salmon smolts from entering hydropower plant turbines with inclined or oriented low bar spacing racks. Ecological Engineering, vol. 122, p. 143-152, DOI:10.1016/j. ecoleng.2018.07.034. [17] Szabo-Meszaros, M., Navaratnam, C.U., Aberle, J., Silva, A.T., Forseth, T., Calles, O., Fjeldstad, H.-P., Alfredsen, K. (2018). Experimental hydraulics on fish-friendly trash-racks: an ecological approach. Ecological Engineering, vol. 113, p. 1120, DOI:10.1016/j.ecoleng.2017.12.032. [18] Rak, G., Hočevar, M., Steinman, F. (2018). Construction of water surface topography using LIDAR data. Strojniški vestnik -Journal of Mechanical Engineering, vol. 64, no. 9, p. 555-565, DOI:10.5545/sv-jme.2017.4619. [19] Klay Systems (2017), from http://www.klay-instruments.com/ hydrobar-i.html, accessed on 2017-04-19. [20] Standards/Manuals/Guidelines series for Small Hydropower Development (2013), 2.2 and 2.3 Civil Works- Hydraulic and Structural Design, Alternate Hydro Energy Center Indian Institute of Technology Roorkee. [21] Hribernik, A. (2017). Optimization of a hydropower plant trash rack cleaning frequency. Proceedings 2017 World Congress on Advances in Nano, Bio, Robotics and Energy (ANBRE17). [22] Hribernik, A. (2016). Influence of debris collected on a trash rack on the economics of a hydropower plant. Proceedings 20th International Research Conference TMT, p. 277-280. 144 Hribernik, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2 Vsebina Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 66, (2020), številka 2 Ljubljana, februar 2020 ISSN 0039-2480 Izhaja mesečno Razširjeni povzetki (extended abstracts) Kamil Urbanowicz, Huan-Feng Duan, Anton Bergant: Prehodni kapljevinski tok v plastičnih ceveh SI 11 Guangjian Wang, Li Su, Shuaidong Zou: Modeliranje dinamike neenakomernih kontaktnih obremenitev in analiza napak prenosa reduktorja 2K-V z ekscentričnim vzbujanjem SI 12 Tomasz Kozior, Al Mamun, Marah Trabelsi, Lilia Sabantina, Andrea Ehrmann: Kakovost površin in mehanske lastnosti preizkušancev, natisnjenih s postopkom FDM, po toplotni in kemični obdelavi SI 13 Peilin Zhang, Bo Li, Xuewen Wang, Chaoyang Liu, Wenjie Bi, Haozhou Ma: Obremenitvene lastnosti razsutega premoga v srednjem žlebu in njegov vpliv na nepremične dele ohišja SI 14 Mihaela Milesan, Constantin Cristinel Gîrdu, Liviu Cîrtînâ, Constanta Râdulescu: Matematični model ter študija hrapavosti in trdote delov iz jekla Hardox 400, odrezanih s CO2-laserjem SI 15 Aleš Hribernik: Ocena izgub zamašenih rešetk hidroelektrarne SI 16 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, SI 11 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2019-09-11 Prejeto popravljeno: 2019-11-07 Odobreno za objavo: 2019-12-03 Prehodni kapljevinski tok v plastičnih ceveh Kamil Urbanowicz1,* - Huan-Feng Duan2 - Anton Bergant3,4 1 Tehnološka univerza Zahodne Pomoranske, Fakulteta za strojništvo in mehatroniko, Poljska 2 Politehniška univerza v Hong Kongu, Fakulteta za gradbeništvo in okoljsko inženirstvo, Kitajska 3 Litostroj Power d.o.o., Slovenija 4 Univerza v Ljubljani, Fakulteta za strojništvo, Slovenija Plastične cevi se v vse večji meri uporabljajo v industrijski praksi, na primer v vodovodnih sistemih. Dinamične analize cevnih sistemov narekujejo uporabo teoretičnih modelov, ki zajemajo viskoelastične lastnosti plastičnih (polimernih) materialov cevi. Viskoelastičnost plastičnih cevi močno poveča dušenje amplitud tlaka med prehodnimi pojavi (vodni udar). Temu ni tako v jeklenih cevovodih z elastičnimi lastnostmi materiala cevi. Teoretični model bazira na enačbah neustaljenega stisljivega kapljevinskega toka v plastičnih ceveh, kontinuitetni in gibalni enačbi. Za viskoelastičen material Hookov zakon ne velja. Plastični material se odzove kombinirano: elastično (trenutno kot na primer jeklo) in zaostalo (vpliv viskoznosti). Vpliv viskoelastičnosti materiala cevi je zajet v kontinuitetni enačbi; elastični del v izrazu za hitrost širjenja tlačnih valov v elastični cevi, viskozni del pa v členu, ki upošteva zaostale napetosti v steni cevi. Obravnavani viskozni člen zapišemo kot konvolucijo časovne spremembe tlaka in viskozne utežne funkcije, ki popiše materialne lastnosti polimerne cevi. Na podoben način je obravnavano neustaljeno stensko trenje v gibalni enači. Popišemo ga z ustaljenim členom (Darcy-Weisbachova enačba) in členom, ki zajema vpliv pospeška kapljevine na trenjske izgube ob steni cevi. V tem prispevku je uporabljen Zielkejev model, ki zajame neustaljenost toka v konvoluciji časovne spremembe pretočne hitrosti in trenjske utežne funkcije. Transformacija postavljenih parcialnih diferencialnih enačb hiperboličnega tipa z uporabo metode karakteristik da navadne diferencijalne enačbe, ki jih rešujemo s pomočjo diferenčne numerične metode. V pravokotno mrežo metode karakteristik sta vgrajena nov poenostavljeni model viskoznega člena polimera in Urbanowiczev model približka neustaljenega kapljevinskega trenjskega člena. Vpeljani aproksimativni utežni funkciji zagotovita računsko učinkovito konvolucijo z uvedbo eksponetnih funkcij na rekruzivni način. Postavljeni numerični model, ki zajema vpliva viskoelastičnosti cevi in neustaljenega stenskega trenja, je validiran z odgovarjajočimi meritvami izvedenimi v evropskih raziskovalnih ustanovah v Londonu (Velika Britanija), Cassinu (Italija) in Lyonu (Francija). Vse tri preizkusne postaje so v osnovi sistemi tipa tlačna posoda-plastični (polimerni) cevovod-ventil. Dolžine (L) in premeri cevovodov (D) so naslednji: (1) London: L = 271,7 m; D = 50,6 mm; (2) Cassino: L = 203,3 m; D = 44 mm; (3) Lyon: L = 43,1 m; D = 41,6 mm. Prehodni pojavi so vzbujeni s hitrim zapiranjem ventila v širokem pasu začetnih Reynoldsovih števil. Primerjave izračuna in meritev so upodobljene v brezdimenzijski obliki. Tak način je posebej prikladen za obravnavo tokov v lamniarnem področju. Za obravnavane primere se rezultati izračuna in meritev opazno razlikujejo, ko v teoretičnem modelu uporabimo ustaljeni model kapljevinskega trenja ob steni cevi. Iskaže se, da ta model nezadostno duši tako nizkofrekvenčne kot visokofrekvenčne tlačne spremembe (pulze). To pomanjklivost smo odpravili z vpeljavo konvolucijskega modela neustaljenega stenskega trenja. V prispevku je pokazan tudi vpliv temperature na viskoelastične lastnosti plastičnih cevi. Temperature v preizkuševališčih, kjer pridobimo viskozne lastnosti materiala cevi, se lahko znatno razlikujejo od temperatur na terenu. Ta vpliv je potrebno ovrednotiti pri preizkusih. Novi model je dva- do trikrat računsko bolj učinkovit od modelov iz literature. Za izračun posebej hitrih prehodov v plastičnih ceveh priporočamo uporabo računsko učinkovitih konvolucijskih modelov, ki zajemajo vpliv viskoelastičnosti materiala cevi (polimeri) in neustaljenega kapljevinskega trenja ob steni cevi. Ključne besede: plastične cevi, viskoelastičnost, vodni udar, prehodni tok, metoda karakteristik, neustaljeno stensko trenje *Naslov avtorja za dopisovanje: Tehnološka univerza Zahodne Pomoranske, Fakulteta za strojništvo in mehatroniko Szczecin, Piastow 19, Szczecin 70-310, Poljska, kamil.urbanowicz@zut.edu.pl SI 11 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, SI 11 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2019-09-11 Prejeto popravljeno: 2019-11-07 Odobreno za objavo: 2019-12-03 Modeliranje dinamike neenakomernih kontaktnih obremenitev in analiza napak prenosa reduktorja 2K-V z ekscentričnim vzbujanjem Guangjian Wang - Li Su - Shuaidong Zou Univerza v Chongqingu, Državni laboratorij za mehanske prenose, Kitajska V praksi manjkajo formule za računanje napak prenosa v prisotnosti ekscentričnosti in neenakomernih obremenitev, malo pa je tudi raziskav natančnosti prenosa za enakomerne kontaktne obremenitve. Neenakomerne kontaktne obremenitve lahko povzročijo težave z vzdržljivostjo (obremenitve na enoto površine so večje od pričakovanih) in dinamična pojave (vibracije zaradi spreminjajočih se obremenitev itd.). Zato obstaja potreba po preučitvi vpliva neenakomernih obremenitev na proces prenosa moči z namenom odprave vibracij in napak prenosa. Za določitev vpliva neenakomernih kontaktnih obremenitev na napako prenosa je bil postavljen dinamični model neenakomernih kontaktnih obremenitev po metodi združenih parametrov, nato pa so bile z modelom analizirane dinamične napake prenosa v pogojih različnih kontaktnih obremenitev zaradi napake ekscentričnosti. Dinamični model je bil preverjen z virtualno simulacijo in eksperimentalno. Teoretična statična napaka prenosa z napako ekscentričnosti je bila izpeljana po metodi koraka ubirnice in metodi združenih parametrov. Model neenakomernih kontaktnih obremenitev je bil določen po metodi združenih parametrov in razrešen z numeričnim integriranjem po metodi Runge-Kutta četrtega reda. Rezultati so bili verificirani s primerjavo dinamične napake prenosa z rezultati simulacije v okolju Recurdyn. Zgrajeno je bilo tudi preizkuševališče napak prenosnikov 2K-V za določanje dejanskih napak prenosa 2K-V in verifikacijo simulacij oz. modelov. Za izboljšanje natančnosti in dinamičnih lastnosti zobniških prenosnikov z majhno razliko v številu zob je bila optimizirana začetna faza ekscentričnosti na podlagi minimalne napake prenosa kot ciljne funkcije. Rezultati raziskave so naslednji: (1) Pri zobniških prenosnikih 2K-V z majhno razliko v številu zob lahko pod vplivom napake ekscentričnosti nastopijo različne kontaktne razmere, vključno z neenakomernimi obremenitvami, trganjem ali povratnim stikom. Razpon dejanskih kontaktnih sil je zato zelo velik. (2) Pri napaki ekscentričnosti je celotna napaka prenosa 2K-V z majhno razliko v številu zob v stanju brez obremenitve enaka maksimalni vrednosti napak prenosa dveh zobniških dvojic z majhno razliko v številu zob. Celotna napaka prenosa v stanju brez obremenitve je zato enaka napaki prenosa zobniške dvojice, ki je trenutno v ubiranju. (3) Pri enaki napaki ekscentričnosti in različnih skupinah faze ekscentričnosti je očitna razlika v krivuljah napake prenosa 2K-V. Napako prenosa je v praksi mogoče izboljšati z optimizacijo začetne faze ekscentričnosti. (4) Krivulja napak prenosa ima nelinearne vrhove, ki so povezani s točkami menjave zobniških dvojic v ubiranju. Zato bi bilo treba preučiti stanje ubiranja v teh točkah. Simulacijski model je bil verificiran s primerjavo simulacijskega diagrama FFT in eksperimentalnih rezultatov. Podatek o začetni fazi napake ekscentričnosti preizkuševališča ni na voljo in simulacijskega modela zato ni mogoče dodatno verificirati s primerjavo diagrama napake prenosa in eksperimentalnih rezultatov. Ta študija je lahko v pomoč pri napovedovanju napake prenosa zaradi napake ekscentričnosti in omogoča kompenzacijo razstopa. Podaja tudi teoretična izhodišča za montažo zobniških dvojic 2K-V z majhno razliko v številu zob. Ključne besede: reduktor 2K-V, nelinearna dinamika, napaka prenosa, neenakomerne kontaktne obremenitve, napaka ekscentričnosti, časovno spremenljiv razstop SI 12 *Naslov avtorja za dopisovanje: Univerza v Chongqingu, Državni laboratorij za mehanske prenose, Chongqing, 400044, Kitajska, gjwang@cqu.edu.cn Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, SI 11 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2019-09-11 Prejeto popravljeno: 2019-11-07 Odobreno za objavo: 2019-12-03 Kakovost površin in mehanske lastnosti preizkušancev, natisnjenih s postopkom FDM, po toplotni in kemični obdelavi Tomasz Kozior1 - Al Mamun2 - Marah Trabelsi2,3 - Lilia Sabantina2 - Andrea Ehrmann2,* 1 Tehniška univerza Kielce, Fakulteta za mehatroniko in strojništvo, Poljska 2 Univerza za aplikativne vede v Bielefeldu, Tehniško-matematična fakulteta, Nemčija 3 Ecole Nationale d'Ingénieurs de Sfax, Tunizija Sodobne dodajalne izdelovalne tehnologije omogočajo izdelavo modelov iz različnih materialov, tehnologija ciljnega nalaganja (FDM) pa spada med pionirske in danes tudi najbolj razširjene postopke. Kakovost površinskega sloja pogosto vpliva na procese obrabljanja in tribološke lastnosti izdelka, toda v literaturi je mogoče najti le podatke o nekaterih izbranih površinskih obdelavah. Avtorji so zato preučili stanje površine preizkušancev, natisnjenih s postopkom FDM, po toplotni oz. kemični obdelavi. Po njihovem vedenju gre za prvo raziskavo na temo različno dolgih obdelav materiala PLA s tekočim acetonom in določitve njihovega vpliva na kakovost površin (hrapavost in valovitost) in na izbrane mehanske lastnosti. Preizkušanci so bili izdelani po postopku FDM, nato pa so bili izpostavljeni toplotni (60 °C do 140 °C) ali kemični obdelavi z acetonom (v tekočem ali plinastem stanju). Površina je bila pred obdelavo in po njej preiskana s konfokalno lasersko mikroskopijo in z infrardečo spektroskopijo s Fourierjevo transformacijo (FTIR). Opravljeni so bili tudi preizkusi za določitev vpliva predstavljene obdelave na natezno trdnost preizkušancev. Parametri površinske hrapavosti Rz, Ra in Rq se niso pomembno spremenili po toplotni obdelavi, zaznana je bila le manjša tendenca k nižjim vrednostim Ra in Rq pri višjih temperaturah. Pri parametrih valovitosti je bila ugotovljena jasna razlika med neobdelanimi preizkušanci in preizkušanci, ki so bili obdelani pri temperaturi 80 °C ali višji. Vsi trije parametri za opis valovitosti Wz, Wa in Wq so se zmanjšali na približno polovico izhodiščne vrednosti. Parametri hrapavosti so se rahlo povišali po najdaljši obdelavi s plinastim acetonom. Pri tekočem acetonu se je vrednost Rz signifikantno spremenila tudi pri najkrajšem času obdelave, medtem ko se ostale vrednosti niso pomembno spremenile. Obdelava z acetonom je signifikantno vplivala na zmanjšanje parametra valovitosti Wz, medtem ko sta se preostala parametra valovitosti le nekoliko zmanjšala. Površina torej postane enakomernejša na večjih skalah, sijaj pa se zmanjša, če površina ni naknadno polirana. Rezultati mehanskih preiskav v praktično nobenem primeru niso pokazali jasnega vpliva toplotne in kemične obdelave na natezno trdnost in na podaljšanje preizkušancev. Le pri 60-minutni obdelavi z acetonom je bilo ugotovljeno signifikantno zmanjšanje napetosti pri porušitvi oz. povečanje raztezka pri porušitvi. Meritve FTIR niso bile konsistentne in morebitnih sprememb povprečne molekulske mase zaradi obdelave z acetonom ni bilo mogoče niti potrditi niti ovreči. Potrebne bodo torej dodatne raziskave z bistveno daljšim trajanjem obdelave s topilom. V splošnem je mogoče privzeti, da obdelava preizkušancev s topilom delno uniči molekularno strukturo PLA in tako oslabi interakcije med molekulskimi verigami, zaradi česar se vedenje materiala spremeni iz krhkega v bolj duktilno. V pričujoči študiji je bil prvič obravnavan vpliv tekočega acetona na 3D-natisnjeni material PLA. Raziskan je bil vpliv kemične in toplotne obdelave na kakovost površine tako na makro- kot na mikroravni (valovitost in hrapavost), s tem pa se odpirajo priložnosti za izbiro vedno optimalne obdelave. Signifikanten vpliv dolgotrajne obdelave z acetonom na mehanske lastnosti omogoča tudi prilagajanje želenih mehanskih lastnosti preizkušancev. Ključne besede: 3D-tiskanje, končna obdelava, FDM, kakovost površine, mehanske lastnosti *Naslov avtorja za dopisovanje: Univerza za aplikativne vede v Bielefeldu, Tehniško-matematična fakulteta, Bielefeld, Nemčija, andrea.ehrmann@fh-bielefeld.de SI 13 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, SI 11 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2019-09-11 Prejeto popravljeno: 2019-11-07 Odobreno za objavo: 2019-12-03 Obremenitvene lastnosti razsutega premoga v srednjem žlebu in njegov vpliv na nepremične dele ohišja Peilin Zhang1,2 - Bo Li1,2 - Xuewen Wang1,2,* - Chaoyang Liu1,2 - Wenjie Bi1,2 - Haozhou Ma1,2 1 Državni laboratorij za mehanizirano premogovniško opremo v provinci Shanxi, Kitajska 2 Tehniška univerza v Taiyuanu, Kolidž za strojništvo in tehniko vozil, Kitajska Zanesljivost obratovanja strgalnih transporterjev je odvisna od obremenitvenih lastnosti razsutega premoga ter od interakcij med razsutim premogom in nepremičnimi deli ohišja. V obstoječih raziskavah so obremenitve nepremičnih delov ohišja, ki jih povzroča premog, modelirane kot konstantne centralne obremenitve v določenih točkah, njihove časovno spremenljive komponente, neenakomerna porazdelitev in nove obremenitve zaradi interakcij pa so prezrte. Zato tudi niso bile razkrite nekatere obremenitvene lastnosti premoga in z njimi povezane težave. V članku je predstavljena sklopitev modela strgalnega transporterja po metodi diskretnih elementov in dinamičnega modela več togih teles za namene simulacij. Zanesljivost sklopljenega modela je bila tudi eksperimentalno preverjena. Preučena je bila porazdelitev gradienta celotne tlačne sile v razsutem premogu v srednjem koritu, nasilne interakcije med razsutim premogom in nepremičnimi deli ohišja ter zastajanje gibanja strgal. Iz rezultatov je razviden gradient tlačne sile v razsutem premogu v koritu ter med dvema strgaloma v smeri y; tlačna sila v prvi polovici razsutega premoga znaša približno 19 % tlačne sile v drugi polovici. Tlačna sila razsutega premoga med dvema verigama v smeri x znaša približno 72 % sile med verigo in robom korita. Strgala rinejo razsuti premog po koritu in potisna sila se prenaša naprej po premogu. Premog v zadnjem delu potiska prednji premog in je zato izpostavljen večjim silam. Gibanje verig vpliva na premog med verigama tako, da se zmanjšujejo tlačne sile v njem. Deli premoga, ki obtičijo med strgali, navpičnimi členi verige in srednjimi ploščami, so v nasilni interakciji z nepremičnimi deli ohišja, povprečna tlačna sila stisnjenih kosov premoga pa je 59-krat večja od tlačne sile v koritu. Kosi premoga z izjemno veliko tlačno silo pospešujejo obrabo srednje plošče, zaradi česar se ta v realnih obratovalnih pogojih najmočneje obrablja pri vodilih za verigo. Strgala zaradi vplivov razsutega premoga zlahka obtičijo na stiku dveh korit, napetost v verigah pred strgali je močno povečana, največja napetost v členih verige pa je približno trikrat večja od največje kontaktne sile med strgalom in koritom. Verige se zaradi razmeroma majhne strukturne trdnosti pogosto pretrgajo. Modelirani strgalni transporter je kratek in v predstavljeni raziskavi niso bile upoštevane nekatere mehanske lastnosti. Modeli strgalnega transporterja v prihodnjih študijah bodo morali upoštevati daljše srednje korito ali kako drugače odpraviti pomanjkljivosti zaradi nezadostne dolžine korita. Sklopljeni model je zanesljiva metoda za preučevanje obremenitvenih lastnosti razsutega premoga ter interakcij med premogom in nepremičnimi deli. Rezultati študije bodo lahko uporabni za izboljšave in optimizacijo konstrukcije strgalnih transporterjev. Ključne besede: strgalni transporter, metoda diskretnih elementov, dinamika več teles, gradient porazdelitve tlačne sile, močna interakcija SI 14 *Naslov avtorja za dopisovanje: Tehniška univerza v Taiyuanu, Kolidž za strojništvo in tehniko vozil, Taiyuan, Kitajska, wxuew@163.com Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, SI 11 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2019-09-11 Prejeto popravljeno: 2019-11-07 Odobreno za objavo: 2019-12-03 Matematični model ter študija hrapavosti in trdote delov iz jekla Hardox 400, odrezanih s CO2-laserjem Mihaela Milesan1 - Constantin Cristinel Girdu1 - Liviu Cirtina2 - Constanta Radulescu2 1 Transilvanska univerza v Brasovu, Fakulteta za tehnološki inženiring in industrijski menedžment, Romunija 2 Univerza Constantin Brancusi, Romunija Tehnologija laserskega rezanja jeklenih materialov se je uveljavila v industrijski proizvodnji na področju natančnega rezanja pločevin, profilov in delov. Kakovost odrezanih površin je bistveno večja kot pri ostalih rezalnih tehnologijah (npr. pri plazemskem ali plamenskem rezanju). V tem kontekstu potekajo prizadevanja za izboljševanje laserskih tehnologij na področju učinkovitosti rezanja, stroškov, produktivnosti in donosnosti. Raziskovalni cilji predstavljene študije so bili: določitev matematičnega izraza za hrapavost Ra v algebrajski in diferencialni obliki; opis rezalnih parametrov CO2-laserskega sistema za izrezanje 45 delov iz plošče materiala Hardox 400; uporaba metodologije odzivne površine in 3D-grafikona za določitev hrapavosti Ra na osnovi vhodnih parametrov, moči in hitrosti; določitev formule za izračun hrapavosti v odvisnosti od moči laserja in rezalne hitrosti po linearnem in kvadratnem modelu, statistični izračuni za linearni in kvadratni model; in analiza variabilnosti trdote površine lasersko odrezanih delov v odvisnosti od hrapavosti. Izveden je bil reduciran faktorski poskus (32) z devetimi laserskimi rezanji materiala HARDOX. Začetni faktorski poskus je bil ponovljen štirikrat v enakih pogojih, pri čemer je bilo izdelanih 45 delov. Na osnovi rezultatov raziskav laserskega rezanja kovinskih materialov je bila sprejeta odločitev za izbiro materiala s posebnimi fizikalno-mehanskimi lastnostmi, o katerem niso na voljo podatki v specializirani literaturi, vendar je močno razširjen v strojegradnji in rudarstvu. Izbrani material je bil HARDOX 400 v obliki 10-milimetrske pločevine. Ta debelina je primerna za preučevanje različnih parametrov, ki karakterizirajo kakovost lasersko odrezanih delov: širina reza, koničnost (nagib bokov), hrapavost bokov in trdota površine. Preliminarne raziskave so identificirale mejne vrednosti izbranih parametrov (moč laserja, tlak kisika kot pomožnega plina, rezalna hitrost), pri katerih poteka rezanje. Izbrane so bile naslednje izhodiščne vrednosti: moč 4200 W, tlak rezalnega plina 0,55 bar, rezalna hitrost 1400 mm/min. Variabilnost moči je znašala 100 W, tlaka pomožnega plina 0,10 bar in rezalne hitrosti 200 mm/min. V preliminarnih preizkusih laserskega rezanja je bil uporabljen kvadratni izrez velikosti 40x40 mm, za glavni eksperiment pa je bil izbran del s tremi ravnimi in eno ukrivljeno stranico. Dimenzije so bile izmerjene na zgornji in spodnji strani izrezanih delov in preostanka plošče. Na vseh 45 delih so bile izmerjene trdota in kemična sestava odrezanih površin ter hrapavost Ra. Površine so bile analizirane po metodi vrtinčnih tokov. Hrapavost je bila izračunana po metodi analize variance ANOVA v programu STATISTICA.7. Predstavljeni rezultati eksperimentov z rezanjem materiala HARDOX 400 in formule za določanje hrapavosti so pokazali, da je najobčutljivejši vhodni parameter moč laserja. Vrednost hrapavosti Ra je manjša pri visokih hitrostih in majhnih močeh, zato je za doseganje majhne vrednosti hrapavosti priporočena minimalna moč laserja. Novost študije je določitev vpliva treh vhodnih parametrov (variabilna moč laserja in rezalna hitrost, konstanten tlak) na hrapavost Ra pri laserskem rezanju delov iz jekla HARDOX 400 ter določitev temperature na površini materiala. Rezultati so uporabni tako za operaterje rezalnih strojev kakor tudi za tehnologe, ki si prizadevajo za dobro hrapavost odrezanih delov. Za doseganje hrapavosti po rezanju, manjše od 3 ^m, je tako npr. treba izbrati moč 4100 W in rezalno hitrost 1400 mm/min pri konstantno nizkem tlaku rezalnega plina. Obstaja množica podobnih raziskav na temo laserskega rezanja jekel in zlitin, le malo podatkov pa je na voljo za obdelavo jekla HARDOX z laserskim žarkom. Rezultat predstavljene študije je matematični opis odvisnosti hrapavosti in kakovosti lasersko izrezanih delov od nastavljenih parametrov laserskega sistema. HARDOX je zahteven jekleni material, ki terja skrbno izbiro vhodnih parametrov. V prihodnjih raziskavah bo treba preučiti še odvisnost parametrov rezalnega procesa od debeline materiala HARDOX 400 za doseganje želene hrapavosti in trdote. Naslednja raziskovalna tema je fokusiranje laserskega žarka v korelaciji z vhodnimi parametri za doseganje optimalne površinske hrapavosti. Ključne besede: lasersko rezanje, hrapavost, trdota, reža, CO2-laser, analiza ANOVA *Naslov avtorja za dopisovanje: Transilvanska univerza v Bra§ovu, Fakulteta za tehnološki inženiring in industrijski menedžment, Bd. Eroilor no. 29, Bratov, Romunija, mihaela.milesan@unitbv.ro SI 15 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, SI 11 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2019-09-11 Prejeto popravljeno: 2019-11-07 Odobreno za objavo: 2019-12-03 Ocena izgub zamašenih rešetk hidroelektrarne Aleš Hribernik* Univerza v Mariboru, Fakulteta za strojništvo, Slovenija Voda, ki žene turbine hidroelektrarne, prinaša s seboj različne plavajoče in lebdeče delce, ki jih imenujemo plavine. Kadar njihova velikost preseže dopustno mejno vrednost, lahko povzročijo poškodbe na turbinskih delih. Zato pred turbinski vtok postavimo rešetko, ki preprečuje vstop večjih plavin. Rešetka je sestavljena iz vzdolžnih palic praviloma pravokotnega preseka, ki jih povezujejo prečno postavljene palice opremljene z distančniki. Tako lahko zgradimo rešetko z izbranim prostim pretočnim presekom. Le-ta je bil sprva izbran tako, da je ščitil turbino pred večjimi plavinami, danes pa so njegove dimenzije pogosto izbrane tako, da preprečuje vstop ribam in jih tako varuje pred poškodbami. Rešetka povzroča energijske izgube, ki so posledica vrtinčenja toka. Le-to povzročajo posamezni elementi rešetke, kot tudi plavine, ki ostajajo ujete na rešetki in jih je zato potrebno redno odstranjevati. Zanimivo je, da lahko v strokovni literaturi zasledimo vrsto raziskav o izgubah rešetk, ki pa se omejujejo predvsem na vpliv geometrije rešetke in oblike preseka vzdolžnih palic in se, z eno samo izjemo, ne dotikajo vpliva plavin, ki mašijo rešetko. Še več, razen površnih ocen ne povedo nič o procesu mašenja rešetk in prispevku tako nastalih izgub k skupni izgubam, ki jih rešetka povzroča v dejanskih obratovalnih razmerah. Zato smo se odločili razmere podrobneje raziskati. Uporabili smo rezultate meritev izgub rešetke agregata 1 HE Vuhred na reki Dravi. Le-te so obsegale podatke o trenutni višini gladine pred in za rešetko ter pretoku, ki so bili izmerjeni s korakom 15 minut v obdobju enega leta. Pri obdelavi rezultatov smo izhajali iz enačbe, ki so jo za napoved izgub rešetke razvili v Inštitutu za hidrologijo na ETH Zurich. Le-to smo preoblikovali tako, da smo lahko na podlagi izmerjenih rezultatov zasledovali spreminjanje stopnje neprepustnosti rešetke, ki jo določata geometrija čiste rešetke in površina s plavinami zamašenega dela rešetke. Ko smo odšteli prispevek geometrije rešetke, smo lahko spremljali časovno odvisno mašenje rešetke. Pri tem smo lahko opazovali intenzivnost rasti količine plavin med posameznimi intervali čiščenja in identificirali obdobja najintenzivnejšega mašenja rešetke, ter ocenili prispevek, ki ga k skupnim izgubam rešetke prispeva mašenje s plavinami in znaša slabih 50 % celoletnih izgub, kljub rednemu, vsakotedenskemu čiščenju rešetke. To nas je spodbudilo k temu, da smo raziskali različne strategije čiščenja rešetke. Za to smo potrebovali model mašenja rešetke, ki lahko napove časovno odvisno rast količine plavin na rešetki med izbranimi intervali čiščenja. V model smo vključili izmerjene rezultate o intenzivnosti mašenja rešetke v posameznih intervalih med dejansko izvedenimi postopki čiščenja. Tako smo lahko z dobro zanesljivostjo napovedali enak potek mašenja rešetke in izgub kot smo ga izmerili v opazovanem obdobju. Prav tako pa smo lahko parametre in pogostost čiščenja rešetke spreminjali, ter napovedovali potek mašenja in izgube rešetke za morebitne drugačne strategije čiščenja. Najprej smo analiziral učinkovitost dejansko opravljenih postopkov čiščenja. Ugotovili smo, da vsa čiščenja v obravnavanem obdobju niso bila povsem učinkovita, in da je bila z izjemo začetnega obdobja takoj po remontu agregata 1 njihova učinkovitost v povprečju 95 %, zaradi česar so se skupne izgube v obravnavanem obdobju povečale za 18 %. V nadaljevanju smo raziskali vpliv pogostosti čiščenja na ekonomsko učinkovitost delovanja agregata. Pri tem smo opazovali razliko med zaslužkom zaradi povečane proizvodnje električne energije, ki jo zagotavlja večja pogostost čiščenja rešetke in povečanimi stroški dela, ki jih to zahteva. Ugotovili smo, da je ne glede na višino cene dela, praksa periodičnega, vsakotedenskega čiščenja rešetke ekonomsko neučinkovita in da je mogoče z načinom čiščenja po potrebi, ko opravimo čiščenje le če stopnja neprepustnosti rešetke preseže 37,5 %, doseči enako proizvodnjo električne energije in znižati število postopkov čiščenja rešetke za 60 %. Pri tem pa se zniža tudi nevarnost mehanskih preobremenitev rešetke zaradi povečane zamašenosti. Ključne besede: hidroelektrarna, rešetka, plavine, energijske izgube, učinkovitost in optimalna strategija čiščenja rešetke SI 16 *Naslov avtorja za dopisovanje: Univerza v Mariboru, Fakulteta za strojništvo, Smetanova 17, 2000 Maribor, Slovenija, ales.hribernik@um.si Information for Authors All manuscripts must be in English. Pages should be numbered sequentially. The manuscript should be composed in accordance with the Article Template given above. The maximum length of contributions is 12 pages (approx. 5000 words). Longer contributions will only be accepted if authors provide justification in a cover letter. For full instructions see the Information for Authors section on the journal's website: http://en.sv-jme.eu . SUBMISSION: Submission to SV-JME is made with the implicit understanding that neither the manuscript nor the essence of its content has been published previously either in whole or in part and that it is not being considered for publication elsewhere. All the listed authors should have agreed on the content and the corresponding (submitting) author is responsible for having ensured that this agreement has been reached. The acceptance of an article is based entirely on its scientific merit, as judged by peer review. Scientific articles comprising simulations only will not be accepted for publication; simulations must be accompanied by experimental results carried out to confirm or deny the accuracy of the simulation. Every manuscript submitted to the SV-JME undergoes a peer-review process. The authors are kindly invited to submit the paper through our web site: http://ojs.sv-jme.eu. The Author is able to track the submission through the editorial process - as well as participate in the copyediting and proofreading of submissions accepted for publication - by logging in, and using the username and password provided. SUBMISSION CONTENT: The typical submission material consists of: - A manuscript (A PDF file, with title, all authors with affiliations, abstract, keywords, highlights, inserted figures and tables and references), - Supplementary files: • a manuscript in a WORD file format • a cover letter (please see instructions for composing the cover letter) • a ZIP file containing figures in high resolution in one of the graphical formats (please see instructions for preparing the figure files) • possible appendicies (optional), cover materials, video materials, etc. Incomplete or improperly prepared submissions will be rejected with explanatory comments provided. In this case we will kindly ask the authors to carefully read the Information for Authors and to resubmit their manuscripts taking into consideration our comments. COVER LETTER INSTRUCTIONS: Please add a cover letter stating the following information about the submitted paper: 1. Paper title, list of authors and their affiliations. 2. Type of paper: original scientific paper (1.01), review scientific paper (1.02) or short scientific paper (1.03). 3. A declaration that neither the manuscript nor the essence of its content has been published in whole or in part previously and that it is not being considered for publication elsewhere. 4. State the value of the paper or its practical, theoretical and scientific implications. What is new in the paper with respect to the state-of-the-art in the published papers? Do not repeat the content of your abstract for this purpose. 5. We kindly ask you to suggest at least two reviewers for your paper and give us their names, their full affiliation and contact information, and their scientific research interest. The suggested reviewers should have at least two relevant references (with an impact factor) to the scientific field concerned; they should not be from the same country as the authors and should have no close connection with the authors. FORMAT OF THE MANUSCRIPT: The manuscript should be composed in accordance with the Article Template. The manuscript should be written in the following format: - A Title that adequately describes the content of the manuscript. - A list of Authors and their affiliations. - An Abstract that should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, as well as the methodology employed. It should summarize the results and state the principal conclusions. - 4 to 6 significant key words should follow the abstract to aid indexing. - 4 to 6 highlights; a short collection of bullet points that convey the core findings and provide readers with a quick textual overview of the article. These four to six bullet points should describe the essence of the research (e.g. results or conclusions) and highlight what is distinctive about it. - An Introduction that should provide a review of recent literature and sufficient background information to allow the results of the article to be understood and evaluated. - A Methods section detailing the theoretical or experimental methods used. - An Experimental section that should provide details of the experimental set-up and the methods used to obtain the results. - A Results section that should clearly and concisely present the data, using figures and tables where appropriate. - A Discussion section that should describe the relationships and generalizations shown by the results and discuss the significance of the results, making comparisons with previously published work. (It may be appropriate to combine the Results and Discussion sections into a single section to improve clarity.) - A Conclusions section that should present one or more conclusions drawn from the results and subsequent discussion and should not duplicate the Abstract. - Acknowledgement (optional) of collaboration or preparation assistance may be included. Please note the source of funding for the research. - Nomenclature (optional). Papers with many symbols should have a nomenclature that defines all symbols with units, inserted above the references. If one is used, it must contain all the symbols used in the manuscript and the definitions should not be repeated in the text. In all cases, identify the symbols used if they are not widely recognized in the profession. Define acronyms in the text, not in the nomenclature. - References must be cited consecutively in the text using square brackets [1] and collected together in a reference list at the end of the manuscript. - Appendix(-icies) if any. SPECIAL NOTES Units: The SI system of units for nomenclature, symbols and abbreviations should be followed closely. Symbols for physical quantities in the text should be written in italics (e.g. v, T, n, etc.). Symbols for units that consist of letters should be in plain text (e.g. ms-1, K, min, mm, etc.). Please also see: http://physics.nist.gov/cuu/pdf/sp811.pdf . Abbreviations should be spelt out in full on first appearance followed by the abbreviation in parentheses, e.g. variable time geometry (VTG). The meaning of symbols and units belonging to symbols should be explained in each case or cited in a nomenclature section at the end of the manuscript before the References. Figures (figures, graphs, illustrations digital images, photographs) must be cited in consecutive numerical order in the text and referred to in both the text and the captions as Fig. 1, Fig. 2, etc. Figures should be prepared without borders and on white grounding and should be sent separately in their original formats. If a figure is composed of several parts, please mark each part with a), b), c), etc. and provide an explanation for each part in Figure caption. The caption should be self-explanatory. Letters and numbers should be readable (Arial or Times New Roman, min 6 pt with equal sizes and fonts in all figures). Graphics (submitted as supplementary files) may be exported in resolution good enough for printing (min. 300 dpi) in any common format, e.g. TIFF, BMP or JPG, PDF and should be named Fig1.jpg, Fig2.tif, etc. However, graphs and line drawings should be prepared as vector images, e.g. CDR, AI. Multi-curve graphs should have individual curves marked with a symbol or otherwise provide distinguishing differences using, for example, different thicknesses or dashing. Tables should carry separate titles and must be numbered in consecutive numerical order in the text and referred to in both the text and the captions as Table 1, Table 2, etc. In addition to the physical quantities, such as t (in italics), the units [s] (normal text) should be added in square brackets. Tables should not duplicate data found elsewhere in the manuscript. Tables should be prepared using a table editor and not inserted as a graphic. REFERENCES: A reference list must be included using the following information as a guide. Only cited text references are to be included. Each reference is to be referred to in the text by a number enclosed in a square bracket (i.e. [3] or [2] to [4] for more references; do not combine more than 3 references, explain each). No reference to the author is necessary. References must be numbered and ordered according to where they are first mentioned in the paper, not alphabetically. All references must be complete and accurate. Please add DOI code when available. Examples follow. Journal Papers: Surname 1, Initials, Surname 2, Initials (year). Title. Journal, volume, number, pages, DOI code. [1] Hackenschmidt, R., Alber-Laukant, B., Rieg, F. (2010). Simulating nonlinear materials under centrifugal forces by using intelligent cross-linked simulations. Strojniški vestnih - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, DOI:10.5545/sv-jme.2011.013. Journal titles should not be abbreviated. Note that journal title is set in italics. Books: Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication. [2] Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken. Note that the title of the book is italicized. Chapters in Books: Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages. [3] Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordič, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553576. Proceedings Papers: Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages. [4] Štefanič, N., Martinčevič-Mikič, S., Tošanovič, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427. Standards: Standard-Code (year). Title. Organisation. Place. [5] ISO/DIS 16000-6.2:2002. Indoor Air - Part 6: Determination of Volatile Organic Compounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva. WWW pages: Surname, Initials or Company name. Title, from http://address, date of access. [6] Rockwell Automation. Arena, from http://www.arenasimulation.com, accessed on 200909-07. EXTENDED ABSTRACT: When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters or approx. 600 words). The instruction for composing the extended abstract are published on-line: http://www.sv-jme. eu/information-for-authors/ . COPYRIGHT: Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to transfer copyright to SV-JME when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Transfer Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on http://en.sv-jme.eu/. PUBLICATION FEE: Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 380 EUR (for articles with maximum of 6 pages), 470 EUR (for articles with maximum of 10 pages), plus 50 EUR for each additional page. The additional cost for a color page is 90.00 EUR (only for a journal hard copy; optional upon author's request). These fees do not include tax. Strojniški vestnik -Journal of Mechanical Engineering Aškerčeva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu http://www.sv-jme.eu Contents Papers Kamil Urbanowicz, Huan-Feng Duan, Anton Bergant: Transient Liquid Flow in Plastic Pipes Guangjian Wang, Li Su, Shuaidong Zou: Uneven Load Contact Dynamic Modelling and Transmission Error Analysis of a 2K-V Reducer with Eccentricity Excitation Tomasz Kozior, Al Mamun, Marah Trabelsi, Lilia Sabantina, Andrea Ehrmann: Quality of the Surface Texture and Mechanical Properties of FDM Printed Samples after Thermal and Chemical Treatment Peilin Zhang, Bo Li, Xuewen Wang, Chaoyang Liu, Wenjie Bi, Haozhou Ma: The Loading Characteristics of Bulk Coal in the Middle Trough and Its Influence on Rigid Body Parts Mihaela MiLesan, Constantin CristineL GTrdu, Liviu CTrtTna, Constanta RaduLescu: Mathematical Modelling Study of Hardox400 Steel Parts' Roughness and Hardness, Cut with CO2 Laser Ales Hribernik: Evaluation of Clogged Hydropower Plant Trash Rack Losses 77 91 105 114 127 9770039248001